冰形预测常用的方法有CFD建模和结冰风洞实验等。CFD建模求解积冰过程一般是采用计算机对合理简化后的数学模型进行求解，包含流场求解、水滴运动求解和热力学过程求解[2-4]。大量文献表明[12-18]，CFD可以在一定程度上较好地模拟积冰过程，为研究结冰提供一定的依据。结冰风洞实验则是把模型安装在实验段内，通过模拟结冰环境来观察积冰增长过程和测量冰形[19-21]。相比较而言，CFD建模具有成本低、数据丰富、研究范围广的优点，但是由于建模过程中模型的误差等原因，对于明冰、过冷大水滴环境等包含明显水滴水膜动力学现象的结冰过程其模拟能力还有进一步提高的空间。另一方面，尽管结冰实验成本较高，参数选取范围很大程度上依赖于结冰风洞的性能，但是它所得的数据是基于物理模拟的，可以较好地反映出结冰问题的本质，同时也能真实地重现结冰过程。当前，受限于结冰风洞的性能和尺寸，全尺寸模型结冰实验是很难实现的，多采用缩比模型进行结冰实验，即采用与真实飞行器几何相似的而实际尺寸较小的模型安装在结冰风洞实验段中进行结冰模拟。
然而，风洞实验作为一种模型实验，必须满足相似性要求才能保证通过风洞实验对真实问题进行准确的物理模拟，否则风洞实验所得的数据是缺乏说服力的[21-24]。国内外对结冰风洞实验的相似准则进行了全面的研究。美国国家航空航天局(NASA)[25]从几何相似、流场相似、水滴运动轨迹相似、水滴收集相似、结冰表面热力学过程相似、结冰表面-水动力学相似这几个方面通过大量的数据分析和方程分析提出了结冰实验的相似准则，并通过结冰实验进行了验证，表明其相似准则具有一定的准确性。在国内，中国空气动力研究与发展中心[16, 26-27]、西北工业大学[28-29]也进行了结冰相似准则的探究，且进行了相应的验证，表明其相似准则的正确性。然而，目前尚无采用完整的相似分析方法导出结冰风洞实验相似准则的报道。
从相似理论出发，导出相似准则的方法有方程分析法和量纲分析法，前者适用于描述物理过程的控制方程已知的情况下对相似准则的推导，当研究对象的物理过程尚不完全清楚的情况下，量纲分析法则更有利于推出相似准则[30]。显然，结冰问题包含了比常规的空气动力学问题更为复杂的物理现象，如水滴撞击、水流相变和积冰的动态增长等，目前对于这些问题的认识还有很多不足，需要进一步加深。因此，结冰问题的相似准则更适合于用量纲分析方法进行导出。
本文在上述大量文献和文献[30]的启发下，采用量纲分析法对结冰问题进行分析，导出其相似准则，并在此基础上提出结冰实验中运行参数的选取方法，最后采用CFD方法进行了验证。
1 结冰相似准则的导出 结冰问题可以表述为：云层中含有不同尺度的过冷水滴，飞行器在穿越云层的过程中这些过冷水滴不断撞击到结冰表面，导致积冰量不断增长，同时积冰反过来不断改变外部的绕流、影响水滴的运动轨迹和撞击位置这样一个动态的过程。
如图 1所示，结冰过程涉及到的具体现象有：气流的流动、水滴的运动、水滴与结冰表面的撞击、大量水滴撞击后在空气动力和重力以及黏附力作用下形成的水膜的运动、水滴和水膜的相变等。可见结冰实验涉及的物理量远多于常规风洞实验中的物理量，可归纳如下：
 |
| 图 1 飞行结冰问题示意图 Fig. 1 Schematic diagram of in-flight icing problem |
| 图选项 |
1) 描述空间的变量：积冰尺度(用积冰的厚度δi(s, t)表示)，过冷水滴的尺度(用过冷水滴的直径δ表示)，过冷水滴的平均间距d, 飞行器的空间尺度(用飞行器特征尺度c表示)，结冰表面的几何特征(用结冰表面的曲线坐标s表示, 坐标原点O为翼形前缘)，水膜尺度(用水膜厚度hfilm表示)，以及水滴的空间位置(用水滴的空间位置矢量X(xi), i=1, 2, 3表示。
2) 描述时间的变量：积冰时间ti。
3) 描述结冰过程中水、空气和积冰以及由它们构成的流动系统性质的变量：黏性(用水和空气的黏性系数νw和νa表示)，密度(水、空气和冰的密度ρw、ρa和ρi)，热传导系数(水、空气和冰的热传导系数kw、ka和ki)，比热容(水、空气和冰的比热容cp, w、cp, a和cp, i)，潜热(水的凝固潜热hf)，对流换热系数h。
4) 描述绕流流场参数的变量：速度(来流速度U∞，流场的当地空气气流速度ua(ua, i), i=1, 2, 3，过冷水滴的运动速度ud(ud, i), i=1, 2, 3，水膜的流动速度uf(uf, i), i=1, 2, 3，压力(空气气流的静压p)，温度(空气流场的温度Ta、水滴的温度Td)。
5) 描述结冰表面性质的变量：温度(结冰表面的温度Ts)，界面张力(水-气界面张力系数σw/a、水-固界面张力系数σw/s)。
6) 其他变量：重力(重力加速度g)，冰点温度Tfre。
因此描述结冰现象的一般关系式可写成
 | (1) |
考虑到结冰问题本质上是一个力学附带热交换问题，那么在国际单位制(SI)下，该问题涉及到4个基本物理量，即长度、质量、时间和温度，它们的量纲可作为基本量纲，分别用符号L、M、T、Θ表示，其他物理量的量纲均可由基本量纲的幂次单项式来表示，这样，在结冰问题中，相关物理量的意义及其量纲可以列在表 1中。
表 1 结冰问题中相关物理量及其量纲 Table 1 Relevant physical variables and their dimensions in icing problem
| 物理量性质 | 参数 | 意义 | 量纲 |
| 空间 | δi(s, t) | 积冰厚度 | L |
| δ | 过冷水滴直径 | L | |
| d | 过冷水滴平均间距 | L | |
| c | 飞行器特征尺度 | L | |
| s | 模型表面曲线坐标 | L | |
| hfilm | 水膜厚度 | L | |
| X(xi), i=1, 2, 3 | 水滴的空间位置矢量 | L | |
| 时间 | ti | 积冰时间 | T |
| 物性参数及系统参数 | νa | 空气的运动黏性系数 | L2T-1 |
| νw | 水的运动黏性系数 | L2T-1 | |
| ρa | 空气的密度 | ML-3 | |
| ρw | 水的密度 | ML-3 | |
| ρi | 冰的密度 | ML-3 | |
| ka | 空气的热传导系数 | MLT-3Θ-1 | |
| kw | 水的热传导系数 | MLT-3Θ-1 | |
| ki | 冰的热传导系数 | MLT-3Θ-1 | |
| cp, a | 空气的比热容 | L2T-2Θ-1 | |
| cp, w | 水的比热容 | L2T-2Θ-1 | |
| cp, i | 冰的比热容 | L2T-2Θ-1 | |
| hf | 水的凝固潜热 | L2T-2 | |
| h | 对流换热系数 | MT-3Θ-1 | |
| 流场参数 | U∞ | 来流速度(飞行速度) | LT-1 |
| ua(ua, i), i=1, 2, 3 | 空气气流速度 | LT-1 | |
| ud(ud, i), i=1, 2, 3 | 过冷水滴速度 | LT-1 | |
| uf(uf, i), i=1, 2, 3 | 水膜速度 | LT-1 | |
| p | 空气气流静压 | ML-1T-2 | |
| Ta | 空气流场的温度 | Θ | |
| Td | 水滴的温度 | Θ | |
| 结冰表面参数 | Ts | 结冰表面温度 | Θ |
| σw/a | 水-气界面张力系数 | MT-2 | |
| σw/s | 水-固表面张力系数 | MT-2 | |
| 其他变量 | g | 重力加速度 | LT-2 |
| Tfre | 冰点温度 | Θ |
表选项
由式(1)或表 1可知，结冰问题涉及33个变量，即总变量个数N=33，而基本量纲数k=4，那么根据Buckingham Π定理，进行量纲分析，式(1)将变成包含N-k=29个无量纲量的方程，即
 | (2) |
式中：29个无量纲参数为
上述29个无量纲参数可按其不同的物理效应范畴分类列于表 2。其中：Π1为无量纲结冰厚度，Π2和Π4分别为表示水滴和模型几何相似的参数，Π3为液态水含量LWC与空气密度的比值，Π5为无量纲水膜厚度，Π6为表示过冷水滴空间位置相似的参数，Π7为无量纲积冰时间，Π8、Π9、Π10、Π11、Π12、Π13、Π15、Π16、Π27、Π29为物性参数的比值，Π11为普朗特数，Π14为马赫数，Π17为雅克伯数，Π18为努塞尔数，Π19为来流雷诺数，Π20、Π21、Π22、Π24、Π25为表示流场几何相似的参数，Π23为欧拉数，Π26为韦伯数，Π28为弗劳德数。
表 2 结冰问题中的无量纲参数 Table 2 Dimensionless parameters in icing problem
| 参数 | 意义 |
 | 无量纲结冰厚度 |
 | 过冷水滴尺寸与模型尺寸之比 |
 | 液态水含量与空气密度的比值 |
 | 模型几何相似 |
 | 无量纲水膜厚度 |
 | 水滴空间位置相似 |
 | 无量纲积冰时间 |
 | 物性参数的比值 |
 | 普朗特数 |
 | 马赫数 |
 | 雅克伯数 |
 | 努赛尔数 |
 | 来流雷诺数 |
 | 速度场相似 |
 | 欧拉数 |
 | 温度场相似 |
 | 韦伯数 |
 | 弗劳德数 |
表选项
在结冰实验中，Π1由其他参数决定，是因变量；Π2由水滴尺寸的选取决定；Π3由LWC的选取决定；当来流温度和压力保持不变时，空气和水的物理性质保持不变，且模型的几何相似是必须满足的，因此这种情形下Π4、Π8、Π9、Π10、Π11、Π12、Π13、Π15、Π16、Π27、Π29自然满足相似性，如果温度和压力也进行对应的缩比，但一般情况下缩比后的温度和压力值并不会发生很大的变化，对空气、水的物理性质(黏性、密度、热传导系数、界面张力系数)的影响是非常有限的，近似仍然可以认为Π4、Π8、Π9、Π10、Π11、Π12、Π13、Π15、Π16、Π27、Π29仍然满足相似性，因此它们可从式(2)中删去；Π5由绕流流场决定，可认为与Π3、Π19相关[19]，因此Π5也可从式(2)中删去；Π6和Π21与水滴的运动过程有关；Π7由实际积冰时间决定；Π14由来流速度和温度决定；Π17与来流温度的选取有关；Π18取决于模型的特征尺度；Π19取决于模型的特征尺度和来流速度的选取；Π20与Π14等价，因此Π20可从式(2)中删去；Π22可认为由Π14、Π19决定，因此Π22也可从式(2)中删去；Π23与Π14、Π19有关，因此Π23可从式式(2)中删去；Π24、Π25由来流温度和模型表面温度以及Π18决定，因此Π24、Π25也可从式(2)中删去；Π26由来流速度和水滴尺寸决定；由于重力加速度可认为为常数，Π28由水滴尺寸和来流速度的选取决定。那么，式(2)可简化为
 | (3) |
式(3)可改写成
 | (4) |
即:无量纲的结冰厚度由无量纲水滴尺寸、液态水含量与空气密度的比值、无量纲水滴空间位置、无量纲积冰时间、马赫数、雅克伯数、努塞尔数、雷诺数、无量纲水滴速度、韦伯数、弗劳德数决定。式(4)便是表示结冰过程的一般无量纲关系式。
2 结冰相似准则的简化 考虑到无量纲的水滴运动方程[2-4, 19]为
 | (5) |
式中：
 | (6) |
 | (7) |
其中：CD为球形水滴在空气中运动的真实阻力系数；Rerel和K分别为基于当地液滴与空气相对速度的相对水滴雷诺数和水滴的无量纲惯性参数；ua为气流的无量纲速度；X为水滴的无量纲位置矢量；t为无量纲时间。由式(5)可知，水滴的运动轨迹可以完全由K和CDRerel确定。如果这2个无量纲参数在缩比模型实验和未缩比情形中保持一致，则可以认为它们的水滴运动轨迹相似。然而在实际应用中，这2个参数常常不能同时匹配。为了克服这一困难，提出了一个综合反映K和CDRerel影响的单一参数K0，只要它匹配，水滴运动就相似。文献[25]给出的K0表达式为
 | (8) |
它只与水滴惯性参数K和水滴雷诺数Reδ有关。水滴雷诺数是基于来流速度的雷诺数，即
 | (9) |
显然，水滴雷诺数Reδ可以写成
 | (10) |
因此式(8)可写成
 | (11) |
而水滴的无量纲惯性参数可重新组合成
 | (12) |
因此可得
 | (13) |
这样，式(4)中Π6和Π21可以用Π2和Π19替换，从而式(4)可改写成
 | (14) |
即
 | (15) |
这就表明影响结冰冰形的主要相似参数有水滴与飞行器尺寸的比值、液态水含量与空气密度的比值、无量纲积冰时间、马赫数、雅克伯数、努赛尔数、雷诺数、韦伯数以及弗劳德数。
3 本文相似准则与传统准则的关系 式(14)尽管从严格的理论上得出了影响结冰冰形的主要无量纲参数，但是在实际应用中，很多无量纲参数是无法同时满足相似性条件的，严格按照式(14)来选取结冰实验的运行参数是很难实现的，比如，模型缩比与水滴缩比如果严格按照同一比例，很可能导致水滴尺寸过小而无法发生与结冰表面的撞击，或马赫数和雷诺数往往无法同时满足相似性要求。
传统的相似准则[16, 25-29]基本上采用分别保证几何相似、流场相似、水滴运动相似、水滴收集相似、水滴-结冰表面动力学相似、结冰表面热力学相似来达到结冰相似的目的。这几方面也被大量实验证实确实能够很好地反映影响结冰过程的几乎所有因素，但是传统相似准则并没有采用无量纲参数来严格描述上述几方面要求。下面根据上述几个方面对所得出的无量纲参数进行分类，找出其作用的结冰相似要求的具体方面，将传统相似准则的要求重新写成本文量纲分析得出的无量纲参数的组合，这样就将每个物理量的作用具体化，也更便于导出结冰实验风洞的运行参数。
3.1 几何相似 几何相似条件对应Π4，只要满足Π4则可保证几何相似。
3.2 流场相似 流场相似需要满足马赫数和来流雷诺数保持不变，即Π14和Π19保持不变。马赫数一般非常容易满足相似性条件，但是雷诺数很难满足，只有特殊的具备降温、增压能力的风洞才能做到马赫数和来流雷诺数同时满足相似性。
对于结冰风洞实验，温度是一项关键运行参数，它的剧烈变化很可能导致结冰的结果发生质变，另外，结冰问题一般属于低速空气动力学的范畴，结冰部位都位于前缘，这些地方边界层厚度很小，因此只要满足马赫数相等即可。
3.3 水滴运动相似 从式(13)可知，要满足水滴运动相似，则需要满足式(16)。当且仅当式(17)和式(18)同时成立，才能保证式(16)成立。
 | (16) |
 | (17) |
 | (18) |
式中：下标“S”表示缩比情形下对应的参数，无下标则表示未缩比情形的参数，下文中也相同。考虑到如果满足水滴尺寸缩比相似，即式(17)成立，则由于风洞实验段参数的限制(U∞和νa的调节会引起其他无量纲量如马赫数、韦伯数等发生变化)必然导致雷诺数减小，使式(18)无法满足，从而即式(16)不能满足，K0不能相等。也就是说，用K0来描述水滴的运动相似尽管精度更高，但是在缩比实验中几乎不可能满足。
现在考察K0与K的关系，即式(11)，也就是要考察函数F(Reδ)的特性。图 2给出了在常规的结冰参数范围内不同水滴直径下F(Reδ)随来流速度的变化曲线。由图 2可见，F(Reδ)的值基本上介于0.01~0.04，并且在实际飞行速度或缩比实验速度范围，F(Reδ)变化平缓，因为缩比模型实验速度和真实飞行速度相差不会很大，所以引起F(Reδ)的差别较小。图 3给出了不同来流速度下F(Reδ)随水滴直径的变化，其变化类似于它随来流速度的变化，即水滴直径的差异不会引起F(Reδ)很大的变化。由此可知，F(Reδ)不会随来流速度或水滴直径发生很大变化。因此可以用K代替K0近似表征水滴运动的相似性，即只要保证K在缩比情形下和未缩比时相等，就近似认为水滴运动相似。
 |
| 图 2 F(Reδ)随来流速度U∞的变化 Fig. 2 F(Reδ) versus freestream velocity U∞ |
| 图选项 |
 |
| 图 3 F(Reδ)随水滴直径δ的变化 Fig. 3 F(Reδ) versus droplet diameter δ |
| 图选项 |
3.4 水滴收集相似 表面水滴收集相似要求结冰表面水滴收集量相似[25]，即
 | (19) |
式中：β和η分别为局部水滴收集系数和冻结比例。考虑到
 | (20) |
那么，只要保证Π3、Π7、Π10、β、η的乘积在缩比条件下与未缩比时相等，就可满足式(19)，其中Π10为常量。
3.5 水滴-结冰表面动力学相似 水滴撞击到结冰表面之后，会发生铺展、融合、飞溅等动力学过程，描述该过程最重要的参数就是韦伯数，即Π26。只要满足
 | (21) |
则可保证水滴-结冰表面动力学相似。
3.6 结冰表面热力学相似 水滴撞结冰表面的热力学过程是一种典型的边界层对流换热现象，根据经典的热力学理论[31]，边界层对流换热现象可表述为
 | (22) |
式中：s为无量纲的物面位置。对应于本文中的无量纲参数可写成
 | (23) |
因此只要保证Π4、Π9、Π11保持不变，即可满足热力学相似。
综上可见，传统的相似准则中的要求都可以表示为本文得出的无量纲参数的函数，它们之间的关系如表 3所示。
表 3 传统相似准则要求与本文相似准则的关系 Table 3 Relationship between traditional similarity laws and proposed similarity law
| 传统相似准则的要求 | 相关的无量纲参数 | 应用中满足相似性要求的条件 |
| 几何相似 |  | 采用缩比模型，自然满足 |
| 流场相似 | Π14=Ma | 选取合理的速度和温度，容易满足 |
| Π19=Re∞ | 往往不容易满足，尽量配合其他条件，使其他条件优先满足 | |
| 水滴运动相似 |  | 令K保持不变 |
| Π19=Re∞ | ||
| 水滴收集相似 |  | 选取合理液态水含量，容易满足 |
 | 取合理的积冰时间，容易满足 | |
| β | 水滴运动相似的前提下即可满足 | |
| η | 对流换热相似，依赖于Π18 | |
| 水滴-结冰表面动力学相似 |  | 明冰条件下非常重要 |
| 结冰表面热力学相似 | Π18=Nu | 很难满足，取决于Π4、Π9、Π11 |
 | 采用缩比模型，自然满足 | |
| Π19=Re∞ | 往往不容易满足，尽量配合其他条件，使其他条件优先满足 | |
| Π11=Pr | 自然满足 |
表选项
4 结冰风洞运行参数选取方法 当模型尺寸缩比时，假定
 | (24) |
式中：k0 < 1为缩比的比例，而当风洞的其他运行参数(风速、温度、压力、水滴直径、液态水含量和积冰时间等)均不发生变化时，式(15)中Π3、Π14、Π17、Π26、Π28均可满足相似性要求，而选取合理的结冰时间也很容易满足
 | (25) |
但式(12)中其余的无量纲量满足以下关系式：
 | (26) |
 | (27) |
 | (28) |
这样，显然只要模型的缩比比例k0≠1，Π2、Π18、Π19是无法满足相似性要求的。而Π2、Π19与水滴运动过程紧密相关，Π18与结冰表面的热力学过程紧密相关，这样就不得不考虑不同结冰条件下结冰问题的具体特征来选取风洞的运行参数。
4.1 霜冰情形下实验运行参数的选取 在霜冰情形下，水滴一撞击到结冰表面就立即冻结成冰，结冰过程显然可以不用考虑水滴的动力学效应和表面水膜的情形，只需要考满足流场相似、水滴运动相似和水滴收集相似即可满足相似性。从表 3可知，水滴运动相似要满足
 | (29) |
即
 | (30) |
考察式(30)，显然取
 | (31) |
 | (32) |
便可满足式(29)。
水滴收集相似要求保证满足式(20)。那么，显然只要满足水滴运动相似，就可以保证
 | (33) |
并且在霜冰情形下，水滴撞击后完全冻结，即
 | (34) |
同时可取积冰时间和液态水含量为
 | (35) |
 | (36) |
便可满足Π3和Π7在缩比条件下与未缩比时相等，且保证式(19)成立。
取
 | (37) |
也可保证
 | (38) |
 | (39) |
这样流场相似也得到了满足。
因此，霜冰情形下对应的风洞运行参数选取为
 | (40) |
4.2 明冰情形下实验运行参数的选取 在明冰情形下，水滴的动力学效应和结冰表面的热力学过程对结冰有很关键的影响，它们的相似性要得到满足，同时也应保证水滴运动轨迹的相似、水滴收集相似、流场相似。可取水滴直径的满足
 | (41) |
式中：λ为待定参数。此时要保证结冰表面水滴动力学相似则要满足
 | (42) |
可得
 | (43) |
这时
 | (44) |
 | (45) |
进一步可得
 | (46) |
要满足水滴运动轨迹相似，显然只有取
 | (47) |
因此可得
 | (48) |
 | (49) |
这样水滴的动力学效应和水滴运动轨迹相似性要求得到了满足。
另外，要保证水滴收集相似要满足式(19)，那么显然与霜冰情形下类似，只要再满足缩比条件下和未缩比时冻结比例和液态水含量相等，并取积冰时间
 | (50) |
就可以满足Π3, S=Π3，Π7, S=Π7且水滴收集相似。同时取
 | (51) |
也可保证
 | (52) |
但是，要保证冻结比例相等就必须要保证结冰表面对流传热相似。表征对流传热强度的无量纲数为努塞尔数Nu，即Π18，显然缩比后Nu减小了。考虑式(23)，缩比后，普朗特数不变，而由式(45)可知雷诺数也减小了，而在结冰部位多在迎风部位，这些地方边界层较薄，当地雷诺数一般也低于临界雷诺数，可认为局部的流动状态仍然为层流状态，缩比后雷诺数进一步减小，并未改变结冰表面边界层的流态，因此可认为缩比后尽管对流换热的强度相对减小了，但是并未引发对流换热过程发生质变。另一方面，结冰表面的水膜厚度是非常小的[25]，一般不超过10 μm量级，其相变的时间尺度也是很小的[32]，在空气动力的剪切作用下的流动速度也很小。因此，尽管缩比后冻结比例有变化，但是相比于水滴动力学效应和运动轨迹对积冰冰形的影响，可以认为冻结比例的改变产生的影响较小，近似可以忽略。此外，按上述方法进行参数选取，有
 | (53) |
 | (54) |
式中：马赫数Π14表征流动的压缩性，在结冰现象对应的来流速度范围内，压缩性要求也可以适当放宽，而Π2与Π28与霜冰情形下类似，可以忽略其不满足相似性要求带来的影响。
综上所述，明冰情形下对应的风洞运行参数选取为
 | (55) |
5 基于CFD的结冰实验相似准则验证 为了对本文提出结冰实验的相似准则进行验证，采用FENSAP-ICE软件对理论分析得出的风洞实验运行参数选取方法进行计算对比，具体计算方法参见文献[12]。首先选取NACA0012翼型分别进行霜冰和明冰情形下积冰冰形的预测并与实验结果对比，验证积冰预测数值方法的有效性，然后采用这种方法预测不同尺寸模型的积冰结果，验证本文提出的积冰实验运行参数选取方法的可行性。为了精确地模拟结冰过程，计算采用多步法，即将结冰的总时间分为若干时间步，在每一个时间步Δt内依次进行流场求解、水滴运动及撞击特性求解、积冰增长求解，然后更新网格，这样在完成所有时间步后得到最终的冰形，如图 4所示。
 |
| 图 4 多步法计算流程 Fig. 4 Calculation flowchart of multi-step method |
| 图选项 |
5.1 NACA0012翼型积冰冰形预测 本文积冰计算的参数选取参照了文献[33]中的结冰风洞实验结果并与其实验参数相同，选取了来流温度较低、水滴尺寸较小、液态水含量较低的霜冰情形以及来流温度较高、水滴尺寸较大、液态水含量较高的明冰情形，具体数值如表 4所示。霜冰的计算时间步数n=3，时间步长Δt=140 s，明冰计算时间步数n=2，时间步长Δt=180 s。
表 4 积冰冰形预测的参数 Table 4 Parameters in icing prediction
| 变量 | 霜冰情形 | 明冰情形 |
| 平均气动弦长c/m | 0.533 4 | 0.533 4 |
| 来流速度U∞/(m·s-1) | 102.8 | 67.1 |
| 来流静温T∞/K | 256.49 | 265.07 |
| 来流静压p∞/Pa | 90 760 | 90 760 |
| 迎角α/(°) | 3.5 | 3.5 |
| 液态水含量LWC/(g·m-3) | 0.55 | 1.30 |
| 平均水滴直径MVD/μm | 20 | 30 |
| 积冰时间ti/s | 420 | 360 |
| 积冰时间步数n | 3 | 2 |
| 积冰时间步长Δt/s | 140 | 180 |
表选项
图 5给出了霜冰的计算结果与实验结果的对比。由图可见，相比于实验结果，本文的计算结果在翼型下表面的积冰有一些微小的差异，但是总体上本文的计算结果与实验结果的积冰范围、积冰厚度、冰形形状都吻合良好。这就表明本文的计算方法可以较好地预测霜冰的冰形。
 |
| 图 5 霜冰冰形计算结果与实验结果 Fig. 5 Computed and experimental icing results and experiment results under rime ice situation |
| 图选项 |
图 6给出了明冰的计算结果与实验结果的对比。由图可见，相比于实验结果，本文的计算结果尽管在冰角的生长方向、冰角的高度上有一定的差异，但是在积冰的范围、积冰的冰角位置和驻点附近的积冰厚度都吻合较好，冰形形状也较为接近。这就表明本文的计算方法可以较好地预测明冰的特征，可用于冰形的预测。
 |
| 图 6 明冰冰形计算结果与实验结果 Fig. 6 Computed and experimental icing results and experiment results under glaze ice situation |
| 图选项 |
综上可见，本文的积冰计算方法可以较好地预测积冰冰形，可用于风洞实验运行参数选取方法的验证。
5.2 风洞实验运行参数选取方法验证 根据式(40)和式(55)，本文选取霜冰情形和明冰情形下的全尺寸模型和缩比模型结冰参数分别如表 5和表 6所示。
表 5 霜冰情形相似性验证计算参数选取 Table 5 Calculation parameter selection in verification of similarity under rime ice situation
| 变量 | 全尺寸模型 | 缩比模型 |
| 缩比比例k0 | 1 | 1/3 |
| 平均气动弦长c/m | 0.533 4 | 0.177 8 |
| 来流速度U∞/(m·s-1) | 102.8 | 102.8 |
| 来流静温T∞/K | 256.49 | 256.49 |
| 来流静压p∞/Pa | 90 760 | 90 760 |
| 迎角α/(°) | 3.5 | 3.5 |
| 液态水含量LWC/(g·m-3) | 0.55 | 0.55 |
| 平均水滴直径MVD/μm | 20 | 11.547 |
| 积冰时间ti/s | 420 | 140 |
| 积冰时间步数n | 3 | 3 |
| 积冰时间步长Δt/s | 140 | 46.667 |
表选项
表 6 明冰情形相似性验证计算参数选取 Table 6 Calculation parameter selection in verification of similarity under glaze ice situation
| 变量 | 全尺寸模型 | 缩比模型 |
| 缩比比例k0 | 1 | 1/3 |
| 平均气动弦长c/m | 0.533 4 | 0.177 8 |
| 来流速度U∞/(m·s-1) | 67.1 | 96.775 |
| 来流静温T∞/K | 265.07 | 265.07 |
| 来流静压p∞/Pa | 90 760 | 90 760 |
| 迎角α/(°) | 3.5 | 3.5 |
| 液态水含量LWC/(g·m-3) | 1.30 | 1.30 |
| 平均水滴直径MVD/μm | 30 | 14.423 |
| 积冰时间ti/s | 360 | 83.2 |
| 积冰时间步数n | 2 | 2 |
| 积冰时间步长Δt/s | 180 | 41.6 |
表选项
图 7给出了按表 5参数计算得出的全尺寸模型和1/3缩比模型的冰形对比结果。由图可见，一方面，采用表 3中的参数计算所得的全尺寸模型和缩比模型的冰形形状基本上完全吻合。另一方面，缩比模型的冰形厚度在每个翼面位置处都略微大于全尺寸模型，这是由于在考虑水滴运动轨迹相似时，只保证了水滴无量纲惯性参数K相等，而水滴相对雷诺数在缩比后变小了，那么从水滴运动方程可知水滴在撞击过程中的加速度就减小了，水滴的惯性相对有所增大，从而导致结冰表面的水滴收集量是略有增大的，但是这种变化影响很小。总体上，根据表 5的参数计算所得的全尺寸模型上的冰形和缩比模型的冰形基本上满足冰形相似。这也就说明根据式(40)选取的缩比模型的参数可以得到相似的冰形。
 |
| 图 7 霜冰冰形相似验证的计算结果 Fig. 7 Similarity verification calculation results of icing under rime ice situation |
| 图选项 |
图 8给出了按表 6参数计算得出的全尺寸模型和1/3缩比模型的冰形对比结果。由图可见，一方面，采用表 6中的参数计算所得的全尺寸模型和缩比模型的冰形形状基本上吻合，冰角的位置和高度基本相同，积冰范围也吻合较好。另一方面，缩比模型驻点附近的积冰厚度要略低于全尺寸模型，而在驻点两侧下游的积冰厚度略高于全尺寸模型，并且上表面的积冰范围比全尺寸模型时偏大，下游形成了冰脊。这是因为模型尺寸缩比后，对流换热的强度降低了，水滴和水膜的冻结比例减小了，在驻点两侧的上下表面的水滴或者水膜在气动力的作用下向下游的流动的比例增大了，导致驻点两侧下游的积冰范围和厚度略大于全尺寸模型，而在驻点附近，由于缩比后来流速度增大了，驻点附近的气动加热效应更强了，导致局部的冻结比例相比于全尺寸模型偏小。总体上，根据表 6的参数计算所得的全尺寸模型上的冰形和缩比模型的冰形基本上满足冰形相似。这也就说明根据式(55)选取的缩比模型的参数可以得到相似的冰形。
 |
| 图 8 明冰冰形相似验证的计算结果 Fig. 8 Similarity verification calculation results of icingunder glaze ice situation |
| 图选项 |
综上所述，式(40)和式(55)表述的选取结冰实验的参数可以得到较好地满足积冰实验相似性的要求。
6 结论 1) 通过严格的量纲分析，得出了描述结冰问题的无量纲关系式以及其中的无量纲参数，并逐个分析了它们的意义，得到了影响结冰的主要无量纲参数：水滴与飞行器尺寸的比值、液态水含量与空气密度的比值、无量纲积冰时间、马赫数、雅克伯数、努赛尔数、雷诺数、韦伯数以及弗劳德数。
2) 传统的结冰相似准则的要求都可以表示成本文推出的各个无量纲参数的组合。
3) 在模型缩比时，通过分别满足传统相似准则的几方面要求，并根据其中各个无量纲参数的关系导出了结冰风洞运行参数的选取方法。
4) 通过FENSAP-ICE软件进行结冰风洞运行参数的选取方法进行验证，表明得出的方法是可行的。
参考文献
| [1] | THOMAS S K, CASSONI R P, MACARTHUR C D. Aircraft anti-icing and de-icing techniques and modeling[J]. Journal of Aircraft, 1996, 33(5): 841-854. DOI:10.2514/3.47027 |
| [2] | 裘燮刚, 韩凤华. 飞机防冰系统[M]. 北京: 国防工业出版社, 2004. QIU X G, HAN F H. Anti-icing systems of aircraft[M]. Beijing: National Defense Industry Press, 2004. (in Chinese) |
| [3] | 朱春玲, 朱程香. 飞机结冰及其防护[M]. 北京: 科学出版社, 2016. ZHU C L, ZHU C X. Aircraft icing and its protection[M]. Beijing: Science Press, 2016. (in Chinese) |
| [4] | 林贵平, 卜雪琴, 申晓斌, 等. 飞机结冰与防冰技术[M]. 北京: 北京航空航天大学出版社, 2016. LIN G P, BU X Q, SHEN X B, et al. Aircraft icing and anti-icing/deicing technologies[M]. Beijing: Beihang University Press, 2016. (in Chinese) |
| [5] | 任鹏飞, 徐宇, 宋娟娟, 等. 霜冰条件对翼型气动性能影响数值研究[J]. 工程热物理学报, 2014, 35(4): 663-668. REN P F, XU Y, SONG J J, et al. Numerical study about the influence of rime ice conditions on airfoil[J]. Journal of Engineering Thermophysics, 2014, 35(4): 663-668. (in Chinese) |
| [6] | CAO Y, HUANG J, XU Z, et al. Insight into rime ice accretion on an aircraft wing and corresponding effects on aerodynamic performance[J]. The Aeronautical Journal, 2016, 120(1229): 1101-1122. DOI:10.1017/aer.2016.49 |
| [7] | BRAGG M B, BROEREN A P, BLUMENTHAL L A. Iced-airfoil aerodynamics[J]. Progress in Aerospace Sciences, 2005, 41(5): 323-362. DOI:10.1016/j.paerosci.2005.07.001 |
| [8] | POURYOUSSEFI S G, MIRZAEI M, NAZEMI M M, et al. Experimental study of ice accretion effects on aerodynamic performance of an NACA 23012 airfoil[J]. Chinese Journal of Aeronautics, 2016, 29(3): 585-595. DOI:10.1016/j.cja.2016.03.002 |
| [9] | 李国知, 曹义华. 旋翼结冰对直升机飞行动力学特性的影响[J]. 航空学报, 2011, 32(2): 187-194. LI G Z, CAO Y H. Effects of rotor icing on helicopterflight dynamic characteristics[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(2): 187-194. (in Chinese) |
| [10] | 张恒, 李杰, 龚志斌. 多段翼型缝翼前缘结冰大迎角分离流动数值模拟[J]. 航空学报, 2017, 38(2): 520733. ZHANG H, LI J, GONG Z B. Numerical simulation of separated flow around a multi-element airfoil at high angle of attack with iced slat[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(2): 520733. (in Chinese) |
| [11] | 周莉, 徐浩军, 龚胜科, 等. 飞机结冰特性及防除冰技术研究[J]. 中国安全科学学报, 2010, 20(6): 105-110. ZHOU L, XU H J, GONG S K, et al. Research of aircraft icing characteristics and anti-icng and deicing technology[J]. China Safty Science Journal, 2010, 20(6): 105-110. DOI:10.3969/j.issn.1003-3033.2010.06.018 (in Chinese) |
| [12] | BEAUGENDRE H, MORENCY F, HABASHI W G. FENSAP-ICE's three-dimensional in-flight ice accretion module:ICE3D[J]. Journal of Aircraft, 2003, 40(2): 239-247. DOI:10.2514/2.3113 |
| [13] | IULIANO E, MINGIONE G, PETROSINO F, et al. Eulerian modeling of large droplet physics toward realistic aircraft icing simulation[J]. Journal of Aircraft, 2011, 48(5): 1621-1632. DOI:10.2514/1.C031326 |
| [14] | JUNG S K, SHIN S, MYONG R S, et al. An efficient CFD-based method for aircraft icing simulation using a reduced order model[J]. Journal of Mechanical Science and Technology, 2011, 25(3): 703-711. DOI:10.1007/s12206-011-0118-4 |
| [15] | PETROSINO F, MINGIONE G, CAROZZA A, et al. Ice accretion model on multi-element airfoil[J]. Journal of Aircraft, 2011, 48(6): 1913-1920. DOI:10.2514/1.C031346 |
| [16] | 易贤.飞机积冰的数值计算与积冰实验相似准则研究[D].绵阳: 中国空气动力研究与发展中心, 2007. YI X.Numerical computation of aircraft icing and study on icing test scaling law[D].Mianyang: Chinese Aerodynamic Research and Development Center, 2007(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-90113-2007189443.htm |
| [17] | 雷梦龙, 常士楠, 杨波. 基于Myers模型的三维结冰数值仿真[J]. 航空学报, 2018, 39(9): 21952. LEI M L, CHANG S N, YANG B. Three-dimensional numerical simulation of icing using Myers model[J]. Acta Aeronautica et Astronautica Sinica, 2018, 39(9): 21952. (in Chinese) |
| [18] | 常士楠, 苏新明, 邱义芬. 三维机翼结冰模拟[J]. 航空学报, 2011, 32(2): 212-222. CHANG S N, SU X M, QIU Y F. Ice accretion simulation on three dimensional wings[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(2): 212-222. (in Chinese) |
| [19] | HUANG X, MYERS B, D'AVIRRO J, et al.Icing wind-tunnel icing test on a contaminated full-scale Wing-Model at takeoff conditions: AIAA-2008-6417[R].Reston: AIAA, 2008. |
| [20] | 战培国. 结冰风洞研究综述[J]. 实验流体力学, 2007, 21(3): 92-96. ZHAN P G. A review of icing wind tunnel[J]. Journal of Experiments in Fluid Mechanics, 2007, 21(3): 92-96. DOI:10.3969/j.issn.1672-9897.2007.03.019 (in Chinese) |
| [21] | 李周复. 风洞特种实验技术[M]. 北京: 航空工业出版社, 2010. LI Z F. Techniques in special wind tunnel tests[M]. Beijing: Aviation Industry Press, 2010. (in Chinese) |
| [22] | 李周复. 风洞试验手册[M]. 北京: 航空工业出版社, 2015. LI Z F. Handbook of wind tunnel test[M]. Beijing: Aviation Industry Press, 2015. (in Chinese) |
| [23] | 赵凯华. 定性与半定量物理学[M]. 北京: 高等教育出版社, 2008. ZHAO K H. Qualitative and semi-quantitative physics[M]. Beijing: Higher Education Press, 2008. (in Chinese) |
| [24] | 李之光. 相似与模化[M]. 北京: 国防工业出版社, 1982. LI Z G. Similarity and modeling[M]. Beijing: National Defense Industry Press, 1982. (in Chinese) |
| [25] | ANDERSON D N.Manual of scaling method: NASA/CR-2004-212875[R].Washington, D.C.: NASA, 2004. |
| [26] | 易贤, 周志宏, 杜雁霞, 等. 考虑相变时间效应的结冰实验相似参数[J]. 实验流体力学, 2016, 30(2): 14-19. YI X, ZHOU Z H, DU Y X, et al. An icing scaling parameter with the effects of phase change time[J]. Journal of Experiments in Fluid Mechanics, 2016, 30(2): 14-19. (in Chinese) |
| [27] | 周志宏, 易贤, 桂业伟, 等. 考虑水滴动力学效应的结冰实验相似准则[J]. 实验流体力学, 2016, 30(2): 20-25. ZHOU Z H, YI X, GUI Y W, et al. Icing scaling law with the dynamics effects of water droplets[J]. Journal of Experiments in Fluid Mechanics, 2016, 30(2): 20-25. (in Chinese) |
| [28] | 李静.发动机进口静止/旋转部件结冰及结冰实验的相似研究[D].西安: 西北工业大学, 2015. LI J.Numerical simulation of ice accretion and study for icing scaling law on aero-engine entry components[D].Xi'an: Northwestern Polytechnical University, 2015(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10699-1016909385.htm |
| [29] | 张丽芬, 张美华, 吴丁毅, 等. 旋转帽罩结冰相似准则的研究[J]. 推进技术, 2015, 36(8): 1164-1169. ZHANG L F, ZHANG M H, WU D Y, et al. Research on icing scaling law for rotating cone[J]. Journal of Propulsion Technology, 2015, 36(8): 1164-1169. (in Chinese) |
| [30] | BILANIN, A J. Proposed modifications to ice accretion/icing scaling theory[J]. Journal of Aircraft, 1991, 28(6): 353-359. DOI:10.2514/3.46034 |
| [31] | 杨世铭, 陶文铨. 传热学[M]. 北京: 高等教育出版社, 2006. YANG S M, TAO W Q. Heat transfer[M]. Beijing: Higher Education Press, 2006. (in Chinese) |
| [32] | 孔维梁.飞机异常结冰的过冷凝固机理及理论研究[D].上海: 上海交通大学, 2015. KONG W L.Mechanism and fundamental theory of supercooled water solidification on the abnormal aircraft icing[D].Shanghai: Shanghai Jiao Tong University, 2015(in Chinese). http://cdmd.cnki.com.cn/Article/CDMD-10248-1016787495.htm |
| [33] | WRIGHT W B, RUTKOWSKI A.Validation results for LEWICE 2.0: NASA/CR-1999-208690[R].Washington, D.C.: NASA, 1999. |
