内曲线径向柱塞马达多作用的特点使之具有功重比大、脉动小等优点，因此被广泛应用于工程机械、军工、船舶等的行走及驱动系统中[1]．它也是目前国内外应用较多的一种低速大扭矩液压马达[2]．其中，凸轮环作为内曲线径向柱塞马达的关键零件，其导轨曲线直接决定着柱塞的运动特性．此外，导轨曲线的设计不合理还会带来冲击和噪声，过大的冲击将加速导轨的磨损．因此，导轨曲线在很大程度上决定着马达的输出特性、噪声和寿命[3]．近年来，国内研究人员对不同的导轨曲线进行了优化，而国外未公开相关的研究成果．其中，文献[4-6]通过对等接触应力曲线进行数值求解，利用多项式对结果进行拟合，进而得到导轨曲线的方程形式．由于他们在最大接触应力的推导过程中忽略了惯性力和摩擦力，因此所设计的曲线并非真正意义的等接触应力曲线．以文献[5]给出的曲线为例，该曲线在与零速区的衔接位置存在硬冲击，进而会产生冲击力．文献[7]最早提出等应力导轨曲线的设想．文献[8]基于该设想对无脉动等应力曲线进行了设计．设计过程中同样忽略了惯性力和摩擦力，且方程组求解困难．文献[1]对带过渡区的等加速曲线进行了设计，等加速曲线由于存在度加速度突变，因此存在多处软冲击．文献[9]提出了无冲击的修正心形线式导轨曲线，并对马达的脉动特性进行了分析，该曲线脉动较大．此外，还有学者对非柱塞式的内曲线马达的导轨曲线进行了分析设计．文献[10-11]对径向双球塞柱塞马达导轨的广义椭圆定子曲线进行设计，并对马达的输出特性进行了分析．文献[12]针对内曲线球塞式水液压马达，给出了14种常用的内曲线方程，并对比分析了14种内曲线的压力角、输出特性和接触应力．文献[13]针对内曲线球塞式水液压马达的导轨曲线进行了分析设计．综上所述，目前内曲线径向柱塞马达导轨曲线仍然存在着冲击或脉动大的问题，因此本研究在考虑惯性力和摩擦力的情况下，推导出轨道与滚柱最大接触应力的数学模型，并以其最小值作为适应度函数，以无冲击、低脉动作为约束条件，并利用遗传算法对高阶导轨曲线的系数进行多目标优化求解．最后以10阶、11阶导轨曲线为例，进行对比分析．1 导轨-滚柱最大接触应力模型最大接触应力与导轨和滚柱的曲率半径及导轨对滚柱的反作用力有关，而在柱塞组件的运动过程中，反作用力会随着度速度和度加速度的变化而变化，因此须要先对柱塞组件进行运动学分析．1.1　柱塞组件运动学分析在柱塞腔液压力的推动下，滚柱沿轨道运动，轨道对柱塞组件的反作用力推动缸体(转子)旋转．轨道曲线示意图如图1所示，图中：滚柱中心轨迹c-c为导轨的理论曲线，它真实反映了柱塞组件的运动规律，因此本研究所提及、求解的导轨曲线均指理论曲线；c-c的外包络线为导轨的实际轮廓曲线a-a；以导轨中心o为极点建立极坐标系，以极点o到滚柱中心o1的距离为极径ρ；ρ与起始位置的夹角为φ；φ0为零速区角度；φ1为工作幅角；γ为压力角；h为柱塞最大行程；Fn为轨道对滚柱的反作用力，Fn分解为切向力Ft和径向力Fr．由此可得柱塞组件的运动学模型．10.13245/j.hust.211005.F001图1凸轮环导轨曲线示意图柱塞组件的度速度为ν1=dρg/dφ=d(ρg+lg)/dφ=dρ/dφ，式中：ρg为柱塞组件质心极径；lg为质心到滚柱中心的距离，即lg=ρ-ρg，为常数．柱塞组件的度加速度为a1=dν1/dφ=d2ρ/dφ2．柱塞组件的相对速度为v2=dρg/dt=(dρ/dφ)(dφ/dt)=ωv1，式中ω为马达缸体(转子)的角速度．柱塞组件的相对加速度为a2=dν2dt=ωdv1dt=ωdv1dφdφdt=ω2a1．柱塞组件的牵连运动加速度及哥氏加速度分别为：a3=ρgω2；a4=2ων2=2ω2ν1．1.2　柱塞组件受力分析柱塞组件受力示意图如图2所示．柱塞组件在运动过程中所受力如下．10.13245/j.hust.211005.F002图2柱塞组件受力分析示意图a. 柱塞腔压力p产生的液压力Fp，有Fp=πd2p/4，其中d为柱塞直径．b. 柱塞组件相对缸体运动的加速度a2所引起的惯性力Fa，有Fa=ma2=mω2a1，其中m为柱塞组件质量．c. 柱塞组件随缸体旋转时产生的离心力F3．F3=ma3d. 马达角加速度ε引起的惯性力F4，有F4=mρgε．e. 牵连运动为缸体的转动，由哥氏加速度a4引起的惯性力F5，有F5=ma4=2mω2ν1．f. 轨道对滚柱的反作用力Fn和摩擦力Ff，Fn可分解为切向力Ft和径向力Fr，有：Ff=fqFn；Ft=Fnsin γ；(1)Fr=Fncos γ，(2)式中fq为滚柱与轨道之间的滚动摩擦系数．g. 柱塞与缸孔内壁的接触力F1和F2，以及F1和F2引起的摩擦力Ff1和Ff2．F1和F2分别是接触应力p1和p2产生的合力，p1max和p2max为p1和p2的最大值．考虑柱塞与缸孔之间的间隙较小，因此p1和p2可以看成是连续的呈直线分布的应力．有：Ff1=fF1；Ff2=fF2，其中f为柱塞与缸体之间的摩擦因数．h. 此外，柱塞组件还受到自身重力G、柱塞缸孔中油液的惯性力、油液黏性摩擦力等(图中未标注)作用．一般情况这些力较小可忽略不计[14]．1.3　求解反作用力Fn根据前述的运动学模型和力学模型求解最大接触应力的关键参数Fn．将所有的力向滚柱中心简化，得到力和力矩平衡方程组，即Ft-Ffcos γ-F4-F5+F2-F1=0；(3)-Fr-Ffsin γ+F3-Fa+Fp-Ff1-Ff2=0；(4)-FfR+F1(l-l0+(l0-l2)/3)-F2(l-l2/3)+Ff1d/2-Ff2d/2+(F4+F5)lg=0, (5)式中：l0为柱塞在缸体内的接触长度，接触长度随柱塞运动而变化，有l0=l00+ρ-ρmin，其中l00为初始接触长度(滚柱处于内死点)，ρmin为导轨曲线基圆半径(最小极径)；l2为p2的作用长度，则p1的作用长度为l0-l2；l为柱塞组件的名义长度；R为滚柱半径．有tan γ=dρ/(ρdφ)=ν1/ρ．由于有4个未知量(F1，F2，Fn，l2)、3个约束方程，因此还须增加一个方程．根据p1和p2分布的三角形相似原理，可得F1/F2=(l0-l2)2/l22．(6)将式(6)代入式(5)求解l2．为简化计算，做以下简化：因为在润滑状态良好的情况下，滚柱滚动时的摩擦阻力很小，且滚柱半径较小，所以FfR项可忽略．当正常运行时，内曲线马达的转速低(额定转速为170 r/min)，转速波动较小，且与极径相比，lg较小，因此(F4+F5)lg项可忽略．进而解得l2=(6l0l-4l02+3fdl0)/(12l+6fd-6l0)．(7)将式(1)、(2)和(6)代入(3)可得：F1=(Fnsin γ-Ffcos γ-F4-F5)(l0-l2)2l02-2l0l2；(8)F2=(Fnsin γ-Ffcos γ-F4-F5)l22l02-2l0l2．(9)将式(8)、(9)代入(4)得Fn=F3-Fa+Fp+f(F4+F5)Φcos γ+fqsin γ+fΦsin γ-ffqΦcos γ, (10)式中Φ为结构参数，有Φ=(l02-2l0l2+2l22)/(l02-2l0l2)．1.4　轨道-滚柱接触应力分析如图3所示，滚柱与轨道的接触为线接触，根据赫兹接触理论[15]，两接触弹性固体在接触区的线接触的最大接触应力σmax的表达式为σmax=2Fn/(πltb)=FnE*/(πR*lt)；(11)1/E*=(1-μ12)/E1+(1-μ22)/E2；(12)R*=RR1/(R1±R)，(13)式中：R*为当量曲率半径；E*为当量弹性模量；μ1，μ2和E1，E2分别为滚柱、导轨的泊松比和弹性模量；lt为线接触长度；R1为轨道接触点处的曲率半径(当轨道与滚柱外切时取正号)，为未知量；b为接触区的半宽度．10.13245/j.hust.211005.F003图3滚柱-导轨线接触应力分布示意图理论曲线上任意一点的曲率半径为R2=(ρ2+ρ'2)3/2/(ρ2+2ρ'2-ρρ″)．(14)当R20时，导轨曲线是凸的，滚柱与导轨内切，此时R1=R2+R，否则R1=-R2-R，代入式(13)可得当量曲率半径R*=-R(-R2-R)/R2=R(R2+R)/R2．(15)为保证滚柱与实际轨道曲线准确接触，理论曲线的曲率半径应满足R20或R2-R．(16)2 高阶导轨曲线求解无冲击要求导轨曲线的三阶导数存在且连续、无突变，而高阶曲线便可以满足这一要求，因此本研究选用高阶曲线作为优化设计的基础，并以上述建立的最大接触应力模型的最小值作为适应度函数(优化目标)，以无冲击、低脉动为约束条件，利用遗传算法对曲线的系数进行求解，最后得到无冲击、低脉动、低应力的高阶导轨曲线．2.1　高阶导轨曲线的一般表达式图4为分段高阶导轨曲线示意图，图中：ϕ为无量纲角度；ρ1(ϕ)为加速区曲线；ρ2(ϕ)为减速区曲线．为了避免柱塞在进油区和回油区之间切换而引起困油现象，如图4所示，在导轨的内死点和外死点位置设置零速区，即此时柱塞的相对速度为v2=0 m/s，柱塞相对缸孔静止．10.13245/j.hust.211005.F004图4分段高阶导轨曲线示意图内死点零速区的导轨曲线表达式为ρ(ϕ)=ρmin，外死点零速区的导轨曲线表达式为ρ(ϕ)=ρmin+h．零速区以外的部分称为工作区．工作区的高阶导轨曲线即为待优化对象．其一般表达式为ρ(ϕ)=ρmin+h(A0+A1ϕ+A2ϕ2+⋯+Anϕn)；ϕ=(φ-φ0)/φ1，式中：n为曲线的最高阶次；φ0=π/180 rad；φ1=2π/x-2φ0；x=8为内曲线作用次数；0≤ϕ≤1；A0，A1，…，An为n阶曲线的各项系数，同时也是待优化求解的系数变量．2.2　约束条件改变曲线的阶次及系数A0，A1，…，An可得到无穷条曲线．根据无冲击、低脉动的约束条件及低应力的优化目标对曲线系数建立约束方程，进而从无穷条曲线中选出符合要求、性能优良的曲线．如图4所示，柱塞组件从内死点零速区运动到外死点零速区的过程中经历了先加速后减速的过程，因此可以把工作区从中心位置分为加速区(图中红色虚线)和减速区(图中绿色虚线)，加速区和减速区分别建立曲线方程和约束方程．设加速区导轨曲线ρ1(ϕ)的各项系数为A1_0，A1_1，…，A1_n；减速区导轨曲线ρ2(ϕ)的各项系数为A2_0，A2_1，…，A2_n．对导轨曲线分段设计求解可以保证加速区和减速区对称，以减小马达脉动．为保证导轨曲线之间的衔接，应满足：ρ1(0)=ρmin;ρ1(0.5)=0.5h+ρmin;ρ2(0.5)=0.5h+ρmin;ρ2(1)=h+ρmin. (17)为消除软冲击、硬冲击、激振冲击，应保证度速度ρ'(ϕ)、度加速度ρ″(ϕ)、度加速度的变化率ρ‴(ϕ)无突变．工作区的起点ϕ=0 rad和终点ϕ=1 rad均与零速区衔接，因此，起点和终点的度速度、度加速度、度加速度的变化率应同时为零，且应保证加速区和减速区衔接处三者对应相等，即：ρ1'(0)=0;ρ1″(0)=0;ρ1‴(0)=0;ρ2'(1)=0;ρ2″(1)=0;ρ2‴(1)=0;ρ1'(0.5)=ρ2'(0.5);ρ1″(0.5)=ρ2″(0.5)=0;ρ1‴(0.5)=ρ2‴(0.5)=0. (18)此外，为保证马达输出无脉动，应使处于进油区所有柱塞组件的度加速度和恒为零[13]．所以，为减小脉动，应尽可能使加速区和减速区的度速度曲线对称．因此，须增加对曲线的约束条件：ρ1″(0.25)=6h/φ12;ρ1‴(0.25)=0;ρ2″(0.75)=-6h/φ12;ρ2‴(0.75)=0. (19)式(17)~(19)共包含19个约束方程，其中9个用以约束加速区导轨曲线ρ1(ϕ)的各项系数，其余10个约束减速区导轨曲线ρ2(ϕ)．为使方程组有解，曲线的最高阶次n最小为9．其中，未知系数有2n+1-19个，加速区未知系数为A1_9，A1_10，A1_11，…，A1_n；减速区未知系数为A2_10，A2_11，A2_12…，A2_n．2.3　利用遗传算法进行优化求解当导轨曲线阶次较高时，未知系数变量较多，且目标函数(最大接触应力的最小值)求解复杂．遗传算法是模拟生物种群优胜劣汰进化法则的一种智能算法，其鲁棒性强，具有天然的并行结构，特别适合于大规模复杂问题的优化求解[16-18]，因此本研究选用遗传算法作为高阶导轨曲线系数的求解算法．以变形区最大接触应力的最小值作为适应度函数，然后利用遗传算法对2n+1-19个未知系数进行求解．在迭代过程中，须计算每一代最优解对应的理论曲线的R2，并判断是否满足式(16)的要求，以保证半径为R的滚柱能够沿着求解的曲线准确运动．3 优化结果对比分析3.1　高阶导轨曲线的特性由于遗传算法自身的特点，虽然求解效率高、求解速度快，但求解结果具有随机性，即不能保证每次求解的结果保持相同，因此本研究以最高阶次n=10和n=11为例进行优化求解，并选取其中一次的求解结果进行对比分析，以验证本研究提出的设计方法的合理性．n=10对应的导轨曲线如图5所示．n=11对应的导轨曲线如图6所示．10.13245/j.hust.211005.F005图510阶导轨曲线的度速度和度加速度曲线10.13245/j.hust.211005.F006图611阶导轨曲线的度速度和度加速度曲线等加速导轨曲线的度速度、度加速度图如参考文献[1，3]中所示．本研究设计的两条导轨曲线方程是10阶、11阶多项式，由于多项式自身的特性，因此其一阶导、二阶导和三阶导存在且连续．这可以从图5和图6中看出，相应的导轨曲线的度速度、度加速度过渡平滑、无突变．3.2　高阶导轨曲线下的最大接触应力两种曲线导轨最大接触应力如图7所示，由图可知：两条曲线的最大接触应力变化规律存在差异，这是由于不同曲线的压力角、度速度、曲率半径不同．根据式(11)，最大接触应力也不相同．此外，在加速阶段起始部位，最大接触应力逐渐增加并达到峰值，10阶导轨曲线的最大接触应力的峰值为3.84×109 Pa，而11阶导轨曲线的峰值为3.63×109 Pa．因此在柱塞组件运动过程中，11阶导轨曲线的最大接触应力峰值比10阶曲线小2.1×108 Pa．10.13245/j.hust.211005.F007图7两种曲线导轨最大接触应力对比3.3　导轨反作用力导轨曲线会影响柱塞组件的度速度、度加速度及压力角，进而影响轨道对滚柱的反作用力，而反作用力又会影响轨道与滚柱的接触应力、滚柱-柱塞摩擦副及柱塞-缸孔摩擦副的摩擦力．反作用力太大会对马达的疲劳、磨损和泄漏产生不利影响，因此选择反作用力作为轨道曲线的一项评价指标．通过式(10)计算导轨对滚柱的反作用力，并与带过渡区的等加速曲线的反作用力对比．不同导轨曲线在进油区导轨反作用力如图8所示，由图8可知：等加速导轨曲线的反作用力存在4个突变点(如图中绿色圆圈标记位置)，以第一个突变点为例，角度变化4.4×10-4 rad，反作用力却减少160 N．反作用力产生突变主要是由于等加速导轨曲线存在度加速度突变．由式(11)、(14)和(15)可知：度加速度的突变还会引起R*的突变，进一步地引起最大接触应力的突变．此外，10阶导轨曲线和11阶导轨曲线在进油区加速阶段的反作用力小于等加速导轨曲线，而在减速阶段大于等加速曲线．10.13245/j.hust.211005.F008图8不同导轨曲线在进油区导轨反作用力对比3.4　马达脉动率马达的脉动会影响执行机构运行的平稳性，当排除泄漏、摩擦、油液压缩等外部因素时，马达的脉动主要取决于其导轨曲线，因此选取马达的脉动作为导轨曲线输出特性的评价指标．马达流量脉动率[19]为δ=2(qshmax-qshmin)/(qshmax+qshmin)，式中：qshmax为最大瞬时流量；qshmin为最小瞬时流量．以单排柱塞、作用次数为8、柱塞数为14(等效柱塞数为7)的马达为例，计算本研究设计的10阶和11阶导轨曲线与文献[9]中设计的无冲击修正心形线导轨曲线的脉动率．10阶导轨曲线下马达的脉动率最小，仅为1.74%；而11阶导轨曲线下马达的脉动率为1.99%；两者均比修正心形线导轨曲线的脉动率2.54%小．4 结论a. 与常用的等加速导轨曲线相比，利用本文的设计方法得到的导轨曲线度速度、度加速度、反作用力及最大接触应力过渡平滑、无冲击．b. 与无冲击的修正心形线导轨曲线相比，利用本文的设计方法得到的导轨曲线马达脉动率较小．c. 对于本文研究的内曲线马达，优化得到的10阶导轨曲线相比于11阶，转速脉动率和反作用力更小，性能更优．d. 对于应用于军工等对冲击、噪声要求较高场合下的内曲线径向柱塞马达，本文的导轨曲线设计方法可为其优化设计提供参考．