冷气式气体发生器爆破片极限破坏载荷的计算

单津晖，魏学哲，王 承

（1. 同济大学汽车学院，上海 201804；2. 均胜汽车安全系统研发（上海）有限公司CAE部，上海 201707）

帘幕式安全气囊可以在汽车发生侧面碰撞或翻滚时和安全带一起挽救司乘人员的生命，是汽车被动安全系统的重要组成部分。冷气式气体发生器将高压冷气存储在发生器腔体内，点火后火药爆炸将起阀门作用的爆破片炸开，高压冷气快速释放并充满气袋。充气过程中没有对人体有害的气体产生，且高压冷气充满气袋后吸热膨胀，对气袋内部气体压力保持十分有利，因此对气体保压时间有要求的防翻滚帘幕气囊多选用冷气式气体发生器。

冷气式气体发生器内存储的是高压缩比的混合气体，爆破片是气体唯一的泄放装置，在未点火时爆破片必须有足够的强度，而在点火时爆破片必须在设计破坏压力载荷下按预设的破坏模式发生破坏，其破坏载荷及破坏方式都有严格的设计要求［1-2］。

在设计开发过程中，为了达到尽可能高的预设压力并加快验证速度，通常会采用水压试验获取爆破片的极限载荷，即通过灌注高压水直至爆破片破坏，试验过程中记录的最高水压值即为极限压力破坏载荷。根据实验室数据，准静态水压和发生器动态点爆加载条件下的爆破片破坏载荷基本一致，因此在设计开发阶段对水压试验爆破片极限载荷的准确预测是发生器结构件设计的关键。

1 极限载荷的分析方法

确定极限载荷的方法主要分为试验法、理论计算法和有限元法三大类，试验法只能用于后期验证，在设计早期主要采用理论计算法和有限元法。当前针对爆破片的极限压力破坏载荷相关计算理论并未统一，理论方法对材料非线性及几何非线性的考虑较少，且无法考虑零件发生破坏前的局部失稳，因此理论计算结果精度尚不能满足工程开发需要。目前基于非线性弹塑性理论的有限元计算正逐渐成为获取极限载荷的主要手段［3-5］。

在有限元分析中，破坏载荷极值的判定方法较多，如塑性功法、2 倍斜率法、双切线相交法、零曲率法、弧长法等［6-11］。其他方法相比零曲率法、弧长法精度不足，在实际中应用也不易掌握，本文也不再讨论。

1.1 零曲率法

极限载荷分析属于塑性力学问题，假定材料为理想弹塑性，结构处于小变形状态，服从Mises 屈服条件及相关流动准则。在上述假定基础上，在有限元分析中，通过多个载荷步、小的载荷增量逐步对结构施加载荷。当载荷增量逐渐趋近于零、即载荷对应位移曲线的曲率趋近于零时，此时结构的承载能力到达极限，再增加载荷会引起极限载荷下的结构失稳，即发生塑性垮塌。这种将载荷位移曲线上零曲率点对应的载荷值定义为极限载荷的方法称为零曲率法。

从力学分析角度看，极限载荷问题求解，就是通过不断求解计入几何非线性和材料非线性的结构平衡方程寻求结构极限荷载的过程。非线性结构的增量平衡方程一般可以表示为［8-9］

式中：K为结构整体切线刚度；ΔX为结构位移增量向量；Δλ为载荷比例系数；F为载荷向量；R为不平衡力向量。

在极限载荷状态下，载荷比例系数Δλ趋近于零，当计算尝试给Δλ一个系统设定的最小增量时，整体刚度矩阵K出现奇异，刚度矩阵变成“病态”矩阵，即系数矩阵的微小误差将导致计算结果的巨大误差，引起计算不收敛而导致计算终止。在有限元计算中，零曲率法将计算最后一步能收敛的载荷定义为极限载荷。

零曲率法基于隐式算法求解，计算位移增量ΔX需要计算结构刚度K的逆矩阵，因此如果结构非线性程度较高，尤其是含复杂接触的大载荷极限工况，该方法收敛难度大。另外零曲率法的计算终止是以载荷增量步为零作为标志，属于计算非正常终止，这要求分析人员有较高的能力对结果做出判断。

1.2 弧长法

弧长法是一种隐式非线性求解的迭代控制方法，不需要事先设定失稳准则或断裂准则，可以在荷载和位移增量均不确定的情况下，追踪结构加载路径，避免传统隐式算法的收敛问题［12-13］。弧长增量法以载荷比例系数和位移增量平方和的开方作为增量参数进行收敛计算分析，其增量参数为

在收敛计算过程中，弧长法以前一步增量计算得到的平衡点作为圆心，基于超平面搜索、求解半径为ΔS的超曲面S与结构平衡路径的交点。在达到极限载荷后，增量弧长法在超曲面S内可沿载荷增量Δλ负曲率方向搜索平衡状态，通过产生负的载荷增量Δλ以匹配结构失稳后负定的刚度矩阵K，继续追踪结构位移ΔX增加对应载荷F的变化。

弧长法将载荷增量为正的最后一个迭代步对应的载荷定义为极限载荷。其极限载荷求解精度与零曲率法基本相当，但相比零曲率法，弧长法可以继续求解结构在失稳状态下刚度负定、载荷增量为负的力学平衡方程。这对研究整个结构在达到极限载荷后应力、应变变化规律具有重要意义。

对每一个计算迭代步，弧长法都要对结构刚度矩阵正定或负定进行判定，并对载荷增量变化在正负曲率变化方向进行搜索，因此相比只考虑载荷增量正曲率变化的零曲率法，弧长法的计算增量步长较小，每一迭代步的计算时间更长。

1.3 针对显式算法的应力矩阵状态分析法

基于隐式算法的零曲率法和弧长法对计算收敛的要求较高，对需要考虑接触的复杂结构，在极限大载荷加载状态下计算通常都很难收敛。对于瞬态高度非线性工况，通常采用显式算法求解。t时刻结构的动力学微分方程可表示为［14］

式中：M为质量矩阵；C为阻尼矩阵；Üt为加速度向量；U̇t为速度向量；Ut为位移向量。

若t时刻节点位移、速度加速度均为已知，使用中心差分法对加速度、速度的导数采用中心差分替代，式（3）可转化为［15］

由式（4）、（5）可知，基于中心差分法的动力学微分方程组求解不需要像式（1）那样需要求解结构刚度K的逆矩阵；另外中心差分法将加速度Ü和速度U̇表示位移U的某种组合，求解t+Δt时刻的瞬时位移向量Ut+Δt时，只需要根据t+Δt时刻以前的状态变量计算出M̑和Ȓ，即可直接计算出Ut+Δt。因此只要设定好计算时间步长Δt，显式算法一般不存在收敛性问题，可有效解决各种复杂非线性问题。

显式算法基于中心差分动力学方程（4），不是力学平衡方程（1）求解，对于极限载荷计算，在结构达到承载极限后，只要有满足收敛的时间步长Δt显式计算就会继续，施加的载荷也会按预先设定继续增加，不会像隐式算法那样出现计算终止或载荷下降的情况。因此对于基于显式算法的极限载荷求解，工程上多采用固定约束位置支撑反力曲线是否出现拐点来判定极限载荷：当结构达到承载极限后发生失稳，其固定约束位置的支撑反力不再随载荷增加而增加，而因结构失去承载能力开始下降，将固定约束位置支撑力反力曲线出现的极值拐点值定义为结构的极限载荷，这种极限载荷定义方法的原理与弧长法类似，计算精度也能满足工程需要。但对于载荷施加位置靠近固定约束位置的工况，就算结构失稳支撑反力也会随载荷的增加而继续增加，这给极限载荷的判定带来了较大难度，也表明该判定方法有相当的局限性。

结构刚度矩阵K由结构所有单元的刚度组合而成，类似式（1）单元刚度矩阵Ke是表征单元节点力向量及节点位移向量的关系矩阵。在结构几何确定的情况下，节点刚度矩阵主要取决于节点应力矩阵σe的状态。

对三维结构上任一节点，其应力矩阵σe为一含9个应力分量的3×3 对称矩阵，为了方便分析，通过矩阵线性变换将σe转化为只含有3个应力分量的特征值矩阵，转化后的σe可表示为［16］

式中：σ1、σ2和σ3分别为节点第一、第二及第三主应力。3 个主应力σ1、σ2和σ3的取值及其方向决定了一点的应力状态。

应力矩阵σe中各应力分量大小取决于材料应力-应变（σ-ε）本构关系。塑性材料σ-ε曲线可以表示为

式中：k为材料应变硬化系数；n为材料应变硬化指数；σT为材料拉伸极限应力；εT为拉伸极限应变。图1为爆破片材料应力-应变曲线，在达到极限应力前，σ-ε按图中曲线实线段变化，在实线末端应力达到极值、材料刚度最大后，若继续增加载荷，对应不同算法：①隐式算法（零曲率法）。σ-ε曲线刚度达到极值，计算因没有大于极限载荷的收敛解而终止。②隐式算法（弧长法）。载荷增量变为负值，应变随位移增加而增加，σ-ε曲线按图1 虚线段变化。③显式算法。载荷继续增加，应变随位移增加继续增加，σ-ε曲线按图1虚线段变化。

对比后2种算法，在达到极限载荷后，由于材料本构相同，其应力-应变对应结构位移变化的趋势接近，关键位置的应力矩阵变化趋势也接近，使用应力矩阵状态分析法判定显式算法的极限载荷，即通过分析、对比隐式弧长法和显式算法2 种算法相同节点在极限状态下的应力矩阵变化规律来判定显式算法的极限载荷。

2 爆破片极限载荷分析

2.1 基于零曲率法计算爆破片极限载荷

使用ABAQUS 软件对爆破片的极限载荷进行计算。爆破片材料牌号为AISI301，材料屈服极限1 100MPa，真实拉伸极限1 482MPa，真实断裂应变0.182，其材料曲线如图1所示，如图中实线部分，先按公式σ=1 711ε0.0831对应变化，σ达到拉伸极限后不再增加，σ-ε曲线按图中虚线段变化。

图1 AISI301材料应力-应变曲线Fig.1 Strain versus stress of AISI301

发生器前端示意图如图2 所示，爆破片为正拱型爆破片，爆破片通过激光焊接与扩散器连接，扩散器和发生器壳体通过摩擦焊焊接成一个整体。爆破片与扩散器焊接后局部截面图如图3所示。爆破片厚度0.203 2±0.012 7mm，仿真取厚度下限0.190 5mm。在实际试验中，除爆破片外，其他零件变形很小，因此为简化仿真，仿真中只考虑爆破片零件。为保证模型精度，对爆破片模型使用六面体实体网格建模。

图2 冷气式气体发生器出气端结构示意Fig.2 Structure of outlet of cold gas inflator

图3 焊接后爆破片区域截面示意Fig.3 Section view of burst disk after welding

为了提高仿真计算精度，在极限载荷计算中计入冲压工艺对爆破片的影响，仿真分2步完成。第1步使用动态显式分析模块将处于自由状态的爆破片圆片冲压成正拱型。冲压前后爆破片的厚度如图4所示。冲压后爆破片厚度沿四周向中心逐渐减薄，中心点位置厚度0.168 4mm，减薄率13.1%。含冲压仿真产生的预应力、厚度不一的爆破片仿真模型是第2 步仿真的起始条件，与无预应力且厚度均一的爆破片模型相比无疑更符合真实情况。

图4 冲压前后爆破片截面形状对比Fig.4 Section of burst disk before and after stamping

第2步仿真在爆破片周边焊接位置添加固定约束，并在爆破片沿内表面法向施加100Mpa 的压力（设计破坏压力大于82.6Mpa）。使用静态隐式非线性分析模块计算爆破片的极限载荷。在第21 个计算迭代步，载荷为83.5MPa 时仿真计算终止，计算过程中爆破片中心N133号节点位移始终最高，其位移对应压力载荷曲线如图5 的虚线所示。由图5 可知，在载荷较低时，载荷增量增加较快，当载荷趋近于83.5Mpa 时，载荷增量逐渐减小并趋近于零。在载荷为83.5Mpa、N133 号节点位移为0.868mm 时载荷增量为零，此时曲线的斜率为零。根据零曲率法准则可以判定爆破片的极限破坏载荷为83.5Mpa。

图5 压力载荷对应爆破片最大变形量变化曲线Fig.5 Curves of pressure versus maximum displacement of burst disk

根据试验室数据，15 次水压试验的平均破坏载荷为86.0MPa，最大值88.4Mpa，最小值83.1MPa。仿真计算获得的极限载荷和试验破坏均值载荷误差为3.0%，与试验最小值误差为0.5%。考虑到仿真中爆破片厚度取公差下限，这个极限载荷计算结果精度十分理想。

试验后爆破片破坏形状如图6a 所示，其破坏从中心位置开始；极限载荷作用下爆破片第一主应变云图如图6b 所示，爆破片中心位置第一主应变最大。对比试验结果，仿真对应变最大位置的预测是准确的。

图6 爆破片试验及仿真结果对比Fig.6 Test and simulation results of burst disk

2.2 基于弧长法计算爆破片极限载荷

2.2.1 使用弧长法计算爆破片极限载荷

使用弧长法计算爆破片的极限载荷，第1 步冲压成形仿真和零曲率法一致，第2 步隐式非线性计算仿真除了零曲率算法已经设定的几何非线性选项外，再增加RIKS（弧长法）算法选项。弧长法基于小载荷增量步计算，相比零曲率法其迭代步数增加较多，对于大极限载荷计算，需要预先设定一个较大的最大迭代步数，本算例设定值为300步。

在弧长法计算过程中，在第183 迭代步载荷达到峰值83.6Mpa，第184 步后载荷增量为负值，载荷也随位移的增加而减小。爆破片中心N133 号节点位移对应压力载荷曲线如图5 实线所示。图5 实线在载荷为83.6MPa 时出现拐点，压力达到峰值，此时对应节点位移0.903mm，在此位移之前结构一直加载，曲线曲率为正值，大于此位移后结构失稳、可承受的最大载荷也相应下降，曲线曲率也变为负值。

将图5实线的极值83.6Mpa定义为弧长法计算获得的极限载荷，其值和试验破坏载荷均值的误差为2.9%，最小值误差为0.6%。其极限载荷计算结果和零曲率法基本一致，只有0.1%的差别。

2.2.2 弧长法爆破片极限载荷计算结果分析

对比分析零曲率法、弧长法计算结果，对图5的2条载荷位移曲线进行对比分析：

（1）在达到极限载荷前，2条曲线走势基本重合，零曲率法载荷增量步较大，直到接近极限载荷时载荷增量才逐渐减小。

（2）从第1个迭代步开始，弧长法每一个迭代步的载荷增量都很小，其曲线变化趋势完全体现出材料弹塑性本构曲线变化和载荷变化的对应关系。

（3）由图5实线可知，在极限载荷作用下结构失稳后，弧长法可以继续计算对应位移增加爆破片结构的应力-应变，这对研究达到极限载荷后爆破片应力-应变变化规律非常有意义。

由式（2）弧长法搜索、计算极限载荷的原理可知，弧长法不依赖于某个强度失效理论，而是基于结构刚度矩阵正定或负定判断结构是否失稳。因此在进行应力-应变结果分析时，要对表征刚度的应力矩阵σe进行分析，即对式（6）中σ1、σ2和σ3及其对应应变变化规律进行考量。

在仿真过程中，爆破片中心位置变形一直最大，对应应力、应变也最大。对爆破片中心N133号节点进行应力-应变分析，其受力状态为两拉一压，沿2个径向拉伸、沿厚度方向压缩。

N133号节点第一、第二主应力对应位移变化曲线如图7 所示，在位移为0.6mm 后，第一、第二主应力增长趋势变得非常缓慢，但始终低于材料极限应力。N133 号节点第三主应力绝对值对应位移变化曲线如图8 所示，在位移为0.903mm 时第三主应力取极大值，后继计算其值随位移增加而减小。

图7 N133号节点第一、第二主应力对应的位移变化曲线Fig.7 Displacement versus 1st and 2nd principle stress of Node 133

图8 N133号节点第三主应力对应的位移变化曲线Fig.8 Displacement versus 3rd principle stress of Node 133

对N133 号节点的应力、应变分析表明，达到极限载荷时，该节点第一、第二主应力接近材料拉伸极限，其位移对应第三主应力及压力载荷2 条曲线同时产生拐点，这说明该节点，即爆破片中心位置沿厚度方向发生了失稳，整个结构的承载能力开始下降，爆破片将从中心位置为起点发生塑性破坏。

N133号节点3个方向主应变对应的位移变化曲线如图9所示。可见，第一、第二主应变变化曲线基本重合，第三主应变绝对值为第一、第二主应变之和；在极限载荷对应位移0.903mm位置，3个主应变变化都比较平缓，没有明显的拐点，这和图1中应变从实线段平缓过渡到虚线段的趋势相匹配。

图9 N133 号节点第一、第二、第三主应变对应的位移变化曲线Fig.9 Displacement versus 1st, 2nd, and 3rd principle strain of Node 133

2.3 基于应力矩阵状态分析法计算爆破片极限载荷

对于本算例，通过零曲率法和弧长法都获得了精度较高的极限载荷值。但受计算方法限制，仿真中未考虑爆破片与下方焊接零件扩散器的接触，由图3 可知，与扩散器的接触限制了爆破片对应位置向下的移动，若忽略此接触，爆破片的变形结果精度会有一定失真。

为考察接触设定对爆破片极限载荷和变形结果的影响，进一步提高仿真精度，使用DYNA 显式分析软件对爆破片的极限载荷进行计算。分析过程同样分为2 步，第1 步为冲压仿真，第2 步在爆破片内表面沿法向施加一个0～100MPa随时间线性增加的压力，同时约束爆破片焊接位置并设置爆破片与扩散器零件的接触。接触类型为面对面接触（surfaceto-surface），由于爆破片与扩散器都是钢材，材料刚度接近，对应接触柔度设定为类型1（Soft设为1）。

第2步仿真计算在载荷达到100MPa后按预先设定终止，在载荷为88.3MPa 时爆破片中心位置单元因变形过大被删除，但在整个仿真过程中，爆破片焊接位置附近区域始终与其下方的扩散器接触并承受不断增加的法向压力，导致爆破片约束位置的支撑反力曲线始终保持持续增加的趋势，没有极值拐点出现。若将仿真中爆破片中心位置出现断裂时对应的内表面法向压力值88.3MPa确定为爆破片的极限载荷，其值与试验结果上限值相当，和试验下限结果83.1MPa误差为6.3%，计算精度低于隐式算法。

使用应力矩阵状态分析法对显式算法的极限载荷进行判定。与弧长法分析过程相同，对爆破片中心N133号节点应力矩阵状态变化进行分析，其3个方向主应力对应位移变化曲线与弧长法对应结果基本一致。限于篇幅，只对出现拐点的第三主应力变化曲线进行讨论。

显式算法计算获得的N133 号节点第三主应力绝对值对应压力载荷及位移变化曲线分别如图10、图11 所示。由图10 可知，当载荷为84.0MPa 时，N133号节点第三主应力达到峰值，并随载荷增加迅速下降。由图11 可知，当N133 节点第三主应力绝对值在节点位移为0.590mm时达到峰值，然后随位移增加其值逐渐减小。

图10 N133号节点第三主应力的应力-位移曲线Fig.10 Pressure versus 3rd principle stress of Node 133

对比图8、图11弧长法和显式算法N133号节点第三主应力变化曲线，可知：

图11 N133号节点第三主应力对应其位移变化曲线Fig.11 Displacement versus 3rd principle stress of Node 133

（1）弧长法和显式算法计算结果曲线都出现极值拐点，表明在厚度方向爆破片中心位置出现了失稳，爆破片结构也将在此位置发生破坏。

（2）在极值点前，2条曲线变化趋势基本相同，在极值点后由于载荷继续增加显式算法曲线下降趋势相比载荷开始下降的弧长法要更快一些。即相比弧长法，显式算法在厚度方向的失稳随位移增加会更迅速。

将显式计算过程中爆破片中心位置失稳即N133 号节点第三主应力取极值时对应的载荷84.0MPa 定义为爆破片的极限破坏载荷，其值与爆破片试验破坏载荷均值误差为2.4%，与试验结果下限误差为1.1%。

上述基于隐式弧长法和显式中心差分算法计算结果对爆破片中心位置关键节点应力矩阵变化规律进行分析、对比，确定显式算法极限载荷的方法即为应力矩阵状态分析法。对比隐式弧长法与针对显式算法的应力矩阵状态分析法计算结果，可知：

（1）使用应力矩阵状态分析法和弧长法判定的极限载荷分别为84.0MPa和83.6MP，2种方法的匹配度达到99.5%，表明应力矩阵状态分析法和弧长法的判定的极限载荷值有极佳的一致性。

（2）显式计算中考虑爆破片下方零件支撑的影响后，极限载荷作用下爆破片最大变形结果由0.902mm 减小到0.590mm，相差34.6%，因此若要准确计算爆破片的变形量，必须使用基于显式算法的应力矩阵状态分析法。

3 结论

首先使用基于隐式算法的零曲率法和弧长法进行计算，确定冷气式气体发生器爆破片水压试验的极限破坏压力载荷，考虑了爆破片的冲压成形工艺后，和试验结果相对比，这2种方法对爆破片极限载荷的计算误差都在3%以内。

对弧长法极限载荷计算过程中爆破片中心节点应力矩阵状态变化规律进行分析，发现在极限载荷作用下，中心位置节点第三主应力对应位移曲线出现极值拐点，该拐点对应载荷步与压力位移曲线极值拐点对应载荷步相同。这表明在达到极限载荷后，爆破片中心位置在厚度方向发生了失稳。

为了弥补隐式算法接触非线性难收敛的不足，基于弧长法应力-应变分析结果提出针对显式算法的极限载荷判定方法，即应力矩阵状态分析法，它将显式计算过程中爆破片中心在厚度方向发生失稳时对应的载荷定义为极限载荷。该方法在利用显式算法在解决接触非线性易收敛的优势同时，避免了结构局部微小失稳时显式算法受力状态不完全平衡带来的误差。最终使用应力矩阵状态分析法判定的爆破片极限载荷值与弧长法结果匹配度达到99.5%，与试验结果相比误差在2.4%以内，且相比隐式算法该方法的变形结果更加准确。

作者贡献声明：

单津晖：理论研究、提炼，论文撰写。

魏学哲：论文统筹、规划。

王 承：仿真，试验协调。