基于分数阶微积分的油气悬架建模与试验分析

孙会来，金 纯，张文明，李 昊，田海勇

（1．北京科技大学 机械工程学院，北京 100083；2．中国北车股份有限公司，北京 100083）

油气悬架具有良好的非线性阻尼特性，结构简单，性能可靠，其特性对车辆行驶平顺性和操纵稳定性有重要影响，在大客车、载重货车以及工程机械车辆上被广泛采用，特别是在载重自卸车领域应用更加普遍。针对油气悬架建模及其运动特性国内外学者做出了大量研究，Concordia大学的研究人员［1］先后建立了单筒单气室油气悬架缸的参数化计算模型和非线性模型，采用基于能量的频率相关等效线性化技术，其等效线性参数化模型被广泛采用；Gao等［2］对重型车辆主动油气悬架进行了建模及仿真；孙涛等［3］建立了油气悬架长通孔紊流阻尼模型；杨波等［4］基于悬架非线性阻尼对双气室悬架进行了建模研究。在上述研究中建立的油气悬架数学模型在考虑油液可压缩的特性时，多认为油液含气以后体积弹性模量发生变化。且同时做了气体的质量保持不变，无泄漏，无形变等多项假设［5－8］。而黄夏旭等［9］研究又表明气体溶解效应和油液可压缩性在油气悬架设计研究过程中不可忽略。

悬架系统的实际工作要更为复杂，不仅与本身参数有关，同时受到各处零部件相对运动位移、温度、磨损情况、油液黏度及其在油气悬架内流动及可压缩特性的影响。以上诸多研究中总是侧重于某几项参数的变化及对悬架系统特性的影响，或者是通过增加计算项来保证计算精度，同时假设过多计算复杂难以从总体实现准确快速描述。

传统整数阶导数模型需要通过引入多个导数项和材料参数，本文提出根据油气悬架多相介质的力学特点，引入分数阶微积分理论对运动微分方程中的阻尼损耗项进行修正建模并分析其等效黏弹性，建立某型工程车辆油气悬架的振动分数阶振动方程并进行数值求解，探讨在该模型及不同工况下悬架的振动特性。

1 油气悬架结构及原理

单气室油气悬架主要由缸筒以及活塞杆和活塞组件组成，油气悬架上端设有气室或储能器。油缸气室与缸筒上端通过浮动活塞隔开，整个悬架缸内形成工作腔和环形腔，活塞杆壁上设有阻尼孔和单向阀。工作腔上部充入惰性气体，下部及环形腔充入油液。图1是根据上述描绘简化后的油气悬架物理模型。

图1 单气室油气悬架原理Fig．1 Single-chamber pneumatic suspension principle

图2 油气悬架实车测试Fig．2 Pneumatic suspension substantial vehicle test

当悬架处在压缩行程时，活塞杆向上运动，工作腔惰性气体体积减小，压力升高，油液受到压缩经单向阀和阻尼孔流入环形腔，使气室容积减小，气压力升高。在这此过程单向阀及阻尼孔会同时连通工作腔与环形腔，油液流量大，产生的油液阻尼力比较小，主要由气室内气体受到压缩产生弹性作用来抑制活塞杆的向上运动。当悬架处在伸张行程时，活塞杆向下运动，工作腔的体积增大进而压力减小，环形腔中的油液流向工作腔，此时单向阀关闭，油液只流经阻尼孔，流速较高，产生较大的阻尼力，抑制了活塞杆的向下运动，产生的阻尼力较大，有较强衰减振动的作用。工程车辆油气悬架安装如图2所示，油气弹簧悬架与轮胎组成振动系统起衰减振动、减缓冲击的作用。

经过如上分析及在试验中发现油液黏度及油液在油气悬架内流动、气体在油液中的溶解率的随温度和压力而变化等这些都不同程度的影响到了油气悬架特性。油气悬架阻尼特性等同于纯黏性体的计算精度有待提高，往往需要涉及过多参数。油气悬架特性和外力大小，温度的改变、力的作用时间及加载历史都有关系，这些符合黏弹性物质的力学特点。黏弹性材料的力学性质既表现出弹性，同时又有黏性特性，介于理想弹性体和牛顿流体之间，其应力－应变响应依赖于时间和应变率有记忆性，与荷载和变形历史有关。而分数阶黏弹性微积分模型仅用较少的参数便能准确地描述大量复杂黏弹性材料广泛频率范围内的动力学特性［10］，能较准确地描述黏弹性材料的动力学特性。

基于以上思想本文将油气悬架等效为黏弹性系统进行建模，引入分数阶模型描述阻尼损耗项对油气悬架系统进行分析，并进一步确定分数阶次，观察在不同激励下分数阶建模与整数阶区别。

2 油气悬架分数阶模型

2．1 分数阶微积分定义

分数阶微积分即非整数阶微积分，数学家从不同的角度入手，给出了分数阶微积分不同的定义。每种定义有其表达式与性质，其中被多数人所采纳的定义形式是 Riemann-Liouville（R-L）分数阶微积分如下［11］：

式中：Γ为伽马函数；n为任意整数；aDqt为分数阶微分算子；q为分数阶阶次；a为积分算子下限；t为积分算子上限。

在工程中应用较多及本文采用的是更为方便的Caputo定义，这种计算形式避免了极限求解，且初始条件以整数阶微积分的形式给出，其定义如下：

（R－L）分数阶微积分定义与Caputo分数阶微积分定义都是在时域范围内对函数f（x）的定义。（RL）的Laplace变换涉及到分数阶积分的初始值与原（R－L）分数微分的初始值，虽可求出含（R－L）分数阶次导数的解，但很难解赋予合理的物理意义。而Caputo分数阶微积分定义的最大优点是其初始值与整数阶微积分一样，有明确的物理意义，具有超奇异性。

2．2 油气悬架的分数阶模型

如前所述，影响悬架输出特性的因素众多，试验条件也无法做到规避所有干扰项，建立数学模型和模拟分析前，都需根据具体情况做出相应假设。本文研究油气悬架阻尼力等效分数阶黏弹性的表现，可借助筒式液力减振器的研究经验［12］，对工作参数做了适当的简化和忽略，做出如下假设：

（1）悬架中活塞、活塞杆与缸筒之间不产生泄漏；

（2）液压油的重力远小于其工作情况下的运动力，故不计工作油液的重力势能；

（3）不计压力变化而引起系统构件的弹性变形；

（4）被研究封闭区域内同一瞬时压力处处相等；

（5）不考虑油气悬架内部温度变化对性能的影响，即假设其工作工程中是恒温过程。

油气悬架的受力分析对象为活塞杆筒组件，活塞杆在内部做往复运动。油气悬架存在压缩运动和拉伸运动，这两种运动中流体的运动方向是相反，因此在对应的两种运动状态下对活塞杆筒组件进行力平衡分析。油气悬架在受到外部激励时，活塞杆及活塞组件和缸筒之间要产生相对运动。从输出力的角度出发整体来分析其振动特点。油气悬架在受到外部激励时，活塞杆及活塞组件和缸筒之间产生相对运动时，其输出作用力主要由气体弹性力、油液阻尼力及活塞和缸筒之间的摩擦力组成。输出合力公式如下：

式中：F为油气悬架缸筒输出力；Fk为气体弹性力项；Fc为油液阻尼力项；Ff为活塞和缸筒之间密封件的摩擦力项

对于弹性力、阻尼力及摩擦力项的计算，先取基本整数阶常用计算方法［13］进行说明，弹性力项为：

式中：P为工作腔气体压力 ；A3为活塞杆面积；M为流体密度；γ为气体多变指数；V0为静平衡位置气体体积；x为活塞和缸筒间相对位移。

油气悬架密封件的摩擦力使得油气悬架的动静太特性呈现出迟滞性，物理特性类似阻尼可以借助液压缸及活塞杆的密封件摩擦阻力的算法来计算：

式中：μ1为密封圈摩擦系数；ΔP为密封圈两侧压力差；D为液压缸缸筒内径d为活塞杆外径；bD，bd为活塞、活塞杆密封圈宽度；kDkd为密封圈摩擦修正系数。

油气悬架阻尼力项根据悬架阻尼孔的节流孔口形式，选取薄壁孔口类型的阻尼流量计算方法如下：

式中：Cz，Az为阻尼孔的流量系数、过流面积；Cd，Ad为工作腔气体压力；ρ为油液密度。

sign为符号函数，根据输入信号的速度方向取值，压缩行程为1，伸张行程－1。

先只考虑压缩行程即sigx＝1。联立方程（4）～（6）可以得到悬架输出力整数阶方程，便于计算进行无量 纲 化 处 理，令 a1＝ （A1－A3）3r／［2n2（CzAz＋CdAd）2］，a2＝Mg·Vγ0，a3＝V0，a4＝A3，得到简化系统方程：

根据图1运行原理及其几何关系可知：

式中：z（t）为悬架上端输出位移；y（t）为悬架下方输入激励位移。

上式描绘了三者位移及加速度关系，在油气悬架的下方施加一个正弦位移激励 y（t）＝A sinωt作为输入，相当于路面激励，A为激励振幅，ω为激励频率。

联立式（3）～（9），在整数阶方程基础上对于耗散阻尼项引入分数阶参数，同时计入自身重力，建立系统等效分数阶 Bagley-Torvik阻尼振动方程［14］，进一步得到油气悬架输出力的非线性分数阶方程：

式中：D为分数阶微分算子；α为分数阶算子阶次。

方程（10）为分数阶非线性常微分方程，弹性力项中变量x位于分母上，阻尼力项亦为平方是非线性项，这使得悬架的振动输出力是非线性的。大部分分数阶微分方程难以求得解析解，此类分数阶微分方程的解析解也只能由比较特殊的函数来表示，针对工程需求一般求取此类方程的数值解。

目前分数阶微积分常用的计算方法包括解析法、数值法和滤波器算法。本文采用改进的Oustaloup滤波器［15］算法来仿真未知信号的分数阶微积分，其原理是用高阶的传递函数近似替代分数阶算子。Oustaloup算法就是在频率段内，用泰勒展开并取一阶近似得到传递函数。另一方面，考虑到滤波器分子和分母阶次的一致性会导致仿真过程出现代数环，故而进行改进，后面再接一个低通滤波器。这样可编写滤波器零极点及增益函数，在Simulink中建立起如图3的分数阶微分器模块。进而搭建悬架系统微分方程，求得数值解。

图3 分数阶微分滤波器Fig．3 Fractional differential filter

3 悬架分数阶模型数值求解和试验

3．1 分数阶求解及试验对比

自行设计油气悬架试验台能获取油气悬架输出力特性的动态数据，方便于对不同工况下的工作特性进行分析和比较，验证油气悬挂分数阶模型的正确性。文中试验使用的试验台其结构图如4所示，主要由驱动装置、悬架系统、传感器、采集系统和簧上负载组成。设计思路主要是参照1／4车的振动模型，模拟车辆实际运行过程中油气悬架的工作状态，测试其受到外部激励和载荷时的输出特性。

整个台架驱动力由下端地基平台上的变频三相异步电动机发出，经减速器后将力传导至与支架接触的偏心轮上。偏心轮与传动梁保持接触，传动梁通过关节轴承与悬架活塞杆铰接。偏心轮随电机匀速旋转同时驱动传动梁做上下方向的摆动，摆动传递至悬架缸活塞杆与传动梁铰接处按照几何关系即为竖直方向的正弦激励。为防止试验时质量块左右晃动，在质量块两侧设计了限位支架。本设计激励频率可通过调节电机频率来随时改变；激励幅值可通过调节偏心轮的偏心距来方便的调整；悬架簧上质量可通过改变质量块数量来调整大小。

对于横梁臂自身质量及质量块质量需通过杠杆比进行等效，即为簧上质量部分，通过气室静态气压推算簧上质量予以校正。试验台系统构件引起的弹性变形以及系统自身各处阻尼较小可以忽略不计。试验台下放为钢结构试验台地基，地基自身振动引发的误差极小也可忽略不计。每次试验前需对偏心距进行校正，保证激励振幅准确性。传动梁及横臂梁的杠杆比问题，需首先参照几何关系进行等效转换和试验校正。

图4 油气悬架试验台Fig．4 Test stand of hydro-pneumatic suspension

为验证分数阶理论的准确性与整数阶对比，完整描述试验过程及结果，本文同时根据试验台结构设计构建了台架三维模型如图5，进而将模型导入多体动力学软件ADAMS／View中，通过设置各零部件的材料属性并添加零部件约束及力关系，与LMSAmesim液压软件联立得到试验台多刚体整数阶仿真模型，目的是作为整数阶结果与论文的分数阶模型进行对比。

试验选用悬架基本参数如表1所示，为某吨位自卸车实车悬架。研究主要对比相同工况下悬架缸的输出特性并针对不同工况展开试验，即簧上载荷的振动加速度、悬架缸输出力幅值和频率、运动行程等。

图5 试验台架软件仿真Fig．5 Software simulation of test stand

表1 试验用悬架基本参数Tab．1 The basic parameters of test suspension

本文针对偏心轮偏心距为5 mm，转动频率为1 Hz正弦激励情况先进行分析，这也是工程车辆常见激励频率。取多个不同分数阶次于振动稳态后在时域范围内进行对比分析发现，与整数阶次结果的周期和相位基本一致，而幅值差异明显，可见分数阶次的选取对振动结果的描述有很大影响。为了最为吻合试验数据，分析前述方程，将最优阶次的选取转化为单变量的最优数学问题［16］，单变量为分数阶次α，根据计算模型与实测数据间误差最小原则进行数值求解选取最优阶次。定义误差评价函数E为计算结果与实测数据的差方均值，这样更为直观方便，取毫米为单位。如下：

式中：E为误差评定；n为对比采样点个数；xli为整数阶结果第i个采样点位移；xi为实测结果第i个采样点位移。

同样x2i为分数阶结果中第i个采样点位移，代入后则为分数阶计算结果与实际测量误差。约束条件的选取设在频段［0．4，1．4］内进行，迭代步长为 0．001，观察最优阶次变化后目标函数的变化如图7。迭代求解得分数阶建模在分数阶次取0．752～1的范围内比整数阶结果误差要小，在1处即为整数阶结果同为1．68的差方均值。迭代结果最小误差即最优值出现在0．912处，其幅值结果最接近实测结果。下文以0．912次分数阶结果、整数阶次仿真结果、实测数据进一步对比分析。

图6 差方均值随迭代阶次变化Fig．6 Error changes with the order

实测数据是悬架稳定振动后测得，为了便于说明对比图的绘制中，实测数据相位滞后或提前90°，取Matlab分数阶数值结果、ADAMS整数阶结果及实测结果合并对比，如图7为悬架缸缸筒相对位移对比，同取在4 s时域范围内对比，后续分析也如是。

分析图7可知缸筒与活塞杆的相对位移中分数阶建模结果更加贴近于实际情况，整数阶结果振幅比试验数据振幅要更大，最大振幅要超过实测值约1．5 mm。在分数阶原理上，对于阻尼项阶数与运动关系可知，随着分数阶次的增大振动会越来越激烈，振动位移幅度也会越大，即前面已有的运动状态对下一时刻影响力减弱。反之越小于1，则增加这种相互关系，使得阻尼力变大，振幅减小。明显分数阶建模在幅值方面更为精确，尤其是在位移下止点位置。受阻尼力影响，活塞压缩和伸张行程是一个不重合过程。在伸张行程，在阻尼的作用下使得其伸张行程要更小，分数阶的引入强调了阻尼的非线性，分数阶模型结果在伸张后半程和压缩前半程的优势更为明显，最大值最小值的出现时间与实测结果更为相近，改变趋势更为吻合。

图7 缸筒相对位移对比Fig．7 data of displacement curve comparison diagram

图8 缸筒加速度对比Fig．8 Accelerated speed of cylinder comparison

图9 悬架缸输出力曲线对比Fig．9 Output force of suspension cylinders comparison

图8 为悬架缸缸筒加速度曲线，与图6相似，频率和幅值都很接近，明显仿真结果的曲线比试验结果更为平滑，这是因为在整数阶仿真模型及软件中把试验台当作刚体处理，没有考虑传动梁等构件自身的阻尼及多处振动，但其对结果影响不大。

图9为悬架缸输出力曲线，动平衡后油气悬架的输出力呈周期性变化，结果曲线幅值变化频率与实测数据相近。油气悬架本身的高度耦合多复杂工况使得黏弹性特性出现且最终趋向稳态。在压缩行程分数阶模型优势并不明显，而在伸张行程分数阶建模远比整数阶结果更贴近实测数据，在伸张行程主要是阻尼力在起作用，分数阶模型的引入是针对阻尼力进行了调整，其优势在伸张行程体现明显，阻尼力的调整同样使得输出力较之于整数阶相位上有延迟，其输出力最大值出现要晚，这与实测数据更为相近。

3．2 多工况下分数阶次与误差变化

为了进一步验证分数阶建模在油气悬架多工况中具有广泛的有效性，可改变激励振幅和频率，观察最优分数阶次的变化，同时对比其与整数阶建模误差变化。本文通过控制偏心轮距和激励电机的转速改变激励频率和振幅，进行多次试验，并参照分数阶结果差方均值最小为目标进行最优阶次的确定。在同一振幅5 mm改变激励频率，在同一频率1 Hz下改变振幅分别进行最优化求解和试验，通过定量计算拟合曲线得到图10最优阶次随激励频率变化以及图11随激励振幅变化。

分析图10和11发现随着频率的改变最优阶次下降并不大，高频激励下最终稳定趋于0．9，而随着激励振幅的变化下降较大趋于稳定0．86，这是因为在高幅快速振动下阻尼力起到了更为明显的作用，其等效黏弹性明显，二者最后都趋于稳定，这说明分数阶建模的取值在高幅阶段可以稳定选取，建模方便直接。

图10 最优阶次随激励频率变化Fig．10 Optimal order changes with the excitation frequency

图11 最优阶次随激励振幅变化Fig．11 Optimal order changes with the excitation amplitude

图12 不同振幅下差方均值对比Fig．12 Error comparison at different amplitude

为了进一步分析和观察分数阶最优阶次确定后与整数阶建模的差别，可对二者按评价指标进行观察。取多组不同振幅下数据进行评价然后拟合绘图得到图12即不同振幅激励下差方均值变化，差方的平均值越小说明精确度越高，可见随着振幅增大二者误差都有进一步的增大，但是整数阶误差增大更为明显，进一步验证了分数阶建模的有效性。另一方面，取多组不同频率下试验数据拟合可得不同频率下差方均值对比，因最优阶次下降不大，二者误差并不明显，但在总体范围内变频和变幅值的分数阶建模误差都要小于整数阶。结合图11可知等效黏弹性的分数阶建模方式能进行有效修正，得到更准确的结果。

图13 阻尼力对比Fig．13 Damping force comparison

分析分数阶建模的特点可以对阻尼力项进行说明，图13为初始试验条件即1 Hz和5 mm激励下二者数值计算的阻尼力项对比，在伸张行程分数阶可以达到稳定的550 N，而整数阶结果只有400 N，分数阶阻尼对系统消耗更大，而二者弹性力项变化不大，原理上这造成了悬架的相对位移非标准正弦波，在阻尼力的非线性作用下使得向下运动的伸张行程偏大，向上运动的压缩行程偏小。

4 结 论

（1）根据油气悬架工作原理，突出油气悬架的多相介质特点，采用等效黏弹性阻尼的分析方法，引入分数阶导数建立了油气悬架的 Bagley-Torvik方程。利用Oustaloup算法对分数阶非线性微分方程进行数值求解得到油气悬架分数阶输出力，位移，加速度等多项数值结果。

（2）将分数阶、整数阶仿真与试验结果进行对比，验证了分数阶模型的正确性。在1 Hz，5 mm激励下测试表明取分数阶次0．912时与实验结果最为相近，该简洁模型能较精确地描述油气悬架的特性。

（3）改变激励振幅和频率，在不同工况下进行对比分析最优阶次的变化，在高频激励下最优阶次趋于0．9，在高幅激励下趋于0．86，在整个测试激励范围内分数阶模型误差都要小于整数阶，这说明了分数阶建模的有效性。

［1］Pazooki A，Cao Dong-pu，Rakheja S．Experimental and analytical evaluations of a torsio-elastic suspension for off-road vehicles［J］．SAE International Journal of Materials and Manufacturing，2010，3（1）：326－338．

［2］Gao B，Darling J，Tilley D，et al．Control of a hydropneumatic active suspension based on a non－linear quarter-car model［J］．Journal of Systems＆Control Engineering，2006，220（1）：15－31．

［3］孙涛，喻凡，张振东．重载车辆油气悬架不确定分析及H∞控制器设计［J］．振动与冲击，2007，26（9）：51－54．SUN Tao，YU Fan，ZHANGZhen-dong．Uncertainty analysis for a hydro-pneumatic suspension system of heavy duty vehicles and H∞controller design［J］．Journal of Vibration and Shock，2007，26（9）：51－54．

［4］杨波，陈思忠，王勋，等．双气室油气悬架特性研究［J］．机械工程学报，2009，45（5）：276－280．YANG Bo，CHEN Sizhong，WANG Xun，et al．Research of Twin accumulator Hydro－pneumatic Suspension［J］．Journal of Mechanical Engineerin，2009，45（5）：276－280．

［5］王增全，申焱华，杨珏．连通式油气悬架数学模型及特性分析［J］．农业工程学报，2012，28（5）：60－65．WANG Zeng-quan，SHEN Yan-hua，YANG Jue．Mathematical model and characteristics analysis of interconnected hydropneumatic suspension［J］．Transactions of the Chinese Society of Agricultural Engineering，2012，28（5）：60－65．

［6］徐道临，张林，周加喜．重型矿用自卸车油气悬架参数优化［J］．振动与冲击，2012，31（24）：98－101．XU Dao-lin， ZHANG Lin， ZHOU Jia-xi． Parametric optimization of hydro pneumatic suspension of a heavy mining dumper［J］．Journal of Vibration and Shock，2012，31（24）：98－101．

［7］黄夏旭，申焱华，杨珏，等．基于集中参数热模型法的自卸车油气悬架系统热分析［J］．农业工程学报，2013，29（10）：64－70．HUANG Xia-xu， SHEN Yan-hua， YANG Jue， et al．Thermal analysis of hydro-pneumatic suspension system for dumper based on a lumped-parameter thermal model［J］．Transactions of the Chinese Society of Agricultural Engineering，2013，29（10）：64－70．

［8］张辉，何景峰，韩俊伟．油气悬架加载试验台正弦波复现控制方法［J］．农业机械学报，2009，40（7）：35－39．ZHANG Hui，HE Jing-feng， HAN Jun-wei． Sine wave replication control methodon hydro pneumatic suspension force loading system［J］．Transactions of the Chinese Society for Agricultural Machinery，2009，40（7）：35－39．

［9］黄夏旭，杨珏，申焱华，等．基于气体溶解与油液可压缩性的油气悬架性能研究［J］．农业机械学报，2013，44（6）：14－18．HANG Xia-xu， YANG Jue， SHEN Yan-hua， et al．Characteristics analysis of hydro pneumatic suspension based on gas dissolution and oil compressibility［J］．Transactions of the Chinese Society for Agricultural Machinery， 2013，44（6）：14－18．

［10］ Song D Y，J TQ．Study on the constitutive equation with fractonal derivative for viscoelastic fluids－modified Jeffreys models and its application［J］．Rhelo Acta 1998，21（37）：512－517．

［11］Connolly J A．The numerical solution of fractional and distributed order differential equations［D］．University of LIverpool，2004．

［12］ Kurfess，Thomas R，Messman，et al．A methodology to integrate a nonlinear shock absorber dynamics into a vehicle model for system identification［J］．SAE International Journal of Materials and Manufacturing，2011，4（1）：527－534．

［13］ Wilkinson，Peter M，Van Dierendonck，et al．Theoretical model for the influence of gas properties and pressure on single-bubble formation at an orifice［J］．Source：Chemical Engineering Science，1994，9（49）：1429－1438．

［14］沈淑君，刘发旺．解分数阶Bagley-Torvik方程的一种计算有效的数值方法［J］．厦门大学学报（自然科学版），2004，7（3）：306－311．SHEN Shu-jun， LIU Fa-wang． A computational efficient numerical method of solving fractional Bagley－Torvik equation［J］．Journal of Xiamen University（Natural Science）2004，7（3）：306－311．

［15］薛定宇．高等应用数学问题的MATLAB求解［M］．北京：清华大学出版社，2008：398－399．

［16］刘鲭洁，陈桂明，杨旗．基于Matlab工具的遗传算法求解有约束最优化问题［J］．兵工自动化，2008，27（11）：43－44．LIU Qing-jie， CHEN Gui-ming， YANG Qi． Solving constrained optimization through genetic algorithm based on matlab toolbox ［J］． Advanced Manufacture and Management，2008，27（11）：43－44．