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

(1. 重慶大學(xué)土木工程學(xué)院，重慶 400045；2. 山地城鎮(zhèn)建設(shè)與新技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(重慶大學(xué))，重慶 400045； 3. 國(guó)立云林科技大學(xué)營(yíng)建工程系，臺(tái)灣 64002)

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

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

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

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

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

1 剛體準(zhǔn)則的基本思想

如圖1 所示剛體單元，在1C狀態(tài)下保持平衡，當(dāng)單元僅發(fā)生剛體位移運(yùn)動(dòng)到2C狀態(tài)時(shí)，由于單元沒有自然變形，不會(huì)產(chǎn)生新的結(jié)點(diǎn)力，仍應(yīng)保持平衡。因此，單元結(jié)點(diǎn)力應(yīng)為伴隨力，即隨單元一起運(yùn)動(dòng)至2C狀態(tài)，大小不變，僅發(fā)生方向的改變。此即“剛體準(zhǔn)則”[7]。

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

對(duì)于發(fā)生大變形或大位移的單元，設(shè)在1C狀態(tài)平衡，由1C狀態(tài)至2C狀態(tài)，單元將經(jīng)歷大變形或大位移，可將這一過程視為如下兩個(gè)過程的疊加：?jiǎn)卧扔?C至2C發(fā)生剛體位移，而后在2C狀態(tài)下發(fā)生有限變形。在剛體位移階段，1C狀態(tài)平衡的單元結(jié)點(diǎn)力同單元一起轉(zhuǎn)動(dòng)而大小不變，在2C狀態(tài)下仍然平衡；而后，再計(jì)算2C狀態(tài)下的有限變形以及由變形引起的單元結(jié)點(diǎn)力增量。

對(duì)于大多數(shù)工程中的大變形和大轉(zhuǎn)動(dòng)問題，當(dāng)增量步較小時(shí)，可以認(rèn)為剛體位移占單元總的變形或位移的絕大部分，而有限變形僅占小部分。這樣處理，物理概念明確清晰，將極大減小2C狀態(tài)下的結(jié)點(diǎn)不平衡力，有效提高迭代計(jì)算的效率，即使產(chǎn)生變形描述的誤差，也可以通過細(xì)分單元或減小步長(zhǎng)來有效降低。

2 三角形空間膜單元

2.1 三角形空間膜單元的構(gòu)造

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

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

圖2 三結(jié)點(diǎn)三角形空間膜單元 Fig.2 The TSME with three nodes

圖3 三角形空間膜單元構(gòu)成 Fig.3 Construction of TSME

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

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

C1狀態(tài)下三角形空間膜單元的結(jié)點(diǎn)力向量｛｝可記作：

圖4 三角形空間膜單元結(jié)點(diǎn)力和桿端力 Fig.4 The nodal and internal forces of TSME

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

三角形空間膜單元與桿單元在各結(jié)點(diǎn)上需滿足平衡方程：

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

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

三角形空間膜單元的結(jié)點(diǎn)力需滿足平衡方程：

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

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

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

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

兩組向量之間有如下坐標(biāo)變換關(guān)系：

其中，坐標(biāo)轉(zhuǎn)換矩陣 3D[ ]T具體如下：

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

將桿單元ij的幾何剛度矩陣 [kg]3Dbar_ij從桿單元局部坐標(biāo)系轉(zhuǎn)換到三角形膜單元的單元坐標(biāo)系，按照單元的結(jié)點(diǎn)自由度編碼，“對(duì)號(hào)入座”疊加得到三角形空間膜單元的幾何剛度矩陣 [kg]TSME：

其中：

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

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

其中：

2.4 剛體準(zhǔn)則檢驗(yàn)

設(shè)1C狀態(tài)下單元處于平衡狀態(tài)，基于UL 格式，建立結(jié)構(gòu)在2C狀態(tài)的增量虛功方程，則單元的增量平衡方程為[8]：

式中：｛u｝ 表示C1狀態(tài)到C2狀態(tài)單元結(jié)點(diǎn)位移增量； {}為以C1狀態(tài)為參考的C2狀態(tài)下的結(jié)點(diǎn)力向量；｛｝為以C1狀態(tài)為參考的C1狀態(tài)下的結(jié)點(diǎn)力向量，為：

假設(shè)膜單元繞z軸逆時(shí)針剛性轉(zhuǎn)動(dòng) 180o到達(dá)C2狀態(tài)，如圖6 中實(shí)線所示，此過程中單元?jiǎng)傮w位移向量{u}r為：

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

易證明彈性剛度矩陣[ke]在剛體位移下不會(huì)產(chǎn)生單元結(jié)點(diǎn)力增量，即：

當(dāng)單元只發(fā)生剛體位移{ }ru時(shí)，式(13)左側(cè)變?yōu)椋?/p>

其中：

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

3 三角形平面膜單元

三角形平面膜單元(Triangular plane membrane element, 簡(jiǎn)稱TPME)可由三角形空間膜單元退化而成，如圖7 所示，設(shè)三角形平面膜單元每個(gè)結(jié)點(diǎn)包含兩個(gè)平面內(nèi)平動(dòng)自由度。

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

圖7 三結(jié)點(diǎn)三角形平面膜單元 Fig.7 The TPME with three nodes

C1狀態(tài)下三角形平面膜單元的結(jié)點(diǎn)力向量｛｝可記作：

圖8 三角形平面膜單元結(jié)點(diǎn)力和桿端結(jié)點(diǎn)力 Fig.8 The nodal and internal force of TPME

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

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

將桿件局部坐標(biāo)系下的幾何剛度矩陣 [kg]2Dbar_ij轉(zhuǎn)換到膜單元坐標(biāo)系下，并按單元結(jié)點(diǎn)自由度編碼對(duì)號(hào)入座，得到三角形平面膜單元的幾何剛度矩陣[kg]TPME：

其中，轉(zhuǎn)換矩陣 2D[ ]T為：

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

其中：

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

4 UL 列式增量迭代法

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

如圖10 所示，增量迭代法分為三個(gè)階段：預(yù)測(cè)階段、結(jié)點(diǎn)力計(jì)算階段和誤差判斷階段。預(yù)測(cè)階段即計(jì)算在給定荷載增量下結(jié)構(gòu)的結(jié)點(diǎn)位移增量：

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

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

單元結(jié)點(diǎn)力計(jì)算階段是以C2狀態(tài)為參考，計(jì)算在C2狀態(tài)時(shí)的單元結(jié)點(diǎn)力。結(jié)點(diǎn)力計(jì)算決定了整個(gè)幾何非線性分析的精度。根據(jù)剛體準(zhǔn)則，由于C1狀態(tài)平衡的單元結(jié)點(diǎn)力，經(jīng)過剛體運(yùn)動(dòng)到C2狀態(tài)時(shí)，僅發(fā)生方向的改變而大小不變，因此可以將C1狀態(tài)的結(jié)點(diǎn)力 {}通過坐標(biāo)變化轉(zhuǎn)動(dòng)到C2狀態(tài)。而后基于小變形線性化理論，采用彈性剛度矩陣計(jì)算單元結(jié)點(diǎn)力增量：

將兩部分相加得到以2C為參考狀態(tài)的2C狀態(tài)下的單元結(jié)點(diǎn)力：

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

以上過程的具體步驟為：

對(duì)于第i增量步、第j迭代步的結(jié)構(gòu)增量平衡方程寫作[15]：

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

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

jλ為：

5 算例與分析

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

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

該算例取自文獻(xiàn)[24]，為驗(yàn)證幾何非線性分析單元及方法的經(jīng)典算例。如圖11，懸臂梁自由端受集中荷載，考慮平面內(nèi)彎曲，彈性模量E=10 MPa，泊松比ν=0.3，懸臂梁長(zhǎng)∶高∶厚=1000 mm∶100 mm∶10 mm。采用本文所建三角形平面膜單元分析，共劃分20 個(gè)單元，如圖12 所示。此外，還采用ANSYS 中的SHELL63 單元進(jìn)行對(duì)比，單元?jiǎng)澐中问酵瑘D12，SHELL63 單元每個(gè)結(jié)點(diǎn)6 個(gè)自由度。僅考慮膜的面內(nèi)剛度，所得懸臂梁端部荷載與豎向位移關(guān)系曲線如圖13 所示。從中可以看出，本文推導(dǎo)的三角形平面膜單元與 ANSYS 中的SHELL63 單元計(jì)算結(jié)果高度吻合，但本文所建單元 每個(gè)結(jié)點(diǎn)僅含兩個(gè)自由度，且?guī)缀蝿偠染仃嚍榫€性舉證，構(gòu)造簡(jiǎn)單。

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

圖12 三角形平面膜單元?jiǎng)澐址桨?Fig.12 Triangular plane membrane element

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

算例2. 空間膜單元算例—內(nèi)壓下的雙拋物線膜片

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

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

圖14 膜片初始未變形狀態(tài) Fig.14 Initial undeformed state of the membrane

圖15 內(nèi)壓為1.5 MPa 時(shí)膜片變形狀態(tài) Fig.15 Deformation state under 1.5 MPa of the membrane

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

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

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

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

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

圖17 初始未變形狀態(tài) Fig.17 Initial undeformed state of the membrane

圖18 內(nèi)壓為3.0 MPa 時(shí)變形狀態(tài) Fig.18 Deformation state under 3.0 MPa of the membrane

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

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

6 結(jié)論

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

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

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