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

2020-08-28 02:29:48曹進軍董凱駿彭福軍惲衛東

工程力學 2020年8期

張亮，曹進軍，董凱駿，彭福軍，惲衛東

(1. 重慶大學航空航天學院工程力學系，非均質材料力學重慶市重點實驗室，重慶 400044；2. 西安交通大學航天航空學院機械結構強度與振動國家重點實驗室，西安 710049；3. 大連理工大學工程力學系，工業裝備結構分析國家重點實驗室，大連 116024；4. 上海宇航系統工程研究所，上海 201109)

褶皺變形是薄膜結構的一種常見的失穩模式，其本質是結構受到局部壓應力作用而發生的屈曲和后屈曲行為。過去30 年，這一問題引起了國內外學者廣泛的研究興趣，研究領域涉及航天工程[1 ? 2]、生物組織工程[3]、微機電系統[4]和軟材料[5 ? 6]等。

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

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

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

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

1 共旋空間膜單元

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

中間局部坐標系 Cr與整體坐標系Cg之間的坐標轉換矩陣為：

單元初始、中間局部坐標系與整體坐標系之間的轉換關系如下：

對上述方程進行變分計算[25]，可得：

式中：下標e 表示單元；P 為消除剛體運動影響的投影矩陣；T 為分塊對角方陣，將整體位移矢量轉換為局部位移矢量。具體表達式如下：

其中：

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

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

其中：

得到單元切線剛度矩陣后，組集得到結構切線剛度矩陣Kt，進而求解有限元增量平衡方程：

式中：Fext和Fint分別為結構等效節點外力和內力向量；Ug為結構在整體坐標系下的節點位移向量。

2 褶皺模型及互補列式

本節基于雙模量材料的應力-應變關系，提出一個消除材料壓縮抗力的褶皺模型。雙模量材料具有拉、壓不對稱的力學行為[15,27]，即在拉伸和壓縮狀態下擁有不同的彈性模量。這意味著，可以通過設置材料的壓縮模量為零，達到消除材料壓縮抗力的目的。本節中，關于雙模量材料的方程推導在單元局部坐標系下進行，這正好借用了共旋坐標法的優勢，即：在單元局部坐標系下，小應變本構關系和線性有限元列式仍然成立。

在局部坐標系Cr中，雙模材料的應力-應變關系建立在主應力空間中，其平面應力問題的本構方程寫作：

其中：

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

對于本構方程，式(33)引入一個非負的參變量(2 階)張量 χ?，即：

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

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

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

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

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

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

其中：

3 算法執行

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

圖3 褶皺判據Fig.3 Wrinkling criterion

4 數值算例

4.1 方形氣囊

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

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

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

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

4.2 十字型氣囊

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