迎风格式在水滴欧拉方程数值求解中的应用

宁义君，顾洪宇，朱东宇，乔龙，李艳亮

1.中国航空工业空气动力研究院 辽宁省飞行器防除冰重点实验室， 辽宁 沈阳 110034

2.中国航空工业空气动力研究院 沈阳市计算流体力学重点实验室， 辽宁 沈阳 110034

数值模拟是研究飞机结冰的重要手段，其中水滴撞击特性是结冰预测以及防/除冰系统设计的必备输入条件。水滴撞击特性计算即通过求解水滴运动方程，模拟水滴运动轨迹，并确定水滴撞击物面的区域以及撞击区域内的水滴收集系数。

对于水滴运动轨迹的数值模拟，最先发展起来的是拉格朗日法[1-2]，国内早期也普遍采用这种方法[3-4]。但模拟复杂三维构型时，该方法释放粒子数、网格搜索和空间插值计算量都大幅增加，且难以精确捕捉几何表面边界并计算水滴撞击特性[5]。而欧拉法将水滴看作空间连续分布的液态相，以水滴经过的空间位置为研究对象建立水滴运动控制方程，引入水滴场的概念，通过设置物面边界条件获取水滴收集系数。相比拉格朗日法，欧拉法在计算复杂外形的水滴撞击特性时更具明显优势。有研究者对比了欧拉法与拉格朗日法的计算结果，表明水滴收集系数分布基本一致[6]。通常认为这两种方法的计算结果是等价的，且都被飞机适航分析人员所接受[7]。

受益于计算流体力学的发展，欧拉法于20世纪80年代从两相流的思想发展而来[8]，此后越来越多的研究者开始采用欧拉法计算水滴撞击特性[9]。当采用欧拉法时，若求解格式不稳定，极易导致非物理解出现[10]。为延缓计算崩溃的出现，很多计算代码都采用一阶格式[11]。Y.Bourgault等[12-13]基于有限单元法实现了水滴欧拉方程二阶精度求解，并将该方法引入FENSAP-ICE[14-15]软件当中（现已被ANSYS 公司收购）。而目前对水滴欧拉方程的求解，绝大多数仍然采用有限体积法[16]。S.K.Jung等[10,17-18]提出了一种水滴欧拉方程改造方法，将其由退化双曲型转化为完全双曲型，并基于有限体积法开发了一种二阶保正性迎风离散格式。

国内在这方面早期以FLUENT等商业软件用户定义函数研究为主[19]。随着数值开发工作的深入，开始编制水滴欧拉方程的数值求解代码。对于对流通量的空间离散，最早出现MUSCL[20-21]、QUICK[22-23]等格式，随后迎风与线性插值相结合的方法[24]、VanLeer 迎风插值方法[25-26]等被提出。近期赵文朝等[27]采用欧拉法对NACA0012 翼型、冰风洞风道、S 形进气道及三段翼的水滴撞击特性进行了网格影响分析，表明当水滴运动没有受到上游部件的导流或者遮挡时，欧拉法计算结果具有网格无关性。

在对流通量离散格式上，一阶格式具备较高的稳定性，但数值耗散大，计算精度较低；而单纯的高阶格式容易引起数值振荡和密度脉冲[10]。有研究者通过添加扩散项以提高稳定性，但人为增加了数值耗散[28-29]。一阶迎风叠加线性重构与限制器函数的二阶离散格式已广泛应用于空气欧拉方程的求解[30]，但鲜有在水滴欧拉方程上的应用。

本文发展了一种基于迎风离散格式的具备二阶精度的水滴欧拉方程数值求解方法，且无须添加扩散项，并已集成至航空工业气动院自主开发的航空飞行器结冰仿真软件穿云（ChuanYun）V1.0 中。通过对NACA0012 翼型水滴撞击特性的模拟，对一阶与二阶精度计算结果进行了分析，并根据SAE ARP5903[31]规范性要求，采用若干特征参数对穿云V1.0计算结果与FENSAP-ICE软件进行量化比较，以验证本文研究方法的有效性与计算精度，最后以相同计算方法应用于NACA23012翼型和NLF0414翼型。

1 数值方法

1.1 控制方程

将水滴看作分布在空气中的连续相，引入液态水含量（LWC）表示单位体积空气中水的质量。水滴在空气中的LWC 通常在10-3量级，因此在建立控制方程时一般忽略水滴对空气的作用，采用松耦合方式进行求解，即先求解空气场，然后基于空气场再求解水滴运动[32]。

为建立水滴运动的数学模型，对其物理过程进行适当的假设[8,32]：（1）水滴尺寸足够小，一般在微米（μm）量级，在空气中均匀分布；（2）水滴为球体，运动过程中不变形、不破碎、不碰撞或凝聚；（3）水滴在运动过程中不发生热交换、不蒸发，水滴温度、密度等物性参数不变；（4）作用于水滴的外力只有重力、空气阻力和浮力；（5）水滴与物面撞击时不反弹、不飞溅。

基于以上假设，根据质量守恒和动量守恒定律，水滴控制方程可表述为[32]

式中，ρ为水滴在空气中的LWC；ρa为空气密度；ρw为水密度；u为水滴速度矢量；ua为空气速度矢量，K为惯性因子；g为重力加速度矢量。ρa和ua都由空气场计算结果中直接提取。

惯性因子K为

式中，μa为空气动力黏度；d为水滴直径；f为阻力函数；CD为水滴阻力系数；Red为水滴雷诺数。

水滴雷诺数Red为

在计算水滴阻力时，采用经典圆球阻力系数计算公式[33]

采用有限体积法时，对于控制体Ω，水滴欧拉方程的积分形式为

式中，W为守恒量；F为对流通量；Q为源项；可表述如下

式中，ua，va，wa分别为空气速度ua的三个分量；u，v，w分别为水滴速度u的三个分量；gx，gy，gz分别为重力加速度g的三个分量。

V为逆变速度，计算公式为

式中,nx，ny，nz分别为面法向矢量n的三个分量。

对于二阶精度有限体积法，物理量在控制体内的分布是均匀的，将式（6）半离散化，表达形式为

式中，m为控制体面的序号；NF为控制体的面数；ΔS为控制体面的面积。

1.2 空间离散格式

基于通量矢量分裂（FVS）方法，将通量矢量F分解为非负和非正两部分，非负部分采用左侧值，用L 来表示；非正部分采用右侧值，用R来表示

由于水滴控制方程雅可比（Jacobi）系数矩阵有多重特征值V，而V为逆变速度，因此根据V与0 的大小关系定义分裂逆变速度，非负和非正逆变速度有

则通量矢量F的具体分裂方案为

1.3 重构与梯度

迎风格式需要确定控制体面的左侧值和右侧值。当假设每个控制体内的物理量为常数，并且左右值分别为相应控制体单元值时，决定了这种空间离散只具备一阶精度。一阶迎风格式虽然具备良好的稳定性和有界性，但会导致过度的数值扩散[30]。

若假设物理量在控制体内是变化的，则左右值插值方式决定了具体的空间精度。对于二阶精度，物理量在控制体内是线性的。为了计算左右值，必须对物理量进行重构。

线性重构是目前最流行的重构方法。假设物理量在控制体内是分段线性分布的，任意物理量ϕ左右值重构如下

式中，下标I和J分别为左右单元编号；r为单元体心到面心距离；Ψ为限制器函数；∇ϕ为物理量ϕ所在单元的梯度。根据格林-高斯（Green-Gauss）定理，任意物理量ϕ的梯度近似为ϕ与向外单位法向量乘积的控制体面积分，计算公式如下

式中，nIJ和SIJ分别为控制体面的单位法矢量和面积；N为与单元I相邻的单元数量。

1.4 限制器

二阶或高阶迎风格式在空间离散时需要使用限制器函数，以防止在大梯度区域（如水滴遮蔽区边界）产生数值振荡和非物理解。因此，所使用限制器至少能够保证单调性，即求解域中物理量的极大值不增加，极小值不减少，并且在时间推进过程中不会产生新的局部极值。

使用两种限制器函数：MinMod[34]和Venkatakrishnan[35-36]，见表1。

表1 限制器函数列表Table 1 List of limiter functions

各限制器函数中，r计算公式为

其中

式中，minJ和maxJ分别为单元I所有相邻单元J的最小值和最大值，ε为机器精度近似值。

单元限制器取面限制器的最小值

1.5 时间推进格式

将水滴欧拉方程的半离散式（8）简化为

式中，R为水滴欧拉方程的残差。

采用多步Runge-Kutta时间推进方法

式中，α1, …,αm为迭代系数；Δt为时间步长。

为加速收敛，采用当地时间步长，计算公式为

式中，CFL为库朗数。

1.6 边界条件

对于初始条件ρ∞= LWC∞，u∞=ua,∞，∞代表自由来流值。对于远场边界条件，如果是流入方向，边界值等于初始条件；如果是流出方向，边界值等于相邻单元值。

对于对称边界条件，速度无穿透，其值为ub-(ubn)n，ub为相邻单元速度矢量；其他物理量边界值等于相邻单元值。

对于物面边界条件，根据物面第一层网格水滴速度物面法向分量ubn，若小于0，代表指向物面内，边界值采用相邻单元值；若大于0，代表指向物面外，边界ρ值取接近于0的小量（如LWC∞×10-7），并代入相应的守恒量ρu中。

1.7 水滴收集系数

水滴撞击特性由水滴收集系数来表征，水滴收集系数计算公式为

式中，U∞为自由来流速度；Un为水滴物面法向速度。

2 计算结果与分析

空气流场基于航空工业气动院自主开发的航空飞行器气动力仿真软件飞廉（FeiLian）V1.0[37]进行计算，采用MENTER_SST 湍流模型、Roe 对流通量离散格式以及Venkatakrishnan限制器。

水滴撞击特性基于航空工业气动院自主开发的航空飞行器结冰仿真软件穿云（ChuanYun）V1.0进行计算，计算结果将与FENSAP-ICE软件进行对比。FENSAP-ICE软件已普遍应用于飞行器防除冰系统设计，具备工程应用可信度，与其对比可以验证本文研究方法的有效性与准确性。

根据SAE ARP5903[31]规范性要求，对比计算结果时采用的水滴收集系数特征参数有：βmax为水滴收集系数最大值误差，(βmax,C-βmax,F)/βmax,F；Sβmax为水滴收集系数最大值所在位置误差，(Sβmax,C-Sβmax,F)/Chord；Sup为上翼面水滴撞击极限误差，(Sup,C-Sup,F)/Chord；Sdw为下翼面水滴撞击极限误差，(Sdw,C-Sdw,F)/Chord，其中下标C 为穿云软件，下标F 代表FENSAP-ICE软件，Chord代表翼型弦长。

2.1 NACA0012翼型

NACA0012 翼型弦长为0.3048m，计算状态见表2，主要从不同限制器函数和不同状态参数两个方面进行对比分析。

表2 NACA0012翼型计算状态Table 2 Calculation condition for NACA0012 airfoil

采用一阶迎风和表1所列两种限制器计算的水滴收集系数与FENSAP-ICE对比如图1所示，图1中横坐标为量纲一化弧长S/c（c为翼型弦长），以翼型前缘点为原点，特征参数对见表3。从表3中可以看出，一阶迎风和两种限制器计算的水滴撞击极限与FENSAP-ICE 一致，误差主要集中在驻点区域。一阶迎风在驻点处的水滴收集系数明显高出其他计算结果，水滴收集系数最大值误差βmax达到6.38%，且出现明显的数据波动，如图1 所示。一阶迎风驻点附近LWC 云图如图2 所示，LWC 在驻点附近出现集中，根据水滴收集系数计算公式（19），该现象导致了驻点处水滴收集系数突然增大。

图1 不同限制器水滴收集系数对比Fig.1 Comparison between collection coefficient for different limiters

图2 一阶迎风LWC云图Fig.2 LWC cloud of first order upwind

表3 不同限制器水滴收集系数特征参数对比Table 3 Comparison between collection coefficient characteristic parameters for different limiters

在两种限制器中，MinMod和Venkatakrishnan计算结果与FENSAP-ICE 基本一致，水滴收集系数最大值βmax误差分别为0.33%和0.34%。说明与一阶迎风相比，MinMod 和Venkatakrishnan两种限制器对水滴收集系数的改进作用是一致的。

与一阶迎风相比，使用线性重构与限制器的主要作用在于对LWC 空间分布间断面进行精确捕捉。一阶迎风和两种限制器、FENSAP-ICE计算的翼型尾迹附近LWC云图如图3所示。从图3可以看出，一阶迎风存在明显的过度耗散，其水滴遮蔽区尺寸明显小于其他限制器和FENSAPICE 计算结果。从翼型尾迹水滴遮蔽区的宽度来看，两种限制器计算的水滴遮蔽区比FENSAP-ICE 计算结果略小；从水滴遮蔽区边界附近的LWC强度来看，两种限制器计算的LWC峰值比FENSAP-ICE计算结果有所减小。

图3 不同限制器LWC云图对比Fig.3 Comparison between LWC cloud of different limiters

根据以上分析，MinMod和Venkatakrishnan两种限制器对水滴遮蔽区的捕捉精度明显高于一阶迎风，但水滴遮蔽区空间分布比FENSAP-ICE 略小，说明在数值耗散上比FENSAP-ICE 略大。并且两种限制器计算的水滴遮蔽区基本无差别，对水滴遮蔽区的改进作用是一致的。针对MinMod和Venkatakrishnan两种限制器，以表1状态为基准，不同水滴直径计算结果对比如图4所示，不同风速计算结果对比如图5所示，特征参数对比分别见表4和表5。从图4和图5、表4和表5可以看出，两种限制器计算的水滴收集系数是一致的，与FENSAP-ICE 软件相比，水滴收集系数最大值误差βmax不超过0.73%。这说明MinMod 和Venkatak rishnan两种限制器对水滴收集系数的计算不受水滴直径和风速等计算状态的影响。

图4 不同水滴直径水滴收集系数对比Fig.4 Comparison between collection coefficient for different droplet diameters

图5 不同风速水滴收集系数对比Fig.5 Comparison between collection coefficient for different velocities

表4 不同水滴直径水滴收集系数特征参数对比Table 4 Comparison between collection coefficient characteristic parameters for different droplet diameters

表5 不同风速水滴收集系数特征参数对比Table 5 Comparison between collection coefficient characteristic parameters for different velocities

2.2 NACA23012和NLF0414翼型

根据NACA0012 翼型计算结果，可以看出MinMod 和Venkatakrishnan两种限制器计算结果彼此之间几乎是一致的，与FENSAP-ICE 软件计算结果也基本一致。将相同计算方法应用于其他翼型以进一步验证该方法的计算精度。

采用NACA23012翼型，弦长为0.4572m，基于表2计算状态，MinMod和Venkatakrishnan两种限制器计算的水滴收集系数与FENSAP-ICE 对比如图6 所示，特征参数对比见表6，水滴收集系数最大值βmax误差分别为0.86% 和0.63%。

图6 NACA23012翼型水滴收集系数对比Fig.6 Comparison between collection coefficient for NACA23012 airfoil

表6 NACA23012翼型水滴收集系数特征参数对比Table 6 Comparison between collection coefficient characteristic parameters for NACA23012 airfoil

采用NLF0414 翼型，弦长为0.9144m，基于表2 计算状态，MinMod和Venkatakrishnan两种限制器计算的水滴收集系数与FENSAP-ICE 对比如图7 所示，特征参数对比见表7，水滴收集系数最大值βmax误差分别为-2.35%和-2.29%。

图7 NLF0414翼型水滴收集系数对比Fig.7 Comparison between collection coefficient for NLF0414 airfoil

表7 NLF0414翼型水滴收集系数特征参数对比Table 7 Comparison between collection coefficient characteristic parameters for NLF0414 airfoil

3 结束语

本文将一种基于逆变速度分裂的迎风格式应用于水滴欧拉方程对流通量的空间离散，通过计算与分析，可以得到如下结论：

（1） 该迎风格式只具备一阶精度，在求解水滴欧拉方程时，驻点处的水滴收集系数出现较大误差，通过叠加线性重构与限制器函数，对水滴欧拉方程的求解可以达到二阶精度。

（2） MinMod 和Venkatakrishnan 限制器计算结果彼此之间基本没有差别，与FENSAP-ICE软件相比，水滴收集系数基本一致。

（3） 由于一阶迎风格式的数值耗散较大，水滴遮蔽区空间分布远小于FENSAP-ICE 软件，通过上述二阶求解方法，水滴遮蔽区计算精度大幅提升，已接近于FENSAP-ICE软件。