CN107092718B - 一种飞行器飞行中遭遇降雨时的数值模拟方法 - Google Patents

Abstract

本发明涉及一种飞行器飞行中遭遇降雨时的数值模拟方法，包括步骤：生成飞行器流场空间计算网格；确定飞行条件和降雨条件；求解空气相流体气动性能数据和空气相流体壁面剪切应力；求解飞行器表面局部雨滴收集率和飞行器表面局部雨滴撞击质量流量；求解雨滴撞击力系数增量；取飞行器流场空间计算网格中的翼型表面贴体网格作为表面水膜控制体，设定翼型表面同时存在气流附着点、气流转捩点和气流再附点，依次求解所有表面水膜控制体的水膜厚度；求解雨滴溅射模型的四个参数；求解等效的水膜层表面粗糙度；根据水膜层厚度更新飞行器几何形状，重构飞行器流场空间计算网格；获得更新空气相流体气动性能数据；获得降雨引起的飞行器气动性能损失数据。

Description

一种飞行器飞行中遭遇降雨时的数值模拟方法

技术领域

本发明属于计算流体力学在航空工程领域的应用，尤其涉及一种飞行器飞行中遭遇降雨时的数值模拟方法。

背景技术

飞行器在空中执行飞行任务期间，经常会穿越云层而遭遇可能的降雨。在降雨环境下飞行时，雨水会与飞行器表面碰撞产生附加力和附加力矩，且雨水会附着于飞行器表面即时形成一层流动的水膜，不但增加了飞行器的质量，而且改变了飞行器的气动外形，对飞行器的空气流场特性造成影响，使飞行器气动性能遭到损失，危害飞行器的平衡性、稳定性、操纵性和飞行品质，影响飞行安全。此外，不容忽视的是，降雨使飞行器舱外视景能见度降低，影响飞行员的视线；雨水渗透可能使机载仪表设备失灵，影响飞行员的判断能力，可能导致仪表飞行困难；雨水被发动机吸入后，可能导致发动机性能恶化，严重时可能导致发动机停车。由于降雨使飞行器气动性能恶化，再加上述可能存在的不利因素，不可避免地使飞行器的飞行安全受到威胁，引发可能的飞行安全事故。

早在上世纪八十年代，为了保证飞行器在降雨环境下的飞行安全性，世界航空强国就已开始进行相关研究。通过对降雨导致的飞行安全事故的调查、分析与研究，人们逐渐意识到气动性能下降是降雨引发飞行器飞行安全事故的首要根源，因此众多研究者对该问题展开了不懈的探索与研究。对降雨致使飞行器气动性能恶化的主要研究手段有飞行试验、风洞试验和计算机数值模拟，其中，飞行试验作为最直接的研究方法，要求在飞行器上安装仪器和设备，以精确记录环境参数，同时飞行员还要在仪表飞行条件下探索安全的飞行包线；因此，这种方法不仅成本高昂、风险性大，而且受自然气候的限制，即使在较长试验期内可能也无法验证完所有自然降雨条件，研究周期长。相比于飞行试验，风洞试验安全系数较大，但试验中较难精确产生飞行中遭遇的雨滴尺寸分布和雨滴运动速度分布，且在试验中不仅要精确模拟降雨，还要确定试验相似准则，并通过技术手段将试验结果应用于全尺寸飞行器；显然，这对试验仪器和试验设备的技术指标以及试验技术手段的要求都相当高，成本依然十分昂贵。计算机数值模拟方法可以利用计算机高级程序语言实现，并通过计算机的运行来模拟飞行器在空中飞行遭遇降雨时的状态，是计算流体力学在航空工程领域的应用；先进的数值模拟技术不但可以减少飞行试验、降低风险，而且还能弥补风洞试验中的误差与缺陷。事实上，在意识到降雨对飞行器恶劣影响的根源后，国际上就有人逐步开始采用数值模拟方法进行研究，不仅在理论上对降雨这一复杂现象有了深刻认识，而且对降雨引起飞行器气动性能损失的机理也有了一定理解，甚至人们开始尝试量化研究降雨引起的飞行器气动性能损失。在研究初期，因在数值模拟技术中仅考虑了雨滴体积力的影响，所得结果与真实情况相差较大；经过数十年的发展，在逐步计入水膜层影响和雨滴溅射影响后，数值模拟结果与真实情况基本符合。

目前，研究者大多基于两相流理论研究降雨对飞行器气动性能的恶劣影响。有人采用欧拉方法描述，将离散的运动雨滴视为与空气相互渗透的连续拟流体，建立雨滴相流体流动控制方程；还有人采用拉格朗日方法描述，直接追踪各个离散雨滴的运动情况，建立雨滴运动方程。对飞行器外部流场的求解方法，有人采用定常流场求解算法，也有人采用非定常流场求解算法。对空气相流体湍流模型的选择上，k-ε湍流模型、k-ω湍流模型和S-A湍流模型等被经常用来作为求解中对空气相流体湍流的模拟。尽管在描述方法、求解算法和湍流模型上有一些差异，但研究降雨对飞行器气动性能恶劣影响的数值模拟方法所采用的流程大致可归纳总结为四个主要功能模块：空气流动模块，用于求解飞行器外部空气相流体流场，即求解空气相流体流动控制方程；雨滴运动模块，用于求解飞行器外部雨滴相流体流场以及雨滴与飞行器表面的碰撞过程，获得附着于飞行器表面的雨水状态(用雨滴收集率表示)，即求解雨滴相流体流动控制方程(或雨滴运动方程)；水膜运动模块，用于求解附着于飞行器表面的水膜生长与发展过程，获得附带水膜的飞行器几何形状；雨滴溅射模块，用于求解雨滴撞击飞行器表面引起的水膜表面粗糙度状态。

如图1所示，研究降雨对飞行器气动性能恶劣影响的数值模拟技术中各模块之间的关系和求解过程为：首先确定飞行条件和生成计算网格，展开空气相流体流场的求解。然后基于降雨条件，结合雨滴相流体流动控制方程(或雨滴运动方程)求解雨滴运动情况。再将获得的数据，如壁面雨滴收集率、壁面雨滴法向撞击速度、壁面空气相流体剪切应力、壁面空气相流体剪切速度等，带入水膜运动模块，根据层流流动的水膜仅受气流摩擦作用的假设条件，利用水膜内水流流动过程中质量守恒原则，在时域内计算飞行器表面各处流动水膜的厚度，获得附带水膜的飞行器几何形状，并围绕这个几何形状重构计算网格。基于雨滴运动模块和水膜运动模块，进一步求解雨滴溅射模块，获得水膜表面粗糙度状态。最后利用更新的计算网格，展开空气相流体绕流流场的更新求解，模拟飞行器在飞行中遭遇降雨时的状态，获得气动性能损失数据。

上述研究降雨对飞行器气动性能恶劣影响的数值模拟方法，在求解水膜运动模块时简单地将壁面空气相流体剪切运动考虑为以气流附着点为分界点，分别沿飞行器上、下翼面作方向一致性的流动，流至飞行器尾部；如此，壁面水膜内的雨水受壁面空气相流体剪切力作用，也以气流附着点为分界点，分别沿飞行器上、下翼面作方向一致性的流动，流至飞行器尾部。然而，这种处理方式仅适用于小迎角飞行条件，在大迎角情况下，如接近失速或失速时，飞行器翼面除了存在气流附着点外，还可能存在气流转捩点和(或)气流再附点。在气流转捩作用下，飞行器上翼面流动水膜内的雨水可能不再紧贴上翼面作方向一致性的流动而流向飞行器尾部，而是在气流转捩点附近被气流卷起，溅至上翼面附近的流场空间中，而转捩点之后的壁面上可能不存在水膜；也可能由于受气流逆流或气流再附的作用，使壁面上收集到的雨水一方面逆流至转捩点附近被气流卷起、溅至上翼面附近的流场空间中，一方面顺流至飞行器上翼面尾部、消散于飞行器尾部流场空间中。因此，上述小迎角假设条件很明显与实际不符。上述研究降雨对飞行器气动性能恶劣影响的数值模拟方法所使用的雨膜厚度计算算法，基于时间域内的水膜厚度增长方程，随着时间的增长，由于质量守恒，水膜的雨水流动必然趋于稳定，水膜厚度也就趋于稳定；但该算法受数值波动影响较大，计算时间较长，计算稳定性较差，尤其受时间步长影响，计算精度低。此外，上述研究降雨对飞行器气动性能恶劣影响的数值模拟方法，在求解雨滴溅射模块时采用的理论计算没有有效考虑雨滴溅射时的能量耗散问题，使得理论值与试验值相差较大，而雨滴溅射试验数据十分有限，又没有进行合理拓展，使得雨滴溅射模块的计算条件相当苛刻，这在数值模拟时必然会产生一定误差，导致结果不精确。

发明内容

针对上述问题，本发明的目的是提供一种飞行器飞行中遭遇降雨时的数值模拟方法，接近真实的降雨条件和飞行条件，兼顾计算精度、效率和功能，探索飞行器表面水膜运动特性及由此引起的气动性能损失的机理，定量分析降雨引起的飞行器气动性能损失。

为实现上述目的，本发明采取以下技术方案：一种飞行器飞行中遭遇降雨时的数值模拟方法，包括以下步骤：

1)根据目标飞行器的几何外形，生成飞行器流场空间计算网格；

2)根据目标飞行器的实际使用条件，确定飞行条件和降雨条件；

3)根据飞行器流场空间计算网格、飞行条件和降雨条件，进行空气相-雨滴相两相流体流动模拟计算，获得流场数据；

4)根据流场数据获得空气相流体气动性能数据和空气相流体壁面剪切应力；

5)利用动量定理建立模拟雨滴对飞行器表面撞击效应的计算式，求解飞行器表面局部雨滴收集率和飞行器表面局部雨滴撞击质量流量；

6)根据飞行器表面局部雨滴撞击质量流量，求解雨滴撞击力系数增量；

7)取飞行器流场空间计算网格中的翼型表面贴体网格作为表面水膜控制体，采用QUICK格式离散步骤4)得到的空气相流体壁面剪切应力，得到各表面水膜控制体边界处的空气相流体壁面剪切应力；

8)设定翼型表面同时存在气流附着点、气流转捩点和气流再附点，分别记气流附着点、气流转捩点和气流再附点所在的表面水膜控制体为i₀控制体、i₁控制体和i₂控制体，并初始化i₀控制体、i₁控制体和i₂控制体的位置；

9)求解i₀控制体水膜厚度；

10)求解i₀控制体左侧界面方向所有表面水膜控制体的水膜厚度；

11)检测i₁控制体的存在性，若存在，更新其位置；否则，不予更新；

12)求解i₀控制体与i₁控制体之间的表面水膜控制体的水膜厚度；

13)检测i₂控制体的存在性，若存在，更新其位置；否则，不予更新；

14)求解i₂控制体水膜厚度；

15)求解i₁控制体与i₂控制体之间表面水膜控制体的水膜厚度；

16)求解i₂控制体右侧边界方向所有表面水膜控制体的水膜厚度；

17)整理计算得到的所有表面水膜控制体的水膜厚度；

18)建立雨滴溅射模型，求解雨滴溅射模型的四个参数：无量纲溅射凹坑高度实际修正值、无量纲溅射凹坑半径、无量纲溅射凹坑形成时间和无量纲凹坑间隔距离；

19)根据雨滴溅射模型的四个参数，求解等效的水膜层表面粗糙度；

20)根据步骤17)得到的所有表面水膜控制体的水膜层厚度更新飞行器几何形状，根据带有附着水膜界面的、新的飞行器几何形状重构飞行器流场空间计算网格；

21)引入步骤19)得到的等效的水膜层表面粗糙度，根据重构的飞行器流场空间计算网格，以及飞行条件和降雨条件，再次进行空气相-雨滴相两相流体流动模拟计算，获得更新流场数据，根据更新流场数据获得更新空气相流体气动性能数据；

22)将步骤21)得到的更新空气相流体气动性能数据与步骤4)得到的空气相流体气动性能数据相减，再和步骤6)得到的雨滴撞击力系数增量数据叠加，获得降雨引起的飞行器气动性能损失数据。

所述步骤2)中的飞行条件包括飞行器飞行高度H、飞行高度H处的空气密度ρ_H；降雨条件包括液态水含量LWC、雨滴体积当量平均直径MVD和雨滴下落的末速V_H；

其中，确定降雨条件的具体方法为：先定义降雨类型和降雨率R，计算液态水含量LWC和雨滴体积当量平均直径MVD；根据飞行器飞行高度H，计算雨滴下落的末速V_H；

或者，直接定义液态水含量LWC和雨滴体积当量平均直径MVD；根据飞行器飞行高度H，计算雨滴下落的末速V_H。

所述步骤3)中采用非定常、不可压缩气体的雷诺平均N-S方程作为空气相流体流动控制方程，采用高雷诺数S-A湍流模型模拟空气相流体湍流，由单个雨滴运动方程推导得到雨滴相流体流动控制方程，采用一阶隐式非定常求解算法进行空气相-雨滴相两相流体流动模拟计算。

所述步骤5)中飞行器表面局部雨滴收集率的计算公式为：

式中，i为翼型表面贴体网格计数，以顺时针方向按序计数，翼型下表面后缘贴体网格记为第1个网格，翼型上表面后缘贴体网格记为第N_max个网格，N_max为翼型表面贴体网格总数目；β_i为第i个贴体网格处的局部雨滴收集率，α_V0＝LWC/ρ_w；α_VS,i和V_⊥S,i分别为第i个贴体网格的雨滴相流体体积分数和雨滴法向撞击速度；V_∞为雨滴相流体的初始速度，V_a为空气相流体的来流速度；ρ_w为雨水密度；

飞行器表面局部雨滴撞击质量流量的计算公式为：

式中，

为第i个贴体网格处的局部雨滴撞击质量流量；A_p,i为第i个贴体网格的壁面表面积。

所述步骤6)中雨滴撞击力系数增量的计算公式为：

式中，ΔC_L,Impinging和ΔC_D,Impinging分别是升力和阻力系数增量；F_i为第i个贴体网格处雨滴对飞行器表面局部的撞击力矢量；F为雨滴对飞行器表面总的撞击力矢量；F_x和F_y分别为撞击力矢量F沿计算网格坐标系x轴和y轴的分量；α为来流迎角。

所述步骤9)中，取第i个表面水膜控制体的空气相流体壁面剪切应力τ_i的符号与前一表面水膜控制体的空气相流体壁面剪切应力τ_i-1的符号刚好不同时的表面水膜控制体为i₀控制体，然后进行一次修正；i₀控制体的水膜厚度

的计算公式为：

式中，

和

分别为i₀控制体左右两个界面处的空气相流体壁面剪切应力，通过QUICK格式离散得到；

为i₀控制体处的局部雨滴撞击质量流量；μ_w为水的粘性系数；

所述步骤11)中，取第i个表面水膜控制体的空气相流体壁面剪切应力τ_i的符号与前一表面水膜控制体的空气相流体壁面剪切应力τ_i-1的符号刚好不同时的表面水膜控制体为i₁控制体，然后进行一次修正；若存在气流转捩点，则更新i₁，否则i₁值不变；

所述步骤13)中，取第i个表面水膜控制体的空气相流体壁面剪切应力τ_i的符号与前一表面水膜控制体的空气相流体壁面剪切应力τ_i-1的符号刚好不同时的表面水膜控制体为i₂控制体，然后进行一次修正；若存在气流再附点，则更新i₂，否则i₂值不变；

所述步骤14)中，i₂控制体的水膜厚度

的计算公式为：

式中，

和

分别为i₂控制体左右两个界面处的空气相流体壁面剪切应力，通过QUICK格式离散得到；

为i₂控制体处的局部雨滴撞击质量流量。

所述步骤10)中，从i₀控制体开始，沿左侧界面方向依次求解第i个表面水膜控制体左侧界面处的质量流量

(i∈[1,i₀-1])，然后依次求解i₀控制体左侧界面方向所有控制体的水膜厚度

其中，任意第i个表面水膜控制体左右两个界面处的质量流量计算公式为：

任意第i个表面水膜控制体左侧界面处的水膜厚度的计算公式为：

所述步骤12)中，从i₀控制体开始，沿右侧界面方向依次求解

(i∈[i₀+1,i₁-1])；然后依次求解i₀控制体与i₁控制体之间的表面水膜控制体的水膜厚度；i₀控制体与i₁控制体之间表面水膜控制体的水膜厚度的计算公式为：

所述步骤15)中，i₁控制体与i₂控制体之间表面水膜控制体的水膜厚度的计算公式为：

所述步骤16)中，从i₂控制体开始，依次展开计算i₂控制体右侧边界方向所有表面水膜控制体的水膜厚度；i₂控制体右侧边界方向所有表面水膜控制体的水膜厚度的计算公式为：

上述各式中，

和

分别为第i个表面水膜控制体左右两个界面处的空气相流体壁面剪切应力，通过QUICK格式离散得到；

和

分别为第i个表面水膜控制体左右两个界面处的水膜厚度；

为第i个表面水膜控制体处的局部雨滴撞击质量流量。

所述步骤18)中求解雨滴溅射模型的四个参数，具体包括以下步骤：

①求解无量纲溅射凹坑高度的理论值

为：

式中，W_b为雨滴韦伯数；F_N为雨滴弗洛德数；h^*为无量纲水膜厚度；

为无量纲溅射凹坑半径；

②求解无量纲溅射凹坑高度的实际修正值

为：

③间接求解无量纲溅射凹坑半径

为：

④求解无量纲溅射凹坑形成时间

为：

⑤求解无量纲凹坑间隔距离

为：

式中，β为飞行器的局部雨滴收集率；I＝nR^-m，R为降雨率，m和n为常数，均与降雨类型相关，为降雨试验测量值；N₀为单位体积空间中的雨滴数密度随雨滴直径大小的分布常数；MVD为雨滴体积当量平均直径。

所述步骤19)中求解等效的水膜层粗糙度，具体包括以下步骤：

a、采用等效沙粒粗糙度模拟雨滴撞击飞行器表面引起的水膜表面粗糙度状态，求解雨滴溅射引起的水膜表面粗糙度k_s，c为：

式中，

系数C与撞击凹坑形状有关；D_c为凹坑间隔距离，

为凹坑当量平均高度，采用M-P雨滴谱分布模拟降雨时有

b、求解水膜层表面水波纹引起的表面粗糙度k_s,w：

k_s，w＝B·h

式中，取B＝1.5；h为水膜厚度，h＝MVD·h^*/2；

c、求解雨滴撞击与溅射引起的等效的水膜层表面粗糙度k_s为：

k_s＝||[k_s,c k_s,w]||。

本发明由于采取以上技术方案，其具有以下优点：1、本发明的一种飞行器飞行中遭遇降雨时的数值模拟方法，根据飞行器壁面附近流场流动特性，从判断翼型表面气流附着点、气流转捩点和气流再附点可能的存在性出发，摒弃传统的时间域内的水膜厚度增长方程，创新性地提出了水膜层生长与发展的算法，不仅能精确模拟飞行器小迎角下遭遇降雨时的状态，更能精确模拟飞行器大迎角条件下遭遇降雨时的状态。2、本发明的一种飞行器飞行中遭遇降雨时的数值模拟方法，基于试验数据，考虑雨滴溅射过程中的能量耗散因素，通过利用数值拟合技术，对试验数据拟合，得到直观、具体的表达式，使计算结果更加符合实际，使传统溅射模型适用性更为广泛。3、本发明的一种飞行器飞行中遭遇降雨时的数值模拟方法，可利用计算机高级程序语言实现，不仅能模拟飞行器正常飞行(小迎角情况)遭遇降雨时的状态，更能模拟飞行器机动飞行(大迎角情况)遭遇降雨时的状态，是对飞行器在不同飞行状态下遭遇降雨过程的认识和分析，能为研究降雨条件下飞行器安全执行各类飞行任务提供重要信息。

附图说明

图1是传统数值模拟计算流程图；

图2是本发明方法的流程图；

图3(a)、(b)分别是仅存在气流附着点时翼型表面附着水膜流动原理图和翼型局部壁面水膜流动平衡原理图；

图4(a)、(b)分别是存在气流附着点和气流转捩点时翼型表面附着水膜流动原理图和翼型局部壁面水膜流动平衡的原理图；

图5(a)、(b)分别是存在气流附着点、气流转捩点和气流再附点时翼型表面附着水膜流动原理图和翼型局部壁面水膜流动平衡的原理图；

图6是雨滴溅射模型示意图；

图7是围绕NACA0012翼型的初始计算网格图；

图8是雨膜厚度放大10倍显示计算得到的翼型表面附着水膜层状态图；

图9(a)、(b)、(c)是雨滴溅射模型试验数据与其数据拟合的对比图；

图10是净翼型与降雨条件下的翼型的升力系数曲线对比图；

图11是净翼型与降雨条件下的翼型的阻力系数曲线对比图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明根据动量定理建立模拟降雨过程中雨滴对飞行器表面撞击效应的算法；综合考虑气流附着点、气流转捩点和气流再附点对壁面空气相流体剪切运动的影响，建立飞行器表面附着水膜层的生长与发展的算法，使得模型能够模拟大迎角条件下的水膜运动；同时摒弃基于时间域内的水膜厚度增长方程，提出一种新的水膜层直接生长法，无需在时间域内进行迭代即得到水膜厚度，避免了数值计算波动，缩短了计算时间，增大了数值计算稳定性；根据雨滴溅射试验数据，采用数值拟合技术，建立雨滴溅射对飞行器表面粗糙度影响的算法，由于实际溅射需合理计入能量耗散，因而拓展了传统溅射模型的计算条件。本发明使用两相流体流动的模拟方法，基于欧拉方法描述，将运动雨滴视为连续拟流体，分别建立空气相流体流动控制方程、雨滴相流体流动控制方程，采用高雷诺数S-A湍流模型模拟空气相流体湍流，利用一阶隐式非定常求解算法展开两相流数值求解。基于计算网格、利用上述模型和算法进行飞行器飞行遭遇降雨时的数值模拟计算流程。完成水膜运动模块的计算后，根据带有附着水膜界面的新的飞行器几何形状，重构计算网格，引入基于雨滴溅射模块得到的水膜层表面粗糙度，展开空气相流体流动模块的更新求解，更新空气相流体气动性能数据，并与之前得到的空气相流体气动性能数据和雨滴撞击力系数增量叠加，获得降雨引起的气动性能损失数据。

基于上述原理，如图2所示，本发明提供的一种飞行器飞行中遭遇降雨时的数值模拟方法，具体包括以下步骤：

1)根据目标飞行器的几何外形，生成飞行器流场空间计算网格。

2)根据目标飞行器的实际使用条件，确定飞行条件和降雨条件。

大气中的液态水均匀分布于大气中，飞行器在飞行中进入液态水含量较高的空域时，易遭遇降雨。据统计，降雨时的雨滴直径较小，一般在2～3mm左右，很少超过6mm，而且雨滴的数密度会随雨滴直径的增大而迅速减小。降雨率R(mm/h)、液态水含量LWC(g/m³)、雨滴直径D(mm)及雨滴下落的末速V_H(m/s)是表征降雨条件的几个重要参数。其中，降雨率R是指单位时间内降落到水平地面上的雨水深度，液态水含量LWC是指单位体积空气中所含的液态水质量。在模拟降雨时，采用经典的M-P雨滴谱建立降雨率R与液态水含量LWC和雨滴直径D之间的量化关系，其中，雨滴谱分布函数为：

式中，I＝nR^-m，N₀(m^-3mm^-1)为单位体积空间中的雨滴数密度随雨滴直径大小的分布常数，m和n也为常数，它们均与降雨类型相关，为降雨试验测量值。

通常可采用如下两种方法确定降雨条件：

方法一：为模拟降雨条件，定义降雨类型和降雨率R，计算液态水含量LWC和雨滴体积当量平均直径MVD；定义飞行器飞行高度H(m)，计算雨滴下落的末速V_H。

方法二：为模拟降雨条件，直接定义液态水含量LWC和雨滴体积当量平均直径MVD；定义飞行器飞行高度H，计算雨滴下落的末速V_H。

其中，液态水含量LWC的计算公式为：

LWC＝ρ_wπN₀I^-4×10^-6 (2)

式中，ρ_w(kg/m³)为雨水密度；π为圆周率；

雨滴体积当量平均直径MVD的计算公式为：

MVD＝4/I (3)

飞行器飞行高度H处的雨滴末速V_H的计算公式为：

式中，ρ₀(kg/m³)为标准海平面处的空气密度；ρ_H(kg/m³)为飞行高度H处的空气密度。

一般采用液态水含量LWC、雨滴体积当量平均直径MVD和雨滴下落的末速V_H作为模拟降雨时的条件输入参数。

3)进行空气相-雨滴相两相流体流动模拟。

采用技术成熟的非定常、不可压缩气体的雷诺平均N-S方程作为空气相流体流动控制方程，采用高雷诺数S-A湍流模型模拟空气相流体湍流，由单个雨滴运动方程推导得到雨滴相流体流动控制方程，采用技术成熟的一阶隐式非定常求解算法进行求解。其中，雨滴相流体流动控制方程的表达式为：

式中，

为雨滴相流体表观密度，

α_v为雨滴相流体体积分数；x和y分别表示计算网格坐标系的x轴和y轴的方向；U(m/s)为雨滴相流体运动速度矢量，U＝[u_x,u_y]^T，u_x和u_y分别为U的沿计算网格坐标系x轴和y轴的分量；F_D和F_G是动量源项，分别为雨滴相流体单位体积质点所受的空气阻力和重力与空气浮力的合力矢量。

4)根据流场数据获得空气相流体气动性能数据和空气相流体壁面剪切应力等相关流场数据。

5)求解局部雨滴收集率和局部雨滴撞击质量流量。

采用飞行器表面局部雨滴收集率描述雨滴对飞行器表面的撞击效应，利用动量定理建立模拟雨滴对飞行器表面撞击效应的计算式。其中，飞行器表面局部雨滴收集率的计算公式为：

式中，i为翼型表面贴体网格计数，以顺时针方向按序计数，翼型下表面后缘贴体网格记为第1个网格，翼型上表面后缘贴体网格记为第N_max个网格，N_max为翼型表面贴体网格总数目；β_i为第i个贴体网格处的局部雨滴收集率，α_V0＝LWC/ρ_w；α_VS,i和V_⊥S,i分别为第i个贴体网格的雨滴相流体体积分数和雨滴法向撞击速度；V_∞为雨滴相流体的初始速度，V_a为空气相流体的来流速度。

飞行器表面局部雨滴撞击质量流量的计算公式为：

式中，

为第i个贴体网格处的局部雨滴撞击质量流量；A_p,i(m²)为第i个贴体网格的壁面表面积。

6)求解雨滴撞击对气动性能的影响，并以系数增量形式表示。

根据动量定理，第i个贴体网格处雨滴对飞行器表面局部的撞击力矢量为：

则雨滴对飞行器表面总的撞击力矢量记为：

将撞击力矢量F沿翼型升、阻力方向分解，以系数增量形式计入其对气动性能的影响如下：

式中，ΔC_L,Impinging和ΔC_D,Impinging分别是升力和阻力系数增量，F_x和F_y分别为撞击力矢量F沿计算网格坐标系x轴和y轴的分量；α为来流迎角。要注意的是，系数增量已带有符号。

7)采用QUICK格式离散步骤4)得到的空气相流体壁面剪切应力，得到各表面水膜控制体边界处的空气相流体壁面剪切应力。

取上述翼型表面贴体网格作为水膜运动状态的求解网格，称为表面水膜控制体。表面水膜控制体网格与翼型表面贴体网格沿壁面的分布尺寸及排序位置都是一致的，但表面水膜控制体的网格厚度取水膜厚度，因此表面水膜控制体与翼型表面贴体网格可以使用相同的标号表示，方便理解。传统的水膜运动模型中，水膜运动采用基于时间域内的水膜厚度增长方程进行模拟，第i个表面水膜控制体的水膜厚度h(i)表示为：

h(i)＝h₀(i)+(1/ρ_w)·(dm_i/dt)·Δt,(1≤i≤N_max) (14)

式中，h₀(i)为第i个表面水膜控制体的初始水膜厚度；dm_i/dt为第i个表面水膜控制体的质量变化率；Δt为时间步长。

随着时间增长，由于质量守恒，根据式(14)求解的水膜厚度必然趋于稳定。但该方法受数值波动影响较大，计算时间较长，计算稳定性较差；尤其受Δt影响，计算精度低。根据式(14)不难看出，一旦dm_i/dt＝0，水膜内雨水流动达到流量平衡，水膜厚度不再随时间变化。以此为突破点，结合翼型表面气流附着点、气流转捩点和气流再附点可能的存在性，分析每个翼型表面水膜控制体的流动情况，可以分三种情况：

第一种情况，如图3(a)、(b)所示，翼型表面仅存在气流附着点，即飞行器在小迎角情况下的飞行。此时，记气流附着点所在的水膜控制体为i₀控制体。根据翼型表面附近空气相流体的流场流动特性，气流分别沿翼型上、下表面流动，则在数值计算中，i₀控制体左右两个界面的水流流动方向均为界面法向朝外。为了达到水膜内雨水流动的流量平衡，对于i₀控制体右侧界面一侧所有控制体，其界面雨水流动方向均为界面法向朝右；对于i₀控制体左侧界面一侧所有控制体，其界面雨水流动方向均为界面法向朝左。最终，水膜内的雨水经流动到达翼型后缘，散于流场空间中。翼型表面水膜控制体的流量平衡方程为：

式中，

和

分别为第i个表面水膜控制体左右两个界面处的质量流量；

为i₀控制体处的局部雨滴撞击质量流量；

和

分别为i₀控制体左右两个界面处的质量流量。

由于任意相邻水膜控制体界面重合，即

根据式(15)和式(16)，有

式中，

为第1个表面水膜控制体左侧界面处的质量流量，

为第N_max个表面水膜控制体右侧界面处的质量流量。

这表明，翼型表面总的雨滴撞击质量流量(即翼型表面收集到的雨水质量流量)等于从翼型后缘处从水膜控制体界面流出的雨水的质量流量，从而实现了整个翼型表面附着水膜层内的雨水流动平衡。

假设水膜层内的雨水处于单向层流流动的状态，且雨水仅受水膜-空气界面的空气相流体剪切应力的作用，则水膜内任意厚度z处的界面法向流动速度为：

u＝τz/μ_w,(0≤z≤h) (18)

式中，μ_w为水的粘性系数；由于水膜层非常薄，τ可取空气相流体壁面剪切应力。

此时，任意第i个表面水膜控制体左右两个界面处的质量流量可记为：

式中，

和

分别为第i个表面水膜控制体左右两个界面处的空气相流体壁面剪切应力，通过QUICK格式离散得到；

和

分别为第i个表面水膜控制体左右两个界面处的水膜厚度。

在第一种情况下，初步取τ_i的符号与前一控制体τ_i-1的符号刚好不同时的控制体为i₀控制体，然后进行一次修正，从而达到寻找i₀控制体的目的。首先求解i₀控制体的水膜厚度

然后从i₀控制体开始，利用式(19)，沿左侧界面方向依次求解

(i∈[1,i₀-1])，再依次求解i₀控制体左侧界面方向所有控制体的水膜厚度；最后从i₀控制体开始，利用式(19)，沿右侧界面方向依次求解

(i∈[i₀+1,N_max])，再依次求解i₀控制体右侧界面方向所有控制体的水膜厚度；如此，即可得到所有水膜控制体左右两个界面所处的水膜厚度值，翼型附着水膜各处厚度得到求解。

第二种情况，如图4(a)、(b)所示，翼型表面同时存在气流附着点和气流转捩点，即飞行器在接近失速或失速情况下的飞行。此时，记气流转捩点所在的水膜控制体为i₁控制体，i₀控制体右侧界面方向的附着水膜内的雨水在运动至i₁控制体后，受壁面空气相流体的逆向剪切应力作用，不再紧贴翼型表面，而是在i₁控制体附近被卷起的气流带起，散于流场空间中。i₁控制体之后的壁面水膜是否还存在，还取决于这部分壁面是否收集到了雨滴。本发明假设在气流转捩点附近，水膜内的雨水均流入i₁控制体，且从i₁控制体所在水膜表面析出。

在第二种情况下，同样取τ_i的符号与前一控制体τ_i-1的符号刚好不同时的控制体为i₁控制体，然后进行一次修正，从而达到寻找i₁控制体的目的。首先求解i₀控制体的水膜厚度

然后从i₀控制体开始，利用式(19)，沿左侧界面方向依次求解

(i∈[1,i₀-1])，再依次求解i₀控制体左侧界面方向所有控制体的水膜厚度；接下来，从i₀控制体开始，利用式(19)，沿右侧界面方向依次求解

(i∈[i₀+1,i₁-1])，再依次求解i₀控制体与i₁控制体之间控制体的水膜厚度；最后，i₁控制体之后水膜厚度可从第N_max个控制体开始依次展开求解(此时，假设

)。

第三种情况，如图5(a)、(b)所示，翼型表面同时存在气流附着点、气流转捩点和气流再附点，即飞行器在完全失速情况下的飞行，此时，记气流再附点所在的水膜控制体为i₂控制体。

在第三种情况下，同样取τ_i的符号与前一控制体τ_i-1的符号刚好不同时的控制体为i₂控制体，然后进行一次修正，从而达到寻找i₂控制体的目的。首先求解i₀控制体的水膜厚度

然后求解i₀控制体左侧界面方向所有控制体的水膜厚度；接下来求解i₀控制体与i₁控制体之间控制体的水膜厚度；再之后，求解i₁控制体与i₂控制体之间控制体的水膜厚度；采用求解i₀控制体水膜厚度涉及到的方法求解i₂控制体的水膜厚度；最后，求解气流再附点右侧界面方向所有控制体的水膜厚度。

鉴于上述三种情况，本发明提出了新的飞行器表面附着水膜层的生长与发展的算法，采用水膜层的“直接生长法”来模拟水膜运动，无需在时间域内进行迭代，即得到水膜厚度；同时在模拟飞行器表面附着水膜层的生长与发展时采取了富有技巧性的处理方法，直接综合处理上述三种情况。在数值计算中，若不存在气流转捩点或气流再附点，则相应计算自然无效，无需作任何处理，在高级编程语言实现上简单易行。

8)设定翼型表面同时存在气流附着点、气流转捩点和气流再附点，并初始化i₀控制体、i₁控制体和i₂控制体的位置，令i₀＝N_max/2，i₁＝N_max+10，i₂＝N_max+20。

9)求解i₀控制体水膜厚度。具体包括以下步骤：

首先取τ_i的符号与前一控制体τ_i-1的符号刚好不同时的控制体为i₀控制体，然后进行一次修正，从而达到寻找i₀控制体的目的。即有：

假设i₀控制体为矩形形状，其它控制体为梯形形状，则有

结合式(15)、式(19)和式(21)，得到i₀控制体的水膜厚度

为：

10)求解i₀控制体左侧界面方向所有控制体的水膜厚度。

从i₀控制体开始，利用式(19)，沿左侧界面方向依次求解

(i∈[1,i₀-1])；然后依次求解i₀控制体左侧界面方向所有控制体的水膜厚度为：

11)检测i₁控制体的存在性，若存在，更新其位置；否则，不予更新。

同样取τ_i的符号与前一控制体τ_i-1的符号刚好不同时的控制体为i₁控制体，然后进行一次修正，从而达到寻找i₁控制体的目的。若存在气流转捩点，则更新i₁，否则i₁值不变。

12)求解i₀控制体与i₁控制体之间的控制体的水膜厚度。

从i₀控制体开始，利用式(19)，沿右侧界面方向依次求解

(i∈[i₀+1,i₁-1])；然后依次求解i₀控制体与i₁控制体之间控制体的水膜厚度为：

13)检测i₂控制体的存在性，若存在，更新其位置；否则，不予更新。

同样取τ_i的符号与前一控制体τ_i-1的符号刚好不同时的控制体为i₂控制体，然后进行一次修正，从而达到寻找i₂控制体的目的。若存在气流再附点，则更新i₂，否则i₂值不变。

14)求解i₂控制体水膜厚度。

假设i₂控制体为矩形形状，其它控制体为梯形形状，则可以结合式(15)和式(19)，并将i₂控制体的相关计算参数替换式(21)和式(22)中的相关计算参数，从而求解i₂控制体的水膜厚度。

15)求解i₁控制体与i₂控制体之间控制体的水膜厚度。

i₁控制体与i₂控制体之间控制体的水膜厚度为：

16)求解i₂控制体右侧边界方向所有控制体的水膜厚度。

从i₂控制体开始，依次展开计算i₂控制体右侧边界方向所有控制体的水膜厚度为：

17)整理计算得到的所有水膜厚度。

18)求解雨滴溅射模型的四个参数。

本发明基于传统雨滴溅射模型的描述方法，如图6所示，采用凹坑高度H_c、凹坑半径R_c、凹坑形成时间t_c和凹坑间隔距离D_c四个参数来描述雨滴撞击壁面水膜后的状态，结合已有的雨滴溅射数据，采用数值拟合技术，建立溅射模型，模拟雨滴溅射。在模拟时，相关长度参数均以雨滴半径进行无量纲化。

求解雨滴溅射模型的四个参数，具体包括以下步骤：

①求解无量纲溅射凹坑高度的理论值

假定运动中的雨滴仅受惯性力、重力和表面张力的作用，基于能量守恒，雨滴碰撞前的总能量等于水膜形态改变为凹坑所需的表面能，可得到无量纲溅射凹坑高度

的理论值为：

式中，W_b为雨滴韦伯数；F_N为雨滴弗洛德数；h^*为无量纲水膜厚度；

为无量纲溅射凹坑半径；

②求解无量纲溅射凹坑高度的实际修正值

充分考虑能量耗散，通过数值拟合技术拟合试验数据，对理论值进行改进，以提高与实际符合的程度，并拓展适用范围，得到表达式：

③间接求解无量纲溅射凹坑半径

试验表明，无量纲溅射凹坑半径

与无量纲溅射凹坑高度

和雨滴韦伯数W_b相关，本发明通过拟合试验数据，得到表达式：

④求解无量纲溅射凹坑形成时间

试验表明，无量纲溅射凹坑形成时间

与无量纲溅射凹坑高度

相关，本发明充分考虑能量耗散，通过拟合试验数据，得到表达式：

⑤求解无量纲凹坑间隔距离

保守估计凹坑存在时间ΔT为凹坑形成时间t_c的两倍。则对于ΔT时间的持续降雨，直径为D的所有雨滴在壁面单位面积内形成的凹坑总数为：

由于采用经典M-P雨滴谱模拟降雨，则降雨在壁面单位面积内形成的凹坑总数为：

于是，壁面单位面积内的无量纲平均凹坑间距

为：

式中，β为飞行器的局部雨滴收集率。

19)求解等效的水膜层粗糙度。包括以下步骤：

a、求解雨滴溅射引起的水膜表面粗糙度k_s,c；

本发明基于上述改进的雨滴溅射模型，采用等效沙粒粗糙度模拟雨滴撞击飞行器表面引起的水膜表面粗糙度状态，表达式为：

式中，

系数C与撞击凹坑形状有关，雨滴撞击凹坑多为圆柱形凹坑，一般取C＝0.64；

为凹坑当量平均高度，采用M-P雨滴谱分布模拟降雨时有：

b、求解水膜层表面水波纹引起的表面粗糙度k_s,w；

试验表明，水膜层表面水波纹引起的表面粗糙度与水膜厚度相关，大致与水膜厚度成正比：

k_s,w＝B·h (36)

式中，根据试验数据，保守取B＝1.5。

c、求解雨滴撞击与溅射引起的等效的水膜层表面粗糙度k_s。

最终，取等效的水膜层表面粗糙度，表达式为：

k_s＝||[k_s,c k_s,w]|| (37)

20)根据步骤17)得到的水膜层厚度更新飞行器几何形状，根据带有附着水膜界面的、新的飞行器几何形状重构计算网格。

21)引入步骤19)基于雨滴溅射模块得到的水膜层表面粗糙度，展开空气相流体流动模块的更新求解，获得更新空气相流体气动性能数据。

22)将步骤21)得到的更新空气相流体气动性能数据与步骤4)得到的空气相流体气动性能数据相减，再和步骤6)得到的雨滴撞击力系数增量数据叠加，获得降雨引起的气动性能损失数据。

下面以采用Fortran90计算机高级程序语言实现、通过计算机模拟NACA0012机翼在空中飞行期间遭遇降雨时的状态及由此引起的气动性能损失情况作为具体实施例，进一步说明本发明的一种飞行器飞行中遭遇降雨时的数值模拟方法，包括以下步骤：

1)根据目标飞行器的几何外形，生成飞行器流场空间计算网格。

在NACA0012翼型周围生成计算网格，这是一个二维非结构化与结构化混合网格，如图7所示，沿翼型壁面法向生成3层结构化网格，壁面首层网格厚度为0.5mm，每层围绕翼型生成4080个网格，网格增长率为1.2，其它非结构化网格约11万。

2)确定数值模拟的飞行条件为：飞行高度0m，来流速度65m/s，大气压力101325Pa，迎角0°、8°、14°、16°，翼型弦长10m，空气密度1.225kg/m³，空气粘性系数1.7894×10⁵Pa·s；按照本发明关于降雨条件的方法一，确定数值模拟的降雨条件为：降雨率2000mm/h，雨水密度1000kg/m³。

3)进行空气相-雨滴相两相流体流动模拟。空气相流体流动控制方程采用技术成熟的非定常、不可压缩气体的雷诺平均N-S方程，空气相流体湍流采用高雷诺数S-A湍流模型模拟；雨滴相流体流动控制方程采用式(6)和式(7)；求解算法采用技术成熟的一阶隐式非定常求解算法。

4)根据流场数据获得空气相流体气动性能数据和空气相流体壁面剪切应力等相关流场数据。

5)根据式(8)求解局部雨滴收集率，根据式(10)求解局部雨滴撞击质量流量。

6)根据式(13)求解雨滴撞击对气动性能的影响，以系数增量形式表示。

7)采用QUICK格式离散步骤4)得到的空气相流体壁面剪切应力，得到各表面水膜控制体边界处的空气相流体壁面剪切应力。

8)设定翼型表面同时存在气流附着点、气流转捩点和气流再附点，令i₀＝N_max/2，i₁＝N_max+10，i₂＝N_max+20；初始化i₀控制体、i₁控制体和i₂控制体的位置。

9)根据式(22)求解i₀控制体水膜厚度。

10)根据式(23)求解i₀控制体左侧界面方向所有控制体的水膜厚度。

11)检测i₁控制体的存在性，若存在，更新其位置；否则，不予更新。

12)根据式(24)求解i₀控制体与i₁控制体之间的控制体的水膜厚度。

13)检测i₂控制体的存在性，若存在，更新其位置；否则，不予更新。

14)根据式(22)的方法求解i₂控制体水膜厚度。

15)根据式(25)求解i₁控制体与i₂控制体之间控制体的水膜厚度。

16)根据式(26)求解i₂控制体右侧边界方向所有控制体的水膜厚度。

17)整理计算得到的所有水膜厚度。

18)根据式(28)、(29)、(30)、(31)求解雨滴溅射模型的四个参数。

19)根据式(32)求解雨滴撞击壁面引起的水膜表面粗糙度；根据式(34)求解水膜表面水波纹引起的表面粗糙度；根据式(35)求解等效的水膜层粗糙度。

20)根据步骤17)得到的水膜层厚度更新飞行器几何形状，根据带有附着水膜界面的、新的飞行器几何形状重构计算网格。

21)引入步骤19)基于雨滴溅射模块得到的水膜层表面粗糙度，展开空气相流体流动模块的更新求解，获得更新空气相流体气动性能数据。

22)将步骤21)得到的更新空气相流体气动性能数据与步骤4)得到的空气相流体气动性能数据相减，再和步骤6)得到的雨滴撞击力系数增量数据叠加，获得降雨引起的气动性能损失数据。

如图8所示，从上到下分别是迎角0°、8°、14°、16°时翼型表面附着水膜层状态，由于雨膜层很薄，为便于观察与研究，将雨膜厚度扩大10倍进行显示。结果表明，飞行器在执行飞行任务过程中遭遇降雨时，雨水撞击飞行器表面，不但在其表面形成一层薄薄的水膜，而且翼型上表面水膜厚度大于翼型下表面水膜厚度。飞行器飞行中，当在小迎角时，雨膜层状态变化不大；随着迎角的继续增大，翼型后缘上表面雨膜厚度增长明显；当迎角达到16°时，由于翼型开始失速，在翼型上表面(转捩位置附近)由于雨水被转捩气流带起，还会促使此处雨膜厚度大大增加，而翼型上表面后段受气流逆压梯度的影响，水膜厚度几乎为零，从而导致翼型气动外形严重受到影响。飞行中飞行器遭遇降雨时在翼型表面所形成的这种水膜层特性必然致使翼型气动性能遭到损失，威胁飞行器的飞行安全。如图9(a)、(b)、(c)所示，本发明考虑能量耗散，将利用数值拟合技术得到的式(28)、式(30)和式(34)以及线性关系式(29)与试验数据进行对比，可以看出，数值拟合与试验数据符合度高，拓展合理，易于程序实现。如图10和图11所示，给出了净翼型与降雨条件下的翼型的气动性能数据对比。

上述各实施例仅用于说明本发明，其中各部件的结构、设置位置及其连接方式等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (7)

1.一种飞行器飞行中遭遇降雨时的数值模拟方法，包括以下步骤：

1)根据目标飞行器的几何外形，生成飞行器流场空间计算网格；

2)根据目标飞行器的实际使用条件，确定飞行条件和降雨条件；

3)根据飞行器流场空间计算网格、飞行条件和降雨条件，进行空气相-雨滴相两相流体流动模拟计算，获得流场数据；

4)根据流场数据获得空气相流体气动性能数据和空气相流体壁面剪切应力；

5)利用动量定理建立模拟雨滴对飞行器表面撞击效应的计算式，求解飞行器表面局部雨滴收集率和飞行器表面局部雨滴撞击质量流量；

6)根据飞行器表面局部雨滴撞击质量流量，求解雨滴撞击力系数增量；

7)取飞行器流场空间计算网格中的翼型表面贴体网格作为表面水膜控制体，采用QUICK格式离散步骤4)得到的空气相流体壁面剪切应力，得到各表面水膜控制体边界处的空气相流体壁面剪切应力；

8)设定翼型表面同时存在气流附着点、气流转捩点和气流再附点，分别记气流附着点、气流转捩点和气流再附点所在的表面水膜控制体为i₀控制体、i₁控制体和i₂控制体，并初始化i₀控制体、i₁控制体和i₂控制体的位置；

9)求解i₀控制体水膜厚度；

10)求解i₀控制体左侧界面方向所有表面水膜控制体的水膜厚度；

11)检测i₁控制体的存在性，若存在，更新其位置；否则，不予更新；

12)求解i₀控制体与i₁控制体之间的表面水膜控制体的水膜厚度；

13)检测i₂控制体的存在性，若存在，更新其位置；否则，不予更新；

14)求解i₂控制体水膜厚度；

15)求解i₁控制体与i₂控制体之间表面水膜控制体的水膜厚度；

16)求解i₂控制体右侧边界方向所有表面水膜控制体的水膜厚度；

17)整理计算得到的所有表面水膜控制体的水膜厚度；

18)建立雨滴溅射模型，求解雨滴溅射模型的四个参数：无量纲溅射凹坑高度实际修正值、无量纲溅射凹坑半径、无量纲溅射凹坑形成时间和无量纲凹坑间隔距离，具体包括以下步骤：

①求解无量纲溅射凹坑高度的理论值

为：

式中，W_b为雨滴韦伯数；F_N为雨滴弗洛德数；h^*为无量纲水膜厚度；

为无量纲溅射凹坑半径；

②求解无量纲溅射凹坑高度的实际修正值

为：

③间接求解无量纲溅射凹坑半径

为：

④求解无量纲溅射凹坑形成时间

为：

⑤求解无量纲凹坑间隔距离

为：

式中，β为飞行器的局部雨滴收集率；I＝nR^-m，R为降雨率，m和n为常数，均与降雨类型相关，为降雨试验测量值；N₀为单位体积空间中的雨滴数密度随雨滴直径大小的分布常数；MVD为雨滴体积当量平均直径；

19)根据雨滴溅射模型的四个参数，求解等效的水膜层表面粗糙度，具体包括以下步骤：

a、采用等效沙粒粗糙度模拟雨滴撞击飞行器表面引起的水膜表面粗糙度状态，求解雨滴溅射引起的水膜表面粗糙度k_s,c为：

式中，

系数C与撞击凹坑形状有关；D_c为凹坑间隔距离，

为凹坑当量平均高度，采用M-P雨滴谱分布模拟降雨时有

b、求解水膜层表面水波纹引起的表面粗糙度k_s,w：

k_s，w＝B·h

式中，取B＝1.5；h为水膜厚度，h＝MVD·h^*/2；

c、求解雨滴撞击与溅射引起的等效的水膜层表面粗糙度k_s为：

k_s＝||[k_s,c k_s,w]||；

20)根据步骤17)得到的所有表面水膜控制体的水膜层厚度更新飞行器几何外形，根据带有附着水膜界面的、新的飞行器几何外形重构飞行器流场空间计算网格；

21)引入步骤19)得到的等效的水膜层表面粗糙度，根据重构的飞行器流场空间计算网格，以及飞行条件和降雨条件，再次进行空气相-雨滴相两相流体流动模拟计算，获得更新流场数据，根据更新流场数据获得更新空气相流体气动性能数据；

22)将步骤21)得到的更新空气相流体气动性能数据与步骤4)得到的空气相流体气动性能数据相减，再和步骤6)得到的雨滴撞击力系数增量数据叠加，获得降雨引起的飞行器气动性能损失数据。

2.如权利要求1所述的一种飞行器飞行中遭遇降雨时的数值模拟方法，其特征在于，所述步骤2)中的飞行条件包括飞行器飞行高度H、飞行高度H处的空气密度ρ_H；降雨条件包括液态水含量LWC、雨滴体积当量平均直径MVD和雨滴下落的末速V_H；

其中，确定降雨条件的具体方法为：先定义降雨类型和降雨率R，计算液态水含量LWC和雨滴体积当量平均直径MVD；根据飞行器飞行高度H，计算雨滴下落的末速V_H；

或者，直接定义液态水含量LWC和雨滴体积当量平均直径MVD；根据飞行器飞行高度H，计算雨滴下落的末速V_H。

3.如权利要求1或2所述的一种飞行器飞行中遭遇降雨时的数值模拟方法，其特征在于，所述步骤3)中采用非定常、不可压缩气体的雷诺平均N-S方程作为空气相流体流动控制方程，采用高雷诺数S-A湍流模型模拟空气相流体湍流，由单个雨滴运动方程推导得到雨滴相流体流动控制方程，采用一阶隐式非定常求解算法进行空气相-雨滴相两相流体流动模拟计算。

4.如权利要求2所述的一种飞行器飞行中遭遇降雨时的数值模拟方法，其特征在于，所述步骤5)中飞行器表面局部雨滴收集率的计算公式为：

式中，i为翼型表面贴体网格计数，以顺时针方向按序计数，翼型下表面后缘贴体网格记为第1个网格，翼型上表面后缘贴体网格记为第N_max个网格，N_max为翼型表面贴体网格总数目；β_i为第i个贴体网格处的局部雨滴收集率，α_V0＝LWC/ρ_w；α_VS,i和V_⊥S,i分别为第i个贴体网格的雨滴相流体体积分数和雨滴法向撞击速度；V_∞为雨滴相流体的初始速度，V_a为空气相流体的来流速度；ρ_w为雨水密度；

飞行器表面局部雨滴撞击质量流量的计算公式为：

式中，

为第i个贴体网格处的局部雨滴撞击质量流量；A_p,i为第i个贴体网格的壁面表面积。

5.如权利要求4所述的一种飞行器飞行中遭遇降雨时的数值模拟方法，其特征在于，所述步骤6)中雨滴撞击力系数增量的计算公式为：

式中，△C_L,Impinging和△C_D,Impinging分别是升力和阻力系数增量；F_i为第i个贴体网格处雨滴对飞行器表面局部的撞击力矢量；F为雨滴对飞行器表面总的撞击力矢量；F_x和F_y分别为撞击力矢量F沿计算网格坐标系x轴和y轴的分量；α为来流迎角。

6.如权利要求4或5所述的一种飞行器飞行中遭遇降雨时的数值模拟方法，其特征在于，所述步骤9)中，取第i个表面水膜控制体的空气相流体壁面剪切应力τ_i的符号与前一表面水膜控制体的空气相流体壁面剪切应力τ_i-1的符号刚好不同时的表面水膜控制体为i₀控制体，然后进行一次修正；i₀控制体的水膜厚度

的计算公式为：

式中，

和

分别为i₀控制体左右两个界面处的空气相流体壁面剪切应力，通过QUICK格式离散得到；

为i₀控制体处的局部雨滴撞击质量流量；μ_w为水的粘性系数；

所述步骤11)中，取第i个表面水膜控制体的空气相流体壁面剪切应力τ_i的符号与前一表面水膜控制体的空气相流体壁面剪切应力τ_i-1的符号刚好不同时的表面水膜控制体为i₁控制体，然后进行一次修正；若存在气流转捩点，则更新i₁，否则i₁值不变；

所述步骤13)中，取第i个表面水膜控制体的空气相流体壁面剪切应力τ_i的符号与前一表面水膜控制体的空气相流体壁面剪切应力τ_i-1的符号刚好不同时的表面水膜控制体为i₂控制体，然后进行一次修正；若存在气流再附点，则更新i₂，否则i₂值不变；

所述步骤14)中，i₂控制体的水膜厚度

的计算公式为：

式中，

和

分别为i₂控制体左右两个界面处的空气相流体壁面剪切应力，通过QUICK格式离散得到；

为i₂控制体处的局部雨滴撞击质量流量。

7.如权利要求6所述的一种飞行器飞行中遭遇降雨时的数值模拟方法，其特征在于，所述步骤10)中，从i₀控制体开始，沿左侧界面方向依次求解第i个表面水膜控制体左侧界面处的质量流量

然后依次求解i₀控制体左侧界面方向所有控制体的水膜厚度

其中，任意第i个表面水膜控制体左右两个界面处的质量流量计算公式为：

任意第i个表面水膜控制体左侧界面处的水膜厚度的计算公式为：

所述步骤12)中，从i₀控制体开始，沿右侧界面方向依次求解

然后依次求解i₀控制体与i₁控制体之间的表面水膜控制体的水膜厚度；i₀控制体与i₁控制体之间表面水膜控制体的水膜厚度的计算公式为：

所述步骤15)中，i₁控制体与i₂控制体之间表面水膜控制体的水膜厚度的计算公式为：

所述步骤16)中，从i₂控制体开始，依次展开计算i₂控制体右侧边界方向所有表面水膜控制体的水膜厚度；i₂控制体右侧边界方向所有表面水膜控制体的水膜厚度的计算公式为：

上述各式中，

和

分别为第i个表面水膜控制体左右两个界面处的空气相流体壁面剪切应力，通过QUICK格式离散得到；

和

分别为第i个表面水膜控制体左右两个界面处的水膜厚度；

为第i个表面水膜控制体处的局部雨滴撞击质量流量。

Images