彈性膜結構幾何非線性分析的剛體準則法

陳朝暉，楊 帥，楊永斌,3

(1. 重慶大學土木工程學院，重慶 400045；2. 山地城鎮建設與新技術教育部重點實驗室(重慶大學)，重慶 400045； 3. 國立云林科技大學營建工程系，臺灣 64002)

膜結構具有輕盈、便于覆蓋大跨屋蓋、造型多變以及可以預制等許多顯著優點，近30 年得到快速發展，在工程實踐中發揮了重要作用。由于膜結構面外剛度很小，即使在較小橫向荷載下，也易發生較大撓曲變形，因而，膜結構通常通過較大變形來抵抗荷載，具有明顯幾何非線性特性。

Masahisa 等[1]提出了張拉結構的非線性分析方法，利用更新拉格朗日列式以考慮膜結構的大變形影響。Tabarrok 和Qin[2]采用非線性有限元方法對膜、索和框架組成的組合張拉結構進行找形與荷載分析，但所推導的膜單元為復雜的高階單元，數值積分時需要較多插值點，效率與精度無法兼顧。Noguchi 等[3]利用無網格伽遼金法(element-free Galerkin method，EFG)對空間膜結構進行幾何非線性分析，突破了傳統有限元分析中單元網格的限制，取得了較為滿意的結果，但在基函數和影響域方面，未明確最佳順序和大小，影響了該方法的穩定性和通用性。Tsiatas 等[4]提出了空間任意形狀薄膜大撓度分析方法，利用變分法推導了新的控制微分方程，并采用模擬方程法(the analogue equation method，AEM)對微分方程進行直接積分求解，即在相同的邊界條件下，用三個非耦合的等效線性平面膜方程替換三個耦合的變系數非線性偏微分方程。該方法精度較高，計算效率高于其他數值方法，一般情況下對于預應力膜，解的收斂性較好，但隨著預應力的減小，收斂性變差。王震[5]和趙陽等[6]基于向量式有限元法(the vector form intrinsic finite element method，VFIFEM)，推導了三角形常應變膜單元和四結點四邊形等參膜單元，詳細描述了變形坐標系下單元結點內力的求解方法，同時對四結點膜單元的位移模式和內力計算的數值積分等問題提出了合理可行的處理方法，有效模擬了空間膜結構大位移和大轉動行為，但該方法計算效率與普通有限元法相比并無優勢。

增量迭代法是結構非線性分析的通行方法，基本思想是將非線性過程在各荷載增量階段線性化，該方法大致可分為三個階段：增量預測階段、單元結點力計算階段和誤差判斷階段。單元結點力與荷載之差通過迭代逐步減小。其中，單元結點力的計算決定了整個計算的精度。對于考慮大變形和大轉動的單元，常見的各類幾何剛度矩陣都存在不同程度的近似，在增量迭代過程中，如果這種近似造成了單元自身不平衡，所產生的誤差是無法通過迭代消除的，甚至可能導致計算發散。解決方法之一是細分單元，但會極大降低計算效率。

為此，Yang 等[7－8]于1987 年提出了幾何非線性分析的“剛體準則”，即發生剛體轉動的單元，在初始狀態下滿足平衡條件的結點力，僅隨單元轉動而不會改變大小。Yang 等[9]運用虛功原理推導了滿足上述剛體準則的桿系結構、板殼結構[10]等結構相應單元幾何剛度矩陣，并基于更新的拉格朗日列式(Updated Languagian，簡記作UL)，建立了與剛體準則相匹配的非線性增量迭代方法。這一思想方法，從根本上解決了大變形過程中單元的平衡問題，有效提高了收斂速度和精度。Yang 等[9]由此建立了一些列滿足剛體準則的桿單元、平面及空間梁單元[8]和殼單元[10]，Yang 和Chen 等[11]結合塑性鉸理論，將彈性幾何非線性桿單元進一步推廣到框架結構的彈塑性幾何非線性分析中。所建立的單元幾何剛度矩陣形式簡潔，幾何非線性分析簡單明了。與現有非線性分析方法相比，計算效率優勢顯著[11－12]。

鑒于前述各膜單元的不足以及剛體準則在幾何非線性分析中簡明優美的出色表現，本文應用剛體準則，建立了一種新的適于工程應用的彈性非線性分析膜單元。該單元由滿足剛體準則的空間桿單元和常規彈性膜單元構成，其中，膜單元的幾何剛度矩陣由桿單元構成，彈性剛度矩陣則由平面應力單元構成。經驗證，由此構成的空間膜單元和平面膜單元滿足剛體準則，即彈性剛度矩陣和幾何剛度矩陣在單元剛體位移上虛功為零。若干經典算例分析以及與已有分析方法的對比，驗證了本文所建單元以及相應分析方法的精度及效率。

1 剛體準則的基本思想

如圖1 所示剛體單元，在1C狀態下保持平衡，當單元僅發生剛體位移運動到2C狀態時，由于單元沒有自然變形，不會產生新的結點力，仍應保持平衡。因此，單元結點力應為伴隨力，即隨單元一起運動至2C狀態，大小不變，僅發生方向的改變。此即“剛體準則”[7]。

圖1 發生剛體位移的初始平衡剛體單元 Fig.1 Rigid body motion of a rigid element with initial forces

對于發生大變形或大位移的單元，設在1C狀態平衡，由1C狀態至2C狀態，單元將經歷大變形或大位移，可將這一過程視為如下兩個過程的疊加：單元先由1C至2C發生剛體位移，而后在2C狀態下發生有限變形。在剛體位移階段，1C狀態平衡的單元結點力同單元一起轉動而大小不變，在2C狀態下仍然平衡；而后，再計算2C狀態下的有限變形以及由變形引起的單元結點力增量。

對于大多數工程中的大變形和大轉動問題，當增量步較小時，可以認為剛體位移占單元總的變形或位移的絕大部分，而有限變形僅占小部分。這樣處理，物理概念明確清晰，將極大減小2C狀態下的結點不平衡力，有效提高迭代計算的效率，即使產生變形描述的誤差，也可以通過細分單元或減小步長來有效降低。

2 三角形空間膜單元

2.1 三角形空間膜單元的構造

如圖 2 所示三結點三角形空間膜單元(Triangular space membrane element，簡稱TSME)，圖中，OXYZ和Cxyz分別為結構整體坐標系，單元坐標系，其中x軸和y軸在單元平面內，z軸垂直于單元平面，坐標原點為單元形心C。設每個結點包含三個自由度，即兩個面內平動自由度和一個面外平動自由度。結點i在單元坐標系下的坐標為(xi,yi,0)，ui、vi和wi分別為結點i沿x軸、y軸和z軸的位移，i=1,2,3。

結構的幾何剛度矩陣本質上體現了由于結構幾何形狀的改變而產生的荷載勢能，是荷載在結構變形上所做的虛功，因此，建立空間膜單元幾何剛度矩陣的關鍵在于對膜單元整體變形的描述。根據前述剛體準則思想，可認為膜單元在變形過程中，其剛體位移對其整體變形的貢獻較大，而單元的彈性變形貢獻較小。進一步地，由三角形穩定性可知，若三角形三條邊的位置和長度確定，三角形的形狀及位置即可唯一確定。由此，可將三角形膜單元視為由三根空間桿件組成鉸接三角形，并在中間張拉薄膜而成，桿件的材料與薄膜相同，如圖3 所示。如此建立的空間膜單元的整體位形與位移即可由桿單元構成的空間鉸接三角形確定，而膜單元的有限彈性變形可由內部張拉的薄膜變形確定。因此，本文所建空間膜單元的幾何剛度矩陣可由桿件鉸接三角形的幾何剛度矩陣得到，彈性剛度矩陣則為內部張拉薄膜的彈性剛度矩陣，后者與通常的線性膜單元相同。以下即由空間桿單元幾何剛度矩陣推導建立空間膜單元的幾何剛度矩陣。

圖2 三結點三角形空間膜單元 Fig.2 The TSME with three nodes

圖3 三角形空間膜單元構成 Fig.3 Construction of TSME

2.2 三角形空間膜單元的幾何剛度矩陣

單元坐標系下三角形空間膜單元的結點力如圖4(a)所示，將其中的鉸接三角形單元拆分為三個空間桿單元，相應桿端力如圖4(b)所示。本文各物理量的左上標表示當前狀態，左下標表示參考狀態，則圖中11xF表示在1C狀態膜單元中結點1 沿x方向的結點力，表示在C1狀態桿單元12 結點1 處沿x方向的桿端力，其余類似。

C1狀態下三角形空間膜單元的結點力向量｛｝可記作：

圖4 三角形空間膜單元結點力和桿端力 Fig.4 The nodal and internal forces of TSME

膜單元與桿件需滿足如下平衡條件：

三角形空間膜單元與桿單元在各結點上需滿足平衡方程：

式中，右上標i,j,k分別用1,2,3 輪換，以下同此。

空間桿單元ij的桿端力需滿足平衡方程：

三角形空間膜單元的結點力需滿足平衡方程：

聯立方程式(2)～式(4)，得到空間桿單元ij的桿端力為：

其中，xf、yf和zf是式(2)～式(4)求解過程中引入的待定常數，可令其為零。

由圖4 可知，式(5)中空間桿單元的桿端力是相對于膜單元的單元坐標系。進一步地，需將其轉換到各桿單元所在的局部坐標系下，如圖5 所示，圖中，、和是桿單元局部坐標軸。用｛｝和分別表示空間桿單元ij在膜單元的單元坐標系與桿單元局部坐標下的桿端力向量：

圖5 局部坐標系下空間桿單元的桿端力 Fig.5 The internal force of 3D bar in local coordinate

兩組向量之間有如下坐標變換關系：

其中，坐標轉換矩陣 3D[ ]T具體如下：

式中，xij=xi-xj，yij=yi-yj，Lij為空間桿單元ij的長度。

將桿單元ij的幾何剛度矩陣 [kg]3Dbar_ij從桿單元局部坐標系轉換到三角形膜單元的單元坐標系，按照單元的結點自由度編碼，“對號入座”疊加得到三角形空間膜單元的幾何剛度矩陣 [kg]TSME：

其中：

2.3 三角形空間膜單元的彈性剛度矩陣

三角形空間膜單元的彈性剛度矩陣可由三角形平面應力單元的彈性剛度矩陣[13]修正得到。將三角形平面應力彈性剛度矩陣擴階，即在平面彈性剛度矩陣的基礎上增加第3、6、9 行和列，并將z坐標所對應元素置0，可得：

其中：

2.4 剛體準則檢驗

設1C狀態下單元處于平衡狀態，基于UL 格式，建立結構在2C狀態的增量虛功方程，則單元的增量平衡方程為[8]：

式中：｛u｝ 表示C1狀態到C2狀態單元結點位移增量； {}為以C1狀態為參考的C2狀態下的結點力向量；｛｝為以C1狀態為參考的C1狀態下的結點力向量，為：

假設膜單元繞z軸逆時針剛性轉動 180o到達C2狀態，如圖6 中實線所示，此過程中單元剛體位移向量{u}r為：

圖6 三角形空間膜單元的剛體位移 Fig.6 Rigid body motion of TSME

易證明彈性剛度矩陣[ke]在剛體位移下不會產生單元結點力增量，即：

當單元只發生剛體位移{ }ru時，式(13)左側變為：

其中：

顯然，單元在發生剛體位移后，仍然滿足平衡條件，即彈性剛度矩陣和幾何剛度矩陣在剛體位移上所作虛功為零，由此證明了本文所推導的空間膜單元幾何剛度矩陣滿足虛功原理。

3 三角形平面膜單元

三角形平面膜單元(Triangular plane membrane element, 簡稱TPME)可由三角形空間膜單元退化而成，如圖7 所示，設三角形平面膜單元每個結點包含兩個平面內平動自由度。

與三角形空間膜單元的構成類似，可將三角形平面膜單元看作是在平面桿單元構成的鉸接三角形內部張拉薄膜而成，如圖8 所示。其單元幾何剛度矩陣即為平面桿單元構成的鉸接三角形的幾何剛度矩陣。單元坐標系下三角形平面膜單元的結點力如圖8(a)所示，平面桿單元的桿端力如圖8(b)所示。

圖7 三結點三角形平面膜單元 Fig.7 The TPME with three nodes

C1狀態下三角形平面膜單元的結點力向量｛｝可記作：

圖8 三角形平面膜單元結點力和桿端結點力 Fig.8 The nodal and internal force of TPME

與空間膜單元幾何剛度陣推導過程類似，先根據平衡條件，建立平面膜單元坐標系下膜單元結點力與各平面桿單元桿端力之間的關系；而后，將各桿單元的桿端力轉換到各桿單元的局部坐標系下，如圖9，并將其代入文獻[9]中的平面桿單元幾何剛度矩陣，得到平面桿單元ij在桿件局部坐標系下的幾何剛度矩陣 [kg]2Dbar_ij：

圖9 平面桿單元局部坐標系下的桿端力 Fig.9 The internal forces of 2D bar in local coordinate

將桿件局部坐標系下的幾何剛度矩陣 [kg]2Dbar_ij轉換到膜單元坐標系下，并按單元結點自由度編碼對號入座，得到三角形平面膜單元的幾何剛度矩陣[kg]TPME：

其中，轉換矩陣 2D[ ]T為：

則三角形平面膜單元的幾何剛度矩陣 [kg]TPME的具體表達式如下：

其中：

三角形平面膜單元的彈性剛度矩陣即為三角形平面常應力單元的彈性剛度矩陣[13]。

4 UL 列式增量迭代法

非線性分析中，UL 列式的增量迭代法包括三種狀態，分別是結構的初始狀態0C、上一步平衡狀態1C和當前狀態2C。非線性分析的增量過程，就是從1C狀態到2C狀態的結構變形過程。每一增量步中，結構的變形很小，但所有增量步累積可以產生很大的變形。根據前述剛體準則，可將每個增量步所產生的位移增量視為剛體位移和自然變形兩部分，其中，剛體位移相對于自然變形大得多。如果在分析的每個階段都充分考慮剛體位移的作用，則可以采用小變形線性化理論來處理自然變形的剩余影響[14]。

如圖10 所示，增量迭代法分為三個階段：預測階段、結點力計算階段和誤差判斷階段。預測階段即計算在給定荷載增量下結構的結點位移增量：

式中， {1P} 和 {2P} 分別表示在C1和C2狀態下的外荷載。由結構位移增量{ }UΔ 可以得到各單元的位移增量{ }uΔ ，并疊加得到結構結點位移。顯然，預測階段主要影響迭代次數和收斂速度[14]，在此，結構整體剛度矩陣[ ]K只需近似即可，但需保障迭代方向的正確。對于非線性不顯著的問題，可以只采用彈性剛度矩陣，對于非線性顯著的結構，結構整體剛度矩陣還需包含結構幾何剛度矩陣。

圖10 增量-迭代法示意圖 Fig.10 Key steps of incremental-iterative analysis method

單元結點力計算階段是以C2狀態為參考，計算在C2狀態時的單元結點力。結點力計算決定了整個幾何非線性分析的精度。根據剛體準則，由于C1狀態平衡的單元結點力，經過剛體運動到C2狀態時，僅發生方向的改變而大小不變，因此可以將C1狀態的結點力 {}通過坐標變化轉動到C2狀態。而后基于小變形線性化理論，采用彈性剛度矩陣計算單元結點力增量：

將兩部分相加得到以2C為參考狀態的2C狀態下的單元結點力：

誤差判斷階段需計算結點不平衡力。由方程(26)得到2C狀態下結構各單元結點力，求和后與總載荷 {2P} 比較，得到C2狀態下結構的結點不平衡力。若不平衡力大于收斂值，則進入下一增量步；否則，迭代以消除不平衡力。

以上過程的具體步驟為：

對于第i增量步、第j迭代步的結構增量平衡方程寫作[15]：

將式(27)進一步分解為：

對于第j次( 2j≥ )迭代計算，荷載增量參數i

jλ為：

5 算例與分析

以下通過若干算例驗證本文所建空間及平面膜單元在幾何非線性分析中的適用性、精度及效率。首先，通過懸臂梁端部受集中荷載作用的算例驗證本文所建平面膜單元的有效性，再通過內壓作用下雙拋物線膜結構、雙拋物面膜結構等的大變形分析，并與Noguchi 等[3]的無網格伽遼金法(EFG)、Tsiatas 等[4]的有限元法(FEM)和模擬方程法(AEM)、王震等[5,25]的向量式有限元法(VFIFEM)以及商業軟件ANSYS 等對比，說明本文所建模型的精度及效率。

算例1. 平面膜單元算例—懸臂梁端部受集中荷載

該算例取自文獻[24]，為驗證幾何非線性分析單元及方法的經典算例。如圖11，懸臂梁自由端受集中荷載，考慮平面內彎曲，彈性模量E=10 MPa，泊松比ν=0.3，懸臂梁長∶高∶厚=1000 mm∶100 mm∶10 mm。采用本文所建三角形平面膜單元分析，共劃分20 個單元，如圖12 所示。此外，還采用ANSYS 中的SHELL63 單元進行對比，單元劃分形式同圖12，SHELL63 單元每個結點6 個自由度。僅考慮膜的面內剛度，所得懸臂梁端部荷載與豎向位移關系曲線如圖13 所示。從中可以看出，本文推導的三角形平面膜單元與 ANSYS 中的SHELL63 單元計算結果高度吻合，但本文所建單元 每個結點僅含兩個自由度，且幾何剛度矩陣為線性舉證，構造簡單。

圖11 受剪懸臂梁 Fig.11 Cantilever beam subjected to transverse end load

圖12 三角形平面膜單元劃分方案 Fig.12 Triangular plane membrane element

圖13 懸臂梁端部荷載位移曲線 Fig.13 Load-deflection curves for cantilever beam

算例2. 空間膜單元算例—內壓下的雙拋物線膜片

如圖14 所示雙拋物線薄膜結構，膜片的投影平面為矩形，其初始形狀由以下雙拋物線曲面函數所確定[4.26]：

設膜片邊長a=b=1000 mm，厚度h=1 mm，中心矢高h=100 mm，彈性模量E=6000 MPa，泊松比ν=0.267，四邊簡支，膜片內表面受均布內壓作用。采用本文所建三角形空間膜單元，共劃分200 個單元，如圖14 所示，采用前述增量迭代法計算，可以觀察到，隨著荷載增大，膜片向外膨脹。圖15為內壓為1.5 MPa 時的變形狀態。雙曲拋物線膜片形心處荷載-豎向位移曲線繪于圖16 中，圖中還給出了有限元法(FEM)、Noguchi 等[3]的無網格伽遼金法(EFG)、Tsiatas 等[4]的模擬方程法(AEM)以及王震[5]的向量式有限元法(VFIFEM)的結果予以比較。各方法在不同荷載下的計算結果及與AEM 方法的相對誤差見表1。與精度較高的模擬方程法(AEM)[4]的比較可知，無網格伽遼金法(EFG)[3]由于采用忽略膜曲率的應變模型和平面彈性應力-應變關系而計算精度最低，向量式有限元法(VFIFEM)[5]通過劃分了800 個三角形常應變薄膜單元而獲得了較高的計算精度。本文方法與Tsiatas 等[4]采用200 個三角形殼單元的有限元法所得結果具有相近的精度，但本文所建單元自由度少，幾何剛度矩陣構造簡單，具有更高的計算效率。

圖14 膜片初始未變形狀態 Fig.14 Initial undeformed state of the membrane

圖15 內壓為1.5 MPa 時膜片變形狀態 Fig.15 Deformation state under 1.5 MPa of the membrane

圖16 雙曲拋物線膜片形心處荷載-豎向位移曲線 Fig.16 Central load-deflection curve of the double parabolic membrane

表1 不同計算方法在不同荷載下的 雙曲拋物線膜形心豎向位移 Table 1 Central deflection of the double parabolic membrane by different method

算例3. 空間膜單元算例—內壓下的雙曲拋物面膜片

如圖17 所示雙曲拋物面薄膜結構，膜片的投影平面形狀為矩形，其初始形狀由以下雙曲拋物面函數所確定[4,27]：

其 中，相 對 矢 高 Δf=C3-C2-C4，C3= 0，C2=C4= 300 mm 。矩形膜片邊長a=b=1000 mm，厚度h=3 mm，彈性模量E=6000 MPa，泊松比ν=3，四邊簡支，膜片內表面受均布內壓作用。采用三角形空間膜單元，共劃分200 個單元。圖18 為內壓達到3 MPa 時膜片的變形狀態。雙曲拋物面膜片形心處荷載-位移曲線繪于圖19 中，圖中還給出了采用模擬方程法(AEM 法)[4]、有限元法(FEM)[4]以及 向量式有限元法(VFIFEM)[5]的計算結果。不同方法在不同荷載下的計算結果與相對誤差見表2。與精度較高的AEM 的比較可知，VFIFEM 法[5]劃分了800 個單元而獲得了較高的計算精度，FEM 法[4]采用了高次插值的三角形膜單元，同樣劃分200 個單元，但計算精度略低于本文方法。表明本文所建空間膜單元對膜結構幾何非線性分析具有較高的精度和效率。

圖17 初始未變形狀態 Fig.17 Initial undeformed state of the membrane

圖18 內壓為3.0 MPa 時變形狀態 Fig.18 Deformation state under 3.0 MPa of the membrane

圖19 雙曲拋物面膜片中心荷載-位移曲線 Fig.19 Central load-deflection curve of the hyperbolic paraboloidal membrane

表2 不同計算方法在不同荷載下的雙曲拋物 面膜形心豎向位移 Table 2 Central deflection of the hyperbolic paraboloidal membrane by different method

6 結論

本文應用剛體準則，基于UL 列式，建立了一種新的彈性非線性分析膜單元及其相應非線性分析增量迭代法。本文所構造的三角形膜單元，由三根空間桿件組成鉸接三角形，并在中間張拉薄膜而成，桿件的材料與薄膜相同，這符合三角形的穩定性，即三角形的形狀及位置由其三條邊的位置和長度可唯一確定。該膜單元的整體位形由桿單元構成的空間鉸接三角形確定，其幾何剛度矩陣由桿件鉸接三角形的幾何剛度矩陣推導得到；而膜單元的彈性變形則由內部張拉的薄膜變形確定，其彈性剛度矩陣即為常規膜單元的彈性剛度矩陣。由此建立的空間膜單元和平面膜單元均通過了剛體準則檢驗。即初始平衡的單元，結點力在剛體位移上的虛功為零。

在此單元基礎上，本文所建立的非線性分析方法，以UL 列式的增量迭代法為依托，在非線性分析的增量迭代法的每個階段充分考慮了剛體轉動的影響，然后利用小變形線性化理論處理剩余的自然變形的影響。在UL 型增量迭代法的預測階段，所用的剛度矩陣[ ]K為彈性剛度矩陣和幾何剛度矩陣的組合，該階段幾何剛度矩陣[kg]不必很準確，可以采用現有軟件通行的剛度矩陣，只要保障增量方向正確即可；在單元結點力計算階段，則需合理計算由單元大轉動大位移引起的結點力增量，為此，首先利用剛體準則更新1C狀態的單元結點力，而后基于小變形線性化理論，僅用彈性剛度矩陣計算結點力增量。

上述基于剛體準則的思想構造的三角形空間膜單元，其幾何剛度矩陣推導過程物理意義明確清晰，且由空間膜單元容易地退化得到平面膜單元。所作經典算例分析以及與已有分析方法、商業軟件的對比表明，所建三角形空間膜單元能有效的模擬空間膜結構的大變形和大轉動行為，且顯示了較高的計算精度和計算效率。