结冰风洞中过冷大水滴云雾演化特性数值研究

郭向东，柳庆林，刘森云，王梓旭，李明

1. 中国空气动力研究与发展中心 空气动力学国家重点实验室，绵阳 621000

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

飞机结冰广泛存在于飞行实践中，并严重威胁飞行安全[1-2]。鉴于结冰的严重危害，适航规章要求制造商标明飞机在结冰气象条件下的适航符合性，确保飞行安全。目前在适航审定中普遍采用的结冰气象条件由适航条例25部附录C给出，该条件下最大液滴直径约为100 μm，但是近年来发生的多起飞行事故表明，大气环境中存在直径超过100 μm的过冷大水滴(Supercooled Large Droplet, SLD)，因此适航条例进一步扩展了结冰气象条件，形成了附录O条款，给出了SLD结冰气象条件[3-4]。该结冰条件由冻细雨和冻雨两种条件组成，其中冻细雨条件下最大液滴直径在100～500 μm，冻雨条件下最大液滴直径超过500 μm。

结冰风洞试验作为重要的结冰适航审定手段，已经广泛应用于传统结冰条件(附录C)下的飞机结冰适航审定[5-6]。具体而言，结冰风洞利用喷雾系统产生水滴云雾，随后水滴云雾由气流携带经过水平收缩段进入试验段，在这一过程中，云雾在水平收缩段内经历了复杂的力学演化过程，最终在试验段内达到了分布均匀且力学平衡(液滴尺寸和浓度分布均匀、速度和温度趋于稳定)的试验状态。但是针对SLD结冰条件，SLD云雾中不仅存在直径小于100 μm的小尺寸液滴，而且包含大尺寸的SLD液滴，而这些SLD液滴存在的重力沉降、变形破碎等力学效应均会影响其在收缩段内的运动、传热和传质过程[7]，进而可能导致SLD经过云雾演化后，在试验段内无法完全达到分布均匀且力学平衡的试验状态。因此，明晰结冰风洞中SLD云雾演化特性，全面考察SLD在试验段内的力学平衡状态，对结冰风洞准确模拟SLD结冰条件具有重要意义。

针对结冰风洞中的云雾演化过程，国外相关研究可以追溯到20世纪70年代，美国阿诺德工程发展中心Willbanks和Schulzt[8]发展了一维两相流动计算方法(AEDC 1DMP)，针对直径小于40 μm的液滴，开展了发动机高空试验台内的液滴热平衡特性研究。加拿大阿尔伯塔大学Gates等[9]基于拉格朗日方法发展了两相流计算模型，模拟了结冰风洞中轴对称收缩段构型内的云雾演化过程，详细研究了直径小于100 μm的球形液滴在试验段气流速度小于50 m/s条件下的力学平衡特性。美国NASA格伦研究中心Miller等[10]采用AEDC 1DMP计算方法，研究了NASA Glenn IRT结冰风洞试验段内的SLD热平衡特性，评估了直径为40、99、160 μm 3种尺寸的液滴，在试验段气流速度195 mph (约87 m/s)条件下的热平衡状态。在接下来的20多年里，IRT结冰风洞开展了大量SLD云雾试验模拟研究[11-13]，逐渐意识到SLD液滴存在的重力沉降、变形破碎等力学特性，会显著影响其在试验段内的力学平衡特性，制约着结冰风洞对SLD结冰条件的准确模拟。最近，加拿大国家研究院Orchard等[14]同样采用AEDC 1DMP计算程序，开展了SLD结冰条件下的结冰风洞水平收缩段构型设计，详细探讨了直径小于500 μm的SLD的力学平衡特性。近年来，随着国内大型结冰风洞——3 m×2 m结冰风洞的建成，一些学者针对结冰风洞中的云雾演化特性开展了初步研究。中国空气动力研究与发展中心易贤等[15]建立了基于拉格朗日法的液滴运动及传质传热数值模拟方法，针对3 m×2 m结冰风洞主试验段水平收缩构型，初步研究了直径小于200 μm液滴在不可压缩气流中的运动传热过程。南京航空航天大学朱春玲等[16]也发展了基于拉格朗日法的气液两相流计算模型，针对直径小于50 μm的液滴，开展了参数影响研究。尽管拉格朗日法可以准确追踪单个液滴的运动、传热和传质过程，但是针对3 m×2 m结冰风洞这种大尺寸计算构型，需要计算大量的液滴运动和传热传质过程，对计算能力提出了严峻的挑战。而欧拉法将云雾液滴视为连续相，这样可以极大地减少计算量，适用于大尺寸计算构型。基于工程应用的需要，本文作者团队[17-18]发展了基于欧拉法的气液两相传质传热耦合流动计算方法，详细研究了3 m×2 m结冰风洞主试验段水平收缩构型内直径小于100 μm的小尺寸液滴的传热传质过程。由此可见，目前国内外研究均基于球形液滴假设、未考虑液滴变形破碎等力学效应，同时研究范围均未覆盖冻雨SLD结冰条件，缺乏对结冰风洞中SLD云雾演化特性的全面认识。

本文在前期计算方法研究基础上，发展了基于欧拉法的SLD液滴运动、传热和传质耦合计算方法，针对3 m×2 m结冰风洞主试验段水平收缩构型，分析了SLD云雾沉降收缩特性、动量平衡特性和热平衡特性，探索了液滴变形破碎的影响，评估了构型出口处SLD液滴动量平衡和热平衡状态，为3 m×2 m结冰风洞准确模拟SLD结冰条件奠定了基础。

1 计算方法

本文基于前期发展的气液两相传质传热耦合流动计算方法，进一步引入了液滴变形和破碎模型，同时考虑了重力作用，发展了基于欧拉法的SLD液滴运动、传热和传质耦合计算方法。为获得简化物理模型，对气液两相进行以下假设[17-19]：① 气相为理想气体，遵循理想气体法则；② 液滴变形过程为准静态过程，不考虑液滴振动过程，且液滴形状从球形变形为扁椭球体(Oblate Spheroid)；③ 液滴体积相对于气相体积可以忽略，不考虑液滴间的碰并作用；④ 液滴内温度均匀分布，热传导仅发生在气液两相间，忽略辐射传热过程，不考虑壁面黏性效应；⑤ 不考虑液滴冻结和水蒸气凝结成核过程。

1.1 物理模型

1.1.1 非球形液滴几何形状参数

许多参数被用于描述非球形液滴的形状特征[20-21]，本文采用以下3个形状参数(见图1)：

图1 非球形液滴几何形状参数描述

1) 液滴体积等效直径d(Droplet Volume Equivalent Diameter)，定义为与非球形液滴体积相同的球形液滴直径：

(1)

式中：V为液滴体积；r为液滴体积等效半径。

2) 液滴纵横比E(Droplet Aspect Ratio)，定义为液滴横向直径b(沿旋转轴)与液滴纵向直径a(垂直于旋转轴，也称为液滴赤道直径)之比：

(2)

根据液滴纵横比，类球形液滴(Spheroids)可以分为扁椭球形液滴(Oblate Spheroids,E<1)和长椭球形液滴(Prolate Spheroids,E>1)。显而易见，当E=1时，液滴为球形，当E→0时，液滴趋于圆盘形，当E→∞时，液滴则趋于针形。

3) 液滴球形度Φ(Droplet Sphericity)，定义为体积等效球形液滴表面积As与液滴实际表面积A之比：

(3)

针对类球形液滴，液滴球形度Φ为液滴纵横比E的函数：

(4)

1.1.2 非球形液滴阻力模型

液滴阻力FD表示为

(5)

式中：CD为液滴阻力系数；μg为气相力学黏性系数；Vg和Vd分别为气相和液滴速度矢量；Red为液滴相对雷诺数，表示为

(6)

式中：ρg为气相密度。针对非球形液滴，液滴阻力系数CD受Red和液滴相对韦伯数Wed的共同影响，Wed表示为

(7)

式中：σST为液滴表面张力。

Luxford等[22-23]基于自由下落液滴阻力实验数据[24]，假设液滴以准静态方式变形为扁椭球体，通过圆球阻力系数CD,Sphere和圆盘阻力系数CD,Disk插值得到了液滴阻力系数CD,equ：

(8)

式中：Red,equ和CD,equ均基于扁椭球体赤道直径a，进一步联立式(1)、式(2)和式(6)，则a和Red,equ表示为

(9)

由于在工程实践中，液滴阻力系数一般基于液滴体积等效直径，因此CD通常表示为

(10)

Luxford 阻力模型利用自由下落液滴阻力实验数据获得插值因子k的函数关系，其中当k→0时，液滴阻力系数趋近于圆盘的阻力系数，而当k→1时，液滴阻力系数趋于圆球的阻力系数。最后，根据文献[20]，Luxford 阻力模型适用于液滴直径在100～1 000 μm范围内的液滴，即适用于本文研究范围。

1.1.3 液滴破碎模型

随着Wed的增加，液滴的变形程度加剧，当Wed超过临界相对韦伯数Wec时，液滴可能发生破碎。文献[25-26]提出了基于欧拉法的液滴数密度方程，并通过方程源项实现了液滴破碎过程的模拟，方程表示为

(11)

Sbu=

(12)

式中：N为液滴数密度；Sbu为液滴数密度源项；m取0.2；dstab为液滴破碎后最大稳定液滴直径，表示为

(13)

其中：水滴的Wec一般取为12；Tbu为无量纲液滴破碎时间，Pilch和Erdman[27]根据实验结果给出了Tbu的经验关系式为

(14)

1.1.4 非球形液滴传热传质模型

(15)

式中：kg为气相热传导系数；Tg和Td分别为气相温度和液滴温度。针对非球形液滴，法国航空航天研究院Villedieu等[29-30]采用雷诺类比法，基于液滴球形度Φ给出了非球形液滴努赛尔数Nu(Φ)的关系式：

(16)

(17)

(18)

式中：DAB为二元质量扩散系数；Rv为水蒸气气体常数；pv,s为液滴表面饱和水气压；S为相对湿度。进一步类比传热过程，Villedieu 模型同时给出了非球形液滴舍伍德数Sh(Φ)关系式，表示为

(19)

(20)

应该指出的是，基于液滴准静态方式变形假设，Villedieu模型中球形度可以通过Luxford模型和式(4)计算得到。

1.2 控制方程组

基于以上假设和物理模型，简化后的控制方程组包括：

1) 气相方程组

(21)

(22)

(23)

(24)

2) 液相方程组

(25)

(26)

(27)

(28)

式中：p为气相压力；气相能量Eg和液相能量Ed分别表示为

(29)

其中：Cv为气相定容比热；Cl为液相比热；σ和α分别为液滴有效密度和水蒸气质量分数，其中σ对应液态水含量；ρl为液滴物理密度；g为重力加速度矢量；qg和ql分别为单位质量传质引起的气相和液相能量变化率，表示为

(30)

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

1.3 数值方法

本文采用文献[17-18]发展的基于有限体积法的数值方法，其中非定常项采用二阶Euler格式离散，空间输运项采用二阶迎风型的Roe格式离散，源项采用隐式格式。

气液两相边界条件采用不同格式处理：气相入口条件采用压力入口，出口条件采用压力出口，壁面采用无黏滑移壁面条件(假设④)；液相入口条件采用速度入口，出口采用出流条件。液相壁面边界条件需要特殊处理[19-20](如图2所示)：当临近壁面处的网格单元内液滴速度的壁面法向分量指向壁面时，壁面处采用出流边界条件；当该速度分量指向流场内部(计算域内)时，为消除这种非物理过程对液相流场的影响，一方面在壁面处将液滴有效密度σ设置为10-9，以保证计算稳定收敛，另一方面采用镜面边界条件处理壁面液相速度。

图2 液相壁面边界条件示意图

本文采用的数值方法在文献[17-18]中得到了验证，此处不再赘述。

2 计算构型和计算工况

3 m×2 m结冰风洞主试验段水平收缩构型如图3所示，图中X方向沿构型水平中心线指向构型出口，Y方向垂直于构型竖直对称面(Symmetry，绿色平面)，Z方向为重力作用的反方向。该构型沿X方向依次由稳定段、收缩段和试验段3部分组成(黄色曲面)，其中构型入口截面(Inlet，粉红色平面)位于稳定段喷雾耙处，出口截面(Outlet，蓝色平面)为试验段中心截面，构型收缩比为14.67。

图3 3 m×2 m结冰风洞主试验段水平收缩构型

为覆盖冻细雨和冻雨两种SLD结冰条件，同时考虑计算模型适用范围，本文仅针对最大液滴直径为1 000 μm的SLD云雾开展研究，计算工况矩阵如表1所示，表中参数根据结冰风洞实际试验条件以及喷嘴测试实际结果确定。应该指出的是文献[18]已经对直径小于100 μm的小尺寸液滴开展了详细研究，因此本文仅考察液滴直径大于100 μm的情况。

表1 计算工况矩阵

3 结果与讨论

3.1 SLD云雾沉降收缩特性

首先，图4给出了构型竖直对称面处液态水含量(LWC)分布云图，其中图4(a) ～ 图4(d)分别为1 000 μm工况下40、80、120、160 m/s对应的计算结果，黑色实线表示液滴运动轨迹线，共选取了7根轨迹线，迹线起始位置等间距分布于构型对称面入口-3 m≤Zi≤3 m范围内。从图中可以看出，云雾首先在重力的影响下整体向构型下部沉降，然后在气动力的影响下表现出显著的收缩汇聚特征，同时液滴轨迹线相对于构型中心线表现出显著的非对称性。文献[17-18]将该构型内液滴运动过程先后分为准一维阶段(Quasi-1D Stage)和三维收缩阶段(Contraction-3D Stage)，本文参考该分类方式：在准一维阶段(0 m≤X≤8 m)， 气流以准一维方式低速稳定流动，此时液滴主要在重力作用下发生沉降；在三维收缩阶段(8 m≤X≤18.33 m)，气流加速收缩汇聚流动，此时液滴主要在气动力的作用下表现出收缩汇聚的运动特征。

图4 构型对称面处液态水含量分布云图(di=1 000 μm)

接着，为揭示液滴尺寸和试验段气流速度Vg,x对SLD液滴沉降收缩特性的影响，图5给出了起始位置位于构型入口中心处(Zi=0 m)的液滴轨迹线，黑色点虚线表示构型中心线处气流速度，背景阴影区表示构型对称面，黄色和青色区域分别对应准一维阶段和三维收缩阶段。从图中可以看出：增大液滴尺寸，会促使液滴轨迹线向构型下部偏折，显著加强SLD液滴的沉降程度；增大试验段气流速度，一方面会减弱准一维阶段内液滴轨迹线的偏折程度，抑制液滴沉降过程，另一方面会增强三维收缩阶段内液滴轨迹线向构型中心线处的收缩汇聚程度，促进液滴收缩汇聚。

图5 构型入口中心处液滴轨迹线

最后，图6给出了构型出口处液态水含量分布。从图中可以看出： SLD液滴在构型出口处会出现沉降现象，液滴尺寸越大，则沉降越显著，尤其在40 m/s工况下，直径超过500 μm的液滴均无法到达构型出口，但增大气流速度会显著减弱液滴沉降程度、增强液滴收缩汇聚趋势，促使液滴趋近于构型中心区域，尤其在160 m/s工况下，出口处液滴分布区域相对于构型中心线表现出显著的中心对称性。

图6 构型出口处液态水含量分布

因此，构型出口SLD云雾(最大液滴直径小于1 000 μm)在重力沉降效应影响下会出现显著的粒径浓度分层现象，云雾上部浓度较低且多为直径小于100 μm的小尺寸液滴，而云雾下部浓度较高且多为大尺寸SLD液滴，但增大气流速度会显著减弱分层程度，尤其在160 m/s工况下，SLD云雾会均匀分布在构型出口中心区域内(-0.75 m

3.2 SLD云雾动量平衡特性

首先，为明晰构型内SLD液滴运动特征，图7给出了构型竖直对称面液滴轨迹线处液滴速度变化曲线，其中液滴轨迹线选取了160 m/s、1 000 μm工况下Zi=2， 1，0，-1，-2 m 5条轨迹线，图7(a)和图7(b)分别对应X方向和Z方向液滴速度分量(Vd,x和Vd,z)变化曲线。从图中可以看出，在X方向上，各轨迹线处液滴速度分量均在准一维阶段不断减小，而在三维收缩阶段则显著增大，可见SLD液滴在X方向上表现出先减速再加速的运动特征。在Z方向上，各迹线处液滴速度分量在准一维阶段均逐渐减小，但在三维收缩阶段内，不同轨迹线间的液滴速度分量却出现了显著差异，液滴迹线越靠近洞壁，对应的液滴速度分量绝对值则越大，同时液滴迹线的非对称性导致了液滴速度分量的非对称分布。由此可见，构型内SLD液滴在Z方向上表现出先加速沉降再非对称收缩汇聚的运动特征。

图7 构型竖直对称面液滴轨迹线处液滴速度变化曲线

然后，为考察构型内SLD液滴动量平衡特性，根据液相动量方程(式(26))，推导出单液滴运动方程：

(31)

可见，当液滴达到动量平衡后(液滴速度时间变化率趋于零)，在X和Y方向上，液滴与气流速度分量趋于一致(Vg,x=Vd,x和Vg,y=Vd,y)，而在Z方向上，在重力的影响下，液滴速度分量小于气流速度分量(Vd,z

接着，为考察液滴尺寸和试验段气流速度对SLD液滴动量平衡特性的影响，图8给出了构型入口中心液滴轨迹线处X方向气液速度分量差(Vg,x-Vd,x)和构型出口竖直对称线处Z方向液滴速度分量与液滴最终自由下落速度差(Vd,z-Vd,ter)变化曲线。从图中可以看出，在X方向上，SLD液滴先经过减速运动后达到了动量平衡状态，然后在加速过程中则显著偏离了动量平衡状态，其中在X=13.5 m处(试验段入口处)出现了速度偏差峰值。增大液滴尺寸和试验段气流速度，不仅会增大准一维阶段内X方向液滴速度平衡距离(液滴X方向速度分量平衡点的X坐标值)，而且会显著增强三维收缩阶段内气液速度偏差的增大趋势，进而显著增大气液速度偏差峰值，最终导致构型出口气液速度偏差不断增大，SLD液滴则不断偏离动量平衡状态。应该指出的是在160 m/s工况下，当液滴直径超过500 μm后，气液速度偏差的增大趋势则趋于一致，出口处气液速度偏差并无显著变化。在Z方向上，构型出口处SLD液滴表现出非对称收缩汇聚运动特征，因此液滴速度分量与液滴最终自由下落速度偏差从出口上部至下部不断增大，进而在最下部达到速度偏差峰值，导致SLD液滴在出口最下部动量平衡状态偏离的最显著；增大液滴尺寸和试验段气流速度会增强出口液滴的收缩汇聚程度，并且增大液滴尺寸还会减小液滴最终自由下落速度，最终导致下边界处的速度偏差峰值不断增大，SLD液滴则不断偏离动量平衡状态。

为探索构型内液滴变形破碎特性，图9给出了构型中心液滴轨迹线处液滴相对韦伯数和液滴形状参数变化曲线，图中仅选取160 m/s工况进行讨论。从图中可以看出： SLD液滴在准一维阶段未发生显著变形，但在三维收缩阶段内，液滴直径超过250 μm后，Wed陡然增大，液滴纵横比和球形度显著减小，液滴则发生了显著形变，其中在X=13.5 m处出现了变形峰值区；当液滴直径超过250 μm后，增大液滴尺寸会显著增强三维收缩阶段内Wed的增大趋势以及液滴纵横比和球形度的减小趋势，进而显著增强了液滴的形变程度，尤其当液滴直径超过750 μm后，Wed超过Wec，液滴尺寸陡然减小，液滴发生显著破碎。

为揭示液滴变形破碎特性对SLD液滴动量平衡特性的影响，根据式(31)可见，液滴加速度正比与液滴阻力系数、反比于液滴直径(dVd/dt～CD/d)，而液滴变形破碎会直接影响液滴阻力系数和直径，因此图10给出了构型入口中心液滴轨迹线处液滴相对雷诺数和液滴阻力系数变化曲线(分别对应图10(a)和图10(c))，图10(b)为Luxford阻力模型给出的CD随Red变化曲线，空心方框和空心圆分别对应构型内Red最小和最大值点，图中仅选取160 m/s工况进行讨论。从图中可以看出：Red在准一维阶段内不断减小，其中最小值点(空心方框)位置对应X方向液滴速度平衡位置(见图8(a))，接着在三维收缩阶段内则陡然增大，并且在X=13.5 m处达到最大值(空心圆)，对应X方向液滴速度分量偏差和液滴变形峰值区(见图8(a)和图9)；增大液滴尺寸会增大构型内的Red。根据Luxford阻力模型，CD随Red的增大表现出先减小后增大的趋势，其中在递减区，CD与球形的阻力系数(点虚线)接近，但在递增区，在液滴变形效应的影响下，CD会显著大于球形的阻力系数；增大液滴尺寸会减小CD的最小值(褐色点线)，进而增大最小值点的Red。因此，在准一维阶段，液滴不发生显著形变，Red处于递减区，CD与球形液滴的阻力系数一致，并在Red最小值处达到峰值(空心方框)。而在三维收缩阶段，当液滴直径超过250 μm后，Red增大至递增区，CD在Red最大值处(X=13.5 m)达到峰值(空心圆)，此时液滴变形效应促使峰值区CD显著大于球形液滴的阻力系数，进而增大了液滴加速度，促使SLD液滴趋近动量平衡状态。增大液滴尺寸，尽管会直接减小液滴加速度，但同时会增强三维收缩阶段内SLD液滴变形程度，增大CD的增大幅度，进而减弱了液滴加速度的减小趋势，抑制了出口处SLD液滴偏离动量平衡状态，尤其在160 m/s工况下，当液滴直径超过500 μm后，X方向气液速度偏差的增大趋势则趋于一致，出口处气液速度偏差并无显著变化(见图8(a))。此外，液滴破碎效应会直接减小液滴直径，增大液滴加速度，同样会抑制SLD液滴动量平衡状态的偏离。

图9 构型入口中心液滴轨迹线处液滴相对韦伯数和液滴形状参数变化曲线

最后，为评估构型出口SLD液滴动量平衡状态，图11给出了构型出口处X方向气液速度分量差(Vg,x-Vd,x)和Z方向液滴速度分量与液滴最终自由下落速度之差的最大绝对值(|Vd,z-Vd,ter|max)，其中图11(b)表格给出了构型出口处液滴最终自由下落速度。从图中可以看出，增大液滴尺寸和试验段气流速度，会同时增大X方向和Z方向上速度偏差，进而增强了SLD液滴动量平衡状态的偏离程度，尤其在160 m/s、1 000 μm工况下，速度偏差将超过18 m/s。这一结论与文献[33] 的试验结果相一致。

图11 构型出口处X方向气液速度分量差(Vg,x- Vd,x)和Z方向液滴速度分量与液滴最终自由下落速度之差最大绝对值(|Vd,z-Vd,ter|max)

因此，构型出口SLD云雾(最大液滴直径小于1 000 μm)会出现显著的动量分层现象，云雾中直径小于100 μm的小尺寸液滴速度快则趋于动量平衡态(与平衡态最大速度偏差约为2 m/s)， 但直径超过500 μm的大尺寸SLD液滴速度慢则显著偏离平衡态，并且增大试验段气流速度会显著增强分层程度，尤其在160 m/s工况下，SLD云雾与其平衡态间的最大速度偏差将会超过18 m/s。

3.3 SLD云雾热平衡特性

首先，根据液相能量方程(式(27))，推导得单液滴能量方程：

(32)

可见，当液滴达到热平衡时(液滴温度时间变化率趋于零)，液滴对流传热率与液滴相变传热率相等，液滴温度趋于平衡温度，则将该平衡温度称为液滴湿球温度Twb。图12给出了构型中心线处气流与液滴温度参数变化曲线，从图中可以看出：在准一维阶段内，气流为欠饱和状态(S<100%)，则Twb100%)，进而导致Twb>Tg；最终，在构型出口处，增大试验段气流速度会显著降低Twb，同时会显著增大相对湿度，进而增大了液滴湿球温度与气流静温差。此处应该指出的是，图12中计算的液滴湿球温度并未考虑液滴直径的影响，这主要是由于在本文计算条件下，Nu与Sh之间的偏差小于5%(Pr与Sc间的偏差约为10%)，因此计算中认为Nu与Sh近似相等，可见液滴湿球温度仅为气流相对湿度和气流静温的函数。

图12 构型中心线处气流和液滴温度参数变化曲线

接着，为揭示液滴尺寸和试验段气流速度对构型内SLD液滴热平衡特性的影响，图13给出了构型入口中心液滴轨迹线处液滴温度与液滴湿球温度差变化曲线，其中图13(a)和图13(b)分别为160 m/s工况和250 μm工况下的计算结果。从图中可以看出：在准一维阶段，温度差不断减小，液滴逐渐趋近热平衡状态，进入三维收缩阶段后，液滴湿球温度陡然下降，温度差则不断增大，液滴逐渐偏离热平衡状态，其中在X=13.5 m处出现了温度差峰值；增大液滴尺寸和试验段气流速度，一方面会减弱准一维阶段内温差的减小趋势，导致准一维阶段出口处(X=8 m) 温差不断增大，另一方面会增强三维收缩阶段内温差的增大趋势，最终导致构型出口处温度差不断增大，液滴则不断偏离热平衡状态。

图13 构型入口中心液滴轨迹线处液滴温度与液滴湿球温度差变化曲线

考虑到SLD液滴温度变化主要受对流传热和相变传热的影响，为分析这两个过程对构型内SLD液滴热平衡特性的影响，将对流传热温度变化率ξconv和相变传热温度变化率ξphase定义为

(33)

图14给出了构型入口中心液滴轨迹线处ξconv和ξphase变化曲线，其中仅选取160 m/s工况进行讨论。从图中可以看出：针对直径小于100 μm的小尺寸液滴，在准一维阶段，液滴处于蒸发状态，则ξconv和ξphase均为负且不断增大，液滴温度陡然下降，当液滴温度趋于液滴湿球温度时，ξconv+ξphase=0，液滴达到稳定蒸发的热平衡状态；紧接着在三维收缩阶段内，气流静温陡然下降，ξconv则陡然降低，促使液滴温度不断降低，同时相对湿度陡然上升，气流从欠饱和状态转换成过饱和状态，液滴则从蒸发状态转换到凝结状态，ξphase不断增大，抑制了液滴温度的降低，这与文献[17]的结论一致。针对SLD液滴，在准一维阶段，增大液滴尺寸会显著减弱对流传热和相变传热强度，抑制液滴趋于热平衡态，导致直径大于500 μm的SLD液滴均无法达到热平衡状态；在三维收缩阶段，增大液滴尺寸一方面会减弱对流传热强度(减弱ξconv的降低趋势)，抑制液滴趋于热平衡态，另一方面会促使液滴进入蒸发状态(直径超过500 μm的液滴均处于蒸发状态)，并不断增强蒸发传热强度(增强ξphase的降低趋势)，促进液滴趋于热平衡态。

图14 构型入口中心液滴轨迹线处ξconv和ξphase变化曲线

为分析液滴变形破碎对SLD液滴热平衡特性的影响，根据式(33)得

(34)

可见，对流传热温度变化率和相变传热温度变化率分别正比于Nu和Sh，而均反比于液滴直径和液滴球形度，同时液滴变形破碎会直接影响Nu、Sh、直径和球形度。图15给出了构型入口中心液滴轨迹线处Nu/Φ和Sh/Φ的变化曲线，空心圆对应构型内Red最大值点(见图10(a))，点虚线为圆球Nu和Sh变化曲线，图中仅选取160 m/s 工况进行讨论。从图中可以看出，Nu/Φ和Sh/Φ随Red增大而不断增大，相较于圆球，液滴变形会增大Nu/Φ和Sh/Φ，并且增大液滴尺寸会减慢Nu/Φ和Sh/Φ的增大趋势。因此，在准一维阶段，液滴不发生显著形变，Nu/Φ和Sh/Φ与球形液滴的一致，但是在三维收缩阶段，Nu/Φ和Sh/Φ在X=13.5 m处达到峰值(空心圆所示)，进而当液滴直径超过500 μm后，液滴变形效应会促使峰值区Nu/Φ和Sh/Φ显著大于球形液滴的，进而增强了对流传热和蒸发传热强度(减小了ξconv和ξphase)，加快了液滴温度的下降，促进液滴趋于热平衡状态。增大液滴尺寸，尽管会直接减弱对流传热和蒸发传热强度(降低了ξconv和ξphase减小趋势)，减小液滴温度下降率，但同时会增强三维收缩阶段内SLD液滴变形程度，增大Nu/Φ和Sh/Φ的增大幅度，降低液滴温度下降率的减小趋势，尤其在160 m/s工况下，当液滴直径超过500 μm后，ξconv和ξphase随液滴尺寸的增大并无显著变化(见图14)，此时液滴温度与液滴湿球温度差的增大幅度则趋于一致(见图13(a))。此外，液滴破碎效应会直接减小液滴直径，进而增强对流传热和蒸发传热强度，同样会抑制SLD液滴热平衡状态的偏离。

图15 构型入口中心液滴轨迹线处Nu/Φ和Sh/Φ变化曲线

图16给出了构型出口处液滴温度与液滴湿球温度差，其中表格给出了出口处气流静温、相对湿度和液滴湿球温度参数值。从图中可以看出，增大液滴尺寸和试验段气流速度会增大液滴温度与液滴湿球温度差，进而显著增强SLD液滴热平衡状态的偏离程度，尤其在160 m/s、1 000 μm工况下，温度偏差将超过20 ℃。

图16 构型出口处液滴温度与液滴湿球温度差

因此，构型出口SLD云雾(最大液滴直径小于1 000 μm)会出现显著的热分层现象，云雾中直径小于100 μm的小尺寸液滴温度低且不断凝结，趋近热平衡态(与平衡态最大温度偏差约为3 ℃)， 但直径超过500 μm的大尺寸SLD液滴温度高且不断蒸发，显著偏离平衡态，并且增大试验段气流速度会显著增强分层程度，尤其在160 m/s工况下，SLD云雾与其平衡态间的最大温差将超过20 ℃。

4 结 论

1) 直径超过250 μm的SLD液滴在构型三维收缩阶段内会发生显著形变，液滴尺寸越大则变形程度越强，尤其在160 m/s工况下，当液滴直径超过750 μm后，液滴甚至会发生破碎。

2) 液滴变形破碎效应会增大构型内SLD液滴阻力系数、增强对流传热和蒸发传热强度，进而增大了液滴加速度和液滴温度下降率，促使SLD液滴趋近动量平衡和热平衡状态。

3) SLD云雾(最大液滴直径小于1 000 μm)在构型出口处会出现显著的粒径浓度分层、动量分层和热分层现象，其中云雾上部浓度较低且多为直径小于100 μm的小尺寸液滴，而云雾下部浓度较高且集中较多大尺寸SLD液滴。直径小于100 μm的小尺寸液滴速度快、温度低且不断凝结，趋于力学平衡态，但直径超过500 μm的大尺寸SLD液滴速度慢、温度高且不断蒸发，显著偏离平衡态。

4) 增大试验段气流速度尽管会减弱SLD云雾粒径浓度分层程度，但是却会显著增强动量分层和热分层程度，尤其在160 m/s工况下，SLD云雾(最大液滴直径小于1 000 μm)会均匀分布在构型出口中心区域内(-0.75 m

综上所述，针对目前结冰风洞构型，SLD云雾在试验段内会出现粒径浓度分层、动量分层和热分层现象，这些现象会导致SLD云雾偏离分布均匀且力学平衡的试验状态，影响SLD云雾的积冰形貌，因此评估这种分层的SLD云雾对模型积冰形貌特征的影响，建立SLD冰形相似准则，是下一步亟待开展的研究。此外，结冰风洞中的SLD云雾演化过程十分复杂，涉及动量传递、质量传递和能量传递等多物理过程的强耦合作用，而本文采用的物理模型和计算方法是基于一定假设发展而来的，缺乏全面的结冰风洞试验验证，因此开展结冰风洞云雾演化特性试验研究也是本文下一步的研究重点。