一种弯曲修正膜单元及对侧翻一步碰撞算法的改进

王 童

（长安大学汽车学院，陕西西安710064）

近年来公共交通的发展呈现迅猛发展态势，客车作为一种公共出行工具，逐步占据社会主流地位。同时，与客车相关的重大交通事故数量也逐步呈现上升趋势，给人类的生命财产带来巨大损失。侧翻是客车最严重的事故类型之一，往往造成群死群伤的恶劣后果［1］。目前，国内外许多学者已经投入到客车侧翻安全性的探索与改进研究中。比如，Keith等［2］发现高强钢顶梁或很便宜的纤维增强复合塑料(FRP）玻璃可以阻止盒状现象出现，并预防客车侧翻的发生.

客车侧翻一步碰撞算法，是笔者参考板料冲压一步成型算法思想［3-5］，提出的一种用于客车侧翻碰撞模拟的新算法。对于板料冲压一步成型算法，王鹏等［6］已经做过了一些研究。在他们的研究中，大部分金属成形过程都通过应用基于塑性变形理论的一步有限元方法进行模拟，它可以快速预测坯料最终构型及应力应变分布。借鉴该算法的核心思想，基于非线性全量理论，客车侧翻一步碰撞算法与现有的LS-DYNA 等增量法软件相比，在略微牺牲一点计算精度的条件下，计算效率可以得到大幅度提升，在车身设计初期，对客车的侧翻安全性能进行快速评价。

算法中单元模型的选择需兼顾精度和计算效率，需在两者间寻找一个平衡点［7］。侧翻过程结构的变形特性主要为弯曲变形，所有膜单元及较简单的三节点三角形壳单元均无法满足要求［8-10］。由于具有低阶、简单、位移连续等优点，初始侧翻一步碰撞算法选择了基于Mindlin 四边形板单元和四节点膜单元组合而成的空间壳单元对结构进行离散。而应用四节点膜单元进行模拟，结构的总自由度会比壳单元减少一半，其计算时间约为原时间的60%左右，对算法模拟更有优势。同时，膜单元弯曲效应的模拟问题也必须解决。

本文针对传统四节点膜单元无法进行弯曲的缺陷，基于单元内力和弯矩的平衡原理，提出一种可进行弯曲修正的新型四节点膜单元模型。将其应用于客车侧翻一步碰撞算法，简便有效地考虑了客车侧翻过程中的单元节点弯曲内力。因其计算量远小于传统壳单元，故在有效保证计算精度的同时，计算效率可以得到显著提升。

1 客车侧翻一步碰撞算法基本理论[11]

客车侧翻一步碰撞算法，基于非线性全量理论和比例加载假定，依据ECE R66 法规，忽略了中间状态和构形的变化，只考虑结构碰撞开始和最大变形两个状态。根据车体侧翻碰撞过程的运动变形特点和能量转换关系，得到满足变形条件的初始解，并采用Newton-Raphson法迭代求解，快速获得结构的最终变形。

将碰撞开始状态的结构作为原始构形X0。此时车体未发生变形，结构动能Ed由式（1）计算：

式中：M为车体质量；Δh为车体重心下降高度；J为车体绕假定转轴的转动惯量；ω为车体角速度。

碰撞开始状态结构各节点的速度v0由式（2）计算：

式中：ri为各节点到侧翻假定转轴距离；n为节点数。

在结构最大变形状态，车体结构产生明显变形。由于侧翻一步碰撞算法中，结构的最大变形是不确定的，因此需假定一个最大变形构型X0。此刻各节点的位移U0由式（3）计算：

式中：x0为结构变形最大构型。

此时车体动能Ed在碰撞中主要转换为结构形变能，结构形变能W由式（4）计算［12］：

式中：ε为单元应变；σ为单元Cauchy应力；Ve为单元体积；N为单元数。

判断此刻的结构形变能W与车体动能Ed是否满足式（5）的能量关系假定：

若不满足，则对节点位移U0进行修正，按照式（4）重新计算结构形变能。将满足式（5）能量关系假定的节点位移U作为Newton-Raphson 迭代初始解。

由于车体结构在空间内变形过程无外力作用，满足能量转换关系的初始解U，节点失衡力R(U)已处于不平衡状态：

应用Newton-Raphson 法，解决节点失衡力的不平衡问题。对初始解U按照式（7）迭代求解，使式（6）达到平衡，得到结构的最终变形。

对于式（7）中切线刚度矩阵的求解，由于需要迭代求解大型线性方程组，算法计算效率大大降低。故在初始算法中，课题组采用改进型拟共轭梯度法代替该切线刚度矩阵的计算，同时引入摄动原理，在无需计算切线刚度矩阵的情况下，对广义失衡力采用一维带宽存储，通过梯度因子快速寻找广义失衡力下降方向，使广义失衡力“2”范数达到极小值，结构碰撞能量达到平衡，获得结构的最终变形。

2 弯曲修正的四节点膜单元模型

本文提出的可进行弯曲修正的新型四节点膜单元模型如图 1 所示。图中i，j，k，l为单元节点在整体坐标系下的坐标；iu，ju，ku，lu为单元上表面 4 个节点在整体坐标系下的坐标；id，jd，kd，ld为单元下表面 4 个节点在整体坐标系下的坐标；ni，nj，nk，nl为单元各节点法线。与传统壳单元模型不同，这种弯曲修正的四节点膜单元节点内力由两部分构成：一部分是由面内变形引起的节点切向力，另一部分是由弯曲变形引起的节点弯矩和法向力，如式（8）所示：

式中：Fint为节点的总内力；Fm为因面内变形而引起的节点内力；Fb为因弯曲变形而引起的节点弯曲内力。Fm和Fb的向量形式由式（9）表示：

面内变形引起的节点内力Fm的计算与普通膜单元相同。为计算单元节点弯曲内力Fb，首先利用四节点膜单元沿节点法线等距的原理，计算弯曲后单元上下表面应力。图1所示为最大变形构形x0下的某一单元A(i，j，k，l)，单元厚度为t。该单元节点在 整 体 坐 标 系 下 的 坐 标 分 别 为i(xi，yi，zi)，j(xj，yj，zj)，k(xk，yk，zk)，l(xl，yl，zl)，各节点法线分别为ni(li，m，iqi)，nj(lj，m，jqj)，nk(lk，m k，qk)，nl(ll，m，lql)，节点法线由节点相邻单元加权平均求得。沿着单元A的4个节点法线方向及负方向平行等距移动半个板厚可得其上、下两个表面Au和Ad，同时分别计算其两个表面的应力。规定上表面Au四 个 节 点 的 坐 标 分 别 为i(xui，yui，zui)，j(xuj，yuj，zuj)，k(xuk，yuk，zuk)，l(xul，yul，zul)，下 表 面Ad四个节点的坐标分别为i(xdi，ydi，zdi)，j(xdj，ydj，zdj)，k(xdk，ydk，zdk)，l(xdl，ydl，zdl)。单元上、下表面节点处的坐标可由式（10）计算：

图1 可进行弯曲修正的四节点膜单元Fig.1 Four-node membrane element considering bending modification

这样，可得到单元上、下表面节点处的位移。以上表面单元Au为例，整体坐标系下各节点的位移可由式（11）计算：

式中：uur表示节点ru(r=i，j，k，l)在整体坐标系中的位移向量；xur表示节点ru(r=i，j，k，l)在最大变形构形x0下的整体坐标；Xur表示节点ru(r=i，j，k，l)在原始构形X0下对应节点的坐标。

根据大变形几何关系可以求得Green 变形张量下的单元主伸长λur(r=1，2，3)，在主应变空间的对数应变εur(r=1，2，3)可由式（12）计算：

考虑材料各向异性的塑性本构关系，在比例加载条件下，上表面单元Au任一点的应力值σu可由式（13）计算［12］：

这样，即可获得上表面单元Au任一点的应力值σu，同理，可获得下表面单元Ad任一点的应力值σd。

根据薄板弯曲问题的有限元理论，薄板的广义内力可由式（14）计算：

式中：Mx为垂直于x轴截面单位长度弯矩；My为垂直于y轴截面单位长度弯矩；Mxy为垂直于x（y）轴截面单位长度扭矩，Mxy=Myx。

由于沿z轴方向的应力成线性分布，薄板内任意一点的应力可由式（15）计算：

对于图1 中的模型，z=t/2。则虚拟板单元某一截面单位长度的弯矩和扭矩可由式（16）计算：

需要注意的是式（15）中，σx，σy，τxy为由板弯曲变形引起的三个应力分量，式（13）给出的是由面内变形和弯曲变形共同引起的应力。故需进一步将应力分量分解为面内变形应力σp和弯曲变形应力σb两部分。根据薄板弯曲理论，弯曲应力应该满足σub=-σdb。令Δ。

因单元内部广义内力相互平衡，故需在单元边界对式（16）进行积分。由弯曲变形引起的合弯矩可由式（17）计算：

式中：s为单元的边界。

考虑到碰撞结束后单元处于平衡状态，如图2所示，图中Fri，Frj，Frk，Frl，r=x，y，z为单元4个节点所受到的法向力，其作用在单元上的合外力为零。由单元平衡条件可得式（18）：

简化式（18）可得到式（19）：

由于式（19）中包含Fzi，Fzj，Fzk，Fzl四个未知参数，但只有三个方程，故求解该方程组仍需要一个补充条件。本文提出一种基于单元节点法向力方差的最小极值法，通过对单元4个节点法向力求方差，对所构成的函数求极小值。结合式（19）求解单元4个节点的法向力，所得各节点法向力的离散程度最小，计算结果最趋于稳定。补充条件可由式（20）表示：

将式（19）代入式（20），可得到一个以Fzk为变量的一元二次方程f(Fzk)。对该方程求极小值，得到满足各节点法向力方差最小的Fzk值，从而得到各节点法向力。应用最终的单元节点内力Fint对客车侧翻一步碰撞算法进行改进，可有效修正弯曲变形的影响，同时提高算法的计算效率。

3 应用实例

本文以某12 m 公路客车车身结构的典型车身段作为分析对象，进行侧翻碰撞试验；将经过弯曲修正的四节点膜单元应用于侧翻一步碰撞算法进行模拟，并与初始算法及LS-DYNA 结果和侧翻试验结果对比，检验所提单元模型的实际工程效果。

图2 单元局部坐标系及节点内力示意图Fig.2 Schematic diagram of element local coordinate system and nodal internal forces

图3 为所选待侧翻试验典型车身段模型。根据标准GB 17578—2013［13］及ECE R66法规［14］，按照车身段实际重量配重，保证重心位置与实车相同，车身段由一个简单的翻滚支架固定，最小离地间隙与实车保持一致。

图3 典型车身段模型Fig.3 Typical bus body section model

将该车身段模型置于图3所示侧翻试验台固定轴边缘位置，车身段处于自由静止状态；侧翻试验台绕着固定轴以0.1(。)·S-1的准静态速度缓慢旋转抬升，车身段试验模型随着试验台的升高旋转，当重心位置越过垂直中心线达到侧翻临界状态后开始自由侧翻下落，车身段骨架撞向地面，试验结束。

图4 所示为图3 所示待试验典型车身段有限元模型，共离散四边形单元259 976 个，节点258 368个。材料性能参数采用Holloman幂次硬化法则强化系数k=710 MPa，应变硬化指数n为0.2，初始应变ε0为0，厚向异性系数r为1.3。

图5 所示为经过改进的客车侧翻一步碰撞算法（简称“改进算法”）模拟的结构最终变形云图，图6为初始算法结果，图7 所示为LS-DYNA 仿真结果，图8为侧翻试验结果。

据图5～图8 对比发现，4 种方式结构最终变形趋势吻合，改进算法计算结果与初始算法、LSDYNA及侧翻碰撞结果趋势基本一致。

为进一步对结构变形进行定量分析，将侧翻试验提供的若干测点作为样本点，在车身段封闭环①和②两侧立柱上各选取11个测点进行数据采集，得到两侧立柱的变形量情况，并作为与改进算法、初始算法和LS-DYNA 仿真分析结构变形的对比数据，测点位置如图9所示。参考侧翻试验中测点选取方式，有限元模型测点的选取位置如图10 所示，两点之间的距离为100 mm。

图4 典型车身段有限元模型Fig.4 Finite element analysis model of typical bus body section

图5 改进算法模拟的结构最终变形云图（单位：m）Fig.5 Final deformation contour plot of improved algorithm(unit:m)

图6 初始算法模拟的结构最终变形云图（单位：m）Fig.6 Final deformation contour plot of original algorithm(unit:m)

图7 LS-DYNA仿真的结构最终变形云图Fig.7 Final deformation contour plot of LS-DYNA

图8 侧翻试验结构最终变形Fig.8 Final deformation of structure of rollover test

图9 侧翻试验测点选取位置Fig.9 Gauging point location for selection of roll over test

改进算法、初始算法、LS-DYNA 仿真及侧翻试验的封闭环①和②两侧立柱的变形量数据如表1所示。由表1可知，部分数据存在微小偏差，是由于受实际侧翻试验制备及测量过程的偶然性影响导致，但误差均在工程允许范围内。

图10 有限元模型测点选取位置Fig.10 Gauging point location for selection of finite element model

表1 各封闭环两侧立柱变形量统计Tab.1 Data of deformations of both sides of each pillar

图11 和图12 所示为典型车身段的封闭环①和②两侧立柱变形量对比柱状图。

图11 封闭环①两侧立柱变形量对比Fig.11 Comparison of deformations of closed-loop ①

图12 封闭环②两侧立柱变形量对比Fig.12 Comparison of deformations of closed-loop ②

从图12可看出，由于受到试验制备及实际测量过程的误差影响，侧翻试验数据在稳定性方面波动相对比较明显，但对改进算法的有效性检验，具有一定工程参考价值。通过对比分析图11 和图12，4 种结果获得的各测点数据走势基本一致，验证了结构最终变形趋势吻合的结论；改进算法与初始算法间的误差平均为2.2%，两者计算精度基本相当；改进算法与LS-DYNA 仿真分析之间误差平均约为4.1%，计算精度基本满足实际需要；改进算法与侧翻试验结果的误差平均约为11.3%，小于有限元工程计算误差的经验值范围15%，计算精度在实际可接受范围内。

接着对比计算效率。算法模拟采用的硬件平台为 DELL 工作站，配置 8 核 Intel i7 处理器、32G 内存及8T固态硬盘。表2为改进算法与初始算法和LSDYNA仿真的模拟效率对比结果，由表2可知，改进算法模拟时间约为初始算法的3/5，改进算法的模拟时间约为LS-DYNA仿真的1/15。

表2 模拟效率对比Tab.2 Comparison of simulation efficiency

综合对比计算精度和计算效率，改进算法虽精度与初始算法基本相当，但计算效率得到了大幅提升，检验了本文所提新型单元在实际工程应用中的有效性，可在客车结构设计初期，在基本保证计算精度的前提下，更快速地获得侧翻碰撞结构的最终变形，对结构侧翻碰撞安全性进行快速评价。

4 结论

本文针对传统四节点膜单元无法进行弯曲的缺陷，基于单元内力和弯矩的平衡原理，提出一种可进行弯曲修正的新型四节点膜单元模型。将其应用于客车侧翻一步碰撞快速模拟算法，简便有效地考虑了客车侧翻过程中由弯曲变形引起的节点内力。通过对典型车身段进行侧翻碰撞模拟，并与初始算法、LS-DYNA及侧翻试验结果对比，4种方法结构最终变形趋势比较吻合，改进算法在略微牺牲一点计算精度的情况下，计算效率得到了显著提高。

本文所提出可进行弯曲修正的四节点膜单元模型，对一般膜单元无法计算单元法向力的缺点进行了修正。整个改进过程依旧采用膜单元理论，节点的整体自由度并未增加，整体刚度矩阵带宽也未改变。既保持了原始膜单元计算效率高的优点，又获得了可观的计算精度，具有一定的实际工程应用价值。

本文改进的客车侧翻一步碰撞算法，在基本保证计算精度的同时，计算效率进一步得到了大幅提高，可在结构设计初期对客车侧翻安全性能进行快速评价，缩短产品开发周期，降低研发成本；同时，为后续开展基于客车侧翻碰撞安全性能的结构尺寸优化、参数优化、拓扑优化及灵敏度分析等算法提供必要支撑条件；并可对非线性隐式理论在结构碰撞模拟中的应用进行探索，为更复杂的客车及乘用车正碰撞快速算法研究奠定坚实的前期基础。