姜 寒，朱春玲

(南京航空航天大学航空学院，南京 210016)

飞机机翼作为提供升力的重要部件，在飞行过程中起着至关重要的作用[1]。当飞机在空中遇到结冰气象条件时，机翼表面会发生结冰现象，从而改变机翼的气动外形，导致升力下降阻力升高，影响飞行安全。因此对机翼结冰的研究就尤为重要，目前对机翼结冰研究的方法主要有结冰风洞实验和数值模拟两种。

结冰风洞实验主要通过将模型固定在结冰风洞内，通过风机和喷雾在制冷条件下模拟结冰气象条件下的飞行过程，观察结冰情况。

数值模拟利用计算机实现不同条件下的流场水滴场以及结冰的计算，因其应用的自由度和灵活性成为结冰研究中的主要方式。国外针对结冰计算开发了多款软件，例如 FENSAP-ICE[2]、LEWICE[3]和 ONERA[4]等。其中一些软件对结冰的模拟主要基于 Messinger 模型[5]。该模型通过求解表面积冰控制体的质量守恒及能量守恒方程，得到控制体的平衡温度，进而求得结冰量。Myers等[6]提出了能在结冰增长过程实现水膜、冰层之间传热传质的计算模型。中国对于结冰计算的研究工作也在发展中，张强等[7]基于欧拉两相流理论进行了三维机翼表面的霜冰模拟；曹广州等[8]根据冰层表面水膜流动特点建立了能模拟水膜流动的三维积冰数学模型；常士楠等[9]利用欧拉法求解空气水滴两相流流场并改进了 Messinger 模型；雷梦龙等[10]通过修改结冰类型判断方法，改进了 Myers 模型；王翊等[11]通过数值模拟的方法研究了二维翼型冰型生长过程；王强等[12]对二维风力机翼型进行了树脂及实验研究。可以看到，近年来中国学者对于机翼结冰的数值及实验研究越来越多，但大多数仍集中于二维翼型，对后掠翼结冰研究较少。

在建立数值计算模型过程中，现考虑水滴撞击在机翼表面先形成一层水膜，基于这层水膜进行结冰增长的计算，实现三维后掠翼的结冰模拟。设计后掠翼三维模型，在结冰风洞进行实验，与计算结果进行对比，验证计算模型的可行性，观察后掠翼积冰外形，得到后掠翼沿展向结冰的规律。

1 后掠翼结冰数值模拟

1.1 三维结冰模型的建立

首先，对模型进行如下假设。

(1)撞击在机翼表面的水滴会先在重力、剪切力及压力梯度的作用下向下游流动形成一层稳态的初始水膜，后结冰。

(2)结冰时忽略物面与气流间的辐射换热。

(3)未冻结的水在向下游流动时会形成连续的水膜且不发生破裂。

基于以上假设，建立初始水膜流动模型。建立贴体坐标系(ξ，η，ζ)，其中ξ、η为曲面上的参考坐标，ζ垂直于ξ、η。

机翼表面初始水膜流动过程中，通过控制体的质量包括：水滴撞击项mimp，水蒸发项mevp，通过控制体前后左右表面的流动项。对控制体应用质量守恒定理得到水膜连续性方程为

(1)

式(1)中：u、v分别为ξ、η方向的速度。

mimp=βV∞CLW

(2)

式(2)中：β为水收集系数；V∞为来流速度，m/s；CLW为液态水含量，kg/m3。

(3)

式(3)中：hc为对流换热系数，W/(m2·K)；ρa、cp、Rv分别为水蒸气密度、定压比热容及气体常数；esur、e∞分别为表面和大气饱和水蒸气分压力，Pa；Tsur、T∞分别为表面和大气温度，K；φ为大气相对湿度。

通过量纲分析法简化水膜动量方程为

(4)

式(4)中：P为空气压力；μw为水动力黏性系数；ρw为水密度；gξ、gη分别为重力在ξ、η方向的分量。

边界采用无滑移边界条件，即

(5)

式(5)中：τa为水膜表面剪切应力；Hw为水膜厚度。

基于上述初始水膜模型，建立结冰增长模型。对每一个控制体建立水的质量守恒方程为

mini+mimp+min-mout-mevp-mice=0

(6)

式(6)中：mini为初始水膜项，mini=ρwHw；min、mout分别为流入、流出项；mice为结冰项。

建立水的能量守恒方程为

qimp+qin+qcnd-qevp-qice-qcnv-qout=0

(7)

式(7)中：qimp为撞击项；qin、qout分别为流入、流出项；qcnd为导热项；qice为结冰项；qcnv为对流换热项；qevp为蒸发项。

各项热流计算式为

(8)

qin=min[ciTm+LF+cw(Tin-Tm)]

(9)

qout=mout[ciTm+LF+cw(Tsur-Tm)]

(10)

qice=miceciTsur

(11)

qcnv=hc(Tsur-T∞)

(12)

qevp=mevp{ciTm+LF+[frzci+(1-frz)cw]×

(Tsur-Tm)+LE}

(13)

式中：Tsur、T∞分别为表面温度和来流温度；Tw、Tin、Tm分别为水膜温度、入流水温及相变平衡温度，K；LF、LE、LS分别为冰的融解潜热、水的蒸发潜热、冰升华潜热，J/kg；cw、ci分别为水和冰的比热；frz为引入的冻结系数，表征冻结水的质量与进入控制体中水质量的比值，计算公式为

frz=mice/(mini+mimp+min)

(14)

1.2 三维结冰模型的求解

求解初始水膜模型时，先利用边界条件式(5)求解水膜动量方程(4)，可得到水膜沿法线方向的速度分布为

(15)

(16)

分别对式(15)、式(16)在η方向积分然后取平均，即可得到水膜沿翼面流动时分别沿两个方向ξ、η的平均速度为

(17)

(18)

将式(17)、式(18)代入式(1)，可得到关于水膜厚度Hw的偏微分方程为

(19)

(1)将所有控制体积中水的流入项赋初值为0。

(2)联立水的质量和能量守恒方程，代入式(14)，得到冻结系数的表达式，从而计算冻结系数frz，可得

frz={qclt+qin+qcnd-qevp-qcnv-(mini+mclt+min-mevp)[ciTm+LF+cw(Tsur-Tm)]}/{(mini+mclt+min)[(ci-cw)(Tsur-Tm)-LF]}

(20)

(3)用式(21)计算流出项mout，并根据水膜模型中计算出水流动速度分量按比例分配各方向的流出分量。

mout=(1-frz)(mini+mclt+min)-mevp

(21)

(4)根据流出分量重新计算各控制体的流入项并计算流出水的残差值e，公式为

e=max[mout(i)-mout(i-1)]

(22)

式(22)中：i为当前时间步；i-1为上一时间步。

(5)判断流出水项残差e是否满足收敛条件，若不收敛，重复上述步骤至收敛。

2 后掠翼结冰试验

2.1 实验模型与测量方法

设计后掠角为 30°的 NACA0012 机翼模型，使其能与结冰风洞支架实现连接，由于试验段宽300 mm，将模型展长定为 340 mm，使其在安装进试验段后与两边玻璃有一点空隙，NACA0012截面弦长定为 300 mm。利用 3D 打印技术加工该模型。进行实验时，将模型通过支架固定在结冰风洞试验段中，结冰试验结束后，将带冰的机翼模型连同支架一同从试验段中取出放在低温环境中以保证冰型不被破坏。结冰试验台如图1所示。

图1 结冰试验台

为准确测得实验所得结冰冰型，需要利用三维扫描仪先对未结冰的干净带支架翼型进行扫描，实验后在低温环境下对结冰带支架翼型进行扫描，将二者根据结冰实验前后完全没发生变化的支架上3个不同方向特征面进行重合处理，以得到二者相差部分即为机翼表面结冰的外形。

2.2 后掠翼结冰试验与结果

为研究机翼表面存在溢流水时的结冰情况，将来流温度控制在-8 ℃，工况如表1所示。

表1 结冰工况

结冰结果如图2所示，4个结冰截面厚度的扫描仪还原模型如图3所示。

图2 结冰试验结果

图3 扫描仪还原冰型

在还原模型上截取4个截面如图4所示，以便后期与结冰计算结果进行对比。图4(a)为截面z=60 mm处的结冰冰型，图4(b)为截面z=110 mm 处的结冰冰型，图4(c)为截面z=190 mm 处的结冰冰型，图4(d)为截面z=240 mm 处的结冰冰型。由图4可知：在攻角为 4°时，下翼面积冰较多且积冰范围较大，冰型轮廓较为光滑，无明显冰角凸起。上翼面积冰较少且积冰范围较小，有冰角，符合明冰结冰特征。由于结冰温度选取了温度较低的结明冰工况，因此冰角并没有典型明冰工况实验结果的冰角突出。图5为由翼根至翼梢各截面的结冰冰型对比。可以看出，沿展向方向由翼根至翼梢结冰冰型轮廓相似但结冰量逐渐增加，这是由于后掠角带来的气流展向流动，使水滴在撞击在机翼上形成的水膜除了弦向的流动外，还会沿展向流动，导致翼梢积冰量大于翼根。

图4 二维冰型截面图

图5 试验结果截面对比图

3 算例验证结果

为验证本文三维结冰模型的正确性，选择实验采用的后掠角为30°的NACA0012机翼作为计算模型，为与实验结果进行对比，对弦长为1 m的标准计算模型进行缩比至实验模型相同尺寸，弦长为0.3 m。其几何外形如图6所示。

图6 NACA0012 30°后掠机翼几何外形

结冰计算状态采用与实验相同的工况，如表1所示。

图7为结冰计算结果的三维视图，图中红色网格为干净翼面。沿翼型展向取4个截面，与FENSAP-ICE软件数值模拟结果与实验结果进行对比，选取的4个截面如图7所示。

图7 三维冰型图

图8为截面z分别为60、110、190、240 mm处的结冰冰型对比。从图8可以看出：本文模型结果与FENSAP-ICE结冰轮廓基本吻合。与实验结果对比可以发现，在靠近翼根的截面z=60 mm 处，本文模型求得的结冰量略小于实验结果，冰型轮廓基本符合，冰角形状不明显。在翼段中部的截面z=110 mm 和截面z=190 mm 处，本文模型计算结果与实验结果较为接近，对于冰角形状的计算略有不足。而在靠近翼梢的截面z=240 mm 处，本文模型求解的结冰量大于实验结果，冰角位置较为相似。

图8 截面冰型对比图

图9为数值计算结果各截面的冰型对比，从图9中可以看出，同实验结果相同，由翼根至翼梢，结冰量逐渐增长，但模型对冰角形状的计算不是很准确，原因可能是在计算结冰增长过程中，采用准稳态方法计算，虽考虑结冰初始水膜流动，但对其在结冰增长过程中由于表面形状变化导致的流动没有考虑，导致了结冰特征外形如冰角形状的不准确。

图9 计算结果截面对比图

4 结论

考虑水滴撞击在机翼表面形成的初始水膜，求解结冰增长过程的质量与能量守恒方程得到能够计算后掠翼结冰的数值方法。在结冰风洞中进行 30°后掠角的 NACA0012 三维模型结冰实验，利用三维扫描仪还原冰型。将计算结果与FENSAP-ICE 结果以及实验结果进行对比，得出以下结论。

(1)本文模型考虑了明冰工况下水滴撞击翼面形成的初始水膜流动，并将算得参数应用在结冰增长过程中，与考虑了结冰过程中水膜流动的FENSAP-ICE 软件计算结果符合良好。

(2)结冰风洞实验结果符合明冰工况结冰冰型，说明实验结果可用来验证计算。将计算结果与之对比可以发现基本吻合，在冰角处形状略有差异，证明本文计算方法可行。

(3)后掠翼模型在结明冰时可以看到由翼根至翼梢结冰厚度逐渐变大，说明后掠翼结冰过程中水沿展向的流动会对结冰外形产生明显的影响。