电磁球型关节全悬浮转子抑振油腔廓形优化

Optimization for Profiles of Vibration-Suppressing Oil Cavities in Electromagnetic Spherical Joint Fully Suspended Rotors

摘要

为满足极端非结构环境对关节变刚度控制性能的特殊需求，解决球型关节普遍存在的电磁驱动结构复杂、耦合磁场建模困难、磁矩奇异性等问题，基于永磁转子轴线始终转向旋转磁场轴线的旋转定轴效应，设计了一种通过旋转磁场轴线控制驱动的两自由度电磁柔性球型关节。由于永磁体转子产生瞬时交变耦合磁矩，针对滚动轴承刚性支撑永磁转子导致的关节输出末端振动的问题，提出一种全悬浮转子多楔形效应动压油膜振动抑制方法。通过建立密封腔油膜流体动态特性方程推导了电磁与转子动压油膜复合变刚度转子系统动力学方程，以转子系统稳定性最优为目标，对密封腔廓形尺寸和数量进行优化设计。仿真和试验结果表明：全悬浮转子系统方案能够有效减小输出端摆动误差，振动幅值最大可减小55.70%，为新型电磁关节抗干扰理论研究奠定了基础。

Abstract

In order to meet the specific requirements of extreme non-structural environment on variable stiffness control performance of the joints， and to address the common problems such as electromagnetic drive structure complexity， coupled magnetic field modelling difficulties， magnetic moment singularity， a two-degree-of-freedom electromagnetic flexible spherical joint， driven through the axis control of the rotating magnetic field， was proposed based on the rotational coaxial effect that the lateral coupling magnetic moment drived the axis of the permanent magnet rotor to turn to the axis of rotating magnetic fields. Since the instantaneous alternating coupling magnetic moment， which generated by the permanent magnet rotor， gave rise to a vibration problem at the output end of the joints if a rigid rolling bearing was used to support the permanent magnet rotor， this paper put forward dynamic pressure oil film vibration suppression method of a fully-suspended rotor based on multi-vee effect. The fluid dynamic characteristic equation of the oil films in sealing cavities was derived， thereafter， rotor system dynamics equation with combined variable stiffness by the electromagnetic stiffness and the stiffness of rotor dynamic pressure oil film was deduced， and the contour size and number of sealing cavities were optimized with the goal of optimal stability of the rotor systems. Simulations and experiments demonstrate that the fully suspended rotor system scheme may reduce the output end’s swing errors effectively， with a maximum reduction in vibration amplitude by 55.70%， laying a foundation for the theoretical research of anti-interference in novel electromagnetic joints.

Graphical abstract

关键词

电磁球型关节 / 旋转定轴效应 / 临界电磁频率 / 廓形优化 / 振动抑制

Key words

electromagnetic spherical joint / rotating magnetic coaxial effect / critical electromagnetic frequency / profile optimization / vibration suppression

0 引言

传动关节是确保机械臂末端执行器实现精确定位和交互性能的关键因素，直接影响机械臂的灵活性与环境适应性。宇航和核电站等环境具有非结构化、干扰性强、安全性差等特征，对传动关节的性能需求颇为严苛，不仅需要结构紧凑、功重比大，还要具有较好的柔顺性和抗干扰能力。关节可变参数繁多，实现其高度集成与高性能设计无疑是一项极具挑战性的任务。

针对极端工作条件，由于传统刚性机械关节难以实时响应外部载荷变化来减缓刚性冲击，故开发刚柔并济的变刚度柔性关节已成为国内外研究焦点。现有柔性关节包括：①液压柔性关节，利用液压缸和液压马达等液压组件驱动^［1-2］，借助油液的流体特征实现关节的柔性控制，受限于高性能液压组件的供给，液压柔性关节普遍存在尺寸及质量偏大、动态柔性不足等问题。②气动柔性关节利用气体的柔性特征来实现关节柔性效果^［3-4］。该关节虽具有结构简单、柔顺性高等优点，但受限于高压气体的非线性特性，其频响特性、刚度特性及动态稳定性还有待提升。③电驱动柔性关节采用电机驱动方式，通过精确控制电机的力矩输出来改变系统的柔性^［5-6］，因结构紧凑、动态性能优越等特点而成为柔性关节研究的主要趋势，但仅依靠力矩输出控制难以满足被动柔顺的需求，因此电驱动柔性关节通常采用串联弹性元件的方式来弥补柔顺不足^［7］。依据不同调节方式，电驱动柔性关节可划分为主动变刚度柔性关节和被动变刚度柔性关节。

相较于传统多自由度电磁驱动球型关节，本文提出的电磁变刚度球型关节基于旋转定轴效应，依靠永磁转子空间磁矩在侧摆和俯仰方向的正交解耦，可确保关节在两自由度方向上的精准运动^［8］。电磁球型关节易于控制、灵活性好、环境适应性强，通过调整输入电流便可实时改变关节输出刚度，应用前景良好，具有重要的研究价值。

研究表明，尽管旋转定轴电磁变刚度球型关节优势明显，但由于永磁体转子绕空间轴线高速旋转，并输出侧摆和俯仰两方向的瞬时空间交变磁矩，如果永磁体转子采用滚动轴承刚性支撑连接，必然产生振动，因此，解决永磁体转子振动问题成为至关重要的任务。尤其是，工作环境对电磁变刚度关节的干扰将加剧振动情形，降低关节操作性能。为此，本文以缓冲转子系统振动和提高变刚度性能为目标，提出了转子系统非接触式驱动模式，利用多楔形效应优化油腔结构承载能力^［9］，拓宽油膜的收敛区间，从而提高电磁球型关节在干扰环境中的工作稳定性。

本文通过构建多楔形油腔流体动压润滑模型，借助油腔端面油膜压力，确定转子轴向偏心最小的端面廓形；并采用小参数法得到油膜侧壁动态特性参数^［10］，结合电磁刚度影响，借助稳定性分析方法，以最大当量刚度为设计标准、临界电磁频率为限制条件，最终得到油腔侧壁廓形优化函数。实验证明，多楔形油腔支撑更具有稳定优越性，有效抑制了瞬时交变磁矩对关节稳定性的影响。实验结果不仅为提升电磁关节转子系统的动态特性、加强整体稳定性提供了有力的理论支撑，还为后续关节的改进与升级奠定了基础。

1 电磁球型关节结构

电磁球型关节结构如图1所示，包括作为定子的三轴正交组合线圈（图1a）和万向随动机构（图1b）两部分。作为定子的三轴正交排布线圈由鞍型线圈和圆形亥姆霍兹线圈组合构成，鞍型组合线圈能量密度高，结构紧凑。作为十字万向节变形结构的万向随动机构通过支撑架与固定端连接，以保证关节输出端不发生自转角度偏转；在传动环外侧两正交支撑轴上分别安装阻尼器与绝对编码器，用于关节的运动实时控制与观测；传动环内侧为转子系统，包括转子和密封油腔，通过不断改变定子生成的旋转磁场轴线来调节随动永磁转子空间姿态从而改变输出轴方位。

随动机构在工作过程中可实现万向灵活转动，但由于转子在空间姿态控制过程中会受到摩擦、负载、外界电磁干扰等未知因素的影响，当采用刚性连接支撑方式时，转子产生的振动将直接作用在输出轴末端，影响关节定位精度和稳定性，故采用动油膜悬浮转子支撑方式，以油膜的动态特性吸收抑制转子振动，实现输出轴末端定位精度和稳定性。两种支撑方式简图见图2。

为了提高油膜支撑能力与稳定性，对密封油腔廓形进行优化设计。受多楔形效应的启发^［11］，由于在流体力学中，二次曲线凸起具有受力均匀、流动阻力小，且在结构设计中具有便于优化的特点，故在侧壁上增设N个圆柱形凸起，上下端面优化廓形设置N_d个二次曲面廓形，由此在供油充分时，油腔内壁与转子表面之间形成了多个收敛的楔形空间，且每一楔形空间均产生稳定油膜动压力，使转子达到动油膜支撑下的全悬浮效果。

2 转子系统动态特性参数计算

为获得转子最优多楔形油腔支撑廓形，需要对电磁球型关节转子系统稳定性进行分析，因此，首先要得到转子系统动态特性参数。采用多楔形油腔支撑方案时，永磁转子受油腔内油膜作用力及定子线圈产生的电磁场共同支撑，故需综合考虑油-电复合作用对转子稳定性的影响。

2.1 油膜压力与油膜动态特性系数

图3中，转子与转子的密封油腔组成了一个转子系统。由于油腔是封闭结构，为完成对油腔内廓形的设计以及对油腔内油膜压力的解算，需要对油腔侧壁及上下端面构建转子系统坐标系。图3a是转子侧表面与油腔侧壁间的坐标系；图3b是图1c中端面A-A区域的局部表示，用于端面廓形的设计与油膜厚度计算，是端面的俯视情况；图3c同样是端面的坐标系，是油腔端面与转子间的侧视情况。当供油足够充分时，随着永磁转子转速与楔形空间内油液流量的增大，油腔内满足形成动压油膜的必要条件。当转速稳定时，转子会悬浮在油腔某一偏心位置，油腔内油膜厚度可以通过几何关系确定。

在图3a油腔侧壁坐标系中，以油腔圆心O₁为原点建立转子初始坐标系O₁X₀Y₀，以转子圆心O₂为原点建立转子稳定点坐标系O₂X₁Y₁。R₁和R₂分别表示转子的最大半径与油腔的内径，故最大径向间隙c=R₂-R₁；R₃表示侧壁优化廓形圆柱凸起半径，凸起圆心O₃距油腔圆心O₁的圆心距为e₁，它满足e₁=R₁+c_min+R₃。当转子随旋转磁场转动且到达动压润滑稳定阶段时，根据流体动压润滑理论，两坐标系间会产生某一偏位角φ，且转子会平衡于油腔某一偏心位置，偏心大小为e₀，故稳定点坐标系由此确定。

选择A-A区域建立上下端面坐标系，如图3b所示，图中R₄代表转子端面半径。当转子稳定转动时，根据流体动压润滑理论，最终转子会稳定在偏离油腔Z轴方向中心位置的某一位置处，如图3c所示，偏心大小为e_d，其中，L为转子侧向有效高度，H为转子最大高度。

忽略曲率半径对油膜压力的影响，对油腔侧壁及上下端面进行展开，如图4所示。

对于油腔侧壁廓形，结合图3a油腔侧壁坐标系，侧壁油膜厚度h如图4a所示，根据油腔侧壁廓形的三角函数关系，在稳定点坐标系O₂X₁Y₁中容易解得：

h=\left\{\begin{array}{l}c+e_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\theta \theta \notin \Omega \\ e_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\theta +e_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{c}}-\\ \sqrt[]{R_{\mathrm{3}}^{\mathrm{2}}-e_{\mathrm{1}}^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{\alpha }_{\mathrm{c}}}-R_{\mathrm{1}}\theta \in \Omega \end{array}\right.（1）

式中：Ω为侧壁凸起区域的集合；α_c为凸起区域动点S和中点S₀分别与转子圆心O₂连线的夹角。

在油腔上下端面中，端面由N_d个二次曲面周向阵列排布，曲面凸起高度与角度θ的关系为

h_{\mathrm{t}}(\theta )=\mathrm{0.002}\mathrm{95}\frac{N_{\mathrm{d}}^{\mathrm{2}}}{{\mathrm{\pi }}^{\mathrm{2}}}(\frac{\mathrm{2}\mathrm{\pi }}{N_{\mathrm{d}}}\theta -{\theta }^{\mathrm{2}})（2）

式中：h_t为随角度θ变化的端面廓形高度。

结合图3b和图3c中坐标系，可得到图4b中端面油膜厚度h_d：

h_{\mathrm{d}}=\frac{h_{\mathrm{z}}-H}{\mathrm{2}}\pm e_{\mathrm{d}}-h_{\mathrm{t}}(\theta )（3）

式中：h_d为上下端面油膜厚度h_s和h_x的统称；h_z为油腔上下两端面最大距离。

因此，根据质点运动方程、连续方程和N-S方程，可分别建立侧壁非定常运动雷诺方程与端面油膜的雷诺方程^［12］：

\frac{\partial }{\partial \theta }(\frac{h^{\mathrm{3}}}{\eta }\frac{\partial p}{R_{\mathrm{1}}^{\mathrm{2}}\partial \theta })+\frac{h^{\mathrm{3}}}{\eta }\frac{{\partial }^{\mathrm{2}}p}{\partial Z^{\mathrm{2}}}=\frac{\mathrm{6}{u}_{\mathrm{c}}\partial h}{R_{\mathrm{1}}\partial \theta }+

\mathrm{12}(\mathrm{\Delta }\dot{Y}\mathrm{c}\mathrm{o}\mathrm{s}\theta +\mathrm{\Delta }\dot{X}\mathrm{s}\mathrm{i}\mathrm{n}\theta )（4）

\frac{\mathrm{1}}{r^{\mathrm{2}}}\frac{\partial }{\partial \theta }(\frac{h_{\mathrm{d}}^{\mathrm{3}}}{\eta }\frac{\partial p_{\mathrm{d}}}{\partial \theta })+\frac{\partial }{\partial r}(\frac{h_{\mathrm{d}}^{\mathrm{3}}}{\eta }\frac{\partial p_{\mathrm{d}}}{\partial r})=\mathrm{6}{u}_{\mathrm{d}}\frac{\partial h_{\mathrm{d}}}{r\partial \theta }（5）

式中：η为油液动力黏度；p为侧壁油膜压力分布，p_d为上下端面油膜压力p_s和p_x的统称；h、h_d分别为侧壁及上下端面油膜厚度；u_c为油腔侧壁与转子间的油膜速度；u_d为端面油膜r轴方向速度分布；\mathrm{\Delta }\dot{X}、\mathrm{\Delta }\dot{Y}分别为O₂X₁Y₁平面内转子X、Y轴方向速度变化量。

对于油膜雷诺方程式（4）及式（5），根据动压油膜特性，借助对称分布关系简化计算，为降低计算变量复杂度，结合电磁频率f，引入量纲一参数：

\begin{array}{c}\overline{Z}=\frac{\mathrm{2}Z}{L}\overline{h}=\frac{h}{c}\overline{p}=\frac{p{c}^{\mathrm{2}}}{\mathrm{4}\mathrm{\pi }\eta f{R}_{\mathrm{1}}^{\mathrm{2}}}\\ \chi ={(\frac{\mathrm{2}{R}_{\mathrm{1}}}{L})}^{\mathrm{2}}{X}^{\text{'}}=\frac{\mathrm{\Delta }\dot{X}}{\mathrm{2}\mathrm{\pi }fc}{Y}^{\text{'}}=\frac{\mathrm{\Delta }\dot{Y}}{\mathrm{2}\mathrm{\pi }fc}\\ \overline{r}=\frac{r}{R_{\mathrm{4}}}{\overline{h}}_{\mathrm{d}}=\frac{\mathrm{2}{h}_{\mathrm{d}}}{h_{\mathrm{z}}-H}{\overline{p}}_{\mathrm{d}}=\frac{p_{\mathrm{d}}(h_z{-H)}^{\mathrm{2}}}{\mathrm{48}\mathrm{\pi }\eta f{R}_{\mathrm{4}}^{\mathrm{2}}}\end{array}

则式（4）、式（5）的量纲一形式为

\frac{\partial }{\partial \theta }({\overline{h}}^{\mathrm{3}}\frac{\partial \overline{p}}{\partial \theta })+\chi {\overline{h}}^{\mathrm{3}}\frac{\partial \overline{p}}{\partial {\overline{Z}}^{\mathrm{2}}}=

\mathrm{3}\frac{\partial \overline{h}}{\partial \theta }+\mathrm{6}(Y^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{s}\theta +X^{\text{'}}\mathrm{s}\mathrm{i}\mathrm{n}\theta )（6）

\frac{\partial }{\partial \theta }({\overline{h}}_{\mathrm{d}}^{\mathrm{3}}\frac{\partial {\overline{p}}_{\mathrm{d}}}{\partial \theta })+{\overline{r}}^{\mathrm{2}}\frac{\partial }{\partial \overline{r}}({\overline{h}}_{\mathrm{d}}^{\mathrm{3}}\frac{\partial {\overline{p}}_{\mathrm{d}}}{\partial \overline{r}})={\overline{r}}^{\mathrm{2}}\frac{\partial {\overline{h}}_{\mathrm{d}}^{}}{\partial \theta }（7）

将最接近实际情况的Reynolds边界条件作为油膜压力分布数值求解边界条件^［13］。侧壁及端面油膜边界条件可分别表示为

\left.\begin{array}{l}\overline{p}{|}_{\overline{Z}=\mathrm{1}}=\mathrm{0}\frac{\partial \overline{p}}{\partial \overline{Z}}{|}_{\overline{Z}=\mathrm{0}}=\mathrm{0}\\ \overline{p}{|}_{\theta =\mathrm{0}}=\overline{p}{|}_{\theta =\mathrm{2}\mathrm{\pi }}=\mathrm{0}\end{array}\right\}（8）

\left.\begin{array}{l}{\overline{p}}_{\mathrm{d}}{|}_{\overline{r}=\mathrm{0}}={\overline{p}}_{\mathrm{d}}{|}_{\overline{r}=\mathrm{1}}=\mathrm{0}{\overline{p}}_{\mathrm{d}}{|}_{\theta ={\theta }_{\mathrm{g}}/\mathrm{2}}={\overline{p}}_{\mathrm{d}}{|}_{\theta =(k-\mathrm{0.5}){\theta }_{\mathrm{g}}}\\ \frac{\partial {\overline{p}}_{\mathrm{d}}}{\partial \theta }{|}_{\theta ={\theta }_{\mathrm{g}}/\mathrm{2}}=\frac{\partial {\overline{p}}_{\mathrm{d}}}{\partial \theta }{|}_{\theta =(k-\mathrm{0.5}){\theta }_{\mathrm{g}}}k=\mathrm{2,3},\cdots ,N_{\mathrm{d}}\end{array}\right\}（9）

根据多楔形油腔内油液动压润滑的支撑方式，由于转子侧向振动对关节输出端转角的影响较大，故忽略转子轴向振动，假设转子轴线始终与油腔轴线平行，利用小参数法，将转子在稳定点坐标系O₂X₁Y₁中受到的微小扰动量代入侧壁油膜厚度和压力的计算公式中，从而得到油膜厚度及油膜压力随扰动的变化公式：

\left.\begin{array}{l}{h}^{\text{'}}=h+\mathrm{\Delta }h=h+\mathrm{\Delta }X\mathrm{s}\mathrm{i}\mathrm{n}\theta +\mathrm{\Delta }Y\mathrm{c}\mathrm{o}\mathrm{s}\theta \\ p^{\text{'}}=p(\theta ,Z;\mathrm{\Delta }X,\mathrm{\Delta }Y,\mathrm{\Delta }\dot{X},\mathrm{\Delta }\dot{Y})=\\ p_{\mathrm{0}}+p_X\mathrm{\Delta }X+p_Y\mathrm{\Delta }Y+p_{\dot{X}}\mathrm{\Delta }\dot{X}+p_{\dot{Y}}\mathrm{\Delta }\dot{Y}\end{array}\right\}（10）

式中：p₀为转子稳定旋转阶段侧壁油膜动压力分布；p_X 、p_Y 、p_{\dot{X}}、p_{\dot{Y}}分别为油膜压力p在平衡位置处对扰动位移和扰动速度的偏导数，即侧壁油膜扰动压力分布。

将式（10）代入式（4），经过量纲一化处理、偏导化简，得到扰动压力方程^［14］：

\frac{\partial }{\partial \theta }({\overline{h}}^{\mathrm{3}}\frac{\partial {\overline{p}}_{\zeta }}{\partial \theta })+\chi {\overline{h}}^{\mathrm{3}}\frac{{\partial }^{\mathrm{2}}{\overline{p}}_{\zeta }}{\partial {\overline{Z}}^{\mathrm{2}}}=

\left\{\begin{array}{ll}\mathrm{3}\frac{\partial \overline{h}}{\partial \theta }& \zeta =\mathrm{0}\\ -\mathrm{9}\frac{\mathrm{s}\mathrm{i}\mathrm{n}\theta }{{\overline{h}}_{\mathrm{0}}}\frac{\partial {\overline{h}}_{\mathrm{0}}}{\partial \theta }-\mathrm{3}{\overline{h}}_{\mathrm{0}}^{\mathrm{3}}\frac{\partial }{\partial \theta }(\frac{\mathrm{s}\mathrm{i}\mathrm{n}\theta }{{\overline{h}}_{\mathrm{0}}})\frac{\partial {\overline{p}}_{\mathrm{0}}}{\partial \theta }+\mathrm{3}\mathrm{c}\mathrm{o}\mathrm{s}\theta & \zeta =X\\ -\mathrm{9}\frac{\mathrm{c}\mathrm{o}\mathrm{s}\theta }{{\overline{h}}_{\mathrm{0}}}\frac{\partial {\overline{h}}_{\mathrm{0}}}{\partial \theta }-\mathrm{3}{\overline{h}}_{\mathrm{0}}^{\mathrm{3}}\frac{\partial }{\partial \theta }(\frac{\mathrm{c}\mathrm{o}\mathrm{s}\theta }{{\overline{h}}_{\mathrm{0}}})\frac{\partial {\overline{p}}_{\mathrm{0}}}{\partial \theta }-\mathrm{3}\mathrm{s}\mathrm{i}\mathrm{n}\theta & \zeta =Y\\ \mathrm{6}\mathrm{s}\mathrm{i}\mathrm{n}\theta & \zeta =\dot{X}\\ \mathrm{6}\mathrm{c}\mathrm{o}\mathrm{s}\theta & \zeta =\dot{Y}\end{array}\right.（11）

借助式（8）、式（9）的Reynolds边界条件，以永磁转子的重力与油膜压力在X₀轴、Y₀轴正交分解的平衡关系作为油膜压力收敛条件，采用有限差分法对式（6）、式（7）进行离散化求解，获得转子侧壁及上下端面油膜动压力分布p₀、p_s和p_x。同时，借助侧壁油膜动压力分布p₀求解结果对式（11）求解，以获得扰动压力分布情况。

永磁转子在扰动源的影响下，侧壁油膜作用力通常不会线性变化，但在扰动足够小时可近似为线性变化，将油膜作用力在平衡位置处展开泰勒级数，并忽略高阶项可得

\left[\begin{array}{l}{F}_X\\ F_Y\end{array}\right]=\left[\begin{array}{l}{F}_{X_{\mathrm{0}}}\\ F_{Y_{\mathrm{0}}}\end{array}\right]+\frac{\partial }{\partial X}\left[\begin{array}{l}{F}_{X_{\mathrm{0}}}\\ F_{Y_{\mathrm{0}}}\end{array}\right]\mathrm{\Delta }X+\frac{\partial }{\partial Y}\left[\begin{array}{l}{F}_{X_{\mathrm{0}}}\\ F_{Y_{\mathrm{0}}}\end{array}\right]\mathrm{\Delta }Y+

\frac{\partial }{\partial \dot{X}}\left[\begin{array}{l}{F}_{X_{\mathrm{0}}}\\ F_{Y_{\mathrm{0}}}\end{array}\right]\mathrm{\Delta }\dot{X}+\frac{\partial }{\partial \dot{Y}}\left[\begin{array}{l}{F}_{X_{\mathrm{0}}}\\ F_{Y_{\mathrm{0}}}\end{array}\right]\mathrm{\Delta }\dot{Y}（12）

式中：X₀、Y₀为转子的平衡位置。

通过积分计算，将油膜作用力在O₂X₁Y₁平面上正交分解后，可得到

\left.\begin{array}{l}{F}_X=-{\int }_{-L/\mathrm{2}}^{L/\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{p}_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\theta R_{\mathrm{1}}\mathrm{d}\theta \mathrm{d}Z\\ F_Y=-{\int }_{-L/\mathrm{2}}^{L/\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{p}_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\theta R_{\mathrm{1}}\mathrm{d}\theta \mathrm{d}Z\end{array}\right\}（13）

式中：F_X 、F_Y 为O₂X₁Y₁平面内油膜正交分解力。

根据刚度和阻尼的定义，推导出多楔形油腔侧壁油膜的8个动态特性系数为^［14］

\mathit{K}=\left[\begin{array}{cc}{k}_{XX}& k_{XY}\\ k_{YX}& k_{YY}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial F_X}{\partial X}& \frac{\partial F_X}{\partial Y}\\ \frac{\partial F_Y}{\partial X}& \frac{\partial F_Y}{\partial Y}\end{array}\right]=

-{\int }_{-L/\mathrm{2}}^{L/\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}\left[\begin{array}{l}{P}_X\\ P_Y\end{array}\right]\times \left[\begin{array}{cc}\mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right]R_{\mathrm{1}}\mathrm{d}\theta \mathrm{d}Z

\mathit{C}=\left[\begin{array}{cc}{c}_{XX}& c_{XY}\\ c_{YX}& c_{YY}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial F_X}{\partial \dot{X}}& \frac{\partial F_X}{\partial \dot{Y}}\\ \frac{\partial F_Y}{\partial \dot{X}}& \frac{\partial F_Y}{\partial \dot{Y}}\end{array}\right]=

-{\int }_{-L/\mathrm{2}}^{L/\mathrm{2}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}\left[\begin{array}{l}{P}_{\dot{X}}\\ P_{\dot{Y}}\end{array}\right]\times \left[\begin{array}{cc}\mathrm{s}\mathrm{i}\mathrm{n}\theta & \mathrm{c}\mathrm{o}\mathrm{s}\theta \end{array}\right]R_{\mathrm{1}}\mathrm{d}\theta \mathrm{d}Z

式中： K 为油膜刚度系数矩阵； C 为油膜阻尼系数矩阵。

2.2 关节电磁刚度

当转子轴线与磁场轴线不重合时，在空间磁矩的作用下便产生了电磁刚度。为了求得该刚度，关节系统受力分析如图5所示。图中，Ox₁y₁z₁为旋转磁场坐标系；Ox₂y₂z₂为关节输出坐标系，它与转子稳定点坐标系O₂X₁Y₁Z重合。根据欧拉旋转定理，关节输出坐标系绕空间某一定轴Ox₃旋转ϕ角后可得到旋转磁场坐标系，因此，令Oz₃与Oz₂重合，从而建立定轴磁矩坐标系Ox₃y₃z₃。

根据耦合磁理论^［8，15］，关节在空间三方向的平均输出磁矩可表示为

\left[\begin{array}{l}{T}_{xav}\\ T_{yav}\\ T_{zav}\end{array}\right]=\frac{m_{\mathrm{0}}{B}_{\mathrm{0}}}{\mathrm{2}}\left[\begin{array}{l}\mathrm{s}\mathrm{i}\mathrm{n}\varphi (r_x\mathrm{c}\mathrm{o}\mathrm{s}(\psi +\delta )+r_y\mathrm{s}\mathrm{i}\mathrm{n}(\psi +\delta ))\\ \mathrm{s}\mathrm{i}\mathrm{n}\varphi (r_y\mathrm{c}\mathrm{o}\mathrm{s}(\psi +\delta )-r_x\mathrm{s}\mathrm{i}\mathrm{n}(\psi +\delta ))\\ \mathrm{s}\mathrm{i}\mathrm{n}\delta (\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\varphi +r_x^{\mathrm{2}}+r_y^{\mathrm{2}}-r_x^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\varphi -r_y^{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\varphi )\end{array}\right]（14）

其中，B_{\mathrm{0}}为旋转磁场强度；m₀为永磁转子磁偶极矩幅值；δ为转差角；ϕ和ψ分别为关节运动过程中随动系统绕定轴磁矩旋转轴Ox₃与Oz₂的偏转角度，具体表示为

\left.\begin{array}{l}\varphi =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{{\mathit{n}}_{\mathit{B}}\cdot {\mathit{n}}_{\mathit{m}}}{\left|{\mathit{n}}_{\mathit{B}}\right|\left|{\mathit{n}}_{\mathit{m}}\right|}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{c}\mathrm{o}\mathrm{s}{\alpha }_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{0}})\\ \psi =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{{\mathit{n}}_{{\mathit{x}}_{\mathrm{2}}}\cdot {\mathit{n}}_{{\mathit{x}}_{\mathrm{3}}}}{\left|{\mathit{n}}_{{\mathit{x}}_{\mathrm{2}}}\right|\left|{\mathit{n}}_{{\mathit{x}}_{\mathrm{3}}}\right|}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}}}{\sqrt[]{\mathrm{1}-\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\alpha }_{\mathrm{0}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\beta }_{\mathrm{0}}}}\end{array}\right\}（15）

根据Rodrigues旋转公式，旋转轴 r 为平行于定轴磁矩旋转轴Ox₃ 的单位矢量，表示为^［8］

\mathit{r}=\left[\begin{array}{l}{r}_x\\ r_y\\ r_z\end{array}\right]=\left[\begin{array}{l}\frac{\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}{\alpha }_{\mathrm{0}}}{\sqrt[]{\mathrm{1}-\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\alpha }_{\mathrm{0}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\beta }_{\mathrm{0}}}}\\ \frac{\mathrm{s}\mathrm{i}\mathrm{n}{\beta }_{\mathrm{0}}}{\sqrt[]{\mathrm{1}-\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\alpha }_{\mathrm{0}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\beta }_{\mathrm{0}}}}\\ \mathrm{0}\end{array}\right]（16）

当关节稳定工作在某一目标位置时，永磁转子会跟随旋转磁场沿其轴线匀速转动，若此时转子受到外界因素干扰，会使关节输出末端偏离稳定工作位置，由于偏转角度ϕ很小，此时转差角对恢复磁矩的影响便可忽视，即δ=0°，联立式（14）～式（16）得输出磁矩为

{\mathit{T}}_{\mathit{a}\mathit{v}}=\left[\begin{array}{l}{T}_x\\ T_y\\ T_z\end{array}\right]=\frac{m_{\mathrm{0}}{B}_{\mathrm{0}}}{\mathrm{2}}\left[\begin{array}{l}\mathrm{s}\mathrm{i}\mathrm{n}\varphi \\ \mathrm{0}\\ \mathrm{0}\end{array}\right]

设输出端沿干扰力F方向的变形量S=l\varphi，其中，l为力F的作用点到电磁变刚度关节中心的距离。干扰力F对电磁变刚度关节产生的力矩应与旋转定轴磁矩T_x 相等，结合刚度定义，并有\mathrm{s}\mathrm{i}\mathrm{n}\varphi /\varphi \approx \mathrm{1}，推导出关节电磁刚度k的表达式为

k=\frac{T_x}{l^{\mathrm{2}}\varphi }=\frac{m_{\mathrm{0}}{B}_{\mathrm{0}}}{\mathrm{2}{l}^{\mathrm{2}}}（17）

2.3 转子系统油-电复合刚度

在扰动过程中，关节永磁转子借助油膜的黏性和惯性吸收耗散扰动能量，同时在电磁及油膜复合刚度支撑下旋转，故在稳定性分析中，结合电磁刚度和油膜刚度的复合计算，考虑油膜阻尼系数的影响，对转子稳定性进行验证。图6为转子系统等效电磁刚度及油-电复合刚度示意图。

计算复合刚度时，可将电磁刚度k等效为两正交轴方向刚度k_x 和k_y，如图6a所示。根据应变能关系^［16］，得到电磁刚度k的等效形式为

k=k_x\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\theta }_k+k_y\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{\theta }_k

由于电磁驱动的物理性质，电磁刚度方向会跟随干扰力方向而变化，且各方向大小相等，故当θ_k =0时，k_x =k；当θ_k =π/2时，k_y =k。故电磁刚度系数矩阵 K_d可表示为

{\mathit{K}}_{\mathrm{d}}=\left[\begin{array}{cc}k& \mathrm{0}\\ \mathrm{0}& k\end{array}\right]

对于油-电复合刚度，如图6b所示，在单方向干扰力F的作用下，力与位移的线性关系为^［17］

\left[\begin{array}{ccc}k& -k& \mathrm{0}\\ -k& k+k_{ij}& -k_{ij}\\ \mathrm{0}& -k_{ij}& k_{ij}\end{array}\right]\left[\begin{array}{c}{u}_{\mathrm{1}}\\ u_{\mathrm{2}}\\ u_{\mathrm{3}}\end{array}\right]=\left[\begin{array}{c}{F}_{R\mathrm{1}}\\ F-F_{I\mathrm{2}}+F_{I\mathrm{2}}\\ F_{R\mathrm{3}}\end{array}\right]

由于位移u₁、u₃为0，故根据两弹性系统的刚度复合方程，容易得到复合刚度矩阵 K'为

{\mathit{K}}^{\text{'}}=\mathit{K}+{\mathit{K}}_{\mathrm{d}}（18）

3 多楔形油腔廓形优化

在油腔廓形优化过程中，对于油腔端面廓形，仅需要考虑转子轴向偏心的影响，因此借助上下端面雷诺方程的求解，端面廓形便可以确定；而侧壁廓形对转子稳定性的影响更大，因此侧壁廓形需借助转子系统动态特性参数对转子稳定性进行分析得到。

3.1 转子系统稳定性分析

以永磁转子为研究对象，结合动态特性参数，建立转子系统的扰动运动方程^［18］：

\mathit{M}\ddot{\mathit{q}}+\mathit{C}\dot{\mathit{q}}+{\mathit{K}}^{\text{'}}\mathit{q}=\mathit{F}（19）

式中： M 为质量矩阵； C 为阻尼系数矩阵； K' 为复合刚度系数矩阵； F 为关节负载力； q 为位移矩阵，\mathit{q}={[\begin{array}{cc}x& y\end{array}]}^{\mathrm{T}}；x、y为转子在O₂X₁Y₁平面内的位移。

故转子系统的扰动运动方程可表示为

\left[\begin{array}{cc}m& \mathrm{0}\\ \mathrm{0}& m\end{array}\right]\left[\begin{array}{c}\ddot{x}\\ \ddot{y}\end{array}\right]+\left[\begin{array}{cc}{c}_{XX}& c_{XY}\\ c_{YX}& c_{YY}\end{array}\right]\left[\begin{array}{c}\dot{x}\\ \dot{y}\end{array}\right]+

\left[\begin{array}{cc}{k}_{XX}+k& k_{XY}\\ k_{YX}& k_{YY}+k\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}{F}_x\\ F_y\end{array}\right]（20）

式中：m为永磁转子质量。

设方程解有指数形式：

\left.\begin{array}{l}x=x_{\mathrm{0}}{\mathrm{e}}^{\lambda t}=x_{\mathrm{0}}{\mathrm{e}}^{(\sigma +{\omega }_{\mathrm{c}}i)t}\\ y=y_{\mathrm{0}}{\mathrm{e}}^{\lambda t}=y_{\mathrm{0}}{\mathrm{e}}^{(\sigma +{\omega }_{\mathrm{c}}i)t}\end{array}\right\}（21）

式中：λ为特征值；x₀、y₀为振幅。

在关节空载条件下，F_x 与F_y 均为0，联立式（20）、式（21），由于e ^λt 恒不为0，则方程有非零解的条件可表示为

\left|\begin{array}{cc}m{\lambda }^{\mathrm{2}}+c_{XX}\lambda +k_{XX}+k& c_{XY}\lambda +k_{XY}\\ c_{YX}\lambda +k_{YX}& m{\lambda }^{\mathrm{2}}+c_{YY}\lambda +k_{YY}+k\end{array}\right|=\mathrm{0}（22）

将其展开，得

a_{\mathrm{0}}{\lambda }^{\mathrm{4}}+a_{\mathrm{1}}{\lambda }^{\mathrm{3}}+a_{\mathrm{2}}{\lambda }^{\mathrm{2}}+a_{\mathrm{3}}\lambda +a_{\mathrm{4}}=\mathrm{0}

a_{\mathrm{0}}=m^{\mathrm{2}}a_{\mathrm{1}}=m(c_{XX}+c_{YY})

a_{\mathrm{2}}=m(k_{XX}+k_{YY}+\mathrm{2}k)+c_{XX}{c}_{YY}-c_{XY}{c}_{YX}

a_{\mathrm{3}}=c_{XX}(k_{YY}+k)+c_{YY}(k_{XX}+k)-c_{XY}{k}_{YX}-c_{YX}{k}_{XY}

a_{\mathrm{4}}=(k_{XX}+k)(k_{YY}+k)-k_{XY}{k}_{YX}

根据Routh-Hurwitz稳定性准则^［19］，得到稳定性充分必要条件：

H_{\mathrm{1}}=a_{\mathrm{1}}>\mathrm{0}H_{\mathrm{2}}=a_{\mathrm{1}}{a}_{\mathrm{2}}-a_{\mathrm{0}}{a}_{\mathrm{3}}>\mathrm{0}

H_{\mathrm{3}}=a_{\mathrm{3}}{H}_{\mathrm{2}}-a_{\mathrm{1}}^{\mathrm{2}}{a}_{\mathrm{4}}>\mathrm{0}H_{\mathrm{4}}=a_{\mathrm{4}}{H}_{\mathrm{3}}>\mathrm{0}

当系统处于失稳边界状态时，根据现代控制理论，若系统稳定，它必满足的充分必要条件为特征方程的根的实部均为负，即σ<0，此时系统受扰动的影响会随着时间推移而衰减。故当σ=0时，转子开始失稳，将λ=ω_ci代入式（22），得

\begin{array}{c}{k}_{\mathrm{e}}=m\frac{a_{\mathrm{3}}}{a_{\mathrm{1}}}=\frac{c_{XX}(k_{YY}+k)+c_{YY}(k_{XX}+k)}{c_{XX}+c_{YY}}-\\ \frac{c_{XY}{k}_{YX}+c_{YX}{k}_{XY}}{c_{XX}+c_{YY}}\\ {\gamma }_{\mathrm{c}}={\omega }_{\mathrm{c}}=\sqrt[]{\frac{k_{\mathrm{e}}}{m}}\end{array}

式中：k_e为当量刚度；γ_c为转子系统临界频率。

因此，关节转子的临界电磁频率f_c为

f_{\mathrm{c}}=\frac{n_{\mathrm{c}}}{\mathrm{60}}=\frac{{\gamma }_{\mathrm{c}}}{\mathrm{2}\mathrm{\pi }}

3.2 多楔形油腔廓形的确定

根据转子动力学理论，当量刚度k_e越大，转子系统稳定性越好^［18］，故以转子系统最大当量刚度k_e为目标，对油腔侧壁廓形参数N、R₃、e₁进行优化。在关节偏转过程中，随动转子向4个工作象限运动的力学性质应趋于一致，因此，设置油腔侧壁的优化廓形凸起数量N为4的整数倍，使支撑保持刚度对称性，便于关节姿态控制。

根据图3a油腔侧壁坐标系，当侧壁凸起数量N及c_min一定时，可推导出油腔最大可设置凸起半径R_3max满足以下三角函数关系：

\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{\pi }}{N}=\frac{R_{\mathrm{2}}^{\mathrm{2}}+(R_{\mathrm{3}\mathrm{m}\mathrm{a}\mathrm{x}}+R_{\mathrm{1}}+c_{\mathrm{m}\mathrm{i}\mathrm{n}}{)}^{\mathrm{2}}-R_{\mathrm{3}\mathrm{m}\mathrm{a}\mathrm{x}}^{\mathrm{2}}}{\mathrm{2}{R}_{\mathrm{2}}(R_{\mathrm{3}\mathrm{m}\mathrm{a}\mathrm{x}}+R_{\mathrm{1}}+c_{\mathrm{m}\mathrm{i}\mathrm{n}})}

因此，得到转子系统油腔侧壁廓形优化函数：

（23）

对于油腔上下端面，在无廓形时，转子旋转过程中不满足油膜收敛条件，无法产生稳定的油膜支撑，因此，在端面设置多组曲面廓形结构，使转子在运动过程中能够带动油液进入多楔形空间，以满足油膜稳定支撑收敛条件。由于定子产生的磁场在中心位置最均匀，故转子偏心程度越小，关节定位精度越高。因此，借助式（7）对上下端面油膜压力p_s和p_x进行求解，以最小轴向偏心e_d为目标，对油腔端面廓形数量N_d进行优化。当N_d<3时，由于廓形数量过少，收敛区间产生的油膜压力不足以保证转子轴向端面稳定支撑，容易引发转子倾斜运动，故油腔端面廓形优化函数可表示为

\left.\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}{e}_{\mathrm{d}}=g(N_{\mathrm{d}})\\ \mathrm{s}.\mathrm{t}.m=\mathrm{0.821}\mathrm{k}\mathrm{g},{\mathrm{0}}^{\circ }\le \varphi \le {\mathrm{18}}^{\circ },f=\mathrm{8}\mathrm{H}\mathrm{z}\\ R_{\mathrm{4}}=\mathrm{22.2}\mathrm{m}\mathrm{m},H=\mathrm{60}\mathrm{m}\mathrm{m},h_{\mathrm{z}}=\mathrm{66}\mathrm{m}\mathrm{m}\\ N_{\mathrm{d}}\in (\mathrm{3,4},\cdots ,n),\eta \le \mathrm{0.192}\mathrm{P}\mathrm{a}\cdot \mathrm{s}\end{array}\right\}（24）

采用粒子群优化算法对优化函数进行计算求解。对于油腔侧壁廓形，在式（23）函数变量N、R₃、e₁的取值范围内进行求解，得到满足优化条件的最大转子系统当量刚度k_e，如图7所示，相较于无优化廓形油腔，当增加凸起数量时，可明显提高转子系统当量刚度，且在N=8、R₃=10 mm、e₁=35.05 mm时，满足条件的当量刚度达到最大值。

在端面廓形优化过程中，根据式（24），优化过程如图8所示，可以明显看出，随着凸起数量N_d的增加，转子偏心逐渐增大。且当N_d=3时，转子在各个工作偏转角度的偏心程度均很小，故可以忽略转子偏心对关节稳定性的影响。

优化后，经仿真计算，关节各工作角度临界电磁频率如图9所示。由于转子旋转时，若驱动频率f恰好与转子系统临界电磁频率f_c相匹配，转子将会产生共振现象，从而导致关节操作性能降低，故优化过程中以f_c为限制条件，确保驱动频率低于转子系统临界电磁频率。由于一般情况下关节驱动永磁转子的电磁频率区间为5~8 Hz，处于亚临界电磁频率区间，故关节在工作区间内的稳定性得以保障。

4 实验验证

4.1 转子系统振动实验

通过实验验证多楔形油腔优化廓形对电磁球型关节稳定性能的提升，搭建电磁球型关节实验平台，如图10所示。

在相同实验参数条件下，分别采用滚动轴承及多楔形油腔两种转子支撑方式进行关节角度偏转实验，并对不同支撑方式下的振动数据进行采集，以验证多楔形油腔支撑方案对关节振动的抑制效果，实验支撑方式如图11所示。

由于关节最大单侧偏转角度ϕ为18°，为验证关节全域偏转角度内多楔形油腔支撑的优越性，以6°为公差跨度设置4组实验偏转角度，在各偏转角度下进行两种支撑方案的振动对比，关节偏转如图12所示。

经多次测量记录，分别以两种支撑方式的最大振动幅值作为评判参数，见表1。由此得出，多楔形油腔对关节的振动抑制效果良好，最大有效振动抑制率可达55.7%，随偏转的振幅增长率最大为19.2%，达到理想预期效果。因此，关节稳定性得到大幅提升。

各实验组随时间变化的振动数据如图13所示。可以清楚看出，多楔形油腔油膜支撑条件下，关节输出末端在各偏转角度的振动情况均优于轴承支撑，且振动幅度随时间的变化维持在±0.5°以内， 呈现出周期性摆动，可以认为此时转子在多楔形油腔中得到了稳定支撑。

4.2 临界电磁频率验证

为得到临界电磁频率，实验中以最大偏转角18°作为实验偏转角度，达到目标方位后不断提高关节工作电磁频率，记录不同频率下的振动角度误差，如图14所示。

由图14可以看出，当关节以8 Hz电磁频率工作时，关节输出末端的振动误差呈现周期性规律振动，且平均振动幅值不超过1°，转子旋转较为稳定；当关节工作电磁频率提升至9 Hz时，振动呈现明显增大现象；当频率提高至10 Hz时，出现发散的振动信号，此时认为转子处于失稳状态；而继续提高电磁频率后，发散信号消失，振动幅度有所改善。由此可确定，临界电磁频率f_c为10 Hz，当量刚度k_e为3237.89 N/m。结果显示，临界电磁频率及当量刚度实际值与模型理论计算值的误差分别为1.9%和3.7%，实验值与理论值差异较小，具有良好的代表性。同时可以进一步说明，当工作电磁频率低于临界电磁频率10 Hz时，手腕工作性能趋于稳定。

5 结论

1）电磁球型关节永磁转子在油-电复合刚度协同支撑作用下转子稳定性更好，关节抗干扰能力与可操作性得到提升，全悬浮转子支撑方案为电磁关节振动抑制问题提供了新思路。

2）在多楔形油腔廓形作用下，转子系统当量刚度大幅提升，转子悬浮定心能力显著增强，滚动轴承刚性支撑引发的电磁球型手腕输出末端振动问题得以解决。

3）通过结构设计提高关节转子系统临界电磁频率，保证关节驱动频率处于亚临界电磁频率区间，避免转子共振降低关节操作性能的可能性，且关节在亚临界电磁频率内工作稳定性表现更优秀，为关节的最优控制提供了理论支撑。

参考文献

[1] YU Bin， LI Huashun， MA Guoliang， et al. Design and Matching Control Strategy of Electro-hydraulic Load-sensitive Hydraulic Power Unit for Legged Robots［J］. Energy， 2024， 313： 133730.

[2] ZONG Huaizhi， AI Jikun， FANG Lizhou， et al. A Novel Hydraulic Swing Actuator with High Torque Density for Legged Robots［J］. Smart Materials and Structures， 2025， 34（1）： 015034.

[3] GARCÍA-SAMARTÍN J F， RIEKER A， BARRIENTOS A.Design， Manufacturing， and Open-loop Control of a Soft Pneumatic Arm［J］. Actuators， 2024， 13（1）： 36.

[4] ZHAO Jinfeng， WU Changqu， WANG Wenbiao， et al. Design and Implementation of Variable Stiffness Rigid-soft Coupling Pneumatic Actuated Joint［C］∥2021 IEEE International Conference on Real-time Computing and Robotics （RCAR）. Xining， 2021： 679-683.

[5] CHIARADIA D， RINALDI G， SOLAZZI M， et al. Design and Control of the Rehab-Exos， a Joint Torque-controlled Upper Limb Exoskeleton［J］. Robotics， 2024， 13（2）： 32.

[6] PRAUTZSCH T， BERNINGER T F C， RIXEN D J. Investigation of Torque Controlled Robots with Flexible Links Using a Flexible Multibody Simulation［C］∥2021 IEEE/ASME International Conference on Advanced Intelligent Mechatronics （AIM）. Delft， 2021： 638-644.

[7] SARKISIAN S V， GABERT L， LENZI T. Series-elastic Actuator with Two Degree-of-freedom PID Control Improves Torque Control in a Powered Knee Exoskeleton［J］. Wearable Technologies， 2023， 4： e25.

[8] 张永顺， 刘高仁， 刘志军， 等. 新型电磁球型手腕解耦驱动机理［J］. 机械工程学报， 2024， 60（19）： 1-10.

[9] ZHANG Yongshun， LIU Gaoren， LIU Zhijun， et al. Decoupling Drive Mechanism of a New Type of Electromagnetic Spherical Wrist［J］. Journal of Mechanical Engineering， 2024， 60（19）： 1-10.

[10] 张永顺， 王履华， 程存欣， 等. 多楔形效应胶囊机器人花瓣廓形优化［J］. 机械工程学报， 2015， 51（7）： 45-52.

[11] ZHANG Yongshun， WANG Lühua， CHENG Cunxin， et al. Petal Shape Optimization of a Capsule Robot with Multiple Wedge Effects［J］. Journal of Mechanical Engineering， 2015， 51（7）： 45-52.

[12] DU Jialei， LIANG Guozhu. Dynamic Coefficients and Stability Analysis of a Water-lubricated Hydrostatic Bearing by Solving the Uncoupled Reynolds Equation［J］. Chinese Journal of Aeronautics， 2020， 33（8）： 2110-2122.

[13] 于英华， 周乐， 阮文新， 等. 特形微织构多油楔滑动轴承特性分析及优化［J］. 西北工业大学学报， 2023， 41（1）： 222-229.

[14] YU Yinghua， ZHOU Le， RUAN Wenxin， et al. Characteristic Analysis and Optimization of Special-shaped Micro-textured Multi-oil Wedge Sliding Bearing［J］. Journal of Northwestern Polytechnical University， 2023， 41（1）： 222-229.

[15] KIM H， JANG G， LEE S. Complete Determination of the Dynamic Coefficients of Coupled Journal and Thrust Bearings Considering Five Degrees of Freedom for a General Rotor-bearing System［J］. Microsystem Technologies， 2011， 17（5）： 749-759.

[16] 赖勇能. 可倾瓦滑动轴承性能分析及稳定性研究［D］. 南京： 东南大学， 2020.

[17] LAI Yongneng. Investgation on Characteristics and Stability of Tilting Pad Journal Bearing［D］. Nanjing： Southeast University， 2020.

[18] MENG Fanming， ZHANG Yifei， SU Linlin， et al. Dynamic Characteristics of Compound Textured Journal Bearing［J］. Proceedings of the Institution of Mechanical Engineers， Part J： Journal of Engineering Tribology， 2021， 235（7）： 1312-1334.

[19] 张永顺， 邢立君， 董海， 等. 双半球胶囊机器人滑滚驱动机理［J］. 机械工程学报， 2023， 59（23）： 87-95.

[20] ZHANG Yongshun， XING Lijun， DONG Hai， et al. Sliding and Rolling Mechanism of a Dual Hemisphere Capsule Robot［J］. Journal of Mechanical Engineering， 2023， 59（23）： 87-95.

[21] 苟文选. 材料力学-Ⅱ［M］. 2版. 北京： 科学出版社， 2010.

[22] GOU Wenxuan. Mechanics of Materials-Ⅱ［M］. 2nd ed. Beijing： Science Press， 2010.

[23] 可帅， 冯治国， 李长虹， 等. 减振器螺旋弹簧组串并联刚度计算与分析［J］. 机械设计与制造， 2018（9）： 25-28.

[24] KE Shuai， FENG Zhiguo， LI Changhong， et al. Shock Absorber Spring Group Series-parallel Stiffness Calculation and Analysis［J］. Machinery Design & Manufacture， 2018（9）： 25-28.

[25] 李争， 陈秋蒴， 王群京. 多自由度永磁同步电机转子动力学特性分析［J］. 电工技术学报， 2019， 34（11）： 2269-2276.

[26] LI Zheng， CHEN Qiushuo， WANG Qunjing. Rotor Dynamics Characteristic Analysis of Multi-degrees-of-freedom Permanent Magnet Synchronous Motor［J］. Transactions of China Electrotechnical Society， 2019， 34（11）： 2269-2276.

[27] WANG Jianlong， ADHIKARI G， KOBAYASHI H， et al. Analysis and Design of Operational Amplifier Stability Using Routh-Hurwitz Stability Criterion［J］. Applied Mechanics and Materials， 2019， 888： 1-10.

基金资助

国家自然科学基金(62173059)