褶皱变形是薄膜结构的一种常见的失稳模式，其本质是结构受到局部压应力作用而发生的屈曲和后屈曲行为。过去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   共旋空间膜单元

工程中常见的充气薄膜结构往往发生大位移、小应变变形，如太阳帆、安全气囊等。在已有的非线性有限元方法中，共旋坐标法是一种适用于大位移、小应变问题的简化方法^[23-25]。共旋坐标法将物体的大变形分解为一个整体坐标系下的刚体运动和一个局部小应变变形。因此，在局部坐标系下，已有的小变形本构关系和线性有限单元均可以直接使用。对于杆、梁、板、壳等结构的大位移分析，共旋坐标法使用起来尤其方便^[25-26]。这里，以空间三节点膜单元为例来说明基于共旋坐标法的大变形运动学描述。如图1所示，共旋框架中有三种构型：初始构型、中间构型和当前构型。同时，建立三个坐标系，即：整体坐标系$C_{}^{\rm{g}}$，以$({{e}}_1^{\rm{g}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_2^{\rm{g}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_3^{\rm{g}})$表示；单元初始局部坐标系$C_{}^0$，以$({{e}}_1^0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_2^0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_3^0)$表示；单元中间局部坐标系$C_{}^{\rm{r}}$，以$({{e}}_1^{\rm{r}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_2^{\rm{r}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_3^{\rm{r}})$表示。单元局部坐标系的原点设置在1号节点。初始局部坐标系$({{e}}_1^0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_2^0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_3^0)$定义为：

图  1  空间三节点膜单元的大变形描述

Figure  1.  Finite deformation of a 3-node spatial membrane element

下载: 全尺寸图片 幻灯片

式中，${{X}}_i^0{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (i = 1,2,3)$表示三个节点在整体坐标系$C_{}^{\rm{g}}$中的坐标。初始局部坐标系$C_{}^0$与整体坐标系$C_{}^{\rm{g}}$之间的坐标转换矩阵为：

类似的，建立在中间构型上的局部坐标系$({{e}}_1^{\rm{r}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_2^{\rm{r}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {{e}}_3^{\rm{r}})$定义为：

中间局部坐标系$C_{}^{\rm{r}}$与整体坐标系$C_{}^{\rm{g}}$之间的坐标转换矩阵为：

以单元内任意一点p为例，整体位移${{U}}_p^{\rm{g}}$和局部纯变形位移${{u}}_p^{\rm{d}}$可分别表示为：

单元初始、中间局部坐标系与整体坐标系之间的转换关系如下：

式中，${{x}}_1^{\rm{g}}$和${{X}}_1^{\rm{g}}$分别为1号节点在整体坐标系中的当前构型和初始构型上的坐标位置矢量。根据式(11)和式(12)，单元在局部坐标系下的纯变形位移可表示为：

对上述方程进行变分计算^[25]，可得：

式中：下标e表示单元；P为消除刚体运动影响的投影矩阵；T为分块对角方阵，将整体位移矢量转换为局部位移矢量。具体表达式如下：

其中：

式中：“spin(x)”表示作用于向量x上的反对称旋转矩阵；${{x}}_{ei}^{\rm{r}}(i = 1,2,3)$表示建立在当前构型上，单元e的第i个节点在局部坐标系$C_{}^{\rm{r}}$中的坐标。

根据虚功原理的恒等性，并考虑式(14)，得到：

式中，${{F}}_e^{\rm{g}}$和${{f}}_e^{\rm{d}}$分别为整体和局部坐标系下的单元节点力向量，之间的转换关系为：

通过式(22)对${{U}}_e^{\rm{g}}$的变分计算可以得到单元一致切线刚度矩阵：

式(24)右端第一项中有${{f}}_e^{\rm{d}} = {{K}}_e {{u}}_e^{\rm{d}}{\kern 1pt}$，因此，变分$\delta {{f}}_e^{\rm{d}}$可写作：

式中，K_e为局部材料刚度矩阵(亦即线性刚度矩阵)。将式(14)和式(25)代入式(24)，化简后得到单元一致切线刚度矩阵：

式中，${{K}}_e^{\rm{m}}$代表材料刚度矩阵，即：

几何刚度矩阵${{K}}_e^{\rm{gr}}$、${{K}}_e^{\rm{gp}}$分别被定义为旋转、平衡投影刚度矩阵，即：

其中：

得到单元切线刚度矩阵后，组集得到结构切线刚度矩阵K^t，进而求解有限元增量平衡方程：

式中：${{{F}}_{{\rm{ext}}}}$和${{{F}}_{{\rm{int}}}}$分别为结构等效节点外力和内力向量；U^g为结构在整体坐标系下的节点位移向量。

2   褶皱模型及互补列式

本节基于双模量材料的应力-应变关系，提出一个消除材料压缩抗力的褶皱模型。双模量材料具有拉、压不对称的力学行为^{[15, 27]}，即在拉伸和压缩状态下拥有不同的弹性模量。这意味着，可以通过设置材料的压缩模量为零，达到消除材料压缩抗力的目的。本节中，关于双模量材料的方程推导在单元局部坐标系下进行，这正好借用了共旋坐标法的优势，即：在单元局部坐标系下，小应变本构关系和线性有限元列式仍然成立。

在局部坐标系$C_{}^{\rm{r}}$中，双模材料的应力-应变关系建立在主应力空间中，其平面应力问题的本构方程写作：

其中：

式中：${{\varepsilon}}_{}^ *$和${{\sigma}}_{}^ *$分别表示主应变张量和主应力张量；${{C}}_{}^ +$和${{C}}_{}^ -$分别为主应力空间中双轴拉伸和双轴压缩的柔度张量。由式(33)不难发现，双模量材料的本构关系是非线性的，其弹性常数由主应力状态确定。已有的研究表明：双模量问题的非线性迭代求解存在收敛困难，其具体原因已在文献[27]中详细分析。本文通过构造一个线性互补问题来克服数值分析的收敛困难。

对于本构方程，式(33)引入一个非负的参变量(2阶)张量${{\;\chi}} _{}^ *$，即：

为了保持与式(32)的等价性，式(36)必须附加一个约束条件，即：

式中，为简便起见，这里只列出该约束条件。参变量${{\\text{χ}} ^*}$的引入及其物理意义的说明请见Zhang等^[15]较早的研究工作。通过在式(37)中引入一个非负松弛变量(2阶)张量${\mu ^{\rm{*}}}$，可以将约束条件转化为标准的线性互补问题，即：

式(33)与式(36)、式(38)之间的等价性可以得到证明^[15]。

在主应力空间中，双模材料的应变能密度可以表示为：

进一步，可在单元局部坐标系下建立双模量材料问题的参变量最小势能原理^[28]：

式中：u^d代表局部坐标系下的位移；b和p分别代表体力和面力。对结构进行有限元离散，并对式(40)进行变分计算^[28]，便得到局部坐标系下的单元平衡方程和互补方程：

其中：

式中：符号“{ }”表示张量的向量形式；而符号“[ ]”表示张量的矩阵形式；I为$2 \times 2$的单位矩阵，B为应变矩阵。通过求解标准线性互补问题，式(43)得到参变量$\{ {\\text{χ}}_e^ * \}$。然后，由式(42)计算局部坐标系下的单元节点内力${{f}}_e^{\rm{d}}$。进而，由式(23)转换得到整体坐标系下的单元节点内力${{F}}_e^{\rm{g}}$。组集所有单元，便得到式(32)中的结构内力向量${{{F}}_{{\rm{int}}}}$。

3   算法执行

在MATLAB语言环境中进行有限元程序实现。程序执行流程如图2所示。本文所有数值算例均采用相对力误差作为收敛准则，即$\dfrac{{\left\| {{{{F}}_{{\rm{ext}}}} - {{{F}}_{{\rm{int}}}}} \right\|}}{{\left\| {{{{F}}_{{\rm{ext}}}}} \right\|}} < \psi$(取$\psi = {10^{ - 4}}$)。采用基于主应力的褶皱判据，如图3所示。

图  2  有限元程序流程图

Figure  2.  Flow chart of finite element procedure

下载: 全尺寸图片 幻灯片

图  3  褶皱判据

Figure  3.  Wrinkling criterion

下载: 全尺寸图片 幻灯片

程序执行过程中，在更新单元的切线刚度矩阵(见式(26))和单元节点内力(见式(42))时，并不需要根据单元所处的主应力状态(张紧、褶皱或松弛)来更新局部材料刚度矩阵K_e。单元所处的主应力状态由参变量${\\text{χ}}_e^ *$自然表征，使得算法能够很好地避免由于单元主应力状态突变可能造成的刚度矩阵奇异和算法振荡问题。这也是本文方法有别于其他材料修正方法的地方。

4   数值算例

4.1   方形气囊

首先分析一个方形气囊的充气膨胀变形。作为标准测试，该算例已被很多学者采用^[17-20]。鉴于结构的对称性，选取单层膜进行分析。如图4所示，对角线AC长为1.2 m；薄膜厚度为$6 \times {10^{ - 4}}$ m。材料的弹性模量为58.8 MPa；泊松比为0.4。边界条件：四条边约束面外位移，面内自由；水平中心线上约束在y向位移，竖直中心线上约束x方向位移。沿x轴、y轴方向的位移分别用u、v表示；垂直于xy平面向外的位移用w表示；A点沿对角线MA方向的位移用r_A表示。

图  4  方形气囊示意图

Figure  4.  Sketch of a square airbag

下载: 全尺寸图片 幻灯片

分别采用128个、200个、512个和800个三角形膜单元对图4所示的结构进行有限元离散。内压$p = 5000{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{Pa}}$以增量的形式逐步施加。表1中列出了A、B和M点的位移，以及M点的最大主应力结果。可以发现，随着网格数目的增加，数值结果逐渐收敛。其中，本文结果与同样采用三角形单元的文献[19]的结果吻合得最好，从而验证了本文方法的正确性。图5呈现了方形薄膜充气膨胀后的变形构型。从图5(a)可以看出，薄膜的四条边向内收缩明显，因而结构在面内是欠约束的。这意味着在充气初始阶段，结构将发生很大的刚体位移，其平衡位置不易找到，这也是该算例收敛困难的原因之一^[17]。图6(a)绘制出了内压$p = 5000{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{Pa}}$时，数值模拟得到的薄膜褶皱区域。其中，灰色代表“张紧”状态；深色代表“褶皱”状态。数值模拟结果与图6(b)所示的实验结果基本吻合。

表  1  方形气囊模拟结果比较

Table  1.  Comparison of simulation results for the square airbag

单元类型	四边形(1/4模型)[17]				四边形(1/4模型)[18]					三角形(1/4模型)[19]				三角形(全模型)(本文)
单元数目	16	25	64		16	64	400	32		50	128	200		128	200	512	800
wM/m	0.209	0.217	0.205		0.212	0.215	0.216	0.221		0.219	0.218	0.218		0.221	0.219	0.217	0.217
rA/m	0.057	0.063	0.047		0.097	0.074	0.069	0.051		0.050	0.050	0.049		0.051	0.050	0.049	0.049
vB/m	0.102	0.110	0.130		0.120	0.122	0.123	0.127		0.125	0.121	0.120		0.128	0.125	0.122	0.121
σM/MPa	3.400	3.500	3.500		3.600	3.500	3.700	3.600		3.900	3.900	3.900		3.260	3.670	3.880	3.910

下载: 导出CSV 

| 显示表格

图  5  p=5000 Pa时的变形

Figure  5.  The deformed shape under p=5000 Pa

下载: 全尺寸图片 幻灯片

图  6  方形气囊的褶皱区域

Figure  6.  Wrinkling regions of the square airbag results of simulation

下载: 全尺寸图片 幻灯片

该算例的数值分析分为三个载荷步，即：I. 面内拉伸；II. 充气膨胀；III. 材料修正，释放面力拉应力。其中，第I载荷步和第II载荷步内均只采用1个增量载荷步，第III载荷步内采用40个增量载荷步。图7给出了程序执行过程中的收敛误差曲线。可以看到，在第I载荷步、第II载荷步内，以及第III载荷步的开始阶段，程序只需1次~2次迭代便可达到收敛容差；而在第III载荷步的中、后阶段，则需要更多的迭代次数。这是因为在第III载荷步的中、后阶段，程序中的材料修正(消除压缩应力)开始生效，材料和几何非线性同时伴随其中，使得问题的非线性更强。从图7中不难看出，收敛误差整体上是逐步减小的，不存在上下跳跃的现象，这说明算法具有较好的稳定性。表2列出了算法收敛所需的迭代次数。算法所需的总迭代次数为133次；当材料修正被激活后(第III步)，算法收敛所需的平均迭代次数仅为3次；单个载荷增量步内所需迭代次数最多为18次。与“拟动态法”^[18]相比，本文方法具有更好的收敛性。

图  7  方形气囊的收敛曲线

Figure  7.  Convergence curve for the square air bag

下载: 全尺寸图片 幻灯片

表  2  方形气囊的迭代次数

Table  2.  Iterations for the square air bag

载荷步	载荷增量步数	总迭代次数	平均迭代次数	最大迭代次数
I	1	1	1	1
II	1	2	2	2
III	40	130	3	18

下载: 导出CSV 

| 显示表格

4.2   十字型气囊

研究一个十字型气囊的充气膨胀，仍然采用单层膜进行模拟分析。整个结构用5120个三角形单元进行有限元离散，总自由度为7971。初始平整的薄膜如图8所示，$L = 0.3\sqrt 3$ m。薄膜的厚度、弹性模量和泊松比均与上一个算例取值相同。边界条件：四周约束面外位移，面内自由；水平中心线上约束y向位移，竖直中心线上约束x向位移。

图  8  十字型气囊示意图

Figure  8.  Sketch of a cross shaped airbag

下载: 全尺寸图片 幻灯片

图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]。

图  9  A、B和M点位移结果

Figure  9.  Results of displacement at the points A, B and M

下载: 全尺寸图片 幻灯片

图  10  $p = 2000{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{Pa}}$时的三维变形

Figure  10.  The three-dimensional deformed shape under $p = 2000{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{Pa}}$

下载: 全尺寸图片 幻灯片

图  11  十字型气囊的褶皱区域

Figure  11.  Wrinkling regions of the cross shaped airbag

下载: 全尺寸图片 幻灯片

5   结论

基于张力场理论，提出了一种适用于充气薄膜结构褶皱分析的互补有限元方法。通过两个典型的数值算例，验证了方法的正确性。该方法能够准确地预测充气薄膜结构的位移、应力水平，以及褶皱区域。具体得到以下两点结论：

(1) 基于共旋坐标方法，推导了一个空间三节点三角形膜单元的切线刚度矩阵，可用于一般充气结构的大位移分析。

(2) 在单元局部坐标系下构造的线性互补列式有效地消除了迭代求解过程中应力重分布导致的算法振荡。算法具有良好的收敛性和稳定性。

在本文方法的基础上，可进一步考虑材料的弯曲刚度，以准确获取褶皱变形的三维形貌。