汽车悬架系统位于车轮和车身之间，除满足汽车的舒适性、操纵稳定性等相关性能要求外，还须满足结构强度和耐久性要求。在汽车实际行驶过程中，地面对轮胎的载荷可以分解为垂直力、侧向力、纵向力和制动力矩，通过车轴作用在悬架上，表现为多轴随机载荷，复杂且严重，而影响悬架耐久性的主要因素正是这4个载荷。

为了确保悬架系统能够满足强度、工作寿命和缩短研发周期的要求，室内道路模拟试验在汽车工业领域被普遍采用，通过电脑控制的多轴耦合疲劳加载试验机复现车辆在试车场测试过程中的应力应变时间历程^[1-3]。此外，主流的伺服液压测试系统供应商还研发了与试验机相匹配的控制软件。1977年，美国MTS公司发布第一版远程参数控制(remote parameter control，RPC)系统^[4]并一直不断更新至今；1979年，德国SCHENCK公司发布了ITFC(iterative transfer function compensation)系统^[5]；20世纪90年代末期，英国INSTRON公司和比利时LMS公司也相继研发推出了类似具有迭代补偿功能的道路模拟系统^[6-8]。

近年来，为了进一步缩短研发周期、降低研发成本以提高市场竞争力，许多车辆生产厂家更青睐于使用计算机虚拟仿真的方法来进行车辆耐久性设计。目前，主流的两种虚拟仿真方法是虚拟试车场(virtual proving ground，VPG)法^[9-10]和混合仿真法^[11-12]，这样就可以在车辆的初样机生产出来之前对车辆结构进行改进和再设计^[13]。

VPG方法的核心是建立车辆整车及轮胎和各种典型路面的有限元模型，进而获得车辆各零部件的应力应变时间历程，从而对车辆的使用寿命和可靠性进行预估，但这种方法一直面临着两大主要挑战：第一，车辆轮胎是一个非线性程度很高的部件，目前在国内外的工程领域内仍无法获得足够准确且适用的轮胎模型；第二，轮胎与地面之间接触关系的仿真计算同样也是一个非线性程度很高的计算过程，非线性计算不仅对有限元模型的网格质量要求很高，而且仿真结果也常常难以收敛，影响仿真的效率和准确性。因此，此方法能够仿真的时长也只有10~20 s，这对于车辆耐久性研究是远远不够的。

混合仿真方法可以理解为VPG仿真方法的一种近似方法，它是把轮胎和车辆悬架等非线性程度很高且难以模拟的部件用实物代替，这些部件被安装在受仿真计算机控制的试验台上，通过这些部件上的传感器传回的响应数据并与计算机仿真环境构成闭环，进而完成一个闭环的混合仿真，这种方法解决了非线性零部件建模困难的问题，但是对计算机的配置需求更高，因为这种仿真方法不是离线计算，而是在线实时计算，这样在成本方面，与实车测试相比优势不大。

在保证模型准确性的前提下，本文用加载杆取代轮胎有限元模型并对板簧结构进行了合理简化，从而将整个驱动后桥模型简化为线性系统，利用频响函数法反求出能够复现测点处应力应变时间历程的载荷，建立了一种确保应力应变时间历程能够精确复现的虚拟多轴向加载疲劳试验的仿真方法。该方法通过线性计算，避免了对计算机配置的苛刻要求、非线性零部件的建模、非线性计算不收敛及仿真成本较高的问题，同时也解决了仿真时长的受限制的问题。

1 有限元模型的建立 1.1 材料参数的确定

合理并正确地建立后驱动桥的有限元模型是确保得到准确仿真结果的基础，有限元模型中主要使用的是实体单元(CTETRA)和壳单元(CTRIA3)。由于驱动桥壳是本文分析的重点，所以该部分的单元尺寸为2 mm，其他部分的尺寸则为5 mm；板簧前后吊耳处的运动副通过RBE2单元建立，满足其运动的协调性并与实际情况一致；两侧加载杆长度为轮胎滚动半径306 mm。该有限元模型中共有2 580 409个单元和627 142个节点，如图 1所示。

模型中各部分的单元尺寸，壳单元厚度，以及对应的材料属性，如表 1所示。

由于设计时车辆是满载状态，因此在做钢板弹簧线性简化时，三片钢板弹簧的两端和中间均使用壳单元模拟钢片直接连接以模拟满载时的受力状态。钢板弹簧满载时的合成刚度为307 N/mm，故钢板弹簧的弹性模量E=100 GPa，Poisson比μ=0.3，密度ρ=7 900 kg/m³；传动轴支撑处的材料为铸铁，其E=170 GPa，μ=0.25，ρ=7 900 kg/m³；橡胶衬套，E=1 MPa，μ=0.48，ρ=1 650 kg/m³；其他部分，E=203 GPa，μ=0.3，ρ=7 900 kg/m³。

1.2 静态标定

进行静态标定的目的是校核和证明在静态载荷作用下本文中所建立的后驱动桥有限元模型的准确性。只有这样才能保证该有限元模型在动态载荷的激励下所得出的响应依然正确，因此，分别对后驱动桥进行了有限元模型静态标定和实验室静态标定试验。测点布置如图 2所示，测点L₁、L₂、L₄、R₁、R₂、R₄和R₅测量的是沿Y轴的正应变，L₃和R₃(各2片)测量的是与轴线成45°的剪切应变。因此，L₁和R₁对纵向力敏感，L₂、L₄、R₂和R₄对垂向力和侧向力敏感，L₃和R₃对制动力矩敏感。

使用Altair HyperWorks 12.0的RADIOSS求解器进行有限元计算，由于有限元的输出结果中没有45°剪切应变输出结果，因此根据广义Hooke定律，沿Y轴的正应变ε_y计算如下：

$ {\varepsilon _y} = \frac{1}{E}\left[{{\sigma _y}-\mu \left( {{\sigma _x} + {\sigma _z}} \right)} \right]. $

(1)

其中：σ_x、σ_y、σ_z分别为沿X、Y、Z轴的正应力。沿与Y轴呈45°的剪切应变ε_45°的计算如下：

$ {\varepsilon _{45{\rm{^\circ }}}} =-\frac{{{\tau _{xy}}}}{E}\left( {1 + \mu } \right). $

(2)

其中：τ_xy是沿与Y轴呈45°的剪切应力。

静态标定试验分为4种加载工况：垂向力(5 kN)、垂向力(5 kN)加侧向力(5 kN)、垂向力(5 kN)加纵向力(5 kN)、垂向力(5 kN)加制动力(3 kN)，其实验结果与有限元静态标定结果对比如表 2和3所示。

由以上静态标定的结果对比可知，有限元静态标定结果与试验结果的误差在1.3%~24.2%，说明后驱动桥的静态有限元模型是准确的。

2 模型的频响识别

在线性系统理论中，一个系统的频率响应函数是输出信号的快速Fourier变换除以对应输入信号的快速Fourier变换，在工程运用中，通常使用白噪声信号作为系统识别的激励信号，但从严格意义上讲，白噪声信号是不能满足采样定理的，因此将白噪声信号进行快速Fourier变换之后就会使频响函数的识别的精度不够，本文则使用“类脉冲”信号作为系统识别的激励信号，定义如下：

$ S\left( t \right) = \frac{A}{N}\sum\limits_{n = 1}^N {\cos } \left( {2\pi {f_n}t} \right). $

(3)

其中：S(t)是类脉冲信号，A是每个余弦信号的幅值，N是初相位为零的余弦信号的个数，f_n是每个余弦信号的频率。这样可以人为控制信号的最高截止频率以保证识别信号满足采样定理，这种信号的快速Fourier变换结果的实部是一条值为A的水平线，虚部是一条值为0的水平线，这是一个非常优良频响识别特性。

后驱动桥动态有限元模型的输入载荷共分为左、右两侧的垂向力、侧向力、纵向力和制动力矩共计8个载荷，依次对模型输入“类脉冲”信号，并求在每个载荷单独作用下8个载荷敏感点的响应信号，根据线性系统理论，该系统的频响函数H(f)就是响应信号的Fourier变换与输入信号的Fourier变换之比，表示如下：

$ \mathit{\boldsymbol{H}}\left( f \right) = \frac{{\mathit{\boldsymbol{Y}}\left( f \right)}}{{\mathit{\boldsymbol{X}}\left( f \right)}}. $

(4)

其中：X(f)为输入(识别激励)信号的快速Fourier变换，Y(f)为输出(响应)信号的快速Fourier变换。

3 目标信号的获取与复现

为了获得目标响应信号，车辆在海南试车场进行了测试并获取了实验数据，因此，根据获得的目标信号和频响函数就可以逆向求得能够复现目标信号的输入信号的快速Fourier变换形式：

$ \mathit{\boldsymbol{X}}\left( f \right) = {\left[{\mathit{\boldsymbol{H}}\left( f \right)} \right]^{ -1}}\mathit{\boldsymbol{Y}}\left( f \right). $

(5)

其中：[H(f)]^-1是H(f)的逆，最后将X(f)通过快速Fourier变换的逆变换变回时域，就能够获得复现目标信号的载荷。

由于本文所研究的后驱动桥的有限元模型是一个标准的线性系统，采用频响识别法可以求出准确的系统传递函数，把求出的载荷输入到模型中可以直接复现目标信号而不需要迭代，复现的综合误差定义如下：

$ \begin{array}{l} {\rm{Error = 0}}{\rm{.8}}\frac{{\sum {{{\left[ {{s_{\rm{d}}}\left( n \right) - s\left( n \right)} \right]}^2}} }}{{\sum {s_{\rm{d}}^2\left( n \right)} }} + \\ \;\;\;\;\;\;0.2\frac{{{\rm{Max}}\left( {\left| {{s_{\rm{d}}}\left( n \right) - s\left( n \right)} \right|} \right)}}{{{\rm{Max}}\left( {\left| {{s_{\rm{d}}}\left( n \right)} \right|} \right)}}. \end{array} $

(6)

其中：s_d(n)是目标信号，s(n)是仿真得到的动态响应信号。

图 3是测点L₂的目标信号与仿真得到的动态响应信号的时域波形对比结果，由于精度非常高，为了更清晰地进行展示，将仿真数据向后错开200个数据点，图 3中仅包含了第40 s至第55 s区间内的部分数据。

复现信号的总仿真时长为540.672 s(270 336个数据点)，图 2中所示的测点L₁、R₁、L₂、R₂、L₃、R₃、L₄和R₄对应的综合复现误差分别为0.7%、0.3%、0.4%、0.3%、1.0%、4.7%、1.1%和0.9%。

测点R₅是一个对结构应力应变分布是否正确的验证测点，其时域波形对比如图 4所示，综合复现误差为18.8%。这充分证明了该方法所求出来的载荷能够保证整个驱动后桥的应力应变分布与实际情况保持一致，这也就为疲劳分析的准确性奠定了基础。

4 疲劳分析 4.1 疲劳分析理论

后驱动桥疲劳强危险点的疲劳寿命预测所采用的方法是局部应力-应变法^[14-16]，其应力-应变方程为：

$ \frac{{\Delta \varepsilon }}{2} = \frac{{\Delta \sigma }}{{2E}} + {\left( {\frac{{\Delta \sigma }}{{2K'}}} \right)^{\frac{1}{{n'}}}}. $

(7)

其中：Δε是应变变程，Δσ是应力变程，K′是循环强度系数，n′是循环应变硬化指数。

应变—寿命方程为：

$ \frac{{\Delta \varepsilon }}{2} = \frac{{{{\sigma '}_f}-{\sigma _0}}}{E}{\left( {2{N_f}} \right)^b} + {\varepsilon _f}{\left( {2{N_f}} \right)^c}. $

(8)

其中：N_f是裂纹形成时材料所承受的应变循环次数，文中对应的意义是在海南试车场循环数，对应的损伤D就是N_f的倒数，σ_f^′是疲劳强度系数，σ₀是平均应力，b是疲劳强度指数，ε_f^′是疲劳延展性系数，c是疲劳延展性指数。

后驱动桥壳所采用的材料是SAE1020钢，这种材料的疲劳性能参数如表 4所示。

4.2 疲劳寿命预测

将在节3求出的能够复现目标信号的载荷输入到后驱动桥的动态有限元模型中进行计算，可以得到模型中各个单元的动态应力应变响应并采用节4.1中的疲劳分析方法计算各个单元对应的疲劳损伤。从整体疲劳损伤的结果情况来看，后驱动桥的右侧的疲劳损伤大于左侧，这与在海南试车场道路实测的情况相关；最终结果显示A、B、C这3处的疲劳损伤最大，如图 5所示。

A、B和C这3处疲劳损伤比较大的区域对应的单元编号以及疲劳损伤如表 5所示。

由表 5的疲劳损伤计算结果可知，这3处均处在后驱动桥截面尺寸发生变化处，这些区域的应力集中比较严重，其中B处的疲劳损伤最大，因此在研发过程中，需要对这3个区域以及所有截面尺寸发生变化的区域的设计重视和改进，以保证后驱动桥有更高的疲劳强度和耐久性。

5 结论

本文仿真时长仅为540.672 s，理论上该仿真方法的时间长度是不受限制的。避免了建立整车和车辆轮胎等复杂模型，避免了接触等非线性计算，极大降低了对计算机性能的苛刻要求，显著地提高了计算速度和计算精度。进行频响识别时所使用的识别信号不是白噪声信号而是“类脉冲”信号，这种信号准确地识别出了系统的频响函数。目标信号的综合复现精度在0.3%~4.7%之间，这在工程应用中，可以认为精度较高，并且经过验证，整个结构的应力应变特征与实际的吻合度也较好。该方法计算出了后驱动桥的疲劳危险点的疲劳损伤，且预报出了驱动后桥的疲劳危险点的具体位置，为CAE设计提供了指导。