迎风格式在水滴欧拉方程数值求解中的应用
迎风格式在水滴欧拉方程数值求解中的应用
Application of Upwind Scheme in Numerical Solution of Euler Equation for Water Droplets
作者: 宁义君中国航空工业空气动力研究院 辽宁省飞行器防除冰重点实验室，辽宁 沈阳 110034 顾洪宇中国航空工业空气动力研究院 辽宁省飞行器防除冰重点实验室，辽宁 沈阳 110034 朱东宇中国航空工业空气动力研究院 辽宁省飞行器防除冰重点实验室，辽宁 沈阳 110034 乔龙中国航空工业空气动力研究院 沈阳市计算流体力学重点实验室，辽宁 沈阳 110034 李艳亮中国航空工业空气动力研究院 沈阳市计算流体力学重点实验室，辽宁 沈阳 110034
Author: Ning YijunLiaoning Provincial Key Laboratory of Aircraft Ice Protection，AVIC Aerodynamics Research Institute，Shenyang 110034，China Gu HongyuLiaoning Provincial Key Laboratory of Aircraft Ice Protection，AVIC Aerodynamics Research Institute，Shenyang 110034，China Zhu DongyuLiaoning Provincial Key Laboratory of Aircraft Ice Protection，AVIC Aerodynamics Research Institute，Shenyang 110034，China Qiao LongShenyang Key Laboratory of Computational Fluid Dynamics，AVIC Aerodynamics Research Institute，Shenyang 110034，China Li YanliangShenyang Key Laboratory of Computational Fluid Dynamics，AVIC Aerodynamics Research Institute，Shenyang 110034，China
摘要：数值模拟是研究飞机结冰的重要手段，为获取水滴撞击特性，通常采用欧拉法求解水滴运动方程。液态水含量在空间分布上存在间断，会导致水滴欧拉方程的求解极易发散；若采用一阶格式，虽然稳定性得到提高，但计算精度难以保证满足工程设计需求；有研究者通过添加扩散项以提升收敛性，但会人为增加数值耗散。针对这一问题，本文采用一阶迎风格式对水滴欧拉方程的对流通量进行离散，并通过线性重构与限制器函数保证计算精度达到二阶，建立一种基于迎风离散格式的水滴欧拉方程二阶求解方法。通过对比NACA0012、NACA23012和NLF0414翼型水滴撞击特性计算结果与FENSAP-ICE软件的结果，验证了该方法的有效性与计算精度。
Abstract：Numerical simulation is an important means to study aircraft icing.In order to obtain water droplet impingement characteristics,Euler method is usually used to solve the control equation of droplet movement.Due to the discontinuity in the spatial distribution of liquid water content,the solution of the Euler equation of water droplet is easily divergent.If the first-order format is adopted,although the stability is improved,the calculation accuracy is difficult to meet the engineering design requirements.Some researchers increase convergence by adding diffusion terms,but artificially increase numerical dissipation.To solve this problem,the convective flux of the water droplet Euler equation is discretized by a first-order upwind scheme,and the calculation accuracy is guaranteed to reach the second-order by linear reconstruction and limiter functions.The effectiveness and accuracy of the method are verified by comparing the calculation results of NACA0012 NACA23012 and NLF0414 airfoils droplet impingement characteristics with FENSAP-ICE software.
关键词：水滴撞击特性;水滴运动轨迹;欧拉法;迎风格式;结冰软件
Keywords：water droplet impingement characteristics;water droplet trajectory;Euler method;upwind scheme;icing software
数值模拟是研究飞机结冰的重要手段，其中水滴撞击特性是结冰预测以及防/除冰系统设计的必备输入条件。水滴撞击特性计算即通过求解水滴运动方程，模拟水滴运动轨迹，并确定水滴撞击物面的区域以及撞击区域内的水滴收集系数。
对于水滴运动轨迹的数值模拟，最先发展起来的是拉格朗日法[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]和Venkatakrish‐nan[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；S up为上翼面水滴撞击极限误差，(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
在两种限制器中，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 不同限制器水滴收集系数特征参数对比Table 3 Comparison between collection coefficient characteristic parameters for different limiters
图2 一阶迎风LWC云图Fig.2 LWC cloud of first order upwind
根据以上分析，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两种限制器对水滴收集系数的计算不受水滴直径和风速等计算状态的影响。
图3 不同限制器LWC云图对比Fig.3 Comparison between LWC cloud of different limiters
图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%。
采用NLF0414翼型，弦长为0.9144m，基于表2计算状态，MinMod和Venkatakrishnan两种限制器计算的水滴收集系数与FENSAP-ICE对比如图7所示，特征参数对比见表7，水滴收集系数最大值βmax误差分别为-2.35%和-2.29%。
图6 NACA23012翼型水滴收集系数对比Fig.6 Comparison between collection coefficient for NACA23012 airfoil
表6 NACA23012翼型水滴收集系数特征参数对比Table 6 Comparison between collection coefficient characteristic parameters for NACA23012 airfoil
图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软件。
参考文献
[1]Ruff G，Berkowitz B.Users manual for the NASA LEWIS ice accretion prediction code（LEWICE)[R].NASA CR-185129,1990.
[2]战培国.美国NASA结冰试验设备体系综述[J].航空科学技术,2021,32(5):1-6.Zhan Peiguo.Review on the system of icing facilities in NASA[J].Aeronautical Science & Technology，2021，32（5）：1-6.（in Chinese）
[3]陈科，曹义华，潘星.改进的翼型积冰数值模拟方法[J].航空动力学报，2007，22（11）：1814-1819.Chen Ke，Cao Yihua，Pan Xing.An improved numerical simulation method for airfoil ice accretion[J].Journal of Aerospace Power，2007，22（11）：1814-1819.（in Chinese）
[4]周志宏，易贤，桂业伟，等.水滴撞击特性的高效计算方法[J].空气动力学学报，2014，32（5）：712-716.Zhou Zhihong，Yi Xian，Gui Yewei，et al.An efficient method to simulate water droplet trajectory and impingement[J].Acta Aerodynamica Sinica，2014，32（5）：712-716.（in Chinese）
[5]Bourgault Y，Boutanios Z，Habashi W G.Three-dimensional Eulerian approach to droplet impingement simulation using FENSAP-ICE，part 1：Model algorithm，and validation[J].Journal of Aircraft，2000，37（1）：95-103.
[6]Luliano E，Brandi V，Mingione G.Water impingement predic‐tion on multi-element airfoils by means of Eulerian and La‐grangian approach with viscous and inviscid air flow[R].AIAA Paper 2006-1270，2006.
[7]王洪伟，李先哲，宋展.通用飞机结冰适航验证关键技术及工程应用[J].航空学报，2016，37（1）：335-350.Wang Hongwei，Li Xianzhe，Song Zhan.Key airworthiness validation technologies for icing of general aviation aircraft and their engineering application[J].Acta Aeronautica et Astronautica Sinica，2016，37（1）：335-350.（in Chinese）
[8]Durst F，Miloievic D，Schönung B.Eulerian and Lagrangian predictions of particulate two-phase flows：A numerical study[J].Applied Mathematical Modelling，1984，8（2）：101-115.
[9]Silva D M D，Bachchan N，Lim K I，et al.Icing collection efficiency prediction using an Eulerian-Eulerian approach[R].AIAA Paper 2014-3073，2014.
[10]Jung S K，Myong R S.A second-order positivity-preserving finite volume upwind scheme for air-mixed droplet flow in atmospheric icing[J].Computers & Fluids，2013，86：459-469.
[11]Kim J W，Dennis P G，Sankar L N，et al.Ice accretion modeling using an Eulerian approach for droplet impingement[R].AIAA Paper 2013-0246，2013.
[12]Bourgault Y，Habashi W G，Dompierre J.An Eulerian approach to supercooled droplets impingement calculations[R].AIAA Paper 1997-0176，1997.
[13]Bourgault Y，Habashi W G，Dopierre J，et al.A finite element method study of Eulerian droplets impingement models[J].International Journal for Numerical Methods in Fluids，1999，29：429-449.
[14]Bourgault Y，Boutanios Z，Habashi W G.Three-dimensional Eulerian approach to droplet impingement simulation using FENSAP-ICE，part 1：Model algorithm and validation[J].Journal of Aircraft，2000，37（1）：95-103.
[15]Morency F，Beaugendre H，Baruzzi G S，et al.FENSAPICE：A comprehensive 3D simulation system for in-flight icing[R].AIAA Paper 2001-2566，2001.
[16]Luliano E.An Eulerian method for the evaluation of ice accretion on airplanes[D].Naples：The University of Naples，2004.
[17]Jung S K，Myong R S，Cho T H.Development of Eulerian droplets impingement model using HLLC Riemann solver and POD-based reduced order method[R].AIAA Paper 2011-3907，2011.
[18]Jung K Y，Jung S K，Myong R S，et al.A three-dimensional unstructured finite volume method for droplet impingement in aircraft icing[R].AIAA Paper 2013-2576，2013.
[19]卜雪琴，林贵平.基于CFD的水收集系数及防冰表面温度预测[J].北京航空航天大学学报，2007，33（10）：1182-1185.Bu Xueqin，Lin Guiping.Predictions of collection efficiency and anti-icing surface temperature[J].Journal of Beijing University of Aeronautics and Astronautics，2007，33（10）：1182-1185.（in Chinese）
[20]张大林，杨曦，昂海松.过冷水滴撞击结冰表面的数值模拟[J].航空动力学报，2003，18（1）：87-91.Zhang Dalin，Yang Xi，Ang Haisong.Numerical simulation of supercooled water droplets impingement on icing surfaces[J].Journal of Aerospace Power，2003，18（1）：87-91.（in Chinese）
[21]陈维建，张大林.瘤状冰结冰过程的数值模拟[J].航空动力学报，2005，20（3）：472-476.Chen Weijian，Zhang Dalin.Numerical simulation of glaze ice accretion process[J].Journal of Aerospace Power，2005，20（3）：472-476.（in Chinese）
[22]张强，曹义华，李栋.采用欧拉两相流法对翼型表面霜冰的数值模拟[J].北京航空航天大学学报，2009，35（3）：351-355.Zhang Qiang，Cao Yihua，Li Dong.Numerical method to simulate rime ice accretions on an airfoil[J].Journal of Beijing University of Aeronautics and Astronautics，2009，35（3）：351-355.（in Chinese）
[23]张强，曹义华，钟国.飞机机翼表面霜冰的三维数值模拟[J].航空动力学报，2010，25（6）：1303-1309.Zhang Qiang，Cao Yihua，Zhong Guo.Three-dimensional numerical simulation of rime ice accretions on an aircraft wing[J].Journal of Aerospace Power，2010，25（6）：1303-1309.（in Chinese）
[24]易贤，王开春，桂业伟，等.结冰面水滴收集率欧拉计算方法研究及应用[J].空气动力学学报，2010，28（5）：596-601.Yi Xian，Wang Kaichun，Gui Yewei，et al.Study on Eulerian method for icing collection efficiency computation and its application[J].Acta Aerodynamica Sinica，2010，28（5）：596601.（in Chinese）
[25]周志宏，李凤蔚，李广宁.基于两相流欧拉方法的翼型结冰数值模拟[J].西北工业大学学报，2010，28（1）：138-142.Zhou Zhihong，Li Fengwei，Li Guangning.Applying Eulerian droplet impingement model to numerically simulating ice accre‐tion but with some improvements[J].Journal of Northwestern Polytechnical University，2010，28（1）：138-142.（in Chinese）
[26]周志宏，易贤，桂业伟，等.复杂三维外形霜状冰数值模拟[J].四川大学学报（工程科学版），2012（S2）：158-162.Zhou Zhihong，Yi Xian，Gui Yewei，et al.Prediction of rime ice accretion on complex configurations[J].Journal of Sichuan University（Engineering Science Edition），2012（S2）：158162.（in Chinese）
[27]赵文朝，宁义君，吴渊，等.飞机结冰中水滴撞击特性欧拉法的网格影响分析[J].南京航空航天大学学报（英文版），2023，40（2）：148-158.Zhao Wenchao，Ning Yijun，Wu Yuan，et al.Mesh impact analysis of Eulerian method for droplet impingement character‐istics under aircraft icing conditions[J].Transactions of Nan‐jing University of Aeronautics and Astronautics，2023，40（2）：148-158.（in Chinese）
[28]Tong X，Luke E A.Robust and accurate Eulerian multiphase simulations of icing collection efficiency using singularity dif‐fusion model[J].Engineering Applications of Computational Fluid Mechanics，2010，4（4）：483-495.
[29]申晓斌，林贵平，杨胜华.三维发动机进气道水滴撞击特性分析[J].北京航空航天大学学报，2011，37（1）：1-5.Shen Xiaobin，Lin Guiping，Yang Shenghua.Analysis on threedimensional water droplets impingement characteristics of engine inlet[J].Journal of Beijing University of Aeronautics and Astronautics，2011，37（1）：1-5.（in Chinese）
[30]Blazek J.Computational fluid dynamics principles and applica‐tions[M].London：Elsevier Science Ltd.，2001.
[31]SAE.Droplet impingement and ice accretion computer codes：ARP5903[S].US：SAE International，2009.
[32]Boutanios Z.An Eulerian 3D analysis of water droplets im‐pingement on a convair-580 nose and cockpit geometry[D].Montreal：Concordia University，1999.
[33]Langmuir I，Blodgett K B.A mathematical investigation of water droplet trajectories[R].US Army Air Forces Technical Report 5418，1946.
[34]Barth T J，Jespersen D C.The design and application of up‐wind schemes on unstructured meshes[R].AIAA Paper 1989 0366，1989.
[35]Venkatakrishnan V.On the accuracy of limiters and convergence to steady state solutions[R].AIAA Paper 19930880，1993.
[36]Venkatakrishnan V.Convergence to steady-state solutions of the Euler equations on unstructured grids with limiters[J].Journal of Computational Physics，1995，30：118-120.
[37]乔龙，李艳亮，杨思源，等.基于面向对象的非结构航空CFD软件体系结构设计[J].航空科学技术，2022，33（7）：66-72.Qiao Long，Li Yanliang，Yang Siyuan，et al.Software archi‐tecture of an unstructured aviation CFD solver using object-ori‐ented techniques[J].Aeronautical Science & Technology，2022，33（7）：66-72.（in Chinese）