翼尖铰接复合飞行器动力学特性研究

刘东旭，谢长川，洪冠新

（北京航空航天大学航空科学与工程学院，北京 100083）

作为一种特定的新概念飞行器系统，翼尖铰接复合飞行技术具体是以单架次分别进行起飞爬升和下降着陆，在巡航过程中通过特定设计的翼尖铰柱接合机构实现飞机复合飞行，这种技术能够有效降低飞机的油耗率，提升飞机的航程和航时等续航性能指标。对于远距离长航时飞机，这种技术的优点具体体现在：可以增加展弦比、翼型弯度，同时飞机还会受到彼此翼尖涡流的上洗抬升作用以增加低速巡航性能；多架飞行器可以在苛刻条件下以单一架次起降，在空中对接后再执行飞行任务和作战任务，这样就降低了飞机的起降环境要求。从经济性方面考虑，气动性能的提升意味着复合飞行系统可以节省燃油的布置空间和携带质量，在总体设计工作流程方面提高飞机的设计灵活性。作为一种提升续航性能的飞行技术，翼尖对接复合技术经常被用来和空中加油技术进行比较。翼尖对接复合技术在一定程度上使作战飞机可以免于空中加油机的辅助燃油保障，从而简化作战出动飞机的机种，改善机队的出动任务规划方法和起降保障作业的流程工作。

翼尖铰接复合飞行［1］的概念最初由德国科学家在二战期间提出，此后美国空军和业内各研究机构制定了Tom-Tom项目、Ficon项目为主的军方试验计划［2］验证这个技术理论，同时设计研发了多个翼尖铰接复合飞行概念验证机型号。曾在美国康奈尔大学和德雷伯实验室工作的Edgar教授通过使用涡格法、CFD法等数值仿真方法［3］和风洞测试模拟仿真方法［4］开展复合机翼的气动特性研究，验证了这类复合飞行技术的气动性能改善效果。美国空军从事的“Ficon”计划设计了一个由战略轰炸机组成的进攻单元，由1架B-29轰炸机和2架F-84 D／E，战斗机以翼尖铰接形式复合飞行；而在风洞测试、地效实验以及几次成功的对接实验等验证工作之后，“Ficon”计划的飞行试验遭遇了严重失败。由于缺少协作飞行策略，其中一架飞机从对接的飞机脱离，该项目组在对接后仅几秒内就坠海失事［2］。此后翼尖铰接复合飞行技术的动力学研究集中在带转动铰约束条件的多刚体系统动力学问题。多刚体系统指的是多个刚体按照确定的方式相互联系所组成的机械系统。多刚体系统动力学是在经典力学基础上产生的新学科分支。处理由多个刚体组成的系统，原则上可以利用传统的经典力学方法，即以Newton-Euler方程为代表的矢量力学方法或以拉格朗日方法为代表的分析力学方法。同时以此为基础，Roberson-Wittenburg方法主要是利用图论的概念及数学工具描述多刚体系统的结构，以邻接刚体之间的相对位移作为广义坐标，导出适合于任意多刚体系统的普遍动力学方程，利用费舍尔的增广体概念对方程的系数矩阵A作出物理解释［5］。在此基础上，佛罗里达大学的Chakravarthy针对变后掠翼飞行器［6］，上海交通大学的安继光［7］、北京航空航天大学的乐挺和王立新教授［8］等针对Z形折叠翼面飞行器［7-8］、非对称斜翼飞行器［9-10］和变后掠翼飞行器［11-12］开展了动力学建模和仿真工作，Thomas［13-14］、Li［15］与Yan［16］等针对类似形式的多种多刚体系统飞行器也开展了建模仿真研究，为多刚体系统飞行器动力学建模和仿真方法提供了参考。对于这种特定的翼尖铰接复合飞行器，Brandon和Montalvo［17］给出了按单体飞机添加约束的传统牛顿欧拉动力学建模方法，同时使用增益调度方法完成了控制特性研究，又对对接飞行轨迹展开了探讨，但对于多个飞行器组成的复合系统基于单个飞机的气动数据和配平方案，没有建立多机复合的方法给出多机的配平方案。

为研究翼尖铰接复合飞行器的动力学特性，本文以2架中空长航时无人机通过翼尖铰接机构组成的复合飞行器系统作为研究对象，考虑复合飞行过程中系统整体和内部的运动特点，利用 Newton-Euler方法、Lagrange方法和Roberson-Wittenburg方法开展系统动力学建模工作，使用CFD方法建立面向复合运动特征参数和复合飞行舵面部件的双机复合飞行器系统的气动力数据库。在此基础上研究复合飞行器的本征动力学特性和动力学响应，研究结果可以为复合飞行器的动力学特性和配平策略的设计提供参考依据。

1 复合动力学建模

1.1 复合飞行器系统

算例采用的复合飞行器系统是由2架常规尾翼布局中空长航时无人机通过翼尖铰接机构实现复合飞行，三维实体设计如图1所示。

翼尖铰接机构主要由铰柱约束机构组成，如图2所示。翼尖铰接约束机构可以限制2架飞机相互之间的3个相对平动自由度和2个相对转动自由度，具体是前后、上下和左右方向上的相对平动运动以及俯仰相对转动、偏航相对转动，同时允许2架飞机之间存在相对滚转运动。

图2 翼尖铰接机构Fig.2 Wingtip joint mechanism

1.2 复合系统动力学和运动学方程

翼尖铰接复合飞行器的研究属于带转动铰约束条件的无根树系统动力学问题，可以根据自由度和约束形式确定列写方程的数量。

在欧美制单刚体机体坐标轴系下，分析飞机的运动学与动力学问题。2架飞机分别各自拥有6个自由度，具体为x、y、z轴的3个平动自由度和3个转动自由度。如图3所示，在翼尖一侧添加铰链约束后，约束了2架飞机的x、y、z轴3个轴向自由度和偏航、俯仰2个转动自由度。因此，2架飞机共有6＋6－3－2＝7种自由度，具体为：整体x、y、z轴方向的3个轴向自由度、1个俯仰转动自由度、1个偏航转动自由度和单独的2个滚转自由度。根据多刚体系统一般动力学原理，建立7自由度力学模型，解决7自由度的动力学问题，需要按照2架飞机的7个自由度建立7个动力学方程才可求解，而每个方向的动力学方程都会需要具体分析才能完整建立。

图3 复合飞行器系统运动分析Fig.3 Composite aircraft system’s motion analysis

推导过程分为2部分：①采用整体法推导复合飞行器系统整体的六自由度动力学方程和运动学方程；②采用分离法推导复合飞行器系统内部的相对转动动力学方程。

1.2.1 系统整体动力学方程和运动学方程

根据复合飞行的运动受力状态分析，重新定义一种适用于复合飞行运动的机体坐标系——复合机体坐标系Oxfyfzf，如图4所示。

图4 翼尖铰接复合飞行器机体坐标系Fig.4 Airframe coordinate system of wingtip-jointed composite aircraft

原点位于两复合飞机机体的铰接点O。Oxf轴在复合飞机对称平面内，平行于各单体飞机机身轴线，指向前；Ozf轴亦在对称平面内，垂直于Oxf轴，指向下；Oyf轴垂直于对称平面，指向右。这样选取对应的坐标系可以直接表示复合飞机整体的运动和姿态，便于飞行器驾驶员或控制系统操纵监视过程；单一架次飞机的气动力和气动力矩也可直接通过矩阵变换以整体形式表示，重力和发动机推力的表示方法也不繁琐。同时，复合飞机整体质心C位于2架飞机质心C1和C2连线的中点，OC可由2飞机铰接角度表出。

参照以上复合机体坐标系对应定义机体复合机体气流坐标系Oxfayfazfa、复合机体航迹坐标系Oxfkyfkzfk和常规大地坐标系Oxeyeze等。

根据一般刚体动力学方法，2架飞机各自在复合机体坐标系Oxfyfzf下的速度可以表示为

vi＝vo＋˙rOCi＋ωi×rOCii＝1，2（1）

式中：vo为O点的速度；rOCi为第i个刚体质心C点相对O点的矢经。

式中：mi为各飞机的质量；vi为各飞机的飞行速度；ω为复合机体坐标系相对于惯性坐标系的转动角速度；ICi为各飞机相对整体质心的转动惯量张量；ωir为各飞机相对复合机体坐标系的转动角速度；S为飞行器系统整体相对质心的静矩。

根据已有数学关系，得到复合飞行器多刚体系统整体动力学方程的标量形式

式中：Lef为复合机体坐标系向大地坐标系转换的转换矩阵；φ、θ、ψ为大地坐标系和复合机体坐标系之间的3个欧拉角；u、v、w为复合机体坐标系上x、y、z三个坐标轴方向上的速度；p、q、r在飞行动力学学科中默认为滚转角速度、俯仰角速度、偏航角速度。

1.2.2 系统内部动力学方程和运动学方程

在飞机复合飞行的过程中，存在多刚体相对转动，如图5所示，φ12为两复合飞机的相对滚转角（右侧飞机滚转角为φ1，左侧飞机滚转角为φ2）。可得系统内部的相对转动动力学方程为

图5 翼尖铰接复合飞行器内部运动分析Fig.5 Internal motion analysis of wingtip-jointed composite aircraft

式中：n为复合的飞机架次；g为重力加速度；mk、Vkx、φk、zk分别为第k架飞机的质量、速度、滚转角、z向坐标。

将动势代入再考虑复合铰接角度和飞行器姿态的变换关系φ1＝φ－φ12后，将式（11）展开得

进一步整理可得采用分离法得到的多刚体内部相对转动参数和复合整体参数之间存在的运动关系表达式：

建立复合飞行多体动力学模型得到最终表达形式，7个自由度坐标方向下的零阶量分别为（xf、yf、zf、φf、θf、ψf、φ12）。根据飞行器动力学基本原理，在复合机体坐标系Oxfyfzf下，得到复合系统整体和内部动力学方程（见式（14）），T为发动机推力。

对应有运动学方程和转动动力学方程以及系统内部运动学方程为

2 CFD方法建立气动力数据库

由于翼尖铰接复合飞行器系统飞行运动区别于单机飞行的特殊运动形式，相关的CFD数值仿真和风洞模拟仿真试验［5］证明复合飞行器系统采用2架飞机以翼尖铰接的形式实现复合飞行时，2架飞机的流场会产生气动耦合效应，从而气动特性相比单机飞行时存在显著变化，同时2架飞机之间的相对滚转造成了飞行器系统的气动参数和气动舵面控制之间高度耦合无法解耦，因此需要使用CFD方法建立复合飞行气动力数据库。

建立气动力数据库时要考虑复合运动特征参数（主要是复合系统内部的相对滚转角φ12）和2架飞机各自的舵面部件偏转角度（δa1，δa2）、（δe1，δe2）和（δr1，δr2）对气动特性的影响。考虑复合运动特征参数的原因是：①不同大小的相对滚转角会对复合飞行器气动参数产生影响；②便于研究关联飞行性能指标的复合运动特征参数最优值。考虑复合飞行时舵面部件偏转角度的原因是：①相比单机常规飞行时，副翼、升降舵和方向舵舵面偏转角度对复合飞行器的影响机理更加复杂；②复合机体配平偏转影响。风洞模拟试验［5］证明复合飞行时机翼存在的内外侧翼面压差使内外侧机翼承受不相等的气动升力，可能要采用不同于单机飞行时的双机舵面差动配平方案，需要研究这种配平舵偏带来的复合飞行气动特性变化；③便于精确研究复合运动特征参数的飞行性能指标最优值。

在建立面向复合运动特征参数和复合飞行舵面部件的气动数据库之前，由多刚体飞行器系统的相关气动耦合效应研究［10-11］仿真结果中可知，复合飞行过程里非定常效应对复合飞行动态特性的影响较小，同时在复合飞机相对滚转角角速度较小时，能够忽略非定常气动力效应，即忽略所有气动力与相对滚转角速度˙φ12的气动导数影响，因此可以利用以上气动力准定常假设来简化计算机翼相对滚转过程中复合飞行器系统的气动力。假设工况动压为k，翼面积为S，气动几何平均弦长为ˉc，复合机翼展长为b，则复合飞行过程中气动力和气动力矩表达式为

式中：CD0为复合机体零升阻力系数；CDα和CDV分别为阻力系数对α和V的无因次导数；CL0为迎角为0时的升力系数；CLV、CLq、CLα分别为升力系数对V、q、α的无因次导数；Cm0为零升俯仰力矩系数；Cmα、CmV分别为俯仰力矩系数对α、V的无因次导数；Cmq和Cm˙α分别为俯仰阻尼导数和洗流时差导数；CDδa1δa2＝fDδa（δa1，δa2）、CDδe1δe2＝fDδe（δe1，δe2）、CDδr1δr2＝fDδr（δr1，δr2）等为复合机体舵面操纵导数，其中（δa1，δa2）、（δe1，δe2）和（δr1，δr2）为复合飞行器系统的舵面偏角二维数组，由双机复合飞行时的舵面配平方案决定；y＝fDδa（x1，x2），y＝fDδe（x1，x2）和y＝fDδr（x1，x2）等函数关系由CFD仿真得到的离散数据点采用二元非线性回归插值获得。

具体地，建立气动数据库的CFD仿真工作内容、气动参数和建模方法具体如表1所示。

表1 气动力数据库获取方法Table 1 Aerodynamic database acquisition method

建立气动力数据库时，复合飞行器系统的静导数和操纵导数计算可以认为是定常流动，采用稳态流动的CFD数值模拟方法，计算过程具体如下：通过复合飞行器系统三维实体建模得到定常流动计算所需的网格文件和边界条件，采用ANSYS R15.0 Fluent软件进行数值模拟迭代求解；复合飞行器巡航速度为声速，因此采用压力基求解器进行稳态流动迭代求解。湍流模型采用航空外流场计算常用的对受到逆压梯度作用下的边界层模拟效果较好的Spalart-Allmaras模型，求解算法使用Simple算法；计算入口边界条件采用压力远场，采用二阶迎风格式对控制方程进行离散推进，并最终通过监测复合机翼上下表面压力平均值的稳定性作为计算结果的收敛判定，复合飞行器系统非结构网格具体如图6所示。建立气动力数据库时，复合飞行器系统的静导数和操纵导数计算可以认为是定常流动，采用稳态流动的CFD数值模拟方法，计算入口边界条件采用压力远场，采用二阶迎风格式对控制方程进行离散推进，并最终通过监测复合机翼上下表面压力平均值的稳定性作为计算结果的收敛判定。动导数CLq、CL˙α、Cmq、Cm˙α等主要由Zaero软件的非定常气动力计算方法获得。

图6 复合飞行器系统ANSYS非结构网格划分Fig.6 Unstructured grid generation of composite aircraft system by ANSYS

经过基本气动参数、静导数、动导数、各舵面偏转角度与复合运动特征参数即相对滚转姿态的CFD仿真计算后，就可以建立得到考虑气动舵面偏转和复合运动特征参数的复合飞行气动数据库。

3 动力学特性

3.1 仿真方法

根据第2节中推导得到的7自由度复合机体动力学方程全量模型，在MATLAB／Simulink平台上搭建复合飞行器系统动力学仿真模块，如图7所示。其中飞机多体系统状态模块主要是计算复合飞行器系统的增广体静矩、增广体张量和增广体质量；飞行器动力学模型计算模块主要负责根据飞机多体状态参数、运动姿态参数和舵面输入来确定动力学方程中力和力矩的计算；多刚体系统非线性方程组求解器主要负责计算下一步计算过程中的运动参数和姿态参数；控制系统反馈模块主要负责根据选取好的复合飞行器配平方案，给定对应舵面偏转参数对应的数值。迭代求解器采用ODE45格式。

图7 翼尖铰接复合飞行器系统动力学仿真模块Fig.7 Dynamic simulation module of wingtip-jointed composite aircraft system

使用全量方程搭建复合飞行器系统的动力学仿真平台，可以开展未配平、准配平和全配平方案下的动力学响应研究。未配平方案指复合飞行器各控制舵面（主要有油门舵偏角、升降舵偏角、方向舵偏角和副翼偏角等）输入量为任意值，不需要使复合飞行器整体运动和内部相对运动动力学方程平衡的舵面控制输入方案；准配平方案指复合飞行器各控制舵面需要满足复合飞行器系统整体运动动力学方程平衡，不需要满足内部相对运动动力学方程平衡的舵面控制输入方案；全配平方案指复合飞行器系统各舵面控制输入量需要同时满足复合飞行器系统整体运动和内部相对运动动力学方程的平衡。

3.2 准配平方案和全配平方案的动力学响应

在算例工况高度5000 m，马赫数0.282（速度95.91 m／s），相对滚转角0°和迎角2.27°条件下，研究准配平方案和全配平方案下的双机复合飞行器系统算例的动力学响应。

准配平方案具体为：为保持复合飞行器整体处于定直平飞状态，使用升降舵偏角、油门杆位移来平衡复合飞行器系统任务飞行高度、速度，同时副翼输入来平衡复合机体系统的滚转力矩，方向舵来配平复合飞行器系统的偏航力矩。按以上策略在算例工况下实现准配平方案如表2所示，保证复合飞行器整体动力学方程中所有初始状态导数（˙xi除外）小于10－8，则检查的所有配置都满足准配平要求；仿真开始后飞机控制舵面在整体运动方程中始终保持系统初始稳定状态的输入。

全配平配平方案具体为：为保持复合飞行器整体处于定直平飞状态，使用升降舵偏角、油门杆位移来平衡复合飞行器系统任务飞行高度、速度。使用副翼输入来平衡复合机体系统的滚转力矩同时平衡系统内部的相对滚转力矩，方向舵来配平复合飞行器系统的偏航力矩。按以上策略在算例工况下实现全配平方案如表2所示，保证复合飞行器整体动力学方程与内部相对转动动力学方程中所有初始状态导数（˙xi除外）小于10－8，则检查的所有配置都满足全配平要求。仿真开始后飞机控制输入在整体运动方程暂不调整，始终保持系统初始稳定状态的输入，仿真时间10 s，仿真步长间隔0.01 s，仿真结果如图8所示。

表2 复合飞行配平方案Table 2 Composite flight trimming strategy

由准配平方案下复合飞行动力学响应仿真结果说明，复合飞行器系统动力学仿真在4.58 s终止，原因是相对滚转角达到180°，即2架复合飞机通过翼尖铰接机构最终在4.58 s以相对滚转运动形式相撞。不同于单机飞行动力学响应，双机复合飞行会存在相对滚转力矩造成内部相对滚转运动，需要通过偏转副翼、升降舵和方向舵实现系统内部动力学方程中的相对滚转力矩。差动双机升降舵会造成多余的俯仰力矩，使得复合飞行器配平问题更加复杂，差动方向舵偏角会造成内部相对偏航力矩，从而对飞行器的翼尖铰接结构强度和机构对接精度产生更高要求。因此，全配平方案主要使用差动副翼舵面偏角来配平内部相对滚转力矩，其中造成运动发散的系统内部相对滚转力矩（l2－l1）可以由面向复合运动特征参数和舵面部件控制参数的气动力数据库得到。仿真结果如图8所示，马赫数、迎角、俯仰角和相对滚转角的发散速度均小于准配平方案仿真结果，同时包括马赫数、迎角、俯仰角和相对滚转角在内的整体和系统内部运动零阶参数在仿真时间10 s内的变化幅值均小于0.3%。

图8 翼尖铰接复合飞行器动力学响应仿真结果Fig.8 Dynamic response simulation results of wingtip-jointed composite aircraft

3.3 本征动力学特性分析

在全配平方案状态下，马赫数、迎角、俯仰角和相对滚转角在内的整体和系统内部运动零阶参数在仿真时间10 s内的数值趋于稳定，以此为基准运动状态，使用小扰动假设［18］下开展全量非解耦的线性化工作，具体包括飞行器系统力和力矩的线化和飞行力学方程组的线化，再将线性化结果整理成状态空间方程有

求解飞行器系统运动状态方程的系统矩阵A的特征值和特征向量可得双机复合运动模态，通过特征值和特征向量的分析，并与单机运动模态对照可得如表3所示的双机复合运动模态分布规律。相比单机运动模态，双机复合运动模态多出了2个发散的相对滚转模态，如表3中命名的“复合运动模态1”和“复合运动模态2”。双机飞行时常规运动模态和单机飞行时的模态特性变化规律可以根据各模态特征值幅值和相位值分析得出为比较双机复合飞行整体系统模态和单机飞行时的变化规律，给出了各运动模态参数指标如表3所示，具体为ωn、ε、ω、T、t1／2、N1／2、t2等特性指标的数值变化；为便于分析2个发散模态的运动特点，在这里分别给出5个模态的特征向量：复合运动模态1的特征矢量：［1.67×10－80.00252.19×10－80.0001－7.88×10－11－1.03×10－10－2.55×10－91.16×10－70.99990.0019］T；复合运动模态2的特征矢量：［0.0024－0.00090.00010.00022.02×10－91.19×10－10－5.6×10－58.42×10－50.15530.9878］T，其中复合运动模态1的特征矢量：忽略数量在10－4以下的参数数值，则相关参数之比α∶θ∶φ12∶˙φ12＝0.0025∶0.0001∶1∶0.0019。复合运动模态2的特征矢量：忽略数量在10－4以下的参数数值，则相关参数之比为V∶α∶q∶θ∶φ12∶˙φ12＝0.0024∶－0.0009∶0.0001∶0.0002∶0.1553∶1。

结合表3、表4中复合运动模态1和复合运动模态2的特征值和模态特性指标参数，可以分析得出复合运动模态1是相对滚转角缓慢发散运动主导的同时伴有微小迎角、俯仰角速度和俯仰角发散的模态，倍幅时间t2为4.813 s；复合运动模态2是相对滚转角速度快速发散运动主导的，同时伴有缓慢的相对滚转角与微小的速度、迎角、俯仰角速度和俯仰角发散运动的模态，倍幅时间t2为0.407 s。复合飞行时各模态指标参数也存在显著变化：长周期的无阻尼自振频率ωn从0.050 rad／s增加到0.054 rad／s，阻尼比ε从0.154减小到0.080，振荡频率从0.0497 rad／s增加到0.0534 rad／s，周期T从126.42 s减小到117.66 s，半幅时间t1／2从90 s提高到161.16 s，半衰期内振荡次数N1／2从0.72次增加到1.37次；纵向模态1的无阻尼自振频率ωn从2.3 rad／s减小到2.27 rad／s，半幅时间t1／2从0.301 s增加到0.305 s；纵向模态2的无阻尼自振频率ωn从0.0552 rad／s减小到0.0486 rad／s，半幅时间t1／2从12.53 s增加到14.26 s；滚转模态的无阻尼自振频率ωn从4.169 rad／s增加到6.658 rad／s，半幅时间t1／2从0.1662 s减小到0.1041 s；螺旋模态的无阻尼自振频率ωn从0.022 rad／s减小到0.008 rad／s，半幅时间 t1／2从31.5 s增加到86.625 s；荷兰滚模态的无阻尼自振频率ωn从0.0169 rad／s增加到0.163 rad／s，阻尼比ε从0.722减小到0.564，振荡频率从0.0118 rad／s增加到0.1346 rad／s，周期T从532.5 s减小到46.7 s，半幅时间t1／2从56.8 s减小到7.541 s，半衰期内振荡次数N1／2从0.1067次增加到0.1615次。

表3 运动模态分析Table 3 Motion modal analysis

表4 模态参数分析Table 4 Modal parameter analysis

由于发散的复合合运动模态2的倍幅时间很短，对飞行稳定性影响很大，飞行器在无控状态下无法稳定飞行，这与传统构型飞行器有很大差别，因此为了保障复合飞行器系统能够长期持续地稳定飞行，需要采用合理的控制策略来保证飞行的稳定性。

4 结 论

1）根据推导得到的7自由度复合机体动力学方程全量模型和建模完成的气动力数据库的动力学仿真结果证明类似单机飞行配平方案的准配平方案无法使复合飞行器系统长期保持飞行稳定，而考虑整体和内部动力学方程组和运动学方程组平衡的全配平方案下复合系统的运动参数在仿真时间10 s内始终趋于稳定，各参数变化幅度均小于5%。

2）在全配平方案下，使用航空飞行器动力学原理的小扰动假设，对复合飞行器动力学方程组和运动学方程组开展线性化工作，整理得到复合系统状态空间方程和系统本征矩阵，分析得到了系统本征矩阵的特征值和特征向量，得到了双机复合飞行模态分布规律，并研究分析了单机飞行模态分布的指标参数数值的变化。

3）研究复合飞行过程出现的2个发散的复杂运动模态，分别由相对滚转角和相对滚转角速度主导，因此需要设计合理的控制方法来保持复合飞行的稳定性。