飞机除冰冲击射流多相流耦合模型研究

龚 淼,申远航

(中国民航大学 航空工程学院, 天津 300300)

0 引言

飞机除冰效果和效率是影响航班起飞安全和准点率的重要环节。近年来我国冬季冰雪天气增多,降雪范围增大,除冰作业压力逐年增大。目前民航局已经在北京大兴机场和首都机场推广了慢车除冰模式,对提升除冰作业效率产生了积极的影响。在慢车除冰作业中,除冰冲击射流在飞机表面形成的速度场与压力场对除冰效果有很大影响,研究冲击射流冲击在飞机表面的压强与流速大小可为缩短慢车除冰融冰时间和优化作业参数提供重要参考。由于飞机除冰研究具有尺寸大、耦合的物理场数目多、涉及研究范围广等难点,很多国内外学者都使用数值模拟法来研究除冰液冲击射流。因此将除冰液冲击射流的多相流模型进行对比分析,选择出适合研究飞机除冰冲射流的多相流模型可以有效提高研究效率,在飞机除冰研究中具有重要的工程意义。

国内外学者用来研究冲击射流的平台包括CFX、Ansys、OpenFOAM、Comsol等。Castillo[1]使用CFX对大坝的溢流瀑布冲击射流进行了研究,得到了大坝自由落体冲击射流的空气夹带程度与射流厚度、破裂形状、冲击底部压力的关系,仿真结果与实验结果具有良好的一致性。此外还发现空气夹带程度对CFX的网格具有较高的敏感性,而改变仿真中的湍流模型只会对冲击射流滞点出的速度压力值有影响。Faiza等[2]使用了CFX模块研究了湍流液体射流流入静止水池中形成的流场结构,通过建立两相流-气泡流模型,使用粒子图像测速法监测水池表面以下区域的时均速度场。此外还用相间动量传递模型预测了下落流体的平均流量,预测结果与实验值具有较高的一致性。Bel等[3]使用能量方法预测气泡爆炸的影响,如果使用微射流来估计气泡爆炸产生的相关压力场,得到的结果更容易被PIF(压力强度函数)和EPF(侵蚀功率函数)响应。Mahdavi等[4]使用Ansys-fluent研究了纳米射流喷射到圆形转盘之后的传热特性与流动状态,通过流体体积模型将两相间张力的作用体现出来,并且发现纳米流体的体积分数越大,传热所需的扭矩就越大。Wang等[5]利用OpenFOAM模拟了深海起重机控制有效系统载荷受水动力冲击的场景,其中6 DOF刚体运动求解器约束条件的引入保证了结果的准确性,研究发现波粒速度和波斜率是圆柱体上压力不对称的基本因素。周洪利等[6]利用OpenFOAM研究了竖直管内复杂的气固液三相流动问题,提出了有效的水合物颗粒聚并破碎模型组合,发现了管内压力损失梯度与流速的变化关系。Islas等[7]使用OpenFOAM模拟了粉尘爆炸中粉尘云的变化情况。该模型预测了粉尘云中双涡流模式的形成,使用较大的颗粒可以增加颗粒惯性,改善粉尘云的均匀性。Elaissi等[8]通过Comsol建立了一个湍流等离子喷射模型用来帮助合成硅纳米粉末。研究考虑了纳米粒子通过对流、扩散和热泳的运输,给出了等离子射流流场以及粒子的轨迹和热传导过程。周忠会等[9]基于Comsol研究了海底沉积物的热流场分布,通过近似原位的电导热加热实验,给沉积物模型加热并且探究了加热后的沉积物热流场分布,这项研究对预测和评估海底天然气水合物资源量有重要意义。Comsol支持多种求解方法包括有限元法、边界元法、VOF法等,可以有效计算多物理场耦合的复杂流动模型,比较适合用来研究飞机除冰过程中的冲击射流流场参数变化。

在冲击射流方面,大量研究结果表明圆形出口的冲击射流满足射流的一般特性。根据冲击射流的流场形状和速度可以将冲击射流分为3个明显的部分:自由射流区、冲击区和壁面射流区。在自由射流区,射流扩展半径随射流中心轴线长度增长而增大。在射流冲击区,射流与冲击壁面接触,发生明显弯曲,存在很大的压力梯度。在壁面射流区,流动速度方向发生显著变化与冲击壁面垂直[10]。贝尔陶斯和拉贾拉南设计了大量实验,结合量纲分析发现了冲击射流滞点压力与出射速度、射流高度、射流半径[11]的定量关系。Zhou等[12]研究了泵射流水动力性能、噪声性能,涉及空化侵蚀和尖端间隙涡流特性,转子与定子及尾流场的相互作用。研究发现,紧凑的结构和复杂的内外部流场使泵获得了更高的性能。

本文以飞机除冰喷洒作业为研究对象,基于二维的N-S方程,根据Comsol中的多相流模型模拟了飞机除冰作业中的除冰冲击射流,得到了冲击射流在冲击壁面上的速度场与压力场分布,并且根据速度与压力参数确定了混合物模型和相传递混合模型更适合做除冰冲击射流的仿真研究,为飞机除冰控制系统的研究工作提供了参考。

1 冲击射流几何及数值模拟

1.1 射流冲击几何模型

本研究以飞机除冰喷洒作业为研究对象,如图1(a)所示。为了更加明显地对比出不同多相流模型对模拟飞机除冰冲击射流的模拟结果,将实际冲击射流模型简化,选取圆柱形喷枪的中心截面为研究平面,模拟垂直冲击角度下所取二维研究平面内气液两相流场的速度与压力分布情况,如图1(b)所示。在Comsol中建立冲击射流平面模型,根据除冰车喷嘴的实际参数确定模型尺寸,其中窄边为除冰枪出口,宽边为除冰枪内壁。为匹配实际除冰过程,选取1 m×1 m的冲击空间,计算域如图1(c)所示。根据行业标准AHM975规定,飞机除冰车喷枪出口温度范围为60～85 ℃,除冰车喷洒压力约为9 bar,除冰作业时喷洒流量约为200 L/min,因此设置喷嘴直径0.012 m,长0.05 m,数值模拟的流场入口为ab,出口选定ce和df,冲击壁面为ef,流体属性设置为标准大气压下85 ℃水的物性参数。

图1 飞机除冰车喷嘴及计算域划分图

1.2 多相流模型的控制方程

本研究以相场模型、相传递混合模型、混合物模型、水平集模型、气泡流模型这5个多相流模型为基础,分别与Realizablek-ε[13]湍流模型耦合,添加稳态研究,探究稳态情形下的冲击射流速度场与压力场分布情况。

Realizablek-ε是标准k-ε[14]模型经过修正之后得出的一种湍流模型,更加适合用来解决流速变化大,流场形状复杂的流动场景。Realizablek-ε的求解基于动量守恒的RANS方程、质量守恒的连续性方程、壁附近的壁函数以及k-ε输运方程。

基于动量守恒的RANS方程:

ρ(u·▽)u= ▽·[ -ρI+K]+F

(1)

连续性方程:

▽(ρu)=0

(2)

壁面无滑移条件:

u·n= 0

(3)

式(1)中ρ是流体密度,u是速度场矢量,▽是二维平面下的拉普拉斯算子,F是包含重力在内的体积力。其中K是粘性应力张量。

K= (μ+μT)(▽u+ (▽u)T)-

(4)

湍流动能k的输运方程:

(5)

湍流动能耗散率ε的输运方程:

(6)

(7)

(8)

在式(4)和式(5)中ep是湍流耗散率;式(7)中μT是湍流粘度;Cμ是随应变率变化的量,简化计算可以取0.09;式(8)中Pk是生成项,Cε1取1.44,Cε2取1.92,σε取1.3,σk取1.0。Realizablek-ε与标准k-ε使用相同的默认初始条件和缩放比例。

相场模型使用式(9)追踪相界面移动情况,其中λ是混合能量密度,γ移动参数。

(9)

相传递混合模型使用式(10)控制粘度,其中μd和μc是连续相和分散相的粘度。

(10)

水平集添加了式(11)跟踪相界面,方程表征的是相界面的移动情况。其中γ是重新初始化参数,一般取1,ε是界面厚度控制参数,φ是水平集变量根据相的不同取1或者-1。

(11)

气泡流模型结合了气液两相的动量方程和连续性方程,其动量方程为:

(12)

1.3 网格划分及边界条件

为了最大程度降低网格对不同流动模型仿真的影响,本研究的5个多相流模型都使用同等级的物理场控制的网格。改变网格最大单元尺寸,对比网格质量如表1所示。

表1 不同最大网格单元大小下的网格参数

由表1可以看出,改变网格尺寸,网格质量变化不大,保持在0.82～0.83,能较好地满足计算精度。最终选取中间值,即总网格数为13 855,平均网格质量为0.828 1的网格划分方案。网格划分结果如图2所示,网格质量分布如图3所示。

图2 计算域网格划分图

在图3中,颜色越趋于蓝色的网格质量越好,可以看出本研究的模型网格具有相对较好的网格质量。网格划分部分使用了最常用的自由三角形网格,能够在满足计算收敛性的前提下减少计算量。为了防止边界处的湍流变量和流体流动变量的不收敛,在湍流模型的壁面边界进行了网格加密。此外,本研究旨在探究不同多相流模型在计算除冰液射流问题中的适用性,对不同模型采用了相同的网格划分办法,因此不再进行网格无关性验证。

为保证计算准确性和稳定性设置了以下边界条件:使用速度条件控制入口,根据行业标准AHM975将入口速度设置为30 m/s,法向流入;左右出口使用一个标准大气压力约束;液气两相流不可压缩且包含重力作用;在喷嘴内壁和其他无热流交换的壁面设置为无滑移条件;使用了P1+P1和一阶线性的离散化方案,即用一阶函数的办法求解速度场和压力场。此外还根据雷诺数指定了统一的中等湍流强度[15]。

2 结果对比分析

2.1 冲击射流速度场结果对比分析

不同多相流模型的冲击射流速度分布云图如图4所示。观察速度云图可以发现喷嘴入口处的速度最大为设置值30 m/s,从喷口处往下,射流在涡旋卷积作用下逐渐扩张,截面逐渐增大。速度沿中心流线逐渐减小,在射流到达冲击壁面时速度急剧减小到0,而后在冲击壁面上形成稳定的壁面射流[16]。

图4 冲击射流多相流模型速度分布云图

然而观察速度云图可以发现相场模型和水平集模型在冲击壁面处发生了偏转,这与实际情况不符。研究2个模型的控制方程可知,相场模型和水平集模型都是基于场的模型,液相与气相之间互不相容,在相邻计算网格之间,相场函数φ[17]会在-1和1之间跳跃,故2种模型对相界面的描述较精确,但是难以准确模拟存在大量流体颗粒或者微团的分散型流场。综上在图4(a)和图4(d)中射流发生偏转是由于射流在破碎时形成了大量流体微团,上述2个模型无法准确描述这种分散相,导致射流两侧流场分布不对称,随着射流的不断发展形成了偏转的现象。

接下来探究冲击壁面上的速度场分布情况,以冲击壁面表面为研究对象,获取不同多相流模型表面的速度值绘制曲线如图5所示。

图5 冲击壁面速度大小分布曲线

不同模型形成的壁面冲击射流具有大致相同的分布模式,在滞点[18]位置速度最小,相场和水平集模型由于其偏转的程度最高,射流冲击到壁面时与壁面存在夹角,动能没有完全转换为压力能,最小速度分别为为0.36 m/s和0.46 m/s。混合物模型、气泡流模型以及相传递混合模型在滞点处的动能都完全转化为压力能,最小速度为0。

壁面射流从滞点位置向两侧发展,压力能逐渐转换为动能,流动速度不断上升。在出口处采用了标准大气压力约束,出口处是自由边界,在壁面流动阻力的作用下,射流速度逐渐降低。不同模型在偏离滞点0.15～0.2 m的位置逐渐获得最大流动速度。混合物模型的壁面最大流速为6.84 m/s,气泡流模型的壁面最大流速为9.05 m/s,水平集模型的壁面最大流速为8.28 m/s,相场模型的壁面最大流速为6.84 m/s,相传递混合模型的壁面最大流速为6.67 m/s。其中气泡流的最高速度比其他模型的最高速度都大,原因是气泡流模型忽略了分散相的动量作用,流动过程中动能耗散更少,因此壁面射流的最高速度也最大。相场模型与水平集模型由于形成了偏转,在滞点撞击壁面时损失的能量更少,所以在壁面射流区获得的最高速度大于相传递混合模型与混合物模型。

除冰作业的实际效果与滞点和最大速度点的位置有密切的关系。在今后的研究中,还将继续探究慢车除冰作业的运动流场中,滞点位置的运动方式对除冰效果的影响。

2.2 冲击射流压力场结果对比分析

不同多相流模型的冲击射流压强分布图如图6所示。观察压强等值线图,发现不同流动模型模拟出的压强分布大致相同。在入口处流体靠压力推动向前喷射,压强达到最大值。随着冲击射流离开喷嘴,射流与大气接触,压力能迅速释放,流体压强降低到大气压力。在滞点处,动能又迅速转换为压力能,滞点处获得冲击壁面的最大压强。

图6 冲击射流多相流模型压力分布图

气泡流的最大压强为5.54×104Pa,明显高于其他4个多相流模型的最大压强。在2.1节中分析得出气泡流模型在模拟冲击射流的过程中忽略了分散相的动量作用,降低了能量损耗,导致其壁面射流最高速度远远高出其他4个模型,所以气泡流模型在动能与压力相互转化的过程中也能获得最大的压强。

在冲击射流发展段的两侧都形成了“眼状”的等值线,从中心处向外逐渐减小到大气压。但是混合物和水平集模型的压强等值线比其余3个模型稀疏,推测可能是其余3个模型模拟的射流左右两侧形成了涡旋,涡旋改变了周围的压强分布梯度。

通过检测射流两侧的涡流强度来验证猜想。取模型空间坐标系坐标为(-0.25 m,-0.5 m)处为监测点,获取该点处不同多相流模型的涡流值。相场模型在该点处的涡流大小为37.742 s-1,相传递混合模型在该点处的涡流大小为6.86 s-1,气泡流模型在该点处的涡流大小为20.521 s-1,水平集和混合物模型在该点未检测到涡流。因此可以认为,在入口速度为30 m/s的情况下,相场、气泡流、相传递混合模型在射流两侧形成了涡旋,造成了更大的压强梯度。

2.3 冲击射流的紊流系数大小对比分析

冲击射流的紊流系数是射流结构的重要表征参数。一般来说出口处的射流紊流系数与紊流强度相关,紊流强度越大,紊流系数越大;一些关于垂直纯水射流的研究发现当紊流系数取0.067时,实验结果与理论计算结果具有较高的一致性[19]。为了进一步对比不同多相流模型对飞机除冰射流的模拟效果,将不同模型研究下的冲击射流紊流强度进行对比。

在自由射流区,射流扩展半径随轴线长度呈线性增长,射流边界与中心轴线形成张角α。一般定义射流张角α与紊流系数a的关系为:

tanα=ka

(13)

k是实验系数,由实验过程中的射流喷口形状确定,本研究中由于出口相比于计算域而言较为狭窄,在工程中一般取2.44。因a是在理想情况下自由冲击射流的紊流系数,本研究考虑了重力作用,故取距离喷口0.2 m以内的射流部分读取其张角α。不同模型的射流边界图如图7所示。

图7 不同模型下的冲击射流边界图

由图7可以读取不同模型下的冲击射流边界夹角,由大到小分别为8.5°、7.7°、7.25°、6.5°、6.25°。再根据式(13)得到紊流系数a,不同模型下的冲击射流紊流系数如表2所示。

从表2中可以看出,混合物模型的紊流系数与实验的紊流系数较为接近,水平集模型、相传递混合模型的计算结果与实验值有一定差距,而相场模型、气泡流模型的计算结果与实验值相差最多,效果最差。

3 结论

1) 冲击射流的截面沿着中心轴线呈锥状发展,在滞点处流速衰减为0。不同模型在偏离滞点0.15～0.2 m的位置逐渐获得最大流动速度。混合物、气泡流、水平集、相场以及相传递混合模型在30 m/s的入射速度下于壁面获得的最大流速分别为6.84、9.05、8.28、6.84、6.67 m/s。壁面的最大流速点到出口之间可以形成比较稳定的壁面射流。喷嘴入口处流体压强最大,随着冲击射流离开喷嘴,流体压强迅速降低到大气压力。在滞点处机械能相互转化获得冲击壁面的最大压强。

2) 相场模型和水平集模型在冲击壁面处发生了偏转,是由于冲击射流在破碎时形成了大量流体微团,2个模型的相场函数φ会在-1和1之间突变,模型不能准确描述发生剧烈拓扑变化的流动,导致射流两侧流场分布不对称,随着射流的不断发展出现了偏转。

3) 通过冲击射流的紊流系数大小对比可以发现,混合物模型、水平集模型、相传递混合模型的紊流系数与实验值更加接近,模拟效果更好。综合考虑模型模拟出的速度、压力分布情况以及射流结构发现:对于飞机除冰射流以及类似的大流速、多流体微团类的流体分析,采用混合物模型和相传递混合模型可达到更加精确的研究结果。