Numerical simulation of droplet deformation in uniform airflow
-
摘要:
飞机在结冰环境下飞行将会遇到过冷水滴,水滴与气流的相互作用会使水滴变形,运动轨迹发生改变。针对水滴在均匀气流中运动变形过程,采用相场方法捕捉了水滴形态变化,结合格子玻尔兹曼通量求解器计算水滴‒空气两相流流场,比较了不同We数和Oh数、不同环境温度下的水滴形态变化及水滴运动加速过程,分析了水滴后缘流动发展特征及水滴变形过程的动力学行为规律;提出了水滴在气流中的运动阻力模型,进行了水滴撞击特性计算,并与实验结果进行对比以验证水滴阻力模型的准确性。研究发现,水滴在气流中变形具有轴向压缩与径向压缩交替进行的特征,轴向压缩过程中的水滴最大纵横比随时间逐渐减小;水滴临近破碎时,边缘较厚,中心逐渐变薄,直至破碎;低温环境影响水滴变形幅度,但对于水滴运动过程的影响较小。本文建立的水滴阻力模型可实现更准确的水滴撞击特性计算,为后续飞机结冰计算提供准确输入。
Abstract:Aircraft flying in icing conditions will encounter supercooled water droplets. These droplets can deform due to the interaction with the airflow, leading to the change of the droplet trajectory. The phase field method is applied to capture the droplet morphology during the droplet deformation and movement in a uniform airflow. Combined with the lattice Boltzmann flux solver, the flow field of the droplet-air two-phase flow is simulated. Time evolutions of the droplet aspect ratio and the accelerations under different We numbers, Oh numbers, and ambient temperatures are compared. Flow development in the wake region of the droplet is analyzed, and the dynamic behavior of the water droplet in the deformation process is revealed. A drag model of the droplet in the airflow is proposed, and the droplet impact characteristics are simulated and compared with experimental data to validate the accuracy of the drag model. It is found that, the droplet has the characteristics of alternating axial compression and radial compression, and the maximum aspect ratio in each axial compression decreases with time. When the water droplet is close to breakup, the edge is thicker while the center gradually becomes thinner until the breakup. Low ambient temperature affects the droplet deformation amplitude but has little effect on the movement of the droplet. Moreover, the proposed droplet drag model can achieve higher accuracy for the droplet impact characteristics calculation, which can provide a more accurate input condition for subsequent aircraft icing prediction.
-
Keywords:
- aircraft icing /
- uniform airflow /
- supercooled droplet /
- numerical simulation
-
0. 引 言
飞机在结冰气象条件下飞行时,空气中的过冷水滴会撞击在机翼表面,造成机翼结冰,严重威胁飞机飞行安全[1-2]。水滴在高空气流作用下,可能会发生变形破碎等动力学行为,这些行为会影响其运动轨迹,撞击在机翼表面不同位置。机翼表面的水滴收集量是飞机结冰量的重要影响因素,因此研究高空气流环境中的水滴变形运动有助于揭示水滴运动力学行为规律,提高飞机结冰问题中水滴撞击特性的预测精度。
在水滴变形破碎过程中,Weber数(We)与Ohnesorge数(Oh)是两个无量纲数,控制着水滴变形破碎机制:在低Oh数(Oh<0.1)下,水滴的变形破碎行为主要受We数控制,随着We数的增加,水滴会发生变形、振动破碎、袋式破碎、雄蕊式破碎等行为。随着Oh数的增加,水滴的变形程度略有降低[3]。
在实验研究方面,一些研究工作探索并提出了不同We数与Oh数下水滴变形过程的转变机制。Kulkarni等[3]重点考虑了单颗低黏度液滴的破碎情况,从理论上分析了从振动破碎到袋式破碎的转变及We数范围,并与实验结果进行了比较,得到了水滴在径向和气流方向上形态变化特征的理论模型。García-Magariño等[4]通过旋转翼型撞击水滴研究了水滴变形与破碎转变过程,在研究中特别考虑了气流加速度对水滴行为的影响,其研究结果指出:基于瞬时速度的We数仅可以定性描述水滴动力学行为机制,这对水滴变形的非定常效应具有重要影响。Soni等[5]实验研究了液滴与倾斜连续气流相互作用时的变形和破碎,发现水滴临界破碎We数随着液滴直径的增加而减小。
水滴-空气大密度比界面的准确捕捉是探索水滴动力学行为的重要内容,在数值研究方面,数值模拟可以提供更加详细的水滴相关行为信息。水气多相流数值模拟方法主要包括Volume of Fluid(VOF)方法[6-8]、Level Set(LS)方法 [9-10]、Front Tracking(FT)方法[11]等。Han等[12]应用FT方法研究了小密度比下二维轴对称液滴,发现不同的临界We数与水滴Reynolds(Re)数呈现负相关关系。Jain等[13]通过实验与VOF方法研究了大We数下的水滴破碎问题,发现了一种多袋式破碎机制(We = 80),并指出水滴破碎机制的转变是由于水滴边缘的动力学行为。为了准确捕捉水气界面,Strotos等[14]应用Fluent采用不同的压力插值格式预测了中等We数下的水滴变形与破碎行为,获得了适用于水滴形态捕捉的最佳格式。Yang等[15]采用带压缩项的VOF方法结合LS方法与理论分析研究了水滴在连续气流中不同破碎行为之间的转变机制(0.001 < Oh < 2),提出了基于Rayleigh–Taylor不稳定性的临界We数预测模型。Yang等[15-16]还采用CLSVOF方法研究了水滴的袋式破碎与雄蕊式破碎,结果表明破碎模式决定于水滴迎风直径和最不稳定Rayleigh-Taylor波数的比值。对于水滴破碎后的水滴尺寸分布情况,Liang等[17]采用VOF方法模拟了不同破碎机制下的水滴尺寸分布,研究所得的对数正态分布近似结果与实验结果[13]大致吻合。对于超声速情形,We数更大,Xiao等[18]采用CLSVOF方法追踪水滴界面演化过程,且气相采用可压缩算法进行求解,研究发现超声速下的水滴比亚声速下的水滴更早开始破碎。Liu等[19]研究指出高速气流与水滴之间的剪切作用是水滴表面不稳定性的原因,而背风面的不稳定性则是源于涡结构。除了We数与Oh数外,水滴与空气的密度比和黏度比也是水滴动力学行为的重要影响因素。Yang等[20]采用CLSVOF方法研究了水滴雾化过程,指出更小的水气密度比会导致更大的水滴变形幅度,在Oh数较小、密度比超过32时,密度比对水滴变形与破碎过程有显著影响。Marcotte与Zaleski[21]分析了密度比为10~2000对于水滴变形破碎机制的影响,研究发现临界We数随密度比的降低而增加。Jain等[22]研究了中等We数下密度比和Re数对水滴变形破碎及运动的影响,其结果表明Re数、RT不稳定波数和黏度比可以更好地代表水滴破碎模式。Xiao等[23]分析了液体黏度对液滴变形破碎的影响,在大Oh数下临界We数与Oh数呈现正相关关系,他们还预测了不同We数与Oh数下的水滴破碎开始时间。
宏观飞机结冰计算中,1 m3的气流中含有数以千万计的过冷水滴,采用水滴阻力模型计算水滴运动轨迹将极大减少计算量,不同的水滴形态具有不同的水滴阻力,其水滴阻力系数有所不同。在微观数值计算层面,Chu等[24]采用耦合LS与VOF方法(CLSVOF)模拟了不同Oh数下的液滴动力学行为,研究发现当气动力占主导时,液滴黏度增加会推迟水滴破碎时间,且他们还基于水滴迎风直径修正了液滴阻力系数公式。Shao等[25]采用守恒LS方法研究了We数对于水滴变形过程中非稳态阻力系数的影响,研究发现低We数下,非稳态阻力系数与稳态阻力系数的差值与非稳态参数呈线性关系。在水滴阻力模型方面,当前飞机结冰广泛采用的是Clift模型[1],这一模型的特点是在圆球阻力系数和圆盘阻力系数间作插值,插值系数决定于We数。相较于水滴破碎,水滴变形与振荡的物理机制相对简单[2],TAB模型(Taylor analogy breakup model)将水滴振荡与变形类比为一个弹簧–质量系统,引入了表面张力、水滴黏性与气动力,通过水滴变形量与圆球阻力系数修正变形水滴的阻力系数。
上述研究工作主要针对水滴破碎行为开展了研究,探索了不同水滴破碎形式,对于低We数下水滴发生变形和振荡破碎两种模式的研究较少,而这一过程的物理机理研究是飞机结冰所需考虑的重要问题之一。本文的主要目的在于针对飞机可能遇到的结冰工况,研究气流引起水滴变形与破碎的过程,建立水滴运动力学模型应用于飞机结冰水滴轨迹预测。本文将对不同We数与Oh数、不同环境温度下的水滴在连续均匀气流中的变形特征与运动加速规律进行研究,探索水滴与空气相互作用下的水滴动力学行为机制,并基于计算结果建立水滴运动阻力模型。本文研究工作可以为飞机结冰问题中的水滴轨迹计算建模提供研究基础,实现机翼表面水收集系数的准确计算。
1. 数值计算方法
为了模拟水滴在连续均匀气流中的变形与运动过程,本文在轴对称坐标系下,应用相场方法捕捉水滴与空气间的界面,采用多相流格子玻尔兹曼通量求解器(multiphase lattice Boltzmann flux solver, MLBFS)[26]计算水滴与空气流场,然后基于所采用的计算方法进行数值模拟,并与Li等[27]的实验结果进行对比,验证数值算法的准确性。
1.1 控制方程
相场方法可用于捕捉大密度比多介质流动的复杂界面变化,相场控制方程为[28]:
$$ \begin{split} &\frac{{\partial C}}{{\partial t}} + \frac{{\partial C{u_r}}}{{\partial r}} + \frac{{\partial C{u_z}}}{{\partial z}} = \\ &\qquad - \frac{{C{u_r}}}{r} + M\left( {\frac{{{\partial ^2}{\mu _C}}}{{\partial {r^2}}} + \frac{{{\partial ^2}{\mu _C}}}{{\partial {z^2}}} + \frac{1}{r}\frac{{\partial {\mu _C}}}{{\partial r}}} \right) \end{split} $$ (1) 其中:C是相变量,表示当前网格单元中水的体积分数,在水滴内部取值为1,在空气中取值为0,水气界面定义为C = 0.5等值线;t为时间;r与z分别为径向与轴向;ur与uz为水滴速度分量;M是迁移率;μC是化学势,其表达式为[28]:
$$ \begin{split} &{\mu _C} = 4{\beta _C}C\left( {C - 1} \right)\left( {C - 0.5} \right) - \\ &\qquad k\left( {\frac{{{\partial ^2}C}}{{\partial {r^2}}} + \frac{{{\partial ^2}C}}{{\partial {z^2}}} + \frac{1}{r}\frac{{\partial C}}{{\partial r}}} \right) \end{split} $$ (2) βC与k是与水气间表面张力σ及数值模拟中水气界面过渡带宽度W有关的参数,可通过如下关系式确定[28]:
$$ \begin{split} {\beta _C} =& \frac{{12\sigma }}{W} \\ k =& 1.5\sigma W \end{split} $$ (3) 相场方法中表面张力为[28]:
$$ {\boldsymbol{F_s}} = - C\nabla {\mu _C} $$ (4) Fs作为外力源项添加至动量方程。
在求解相场控制方程时,为了保证水气界面的准确捕捉,提高计算精度,对流项采用了有限差分型五阶WENO格式进行离散[29]。此外,计算过程中的质量守恒是数值解的一个重要性质,本文在后续数值模拟计算中采用质量守恒算法[28],以保证连续气流场中变形水滴总体质量守恒。
除了相场控制方程外,水滴与空气流场计算采用了轴对称坐标系下的MLBFS方法[26],该算法将有限体积法与格子玻尔兹曼方法(lattice Boltzmann method, LBM)相结合用于计算网格单元的界面通量。基于MLBFS方法的宏观控制方程为[28,30]:
$$ \begin{split} &\frac{{\partial p}}{{\partial t}} + \frac{{\partial \rho {u_r}c_s^2}}{{\partial r}} + \frac{{\partial \rho {u_z}c_s^2}}{{\partial z}} = \\ &\qquad \left( {{u_r}\frac{{\partial \rho }}{{\partial r}} + {u_z}\frac{{\partial \rho }}{{\partial z}}} \right)c_s^2 - \frac{{\rho {u_r}c_s^2}}{r} \end{split} $$ (5) $$ \begin{split} &\frac{{\partial \rho {u_r}}}{{\partial t}} + \frac{{\partial \rho u_r^2}}{{\partial r}} + \frac{{\partial \rho {u_r}{u_z}}}{{\partial z}} = \frac{{\partial p}}{{\partial r}} + \\ &\qquad \frac{\partial }{{\partial z}}\left[ {\mu \left( {\frac{{\partial {u_r}}}{{\partial z}} + \frac{{\partial {u_z}}}{{\partial r}}} \right)} \right] + \frac{\partial }{{\partial r}}\left( {2\mu \frac{{\partial {u_r}}}{{\partial r}}} \right) + \\ &\qquad \frac{{2\mu }}{r}\left( {\frac{{\partial {u_r}}}{{\partial r}} - \frac{{{u_r}}}{r}} \right) - \frac{{\rho u_r^2}}{r} + {F_{sr}} \end{split} $$ (6) $$ \begin{split} &\frac{{\partial \rho {u_z}}}{{\partial t}} + \frac{{\partial \rho u_z^2}}{{\partial z}} + \frac{{\partial \rho {u_r}{u_z}}}{{\partial r}} = \frac{{\partial p}}{{\partial z}} + \\ &\qquad \frac{\partial }{{\partial r}}\left[ {\mu \left( {\frac{{\partial {u_r}}}{{\partial z}} + \frac{{\partial {u_z}}}{{\partial r}}} \right)} \right] + \frac{\partial }{{\partial z}}\left( {2\mu \frac{{\partial {u_z}}}{{\partial z}}} \right) + \\ &\qquad \frac{\mu }{r}\left( {\frac{{\partial {u_r}}}{{\partial z}} + \frac{{\partial {u_z}}}{{\partial r}}} \right) - \frac{{\rho {u_r}{u_z}}}{r} + {F_{sz}} \end{split} $$ (7) 其中:p为压力,ρ为密度,cs是声速,μ是动力黏度,Fsr与Fsz是表面张力分量。
本文求解相场方程与流动控制方程的非稳态项采用三阶TVD Runge-Kutta时间推进格式[31]。为了保证空气流经水滴后的边界条件不会影响水滴动力学行为,需要使水滴始终处于计算域中,为此本文采用了移动计算域的方法[32],该方法能够以较小的计算量模拟水滴更长时间的变形与运动情况。
1.2 计算设置
水滴在连续气流中运动,空气流经水滴背风面后为水滴后缘区域,考虑到此区域内可能出现复杂的流动结构,故此区域应布置足够的网格。本文采用均匀笛卡尔网格进行计算,计算域设置如图1所示。
在5D×15D的轴对称计算域内(D为水滴直径),初始静止的水滴置于距离速度入口边界2.5D位置处,水滴径向边界与另一轴向边界为自由出口边界。在计算中将Navier-Stokes(N-S)方程无量纲化后进行求解,空气与水滴的物理性质为(T1 = 25℃):水滴密度ρl = 997.08 kg/m3,水滴动力黏度μl = 8.937×10−4 Pa∙s,空气密度ρg = 1.185 kg/m3,空气动力黏度μg =
1.8365 ×10−5 Pa∙s,表面张力σ =0.07197 N/m。1.3 算法验证
应用上述介绍的数值模拟算法,定义计算采用的均匀网格单元尺度为h,分别采用粗网格(D/h = 50)、中等网格(D/h = 100)和细网格(D/h = 150)对直径D = 2.2 mm的水滴在气流速度Ug = 13.2 m/s的空气中的运动变形情况进行数值计算。提取相变量C = 0.5等值线作为水滴边界,不同网格下水滴形态的计算结果如图2所示,其中无量纲时间t* = Ugt/D。结果表明,中网格和细网格的计算结果差异很小,本文后续采用中等网格进行数值计算与分析。
将水滴运动过程中形态变化的数值计算结果与Li等[27]的实验结果进行对比,如图3所示。可以看出,水滴在来流方向上逐渐压缩至椭球状(t* = 24),随后逐渐恢复至圆球状。定义变形的水滴沿来流方向的直径为Dz(即水滴轴向直径),垂直于来流方向的直径为Dr(即水滴径向直径),Dr的计算结果与实验结果对比如表1所示。计算中的空气为均匀来流,水滴形态是对称的,而实验中的水滴是滴落进空气流场中的,引入了非对称初始扰动,因此本文计算结果与文献实验结果相比存在一定误差,但平均误差小于7.5%,总体来说,本文的计算方法可以较为准确地捕捉水滴在气流中的变形过程。
时刻 Dr计算值/mm Dr实验值/mm t* = 1.2 2.33 2.41 t* = 6 2.64 2.44 t* = 24 3.40 2.85 t* = 30 3.07 2.81 t* = 36 2.66 2.66 t* = 42 2.38 2.50 2. 数值计算结果与分析
水滴在连续气流中的运动变形过程受水滴尺寸、空气速度、环境温度等因素的影响,本文在不同工况下进行数值计算,以分析水滴与空气相互作用导致水滴变形的物理机制以及不同影响因素对水滴运动过程的影响规律。飞机结冰情况下,环境温度通常在−40℃以上,液态水含量通常为1.0 g/m3左右,水滴平均直径为15~40 μm,且近年来研究表明50 μm以上的大水滴也会造成飞机结冰。飞机结冰事故统计结果表明,飞机结冰时,空气速度通常在100 m/s以下,飞机在云层中遭遇过冷水滴发生结冰的时长不超过45 min[1]。因此,本文基于以上数据信息选取了数值模拟计算工况,具体如表2所示,定义We = ρgUg2D/σ,Oh = μl/(ρlDσ)1/2,T1 = 25℃,T2 = 0℃,T3 = −10℃,T4 = −20℃。
编号 水滴直径
D/μm空气速度
Ug/(m·s−1)温度
T/℃We数 Oh数 1 50 30 25 0.74 0.015 2 50 40 25 1.32 0.015 3 50 50 25 2.06 0.015 4 50 60 25 2.96 0.015 5 70 30 25 1.04 0.013 6 70 40 25 1.84 0.013 7 70 50 25 2.88 0.013 8 70 60 25 4.15 0.013 9 100 30 25 1.48 0.011 10 100 40 25 2.63 0.011 11 100 50 25 4.12 0.011 12 100 60 25 5.93 0.011 13 160 30 25 2.37 0.008 14 160 40 25 4.22 0.008 15 160 50 25 6.60 0.008 16 160 60 25 9.50 0.008 17 30 80 0 3.29 0.038 18 30 80 −10 3.34 0.055 19 30 80 −20 3.42 0.09 20 40 80 0 4.38 0.033 21 40 80 −10 4.46 0.047 22 40 80 −20 4.55 0.078 23 80 60 0 4.93 0.023 24 80 60 −10 5.01 0.034 25 80 60 −20 5.12 0.055 26 100 50 0 4.28 0.021 27 100 50 −10 4.35 0.03 28 100 50 −20 4.45 0.049 29 140 40 0 3.84 0.017 30 140 40 −10 3.9 0.025 31 140 40 −20 3.99 0.041 2.1 We数与Oh数的影响
水滴尺寸D和气流速度Ug是影响水滴在连续气流中变形及运动的重要因素,其影响可以通过We数与Oh数进行分析。图4给出了表2中算例4、算例8和算例12的水滴形态随时间的变化过程。从图4可以看出,不同直径的水滴在气流中均表现出轴向压缩与径向压缩交替进行的特征:气流中运动的水滴受气动力和表面张力的综合作用,初始时刻因气动力的影响,表面张力不足以维持水滴球形形态,水滴发生轴向压缩,水滴变形量增加;待水滴径向直径达到最大时,在表面张力的主导作用下,水滴向球形形态恢复。水滴变形过程是气动力和表面张力综合作用下水滴动力学过程的动态表现。图中红色分隔表示水滴完成了一次轴向压缩与径向压缩(下文同),这与Li等[27]和Rimbert等[33]所得出的水滴周期性震荡的研究结论一致。从图4中水滴形态变化可知,随着水滴尺寸的增加,We数增加,Oh数变化较小,水滴的轴向压缩程度增加。We = 2.96时,水滴变形特征不明显,但当We数增加到5.93时,水滴首先在轴向压缩阶段向后缘方向弯曲(t* = 35),然后在表面张力与气流持续作用下反向弯曲(t* = 45),接着恢复成圆盘状(t* = 60),最后因表面张力的作用逐渐恢复成球形(t* = 85)。此后,水滴将进入新一轮的压缩过程。值得注意的是,由于水滴向后运动,水滴与空气间的速度差减小,水滴在第二次压缩过程中表现为圆盘形,并未表现出第一次压缩过程中的向后缘方向弯曲的形态(t* = 135),水滴随后表现出反向弯曲形态(t* = 145)。
为量化分析水滴变形过程,这里引入水滴纵横比Ar用于表示水滴变形程度,Ar定义为水滴径向直径Dr与轴向直径Dz之比。图5给出了不同We数下水滴纵横比Ar的变化过程,可以得知:水滴在变形过程中受到气动力和表面张力的影响,轴向压缩与径向压缩交替进行,首次轴向压缩达到的最大纵横比最大;随着We数增大,达到首次轴向压缩的最大纵横比所需的时间在逐渐增加;水滴每一次压缩过程中的最大纵横比在逐渐降低,这是因为水滴与空气之间的速度差在逐渐减小,水滴受力也随之减小。此外,水滴纵横比曲线还表明:当We = 5.93时,水滴变形更加剧烈,水滴在首次恢复过程中表现出微小的二次轴向压缩特征(t* = 49到t* = 61),这是水滴由向前缘弯曲过渡至圆盘状过程中,轴向直径减少所导致。
在气流作用下的形态变化过程中,水滴后缘空气流场结构也会发生相应变化,图6给出了We = 5.93、Oh = 0.011的水滴后缘空气流场结构变化及来流方向无量纲空气速度Ua*的变化。空气速度Ua通过空气来流速度Ug进行无量纲化,即Ua* = Ua/Ug。从图6可知,空气流场速度在来流方向上可以增加至初始来流速度的1.5倍左右。在变形初期,水滴被压缩为圆盘状,后缘流场发展出涡结构,水滴纵横比大,涡距离水滴后缘较近(t* = 35);在水滴第一次恢复成球形的过程中,水滴后缘的涡向后移动并继续发展,涡结构沿轴向拉长(t* = 85);在之后的轴向压缩过程中,水滴变形程度增加、径向直径增加,水滴后缘涡结构持续发展,沿径向增加、沿轴向减小(t* = 145)。水滴后缘涡结构的发生发展是水滴与持续来流相互作用的结果。
水滴在空气流场中受到气动力作用,运动速度逐渐增加。图7给出了算例4、8、12的无量纲水滴速度U*的变化。在第一次轴向压缩过程中,We = 5.93时的无量纲水滴速度相较于We = 2.96、4.15时,表现出增速先增大随后降低的特性,这是由于初始时水滴变形剧烈,水滴迎风面积增加,所受气动力增大,但随着变形程度降低,无量纲水滴速度增速逐渐下降。较低We数下的无量纲水滴速度变化如图8所示,由图可知,在运动过程中,水滴增速逐渐减缓,且We数更低时无量纲水滴速度更大,这是由于低We数下的水滴变形程度均较小,质量更小的水滴随气流运动更易被增加速度。
当We数继续增加,水滴的变形程度更大。如图9所示的算例14、15、16的水滴变形情况:当We数为4.22时,水滴形态呈现椭球形-球形交替转换状态,而当We数增至9.5时,水滴的表面张力不能使其恢复至球形,水滴变形程度增加;水滴形态变化表现为先向后缘弯曲(t* = 25),然后呈圆盘状,接着在气流作用下,水滴表现出反向弯曲,其径向直径持续增加,边缘增厚,中心逐渐变薄,直至气流吹破水滴。
算例14、15、16相应的水滴纵横比Ar随时间变化过程如图10所示。We数较低时,水滴呈现轴向压缩和径向压缩交替进行的特性。当We数增加,水滴最大纵横比变大,达到最大纵横比所需时间变长,且由于表面张力的恢复作用,最大纵横比随时间逐渐下降。当We数增至9.5时,水滴变形过程表现为轴向压缩至水滴中心破碎,纵横比随时间先增加后减小。其中,纵横比增加是因为轴向压缩使得水滴中心厚度变薄、径向直径Dr增加,但随后的纵横比减小是由于水滴反向弯曲,中间形成了空腔,在三维立体视角下,水滴的中心厚度“增加”了,因此水滴径向直径持续增加,但水滴纵横比反而降低了。
水滴后缘的流场结构变化与气流中的水滴形态有关,图11给出了We = 9.5(算例16)时,水滴后缘空气流场结构变化及来流方向无量纲空气速度变化。图11表明空气流场速度在来流方向上可以增加至初始来流速度的1.7倍左右。在变形初期,水滴前缘在强气流作用下向后移动,发生剧烈的轴向压缩,并在水滴径向边缘形成尖角(t* = 15),水滴后缘逐渐发展出涡结构;当t* = 35时,涡结构在径向和轴向上进一步发展,水滴中心位置开始变薄;当水滴呈现为反向弯曲时(t* = 50),水滴已经临近破碎,此时的水滴后缘涡较t* = 35时涡核向后移动。
Oh数相同,不同We数(算例14、15、16)下的水滴运动速度随时间变化情况如图12所示。We = 4.22时,水滴速度增加较为平缓,而随着We数的增加,由于气流作用下的水滴迎风面积增加,所受气动力更大,无量纲水滴速度迅速增加。但特别的是,Oh数相同、We数较小时,不同We数下的无量纲水滴速度发展曲线基本重合(如图13所示),可以归一化。这一结果表明,对于非剧烈变形的水滴,可以通过分析其速度发展曲线,建立水滴运动力学模型,从而实现飞机结冰环境下的水滴轨迹精确预测。
2.2 环境温度的影响
飞机结冰条件下,水滴周围环境温度均为低温,本小节针对低温下的水滴变形与运动过程进行了数值计算。低温下的水滴与空气物性参数如表3所示,不考虑水滴随气流运动过程中的结冰相变问题,因而环境温度的影响可归结为水滴与空气物性参数的影响。算例17~31采用了表3所示的水滴与空气物性参数[34]。
环境温度
T/℃水滴密度
ρl/kg·m−3空气密度
ρg/kg·m−3表面张力
σ/N·m水滴动力黏度
μl/10−3Pa·s空气动力黏度
μg/10−5Pa·s0 1000.00 1.293 0.07550 1.7921 1.72 −10 998.15 1.342 0.07710 2.6320 1.67 −20 993.60 1.395 0.07841 4.3300 1.62 图14给出了环境温度为−20℃、Oh = 0.041、We = 3.99(算例31)时的水滴形态变化,环境温度为0℃和−10℃下的水滴形态变化与其类似,均形成轴向压缩的椭球形水滴。不同温度下(算例29、30、31)水滴变形程度如图15所示,可见:环境温度越低,水滴表面张力越大,首次轴向压缩过程中所达到的最大纵横比越小;另外,每次轴向压缩与径向压缩过程中,水滴达到最大纵横比和最小纵横比的无量纲时间基本相同。
为了研究环境温度对水滴运动过程的影响,图16给出了算例29、30、31的无量纲水滴速度变化过程。在不同的低温环境下,无量纲水滴速度变化曲线基本重合,且由表3可知,不同环境温度下水滴的动力黏度变化最大,这表明水滴动力黏度对水滴运动过程影响较小。
2.3 水滴运动阻力模型
水滴在气流中的运动轨迹对于研究飞机结冰问题至关重要。由于水滴(~101 μm)和机翼(~101 m)的尺度差异巨大,采用同一套网格计算翼型的气动特性和微米级别的水滴形态变化及运动轨迹,其计算量极大,几乎是不可能的。可行的方法是建立水滴阻力模型,在宏观水平上,应用此模型计算机翼表面水滴收集系数。宏观飞机结冰数值计算首先需计算空气流场[1-2],再基于流场信息,求解气流中的水滴运动控制方程,确定水滴撞击在机翼表面的位置,并统计水滴收集量,获得水滴撞击特性。准确的水滴阻力模型可以为后续结冰相变计算提供准确的水滴收集系数作为输入条件,降低冻结相变过程的模拟误差,提升飞机结冰全过程的预测精度。
基于本文数值计算结果,水滴在气流中非剧烈变形时的无量纲水滴速度发展过程表明:水滴速度增加,水滴增速减缓;水滴速度从0逐渐增加,经过相当长时间的空气加速后,其速度会逐渐趋于气流速度,即无量纲水滴速度趋于1。采用如下公式通过最小二乘法可拟合无量纲水滴速度发展过程:
$$ {U^*} = 1 - {{\rm{e}}^{ - B{t^*}}} $$ (8) 以算例4和算例9为例,计算所得的无量纲水滴速度发展过程及拟合结果如图17所示,可见两者吻合良好。
基于公式(8),可得无量纲水滴加速度为:
$$ {a^*} = \frac{{d{U^*}}}{{d{t^*}}} = B{{\rm{e}}^{ - B{t^*}}} $$ (9) 水滴与气流间的速度差∆U = Ug–U,无量纲化的水滴与气流间的速度差∆U*为:
$$ \Delta {U^*} = \frac{{{U_g} - U}}{{{U_g}}} = 1 - {U^*} $$ (10) 结合式(8)和式(9),可以得出无量纲水滴加速度a*和水滴与气流间的速度差∆U*的关系为:
$$ {a^*} = B\Delta {U^*} $$ (11) 通过数据拟合可以获得各个工况条件下的拟合参数B,B与Oh数的关系如图18所示,图中数据表明拟合参数B大致随Oh数呈线性增长。
为了考虑不同温度下的水滴加速度,这里引入无量纲温度θ与拟合参数p1、p2、p3、p4,B与Oh、θ的关系为:
$$ \begin{split} &{B = \left( {{p_1}\theta + {p_2}} \right)Oh + {p_3}\theta + {p_4}} \\ &\qquad {Oh = \dfrac{{{\mu _l}}}{{\sqrt {{\rho _l}D\sigma } }},}\;\;{\theta = \dfrac{{T - {T_2}}}{{{T_1} - {T_2}}}} \end{split} $$ (12) 其中:p1 =
0.01762 ,p2 =0.02304 ,p3 =−8.07×10−5,p4 = 6.98×10−4。在飞机结冰数值计算中,水滴撞击特性的计算主要采用Euler法或Lagrange法[2],气流中的水滴运动遵循牛顿第二定律,在拉格朗日体系下,水滴轨迹的预测广泛采用Clift阻力模型[35],基于Clift模型的水滴运动控制方程为[36]:
$$ \begin{split} &\frac{{{\rm{d}}{\boldsymbol{U}}}}{{{\rm{d}}t}} = K\left( {{{\boldsymbol{U_g}}} - {\boldsymbol{U}}} \right) + \left( {1 - \frac{{{\rho _g}}}{{{\rho _l}}}} \right){\boldsymbol{g}} \\ &\qquad K = \frac{{18{\mu _g}}}{{{\rho _l}{D^2}}}\frac{{{C_d}Re}}{{24}} \\ &\qquad Re = \frac{{{\rho _g}\left\| {{\boldsymbol{U_g}} - {\boldsymbol{U}}} \right\|D}}{{{\mu _g}}} \end{split} $$ (13) 其中:U为水滴速度,Ug为空气速度,g为重力加速度,ρg和ρl分别为空气和水滴的密度,μg是空气动力黏度,Re为水滴相对雷诺数,Cd是水滴阻力系数,其表达式为[35-36]:
$$ \begin{split} & {C_d} = \varphi {C_{d,{\rm{disk}}}} + \left( {1 - \varphi } \right){C_{d,{\rm{sphere}}}} \\ &\qquad {C_{d,{\rm{disk}}}} = 1.1 + \frac{{64}}{{\text{π} Re}} \\ &\qquad {C_{d,{\rm{sphere}}}} = 0.36 + 5.48R{e^{ - 0.573}} + \frac{{24}}{{Re}} \end{split} $$ (14) Clift模型[35]考虑了水滴变形效应,认为水滴阻力系数Cd是圆盘形阻力系数Cd, disk和球形阻力系数Cd, sphere的组合,参数φ表示水滴接近圆盘形的程度,φ与相对We数有关,其表达式如下[36]:
$$ \begin{split} & \varphi = 1 - {\left( {1 + 0.007\sqrt {We} } \right)^{ - 6}} \\ &\qquad We = \frac{{{\rho _g}{{\left\| {{\boldsymbol{U_g}} - {\boldsymbol{U}}} \right\|}^2}D}}{\sigma } \end{split} $$ (15) Clift模型中的水滴阻力系数与We数直接相关,根据本文研究结果可知:在We数较低、水滴非剧烈变形情况下,水滴直径不同、空气速度相同时,不同直径的水滴速度变化具有明显差异性,与We数变化有关(如算例1、5和9);当水滴直径相同、空气速度不同时,水滴速度增长过程差异性较小,与We数变化无关(如算例2、3和4)。本文建立水滴非剧烈变形情况下的阻力模型,此模型不直接计算水滴形态的动态变化,而是将无量纲水滴加速度a*与水滴和气流间的速度差∆U*相关联,不与We数直接联系,如公式(11)所示。
气流中运动水滴所受的单位质量气动力为a,通过公式(9)与无量纲水滴速度U*和无量纲时间t*的定义,可得a与无量纲水滴加速度a*的关系为:
$$ {a^*} = \frac{{Da}}{{U_g^2}} $$ (16) 结合公式(10)和公式(11),可得水滴所受单位质量气动力为:
$$ {\boldsymbol{a}} = \frac{{B{U_g}}}{D}\left( {{{\boldsymbol{U_g}}} - {\boldsymbol{U}}} \right) $$ (17) 在拉格朗日体系下,基于此模型的水滴轨迹控制方程为:
$$ \frac{{{\rm{d}}{\boldsymbol{U}}}}{{{\rm{d}}t}} = \frac{{B{{U_g}}}}{D}\left( {{{\boldsymbol{U_g}}} - {\boldsymbol{U}}} \right) + \left( {1 - \frac{{{\rho _g}}}{{{\rho _l}}}} \right){\boldsymbol{g}} $$ (18) 其中B由公式(12)确定。
公式(13)和公式(18)通过四步Runge-Kutta法进行时间推进[2],确定新时刻的水滴位置及速度,最终获得水滴在机翼表面的撞击位置。
应用飞机结冰研究中水滴轨迹广泛采用的Clift模型[35]和本文提出的水滴运动阻力模型计算飞机结冰工况下机翼表面的水滴撞击特性,并与实验结果进行对比[37-38]。实验翼型为NACA23012,弦长c = 0.914 m,攻角α = 2.5 °,速度V = 78.2 m/s,水滴中值体积直径(median volume diameter, MVD)与液态水含量(liquid water content, LWC)如表4所示。
工况编号 a b c d MVD/μm 52 111 154 236 LWC/g·m−3 0.4 0.73 1.44 1.89 本文数值计算采用SST湍流模型、贴体网格求解空气流场及水滴轨迹控制方程,为保证网格无关性,分别采用了细网格(网格量
44000 )、中等网格(网格量22000 )和粗网格(网格量11000 )计算翼型表面的压力系数Cp,并与实验结果进行对比,对比结果如图19所示。计算结果表明,采用中等网格的计算结果与实验值已经吻合较好,因此本文基于中等网格进行计算与结果分析,中等网格翼型前缘如图20所示。
离散求解水滴轨迹控制方程,水滴收集系数β定义为机翼表面过冷水滴撞击量与可能的最大过冷水滴撞击量之比。计算结果与实验数据对比如图21所示。整体而言,在上翼面,采用Clift模型[35]和本文所提出模型所得的水滴收集系数均与实验结果吻合较好;在下翼面,基于Clift模型[35]的计算结果与实验值有一定误差,而基于本文模型所得的水滴收集系数与实验值更吻合,具有更高的预测精度。这表明本文所提出的水滴运动阻力模型在飞机结冰水滴撞击特性计算中具有更高的准确性,这将为后续结冰相变过程数值计算提供更准确的输入,有助于提升飞机结冰全过程的预测精度。
3. 结 论
本文针对水滴在连续气流中的变形及运动过程,基于相场方法捕捉水滴–空气界面,应用MLBFS通量求解器解算流场,分析数值计算结果,讨论了不同We数与Oh数、不同环境温度对水滴的变形及速度变化的影响,建立了水滴运动阻力模型,并进行了水滴撞击特性计算,完成了与实验结果的对比验证。所得结论如下:
1) 气流中的水滴变形表现出轴向压缩与径向压缩交替进行的特征,水滴在第一次轴向压缩时达到的最大纵横比最大,We数越大,则达到首次轴向压缩的最大纵横比所需的时间越长,轴向压缩的最大纵横比随时间逐渐减小。当水滴临近破碎时,水滴径向直径随时间增加,边缘较厚,水滴中心在气流作用下逐渐变薄,直至破碎。
2) 低温环境对水滴变形幅度有影响,对每次轴向压缩–径向压缩过程中水滴达到最大纵横比和最小纵横比的无量纲时间基本无影响,对气流中水滴加速过程的影响较小。
3) 相较于Clift模型,基于本文建立的水滴运动阻力模型计算的飞机结冰水滴撞击特性结果与实验值更为吻合,可获得更高精度的水滴收集系数,为飞机结冰全过程的准确预测提供了基础。
-
表 1 水滴径向直径Dr计算结果与实验结果[27]对比
Table 1 Comparison of the droplet radial diameter Dr between the numerical and experimental results[27]
时刻 Dr计算值/mm Dr实验值/mm t* = 1.2 2.33 2.41 t* = 6 2.64 2.44 t* = 24 3.40 2.85 t* = 30 3.07 2.81 t* = 36 2.66 2.66 t* = 42 2.38 2.50 表 2 数值模拟计算工况
Table 2 Computational conditions for numerical simulations
编号 水滴直径
D/μm空气速度
Ug/(m·s−1)温度
T/℃We数 Oh数 1 50 30 25 0.74 0.015 2 50 40 25 1.32 0.015 3 50 50 25 2.06 0.015 4 50 60 25 2.96 0.015 5 70 30 25 1.04 0.013 6 70 40 25 1.84 0.013 7 70 50 25 2.88 0.013 8 70 60 25 4.15 0.013 9 100 30 25 1.48 0.011 10 100 40 25 2.63 0.011 11 100 50 25 4.12 0.011 12 100 60 25 5.93 0.011 13 160 30 25 2.37 0.008 14 160 40 25 4.22 0.008 15 160 50 25 6.60 0.008 16 160 60 25 9.50 0.008 17 30 80 0 3.29 0.038 18 30 80 −10 3.34 0.055 19 30 80 −20 3.42 0.09 20 40 80 0 4.38 0.033 21 40 80 −10 4.46 0.047 22 40 80 −20 4.55 0.078 23 80 60 0 4.93 0.023 24 80 60 −10 5.01 0.034 25 80 60 −20 5.12 0.055 26 100 50 0 4.28 0.021 27 100 50 −10 4.35 0.03 28 100 50 −20 4.45 0.049 29 140 40 0 3.84 0.017 30 140 40 −10 3.9 0.025 31 140 40 −20 3.99 0.041 环境温度
T/℃水滴密度
ρl/kg·m−3空气密度
ρg/kg·m−3表面张力
σ/N·m水滴动力黏度
μl/10−3Pa·s空气动力黏度
μg/10−5Pa·s0 1000.00 1.293 0.07550 1.7921 1.72 −10 998.15 1.342 0.07710 2.6320 1.67 −20 993.60 1.395 0.07841 4.3300 1.62 -
[1] CAO Y H, TAN W Y, WU Z L . Aircraft icing: an ongoing threat to aviation safety[J]. Aerospace Science and Technology,2018 ,75 :353 −385 . doi: 10.1016/j.ast.2017.12.028[2] CAO Y H, XIN M . Numerical simulation of supercooled large droplet icing phenomenon: a review[J]. Archives of Computational Methods in Engineering,2020 ,27 (4 ):1231 −1265 . doi: 10.1007/s11831-019-09349-5[3] KULKARNI V, SOJKA P E . Bag breakup of low viscosity drops in the presence of a continuous air jet[J]. Physics of Fluids,2014 ,26 (7 ):072103 . doi: 10.1063/1.4887817[4] GARCÍA-MAGARIÑO A, SOR S, VELAZQUEZ A . Experimental characterization of water droplet deformation and breakup in the vicinity of a moving airfoil[J]. Aerospace Science and Technology,2015 ,45 :490 −500 . doi: 10.1016/j.ast.2015.06.025[5] SONI S K, KIRAR P K, KOLHE P, et al . Deformation and breakup of droplets in an oblique continuous air stream[J]. International Journal of Multiphase Flow,2020 ,122 :103141 . doi: 10.1016/j.ijmultiphaseflow.2019.103141[6] JALAAL M, MEHRAVARAN K . Fragmentation of falling liquid droplets in bag breakup mode[J]. International Journal of Multiphase Flow,2012 ,47 :115 −132 . doi: 10.1016/j.ijmultiphaseflow.2012.07.011[7] JALAAL M, MEHRAVARAN K . Transient growth of droplet instabilities in a stream[J]. Physics of Fluids,2014 ,26 (1 ):012101 . doi: 10.1063/1.4851056[8] KÉKESI T, AMBERG G, PRAHL WITTBERG L . Drop deformation and breakup[J]. International Journal of Multiphase Flow,2014 ,66 :1 −10 . doi: 10.1016/j.ijmultiphaseflow.2014.06.006[9] AALBURG C, VAN LEER B, FAETH G M . Deformation and drag properties of round drops subjected to shock-wave disturbances[J]. AIAA Journal,2003 ,41 (12 ):2371 −2378 . doi: 10.2514/2.6862[10] NI M J, KOMORI S, MORLEY N B . Direct simulation of falling droplet in a closed channel[J]. International Journal of Heat and Mass Transfer,2006 ,49 (1-2 ):366 −376 . doi: 10.1016/j.ijheatmasstransfer.2005.03.025[11] HAN J, TRYGGVASON G . Secondary breakup of axisymmetric liquid drops. I. Acceleration by a constant body force[J]. Physics of Fluids,1999 ,11 (12 ):3650 −3667 . doi: 10.1063/1.870229[12] HAN J, TRYGGVASON G . Secondary breakup of axisymmetric liquid drops. II. Impulsive acceleration[J]. Physics of Fluids,2001 ,13 (6 ):1554 −1565 . doi: 10.1063/1.1370389[13] JAIN M, PRAKASH R S, TOMAR G, et al. Secondary breakup of a drop at moderate Weber numbers[J]. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2015, 471(2177): 20140930. DOI: 10.1098/rspa.2014.0930
[14] STROTOS G, MALGARINOS I, NIKOLOPOULOS N, et al . Predicting droplet deformation and breakup for moderate Weber numbers[J]. International Journal of Multiphase Flow,2016 ,85 :96 −109 . doi: 10.1016/j.ijmultiphaseflow.2016.06.001[15] YANG W, JIA M, CHE Z Z, et al . Transitions of deformation to bag breakup and bag to bag-stamen breakup for droplets subjected to a continuous gas flow[J]. International Journal of Heat and Mass Transfer,2017 ,111 :884 −894 . doi: 10.1016/j.ijheatmasstransfer.2017.04.012[16] 杨威, 贾明, 孙凯, 等 . 液滴变形-袋式-多模式破碎转换研究[J]. 工程热物理学报,2017 ,38 (02 ):416 −420 .YANG W, JIA M, SUN K, et al . Investigation on transitions of deformation-bag-multimode breakup for liquid droplets[J]. Journal of Engineering Thermophysics,2017 ,38 (02 ):416 −420 . (in Chinese)[17] LIANG C, FEIGL K A, TANNER F X, et al . Three-dimensional simulations of drop deformation and breakup in air flow and comparisons with experimental observations[J]. Atomization and Sprays,2018 ,28 (7 ):621 −641 . doi: 10.1615/atomizspr.2018025948[18] XIAO F, WANG Z G, SUN M B, et al . Simulation of drop deformation and breakup in supersonic flow[J]. Proceedings of the Combustion Institute,2017 ,36 (2 ):2417 −2424 . doi: 10.1016/j.proci.2016.09.016[19] LIU N, WANG Z G, SUN M B, et al . Numerical simulation of liquid droplet breakup in supersonic flows[J]. Acta Astronautica,2018 ,145 :116 −130 . doi: 10.1016/j.actaastro.2018.01.010[20] YANG W, JIA M, SUN K, et al . Influence of density ratio on the secondary atomization of liquid droplets under highly unstable conditions[J]. Fuel,2016 ,174 :25 −35 . doi: 10.1016/j.fuel.2016.01.078[21] MARCOTTE F, ZALESKI S . Density contrast matters for drop fragmentation thresholds at low Ohnesorge number[J]. Physical Review Fluids,2019 ,4 (10 ):103604 . doi: 10.1103/physrevfluids.4.103604[22] JAIN S S, TYAGI N, PRAKASH R S, et al . Secondary breakup of drops at moderate Weber numbers: effect of Density ratio and Reynolds number[J]. International Journal of Multiphase Flow,2019 ,117 :25 −41 . doi: 10.1016/j.ijmultiphaseflow.2019.04.026[23] XIAO F, DIANAT M, MCGUIRK J J . A robust interface method for drop formation and breakup simulation at high density ratio using an extrapolated liquid velocity[J]. Computers & Fluids,2016 ,136 :402 −420 . doi: 10.1016/j.compfluid.2016.06.021[24] CHU G D, QIAN L J, ZHONG X K, et al . A numerical investigation on droplet bag breakup behavior of polymer solution[J]. Polymers,2020 ,12 (10 ):2172 . doi: 10.3390/polym12102172[25] SHAO C X, LUO K, FAN J R . Detailed numerical simulation of unsteady drag coefficient of deformable droplet[J]. Chemical Engineering Journal,2017 ,308 :619 −631 . doi: 10.1016/j.cej.2016.09.062[26] WANG Y, SHU C, YANG L M . An improved multiphase lattice Boltzmann flux solver for three-dimensional flows with large density ratio and high Reynolds number[J]. Journal of Computational Physics,2015 ,302 :41 −58 . doi: 10.1016/j.jcp.2015.08.049[27] LI W, WANG J X, ZHU C L, et al . Deformation and acceleration of water droplet in continuous airflow[J]. Physics of Fluids,2022 ,34 (3 ):033313 . doi: 10.1063/5.0085210[28] YANG L M, SHU C, YU Y, et al . A mass-conserved fractional step axisymmetric lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio[J]. Physics of Fluids,2020 ,32 (10 ):103308 . doi: 10.1063/5.0022050[29] LIU X D, OSHER S, CHAN T . Weighted essentially non-oscillatory schemes[J]. Journal of Computational Physics,1994 ,115 (1 ):200 −212 . doi: 10.1006/jcph.1994.1187[30] BIAN Q Y, ZHU C X, WANG J X, et al . Numerical investigation on the characteristics of single droplet deformation in the airflow at different temperatures[J]. Physics of Fluids,2022 ,34 (7 ):073307 . doi: 10.1063/5.0094748[31] SHU C W, OSHER S . Efficient implementation of essentially non-oscillatory shock-capturing schemes[J]. Journal of Computational Physics,1988 ,77 (2 ):439 −471 . doi: 10.1016/0021-9991(88)90177-5[32] LIANG C, FEIGL K A, TANNER F X . Axisymmetric simulations of drop deformation and breakup for validation of a modified Taylor analogy breakup model[J]. Atomization and Sprays,2017 ,27 (5 ):439 −455 . doi: 10.1615/atomizspr.2017017091[33] RIMBERT N, CASTRILLON ESCOBAR S, MEIGNEN R, et al . Spheroidal droplet deformation, oscillation and breakup in uniform outer flow[J]. Journal of Fluid Mechanics,2020 ,904 :A15 . doi: 10.1017/jfm.2020.675[34] ZHANG X, LIU X, WU X M, et al . Impacting-freezing dynamics of a supercooled water droplet on a cold surface: rebound and adhesion[J]. International Journal of Heat and Mass Transfer,2020 ,158 :119997 . doi: 10.1016/j.ijheatmasstransfer.2020.119997[35] CLIFT R, GRACE J R, WEBER M E. Bubbles, drops, and particles[M]. New York: Academic Press, 1978.
[36] WANG C, CHANG S, LENG M, et al . A two-dimensional splashing model for investigating impingement characteristics of supercooled large droplets[J]. International Journal of Multiphase Flow,2016 ,80 :131 −149 . doi: 10.1016/j.ijmultiphaseflow.2015.12.005[37] PAPADAKIS M, RACHMAN A, WONG S C, et al. Water impingement experiments on a NACA 23012 airfoil with simulated glaze ice shapes[C]//Proceedings of the 42nd AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada. Reston, Virigina: AIAA, 2004: AIAA2004-565. DOI: 10.2514/6.2004-565
[38] PAPADAKIS M, RACHMAN A, WONG S C, et al. Water droplet impingement on simulated glaze, mixed, and rime ice accretions, E-15277[R]. Washington, D. C. : NASA, 2007.