改進(jìn)共旋坐標(biāo)法的Timoshenko 梁單元非線性分析

李東升，高嚴(yán)培，郭 鑫

(1. 汕頭大學(xué)土木與環(huán)境工程系，廣東省結(jié)構(gòu)安全與監(jiān)測工程技術(shù)研究中心，汕頭 515063；2. 大連理工大學(xué)海岸及近海工程國家重點(diǎn)實(shí)驗(yàn)室，大連 116024)

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

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

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

1 空間梁單元

使用共旋坐標(biāo)法求解非線性問題時(shí)，首先需要求解單元的線剛度矩陣，分析時(shí)通常采用插值形函數(shù)法構(gòu)造梁單元，但這些插值函數(shù)大多為位移的近似方程，計(jì)算往往存在截?cái)嗾`差，導(dǎo)致精度較低。本文采用一種改進(jìn)的空間Timoshenko 梁單元[16]求解單元線剛度矩陣，具體推導(dǎo)過程如下。

1.1 梁單元橫向位移的形函數(shù)

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

式中：v、w分別為梁兩個(gè)方向總的橫向位移；下標(biāo)b為彎曲變形產(chǎn)生的位移；下標(biāo)s為剪切變形產(chǎn)生的位移。

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

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

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

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

1.2 梁單元?jiǎng)偠染仃?/h3>

與傳統(tǒng)Timoshenko 梁單元一樣，梁軸方向的位移u及轉(zhuǎn)角 θx的形函數(shù)假設(shè)為線性插值，而彎曲位移的插值函數(shù)為三次多項(xiàng)式，表達(dá)式如下：

通常，Timoshenko 梁單元的應(yīng)變矢量可以表示為：

根據(jù)空間梁單元在x=0,L上的邊界條件，可以得到形函數(shù)的系數(shù)與節(jié)點(diǎn)位移之間的關(guān)系，其矩陣表達(dá)式為：

式中：E(x)為形函數(shù)系數(shù)向量的系數(shù)矩陣；d為節(jié)點(diǎn)位移，其表達(dá)式為：

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

最后通過對梁長L的數(shù)值積分，得到Timoshenko梁單元的單元?jiǎng)偠染仃嘖e的表達(dá)式為：

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

2 改進(jìn)共旋坐標(biāo)法

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

對于經(jīng)典共旋坐標(biāo)法，一般通過線性插值構(gòu)造剛體旋轉(zhuǎn)矩陣[19]，當(dāng)變形較小時(shí)，計(jì)算結(jié)果較為精確，但考慮大應(yīng)變和大轉(zhuǎn)動問題時(shí)則不可忽略計(jì)算產(chǎn)生的虛假剛體轉(zhuǎn)動，本文通過使用逆向運(yùn)動改進(jìn)共旋坐標(biāo)法從而減小這一部分所產(chǎn)生的誤差。

2.1 空間梁單元的參考坐標(biāo)系定義與轉(zhuǎn)換

首先定義梁單元相關(guān)參考坐標(biāo)系。

將梁單元的逆向運(yùn)動分解為兩步：一是剛體旋轉(zhuǎn)前后梁單元軸線的轉(zhuǎn)動；二是剛體旋轉(zhuǎn)前后繞著梁軸線的定軸轉(zhuǎn)動。首先求得剛體旋轉(zhuǎn)前后梁的軸向坐標(biāo)：

2.2 局部坐標(biāo)系與整體坐標(biāo)系之間位移向量的轉(zhuǎn)化

定義梁單元整體位移向量為Pgg，去除剛體變形后局部坐標(biāo)系下位移向量為Pl。通過2.1 節(jié)中介紹的旋轉(zhuǎn)框架，從總位移Pgg中剔除剛體位移來計(jì)算局部位移Pl。通過兩者的轉(zhuǎn)換關(guān)系計(jì)算局部坐標(biāo)系下的內(nèi)力向量fl和切線剛度矩陣Kl，而內(nèi)力向量在整體坐標(biāo)系Fg中的表達(dá)式，可以通過平衡整體與局部系統(tǒng)中的內(nèi)部虛功來得到：

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

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

3 數(shù)值算例

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

3.1 算例1

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

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

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

3.2 算例2

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

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

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

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

3.3 算例3

考慮耦合效應(yīng)的懸臂梁結(jié)構(gòu)如圖7 所示，在自由端作用一集中荷載F3=150 N，其耦合的截面剛度矩陣參數(shù)詳見文獻(xiàn)[23]。同樣將梁劃分為10 個(gè)單元，圖8 給出了沿梁軸的每個(gè)節(jié)點(diǎn)的位移和轉(zhuǎn)角所得計(jì)算結(jié)果。表3 列出了自由端3 個(gè)方向的位移和轉(zhuǎn)角與參考文獻(xiàn)[7]經(jīng)典共旋坐標(biāo)法和HAWC2[25]方法計(jì)算結(jié)果的對比。

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

表3 自由端位移和轉(zhuǎn)角參數(shù)比較Table 3 Comparison of tip displacements and rotations

3.4 算例4

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

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

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

4 結(jié)論

本文基于共旋坐標(biāo)法和Timoshenko 梁理論，使用新型空間Timoshenko 梁單元計(jì)算單元線剛度，使用逆向運(yùn)動方法得到局部坐標(biāo)系下的剛體轉(zhuǎn)換矩陣，進(jìn)而求得整體坐標(biāo)系下的單元切線剛度。通過四個(gè)典型算例的比較，得到如下結(jié)論：

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