基于位移控制新法的結構非線性有限元數值模擬

陳夢成，李嘉鈺

（華東交通大學土木建筑學院,江西 南昌330013）

結構全過程分析是非線性分析理論的基礎，它常常應用于大跨度空間結構的設計之中。 在結構非線性全過程分析中，剛度矩陣的組裝、計算方法及迭代方式的選取，都將對結果產生重大的影響。 近年來，人們一直在對結構非線性全過程分析中的迭代控制算法進行研究，提出了壓縮平衡迭代法、人工彈簧法、位移控制法、弧長法等多種方法。 王志軍等[1]在常規的位移控制法和弧長法中引入了簡單適用的考慮路徑搜尋和單因子加速的方法，能明顯改善收斂情況；鄭宏等[2]利用Sherman-Morrison 定理，對經典的荷載增量控制法程序做了少量的修補，構建了一個有良好數值特性且僅需很少運算量的位移控制法程序；杜柏松等[3]引入廣義位移作為約束方程來確定荷載增量因子，從而控制結構荷載加載的全過程，形成了廣義的位移控制法；許茂等[4]針對負剛度結構，采用了附加正剛度彈簧的方法進行了分析，解決了負剛度切線矩陣在有限元分析中造成的迭代求解不收斂的問題；李鋼等[5]證明了擬力法在計算結構負剛度時的優勢，它是一種高效的非線性分析方法；曲激婷等[6]基于向量式有限元的基本思想，提出了力-位移控制的復合模擬方法，此方法可以很好地求解結構的大變形問題；AL－Aukaily 等[7]采用了直接微分法和向前差分法，同時考慮結構的幾何非線性和材料非線性，對求解結構響應所用位移控制法的相關參數進行了敏感度分析。

利用杜修力等提出的位移控制新方法[8]，編寫相關計算程序，并結合實際算例，對結構的非線性全過程進行有限元分析，驗證位移控制新方法的有效性、正確性和可行性。

1 U.L.列式下增量有限元方程

根據連續性介質力學和虛位移原理[9]，并以時刻結構的狀態作為基準參考構型，應變、應力以及位移的增量關系為[10]

式中：KL為整體坐標下的彈性剛度矩陣；KG為整體坐標下的幾何剛度矩陣，兩者之和為單元切線剛度矩陣；fe為單元抵抗力向量。

2 位移控制新方法

2.1 傳統位移控制法

傳統的位移控制法（displacement control method，DCM）是求解材料、構件及結構包含極限強度后軟化下降段的荷載-位移關系非線性全過程曲線的常用方法之一，它是通過引入給定力作用節點（位移控制點）的位移增量來反求荷載增量和其它未知節點的位移增量。 其一般迭代方程為[11]

式中：[KT]為結構的切線剛度矩陣；Δλ 為荷載因子增量；{P}為參考荷載向量；[R]為不平衡力向量；j 為增量步數；i 為迭代次數。 接下來根據需要控制的位移向量對結構的迭代平衡方程式（22）進行重新排序，即對結構的切線剛度矩陣、荷載向量以及不平衡力向量進行如下分塊

至此，可以計算出每一步迭代中的荷載因子增量Δλ，最終可計算其余點位移分量。

2.2 位移控制新方法

傳統位移控制法的求解過程需將有限元方程的剛度矩陣進行分塊并重新排列，從而得到新的方程和系數矩陣，但新的系數矩陣將不再具有原剛度矩陣的對稱性及帶狀性，會使計算量顯著增加，且計算過程較為復雜，最終帶來諸多計算上的不便。

一般地，有限元增量方程可寫為

3 算例分析

3.1 幾何非線性

選取William 肘式框架[12]進行計算分析，William 肘式框架是眾多學者檢驗其所建立理論有效性及正確性的典型算例，其原因在于該結構具有顯著的幾何非線性特征，同時又具備可以用來與理論值對比的試驗結果，如圖1 所示。 該結構兩端固結，水平方向跨度為657.6 mm， 豎直方向高為9.8 mm，桿件為細長矩形等截面直桿，截面寬19.1 mm，高6.17 mm，彈性模量為71 018.5 MPa，框架頂端作用一豎向集中力荷載。

在分析時， 將該結構離散成20 個Euler-Bernoulli[13]梁單元進行模擬，取控制點的位移增量步步長u=0.1 mm，總位移增量步n=180 即分析結構位移控制點0 至18 mm 的荷載-位移關系曲線，允許位移收斂誤差δ=10-4，在U.L.列式下，使用MATLAB 進行編程計算。

圖2 給出了運用傳統位移控制法、 位移控制新方法得到的數值分析結果與William 試驗結果之間的對比；圖3 給出了運用傳統位移控制法與位移控制新方法得到的結構最終變形曲線。可見，位移控制新方法能跨越荷載-位移關系中的極值點以及結構極限強度后的軟化下降段模擬，且能準確地計算出結構的荷載-位移全過程曲線。兩種方法所得到的結構最終變形曲線高度吻合，證明了位移控制新方法的可行性和正確性。

圖2 William 肘式框架荷載-位移曲線Fig.2 The load-displacement curve for William frame

圖3 William 肘式框架變形曲線Fig.3 The deformation curve for William frame

3.2 材料非線性

一混凝土棱柱體試件[14]尺寸為150 mm×150 mm×300 mm，混凝土等級為C30，考慮材料非線性，對該構件進行非線性有限元分析。

計算分析時，采用三維纖維梁單元，沿長度方向離散成3 個單元， 截面劃分成15×15 的纖維網格，取控制點的位移增量步步長u=0.1 mm，總位移增量步n=24，即分析構件應變為0～0.008 時對應的應力大小，允許位移收斂誤差δ=10-4。混凝土的本構模型根據我國混凝土設計規范[15]，采用單軸受壓應力-應變曲線。

圖4 混凝土柱應力-應變曲線Fig.4 Stress-strain curve of concrete column

將運用傳統位移控制法和位移控制新方法得到的數值分析結果與文獻[14]結果進行對比，如圖4 所示。可見，兩種方法計算出的曲線與文獻中曲線的吻合程度良好，且兩種方法的結果高度吻合，證明了位移控制新方法的準確性和可行性；除此之外，本文在利用MATLAB 進行編程計算時，對兩種方法的計算完成時間進行了測試，結果發現，傳統位移控制法所需時間約為333 s，而位移控制新方法所需時間約為202 s，新方法近節約了40%的時間，優化了計算過程。

4 結論

通過以上理論分析與數值模擬計算，可以得到以下結論：

1） 位移控制法可以模擬出結構的荷載-位移全過程非線性曲線，尤其是可以模擬出曲線含有極值點或軟化下降段的情況；

2） 位移控制新方法不會破壞有限元剛度矩陣的對稱性和帶狀性，且能得到穩定、精確的數值解；

3） 位移控制新方法較傳統位移控制法更為簡便實用，求解所需的計算步驟和計算時間都可以得到一定程度上的優化，節約了計算成本。