伍强，徐浩军，魏扬，裴彬彬, 3，薛源

1. 空军工程大学 航空工程学院，西安 710038 2. 中国人民解放军 94639部队，南京 211500 3. 西安电子科技大学 电子工程学院, 西安 710071

波音公司统计结果表明，飞机结冰是诱发飞行失控严重事故的三大环境因素之一。结冰会破坏飞机原有的气动外形，结冰后的流动表现出强烈的非定常、非线性和随机性，如边界层转捩位置的改变、大范围分离流动、不同尺度旋涡的产生及脱落等。这种复杂流动使得飞机部件表面的气动载荷处于不规则的变化之中，进而产生异常的气动力和力矩。异常气动特性作用下的飞机，其操纵性、稳定性等飞行力学特性被破坏，飞机的姿态和飞行状态往往会产生不可控的变化，这种变化又反作用于飞机的流场和气动特性，从而改变冰的演变过程。同时，结冰的动态增长又将导致飞机纵横向气动特性的耦合程度加剧、飞行模态改变，甚至舵面操纵困难。现有研究工作通常是把结冰之后的空气动力学和飞行力学解耦，即是基于气动力数据库的飞行动力学仿真，飞行过程中的气动力及力矩系数是由静动态气动数据库插值得到。该方法虽然在气动模型中加入了动导数项，一定程度上考虑了非定常效应，但本质上仍是一种基于线性叠加原理的准定常方法。该方法无法适用于飞行器处于强非线性、非定常流场的情形，且未深入考虑非定常空气动力学和飞行动力学之间的相互耦合作用。

通常研究飞行器气动/运动耦合特性的方法有以下3种：① 虚拟飞行试验；② 风洞自由飞行试验；③ 基于CFD的气动/运动等多学科耦合的一体化数值模拟，即“数值虚拟飞行”技术。前两种方法通过试验手段能较为真实地模拟飞行器机动飞行过程，得到更为准确的动态气动特性，然而试验技术十分复杂，难度较大。随着CFD基础理论和算法的发展，数值虚拟飞行已成为当前CFD研究领域的热点。2007年，美国国防部投资3.6亿美元开展了计算研究和工程采办工具与环境项目(CREATE)。研究人员通过该项目开发了固定翼飞机虚拟模拟工具Kestrel软件，实现了固定翼飞机空气动力学、结构动力学、飞行动力学、发动机推进等多学科的一体化耦合模拟，基于该软件对不考虑舵偏的F-22飞机做快速拉升及失速过程进行了模拟，通过与飞行数据的对比，发现模拟结果准确度较高，英国的Allan等对某标模在给定舵偏规律下的自由俯仰运动进行了模拟。在非定常气动力计算和飞行动力学方程耦合求解的研究方面，国防科技大学的杨小亮对飞行器的多自由度耦合摇滚运动进行了模拟；北京航空航天大学的阎超等对返回舱、带翼导弹的动态特性进行了数值模拟；中国航天空气动力研究院的杨云军等基于数值虚拟飞行技术研究了三角翼在滚转和侧滑两自由度下的耦合运动特征。中国空气动力研究与发展中心的张来平等系统论述了虚拟飞行中气动、运动和控制耦合的数值模拟技术，指出动态网格的生成、高精度湍流模拟方法及耦合求解算法是该技术的难点，并给出了典型的气动/运动/控制耦合计算流程。

目前国内外没有专门针对结冰条件下飞机气动/运动耦合特性展开研究。国内已有的研究主要是对飞行器机动飞行过程的模拟，且大多是对飞行器单自由度或双自由度自激运动的简单模拟。本文研究了机翼前缘结冰条件下飞机气动/运动的耦合特性，介绍了松耦合气动力与飞行动力学耦合求解策略。通过NACA0012翼型强迫俯仰振荡标准算例验证了开发的非定常CFD求解程序及耦合求解算法的有效性和准确性。在此基础上开发了高精度非定常流场求解器与飞行动力学的双向耦合仿真平台，对流场控制方程、飞行动力学方程和运动学方程组在时间尺度上耦合推进求解，能够刻画结冰增长特性-分离流动非定常特性-飞行动力学特性三者的相互关系，对于研究结冰条件下的飞行安全问题奠定了基础，具有重要意义。

1 控制方程及其离散格式

1.1 耦合控制方程组

建立描述飞机运动与非定常流场耦合的控制方程组：

(1)

式中：为Navier-Stokes方程中流体的守恒变量;、和分别为直角坐标系下、、这3个方向上的无黏通量；、、分别对应3个方向黏性通量；为动力学方程状态变量，包含(飞行速度)、(迎角)、(侧滑角)、(滚转角速度)、(俯仰角速度)和(偏航角速度)；为运动学方程的状态变量，包含(滚转角)、(俯仰角)、(偏航角)和位移、、。

从式(1)可以看出，非定常Navier-Stokes方程与飞行动力学方程组互相耦合，相互影响，互为输入。具体来说就是：由求解非定常Navier-Stokes方程计算出飞行器的气动力及力矩系数，传递给飞行动力学方程，再由飞行动力学方程组解算得到姿态、位移等物理量传递给非定常Navier-Stokes方程组，作为流场初始状态的输入条件，然后更新网格进行解算。

在耦合计算的过程中，飞机所受气动力及气动力矩系数可由Navier-Stokes方程计算得到。根据有限体积法的思想，飞行器受到的气动力/力矩为物面所有微元所受气动力/力矩的总和，而每个微元所受气动力为该处应力张量和外法向矢量的点积，所受气动力矩为矢径与气动力的叉积，具体可表示为

(2)

式中：表示一个微元；为外法向单位向量；为表面切应力；为表面压应力；(，，)为飞机的气动中心位置。根据式(2)可计算得到飞行器的气动力及力矩系数，如式(3)所示：

(3)

1.2 离散格式

1.2.1 非定常Navier-Stokes方程

采用Jameson有限体积法对三维可压缩非定常Navier-Stokes方程进行空间离散，对无黏项(对流和压力项)的空间离散采用Roe-FDS通量差分分裂格式，对网格界面通量采用三阶MUSCL格式，以提高空间分辨率，对黏性项(剪切应力和热传导项)的空间离散采用二阶中心格式。采用隐式LU-SGS双时间步法进行时间推进，该方法在经典隐式LU-SGS格式中引入伪时间导数项，借助伪时间方向的子迭代技术(-TS 迭代)，使得时间离散精度达到二阶。

1.2.2 飞行动力学方程

飞机六自由度动力学和运动学方程组是典型的非线性常微分方程组，最常用的是显式的四阶Runge-Kutta数值解法。其计算公式可描述为

(4)

式中：Δ为时间步长；，，，表示为

(5)

参考文献[17]中的方法，对飞行器刚体六自由度动力学方程采用与Navier-Stokes方程相同的时间离散格式，根据线性多步法松弛法迭代求解：

(6)

式中:为六自由度动力学/运动学方程中的变量、；为对应的六自由度动力学/运动学变量的通量(导数)；为松弛因子；、、、和为系数，通过改变系数值可以选取时间推进精度，本文依次取值为1、0、0、1、0。

2 湍流模拟方法

针对飞行器整机级高雷诺数湍流流动的数值模拟，目前应用最多的还是采用湍流模型的RANS方法，该方法模化了所有尺度的湍流结构，仅对流场平均量求解，计算量相对较小，且能够保证一定的精度。本文采用工程上常用的湍流模型Spalart-Allmaras(SA)一方程湍流模型模拟飞行器的湍流效应。

3 气动/运动一体化松耦合计算方法

本文采用松耦合气动力与飞行动力学耦合求解策略进行气动/运动耦合计算。松耦合是将流体控制方程Navier-Stokes方程和飞机动力学方程组在各自的时间域上独立求解，交错时间推进获得系统耦合解的算法。其中飞机刚体六自由度运动/动力学方程采用显式推进，在物理时间推进的每一步内，飞行动力学方程与Navier-Stokes方程之间只交互一次数据。具体的算法步骤如下：

1) 将第物理时刻的飞行器位移、姿态等参数(主要是当前状态下的高度、马赫数、姿态角(如迎角、侧滑角和滚转角)等)传递给流体控制方程。

2) 根据第时刻的位移、姿态等参数更新流场网格和边界条件设置，求解Navier-Stokes方程并获得第物理时刻的气动力和气动力矩。

3) 根据第物理时刻的气动力和气动力矩，飞行动力学方程进行显式推进求解，获得第+1物理时刻飞行器位移、姿态等参数，并传递给Navier-Stokes方程。

4) 进行第+1物理时刻的Navier-Stokes方程双时间步的隐式推进求解，获得第+1物理时刻飞行器的气动力和气动力矩，并传递给六自由度全量飞行动力学方程。

5) 按照上述步骤继续进行迭代，直到达到总的物理时间，仿真结束。

松耦合计算过程中的参数传递过程如图1所示。可以看到，松耦合的求解策略很大程度上保留了CFD程序和飞行仿真程序的独立性和模块化，只需编写数据交互的程序即可，避免了联立方程面临的求解难度。但由于两种方程积分时间的不同步，使得两者之间始终存在一个物理时间步长的延迟，需要选择足够小的时间步长才能保持整个耦合求解的稳定性。

图1 松耦合的参数传递与计算流程Fig.1 Parameter transfer and calculation flow of loose coupling

4 算例验证

为了考核CFD计算程序模拟动姿态非定常流动的能力，以及刚体动力学/运动学方程与流动控制方程耦合求解算法的有效性，采用二维NACA0012翼型强迫俯仰振荡的标准算例来进行验证。

4.1 强迫俯仰振荡动力学方程

如图2所示，在亚声速来流情况下模拟二维NACA0012翼型的强迫俯仰振荡运动。常规强迫俯仰振荡运动采用解析形式，定义为翼型迎角随时间变化的正弦函数：

()=+sin()=+sin(2)

(7)

图2 NACA0012翼型强迫俯仰振荡示意图Fig.2 Forced pitching oscillation diagram of NACA0012 airfoil

在本算例中，有黏、可压缩自由来流的马赫数为=0.755，雷诺数为=5.5×10。设定翼型俯仰振荡的轴心处于距翼型前缘点0.25倍弦长的弦线上，翼型在=0016°初始迎角自由来流作用下开始绕俯仰轴做强迫俯仰振荡运动，俯仰振荡振幅=251°。

为了考核刚体动力学/运动学方程与流动控制方程耦合求解方法的有效性，本文将强迫俯仰振荡运动采用等价的微分方程形式给出：

(8)

通过联立流动控制方程与式(6)即可实现俯仰振荡计算。设置物理时间步长为0.001 s，非定常子迭代步数设为100步。

4.2 计算条件设置

如图3所示，为NACA0012翼型的C型拓扑网格，规模为500×411(流向×法向)，物面第一层网格节点距壁面的无量纲距离满足<1。远场边界距翼型为20倍的弦长，弦长为1 m。边界条件设置：计算域采用远场自由来流条件，翼型物面采用无滑移壁面边界条件。采用数值模拟手段

图3 NACA0012翼型C型网格Fig.3 NACA0012 airfoil C-grid

来模拟这一非定常流动现象，通过与实验数据的对比来验证耦合算法和计算程序模拟非定常流动的准确性。

4.3 计算结果分析

为了提高计算效率，这里先在初始迎角的自由来流条件下计算定常流动，并作为非定常计算的稳定初始流场。计算得到的稳定初场压力等值线分布如图4所示。从图4(a)中可以看到，由于初始迎角较小，同时NACA0012翼型又是对称翼型，压力分布大致呈现上下对称。这与图4(b)中文献[13]的结果较为吻合，说明了本文CFD程序具备精确计算定常流动的能力。

图4 NACA0012翼型压力等值线分布对比Fig.4 Comparison of pressure isoline distribution of NACA0012 airfoil

在得到稳定的初始流场后，开始进行非定常的计算。采用松耦合算法进行了耦合计算，将计算得到的升力系数及俯仰力矩系数随迎角的变化曲线与AGARD的实验数据及Batina的计算结果进行了比较，如图5所示。图中：为升力系数；为俯仰力矩系数。

从图5中可以看到，升力系数的数值模拟结果与Batina的计算结果吻合度较高，与实验值也较为吻合，只存在较小的偏差；俯仰力矩系数的计算结果总体上与实验基本吻合，存在小量偏差的原因可能是：翼型绕0.25弦长点进行俯仰振荡对流场模拟的精度较为敏感，力矩特性的数值模拟难度较大，因此计算结果吻合度不如升力系数吻合度高。

图5 NACA0012翼型强迫俯仰振荡过程中的升力系数和俯仰力矩系数Fig.5 Lift coefficient and pitch moment coefficient of NACA0012 airfoil during forced pitching oscillation

上述案例验证了本文的CFD数值模拟程序能够较为精确的模拟非定常流动，同时说明了本文提出的耦合算法是可行的。

5 机翼前缘结冰条件下飞机的气动/运动耦合特性

由于结冰动态增长的时间尺度、Navier-Stokes方程解算的时间尺度和飞行动力学方程解算的时间步长各不相同，且差距较大，因此选取合适的物理时间步长至关重要。由于冰形增长网格重构的复杂性及背景飞机整机带冰构型巨量的网格将对计算效率造成不可忽视的影响，为了有效解决该问题，做出以下假设：

1) 由于冰的增长速率与形状会受到外部气象环境参数、飞行速度、飞行姿态等时变因素的影响，需要在每一物理时间步长进行积冰冰形计算。这里简化了该项步骤，根据NASA Lewis中心IRT(Icing Research Tunnel)冰风洞公开实验数据直接给定不同时刻的冰形，同时不考虑冰形沿飞机机翼展向的变化特征。

2) 根据NASA Lewis中心IRT冰风洞公开实验数据和本文结冰实验数据表明：冰形的变化速率虽然受到不同的来流条件和飞行参数的影响而发生变化，但总体来看，冰形的增长是缓慢的，通常长达几十分钟才能有较为明显的变化特征。因此，为了兼顾计算效率，在计算机翼积冰动态增长条件下气动/运动耦合特性时，每20 s更新一次结冰冰形网格，在20 s周期内假定冰形不发生变化。

3) 本文着重研究机翼前缘动态结冰条件下的气动/运动耦合特性，这里设定平尾前缘结冰冰形保持不变，不发生增长，发动机吊舱唇口处和垂尾前缘无结冰。

5.1 机翼带冰构型全机三维数模及结构化网格生成方法

结合NASA Lewis中心IRT冰风洞公开数据，形成了一组背景飞机在不同时刻的机翼二维结冰冰形的输入如图6所示。根据每一时刻的二维冰形，基于该数模能够快速建立沿机翼展向的三维冰模如图7所示。机翼带冰构型三维数模能够反映不同时刻二维冰形的几何细节特征。建立了适合背景飞机带冰构型气动特性-飞行力学特性在线耦合分析的机翼带冰构型计算网格拓扑如图8所示。本文采用网格生成工具ICEM生成点对接多块网格，针对不同时刻的机翼带冰构型，能够在同一套网格拓扑下生成相应的全机三维结构化网格。

图6 机翼带冰构型数模及不同时刻二维冰形输入Fig.6 Numerical model of wing icing configuration and two-dimensional ice shape input at different time

图7 基于二维冰形的机翼带冰三维数模Fig.7 Three-dimensional numerical simulation of wing with ice based on two-dimensional ice shape

由于网格质量直接影响数值模拟结果的好坏，这里充分考虑到冰增长带来的网格重构的复杂性，在生成三维结构化网格时均进行质量检查，对于几何外形变化较大的区域进行局部加密处理，以保证数值模拟计算的精度。

图8 机翼结冰构型三维结构化计算网格拓扑与三维结构网格的生成Fig.8 Grid topology and generation of three-dimensional structural grid for wing icing configuration

5.2 非定常气动力求解器与飞行动力学耦合仿真平台的开发

结合多块结构化网格高质量计算网格生成技术，开发了求解非定常Navier-Stokes方程的CFD程序。该程序具备多核并行计算、计算条件输入和流场后处理的功能，并能实现与飞行动力学仿真程序的数据交互。

图9给出了进行耦合计算的流程。在非定常气动力计算模块中，基于上一时间段的飞行力学数据输入，将ICEM提前生成好的当前时刻结冰冰形的多块结构化网格输入到NASA著名的开源计算流体力学软件CFL3D进行流场计算，得到本时间段的气动力数据传递给飞行动力学特性计算分析模块进行飞行仿真计算。进行精确的耦合分析时，时间段为每一物理时间步长，通常设置在10量级。如果兼顾计算效率，在结冰缓慢增长带来气动力变化幅度较小的前提下，时间段可定义为结冰增长的时间尺度，本文设为20 s，即每20 s更新一次结冰冰形的网格，在此期间与飞行力学交互时冰形的网格不做更新，只改变来流条件。

图9 结冰条件下背景飞机气动/运动耦合计算流程Fig.9 Aerodynamic and flight dynamics coupling calculation process of background aircraft under icing condition

本文在进行流场计算时，首先基于RANS方法获得充分发展的初始流场，在稳定初场的基础上进行后续非定常计算至非定常流场基本稳定，在进行动态结冰条件下的流场计算时，流场在新结冰外形下从本时间段的稳定流场开始计算至非定常趋势明显时结束(具体根据气动力系数在稳定值上下波动的幅值)，为该时间段的非定常气动力计算全过程，得到较为精确的大飞机结冰复杂外形下流动的非定常效应。然后进行耦合交互计算，实现非定常气动力与飞行动力学的耦合分析。

5.3 算例分析

5.3.1 三维整机固定冰形条件下的气动/运动耦合特征

对10 min结冰构型(中度结冰)下的背景飞机气动/运动耦合特性进行了研究。背景飞机10 min结冰构型下的初始飞行状态设定为：以高度6 000 m， 飞行速度为150 m/s(=0.473 5)定直平飞，飞机空重72 000 kg。假设初始时刻飞机迎角和舵偏角为未结冰时该飞行状态的配平量，配平迎角为6.26°，升降舵偏角为=-6.71°，油门偏度为=34.91%。驾驶员在=0.2 s时感知飞机结冰带来的低头力矩，操纵舵面使飞机抬头，此时=-14.11°。同时本仿真条件中不考虑非对称结冰情形，通过耦合计算分析在2 s内飞行器姿态及流场的非定常变化特征。设置物理时间步长为0.01 s，进行非定常计算时子迭代步数设为25步，采用松耦合计算策略，因此在2 s的总物理时间域内耦合交互次数为200次。由于在此期间冰形不发生变化，因此不考虑网格的重构，只需改变来流条件。

图10给出了通过耦合计算得到的该状态条件下飞行姿态的变化，并对迎角响应过程选取了6个时刻点进行流场分析。可以看到主要是短周期参数迎角和俯仰角速度变化较为剧烈，飞行速度和高度虽然总体变化不大，但均呈现下降趋势，且下降率越来越大，这是由于迎角超过失速迎角后升力系数陡降导致的。此时飞机处于失速状态，如果驾驶员不能及时改出，可以预见，飞行高度下降率会急剧增大，飞机将在较短时间内损失较多的高度，对飞行安全产生威胁。图11为耦合计算过程中升力系数、阻力系数及俯仰力矩系数的变化曲线。可以看到，随着迎角不断增大，气动力及力矩系数的脉动越来越大，非定常效应越来越显著，导致俯仰角速度发生振荡。当迎角超过12°左右时，由于失速，升力系数逐渐降低。

图10 固定冰形条件下耦合计算过程中飞行参数 变化曲线Fig.10 Change curves of flight parameters in coupling calculation under fixed ice shape conditions

图11 固定冰形条件下耦合计算过程中气动力及力矩系数变化曲线Fig.11 Change curves of aerodynamic force and moment coefficients in coupling calculation under fixed ice shape conditions

图12给出了背景飞机在6个时刻下的翼面分离流动形态(左图)及空间流场结构(右图)。在右图中截取了位于机翼9个不同位置的横截面，画出了不同时刻流向速度的流线分布，而整机是压力系数的分布云图。由图可以看出，在10 min结冰构型、中等迎角状态(迎角为6.26°时)时，内外翼分离同时出现，始发位置均位于机翼前缘，平尾处流动无分离。表明结冰使得分离流动始发迎角大幅提前，在中等迎角状态下就表现出一定的翼面分离流动趋势。随着迎角增加，在迎角为8°时就已经出现了明显的翼面分离流动。在迎角为12°时，空间流动形态发生较大变化，流动分离效应较为显著，内翼和外翼几个截面处的分离泡完全破碎，再附效应消失，产生大范围的强回流区域。平尾处开始发生分离。随着迎角的继续增大，机翼处流动分离效应则更为明显，而平尾处分离流动形态未发生较大变化。

图12 固定冰形条件下耦合计算过程中不同时刻飞机翼面分离流动形态及空间流场结构Fig.12 Separation flow pattern and spatial flow field structure of aircraft airfoil at different time in coupling calculation process under fixed ice shape conditions

图13给出了采用准则计算得到的在=0 s迎角为6.26°及=1.98 s、迎角为16°时流动的涡结构。可以看到，在大迎角情形下机翼及平尾处的分离涡较为明显，而小迎角条件下未发现明显的分离涡。由于是采用非定常RANS方法和S-A模型对流场进行求解，该方法对强非定常涡运动的模拟无法达到较高精度。因此计算得到的分离涡只能作定性的分析，下一步需要精细化分析流动中涡脱落与非定常特征还需采用RANS/LES混合方法等。

图13 固定冰形条件下耦合计算过程中不同迎角下的分离涡结构Fig.13 Separated vortex structure at different angles of attack in coupled calculation under fixed ice shape conditions

5.3.2 三维整机动态结冰条件下的气动/运动耦合特征

采用相同的方法对整机动态结冰条件下的气动/运动耦合特征进行分析。在该模拟情形下，冰形是不断增长变化的。仿真模拟的条件设置为：初始飞行高度6 000 m，飞行速度150 m/s，飞机空重42 000 kg，初始航迹俯仰角为-2.5°，并做直线下滑。1 min后穿越含大量过冷水滴的云层，开始逐渐积冰。在整个飞行过程中舵偏量及油门开度保持为初始配平值。假定飞机持续结冰5 min，而后保持不变。总的物理时间设定为360 s。

首先根据设置的条件计算配平量，配平迎角为2.98°，配平升降舵偏角为-4.59°，配平油门偏角为12.41(大小范围为0～100)，侧滑角、副翼偏角及方向舵偏角均为0°。由于在飞行过程中飞机逐渐结冰，导致飞行姿态发生变化。且结冰增长的时间较长，为了减小耦合计算的工作量，采用变步长的思想：在前1 min未结冰期间内，时间步长设为0.1 s，且不进行非定常气动力的解算，只解算定常的气动力，这是由于迎角较小且飞机未结冰，非定常效应较为微弱可忽略不计；在5 min结冰增长期间，视结冰气动力系数变化量的大小来改变时间步长，当飞行姿态改变较小导致计算得到的气动力变化不大时，时间步长可取为0.1 s，变化较大时可取为0.02 s，以实现减小计算量的同时尽量不影响计算精度。

图14为通过耦合计算得到的飞行参数的变化曲线。图15为计算得到的升力系数、阻力系数及俯仰力矩系数变化曲线。由图可以看出，在动态结冰的过程中，飞机的阻力系数不断增大，俯仰力矩系数不断减小。在配平状态下，气动力产生的俯仰力矩与发动机推力产生的俯仰力矩平衡；在动态结冰状态时，由于结冰导致气动力产生的俯仰力矩出现变化，与发动机推力产生的俯仰力矩不能平衡，导致迎角不断增大。而升力系数随着迎角的增大而增大，直到飞机迎角超过失速迎角后，升力系数逐渐降低。随着迎角的增大，非定常效应越来越显著，气动力及力矩系数脉动现象逐渐加强。由于总的俯仰力矩系数存在上下振荡，导致俯仰角速度的波动较为明显，从而使飞机的俯仰姿态也出现了一定的振荡。在整个飞行过程中控制舵面和油门开度没有变化，保持在初始配平位置，但阻力的增大，高度和速度整体呈现下降趋势。

图14 动态结冰条件下耦合计算过程中 飞行参数变化曲线Fig.14 Change curves of flight parameters in coupling calculation under dynamic icing conditions

图15 动态结冰条件下耦合计算过程中气动力及力矩系数变化曲线Fig.15 Change curves of aerodynamic force and moment coefficient in coupling calculation under dynamic icing conditions

图16给出了不同时刻飞机翼面分离流动形态及空间流场结构特征。由图可以看出，在未结冰、小迎角情形下飞机翼面流动较为平稳；随着结冰的增长，2 min结冰构型在迎角为5.62°时，翼面沿展向已出现分离特征，从空间流场结构来看未见较大的分离流动；3 min结冰构型的流动情况与之类似，在内翼面后侧出现一定的回流区域，同时外翼的展向横流效应均更加明显；当冰形增长到4 min结冰构型时，在迎角为7.21°时内翼面后侧出现了显著的分离，从空间流场分布来看，内翼后缘出现再附效应，并显示了较高的回流强度，从表面压力分布来看，分离区翼面压力增加，形成了分离泡结构。随着迎角进一步增加和冰形的增长，在5 min结冰构型时，分离泡体积迅速膨胀，导致内翼面形成了大范围回流区域。总体来看，外翼处的分离流动趋势要显著弱于内翼面，这是由于机翼/短舱等部件之间的气动干扰削弱了分离泡的弦向扩展过程，同时抑制了横向流动，降低了当地流动分离的趋势。

图16 动态结冰条件下耦合计算过程中不同时刻飞机翼面分离流动形态及空间流场结构Fig.16 Separation flow pattern and spatial flow field structure of aircraft wing at different time in coupling calculation process under dynamic icing conditions

6 结 论

1) 本文针对机翼前缘结冰条件下空气动力学与飞行力学耦合特性进行了分析研究。通过对三维整机固定冰形、动态结冰条件下的气动/运动耦合特性分析，算例计算结果表明：结冰条件下非定常效应的影响在中等迎角阶段就开始显著，随着结冰的动态增长和迎角的增大，内翼分离流动趋势进一步增强，表现为分离泡体积的迅速膨胀。由于中等迎角条件下显著的非定常效应将导致飞机出现一定程度的振荡，影响飞行稳定性，此时驾驶员应该及时介入，解除危险的飞行状态。

2) 考虑到计算量的限制，本文没有考虑冰形沿机翼展向的变化特征。同时在进行动态结冰的流场计算时，网格在新结冰外形下是重新生成的。在后续的研究工作中需要减小这两处的模拟误差。