充氣薄膜褶皺分析的高效互補(bǔ)有限元列式

張 亮，曹進(jìn)軍，董凱駿，彭福軍，惲衛(wèi)東

(1. 重慶大學(xué)航空航天學(xué)院工程力學(xué)系，非均質(zhì)材料力學(xué)重慶市重點(diǎn)實(shí)驗(yàn)室，重慶 400044；2. 西安交通大學(xué)航天航空學(xué)院機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710049；3. 大連理工大學(xué)工程力學(xué)系，工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024；4. 上海宇航系統(tǒng)工程研究所，上海 201109)

褶皺變形是薄膜結(jié)構(gòu)的一種常見的失穩(wěn)模式，其本質(zhì)是結(jié)構(gòu)受到局部壓應(yīng)力作用而發(fā)生的屈曲和后屈曲行為。過去30 年，這一問題引起了國(guó)內(nèi)外學(xué)者廣泛的研究興趣，研究領(lǐng)域涉及航天工程[1 ? 2]、生物組織工程[3]、微機(jī)電系統(tǒng)[4]和軟材料[5 ? 6]等。

除實(shí)驗(yàn)和理論研究外，作為工程結(jié)構(gòu)分析的一般性手段，褶皺變形問題的數(shù)值模擬也已大量涌現(xiàn)。已有的薄膜褶皺分析數(shù)值模擬方法可分為兩大類：1) 薄殼有限元后屈曲分析[7 ? 12]；2) 基于張力場(chǎng)理論的材料修正方法[13 ? 21]。這兩種方法各有優(yōu)、缺點(diǎn)。

在薄殼后屈曲分析中，材料既具有面內(nèi)的拉伸剛度，又具有面外的彎曲剛度，可以定量地獲取褶皺的細(xì)節(jié)信息，如波長(zhǎng)和幅值。然而，數(shù)值模擬結(jié)果對(duì)薄殼單元的尺寸、類型，以及初始缺陷等參數(shù)都有著很強(qiáng)的依賴性[7,12]?；趶埩?chǎng)理論的材料修正方法是研究薄膜褶皺行為的一種簡(jiǎn)化辦法，它假設(shè)材料的壓縮和彎曲抗彎剛度均可忽略，薄膜一旦處于壓應(yīng)力狀態(tài)便會(huì)產(chǎn)生褶皺或松弛[22]。國(guó)內(nèi)外學(xué)者在這方面做了大量的研究工作。Miller 和Hedgepeth[13]提出了一種修正材料彈性矩陣的方法來(lái)表征薄膜單元所處的應(yīng)力狀態(tài)，并以此預(yù)測(cè)結(jié)構(gòu)中的張緊、褶皺和松弛區(qū)域。Ding 和Yang[14]指出許多應(yīng)力迭代方法本質(zhì)上是“試錯(cuò)法”(try-error)，對(duì)于存在松弛區(qū)域的薄膜，不同迭代方法得到的結(jié)果可能差異很大。Ding和Yang[14]推導(dǎo)了2-VP 模型的參變量變分原理，該方法計(jì)算效率高，并且剛度矩陣不依賴于薄膜應(yīng)力狀態(tài)的迭代而更新；Zhang 等[15]提出了一種用于雙模材料非線性分析的二次規(guī)劃算法，并通過消除材料的壓縮抗力來(lái)預(yù)測(cè)薄膜的褶皺區(qū)域；Ding 和Yang[14]、Zhang 等[15]提出的方法具有良好的算法穩(wěn)定性，但都局限于幾何小變形情況；Zhang等[16]后來(lái)將其提出的方法擴(kuò)展到幾何大變形，但有限元列式仍局限于二維平面問題。

對(duì)于充氣薄膜結(jié)構(gòu)，尤其是具有較大剛體位移和嚴(yán)重褶皺變形的情況，直接對(duì)材料進(jìn)行修正以消除其壓縮抗力的辦法往往會(huì)導(dǎo)致算法無(wú)法收斂。其原因在于Newton 類方法只能在平衡構(gòu)型附近進(jìn)行有效的迭代求解。然而，一方面，在充氣的初始階段，薄膜結(jié)構(gòu)的承載能力非常有限，平衡構(gòu)型遠(yuǎn)離初始構(gòu)型，給迭代求解帶來(lái)障礙；另一方面，材料修正容易引起內(nèi)力的突變和局部區(qū)域剛度矩陣的奇異，進(jìn)而造成算法振蕩，甚至發(fā)散。一個(gè)欠約束的充氣氣囊就是很好的例子，它被很多學(xué)者作為標(biāo)準(zhǔn)測(cè)試算例來(lái)研究。例如，Contri 和Schrefler[17]提出了一個(gè)“兩步求解法”有限元程序來(lái)處理充氣氣囊的褶皺問題，其本質(zhì)思想類似于線搜索技術(shù)，擴(kuò)大了Newton 迭代的收斂半徑，這種求解方法的經(jīng)驗(yàn)性很強(qiáng)，對(duì)于不同的問題往往需要反復(fù)的嘗試；Lee 和Youn[18]采用“擬動(dòng)態(tài)法”對(duì)Newton 迭代進(jìn)行初始化，即先求解一個(gè)虛假的動(dòng)態(tài)平衡問題，得到一個(gè)近似平衡構(gòu)型，在此基礎(chǔ)上，再進(jìn)行常規(guī)的增量迭代求解。需要指出的是，“擬動(dòng)態(tài)法”本質(zhì)上屬于顯式動(dòng)力學(xué)方法，而Newton 迭代是隱式求解方法，二者的迭代格式截然不同。在兩種不同的方案之間切換，將增加程序?qū)崿F(xiàn)的復(fù)雜程度；Jarasjarungkiat 等[19]類比塑性問題，引入了“褶皺應(yīng)變”的概念，并對(duì)變形梯度張量進(jìn)行了修正，以消除薄膜的壓縮抗力。文中指出，這種修正不可避免地導(dǎo)致由于應(yīng)力重分布而引起的算法振蕩。算法收斂性可以通過引入“懲罰因子”的辦法得到改善[19]。然而，“懲罰因子”的確定具有一定的經(jīng)驗(yàn)性。

綜上所述，鑒于充氣薄膜結(jié)構(gòu)的強(qiáng)非線性特征和現(xiàn)有求解方法的不足，針對(duì)空間充氣薄膜結(jié)構(gòu)褶皺分析的高效數(shù)值方法仍然值得研究。本文針對(duì)大變形(大位移、小應(yīng)變)充氣薄膜，提出了一種能夠準(zhǔn)確預(yù)測(cè)結(jié)構(gòu)位移、應(yīng)力和褶皺區(qū)域的互補(bǔ)有限元方法。首先，基于共旋坐標(biāo)法，將薄膜的大變形分解為整體坐標(biāo)系下的剛體運(yùn)動(dòng)和局部坐標(biāo)系下的小應(yīng)變變形；其次，在單元局部坐標(biāo)系下，基于張力場(chǎng)理論構(gòu)造了一個(gè)褶皺模型及相應(yīng)的線性互補(bǔ)問題，用于計(jì)算單元節(jié)點(diǎn)內(nèi)力。由于在迭代求解過程中，并不需要依據(jù)薄膜所處的狀態(tài)(張緊、褶皺或松弛)來(lái)更新單元的材料剛度矩陣，該方法能夠有效地消除迭代求解過程中的內(nèi)力振蕩，具有良好的收斂性和穩(wěn)定性。

1 共旋空間膜單元

圖1 空間三節(jié)點(diǎn)膜單元的大變形描述Fig.1 Finite deformation of a 3-node spatial membrane element

中間局部坐標(biāo)系 Cr與整體坐標(biāo)系Cg之間的坐標(biāo)轉(zhuǎn)換矩陣為：

單元初始、中間局部坐標(biāo)系與整體坐標(biāo)系之間的轉(zhuǎn)換關(guān)系如下：

對(duì)上述方程進(jìn)行變分計(jì)算[25]，可得：

式中：下標(biāo)e 表示單元；P 為消除剛體運(yùn)動(dòng)影響的投影矩陣；T 為分塊對(duì)角方陣，將整體位移矢量轉(zhuǎn)換為局部位移矢量。具體表達(dá)式如下：

其中：

根據(jù)虛功原理的恒等性，并考慮式(14)，得到：

式中，Ke為局部材料剛度矩陣(亦即線性剛度矩陣)。將式(14)和式(25)代入式(24)，化簡(jiǎn)后得到單元一致切線剛度矩陣：

其中：

得到單元切線剛度矩陣后，組集得到結(jié)構(gòu)切線剛度矩陣Kt，進(jìn)而求解有限元增量平衡方程：

式中：Fext和Fint分別為結(jié)構(gòu)等效節(jié)點(diǎn)外力和內(nèi)力向量；Ug為結(jié)構(gòu)在整體坐標(biāo)系下的節(jié)點(diǎn)位移向量。

2 褶皺模型及互補(bǔ)列式

本節(jié)基于雙模量材料的應(yīng)力-應(yīng)變關(guān)系，提出一個(gè)消除材料壓縮抗力的褶皺模型。雙模量材料具有拉、壓不對(duì)稱的力學(xué)行為[15,27]，即在拉伸和壓縮狀態(tài)下?lián)碛胁煌膹椥阅Ａ?。這意味著，可以通過設(shè)置材料的壓縮模量為零，達(dá)到消除材料壓縮抗力的目的。本節(jié)中，關(guān)于雙模量材料的方程推導(dǎo)在單元局部坐標(biāo)系下進(jìn)行，這正好借用了共旋坐標(biāo)法的優(yōu)勢(shì)，即：在單元局部坐標(biāo)系下，小應(yīng)變本構(gòu)關(guān)系和線性有限元列式仍然成立。

在局部坐標(biāo)系Cr中，雙模材料的應(yīng)力-應(yīng)變關(guān)系建立在主應(yīng)力空間中，其平面應(yīng)力問題的本構(gòu)方程寫作：

其中：

式中： ε?和 σ?分別表示主應(yīng)變張量和主應(yīng)力張量； C+和 C?分別為主應(yīng)力空間中雙軸拉伸和雙軸壓縮的柔度張量。由式(33)不難發(fā)現(xiàn)，雙模量材料的本構(gòu)關(guān)系是非線性的，其彈性常數(shù)由主應(yīng)力狀態(tài)確定。已有的研究表明：雙模量問題的非線性迭代求解存在收斂困難，其具體原因已在文獻(xiàn)[27]中詳細(xì)分析。本文通過構(gòu)造一個(gè)線性互補(bǔ)問題來(lái)克服數(shù)值分析的收斂困難。

對(duì)于本構(gòu)方程，式(33)引入一個(gè)非負(fù)的參變量(2 階)張量 χ?，即：

為了保持與式(32)的等價(jià)性，式(36)必須附加一個(gè)約束條件，即：

式中，為簡(jiǎn)便起見，這里只列出該約束條件。參變量 χ?的引入及其物理意義的說明請(qǐng)見Zhang等[15]較早的研究工作。通過在式(37)中引入一個(gè)非負(fù)松弛變量(2 階)張量 μ?，可以將約束條件轉(zhuǎn)化為標(biāo)準(zhǔn)的線性互補(bǔ)問題，即：

式(33)與式(36)、式(38)之間的等價(jià)性可以得到證明[15]。

在主應(yīng)力空間中，雙模材料的應(yīng)變能密度可以表示為：

進(jìn)一步，可在單元局部坐標(biāo)系下建立雙模量材料問題的參變量最小勢(shì)能原理[28]：

式中：ud代表局部坐標(biāo)系下的位移；b 和p 分別代表體力和面力。對(duì)結(jié)構(gòu)進(jìn)行有限元離散，并對(duì)式(40)進(jìn)行變分計(jì)算[28]，便得到局部坐標(biāo)系下的單元平衡方程和互補(bǔ)方程：

其中：

3 算法執(zhí)行

圖2 有限元程序流程圖Fig.2 Flow chart of finite element procedure

圖3 褶皺判據(jù)Fig.3 Wrinkling criterion

4 數(shù)值算例

4.1 方形氣囊

首先分析一個(gè)方形氣囊的充氣膨脹變形。作為標(biāo)準(zhǔn)測(cè)試，該算例已被很多學(xué)者采用[17 ? 20]。鑒于結(jié)構(gòu)的對(duì)稱性，選取單層膜進(jìn)行分析。如圖4 所示，對(duì)角線AC 長(zhǎng)為1.2 m；薄膜厚度為6×10?4m。材料的彈性模量為58.8 MPa；泊松比為0.4。邊界條件：四條邊約束面外位移，面內(nèi)自由；水平中心線上約束在y 向位移，豎直中心線上約束x 方向位移。沿x 軸、y 軸方向的位移分別用u、v 表示；垂直于xy 平面向外的位移用w 表示；A 點(diǎn)沿對(duì)角線MA 方向的位移用rA表示。

分別采用128 個(gè)、200 個(gè)、512 個(gè)和800 個(gè)三角形膜單元對(duì)圖4 所示的結(jié)構(gòu)進(jìn)行有限元離散。內(nèi)壓p=5000 Pa以增量的形式逐步施加。表1 中列出了A、B 和M 點(diǎn)的位移，以及M 點(diǎn)的最大主應(yīng)力結(jié)果?？梢园l(fā)現(xiàn)，隨著網(wǎng)格數(shù)目的增加，數(shù)值結(jié)果逐漸收斂。其中，本文結(jié)果與同樣采用三角形單元的文獻(xiàn)[19]的結(jié)果吻合得最好，從而驗(yàn)證了本文方法的正確性。圖5 呈現(xiàn)了方形薄膜充氣膨脹后的變形構(gòu)型。從圖5(a)可以看出，薄膜的四條邊向內(nèi)收縮明顯，因而結(jié)構(gòu)在面內(nèi)是欠約束的。這意味著在充氣初始階段，結(jié)構(gòu)將發(fā)生很大的剛體位移，其平衡位置不易找到，這也是該算例收斂困難的原因之一[17]。圖6(a)繪制出了內(nèi)壓p=5000 Pa時(shí)，數(shù)值模擬得到的薄膜褶皺區(qū)域。其中，灰色代表“張緊”狀態(tài)；深色代表“褶皺”狀態(tài)。數(shù)值模擬結(jié)果與圖6(b)所示的實(shí)驗(yàn)結(jié)果基本吻合。

圖4 方形氣囊示意圖Fig.4 Sketch of a square airbag

該算例的數(shù)值分析分為三個(gè)載荷步，即：I. 面內(nèi)拉伸；II. 充氣膨脹；III. 材料修正，釋放面力拉應(yīng)力。其中，第I 載荷步和第II 載荷步內(nèi)均只采用1 個(gè)增量載荷步，第III 載荷步內(nèi)采用40 個(gè)增量載荷步。圖7 給出了程序執(zhí)行過程中的收斂誤差曲線?？梢钥吹?，在第I 載荷步、第II 載荷步內(nèi)，以及第III 載荷步的開始階段，程序只需1 次～2 次迭代便可達(dá)到收斂容差；而在第III 載荷步的中、后階段，則需要更多的迭代次數(shù)。這是因?yàn)樵诘贗II 載荷步的中、后階段，程序中的材料修正(消除壓縮應(yīng)力)開始生效，材料和幾何非線性同時(shí)伴隨其中，使得問題的非線性更強(qiáng)。從圖7中不難看出，收斂誤差整體上是逐步減小的，不存在上下跳躍的現(xiàn)象，這說明算法具有較好的穩(wěn)定性。表2 列出了算法收斂所需的迭代次數(shù)。算法所需的總迭代次數(shù)為133 次；當(dāng)材料修正被激活后(第III 步)，算法收斂所需的平均迭代次數(shù)僅為3 次；單個(gè)載荷增量步內(nèi)所需迭代次數(shù)最多為18 次。與“擬動(dòng)態(tài)法”[18]相比，本文方法具有更好的收斂性。

4.2 十字型氣囊

圖9 繪制了M、B 兩點(diǎn)的面外位移，以及A 點(diǎn)的y 向位移隨內(nèi)壓的變化關(guān)系曲線。其中，B 點(diǎn)和M 點(diǎn)的面外位移結(jié)果與文獻(xiàn)的結(jié)果吻合良好，而在充氣的初始階段，A 點(diǎn)的y 向位移則存在一定的誤差[20]。這是由于本文與參考文獻(xiàn)采用了不同的單元類型和網(wǎng)格密度。在“十字型”交叉處，薄膜單元嚴(yán)重收縮，甚至可能發(fā)生局部重疊。該問題具有很強(qiáng)的網(wǎng)格依賴性。在充氣初期，薄膜變形以向內(nèi)收縮的剛體位移為主(應(yīng)變能很小)，數(shù)值結(jié)果對(duì)網(wǎng)格的依賴性更為明顯。從力學(xué)本質(zhì)上講，向內(nèi)收縮的剛體位移可能存在多解的情況，采用不同的單元類型、不同的網(wǎng)格密度可能得到不同的位移解。隨著內(nèi)壓逐漸增大，薄膜的應(yīng)變能會(huì)逐漸增大，剛體位移所占成分會(huì)逐漸減小，由不同網(wǎng)格造成的位移誤差也會(huì)隨之減小。內(nèi)壓為p=2000 Pa 時(shí)的三維變形如圖10 所示。可見，薄膜四邊在面內(nèi)收縮，并在中心交叉位置出現(xiàn)了局部重疊。圖11 呈現(xiàn)了數(shù)值模擬所得到的褶皺分布區(qū)域。圖11(a)表明：內(nèi)壓p=2000 Pa 時(shí)，褶皺出現(xiàn)在兩條交叉氣柱的端部和中心交叉區(qū)域；當(dāng)內(nèi)壓增大為p=150 kPa 時(shí)，氣囊進(jìn)一步膨脹至完全張緊狀態(tài)，褶皺區(qū)域消失，如圖11(b)所示。本文模擬結(jié)果與文獻(xiàn)所呈現(xiàn)的現(xiàn)象相一致[20]。

表1 方形氣囊模擬結(jié)果比較Table1 Comparison of simulation results for the square airbag

圖5 p=5000 Pa 時(shí)的變形Fig.5 The deformed shape under p=5000 Pa

圖6 方形氣囊的褶皺區(qū)域Fig.6 Wrinkling regions of the square airbag results of simulation

圖7 方形氣囊的收斂曲線Fig.7 Convergence curve for the square air bag

表2 方形氣囊的迭代次數(shù)Table2 Iterations for the square air bag

圖8 十字型氣囊示意圖Fig.8 Sketch of a cross shaped airbag

圖9 A、B 和M 點(diǎn)位移結(jié)果Fig.9 Results of displacement at the points A, B and M

圖10 p=2000 Pa時(shí)的三維變形Fig.10 The three-dimensional deformed shape under p=2000 Pa

圖11 十字型氣囊的褶皺區(qū)域Fig.11 Wrinkling regions of the cross shaped airbag

5 結(jié)論

基于張力場(chǎng)理論，提出了一種適用于充氣薄膜結(jié)構(gòu)褶皺分析的互補(bǔ)有限元方法。通過兩個(gè)典型的數(shù)值算例，驗(yàn)證了方法的正確性。該方法能夠準(zhǔn)確地預(yù)測(cè)充氣薄膜結(jié)構(gòu)的位移、應(yīng)力水平，以及褶皺區(qū)域。具體得到以下兩點(diǎn)結(jié)論：

(1) 基于共旋坐標(biāo)方法，推導(dǎo)了一個(gè)空間三節(jié)點(diǎn)三角形膜單元的切線剛度矩陣，可用于一般充氣結(jié)構(gòu)的大位移分析。

(2) 在單元局部坐標(biāo)系下構(gòu)造的線性互補(bǔ)列式有效地消除了迭代求解過程中應(yīng)力重分布導(dǎo)致的算法振蕩。算法具有良好的收斂性和穩(wěn)定性。

在本文方法的基礎(chǔ)上，可進(jìn)一步考慮材料的彎曲剛度，以準(zhǔn)確獲取褶皺變形的三維形貌。