基于非线性有限元的空气弹簧垂向刚度分析

邓长喜

【摘 要】针对传统计算空气弹簧垂向刚度方法的局限性，本文利用有限元法分析了空气弹簧的垂向刚度。以某型空气弹簧为例，进行了充气过程模拟和刚度计算，并对比了仿真分析和试验所得刚度结果，结果表明仿真分析可对空气弹簧的充气过程以及刚度计算进行有效模拟。

【关键词】空气弹簧；垂向刚度；有限元法

Analysis of Air-spring Stiffness by Nonlinear Finite Element Method

DENG Chang-xi

（Avic Xian Aircraft Industry （Group） Company， Yanliang Shaanxi 710089， China）

【Abstract】Finite element method（FEM） was introduced to analyze vertical stiffness of air spring because of traditional methods limitation. Taking a rubber air spring as an example， inflation process and stiffness calculation were performed. And according to comparison of results between numerical simulation and experimental method， it shows that the numerical simulation method is feasible and effective for simulating the inflation process of air spring， and the vertical stiffness calculation.

【Key words】Air spring； Vertical stiffness； Finite element method

0 引言

空气弹簧诞生于19世纪中叶，有专利记载在1847年John Lewis申请了空气弹簧的发明专利[1]，同年在《美国科学》的创刊号上提出了“ride on air”的概念。空气弹簧以其固有频率低、水平稳定性好、承受载荷能力范围大、无需另外设置阻尼器等优点，作为隔振器在载重汽车、有轨机车、超精密机床等领域内被广泛采用[2-3]。然而，在实际工程应用中空气弹簧的相关计算较繁琐，且忽略了众多影响因素，设计时对影响空气弹簧主要性能的工作压力、附加气室以及阻尼的选取原则不太明确，使其不能发挥应有的最大作用。

目前，确定空气弹簧刚度参数主要有理论分析和实验测定两种途径。对于形状简单的空气弹簧，理论分析可以用解析法或图解法计算出近似解，但对于形状复杂的空气弹簧，要计算其刚度就比较困难。而且，理论分析得出的解析解并没有考虑空气弹簧囊壁材料对刚度的影响，实验测定虽可以得出实际刚度大小，但试验过程需要有样品，而且试验过程参数的调节比较困难，这就给空气弹簧的设计带来了较大的难度，不仅延长了设计周期，而且增加了研发费用[4]。所幸随着计算机辅助工程的快速发展和应用，利用仿真分析进行空气弹簧设计和优化计算逐步成为不可或缺的有效手段之一，用有限元方法分析空气弹簧的刚度，不仅考虑了囊壁材料对刚度的影响，而且对于形状复杂的空气弹簧，也能较好地计算其刚度。自20世纪90年代以来，空气弹簧系统的仿真模拟研究成为热点，国外众多学者利用有限元对空气弹簧进行了一系列模拟研究得到了很多重要结论[5-6]。国内空气弹簧的研究始于1957年，当时只是局限于车辆用的空气弹簧，对空气弹簧的理论进行了初步研究，在理论和试验方面，积累了一些经验[7]。

本文分析了工程中对空气弹簧垂向刚度的简化计算，随后利用仿真分析手段充分考虑材料、几何以及接触非线性的影响，对某型空气弹簧进行非线性有限元分析，计算其垂向刚度，通过与试验结果进行对比，表明仿真分析结果具有较高精度，可用于指导空气弹簧工程设计和优化。

1 空气弹簧的刚度特性

空气弹簧是一种由有帘线的橡胶囊和充入其内腔的压缩空气所组成的非金属弹性元件，其隔振性能主要受到刚度特性、频率特性和阻尼特性的影响。尤其是其刚度特性，一致是研究热点之一，结构设计中，往往先根据初步的刚度需求进行设计。

设空气弹簧胶囊内的气体为理想气体，空气弹簧受到垂直载荷F的作用，此时内部的绝对气压为P，外部气压为Pa，则有以下关系式：

F=（P-Pa）Ae（1）

PVn=P0Vn 0=const（2）

式中：Ae——空气弹簧的有效面积；

n——多变指数，与气体的流动速度相关，缓慢流动时相当于等温过程，n=1；快速流动时接近绝热过程，n=1.4，一般情况下，n=1.3～1.38；

P0、V0——为空气弹簧处于静平衡位置的内部气体的绝对压力和容积。

由式（1）、式（2）可得：

F=（P0Vn 0-Pa）Ae（3）

将式（3）对空气弹簧的垂直位移S求导数，其中Ae=-dV/dS，负号表示空气弹簧受载荷压缩时容积减少，则可得空气弹簧的垂直刚度公式：

K=（P0Vn 0Nn-Pa）（dAe/dS）-Aen P0Vn 0Nn+1（dV/dS）（4）

由上述可知空气弹簧刚度的影响因素有两个，容积变化率和有效面积变化率。通过添加附加气室可以改变容积变化率；通过改变空气弹簧的结构可以控制有效面积变化率，则可以进一步降低空气弹簧的刚度。显然，常规空气弹簧刚度计算中，将问题进行了诸多简化，可用于初期的空气弹簧设计，但并未考虑实际结构和工况中的诸多因素影响，与实际情况有一定偏差，利用有限元分析进行刚度求解是十分必要的。

2 有限元模型

2.1 空气弹簧刚度分析中的非线性

空气弹簧的性能分析过程中，包含有几何非线性、材料非线性、边界非线性三大非线性问题，即为一个三重耦合的强非线性问题。空气弹簧的气囊和橡胶堆的材料是典型的超弹性材料，即橡胶，其计算就涉及到了材料的非线性-超弹性材料本构，本文中空气弹簧的材料模型选择M-R模型[8]：

W=C1（I1-3）+C2（I2-3）（5）

其中C1，C2的取值依赖于实验数据而得。对此模型，既可运用大变形非线性弹性有限元法中的全拉格朗日模式来求解，也可用更新拉格朗日模式来求解。气囊和橡胶堆在模拟的过程中变形是很大的，传统有限元计算的小变形情况的理论就不再适用了，因而必须考虑有限元计算中的几何非线性问题。而空气弹簧的接触非线性问题则是由上下盖板与气囊相接触引起的，当气囊与盖板相接触，它们相互之间会产生法向压力和摩擦剪切力，在接触有限元法中可用罚单元法来求解。本文主要的分析目的就是确定空气弹簧对应的载荷-位移响应，故本文要求解的方程组属求解非线性结构的响应类型。一般地应用于非线性有限元方程的非线性结构问题可表示为：

K（q）q=g（6）

其中q是未知位移的矢量，g是施加的结点载荷矢量，三大非线性因素则隐含于内力K（q）q中，而内力由内应力σ的空间积分得到：

其中B是结构刚度矩阵。一般的几何非线性有限元法中的应变与位移的非线性关系，材料非线性有限元法中的应力与应变的非线性关系接触非线性问题中，因载荷与节点位移有关，最后均会导致整体刚度方程为非线性方程。求解非线性方程一般都采用线性化方法，即把非线性线性问题转化为一序列线性问题求解。本文分析过程中是通过逐步施加给定的位移，以增量形式趋于最终解而得到结果，因此ABAQUS将模拟计算分为许多位移增量步，并在每个位移增量步结束时确定近似的平衡，所有增量响应和就是非线性分析的近似解[8]。

2.2 模型及载荷边界条件

根据空气弹簧的轴对称特性，建立如图1所示几何模型，考虑帘线等结构对气囊的增强。

Fig.1 Axisymmetric model of air spring

结构网格模型如图2所示，并考虑接触分析特点对接触区网格进行了网格设计，气囊内表面设置静水流体单元，如图3所示（为便于观察，图中给出的是轴对称模型旋转180°后的效果），用于表征内部气体域边界并实现内部气压控制。

载荷及边界条件设置为上板顶面设置强制位移（用于计算垂向刚度），空气弹簧底座下端面固定，左侧边界为轴对称条件，在模型可能接触的区域均设置接触关系。分析工况为，首先固定上板顶面和底座，给气囊内冲压0.2MPa，然后封闭空腔，上板下移40mm，随后继续充压，使上板反力达到要求的值，为计算此时的刚度，上板继续下移10mm。

材料参数如表1所示：

2.3 气囊充压及上板下移过程的实现

为实现气囊充压过程以及垂向刚度的模拟仿真，在气囊内表面，利用共节点技术生成静水流体单元，如图3所示，单元类型为FAX 2，将这些节点利用耦合的方式，耦合到一个参考点，对于轴对称模型，通常选择对称轴上的一个点为参考点。由于设置了耦合，这些静水流体单元的表征就通过唯一的参考点实现，参考点的自由度8即为气囊的内压，通过设定该值即可实现气囊充压的模拟，同时，由于采用了共节点技术，气囊结构就会随着内压的变化而变化。于此类似，当上板下移时，由于气囊结构会发生变化，内压也会随之发生变化，以此实现垂向刚度的求解。

3 结果

计算刚度前的预充压状态应力云图如图4所示，最大Mises等效应力74.4MPa，上板进一步下移10mm后等效应力79.5MPa，如图5所示，且空气弹簧下端的橡胶堆进一步受压，显然，随着上板的下移，结构应力增大，变形增大。

Fig.4 Von Mises stress before stiffness calculation

Fig.5 Von Mises stress after stiffness calculation

图6为空气弹簧上板下移10mm过程中气囊内压的变化曲线，内压由0.407MPa，增大到0.424MPa，内压的变化过程与上板位移近似呈线性关系。

Fig.6 Von Mises stress after stiffness calculation

图为空气弹簧上板下移10mm过程中气囊容积的变化曲线，容积由3.69E7mm3，减小到3.609E7mm3，容积的变化过程基本与上板下移同样近似呈线性关系。

Fig.7 Von Mises stress after stiffness calculation

图8为仿真分析和实验所得空气弹簧垂向刚度曲线（上板反力-位移曲线），对比实验和仿真分析结果可知，仿真所得刚度反力略低于实验所得反力，这主要是仿真分析中接触摩擦参数以及部分模型简化导致的，但整体来看仿真分析和实验所得刚度符合较好。

Fig.8 Displacement-reaction force relationship for air spring

4 结论

1）传统工程算法计算空气弹簧刚度方法简单，可用于初步的空气弹簧设计；

2）空气弹簧等效应力和内压随上板下移进一步增大，气囊容积随上板下移进一步减小，内压增大和容积减小与上板下移位移近似呈现线性关系；

3）仿真所得空气弹簧刚度与试验所得结果匹配较好，表明利用仿真分析空气弹簧刚度可用于空气弹簧的精细设计。

【参考文献】

[1]John Lewis. Car spring[P]. US Patent. No. 4965， 1847-02-10.

[2]李蒂，付茂海，黄运华.空气弹簧动力学特性参数分析[J].西南交通大学学报， 2003，38（3）：276-281.

[3]盖玉先，董申，李旦，等.超精密机床的振动混合控制[J].中国机械工程，2000，11（3）：289-291.

[4]吴善跃，黄映云.空气弹簧刚度的有限元分析[J].海军工程大学学报，2001，13（6）：94-98.

[5]Giuseppe Quaglia， Massimo Rorli. Air suspension dimensionless analysis and design procedure[J].Vehicle System Dynamic，2001，35：443-475.

[6]Takuye Yuasa. The application of CAE in the development of air suspension beam[Z].SAE Inc， 1997. SAE 840342.

[7]张利国，张嘉钟，贾力萍，等.空气弹簧的现状及其发展[J].振动与冲击，2007， 26（2）：146-151.

[8]卡德斯图赛H.有限元法手册[M].北京：科学出版社，1996.