0 引言

飞机/发动机结冰会严重影响飞机性能和操稳特性，甚至危害飞行安全^[1-3]，因此防除冰系统设计非常重要。水滴撞击特性是部件防除冰系统设计的基础和关键^[4-5]，主要计算研究方法包括拉格朗日法^[6-8]和欧拉法^[9-11]。拉格朗日法是将过冷水滴当作离散系，对单个水滴建立控制方程，追踪水滴运动轨迹。欧拉法则是把过冷水滴视作连续流体，在欧拉坐标系下建立水滴相控制方程，使用同一套网格求解过冷水滴流场和空气流场，避免了拉格朗日法的大量插值问题，可获得各网格单元内水滴体积份额和表面撞击水滴的质量。对比研究^[12-14]表明二者都能在特定情况下很好分析物面水滴收集规律。

借助数值模拟工具，国内外对常规防冰系统和部件水滴收集规律开展了大量研究。赵勇等^[15]发现飞机三维水平尾翼从翼根到翼尖局部水滴撞击系数逐渐增加；易贤等^[16]发现三段翼最大收集率会出现在缝翼，而水滴撞击范围最大的是襟翼，T型尾翼平尾上水滴收集率要整体大于立尾。杨胜华等^[17]发现在多段翼后方存在水滴无法撞击的遮蔽区。杨军等^[18]发现随着飞行高度、来流速度及水滴直径的增加，水滴撞击极限、总收集系数和局部收集系均增加。申晓斌等^[19]发现进气道唇口附近水滴撞击区和撞击量都受进气流量的影响，水滴撞击特性受到整流罩影响较小。

鉴于冲击换热具有高的对流传热系数，常规热气防冰系统中通常利用热气射流冲击防冰部件内表面来提升加热效率，近年来也有研究人员直接利用热气射流除冰和防冰。例如，罗振兵等^[20-21]开展了合成热射流机翼除冰实验，发现合成热射流能有效促进热能的扩散，缩短除冰时间，出口通道开缝角度与表面积冰越垂直，加速除冰效果越明显。柯鹏等^[22]发展了一种针对复合材料帽罩的热气膜防冰系统，采用数值模拟和实验方法研究了其内部冲击换热和外部气膜加热效果，发现射流气膜能大幅提高缝后外壁面温度，在前缘开设气膜缝还会提升前缘内部的冲击换热效果。

外部射流直接加热待防护表面的同时还会影响近壁区的外部流场^[23]，改变物面水滴撞击特性。例如，柯鹏等^[24]计算发现气膜缝位置和热气吹风比等都影响开缝支板表面局部水滴撞击系数分布，吹风比越大，水滴撞击系数越小。刘华等^[25]计算了不同气膜缝出流位置、射流角度等结构参数下，支板表面水滴撞击特性，发现气膜缝位置会影响局部水滴撞击系数和总撞击效率。但目前射流对结构水滴撞击特性影响的研究仍十分有限，尤其是对于飞机/航空发动机的典型防冰部件，如机翼、进气道、帽罩和支板等。

本文选取了双开孔的三维圆柱作为基本构型，采用基于欧拉-欧拉方法的数值模拟算法对不同吹风比和来流速度下的水滴撞击系数进行了数值模拟，揭示射流对开孔圆柱表面不同位置水滴撞击系数分布、最大水滴收集系数和撞击极限的影响规律，并分析原因，为结冰预测与防冰设计提供参考。

1 考虑数值扩散的水滴相模型

一般防冰工况下可认为空气和水滴之间单向耦合，只考虑空气流场对水滴运动的影响。基本假设如下：

(1) 水滴运动过程保持常物性、无变形、无破碎、无碰撞；

(2) 水滴撞击壁面后无飞溅；

(3) 忽略水滴与空气之间的传热传质；

(4) 只考虑作用在水滴上的空气阻力和重力。

1.1 水滴相控制方程

Slater^[26]指出速度方向不同的水滴相交后形成一个无限大的密度脉冲，会导致欧拉法计算发散，为此建议在水滴连续性方程中加入数值扩散项，详细分析参见文献[27-28]等。

考虑数值扩散后水滴相控制方程如下：

$ \nabla \cdot\left(\rho_{w} \alpha \boldsymbol{V}\right)=b \nabla^{2} \alpha \rho_{w} $

(1)

$ \nabla \cdot\left(\rho_{w} \alpha \boldsymbol{V} \boldsymbol{V}\right)=\rho_{w} \alpha K\left(\boldsymbol{V}_{a}-\boldsymbol{V}\right)+\rho_{w} \alpha \boldsymbol{G} $

(2)

其中b为数值扩散系数，ρ_w为水滴密度，α为局部水滴体积分数，V为水滴速度矢量，V_a为气流运动速度矢量，G是重力加速度矢量。K为空气-水滴作用因子，

$ K=\frac{18 \mu_{a} f}{\rho_{w} D_{w}^{2}} $

(3)

其中μ_a是空气动力黏度，D_w是水滴直径，f为水滴阻力系数。

1.2 边界条件

基于假设(2)，撞击到壁面的水滴应从计算域中排除。设n为撞击壁面的单位法向量，n · V < 0即表明水滴撞击到了壁面。壁面水滴容积分数和速度都取紧贴壁面网格的中心点的值。

1.3 局部水滴撞击系数计算

局部水滴撞击系数β定义为微元表面上的实际水收集量与该微元表面上最大可能的收集量之比，计算式为：

$ \beta=\alpha(\boldsymbol{V} \cdot \boldsymbol{n}) /\left(\alpha_{\infty}\left|\boldsymbol{V}_{\infty}\right|\right) $

(4)

其中α_∞为来流体积分数，V_∞为来流速度。

定义吹风比M来表征射流强度，

$ M=\rho_{j} \boldsymbol{V}_{j} / \rho_{\infty} \boldsymbol{V}_{\infty} $

(5)

其中ρ_j为射流气体密度，V_j为射流气体速度。

1.4 算法验证

选取文献[29]中三维圆柱算例，圆柱直径101.6 mm，来流速度80 m/s，压力89867 Pa，水滴直径16 μm。计算边界采用速度入口，压力出口。图 1给出了中心截面上水滴撞击系数对比结果。可以发现，本文计算结果中水滴撞击系数分布与文献中数值模拟和实验结果吻合良好，撞击极限和实验结果相比误差大约为6.5%，最大水滴撞击系数和实验结果相比误差大约为1%。更多验证结果见文献[14, 28]。

2 射流对水滴撞击特性分布的影响 2.1 几何模型及网格

研究构型为开直孔的三维圆柱，几何模型和计算域如图 2所示，直径D取76 mm。射流孔为上下对称分布，水平方向，直径d为7 mm，孔间距为15.4 mm。计算域为前场8D，后场14D，上下各7D。其中图 2(d)为局部放大图，图 2(a)、图 2(b)和图 2(c)分别图 2(d)的正视图、左视图和俯视图。

计算网格为三维结构化网格，通过网格无关性验证后，最终网格数约300万，如图 3所示，其中近壁网格进行了加密，y⁺≈1。

2.2 边界条件及求解器设置

来流入口采用速度入口边界条件，出口采用压力出口边界条件，射流采用速度入口边界条件，展向两侧采用周期性边界条件，圆柱表面采用1.2节定义的边界条件，其他壁面为无滑移壁面边界。

利用Fluent软件求解空气流场，水滴流场通过用户自定义函数求解。采用标准k-ε湍流模型和增强型壁面函数，使用基于压力的分离式求解器，各项参数均采用二阶精度的空间离散方法。

2.3 计算结果及分析

选取典型水滴大小(单一直径20 μm)在来流速度(30 m/s)和射流吹风比(M=1.0)条件下进行计算，发现射流对水滴撞击特性存在显著影响，无射流孔圆柱和有射流孔圆柱的表面水滴撞击系数分布对比如图 4所示。

对比图 4(a)、图 4(b)可见，射流孔出口下游水滴被完全吹散，小孔附近则随表面位置不同而呈现出不同变化规律。为此，进一步选取几个特征截面来分析不同位置水滴撞击特性。各截面对应位置关系如图 5所示，展向从左往右依次为截面L1~L6，中心截面L6通过两个射流孔中线，弦向中心截面H1则为二个射流孔上下方向的正中间截面。各个截面的无量纲位置坐标如表 1所示，其中，s表示截面到二孔中心线的距离，r为射流孔半径。

各弦向截面在吹风比为1.0时的水滴撞击系数分布如图 6所示, 射流对水滴撞击系数分布规律影响较为复杂。例如，中心截面L6的水滴撞击系数显著小于无气膜情况，撞击极限终止于孔边缘，最大水滴撞击系数降低了6.88%。近孔截面L5的水滴撞击系数分布曲线呈现显著不同的凹曲线，在前缘的一小段范围内曲线保持水平，之后呈现为快速下降的凹曲线，到撞击边缘附近时降幅趋缓。水滴撞击系数变化最剧烈的区域(距离中心-0.015~-0.008 m和0.008~0.015 m)正是射流吹袭作用最强的区域。最大水滴撞击系数相比于无气膜时减少了5%，撞击极限减少了8%。中间截面L4上最大水滴撞击系数和撞击极限都略小于无气膜情况，但是水滴撞击系数曲线从撞击极限到驻点前缘(距离中心-0.030~-0.012 m和0.012~0.030 m)几乎呈线性分布，与无射流情况下差异不大。在前缘一小段范围内(距离中心-0.012~0.012 m)曲线几乎保持水平。

为分析水滴撞击系数变化的原因，将图 6中圆柱表面沿弦向划分为三个区域，Ⅰ区、Ⅱ区和Ⅲ区，对应表面的孔间、孔两侧和孔后区域，图 7进一步给出与三个区相应的水滴速度矢量图。

Ⅰ区所有截面的水滴撞击系数均小于无气膜作用情况，从图 7中Ⅰ区水滴速度矢量图可以明显看出，水滴有向展向和弦向外侧运动的趋势，这是在射流展向扩散和弦向吹袭共同作用下的结果。

1) 展向射流扩散作用可以从图 8的H1剖面上来流水滴轨迹和流场y方向涡量场可以明显看出，两射流孔间的水滴被吹向展向外侧，造成Ⅰ区的来流水量减少。

2) 弦向射流吹袭作用，由图 9中L6和L1截面吹风比为1.0时空气运动轨迹和水滴运动轨迹可以看出，水滴轨迹有明显的上扬，被吹向弦向后侧，而且对比图 9(b)、图 9(c)发现，近流孔截面受到射流吹袭作用更强，所以L6截面上水滴撞击特性变化更为剧烈。

Ⅱ区为射流孔侧区域，在气流吹袭和扩散作用下(如图 7中Ⅱ区水滴速度矢量图)下，水滴撞击系数会下降。近孔截面L5和L6则受射流吹袭作用影响更大，下降最为显著，而外侧截面L1和L2则影响不大；但是，对于L3和L4截面，由于射流扩散作用将Ⅰ区的水滴吹过来，导致来流水量增加，所以水滴撞击系数反而会增大。

Ⅲ区为射流孔后的区域，与Ⅱ区受到的作用力类似，即在射流吹袭和射流扩散的共同作用下，水滴向孔后、向外侧方向运动(如图 7中Ⅲ区水滴速度矢量图)，射流影响范围仅局限在近孔截面L4~L6。

3 射流强度影响分析

射流强度会影响圆柱表面的水滴撞击特性分布规律，为此，进一步分析了吹风比(0.5和1.5)对典型截面位置的水滴撞击系数的分布特征和最大水滴撞击系数的影响。

图 10给出了不同吹风比条件下表面水滴撞击系数的分布规律。总体来说，吹风比对不同位置的水滴撞击系数的影响规律类似，孔后基本无水滴，孔间和近孔位置受影响最为显著，差异也很明显。

当吹风比减小到0.5时，由图 10(a)可以看到，总体分布规律与吹风比为1.0时基本一致，是射流吹袭和射流扩散共同作用下的结果。由于吹风比小，射流孔气流吹袭影响仅限于近孔侧截面L5和孔间截面L6。对其他截面的影响较小，基本限于孔间Ⅰ区范围。

当吹风比增大到1.5时，图 10(b)显示了射流的显著影响。外侧截面L1和L2的水滴撞击系数都明显高于无气膜情况，其他截面变化与吹风比1.0的情况类似，但是孔间Ⅰ区的水滴撞击系数也明显增加，高于无气膜状态。原因在于当吹风比大于1.0时，射流吹袭和扩散作用更加强烈，来自两孔间(Ⅰ区)的水滴沿展向被吹向外侧，造成近孔截面L6水滴量减少，小于无气膜结果，而远孔各截面上增加，在不同程度上高于无气膜情况。

各截面水滴撞击系数变化规律的差异，主要来源于物面局部流场的巨大差异。图 11给出了吹风比为1.5时的物面涡量分布情况，近壁区域存在射流诱发的复杂涡流场，会对水滴撞击轨迹产生显著影响。接下来针对不同截面来具体分析。

图 12为吹风比为1.5时L6和L1截面水滴运动轨迹。对比图 12(b)和图 9(b)可明显看出吹风比增大时，中心截面L6处射流吹袭作用加强，减小了局部水滴撞击系数。综上，在近孔截面射流吹袭作用和远孔截面射流扩散作用的综合影响下，会出现图 10(b)中的水滴撞击系数分布规律。

展向位置上的最大水滴撞击系数的分布结果能更好的显示射流在不同区域的综合作用。图 13给出了三个吹风比下, 截面H1的水滴撞击系数。吹风比较小时(0.5和1.0)，总体水滴撞击系数分布低于无气膜工况，孔间区域水滴撞击系数相比于无气膜工况最大下降了7.51%，吹风比增大到1.5时，射流动量大大增强，对孔间主气流吹袭作用更为强烈，导致孔间水滴撞击系数比无气膜工况最大下降了8.95%。但水滴撞击系数从展向外侧截面向孔间截面呈下降趋势，最大降幅约17.65%。

4 结论

采用基于欧拉-欧拉模型的水滴撞击算法，针对典型水滴直径，分析了射流对三维带双开孔射流的圆柱表面水滴撞击特性的影响规律和内在机理，得出以下结论：

1) 射流与主流之间存在的强耦合作用会显著改变来流的流动规律，进而会影响开孔圆柱物面的水滴撞击特性。在物面的不同区域，射流与水滴之间存在吹袭作用和扩散作用等不同的作用关系，且其强度随物面区域变化而差异很大。例如，吹风比为1.0时，弦向孔后区域影响最大，完全没有水滴撞击，展向向外区域影响逐渐减小。

2) 射流和来流水滴间的吹袭和扩散作用不仅与物面位置有关，与射流强度也密切相关。水滴直径为20 μm，在小吹风比(0.5和1.0)时，射流仅对展向近孔截面有影响，孔间水滴撞击系数相比于无气膜工况最大下降7.51%；吹风比为1.5时，射流吹袭和扩散作用进一步加强，孔间的水滴沿展向被吹向外侧，造成近孔截面水滴量减少而远孔截面增加，孔间截面水滴撞击系数比无气膜工况最大下降了8.95%，水滴撞击系数从展向外侧向孔间最大降幅约17.65%。

本文初步研究了射流对三维开孔圆柱物面不同区域水滴撞击特性的影响规律，分析了其内在机理。虽然所得规律在相同工况下对多种结构具有一定普遍性，但仅针对了特定的水滴直径(20 μm)和开孔布置方式，而射流对水滴撞击特性的影响程度会随水滴直径减小进一步增强，具体规律还有待后续进一步研究。