甘华权，李海艳

（广东工业大学机电工程学院，广州 510006）

0 引言

随着现代科学技术的不断发展和生活水平的不断提高，人们对汽车的功能需求越来越趋于多样化和多元化。除了作为日常代步行驶之外，人们对汽车舒适性和平稳性的要求也越来越高。由于汽车悬挂系统综合多种作用力，对车辆行驶的振动影响较大，进而汽车振动直接影响了汽车的稳定性、舒适性以及安全性，所以研究汽车悬挂系统的振动变形对改进汽车的稳定性、舒适性和安全性有着重大的意义[1]。

汽车悬挂系统属于柔性多体系统的范畴[2]，其在运动时会产生刚性与弹性变形的耦合运动[3]。所以汽车悬挂系统动力学建模需要考虑弹性变形，但是弹性变形有无限维自由度的特点，所以无法求得其精确解。目前比较常用的是通过有限单元法进行离散化处理进而求取其近似解[4]。例如，Patil等[5]利用有限元方法来求解柔性臂运动到特定位置情况下的挠度。许志华[6]使用有限元方法建立自卸车悬挂系统非线性有限元接触模型。Nabawy等[7]使用平面杆单元建立了一种综合的有限元模型来研究双叉臂悬架系统的动态响应。然而为了方便模型求解，很多学者使用有限单元法对汽车悬挂系统动力学建模中一般会选择简化悬挂系统模型同时忽略一部分小变形，虽然减少了一些计算量，但也降低了一些动力学模型的精度，而且有限单元法的精度和网格单元划分的大小相关，对于像汽车悬挂系统这样的复杂系统，即使在简化了模型的情况下，为了保证计算精度的网格单元会划分的比较小，计算量也仍然比较大。

因此，本文以汽车悬挂系统为研究对象，提出了一种利用模态方法进行降维求解的方法，该方法首先通过拉格朗日方程法建立准确的汽车悬挂系统的非线性柔性多体动力学模型。然后对悬挂系统的弹性变形进行离散化处理，并通过有限阶级数展开式精确表示该弹性变形。最后使用广义α方法迭代求解该模型，并与有限单元法的柔性多体动力学建模求解结果进行对比分析，得出该方法应用于柔性多体动力学建模的优势。

1 汽车悬挂系统动力学建模

1.1 建立参考坐标系

如图1所示，1、3、5、7为车轮，2、4、6、8为柔性臂，9为车架，10～13为减震器。分析整车悬挂系统需要建立一系列的坐标系来表示其位置，其中oxyz为整体坐标系，然后以质量块i(i=1,2,…,9)的重心oi为局部坐标系原点，建立与整体坐标系平行的初始状态的局部坐标系oixiyizi。

图1 汽车悬挂系统

整车悬挂系统根据解耦控制方法，可以由4个1/4悬挂系统和一个为车架与受到的外力组成的系统组成。在整车悬挂系统中，所有的运动部件全部考虑为均质质量块[8]。

如图2所示，以1/4悬挂系统为例进行分析。首先，1为车轮，2为柔性臂，10为减震器，可替换为弹簧阻尼器，质量忽略不计。整体坐标系为oxyz，车轮1的局部坐标系为o1x1y1z1，柔性臂2的局部坐标系为o2x2y2z2（图中未显示的x、x1与x2轴为由纸面向里）。由于柔性臂2的长度远大于其截面的直径，柔性臂的截面方向的变形变化率非常小，所以将柔性臂简化为空间梁结构。

图2 1/4悬挂系统示意图

假设在简化的空间梁2上存在某一点p，则该点p在局部坐标系o2x2y2z2上的位移矢量为：

式中：R2为局部坐标系o1x1y1z1的原点在整体坐标系中的位置；A2为将柔性臂局部坐标系位移转换到整体坐标系的旋转矩阵；u2r和u2e分别为p点在局部坐标系o2x2y2z2上的刚性和弹性变形位移。

1.2 动力学模型建立

考虑到汽车悬挂系统动力学方程非常复杂，而使用拉格朗日方程法是直接建立外力与运动部件变量的关系，不涉及运动部件之间的各种约束力，原理相对简单，所以本文利用拉格朗日方程法建立汽车悬挂系统的动力学模型。下面以1/4汽车悬挂系统为例，其拉格朗日方程为[9]：

式中：Qi为运动部件i广义坐标的广义外力；λ为拉格朗日乘子向量；Cqi为运动部件i的约束雅克比矩阵；qi为运动部件i的广义坐标；Ti为运动部件i的动能。

Ti计算如下：

式中：q̇i为qi对时间t的一阶微分；Mi为运动部件i的质量矩阵。

根据式（3）、虚功原理以及边界条件可求得1/4悬挂系统的动力学方程为：

式中：q̈ii为qi对时间t的二阶微分；Ki为运动部件i的刚度矩阵；Qvi为运动部件i二次速度矢量；Qci为运动部件i的约束雅克比矩阵余量。

其中，可根据以下式计算运动部件i二次速度矢量[10]：

进一步的，可以将式（4）写成以下的矩阵形式：

最后同理，可得整车悬挂系统动力学模型为：

2 动力学模型求解

2.1 动力学方程降维

由于拉格朗日方程法推导得到的汽车悬挂系统柔性多体动力学模型的方程是一组时变、非线性和强耦合的偏微分-积分方程，所以求汽车悬挂系统动力学模型的解析解非常困难。在这种情况下，一般是将无限维的动力学模型通过离散的方法进行降维，得到有限维的动力学模型。目前柔性多体动力学模型的常用的降维处理方法主要是有限单元法，该方法通过划分单元的大小保证其求解精度。单元划分越细时，离散后的动力学模型维度就越高，计算量就越大，计算效率就越低。而本文通过将弹性变形离散成基函数与对应时间系数乘积的级数展开式进而进行降维处理，将其代入前面求得的动力学模型即可得到降维后的低维近似模型。

首先根据柔性臂j的局部坐标系ojxjyjzj和柔性臂j的振型函数取特殊组合的三角函数为基函数：

式中：yj为柔性臂j在y轴上的位置表示；Lj为柔性臂j的长度，(i=1,2,…,N)，ωi为柔性臂的固有频率。

然后，将柔性臂j的弹性变形离散成式（8）的基函数与对应时间系数乘积的有限级数展开式表示：

式中：axi(t)、ayi(t)、aji(t)为待求的与时间相关的权函数。

根据柔性臂两端变形为零的边界条件，可得：

进而计算求得参数dji为：

最后，本文取前2阶级数展开式精确表示其柔性臂j的弹性变形：

将上述离散的弹性变形代入前面求得的汽车悬挂系统柔性多体动力学模型中，得到降维后的低维柔性多体动力学模型。

根据上面弹性变形的表示形式可以得到以下刚度矩阵[11]：

由上式可求得柔性臂j的形函数Sj、几何矩阵Bj和刚度矩阵Kj，式中qjf为柔性臂j变形对应的广义坐标，E为弹性模量。则有：

式中：Smj为质量块j的截面面积；Igxj为质量块j关于x方向的惯性矩；Igzj为质量块j关于z方向的惯性矩。

2.2 动力学方程求解

如图3所示，通过拉格朗日方程法可建立的汽车悬挂系统柔性多体动力学模型，通过求逆矩阵的方法求得初始加速度矢量和拉格朗日乘子。然后将初始时刻t0时刻的位移、速度、加速度等矢量作为广义α方法[12]的初始迭代值，通过广义α方法迭代求取t1时刻的相应位移、速度、加速度等，然后以t1时刻的相应位移、速度、加速度为作为t2时刻广义α方法迭代的初始值求t2时刻的相应位移、速度、加速度等矢量。然后依次类推，依次求解，直到误差小于预设误差值且执行到未尾时刻tn则程序迭代结束，并输出最后tn时刻相应的位移、速度、加速度矢量等。

图3 动力学模型求解流程

3 仿真与结果分析

通过具体参数进行数值仿真，车前端、后端的柔性臂长度分别为0.38 m、0.35 m，车轮、车架质量分别为20 kg、300 kg，弹性模量为200 GPa，密度为7800 kg/m3，弹簧弹性系数为10000 N/m，车架质心( 0.25,0,0.36)m，柔性臂截面内、外圆半径分别为0.015 m，0.01 m。

图4所示为沿x轴方向的弹性位移的比较，实线为本文方法求解得到的x轴弹性位移，虚线为有限元方法求解得到的x轴弹性模型，第一到第四行分别为柔性臂2、4、6、8的x轴弹性位移对比。

图4 不同时刻柔性臂x轴变形位移

图5所示为柔性臂中点沿z轴方向的弹性位移的比较，其中实线为本文方法表示的动力学模型求解得到的z轴弹性位移，虚线为有限元方法表示的动力学模型求解得到的z轴弹性模型，图5（a）～（d）分别为柔性臂2、4、6、8的z轴弹性位移对比。

图5 不同柔性臂中间节点变形位移

从图4、图5可知：首先，将弹性变形降维到64维的动力学建模求解的有限元方法与将弹性变形降维为二维的动力学模型求解的本文方法计算得到的不同时刻柔性臂x轴变形位移及不同柔性臂中间节点变形位移趋势一致，且结果相差无几，即精度相当；其次，本文方法求解所用时间仅为0.79 s，而有限元方法求解所用时间为2.3 s，其所用时间是本文方法的4倍多，即本文方法的计算效率远高于有限单元方法；最后，如果进一步提升求解动力学模型的精度，有限单元法的网格划分将会更多，计算量将大大增大，而本文则是只需增加一两阶基函数即可，所以在相同精度下，本文方法在求解效率方面优势非常显著。

4 结束语

本文通过以整车汽车悬挂系统作为研究对象，考虑弹性变形对汽车悬挂系统的影响，将柔性臂简化成柔性空间梁结构，不忽略柔性臂的轴向和截面变形，使建立的汽车悬挂系统的柔性多体动力学模型比以往的将柔性臂简化成平面梁表示动力学模型更完整。通过对比整车汽车悬挂系统动力学模型有限单元法降维求解和本文所提方法可得，两种方法都可以保证整车汽车悬挂系统动力学建模的精度，在求解精度相差不大的情况下，本文所提方法的计算效率明显要高得多。

为了进一步研究汽车悬挂系统运动中产生的弹性变形对汽车悬挂系统振动的影响，以后的研究通过建立更高维的柔性多体动力学模型和研究其振动频率进行分析。