考虑关节间隙的可调矢量喷管动力学特性

Dynamic Characteristics of Adjustable Vector Nozzles Considering Joint Clearance

摘要

针对航空发动机可调矢量喷管机构运动稳定性差的问题，开展了间隙对其动力学特性影响的研究.首先，采用第一类Lagrange方程建立可调矢量喷管单链机构动力学模型，得到单链机构驱动力变化曲线.然后采用Lankarani-Nikravesh接触力模型和修正的Coulomb摩擦模型建立运动副间隙模型，并利用函数在动力学仿真软件中构建间隙模型.最后基于虚拟样机技术建立含关节间隙的单链机构动力学模型，仿真分析了间隙位置、间隙尺寸对可调矢量喷管机构动力学特性的影响.结果表明，间隙与气动力的力臂越大，系统动力学特性波动越显著，小间隙范围内关节间隙尺寸增大时，可调矢量喷管机构动力学特性的波动幅值先减小后增大，当三角拉杆与转向控制环之间的间隙为0.3 mm时，系统动力学最稳定.

Abstract

Aiming at the problem of poor kinematic stability of the adjustable vector nozzle mechanism used in the aero‑engine system， the effect of clearance on its dynamic characteristics was studied. Firstly， the dynamic model of the adjustable vector nozzle single‑chain mechanism was established using the first type of Lagrange equation， and the driving force variation curve of the single‑chain mechanism was obtained. Then， the contact force model proposed by Lankarani and Nikravesh and the modified Coulomb friction model were used to establish the moving pair clearance model， which was constructed in the dynamic simulation software by using the function. Finally， a dynamic model of the single‑chain mechanism with joint clearance was established based on the virtual prototyping technology. The effects of clearance position and clearance size on the dynamic characteristics of the adjustable vectoring nozzle mechanism were simulated and analyzed. The results showed that the larger the clearance and aerodynamic force arm， the more significant the fluctuation of the system dynamics characteristics. As the joint clearance size increases within the small clearance range， the fluctuation amplitude of the dynamic characteristics of the adjustable vector nozzle mechanism decreases first and then increases. When the clearance between the triangular rod and the steering control ring is 0.3 mm， the system dynamics become the most stable.

Graphical abstract

关键词

航空发动机 / 可调矢量喷管 / 动力学建模 / 关节间隙 / 动力学特性

Key words

aircraft engine / adjustable vector nozzle / dynamic modeling / joint clearance / dynamic characteristics

可调矢量喷管机构工作在高温重载等复杂环境下，温度载荷和气动载荷等会造成可调矢量喷管的变形和磨损^［1］，使运动关节产生较大间隙，影响喷管运动精度.间隙关节在外载荷作用下会产生较大的接触碰撞，导致喷管机构在运动过程中出现由冲击诱发的结构振动，降低机构运动精度的同时加剧了关节磨损，降低了机构寿命^［2-3］.对含间隙机构进行动力学建模分析，一方面可以得到构件受力情况，以此进行构件强度分析、疲劳寿命预测；另一方面可以通过改变间隙值等因素探究间隙对机构动力学特性的影响规律，为可调矢量喷管机构结构优化设计和失效评估标准提供理论依据.

关节间隙的冲击载荷会导致机构运动产生振动，严重影响系统的动态响应^［4-5］.为此，部分学者在理论研究方面分析了间隙对运动机构动态特性的影响，Flores等^［6］采用Lankarani-Nikravesh（L-N）连续接触力模型研究了多间隙多体系统的动态响应，指出了含转动关节的机械系统动力学输出对间隙大小和操作条件敏感；Marques等^［7］采用一种新的混合接触力模型分析了空间旋转关节的动力学特性，通过应用实例验证了新模型的正确性和有效性.还有部分学者通过软件仿真和实验，验证了前人的研究结果，肖安建等^［8］基于正交仿真试验分析运动副间隙对塑壳断路器操作机构动力学特性的影响，指出了忽略间隙的仿真所得各指标平均偏差至少在5%以上，最大偏差高达49.45%；Bai等^［9］采用理论分析和实验验证相结合的方法分析了轴向和径向间隙对平面曲柄-滑块机构动态响应的影响，指出混合间隙机构的动态响应具有显著的抖动特征，并出现较大峰值.

虽然众多学者对含间隙系统进行了动态特性研究，但其研究对象主要是一些简单常规系统，对含有多个关节间隙的可调矢量喷管机构研究较少^［10-11］.目前针对可调矢量喷管机构的研究主要集中于矢量喷管喉道面积精确建模^［12］、运动轨迹控制^［13］等运动学问题，闫世洲等^［14］虽采用Lagrange方法推导了喷管动力学方程并得到了不同姿态下驱动力的变化规律，但并未考虑关节间隙对可调矢量喷管动力学特性的影响.因此，本文以含间隙可调矢量喷管机构为研究对象，采用L-N接触力模型和修正Coulomb摩擦模型构建间隙关节处碰撞力和摩擦力模型，并基于虚拟样机和动力学仿真技术，模拟间隙的接触碰撞，分析间隙位置、间隙值对可调矢量喷管机构动力学特性的影响规律.

1 可调矢量喷管单链动力学模型

以某型号可调矢量喷管（AVEN）机构为研究对象，其由多条图1所示的运动链并联而成，众多运动链的运动规律基本相同，因此以单运动链为对象展开AVEN机构的动力学分析. AVEN机构运动由收敛调节环和转向控制环同时控制，使喷管喉道面积满足最佳推重比^［1］，可调矢量喷管单链机构可简化为凸轮-滚子-连杆和曲柄-滑块相组合的机构，其运动简图如图2所示.

采用式（1）所示的第一类Lagrange方程对单链机构进行动力学建模.

\left[\begin{array}{cc}\mathit{M}& {\mathit{\Phi }}_{\mathit{q}}^{\mathrm{T}}\\ {\mathit{\Phi }}_{\mathit{q}}& \mathrm{0}\end{array}\right]\left[\begin{array}{l}\ddot{\mathit{q}}\\ \mathit{\lambda }\end{array}\right]=\left[\begin{array}{l}\mathit{Q}\\ \mathit{\gamma }\end{array}\right].（1）

式中： q 为广义坐标向量；\ddot{\mathit{q}}为广义坐标加速度向量； M 为各构件质量或转动惯量组成的对角矩阵； Q 为系统受到的对应于广义坐标上的作用力向量； Φ_q 为机构约束方程对广义坐标求导的雅可比矩阵；λ为Lagrange乘子； γ 为速度平方项向量.

约束方程是利用结构稳定性约束建立的构件之间的关联方程^［15］，将体坐标系建立在构件的质心处，根据构件之间的空间位置关系，运动链的约束方程如式（2）所示.

\mathit{\Phi }(\mathit{q},t)=\left[\begin{array}{c}{y}_{\mathrm{1}}-h_{\mathrm{1}}\\ x_{\mathrm{1}}+(R_{\mathrm{1}}+R_{\mathrm{2}})\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{1}}-l_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}+l_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}\\ y_{\mathrm{1}}-(R_{\mathrm{1}}+R_{\mathrm{2}})\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{1}}+l_{\mathrm{1}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}+l_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}\\ x_{\mathrm{2}}-l_{\mathrm{4}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}-l_{\mathrm{5}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}\\ y_{\mathrm{2}}+l_{\mathrm{4}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}-l_{\mathrm{5}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}\\ x_{\mathrm{3}}-l_{\mathrm{6}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{3}}-l_{\mathrm{7}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{3}}-l_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{2}}\\ y_{\mathrm{3}}+l_{\mathrm{6}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{3}}-l_{\mathrm{7}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{3}}+l_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{2}}\\ x_{\mathrm{4}}+l_{\mathrm{10}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{4}}-x_{\mathrm{3}}-l_{\mathrm{8}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{3}}-l_{\mathrm{9}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{3}}\\ y_{\mathrm{4}}-l_{\mathrm{10}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{4}}-y_{\mathrm{3}}+l_{\mathrm{8}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{3}}-l_{\mathrm{9}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{3}}\\ x_{\mathrm{5}}+l_{\mathrm{12}}-x_{\mathrm{4}}+l_{\mathrm{11}}\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{4}}\\ y_{\mathrm{5}}-y_{\mathrm{4}}-l_{\mathrm{11}}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{4}}\\ y_{\mathrm{5}}-h_{\mathrm{2}}\\ x_{\mathrm{1}}-f(t)\\ x_{\mathrm{5}}-g(t)\end{array}\right].（2）

式中： q =［x_i，y_i，{\phi }_i］，它唯一地定义了可调矢量喷管机构中所有构件的位置和方向；t为时间；h₁，h₂，R₁，R₂，l₁，l₂，…，l₁₁为结构参数；f（t），g（t）分别为收敛调节环和转向控制环的驱动方程.

2 考虑关节间隙的动力学建模

2.1 接触力模型

由于间隙的存在，轴销和轴孔之间存在3种状态^［15］：自由飞行状态、连续接触状态和碰撞穿透状态.如图3所示，在关节转动副中P_i 和P_j 分别表示相邻构件上部件i和部件j的质心，O_i 和O_j 分别表示相邻构件部件i和部件j的中心位置，定义间隙半径为

c=R_j-R_i.（3）

式中，R_i 和R_j 分别表示相邻构件部件i和部件j的半径.法向穿透深度表示为

\delta =\left|\mathit{e}\right|-c.（4）

式中， e 表示轴销和轴孔中心偏心矢量，\left|\mathit{e}\right|=\sqrt[]{{\mathit{e}}^{\mathrm{T}}\mathit{e}}.

L-N模型将滞回阻尼因子用内部阻尼引起的动能损失作为评估标准，并充分利用了间隙关节相互碰撞的构件的接触面形状特征、材料性能和运动参数等因素. Lankarani等^［16］指出，相互接触的物体接触指数在1.5附近时，L-N模型能够准确描述运动副的接触力状况，该模型的表达式为

F_{\mathrm{n}}=K{\delta }^n(\mathrm{1}+\frac{\mathrm{3}}{\mathrm{4}}\frac{\dot{\delta }}{{\dot{\delta }}_{\mathrm{0}}}(\mathrm{1}-c_{\mathrm{r}}^{\mathrm{2}})).（5）

式中：F_n表示关节法向接触力；n表示材料接触指数，金属接触时一般取1.5；\delta最大值取0.01 mm；{\dot{\delta }}_{\mathrm{0}}和c_r分别表示初始冲击速度和恢复系数；K表示广义参数，其数值表达了接触构件的材料性能和接触面形状特征，表达式为

K=\frac{\mathrm{4}}{\mathrm{3}\mathrm{\pi }({\sigma }_{\mathrm{1}}+{\sigma }_{\mathrm{2}})}{(\frac{R_i{R}_j}{R_i+R_j})}^{\frac{\mathrm{1}}{\mathrm{2}}}.（6）

{\sigma }_i=\frac{\mathrm{1}-{\nu }_i^{\mathrm{2}}}{\mathrm{\pi }{E}_i}(i=\mathrm{1,2}).（7）

式中：E_i 为弹性模量；ν_i 为泊松比.

2.2 摩擦力模型

可调矢量喷管机构运动过程中，间隙关节的轴销和轴孔接触碰撞时接触点处存在相互滑移，产生了滑动摩擦力. Coulomb摩擦定律是干接触面上最基本、最简单的摩擦模型. Coulomb摩擦定律认为摩擦系数与相对切向速度无关.但接触物体表面存在不同的摩擦状态，如滑动和黏着等，为了更加接近实际，采用修正的Coulomb摩擦模型^［17］为

F_{\mathrm{t}}=-\mu (v_{\mathrm{t}})F_{\mathrm{n}}.（8）

式中：F_t表示关节摩擦力；v_t表示相对切向速度；\mu（v_t）表示摩擦系数是相对切向速度的函数，其表达式为

\mu (v_{\mathrm{t}})=\left\{\begin{array}{cc}{\mu }_{\mathrm{s}}-\mathrm{2}{\mu }_{\mathrm{s}}{(\frac{v_{\mathrm{t}}+v_{\mathrm{s}}}{\mathrm{2}{v}_{\mathrm{s}}})}^{\mathrm{2}}(\mathrm{3}-\frac{v_{\mathrm{t}}+v_{\mathrm{s}}}{v_{\mathrm{s}}})，& ‖v_{\mathrm{t}}‖\le v_{\mathrm{s}}；\\ -\left\{{\mu }_{\mathrm{d}}+({\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}){(\frac{v_{\mathrm{t}}-v_{\mathrm{s}}}{v_{\mathrm{d}}-v_{\mathrm{s}}})}^{\mathrm{2}}[\mathrm{3}-\mathrm{2}(\frac{v_{\mathrm{t}}-v_{\mathrm{s}}}{v_{\mathrm{d}}-v_{\mathrm{s}}})]\right\}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(v_{\mathrm{t}})，& v_{\mathrm{s}}\le ‖v_{\mathrm{t}}‖\le v_{\mathrm{d}}；\\ -{\mu }_{\mathrm{d}}\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(v_{\mathrm{t}}),& ‖v_{\mathrm{t}}‖\ge v_{\mathrm{d}}.\end{array}\right.（9）

式中：{\mu }_{\mathrm{s}}和{\mu }_{\mathrm{d}}分别表示接触中的静摩擦系数和动摩擦系数；v_s和v_d分别表示最大静摩擦系数对应的相对滑动速度和稳定动摩擦系数对应的最小相对滑动速度.摩擦系数与相对滑动速度的关系曲线如图4所示.

2.3 考虑间隙的动力学模型

在对间隙关节建模时，众多学者普遍将间隙关节处的碰撞接触力视为法向接触力和切向摩擦力的合力^［18］，其表达式为

\mathit{F}=F_{\mathrm{n}}{\mathit{I}}_{\mathrm{n}}+F_{\mathrm{t}}{\mathit{I}}_{\mathrm{t}}.（10）

式中： I_n为碰撞穿透方向上的单位向量，即接触法向单位向量； I_t为接触切向单位向量，其方向为 I_n逆时针旋转90°.

当关节处存在间隙时，该位置的结构约束不再起作用，\mathit{\Phi }（ q，t）中以间隙关节为关联的方程失效，此时需要力约束替换失效的结构约束，即将式（10）所示的碰撞力 F 引入式（1）作用力 Q 中，此时系统动力学方程为

\left[\begin{array}{cc}\mathit{M}& {\mathit{\Phi }}_{\mathit{q}}^{\mathrm{T}}\\ {\mathit{\Phi }}_{\mathit{q}}& \mathrm{0}\end{array}\right]\left[\begin{array}{l}\ddot{\mathit{q}}\\ \mathit{\lambda }\end{array}\right]=\left[\begin{array}{c}\mathit{Q}+\mathit{F}\\ \mathit{\gamma }\end{array}\right].（11）

3 虚拟样机与模型仿真

3.1 含间隙单链机构虚拟样机

基于虚拟样机技术建立如图5所示的含间隙可调矢量喷管单链机动力学仿真模型，各构件关系及位置如图1所示.首先在三维建模软件中构建喷管模型，保存成工业标准格式文件导入ADAMS中，在ADAMS中进行约束、负载、接触、求解器等参数的设置.使用ADAMS提供的约束库，在构件之间添加运动副约束，运动副关系如表1所示.

利用ADAMS中的IMPACT冲击函数和STEP函数，并将2.1节中的碰撞接触力模型嵌入到ADAMS软件中，IMPACT冲击函数将接触力分解为弹性变形力和阻尼力两部分进行设置，其表达式为

（12）

（13）

式中：广义参数K由式（6）计算得到；step为STEP函数；δ_max为最大穿透深度；C_max为最大阻尼系数；D_max为接触碰撞过程中最大非线性阻尼系数，其计算公式为

D_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{\mathrm{3}}{\mathrm{4}}\frac{(\mathrm{1}-c_{\mathrm{r}}^{\mathrm{2}})}{{\dot{\delta }}_{\mathrm{0}}}\dot{\delta }.（14）

3.2 模型仿真

结合喉道面积的调节范围，式（2）中收敛调节环和转向控制环的驱动方程为

\left.\begin{array}{l}f(t)=\mathrm{115.35}-\mathrm{30}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{\pi }t)，\\ g(t)=\mathrm{46.97}-\mathrm{40}\mathrm{c}\mathrm{o}\mathrm{s}(\mathrm{2}\mathrm{\pi }t).\end{array}\right\}（15）

采用MATLAB软件对式（1）进行数值求解，得到关节5和关节11分别驱动时的驱动力随时间变化的曲线，如图6所示.将图6所示曲线用三角函数拟合后，利用STEP函数构建驱动力方程并添加为ADAMS仿真的驱动.虚拟样机中主要运动构件的材料属性如表2所示，其他构件材料定义为结构钢.仿真分析时根据表2计算式（12）和式（13）中的参数.ADAMS软件中已集成了修正的Coulomb摩擦模型，根据材料属性，设定静摩擦系数为0.49，动摩擦系数为0.46，静平移速度和摩擦平移速度采用默认值.

由于收敛调节环和转向控制环的运动相互独立，为了提高求解效率，仿真时采用单驱动方式，在初始时刻，收敛骨架与水平方向夹角为45°，扩张调节片处于水平位置，各个构件速度均为零.可调矢量喷管机构运动过程中调节片所受气动力实时变化，其变化规律尚不明确，因此在可调矢量喷管机构运动过程中认为气动力大小保持不变，并以集中力的形式添加到构件上.

4 仿真结果与分析

4.1 间隙位置对动力学特性的影响

可调矢量喷管机构连接位置均有间隙存在，为了确定间隙关节的重要性，有必要探究间隙位置对系统动力学特性的影响程度.考虑到可调矢量喷管机构工作在高温环境中，关节尺寸会因为热效应而发生改变，并根据可调矢量喷管机构实际尺寸，仿真时在模型中的主要关节位置分别设置0.5 mm间隙，间隙位置在构件1和构件2之间的关节处，如表3所示.在实际工作中，转向控制环并不会影响收敛骨架的运动，为此，不分析间隙位置1对动力学特性的影响.图7为转向控制环驱动下不同间隙位置的扩张调节片偏转角度、角速度和系统驱动力随时间变化的曲线.图8为收敛调节环驱动下不同间隙位置的收敛骨架偏转角度、角速度和系统驱动力随时间变化的曲线.

由图7a和图8a可知，间隙对扩张调节片和收敛骨架偏转角度的影响很小，即间隙不会使可调矢量喷管运动产生较大误差.在图7b和图8b中，含间隙系统的角速度出现了阶梯状波动，并且伴有脉冲特性，说明含间隙系统在运动过程中产生了冲击碰撞，且碰撞时间极短，间隙位置4导致的波动幅值相较于其他间隙位置更为显著，对系统的运动影响程度最大.在图7c和图8c中观察到驱动力的变化曲线具有相同的变化规律，其原因是气动力对间隙关节产生了额外的力矩，影响构件的偏转运动，气动力对间隙位置4的力臂最大，因此对机构运动的影响最明显.

4.2 间隙尺寸对动力学特性的影响

高温高压使得关节处的变形磨损较为显著，关节间隙也会逐渐增大，为了确定间隙的有效范围，进一步仿真分析间隙尺寸对可调矢量喷管机构动力学特性的影响规律，在对机构动力学特性影响显著的关节间隙位置（间隙位置4）设置不同的间隙尺寸.

图9给出了关节10碰撞力随时间变化的曲线.从图9中可以看出，碰撞力具有先减小后增大的变化规律，这是关节变形和间隙尺寸共同影响的结果，在间隙尺寸较小时，关节处接触变形挤压非接触区域使之发生接触，即增大了接触面积，从而导致关节接触力的切向分量增大；随着间隙尺寸的增大，有充足的空间容纳关节接触变形，额外产生的接触减少，导致接触力的切向分量减小；间隙尺寸持续增大时，轴径和孔发生碰撞的初速度增大，接触力的法向分量增大.

图10和图11给出了收敛调节环驱动下的间隙在0.2~0.5 mm变化时收敛骨架偏转角速度和关节3接触力随时间变化的曲线.由图10和图11可知，随着间隙值的增加，收敛骨架偏转的角速度和关节3接触力的振荡幅值均有先减小再增大的变化趋势，与间隙关节碰撞力变化趋势一致，说明系统的动力学特性的振荡是由间隙关节碰撞力的振荡引起的；波动频率随着间隙值的增加逐渐减小，说明间隙尺寸越大，相同时间内接触碰撞次数越少.

5 结 论

1） 间隙关节的轴和孔在运动过程中存在高频碰撞，使机构动力学特性产生波动，导致机构运动稳定性劣化.

2） 小间隙不会导致机构运动位置产生较大偏差，但会使机构在运动过程中产生冲击振动，与气动力作用方向的力臂越大的间隙位置对系统动力学特性的影响越显著.

3） 随着间隙尺寸的增大，收敛骨架的偏转角速度和滚子-骨架接触力的波动幅值先减小后增大，影响间隙关节接触力的主导因素逐渐由接触变形转变为碰撞速度，且间隙越大，碰撞速度所占的权重越高.

4） 研究发现三角拉杆与转向控制环之间的间隙对系统动力学特性的影响最为明显，且当间隙为0.3 mm时，其波动幅值最小，系统动力学最稳定.本文为改善可调矢量喷管机构的运动稳定性提供理论基础和技术支撑.

参考文献

[1] 王玉新.喷气发动机轴对称推力矢量喷管［M］.北京：国防工业出版社，2006：122-130.

[2] Wang Yu‑xin.Axial‑symmetric vectoring exhaust nozzle for jet‑thrust‑aircraft［M］.Beijing：National Defense Industry Press，2006：122-130.

[3] Erkaya S.Prediction of vibration characteristics of a planar mechanism having imperfect joints using neural network［J］.Journal of Mechanical Science and Technology，2012，26（5）：1419-1430.

[4] Ordiz M， Cuadrado J， Cabello M，et al.Prediction of fatigue life in multibody systems considering the increase of dynamic loads due to wear in clearances［J］.Mechanism and Machine Theory，2021，160：104293.

[5] Erkaya S， Uzmay I.Experimental investigation of joint clearance effects on the dynamics of a slider‑crank mechanism［J］.Multibody System Dynamics，2010，24（1）：81-102.

[6] Erkaya S.Clearance‑induced vibration responses of mechanical systems： computational and experimental investigations［J］.Journal of the Brazilian Society of Mechanical Sciences and Engineering，2018，40（2）：1-12.

[7] Flores P， Lankarani H M.Dynamic response of multibody systems with multiple clearance joints［J］.Journal of Computational and Nonlinear Dynamics，2012，7（3）：031003.

[8] Marques F， Isaac F， Dourado N，et al.An enhanced formulation to model spatial revolute joints with radial and axial clearances［J］.Mechanism and Machine Theory，2017，116：123-144.

[9] 肖安建，李红军，戴水东，等.运动副间隙对MCCB操作机构动力学特性的影响研究［J］.机械科学与技术，2024，43（10）：1654-1663.

[10] Xiao An‑jian， Li Hong‑jun， Dai Shui‑dong，et al.Study on the influence of kinematic pair clearance on the dynamic characteristics of MCCB manipulator［J］.Mechanical Science and Technology for Aerospace Engineering，2024，43（10）：1654-1663.

[11] Bai Z F， Liu T X， Li J Y，et al.Numerical and experimental study on dynamic characteristics of planar mechanism with mixed clearances［J］.Mechanics Based Design of Structures and Machines，2023，51（11）：6142-6165.

[12] Yao S B， Luo Z， Wei K，et al.Analysis of fluid‑solid‑thermal coupling characteristics of axial‑symmetric vectoring exhaust nozzle［J］.Proceedings of the Institution of Mechanical Engineers，Part C：Journal of Mechanical Engineering Science，2022，236（17）：9472-9484.

[13] Wang X Y， Wang H F， Zhao J C，et al.Rigid‑flexible coupling dynamics modeling of spatial crank‑slider mechanism based on absolute node coordinate formulation［J］.Mathematics，2022，10（6）：1-13.

[14] 刘洋，赵志刚，李维维，等.轴对称推力矢量喷管喉道面积的精度分析与补偿［J］.机械科学与技术，2018，37（7）：1137-1142.

[15] Liu Yang， Zhao Zhi‑gang， Li Wei‑wei，et al.Precision compensation of throat area for axial‑symmetric vectoring exhaust nozzle［J］.Mechanical Science and Technology for Aerospace Engineering，2018，37（7）：1137-1142.

[16] 柳亚冰，符大伟，蔡常鹏，等.轴对称矢量喷管空间运动学建模仿真［J］.航空发动机，2020，46（6）：34-40.

[17] Liu Ya‑bing， Fu Da‑wei， Cai Chang‑peng，et al.Modeling and simulation of spatial kinematics of axisymmetric vectoring nozzle［J］.Aeroengine，2020，46（6）：34-40.

[18] 闫世洲，赵志刚，霍树林.轴对称矢量喷管喷口Lagrange动力学分析［J］.机械科学与技术，2019，38（5）：809-815.

[19] Yan Shi‑zhou， Zhao Zhi‑gang， Huo Shu‑lin.Lagrange dynamics analysis of axisymmetric vector nozzle［J］.Mechanical Science and Technology for Aerospace Engineering，2019，38（5）：809-815.

[20] Marques F， Roupa I， Silva M T，et al.Examination and comparison of different methods to model closed loop kinematic chains using Lagrangian formulation with cut joint，clearance joint constraint and elastic joint approaches［J］.Mechanism and Machine Theory，2021，160：104294.

[21] Lankarani H M， Nikravesh P E.A contact force model with hysteresis damping for impact analysis of multibody systems［J］.Journal of Mechanical Design，1990，112（3）：369-376.

[22] Marques F， Flores P， Claro J C P，et al.A survey and comparison of several friction force models for dynamic analysis of multibody mechanical systems［J］.Nonlinear Dynamics，2016，86（3）：1407-1443.

[23] Skrinjar L， Slavič J， Boltežar M.A review of continuous contact‑force models in multibody dynamic［J］.International Journal of Mechanical Sciences，2018，145：171-187.

基金资助

国家科技重大专项(J2019-IV-0002-0069)

国家自然科学基金资助项目(12272089)

国家自然科学基金资助项目(U1908217)

广东省基础与应用基础研究联合基金资助项目(2020B1515120015)

中央高校基本科研业务费专项资金资助项目(N2224001-4)

中央高校基本科研业务费专项资金资助项目(N2003013)

中央高校基本科研业务费专项资金资助项目(N2303002)