类X-51A飞行器纵向机动数值虚拟飞行仿真

王胜，王强，林博希，阎超，*

（1.北京航空航天大学 航空科学与工程学院，北京100083； 2.中国航天空气动力技术研究院，北京100074）

吸气式高超声速飞行器由机体和超燃冲压发动机一体化组成，其构型往往较为复杂［1-2］。飞行器高速飞行时，经常需要快速机动到一个较大的迎角，以实现快速的爬升。在迎角快速拉升的过程中，会产生大范围的流动分离，使得飞行器的气动特性出现强烈的非定常特性，气动力迟滞效应明显，这对控制系统的设计带来较大挑战。传统基于数据库或气动力模型的飞行仿真，割裂了气动、运动和控制之间的耦合关系，不能准确描述飞行器机动过程中复杂的气动特性和运动规律，以此设计出的控制律可能难以取得令人满意的效果。因此，需要发展一种更为先进的方法以精确模拟吸气式高超声速飞行器的快速机动过程。

近年来，随着计算机技术的快速发展，数值虚拟飞行技术成为模拟飞行器机动飞行的一种新选择［3］。数值虚拟飞行技术是一种将计算流体力学（CFD）、刚体动力学（RBD）和飞行控制系统（FCS）耦合在一起的高保真计算方法，该方法考虑了流场的非定常特性以及气动特性和运动特性的耦合效应，可以更为精确地获得飞行器的闭环响应特性，受到越来越多的关注［4-10］。

数值虚拟飞行的关键技术主要有3点：一是耦合求解CFD/RBD方程，软件需要具备处理网格运动和变形的能力；二是姿态控制律的设计；三是控制系统与CFD系统的耦合求解方法。这3个关键技术中，耦合求解CFD/RBD方程、处理网格运动和变形的技术在过去几十年中得到了长足的发展，目前已经比较成熟［11-13］。姿态控制律的设计在导航控制领域也研究较多，方法相对成熟［14-15］。而控制系统与CFD系统的耦合求解方法在国内外研究比较少，已有的关于虚拟飞行技术的研究大都采用简化的控制律设计控制器，并将设计好的控制器以CFD 代码的形式嵌入到CFD模块中，达到控制目的［16］。这种做法能够实现运动/流动/控制三者之间的耦合仿真，从而获得比传统方法更精确的闭环响应特性。但是也存在一些不足，比如在CFD软件中实现控制功能编程比较复杂；每次调试更改控制器都要重新编译CFD软件，过程繁琐；通常只能采用简化的控制律函数设计控制器，只能实现较为简单的控制功能等。

为了降低控制系统与CFD系统的耦合难度，本文基于现代软件分布式、模块化的发展趋势，使用在航空航天器导航与姿态控制等领域得到广泛应用的商业软件Simulink实现姿态控制功能，耦合自研CFD软件MICFD［17-18］，建立了Simulink/MICFD数值虚拟飞行仿真平台。利用该仿真平台，对类X-51A外形吸气式高超声速飞行器进行了纵向机动闭环数值仿真，通过与工程仿真结果对比，研究了运动和气动耦合情况下非定常效应对飞行器控制响应的影响。此外，利用该平台，还进行了一些简单的应用，研究了纵向机动过程中舵回路时间常数对控制性能的影响。

1 数值虚拟飞行仿真平台的实现

为了建立Simulink/MICFD数值虚拟飞行仿真平台，首先要解决2个软件之间的数据传输问题，在保证Simulink和MICFD同步运行的前提下，实现两者之间稳定高效的数据传输。

Simulink可以利用MATLAB的各种命令和库函数，实现复杂的控制仿真任务，再利用MATLAB的RTW 模块可以将生成的Simulink模型转变为可以直接运行的C++程序，这就为Simulink与其他应用程序的耦合提供了技术途径。

远程过程调用（Remote Procedure Call，RPC）［19］是一种通过网络从远程计算机程序上请求服务，而无需了解底层网络技术的协议。RPC采用客户机/服务器模式，实现进程间的同步机制，为用户提供请求/应答的通信方式。发出请求的程序称为客户机，提供服务的程序称为服务器。gRPC是由Google公司主导开发的一款语言中立、平台中立、开源的RPC框架，支持C、C++、Python、java等多种语言版本［20］。gRPC客户端和服务端可以在多种环境中运行和交互，由于其跨平台、跨语言的特点，因此能够实现不同软件系统之间的通信，如图1所示。

图1 不同环境下gRPC实现过程Fig.1 Schematic diagram of gRPC implementation process in different environments

根据上述分析，可以基于gRPC网络通信原理，在Simulink和MICFD软件中分别建立服务器和客户端，通过远程调用实现2个软件之间高效稳定的数据通信，并在此基础上建立Simulink/MICFD数值虚拟飞行仿真平台。如图2所示，其具体实现步骤如下：

步骤1启动Simulink，根据需要设计控制器，建立仿真程序，并生成可独立运行的C++程序。

步骤2启动MICFD，读入网格和初始计算条件，进行定常计算，得到初始的计算流场。

步骤3进行非定常求解，计算当前时刻的力和力矩。

步骤4RBD模块通过求解飞行力学方程，计算飞行器的姿态，将当前的飞行姿态与期望的姿态进行比较，得到姿态偏差，并将偏差信息传递到gRPC客户端。

步骤5gRPC客户端通过网络通信向gRPC服务器发送请求，将偏差信息传递到gRPC服务器。

步骤6gRPC服务器接收到请求之后，调用步骤1生成的C++project，根据设计的控制器，求出每个控制面的偏转指令，然后将舵偏指令通过网络通信传递到gRPC客户端。

步骤7CFD模块获得gRPC客户端得到的舵偏指令，并根据当前时刻的姿态进行网格移动和网格重叠，准备下一物理时间步的非定常计算。

步骤8重复步骤3～步骤7，直到物理时间推进结束。

图2 Simulink/MICFD数值虚拟飞行仿真平台流程图Fig.2 Flowchart of Simulink/M ICFD numerical virtual flight simulation platform

2 数值虚拟飞行计算方法

2.1 流动控制方程与求解方法

流动控制方程为三维可压缩Navier-Stokes方程，无量纲的守恒形式为

式中：Q为守恒变量；F、G和H 为对流通量；FV、GV和HV为黏性通量；ξ、η和ζ为3个贴体坐标系方向；t为无量纲时间；Re∞为自由来流雷诺数。

流场求解采用基于结构重叠网格的有限体积法。空间离散上，无黏通量使用Roe通量差分分裂格式求解，黏性通量使用二阶中心差分格式进行离散；时间推进上，采用双时间步LU-SGS（Lower-Upper Symmetric Gauss-Seidel）方法，湍流模型采用两方程的SST模型假设。在涉及到动态网格计算时，由于离散后的方程中含有网格体积对时间的导数项，为避免网格变化引入的额外误差，满足几何守恒律十分重要［21］。由于本文使用的是刚性重叠网格技术，在计算过程中不涉及网格体积的变化，并且本文采用的空间离散方式与坐标变换格式相匹配，极大地缓解了几何守恒律的影响，因此在本文的研究中未对几何守恒律做单独处理。

2.2 行力学方程求解

飞行器的运动可以分解为两部分：质心的平动和绕质心的转动。惯性坐标系下的质心平动动力学方程组可以表示为

式中：m为飞行器的质量；V为速度矢量；Fa为作用在飞行器上的空气动力矢量；Fe为外力矢量，如推力等；Fg为重力矢量。体轴坐标系下，绕质心转动的动力学方程可以表示为

式中：Hc为飞行器相对于质心的动量矩矢量；ω为角速度矢量，其分量为（ωx，ωy，ωz）；M 为力矩矢量。体轴坐标系下，绕质心转动的运动学方程可以表示为

式中：（φ，θ，γ）为飞行器3个方向的姿态角。按式（4）方程可以用龙格-库塔方法进行积分求解，并采用双欧法克服方程的奇异性［22］。

2.3 姿态控制系统设计

姿态控制系统的功能是根据指令自动调整飞行器的飞行姿态，使其能够保持稳定或者按照预定指令绕其质心转动［23］。本文所研究的纵向姿态控制的实现如图3所示。图中：αc为目标迎角；α为当前时刻迎角；e为姿态角误差；δe为控制器输出的舵偏指令；δc为舵机输出的舵偏角。图3（a）为工程上常用的基于静态数据库的方法，该方法根据已有的静态数据库通过插值的方法获得不同飞行姿态下的气动力数据，以此气动力数据进行飞行力学求解以获得新的飞行姿态。图3（b）为本文所采用的耦合CFD和控制系统的一体化模拟方法，该方法通过求解飞行器运动过程中的实时流场以获取更为精确的非定常气动力来代替工程方法中插值得到的气动力，由于考虑了流场的非定常特性，其计算结果更为真实、可靠。

理论上，在数值虚拟飞行技术中可以使用任意的控制器以满足不同的控制需求，复杂控制器的实现以及性能考察不是本文的研究重点。因此，简单起见，在本文的研究中，控制器采用常用的PID控制器，其输入输出之间的关系可以表示为

图3 纵向姿态控制Fig.3 Longitudinal attitude control

式中：e（t）为t时刻的误差；Kp、Ki和Kd分别为比例系数、积分系数和微分系数。

舵回路采用基于位置反馈的硬反馈式回路，其近似模型可用一阶惯性环节表示［24］。

式中：K为反馈增益，一般取值为1；δc（s）为δc（t）在变换后空间的表述，δc（t）为t时刻控制舵偏指令值；Td为舵回路的时间常数，反应了舵机的迟滞特性，其值的大小受舵机的功率等自身因素影响。

3 算例验证

为了考核MICFD软件中CFD/RBD耦合计算能力，选取三维机翼挂载分离模型，对程序进行验证。该模型由3部分组成，即机翼、挂架以及外挂物，其几何外形具体参数详见文献［25］。计算条件为：Ma=0.95、H=7.92 km、α=0°。

图4和图5分别为分离过程中挂载物的气动力、力矩以及质心位移和角位移随时间的变化曲线。可以看出，各曲线的计算结果与实验结果［25］吻合良好，验证了本文CFD/RBD耦合计算方法的准确性，表明软件具备良好的非定常多体相对运动数值模拟能力。

图4 挂载物气动力系数、气动力矩系数随时间变化曲线Fig.4 Time history of aerodynam ic coefficients and aerodynamic moment coefficients for store

图5 挂载物线位移、角位移随时间变化曲线Fig.5 Time history of linear and angular displacement for store

4 类X-51A飞行器纵向机动过程闭环模拟

4.1 计算模型和网格

本文首先基于X-51A飞行器的相关资料［1，26］，获得其几何外形的主要参数，采用反向建模技术建立了一种与X-51A外形相似的高超声速飞行器模型，如图6所示。本文主要研究飞行器纵向机动过程，模型左右对称且不考虑侧滑，为节省计算量采用半模计算。网格生成时，飞行器本体网格和控制舵面网格分别独立生成，并采用动态网格重叠的方式实现两者的相对运动。网格总量约1 600万。机动过程中不同时刻控制舵与机身的动态重叠边界示意图如图7所示。

图6 计算模型及对称面网格Fig.6 Computationalmodel and symmetry plane grids

图7 不同时刻控制舵与机身动态重叠边界示意图Fig.7 Schematic diagram of dynamic overlapping boundary of control rudder and airframe at differentmoments

4.2 静态气动特性

选取的计算状态为典型的进气道起动状态。来流马赫数为4.8，来流单位雷诺数为1.2×107。力矩计算时，为方便起见，参考长度取Lref=1m，参考面积取Sref=1m2。坐标原点设置在飞行器前缘纵向对称面处。力矩参考点（x，y，z）=（1.8，-0.5，0）m。

图8显示了不同迎角和舵偏角下飞行器的俯仰力矩，其中点画线代表的是全机俯仰力矩，实线是控制舵产生的俯仰力矩。可以看出，在所计算的迎角范围内，除20°舵偏角之外，其余舵偏下俯仰力矩的线性度比较好。而全机俯仰力矩具有较明显的非线性。图9给出了各迎角和舵偏角下除控制舵以外的机身部分的俯仰力矩。可以看出，与传统的轴对称外形或者升力式外形飞行器不同，该外形机身部分的俯仰力矩呈现明显的非线性，说明由于进气道和内流道的存在，迎角的改变对气动特性的影响更加复杂。

图8 全机及控制舵俯仰力矩Fig.8 Pitching moment of whole aircraft and control rudder

图9 除控制舵以外机身的俯仰力矩Fig.9 Pitching moment of airframe except control rudder

4.3 纵向机动过程模拟

本节对建立的类X-51A模型进行纵向机动仿真。计算的来流条件与静态计算时相同，初始的来流迎角为0°，控制指令为：以阶跃响应的方式机动到10°迎角，并保持稳定。分别采用图3（a）介绍的工程方法和图3（b）介绍的耦合方法进行计算。由4.2节计算可知，初始状态下控制舵面的配平偏转角为-11°，闭环模拟时以该舵偏角作为初始舵偏，即假定开始机动时，飞行器处于平衡状态（后文中的舵偏角均指在此基础上的相对舵偏角）。

通过试凑，选取一组参数：Kp=2，Ki=1.5，Kd=0.6作为控制器的增益。舵回路中，反馈增益K=1，舵回路时间常数Td=0.1。采用工程方法计算时，飞行器不同迎角和舵偏角下的气动力由静态数据库通过插值得到。非定常计算时，物理时间步长设为1ms，子迭代残差指标为0.01，同时限定最大子迭代步数为50步，子迭代CFL数为3.0。

图10给出了采用耦合方法模拟时不同时刻飞行器周围的流场结构，包括壁面压力分布、轴向不同截面的马赫数等值线以及控制舵附近的流线等。可以看出，相比于图10（a）的初始定常状态，飞行器在拉起过程中迎风面和背风面出现明显差异，控制舵附近流线偏转，出现了较强的非定常特征。

图10 机动过程中典型时刻的流场结构（T d=0.1）Fig.10 Typical-moment flow field structure during maneuvering process（T d=0.1）

图11～图13分别表示了2种方法得到的迎角、舵偏角以及俯仰力矩的响应过程。可知，2种方法均可以实现控制指令要求的机动过程，最终的配平舵偏角均保持在-5.6°附近，但是两者之间也存在一些差异。从迎角响应过程来看，耦合方法计算得到的超调量（16%）略高于工程方法（13.4%），耦合方法计算得到的调节时间为2.01 s，远高于工程方法得到的1.48 s，说明工程算法给出的结果可能低估了控制系统的迟滞特性。从舵偏角响应过程来看，尽管最终的配平舵偏相同，但是两者在峰值位置存在较大差异，这正反映了传统的基于静态气动力数据库的飞行仿真在模拟快速机动过程中的缺陷，体现了流场非定常效应对控制过程的影响，这一点通过图13中俯仰力矩在峰值处的差异也可以看出。因此，对于类似X-51A、外形复杂的高超声速飞行器，由于非定常效应明显，很有必要采用基于耦合方法的数值虚拟飞行技术来研究飞行器的机动过程、评估飞行控制律。

图11 迎角响应过程Fig.11 Response process of angle of attack

图12 舵偏角响应过程Fig.12 Response process of rudder deflection angle

图13 俯仰力矩时间历程Fig.13 Time history of pitching moment

4.4 舵回路时间常数的影响

舵回路系统也称伺服系统，是飞行控制系统中的重要组成部分。舵机响应的迟滞特性对飞行器的控制系统有着显著影响。因此，对控制系统在舵机不同响应特性下的控制性能进行评估和分析，为控制律的设计提供参考具有重要的意义。在本节中，通过改变舵回路系统的时间常数Td进行数值虚拟飞行仿真，以评估其对飞行操纵过程的影响。其他仿真条件与4.3节一致。

图14～图16给出了迎角、舵偏角以及俯仰力矩在不同舵回路时间常数下的响应过程。可以看出，不考虑舵回路（即Td=0）时，舵偏操纵是理想的，即δc=δe。控制舵一开始就存在一个较大的偏转角（约-20°）；而在考虑舵回路时，由于舵机的迟滞特性，舵偏角是从0°开始逐渐增大的，这也导致2种情况下飞行器的机动过程出现较大差异。表1列出了不同时间常数下控制系统的性能指标，包括系统的延迟时间（td）、上升时间（tr）、峰值时间（tp）、调节时间（ts）以及超调量（σ%）。可以看出，随着Td的增加，系统的调节时间和超调量显著增加，动态性能变差。值得注意的是，当Td=0.2时，超调量达到31.6%，并且纵向姿态在数值模拟的时间历程中未能收敛到指定的稳定姿态。可以预见，继续增大Td可能会使系统难以收敛甚至发散。

从以上分析可以看出，舵回路的时间常数对控制过程有重要影响。在满足舵机功率限制的前提下，通过设计适当的舵回路来减小Td，可以降低超调量，提高响应速度，改善系统的动态性能。

图14 不同时间常数下迎角响应过程Fig.14 Response process of angle of attack under different time constants

图15 不同时间常数下舵偏角响应过程Fig.15 Response process of rudder deflection angle under different time constants

图16 不同时间常数下俯仰力矩时间历程Fig.16 Time history of pitching moment under different time constants

表1 不同T d 下控制系统的性能指标Tab le 1 Perform ance indexes of control system under d ifferent T d

5 结 论

本文从现代软件分布式、模块化的发展趋势出发，基于 gRPC 网络通信建立了 Simulink/MICFD数值虚拟飞行仿真平台，有效降低了控制模块和CFD模块的耦合难度。利用该仿真平台，对类X-51A外形的吸气式高超声速飞行器进行了纵向机动闭环数值仿真。得到如下结论：

1）本文采用的远程过程调用的方法在形式上完全分割了控制和CFD这2个模块，但是在逻辑上又将2者紧密耦合在一起。相比代码级耦合，基于此方法建立的数值虚拟飞行仿真平台能够充分发挥控制和CFD软件各自的优点，显著降低了多学科耦合的难度。

2）仿真结果表明，该数值平台具备针对复杂飞行力学行为和控制响应特性的先进数值模拟能力，能够实现对飞行器有控机动过程的精细化模拟。

3）对于类X-51A外形的吸气式高超声速飞行器，在纵向拉起时，工程算法给出的结果可能不能完全反映非定常效应的影响。有必要采用基于耦合方法的数值虚拟飞行技术来研究飞行器的机动过程、评估飞行控制律。

4）舵回路的时间常数对控制系统的性能有重要影响，减小时间常数可以降低超调量，提高响应速度，改善系统的动态性能。