共轴式双旋翼构型以其特有的优势被认为是发展高速直升机的一个主要方向^[1]，但近年来出现了多次共轴式直升机舰面共振事故，使共轴式直升机动不稳定问题成为迫切需要解决的问题。对于直升机旋翼/机体耦合动不稳定现象的研究，最经典的是采用旋翼多桨叶坐标变换和系统特征值分析直升机的自激振动^[2]。该方法将直升机动力学系统简化成线性系统，但桨叶减摆器和起落架的刚度和阻尼的非线性特性决定了直升机动力学系统事实上是一个旋翼/机体耦合的非线性系统。对直升机这类较复杂的非线性动力学系统现多采用数值仿真和时域参数识别，不仅能够模拟缩比模型的试验^[3]，更可以用于分析减摆器非线性^[4]、动态失速^[5]和减摆器故障^[6]等因素对直升机非线性动力学行为和动稳定性的影响。但是，现有大部分研究不管是地面共振机理分析^[7]还是设计参数影响分析^[8-10]都是针对单旋翼带尾桨直升机，而对共轴式直升机旋翼/机体耦合动稳定性的研究比较少，关于共轴式直升机旋翼/机体耦合非线性动稳定性的研究则更少。文献^[11]对共轴式直升机地面共振进行了理论和试验研究，地面共振的数学模型考虑上/下旋翼桨叶的摆振运动、上旋翼轴的弹性弯曲和机体6自由度运动，起落架用线性弹簧-线性阻尼器代替，桨叶摆振阻尼也等效成线性阻尼。文献^[12]对两片桨叶的单轴或共轴双旋翼系统在各向异性柔性支座上振动的稳定性问题，提出了一种数值仿真的分析方法。文献^[13]提出了一种共轴式直升机地面共振的时-频分析方法，通过各模态的特征值和时域响应分析，对上、下旋翼与机体之间的相互作用进行了解释。

本文建立了共轴式直升机旋翼/机体耦合的非线性动力学模型。在不考虑桨叶减摆器和起落架非线性因素的情况下，对某模型直升机进行了动态仿真，并通过对比线性化模型计算结果分析了非线性因素的影响。该非线性动力学模型也可进一步用于分析减摆器和起落架非线性因素对共轴式直升机地面共振的影响。

1 非线性动力学模型

共轴式直升机有上、下两副旋翼，各N片桨叶，以转速Ω反向旋转，如图 1所示。

由图 1可见，下旋翼俯视逆时针旋转，o₁x₁y₁z₁为下旋翼桨毂固定坐标系，跟随桨毂运动，原点o₁位于桨毂中心，o₁x₁指向机身尾部；上旋翼俯视顺时针旋转，o₂x₂y₂z₂为上旋翼桨毂固定坐标系，原点o₂位于桨毂中心，o₂x₂也指向机身尾部；摆振铰外伸量均为e，桨叶的摆振角分别为ζ₁和ζ₂；桨叶方位角分别为ψ₁和ψ₂；oxyz为固定坐标系，与机体静止时的坐标系重合。桨毂中心将随机体绕瞬心转动，绕x轴的滚转角φ_x，其瞬心位于下旋翼桨毂中心以下h_x处，绕y轴的俯仰角φ_y，其瞬心位于下旋翼桨毂中心以下h_y处，上、下旋翼之间的距离为h。

通过坐标变换的方法可以求得下旋翼和上旋翼任意剖面(距摆振铰r处)在固定坐标系中的坐标x₁，y₁和x₂，y₂，对时间进行二次求导可得加速度$\ddot{x}$₁，${\ddot{y}}$₁和$\ddot{x}$₂，${\ddot{y}}$₂，对r进行积分后可以写成如下形式

$\left\{ \begin{align} & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{x}}}_{1}}\bar{m}rdr={{A}_{1}}{{{\ddot{\zeta }}}_{1}}+{{A}_{2}}{{{\ddot{\varphi }}}_{x}}+{{A}_{3}}{{{\ddot{\varphi }}}_{y}}+{{A}_{4}} \\ & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{y}}}_{1}}\bar{m}rdr={{A}_{5}}{{{\ddot{\zeta }}}_{1}}+{{A}_{6}}{{{\ddot{\varphi }}}_{x}}+{{A}_{7}} \\ & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{x}}}_{2}}\bar{m}rdr={{B}_{1}}{{{\ddot{\zeta }}}_{2}}+{{B}_{2}}{{{\ddot{\varphi }}}_{x}}+{{B}_{3}}{{{\ddot{\varphi }}}_{y}}+{{B}_{4}} \\ & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{y}}}_{2}}\bar{m}rdr={{B}_{5}}{{{\ddot{\zeta }}}_{1}}+{{B}_{6}}{{{\ddot{\varphi }}}_{x}}+{{B}_{7}} \\ \end{align} \right.$

(1)

$\left\{ \begin{align} & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{x}}}_{1}}\bar{m}rdr={{a}_{1}}{{{\ddot{\zeta }}}_{1}}+{{a}_{2}}{{{\ddot{\varphi }}}_{x}}+{{a}_{3}}{{{\ddot{\varphi }}}_{y}}+{{a}_{4}} \\ & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{y}}}_{1}}\bar{m}rdr={{a}_{5}}{{{\ddot{\zeta }}}_{1}}+{{a}_{6}}{{{\ddot{\varphi }}}_{x}}+{{a}_{7}} \\ & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{x}}}_{2}}\bar{m}rdr={{b}_{1}}{{{\ddot{\zeta }}}_{2}}+{{b}_{2}}{{{\ddot{\varphi }}}_{x}}+{{b}_{3}}{{{\ddot{\varphi }}}_{y}}+{{b}_{4}} \\ & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{y}}}_{2}}\bar{m}rdr={{b}_{5}}{{{\ddot{\zeta }}}_{1}}+{{b}_{6}}{{{\ddot{\varphi }}}_{x}}+{{b}_{7}} \\ \end{align} \right.$

(2)

式中：A₁~A₇，B₁~B₇，a₁~a₇,b₁~b₇均不含变量的二阶导数，这样处理便于下文进行数值分析。

下旋翼和上旋翼桨叶摆振运动方程为

$\left\{ \begin{align} & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,\bar{m}{{{\ddot{x}}}_{1}}\sin {{{\bar{\varphi }}}_{1}}rdr-\overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{y}}}_{1}}\cos {{{\bar{\varphi }}}_{1}}rdr+ \\ & {{M}_{C\zeta }}\left( {{{\dot{\zeta }}}_{1}} \right)+{{M}_{K\zeta }}\left( {{\zeta }_{1}} \right)=0 \\ & \overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,\bar{m}{{{\ddot{x}}}_{2}}\sin {{{\bar{\varphi }}}_{2}}rdr-\overset{R-e}{\mathop{\underset{0}{\mathop{\int }}\,}}\,{{{\ddot{y}}}_{2}}\cos {{{\bar{\varphi }}}_{2}}rdr+ \\ & {{M}_{C\zeta }}\left( {{{\dot{\zeta }}}_{2}} \right)+{{M}_{K\zeta }}\left( {{\zeta }_{2}} \right)=0 \\ \end{align} \right.$

(3)

式中：m=m/R，m为单片桨叶质量，R为旋翼半径；ψ=ψ-ζ(桨叶后摆为正)；M_Cζ和M_Kζ分别为摆振阻尼力矩和弹性力矩。

式(1) 代入式(3) 可得下旋翼和上旋翼第k片桨叶摆振运动方程为

$\left\{ \begin{align} & ({{A}_{1k}}si{{n}_{1k}}-{{A}_{5k}}co{{s}_{1k}}){{{\ddot{\zeta }}}_{1k}}+ \\ & ({{A}_{2k}}si{{n}_{1k}}-{{A}_{6k}}co{{s}_{1k}}){{{\ddot{\varphi }}}_{x}}+ \\ & {{A}_{3k}}si{{n}_{1k}}{{{\ddot{\varphi }}}_{y}}+{{A}_{4k}}si{{n}_{1k}}-{{A}_{7k}}co{{s}_{1k}}+ \\ & {{M}_{C\zeta }}\left( {{{\dot{\zeta }}}_{1k}} \right)+{{M}_{K\zeta }}\left( {{\zeta }_{1k}} \right)=0 \\ & ({{B}_{1k}}si{{n}_{2k}}+{{B}_{5k}}co{{s}_{2k}}){{{\ddot{\zeta }}}_{2k}}+ \\ & ({{B}_{2k}}si{{n}_{2k}}+{{B}_{6k}}co{{s}_{2k}}){{{\ddot{\varphi }}}_{x}}+ \\ & {{B}_{3k}}si{{n}_{2k}}{{{\ddot{\varphi }}}_{y}}+{{B}_{4k}}si{{n}_{2k}}+{{B}_{7k}}co{{s}_{2k}}+ \\ & {{M}_{C\zeta }}\left( {{{\dot{\zeta }}}_{2k}} \right)+{{M}_{K\zeta }}\left( {{\zeta }_{2k}} \right)=0 \\ \end{align} \right.$

(4)

机体的滚转和俯仰运动方程为

$\left\{ \begin{align} & \sum\limits_{k=1}^{N}{\int\limits_{0}^{R-e}{{{h}_{x}}{{{\ddot{y}}}_{1k}}\bar{m}dr}}+\sum\limits_{k=1}^{N}{\int\limits_{0}^{R-e}{{{{\bar{h}}}_{x}}{{{\ddot{y}}}_{1k}}\bar{m}dr}}= \\ & {{I}_{fx}}{{{\ddot{\varphi }}}_{x}}+{{M}_{Cx}}\left( {{{\dot{\varphi }}}_{x}} \right)+{{M}_{Kx}}\left( {{\varphi }_{x}} \right) \\ & -(\sum\limits_{k=1}^{N}{\int\limits_{0}^{R-e}{{{h}_{x}}{{{\ddot{y}}}_{1k}}\bar{m}dr}}+\sum\limits_{k=1}^{N}{\int\limits_{0}^{R-e}{{{{\bar{h}}}_{x}}{{{\ddot{y}}}_{2k}}\bar{m}dr}}= \\ & {{I}_{fy}}{{{\ddot{\varphi }}}_{y}}+{{M}_{Cy}}\left( {{{\dot{\varphi }}}_{y}} \right)+{{M}_{Ky}}\left( {{\varphi }_{y}} \right) \\ \end{align} \right.$

(5)

式中：h_x=h+h_x，h_y=h+h_y；I_fx，I_fy分别表示机体绕x和y轴的惯性矩;M_Cx,M_Kx分别为机体滚转阻尼力矩和弹性力矩;M_Cy,M_Ky分别为机体俯仰阻尼力矩和弹性力矩。

式(2) 代入式(5) 后可得机体滚转和俯仰运动方程为

$\left\{ \begin{align} & ~-\sum\limits_{k=1}^{N}{{{a}_{5k}}}{{h}_{x}}{{{\ddot{\zeta }}}_{1k}}-\sum\limits_{k=1}^{N}{{{b}_{5k}}}{{{\bar{h}}}_{x}}{{{\ddot{\zeta }}}_{2k}}+ \\ & \begin{array}{*{35}{l}} \left[ -\sum\limits_{k=1}^{N}{\left( {{a}_{6k}}{{h}_{x}}+{{b}_{6k}}{{{\bar{h}}}_{x}} \right)}+{{I}_{fx}} \right]{{{\ddot{\varphi }}}_{x}}- \\ \sum\limits_{k=1}^{N}{\left( {{a}_{7k}}{{h}_{x}}+{{b}_{7k}}{{{\bar{h}}}_{x}} \right)}+{{M}_{Cx}}\left( {{{\bar{h}}}_{x}} \right)+{{M}_{Kx}}\left( {{{\dot{\varphi }}}_{x}} \right)=0 \\ \end{array} \\ & \sum\limits_{k=1}^{N}{{{a}_{1k}}}{{h}_{y}}{{{\ddot{\zeta }}}_{1k}}+\sum\limits_{k=1}^{N}{{{b}_{1k}}}{{{\bar{h}}}_{y}}{{{\ddot{\zeta }}}_{2k}}+\sum\limits_{k=1}^{N}{({{a}_{2k}}{{h}_{y}}+{{b}_{2k}}{{{\bar{h}}}_{y}}){{{\ddot{\varphi }}}_{x}}}+ \\ & \left[ \sum\limits_{k=1}^{N}{\left( {{a}_{3k}}{{h}_{y}}+{{b}_{3k}}{{{\bar{h}}}_{y}} \right)}+{{I}_{fy}} \right]{{{\ddot{\varphi }}}_{y}}+ \\ & \sum\limits_{k=1}^{N}{\left( {{a}_{4k}}{{h}_{y}}+{{b}_{4k}}{{{\bar{h}}}_{y}} \right)}+{{M}_{Cy}}\left( {{{\dot{\varphi }}}_{y}} \right)+{{M}_{Ky}}\left( {{\varphi }_{y}} \right)=0 \\ \end{align} \right.$

(6)

设p=[ξ₁₁,…,ξ_1N,ξ₂₁,…,ξ_2N,φ_x,φ_y]^T，则共轴式直升机旋翼/机体耦合运动方程组可写为

$M\left( p,\dot{p} \right)\ddot{p}=F\left( p,\dot{p} \right)$

(7)

令${\dot{p}}$=q，则变为一阶微分方程组

$\left[ \begin{gathered} {\dot q} \hfill \\ {\dot p} \hfill \\ \end{gathered} \right] = \left[ \begin{gathered} {M^{ - 1}}\left( {p,q} \right)F\left( {p,q} \right){\text{ }} \hfill \\ q \hfill \\ \end{gathered} \right]$

(8)

M,F均不含${\dot{q}}$，则式(8) 可以采用龙格库塔法对桨叶摆振和机体运动的时域响应进行求解。在桨叶摆振以及机体运动角度较小时，作线性化处理(sinζ=ζ，cosζ=1，sinφ=φ，cosφ=1) ，且略去二阶小量，则下旋翼和上旋翼摆振运动方程(式(4) )变为

$\left\{ \begin{align} & {{I}_{b}}{{{\ddot{\zeta }}}_{1k}}+{{S}_{b}}{{h}_{x}}cos{{\psi }_{1k}}{{{\ddot{\varphi }}}_{x}}+{{S}_{b}}{{h}_{y}}sin{{\psi }_{1k}}{{{\ddot{\varphi }}}_{y}}+ \\ & {{C}_{\zeta }}{{{\dot{\zeta }}}_{1k}}+\left( e{{\Omega }^{2}}{{S}_{b}}+{{K}_{\zeta }} \right){{\zeta }_{1k}}=0 \\ & {{I}_{b}}{{{\ddot{\zeta }}}_{2k}}-{{S}_{b}}{{{\bar{h}}}_{x}}cos{{\psi }_{2k}}{{{\ddot{\varphi }}}_{x}}+{{S}_{b}}{{{\bar{h}}}_{y}}sin{{\psi }_{2k}}{{{\ddot{\varphi }}}_{y}}+ \\ & {{C}_{\zeta }}{{{\ddot{\zeta }}}_{2k}}+\left( e{{\Omega }^{2}}{{S}_{b}}+{{K}_{\zeta }} \right){{\zeta }_{2k}}=0 \\ \end{align} \right.$

(9)

式中：S_b,I_b为桨叶的静矩和惯性矩；K_ζ,C_ζ为桨叶摆振的等效线性刚度和阻尼。机体的滚转和俯仰运动方程(式(6) )变为

$\left\{ \begin{align} & \sum\limits_{k=1}^{N}{{{S}_{b}}}{{h}_{x}}cos{{\psi }_{1k}}{{{\ddot{\zeta }}}_{1k}}-{{\sum\limits_{k=1}^{N}{{{S}_{b}}}}_{x}}{{{\bar{h}}}_{x}}cos{{\psi }_{2k}}{{{\ddot{\zeta }}}_{2k}}- \\ & \sum\limits_{k=1}^{N}{{{S}_{b}}}{{h}_{x}}\Omega (\Omega cos{{\psi }_{1k}}{{\zeta }_{1k}}+2sin{{\psi }_{1k}}{{{\dot{\zeta }}}_{1k}})+ \\ & \sum\limits_{k=1}^{N}{{{S}_{b}}}{{{\bar{h}}}_{x}}\Omega (\Omega cos{{\psi }_{2k}}{{\zeta }_{2k}}+2sin{{\psi }_{2k}}{{{\dot{\zeta }}}_{2k}})+ \\ & {{I}_{x}}{{{\ddot{\varphi }}}_{x}}+{{C}_{x}}{{{\dot{\varphi }}}_{x}}+{{K}_{x}}{{\varphi }_{x}}=0 \\ & \sum\limits_{k=1}^{N}{{{S}_{b}}}{{h}_{y}}sin{{\psi }_{1k}}{{{\ddot{\zeta }}}_{1k}}+\sum\limits_{k=1}^{N}{{{S}_{b}}}{{{\bar{h}}}_{y}}sin{{\psi }_{2k}}{{{\ddot{\zeta }}}_{2k}}- \\ & \sum\limits_{k=1}^{N}{{{S}_{b}}}{{h}_{y}}\Omega (\Omega sin{{\psi }_{1k}}{{\zeta }_{1k}}-2cos{{\psi }_{1k}}{{{\dot{\zeta }}}_{1k}})- \\ & \sum\limits_{k=1}^{N}{{{S}_{b}}}{{{\bar{h}}}_{y}}\Omega (\Omega sin{{\psi }_{2k}}{{\zeta }_{2k}}-2cos{{\psi }_{2k}}{{{\dot{\zeta }}}_{2k}})+ \\ & {{I}_{y}}{{{\ddot{\varphi }}}_{y}}+{{C}_{y}}{{{\dot{\varphi }}}_{y}}+{{K}_{y}}{{\varphi }_{y}}=0 \\ \end{align} \right.$

(10)

式中：I_x=I_fx+Nm(h_x²+h_x²)为全机绕x轴的惯性矩，I_y=I_fy+Nm(h_y²+h_y²)为全机绕y轴的惯性矩;K_x,C_x分别为机体滚转的等效线性刚度和阻尼;K_y,C_y分别为俯仰的等效线性刚度和阻尼。则共轴式直升机旋翼/机体耦合运动方程组可写为

${\rm{M\ddot p + C\dot p + Kp = 0}}$

(11)

该线性系统可采用特征值法对共轴式直升机动稳定性进行分析，若进行时域分析则可写成式(8) 的形式。

2 实例仿真

本文以铰接式旋翼为例，参数主要以NASA Ames研究中心采用的模型为主^[14]，该模型是无铰旋翼，摆振结构阻尼比为0.52%，铰接式旋翼摆振结构阻尼更小，一般设置有减摆器增大阻尼，且本文为观察不稳定区的摆振响应，人为放大摆振阻尼，延缓摆振响应发散；另外，由于模型旋翼是单旋翼，只给出转动中心至旋翼的距离，一般卡型机上下旋翼之间的距离与转动中心至下旋翼的距离相当，故本文假定二者相同。

文献^[12]对共轴式直升机旋翼/机体耦合线性化模型进行特征值分析，根据特征值计算结果绘制出模型直升机的模态阻尼曲线如图 2所示。

图 2中，LA2,LA1分别表示上、下旋翼摆振前进型模态;LR2,LR1分别表示上、下旋翼摆振后退型模态;BX,BY分别表示机体滚转和俯仰模态。由图 2可知，对线性化模型进行特征值分析确定的共振转速范围为243~327 r/min。根据时域响应是否发散，对非线性模型进行时域分析确定的共振转速范围为244~326 r/min。从共振转速范围来看，所建立的计算模型是可靠的，而且可以看出在不考虑桨叶减摆器和起落架非线性因素的前提下，线性化处理对共振转速范围影响并不大。分别选取稳定区转速233 r/min和不稳定区转速253 r/min对共轴式直升机桨叶摆振和机体运动的时域响应进行计算并分析，选取激起摆振的诱因为0.1 rad的机体初始滚转角，计算结果如图 3~6所示。

由图 3,4可见，在稳定区内，桨叶摆振和机体运动呈衰减趋势。由于桨叶摆振和机体运动是以小角度开始衰减，运动一直处于小角度范围内，且未考虑桨叶减摆器和起落架的非线性因素，此时线性化处理对桨叶摆振和机体滚转的时域响应影响并不大。

由图 5,6可见，在不稳定区内，非线性系统和线性化系统的响应特性呈现显著差异，非线性系统的桨叶摆振和机体运动发散到一定程度后会出现极限环现象。

3 结论

本文建立了共轴式直升机旋翼/机体耦合的非线性动力学模型。在不考虑桨叶减摆器和起落架非线性因素的前提下，通过求解某模型直升机不同转速时桨叶摆振和机体运动的时域响应，确定了该直升机发生“地面共振”的转速范围，并与特征值分析确定的共振转速范围进行了对比验证，证实了该动力学模型的可靠性。最后，分别选取稳定区转速233 r/min和不稳定区转速253 r/min对共轴式直升机桨叶摆振和机体运动的时域响应进行计算，从结果可知：

(1) 在稳定区内，未考虑桨叶减摆器和起落架非线性因素时，系统非线性因素的影响不大；

(2) 在不稳定区，非线性系统与线性化系统的响应特性呈现显著差异，且非线性系统将出现极限环现象。

本文所建立的共轴式直升机旋翼/机体耦合非线性动力学模型，也可用于分析减摆器和起落架非线性因素对共轴式直升机地面共振的影响。