结冰风洞中液滴相变效应数值模拟

郭向东，王梓旭，李明，刘蓓

1.中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000 2. 中国空气动力研究与发展中心 飞行器结冰与防除冰重点实验室，绵阳 621000

飞机结冰广泛存在于飞行实践中，并严重威胁飞行安全[1-2]。研究飞机结冰的主要途径包括数值模拟、结冰风洞试验和飞行试验3种，其中结冰风洞试验顾名思义是利用结冰风洞在地面模拟飞机结冰过程，验证飞机防除冰系统，开展飞机结冰问题研究[3]。相对于数值模拟和飞行试验，结冰风洞试验具有结冰条件易控、试验成本相对较低、结果可靠等优点，是目前研究飞机结冰的主要手段。为模拟高空低温低压的云雾环境，结冰风洞主要包括制冷系统、高度模拟系统和喷雾系统三大部分，其中喷雾系统利用喷嘴在风洞稳定段内生成大量悬浮液滴，随后液滴由低温气流携带进入试验段，形成试验所需的过冷水滴云雾[4]。结冰试验中，气流静温低于冰点，为防止喷雾系统中液态水结冰，须采用加热方式提高液态水温度，这样便导致稳定段内生成的液滴温度高于气流静温，通常情况下，液滴会与气流发生充分的热交换，进而在试验段内达到过冷状态(液滴温度与试验段静温一致)，形成试验所需的过冷水滴[5]。但是，由于水滴与气流间的热交换过程十分复杂，涉及对流传热、相变传质和三维收缩效应等多物理过程的耦合作用，试验参数的变化可能导致试验段内液滴偏离过冷状态，进而改变积冰生长过程和冰形特征，影响结冰风洞试验结果，因此探索结冰风洞中液滴传热过程，揭示试验参数的影响规律，评估液滴过冷状态，对结冰风洞试验准确模拟结冰过程具有重要意义。

针对结冰风洞液滴过冷问题，Willbanks和Schulzt[6]率先发展了基于拉格朗日法的运动液滴传质传热计算方法，并对发动机高空试验台内的液滴传热过程开展了研究，随后Miller等[7]采用这一方法，研究了NASA Glenn IRT结冰风洞内的液滴过冷问题，发现初始液滴温度对液滴过冷状态无显著影响，此外Bellucci等[8]采用试验方法对CIRA结冰风洞高速试验段和低速试验段构型的液滴过冷状态进行了评估。但是这些研究主要关注工程层面上液滴过冷状态的评估，缺乏对结冰风洞中液滴传热传质过程的规律性认识。

近年来，随着国内大型结冰风洞——3 m×2 m结冰风洞的建成，液滴过冷问题逐渐引起了一些国内学者的关注[9-10]。其中郭向东等[10]发展了基于Euler法的气液两相耦合流动计算方法，研究了3 m×2 m结冰风洞主试验段构型液滴传热过程，并评估了液滴过冷状态，研究结果表明结冰风洞中液滴传热过程可以分为准一维传热和三维收缩传热两个阶段，其中三维收缩效应对液滴过冷状态起决定性作用，但是他们的研究中未考虑液滴相变传质对液滴传热过程的影响，导致对结冰风洞液滴过冷状态的评估结果较为保守。

因此，本文在郭向东等的计算方法基础上，进一步引入Hill矩方法[11]，发展了基于Euler法的气液两相传质传热耦合流动计算方法，实现了液滴相变传质过程的模拟。进而利用该方法，针对3 m×2 m结冰风洞主试验段构型，考察了液滴相变效应对液滴传热过程的影响，揭示了相变效应的影响规律，开展了参数影响研究，分析了相对湿度、试验段气流速度和液滴尺寸对液滴传质传热过程的影响，评估了液滴过冷状态，为国内飞机结冰地面模拟试验数据应用评估提供支撑。

1 计算方法

1.1 控制方程和物理模型

本文借鉴了文献[10]中基于Euler法的气液两相耦合流动计算方法，进一步耦合Hill矩函数方程组[12-14]，实现了相变传质过程的模拟。为获得简化物理模型，对气液两相进行以下假设[15]：①气相为理想气体，遵循理想气体法则；②液滴为球形，不发生变形和破碎，液滴间不发生碰撞；③液滴内温度均匀分布；④忽略重力效应；⑤ 黏性和热传导效应仅发生在气液两相间；⑥不考虑液滴冻结过程以及水蒸气凝结成核过程。

基于以上假设，简化后的控制方程组为

1) 气相方程组

(1)

(2)

(3)

(4)

2) 液相方程组

(5)

(6)

(7)

3) 矩函数方程组

(8)

(9)

(10)

式中：下标g和d分别代表气相和液相;ρ、p、V和E分别为密度、压力、速度和单位质量总能量；气相能量Eg和液相能量Ed的表达式分别为

(11)

式中：cv为气相定容比热；cl为液滴比热；T为温度。

σ和α分别为液滴有效密度和水蒸气质量分数，其中σ对应液态水含量(LWC)；n阶矩函数定义为

(12)

其中：r为液滴半径；f为粒子尺寸分布函数；Q0、Q1和Q2为零阶、一阶和二阶矩函数，根据定义可见Q0表示单位体积内液滴总个数。

S0为质量源项，表示为

(13)

式中：ρl为液滴物理密度；dr/dt为液滴半径变化率。

SD和SM分别为两相动量交换源项和两相间传质引起的动量传递源项，表示为

SD=Q0FD,SM=S0Vd

(14)

式中：FD为单个液滴受到的气相作用力。

SEg和SEd分别为气相和液相能量源项，表示为

(15)

式中：qheat与qrad分别为单个液滴与气相间的对流传热率和辐射传热率；qgas,mass和qliquid,mass分别为单位质量传质引起的气相和液相能量变化率，表示为[16-17]

(16)

其中：cp为水蒸气定压比热；L为汽化潜热。

S1和S2分别为一阶和二阶矩函数方程源项，表示为

(17)

单个液滴受到的气相作用力FD和单个液滴与气相间的对流传热率qheat表示为[18]

(18)

qheat=Nukgπd(Tg-Td)

(19)

式中：d为液滴直径；μg为气相黏性系数；kg为气相热传导系数；Red为液滴相对雷诺数，表示为

(20)

CD为液滴阻力系数，针对球型液滴，采用Schiller-Naumann圆球阻力模型，表示为

(21)

针对球型液滴，Nu采用Ranz-Marshall半经验关系式，表示为

(22)

式中：Pr为Prandtl数。

单个液滴与气流间的辐射传热率qrad表示为[19]

(23)

式中：ε为热辐射发射率，根据文献[19]，ε取0.95;σSB为Stefan-Boltzmann常数。

(24)

式中：DAB为二元质量扩散系数；Rv为水蒸气气体常数；pv,s为液滴表面饱和水气压；Hu为相对湿度，表示为

(25)

式中：pv,∞为气相环境中水蒸气分压。

Sh为Sherwood数，针对球型液滴，Sh采用Ranz-Marshall半经验关系式，表示为

(26)

式中：Sc为Schmidt数。根据文献[16，21]，Ranz-Marshall半经验关系式适用于Red<137、温度低于0 ℃时球型液滴的传质传热过程，匹配本文的研究范围(尤其对于高风速工况下依然成立)，因此本文计算模型中Nu和Sh采用该关系式。

1.2 数值方法

本文采用与文献[10]相同的数值方法，控制方程组采用有限体积方法离散，非定常项采用二阶Euler格式离散，空间输运项采用二阶迎风型的Roe格式离散，源项采用隐式格式[10]。

气液两相边界条件采用不同形式处理：气相入口条件采用压力入口，出口条件采用压力出口，壁面采用无黏滑移壁面条件(1.1节假设⑤)；液相入口条件采用速度入口，出口采用出流条件，壁面处认为液滴可以直接穿过，不考虑壁面对液相的作用。

1.3 验证算例

选取文献[21]中的试验结果验证本文计算方法。文献[21]采用悬挂液滴试验方法研究了液滴冻结过程中的液滴温度变化特性，文中将液滴冻结过程分为过冷(Supercooling Stage)、复辉(Recalescence Stage)、冻结(Freezing Stage)和冷却(Cooling Stage)4个阶段，其中过冷阶段对应本文研究的内容，因此选取文献中过冷阶段的试验结果验证本文计算方法。验证结果如图1所示，图中：Va和Ta分别表示气流速率和气流静温。从图中可以看出：本文计算结果与试验结果吻合得较好，不同工况下液滴温度误差均小于1%，因此本文计算方法可以合理模拟液滴过冷阶段的传热过程。针对以上的验证过程，应该指出的是考虑到悬挂液滴法对液滴尺寸和试验风速的限制，该方法仅获得了大液滴(毫米量级)低风速(小于1 m/s)下的试验结果，而缺乏小液滴(微米量级)高风速的试验结果。进一步，文献[21]指出，尽管缺乏真实喷雾的试验数据，但是试验数据中液滴相对雷诺数范围(36

2 结果与分析

本节针对3 m×2 m结冰风洞主试验段收缩构型，采用本文发展的计算方法，考察了基准工况下液滴相变效应对液滴传热过程的影响，揭示了相变效应的影响规律，开展了参数影响研究，分析了相对湿度、试验段气流速度和液滴尺寸对液滴传质传热过程的影响，评估了试验段中心处液滴过冷状态。

3 m×2 m结冰风洞主试验段构型如图2所示，该构型包括稳定段、收缩段和试验段3部分，其中构型入口截面为高8 m×宽11 m，位于稳定段喷雾耙处，出口截面为高2 m×宽3 m，位于试验段中心处，构型全长16 m，收缩比为14.67。

根据文献[10]的研究结果可知，试验段气流速度和液滴尺寸是影响液滴过冷状态的主要参数，因此本文选取相对湿度、试验段气流速度和液滴尺寸3个参数开展参数研究，计算工况矩阵如表1所示，表中LWCi、di、Tdi、T0i、Udi、UTS和Hui分别代表初始液态水含量、初始液滴直径、初始液滴温度、初始气流总温、初始液滴速率、试验段气流速度、初始相对湿度。Case 1为基准工况，考察典型工况下结冰风洞中相变对液滴传热过程的影响；Case 2选取了70%、80%、90%和100%这4个典型的相对湿度状态，考察相对湿度的影响；Case 3选取了55 、77、130、164 m/s这4个典型试验段风速，考察气流速度的影响，其中各风速对应的入口总压分别为2 000、4 000、12 000、20 000 Pa，参考压力为105Pa；Case 4选取了40、60、80、100 μm这4个典型液滴直径，考察液滴尺寸的影响，值得指出的是由于假设液滴为球型，不考虑液滴变形和破碎等动力学特性，因此仅选取了小尺寸液滴开展研究。

表1 计算工况矩阵Table 1 Calculation condition matrix

2.1 基准工况

基准工况选取FAR Part25 附录C中典型结冰条件[22]，为考察相变对液滴传热过程的影响，下文通过对比构型内液滴和颗粒的传热过程，进而揭示液滴相变效应的影响规律，其中针对颗粒传热过程的模拟，本文采用文献[10]中的计算方法来实现。

图3给出了基准工况下构型中心线处气流静温和液滴温度变化曲线，从图中可以看出，类似于文献[10]中颗粒传热过程，整个液滴传热过程近似以构型6 m处截面为界分成准一维传热和三维收缩传热两个阶段。在准一维传热阶段内(入口至6 m处截面，包括稳定段和收缩段入口区域)，风洞型面变化较缓，气流以准一维方式低速流动，气流速度和静温变化较小，接近入口流速和总温，液滴以准一维稳态方式传热；在三维收缩传热阶段内(6 m处截面至出口，包括收缩段后部区域和试验段)，风洞型面变化较大，三维收缩效应较强，气流场在此阶段内发生了显著变化，气流速度加速增大，气流静温则等熵下降，出口气流速度越大，出口气流静温越低，三维收缩效应(气流可压缩性)显著影响液滴传热过程，液滴以三维动态方式传热。对比液滴与颗粒的温度变化曲线可见：在准一维传热阶段，液滴温度下降趋势快于颗粒温度，当液滴温度和颗粒温度趋于稳定时(约在0.6 m处)，液滴温度(约-5.4 ℃)低于颗粒温度(-5 ℃)，此时液滴温度称为液滴湿球温度[16]；在三维收缩传热阶段，液滴温度下降趋势慢于颗粒，最终在构型出口处(试验段中心)，液滴温度(-7.41 ℃)高于颗粒温度(-7.77 ℃)。

图4给出了构型中心线处液滴直径和相对湿度变化曲线，从图中可以看出：结冰风洞中液滴经历了先蒸发后凝结两个阶段，在准一维传热阶段，气流中水蒸气未达到饱和状态，液滴处于蒸发状态，液滴直径不断减小，同时蒸发产生的水蒸气则使气流的相对湿度不断增大；在三维收缩阶段，气流静温迅速下降，气流饱和水汽压降低，而水蒸气分压变化较小(适用于低速流动)，导致气流的相对湿度迅速增大(对应式(25))，气流中水蒸气从未饱和状态转换为过饱和状态，液滴则从蒸发状态转换为凝结状态。

进一步理论分析相变效应对液滴传热的影响。首先液相能量方程式(7)可以化简为

(27)

式中:由于液滴过冷阶段辐射传热远小于对流传热和相变传热，因此忽略辐射传热项，同时联立式(19)和式(24)，则式(27)可化简为

(28)

ξheat+ξmass

(29)

式中：ξheat为对流传热引起的液滴温度空间变化率；ξmass为液滴相变引起的液滴温度空间变化率。需要指出的是，式(29)是在一维定常流动条件下得到的，但是对于三维收缩阶段，构型中心线处的气流沿流向加速流动，仍可以近似认为气流运动为准一维，因此式(29)适用于构型中心线处整个液滴传热过程的分析。

图5给出了构型中心线处ξmass和ξheat的变化曲线，从图中可以看出，在准一维传热阶段，当温度较高的液滴进入未饱和的低温气流(Hui<100%)中时，液滴与气流间发生对流传热，液滴向气流传递热能(ξheat<0)，同时蒸发效应从液滴吸收所需的汽化潜热(ξmass<0)，而颗粒与气流间仅发生对流传热过程，因此液滴温度下降速度快于颗粒的速度(0～0.6 m)；当液滴温度下降至气流静温时[23]，由于蒸发效应仍持续从液滴吸收热量，因此液滴温度会持续下降，但对于颗粒而言，其温度此时将会与气流静温保持一致，颗粒与气流间的对流传热过程停止，颗粒达到稳定状态；进一步，当液滴温度下降至液滴湿球温度时，蒸发效应所需的汽化潜热(ξmass<0)与气流向液滴传递的热能(ξheat>0)达到平衡(|ξmass|=|ξheat|)，液滴温度趋于稳定，此时ξmass与ξheat沿流向(x轴)对称分布(0.6～6 m)，液滴进入稳定蒸发阶段。在三维收缩阶段，三维收缩效应导致气流静温迅速下降，液滴稳定蒸发过程被破坏，ξmass与ξheat不再对称分布，液滴温度下降(ξmass+ξheat<0)；此阶段内对流传热的能量传递方向发生了转换，即ξheat>0转换为ξheat<0，促进了液滴温度下降，同时气流中水蒸气从未饱和状态转换为过饱和状态，液滴则从蒸发状态转换为凝结状态，即ξmass<0转换为ξmass>0，则抑制了液滴温度下降，而颗粒与气流间仅存在对流换热过程，缺少了凝结过程的抑制效应，因此可见凝结效应是三维收缩阶段内液滴温度下降趋势慢于颗粒的主要原因，最终导致构型出口处(试验段中心)，液滴温度高于颗粒温度。

2.2 相对湿度影响

首先考察相对湿度的影响，在结冰风洞试验中，为避免气流湿度过低导致液滴大量蒸发，会采用加湿器提高风洞回路内气流湿度，因此此处选取较高的相对湿度参数开展研究。

图6(a)和图6(b)给出了不同初始相对湿度下构型中心线处液滴直径和相对湿度变化曲线，从图中可以看出：在准一维传热阶段，随着初始相对湿度的增加，液滴尺寸的减小趋势减慢，表明液滴蒸发效应减弱，尤其当湿度100%时，液滴尺寸在温度稳定阶段不再发生变化，蒸发效应消失；与此相反，在三维收缩阶段，随着初始相对湿度的增加，水蒸气过饱和程度增强，液滴尺寸的增大趋势加快，表明凝结效应增强。

图6(c)给出了不同初始相对湿度下构型中心线处液滴温度变化曲线，从图中可以看出：在准一维传热阶段，随着初始相对湿度增加，液滴蒸发效应减弱，液滴温度下降趋势减慢(如放大图所示)，当液滴趋于稳定蒸发状态时(约在0.6 m处)，初始相对湿度越大则液滴湿球温度越高，尤其当相对湿度为100%时，蒸发效应消失，液滴温度与颗粒温度趋于一致；在三维收缩传热阶段，随着初始相对湿度增加，液滴凝结效应增强，导致构型出口处(试验段中心)液滴温度不断增大，尤其当相对湿度100%时，液滴凝结效应最强，出口液滴温度最高，比颗粒温度高近1 ℃。

图7给出了不同初始相对湿度下构型出口中心处液气温度差，从图中可以看出：随着初始相对湿度增加，液气温差不断增大，其中初始相对湿度为70%时，液滴温差与颗粒的一致，约为0.1 ℃，而当初始相对湿度增加至100%时，液滴温差增至约0.87 ℃，比颗粒温度差高近0.8 ℃。由此可见，相对湿度会影响试验段内液滴过冷状态，但影响程度较弱，本文计算工况下最大液气温差(相对湿度为100%)小于1 ℃，液滴仍认为处于过冷状态。

2.3 试验段气流速度影响

考察试验段气流速度的影响，根据3 m×2 m结冰风洞典型试验工况，选取55、77、130、164 m/s这4个试验段气流速度开展研究，其中各试验段气流速度(构型出口气流速度)对应的稳定段气流速度(构型入口气流速度)分别为4.11、5.70、9.15及11.00 m/s。

图8给出了不同试验段气流速度下构型中心线处液滴直径和相对湿度变化曲线，从图中可以看出：在准一维传热阶段，随着试验段气流速度的增加，稳定段气流速度增大，液滴相变时间减少，液滴尺寸的减小趋势减慢，表明蒸发效应减弱；在三维收缩阶段，随着试验段气流速度的增加，尽管液滴相变时间不断减少，但是在三维收缩效应的影响下，气流静温显著降低，水蒸气过饱和度显著增大，导致液滴尺寸的增大趋势加快，进而表明凝结效应增强。

图9给出了不同试验段气流速度下构型中心线处气流静温和液滴温度变化曲线，从图中可以看出，在准一维传热阶段，随着试验段气流速度的增加，稳定段气流速度增大，但是由于稳定段速度增大幅度较小(从约4 m/s增至11 m/s)，气流可压缩性较弱，因此各工况下的气流静温和液滴温度变化趋势无显著差异，气流静温接近气流总温，并且液滴温度稳定后，液滴湿球温度也基本一致。在三维收缩阶段，随着试验段气流速度的增加，气流可压缩性显著增强，试验段气流静温显著下降，具体的当试验段气流速度从55 m/s增至164 m/s时，试验段气流静温约从-7 ℃降至-18 ℃，因此在结冰风洞试验中，气流可压缩性导致的气流静温变化是不可忽略的；同时，在这一阶段随着试验段气流速度的增加，液滴温度的下降趋势加快，但是相对于颗粒，由于凝结效应的不断增强，抑制了液滴温度的下降趋势，导致出口处(试验段中心)液滴与颗粒温差不断增大。

图10给出了不同试验段气流速度下构型出口中心处液气温度差，从图中可以看出：随着试验段气流速度的增加，液气温差不断增大，并且在相变效应的影响下，其增大趋势强于颗粒，其中试验段气流速度为55 m/s时，液滴温差约为0.3 ℃，而颗粒的温差则趋近于0 ℃，当试验段气流速度增加至164 m/s时，液滴温差增至约3.3 ℃，比颗粒的温差高近2 ℃。由此可见，气流速度会显著影响试验段内液滴过冷状态，同时相变效应则会显著增强影响程度，其中小尺寸液滴(40 μm≤di≤100 μm)在高风速时(UTS≥164 m/s)将偏离过冷状态(液气温差超过3 ℃)。

2.4 液滴尺寸影响

最后，考察液滴尺寸的影响，由于采用球型液滴阻力模型，不考虑液滴变形和破碎效应，因此在We数的约束下[24]，仅选择尺寸较小的液滴，以较好地满足球型液滴假设。

图12给出了不同初始液滴尺寸下构型中心线处液滴温度变化曲线，从图中可以看出：在准一维传热阶段，随着初始液滴尺寸增加，液滴热容量增大，液滴温度下降趋势减慢，液滴温度稳定后，液滴湿球温度无显著变化，液滴与颗粒温差基本一致；在三维收缩阶段，随着初始液滴尺寸增加，液滴温度下降趋势减慢(与文献[10]的结论一致)，但是相对于颗粒，凝结效应不断减弱，促进了液滴温度的下降趋势，导致出口处(试验段中心)液滴与颗粒温差不断减小。

图13给出了不同初始液滴尺寸下构型出口中心处液气温度差，从图中可以看出：随着初始液滴尺寸增加，液气温差不断增大，但是在相变效应的影响下，其增大趋势弱于颗粒的温度差，其中液滴直径为40 μm时，液滴温度差约为0.6 ℃，而颗粒温度差约为0.1 ℃，当液滴直径增加至100 μm时，液滴温度差增至约0.85 ℃，比颗粒的温差低近0.05 ℃。由此可见，液滴尺寸会影响试验段液滴过冷状态，但相变效应却会减弱影响程差，其中在低风速(UTS≤77 m/s)试验工况下，小尺寸液滴(di≤100 μm)仍可以认为处于过冷状态(液气温差小于1 ℃)，这与文献[10]的评估结果一致。

3 结 论

本文发展了基于Euler法的气液两相传质传热耦合流动计算方法，模拟了结冰风洞中气液两相传质传热耦合流动过程，研究了结冰风洞中液滴相变效应。

1) 结冰风洞中液滴经历了先蒸发后凝结两个阶段。在准一维传热阶段，液滴为蒸发状态，蒸发效应加快了液滴温度下降趋势，并且促使液滴温度趋于湿球温度；在三维收缩阶段，液滴则从蒸发状态转换为凝结状态，凝结效应减慢了液滴温度下降趋势，导致构型出口处(试验段中心)液气温差增大，影响液滴过冷状态。

2) 初始相对湿度增加，会减弱蒸发效应而增强凝结效应，进而影响试验段内液滴过冷状态。但其影响程度较弱，本文计算工况下最大液气温差(相对湿度为100%)小于1 ℃，液滴仍认为处于过冷状态。

3) 气流速度会显著影响试验段内液滴过冷状态，随着气流速度的增加，蒸发效应减弱而凝结效应增强，进而显著增强了影响程度，其中小尺寸液滴(40 μm≤di≤100 μm)在高风速时(UTS≥164 m/s)将偏离过冷状态(液气温差超过3 ℃)。

4) 初始液滴尺寸增加，会导致蒸发效应和凝结效应均减弱，进而减弱了液滴尺寸对试验段内液滴过冷状态的影响程度，对于小尺寸液滴(di≤100 μm)在低风速(UTS≤77 m/s)试验工况下，仍可以认为处于过冷状态(液气温差小于1 ℃)。

需要指出的是，实际结冰风洞中，气液两相流动过程十分复杂，涉及动量传递、质量传递和能量传递等多物理过程的强耦合作用，而本文发展的基于Euler法的气液两相传质传热耦合流动计算方法是基于一定假设下的简化分析方法，因此开展计算方法验证试验则是下一步急需开展的工作。

此外，针对液滴直径超过100 μm的大尺寸液滴，由于其存在变形和破碎等动力学效应，超出了本文采用的球形液滴数学模型适用范围，因此建立大尺寸液滴传质传热模型，则是今后需要开展的研究。

参 考 文 献

[1] 杜雁霞， 李明， 桂业伟， 等. 飞机结冰热力学行为研究综述[J]. 航空学报， 2017， 38(2)： 520706.

DU Y X，LI M, GUI Y W, et al. Review of thermodynamic behaviors in aircraft icing process[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 520706 (in Chinese).

[2] 桂业伟， 周志宏， 李颖晖， 等. 关于飞机结冰的多重安全边界问题[J]. 航空学报， 2017， 38(2): 520723.

GUI Y W, ZHOU Z H, LI Y H, et al. Multiple safety boundaries protection on aircraft icing[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 520723 (in Chinese).

[3] 易贤， 王斌， 李伟斌， 等. 飞机结冰冰形测量方法研究进展[J]. 航空学报， 2017, 38(2): 520700.

YI X, WANG B, LI W B, et al. Research progress on ice shape measurement approaches for aircraft icing[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 520700 (in Chinese).

[4] RAGNI A, ESPOSITO B, MARRAZZO M, et al. Calibration of the CIRA IWT in the high speed configuration[C]//43th AIAA Aerospace Science Meeting and Exhibit. Reston, VA: AIAA, 2005.

[5] STEEN L E, IDE R F, van ZANTE J F, et al. NASA Glenn icing research tunnel: 2014 and 2015 cloud calibration procedures and results: NASA/TM-2015-218758 [R]. Washington, D. C. : NASA, 2015.

[6] WILLBANKS C E, SCHULZT R J. Analytical study of icing simulation for turbine engines in altitude test cells[J]. Journal of Aircraft, 1975, 12(12): 960-967.

[7] MILLER D R, ADDY H E, IDE R F. A study of large droplet ice accretions in the nasa glenn irt at near-freezing conditions: NASA TM-1996-107142-REV1 [R]. Washington, D. C. : NASA, 1996.

[8] BELLUCCI M, ESPOSITO B M, MARRAZZO M, et al. Calibration of the CIRA IWT in the low speed configuration[C]∥45th AIAA Aerospace Science Meeting and Exhibit. Reston, VA： AIAA，2007.

[9] 易贤, 马洪林, 王开春, 等. 结冰风洞液滴运动及传质传热特性分析[J]. 四川大学学报, 2012, 44(2): 132-135.

YI X, MA H L, WANG K C, et al. Analysis of water droplet movement and heat / mass transfer in an icing wind tunnel[J]. Journal of Sichuan University, 2012, 44(2): 132-135 (in Chinese).

[10] 郭向东，王梓旭，李明，等．结冰风洞中液滴过冷特性数值研究[J]．航空学报, 2018, 39(1): 121254.

GUO X D, WANG Z X, LI M, et al. Numerical study of droplet supercooling in an icing wind tunnel [J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(1): 121254 (in Chinese).

[11] HILL P G. Condensation of water vapor during supersonic expansion in nozzles[J]. Journal of Fluid Mechanics, 1966, 25(3): 593-620.

[12] CHENG W, LUO X S, van DONGEN M E H. On condensation-induced waves[J]. Journal of Fluid Mechanics, 2010, 651(1): 145-164.

[13] 郭向东, 黄生洪, 吴颖川, 等. 氢燃料燃烧加热风洞中水蒸气相变效应的数值研究[J]. 推进技术, 2017, 38(4):932-941.

GUO X D, HUANG S H, WU Y C, et al. Numerical evaluation of water vapor phase transition effects in a hydrogen-fueled combustion-heated wind tunnel[J]. Journal of Propulsion Technology, 2017, 38(4): 932-941 (in Chinese).

[14] CAO Y, LUO X S. Numerical study on Homogeneous condensation in a vortex[J]. Procedia Engineering, 2015, 126: 607-611.

[15] LUO X S, WANG G, OLIVIER H. Parametric investigation of particle acceleration in high enthalpy conical nozzle flows for coating applications[J]. Shock Waves, 2008, 17(5): 351-362.

[16] KERSEY J, LOTH E, LANKFORD D. Effects of Evaporating Droplets on Shock Waves[J]. AIAA Journal, 2010, 48(9): 1975-1986.

[17] RUSSO E, KUERTEN J G M, van der Geld C W M, et al. Water droplet condensation and evaporation in turbulent channel flow[J]. Journal of Fluid Mechanics, 2014, 749: 666-700.

[18] SAITO T. Numerical analysis of dusty-gas flows[J]. Journal of Computational Physics, 2002, 176(1): 129-144.

[19] FUJITA A, KUROSE R, KOMORI S. Experimental study on effect of relative humidity on heat transfer of an evaporating water droplet in air flow[J]. International Journal of Heat and Mass Transfer, 2010, 36(3): 244-247.

[20] MILLER R S, HARSTAD K, BELLAN J. Evaluation of equilibrium and non-equilibrium evaporation models for many-droplet gas-liquid flow simulations[J]. International Journal of Multiphase Flow, 1998, 24(6): 1025-1055.

[21] HINDMARSH J P, RUSSEL A B, CHEN X D. Experimental and numerical analysis of the temperature transition of a suspended freezing water droplet[J]. International Journal of Heat and Mass Transfer, 2003, 46(7): 1199-1213.

[22] COBER S, BERNSTEIN B, JECK R, et al. Data and analysis for the development of an engineering standard for supercooled large drop conditions: DOT/FAA/AR-09/10[R]. Washington, D. C.: FAA, 2009.

[23] KUROSE R, FUJITA A, KOMORI S. Effect of relative humidity on heat transfer across the surface of an evaporating water droplet in air flow[J]. Journal of Fluid Mechanics, 2009, 624: 57-67.

[24] LUXFORD G. Experimental and modelling investigation of the deformation, drag and break-up of drizzle droplets subjected to strong aerodynamic forces in relation to SLD aircraft icing[D]. Bedfordshire: Cranfield University, 2005: 5-7 .