【文章內(nèi)容簡介】 元 等效結點恢復力的計算式為： ? ? ? ? ? ? ? ?? ?11nn Te ijijF H H B t J???? ?? 其中， ??B 為應變矩陣， J 為雅克比行列式。 ??B 和 J 的值與高斯積分點相對應。 針對不同的單元，一般通過經(jīng)驗和試算來確定最佳的積分點數(shù)。對于四邊形單元而言， 22? 的高斯積分點就可以較好地滿足精度要求 。表 中所示的是各個高斯積分點的坐標和加權系數(shù)，圖 所示的是 22? 高斯積分點的分布情況。 表 高斯積分中的 i? 和 iH Table i? and iH in Gauss Integration 點數(shù)n 積分點坐標 i? 加權系數(shù) iH 點數(shù)n 積分點坐標 i? 加權系數(shù) iH 1 1? ＝ 0 1H ＝ 2 4 1? ＝ 2? ＝ 3? ＝ 2?? 4? ＝ 1?? 1H ＝ 2H ＝ 3H ＝ 2H 4H ＝ 1H 2 1? ＝ 1? ＝ 2?? 1H ＝ 1 2H ＝ 1 3 1? ＝ 2? ＝ 0 3? ＝ 1?? 1H ＝ 2H ＝ 3H ＝ 1H 5 1? ＝ 2? ＝ 3? ＝ 0 4? ＝ 2?? 5? ＝ 1?? 1H ＝ 2H ＝ 3H ＝ 4H ＝ 2H 5H ＝ 1H ζη 圖 高斯點的位置 Fig The location of Gauss Points NewtonRaphson 迭代法 基于切線剛度的 NewtonRapson 迭代法如圖 。 圖 基于切線剛度的 NewtonRaphson 算法 NewtonRaphson algorithm based on tangent stiffness NewtonRaphson 是一種變剛度 的迭代法。 首先 ， 取 單元的 初始剛度 為 ? ?0K ，根據(jù)所施加的外荷載等效結點力 ??P 可以求出 位移的第一次近似值 為 ： ? ? ? ? ? ?110KP? ?? () 由 得到的 初始位移可 以計算 單元 的 應變 值，再結合本構模型求到單元應力，進而得到相應的結點恢復力 ??1P 。 然后，用外荷載等效結點力 ??P 和結構恢復力 ??1P 求出殘余力 ? ?1P? ，并連同 相應 于初始位移 ??1? 的即時切線 剛度 ? ?1K ， 求出 位移增量? ?2?? ，即 ? ? ? ? ? ?11P P P? ? ? () ? ? ? ? ? ?12 1 1KP? ?? ? ? () 得到 位移的第二次近似值 為： ? ? ? ? ? ?2 1 2? ? ?? ? ? () 以此類推，第 k 步的計算過程為： P 0 1?0 2?0 3?0 1P?0 1P 2P?0 2P 1??0 2??0 3??0 ?0 ? ? ? ? ? ?kkP P P? ? ? () ? ? ? ? ? ?11k k kKP? ??? ? ? () ? ? ? ? ? ?11k k k? ? ???? ? ? () 當 ? ?1k?? 與 ??k? 的誤差達到所規(guī)定的限值 ，或 者 當 ? ?kP? 足夠小的時候，迭代結束 。 弧長控制法 弧長控制法的基本原理是，通過 在結構分析的一般方程 ? ?iiF x f? ?? 中 添入一個 新的 變量 ? ， 來實現(xiàn)對整個計算過程的 自動控制 。 加入變量 ? 后， 原方程組 將 變?yōu)?： ? ?iiF x f?? ? ? （ ） 上述 方程組 一 共有 n+1 個未知量： nxR?? 、 ? ， 但是 方程 卻 只有 n 個 。所以，需要補充如下條件來完成方程組的求解： 2222f x S? ? ? ? ? （ ） S 稱 為 “ 弧長 ”，是 附加的控制變量 。這時，方程組變?yōu)?? ? ? ?11nnFR? ? ??? ?；￠L法的迭代過程如圖 所示。 圖 弧長法迭代過程示意圖 Fig The Arclength method 收斂準則 基于切線剛度的 NewtonRaphson迭代法的收斂準則有三種：力準則、位移準則以及能量準則。選取合適的收斂準則，有利于提高迭代法的計算效率，是增量求解策略中一個重要的組成部分。 當?shù)Y束時，結構應當處于一個穩(wěn)定的力平衡狀態(tài)，力準則就是要求結構的恢復力和外力 之差等于零。但是在數(shù)值計算中，不可能實現(xiàn)不平衡力完全為零，這就需要設定一個允許值。力準則的表達式為： ? ? ? ? ? ? ? ?1 1 1 1 22im m m mFR F R F?? ? ? ?? ? ? 其中， F? 就是給殘余力設置的一個允許值。 在進行位移分析時，計算位移應該與真實值接近。位移準則的表達式為： ? ? ? ? ? ? ? ? ? ?1122iimmDU U U? ??? ? ? 其中，2是 矢量的歐幾里德范數(shù)， ? ???iU? 是第 i 步迭代的位移增量， D? 為 設置的允許值。 能量準則是用來度量力和位移同它們的平衡值之間的接近程度的，其表達式為： ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ?? ?1 1 1i T i i Tm m m mEU R F U R F?? ? ?? ? ? ? ? （ ) 其中，不等式 的 左邊 表示殘余力 在位移增量上所做的功 ，不等式的 右邊是 殘余力所 做 的 功的初始值 ， E? 是 為能量設置 的 允許值。 鄧興龍和蔣曉華等編制的 FEAPpv程序采用的是能 量收斂準則。 基于 MCFT 理論的 割線剛度迭代法 F. J. Vecchio 于 1990 年在 ACI 發(fā)表的論文中采用了基于 MCFT 理論的割線剛度迭代法來進行鋼筋混凝土的有限元分析。通過比較割線剛度迭代法與基于切線剛度的 NewtonRaphson 法，可以得出以下幾點結論： 1. 在結構的應力 應變曲線進入到下降段后，對于切線剛度法，剛度會出現(xiàn)負值，從而給計算處理帶來難度。而割線剛度迭代法則沒有這種負剛度問題，但是割線剛度迭代法會在卸載曲線穿越坐標軸時出現(xiàn)數(shù)值問題，這是由于此時割線剛度極小造成的。 2. 以往采用切線剛度法的程序，在更新應變時，對于一個單元采用了 4 個高斯積分點，每次迭代得到了一個新的結構位移，就需要在 4 個積分點上分別求出它們的應變，再以此計算出 4 個新的割線剛度，用以判斷收斂與否。而 F. J. Vecchio 采用的割線剛度迭代法只涉及了 1 個積分點，即中心積分點。 根據(jù)附錄 1 可知，單元應變和位移的關系為： jpimx i j m pNNNNu u u u ux x x x x? ?????? ? ? ? ?? ? ? ? ? jpimy i j m pNNNNv v v v vy y y y y? ?????? ? ? ? ?? ? ? ? ? j p j pi m i mx y i j m p i j m pN N N NN N N Nuv u u u u v v v vy x y y y y x x x x? ? ? ? ?? ? ? ???? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ?其中： 1 (1 )(1 )4i xyN ab? ? ? 1 (1 )(1 )4j xyN ab? ? ? 1 (1 )(1 )4m xyN ab? ? ? 1 (1 )(1 )4p xyN ab? ? ? 令 2ax? ， 2by? ，則 可以得到割線剛度迭代法中心積分點的應變計算式為： 2j m i px u u u ua? ? ? ?? 2m p i jy v v v vb? ? ? ?? 22m p i j j m i pxy u u u u v v v vba? ? ? ? ? ? ??? 每次循環(huán)迭代出的新位移，僅需要換算成中心積分點上的應變，通過一次計算便可以得出割線剛度用于收斂判斷。相比于切線剛度法，割線剛度迭代法簡化了計算過程。 3. 切線剛度法的過程，需要通過外荷載節(jié)點力和結構恢復力求解出一個殘余力，再以此計算出增量位移來判斷是否滿足收斂條件，若不收斂，便更新結構位移并重新計算結構的切線剛度矩陣。對于割線剛度迭代法，并不需要求解殘余力，而是直接利用外荷載節(jié)點力和結構割線剛度來求解出一個全量位移，之后通過計算出的應變來重新計算結構的割線剛度矩陣。二者 相比起來，割線剛度迭代法不需要啟用一個新的數(shù)組來儲存殘余力值，降低了對程序的存儲空間要求并提高了程序的運算速度。 從上面三點可以看出，切線剛度法和割線剛度迭代法各有特點和不足，在某些方面，割線剛度迭代法還體現(xiàn)出了一定的優(yōu)勢： 4. 割線剛度迭代法相比于切線剛度法的收斂準則更為簡單。 5. 割線剛度迭代法相比于切線剛度法的有限元模型更為簡單。 6. 割線剛度迭代法的穩(wěn)定性較好。 所以對割線剛度迭代法進行研究佐證是十分有必要的。 算法流程 以圖 所示的單元為例， Vecchio 在其論文中所述的割線剛度迭代法的具體計算 過程如下： 圖 單元尺寸 Element dimension 第 1 步：輸入結構和材料的屬性數(shù)據(jù) 第 2 步：計算外荷載節(jié)點力 ??F 第 3 步：計算材料的剛度矩陣 ： 1) 確定初始應變： x? ， y? ， xy? ； 2) 計算主方向應變 1c? ， 2c? ， sx? ， sy? ，以及主拉應變角 c? ； 3) 根據(jù)混凝土的本構關系，確定混凝土主拉和主壓應力 1cf ， 2cf 4) 根據(jù)鋼筋的本構關系，確定鋼筋的主拉和主壓應力 sxf ， syf ； 5) 計算混凝土和鋼筋的割線剛度元素值： 222cccfE ?? 111cccfE ?? 12ccccEEG EE? ? sysysyfE ?? sxsxsxfE ?? 第 4 步： 計算混凝土和鋼筋的材料割線剛度矩陣 1) 計算混凝土主方向上的材料割線剛度矩陣 ? ?39。cD ： ? ? 239。1000000cccEDEG??????? (32) 2) 計算轉換矩陣 ??cT ? ?2222c o s sin c o s sinsin c o s c o s sin2 c o s sin 2 c o s sin c o s sinT? ? ? ?? ? ? ?? ? ? ? ? ??????????? 3) 計算整體坐標系（ XY 坐標）下的混凝土割線剛度矩陣 ? ? ? ? ? ? ? ?39。TccD T D T? (33) 4) 計算整體坐標系 (XY 坐標系 )下的鋼筋割線剛度矩陣： ? ? 00000 0 0x sxs y syEDE????????? (34) 第 5 步：計算組合材料剛度 ??D ： ? ? ? ? ? ?csD D D?? (35) 第 6 步：計算單元剛度矩陣 ??k ： ? ? ? ? ? ? ? ?Tk B D B d x d y d t? ??? (36) 第 7 步 ： 計算單元預應變產(chǎn)生的位移 00,csrr? ? ? ?? ? ? ? 第 8 步：計算單元預應變產(chǎn)生的節(jié)點力 *F???? 第 9 步：計算單元節(jié)點力： 39。*F F F? ? ? ? ? ???? ? ? ? ? ? 第 10 步：計算結構總剛矩陣 ??K ： ? ? ? ?Kk?? (37) 第 11 步：計算結構總剛矩陣的逆矩陣 ? ?1K? ： 第 12 步：計算節(jié)點位移 ? ? ? ? 1 39。r K F? ??? ?? (38) 第 13