鄧?yán)^華，譚建平，譚 平，田仲初

(1. 長沙理工大學(xué)土木工程學(xué)院，長沙 410076；2. 廣州大學(xué)工程抗震研究中心，廣州 510405)

隨著計算機軟硬件技術(shù)的發(fā)展，基于非線性有限元開展結(jié)構(gòu)破壞乃至失效倒塌分析已經(jīng)成為可能[1 ? 2]，這其中研究開發(fā)出各種高效精確的非線性單元模型和算法是基礎(chǔ)和前提[3 ? 5]。對于幾何非線性梁單元，許多學(xué)者從考慮應(yīng)變的高階項推導(dǎo)精確的單元切線剛度矩陣、桿端力與變形后單元坐標(biāo)轉(zhuǎn)換矩陣的精確計算等入手，進(jìn)行了卓有成效的研究[4, 6]。在各種幾何非線性有限元列式中，共旋坐標(biāo)(CR)法認(rèn)為結(jié)構(gòu)的幾何非線性主要由結(jié)構(gòu)大位移引起，因此在分析中每一個單元都建立一個共旋坐標(biāo)系(局部坐標(biāo)系)，該坐標(biāo)系隨單元位移發(fā)生平移和轉(zhuǎn)動，該坐標(biāo)系下單元的純變形可以通過扣除單元剛體平移和轉(zhuǎn)動而得到，相應(yīng)可求得單元切線剛度矩陣和抗力，再通過變換的坐標(biāo)轉(zhuǎn)換矩陣將其轉(zhuǎn)換到結(jié)構(gòu)坐標(biāo)系下，此過程中幾何非線性效應(yīng)被自動計入[7 ? 9]。研究表明，基于共旋坐標(biāo)法研究幾何非線性，不僅具有列式簡單、計算準(zhǔn)確，以及適用于桿、梁、平面、殼、實體單元等各類型單元的優(yōu)點[10 ? 14]，而且能將幾何與材料非線性脫藕，即材料非線性在局部坐標(biāo)系內(nèi)考慮，幾何非線性則通過局部坐標(biāo)系與結(jié)構(gòu)坐標(biāo)系之間的轉(zhuǎn)換予以實現(xiàn)[15?16]。對于梁單元而言，除了大位移能引起幾何非線性以外，梁-柱效應(yīng)也是引起幾何非線性的一個重要因素[17]，這其中采用能考慮構(gòu)件軸力與彎矩相互作用的穩(wěn)定函數(shù)[18 ? 19]來考慮梁柱效應(yīng)是一種最精確的方法[20 ? 21]，實際工程分析中常采用的幾何剛度矩陣實質(zhì)是將穩(wěn)定函數(shù)法表達(dá)的單元剛度矩陣用級數(shù)展開，取其一階近似而得到[22]。但應(yīng)指出的是，以往基于穩(wěn)定函數(shù)法的研究采用的有限元列式一般為T.L或U.L 列式，筆者尚未見到將穩(wěn)定函數(shù)法與共旋坐標(biāo)(CR)法結(jié)合起來進(jìn)行幾何非線性平面梁單元的研究。因此，本文基于共旋坐標(biāo)法，在局部坐標(biāo)系下分離單元剛體位移和變形，采用穩(wěn)定函數(shù)考慮梁柱效應(yīng)；再基于結(jié)構(gòu)坐標(biāo)系與局部坐標(biāo)系下節(jié)點力之間及節(jié)點位移之間的總量關(guān)系及微分導(dǎo)出的增量關(guān)系，獲得平面梁單元在結(jié)構(gòu)坐標(biāo)系下的全量平衡方程和切線剛度矩陣；在此基礎(chǔ)上根據(jù)帶鉸梁端彎矩為零的受力特征，導(dǎo)出了能考慮梁端帶鉸的單元切線剛度矩陣。最后通過數(shù)個經(jīng)典算例對本文算法及程序進(jìn)行了較全面和嚴(yán)格的檢驗。

1 平面梁單元切線剛度矩陣

如圖1(a)、圖1(b)所示分別為初始時刻和計算t時刻的平面梁單元，XY為結(jié)構(gòu)坐標(biāo)系，xy為單元的共旋坐標(biāo)系，它具有隨單元變形而轉(zhuǎn)動的特點，在運動中始終以節(jié)點i為原點，以節(jié)點i到j(luò)的連線方向為x軸，y軸則由x軸按逆時鐘方向旋轉(zhuǎn)90°形成。

圖1 變形前后平面梁單元Fig. 1 Plane beam element before and after deformation

初始時刻設(shè)單元在結(jié)構(gòu)坐標(biāo)系下水平及豎向的投影長度為x0和y0，在計算t時刻結(jié)構(gòu)坐標(biāo)系中的位移向量d=[ui viθi uj vjθj]T，在共旋坐標(biāo)系中的位移向量dL=[uLi vLiθiLuLj vLjθLj]T，聯(lián)立圖1(a)、圖1(b)，有：

平面梁單元在共旋坐標(biāo)系中的節(jié)點位移為：

其中：

式中：si=sinαt,co=cosαt。

圖2(a)、圖2(b)所示分別為梁單元及微段的受力平衡情況。

圖2 梁單元及微段Fig. 2 beam element and micro-segment

對式(11)進(jìn)行與單元受拉時同樣的推導(dǎo)，可得到與式(9)相同的表達(dá)式，只是各參數(shù)s、c、k、ω的計算式變成以下形式：

kT的具體表達(dá)式較繁瑣，限于篇幅此處未列出。

聯(lián)立式(19)和式(20)，可將式(17)寫成：

2 梁端帶鉸時的單元切線剛度矩陣

3 非線性方程組求解和計算流程

非線性方程組的增量法求解主要有荷載增量法和位移增量法[23]。理論上而言，荷載增量法只能計算出荷載-位移曲線的上升段，而位移增量法能計算荷載-位移曲線的下降段，從而能獲得極限荷載值與位移值；但在實際計算中發(fā)現(xiàn)，由于在荷載-位移曲線的頂點附近結(jié)構(gòu)切線剛度接近奇異，即使采用位移增量法，非線性計算也很難收斂。鑒于此，很多文獻(xiàn)都大力推薦弧長增量法，文獻(xiàn)[24]對傳統(tǒng)弧長增量法[25]研究后認(rèn)為，該方法能順利進(jìn)行的前提是必須對結(jié)構(gòu)的特性有深刻的了解并進(jìn)行復(fù)雜的試算，顯然這會給分析工作帶來極大困難和不便；因此，該文獻(xiàn)在傳統(tǒng)弧長增量法的基礎(chǔ)上提出一種如圖3 所示的基于牛頓-拉菲遜法的投影增量法，該方法能有效克服傳統(tǒng)弧長增量法的上述缺點，且該方法收斂速度要稍快，計算量也略小，本文采用該方法進(jìn)行非線性方程組的求解。

本文非線性計算流程如下：

1)根據(jù)式(2)由上一荷載步末單元i、j節(jié)點在結(jié)構(gòu)坐標(biāo)系下的總位移向量d求出共旋坐標(biāo)系下的節(jié)點位移向量dL；

2)由式(13)中最后一行單元軸力N的計算式得到N，再根據(jù)軸力N是受拉還是受壓分別基于式(10)或式(12)計算參數(shù)s、c、k、ω，進(jìn)而根據(jù)式(9)得到單元彎矩Mi和Mj；

3)根據(jù)式(13)可得到單元在共旋坐標(biāo)系下的桿端力向量為f；

4)由式(4)和式(16)分別計算矩陣T和t；

5)由式(19)計算Kn，再由式(22)得到單元在結(jié)構(gòu)坐標(biāo)系下的切線剛度矩陣K；

6)由式(16)得到單元在結(jié)構(gòu)坐標(biāo)系下的節(jié)點抗力向量F；

7)對所有單元重復(fù)步驟1)～步驟6)，基于常規(guī)的“對號入座”原則疊加形成結(jié)構(gòu)總剛矩陣和總抗力矩陣；

8)由總抗力矩陣與總荷載矩陣之差形成不平衡力矩陣作為非線性平衡方程的右端項，解方程得到增量位移矩陣 ?d，與前述總位移向量d疊加形成新的總位移向量d；

9)判斷是否收斂，如是，轉(zhuǎn)至下一個荷載步；如否，進(jìn)行下一輪迭代計算。

圖3 投影增量法示意Fig. 3 Projection increment method

4 算例分析

圖4 Lee’s 框架及荷載 /cm Fig. 4 Lee’s frame geometry and loading

圖5 Lee’s 框架Fig. 5 Lee’s frame

采用3 個經(jīng)典算例對本文算法及程序從正確性、計算效率及精度方面進(jìn)行較全面和嚴(yán)格的檢驗。算例1. 如圖4 所示為Lee’s 框架[26]，材料及幾何參數(shù)如下：線彈性模量E=7060.8 kN/m2，梁和柱具有同樣的截面，且截面面積為6 cm2，慣性矩為2 cm4，對結(jié)構(gòu)進(jìn)行幾何非線性分析，分析中柱均分為8 個單元，荷載作用點左邊梁均分為2 個單元，右邊梁均分為6 個單元，圖5(a)所示為框架的初始構(gòu)型以及不同變形時刻的構(gòu)型，圖5(b)、圖5(c)所示分別為荷載作用點處荷載-水平位移及荷載-豎向位移曲線，可看出與相關(guān)文獻(xiàn)[27]的計算結(jié)果很吻合(該文獻(xiàn)未提供荷載-水平位移曲線)。算例2 . 如圖6 所示懸臂梁，集中荷載P=35 kN作用在自由端，梁跨為L=100 m、彎曲剛度為EI=3.5×104kN·m2、梁端水平位移為u、垂直位移為w、轉(zhuǎn)角為θ，對本算例按以下兩種方法求解：

圖6 端部承受集中荷載的懸臂梁Fig. 6 Cantilever subjected to concentrated load at free end

方法1. 采用本文算法研制的程序進(jìn)行計算；

方法2. 按文獻(xiàn)[11]和文獻(xiàn)[28](文獻(xiàn)[28]就是將文獻(xiàn)[11]的算法從平面梁延伸到空間梁)的方法進(jìn)行計算，即共旋坐標(biāo)系下的桿端力就由節(jié)點位移(去除了剛體位移)與小位移線彈性剛度矩陣相乘得到。

為使計算結(jié)果具有可比性，兩個程序采用的方程組求解方法、數(shù)據(jù)精度及收斂精度標(biāo)準(zhǔn)等都完全相同，方法2 中桿端力也是基于全量計算，詳見文獻(xiàn)[11]及文獻(xiàn)[28]。

兩種方法得到的計算結(jié)果及比較見表1 和表2，解析解見文獻(xiàn)[29]。可以看出，無論是從非線性計算能力、精度還是穩(wěn)定性方面，方法1 均優(yōu)于方法2(方法2 在劃分為1 個單元及4 個荷載步情況下還出現(xiàn)錯誤的收斂解)；還可以看出，由于方法1 與方法2 中桿端力的計算均是基于全量計算，所以在單元數(shù)劃分固定的前提下，計算精度與荷載步數(shù)無關(guān)。

表1 懸臂梁自由端的位移Table 1 Displacement of cantilever at free end

表2 計算誤差 /(%)Table 2 Calculation error

算例3. 如圖7(a)所示，一對拉力2P作用于鉸接方棱形框架的對角點，各項幾何、材料參數(shù)及位移變量詳見圖7，EI為桿件的抗彎剛度。

圖7 對角點受拉力作用時的鉸接方棱形框架Fig. 7 Diamond-shaped frame under a pare of opposite concentrated tensions in opposite angles

為驗證本文所導(dǎo)出的帶鉸梁單元模型，基于荷載及結(jié)構(gòu)對稱的特征，取圖7(b)所示的1/2 結(jié)構(gòu)進(jìn)行計算，按每根桿件分成1、2、3 和4 個單元(中間鉸左、右兩個相鄰單元類型選用本文推導(dǎo)的帶鉸單元)，在荷載參數(shù)PL2/EI=10 時所得的計算結(jié)果如表3 所示，其中解析解見文獻(xiàn)[29]。從表3可看出，計算結(jié)果收斂非常好，如按工程精度誤差不超過5%的要求，則每根桿件劃分成3 個單元就已經(jīng)足夠。

表3 鉸接方棱形框架位移Table 3 Displacement of diamond-shaped frame

5 結(jié)論

本文在局部坐標(biāo)系下采用穩(wěn)定函數(shù)考慮P-δ效應(yīng)，再基于結(jié)構(gòu)坐標(biāo)系與局部坐標(biāo)系下節(jié)點力之間及節(jié)點位移之間的總量關(guān)系及微分導(dǎo)出的增量關(guān)系，獲得平面梁單元在結(jié)構(gòu)坐標(biāo)系下的切線剛度矩陣，同時還獲得單元桿端力的全量算法；在此基礎(chǔ)上，根據(jù)帶鉸梁端彎矩為零的受力特征，導(dǎo)出了能考慮梁端帶鉸的單元切線剛度矩陣表達(dá)式。多個算例表明本文算法是正確的，相對于已有文獻(xiàn)算法，本文算法在非線性計算能力、精度以及穩(wěn)定性方面均占有一定優(yōu)勢，可用于平面桿系結(jié)構(gòu)的幾何非線性分析。