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

王動，鄧繼華

（1.淮北市交通投資控股集團有限公司，安徽 淮北 235000；2．長沙理工大學 土木工程學院，湖南 長沙 410114）

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

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

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

1 共旋法有限元列式

1.1 基本圖式

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

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

在圖1 中，XOY為固定不變的結構坐標系，xoy為隨單元平移和轉動的共旋坐標系，取三角形3 個頂點M1，M2，M3坐標分量的平均值作為點C在結構坐標系XY下的坐標分量值。在初始時刻，x軸方向平行于X軸；在任意計算時刻，x軸的方向由剛性轉角θ確定，y軸的方向則由x軸逆時針旋轉90°得到。

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

對圖1所示三角形平面單元，設

式中：

為初始時刻的節點單元；

Ri為計算時刻單元在結構坐標系下的節點坐標；

ui為單元在結構坐標系下的節點位移；

為初始時刻共旋坐標系xy的原點C在結構坐標系下的坐標；

rC計算時刻共旋坐標系xy的原點C在結構坐標系下的坐標。故有

在式（2）中，uc=()UC，VC，是在計算時刻，共旋坐標系原點C在結構坐標系下的位移，故有

同時，Ri=(Xi，Yi)與之間有關系：

設=(xi，yi)，i= 1，2，3 為初始時刻單元節點在共旋坐標系下的坐標，顯然有關系：

假定計算時刻單元發生節點位移ui后，共旋坐標系x軸與結構坐標系X軸的夾角θ已求出，則計算時刻單元在共旋坐標系下的節點位移= 1，2，3的計算式為

對于計算時刻x軸與X軸夾角θ的確定，可參照文獻［20］中的方法。夾角θ的選取的原則是：選取能使得值最小的θ，故可將其帶入式（11），這時式（11）變為θ的函數，對其進行關于θ的求導，并令該導數為0。可得

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

為描述方便，令

設共旋坐標系與結構坐標系下與Pl和Pg分別對應的節點力為fl和fg，設三角形平面單元在共旋坐標系下的剛度矩陣為kl，根據前面所述共旋坐標法的性質，有

由共旋坐標系與結構坐標系下虛功相等的原則［18］，有

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

先對式（12）進行微分，得到δθ與δPg之間的關系，再對式（11）進行微分，可得到δPl與δPg的關系，其表達式為

聯立式（14）～（15），可得到fl和fg之間的關系，其表達式為

對式（16）進行微分，并聯立公式（13），有

設三角形平面單元在結構坐標系下的幾何非線性切線剛度矩陣為kg，即：

聯立式（17）、（18），則有

且

顯然，要得到kh的具體表達式，關鍵是求出δBT。因此，對式（11）進行變分，可得

式中：xdi=+xi；

ydi=+yi。

為求出式（21）中的δθ，對式（12）進行微分，則有關系：

聯立式（21）～（22）及式（15），可得到B的表達式：

式中：P=I-AG，且有

其中，

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

其中，有

將式（24）代入式（19），可得到三角形平面單元在結構坐標系下的幾何非線性切線剛度矩陣為kg，其表達式為

考慮到式（25）中等號右邊各項矩陣具有的特殊形式，不難證明kg是對稱的。

2 非線性計算流程及方程組求解

先輸入結構的初始幾何、材料、荷載及邊界約束等參數，再結合單元網格劃分得到每個單元初始時刻在共旋坐標系下的節點坐標、原點坐標、單元剛度矩陣及總荷載矩陣，保留這些數據。在每一級荷載增量步中，迭代的具體步驟是：

步驟1：由上一級荷載增量步末的總節點位移確定當前單元的節點位移ui=()Ui，Vi，i= 1，2，3；

步驟2：由式（8）得到當前單元共旋坐標系原點C在結構坐標系下的位移(UC，VC)；

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

步驟4：由式（11）得到當前單元在共旋坐標系下的節點位移

步驟5：由式（13）得到當前單元在共旋坐標系下的節點力fl；

步驟6：對式（11）微分，基于δPl與δPg的關系，得到矩陣B，進而由式（16）得到當前單元在結構坐標系下的節點力fg；

步驟7：聯立式（23）～（24）與（25）得到當前單元在結構坐標系下的切線剛度矩陣kg；

步驟8：對所有單元進行步驟1～7的操作，將每個單元在結構坐標系下的切線剛度矩陣kg進行疊加，得到結構總剛矩陣與總抗力矩陣；

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

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

步驟11：進行收斂判斷，如收斂，轉至下一個荷載增量步；如不收斂，進入下一個迭代步。

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

3 算例分析

3.1 算例1

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

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

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

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

3.2 算例2

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

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

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

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

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

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

為了比較本研究的推導算法與ANSYS 商業數值軟件算法的優劣，本研究也采用ANSYS 商業數值軟件對該算例進行了仿真與計算。其單元網格劃分與本研究的第二種網格劃分（2 000 個單元）完全相同，也取10個荷載步進行計算。為便于與本研究所導出的單元進行比較，采用由 Plane 42 單元退化的三角形單元。當α為10時，由ANSYS商業數值軟件計算得到的懸臂梁變形及初始狀態如圖6所示。

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

從圖6 可以看出，當α為10 時，ANSYS 商業數值軟件計算的節點A的垂直位移w為8.180 8 m，該值與線性方程組解析解的相對誤差為0.92%，而本研究算法僅與解析解的相對誤差僅為0.39%。因此，本研究算法的精度略高于ANSYS 商業數值軟件的計算精度。

4 結論

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

2） 基于本研究推導公式及算法程序得到的數值解與已有文獻數值解及經典解析解均吻合良好。且在同等條件下，本研究對經典算例的計算精度略高于ANSYS 商業數值軟件的計算精度。為提高結構非線性分析的精度和效率，可將本研究方法推廣至高階單元，如六節點的三角形單元。

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