膜材是一种得到广泛使用的柔性材料, 仅有很小的抗压或抗弯刚度, 面外荷载作用下容易产生大变位大转动。碰撞接触是膜结构分析中的一类重要问题, 往往同时涉及大变形引起的几何、材料非线性以及接触界面的非线性, 具有极强的非线性效应。

膜材的碰撞接触问题包括膜材与刚体、膜材和膜材的碰撞两大类, 主要百富策略白菜网在建筑膜结构、空间充气可展结构、安全气囊和布料运动仿真等领域。文献[1]利用LS-DYNA对具有初速度的刚性小球撞击薄膜结构动力过程进行数值模拟。文献[2]采用气体动力学和多刚体动力学对卷曲折叠管进行了展开动力学研究, 展开过程中存在膜材自身的碰撞接触。文献[3]基于控制体积法采用LS-DYNA软件对空间充气天线的整体模型进行展开仿真分析。文献[4, 5]对汽车安全气囊的充气展开及与假人模型的碰撞过程进行了模拟分析及研究。文献[6]采用LS-DYNA对一种新型着陆缓冲气囊进行了仿真分析, 考虑了气囊与着陆面、气囊和气囊两类接触行为。文献[7, 8]对布料运动中的碰撞及自碰撞、接触和摩擦问题进行了研究并提出相应处理方法。文献[9]基于碰撞历史和精确摩擦, 采用弹簧质点模型对高分辨率的布料折叠和褶皱现象进行了数值仿真模拟。

不同于传统有限元的分析模式, 向量式有限元[10,11,12]是以向量力学理论为分析基础并基于点值描述来获得结构体系行为的新型分析方法。质点的运动方程求解采用逐点逐步循环计算, 不需求解非线性联立方程组和刚度矩阵, 即不存在迭代不收敛和刚度矩阵奇异问题。通过引入逆向运动和变形坐标系, 可消除刚体位移带来的数值误差。因此, 向量式有限元非常适合于结构大变位大转动、碰撞和断裂等复杂行为问题的分析。

本文首先给出了三角形膜单元的向量式有限元基本理论和求解过程, 然后针对膜材与刚体、膜材与膜材两类碰撞接触问题, 提出了碰撞检测和碰撞响应的处理方法。在此基础上编制了向量式有限元膜单元的碰撞接触分析程序, 通过算例分析验证了理论推导和编制程序的可靠性和计算稳定性。

向量式有限元是将结构离散为有质量的质点和质点间无质量的单元, 本文所用为常应变三角形单元。质点a的运动满足质点平动微分方程 (α=0时, 即为牛顿第二定律) :

式中ma为质点a的质量, 为质点的速度 (加速度) 向量, α为阻尼参数, Fa=faext+faint为质点受到的合力向量, 其中faext为质点的外力向量, faint是膜单元传递给质点的内力向量。

利用中央差分公式数值求解运动方程式 (1) 时, 需将其转化为差分形式, 无初始条件 (连续) 和有初始条件 (不连续) 时的差分式如式 (2) 所示 (省略下标a) :

式中Δt为时间步长, ;发生碰撞时刻前后属于不连续行为, 需使用式 (2b) 。

向量式有限元三角形膜单元的求解过程包括单元节点位移计算、单元节点内力计算和应力转换计算三个部分, 以下给出主要求解过程, 详细推导过程见文献[13]。

向量式有限元采用逆向运动方法从单元节点全位移中扣除刚体位移, 以获得单元节点的纯变形位移量。刚体位移包括刚体平移和刚体转动, 前者可以三角形单元的一个节点为参考点取其总位移向量;后者包括平面外和平面内转动位移, 可通过单元平面外 (内) 刚体转动向量和对应平面外 (内) 逆向转动矩阵换算获得。

在获得单元节点的纯变形位移量后, 通过单元变形满足的虚功方程来求解单元节点内力。首先引入单元变形坐标系, 将空间膜元问题转化到平面上;接着采用传统有限元的常应变三角形单元形函数分布, 获得变形坐标系下单元节点内力向量;然后, 利用坐标转换矩阵和平面外 (内) 正向转动矩阵, 将变形坐标系下单元节点内力向量转换回到整体坐标系下分量;最后, 反向作用于质点上并进行集成得到三角形膜单元传给质点的内力向量。

由于应力增量是通过虚拟位置时变形坐标系下的纯变形位移量得到的, 对于每一个时间子步均有不同的变形坐标系。因而上一时刻末变形坐标系下的应力需要利用坐标转换矩阵和平面外 (内) 正向转动矩阵, 首先转换回到整体坐标系下分量, 然后再转换到下一时刻初的变形坐标系下, 以实现应力的逐步循环计算。

膜材的碰撞接触行为包括膜材与刚体、膜材与膜材的碰撞两种情况, 本文研究这两种情况下膜材的碰撞检测和碰撞响应的处理, 不涉及碰撞检测中的优化加速问题。

刚体是指刚度相对膜材大得多或其变形可忽略的物体, 本节方法用于膜材与刚体之间碰撞接触的处理。

膜材与刚体的碰撞检测有刚体边界条件法和刚体通用网格法, 前者仅适用于几何形状比较简单的特定刚体 (如刚性球、刚平面等) ;后者, 则适用于任意几何边界形状的刚体。本文采用刚体通用网格法进行膜材碰撞检测处理, 其原理是将刚体表面进行三角形网格划分, 即将刚体表面离散为三角形平面单元。

本算法仅考虑膜材质点与刚体三角形网格面之间的单向碰撞检测, 即“点-三角形”检测[7,14]。算法中刚体网格实际所用变量仅为各三角形单元的节点坐标和法线向量。

若第i个三角形刚体单元的三个节点坐标分别记为XR_1、XR_2和XR_3, 则第i刚体单元的单位法线向量为

膜材在第n时间子步的运动过程中, 质点j的坐标位置从Xn运动到Xn+1, 则n和n+1时刻质点j到刚体单元i的距离 (含符号) 为

对应的距离向量为dn=dnn0和dn+1=dn+1n0。

若sign (dn·dn+1) ≤0, 则质点j与第i单元所在的平面αi发生碰撞并穿透。以下还需判断该碰撞点是否位于第i单元内部, 才能说明质点j与第i单元是否实际发生了碰撞。

由于时间子步比较短, 质点的位移变化量相对刚体单元为小量值, 因而可用n时刻质点j在第i单元上的垂足点Xn_p来代替实际碰撞点, 有

利用下式计算质点j的垂足点Xn_p在第i单元上的面积坐标L:

式中L={L1L2L3}T

若满足0≤Li≤1 (i=1, 2, 3) , 则垂足点Xn_p在第i刚体单元内部, 即质点j与第i刚体三角形单元实际发生碰撞, 需进行碰撞响应处理。

对于膜材质点j, 与其同时发生碰撞的刚体单元数r可为1、2或2以上。当r=1时, 即质点j的碰撞点仅落在1个刚体单元内;当r≥2时, 即质点j的碰撞点可能落在2个刚体单元交界线或多个单元交界节点上。处理碰撞响应时, 质点j的碰撞单位法向量可取所有发生碰撞单元法向量的平均值, 即

式中nj0为膜材质点j的碰撞单位法向量, n0 (i) 为与质点j发生碰撞的第i刚体单元的单位法向量。

当检测到膜材质点与刚体发生碰撞时, 需要对质点进行碰撞响应处理。常用的碰撞处理方法有初始条件更新法和罚函数法[15,16], 后者由于逻辑简单易于实现, 在碰撞领域获得了广泛百富策略白菜网。本文采用罚函数法进行膜材碰撞响应处理, 仅考虑光滑碰撞。

罚函数法的原理是在发生穿透的从动节点和主动面之间引入一个较大的界面接触力, 以限制从动节点对主动面的穿透。光滑碰撞时在物理上相当于放置一法向弹簧, 法向接触力 (即罚函数值) 的大小与穿透深度、主动片的刚度成正比。此方法无需引入碰撞和释放条件, 逻辑简单易于实现, 在接触-碰撞计算中已广泛百富策略白菜网;其缺点是罚参数较难选取, 目前尚没有明确的规则。

本文结合罚函数法与中央差分位移式, 提出基于中央差分式的罚接触力响应方法;同时, 赋予罚参数的选取规则, 以实现膜材与刚体的碰撞响应处理。

若在第n时间子步运动中, 质点j与刚体单元i (i=1, 2, …, r) 发生了碰撞穿透, 如图1所示。

质点j的碰撞单位法向量为nj0, 法向穿透向量δ=δnj0, 则附加法向罚接触力取为FN0=-kδnj0, 其中k是刚度因子 (即罚参数) 。质点j碰撞后初始时间子步内运动, 需采用含初始条件 (不连续) 的中央差分位移式 (2b) 进行计算。由式 (2b) 可知, 附加法向罚接触力FN0引起的质点j的附加位移为

因而, 当满足时, 即满足,

才能使质点j与刚体单元不会发生穿透。

此时罚参数为k=β (2m/Δt2) , 通过选取β的值来确定罚参数k的大小, 对于非弹性碰撞可取β=1, 弹性碰撞则可取较大的β值来计算k值。

本节方法用于膜材与膜材之间碰撞接触的处理, 对于同一膜材不同区域之间存在相互碰撞的行为, 也可通过膜材分区定义相互接触来实现。

同样, 本算法仅考虑“点-三角形”检测, 即膜材A的质点与膜材B的三角形网格面之间的单向碰撞检测。和膜材与刚体的碰撞问题不同, 由于被碰撞膜材B会发生变形, 故膜材B单元节点和法向向量均会发生改变。在第n时间子步的运动过程中, 膜材B第i单元的三个节点坐标位置从XBn_k (k=1, 2, 3) 运动到XBn+1_k, 则n和n+1时刻膜材B第i单元的单位法线向量n0n和n0n+1可由式 (3) 中的XR_k分别替换为XBn_k和XBn+1_k来计算获得。

膜材A质点j的坐标位置从XnA运动到XAn+1, 则n和n+1时刻膜材A质点j到膜材B第i单元的距离 (含符号) 为

对应的距离向量为dn=dnnn0和dn+1=dn+1n0n+1。

若sign (dn·dn+1) ≤0, 则膜材A质点j与膜材B第i单元所在的平面αi发生碰撞并穿透;同样, 还需判断该碰撞点是否位于膜材B的第i单元内部。

由于时间步长比较短, 质点的位移变化量相对于单元尺寸为小量值, 在计算碰撞点位置和附加罚接触力时均可以n+1时刻膜材B第i单元位置作为计算参考位置。

取n时刻膜材A质点j在n+1时刻膜材B第i单元位置上的垂足点Xn_p来代替实际碰撞点, 有

式中dnp= (XnA-XBn+1_1) ·n0n+1是n时刻膜材A质点j到n+1时刻膜材B第i单元的距离。

将式 (6) 中XR_k替换为XBn+1_k, 可计算得到膜材A质点j的垂足点Xn_p在膜材B第i单元上面积坐标L。

若满足0≤Li≤1 (i=1, 2, 3) , 则垂足点Xn_p在膜材B第i单元内部, 即膜材A质点j与膜材B第i单元实际发生碰撞, 需进行碰撞响应处理。

对于膜材A质点j, 与其同时发生碰撞的膜材B单元数r可为1、2或2以上。处理碰撞响应时, 膜材A质点j的碰撞单位法向量可由式 (7) 的n0 (i) 替换为n0n+1 (i) 来计算获得。

同膜材与刚体的碰撞响应处理, 采用罚函数法进行膜材与膜材的碰撞响应处理。

若在第n时间子步运动中, 膜材A质点j与膜材B第i单元 (i=1, 2, …, r) 发生了碰撞穿透如图1所示。此时, 膜材B第i单元取n+1时刻的位置作为参考位置来计算附加罚接触力。

膜材A质点j的附加法向罚接触力为FN0=-kδnj0。膜材A质点j碰撞后初始时间子步内运动, 需采用中央差分位移式 (2b) 进行计算。附加位移δ0和β满足的条件同式 (8, 9) 。

对于膜材与膜材的碰撞, 对应罚接触力FN0的膜材B上总附加罚接触反力-FN0;先将-FN0均分到膜材B第i单元 (i=1, 2, …, r) 的碰撞点 (垂足点) Xn_p (i) 处, 即

再由等效分配原则, 利用碰撞点Xn_p (i) 在膜材B第i单元上面积坐标L (i) ={L1 (i) L2 (i) L3 (i) }T, 将F′N0分配到膜材B第i单元的三个节点上, 即

式中k=1, 2, 3为膜材B第i单元的三个节点。最后进行集成以获得质点上的附加罚接触力。

基于膜单元基本理论及碰撞检测与碰撞响应的处理方法, 本文采用Matlab编制了向量式有限元膜单元的碰撞接触分析程序。下面进行算例验证。

以图2所示的方形膜片覆盖刚性球问题为例进行分析。膜片边长L=0.5 m, 厚度ta=1.0e-3m, 初始状态为水平放置且边界无约束, 刚性球半径R=0.125m, 位于膜片正下方;膜片仅在重力作用下做下落运动, 与刚性球发生碰撞接触并覆盖刚性球。膜材弹性模量E=7.0e-4Pa, 泊松比ν=0.3, 质量密度ρ=27kg/m3。采用三角形单元分别对膜片和刚性球进行网格结构化划分, 膜片和刚性球分别有1800个、1200个三角形薄膜单元和961个、602个节点如图2所示, 分析时其时间步长取h=2.0e-4s, 阻尼参数取α=20。分析中不考虑膜片自身的碰撞接触行为。

图3给出了膜片在重力作用下覆盖刚球过程中几个典型时刻的变形图。可以看出, 随着膜片在重力作用下的下落运动, 首先与刚性球发生碰撞接触, 在t=0.2s时, 沿对角线开始出现较为明显的折痕;之后, 未碰撞的膜片外围区域继续下落而呈现大变形 (如t=0.4s时) , 在t=0.6s时, 下落变形已扩展至整个膜片区域, 并出现了大量的细小折痕效应;在t=1.2s时, 膜片四个角端已下落至最低位置;之后, 由于动力振荡效应膜片继续发生细微的摆动;在t=1.6s时, 达到最终覆盖刚性球的平衡状态。本算例体现了本文膜材与刚体碰撞接触算法的可靠性和计算稳定性。

如图4所示为四边形折叠气囊充气模型, 由上、下两片四边形薄膜封闭并对折而成如图4 (b) 所示, 边长L=0.5 m, 厚度ta=1.0e-3m, 初始状态为水平放置且边界无约束, 仅在一连续均匀填充气压q作用下做充气膨胀变形, 膨胀时存在膜材与膜材的碰撞接触行为。气压q采用斜坡-平台方式缓慢施加, 在t0=0.3s时达到最大值q0=60Pa。材料弹性模量E=7.0e-4Pa, 泊松比ν=0.3, 质量密度ρ=27kg/m3。采用三角形单元进行网格结构化划分, 共1600个三角形薄膜单元和802个节点 (图4) , 分析时其时间步长取h=2.0e-4s, 阻尼参数取α=100。

由于本算例考察的是膜材自身不同区域间的碰撞接触行为, 需将折叠气囊分成上、下两折叠部分区域, 并定义相互碰撞接触算法, 以实现其充气展开过程的模拟。图5给出了折叠气囊充气过程中几个典型时刻的变形图。可以看出, 随着气压q的增大, 折叠气囊在t=0.05s前上、下两折叠部分各自逐渐膨胀变形直至开始出现相互碰撞接触;在t=0.1s时, 由于碰撞接触约束的作用上、下两折叠部分呈现向两边展开变形;之后, 展开变形继续增大并扩展至整个气囊而呈现大变形 (如t=0.15s时) ;在t=0.2s时气囊大部分区域已展开而仅在折边处尚存在局部折痕;最后, 由于气囊内压力达到饱和状态;在t=0.3s时, 气囊已完全展开, 并在周围边界上产生了皱折效应。计算结果符合折叠气囊充气展开的实际变形过程。本算例体现了本文膜材与膜材碰撞接触算法的可靠性和计算稳定性。

将向量式有限元引入膜材的碰撞接触分析, 针对膜材与刚体、膜材与膜材两类碰撞接触问题, 采用“点-三角形”检测和基于中央差分式的罚接触力方法处理膜材的碰撞检测和碰撞响应问题, 编制了向量式有限元膜单元的碰撞接触分析计算程序。

通过膜片覆盖刚球、折叠气囊展开等典型算例考察了本文方法在膜材碰撞接触问题中的百富策略白菜网, 跟踪获得从接触出现、接触大范围产生直至处于最终位置状态的全过程, 验证了理论推导和编制程序的可靠性和计算稳定性;同时, 体现出向量式有限元方法进行膜材碰撞接触分析的优势。