0 引言

目前飞行力学分析和飞行控制系统设计时主要采用线性叠加方法建立气动力模型，构建气动力模型的数据主要通过静态测力试验和小振幅动导数试验获得。在线性迎角范围内，动导数与迎角、振幅及频率等关系不大；而在大迎角区域，迎角、振幅及频率对动导数的非线性影响明显加大^[1]。现有设计流程没有考虑这些非线性气动力，因此需要通过原型机飞行试验对气动力和飞行控制系统的性能进行验证。目前获取动态气动力的试验装置存在一定的局限性：动导数装置仅能在单个自由度上强迫振动，其幅值小、频率调节范围有限；旋转天平装置仅能实现绕速度矢作锥形运动。而研究大迎角飞行出现的气动力问题，如机翼摇滚、横向偏离、上仰等，就需要把试验装置的某些自由度放开来复现某些复杂飞行现象^[2]。比如自由摇滚、自由偏航等试验装置较好地解决了摇晃、偏离等非线性飞行现象的模拟和研究。大迎角气动问题具有非线性、纵横向耦合强等特点，因此为了尽量逼真地模拟飞行器飞行机动动作，新的多自由度试验装置的研究一直在进行之中。近年来发展的一些两自由度装置模拟了两自由度耦合振荡运动、旋转流场下的振荡运动状态^[3]，多自由度机动模拟装置在风洞中实现了对大飞行包线机动动作的模拟和多自由度耦合运动的研究^[4-5]。

另外，近年来发展的风洞虚拟飞行试验技术^[6-8]不仅要求多自由度耦合、模型自由转动、流场环境逼真；还要求在模型上集成飞行控制设备，使模型运动具备实时可控的能力^[9-10]。虚拟飞行试验技术中模型支撑机构设计是其中的关键技术之一。支撑机构要求气动干扰小、摩擦影响小、角度范围大、模拟更逼真。美国Magill等采用6根张线和球铰的组合方式支撑BOA导弹进行试验，结果与飞行试验结果吻合较好^[11-12]。英国Lowenberg等发展了单自由度、两自由度、三自由度和五自由度等机动试验装置^[13-15]，开展了支杆-飞机多体动力学建模、非线性气动力、参数辨识等方面的研究，获得了试验装置的动力学模型描述，以及Hawk模型和M2370模型三轴的静、动导数数据^[16-17]。俄罗斯Sohi采用腹撑三自由度机构在水平风洞研究稳定尾旋性能，与立式风洞试验结果较为吻合^[18]。TsAGI采用背撑三自由度机构开展了机翼摇滚问题的研究，通过H∞控制技术抑制了试验模型的机翼摇滚现象^[19-20]。中国空气动力研究与发展中心在2.4m跨声速风洞采用纵向单吊臂的支撑方式开展了某导弹虚拟飞行试验，验证了俯仰/滚转耦合运动的解耦控制方法^[21]。

综合分析以上动态试验的支撑方式可以看出，采用直支杆腹撑方式的三自由度装置迎角模拟范围、滚转运动范围有限，难以满足大迎角气动/运动耦合特性的研究。本文提出了一种三自由度装置，能够实现大角度范围的机动动作模拟。针对此支撑装置，采用拉格朗日乘子法导出该虚拟飞行试验装置的动力学数学模型，并进行典型激励下纵向、横航向动力学响应以及大迎角混合操纵响应的非线性仿真与分析。

1 试验装置 1.1 装置介绍

图 1所示的三自由度虚拟飞行试验装置通过置于模型内部的万向铰实现飞机模型的俯仰和偏航运动，通过与旋转曲杆一起转动实现飞机模型的滚转运动。旋转曲杆经过三次预弯，一是为了扩展模型的运动范围，二是为了减小曲杆惯量对系统的影响。旋转曲杆与模型腹部连接的部分进行了第三次预弯，其偏斜角即该部分曲杆轴线(图 1中虚线所示的偏航轴)与水平面的夹角为45°，这样可以减小大迎角状态下模型的腹部开孔。旋转曲杆上方设计了具有流线外形的配重，其高度可调节，可使曲杆重心在曲杆的旋转轴上以保持动平衡。试验模型的三轴姿态及角速率可以通过机载的航姿参考系统进行测量，迎角、侧滑角可利用风洞来流方向不变的特点由三轴姿态推算导出。模型的升降舵、副翼、方向舵等舵面与舵机连接采用平行四边形传动方式；平行于舵面转轴在机身或机翼内部布置舵机转轴，舵面和舵机转轴上设计摇臂，通过推杆连接摇臂推动舵面转动。该装置能实现接近90°迎角的极限飞行动作模拟，可开展失速偏离、尾旋初始阶段和稳定尾旋的研究。该装置的主要的几何尺寸见图 2。

1.2 技术分析

飞机的三自由度角运动是快变化过程，相对慢变化的线运动，其与飞行安全的关系更加紧密。现目前立式风洞尾旋试验、水平风洞自由飞试验中模型的线位移十分有限，其实质也主要为三自由度角运动。因此，三自由度动态试验装置可以捕获飞机三轴转动运动的主要特征，以达到开展气动/运动耦合、多轴运动耦合研究的目的。

采用常规直支杆腹撑或背撑方式实现大角度范围的运动模拟，需要在模型表面开设较大的孔，这样将严重影响模型的气动特性。而采取翼型主支杆和预弯曲杆的方式，既大幅减小了对后方流场的干扰，又可以很好控制对模型外形的破坏。当然，由于支撑的存在会对模型的气动力产生一定的干扰，比如支撑会在模型区诱导一定的上洗气流，从而对模型C_m0的产生一个平移量等。支架干扰的贡献量可以通过CFD计算模拟或风洞支架干扰试验等手段获取，并在气动力建模中予以考虑。

绕速度矢横滚是飞机最基本的机动动作之一。常规腹撑或背撑方式无法实现360°连续滚转，但利用风洞的气流方向保持不变的特点，通过一个小巧的曲杆与模型同步滚转，就可以巧妙地解决绕速度矢旋转的问题。另外，置于模型内部的万向铰可以实现模型迎角、侧滑角的连续变化。这样，在可控的风洞流场环境下就实现了物理意义明晰的三自由度运动。

2 数学模型 2.1 前提和假设

风洞试验时，模型质心位于风洞试验段中心保持不变，但模型的姿态可绕三轴转动，因此，在动力学建模时可作如下假设：

(1) 假设飞机模型及试验装置均为刚体，不存在变形情况；

(2) 假定曲杆所受气动力为零，曲杆的旋转轴平行于风洞来流方向；

(3) 模型的转动中心为其重心，且模型和曲杆的重心在曲杆的旋转轴上；

(4) 试验模型的惯量和旋转速率较小，可忽略陀螺力矩的影响。

2.2 轴系定义

在系统动力学分析时需涉及三种坐标系(如图 3所示)，分别为：

地轴系ox_gy_gz_g，简称S_g，其oz_g沿铅垂方向向下，ox_g在水平面内，并与风洞轴线平行，与来流方向相反，oy_g与平面ox_gz_g垂直，指向右。

模型体轴系ox_ay_az_a，简称S_a，原点位于模型重心，ox_a轴在飞行器对称平面内，平行于机身轴线或机翼的平均气动弦线，指向前；oz_a位于对称面内，垂直于ox_a，指向下；oy_a垂直于对称平面，指向右。

曲杆体轴系ox_ry_rz_r，简称S_r，原点位于曲杆重心，轴系定义与S_a类似。

用φ_a、θ_a、ψ_a分别表示模型滚转角、俯仰角、偏航角，导出S_a与S_g的转换矩阵T_ga、T_ag，且有T_ag=T_ga^T。

用φ_r、θ_r、ψ_r分别表示曲杆滚转角、俯仰角、偏航角，导出S_r与S_g的转换矩阵T_gr、T_rg。

2.3 基于绝对坐标方法的动力学方程 2.3.1 绝对坐标

利用绝对坐标方法进行动力学分析，首先假设系统内所有刚体为无约束的自由刚体，以各刚体的质心笛卡尔坐标和绕质心转动的角度坐标或欧拉参数作为系统的绝对坐标，对各刚体建立无约束的动力学方程，再利用拉格朗日乘子法将其与系统铰约束方程联立，构成完整的动力学方程。

系统的绝对坐标为：

$ \mathit{\boldsymbol{q}} = {(\mathit{\boldsymbol{q}}_{\rm{r}}^{\rm{T}}\;\;\;\mathit{\boldsymbol{q}}_a^{\rm{T}})^{\rm{T}}} = {({\phi _{\rm{r}}}\;\;\;{\theta _{\rm{r}}}\;\;\;{\psi _{\rm{r}}}\;\;\;{\theta _{\rm{a}}}\;\;\;{\psi _{\rm{a}}})^{\rm{T}}} $

刚体基于自身体轴的角速度向量为：

$ {\mathit{\boldsymbol{\omega }}_i} = {({\mathit{p}_i}\;\;\;{\mathit{q}_i}\;\;\;{\mathit{r}_i})^{\rm{T}}},i = {\rm{r}}\mathit{,}{\rm{a}} $

2.3.2 动力学方程的一般形式

确定坐标后，根据拉格朗日乘子法，可以给出受约束的动力学方程一般形式如下^[22]：

$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{A}}&{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q^{\rm{T}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{\rm{q}}}}&{\bf{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\bf{\ddot q}}}\\ {\bf{ \pmb{\mathit{ λ}} }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{B}}\\ \zeta \end{array}} \right] $

(1)

其中，A是与无约束动力方程有关的系数矩阵，B是该方程组右侧向量，Φ_q为约束方程Φ(q)的Jacobian矩阵，λ为拉格朗日乘子，ζ为加速度约束方程组右侧向量。下面将给出该动力学方程的推导过程。

2.3.3 无约束的动力学方程

ω_i是绕体轴转动的角速度，由于体轴是动系，${{{\mathit{\boldsymbol{\dot{q}}}}}_{i}}$与ω_i并不等价，而是有如下关系：

$ {\mathit{\boldsymbol{\omega }}_i} = {\mathit{\boldsymbol{D}}_i}{{\mathit{\boldsymbol{\dot q}}}_i} $

(2)

上式中D_i是与${{{\mathit{\boldsymbol{\dot{q}}}}}_{i}}$相关的矩阵，具体为：

$ {\mathit{\boldsymbol{D}}_i} = \left[ {\begin{array}{*{20}{c}} 1&0&{ - \sin {\theta _i}}\\ 0&{\cos {\phi _i}}&{\sin {\phi _i}\cos {\theta _i}}\\ 0&{ - \sin {\phi _i}}&{\cos {\phi _i}\cos {\theta _i}} \end{array}} \right] $

${\mathit{\boldsymbol{\dot \omega }}_i}$与${{{\mathit{\boldsymbol{\ddot{q}}}}}_{i}}$的关系如下：

$ {\mathit{\boldsymbol{\dot \omega }}_i} = {{\mathit{\boldsymbol{\dot D}}}_i}{{\mathit{\boldsymbol{\dot q}}}_i} + {\mathit{\boldsymbol{D}}_i}{{\mathit{\boldsymbol{\ddot q}}}_i} $

(3)

$ {{\mathit{\boldsymbol{\dot D}}}_i} = \left[ {\begin{array}{*{20}{c}} 0&0&{ - {{\dot \theta }_i}\cos {\theta _i}}\\ 0&{ - {{\dot \phi }_i}\sin {\phi _i}}&{{{\dot \phi }_i}\cos {\phi _i}\cos {\theta _i} - {{\dot \theta }_i}\sin {\phi _i}\sin {\theta _i}}\\ 0&{ - {{\dot \phi }_i}\cos {\phi _i}}&{ - {{\dot \phi }_i}\sin {\phi _i}\cos {\theta _i} - {{\dot \theta }_i}\cos {\phi _i}\sin {\theta _i}} \end{array}} \right] $

一般情况下，刚体R_i绕质心转动的运动方程为如下形式：

$ {\mathit{\boldsymbol{J}}_i} \cdot {{\mathit{\boldsymbol{\dot \omega }}}_{{\rm{ab}}i}} + {\mathit{\boldsymbol{\omega }}_{{\rm{ab}}i}} \times ({\mathit{\boldsymbol{J}}_i} \cdot {\mathit{\boldsymbol{\omega }}_{{\rm{ab}}i}}) = {\mathit{\boldsymbol{M}}_i} $

(4)

式中ω_abi为刚体R_i的绝对角速度，其可以表示为：

$ {\mathit{\boldsymbol{\omega }}_{{\rm{ab}}i}} = {\mathit{\boldsymbol{\omega }}_0} + {\mathit{\boldsymbol{\omega }}_i} $

(5)

式中ω₀为多体系统零刚体的角速度，系统选定支撑座为零刚体，因此ω₀=0，式(5) 可简化为：

$ {\mathit{\boldsymbol{\omega }}_{{\rm{ab}}i}} = {\mathit{\boldsymbol{\omega }}_i} = {\mathit{\boldsymbol{D}}_i}{{\mathit{\boldsymbol{\dot q}}}_i} $

(6)

将上式代入式(4) 后得到：

$ {\mathit{\boldsymbol{J}}_i} \cdot {\mathit{\boldsymbol{D}}_i}{{\mathit{\boldsymbol{\ddot q}}}_i} = {\mathit{\boldsymbol{M}}_i} - {\mathit{\boldsymbol{J}}_i} \cdot {{\mathit{\boldsymbol{\dot D}}}_i}{{\mathit{\boldsymbol{\dot q}}}_i} - {{\mathit{\boldsymbol{\dot D}}}_i}{{\mathit{\boldsymbol{\dot q}}}_i} \times ({\mathit{\boldsymbol{J}}_i} \cdot {\mathit{\boldsymbol{D}}_i}{{\mathit{\boldsymbol{\dot q}}}_i}) $

(7)

此处引入反对称阵的概念：设有向量a=(a₁  a₂  a₃)^T, 则a的反对称阵为：

$ \mathit{\boldsymbol{\tilde a = }}{\left[ {\begin{array}{*{20}{c}} 0&{ - {a_3}}&{{a_2}}\\ {{a_3}}&0&{ - {a_1}}\\ { - {a_2}}&{{a_1}}&0 \end{array}} \right]^{\rm{T}}} $

则式(7) 可以整理为：

$ {\mathit{\boldsymbol{J}}_i} \cdot {\mathit{\boldsymbol{D}}_i}{{\mathit{\boldsymbol{\ddot q}}}_i} = {\mathit{\boldsymbol{M}}_i} - ({\mathit{\boldsymbol{J}}_i}{{\mathit{\boldsymbol{\dot D}}}_i} + \widetilde {{{\mathit{\boldsymbol{\dot D}}}_i}{{\mathit{\boldsymbol{\dot q}}}_i}}{\mathit{\boldsymbol{J}}_i}{\mathit{\boldsymbol{D}}_i}){{\mathit{\boldsymbol{\dot q}}}_i} $

(8)

令：${{\mathit{\boldsymbol{A}}}_{i}}={{\mathit{\boldsymbol{J}}}_{i}}{{\mathit{\boldsymbol{D}}}_{i}},{{\mathit{\boldsymbol{B}}}_{i}}={{\mathit{\boldsymbol{M}}}_{i}}-\left( {{\mathit{\boldsymbol{J}}}_{i}}{{{\mathit{\boldsymbol{\dot{D}}}}}_{i}}+\widetilde{{{\mathit{\boldsymbol{D}}}_{i}}{{{\mathit{\boldsymbol{\dot{q}}}}}_{i}}}{{\mathit{\boldsymbol{J}}}_{i}}{{\mathit{\boldsymbol{D}}}_{i}} \right){{\mathit{\boldsymbol{\dot{q}}}}_{i}},$ i=r, a，定义方阵A=diag(A_r   A_a)和列阵B=(B_r^T  B_a^T)^T，则得到系统的无约束运动学方程组：

$ \mathit{\boldsymbol{A\ddot q}} - \mathit{\boldsymbol{B = }}0 $

(9)

2.3.4 约束条件

系统中存在两种约束：第一种是约束曲杆仅有一个旋转自由度的旋转铰约束；第二种是连接曲杆和模型的十字万向铰约束，即万向铰的偏航旋转轴单位向量e_z^c与俯仰旋转轴单位向量e_y^c始终垂直，而e_y^c与模型体轴系S_a的轴Oy_a重合，设e_y^a为轴Oy_a上的单位向量，则有e_y^a·e_z^c=0。

因此，本系统的几何约束方程组Φ(q)为：

$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(\mathit{\boldsymbol{q}}) = \left[ \begin{array}{l} {\theta _{\rm{r}}}\\ {\psi _{\rm{r}}}\\ \mathit{\boldsymbol{e}}_y^{\rm{a}} \cdot \mathit{\boldsymbol{e}}_z^{\rm{c}} \end{array} \right],\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}(\mathit{\boldsymbol{q}}) = {\bf{0}} $

(10)

上述约束条件为不显含时间的定常约束，上式对t求导即可得到速度约束方程：

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}\mathit{\boldsymbol{\dot {q} }} = {\bf{0}} $

(11)

式中Φ_q为约束方程Φ(q)对坐标q的雅可比矩阵。

速度约束方程再次对时间求导，得到加速度约束方程：

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}\mathit{\boldsymbol{\ddot q = }}{({\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}\mathit{\boldsymbol{\dot q}})_q}\mathit{\boldsymbol{\dot q}} = {\bf{0}} $

(12)

令ζ=-(Φ_q$\mathit{\boldsymbol{\dot{q}}}$)_q$\mathit{\boldsymbol{\dot{q}}}$，上式整理为：

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}\mathit{\boldsymbol{\ddot q = }}\zeta $

(13)

2.3.5 拉格朗日乘子法

将无约束的运动方程与坐标的高斯加速度变分δ${\mathit{\boldsymbol{\ddot{q}}}}$相乘，得到：

$ {(\mathit{\boldsymbol{A\ddot q}} - \mathit{\boldsymbol{B}})^{\rm{T}}}\mathit{\delta }\mathit{\boldsymbol{\ddot q}} = {\bf{0}} $

(14)

其中的变分δ${\mathit{\boldsymbol{\ddot{q}}}}$并非独立变量，必须满足3个加速度约束方程的限制，即：

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}\mathit{\delta }\mathit{\boldsymbol{\ddot q}} = {\bf{0}} $

(15)

引入3个拉格朗日乘子λ=(λ₁  λ₂  λ₃)^T，将上述约束方程与相同标号的拉格朗日乘子相乘，加入动力学方程中，得到下式：

$ {(\mathit{\boldsymbol{A\ddot q}} - \mathit{\boldsymbol{B + \boldsymbol{\varPhi} }}_q^{\rm{T}}\mathit{\boldsymbol{\lambda }})^{\rm{T}}}\mathit{\delta }\mathit{\boldsymbol{\ddot q}} = {\bf{0}} $

(16)

令δ$\mathit{\boldsymbol{\ddot q}}$中3个不独立的变分的系数为零，上式中仅剩余与3个独立变分有关的和式。方程成立的充要条件是独立变分的系数为零，即括弧内所有元素均为零，这样导出第一类拉格朗日方程如下：

$ \mathit{\boldsymbol{A\ddot q}} = \mathit{\boldsymbol{B - \boldsymbol{\varPhi} }}_q^{\rm{T}}\mathit{\boldsymbol{\lambda }} $

(17)

与加速度约束方程联立，得到受约束的动力学方程：

$ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{A}}&{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q^{\rm{T}}}\\ {{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_q}}&0 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{\ddot q}}}\\ \mathit{\boldsymbol{\lambda }} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{B}}\\ \zeta \end{array}} \right] $

2.4 力和力矩 2.4.1 气动力模型

本文基于常规风洞试验得到的基本气动力、舵效及动导数等数据建立了如下的气动力模型，各系数均在体轴系下描述，气动力模型中支架干扰的贡献量可通过后续补充风洞试验获得。该飞机模型舵面全零状态下的基本气动力矩见图 4所示。

$ ({\mathit{\boldsymbol{F}}_A},{\mathit{\boldsymbol{M}}_A}) = f(\alpha ,\beta ,{\mathit{\boldsymbol{\omega }}_{\rm{a}}},\mathit{V,}{\delta _{\rm{e}}},{\delta _{\rm{r}}},{\delta _{\rm{f}}},{\delta _{\rm{a}}}) $

(18)

式中V为来流速度，δ_e、δ_r、δ_f、δ_a分别表示模型升降舵、方向舵、襟翼和副翼角度，各气动系数按照基本量叠加增量的方式计算。以俯仰力矩系数C_m为例，按照下式进行计算：

$ \begin{array}{l} {C_m} = {C_{m{\rm{basic}}}}\left( {\alpha + {\delta _{\rm{f}}}} \right) + \Delta {C_{m\beta }}(\alpha ,\beta ,{\delta _{\rm{f}}}) + \Delta {C_{m\sup }}\left( {\alpha ,\beta } \right) + \\ \;\;\;\;\;\;\;\;\Delta {C_{m\dot \alpha }}\left( \alpha \right)\dot {\alpha} + {\Delta} {C_{mq}} + \Delta {C_{m{\delta _{\rm{e}}}}}(\alpha ,\beta ,{\delta _{\rm{e}}}) \end{array} $

(19)

式中，C_mbasic为与迎角α、襟翼δ_f有关的基本量，ΔC_mβ是由侧滑角β引起的增量，ΔC_msup是支撑装置带来的支架干扰量，ΔC_{m${\mathit{\boldsymbol{\dot{a}}}}$}和ΔC_mq是升沉运动和俯仰运动引起的增量，ΔC_mδe是升降舵带来的增量。偏航力矩系数C_n、滚转力矩系数C_l的计算与C_m类似。

2.4.2 摩擦力矩计算

常用的动态摩擦模型有Dahl模型、Bristle模型、Bliman and Sorine模型、LuGre模型和Leuven模型^[23-24]。其中LuGre模型能精确的描述摩擦静态和动态特性，因此选择LuGre模型的一种简化形式对装置中两组铰摩擦力矩M_cafⁱ、M_Rrfⁱ进行计算，前者为万向铰的摩擦力矩，后者为旋转铰的摩擦力矩。LuGre模型如下：

$ \frac{{{\rm{d}}{Z_j}}}{{{\rm{d}}t}} = {{\dot q}_j} - \frac{{{\sigma _0}\left| {{{\dot q}_j}} \right|}}{{S\left( {{{\dot q}_j}} \right)}}Z $

(20)

$ S\left( {{{\dot q}_j}} \right) = {T_{\rm{c}}} + ({T_{\rm{s}}} - {T_{\rm{c}}}){e^{ - {\alpha _{\rm{f}}}\left| {{{\dot q}_j}} \right|}} $

(21)

$ {T_{\rm{f}}} = {\sigma _0}Z + {\sigma _1}{{\rm{e}}^{ - {{({{\dot q}_j}/{\mathit{v}_d})}^2}}}\frac{{{\rm{d}}{\mathit{Z}_j}}}{{{\rm{d}}t}} + {\sigma _2}{{\dot q}_j} $

(22)

式中，σ₀为刚性系数，N·m/rad；σ₁为阻尼系数，σ₂为粘性系数，N·m·s/rad；Z为刚毛平均形变，T_s为静摩擦力矩，T_c为库仑摩擦力矩，N·m；参数α_f的引入使该模型符合Stribeck效应，单位为s/m；参数ν_d用以减弱高转速下刚毛偏转的影响，并且保证了摩擦力矩模型的钝性，单位为rad/s。

当系统处于稳态时，dZ_j/dt=0，此时摩擦力矩模型可以简化为式(23)。图 5给出了该装置俯仰运动方向摩擦特性的地面测试结果。

$ {T_{\rm{f}}} = {\mathop{\rm sgn}} ({{\dot q}_j})({\mathit{T}_c} + ({\mathit{T}_{\rm{s}}} - {\mathit{T}_{\rm{c}}}){e^{ - {\alpha _{\rm{f}}}\left| {{{\dot q}_j}} \right|}}) + {\sigma _2}{{\dot q}_j} $

(23)

3 仿真与结果分析 3.1 仿真程序

仿真程序在Matlab环境下编写，采用Nelder-Mead算法在给定初始状态下配平飞机；导出约束方程的Jacobian矩阵的解析形式，与无约束的动力学方程联立构成动力学模型，然后采用无违约算法进行求解^[25-27]；采用多维线性插值方法获取气动力，根据LuGre模型解算摩擦力矩。图 6给出了计算程序的相关模块和主要流程。

为了分析曲杆、摩擦力矩等因素对飞机模型响应特性的影响规律，选取三种工况进行对比研究，即Ⅰ.仅试验模型的三自由度运动(无摩擦), Ⅱ.试验模型和曲杆的三自由度运动(无摩擦), Ⅲ.试验模型和曲杆的三自由度运动(有摩擦)。通过Ⅰ和Ⅱ的对比，可以获得曲杆对系统的影响规律；通过Ⅱ和Ⅲ的对比，可以获得摩擦对系统的影响规律。

3.2 算法验证

Adams是美国MSC公司开发的可用于开展多体动力学仿真的成熟软件。为验证前述动力学模型与算法的正确性，利用Adams-Matlab联合仿真功能进行确认。采用较为常用的斜拉杆操纵进行对比(升降舵-10°，副翼-3°)，两种方法的对比验证结果见图 7，图中Rig为本文算法的结果。不同算法下，纵向响应的一致性很好；横航向响应基本吻合，曲线形态相似，只是其动态振荡过程略有差异。总的来说，不同算法的差异可以接受，达到了算法验证的目的。后文在此基础上利用本文的方法对试验装置进行详细分析。

3.3 典型操纵

为了便于分析，针对纵、横、航向分别采用典型的阶跃信号激励飞机的动态响应；纵向、横向及航向的操纵信号分别为1°升降舵阶跃、1°副翼阶跃和1°方向舵阶跃信号。同时给出一组大迎角混合操纵的响应结果，各舵面输入信号为：在1s时输入-30°升降舵阶跃，4s时输入6°方向舵阶跃，12s时输入10°副翼阶跃。飞机模型的机翼后掠角约30°，其主要参数见表 1。仿真初始状态的主要参数见表 2。图 8~图 11给出了相应的动力学响应仿真结果。图中标注3-DOF为基于飞机模型三自由度转动方程无杆状态的仿真结果，Rig为飞机模型在虚拟飞行试验装置上的仿真结果，Rig+Friction为在该装置基础上增加摩擦力影响的仿真结果。

3.4 纵向阶跃操纵响应分析

升降舵操纵的响应结果如图 8所示，有杆、无杆仿真结果吻合较好。从理论分析，纯纵向运动中曲杆不会对俯仰运动的响应产生影响，因此模型与机构均无滚转、偏航等横航向运动。

但摩擦力矩对俯仰角、俯仰速率均有一定影响。对俯仰速率的影响主要体现在阻尼作用和死区现象。阻尼作用体现在俯仰速率峰值变小，在第一、二个峰值分别比无摩擦仿真结果减少0.11°/s和0.25°/s；死区现象表现为俯仰速率曲线在第二次过零时(约1.8s)不再振荡衰减，而是自此时以后一直为零。与之相应的，俯仰角响应峰值略有减小，并在1.8s后不再变化，提前0.7s进入稳定状态，其稳态值与无摩擦状态相比约有0.02°的微小差量。

3.5 横向阶跃操纵响应分析

副翼操纵的响应结果如图 9所示，曲杆使模型滚转速率、偏航速率响应变慢，速率响应振荡稍平缓，稳态值差异较小；曲杆带来模型滚转角响应差量约0.8°(仿真结束时)。模型侧滑角、偏航角响应差异不大。曲杆滚转角与模型滚转角的响应基本同步一致，最大差量约0.3°。

摩擦力矩对速率响应的动态过程和稳态值均有一定影响。模型滚转速率峰值较无摩擦结果小0.4°/s，稳态值存在0.4°/s的差量；偏航速率峰值较无摩擦结果相差约0.13°/s，稳态值差量为0.04°/s。由于模型滚转角由速率积分决定，有无摩擦两种状态滚转角响应出现明显的斜率差异，滚转角最大差量达3.4°；模型偏航角最大相差约0.18°，侧滑角稳态值相差约0.05°。

3.6 航向阶跃操纵响应分析

方向舵操纵的响应结果如图 10所示，三种状态下总体响应趋势一致，但也存在一定差异。曲杆使偏航速率振荡峰值先减小后增加，峰值时间略有推迟，稳态值变化不大；与之对应，曲杆使侧滑角、偏航角响应的振荡幅值也略增大。此种差异主要来自于万向铰偏航旋转轴的偏斜，曲杆惯量的影响次之。

曲杆使模型滚转速率峰值减小，峰值时间稍推迟，稳态值基本不变；但操纵初始阶段响应曲线形态差异明显。在操纵后0.2s内无杆状态下模型先正滚转后减速变为负滚转，而有杆状态下模型直接表现为负滚转。模型滚转角响应的差异与速率差异吻合，二者在仿真结束时存在0.26°的差量。另一方面模型滚转角与曲杆滚转角在操纵初始时运动趋势相反，曲杆滚转角正滚转最大达0.8°。结合加速度曲线分析发现，曲杆的存在使滚转运动对应惯量增大，导致滚转模态时间常数变大，使得滚转角速率响应变慢，加上万向铰偏航旋转轴偏斜的影响，导致了滚转速率在1~5s的较大差异；而万向铰偏航旋转轴的偏斜使偏航运动分解为绕该轴的旋转和绕速度矢的旋转，后者造成了曲杆与模型在初始响应时滚转方向不一致。

摩擦力矩的影响主要表现为阻尼作用，模型的偏航速率和滚转速率收敛更快，其中偏航速率在3.6~3.8s时呈现了死区现象，偏航速率稳态值差量约0.04°/s，使偏航角斜率略有差异，在仿真结束10s时差量约0.3°，侧滑角差量约0.05°；滚转速率稳态值与无摩擦结果存在约0.4°/s的差量，使滚转角响应的斜率呈现明显差别，仿真结束时滚转角相差约3.2°。

3.7 大迎角混合操纵响应分析

如图 11所示，输入升降舵信号后，模型快速抬头，迎角在1.3s时达到峰值，其后迅速收敛稳定在28.6°；输入方向舵信号后，模型4.6s时达到侧滑角峰值后略有收敛但振荡剧烈，均值在8°左右，同时产生约-49°/s的滚转速率和-31°/s的偏航速率，模型开始绕速度矢连续滚转；输入副翼信号后，滚转速率和偏航速率继续负增长，即模型绕速度矢加速滚转。

与前面小迎角状态相比，模型的响应特性有所不同。主要体现在：一是纵向和横航向之间的耦合更为严重。操纵方向舵后，模型迎角、俯仰速率出现小幅振荡，且振荡形态与侧滑角类似；这主要是由于模型随侧滑变大在俯仰方向会产生一定的上仰力矩(见图 4、图 11)。二是典型模态的响应特性和操纵性发生变化。纵向运动方面，小迎角下单位舵偏产生的迎角稳态响应为1.1°，而大迎角状态下为0.73°，与升降舵效率随迎角增大而降低的规律吻合。航向运动方面，小迎角下单位舵偏的侧滑稳态响应为1.4°，大迎角状态下为1.2°，量值相当；但大迎角状态下模型荷兰滚模态响应特性恶化、振荡加剧，这是由于随迎角增大航向稳定性降低、动态阻尼下降所致。横向运动方面，小迎角下单位舵偏滚转速率的稳态响应为-8.7°/s，而大迎角状态下为-5.4°/s，与副翼舵效、滚转阻尼随迎角和侧滑角增加而显著降低的规律相符。

大迎角状态下，机构摩擦对响应结果的稳态值有微小影响，其影响规律与小迎角状态响应结果一致。从图中还发现，大迎角状态下曲杆对模型航向操纵响应的影响明显增大，其差异仍来自于万向铰偏航旋转轴的偏斜以及曲杆惯量，但与小迎角状态不同的是，此时曲杆惯量为主要影响因素，偏航旋转轴的偏斜次之。这是因为：偏航旋转轴的偏斜主要影响偏航运动，曲杆惯量主要影响滚转运动，并且随着模型迎角增大，曲杆惯量对于滚转运动的影响加大。而在航向操纵响应中偏航运动和滚转运动又是相互耦合的，最终就导致了侧滑响应的明显差异。目前设计的曲杆惯量约占飞机滚转惯量的20%，小迎角状态下差异不大；但在大迎角状态，侧滑响应峰值相差约30%。因此曲杆惯量最好在20%以内甚至更小，以保证试验模拟的合理性和相似性。

4 结论

本文基于绝对坐标方法建立了一种虚拟飞行试验装置的动力学模型；基于该模型所获得的典型操纵响应结果趋势合理，量值正确，表明该动力学模型可从理论上有效指导该类试验装置的设计与实现。通过仿真分析得到：

1) 曲杆对试验模型滚转角速率响应变慢；对偏航运动的振荡过程有一定的增强效果，该影响随模型迎角增大而变大；但曲杆对稳态响应影响不大。因此建议旋转曲杆应尽量采用轻质材料制作，其惯量最好控制在飞机滚转惯量的20%以内，以降低对模型操纵响应的影响。

2) 摩擦力矩对短周期、荷兰滚等模态的动态过程和各特征量的稳态值有一定影响，在低旋转速度下可能产生死区现象。因此该装置应采用摩擦贡献较小的精密轴承，以尽量减小摩擦力矩对模型操纵响应的影响。

本文提出的装置尚未开展风洞试验；后续将对现有试验装置和模型进行适应性改造，以开展大迎角风洞虚拟飞行试验进一步验证数学模型和算法。