基于UDS框架的水滴撞击特性数值计算方法

毛晓东，卜雪琴，赵国昌，王建明

(1.沈阳航空航天大学 航空航天工程学部(院)，沈阳 110136； 2.北京航空航天大学 航空科学与工程学院，北京 100191)



基于UDS框架的水滴撞击特性数值计算方法

毛晓东1，卜雪琴2，赵国昌1，王建明1

(1.沈阳航空航天大学 航空航天工程学部(院)，沈阳 110136； 2.北京航空航天大学 航空科学与工程学院，北京 100191)

摘要：为了预测二维/三维复杂表面在可压/不可压流动状态下的水滴撞击特性，提出了一种基于UDS框架的水滴撞击特性数值计算方法。水滴相和空气相认为单向耦合，空气相流场由N-S方程预先求得。水滴相控制方程基于FLUENT中UDS输运方程框架，自定义编程迭代求解。为了排除壁面处水滴，给出一种水滴撞击壁面边界条件处理方法。引入具有自适应性的人工数值扩散项以增加计算稳定性，避免局部水滴容积分数异常大导致解的奇异。验证计算了二维圆柱和三维球面，通过与公开文献中试验数据的对比，证明了本文计算方法的准确性和有效性，能很好的解决复杂表面在各种流动状态下的水滴撞击特性预测。

关键词：UDS框架；水滴撞击；结冰；二相流；数值计算

当飞机在某些特殊的气象条件下飞行时，飞机表面以及航空发动机相关部件表面，均有可能发生结冰现象。结冰将改变飞机绕流流场以及发动机进气道流场特性，影响飞机的操纵性和稳定性以及发动机正常工作特性，危害飞行安全。如果冰层脱落并吸入发动机进气道，可能直接导致发动机停机，最终发生机毁人亡的严重事故[1]。水滴撞击特性计算，是飞机及发动机结冰研究的首要内容，是结冰冰形预测及防除冰系统设计的前提和基础。

水滴撞击特性计算方法主要分为拉格朗日法和欧拉法。拉格朗日法把水滴看作离散相，在求解空气相绕流流场的基础上，采用有限差分法追踪水滴在流场中的运动轨迹，从而得到水滴撞击极限及局部水滴收集系数等参数[2-3]。该方法对二维或几何形状简单的表面比较简便，但对于三维或复杂外形，由于粒子释放位置较难确定，处理起来比较复杂。欧拉法把水滴看作连续相介质，通过求解水滴相的连续方程和动量方程，得到计算域各网格节点处的水滴容积分数和速度矢量。相对于拉格朗日法，欧拉法无需确定粒子的释放位置，特别适合如三维机翼[4-5]、多段翼[6]、发动机进气道[7]、发动机整流帽罩[8]和发动机整流支板[9]等处的三维复杂表面的水滴撞击特性计算。目前FLUENT商用软件中的欧拉模型求解器被广泛应用于水滴撞击特性的数值求解。欧拉模型是非常复杂的多相流模型，采用完全耦合，计算量大，收敛速度慢，容易发散。而最重要的是它不支持模拟可压缩流动所需的压力远场边界条件，因此直接影响到结冰模拟中马赫数超过0.3的情况。本文提出了一种基于FLUENT中UDS框架的水滴撞击特性计算方法，该方法计算快速、准确，并且普遍适用于可压缩流及不可压缩流。

1数学模型

1.1模型假设

对于空气相和过冷水滴相所组成的两项流动，空气相通过阻力、湍流效应等方式影响水滴相。但由于水滴相的微粒浓度很小，故认为水滴对空气相无影响，即单向耦合[7]。在建立水滴相控制方程前，先作如下假设：(1)水滴为大小统一的球形刚体，忽略变形和破裂；(2)水滴之间没有碰撞和融合，在撞击到壁面后不发生反弹和飞溅；(3)忽略水滴与空气间的传热和传质；(4)忽略水滴相动量方程中的粘性项和压力项；(5)仅考虑水滴所受到的阻力和重力，忽略非稳态力；(6)阻力认为是稳态的；(7)不考虑空气对水滴的湍流作用。

1.2控制方程

由于单向耦合，故空气流场仍由Navier-Stokes方程描述，且可以预先计算收敛，作为水滴相控制方程求解的已知条件。水滴的连续性方程和动量方程如式(1)和式(2)所示。

(1)

(2)

式中，ρ为水滴密度，α为水滴容积分数，u为水滴速度矢量，ua为空气速度矢量，F为作用于水滴上的除阻力外的外力，K是空气和水滴的动量交换系数，由式(3)定义。

(3)

其中，μa为空气动力粘度，dp为水滴直径，f是阻力函数，由于假设水滴为球形刚体，阻力函数如式(4)。

(4)

其中，CD为阻力系数，Re为相对雷诺数，分别由式(5)和式(6)给出。

(5)

(6)

式中，ρa为空气密度。

2基于FLUENT的UDS框架

2.1UDS输运方程

UDS(User Defined Scalar)，即用户自定义标量。FLUENT商用软件提供了求解UDS输运方程的相关功能，只需待求解方程满足输运方程一般形式，即可通过宏定义进行迭代求解。FLUENT中UDS输运方程的一般形式如式(7)。

(7)

式中，φk为用户自定义的一个标量，Fi，Гk，Sφk分别为输运方程的对流项、扩散项和源项系数。以三维问题为例，根据水滴相控制方程，需要定义4个UDS，分别是水滴容积分数α，水滴速度分量u，v，w。对照式(1)、式(2)和式(7)，即可得到4个UDS输运方程所对应的系数如表1。

表1　输运方程系数表

根据以上定义，利用FLUENT宏函数进行编程和调用，对流项、扩散项和源项的定义宏分别为DEFINE_UDS_FLUX、DEFINE_DIFFUSIVITY、DEFINE_SOURCE。以预先计算收敛的空气相流场作为已知条件，耦合水滴相边界条件，进行迭代求解，即可获得各网格点的水滴容积分数及水滴速度矢量。

2.2边界条件处理

飞机和发动机结冰模拟中，一般为绕流流动，因此边界条件有远场边界条件及壁面边界条件。对于远场边界，水滴容积分数α∞可由液态水含量计算得到，如式(8)。

(8)

远场处认为水滴根随空气一同运动，无相对速度，故水滴速度即为空气速度。

u∞=ua，∞

(9)

对于壁面边界，当过冷水滴撞击到壁面处时，将不会继续以水滴形态留存于流体计算域中，而是会发生结冰，或者变成水膜、水珠等方式附着在壁面处。由假设忽略水滴撞击壁面处的反弹和飞溅，则固壁面对于水滴而言，相当于一个单向出口，只许流出计算域，而不能流入计算域。因此必需将撞击到壁面处的水滴予以排除，否则会产生非物理解。

水滴在物体表面的绕流，根据水滴速度矢量方向与表面法线向量的关系，可将壁面分为撞击区和非撞击区，如图1所示。对于撞击区的面网格，其UDS输运方程的对流项按表1所示进行计算，即允许水滴流出计算域。而对于非撞击区的面网格，对流项设定为返回0值，即表示此处没有水滴流入计算域。

图1　物体表面水滴速度示意图

2.3数值扩散项

计算中发现，在某些局部区域，一般在绕流物体后方上下两股气流交汇处，会出现水滴容积分数异常增大的情况，体现在某个网格节点的水滴容积分数值可能比它周围网格值大100倍以上。Tony把这种现象解释为两个沿不同方向运动的水滴颗粒交汇后，变成以同一个速度沿同一方向运动了，如图2所示[8]。这使得交汇点的密度出现了一个无限大的阶跃，这个无限大的密度脉冲，就可能带来解的奇异。

图2　水滴交汇后的情况

为了及时清除水滴容积分数异常的区域，在水滴连续性方程(1)中加入人工数值扩散项。加入数值扩散项后的连续性方程变为式(10)。

(10)

数值扩散系数Г的选取具有一定的自适应性，在正常区域，Г应趋于0，以保证解的精度；在水滴容积分数异常区域，Г应足够大，使之可以迅速清除奇异区。基于以上要求，本文构造数值扩散系数如式(11)、式(12)所示。

Γ=a(eb-1.0)

(11)

(12)

式中，下标P代表当前网格，Ni代表与之相邻网格，b表征周围网格和当前网格水滴容积分数的偏离程度，a是比例因子，可根据经验进行调整。由式(11)可知，当前网格容积分数正常时，b值趋于0，扩散系数也趋于0。而对于异常区域，偏离程度越大，b值越大，指数型函数也可以确保扩散系数足够大，从而达到快速清除奇异区的目的。综上所述，2.1节中表1所列UDS(α)方程的扩散项系数Гk应变为ρГ。

3算例验证

为了验证建立模型的准确性，与相关试验结果进行对比。由于本文提出的方法同时适用于二维及三维物体，故分别选取二维圆柱及三维球面作为对比算例。对比参数为局部水滴收集系数β，其计算式为式(13)。

(13)

式中，u为壁面处水滴速度矢量，n为壁面法向向量，V∞为远场速度，αn为水滴的相对容积分数，定义为壁面处水滴容积分数与远场水滴容积分数之比。由模型假设可知，它基于单一水滴直径，实际上，云层中的水滴不可能是同一种直径大小。为了提高模拟精度，使用呈某种分布的若干种直径来代替实际情况，最常见的是Langmuir-D分布[9]，由7种水滴直径按一定比例构成。计算中首先针对每种水滴直径(dp)i求解一次，得到相应水滴直径的局部水滴收集系数βi，最后由式(14)计算得到综合收集系数。

(14)

式中pi为第i种直径水滴所占液态水含量比重。

3.1二维圆柱

二维圆柱作为典型的二维情况，在多篇文献中都作为验证模型使用[10-11]。本文计算选取与试验相同的条件：圆柱直径10.16 cm，自由流速80 m/s，入口空气密度1.097 kg/m3，环境压力89 867 Pa。计算7种水滴直径下的收集系数，平均容积直径(MVD)为16 μm，满足Langmuir-D分布。图3显示了7种水滴直径下的局部水滴收集系数以及由式(14)计算得到的Langmuir-D分布下的综合收集系数，可以看出水滴直径越大，局部水滴收集系数峰值和撞击极限均增加。这是由于水滴直径越大，质量越大，惯性力越大，受气流影响越小，因此更容易撞击到物体表面。

图3　二维圆柱不同直径水滴局部收集系数

图4　二维圆柱水滴收集系数预测值与试验结果对比

图4显示了二维圆柱局部水滴收集系数的计算结果与试验值的比较。图中红色圆圈表示使用单一平均容积直径MVD的计算结果，蓝色三角表示采用相应Langmuir-D分布的7种水滴直径的综合收集系数，阴影部分为试验重复性范围。由图4可知，2种计算结果与试验值均吻合较好，单一MVD计算得到的局部水滴收集系数峰值较Langmuir-D分布的综合收集系数略大，但基本都处于试验范围之内。而对于撞击极限预测，使用Langmuir-D分布的预测结果更接近试验值，使用单一MVD水滴直径的撞击极限预测值偏小。

3.2三维球面

第二个验证算例选自Bidwell等人[12]对三维球面所做的水滴撞击试验。试验对象为直径为15.04 cm的球体，自由流速度75 m/s，环境静温度7 ℃，环境压力95 840 Pa，水滴平均容积直径18.6 μm，满足Langmuir-D分布。图5显示了不同直径下的局部水滴收集系数，以及据此计算出的Langmuir-D分布下综合收集系数。图6给出了预测值与试验数据的对比，可以看出本文的预测结果与试验数据非常吻合，验证了本文方法在三维情况下的有效性。

图5　三维球面不同直径水滴局部收集系数

图6　三维球面水滴收集系数预测值与试验结果对比

4结论

本文提出了一种基于UDS框架的水滴撞击特性计算方法，通过对水滴撞击边界进行处理从而将撞击水从计算域中排除，引入数值扩散项避免了解的奇异性。通过二维及三维算例的试验对比，验证了本文方法的准确性和有效性，能够满足工程应用需要。

参考文献(References)：

[1]VENKATARAMANI K S,PLYBON R C,HOLM R G,et al.Aircraft engine icing model[R].AIAA 2008-440.

[2]杨倩,常士楠,袁修干.水滴撞击特性的数值计算方法研究[J].航空学报,2002,23(2):173-176.

[3]杨倩,常士楠,袁修干.发动机进气道水滴撞击特性分析[J].北京航空航天大学学报,2002,28(3):362-365.

[4]韩凤华,张朝民,王跃欣.飞机机翼表面水滴撞击特性计算[J].北京航空航天大学学报,1995，21(3)：16-21.

[5]张强,胡利,曹义华.过冷水滴撞击三维机翼的数值模拟[J].航空动力学报,2009,24(6):1345-1350.

[6]杨胜华,林贵平,申晓斌.三维复杂表面水滴撞击特性计算[J].航空动力学报,2010,25(2):284-290.

[7]申晓斌,林贵平,杨胜华.三维发动机进气道水滴撞击特性分析[J].北京航空航天大学学报,2011,37(1):1-5.

[8]吴孟龙,常士楠,冷梦尧,等.基于欧拉法模拟旋转帽罩水滴撞击特性[J].北京航空航天大学学报,2014,40(9):1263-1267.

[9]胡剑平,刘振侠,张丽芬.发动机整流支板大尺寸过冷水滴撞击特性[J].航空学报,2011,32(10):1778-1785.

[10]CROWE C T.Review-Numerical models for dilute gas-particle[J].Journal of Fluids Engineering,1982(104):297-303.

[11]TONG X L,LUKE E A.Eulerian simulations of icing collection efficiency using a singularity diffusion model[R].AIAA 2005-1246.

[12]PAPAKAKIS M,WONG S,RACHMAN A A,et al.Large and small droplet impingement data on airfoils and two simulated ice shapes[R].NASA/TM-2007-213959.

[13]BOURGAULT Y,BOUTANIOS Z,HABASHI W G.Three-dimensional eulerian approach to droplet impingement simulation using fensap-ice,part1:model,algorithm,and validation[J].Journal of Aircraft,2000,37(1):95-103.

[14]WIROGO S,SRIRAMBHATLA S.An eulerian method to calculate the collection efficiency on two and three dimensional bodied[R].AIAA 2003-1073.

[15]BIDWELL C S,MOHLER S R.Collection efficiency and ice accretion calculations for a sphere,a swept MS(1)-317 wing,a swept NACA-0012 wing tip,and axisymmetric inlet,and a Boeing 737-300[R].AIAA-95-0755.

(责任编辑：宋丽萍英文审校：隋华)

Calculation of water droplet impingement characteristics based on UDS frame

MAO Xiao-dong1，BU Xue-qin2，Zhao Guo-chang1，Wang Jian-ming1

(1.Faculty of Aerospace Engineering,Shenyang Aerospace University,Shenyang 110136，China;2.School of Aeronautic Science and Engineering,Beijing University of Aeronautics and Astronautics,Beijing 100191,China)

Abstract：In order to predict the water droplet impingement characteristics of 2-dimensional and 3-dimensional surfaces under compressible and uncompressible flow conditions,a calculation method based on UDS frame was approached.The air-droplet two-phase flow was considered as one-way coupled.The airflow field could be solved previously and independently using N-S equations.The droplet control equations were user-programed and iterated under fluent UDS frame.To remove the droplet from wall boundaries,a treatment of boundary conditions on the impingement surfaces was given.A self-adaptation artificial numerical diffusion model was proposed to enhance the iteration stability and avoid the singularity of droplet volume fraction in localized region.Validation cases were performed for typical 2-dimensional cylinder and 3-dimensional sphere,through the comparison with the published experimental results,the approached calculation method was proved to be accurate and effective in predicting the droplet impingement characteristics on complex surfaces under various flow conditions.

Key words：UDS frame;water droplet impingement;icing;two-phase flow;numerical calculation

doi:10.3969/j.issn.2095-1248.2016.01.002

中图分类号：V211.3

文献标志码：A

文章编号：2095-1248(2016)01-0008-05

作者简介：毛晓东(1984-)，男，辽宁大连人，副教授，主要研究方向：飞机与发动机防除冰，E-mail：mxdbh@163.com。

基金项目：国家自然科学基金(项目编号：51206008)

收稿日期：2015-09-17