谐波减速器在军事领域、通用机械领域、机器人领域以及航空航天领域中都具有广泛的运用，谐波减速器设计追求小体积、高寿命与高刚度，分析影响柔轮应力和刚度的结构参数对谐波减速器的设计尤为重要.

目前，国内该领域的研究趋势为短筒谐波减速器，不少学者在齿形上进行了大量研究. 如王家序等^[1-2]提出三维双圆弧空间齿廓，以提高谐波传动装置的啮合性能. 吴伟国等^[3]提出有一定倾角的刚轮齿形，提高了传动刚度；柴文杰等^[4-8]研究柔轮的材料^[4]、弯曲特性^[5]、动态特性^[6-8]等. 国内外广泛运用有限元法^[9-10]，对谐波减速器进行应力分析和齿形优化设计^[11]. Tjahjowidodo等^[12]为分析谐波传动不同负载下的扭矩情况，设计测试装置. Timofeev等^[13]针对高刚度谐波减速器，提出新的计算模型. 国内外学者为了提高柔轮的寿命和承载能力，对柔轮设计作了很多工作，但针对短筒（长径比为0.3~0.6）、超薄壁厚（壁厚区间小于0.6 mm）型柔轮的结构力学特性分析较少. 特别是针对柔轮小体积、高寿命、高刚度设计中尺寸的选取依据分析不够完善，缺少尺寸间的合理搭配与补偿方法，本文针对该问题开展进一步研究.

本文以工业机器人关节谐波减速器为研究对象，基于响应面理论和中心复合设计采样法，采用参数化建模的有限元分析方法，分析柔轮结构参数对应力和刚度的影响，得出结构参数对应力和刚度的影响规律. 通过研究各参数灵敏度，提出柔轮设计中各结构参数的补偿方法，拟合出响应曲线和回归方程.

为了探究各结构参数对柔轮应力和刚度的影响规律，需要确定影响应力和刚度的主要结构参数. 如图1所示为柔轮的结构简图. 柔轮各部分的几何结构参数如下：柔轮筒壁厚度为t，光滑筒体壁厚为t₁，筒长为L，内壁直径为d，变形前柔轮中线半径为r_m. 根据柔轮疲劳强度经验计算公式和圆柱壳体理论^[14]可得，影响柔轮应力的主要几何参数如下：筒壁厚度为t，筒体长度为L，柔轮半径为r，径向变形量为w₀. 为了增强参数的普适性，在研究分析中定义长径比i=L/(2r)代替参数L，定义厚径比j=t/(2r)代替t.

由于柔轮齿圈部位在装配与啮合时，同时受到柔轮联接边界效应的影响、参与啮合的齿对数、齿间载荷分布规律以及轮齿对柔轮应力分布的影响，使得柔轮应力的计算模型十分复杂. 为了简化问题，在开展齿圈建模时作出如下假设^[15].

1）考虑到齿圈的柔性，可以在满足不可延伸的条件下将齿圈部位简化为矩形截面弹性圆环.

2）弹性齿圈在装配时仅产生周向和径向2个方向的位移.

如图2所示为柔轮弹性齿圈微元段受力模型. 图中，u、w分别为周向位移和径向位移；M₁、N₁、Q₁分别为微元段界面受到的弯矩、轴向力和切向力；F_r、F_t分别为负载对微元段产生的径向力和切向力. 列出该截面受力平衡方程：

(1) $\left. \begin{split} & \rho {r_{\rm{m}}}\ddot u{\rm{d}}\theta - \frac{{\partial {N_1}}}{{\partial \theta }}{\rm{d}}\theta - {Q_1}{\rm{d}}\theta - {F_{\rm{t}}}{r_{\rm{m}}}{\rm{d}}\theta = 0,\\ & \rho {r_{\rm{m}}}\ddot w{\rm{d}}\theta - \frac{{\partial {Q_1}}}{{\partial \theta }}{\rm{d}}\theta + {N_1}{\rm{d}}\theta - {F_{\rm{r}}}{r_{\rm{m}}}{\rm{d}}\theta = 0. \end{split} \right\}$

(2) ${Q_1} = \frac{{\partial {M_1}}}{{{r_{\text{m}}}\partial \theta }},\;\;{M_1} = \frac{{EJ}}{{1 - {{{\nu }}^2}}}\left( {\frac{{\partial u}}{{r_m^2\partial \theta }} - \frac{{{\partial ^2}w}}{{r_m^2\partial {\theta ^2}}}} \right).$

式中：ρ为柔轮齿圈的线密度，J为齿圈截面系数，E为材料弹性模量，ν为材料泊松比. 由柔轮齿圈不可延伸条件，可知

(3) $\frac{{\partial u}}{{\partial \theta }} + w = 0.$

综合式（1）~（3），可得柔轮齿圈运动微分方程：

(4) $\left. \begin{aligned} & A\frac{{{\partial ^2}u}}{{\partial {t^2}}} - B\left( {\frac{{{\partial ^6}}}{{\partial {\theta ^6}}} + 2\frac{{{\partial ^4}}}{{\partial {\theta ^4}}} + \frac{{{\partial ^2}}}{{\partial {\theta ^2}}}} \right)u = {F_{\rm{t}}} + \frac{{\partial {F_{\rm{r}}}}}{{\partial \theta }},\\ & A = \rho \left( {1 - \frac{{{\partial ^2}}}{{\partial {\theta ^2}}}} \right),\;\;\;\;B = \frac{{EJ}}{{r_{\rm m}^4\left( {1 - {v^4}} \right)}}. \end{aligned}\right\}$

式中：F_r为产生的径向力. 考虑到柔轮在装配和啮合过程中运动的复杂性，对该微分方程引入修正系数k进行修正，得到以下微分运动方程^[16]：

(5) $A\frac{{{\partial ^2}u}}{{\partial {t^2}}} - B\left( {\frac{{{\partial ^6}}}{{\partial {\theta ^6}}} + 2\frac{{{\partial ^4}}}{{\partial {\theta ^4}}} + \frac{{{\partial ^2}}}{{\partial {\theta ^2}}}} \right)u = k\left( {{F_{\rm{t}}} + \frac{{\partial {F_{\rm{r}}}}}{{\partial \theta }}} \right).$

式中：k由齿圈单位长度上的径向约束刚度和切向约束刚度决定.

目前，分析柔轮应力和刚度常采用有限元结合正交试验的方法，但正交试验取样本点进行分析需要做大量的工作，试验周期长，后期对真实曲线的拟合误差较大. 近年来，基于响应面法的多目标优化设计在机械设计领域中取得重要进展，且模型函数精度高，计算效率高的优势逐渐凸显^[17-18].

响应面法的思想是对指定区域内的一系列样本点进行设计试验，通过试验结果拟合出一个响应曲面函数，近似代替真实响应曲面函数. 假设用式（6）描述响应结构Z与系统随机参量 ${{Q}} = \left[ {{Q_1},{Q_2},\cdots,{Q_{\rm{r}}}} \right]$ 之间的关系，采用一定的抽样法获得N个样本值，对相应样本值进行计算获得系统响应的样本值 $\left( {{z_1},{z_2},\cdots,{z_{\rm{s}}}} \right)$，利用最小二乘法可以得到系统响应函数. 通过响应面方程代替有限元模型进行分析，可以极大地减少计算时间^[19].

(6) ${Z_1} = {a_0} + \sum\limits_{i = 1}^R {{a_i}{Q_i}} + \sum\limits_{i = 1}^R {{a_{ii}}Q_i^2} + \sum\limits_{i < j}^{} {\sum\limits_{}^{} {{a_{ij}}{Q_i}} } {Q_j} + \varepsilon .$

式中： ${a_0}\text{、}{a_i}\text{、}{a_{ii}}\text{、}{a_{ij}}\left( {i = 1 \sim R,j = i \sim R} \right)$ 为待定系数，ε为Z₁的误差.

常用的响应面试验设计方法为Central Composite Design（CCD）和Box-Behnken Design（BBD），CCD称为中心复合设计方法；由于CCD法能够比BBD法更好地拟合响应曲面（线），采用CCD法进行采样.

通过有限元分析软件ANSYS的优化设计模块，计算多参数的应力和刚度. 采取响应曲面的功能进行分析计算，该功能通过图表形式能够动态地显示输入与输出参数之间的关系. 采用该方法需要对柔轮结构进行参数化建模，设定i、j、r、w₀为输入参数.

采用某工业机器人关节谐波减速器为建模基础. 该型号柔轮属于短筒型柔轮，基本参数如下：传动比e=80，齿数z₁=160，模数m=0.396 mm，d=61.3 mm，i=0.52，j=0.012，w₀=0.414 mm. 柔轮的材料为30CrMnSiA，弹性模量为206 GPa，泊松比为0.3. 该谐波减速器采用标准椭圆凸轮波发生器. 采用省略柔性轴承的简化方法，凸轮椭圆轮廓设计尺寸将增加柔性轴承厚度. 波发生器材料采用45钢，弹性模量为200 GPa，泊松比为0.277.

以该型号减速器为基础，开展小体积优化设计，设定i、j、r、w₀为输入参数. i的取值为 $\left[ {0.3,0.6} \right]$，j的取值为 $\left[ {{\rm{0}}{\rm{.008,0}}{\rm{.015}}} \right]$，r的取值为 $\left[ {{\rm{15,35}}} \right]$，w₀的取值为 $\left[ {0.352,0.476} \right]$.

在Workbench下的DM中，分别对柔轮和波发生器进行参数化建模. 建模时，将柔轮齿圈部分的齿简化为当量齿圈，当量厚度 $h = \sqrt[3]{{1.67}}t$. 在DM中开展柔轮和波发生器的装配，装配前先将波发生器沿短轴所在平面切成两半，然后两半波发生器分别沿长轴朝中心轴位置移动1 mm，使得柔轮与波发生器之间不发生干涉. 装配时，两半波发生器分别沿长轴朝远离中心轴方向运动1 mm，回到原始位置，便可将柔轮撑开. 柔轮参数化模型及装配图如图3所示.

波发生器与柔轮之间的接触形式采用面-面接触、刚-柔分析模式，将凸轮轮廓面定义为刚性目标面. 设置摩擦系数为0.15，接触刚度系数为0.1. 在柔轮底部施加固定约束，给两半凸轮分别设置位移约束. 在柔轮齿圈部分与刚轮拟啮合的区域的曲面上施加转矩T=63 N·m，该啮合区角度取46°.

该模型是一个几何非线性模型，根据模型的几何结构，选用线性单元和较大的网格密度，使用扫掠法划分网格以获得质量较好的六面体网格.

在ANSYS的优化设计中，需要将输出结果进行参数化，可得设计参数与目标参数之间的关系. 选取柔轮内壁、齿圈等效应力和刚度为目标参数，将这3个参数在ANSYS中参数化.

在DE模块下选取Response Surface，开展谐波减速器的优化设计分析. 设置设计变量类型为Continuous，实验设计类型为中心组合设计方法（CCD）. 按照输入参数范围，设置参数上、下界，提交进行有限元分析.

对未进行参数变动的柔轮进行静力学分析，分别建立带有轮齿形状的模型和简化后的齿圈模型. 设定相同的有限元前处理条件，求解后得到的等效应力云图如图4所示.

由图4可知，带有轮齿的模型和齿圈等效模型柔轮的最大应力均出现在齿圈后端与筒体过渡的位置，此处是由于应力集中. 另一应力较大的位置在柔轮内壁与凸轮短轴接触的位置，此处是由于径向变形较大而产生较大的应力；选取内壁应力最大点与齿圈应力最大点两处为观测点.

2种模型应力结果显示，在有轮齿的模型中，齿圈应力最大值为432.49 MPa，内壁应力最大值为305.92 MPa；齿圈等效模型齿圈处应力最大值为404.1 MPa，内壁应力最大值为302.3 MPa. 相对误差较小，且应力分布位置及云图类似，验证了齿圈简化后的等效模型的正确性.

柔轮各结构参数对应力的灵敏度定义为应力对结构参数的变化率，可以利用式（7）对灵敏度进行数学定义. 假设影响应力的各结构参数组成为 ${{X}} = {\left[ {{x_1},{x_2}, \cdots, {x_k}} \right]^{\rm{T}}}$，柔轮应力为 $\sigma = f\left( X \right)$，则柔轮的应力灵敏度 ${ s}$ 定义为

(7) ${{s}} = \nabla f = {\left[ {\frac{{\partial f}}{{\partial {x_1}}},\frac{{\partial f}}{{\partial {x_2}}}, \cdot \cdot \cdot ,\frac{{\partial f}}{{\partial {x_k}}}} \right]^{\rm{T}}} = {\left[ {{s_1},{s_2}, \cdot \cdot \cdot ,{s_k}} \right]^{\rm{T}}}.$

式中： ${s_i}(i = 1,2, \cdot \cdot \cdot ,k)$ 表示柔轮各参数对应力的灵敏度.

采用中心复合设计方法，采样得到不同参数不同水平对应的设计点；通过响应曲面分析功能，可以分析应力对各结构参数的灵敏度. i、j、r、w₀对柔轮应力的局部灵敏度图如图5所示.

从图5可知，w₀和j的分量为正值，r和i的分量为负值，这种现象表明应力对径向变形量和壁厚正响应，应力对半径和筒长负响应. 对比得出灵敏度较强的是w₀和r. 在进行柔轮小体积设计时，可以减小径向变形量来抵消由于减小半径而增加的应力. 减小长径比会增大应力，此时可以用厚径比进行补偿.

基于柔轮小体积设计时的应力补偿，应结合刚度同时开展分析，在增加柔轮寿命的同时增大柔轮的刚度. 在响应面优化设计功能中，采用CCD进行采样，可以得到柔轮参数r、w₀耦合的试验设计点和计算结果，如表1所示. 表中，p₁为内壁应力，p₂为齿圈应力，k为扭转刚度. 通过响应面拟合运算结果的变化曲线如图6、7所示.

如图6所示为柔轮应力和扭转刚度对内径r的变化曲线. 可以看出，整体上应力随着r的增大而减小，扭转刚度随着半径的增大而增大；当r<19 mm时，齿圈应力大于内壁应力，齿圈后端存在应力集中现象，则柔轮设计半径越小，越容易在齿圈后端发生疲劳破坏. 在柔轮小体积设计时，在保证小尺寸要求的情况下，应适度增大半径以减小应力，提升承载能力.

如图7所示为柔轮应力和扭转刚度对径向变形量w₀的变化曲线. 可以看出，随着w₀的增加，柔轮径向变形程度增大，整体上柔轮应力增大；柔轮内壁应力的增大速度较齿圈应力的增大速度更快，说明径向变形量对与柔性轴承直接接触的柔轮内壁的影响更明显；柔轮扭转刚度随着径向变形量的增大无明显区别，说明径向变形量对柔轮的扭转刚度影响微小. 在柔轮进行小体积设计时，应在保证轮齿啮合状态和承载能力的条件下，尽量减小径向变形量，同时加强对柔轮润滑的研究以减轻复杂应力对柔轮内壁的影响.

用CCD对柔轮的i和j进行采样，得到对应的设计点和计算结果，如表2所示. 通过响应面拟合运算结果的变化曲线，如图8、9所示.

如图8所示为柔轮应力和扭转刚度对i的变化曲线. 可以看出，随着i的增加，柔轮筒体挠性变好，应力逐渐减小，同时扭转刚度随着长径比的增加而逐渐减小. 当i<0.42时，长径比减小柔轮应力增加不明显，说明此时长径比已不是影响柔轮应力的主要因素；当i> 0.42时，齿圈应力的变化速度较内壁的变化速度更快，说明此时长径比对齿圈后端的应力影响非常明显，且当i<0.49时，柔轮齿圈应力大于内壁应力，表明长径比减小会加强齿圈后端的应力集中，使得柔轮的啮合和接触状况恶化；在柔轮进行小体积设计时，应在保证啮合状态的情况下减小筒长，提升柔轮承载能力. 同时应特别注意减轻柔轮齿圈后端的应力集中现象.

如图9所示为柔轮应力和扭转刚度对j的变化曲线. 可以看出，随着j的增大，柔轮筒体挠性降低，使得与柔性轴承直接接触的柔轮内壁应力增大；厚径比增大使得齿圈应力有减小趋势，且当j< 0.085时，齿圈应力大于内壁应力，说明增大厚径比能够在一定程度上降低齿圈后端的应力集中现象；整体上，齿圈应力的变化范围大于内壁应力，说明厚径比对齿圈应力的影响更明显；扭转刚度随着j的增大而增大. 在开展柔轮的小体积设计时，应在满足柔轮内壁接触状况的条件下尽量增大厚径比，以减小齿圈应力集中，提升啮合性能和承载能力.

柔轮设计追求高寿命、高刚度、小体积. 在柔轮设计过程中，需要在减小半径和长度的同时尽量减小应力，提高刚度. 在ANSYS中得到有限元静力学分析的数据，将模拟试验结果输入曲面响应试验分析处理软件Design Expert，可得各参数对应力和刚度的回归曲线. 该方法基于响应面模型所选取的试验骨架点，采用最小二乘法对数据结果进行回归拟合，通过回归方程寻求最优的结构参数，解决多变量问题. 该方法能够针对多目标参数分析，寻找出最优区域.

通过对灵敏度的分析可得，对柔轮应力影响最大的参数是w₀和r，其次是i和j，故在进行柔轮小体积设计时，可以将w₀和r组合考虑，以求通过减小径向变形量来补偿由于减小柔轮半径而增加的柔轮应力. 同理，可以将i和j组合考虑.

在柔轮设计中，首先根据负载确定半径，用径向变形量对应力和扭转刚度进行补偿. 对r和w₀的模拟结果进行分析，拟合出回归曲线：

(8) $\begin{split} {p_{r{w_0}}} = & 1.46{r^2} - 84.8r + 200w_0^2 + 2.72 \times {10^3}{w_0} - \\ & 51.4r{w_0} + 1.22 \times {10^3} , \end{split}$

(9) $\begin{split} {k_{r{w_0}}} = & 200{r^2} + 688r - 6.544 \times {10^4}w_0^2 + \\ & 6.712 \times {10^4}{w_0} + 1.76 \times {10^3}r{w_0} + 3.3 \times {10^4} . \end{split}$

经拟合得到的回归曲线p_rw的决定系数R²=0.999 9，校正的决定系数 $R_{{\rm{adj}}}^2$=0.999 7. 回归曲线k_rw的决定系数R²=1.0， $R_{{\rm{adj}}}^2$=1.0. 说明应力和扭转刚度的回归曲线拟合模型较好.

对柔轮半径进行设计选择时，可以根据式（8）、（9）在减小r的同时利用w₀进行补偿，以取得最优应力和刚度.

在设计中，通过减小长径比来减小柔轮尺寸，利用厚径比对应力和扭转刚度进行补偿. 对i和j的模拟结果进行分析可得，拟合出的回归曲线方程如下：

(10) $\begin{split} {p_{ij}} = & - 1.98 \times {10^3}{i^2} + 1.9 \times {10^3}i - 1.17 \times {10^6}{j^2} + \\ & 3.7 \times {10^4}j - 765.7ij - 181.5 , \end{split}$

(11) $\begin{split} {k_{ij}} = & - 2.76 \times {10^5}{i^2} + 1.56 \times {10^5}i - 1.92 \times {10^8}{j^2} + \\ & 2.56 \times {10^7}j - 1.13 \times {10^7}ij - 569 . \end{split}$

经拟合得到的回归曲线p_ij的决定系数R²=0.942 6， $R_{{\rm{adj}}}^2$=0.847 1. 回归曲线k_ij的决定系数R²=0.999 7， $R_{{\rm{adj}}}^2$=0.999 2. 说明应力和扭转刚度的回归曲线拟合模型较好.

对柔轮筒长进行设计选择时，可以根据式（10）、（11）在减小i的同时利用j进行补偿，以取得最优应力和刚度.

（1）i、j、r、w₀ 4种因素对柔轮应力的影响程度如下：影响最明显的是径向变形量和半径，其次是长径比和厚径比. 在柔轮设计时，应重点关注柔轮半径和径向变形量的选取.

（2）在柔轮的小体积设计中，可以减小径向变形量以抵消半径减小时增加的应力，减小壁厚以抵消筒长减小时内壁产生的应力. 通过响应面法可得关于柔轮应力和刚度的回归曲线，可以用于设计时获取最优值.

（3）柔轮的小体积设计减弱了柔轮的承载能力，加剧了齿圈后端的应力集中，容易在齿圈后端发生疲劳破坏. 在保证啮合状态和承载能力的情况下可以增大厚径比，减小径向变形量以减轻齿圈后端的应力集中现象.