凸轮机构广泛应用于各类型数控机床及其功能部件中，其中包括回转工作台的分度机构、动力伺服刀塔的转位机构、铲齿车床的铲削机构等。经表面处理后，凸轮表面的高硬度和耐磨性，使其具有良好的精度保持性。用于中高档机床产品的精密凸轮机构对设计、制造水平要求较高。铲齿凸轮是其中的核心部件，如果设计不当会产生很大的运动冲击，并给铲削机构带来振动噪声，直接降低工件加工精度，因此，研究铲齿凸轮的优化设计问题很有意义。

几十年来众多学者对凸轮综合设计问题进行了深入研究。李文辉^[1]以盘型凸轮机构在平面内占有的面积为目标函数，分别分析并建立了滚子摆动从动件盘型凸轮和平底直动从动件盘型凸轮机构的优化设计数学模型，但文中并未结合设计实例加以阐述。Angeles等^[2]给出了保证传动压力角的凸轮尺寸优化设计的最大压力角约束设计法，该方法计算简便，但当凸轮基圆半径较小时，将不满足轮廓曲率半径的约束。Navarro等^[3]提出基于最大压力角约束法的移动和摆动从动件凸轮机构的尺寸优化设计方法，并将其应用于凸轮急回机构、凸轮替代椭圆齿轮-四杆机构的设计。Flores^[4-5]建立了以凸轮基圆半径为优化目标，以许用压力角和最小曲率半径为约束的移动滚子从动件凸轮优化设计问题的数学模型，并利用MATLAB内置函数fmincon求解，得到了良好的设计结果。Bravo和Flocker^[6]采用PSO算法研究了带有匀速部分的凸轮回程时间最小化问题，并得到了满意的结果。Flocker等^[7-8]研究了关于多休止凸轮设计中的回程时间最小化问题。Qiu等^[9]给出用B样条曲线描述凸轮规律曲线的方法，用于凸轮廓线的优化设计，并将此方法用于分度凸轮机构的残余振动控制技术中。Fissette等^[10]利用多体力学方法建立了凸轮机构的动力学模型。

凸轮机构的综合设计问题包含系统分析和优化两个层次，现有文献主要从机构运动学方向和动力学分析方向进行研究。机构运动学角度的研究成果较为丰富，理论较为深入。这类凸轮适用于低速、间歇运动的场合，如分度凸轮机构。从凸轮机构动力学角度展开的研究，主要针对高速连续运动的场合，如发动机配气机构。凸轮动力学理论研究尚不完善，与实际应用还有很大差距。目前的凸轮综合设计问题研究多数处于理论分析层次，鲜有结合某类特定应用场合的凸轮设计方法并能够对方法加以深入讨论的研究文献。本文针对铲齿加工的凸轮机构设计场景，提出H型铲齿凸轮优化设计问题的标准解决方案，并结合实例深入讨论优化设计问题。

本文提出一种多升程H型铲齿凸轮优化设计问题的标准解决方案。首先，针对传统铲齿凸轮存在的过渡点冲击问题，提出基于H型从动件运动规律的铲齿凸轮设计方案，并建立多升程H型铲齿凸轮的优化设计模型。其次，针对上述优化设计模型的特点，将多项式变异算子和粒子群优化结合，提出多项式变异粒子群优化 (PMOPSO) 方法。最后，阐述了将PMOPSO方法应用于三升程和四升程的H型铲齿凸轮优化设计问题的计算流程。

1 H型铲齿凸轮的数学模型

铲齿凸轮的基本工作原理如图 1所示。铲刀的前刀面与滚刀轴线中心平面对齐。滚刀匀速转动，铲刀在凸轮控制下实现匀速进给，凸轮转过1周完成1个齿背的加工，凸轮铲背量为K^[11-12]。

1.1 H型从动件运动规律设计

文献[13]提出Hermite型回程曲线的铲齿凸轮设计方案并结合设计实例分析，阐明Hermite型回程曲线的设计优势。如图 2所示，H型凸轮从动件的位移规律曲线由直线和三次多项式曲线构成：

$ {s_2}\left( \theta \right) = \left\{ {\begin{array}{*{20}{l}} {\frac{K}{{{\theta _1}}}\theta, }&{\theta \in [0, {\theta _0}]}\\ {{H_3}\left( \theta \right), }&{\theta \in [{\theta _0}, {\theta _1}]} \end{array}} \right. $

(1)

式中H₃(θ) 是三次Hermite插值函数，

$ \begin{array}{l} {H_3}\left( \theta \right) = {f_0}\cdot{\alpha _0}\left( \theta \right) + {f_1}\cdot{\alpha _1}\left( \theta \right) + f{\prime _0}\cdot{\beta _0}\left( \theta \right) + \\ f{\prime _1}\cdot{\beta _1}\left( \theta \right) \end{array} $

(2)

其中：${f_0} = \frac{{K{\theta _0}}}{{{\theta _1}}}, {f_1} = 0, f{\prime _0} = f{\prime _1} = \frac{K}{{{\theta _1}}}$；基函数是

$ \begin{array}{c} {\alpha _0}\left( \theta \right) = \left( {1 + 2\frac{{\theta-{\theta _0}}}{{{\theta _1}-{\theta _0}}}} \right){\left( {\frac{{\theta-{\theta _1}}}{{{\theta _0} - {\theta _1}}}} \right)^2}\\ {\alpha _1}\left( \theta \right) = \left( {1 + 2\frac{{\theta - {\theta _1}}}{{{\theta _0} - {\theta _1}}}} \right){\left( {\frac{{\theta - {\theta _0}}}{{{\theta _1} - {\theta _0}}}} \right)^2}\\ {\beta _0}\left( \theta \right) = (\theta - {\theta _0}){\left( {\frac{{\theta - {\theta _1}}}{{{\theta _0} - {\theta _1}}}} \right)^2}\\ {\beta _1}\left( \theta \right) = (\theta - {\theta _1}){\frac{{\theta - {\theta _1}}}{{{\theta _1} - {\theta _0}}}^2} \end{array} $

令参数u=(θ-θ₀)/h，u∈[0,1], h=θ₁-θ₀, 则式 (2) 为

$ \left\{ \begin{array}{l} {s_2} = {H_3}\left( u \right) = {f_0}\cdot{\alpha _k}\left( u \right) + {f_1}\cdot{\alpha _{k + 1}}\left( u \right) + \\ \;\;\;\;\;\;f{\prime _0}\cdot{\beta _k}\left( u \right) + f{\prime _1}\cdot{\beta _{k + 1}}\left( u \right)\\ \theta = \theta \left( u \right) = {\theta _0} + hu \end{array} \right. $

(3)

图 2是H型凸轮的从动件位移和速度特性曲线。H型回程曲线可以保证过渡点处的速度连续，理论上不存在刚性冲击。但H型回程曲线在过渡点处依然存在加速度突变的柔性冲击。为了保证滚刀刀齿廓形的重磨不变性，工作时从动件应做严格的匀速运动，故首位过渡点处的加速度为零。理论上只有直线能满足过渡点处加速度连续，即无柔性冲击。实际中，直线显然不合理。因此，在此种工况下，柔性冲击和刚性冲击都不存在的回程位移曲线是不存在的。H型回程曲线设计问题是一个考虑区间一阶导数边界条件的样条曲线插值问题。H型回程曲线是一条高阶多项式曲线。可通过增加区间内设计点来改善H型回程曲线凸轮的动力学性能，但这会增加设计和制造的难度。采用三次Hermite插值曲线可平衡设计性能和制造成本。

1.2 H型凸轮廓线设计

如图 3所示，设凸轮内部区域D由理论廓线η围成，η是D的负向边界曲线。设凸轮逆时针旋转 (ω>0), θ∈[0, 2π], 偏置量是e。P点是从动件与凸轮的实际接触点，P′点是凸轮轮廓η上与P对应的点。当凸轮绕回转中心O点旋转θ(θ>0) 角后，P点坐标为 (s₀+s, -e)。理论廓线η的参数方程是：

$ \left\{ \begin{array}{l} x = ({s_0} + s){\rm{cos}}\theta-e{\rm{sin}}\theta \\ y =-({s_0} + s){\rm{sin}}\theta-e{\rm{cos}}\theta \end{array} \right. $

(4)

式中：${s_0} = \sqrt {{r^2}_0- {e^2}}, \theta \in [0, {\theta _1}]$。

工作廓线是η_K，P′点的切矢是τ，且τ与η的方向一致，则：

$ \boldsymbol{\tau} = {({\tau _x}, {\tau _y})^{\rm{T}}} = \frac{{\rm{d}}}{{{\rm{d}}\theta }}{\left( {x, y} \right)^{\rm{T}}} $

其中θ∈[0, θ₁]。令切矢τ的方向角为ψ_τ，有：

$ {\rm{cos}}{\psi _\tau } = \frac{{{\tau _x}}}{{\left| \boldsymbol{\tau} \right|}}, {\rm{sin}}{\psi _\tau } = \frac{{{\tau _y}}}{{|\boldsymbol{\tau} |}} $

(5)

工作廓线η_K的参数方程是：

$ \left\{ \begin{array}{l} {x_{\rm{K}}} = {x_{\rm{P}}} + {R_{\tau {\rm{sin}}{\psi _\tau }}}\\ {y_{\rm{K}}} = {y_{\rm{P}}}-{R_{\tau {\rm{cos}}{\psi _\tau }}} \end{array} \right. $

(6)

式中：(x_P，y_P) 是理论廓线上的点，由式 (4) 确定；滚子半径R_τ由设计计算确定；sinψ_τ，cos ψ_τ由式 (5) 确定。η_τ亦是凸轮内部区域D的负向边界曲线。

P′点处的外法矢n₀的方向角为ϕ_no。于是有${\phi _{{n_o}}} = {\psi _\tau } + \frac{{\rm{\pi }}}{2}$。曲线η_g是η_K的外等距线，则刀具中心轨迹η_g的参数方程是：

$ \left\{ \begin{array}{l} {x_C} = {x_{\rm{K}}}-{R_{\rm{g}}}{\rm{sin}}{\psi _\tau }\\ {y_C} = {y_{\rm{K}}} + {R_{\rm{g}}}{\rm{cos}}{\psi _\tau } \end{array} \right. $

(7)

式中：(x_K，y_K) 由式 (6) 确定；R_g由机床刀具确定；sin ψ_τ，cos ψ_τ由式 (5) 确定。η_g亦是凸轮内部区域D的负向边界曲线。

铲齿凸轮面积计算公式为：

$ \begin{array}{l} A = \frac{1}{2}\int_0^{2{\rm{\pi }}} {|\rho \left( \theta \right)-{R_\tau }{|^2}{\rm{d}}\theta } \approx \\ \;\;\;\;\frac{{\Delta \theta }}{2}\sum\limits_{i = 0}^{N-1} {|\rho \left( {{\theta _i}} \right)-{R_\tau }{|^2}} \end{array} $

(8)

其中：$\rho \left( \theta \right) = \sqrt {{{\left( {{s_0} + s} \right)}^2} + {e^2}}, {\rm{d}}\theta \approx \Delta \theta = \frac{{2{\rm{\pi }}}}{N}, N$是积分区间数量。

1.3 H型凸轮的理论廓线曲率和压力角分布

令有向曲率是κ(θ)，曲率半径是R(θ)，则：

$ \begin{array}{l} \kappa \left( \theta \right) = \frac{{x'y''- y'x''}}{{{{(x{'^2} + y{'^2})}^{3/2}}}}\\ R\left( \theta \right) = \frac{1}{{\left| {\kappa \left( \theta \right)} \right|}}, \theta \in [0, {\theta _1}] \end{array} $

(9)

由于边界曲线方向为负，凸轮表面外凸处有κ(θ) < 0，内凸处有κ(θ)>0。

如图 4所示，凸轮机构压力角定义为凸轮从动件与凸轮表面的外公法线方向与从动件位移方向的夹角。压力角取值范围是$\left[{-\frac{{\rm{\pi }}}{2}, \frac{{\rm{\pi }}}{2}} \right]$，其计算公式是：

$ {\rm{tan}}\phi = \frac{{{\tau _x}{\rm{cos}}\theta-{\tau _y}{\rm{sin}}\theta }}{{{\tau _y}{\rm{cos}}\theta + {\tau _x}{\rm{sin}}\theta }} = \frac{{s'-e}}{{{s_0} + s}} $

(10)

图 4是由式 (10) 得到H型铲齿凸轮的压力角分布曲线。定义压力角ϕ∈[ϕ_m, ϕ_M]，对应的凸轮角位移θ∈[θ_max, θ_min]。最大压力角的位置θ_max=0，代入式 (10) 可得最大压力角：

$ {\rm{tan}}{\phi _M} = \frac{{s\prime ({\theta _{{\rm{max}}}})-e}}{{s({\theta _{{\rm{max}}}}) + {s_0}}} $

由于θ_min∈[θ₀, θ₁), 用迭代法可得θ_min及tan ϕ_m。

1.4 H型铲齿凸轮优化设计模型

以凸轮面积为设计目标函数，选取从动件初始位移和偏置量为设计变量，考虑机构压力角分布和理论廓线曲率范围的约束条件，并结合凸轮机构的空间布局限制，建立H型铲齿凸轮的优化设计模型，如式 (11) 所示。表 1是式 (11) 的参数列表。

$ \begin{array}{l} \mathop {{\rm{min}}}\limits_\mathit{\boldsymbol{x}} \left( \mathit{\boldsymbol{x}} \right), \\ {\rm{s}}{\rm{.t}}{\rm{.}}\\ \left\{ \begin{array}{l} {g_1}\left( \mathit{\boldsymbol{x}} \right) = {\kappa _{{\rm{max}}}}\left( \mathit{\boldsymbol{x}} \right)- \frac{1}{{{R_{\rm{r}}} + {R_{\rm{g}}}}} \le 0\\ {g_2}\left( \mathit{\boldsymbol{x}} \right) = \frac{{- 1}}{{\eta {R_{\rm{r}}}}}- {\kappa _{{\rm{min}}}}\left( \mathit{\boldsymbol{x}} \right) \le 0\\ {g_3}\left( \mathit{\boldsymbol{x}} \right) = {\phi _{{\rm{max}}}}\left( \mathit{\boldsymbol{x}} \right) - {\phi _u} \le 0\\ {g_4}\left( \mathit{\boldsymbol{x}} \right) = {\phi _{\rm{l}}} - {\phi _{{\rm{min}}}}\left( \mathit{\boldsymbol{x}} \right) \le 0\\ {x_1} \in \left[{{s_{{\rm{0l}}}}, {s_{{\rm{0u}}}}} \right]\\ {x_2} \in [{e_{\rm{l}}}, {e_{\rm{u}}}] \end{array} \right. \end{array} $

(11)

设计变量x=(x₁, x₂)^T。目标函数f(x) 是H型铲齿凸轮的面积，其由数值积分公式计算得到:

$ f\left( \boldsymbol{x} \right) = \frac{1}{2}\int_0^{2{\rm{\pi }}} {{{\left[{\sqrt {{{({x_1} + s\left( \theta \right))}^2} + x_2^2}-{R_r}} \right]}^2}{\rm{d}}\theta } $

g₁(x) 和g₂(x) 是由凸轮理论廓线曲率确定的约束函数。如图 5(a)所示，κ(x) 是曲率分布函数，由式 (9) 确定，其范围由滚子半径R_τ和刀具半径R_g确定。当κ>0时凸轮轮廓内凸，当κ < 0时凸轮轮廓外凸。该模型中还引入了曲率安全系数η。

g₃(x) 和g₄(x) 是由凸轮机构压力角确定的约束函数。如图 5(b)所示，ϕ(x) 是压力角分布计算函数，由式 (10) 确定。凸轮机构在运行过程中，移动从动件需承受工作载荷，同时需要保证较高的零件制造精度。这里应选择尽量小的机构压力角，尤其是升程压力角，这样可以显著改善凸轮机构的传力状态，亦能防止细长移动部件因压力角过大而自锁。

凸轮机构的空间布局是有一定限制的，为保证设计结果的合理性，这里需要根据设计要求确定从动件初始位移和偏置量的取值范围。式 (11) 中考虑了设计变量的取值范围，如表 1所示。

2 PMOPSO方法

多项式变异粒子群优化 (PMOPSO) 方法在粒子群演化过程中，利用多项式变异算法扰动当前粒子，并以一定概率接受劣化解。随着时间的推移，劣化解接受概率逐渐趋于零，算法也收敛于最优解。在搜索后期，多项式变异算子变异产生的新点逐渐逼近当前点，因而可提升算法的局部搜索精度。PMOPSO方法是对标准粒子群优化^[14]方法的改进。

这里采用带有惯性权重w_t和收缩因子χ₀的粒子群进化方程:

$ \begin{array}{l} {v^t}_j = {w_t}\cdot{v^t}_j + {c_1}{r_1}({p_j}-{x^t}_j) + {c_2}{r_2}({N_j}-{x^t}_j)\\ v_{_j}^{t + 1} = {\chi _0} \times {v^t}_j\\ x_{_j}^{t + 1} = {x^t}_j + v_{_j}^{t + 1} \end{array} $

其中，惯性权重采用线性递减的计算公式：

$ {w_t} = {w_{{\rm{max}}}}-({w_{{\rm{max}}}}-{w_{{\rm{min}}}}) \times \frac{t}{{{t_{{\rm{max}}}}}} $

引入含固定边界约束的多项式变异算子^[15]。当前点p和变异新点x满足x，p∈[x_l, x_u], 定义扰动因子为：

$ \delta = \frac{{x-p}}{{{x_{\rm{u}}}-{x_{\rm{l}}}}} $

于是得到扰动因子δ的计算公式是：

$ \delta = \left\{ {\begin{array}{*{20}{l}} {{{\left[{2u\left( {1-2v} \right) + 2v} \right]}^{1/(\eta + 1)}} - 1, }&{u \in \left[{0, 0.5} \right]}\\ {1 - {{\left[{2u\left( {2v-1} \right) + 2\left( {1-v} \right)} \right]}^{1/(\eta + 1)}}, }&{u \in (0.5, 1]} \end{array}} \right.{\rm{ }} $

其中，随机数u满足[0,1]区间的均匀分布并记作u~U(0, 1), v=0.5(1-β)^η+1, 常数β的计算公式是:$\beta = \frac{{{\rm{min}}(p-{x_{\rm{l}}}, {x_{\rm{u}}}-p)}}{{{x_{\rm{u}}}-{x_{\rm{l}}}}}$。η是分布常数, $\eta = \frac{t}{{{t_{{\rm{max}}}}}} \times {\eta _{{\rm{max}}}}$。

一般取η_max=30.0~50.0。于是得到变异新点的计算公式是:x=p+δ(x_u-x_l)。图 6是含有固定边界的多项式变异算子产生变异新点相对于当前点的概率密度函数，图中给出了样本量为10⁴的多项式变异算子频率抽样的结果。

对于式 (11)，通过构造罚函数来处理不等式约束函数，罚函数的构造形式是：

$ F\left( {\boldsymbol{x},\pi } \right) = f\left( \boldsymbol{x} \right) + \sum\limits_{i = 1}^n {\left\{ {\begin{array}{*{20}{l}} {{\pi _i} + {{[{g_i}\left( \boldsymbol{x} \right)]}^2},}&{{g_i}\left( \boldsymbol{x} \right) > 0}\\ {0,}&{{g_i}\left( \boldsymbol{x} \right) \le 0} \end{array}} \right.} $

(12)

铲齿凸轮优化设计模型中的不等式约束函数g_i(x) 的值较小，若采用乘积形式的罚函数则要求罚因子的值较大，可能有数值计算困难。采用求和罚函数形式，便于罚因子的选取。罚因子π_i(i=1, …, n) 根据目标函数最优值的估计值选取。

PMOPSO方法中的算法参数对凸轮优化设计问题的求解精度和计算复杂度的影响很大。在满足计算精度要求的条件下，尽量降低种群规模和最大进化次数，以降低计算复杂度。采用环形 (ring) 邻域拓扑结构，个体邻域半径约为种群规模的1/5。

3 设计实例

三升程和四升程铲齿凸轮的初始设计方案的铲背量相同 (K=3.0 mm)，从动件偏距均为零。升程数的差异使2个方案的从动件初始位移值不同，四升程方案的从动件初始位移比较大。初始设计方案并未考虑凸轮机构压力角分布和理论廓线曲率范围对目标函数的影响，因而凸轮轮廓的初始设计有进一步优化的空间。表 2是三升程和四升程铲齿凸轮的初始设计参数表。2个方案中的滚子半径和刀具半径是相同的，用以确定凸轮廓线曲率的分布范围。

3.1 三升程凸轮设计实例

本小节以三升程H型铲齿凸轮为优化设计对象。如图 7所示，三升程H型铲齿凸轮的从动件在1个回转周期内有3个上升行程。与传统的单升程铲齿凸轮不同，三升程铲齿凸轮的工作轮廓是旋转对称，其几何中心与回转中心重合。在高速转动时，三升程凸轮的心轴负荷较小。三升程H型铲齿凸轮的设计约束列于表 3，其中凸轮理论廓线曲率、压力角分布、从动件初始位移和偏置量的取值范围保证了设计的合理性。

由式 (11) 建立三升程H型铲齿凸轮的优化设计模型，并采用PMOPSO方法求解，得到的优化设计结果列于表 4。表 4中由于初始设计偏于保守，有较大的改进空间。表 5是三升程H型凸轮设计实例中的PMOPSO算法参数表，其中最优目标函数值的估计值是1 100，故选取罚因子π=1 500。

图 8是优化后铲齿凸轮的理论廓线和工作廓线。图 9是优化后铲齿凸轮的理论廓线曲率范围和压力角分布情况。图 10是优化设计前后凸轮实际工作廓线的变化情况。由表 4可知，优化设计使H型铲齿凸轮面积降低约75%，铲削机构更加紧凑。

3.2 四升程凸轮设计实例

本小节以四升程H型铲齿凸轮为优化设计对象。如图 11所示，该凸轮有4个从动件上升行程。与三升程凸轮类似，四升程凸轮的工作轮廓是旋转对称的，其几何中心与回转中心重合，在高速场合下惯性力较小。四升程凸轮在每个旋转周期中要完成4次进给运动。四升程铲齿凸轮的设计约束列于表 6，凸轮理论廓线曲率、压力角、从动件初始位移和偏置量的取值范围保证了设计的合理性。四升程凸轮面积大于三升程凸轮，在铲削机构尺寸一定的条件下，四升程方案从动件初始位移的上限略低。

表 7中初始设计解能够满足公式 (11) 中的约束条件，但其较为保守，凸轮面积存在较大的设计优化空间。由式 (11) 建立四升程H型凸轮的优化设计模型，并采用PMOPSO方法获得优化设计解。表 8是PMOPSO的算法参数列表。罚因子需要根据最优目标函数的估计值选取，这里最优目标函数的估计值约为2 800，故可选取罚因子π=3 000。

图 12是PMOPSO方法中当前种群最优解的位置及其目标函数值随时间的变化曲线图。在种群进化的初始阶段，虽然目标函数的值较大，但其下降速度较大。当种群进化时间t>40时，目标函数f(x) 的值下降速度趋缓，当前种群最优个体逐渐逼近理论最优点，并逐渐收敛于理论最优解f(x^*)。

对比三升程H型铲齿凸轮 (图 8) 和四升程H型铲齿凸轮 (图 13) 的优化设计结果。由于升程数的增加，四升程凸轮轮廓的整体尺寸略大于三升程凸轮轮廓，其理论曲率也略高于三升程凸轮。由从动件的压力角分布曲线可知，四升程凸轮的回程段压力角略大于三升程凸轮，升程段的最大压力角是一致的。图 14是优化前后的凸轮工作廓线的对比图。图 15是优化后的理论廓线曲率分布和压力角分布曲线。在满足公式 (11) 中的设计约束条件下，优化设计后，四升程铲齿凸轮的面积降低约50%，这使得铲削机构更紧凑。

4 分析与讨论

本文给出一种H型铲齿凸轮的标准设计方法，可用于解决多升程H型铲齿凸轮的优化设计问题。以H型从动件运动规律为基础，建立了H型铲齿凸轮的优化设计模型，并提出多项式变异粒子群优化 (PMOPSO) 方法，将PMOPSO方法应用于三升程和四升程H型铲齿凸轮的优化设计问题。

计算结果表明，本文提出的设计方法在满足各项设计要求的前提下，显著降低多升程H型铲齿凸轮的工作轮廓面积并使铲齿机构更加紧凑。四升程方案中的从动件初始位移和偏置量均大于三升程方案，因而四升程H型铲齿凸轮的铲削机构体积大于三升程H型铲齿凸轮。然而，四升程H型铲齿凸轮的理论廓线曲率范围和压力角分布范围均小于三升程H型铲齿凸轮。综合而言：从铲削机构尺寸角度分析，三升程方案具有优势；从机构受力情况角度分析，四升程方案较佳。