改進共旋坐標法的Timoshenko 梁單元非線性分析

李東升，高嚴培，郭 鑫

(1. 汕頭大學土木與環境工程系，廣東省結構安全與監測工程技術研究中心，汕頭 515063；2. 大連理工大學海岸及近海工程國家重點實驗室，大連 116024)

隨著科學技術的不斷進步，現代工程結構正在向大型化、復雜化和智能化方向發展，如大跨度橋梁、高層建筑和大型風力機等。大型化將導致結構更加趨于柔性，基于傳統小變形假設的有限元分析方法有時已不能滿足計算需求，因此對結構進行幾何非線性分析成為近年研究的熱點問題，也是難點問題。關于幾何非線性的有限元分析方法，最常用的是以結構變形前為參考建立平衡方程的完全Lagrange 列式法(TL 列式)和以結構變形后為參考建立平衡方程的更新Lagrange 列式法(UL 列式)[1]。其中，TL 列式法在求解大轉動問題時，會因參考系的轉動引起不可忽略的誤差，甚至出現錯誤的解，而UL 列式法雖然可以求解大轉角問題，但非線性應變關系的計算非常復雜，導致計算量也隨之偏大[2]。近些年發展起來的共旋坐標法是對上述兩種方法的一種改進，它具有建模簡單，意義明確，并且能借用現有的有限元理論，求解效率和精度相對較高[3]等特點，成為幾何非線性研究的熱點。

共旋坐標法的概念最早由WEMPNER 等[4]提出，即通過剔除結構的剛體位移來得到單元節點的真實位移。CRISFIELD[5]使用共旋坐標法求解各種單元的幾何非線性，提出了一致的單元平衡方程計算方法。AREIAS 等[6]將共旋坐標法用于分析材料失效破壞及斷裂，證明共旋坐標法可以對更復雜的問題進行研究。ST?BLEIN 等[7]基于共旋坐標法，提出了各向異性的Timoshenko 梁單元用于梁單元的幾何非線性分析。鄧繼華等[8]基于共旋坐標法和穩定函數提出了一種幾何非線性平面梁單元，楊浩文等[9]將能量守恒和共旋坐標法相結合對二維Timoshenko 梁單元進行動力學分析，杜柯等[10]利用共旋坐標法模擬框架結構的連續倒塌，ZIYUN 等[11]簡化并開發了一種新的共旋分析框架，來避免非線性投影矩陣的復雜推導。在共旋坐標法中，材料非線性和結構純變形等位于局部坐標系內，因此局部坐標系的選取對于精度的計算有顯著的影響。通常梁單元局部坐標系的選擇有以下幾種方法：RANKIN 法[12]、BATTINI法[13]、CRISFIELD 法[14]和HSIAO 法[15]。前三種方法都是通過左右兩端總轉角值線性插值構造剛體旋轉矩陣，當變形較小時，計算結果較為精確；而當變形較大時，會產生較大的虛假應變。對于HSIAO 法，可以得到比前面幾種方法精確的剛體旋轉矩陣，但對于大應變和大轉動問題的誤差仍然較大，并且計算更加復雜。同時Timoshenko梁單元的剪切自鎖現象也會使得部分變形被放大，從而導致虛假的剛體轉動被放大，誤差較大。

為了解決上述方法中存在計算精度和效率較低的問題，需對現有的共旋坐標法進行改進。因此，本文在共旋坐標法的框架上，基于改進的空間Timoshenko 梁單元[16]，使用逆向運動方法，剔除結構剛體運動部分，再結合梁單元形函數和截面剛度矩陣得到整體坐標系下的單元切線剛度矩陣。最后，通過多種類型的算例對本文提出的方法和程序進行了較全面和嚴格的檢驗。

1 空間梁單元

使用共旋坐標法求解非線性問題時，首先需要求解單元的線剛度矩陣，分析時通常采用插值形函數法構造梁單元，但這些插值函數大多為位移的近似方程，計算往往存在截斷誤差，導致精度較低。本文采用一種改進的空間Timoshenko 梁單元[16]求解單元線剛度矩陣，具體推導過程如下。

1.1 梁單元橫向位移的形函數

對于空間Timoshenko 梁，截面橫向位移包含彎曲位移和剪切位移兩部分，可以表示為：

式中：v、w分別為梁兩個方向總的橫向位移；下標b為彎曲變形產生的位移；下標s為剪切變形產生的位移。

對式(1)求一階導，可得：

式中：E為彈性模量；Iy為截面關于y軸的慣性矩；G為材料剪切模量；A為截面面積；kz為z方向的截面不均勻系數。

將式(3)對x求一階導，化簡之后可得：

由式(5)可知，滿足彎曲位移wb的解析解為三次多項式，同理可得到滿足vb的解析解也為三次多項式。

1.2 梁單元剛度矩陣

與傳統Timoshenko 梁單元一樣，梁軸方向的位移u及轉角 θx的形函數假設為線性插值，而彎曲位移的插值函數為三次多項式，表達式如下：

通常，Timoshenko 梁單元的應變矢量可以表示為：

根據空間梁單元在x=0,L上的邊界條件，可以得到形函數的系數與節點位移之間的關系，其矩陣表達式為：

式中：E(x)為形函數系數向量的系數矩陣；d為節點位移，其表達式為：

把式(13)代入式(10)和式(11)，可得：

最后通過對梁長L的數值積分，得到Timoshenko梁單元的單元剛度矩陣Ke的表達式為：

式中：K11=EA,K22=kyGA，其余參數具體物理意義在HODGES 等[17]和GIAVOTTO 等[18]論文中做出了具體解釋。

2 改進共旋坐標法

共旋坐標法的獨特之處在于從總體位移中提取彈性變形位移來預先定義投影關系。將梁單元運動從初始狀態到最終變形狀態分解為剛體運動部分和純變形部分，其中剛體運動部分包括局部參考坐標系的剛體平移和旋轉。因此共旋坐標法的核心問題就是處理不同坐標之間的轉換，從而建立純變形與整體變形之間關系。

對于經典共旋坐標法，一般通過線性插值構造剛體旋轉矩陣[19]，當變形較小時，計算結果較為精確，但考慮大應變和大轉動問題時則不可忽略計算產生的虛假剛體轉動，本文通過使用逆向運動改進共旋坐標法從而減小這一部分所產生的誤差。

2.1 空間梁單元的參考坐標系定義與轉換

首先定義梁單元相關參考坐標系。

將梁單元的逆向運動分解為兩步：一是剛體旋轉前后梁單元軸線的轉動；二是剛體旋轉前后繞著梁軸線的定軸轉動。首先求得剛體旋轉前后梁的軸向坐標：

2.2 局部坐標系與整體坐標系之間位移向量的轉化

定義梁單元整體位移向量為Pgg，去除剛體變形后局部坐標系下位移向量為Pl。通過2.1 節中介紹的旋轉框架，從總位移Pgg中剔除剛體位移來計算局部位移Pl。通過兩者的轉換關系計算局部坐標系下的內力向量fl和切線剛度矩陣Kl，而內力向量在整體坐標系Fg中的表達式，可以通過平衡整體與局部系統中的內部虛功來得到：

式中，r1、r2、r3和 δr1都容易求得，對于 δr3，由式(36)可知與 δq有關。對于共旋法，其輔助向量的定義假設了小的純變形及其變化且認為輔助向量和r1、r2在同一個平面內，但在求解過程中使用了簡單的線性插值函數來求得q。在使用逆向運動求解輔助向量后，輔助向量和r1、r2在同一個平面，求得的剛體轉換矩陣更加精確，并且這里不假設小的純變形，僅僅假設純變形的增量較小，從而可以求得輔助向量的變化：e

利用求得的切線剛度矩陣計算整體力向量的差值，通過迭代使結果逐漸收斂于精確解，計算流程圖如圖3 所示。

3 數值算例

為檢驗本文所提方法的合理性，本文使用MATLAB 軟件，編制相應的有限元程序。通過算例1 驗證該方法計算平面梁單元所得結果與解析解相比足夠接近；算例2 驗證該方法可以用于計算不考慮耦合效應的空間梁單元非線性分析；算例3 驗證該方法可以用于計算空間耦合梁單元非線性分析；算例4 驗證該方法可以用于空間預彎梁單元非線性分析。

3.1 算例1

如圖4 所示，懸臂梁在自由端受到橫向集中力P的作用，將梁劃分成10 個單元，迭代容許誤差設定為 10-6，表1 給出了本文的計算豎直位移w、水平位移u及 轉角 θ0的結果與文獻[21]的解析解對比。

表1 懸臂梁承受集中荷載時的大撓度變形Table 1 Cantilever's large deformation under concentrated load

從表1 結果對比可知，采用改進共旋坐標法計算所得結果與解析解的結果相比，誤差基本都小于0.5%，能較為精確的計算懸臂梁在集中荷載作用下的大撓度變形。驗證了該方法用于計算平面梁單元非線性分析的正確性。

3.2 算例2

如圖5 所示的空間懸臂梁，彎矩M2作用在自由端，大小為：

式中：參數k在0 到2 之間取值。將懸臂梁劃分為10 個單元，梁端水平位移為u，豎向位移為v。所得不同彎矩大小作用下圖5 所示懸臂梁的變形結果示于圖6；表2 給出了本文算法計算的梁端水平位移和豎向位移結果，與解析解[22-24]和文獻[25]中HAWC2 方法計算結果進行對比。

從表2 可知，對于懸臂梁承受彎矩作用發生大變形，當k取值從0.4 變化到2.0(即彎矩逐漸變大)，改進共旋坐標法的水平位移計算誤差均在0.20%之內，豎直位移計算誤差均在1.10%之內，而HAWC2-30 水平位移計算誤差多數超過了1.0%，豎直位移計算誤差最大達到了22.7%，HAWC2-50 水平位移計算誤差多數超過0.4%，豎直位移計算誤差最大達到了9.7%，說明本文方法與HAWC2 軟件計算結果相比，尤其是大變形情況的位移轉角計算，本文方法要更加精確，驗證了該方法適用于計算不考慮耦合效應的空間梁單元大轉動變形。

表2 懸臂梁承受彎矩作用的大變形Table 2 Cantilever’s large deformation under bending moment

3.3 算例3

考慮耦合效應的懸臂梁結構如圖7 所示，在自由端作用一集中荷載F3=150 N，其耦合的截面剛度矩陣參數詳見文獻[23]。同樣將梁劃分為10 個單元，圖8 給出了沿梁軸的每個節點的位移和轉角所得計算結果。表3 列出了自由端3 個方向的位移和轉角與參考文獻[7]經典共旋坐標法和HAWC2[25]方法計算結果的對比。

從圖8 可知，當耦合梁單元作用一面內荷載時，本文方法能計算出其他方向因耦合作用產生的位移和旋轉。從表3 可知，對于考慮耦合效應的空間Timoshenko 梁單元，改進共旋坐標法的計算結果與經典共旋坐標法和HAWC2 軟件計算結果相比，誤差均有所減小。尤其是對于經典共旋坐標法在計算時因虛假剛體轉動產生的較大誤差(如本算例中的y方向位移和z方向轉角)，使用改進共旋坐標法后這一部分的誤差分別減小了2.4%和4.5%，與HAWC2 計算結果相比除z軸誤差擴大了0.36%外，其他位移和轉角結果均較好，尤其是HAWC2 方法計算誤差較大的y方向位移和x方向轉角，誤差分別減小了約3.6%和2.3%，驗證了該方法適用于計算考慮耦合效應的空間Timoshenko 梁單元的變形。

表3 自由端位移和轉角參數比較Table 3 Comparison of tip displacements and rotations

3.4 算例4

如圖9 所示，一個 45°的懸臂圓弧梁，圓弧半徑R=100 m，在自由端受到豎向集中荷載F=600 N的作用，截面剛度矩陣詳細參數見文獻[26]所示，表4 列出了本文計算的梁自由端三個方向位移x、y、z計算結果與解析解[26]、經典共旋坐標法的計算結果[7]以及商用軟件HAWC2[25]的計算結果的對比分析。

從表4 可知，對于預彎的空間Timoshenko梁，改進共旋坐標法的計算結果與經典共旋坐標法相比，三個方向位移的誤差分別減小了0.2%、0.9%、0.1%，與HAWC2 計算誤差較小的結果相比本文方法也能較為精確的求得，而HAWC2 計算誤差較大的x和y方向位移，本文方法誤差分別減小了1.9%和1.8%，驗證了該方法適用于計算預彎式懸臂梁的大變形，也為后續研究預彎風力機葉片結構提供了研究基礎。

表4 自由端受力下曲梁端部位移比較Table 4 Comparison of the curved beam tip displacements under a force applied at the free end

4 結論

本文基于共旋坐標法和Timoshenko 梁理論，使用新型空間Timoshenko 梁單元計算單元線剛度，使用逆向運動方法得到局部坐標系下的剛體轉換矩陣，進而求得整體坐標系下的單元切線剛度。通過四個典型算例的比較，得到如下結論：

本文提出的改進共旋坐標法適用于不考慮耦合效應、考慮耦合效應及預彎等多種情況下的Timoshenko 梁單元非線性分析，且計算精度與經典共旋理論和HAWC2 方法相比均有所提高。