结冰风洞试验混合相似转换方法及其验证

刘 宇，易 贤，王 强，李维浩

(中国空气动力研究与发展中心，绵阳 621000)

0 引 言

飞机结冰是飞行安全的主要威胁之一[1]，为了评估结冰对飞行性能的影响，需要对飞机进行结冰风洞试验。若结冰风洞受试验段尺寸限制无法进行全尺寸模型试验，或参考工况超过结冰风洞的模拟能力范围，需要对试验工况进行相似转换。相似参数是描述结冰过程的重要无量纲数，相似转换是指利用相似参数计算合适的试验工况，使试验可在风洞能力范围内进行，并获得与参考工况相同的冰形。但是在实际试验过程中发现，现有相似转换方法的结果很多时候超出风洞的实际能力范围：在参考工况速度较高或液态水含量较小时，往往转换结果的风速更高或液态水含量更小，导致转换后的参数在风洞中难以实现。因此需要开展更具实用性的相似转换方法研究。

国外的相似转换研究始于20世纪50年代，对相似参数的研究已经相对完备，通过不同相似参数的组合形成了多种相似转换方法。目前常用的转换方法为Ruff方法[2]（或称AEDC方法），以及Ruff与Anderson提出的改进Ruff方法[3]。法国宇航局的Charpin等提出了ONERA方法[4]。还有LWC×Time方法[5]和Olsen方法[6]，可以有效对液态水含量（LWC）进行转换。单独对温度、速度和压力进行的转换方法[7-9]由于没有试验验证很少被采用。国内对相似转换的研究多在已有方法上拓展现应用场景，如针对过冷大水滴[10-11]转部件[12-13]的相似转换方法，对转换方法的优化和改进较少。

国外现有的相似转换方法对霜冰条件的转换已经成熟，明冰条件的结冰机理还缺乏理解[14]，转换效果有限。对于转换速度如何确定还存在一些讨论，主要集中在韦伯数We特征长度的定义上，可能作为特征长度的有：翼型前缘半径[15]、水滴直径[16]、水膜厚度[17]等，对于如何合理确定转换速度仍未得到合理解释。其中Kind[17]和Feo[18]等从表面水膜流动和水滴表面张力等方面对Ruff方法进行了修正。在实际应用上国外以Ruff方法为主，近年关于ONERA方法的改进和验证的文献相对较少。

国外对结冰相似准则的研究随应用需求转入过冷大水滴[19]、后掠翼[20-21]、发动机冰晶结冰[22]和防冰试验相似转换[23]等研究领域，对传统的相似转换理论研究逐渐降温。

为拓展相似转换方法的适用范围，提高其工程实用价值，发展适合我国大型结冰风洞的相似转换方法，本文对相似参数及其变化规律进行了深入分析，根据相似参数和来流条件的关联性，提出了一种混合相似转换方法。该方法根据满足的相似参数不同分为四种模式，每种模式的转换参数选择范围不同。使用混合转换方法对某机翼结冰试验的工况进行了相似转换，并通过数值模拟[24]计算了试验工况、参考工况以及被普遍使用的Ruff方法转换工况的冰形，并对试验和计算结果进行了对比。

1 结冰试验相似转换理论

1.1 相似参数

为保证水收集系数相似，需要满足水滴运动轨迹相似。在流场相似的前提下满足水滴惯性系数K0相等即可认为满足水滴轨迹相似。其定义如下[25]：

其中，λ/λstokes是平均阻力比，K是未做修正的惯性系数[26]。

模型表面水收集总量受水收集系数、LWC和结冰时间影响，用积聚系数Ac表示：

其中，τ为结冰时间，ρi为冰密度，D为前缘直径，V为来流速度。

根据Messinger模型[27]的能量守恒方程，可以推导出模型表面冻结系数：

其中，Λf为结冰相变潜热，cp,ws为水的定压比热容。另外有描述热平衡过程的三个参数：相对热系数b、水滴能量传递势φ和空气能量传递势θ。定义如下：

其中，ts为壁面温度，tf为水膜温度，pw为水蒸汽分压力，pww为饱和水蒸气压力，Λf为水的蒸发潜热，hc为对流换热系数，hG为对流传质系数。

翼型驻点位置冻结系数n0也可通过公式（3）计算。通常认为在流场和水收集系数相似的前提下，满足n0相等则其余部分的冰生长过程也相似。

1.2 常见相似转换方法

常见的相似转换方法有ONERA方法、Ruff方法和Olsen方法。以上方法均有一定试验数据进行支撑，但每种方法均有一定局限性。其中ONERA方法根据Modane风洞的硬件条件设计，没有将压力作为可调整变量[28]纳入转换方法中。Ruff方法建议使用WeL,s=WeL,r条件[19-20]确定Vs， 使转换后速度增加，不利于开展试验。WeL的定义为：

LWC×Time方法仅适用于温度较低的霜冰工况，Olsen方法的部分试验结果显示其转换结果欠佳[29]。

2 相似参数的变化规律

2.1 相似参数的影响因素

通过上节中相似参数定义分析，可以定性的得出以下推论：1）K0主要受 δ和V影响，温度和压力变化对K0影响相对较小；2）当K0,s=K0,r条件满足时，n0仅受到LWC和Tst影响，而三个热平衡参数相等不是n0,s=n0,r的必要条件；3）相对热系数b是LWC和V的函数；4）水滴能量传递势φ是Tst和V的函数。

由以上推论可将相似参数分为三组：一是K0，通过相等条件K0,s=K0,r限制δs和Vs；二是n0，通过n0,s=n0,r限制LWCs和Tst,s；三是b、φ、θ，对来流条件的限制较弱，即使这三个参数均不相等或部分相等，相似转换依然有效。以下将针对这三组参数进行分析。

2.2 惯性系数的变化规律

惯性系数K0的经验计算公式如下，该公式在同类型研究中得到广泛应用：

其中，k和kl均为常数，取值方法可参见文献[2]。实际试验以水滴中值体积直径（MVD）代替公式中的单粒径水滴直径δ。忽略温度和压力对空气粘性和密度的影响，根据公式（8）绘制在满足K0,s=K0,r条件下，K0,s关于V和MVD的变化曲面如图1所示，参考工况已在图中标注。从曲线变化趋势可以看出，MVD随V增大而减小，在K0,s取不同值时曲线变化趋势一致。

图1 K0,s 关于速度和MVD的曲面关系Fig. 1 Curvature plane of correlation between velocity and MVD in constantK0,s

2.3 冻结系数的变化规律

根据3.1节的推论（2）可知，n0,s=n0,r是 较弱的约束条件，对任意状态，理论上有无数LWC和Tst组 合满足该条件。将公式（4 ～ 6）带入到公式（3）并化简，可以得到：

由于计算n0的公式展开后有若干物性参数，使用A1、A2、A3简化展开公式产生的常数项，其中部分如空气密度和静压与温度直接相关的物性参数在简化过程被作为常数处理。

根据公式（9）绘制在满足n0,s=n0,r条件下，n0,s关于LWC和Tst变化的曲面如图2所示，参考工况已在图中标注，其中n0,s坐标轴取对数。从曲线的趋势可以看出，Tst随LWC增加而减少，n0越大可供选择的LWC范围越小。

图2 n0,s 关于LWC和静温的曲面关系Fig. 2 Curvature plane of correlation between static temperature and LWC in constantn0,s

2.4 热平衡参数和水积聚系数的变化规律

为了保证相似转换前后的冰形尽可能相似，Ruff方法和ONERA方法均在n0相等的基础上约束了φ或b相等，以期获得与参考工况相似的热平衡过程。为了进一步分析n0关于来流条件的变化规律，需要分别对b、φ、θ进行分析。

相对热系数b是在n0,s=n0,r条件的基础上对LWC的进一步约束，展开对流换热系数，公式（5）可以写成：

其中，ka为空气导热系数，Pra为普朗特数，通常取0.7。满足了水滴惯性系数相等条件时式中β0可作为常数，即式中除LWC和V以外均为常数。

依据公式（10）绘制在满足bs=br条件下，bs关于LWC和V的空间曲面如图3所示，参考工况条件已在图中标注。从曲线的规律可以看出，LWC随V的增大而单调减小。

图3 bs 关于LWC和速度的曲面关系Fig. 3 Curvature plane of correlation between velocity and LWC in constantbs

同理根据公式（5）做出φs关于Tst和V的空间曲面如图4所示，参考工况已在图中标注。从方程和投影曲线均可以看出，Tst和V是单调的线性关系，斜率与φs的值无关。

图4 φs 关于静温和速度的曲面关系Fig. 4 Curvature plane of correlation of static temperature and velocity in constantφs

由以上分析可以知道，在相似转换方法仅满足ns=nr的弱约束条件时，LWC和Tst存在非线性的函数关系；在引入bs=br或φs=φr条件后，LWC和Tst实际上是两个完全不相关的参数。因此在速度确定的情况下，可计算出同时满足bs=br和φs=φr的LWCs和Tst,s。此时θs=θr成了满足冻结系数相似的充要条件。

在没有其他约束条件作为前提时，θ受到多个来流条件的影响，无法通过 θ约束单一来流条件。但在其他参数已经确定的情况，可求解Pst使 θ满足相似条件，分析 θ关于Pst的变化规律具有一定的实用价值。

根据公式（6），当其余来流条件确定时，θ关于Pst的变化规律如图5所示。需要额外注意的是，根据Ruff对 θ的分析，当Vs大于Vr时，满足θs=θr的静压Pst,s小于Pst,r，反之亦然。这个规律意味着当参考工况的静压指定为当地压力时，更小Vs会导致相似转换后的试验压力高于当地压力。

图5 不同条件下θ 关于静压的变化曲线Fig. 5 Curve of correlation between static pressure and θ in different condition

经过对水滴惯性系数、冻结系数以及三个热平衡参数的分析，已经确定了来流的主要参数，仅剩结冰时间。根据Acs=Acr计算相应的结冰时间即可满足相似。

2.5 小 结

相似参数可分为三组：约束 MVDs和Vs的水滴惯性系数K0；约束LWCs和Tst,s的驻点冻结系数n0；最后一组为热平衡参数b、φ和θ，对试验参数没有直接约束，即使转换前后没有满足相等关系也可获得相似的冰形。分析了三组相似参数关于来流参数的变化规律，并绘制了当满足相似参数相等时来流参数之间的关系曲线。

3 混合相似转换方法

3.1 实施流程

一种混合相似转换方法的流程如图6所示，具体实施流程为：

1）指定速度Vs，根据K0条件计算MVDs，Vs可通过WeL条件计算获得，若计算所得MVDs超出风洞能力包线，可参考图1中MVD-V的曲线，重新调整Vs；

2）根据φ和b的相似条件，分别计算Tst,s和LWCS，当计算结果超出风洞能力包线，可参考图3和图4中两个参数分别关于V的曲线关系，重新调整Vs并重复步骤1；

图6 混合转换方法流程图Fig. 6 Flow chart for hybrid scaling method

3）a. 如果风洞具备连续调节压力的能力，可根据公式（7）计算Pst,s，使之满足 θ条件。由于Pst,s改变导致K0,s出现变化，进一步导致bs不相等，可再次调整MVDs和LWCs使之重新满足。此时为混和相似转换方法模式1；

b. 如果风洞无法连续调节压力，则在b和φ条件中选择一个，并通过n0条件计算LWCs或Tst,s。此时为混合相似转换方法模式2或模式3；

c. 当模式2或模式3所计算的参数超出风洞能力范围时改用Olsen方法，指定LWCs或Tst,s并通过n0条件计算另一个。此时为模式4；

d. 根据Ac条件计算结冰时间τs，并再次检查各参数是否超出风洞能力包线，转换结束。

3.2 四种模式

根据满足的相似参数不同，混合相似转换方法又区分为四种模式。混合相似转换方法与常见方法的对比见表1。

表1 相似转换方法满足相似参数表Table 1 Scaling parameters fitted by different scaling method

模式1同时满足了n0相关的三个热平衡参数，约束最多，理论上可以获得最接近参考冰形的转换结果；模式2可视为改良的ONERA方法，模式3与修正的Ruff方法类似；模式4直接使用了Olsen方法，在模式4下参数选择范围最大，任意参考工况均能找到试验能力范围内的转换结果。

根据需要转换的参考工况和试验设备的能力在四种模式中进行选择：试验设备具有高度模拟的能力可以考虑使用模式1；使用模式2或模式3的参数选择范围相对自由，同时也有大量的文献和试验结果支撑其转换效果，可信度最高；模式4最具泛用性，考虑到有公开文献验证Olsen方法的转换结果并不总是令人满意，特定情况下模式4下转换效果可能较差。

4 方法验证

4.1 数值验证

分别使用混合相似转换方法的四种模式，对文献[30]中的试验工况进行了相似转换，使用数值模拟方法计算了转换工况的冰形，并与文献中的试验结果进行了对比。参考状态使用弦长0.533 4 m的NACA0012翼型，转换状态的试验和计算均使用0.266 7 m的缩比模型。算例如表2所示。

图7展示了两组不同速度的计算结果与试验的对比结果，试验与计算结果基本吻合，主要差异在于：转换工况的驻点冰厚较试验更低；下冰角角度和厚度基本一致，但上冰角吻合程度稍差；结冰极限相对试验更靠近前缘。考虑到数值模拟采用单粒径水滴和单步法计算结冰，与试验结果的误差在可接受范围内。

从验证结果来看，混合相似转换方法的四种模式均能有效获得与参考工况相近的转换结果。速度不同时转换结果的冰形存在一定差异，但区别并不明显。

4.2 试验验证

试验在中国空气动力研究与发展中心（CARDC）的3 m×2 m结冰风洞进行。该风洞试验能力覆盖液态水含量0.2～2.0 g/m3，水滴中值体积直径15～50 μm的云雾参数范围，主试验段最大风速210 m/s。

使用混合转换方法，对某型飞机的结冰试验工况进行了相似转换，具体工况见表3至表5。其中工况1是因模型存在缩比而进行相似转换；工况2和工况3是因原定试验参数超出结冰风洞试验能力而进行相似转换。试验工况均采用混合转换方法所得，在表中以#号标出。通过热刀法截取的翼型中截面冰形。对试验工况、原始工况和其他相似转换工况进行了数值模拟，并对比了试验结果和计算结果。

表2 混合相似转换方法验证算例Table 2 Validation cases by hybrid scaling method

工况1中对比了Ruff方法和混合相似准换方法，如图8所示。从计算所得冰形结果分析，Ruff方法和混合转换方法均能获得与参考冰形类似的结果，且均与试验结果对比较好。但Ruff方法的转换结果风速高，液态水含量小，超过了试验设备的能力范围。

图7 混合相似转换方法计算验证结果图Fig. 7 Ice shape results comparison with validation cases by hybrid scaling method

表3 工况1参数表Table 3 Reference and scaled conditions for test 1

表4 工况2参数表Table 4 The reference and scaled conditions for test 2

表5 工况3参数表Table 5 Reference and scaled conditions for test 3

工况2和工况3来流条件相同，但所用翼型不同，对比结果分别如图9、图10所示。转换前后的冰形基本相同，与试验结果对比，冰形特征基本吻合。其中工况3计算冰形的冰角角度和下冰角厚度与试验结果有表明显的差异，其来源可能是没有对冰增长过程进行多步法求解。

图8 工况1试验与计算结果冰形对比Fig. 8 Ice shape results comparison with validation cases for test 1

图9 工况2试验与计算结果冰形对比Fig. 9 Ice shape results comparison with validation cases for test 2

图10 工况3试验与计算结果冰形对比Fig. 10 Ice shape results comparison with validation cases for test 3

5 结 论

本文分析了相似参数随试验参数的变化规律，提出了一种混合相似转换方法，通过试验和数值模拟进行了对比验证，得到以下结论：

1）使用数值模拟方法验证了混合相似转换方法的四种模式，对比文献试验结果，四种模式均与试验吻合良好，验证了本文方法的有效性。又以混合相似转换方法对某机翼的结冰试验进行了相似转换，对比了试验结果和计算结果，初步验证了本方法的实用性；

2）混合相似转换方法相比传统的相似转换方法，提供了更宽的适用范围，并根据参考状态的需求确定不同的初始转换条件，选择对应的转换模式，能够更好的解决实际工程问题；

3）混合相似转换方法的四种模式中，模式1创新的提出了一种使三个热平衡参数同时满足相等的转换方法，相比经典的Ruff方法和ONERA方法只能满足三个热平衡参数中的一个，其转换结果在理论上更接近参考状态的结冰热力学过程；

4）后续工作一方面会继续在理论上拓展相似转换方法，通过考虑水膜流动对结冰过程的影响，对来流风速的选择缺乏有力的理论约束；另一方面将积极开展结冰风洞标模试验，进一步验证混合相似转换方法的有效性。