0 引 言
飞行稳定性与控制都是必须集智攻关的核心技术问题。当前对飞行器的飞行品质评估和飞行安全分析通常采用传统的飞行仿真，即通过计算或实验获得飞行器在不同攻角、侧滑角和舵面偏转等不同状态的定常气动特性和动导数[1, 2, 3]，构建数据库和数学模型来分析飞行器的飞行性能。
随着CFD数值模拟技术的发展以及硬件水平的快速提高，当前CFD已经能够产生用于飞行器载荷计算所需的气动特性输入，因此通过与控制飞行器运动的动力学方程相耦合，可以直接模拟机动状态下飞行器的气动响应。耦合刚体运动的CFD已经研究多年，但主要局限于飞行器动导数的计算、二维验证算例以及简单外形的三维算例[4, 5, 6]。在最近十年，对于复杂几何外形的计算流体力学/飞行动力学耦合的研究开始增多，已经不再满足于无控开环状态下飞行器运动特性模拟和分析，结合操纵面响应的闭环飞行器控制机动行为的模拟也逐渐吸引了越来越多研究者的关注[7, 8, 9]。
飞行器数值虚拟飞行(飞行器连续偏舵闭环控制的流动/运动/控制一体化数值仿真)的难点主要有两个方面：1) 在计算流体力学CFD中，操纵面每时每刻都处在不同的位置，因此需要数值模拟软件具备处理网格运动和变形的能力；2) 舵面控制律的实现以及与数值模拟软件的结合。
总的来说，对飞行器闭环机动的数值模拟研究是一个尚未深入的领域，有许多尚待解决的问题。尽管由于费用等因素这种研究当前不大可能普遍用于飞行力学的常规研究，但对调查和理解飞行器飞行中的潜在问题能提供重要的帮助。
本文采用基于非结构重叠网格的气动舵连续操纵，耦合刚体运动方程，非定常模拟高超声速飞行器的纵向机动过程，其目标是发展能够用于研究飞行器复杂飞行力学行为和控制响应的先进数值模拟能力和技术。 1 数值虚拟飞行的关键技术 1.1 非结构重叠网格技术
网格运动和变形的处理方法通常分为两大类，分别是重叠网格[10, 11, 12]和网格变形[13, 14]。对刚性大位移运动，如果采用网格变形的方法，大位移时势必要对网格进行重构，增加了解决问题的难度，且重构不可避免丧失部分精度。而重叠网格技术中各个独立网格的拓扑结构不会随着物体的相对运动而改变，因此基于重叠网格技术的数值模拟方法特别适合于处理刚性物体之间的相对运动。最早的重叠网格是基于结构网格的，目的是用于航天飞机之类复杂外形的网格生成。基于结构网格的重叠网格方法虽然部分改善了网格生成的难度，但仍然需要使用者大量的人工干预和较为丰富的使用经验，极大地限制并抵消了重叠网格所带来的优点。非结构重叠网格技术在最近10年得到快速发展，它结合了非结构网格和重叠网格的优点，得到了越来越多研究者的重视。
本文采用非结构重叠网格方法处理操纵面的偏转。非结构重叠网格技术主要包括以下四个方面：(1) 洞边界的构造；(2) 边缘单元的确定；(3) 贡献单元的查找；(4) 网格之间的插值。采用高效的ADT数据结构快速地实现鲁棒性的洞边界的构造，准一维的快速搜索方法得到贡献单元，二阶精度的插值计算以保证和流场解的空间精度相同。
图 1为使用非结构重叠网格方法后弹体网格在舵面上形成的洞边界，图 2为舵网格在弹体上形成的洞边界。弹体和舵面各自的独立网格由笛卡尔网格方法生成[15]。从图 2中可见，在舵缝位置处弹体网格和控制舵网格都进行了加密，避免了孤立单元的出现，同时保证边缘单元能够得到足够的插值模板。
 |
| 图 1 弹体网格在舵上形成的洞边界Fig. 1 Hole boundary of body grid around the elevator |
 |
| 图 2 舵网格在弹体上形成的洞边界Fig. 2 Hole boundary of elevator grid around the body |
图 3给出了高速飞行器(Ma=6.0，H=30km)采用重叠网格和单块网格计算得到的气动特性随攻角变化的曲线，从图 3中可见，重叠网格获得的气动特性和单块网格的结果吻合得非常好。重叠网格和单一网格在静态计算的比较结果表明，当前所发展的重叠网格技术在精度上和单一网格的计算精度相当。
 |
| 图 3 重叠网格和单块网格的气动特性计算结果比较Fig. 3 Comparison of aerodynamics between overset grid and single grid |
本文推广应用了亚迭代思想，采用全局同步亚迭代的求解策略来获得高的耦合时间精度(图 4)。在构造流体力学方程和刚体动力学方程亚迭代形式的基础上，同步推进流体力学方程和刚体动力学方程，在每一子迭代的过程中，及时更新刚体动力学方程中的气动力(矩)源项且及时更新流体动力学方程中的网格运动速度与飞行姿态；在亚迭代过程中实现子系统间的信息互换；当亚迭代残差收敛于0时，可以发现无论流体动力学方程还是刚体动力学方程都同步达到二阶精度，即耦合推进的时间精度为二阶；该方法避免了通常解耦计算的一阶时间滞后。通过全局亚迭代策略可以获得物理时间精确的非定常CFD/RBD耦合方法。本质上讲，全局亚迭代耦合求解方法是将源于非定常流体力学计算的亚迭代思想推广应用到多系统动力学耦合问题。
 |
| 图 4 全局亚迭代耦合方法的流程图Fig. 4 Flow chart of global sub-iteration coupling method |
飞行控制律设计及飞行仿真，一方面为背景飞行器控制飞行数值仿真提供舵面输入或控制规律，另一方面控制律设计后的六自由度飞行仿真也为流动/运动/控制一体化数值虚拟飞行仿真提供比较对象。
针对静不稳定布局开展稳定控制飞行，研究设计一个简单的纵向控制律，实现变攻角调姿的同时进行静不稳定控制。控制律设计的数据源为定常数值计算结果和动导数计算结果，非计算点的结果通过线性插值得到。
经过稳定性分析，选择合适的反馈增益，考虑舵面限制和仿真任务，在理想舵面偏转情况下，设计仿真模型的控制律与控制指令为(各变量单位均为度)：
其中指令变量：
在上述控制系统和控制操纵下，仿真10s所得结果如图 5所示，攻角可以实现不同攻角位置的平衡状态。当攻角从第一个平衡位置改变到第二个平衡位置，舵偏首先减小(在这里由于舵偏限制而保持0°)，飞行器抬头，攻角和俯仰角速度增大，到达某个位置舵偏角也开始逐渐增大，飞行器在控制系统作用下逐渐平衡；从第二个平衡状态到第三个平衡状态也有相同的变化过程。
 |
| 图 5 控制律仿真Fig. 5 Flight simulation according to the designed control law |
为了检验飞行器部件运动条件下非结构重叠网格适应性，在配平状态下研究了周期性俯仰舵偏对高超声速飞行器气动特性的影响。计算网格和计算条件与前述定常计算所用的非结构重叠网格相同。飞行器俯仰舵偏δm0为1°时对应的配平攻角为12.21°；在此固定姿态下，舵面围绕平衡舵偏进行谐波形式的振荡：
式中,A0为舵面的振幅，f是振动的频率。在本文研究中，选取的振动频率有2个，分别为20Hz和40Hz。计算中每一个振荡周期物理时间划分2 000步以满足非定常计算的要求。从图 6的物理时间步内全弹气动力系数的内迭代收敛情况看，每一物理时间步的内迭代收敛都得到保证。
 |
| 图 6 内迭代收敛情况Fig. 6 Convergence of sub-iteration |
从图 7和图 8可发现，无论是全弹气动特性还是舵面气动特性，都对俯仰舵偏呈现良好的线性。舵面振荡频率较高时，气动特性的迟滞效应略明显一些，但仍与定常状态获得的结果非常一致，表明实际的气动迟滞效应并不显著。从图 7中还可见，在此条件下，俯仰舵效不随偏转频率的不同而发生改变，而且和定常结果一致，表明从定常计算结果得到的舵效是可靠的。
 |
| 图 7 周期舵偏运动过程的全机质心力矩系数Fig. 7 Vehicle moment coefficient around gravity during periodic elevator deflection |
 |
| 图 8 周期舵偏运动过程的舵面质心力矩系数Fig. 8 Elevator moment coefficient around gravity during periodic elevator deflection |
图 9给出了俯仰舵偏在0°~2°振幅区间和平衡位置时中间截面的压力等值线云图。随着舵偏的增大，舵下表面的压力显著增加。由于离弹体较远，舵面后部压力的增加更明显。
 |
| 图 9 俯仰舵偏不同位置时的中间截面压力云图Fig. 9 Pressure contour of medial section during periodic elevator deflection |
针对高超声速布局，采用上述的简单控制律进行纵向闭环数值虚拟飞行模拟。计算网格与前述定常计算所用的非结构重叠网格相同。计算初始条件是H=30km，Ma∞=6.0，俯仰平衡舵偏δm=1°，飞行器在该平衡舵偏下的配平攻角为12.21°。采用1.3节的控制律和控制指令。
图 10给出了飞行器在带控制飞行过程中攻角随时间的变化、舵偏角随时间的变化情况，并和仿真结果进行了比较。从图中可见，带控制闭环数值虚拟飞行的结果和采用控制律仿真结果的趋势吻合得很好。在CFD/RBD模拟过程 中加入设计的控制律，正确模拟了飞行器从较低攻角拉起然后在新攻角配平的过程。
 |
| 图 10 飞行器姿态和舵偏随时间的变化Fig. 10 Variations of state and elevator deflection |
从图 10中可见，为了拉起飞行器，俯仰舵偏从初始的1°回到0°，维持一段时间后再打正舵偏。初始俯仰舵偏变小的原因是飞行器是纵向静不稳定的。因此要使飞行器抬头必须减小俯仰舵偏，在获得足够的向上速度和攻角后再增加俯仰舵偏角最终使得飞行器在新的攻角配平。
图 10中数值虚拟飞行的结果与仿真结果相比，攻角基本吻合，第一段拉起配平过程中的最大舵偏角和最终配平舵偏角稍大。原因主要是由于仿真和数值模拟两者对应的配平舵偏有差异所造成的，因此将仿真分析获得的控制律用于数值模拟必然会导致舵偏角的偏差。造成仿真和数值模拟两者对应的配平舵偏差异的主要原因有：① 用于仿真分析的气动力数据计算时采用的模型和网格与数值飞行时有差别(图 11)，从而导致同一个状态的气动力数据不同。② 仿真时用到的气动力数据较稀，非计算点的数据采用线性插值得到，而数值模拟时的气动力实时计算连续获得。后续将发展更合适的气动模型用于仿真分析。③ 当前的舵面控制律是无舵机函数的控制律，舵偏减小的过程是阶跃性的，在数值模拟中可能会导致计算的偏差。
 |
| 图 11 静态仿真和数值虚拟飞行计算网格的比较Fig. 11 Surface mesh of single grid and overset grid |
图 12给出了弹道和舵偏典型点的飞行器姿态和舵偏位置。流场显示表明两次拉起变攻角过程，下表面压力屡次增加。舵偏先收回再打到14.5°攻角配平所需要的2.8°正舵偏，舵偏第二次收回再打到17.5°攻角配平所需要的3.6°正舵偏。上述流动和舵偏运动与纵向变攻角操纵的简单控制系统设计仿真历程吻合。
 |
| 图 12 典型时刻飞行器流场、姿态和舵偏位置Fig. 12 Flow field,state and elevator deflection at the typical moments |
本文进行了飞行器连续偏舵变攻角闭环控制的流动/运动/控制一体化数值仿真(数值虚拟飞行)。研究获得以下启示：
(1) 依据大型计算系统和先进的非定常计算软件实现流动/运动/控制一体化数值仿真是现实可行的。
(2) 数值虚拟飞行的关键技术主要有：①跨学科非定常一体化耦合模拟技术；②保持时空精度的运动网格技术(非结构重叠网格技术)。
(3) 数值虚拟飞行与依据气动数据模型的理想控制仿真基本吻合，但过程偏差是存在的。
(4) 在复杂飞行包络的典型状态和极限状态，数值虚拟飞行具有更突出的必要和研究价值。
| [1] | Frink N T. Stability and control CFD investigations of a generic 53° swept UCAV configuration[R]. AIAA 2014-2133, 2014. |
| [2] | Yuan Xianxu, Zhang Hanxin, Xie Yufei. The pitching static/dynamic derivatives computation based on CFD methods[J]. Acta Aerodynamica Sinica, 2005, 23(4):458-463. (in Chinese)袁先旭, 张涵信, 谢昱飞. 基于CFD方法的俯仰静、动导数数值计算[J]. 空气动力学学报, 2005, 23(4):458-463. |
| [3] | Schoenenberger M. Supersonic pitch damping predictions of blunt entry vehicles from static CFD solutions[R]. AIAA 2013-0356, 2013. |
| [4] | Zhou Weijiang. Numerical simulation of unsteady transonic flowaround a free oscillating reentry vehicle[J]. Acta Aerodynamica Sinica, 2000, 18(1):46-51. (in Chinese)周伟江. 返回舱自由振动跨声速非定常流场数值模拟[J]. 空气动力学学报, 2000, 18(1):46-51. |
| [5] | Sahu J. Time-accurate numerical prediction of free flight aerodynamics of a finned projectile[R]. AIAA 2005-5817, 2005. |
| [6] | Ronch A D, Vallespin D, Ghoreyshi M, et al. Evaluation of dynamic derivatives using computational fluid dynamics[J]. AIAA Journal, 2012, 50(2):470-484. |
| [7] | Clifton J D, Ratcliff C J, Bodkin D J, et al. Determining the stabilityand control characteristics of high-performance maneuvering aircraft using high-resolution CFD simulation with and without moving control surfaces[R]. AIAA 2013-0972, 2013. |
| [8] | SchüTte A, Einarsson G, Raichle A, et al. Numerical simulation of maneuvering aircraft by aerodynamic, flight-mechanics, and structural-mechanics coupling[J]. Journal of Aircraft, 2009, 46(1):54-64. |
| [9] | Tao Yang, Fan Zhaolin, Wu Jifei. CFD based virtual flight simulation of square cross-section missile with control in longitudinal flight[J]. Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(2):169-176. (in Chinese)陶洋, 范召林, 吴继飞. 基于CFD的方形截面导弹纵向虚拟飞行模拟[J]. 力学学报, 2010, 42(2):169-176. |
| [10] | Noack R W, Boger D A, Kunz R F, et al. Suggar++:An improved general overset grid assembly capability[R]. AIAA 2009-3992, 2009. |
| [11] | Biedron R T, Thomas J L. Recent enhancements to the FUN3D flow solver for moving-mesh applications[R]. AIAA 2009-1360, 2009. |
| [12] | Zhang S J, Zhao X. Computational studies of stage separation with an unstructured chimera grid method[R]. AIAA 2007-5409, 2007. |
| [13] | Biedron R T, Vatsa V N, Atkins H L. Simulation of unsteady flows using an unstructured Navier-Stokes solver on moving and stationary grids[R]. AIAA 2005-5093, 2005. |
| [14] | Sitaraman J, Baeder J D. Field velocity approach and geometric conservation law for unsteady flow simulations[J]. AIAA Journal, 2006, 44(9):2084-2094. |
| [15] | Liu Zhou, Zhou Weijiang. Adaptive viscous cartesian grid generation and application[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(12):2280-2287. (in Chinese)刘周, 周伟江. 适于黏性计算的自适应笛卡儿网格生成及其应用[J]. 航空学报, 2009, 30(12):2280-2287. |