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

摘要

针对传统四节点膜单元无法进行弯曲的缺陷，基于单元内力和弯矩的平衡原理，提出一种可进行弯曲修正的新型四节点膜单元模型，将其应用于客车侧翻一步碰撞算法，对初始算法进行改进。通过对典型车身段模型进行模拟，并与客车侧翻一步碰撞初始算法、LS-DYNA仿真及侧翻试验结果在计算精度和计算效率方面分别进行对比，检验了所提新型单元模型的有效性和实际工程应用价值。

关键词

四节点膜单元; 平衡原理; 弯曲修正; 客车侧翻; 一步碰撞算法

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

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

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

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

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

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

将碰撞开始状态的结构作为原始构形${\mathit{X}}_{\mathrm{0}}$。此时车体未发生变形，结构动能$E_{\mathrm{d}}$由式（1）计算：

$$E_{\mathrm{d}}=\frac{\mathrm{1}}{\mathrm{2}}J{\omega }^{\mathrm{2}}=Mg\mathrm{\Delta }h$$

（1）

式中：$M$为车体质量；$\mathrm{\Delta }h$为车体重心下降高度；$J$为车体绕假定转轴的转动惯量；$\omega$为车体角速度。

碰撞开始状态结构各节点的速度$v_{\mathrm{0}}$由式（2）计算：

$$v_{\mathrm{0}i}=r_i\omega ,i=\mathrm{1},...,n$$

（2）

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

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

$${\mathit{U}}_{\mathrm{0}}={\mathit{x}}_{\mathrm{0}}-{\mathit{X}}_{\mathrm{0}}$$

（3）

式中：${\mathit{x}}_{\mathrm{0}}$为结构变形最大构型。

此时车体动能$E_{\mathrm{d}}$在碰撞中主要转换为结构形变能，结构形变能$W$由式（4）计算^［12］：

$$W={\displaystyle \sum _{\mathrm{i}=\mathrm{1}}^{\mathrm{N}}}\left({{\displaystyle \underset{{\mathrm{V}}^{\mathrm{e}}}{\overset{}{\iiint }}}\mathit{\epsilon }}^{\mathrm{T}}\mathit{\sigma }\mathrm{d}{\mathrm{V}}^{\mathrm{e}}\right)$$

（4）

式中：$\mathit{\epsilon }$为单元应变；$\mathit{\sigma }$为单元Cauchy应力；${\mathrm{V}}^{\mathrm{e}}$为单元体积；$N$为单元数。

判断此刻的结构形变能$W$与车体动能$E_{\mathrm{d}}$是否满足式（5）的能量关系假定：

$$W=E_{\mathrm{d}}$$

（5）

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

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

$$\mathit{R}\left(U_{{}^i}\right)={\mathit{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\left(U_{{}^i}\right)-{\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}\left(U_{{}^i}\right)=\mathrm{0}-{\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}\left(U_{{}^i}\right)\ne \mathrm{0}$$

（6）式中：${\mathit{F}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\left(U_{{}^i}\right)$为结构所受外力；${\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}\left(U_{{}^i}\right)$为结构所受内力。

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

$$\begin{array}{l}{\mathit{K}}_{\mathrm{T}}\left(U_{{}^i}\right)\mathrm{\Delta }{\mathit{U}}_{{}^i}={\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}\left(U_{{}^i}\right)(\mathrm{7}\mathrm{a})\\ {\mathit{U}}_{{}^{i+\mathrm{1}}}={\mathit{U}}_{{}^i}+\eta \mathrm{\Delta }{\mathit{U}}_{{}^i}(\mathrm{7}\mathrm{b})\end{array}$$

式中：切线刚度阵${\mathit{K}}_{\mathrm{T}}\left(U_{{}^i}\right)={\left[-\frac{\partial {\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}\left(U\right)}{\partial U}\right]}_{U=U_{{}^i}}$；$\eta$为松弛因子，取值在0~1之间。

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

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

本文提出的可进行弯曲修正的新型四节点膜单元模型如图1所示。图中$i$，$j$，$k$，$l$为单元节点在整体坐标系下的坐标；$i_{\mathrm{u}}$，$j_{\mathrm{u}}$，$k_{\mathrm{u}}$，$l_{\mathrm{u}}$为单元上表面4个节点在整体坐标系下的坐标；$i_{\mathrm{d}}$，$j_{\mathrm{d}}$，$k_{\mathrm{d}}$，$l_{\mathrm{d}}$为单元下表面4个节点在整体坐标系下的坐标；$n_{\mathrm{i}}$，$n_{\mathrm{j}}$，$n_{\mathrm{k}}$，$n_{\mathrm{l}}$为单元各节点法线。与传统壳单元模型不同，这种弯曲修正的四节点膜单元节点内力由两部分构成：一部分是由面内变形引起的节点切向力，另一部分是由弯曲变形引起的节点弯矩和法向力，如式（8）所示：

$${\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}={\mathit{F}}_{\mathrm{m}}+{\mathit{F}}_{\mathrm{b}}$$

（8）

式中：${\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}$为节点的总内力；${\mathit{F}}_{\mathrm{m}}$为因面内变形而引起的节点内力；${\mathit{F}}_{\mathrm{b}}$为因弯曲变形而引起的节点弯曲内力。${\mathit{F}}_{\mathrm{m}}$和${\mathit{F}}_{\mathrm{b}}$的向量形式由式（9）表示：

$$\begin{array}{l}{\mathit{F}}_{\mathrm{m}}={\left\{{\mathit{F}}_{xi}{\mathit{F}}_{yi}\mathrm{0}{\mathit{F}}_{xj}{\mathit{F}}_{yj}\mathrm{0}{\mathit{F}}_{xk}{\mathit{F}}_{yk}\mathrm{0}{\mathit{F}}_{xl}{\mathit{F}}_{yl}\mathrm{0}\right\}}^{\mathrm{T}}\\ {\mathit{F}}_{\mathrm{b}}={\left\{\mathrm{0}\mathrm{0}{\mathit{F}}_{zi}\mathrm{0}\mathrm{0}{\mathit{F}}_{zj}\mathrm{0}\mathrm{0}{\mathit{F}}_{zk}\mathrm{0}\mathrm{0}{\mathit{F}}_{zl}\right\}}^{\mathrm{T}}\end{array}$$

（9）

面内变形引起的节点内力${\mathit{F}}_{\mathrm{m}}$的计算与普通膜单元相同。为计算单元节点弯曲内力${\mathit{F}}_{\mathrm{b}}$，首先利用四节点膜单元沿节点法线等距的原理，计算弯曲后单元上下表面应力。图1所示为最大变形构形${\mathit{x}}_{\mathrm{0}}$下的某一单元$A(i,j,k,l)$，单元厚度为$t$。该单元节点在整体坐标系下的坐标分别为$i(x_i,y_i,z_i)$，$j(x_j,y_j,z_j)$，$k(x_k,y_k,z_k)$，$l(x_l,y_l,z_l)$，各节点法线分别为$n_i(l_i,m{}_i{}^{},q_i)$，$n_j(l_j,m{}_j{}^{},q_j)$，$n_k(l_k,m{}_k{}^{},q_k)$，$n_l(l_l,m{}_l{}^{},q_l)$，节点法线由节点相邻单元加权平均求得。沿着单元$A$的4个节点法线方向及负方向平行等距移动半个板厚$\frac{t}{\mathrm{2}}$，可得其上、下两个表面$A_{\mathrm{u}}$和$A_{\mathrm{d}}$，同时分别计算其两个表面的应力。规定上表面$A_{\mathrm{u}}$四个节点的坐标分别为$i(x_{\mathrm{u}i},y_{\mathrm{u}i},z_{\mathrm{u}i})$，$j(x_{\mathrm{u}j},y_{\mathrm{u}j},z_{\mathrm{u}j})$，$k(x_{\mathrm{u}k},y_{\mathrm{u}k},z_{\mathrm{u}k})$，$l(x_{\mathrm{u}l},y_{\mathrm{u}l},z_{\mathrm{u}l})$，下表面$A_d$四个节点的坐标分别为$i(x_{\mathrm{d}i},y_{\mathrm{d}i},z_{\mathrm{d}i})$，$j(x_{\mathrm{d}j},y_{\mathrm{d}j},z_{\mathrm{d}j})$，$k(x_{\mathrm{d}k},y_{\mathrm{d}k},z_{\mathrm{d}k})$，$l(x_{\mathrm{d}l},y_{\mathrm{d}l},z_{\mathrm{d}l})$。单元上、下表面节点处的坐标可由式（10）计算：

$$A_{\mathrm{u}}:\left\{\begin{array}{c}{x}_{\mathrm{u}r}=\frac{t}{\mathrm{2}}{l}_r+x_r\\ y_{\mathrm{u}r}=\frac{t}{\mathrm{2}}{m}_r+y_r\\ z_{\mathrm{u}r}=\frac{t}{\mathrm{2}}{q}_r+z_r\end{array}\right.r=i,j,k,l$$

（10a）

$$A_{\mathrm{d}}:\left\{\begin{array}{c}{x}_{\mathrm{d}r}=-\frac{t}{\mathrm{2}}{l}_r+x_r\\ y_{\mathrm{d}r}=-\frac{t}{\mathrm{2}}{m}_r+y_r\\ z_{\mathrm{d}r}=-\frac{t}{\mathrm{2}}{q}_r+z_r\end{array}\right.r=i,j,k,l$$

（10b）

这样，可得到单元上、下表面节点处的位移。以上表面单元$A_{\mathrm{u}}$为例，整体坐标系下各节点的位移${\mathit{Q}}_{\mathrm{u}}^{\mathrm{e}}$可由式（11）计算：

$${\mathit{Q}}_{\mathrm{u}}^{\mathrm{e}}=\left\{\begin{array}{c}{u}_{\mathrm{u}i}\\ u_{\mathrm{u}j}\\ u_{\mathrm{u}k}\\ u_{\mathrm{u}l}\end{array}\right\}=\left\{\begin{array}{c}{x}_{\mathrm{u}i}-X_{\mathrm{u}i}\\ x_{\mathrm{u}j}-X_{\mathrm{u}j}\\ x_{\mathrm{u}k}-X_{\mathrm{u}k}\\ x_{\mathrm{u}l}-X_{\mathrm{u}l}\end{array}\right\}$$

（11）

式中：${\mathit{u}}_{\mathrm{u}r}$表示节点$r_{\mathrm{u}}(r=i,j,k,l)$在整体坐标系中的位移向量；${\mathit{x}}_{\mathrm{u}r}$表示节点$r_{\mathrm{u}}(r=i,j,k,l)$在最大变形构形$x_{\mathrm{0}}$下的整体坐标；${\mathit{X}}_{\mathrm{u}r}$表示节点$r_{\mathrm{u}}(r=i,j,k,l)$在原始构形$X_{\mathrm{0}}$下对应节点的坐标。

图1 可进行弯曲修正的四节点膜单元

Fig. 1 Four-node membrane element considering bending modification

根据大变形几何关系可以求得Green变形张量下的单元主伸长${\lambda }_{\mathrm{u}r}(r=\mathrm{1,2},\mathrm{3})$，在主应变空间的对数应变${\epsilon }_{\mathrm{u}r}(r=\mathrm{1,2},\mathrm{3})$可由式（12）计算：

$${\mathit{\epsilon }}_{\mathrm{u}r}={\left\{\begin{array}{ccc}{\epsilon }_{\mathrm{u}\mathrm{1}}& {\epsilon }_{\mathrm{u}\mathrm{2}}& {\epsilon }_{\mathrm{u}\mathrm{3}}\end{array}\right\}}^{\mathrm{T}}={\left\{\begin{array}{ccc}\mathrm{l}\mathrm{n}{\lambda }_{\mathrm{u}\mathrm{1}}& \mathrm{l}\mathrm{n}{\lambda }_{\mathrm{u}\mathrm{2}}& \mathrm{l}\mathrm{n}{\lambda }_{\mathrm{u}\mathrm{3}}\end{array}\right\}}^{\mathrm{T}}$$

（12）

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

$${\mathit{\sigma }}_{\mathrm{u}}=\frac{\overline{\sigma }}{\overline{\epsilon }}{\mathit{P}}^{-\mathrm{1}}{\mathit{\epsilon }}_{\mathrm{u}}$$

（13）

式中：$\overline{\sigma }$为等效应力；且$\overline{\sigma }=k{\left(\overline{\epsilon }+{\epsilon }_{\mathrm{0}}\right)}^n$；$k$为强化系数；$n$为应变硬化指数；${\epsilon }_{\mathrm{0}}$为初始应变；$\overline{\epsilon }$为等效应变；${\mathit{P}}^{-\mathrm{1}}$为各向异性系数矩阵。

这样，即可获得上表面单元$A_{\mathrm{u}}$任一点的应力值${\mathit{\sigma }}_{\mathrm{u}}$，同理，可获得下表面单元$A_{\mathrm{d}}$任一点的应力值${\sigma }_{\mathrm{d}}$。

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

$$M=\left\{\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right\}$$

（14）

式中：$M_x$为垂直于x轴截面单位长度弯矩；$M_y$为垂直于y轴截面单位长度弯矩；$M_{xy}$为垂直于x（y）轴截面单位长度扭矩，$M_{xy}=M_{yx}$。

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

$${\sigma }_x=\frac{\mathrm{12}{M}_x}{t^{\mathrm{3}}}z$$

（15a）

$${\sigma }_y=\frac{\mathrm{12}{M}_y}{t^{\mathrm{3}}}z$$

（15b）

$${\tau }_{xy}={\tau }_{yx}=\frac{\mathrm{12}{M}_{xy}}{t^{\mathrm{3}}}z$$

（15c）

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

$$M_x=\frac{{\sigma }_x{t}^{\mathrm{2}}}{\mathrm{6}}$$

（16a）

$$M_y=\frac{{\sigma }_y{t}^{\mathrm{2}}}{\mathrm{6}}$$

（16b）

$$M_{xy}=M_{yx}=\frac{{\tau }_{xy}{t}^{\mathrm{2}}}{\mathrm{6}}$$

（16c）

需要注意的是式（15）中，${\sigma }_x,{\sigma }_y,{\tau }_{xy}$为由板弯曲变形引起的三个应力分量，式（13）给出的是由面内变形和弯曲变形共同引起的应力。故需进一步将应力分量分解为面内变形应力${\sigma }_{{}^{\mathrm{p}}}$和弯曲变形应力${\sigma }_{{}^{\mathrm{b}}}$两部分。根据薄板弯曲理论，弯曲应力应该满足${\sigma }_{\mathrm{u}\mathrm{b}}=-{\sigma }_{\mathrm{d}\mathrm{b}}$。令$\Delta =\frac{{\sigma }_{\mathrm{u}}+{\sigma }_{\mathrm{d}}}{\mathrm{2}}$，则${{\sigma }_{\mathrm{u}\mathrm{b}}}^{\mathrm{*}}={\sigma }_{\mathrm{u}}-\Delta$，${{\sigma }_{\mathrm{d}\mathrm{b}}}^{\mathrm{*}}={\sigma }_{\mathrm{d}}-\Delta$。

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

$$\begin{array}{l}{M}_x^{\mathrm{*}}={\int }_{\mathrm{s}}(M_x+M_{xy})\mathrm{d}\mathrm{s}\\ M_y^{\mathrm{*}}={\int }_{\mathrm{s}}(M_y+M_{yx})\mathrm{d}\mathrm{s}\end{array}$$

（17）

式中：s为单元的边界。

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

$$\left\{\begin{array}{c}{F}_{zi}\cdot y_i+F_{zj}\cdot y_j+M_x^{\mathrm{*}}=\mathrm{0}\\ -F_{zi}\cdot x_i-F_{zj}\cdot x_j-F_{zk}\cdot x_k+M_y^{\mathrm{*}}=\mathrm{0}\\ F_{zi}+F_{zj}+F_{zk}+F_{zl}=\mathrm{0}\end{array}\right.$$

（18）

图2 单元局部坐标系及节点内力示意图

Fig. 2 Schematic diagram of element local coordinate system and nodal internal forces

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

$$\left\{\begin{array}{c}{F}_{zi}=F_{zk}\cdot \frac{-x_k{y}_i{y}_j}{x_i{y}_i{y}_j-x_j{y}_i^{\mathrm{2}}}+\frac{x_j{y}_i{M}_x^{\mathrm{*}}+y_i{y}_j{M}_y^{\mathrm{*}}}{x_i{y}_i{y}_j-x_j{y}_i^{\mathrm{2}}}\\ F_{zj}=F_{zk}\cdot \frac{x_k{y}_i}{x_i{y}_j-x_j{y}_i}-\frac{x_i{M}_x^{\mathrm{*}}+y_i{M}_y^{\mathrm{*}}}{x_i{y}_j-x_j{y}_i}\\ F_{zl}=F_{zk}\cdot \frac{x_k{y}_i{y}_j-x_i{y}_i{y}_j-x_k{y}_i^{\mathrm{2}}+x_j{y}_i^{\mathrm{2}}}{x_i{y}_i{y}_j-x_j{y}_i^{\mathrm{2}}}+\frac{(x_i{y}_i-x_j{y}_i)M_x^{\mathrm{*}}+(y_i^{\mathrm{2}}-y_i{y}_j)M_y^{\mathrm{*}}}{x_i{y}_i{y}_j-x_j{y}_i^{\mathrm{2}}}\end{array}\right.$$

（19）

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

$$\begin{array}{l}f(F_{zk})=\frac{\mathrm{1}}{\mathrm{4}}\left[(F_{zi}{-\overline{F})}^{\mathrm{2}}+(F_{zj}{-\overline{F})}^{\mathrm{2}}+(F_{zk}{-\overline{F})}^{\mathrm{2}}+(F_{zl}{-\overline{F})}^{\mathrm{2}}\right]\\ \overline{F}=\frac{F_{zi}+F_{zj}+F_{zk}+F_{zl}}{\mathrm{4}}=\mathrm{0}\end{array}$$

（20）

式中：$\overline{F}$为4个节点法向力的平均值。

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

3 应用实例

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

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

图3 典型车身段模型

Fig. 3 Typical bus body section model

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

图4所示为图3所示待试验典型车身段有限元模型，共离散四边形单元259 976个，节点258 368个。材料性能参数采用Holloman幂次硬化法则$\overline{\sigma }=k{\left(\overline{\epsilon }+{\epsilon }_{\mathrm{0}}\right)}^n$，强化系数$k=\mathrm{710}\mathrm{M}\mathrm{P}\mathrm{a}$，应变硬化指数$n$为0.2，初始应变${\epsilon }_{\mathrm{0}}$为0，厚向异性系数$r$为1.3。

图4 典型车身段有限元模型

Fig. 4 Finite element analysis model of typical bus body section

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

图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

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

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

图9 侧翻试验测点选取位置

Fig. 9 Gauging point location for selection of roll over test

图10 有限元模型测点选取位置

Fig. 10 Gauging point location for selection of finite element model

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

表1 各封闭环两侧立柱变形量统计

Tab. 1 Data of deformations of both sides of each pillar

测点 序号	改进算法/mm		初始算法/mm		LS-DYNA /mm		侧翻试验/mm
	封闭环①	封闭环②	封闭环①	封闭环②	封闭环①	封闭环②	封闭环①	封闭环②
1	59	61	61	65	65	70	80	80
2	72	78	76	81	78	83	90	90
3	83	92	88	94	90	94	105	100
4	97	101	99	104	105	108	115	110
5	119	118	122	124	125	127	125	120
6	134	137	137	141	142	143	130	130
7	124	127	126	127	126	129	120	145
8	101	106	106	111	108	112	110	130
9	86	95	91	97	93	99	100	120
10	77	84	80	84	82	88	85	105
11	66	71	70	76	74	80	70	85

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

图11 封闭环①两侧立柱变形量对比

Fig. 11 Comparison of deformations of closed-loop ①

图12 封闭环②两侧立柱变形量对比

Fig. 12 Comparison of deformations of closed-loop ②

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

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

表2 模拟效率对比

Tab. 2 Comparison of simulation efficiency

时间	改进算法	初始算法	LS-DYNA
T/min	12	20	180

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

4 结论

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

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

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

作者贡献声明

王童：负责整个论文的思路整理、方法研究、理论推导、方案验证、试验及数据整理、论文撰写及修改等全部工作，是此文章的主要贡献者。