對(duì)稱幾何非線性三角形平面單元研究

王動(dòng)，鄧?yán)^華

（1.淮北市交通投資控股集團(tuán)有限公司，安徽 淮北 235000；2．長(zhǎng)沙理工大學(xué) 土木工程學(xué)院，湖南 長(zhǎng)沙 410114）

隨著高強(qiáng)度、輕質(zhì)材料的廣泛應(yīng)用，工程結(jié)構(gòu)不斷向大跨、柔性的方向發(fā)展，結(jié)構(gòu)的非線性效應(yīng)也愈加明顯，對(duì)結(jié)構(gòu)進(jìn)行整體或局部的動(dòng)力、靜力分析時(shí)必須考慮非線性效應(yīng)的影響，才能滿足安全、經(jīng)濟(jì)的設(shè)計(jì)要求。在這一過程中，建立各種精確、高效的單元模型，并考慮其所面對(duì)的非線性效應(yīng)是關(guān)鍵［1-2］。目前，采用諸如ANSYS 等大型商業(yè)數(shù)值軟件對(duì)結(jié)構(gòu)進(jìn)行非線性分析是主流的研究方法。但事實(shí)上，完全拉格朗日法（total lagrange，TL）、修正的拉格朗日法（updated lagrange， UL）以及共旋坐標(biāo)法（corotational， CR 法）才是各種非線性單元分析的主要算法。相對(duì)于TL 法和UL 法而言，CR 法出現(xiàn)的時(shí)間較晚，且相關(guān)研究較少。因此，TL 與UL 這兩種算法被ANSYS 等商業(yè)數(shù)值軟件廣泛采用，而到目前為止，CR法尚未被商業(yè)數(shù)值軟件采用。

FELIPPA 等［3］認(rèn)為CR 法僅適用于大轉(zhuǎn)動(dòng)、小應(yīng)變的運(yùn)動(dòng)條件，苛刻的應(yīng)用條件使其不能滿足大型商業(yè)工程仿真數(shù)值軟件的算法適用范圍必須很寬的要求，這也是其被眾多商業(yè)工程仿真數(shù)值軟件棄用的主要原因。但RANKIN 等［4-6］認(rèn)為在物體發(fā)生大轉(zhuǎn)動(dòng)、小應(yīng)變的運(yùn)動(dòng)條件下，CR 法相對(duì)于TL 法和UL 法具有兩個(gè)非常明顯的優(yōu)勢(shì)：一是在進(jìn)行非線性單元研究時(shí)，CR 法能非常方便地引入各種先進(jìn)線性單元，而無須考慮推導(dǎo)這些先進(jìn)線性單元時(shí)所采用的各種假設(shè)條件；二是在非線性列式中，幾何非線性與材料非線性是解耦的，即材料非線性在共旋坐標(biāo)系下進(jìn)行計(jì)算，而幾何非線性則通過變量在共旋坐標(biāo)系與結(jié)構(gòu)坐標(biāo)系之間的轉(zhuǎn)換矩陣被自動(dòng)計(jì)入。由于工程結(jié)構(gòu)中的非線性問題絕大部分都屬于大轉(zhuǎn)動(dòng)、小應(yīng)變的范圍［7］，因此，基于CR 法研究各種非線性單元，對(duì)促進(jìn)結(jié)構(gòu)非線性效應(yīng)研究具有十分重要的意義。楊浩文等［8］采用共旋坐標(biāo)法建立了考慮幾何非線性的Timoshenko 梁?jiǎn)卧?，并基于該單元建立了用于非線性動(dòng)力分析的能量守恒逐步積分算法。馮曉東等［9］基于共旋坐標(biāo)法，推導(dǎo)了大位移空間桿單元的切線剛度矩陣，并將其結(jié)論運(yùn)用在一個(gè)對(duì)四桿張拉整體結(jié)構(gòu)單元體的幾何非線性彈塑性特性進(jìn)行分析的研究上，該研究結(jié)果表明，相比于傳統(tǒng)的TL 法和UL法，CL法具有更高的計(jì)算效率。鄧?yán)^華等［10］將共旋坐標(biāo)法和穩(wěn)定函數(shù)結(jié)合起來，利用共旋坐標(biāo)法考慮大位移，利用穩(wěn)定函數(shù)考慮P-δ效應(yīng)，多個(gè)經(jīng)典算例分析均表明該非線性平面梁?jiǎn)卧哂休^高的計(jì)算精度和效率。張亮等［11］則采用共旋坐標(biāo)法，推導(dǎo)了一個(gè)適用于充氣薄膜結(jié)構(gòu)褶皺分析的空間三節(jié)點(diǎn)的三角形膜單元。目前，也有一些學(xué)者基于共旋法對(duì)平面單元方面進(jìn)行了研究。蔡松柏等［12］基于共旋法與場(chǎng)一致性原則，分別建立了四節(jié)點(diǎn)四邊形和八節(jié)點(diǎn)四邊形幾何非線性平面單元。鄧?yán)^華等［13-14］還將僅限于幾何非線性的研究拓展到幾何與材料雙非線性領(lǐng)域。相對(duì)于四邊形平面單元而言，基于共旋法分析幾何非線性三角形平面單元的研究很少，筆者僅查閱到文獻(xiàn)［15］，該文獻(xiàn)基于共旋法與場(chǎng)一致性原則，對(duì)幾何非線性三角形平面單元進(jìn)行了研究。不過需要指出的是，鄧?yán)^華等［15］與蔡松柏等［12-14］學(xué)者一樣，其推導(dǎo)出的單元切線剛度矩陣均是非對(duì)稱的矩陣，這將給非線性方程組的求解帶來一定的困難。眾所周知，平面三角形單元作為常應(yīng)變單元，當(dāng)網(wǎng)格化較為稀疏時(shí)，其計(jì)算精度較低；但隨著單元網(wǎng)格數(shù)量的增加，其計(jì)算精度會(huì)迅速提高。且三角形平面單元可適用于復(fù)雜和不規(guī)則邊界形狀的分析［16］。因此，基于CR 法，建立具有對(duì)稱切線剛度矩陣的幾何非線性三角形平面單元是十分必要的。

綜上所述，本研究在這些研究的基礎(chǔ)上，針對(duì)三角形平面單元，借鑒YAW 等［17］的研究方法，采用具有對(duì)稱切線剛度矩陣的幾何非線性梁?jiǎn)卧?，合理選擇三角形平面單元共旋坐標(biāo)系的原點(diǎn)及坐標(biāo)軸方向，基于與虛功原理相結(jié)合的幾何一致性原則［18］，建立三角形平面單元在大轉(zhuǎn)動(dòng)、小應(yīng)變條件下，具有對(duì)稱特征的幾何非線性單元切線剛度矩陣和相應(yīng)的節(jié)點(diǎn)抗力算法；對(duì)于該算法中待求解的非線性方程組，采用一種能將位移增量法與荷載增量法形成統(tǒng)一迭代格式的新方法［19］，開發(fā)了相應(yīng)的計(jì)算程序，并將該計(jì)算程序應(yīng)用于大變形懸臂梁計(jì)算的兩個(gè)經(jīng)典算例，將該算法對(duì)這兩個(gè)算例的計(jì)算結(jié)果與ANSYS 商業(yè)數(shù)值軟件的計(jì)算結(jié)果進(jìn)行了對(duì)比，驗(yàn)證了本研究建立的有限元列式和算法的正確性，且算例2的計(jì)算比較結(jié)果也表明該算法的計(jì)算精度略高于ANSYS商業(yè)數(shù)值軟件的計(jì)算精度。

1 共旋法有限元列式

1.1 基本圖式

變形前、后的三角形平面單元如圖1所示。

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

在圖1 中，XOY為固定不變的結(jié)構(gòu)坐標(biāo)系，xoy為隨單元平移和轉(zhuǎn)動(dòng)的共旋坐標(biāo)系，取三角形3 個(gè)頂點(diǎn)M1，M2，M3坐標(biāo)分量的平均值作為點(diǎn)C在結(jié)構(gòu)坐標(biāo)系XY下的坐標(biāo)分量值。在初始時(shí)刻，x軸方向平行于X軸；在任意計(jì)算時(shí)刻，x軸的方向由剛性轉(zhuǎn)角θ確定，y軸的方向則由x軸逆時(shí)針旋轉(zhuǎn)90°得到。

1.2 平面單元的切線剛度和抗力矩陣

對(duì)圖1所示三角形平面單元，設(shè)

式中：

為初始時(shí)刻的節(jié)點(diǎn)單元；

Ri為計(jì)算時(shí)刻單元在結(jié)構(gòu)坐標(biāo)系下的節(jié)點(diǎn)坐標(biāo)；

ui為單元在結(jié)構(gòu)坐標(biāo)系下的節(jié)點(diǎn)位移；

為初始時(shí)刻共旋坐標(biāo)系xy的原點(diǎn)C在結(jié)構(gòu)坐標(biāo)系下的坐標(biāo)；

rC計(jì)算時(shí)刻共旋坐標(biāo)系xy的原點(diǎn)C在結(jié)構(gòu)坐標(biāo)系下的坐標(biāo)。故有

在式（2）中，uc=()UC，VC，是在計(jì)算時(shí)刻，共旋坐標(biāo)系原點(diǎn)C在結(jié)構(gòu)坐標(biāo)系下的位移，故有

同時(shí)，Ri=(Xi，Yi)與之間有關(guān)系：

設(shè)=(xi，yi)，i= 1，2，3 為初始時(shí)刻單元節(jié)點(diǎn)在共旋坐標(biāo)系下的坐標(biāo)，顯然有關(guān)系：

假定計(jì)算時(shí)刻單元發(fā)生節(jié)點(diǎn)位移ui后，共旋坐標(biāo)系x軸與結(jié)構(gòu)坐標(biāo)系X軸的夾角θ已求出，則計(jì)算時(shí)刻單元在共旋坐標(biāo)系下的節(jié)點(diǎn)位移= 1，2，3的計(jì)算式為

對(duì)于計(jì)算時(shí)刻x軸與X軸夾角θ的確定，可參照文獻(xiàn)［20］中的方法。夾角θ的選取的原則是：選取能使得值最小的θ，故可將其帶入式（11），這時(shí)式（11）變?yōu)棣鹊暮瘮?shù)，對(duì)其進(jìn)行關(guān)于θ的求導(dǎo)，并令該導(dǎo)數(shù)為0。可得

顯然，由式（12）可得到兩個(gè)解，分別為θ和θ+π，將這兩個(gè)解代入式（11），選取使較小的那個(gè)作為最終解。

為描述方便，令

設(shè)共旋坐標(biāo)系與結(jié)構(gòu)坐標(biāo)系下與Pl和Pg分別對(duì)應(yīng)的節(jié)點(diǎn)力為fl和fg，設(shè)三角形平面單元在共旋坐標(biāo)系下的剛度矩陣為kl，根據(jù)前面所述共旋坐標(biāo)法的性質(zhì)，有

由共旋坐標(biāo)系與結(jié)構(gòu)坐標(biāo)系下虛功相等的原則［18］，有

其中，位移微分的下標(biāo)v表示該變量為虛擬量。

先對(duì)式（12）進(jìn)行微分，得到δθ與δPg之間的關(guān)系，再對(duì)式（11）進(jìn)行微分，可得到δPl與δPg的關(guān)系，其表達(dá)式為

聯(lián)立式（14）～（15），可得到fl和fg之間的關(guān)系，其表達(dá)式為

對(duì)式（16）進(jìn)行微分，并聯(lián)立公式（13），有

設(shè)三角形平面單元在結(jié)構(gòu)坐標(biāo)系下的幾何非線性切線剛度矩陣為kg，即：

聯(lián)立式（17）、（18），則有

且

顯然，要得到kh的具體表達(dá)式，關(guān)鍵是求出δBT。因此，對(duì)式（11）進(jìn)行變分，可得

式中：xdi=+xi；

ydi=+yi。

為求出式（21）中的δθ，對(duì)式（12）進(jìn)行微分，則有關(guān)系：

聯(lián)立式（21）～（22）及式（15），可得到B的表達(dá)式：

式中：P=I-AG，且有

其中，

至此，聯(lián)立式（20）、（23），并考慮ATGT=I，可得到kh的具體表達(dá)式：

其中，有

將式（24）代入式（19），可得到三角形平面單元在結(jié)構(gòu)坐標(biāo)系下的幾何非線性切線剛度矩陣為kg，其表達(dá)式為

考慮到式（25）中等號(hào)右邊各項(xiàng)矩陣具有的特殊形式，不難證明kg是對(duì)稱的。

2 非線性計(jì)算流程及方程組求解

先輸入結(jié)構(gòu)的初始幾何、材料、荷載及邊界約束等參數(shù)，再結(jié)合單元網(wǎng)格劃分得到每個(gè)單元初始時(shí)刻在共旋坐標(biāo)系下的節(jié)點(diǎn)坐標(biāo)、原點(diǎn)坐標(biāo)、單元?jiǎng)偠染仃嚰翱偤奢d矩陣，保留這些數(shù)據(jù)。在每一級(jí)荷載增量步中，迭代的具體步驟是：

步驟1：由上一級(jí)荷載增量步末的總節(jié)點(diǎn)位移確定當(dāng)前單元的節(jié)點(diǎn)位移ui=()Ui，Vi，i= 1，2，3；

步驟2：由式（8）得到當(dāng)前單元共旋坐標(biāo)系原點(diǎn)C在結(jié)構(gòu)坐標(biāo)系下的位移(UC，VC)；

步驟3：由式（12）并結(jié)合所得θ值需使較小的原則確定單元的剛體轉(zhuǎn)角值θ；

步驟4：由式（11）得到當(dāng)前單元在共旋坐標(biāo)系下的節(jié)點(diǎn)位移

步驟5：由式（13）得到當(dāng)前單元在共旋坐標(biāo)系下的節(jié)點(diǎn)力fl；

步驟6：對(duì)式（11）微分，基于δPl與δPg的關(guān)系，得到矩陣B，進(jìn)而由式（16）得到當(dāng)前單元在結(jié)構(gòu)坐標(biāo)系下的節(jié)點(diǎn)力fg；

步驟7：聯(lián)立式（23）～（24）與（25）得到當(dāng)前單元在結(jié)構(gòu)坐標(biāo)系下的切線剛度矩陣kg；

步驟8：對(duì)所有單元進(jìn)行步驟1～7的操作，將每個(gè)單元在結(jié)構(gòu)坐標(biāo)系下的切線剛度矩陣kg進(jìn)行疊加，得到結(jié)構(gòu)總剛矩陣與總抗力矩陣；

步驟9：結(jié)構(gòu)非線性平衡方程組右端的不平衡力矩陣則由總荷載矩陣減去總抗力矩陣形成；

步驟10：引入邊界條件后，進(jìn)行非線性平衡方程組的求解，得到當(dāng)前迭代步的增量節(jié)點(diǎn)位移，將其與之前的總節(jié)點(diǎn)位移疊加，形成新的總節(jié)點(diǎn)位移；

步驟11：進(jìn)行收斂判斷，如收斂，轉(zhuǎn)至下一個(gè)荷載增量步；如不收斂，進(jìn)入下一個(gè)迭代步。

最常用的非線性方程組求解方法是荷載增量法和位移增量法，這兩種方法均存在一些各自的優(yōu)缺點(diǎn)。為綜合利用這兩種方法的優(yōu)點(diǎn)而又能避免各自的缺點(diǎn)，本研究采用已有的非線性方程組求解方法［19］，該方法將荷載因子λ視為第n+ 1 個(gè)變量，從而將非線性方程組的荷載-位移求解空間由Rn擴(kuò)展到Rn+1，推導(dǎo)出一種能集荷載增量法和位移增量法于一體的統(tǒng)一迭代格式，該方法在實(shí)際研究中得到了較廣泛的應(yīng)用［12，21］。

3 算例分析

3.1 算例1

底端固結(jié)上端自由的柔性立柱，該柔性立柱高為9 m，寬為3 m，厚為1 m，具體如圖2 所示。在該柔性立柱中，自由端截面中心承受的水平力P的數(shù)值為2.777 8×109kN，所使用材料的彈性模量E為1.0×1010kN/ m2，泊松比μ為0。文獻(xiàn)［15］基于其推導(dǎo)的具有非對(duì)稱切線剛度矩陣的三角形平面單元對(duì)該算例進(jìn)行了計(jì)算，并用ANSYS 商業(yè)數(shù)值軟件進(jìn)行了驗(yàn)證。為驗(yàn)證本研究推導(dǎo)公式和算法的合理性與有效性，將本研究的計(jì)算結(jié)果及與文獻(xiàn)［15］的計(jì)算結(jié)果進(jìn)行對(duì)比（這兩種計(jì)算都是基于將柔性立柱在X軸向與Y軸方向，并按0.25 m 的間距對(duì)其進(jìn)行單元網(wǎng)格劃分的，單元網(wǎng)格劃分均產(chǎn)生了481 個(gè)節(jié)點(diǎn)和864個(gè)單元），結(jié)果見表1。

表1 柔性立柱自由端的水平位移Table 1 Horizontal displacement at free end m

圖2 承受水平集中荷載的柔性立柱（單位：m）Fig. 2 Flexible column under horizontal concentrated load（unit： m）

由表1 可知，本研究得到的解與文獻(xiàn)［15］的解吻合良好。這是因?yàn)樵谇蠼饨Y(jié)構(gòu)坐標(biāo)系下的單元節(jié)點(diǎn)抗力時(shí)，本研究與文獻(xiàn)［15］都是基于先在共旋坐標(biāo)系，再進(jìn)行轉(zhuǎn)換得到的。具體都是先用單元線性剛度矩陣乘以去除剛體后的純變形，得到共旋坐標(biāo)系下的單元節(jié)點(diǎn)抗力，再將其轉(zhuǎn)換到結(jié)構(gòu)坐標(biāo)系下進(jìn)行計(jì)算，計(jì)算得到的矩陣與對(duì)應(yīng)階段的荷載矩陣之差即為不平衡力矩陣。CRISFIELD［18］也指出，結(jié)構(gòu)非線性有限元計(jì)算的精度主要由不平衡力的計(jì)算決定，單元切線剛度矩陣的計(jì)算僅決定非線性有限元計(jì)算的效率。因此，在按5 個(gè)等荷載步進(jìn)行加載時(shí)，在同樣的收斂精度要求下，若每一步迭代采用本研究的迭代步驟，則僅需3個(gè)迭代步；而若采用文獻(xiàn)［15］的迭代步驟，則需要不少于5個(gè)迭代步。這也表明了本研究方法在保證計(jì)算精度的前提下，能較好地提高算法的計(jì)算效率。

3.2 算例2

在本算例中，假設(shè)當(dāng)自由端承受水平荷載P（大小為287 500 kN）時(shí)，集中荷載作用的懸臂梁的梁長(zhǎng)L為10 m、彈性模量E為3.45×107kN m2，具體如圖3所示。梁的高度與厚度均為1 m，泊松比μ取為0，并令α=PL2/EI，P是外荷載，L是梁的長(zhǎng)度，E為彈性模量，I為截面抗彎慣性矩。沿懸臂梁的高度方向，將其分別均勻劃分成10、20層，沿其長(zhǎng)度方向均分成50段。

圖3 端部承受集中荷載的懸臂梁（單位：m）Fig. 3 Cantilever beam subjected to concentrated load at free end（unit： m）

在這兩種劃分方式下，分別形成了561個(gè)節(jié)點(diǎn)、1 000個(gè)單元及1 071個(gè)節(jié)點(diǎn)、2 000個(gè)單元的兩種有限元網(wǎng)格。本算例將荷載P分成10 級(jí)增量，采取逐漸施加的方式，并采用荷載增量法對(duì)其結(jié)構(gòu)非線性方程組進(jìn)行求解。對(duì)于梁端截面1/2 高度處節(jié)點(diǎn)A的垂直位移w，本研究計(jì)算值和線型方程的解析解［22］進(jìn)行了對(duì)比，結(jié)果如圖4 所示。在圖4 中，計(jì)算值1 是將懸臂梁劃分為1 000 個(gè)單元的網(wǎng)格計(jì)算得到的，計(jì)算值2則是將懸臂梁劃分為2 000個(gè)單元的網(wǎng)格計(jì)算得到的。

圖4 荷載-豎向位移曲線Fig. 4 Load-vertical displacement diagram

從圖4 中可以看出，第一種網(wǎng)格劃分（1 000 個(gè)單元）下的計(jì)算精度較差，其中，當(dāng)α為1 時(shí)，解析解和數(shù)值解分別為3.017 2 m 和5.375 8 m，相對(duì)誤差為78.2%。當(dāng)α為10時(shí)，解析解和數(shù)值解分別為8.106 1 m和8.328 4 m，相對(duì)誤差為2.74%；但將該懸梁臂進(jìn)行更細(xì)致地劃分，分隔成第二種網(wǎng)格（2 000 個(gè)單元）后，計(jì)算精度得到了提高。其中，當(dāng)α為1 時(shí)，解析解和數(shù)值解分別為3.017 2 m 和3.322 2 m，相對(duì)誤差為10.10%，當(dāng)α為10 時(shí)，解析解和數(shù)值解分別為8.106 1 m 和8.137 3 m，相對(duì)誤差為0.38%。當(dāng)α分別為0（初始狀態(tài)）、1、5及10時(shí)，懸臂梁的變形如圖5所示。

圖5 懸臂梁的變形Fig. 5 Deformation configuration of cantilever beam

為了比較本研究的推導(dǎo)算法與ANSYS 商業(yè)數(shù)值軟件算法的優(yōu)劣，本研究也采用ANSYS 商業(yè)數(shù)值軟件對(duì)該算例進(jìn)行了仿真與計(jì)算。其單元網(wǎng)格劃分與本研究的第二種網(wǎng)格劃分（2 000 個(gè)單元）完全相同，也取10個(gè)荷載步進(jìn)行計(jì)算。為便于與本研究所導(dǎo)出的單元進(jìn)行比較，采用由 Plane 42 單元退化的三角形單元。當(dāng)α為10時(shí)，由ANSYS商業(yè)數(shù)值軟件計(jì)算得到的懸臂梁變形及初始狀態(tài)如圖6所示。

圖6 懸臂梁的變形Fig. 6 Deformation configuration of cantilever beam

從圖6 可以看出，當(dāng)α為10 時(shí)，ANSYS 商業(yè)數(shù)值軟件計(jì)算的節(jié)點(diǎn)A的垂直位移w為8.180 8 m，該值與線性方程組解析解的相對(duì)誤差為0.92%，而本研究算法僅與解析解的相對(duì)誤差僅為0.39%。因此，本研究算法的精度略高于ANSYS 商業(yè)數(shù)值軟件的計(jì)算精度。

4 結(jié)論

1） 通過合理確定共旋坐標(biāo)系原點(diǎn)及坐標(biāo)軸方向，結(jié)合虛功原理與幾何一致性原則，建立基于共旋法的具有對(duì)稱切線剛度矩陣的幾何非線性三角形平面單元，從而避免切線剛度矩陣不對(duì)稱給非線性方程組求解帶來的麻煩。

2） 基于本研究推導(dǎo)公式及算法程序得到的數(shù)值解與已有文獻(xiàn)數(shù)值解及經(jīng)典解析解均吻合良好。且在同等條件下，本研究對(duì)經(jīng)典算例的計(jì)算精度略高于ANSYS 商業(yè)數(shù)值軟件的計(jì)算精度。為提高結(jié)構(gòu)非線性分析的精度和效率，可將本研究方法推廣至高階單元，如六節(jié)點(diǎn)的三角形單元。

3） 利用共旋法特有的幾何非線性與材料非線性解耦的特點(diǎn)和優(yōu)點(diǎn)，可將本研究方法拓展到鋼筋混凝土結(jié)構(gòu)的幾何與材料的雙非線性分析研究中，這也是下一階段工作將重點(diǎn)研究的內(nèi)容。