0 引言

积冰是生活中常常遇到且涉及安全的一个问题，尤其是一些生活生产设备，当遇到积冰时会降低效率甚至发生安全事故。例如，飞机飞行过程中常常会遇到低温大气环境，这种低温大气环境往往会使得飞机的一些关键部位(例如发动机短舱、机翼前缘等)结冰，这些关键部位一旦结冰，会导致飞机的阻力增加、升力减小、操控性能下降以及燃油消耗增加，甚至引起发动机故障等情况，从而使得飞行器的飞行安全裕度大幅度降低，进而造成一些飞行安全事故。风力发电机在寒冷气候中会在叶片上产生积雪或者积冰，当积雪或者积冰量达到一定程度后，风力发电机的效率会大大降低甚至造成风力发电机损毁，其叶片上甩脱的冰块会对周边风力发电机造成损伤，甚至对人造成伤害^[1]。据调查，每年因积冰问题引发的飞机安全事故以及风力机效率降低甚至损毁就有近几十起，所以，积冰是一个重要的安全问题，急需加以解决。为了解决积冰问题，需要对积冰的原理、防冰及除冰方法进行研究。在当前风力机和飞机除冰的相关研究中，冰的融化过程是非常重要的一方面。

冰的融化过程比较复杂，涉及到多个问题，如：积冰密度不同导致融化过程不同、融化区域的对流换热问题、冰内部融化过程、冰外部融化过程、冰脱落过程、水冻结成冰的过程、冰融化过程中的传热特性、冰孔隙率问题、固-液交界面问题等等。针对这些问题，国内外学者开展了不少相关的研究工作。

关于积冰密度对热物性的影响，Aguirre varela等人^[2]研究了积冰密度对导热系数的影响，研究了水滴直径对努森数Nu的影响。

关于融化区域内的对流扩散问题，Voller等人^[3]使用基于焓法的固定网格方法对相变问题中混合区域的对流扩散进行了数值模拟研究，这种方法的基本特征在于通过选择适当的源项后能够表现潜热及固-液混合区域的演变过程。他们^[4]还开发了一种用于求解相变过程中对流/扩散问题的焓法，这种方法的基本特征是潜热变化的影响与源项相互隔离，这种方法可适用于常规的对流/扩散相变问题，例如在等温或者一定温度范围内，可以有效地求解潜热变化。

关于冰块浸入水中的融化过程，Scanlon等人^[5]对圆柱形冰柱在水中的融化过程进行了试验及数值研究，使用PIV系统观察冰柱的融化过程，并且得到了冰柱周围的水流方向，最后通过数值方法对冰柱的融化过程进行了研究，并和试验值进行了对比，对比结果显示从定性方面PIV测试结果和CFD计算结果具有良好的一致性。

关于球内的融合过程，Tan等人^[6]对圆球内部相变材料的融化过程进行了数值和试验研究，基于可迭代的有限容积法对融化过程进行数值模拟，并且使用单域焓法模拟相变过程，最后对试验与数值的结果进行了分析对比。

关于冰颗粒的融化过程，Hauk等人^[7]对不规则形状冰粒的融化特性进行了试验研究，并使用球形冰粒融化模型对冰粒的融化时间进行了理论计算，试验值和理论计算时间吻合良好。

关于冰脱落过程的传热特性问题，肖春华^[8]等人基于热-焓值方法研究了不同加热模式和加热单元间隔情况下冰融化脱落对电加热除冰过程中传热特性的影响。

关于水冻结成冰的过程，王凯旋等人^[9]对有限厚度平板的冻冰过程中固-液相变交界面进行了研究，采用精确解与积分解相结合的方法，分析相变过程中的固-液交界面、边界热流密度以及固、液能量随时间的变化规律。

关于冰水相变过程中的传热特性问题，管建春等人^[10]对圆球形相变储冷器的相变传热过程进行了研究，采用焓法对水在球内的凝结相变过程进行了数值模拟，并且得到了冰水两相交界面位置与球的直径及温差的关系。张敏等人^[11]使用移动热源法对冰水相变传热过程进行了数值模拟研究，通过导热方程源项的线性离散和数值计算，得到了冰水区域内的温度分布。鄢新^[12]对冰水相变传热过程进行了数值仿真模拟，研究了冰水相变的动态过程、温度场对冰融化速率的影响、冰水相变界面传热问题等，并且设计了冰水相变的试验装置进行了试验研究。雷桂林等人^[13-14]基于融化模型对电热除冰系统中冰的相变过程进行了动态数值模拟，研究了冰相变过程中的传热特性以及融化区域厚度随着加热时间的变化规律。

关于把冰看作为多孔介质进行相变研究，黄晓明^[15]发展了一种新的多孔介质相变传热模型，他首先分析了包含湿润的多孔介质内部湿润部分和热量的传递机理，然后使用连续的介质力学以及局部体积平均技术，发展了一种能够比较全面的处理包含湿饱和多孔介质内部传热问题的模型；朱杰^[16]研究了多孔介质内湿饱和流体的凝固相变过程及其传热特性，使用Whitaker提出的体积平均理论将多孔介质近似的看做一种拟连续介质，并且使用了焓法求解流体凝固的相变过程。崔海亭^[17]等人研究了填充相变材料的泡沫金属的相变过程，使用FLUENT软件的凝固/熔化模型对相变材料的相变过程进行了数值模拟研究，分析了材料的温度场及相界面变动规律。

对于固-液交界面问题的研究，Du^[18]等人在方形的封闭腔中对乙醇胺和水的二元混合物的融化特性进行了试验研究，追踪了液-固交界面并测量了瞬态的液相分数、研究了自然对流对熔化分数及固-液交界面形状的影响。

冰的物理特性对融化过程影响十分明显。而积冰具有不同的类型，例如霜冰、明冰及混合冰等等。不同类型的积冰主要是由不同的气候条件造成的。影响结冰类型的气候参数主要有：温度、风速、液态水含量以及水滴大小等, 这些参数的变化给冰融化过程带来的影响有些会相互抵消，即来流条件与冰融化过程并非一一映射关系，因此需要寻找一个更简洁的参数与冰融化过程进行关联。在不同来流条件下形成的霜冰、明冰以及混合冰的最主要区别是其密度不同，因此可以研究积冰密度与冰融化过程的关系，从而得出积冰密度与除冰过程的关系。目前，关于密度对冰融化过程影响的相关研究较少。本文针对机翼上不同类型的积冰，采用数值方法研究了不同冰密度对于电加热除冰过程中热传导及冰融化过程的影响。

1 融化模型 1.1 模型假设

融化模型把积冰看作是含有不同孔隙率的多孔介质。其中，介质本体为积冰，作为固相；积冰中含有大小一致均匀分布的孔隙，内部充满空气，作为气相。孔隙总体积占积冰层总体积的比值为孔隙率。融化相变过程中只考虑积冰-水的相变过程，忽略孔隙中空气的相变。在电加热除冰过程中，使用融化模型来求解积冰层的融化相变过程，得出不同时刻积冰层的温度分布以及不同相的交界面。实际的冰-水相变物理过程非常复杂，为了能够在突出冰-水相变规律的同时有利于数值模拟的实施，采用了如下假设条件：1)材料的物性参数为常数；2)来流条件稳定不变，机翼蒙皮表面积冰融化过程中积冰的外形保持不变；3)电加热除冰模型中各层材料之间接触良好，不会阻碍传热；4)积冰层为纯净物，除孔隙中的空气外不含有任何杂质；5)相变过程发生在一个很小的温度范围内。

1.2 焓值-多孔介质法

采用的融化模型中使用焓值-多孔介质法替代常规的显式追踪法来确定相变过程中不同相的交界面。焓值-多孔介质法最早由Voller提出^{[3, 19-21]}并用于固化/融化过程的模拟，其中，将控制体内液态相所占体积与总控制体体积的比值定义为液相分数，用β表示。在迭代计算过程中，每次迭代都会基于焓值守恒原理对液相分数进行计算。如果液相分数值在0~1之间，则此区域被称之为混合区域。在全部融化的区域其液相分数为1；反之，在全部固化的区域其液相分数为0。对于积冰，当冰的温度值低于融点时为固态，其液相分数为0；当冰的温度值高于融点时，冰融化成水为液体，其液相分数为1；当冰的温度处于融点时，为冰水混合状态，其液相分数介于0~1之间。用β表示冰中的液相分数，其表达式如下：

$ \beta = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 0,\\ 0 \sim 1,\\ 1, \end{array}&\begin{array}{l} T < {T_{{\rm{melt}}}}\\ T = {T_{{\rm{melt}}}}\\ T > {T_{{\rm{melt}}}} \end{array} \end{array}} \right. $

(1)

式中，T_melt为融化温度。

1.3 控制方程

电加热除冰模型中，含内热源的非稳态导热方程为：

$ {\rho _j}{C_{pj}}\frac{{\partial {T_j}}}{{\partial t}} = {k_{xj}}\frac{{{\partial ^2}{T_j}}}{{\partial {x^2}}} + {k_{yj}}\frac{{{\partial ^2}{T_j}}}{{\partial {y^2}}} + {k_{zj}}\frac{{{\partial ^2}{T_j}}}{{\partial {z^2}}} + {Q_j} $

(2)

式中j表示电加热除冰模型中不同的层。

在冰层中，能量控制方程可表示为：

$ \frac{{\partial {H_{{\rm{ice}}}}}}{{\partial t}} = \frac{\partial }{{\partial x}}\left( {{k_x}\frac{{\partial T}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {{k_y}\frac{{\partial T}}{{\partial y}}} \right) + \frac{\partial }{{\partial z}}\left( {{k_z}\frac{{\partial T}}{{\partial z}}} \right) $

(3)

式中H_ice为冰层或者融化水的焓值，其温度与焓值的关系如下：

$ {H_s} = {C_{ps}}T,\;\;\;\;T < {T_{{\rm{melt}}}} $

(4)

$ {H_l} = {C_{pl}}\left( {T - {T_{{\rm{melt}}}}} \right) + {C_{ps}}{T_{{\rm{melt}}}} + {L_f} $

(5)

式中，H_s为冰层融化前的焓值，H_l为冰层融化成水后的焓值，C_ps为冰层的热容，C_pl为液态水的热容，L_f为冰融化成水所需要的潜热。

对式(4)、式(5)进行变换，得到温度的表达式：

$ T = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} H/{C_{ps}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;,\\ {T_{{\rm{melt}}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;,\\ {T_{{\rm{melt}}}} + \left( {H - {H_{lm}}} \right)/{C_{pl}}, \end{array}&\begin{array}{l} H \le {H_{sm}}\\ {H_{sm}} < H < {H_{lm}}\\ H \ge {H_{lm}} \end{array} \end{array}} \right. $

(6)

式中，H_sm、H_lm分别为固体材料与液体材料在温度值为融点时的焓值，其表达式为：

$ {H_{sm}} = {C_{ps}}{T_{{\rm{melt}}}} $

(7)

$ {H_{lm}} = {C_{ps}}{T_{{\rm{melt}}}} + {L_f} $

(8)

相变过程中材料的温度与焓值的关系如图 1所示。

材料的焓值可以由材料的显焓h与潜热ΔH之和表示：

$ H = h + \Delta H $

(9)

式中，

$ h = {h_{{\rm{ref}}}} + \int_{{T_{{\rm{ref}}}}}^T {{C_p}{\rm{d}}T} $

(10)

式(10)中，h_ref为参考焓值，T_ref为参考温度。

确定了液相比后，就可以把式(9)中的潜热表示成材料潜热L的一部分：

$ \Delta H = \beta L $

(11)

式中L表示材料的潜热。由式(11)可以看出，潜热值的范围是从0(固体状态)到L(液体状态)。于是，对于融化问题，能量方程可以写成以下形式：

$ \frac{\partial }{{\partial t}}\left( {\rho H} \right) + \nabla \cdot \left( {\rho \mathit{\boldsymbol{v}}H} \right) = \nabla \cdot \left( {k\nabla T} \right) + S $

(12)

式中，H为总焓，ρ为密度，v为流体速度，S为源项。

2 积冰参数 2.1 积冰密度

飞机或发动机上产生的积冰多为霜冰、明冰或者混合冰，主要的影响因素为结冰环境。根据Jones等人的研究结果^[22-27]可知，积冰的密度与来流中水滴直径、速度、液态水含量以及机翼前缘直径及其温度的关系式为：

$ \begin{array}{l} \rho = 1000\exp \left[ { - 0.15\left( {1 + \frac{{6043}}{{{S^{2.65}}}}} \right)} \right]\\ S = \frac{{MV{D^{0.82}}{V^{0.59}}LW{C^{0.21}}}}{{{D^{0.48}}{{\left( { - {T_{{\rm{wall}}}}} \right)}^{0.23}}}} \end{array} $

(13)

式中，MVD为水滴平均直径，单位为μm；V为来流速度，单位为m/s；LWC为液态水含量，单位为g/m³；D为机翼前缘直径，单位为cm；T_wall为机翼表面温度。由式(13)可知，当S≤6时，ρ接近于0；当7≤S≤60时，ρ随着S的增大呈非线性增大；当S>60时，ρ增长缓慢；当S趋于无穷大时，ρ≈860.7。由此可见，对于密度大于860.7 kg/m³的积冰，式(13)并不适用。

当积冰为霜冰或者混合冰时，其密度可采用式(13)估算；当积冰为明冰时，积冰密度达到最大值，约为920.0 kg/m³。

2.2 冰层热物性

霜冰与明冰的主要差异是冰中所含气泡量不同。当积冰为明冰时，冰中气泡含量最少，一般可认为气泡含量为0。当积冰为霜冰或者混合冰时，冰中含有一定量的气泡，冰中气泡占据了一定的体积，气泡的总体积与冰块(包含气泡与冰骨架)的体积比即为冰块的孔隙率。通过上述关系推算，可知积冰的孔隙率可以通过空气密度、积冰(霜冰或者混合冰)密度及明冰的密度表示，关系式为：

$ \gamma = \left( {\frac{{{\rho _{\max }} - \rho }}{{{\rho _{\max }} - {\rho _{{\rm{air}}}}}}} \right) \times 100\% $

(14)

式中，γ为积冰的孔隙率；ρ_air为空气密度，取值为0 ℃时的密度值1.293 kg/m³；ρ_max为明冰的密度，取值为920 kg/m³；ρ为霜冰或者混合冰的密度。

积冰被当作多孔介质时，其内部孔隙中的填充介质为空气。假设通过积冰、空气的热传导同时发生，而在积冰与空气的交界面上不发生热量交换，且满足平行热传导模型^[28]，则多孔介质积冰的当量导热系数可以通过如下公式得出：

$ {k_{{\rm{eff}}}} = \gamma {k_{{\rm{air}}}} + \left( {1 - \gamma } \right){k_{{\rm{ice}}}} $

(15)

式中，k_eff为多孔介质积冰的当量导热系数，k_air为空气的导热系数，k_ice为冰处于最大密度时的导热系数。

多孔介质积冰的当量热容可以通过如下公式得出：

$ {C_{p{\rm{eff}}}} = \gamma {C_{p{\rm{air}}}} + \left( {1 - \gamma } \right){C_{p{\rm{ice}}}} $

(16)

式中C_peff为多孔介质积冰的当量定压热容；C_pair为空气的定压热容；C_pice为冰处于最大密度时的定压热容。

为了研究密度对积冰融化相变过程的影响，选取了典型结冰气候条件下形成的霜冰以及积冰具有最大密度时的明冰作为电加热除冰系统的冰层，在相同的加热条件下进行除冰分析，比较不同密度情况下的电加热除冰过程中积冰层的温度变化及融化相变过程。在进行除冰计算时，根据常见的结冰条件确定来流速度，从而确定霜冰的密度。除冰过程中计算的时间较短，冰层的融化的厚度较小，还没有发生明显的位置或者形状变化，且冰层内部环境仍然保持相对不变，因此忽略了冰层形状变化对外部流体流动的影响。霜冰的物性参数由式(13)根据常见结冰条件确定。取两个不同的来流速度，分别为60 m/s和180 m/s，来流温度均为-12.2 ℃，平均水滴直径均为20 μm，液态水含量均为1 g/m³，机翼前缘直径为5 cm。由式(13)可得出积冰的密度分别为794.8 kg/m³和848.5 kg/m³。积冰的相关参数如表 1所示。

3 物理模型

计算分析中采用的三维电热除冰分析模型如图 2所示。模型共分为六层，由下至上分别是机翼结构、内部绝热层、电加热片、外部绝热层、机翼蒙皮及冰层。此模型与实际机翼的铺层相似，采用此模型进行电热除冰分析能取得与真实情况较为接近的结果。

图 2所示计算分析模型的长×宽×高尺寸为6.35 mm ×6.35 mm×10.489 mm。各层材料的厚度及物性参数如表 2所示，表 2中d为厚度，k为导热系数，C_p为热容。冰层初始厚度为6.35 mm，其导热系数及热容根据表 1确定，随着积冰融化过程的进行，冰层内部会逐渐产生一层液态水，其厚度随着积冰融化过程发生变化，其导热系数为0.5538 W/(m·K)，热容为4186.553 J/(kg·K)。

分析模型中，电加热片在长宽方向的尺寸为模型的一半，热流密度为46 500 W/m²。整个模型的初始温度与环境温度相同，均为-12.2 ℃。模型的四个垂直侧面均为中间对称面，底部平面为机翼内部表面，顶部平面为积冰层与来流空气的接触面，故在数值计算时，四个侧面均设置为绝热壁面。在除冰过程中，只有机翼蒙皮表面很薄的一层冰产生了融化现象，此时热源的热量对积冰外表面温度影响很小，因此假定积冰外表面温度与来流空气的温度相同。底部平面与内部空气之间的对流换热系数为常值5.68 W/(m²·K)。

为了分析电热除冰效果，取冰层与蒙皮交界面作为除冰效果的监测面，如图 3所示。图中线Aa为加热片的中心线，垂直于电加热片并过中心；线Bb及线Dd位于两个相邻电加热片的中间位置，与交界面垂直；线Cc距离电加热片最远，同样垂直于交界面。图 2所示的三维电热除冰模型计算域在除冰系统中的相对位置如图 4所示。图 3所示各点与电加热片的相对位置如图 4所示。由图 4可知，在除冰系统中电加热片等间距布置，计算域在长宽方向上取了电加热片及中间隔热材料的各一半。加热片长宽尺寸及中间隔热材料尺寸均为6.35 mm。由图中各点与电加热片的距离可知在除冰过程中，A点位置的积冰最先开始融化，B点及D点位置附近的积冰随后同时开始融化，C点位置附近的积冰最晚开始融化。因此，C点是判断交界面积冰是否完全融化的关键位置，只有C点位置附近的积冰全部融化之后，才能基本判断交界面上的积冰已经完全融化。

4 结果与讨论

通过不同来流条件得到不同的积冰密度，在电加热除冰模型中研究了在同样的加热条件下积冰密度对除冰效果的影响。主要研究了密度对积冰温度分布、初始融化时间以及除冰时间的影响。针对表 1中的霜冰和明冰开展了电加热除冰及冰融化相变过程的研究。

4.1 算法验证

为了验证数学物理模型计算结果的可靠性，将积冰为明冰时除冰过程中A点位置的温度计算结果与Alan D. Yaslik的研究结果^[29]进行了对比，如图 5所示。结果显示，所用的计算方法及结果是可靠的。

4.2 密度对温度的影响

交界面的温度分布情况决定了积冰的融化与脱落情况。图 6显示了不同密度的积冰在同样的电加热条件下在不同时刻直线AC处的温度曲线。由图可知，A点附近的温度最高，C点附近的温度最低；在同样的电热功率情况下，积冰密度低的情况直线AC上的温度高于积冰密度高的情况；直线AC上发生融化的区域随着加热时间的增长由A点逐渐向C点移动。

4.3 密度对积冰初始融化时间的影响

由图 3及图 6可知，在电加热除冰系统工作时，冰层与蒙皮交界面上A点位置的温度最高，即A点位置的积冰最先开始融化，因此可以把A点位置温度大于0 ℃的时刻看作为积冰初始融化的时刻。在同样的电加热及结冰条件下，A点位置的瞬态温度曲线如图 7所示。图中黑色实线表示积冰为明冰时A点位置的温度随加热时间的变化情况，红色线段和蓝色点画线表示积冰为不同密度的霜冰时A点位置的温度随加热时间的变化情况。由图可知，当积冰密度为920 kg/m³时，电加热除冰系统工作2.7 s时，A点位置的积冰开始融化，加热至4.6 s时，A点位置的积冰完全融化，这个融化过程持续约1.9 s；当积冰密度为848.5 kg/m³时，电加热除冰系统工作2.5 s时，A点位置的积冰开始融化，加热至4.2 s时，A点位置的积冰完全融化，整个融化过程持续约1.7 s；当积冰密度为794.8 kg/m³时，电加热除冰系统工作2.4 s时，A点位置的积冰开始融化，加热至3.9 s时，A点位置的积冰完全融化，整个融化过程持续约1.5 s。由此可知，当积冰的密度降低时，在同样的加热条件下，A点位置温度升高的速度会加快，初始融化的时间会提前，整个融化过程会缩短。

4.4 密度对除冰时间的影响

由图 3及图 6可知，在电加热除冰系统工作时，冰层与蒙皮交界面上C点位置的温度最低，即C点是判断交界面上积冰是否完全融合的关键位置。因此，假设C点位置的温度大于0 ℃时，交界面上的积冰已经完全融合。C点位置的温度随时间变化的曲线如图 8所示。图中黑色实线表示积冰为明冰时C点位置的温度随加热时间的变化情况，红色线段和蓝色点画线表示积冰为不同密度的霜冰时C点位置的温度随加热时间的变化情况。由图可知，当积冰密度为920.0 kg/m³时，即为明冰时，C点位置于除冰系统工作42.5 s后开始融化，并于57.5 s完全融化，整个融化过程持续大约15 s；当积冰密度为848.5 kg/m³时，即为霜冰时，C点位置于除冰系统工作39 s后开始融化，并于53 s完全融化，整个融化过程持续大约14 s；当积冰密度为794.8 kg/m³时，即为霜冰时，C点位置于除冰系统工作34.5 s后开始融化，并于47.5 s完全融化，整个融化过程持续大约13 s。由此可知，当积冰密度降低时，在同样的加热条件下，C点位置的温度会升高，融化起始时刻会提前，整个融化周期会变短，即冰层与蒙皮交界面处的积冰完全融化的时间会缩短。

5 结论

针对不同来流条件下形成的不同密度的积冰，采用数值方法研究了积冰密度对于电加热除冰过程中冰融化过程及冰层与蒙皮交界面温度的影响。选取了两个典型来流条件下的霜冰与明冰共三种不同类型的积冰作为分析对象，通过比较霜冰与明冰在融化过程中交界面的温度变化情况，分析了密度对冰融化过程的影响。计算结果表明，冰密度对电热除冰过程具有明显的影响，密度降低会提高积冰温度升高的速度，使得冰初始融化的时刻提前，整个融化过程缩短。

致谢: 本文的研究得到了上海交通大学机械与动力工程学院工程热物理研究所陈勇老师的支持和帮助，在此表示感谢。