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

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

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

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

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

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

初始時(shí)刻設(shè)單元在結(jié)構(gòu)坐標(biāo)系下水平及豎向的投影長度為x0和y0，在計(jì)算t時(shí)刻結(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é)點(diǎn)位移為：

其中：

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

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

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

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

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

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

2 梁端帶鉸時(shí)的單元切線剛度矩陣

3 非線性方程組求解和計(jì)算流程

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

本文非線性計(jì)算流程如下：

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

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

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

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

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

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

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

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

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

圖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 個(gè)經(jīng)典算例對本文算法及程序從正確性、計(jì)算效率及精度方面進(jìn)行較全面和嚴(yán)格的檢驗(yàn)。算例1. 如圖4 所示為Lee’s 框架[26]，材料及幾何參數(shù)如下：線彈性模量E=7060.8 kN/m2，梁和柱具有同樣的截面，且截面面積為6 cm2，慣性矩為2 cm4，對結(jié)構(gòu)進(jìn)行幾何非線性分析，分析中柱均分為8 個(gè)單元，荷載作用點(diǎn)左邊梁均分為2 個(gè)單元，右邊梁均分為6 個(gè)單元，圖5(a)所示為框架的初始構(gòu)型以及不同變形時(shí)刻的構(gòu)型，圖5(b)、圖5(c)所示分別為荷載作用點(diǎn)處荷載-水平位移及荷載-豎向位移曲線，可看出與相關(guān)文獻(xiàn)[27]的計(jì)算結(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)行計(jì)算；

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

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

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

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

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

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

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

為驗(yàn)證本文所導(dǎo)出的帶鉸梁單元模型，基于荷載及結(jié)構(gòu)對稱的特征，取圖7(b)所示的1/2 結(jié)構(gòu)進(jìn)行計(jì)算，按每根桿件分成1、2、3 和4 個(gè)單元(中間鉸左、右兩個(gè)相鄰單元類型選用本文推導(dǎo)的帶鉸單元)，在荷載參數(shù)PL2/EI=10 時(shí)所得的計(jì)算結(jié)果如表3 所示，其中解析解見文獻(xiàn)[29]。從表3可看出，計(jì)算結(jié)果收斂非常好，如按工程精度誤差不超過5%的要求，則每根桿件劃分成3 個(gè)單元就已經(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é)點(diǎn)力之間及節(jié)點(diǎn)位移之間的總量關(guān)系及微分導(dǎo)出的增量關(guān)系，獲得平面梁單元在結(jié)構(gòu)坐標(biāo)系下的切線剛度矩陣，同時(shí)還獲得單元桿端力的全量算法；在此基礎(chǔ)上，根據(jù)帶鉸梁端彎矩為零的受力特征，導(dǎo)出了能考慮梁端帶鉸的單元切線剛度矩陣表達(dá)式。多個(gè)算例表明本文算法是正確的，相對于已有文獻(xiàn)算法，本文算法在非線性計(jì)算能力、精度以及穩(wěn)定性方面均占有一定優(yōu)勢，可用于平面桿系結(jié)構(gòu)的幾何非線性分析。