充气薄膜褶皱分析的高效互补有限元列式

2020-08-28曹进军董凯骏彭福军恽卫东

工程力学 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]。