颤振是一种由气动弹性引起的复杂的、自激振动的不稳定现象，会对飞行安全构成极大的威胁，是飞行器设计和研制过程中需要重点解决的关键问题^[1]。已有研究表明：尤其在跨声速区域颤振临界速度会急剧下降，存在明显的“跨声速凹坑”现象^[2-4]，此问题更值得关注。而由于目前颤振建模和工程分析方法^[5]还无法准确考虑跨声速范围内的激波运动和激波诱导分离等非线性因素，因此风洞颤振试验成为研究跨声速颤振机理和预测颤振边界的重要手段。

风洞颤振试验系统包括单独部件、组合体部件及全模颤振试验结构。随着飞行器各部件的刚度越来越接近，部件之间的耦合作用无法忽视，因而全模颤振成为一项不可或缺的试验技术。开展全模颤振试验不仅要求支撑系统具有较低的刚体模态频率，即典型颤振模型结构最低阶弹性一阶模态频率的1/10~1/3^[6]以减少支撑对全机颤振特性的影响，且还应实现能实时调整模型姿态的功能，以确保模型安全。因此，设计具有稳定性、符合要求的全模颤振支撑系统一直是研究的重点。

在全模颤振试验中，为尽量减小支撑刚度对飞行器模态参数测量的影响，需要尽可能使用柔软的弹性支撑方式。目前，悬浮或悬挂式绳系支撑是全模颤振试验主要采用的形式。典型的有美国航空航天局(National Aeronautics and Space Administration，NASA)兰利中心于1960年在跨声速动力学风洞(transonic dynamic tunnel, TDT)中建成的专门用于研究飞行器气动弹性的双索悬挂系统^[7-9]。该支撑系统具有准确度高、系统误差较小等优点，且可通过控制舵面动态调整模型位姿，还具备对静不稳定模型的支撑能力。但目前仅能在公开发表的文献中了解此类悬挂系统的基本原理、结构组成及稳定性分析等^[10]，还无法完全掌握其详细结构和运行过程。此外，俄罗斯中央流体研究院(Central Aerohydrodynamic Institute, TsAGI)采用不同支撑方式，在T-128风洞中进行了跨声速气动特性分析试验，探究了不同支撑方式对试验结果的影响^[11]。国内，多个研究单位开展了低速和跨声速全模颤振支撑的试验研究，如中国空气动力研究与发展中心分别在直径为3.2 m的低速风洞和直径为2.4 m的暂冲式亚跨超风洞FL-26中建立了一套全模颤振悬浮支撑系统(float suspension system，FSS)。唐建平等^[12]在不同低速风洞中证明了将FSS用于低速全模颤振试验合理可行。郭洪涛等^[13]采用FSS支撑方式验证了这种悬浮支撑系统可使颤振模型具有除轴向外的5个方向的运动自由度，并具有较好的稳定性和安全性，且研究了跨声速时的典型颤振特性。钱卫等^[14]针对某全机结构相似跨声速颤振模型，重点进行了跨声速颤振特性仿真，并与基于FSS技术的试验结果进行对比，验证了一种全机复杂耦合的颤振形式。文[15-16]分别提供了一种全机低速颤振模型悬挂系统。赵振军等^[17]设计了一种三索悬挂系统，通过弹簧刚度选择和钢绳张力控制实现了较低的振动频率，并通过建立包含钢绳、气动力及姿态控制系统等的多体动力学模型，进一步仿真分析了系统姿态调整能力。厦门大学课题组前期基于双索悬挂系统研究了在低速情况下，绳系结构和绳索预紧力对该全模颤振悬挂系统刚体模态频率的影响^[18]。

综上所述，在跨声速全模颤振试验方面，国内主要采用FSS、双索和三索等悬挂支撑方式，这3种支撑方式均是可以释放五自由度的软支撑形式。FSS主要用于暂冲式跨声速全模颤振风洞试验，若要应用于连续式跨声速风洞，则该支撑的试验精准度还需提高；而双索和三索支撑的形式还有待进一步试验验证。此外，国内还没有相关文献从理论上深入研究或证明绳系支撑系统稳定，尤其是模型姿态调整时还未考虑舵面作用，尽管这会涉及更复杂的气动伺服弹性问题。发展新型的全模颤振支撑方式，进行详尽的支撑系统稳定性分析(颤振试验中的关键)，均是亟待解决的问题。

本文以连续式跨声速风洞全模颤振试验为背景，在双索支撑研究的基础上，探索一种基于绳系并联机器人技术的四绳支撑方式。首先，建立一种四绳-三弹簧支撑系统，基于运动学分析，推导其刚度表达式；其次，根据系统动力学方程与虚功原理，开展稳定性理论研究；最后，从支撑系统的固有频率、运动空间及动态响应等方面进行分析，验证所提四绳并联支撑方式的有效性和可行性。

1 四绳并联支撑系统建模

根据全模颤振试验的要求，支撑系统应具有足够的刚度和稳定性，能够为刚体，即模型的运动提供滚转、俯仰、偏航、侧滑和升沉5个方向上的自由度，或在软支撑条件下，允许模型在多个自由度上实现一定幅度的运动，同时支撑系统的固有频率应较低，尽可能不影响模型的结构弹性振动模态和刚体运动模态。

考虑上述因素，本文提出一种四绳-三弹簧的支撑机构，如图 1所示，保证系统安全性的辅助绳索不再描述。四绳-三弹簧支撑系统主要由2根竖直绳索(B₁P₁、B₂P₂)、2根水平绳索(B₃P₃、B₄P₄)、3根弹簧、传动装置、预紧装置和驱动控制单元组成。其中弹簧与绳索串联，绳索一端固接在飞机模型上，另一端通过滑轮直接或经由弹簧间接连接到由驱动组件控制的绳索传动装置上。弹簧的作用是：1) 维持预紧装置加载到绳索上的预紧力；2) 在受到突发外部载荷的情况下，依靠弹簧实现机构的被动调节，起到安全保护作用。四绳-三弹簧悬挂系统通过控制绳索拉力和绳索长度保持模型平衡。

1.1 系统运动学模型

系统运动学几何关系如图 2所示，其中黑色“十”字表示飞机模型。以飞机模型质心为原点建立动坐标系O_bX_bY_bZ_b，即机体坐标系，X_b轴正方向指向机头方向，Z_b轴位于模型对称面内，垂直于X_b轴，正方向指向机身下方，Y_b轴与模型对称面垂直，正方向指向模型右舷方向；以机架底部中心为原点建立静坐标系O_gX_gY_gZ_g，Z_g轴正方向竖直向下，X_g正方向指向风洞来流方向；Y_g轴垂直于X_g轴，正方向遵循右手定则，如图 1所示。飞机模型相对静坐标系的姿态角分别用滚转角ϕ、俯仰角θ和偏航角ψ表示。B_i为固定于机架上的第i个定滑轮点(忽略滑轮直径)，P_i为位于飞机模型上的第i个牵引点，B_i和P_i是在静坐标系下的向量；L_i为P_i点到B_i点的绳长矢量；第i根绳的单位方向矢量u_i=L_i/L_i；机体坐标系原点在静坐标下的位置矢量可表示为O_gb=O_gO_b。

由矢量闭合原理可知，L_i满足：

$ \boldsymbol{L}_i=\boldsymbol{B}_i-\boldsymbol{O}_{\mathrm{gb}}-\boldsymbol{R} \boldsymbol{r}_i . $

(1)

其中：r_i为在机体坐标系中P_i相对于原点O_b的矢量，R为由机体坐标系到静坐标系的转换矩阵。

第i根绳的绳长表示为

$ L_i=\sqrt{\boldsymbol{L}_i^{\mathrm{T}} \boldsymbol{L}_i} . $

(2)

定义飞机模型的位姿参数为$ \boldsymbol{H} = (x, y, z, \phi, \theta, \psi)^{\mathrm{T}}$，运动速度矢量$ \dot{\boldsymbol{H}} = \left(\dot{x}, \dot{y}, \dot{z}, \omega_X, \omega_Y, \omega_Z\right)^{\mathrm{T}}$。x、y和z分别为在静坐标系下的坐标参数，φ、θ和ψ分别为动坐标系下绕3条轴线的滚转角、俯仰角和偏航角，$\dot{x}, \dot{y}, \dot{z} $分别为在静坐标系下的线速度，ω_X、ω_Y和ω_Z分别为飞机绕静坐标系的3条轴线的角速度。姿态角速度矢量$ (\dot{\phi}, \dot{\theta}, \dot{\psi})^{\mathrm{T}}$与角速度矢量$\boldsymbol{\omega}=\left(\omega_X, \omega_Y, \omega_Z\right)^{\mathrm{T}} $存在如下关系：

$ \boldsymbol{\omega}=\left[\begin{array}{l} \omega_X \\ \omega_Y \\ \omega_Z \end{array}\right]=\left[\begin{array}{ccc} \cos \theta \cos \psi & -\sin \psi & 0 \\ \cos \theta \sin \psi & \cos \psi & 0 \\ -\sin \theta & 0 & 1 \end{array}\right]\left[\begin{array}{c} \dot{\phi} \\ \dot{\theta} \\ \dot{\psi} \end{array}\right] . $

(3)

其中：$\dot{\phi} 、\dot{\theta} \text { 和 } \dot{\psi} $分别为动坐标系下绕3条轴线的角速度。

根据并联机器人微分运动学^[19]，可得绳长变化速度矢量$\dot{\boldsymbol{L}}_i $与飞机模型运动速度矢量$ \dot{\boldsymbol{H}}$之间的关系表示如下：

$ \dot{\boldsymbol{L}}_i=\boldsymbol{J} \dot{\boldsymbol{H}}, $

(4)

$ \boldsymbol{J}=\left[\begin{array}{cccc} \boldsymbol{u}_1 & \boldsymbol{u}_2 & \cdots & \boldsymbol{u}_4 \\ \boldsymbol{r}_1 \times \boldsymbol{u}_1 & \boldsymbol{r}_2 \times \boldsymbol{u}_2 & \cdots & \boldsymbol{r}_4 \times \boldsymbol{u}_4 \end{array}\right]^{\mathrm{T}}. $

(5)

其中：J为并联支撑系统的Jacobi矩阵(4×6)，与机构布局参数有关。

由于有3根绳索(B₂P₂、B₃P₃和B₄P₄)通过弹簧与传动装置相连，在绳索运动过程中，将拉力传导至弹簧并在绳索伸长方向产生位移，因此需要进行补偿，弹簧绳受力后的伸长量表示如下：

$ \Delta \boldsymbol{L}_i=\left(\boldsymbol{T}_i-\boldsymbol{T}_{i 0}\right) / k_{\mathrm{s}} . $

(6)

其中：T_i为第i根绳的拉力，T_i0为第i根绳的初始预紧力，k_s为弹簧刚度系数。

1.2 系统静力学模型

对于完全由绳索组成的绳系并联支撑机构，飞机模型受到绳拉力、外力(外力矩)及重力的作用，当其处于平衡状态时，满足静力平衡条件：

$ \boldsymbol{J}^{\mathrm{T}} \boldsymbol{T}+\boldsymbol{W}=\bf{0} . $

(7)

其中：T为绳拉力矢量(4×1)；W为飞机模型所受的外力和重力矢量。对于由绳索和弹簧组成的绳系并联支撑机构，可把串联弹簧的绳索看作1根刚度较低的绳索，本质上对Jacobi矩阵并无影响。

由于绳索只能承受拉力而不能受压，而且应该满足一定的强度要求，因此绳索i承受的拉力值T_i应满足：

$ T_{\max } \geqslant T_i \geqslant T_{\min } \geqslant 0, \quad i=1, 2, 3, 4 . $

(8)

其中：T_max与T_min分别是绳索在安全范围内所能承受拉力的最大值与防止发生虚牵现象的最小值，与绳索材料相关。

1.3 系统刚度模型

支撑系统的刚度定义为飞机模型产生单位位姿变化所需的外力或力矩，根据受力平衡方程，该四绳-三弹簧并联机构的刚度K表示为

$ \boldsymbol{K}=-\frac{\partial \boldsymbol{W}}{\partial \boldsymbol{H}}=\frac{\partial\left(\boldsymbol{J}^{\mathrm{T}} \boldsymbol{T}\right)}{\partial \boldsymbol{H}}=\frac{\partial \boldsymbol{J}^{\mathrm{T}}}{\partial \boldsymbol{H}} \boldsymbol{T}+\boldsymbol{J}^{\mathrm{T}} \frac{\partial \boldsymbol{T}}{\partial \boldsymbol{H}}. $

(9)

其中：$\frac{\partial \boldsymbol{J}^{\mathrm{T}}}{\partial \boldsymbol{H}} \boldsymbol{T} $表示主动刚度项，用K₁表示，与绳拉力和支撑系统绳系布局相关；$\frac{\partial \boldsymbol{J}^{\mathrm{T}}}{\partial \boldsymbol{H}} $表示Jacobi矩阵相对位姿矢量的变化，又称为系统Hessian矩阵^[20]；$\boldsymbol{J}^{\mathrm{T}} \frac{\partial \boldsymbol{T}}{\partial \boldsymbol{H}} $为被动刚度项，用K₂表示，主要与Jacobi矩阵、绳索材料和弹簧参数等相关，根据绳拉力与绳长关系，进而可表示为

$ \begin{aligned} \boldsymbol{K}_2 & =\boldsymbol{J}^{\mathrm{T}} \frac{\partial \boldsymbol{T}}{\partial \boldsymbol{L}_i} \frac{\partial \boldsymbol{L}_i}{\partial \boldsymbol{H}}, \\ k_{\mathrm{c} i} & =\frac{\partial T_i}{\partial L_i}=\frac{E A}{L_i}, \\ k_{\mathrm{L} i} & =\frac{\left(k_{\mathrm{s} i} E A\right)}{\left(E A+k_{\mathrm{s} i} L_{\mathrm{s} i}\right)} . \end{aligned} $

(10)

其中：k_ci为绳索的刚度系数，E为绳索的弹性模量，A为绳索横截面积，k_Li为串联弹簧的绳索刚度，k_si为第i根弹簧的刚度系数，L_si为第i根串联弹簧的绳索变形前的长度。当弹簧刚度$ k_{\mathrm{s} i}=\frac{E A}{L_{\mathrm{s} i}}$时，k_Li≈k_si，根据本文弹簧布置方式，可得$\frac{\partial \boldsymbol{T}}{\partial \boldsymbol{L}_i+\Delta \boldsymbol{L}_i}=\operatorname{diag}\left(k_{\mathrm{c} 1}, k_{\mathrm{s} 2}, k_{\mathrm{s} 3}, k_{\mathrm{s} 4}\right) $；又根据式(4)，可得被动刚度K₂，表示如下：

$ \boldsymbol{K}_2=\boldsymbol{J}^T \operatorname{diag}\left(k_{\mathrm{c} 1}, k_{\mathrm{s} 2}, k_{\mathrm{s} 3}, k_{\mathrm{s} 4}\right) \boldsymbol{J} . $

(11)

经过推导，支撑系统刚度K₁和K₂可进一步表示为：

$ \begin{array}{*{20}{c}} {\left\{ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{K}}_1} = \sum\limits_{i = 1}^4 {\frac{{{\mathit{\boldsymbol{T}}_i}}}{{{\mathit{\boldsymbol{L}}_i}}}} \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}}&{{{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]}^{\rm{T}}} - {\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}{{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]}^{\rm{T}}}}\\ {\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right] - \left[ {{{\bf{r}}_i} \times } \right]{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}}&{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]{{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]}^{\rm{T}}} - \left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}{{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]}^{\rm{T}}}} \end{array}} \right]{\rm{ - }}\sum\limits_{i = 1}^4 {{\mathit{\boldsymbol{T}}_i}} \left[ {\begin{array}{*{20}{c}} 0&0\\ 0&{\left[ {{\mathit{\boldsymbol{u}}_i} \times } \right]\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]} \end{array}} \right],}\\ {{\mathit{\boldsymbol{K}}_2} = {k_{{\rm{cl}}}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{u}}_1}\mathit{\boldsymbol{u}}_1^{\rm{T}}}&{{\mathit{\boldsymbol{u}}_1}\mathit{\boldsymbol{u}}_1^{\rm{T}}{{\left[ {{\mathit{\boldsymbol{r}}_{i = 1}} \times } \right]}^{\rm{T}}}}\\ {\left[ {{\mathit{\boldsymbol{r}}_{i = 1}} \times } \right]{\mathit{\boldsymbol{u}}_1}\mathit{\boldsymbol{u}}_1^{\rm{T}}}&{\left[ {{\mathit{\boldsymbol{r}}_{i = 1}} \times } \right]{\mathit{\boldsymbol{u}}_1}\mathit{\boldsymbol{u}}_1^{\rm{T}}{{\left[ {{\mathit{\boldsymbol{r}}_{i = 1}} \times } \right]}^{\rm{T}}}} \end{array}} \right] + \sum\limits_{i = 2}^4 {{k_{{\rm{s}}i}}} \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}}&{{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}{{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]}^{\rm{T}}}}\\ {\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}}&{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]{\mathit{\boldsymbol{u}}_i}\mathit{\boldsymbol{u}}_i^{\rm{T}}{{\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right]}^{\rm{T}}}} \end{array}} \right]{\rm{, }}} \end{array}} \right.}\\ {\left[ {{\mathit{\boldsymbol{r}}_i} \times } \right] = \left[ {\begin{array}{*{20}{c}} 0&{ - {r_Z}}&{{r_Y}}\\ {{r_Z}}&0&{ - {r_X}}\\ { - {r_Y}}&{{r_X}}&0 \end{array}} \right].} \end{array} $

(12)

其中：I为3阶单位矩阵; [r_i×]为向量r_i的叉乘运算, r_X、r_Y和r_Z分别为P_i在X_b、Y_b和Z_b轴上相对于原点O_b的距离。

由式(12)可知，相较于无弹簧的四绳支撑机构，四绳-三弹簧支撑系统的Hessian矩阵没有变化，但引入弹簧会改变绳拉力T，进而影响主动刚度项K₁；同样，弹簧引入也会影响被动刚度项K₂。进一步分析刚度表达式可知：主动刚度中第1项为正定矩阵，第2项子矩阵[u_i×] [r_i×]的特征值为－u_ir_i，是有限值。另外，由式(11)可知被动刚度也是正定矩阵，其中绳索刚度k_ci一般具有较大的数量级，因此K₂在整个刚度中起主导作用，进而保证了K的每个特征值均大于0。

由支撑系统的刚度可进一步得到系统的刚体模态频率，即支撑无阻尼固有圆频率f_cir和无阻尼固有频率f_n，分别表示如下：

$ \left|\boldsymbol{K}-f_{\text {cir }} \boldsymbol{M}\right|=0, $

(13)

$ f_{\mathrm{n}}=\frac{f_{\text {cir }}}{2 \pi} . $

(14)

其中M为飞机模型的惯量矩阵。但值得注意的是，由于K并非严格的对角矩阵，不同方向的刚度存在一定程度的耦合，在求解固有频率时还需结合特征向量具体分析。

2 支撑系统稳定性

支撑系统的稳定性分析是全模颤振技术研究的关键，本节基于动力学方程和稳定性理论，推导和建立系统稳定的条件。

2.1 系统动力学模型

考虑风洞试验实际应用条件，绳索长度和直径相对较小，忽略绳索垂度及动力学特性，且绳索变形影响有限，因此仅将其视为单向受力的线性模型。采用Newton-Euler法，建立飞机模型的动力学方程，其矩阵形式可表示为

$ \boldsymbol{M}(\boldsymbol{H}) \ddot{\boldsymbol{H}}+\boldsymbol{N}(\boldsymbol{H}, \dot{\boldsymbol{H}})-\boldsymbol{w}_{\mathrm{gra}}-\boldsymbol{w}_{\mathrm{e}}=-\boldsymbol{J}^{\mathrm{T}} \boldsymbol{T} . $

(15)

其中：$ \ddot{\boldsymbol{H}}$为飞机运动模型的加速度矢量；$\boldsymbol{N}(\boldsymbol{H}, \dot{\boldsymbol{H}}) $为非线性Coriolis力；w_gra为重力矢量；w_e为气动力或力矩矢量，是状态变量H和dH的函数；在气流坐标系中，气动力一般表示为$Q S\left[-C_{\mathrm{drg}} \quad C_{\mathrm{sid}} \quad -C_{\text {lif }}]^{\mathrm{T}}\right.$；气动力矩表示为$ Q S\left[\begin{array}{ll} l_{\text {win }} C_{\text {rol }} & \bar{l}_{\text {cho }} C_{\text {pit }} & l_{\text {win }} C_{\text {yaw }}\end{array}\right]^{\mathrm{T}}$，通常用动导数气动模型表示飞机模型的气动力和气动力矩；Q为动压，S为参考机翼面积；总气动力和气动力矩系数，即阻力系数C_drg、侧向力系数C_lif、升力系数C_sid、滚转力矩系数C_rol、俯仰力矩系数C_pit和偏航力矩系数C_yaw表示如下：

$ \begin{gathered} C_{\mathrm{drg}}=C_{\mathrm{drg} 0}+C_{\mathrm{drg}}^\alpha \alpha, \\ C_{\mathrm{sid}}=C_{\mathrm{sid}}^\beta \beta, \\ C_{\mathrm{lif}}=C_{\mathrm{lif} 0}+C_{\mathrm{lif}}^\alpha \alpha+C_{\mathrm{lif}}^{\delta_{\mathrm{e}}} \delta_{\mathrm{e}}, \end{gathered} $

(16)

$ \begin{gathered} C_{\mathrm{rol}}=C_{\mathrm{rol}}^\beta \beta+\frac{l_{\mathrm{win}}}{2 V}\left(C_{\mathrm{rol}}^{\omega_X} \omega_X+C_{\mathrm{rol}}^{\omega_Z} \omega_Z\right)+C_{\mathrm{rol}}^{\delta_{\mathrm{a}}} \delta_a+C_{\mathrm{rol}}^{\delta_{\mathrm{r}}} \delta_{\mathrm{r}}, \\ C_{\mathrm{pit}}=C_{\mathrm{pit0}}+C_{\mathrm{pit}}^\alpha \alpha+\frac{l_{\mathrm{cho}}}{2 V}\left(C_{\mathrm{pit}}^{\omega_Y} \omega_Y+C_{\mathrm{pit}}^{\dot{\alpha}} \dot{\alpha}\right)+C_{\mathrm{pit}}^{\delta_{\mathrm{e}}} \delta_{\mathrm{e}}, \\ C_{\mathrm{yaw}}=C_{\mathrm{yaw}}^\beta \beta+\frac{l_{\mathrm{win}}}{2 V}\left(C_{\mathrm{yaw}}^{\omega_X} \omega_X+C_{\mathrm{yaw}}^{\omega_Z} \omega_Z\right)+ \\ C_{\mathrm{yaw}}^{\delta_{\mathrm{a}}} \delta_{\mathrm{a}}+C_{\mathrm{yaw}}^{\delta_{\mathrm{r}}} \delta_{\mathrm{r}} . \end{gathered} $

(17)

其中：α为迎角，$ \dot{\alpha}$为无量纲迎角；β为侧滑角；l_win为机翼展长；$ \bar{l}_{\text {cho }}$为平均气动弦长; V为气流速度：δ_e、δ_r和δ_a分别为升降舵、方向舵和副翼操纵效率；C_drg0、C_lif0和C_pit0分别为零姿态角下阻力、侧向力和升力系数；C_drg^α、C_lif^α和C_pit^α分别为阻力、升力和俯仰力矩系数对迎角的气动静导数；C_sid^β、C_rol^β和C_yaw^β分别为侧向力、滚动力矩和偏航力矩系数对侧滑角的气动静导数；$ C_{\mathrm{pit}}^{\dot{\alpha}}$为俯仰力矩系数对迎角变化率的导数；C_rol^ω_X和C_yaw^ω_X分别为滚转力矩和偏航力矩系数对滚转角速度的动导数系数；C_pit^ω_Y为俯仰力矩系数对俯仰角速度的动导数系数；C_rol^ω_Z和C_yaw^ω_Z分别为滚转力矩和偏航力矩系数对偏航角速度的导数系数；C_lif^δ_e和C_pit^δ_e分别为升力和俯仰力矩系数对升降舵偏角的导数；C_rol^δ_a和C_yaw^δ_a分别为滚转力矩和偏航力矩系数对副翼偏角的导数；C_rol^δ_r和C_yaw^δ_r分别为滚转力矩和偏航力矩系数对方向舵偏角的导数，可根据典型的小振幅强迫振荡试验数据确定。以飞机质心为原点，与来流风向相反的方向为X_air轴的正方向创建气流坐标系，考虑在风洞试验中，气流坐标系X_air轴与静坐标系的X_g轴的方向一致，因此经过转换矩阵推导，可得气流角与飞机模型姿态角之间的关系为：

$ \tan \alpha=\tan \theta \cos \phi+\sin \phi \tan \psi / \cos \theta, $

(18)

$ \sin \beta=\sin \theta \sin \phi \cos \psi-\sin \psi \cos \phi. $

(19)

2.2 系统稳定性条件

由虚位移原理可知：当系统处于平衡位置时，式(15)中前两项为0，所有作用力产生的虚功和为0，即

$ \delta U=-\delta \boldsymbol{L}_i^{\mathrm{T}} \boldsymbol{T}+\delta \boldsymbol{H}^{\mathrm{T}} \boldsymbol{W}=0. $

(20)

其中：δU表示Lagrange能量函数U的一阶变分，δL_i和δH分别表示绳拉力和外力作用产生的虚位移。

暂不考虑对气动力项的变分，函数U的二阶变分可进一步表示为

$ \delta^2 U=-\delta^2 \boldsymbol{L}_i^{\mathrm{T}} \boldsymbol{T}-\delta \boldsymbol{L}_i^{\mathrm{T}} \delta \boldsymbol{T}+\delta^2 \boldsymbol{H}_i^{\mathrm{T}} \boldsymbol{W}. $

(21)

由绳长变化率与Jacobi矩阵关系式δL_i=JδH，可进一步得到绳长的二阶变化量为

$ \delta^2 \boldsymbol{L}_i=\delta \boldsymbol{J} \delta \boldsymbol{H}+\boldsymbol{J} \delta^2 \boldsymbol{H}. $

(22)

将上式代入式(21)，并根据静力平衡关系式，可得：

$ \begin{gathered} \delta^2 U=-\left(\boldsymbol{J} \delta^2 \boldsymbol{H}\right)^{\mathrm{T}} \boldsymbol{T}-(\delta \boldsymbol{J} \delta \boldsymbol{H})^{\mathrm{T}} \boldsymbol{T}- \\ \delta \boldsymbol{L}_i^{\mathrm{T}} \delta \boldsymbol{T}+\delta^2 \boldsymbol{H}^{\mathrm{T}} \boldsymbol{W}=-\delta^2 \boldsymbol{H}^{\mathrm{T}}\left(\boldsymbol{J}^{\mathrm{T}} \boldsymbol{T}-\boldsymbol{W}\right)- \\ \delta \boldsymbol{H}^{\mathrm{T}} \delta \boldsymbol{J}^{\mathrm{T}} \boldsymbol{T}-\delta \boldsymbol{L}_i^{\mathrm{T}} \delta \boldsymbol{T}=-\delta \boldsymbol{H}^{\mathrm{T}} \delta \boldsymbol{J}^{\mathrm{T}} \boldsymbol{T}-\delta \boldsymbol{L}_i^{\mathrm{T}} \delta \boldsymbol{T} . \end{gathered} $

(23)

根据2.1节刚度推导过程可知，$\delta \boldsymbol{J}^{\mathrm{T}} \boldsymbol{T}=\boldsymbol{K}_1 \delta \boldsymbol{H}, \delta \boldsymbol{T}=\boldsymbol{K}_2 \delta \boldsymbol{L}_i=\boldsymbol{K}_2 \boldsymbol{J} \delta \boldsymbol{H} $，则式(23)可进一步表示为

$ \delta^2 U=-\delta \boldsymbol{H}^{\mathrm{T}}\left(\boldsymbol{K}_1+\boldsymbol{J}^{\mathrm{T}} \boldsymbol{K}_2 \boldsymbol{J}\right) \delta \boldsymbol{H}=-\delta \boldsymbol{H}^{\mathrm{T}} \boldsymbol{K} \delta \boldsymbol{H} . $

(24)

由此可知，二次虚功实质上可表示为关于负刚度矩阵的二次型，而由于刚度矩阵又可表示为正定矩阵与一般有限值矩阵的线性组合，且被动刚度起主要作用，J^TK₂J的特征值始终为正，且远大于K₁的特征值，并高出几个量级，因此可保证K是正定的，即式(24)在非零状态始终为负。根据二次虚位移原理，该二次型恒负是系统平衡且稳定的充分必要条件。因此，在四绳并联支撑作用下，系统在平衡点是稳定的。

进一步考虑对气动力或力矩项的变分，式(24)可表示为

$\begin{aligned} \delta^2 U & =-\delta \boldsymbol{H}^{\mathrm{T}} \boldsymbol{K} \delta \boldsymbol{X}+\delta \boldsymbol{H}^{\mathrm{T}} \delta \boldsymbol{W}= \\ & -\delta \boldsymbol{H}^{\mathrm{T}}\left(\boldsymbol{K}-\frac{\delta \boldsymbol{W}}{\delta \boldsymbol{H}}\right) \delta \boldsymbol{H} . \end{aligned} $

(25)

其中$ \frac{\delta \boldsymbol{W}}{\delta \boldsymbol{H}}$可等效为气动静导数矩阵。

显然，当飞机模型为静稳定状态，即C_pit^α < 0，C_rol^β < 0时，系统稳定性得到进一步增强；即使不满足静稳定条件，即C_pit^α>0，C_rol^β>0时，只要刚度矩阵的特征值大于0，仍然能保证系统的稳定性。显然对于绳系并联支撑而言，被动刚度J^TK₂J起主要作用，只需保证其特征值高于气动导数矩阵的特征值，则全模颤振支撑系统可保持稳定，即

$ \boldsymbol{K}-\frac{\delta \boldsymbol{W}}{\delta \boldsymbol{H}}>0. $

(26)

3 仿真分析与验证 3.1 已知条件

本节针对所提出的四绳-三弹簧机构进行运动特性仿真与分析，以0.6 m连续式跨声速风洞为应用背景，动压范围为0.6×10⁴~8.5×10⁴ Pa，设计支撑布局与模型参数，具体参数设置如下：取类F-22飞机模型(见图 3)，质量为2.5 kg，等效面积为0.085 m²，弦长为0.2 m，模型关于质心的主惯量矩为[0.137，0.882，0.994] kg·m²。绳索材料根据受力情况选择常用于风洞试验支撑的钢丝绳，绳索直径取4 mm，弹性模量为206 GPa；根据固有频率估算，3根弹簧刚度系数分别取k_s2=0.100k_c=393 N/mm，k_s3=k_s4=0.008k_c=31 N/mm，k_c为绳索对应刚度。根据飞机模型特点，采用前2根绳索水平、后2根绳索竖直的方式布置。

四绳-三弹簧支撑系统与机构连接点坐标如表 1所示，可知初始绳长为[0.659 4，0.659 4，0.672 2，0.672 2] m。

在初始位置条件下，由静力平衡式(7)，可得初始预紧力为[1 086.8，864.6，891.6，891.6] N。根据上述参数，可进一步得到绳索变形为亚毫米量级，因此忽略形变影响是合理的。

3.2 稳定性分析与气动响应特性仿真

首先根据系统稳定性判断条件式(26)，计算刚度矩阵和气动导数矩阵的特征值，进行理论判断；进一步开展在气动力作用下的系统动态响应仿真与分析。

在初始平衡状态，系统主动刚度K₁与总刚度矩阵K分别为：

$ \boldsymbol{K}_1=\left[\begin{array}{cccccc} 1002 & 0 & -50 & 0 & -25 & 0 \\ 0 & 5158 & 0 & 9 & 0 & 203 \\ -50 & 0 & 5064 & 0 & -165 & 0 \\ 0 & 9 & 0 & 35 & 0 & 5 \\ -25 & 0 & -165 & 0 & 340 & 0 \\ 0 & 203 & 0 & 5 & 0 & 343 \end{array}\right], $

$ \boldsymbol{K}=\left[\begin{array}{cccccc} 359150 & 0 & -135420 & 0 & 9580 & 0 \\ 0 & 1570 & 0 & 0 & 0 & 20 \\ -135420 & 0 & 78400 & 0 & -5490 & 0 \\ 0 & 0 & 0 & 1 & 0 & 3 \\ 9580 & 0 & -5490 & 0 & 420 & 0 \\ 0 & 20 & 0 & 3 & 0 & 30 \end{array}\right] \times 10. $

支撑系统总刚度矩阵对角线元素分别对应模型的平动刚度(N/m)与转动刚度(N·m/rad)。从数值来看，在此构型下，飞机模型沿轴的平动刚度明显高于其他2个方向的刚度；滚转方向的转动刚度相对较弱，这与全模颤振试验对模型运动自由度的约束要求相一致，即支撑系统提供模型除来流方向外的其他5个自由度。进一步，可求得总刚度矩阵K的特征值为：[414 110，23 820，1 569，34，34，4]×10。

由文[21]提供的标模气动试验数据^[21]可知，该飞机模型的静导数参数如表 2所示。由表 2可知，C_pit^α为正，说明纵向静不稳定；C_rol^β为负，C_yaw^β为正，说明横滚和航向静稳定。

代入动压、平均气动弦长等参数，进一步可得气动参数对应的矩阵为

$ \frac{\delta \boldsymbol{W}}{\delta \boldsymbol{H}}=\left[\begin{array}{cccccc} 0 & 0 & 0 & 49 & -1129 & 49 \\ 0 & 0 & 0 & 2868 & 54 & 2868 \\ 0 & 0 & 0 & -3724 & 1 & -3724 \\ 0 & 0 & 0 & 7 & 259 & 7 \\ 0 & 0 & 0 & 91 & 0 & 91 \end{array}\right]. $

进一步可得$ \boldsymbol{K}-\frac{\delta \boldsymbol{W}}{\delta \boldsymbol{H}}$的特征值为[414 090，23 730，1 566，117，371，34]×10，均为正，且特征值远大于0，证明即使在飞机模型静不稳定的情况下，采用该绳系支撑系统，仍然能满足稳定性要求。

需要强调的是，尽管本文针对给定的绳索布局和飞行器模型参数，证明了该四绳-三弹簧支撑系统是稳定的，但从刚度矩阵特性分析可知，被动刚度的特征值起主要作用，因此对于静稳定状态下的飞机，支撑系统是稳定的；对于静不稳定状态下的飞机，需根据稳定性判断条件式(26)，并进一步计算和判断，式(26)能为支撑系统设计提供指导。

此外，基于系统动力学方程式(15)，考虑气动力项，开展动力学响应仿真，飞行器模型质心位移、姿态角变化及支撑系统绳拉力变化如图 4—6所示。由图 4和5可知，在高速来流作用下，四绳-三弹簧并联支撑的飞行器模型Y_b轴方向质心位移最大，但小于0.04 m；滚转角和偏航角变化幅值不超过15.0°，俯仰角变化范围不超过5.0°，进一步验证了该绳系支撑系统具有良好的稳定性。

此外，由绳拉力变化可知，在振动过程中，竖直向上的绳拉力T₁、水平牵引的绳拉力T₃和T₄始终处于张紧状态，且最大值在安全范围内；在飞机模型处于某些极值位姿时，竖直向下的绳拉力T₂会处于瞬时虚牵状态，但由位姿变化可知，这并不影响系统的稳定性和安全性，同样更有利于减小支撑对模型的约束，不过还是可以通过增大初始预紧力避免此现象。需要指出的是，系统在理论建模时，并未考虑绳索和弹簧的阻尼特性，因此飞机模型在有界范围内振荡，且振荡没有明显衰减。

3.3 刚体模态频率及影响因素

支撑系统刚体模态频率是颤振试验首要关注的指标。为避免支撑系统影响模型结构弹性模态，通常要求支撑系统的刚体模态频率低于模型第1阶弹性模态固有频率的1/3，即11.5 Hz。

将SOLIDWORKS三维模型导入ANSYS Workbench软件中，材料设置为结构钢，密度为7 850 kg/m，Young模量为200 GPa，Poisson比为0.3，在Mesh模块中采用四面体单元划分飞机模型网格，设置体划分尺寸为40 mm，部分面尺寸为1 mm，使用自适应尺寸调整，物理场设置为“Mechanical”，求解器类型选择“subspace”，计算前10阶固有频率，分析本文飞机模型进行自由-自由状态下的模态，结果表明前6阶振型是无约束状态下沿X_b、Y_b和Z_b轴的平移和旋转刚体运动，应当舍去，因此可得飞机模型一阶弹性模态固有频率为34.6 Hz(见图 7)，这意味着支撑系统的刚体模态频率需限制在11.5 Hz以内。

模型初始平衡后，在时间t=1 s时对不同方向分别施加外部激励，可得系统的响应曲线并进一步做频域转换，如图 8和9所示。

由图 9可知，在该构型下转动方向，即滚转、俯仰和偏航方向的固有频率分别为0.8、1.0和0.9 Hz。同样地，根据位置冲击响应频率，可以得到3个平动方向频率，分别是升沉方向1.4 Hz，侧向9.6 Hz，来流方向44.6 Hz，如图 10所示。除来流方向外，均能满足全模颤振对支撑系统刚体模态频率的要求。

为进一步验证上述结果的正确性，通过Adams软件建立四绳-三弹簧支撑系统模型(见图 11)，参数设置参考3.1节，进行系统振动频率分析。

Adams中的Vibration模块可用于分析支撑系统的振动性能。在飞机模型质心位置施加一定幅值、相位角为0°且频率为0.01~100.00 Hz的正弦扫频信号，然后对支撑系统进行强迫振动分析，可得四绳-三弹簧支撑系统的刚体模态频率，结果显示支撑系统的滚转方向频率为0.81 Hz，偏航方向频率为0.95 Hz，俯仰方向频率为1.03 Hz，与数值计算结果吻合，以俯仰方向为例，振动时序图如图 12所示。

与经典的两绳、三绳支撑方式不同，对于四绳并联支撑，在平衡位姿情况下，绳拉力不存在多解问题，即如果改变绳拉力，飞机模型的平衡位姿将发生变化，或处于另一新的平衡状态，故不再单独讨论绳拉力对系统刚体模态频率的影响，而是结合绳索布局参数及模型位姿等，进行综合分析。

1) 不同牵引位置的影响。

分别改变绳索与飞机模型连接点P_i位置参数中的x_Pi和y_Pi坐标，以考察其对支撑系统俯仰和滚转方向固有频率的影响。本文主要改变P₃和P₄，例如原构型x_P3=-50，x_P4=50，新增2个算例，即x_P3=-100, -150, x_P4=100, 150，其他坐标不变；类似原构型y_P3=-22，y_P4=22，新增2个算例，即y_P3=-62, -102，y_P4=62, 102，单位为mm，其他坐标不变。通过施加外部激励，对响应曲线做频域转换，可得不同位置对应的固有频率，如表 3所示。

由表 3可知，当绳索连接点在飞机模型沿X_b方向增大间距时，俯仰方向对应的固有频率适当增大，意味着绳索对产生的俯仰约束力矩变大。但当(x_P3, x_P4)分别为(-100, 100)和(-150, 150)时，对应的频率并未发生变化，这是由于虽然力臂增加，但对应的平衡绳拉力会一定程度减小。其中，前者平衡预紧力为[78.8, 45.3, 62.9, 62.9] N；后者为[63.1, 24.2, 45.2, 45.2] N。同样，当绳索连接点在飞机模型沿Y_b轴方向增大间距时，可得绳索产生的滚转力矩变大，使该方向的固有频率逐渐增大。但整体而言，俯仰和滚转方向固有频率均小于1.4 Hz, 支撑系统的固有频率都处于全模颤振试验的要求范围内。

2) 不同弹簧数量(刚度)的影响。

弹簧数量(刚度)将会改变支撑系统的频率，但由于弹簧与绳索串联，弹簧数量不同，本质上相当于弹簧刚度发生变化。

在本文所采用的四绳-三弹簧基础上，本节又计算了四绳、四绳-一弹簧(k_s2或k_s4=0.100k_c)和四绳-二弹簧(k_s2=k_s4=0.100k_c)算例，计算得到其对应的固有频率，并进行对比分析，结果如表 4所示。

由表 4可知，对四绳并联支撑系统而言，其转动方向固有频率小于3.0 Hz，是否采用弹簧对滚动和升沉方向的频率影响相对不大，均满足试验要求。但对侧向方向固有频率的影响很大，与无弹簧的四绳构型相比，三弹簧构型极大地降低了侧向的频率，满足固有模态低于11.5 Hz的要求；当只加入一根弹簧时，采用水平放置方式比竖直方式放置更有效。综合而言，采用弹簧不仅是为了满足系统安全要求，即提供一定的缓冲，而且更为重要的目的是须满足支撑系统频率的要求。

3.4 姿态调整范围

全模颤振试验中，当飞机模型受到来流冲击时，位姿会出现一定程度的偏离。当位姿变化较大时，必须及时调整，使其能快速恢复至安全范围内，避免对模型和风洞的破坏。针对该四绳支撑系统，可采用控制绳长或操纵舵面调姿。分别阐述如下：

1) 基于绳长变化的位姿调整方式。

该绳系支撑方式属于并联机构，根据绳长与模型位姿的关系式及约束条件，可计算出模型运动范围，具体求解步骤和流程参考文[22]。由3.2节气动响应特性分析可知，四绳-三弹簧支撑系统受到来流冲击力后，飞机模型在3个平动方向上的位移较小，姿态角变化较大，因此本节重点分析支撑系统姿态角的调整范围。考虑到避免绳索与机体结构产生碰撞干涉及全模颤振试验要求具备一定的姿态角调整能力^[13]，因此，设置姿态角的搜寻范围为：滚转角[-40.0°, 40.0°]，俯仰角[-40.0°, 40.0°]，偏航角[-30.0°, 30.0°]；令绳索拉力范围为0~200 N，以保持绳索张紧，其他系统参数与3.1节一致，对离散点进行遍历计算得出姿态调整范围，仿真结果如图 13所示。

图 13中坐标点是能够读取到的边界点，由边界点数值可知，在该B_i、P_i布局下，滚转角调整范围约为[-40.0°, 40.0°]；俯仰角调整范围约为[-9.5°, 25.0°]，这是由俯仰方向绳索-弹簧支链刚度差别导致的；偏航角调整范围为[-22.0°, 22.0°]。在Adams中通过调整绳长和施加脉冲力的方式也可得到模型各姿态角调整范围，图 14为动态试验中飞机模型姿态角调整范围，测量结果表明：四绳-三弹簧支撑系统飞机模型的滚转角调整范围超过20.0°，俯仰角最大可达23.0°，偏航角最大可达22.0°。

综上所述，通过MATLAB的数学建模理论计算和基于Adams的软件仿真结果对比可以看出，动态试验中的四绳-三弹簧支撑系统姿态角调整范围能满足全模颤振试验的姿态调整范围要求。

2) 基于舵面偏转角的姿态控制方式。

由3.2节可知，在吹风过程中当发生较大角度偏移时，可通过调节绳长，恢复至安全范围，同样加入舵偏指令，也可以调节飞机模型的姿态角。例如针对图 5中无舵偏时滚转角和偏航角较大的现象，可将副翼调整-0.3°，方向舵调整-1.0°，升降舵调整-2.0°，仿真结果如图 15所示，可知滚转角和偏航角变化降低幅度大，最大值约为1.0°，俯仰角最大约为3.2°，有效实现了试验过程中对飞机模型位姿的调控。

3.5 试验验证

考虑到实际风洞试验条件下的弹簧刚度和绳索直径较大，需特别订制。为进一步验证理论分析结果，本文仅在实验室搭建简易的四绳-三弹簧支撑地面样机，如图 16所示，采用60 mm×60 mm×1 300 mm的铝合金型材搭建机架外框(1.3 m×1.3 m×1.3 m)以模拟风洞试验段，类F-22飞机模型质量参数如3.1节所述，弹簧刚度系数分别取k_s2=2 777 N/m、k_s3=k_s4=925 N/m；用红色加粗线段示意绳索，直径为1 mm，采用Kevlar绳，预紧力为[81.6, 61.5, 71.3, 71.3] N。

采用锤击法进行模态频率测试试验，在飞机模型静止后，分别沿滚转、俯仰、偏航3个方向敲击，激励模型自由振荡，此时飞机位姿发生变化，使用AHRS位姿传感器采集姿态角响应信号，经过Fourier变换进一步得到支撑系统3个转动方向的固有频率，结果如图 17所示。

结果表明：飞机模型分别受滚转、俯仰和偏航方向激励时，滚转方向模态频率约为1.5 Hz，俯仰方向模态频率约为2.3 Hz，偏航方向模态频率约为1.6 Hz。试验结果与仿真结果稍有差别，一是由于支撑系统实际的B_i、P_i与仿真参数有所差异，二是由于实际试验中的弹簧刚度与绳索弹性模量参数不同导致固有频率有所偏差。但总体而言，仍然可以验证所提支撑方式能够满足全模颤振固有频率的要求。

4 结论

1) 针对跨声速全模颤振支撑技术需求，提出一种基于并联机构的四绳-三弹簧牵引支撑方式，推导了该系支撑系统的稳定性条件，证明了该四绳-三弹簧支撑适用于静不稳定飞机模型。

2) 由数值模拟与Adams仿真可知，通过调整绳索连接点及结构参数，四绳-三弹簧支撑系统固有频率可控制在11.5 Hz以下，其中转动方向固有频率小于3.0 Hz；姿态角调整范围超过20.0°；经由舵面也可以有效控制飞机模型姿态，满足全模颤振试验的要求；此外，通过简易样机实验验证了理论方法可行。

3) 所提四绳-三弹簧并联支撑方式为全模颤振试验提供了一种新的手段，可为下一步风洞试验研究提供技术支持。