单股射流倾斜撞壁影响因素的数值模拟

戴青雯,杨建文,杨 斌,王 莹

(1.上海理工大学 能源与动力工程学院 上海市动力工程多相流动与传热重点实验室,上海 200093;2.西安航天动力研究所,陕西 西安 710100)

0 引言

倾斜射流撞击固体表面设计简单,并且生产成本低,广泛应用于燃料雾化[1]、液体火箭发动机[2]的液膜冷却[3-4]、射流清洗[5-7]等领域。在液体火箭发动机工作过程中,核心区域温度高达3 000 ℃,为避免燃烧室壁面被高温热气流烧蚀,常采用液膜冷却的方式来降低传递到壁面的热量,达到保护壁面的作用[8]。在液膜冷却设计中,发动机推力室内冷却液膜的覆盖面积是影响冷却性能的重要指标,因此有必要对液膜形态的影响因素开展深入的研究。

文献[9]的研究表明,射流角度、射流速度、流体物性及壁面接触角等因素皆会影响射流撞壁后的液膜形态,例如液膜呈现圆形水力跳跃、辫状流、撞壁后反弹及先铺展后破碎等形态。针对倾斜射流撞壁后的液膜流动这一典型气液两相流问题,国内外很多学者通过开展一系列的实验来研究液膜的形态变化规律。文献[10]针对疏水和超疏水壁面开展了射流撞壁的实验研究,发现接触角增大,而液膜铺展面积减小,同时还对比分析了壁面接触角、射流倾角及韦伯数对射流反射角的影响。Good等探究了不同射流倾角对撞击后亲水表面液膜分布的影响,揭示了射流速度与液膜最大宽度和壁面液膜溅射量之间的关系[11]。唐亮等在开展单股圆柱射流撞壁的实验研究的基础上,得出了射流撞壁后液膜最大宽度、液膜长度等几何参数的表达式[12]。林庆国等开展了射流撞击曲壁的实验研究,采用探针法测量液膜厚度,探究了射流角度和壁面曲率半径对液膜厚度和液膜形态的影响[13]。然而,射流撞壁后液膜的铺展过程复杂且实验获取液膜几何参数的成本较高。

在实验研究的基础上,国内外很多学者通过建立半经验模型开展液膜变化规律的理论研究。Wang等建立了半经验模型来分析射流速度、射流角度、流体黏度和表面张力等因素对液膜关键特征的影响,并对射流撞壁开展了一系列实验,得出的半经验模型计算结果与实验结果相符合[14]。Hasson等提出两股射流对撞的流动滞止点位于对撞后形成的椭圆形液体薄片的一个焦点上,并在不考虑边界层理论和黏性作用的情况下建立了液膜厚度分布的理论模型[15]。Watson考虑边界层理论,建立了射流垂直撞击壁面后层流和湍流两种情况下的液膜厚度分布模型[16]。更进一步地,Inamura等[17]考虑到层流边界层的发展,对Hasson等[15]的模型进行修正后建立了适用于射流撞击曲壁的液膜厚度分布模型。文献[18]在考虑黏性作用的基础上建立了单股圆柱射流撞击壁面后的液膜铺展模型,在实验验证模型可靠性的基础上研究了表面亲疏水性对液膜形态的影响。为进一步深入研究液膜的形态变化,文献[19]还建立了两股流体射流的理论分析模型,研究了射流角度、韦伯数和雷诺数对液膜形态的影响。然而,传统经验模型目前还无法开展下游流动复杂处的壁面压力和液膜变化规律的研究。

相比于实验和理论模型预测,倾斜射流撞壁在数值仿真方面的研究相对较少。对于垂直射流撞击壁面,Gradeck等开展了不同射流速度下撞击运动表面的实验和数值仿真研究,通过对比数值仿真结果与实验结果得到的水跃位置,验证了数值模型的准确性[20]。在倾斜射流撞击壁面方面,Fard等运用三维数值模拟研究了工业喷嘴射流撞壁后液膜的形成及破碎成液滴过程,并将仿真得到的液膜厚度与实验结果进行对比,验证了数值模型的可靠性[21]。林庆国通过数值仿真得到壁面上的液膜分布,研究了不同壁面曲率和不同压降对液膜表面形态的影响,且只定性研究了液膜的基本形态,并未定量对比数值仿真结果与实验结果[22]。邱添开展了射流冲击平板的全三维和准三维数值模拟研究,由于其全三维模拟研究的网格不够细密,使得模拟得到的液膜形态存在失真现象[23]。上述文献并未对不同影响因素,例如不同射流速度下壁面的压力变化和不同孔径下液膜表面速度变化开展相关研究。

对于单股射流撞壁后形成的液膜,可以用液膜长度、宽度和厚度等关键几何参数来描述。但当前在液膜几何参数方面的定量研究较少,一是实验获取液膜几何参数的耗时长、成本高,二是传统经验模型会忽略黏性、重力等因素对液膜的影响导致不能精准捕捉液膜的基本形态。因此,本文采用数值模拟方法定量研究射流撞壁后的流场区域。首先,比较实验和数值模拟所得的液膜宽度,验证了本文采用的VOF多相流模型及数值模拟方法的可靠性;其次,分析液膜形成过程中的关键特征;最后,定量研究射流倾角、射流速度及射流孔径对液膜形态的影响。

1 数值求解方法

1.1 控制方程

在本文射流撞壁数值仿真的研究中,研究对象为水和空气。在计算中,流体的物性设为常数,因此连续性方程和动量方程为

(1)

(2)

式中:ρ为流体密度;v为流体速度;t为时间;p为压强;μ为动力黏度;g为重力加速度;f为动量源项。

在数值模拟过程中,气液两相的体积分数相加和为1,引入流体体积比函数φ的概念,φ表示流体体积与网格单元体积的比值[24]。φ=1时,表征区域内只有液态水;φ=0时,表征区域内只有空气;0<φ<1时,表征区域内既有液态水也有空气。

因此,密度和动力黏度可用流体体积比函数φ进行插值计算,即

(3)

(4)

式中:ρl和ρg分别为液体和气体的密度;μl和μg分别为液体和气体的动力黏度。

1.2 几何模型

在倾斜射流撞壁的过程中,液态水从圆形喷孔以一定的角度射出,到达壁面后液态水先在平板表面向四周铺展并继续向出口处流动,然后在表面张力的作用下收缩,形成一个类似椭圆形的液相区域,最后汇聚成一股流动从出口流出,射流撞壁后形成的液膜如图1所示[14]。

图1 液膜图像Fig.1 Picture of the liquid sheet

由于液态水射流撞壁铺展成液膜的全过程是关于喷孔中沿轴线对称的,计算域模型采用半模型,可以减小计算量;同时液膜厚度很小,在不影响液膜形态的前提下可以减小计算域模型在液膜厚度方向上的尺寸。为了充分展示流体从射流入口以一定的速度v沿着射流倾角θ方向流出后液膜的铺展情况,故计算域大小为9 mm×51 mm×6 mm,如图2所示。图2中,y轴正方向为重力方向,射流倾角θ为速度与y轴之间的夹角,同时液膜主流方向与y轴正方向保持一致。

图2 计算域三维示意图(单位:mm)Fig.2 Three-dimensional schematic chart of computational domain(unit:mm)

1.3 计算方法

选择基于有限体积法的Fluent求解器计算,采用的VOF模型是一种基于欧拉法的多相流表面跟踪方法[25],具有易守恒和占用内存小的优点,适用于相界面变化复杂的液体流动计算。使用VOF模型计算倾斜射流撞壁时,在多相流模型设置中将空气和液态水分别设置为第一和第二相,初始计算域内均为空气。

将射流入口设置为速度入口,壁面设置为无滑移壁面,与壁面的接触角设为73.3°,其他各个面都设置为压力出口。物性参数设置参考水的物性参数,表面张力系数设为72.8 mN/m,液体的密度设为998.2 kg/m3,黏性系数设为0.001 Pa·s。湍流模型选择RNGk-ε模型,采用Menter-Lechner壁面函数处理近壁面流动,压力离散选用PRESTO算法,压力与速度耦合关系选用SIMPLE算法,时间步长为10 μs,每时间步长最大迭代步数为40。

1.4 网格划分和无关性验证

采用多面体网格划分方法,并且为了更好地捕捉液膜形态,在壁面附近设置边界层,如图3所示。通过4套不同尺寸的网格进行网格无关性验证,网格数量分别为54万、99万、115.5万和175万。参考文献[26]取射流孔径为0.8 mm,设置进口速度为5 m/s,射流倾角为20°,以y=6 mm、y=12 mm和y=20 mm这3个截面处的液膜宽度为判断依据。

图3 计算域网格及壁面边界层网格局部放大图Fig.3 Computational domain mesh and partial magnified drawing of wall boundary layer mesh

图4为在射流孔径为0.8 mm、速度为20 m/s和倾角为20°的工况下,射流撞击固体表面后形成的液膜形态及液膜关键几何参数示意图。射流撞壁后可分为自由射流区、薄膜区、边缘水跃区及液膜汇集区。

图4 射流撞壁后的液膜形态及液膜关键几何参数示意图Fig.4 Diagram of liquid film shape after jet impingement and key geometric parameters of liquid film

图4还展示了液膜的关键几何参数液膜长度l和液膜宽度w的提取方式。液膜长度为沿着流动方向从壁面上开始出现液膜到液膜第一次汇集成单股流动的位置;液膜宽度为液膜两侧边缘在某一截面上的铺展距离。

网格无关性验证结果如图5所示。当网格数量为54万时,这3个截面上的液膜宽度与其他3套网格液膜宽度的相对误差超过5%,而99万、115.5万和175万3套网格之间计算出来的液膜宽度相对误差均低于5%。为兼顾计算精度和成本,最终应用115.5万的网格数量进行后续的仿真计算。

图5 网格无关性验证结果Fig.5 The result of mesh independence verification

1.5 模型验证

当液态水以射流的形式撞击固体壁面,由于惯性力的作用,液态水以撞击点为中心向周围扩散,并形成两个突出的边缘,即边缘水跃区;在边缘水跃区之间形成的一层薄液体层称为薄膜区;在液膜达到最大宽度后,由于表面张力克服惯性力的作用,边缘水跃区开始向内收缩流动,最终在液膜汇集点变成一股流体继续向下游流动;在液膜汇集点之后的区域称为液膜汇集区,至此,液膜汇集区上游的液膜稳定发展,液膜铺展的关键特征不再发生变化。图4中观察到的现象与文献[14]实验结果相似,证明了本文模型设置的准确性。

为了进一步验证数值仿真结果的可靠性,将数值仿真结果与实验结果进行定量对比。当液膜稳定后,在y=6 mm、y=12 mm和y=20 mm这3个截面上提取液膜宽度,并与文献[26]实验得到的液膜宽度进行对比,如表1所示。水射流撞击平板后形成液膜,液膜与壁面之间存在一定的相互作用,而在本文数值仿真过程中简化了液膜与壁面之间温度变化和壁面接触角变化的相互作用关系,这使得仿真得到的液膜宽度略小于实验值。但数值模拟结果与实验结果之间的误差均在20%以内,低于文献[26]仿真结果的误差。因此认为本文的数值仿真结果可靠,可采用本文的数值模拟方法进行射流撞壁的相关研究。

表1 不同截面的液膜宽度对比Tab.1 Width comparison of liquid film in different cross sections

2 结果与讨论

2.1 射流倾角对液膜的影响

图6为倾斜射流撞壁原理的局部示意图。

图6 倾斜射流撞壁原理的局部示意图Fig.6 Partial schematic diagram of oblique jet impingement principle

Kate等通过理论预测和实验验证分析了射流撞壁后的流场变化,推导得出射流撞壁形成液膜的滞止点位置[27]。液态水从喷孔射出后,在与壁面平行的自由射流区截面A-A上呈椭圆形,椭圆的焦点S1在分离流线上,分离流线与射流中心线平行,且在分离流线右侧射流向液膜汇集区方向铺展,左侧射流向反方向铺展,将分离流线与壁面的交点称为滞止点,如图6中S所示。

在不同的射流倾角下,液膜的铺展方式相同。首先液膜在铺展过程中,沿着液膜铺展方向,液膜宽度先增大后减小,但由于不同射流倾角下的液膜在各个方向上的流动分速度不同,使得液膜宽度和液膜厚度有所差异。

取滞止点所在截面为y=0 mm截面,图7是射流倾角θ分别为20°、15°和10°时沿y方向不同截面位置的液相体积分数分布云图。由图7可知,射流倾角越大,液膜的宽度越大。虽然液膜倾角越大,在各个截面上的铺展面积明显增大,但是倾角越小,液膜的能量和受力越集中,在下游越不容易发生破碎。

图7 沿y轴不同截面的液相体积分数分布云图Fig.7 Liquid phase volume fraction in the y-direction at different cross sections

图8对比分析了不同射流倾角下y=4 mm、y=10 mm、y=16 mm及y=22 mm这4个截面上的液膜厚度变化。射流倾角越大,不同横截面上薄膜区与水跃区之间的液膜厚度变化越剧烈。这是由于射流倾角越大,液膜轴向分速度越小,流动越慢,横向分速度越大,横向扩展面积越大。此外,由于在不改变射流流量情况下改变倾角,液膜在薄膜区与边缘水跃区厚度变化越剧烈,液膜稳定性越差,且过大的射流倾角会使边缘水跃区在横向不断向外扩展,在轴向不断向内收缩,从而可能会导致后续液膜两侧破碎现象的产生。

图8 不同射流倾角下的液膜厚度变化Fig.8 Thickness variation of liquid film at different inclination angles

结合图7和图8分析可知:当液膜完全铺展并达到稳定后,在薄膜区,液膜厚度从中心轴线向两侧减小;当液膜铺展到边缘水跃区时,厚度突然增大。因此,可认为在液膜汇集区上游各个截面上液膜整体呈中间薄两侧厚的分布。

射流倾角越大,薄膜区铺展面积越大,边缘水跃区的厚度也越大;同时液膜的宽度越大,撞击点与液膜汇集点之间的距离越长。文献[23]中也观察到同样的现象,这是因为不同射流角度下流量在各个方向上分配不一致导致的,射流倾角越大,y方向流量越大,这也解释了不同入射条件下撞壁后液膜轮廓不同的原因。

2.2 射流速度对液膜的影响

射流孔径和射流倾角一定时,不同射流速度会影响液膜在流动方向上的延伸距离和液膜横向铺展方向上宽度的变化。表2是不同射流速度下液膜的长度和最大宽度的对比。

表2 不同射流速度下液膜最大宽度和液膜长度Tab.2 Maximum width and length of liquid film at different inclination velocity

图9是不同射流速度下的液膜形态,虽然射流速度不同,但描述液膜关键特征的自由射流区、薄膜区、边缘水跃区及液膜汇集区的轮廓和边界线都清晰可见。结合表2和图9分析可得,射流速度越大,液膜的最大宽度越大,长度也越大。这是因为射流速度越大会导致射流总流量和动量越大,液膜的铺展面积越大。

图9 不同射流速度下的液膜形态Fig.9 Liquid film shape at different inclination velocities

液态水撞击壁面后铺展时,壁面会受到水的压力。忽略壁面粗糙度、黏性等参数的影响,由于力的作用是相互的,该壁面压力值也能反映液膜的压力分布情况。图10是不同射流速度下在0～150 Pa范围内壁面的压力分布。红色区域表示此处壁面的压力高于150 Pa,蓝色区域表示此处壁面的压力低于0 Pa。

图10 不同射流速度下壁面的压力分布Fig.10 Wall pressure distribution at different inclination velocities

从图10中可以看出,液态水在刚撞击壁面时的压力较高,之后在薄膜区壁面压力迅速降低,趋近于环境压力;边缘水跃区则是外边界压力最大,从水跃区外边界到内边界压力递减,外侧压力大能抑制薄膜区的持续铺展,推动液膜向内侧收缩;而在薄膜区和边缘水跃区内边界之间的分界处,压力值小于0 Pa,认为该处压力低于环境压力。同时,在液膜汇集区之后的壁面会出现高压,文献[9]在进行超疏水壁面倾斜射流时的研究也发现液膜汇集区之后的壁面存在高压区。因此,推测本文数值模拟的亲水表面出现这一现象也是因为倾斜射流是由固体壁面支撑的,壁面压力与射流的切向动量密切相关,液膜在汇集区发生堆积,之后在重力作用下回落到壁面并继续向下游流动,这就使得汇集区之后的切向动量增大,从而汇集区之后的壁面出现高压区。

图10对比分析了不同射流速度时壁面的压力分布,速度越大,撞壁后压力高于150 Pa的面积越大,可认为射流撞壁区域越大;且速度越大,撞壁处液膜的惯性力越大,使得液膜在薄膜区铺展面积越大,液膜宽度越大;当速度越小时,靠近液膜汇集区的水跃区压力越大,且速度越小惯性力越小,表面张力越大,较高的表面张力推动着液膜收缩,使得液膜长度越短,但液膜的稳定性会有所提高。

2.3 射流孔径对液膜的影响

当保持射流倾角θ=20°和射流速度v=5 m/s时,射流孔径越大,液态水的质量流量也越大,质量流量的差异导致液膜的长度、宽度和厚度发生改变。图11是不同射流孔径下液膜表面的速度分布,在撞击点附近流速最大,使得液膜向周围铺展,而随着液膜向下游运动,速度不断减小,在水跃区汇集点附近速度达到最小。此时由于重力和气体回流的影响,液膜下游初步呈现向内侧聚拢的现象。

图11 不同射流孔径下液膜表面速度分布Fig.11 The velocity distribution of liquid film surface at different inclination diameter

不同射流孔径下液膜宽度和长度变化如表3所示。由表3数据分析可知,射流孔径越大,质量流量越大,液膜铺展面积越广,液膜的宽度和长度也越大。

表3 不同射流孔径下液膜最大宽度和液膜长度Tab.3 Maximum width and length of liquid film at different inclination diameters 单位:mm

图12是不同射流孔径下中心轴线上的液膜厚度分布。从图12中可以看出,从撞击点处开始随着液态水不断地向下游流动,液膜厚度先逐渐减小再逐渐增大。发生这一现象的主要原因是,一开始射流撞击壁面,较大的速度使得液膜迅速在壁面铺展,液膜在撞击点附近由厚变薄,之后由于壁面粗糙度及回流等因素的影响,流动速度变得缓慢,液膜逐渐开始收缩堆积,最后形成液膜汇集区。在这一过程中,中心轴线处的液膜厚度不断增加,甚至在液膜汇集区的最大厚度与薄膜区的最小厚度相差一个量级。结合表3和图12分析可知,射流孔径越大,液膜最大宽度越大,液膜长度越长,液膜厚度越大。这与文献[28]中理论模型计算得到液膜轮廓结论一致。表3中射流孔径从0.6 mm增大至0.8 mm时,液膜最大宽度和长度的增幅分别为37.6%和37.7%,而孔径变化对射流速度在各个方向的分量影响几乎可忽略,可推测孔径变化对液膜稳定性的影响较小。

图12 不同射流孔径下中心轴线上液膜厚度变化Fig.12 Thickness variation of liquid film on central axis at different inclination diameters

3 结论

本文使用VOF多相流模型并结合多面体混合网格,对多工况下倾斜射流撞壁进行了数值模拟,并得出以下结论。

1) 本文数值仿真得到的液膜宽度与实验值变化趋势保持一致,且两者之间误差均在20%以内,表明本文构建的数值模型适用于倾斜射流撞壁的数值模拟研究。

2)在液膜宽度方向,液膜厚度整体呈中间薄两侧厚的分布;在液膜主流方向,从撞击点处开始,中心轴线上液膜厚度先逐渐减小再迅速增大;射流孔径越大,液膜的中心轴线上最大厚度也越大,当射流孔径为0.6 mm时,液膜最大厚度约为1.3 mm,当射流孔径为0.8 mm时,液膜的最大厚度超过2 mm。随着射流倾角增大,液膜的宽度增大而长度减小;随着射流速度和孔径的增大,液膜的宽度和长度都增大。

3) 射流撞壁后,壁面压力从液膜的边缘水跃区到中心轴线先逐渐减小后增大再减小,且在薄膜区和边缘水跃区内边界之间的分界处,压力值小于0 Pa。在撞击点附近,射流速度越大,液膜的惯性力越大,使得壁面受到的压力越大,而在液膜汇集区,表面张力起主导作用,表面张力增大时惯性力减小,因而壁面受到的压力越小。

4) 在液膜汇集点之前,沿着液膜长度方向,液膜的流动速度先增大后减小;液膜在撞击点附近速度最大,近似等于射流速度,在水跃区汇集点附近速度最小;当射流孔径从0.6 mm增大到0.8 mm时,液膜的最大质量流量增大,体现在液膜最大宽度从5.24 mm增大到8.4 mm,液膜长度从19.55 mm增大到31.55 mm,液膜铺展面积越广,液膜最大速度和最小速度区域的面积越大。