不同浸润性冷表面上水滴碰撞结冰的数值模拟

冷梦尧，常士楠，丁　亮



研究简报

不同浸润性冷表面上水滴碰撞结冰的数值模拟

冷梦尧，常士楠，丁亮

（ 北京航空航天大学航空科学与工程学院，北京 100191）

摘要:对冷水滴撞击不同表面时的动力学行为和相变过程进行了模拟。通过耦合 VOF和Level-set方法追踪气液自由界面，结合焓-孔隙度相变模型，模拟水滴撞击冷表面的动力学行为及相变特征。选取亲水(接触角30°)、疏水(接触角114°)和超疏水(接触角163°) 3种典型浸润性的表面，计算了多种壁温条件下的水滴撞击结冰过程。结果表明提高表面疏水性，将减小水滴与冷表面的接触时间和接触面积，降低水滴内的相变速率，延缓水滴结冰的时间。在表面温度高于-15℃时，超疏水表面可以避免冷水滴的冻结黏附，保持表面洁净。将模拟得到的最大铺展直径、回缩速率以及冻结情况，与已有实验结果进行对比验证，表明了模拟方法的有效性和准确性。

关键词:表面；液滴碰撞; 动力学; 多相流；相变；数值模拟

引 言

水滴撞击冷表面发生结冰的现象，属于气-液-固三相流的流动传热问题，包含水滴撞击后的动力学行为以及水的凝固相变过程；而这二者之间的相互影响，使得水滴碰撞表面结冰往往会导致严重的后果且难以准确预测[1-2]。近几年研究者们通过实验发现，提高表面疏水性能够获得减少结冰量[3]，降低防冰能耗[4]，削弱冰的黏附力[5]等效果。然而，传统的结冰理论模型[6-7]主要从水滴撞击时的各项热流出发，通过热力学平衡判断水滴发生冻结与否，无法解释表面特性对水滴碰撞后的相变行为的影响。其原因在于计算时假设相变过程的时间极短而将其忽略，缺少对水滴撞击表面时其动力学行为与相变过程相互影响的考虑。

关于过冷水滴静止状态下的结冰过程已有相关文献研究[8-11]。实验表明，水滴的完全冻结过程平均长达几十秒[12]；在超疏水表面上，水滴冻结的时间会进一步被延长[13]，而碰撞过程的时间则要短得多。对于冷水滴碰撞的动态过程，Mishchenko等[14]和 Li等[15]的研究均揭示出超疏水表面可以显著改变结冰情况，有效防止水滴的冻结黏附。但是受实验手段的限制，运动水滴内部的相变详情难以体现。作为实验观测的补充手段，研究者们通过数值模拟研究水滴动力学行为，如VOF和Level-set相结合的方法，获得了与实验结果符合良好的结果[16-17]。数值模拟能够获取更为全面的水滴内部情况，通过分析水滴内的压力场[18]、速度矢量[19-20]和温度分布[21-22]等，能够进一步解释水滴动力学行为和传热机理。

本工作的目标是建立能够准确地同时追踪气液界面和揭示液固相变的数值方法，获得表面浸润性对水滴碰撞冻结的影响。针对水滴碰撞不同浸润性的冷表面，将多相流的界面追踪方法（耦合Level-set和VOF方法）和相变模型（焓-孔隙度模型）相结合，模拟了冷水滴撞击到不同表面上的铺展-回缩行为，与实验得到的湿润因子对比验证本方法的准确性，并进一步揭示水滴内部相变区域的发展过程以及相变与回缩行为的相互作用，研究了表面特性对水滴内部相变动态过程的影响，以期为表面结冰的界面现象和防冰/除冰机理研究提供参考。

1　数学模型

1.1 自由界面追踪模型

本文采用VOF耦合Level-set方法[23]追踪气液自由界面，流场求解以 VOF的思想为主，在考虑界面效应时，以Level-set方法为辅助。VOF方法能够保持追踪流体的“质量”，并且能够容易地追踪拓扑结构变化的流体界面。在 VOF方法中，相与相之间的界面追踪是通过计算次相的体积分数连续方程完成的。由于本文的研究对象为两相流，主相为空气，次相为水，故而只需要求解次相体积分数的连续方程

主相的体积分数则通过总体积分数之和为1计算得到，αP=1-αS。动量方程和能量方程中的物性参数（如密度、黏度等）则通过各相体积分数进行加权计算。

VOF方法通过几何重构法对界面进行插值构造，即基于部分填充单元内的体积分数及其梯度，用分段线性函数构造相界面；然后根据通过界面的流量，更新单元内体积分数。但是分段线性界面不符合光滑和连续的要求，无法提供准确的界面方向和曲率，故而在此引入 Level-set方法来解决这个问题。

Level-set方法是将气液相界面用一个高阶函数，即Level-set函数ϕ的零值点来表示，由ϕ的代数值来区分计算区域中的各相。Level-set函数ϕ为一距离函数，其定义如下

但是由于数值计算过程的内在效应，Level-set函数在经过几个时层的计算后，将不再是距离函数，需要进行距离函数的初始化运算，使得Level-set函数在不改变其0等值面位置的前提下重新构成为距离函数。另外，Level-set方法的数值求解过程中，由于数值耗散和数值误差等原因可能造成封闭界面内流体质量不守恒[24]。

耦合VOF和Level-set的自由界面追踪方法将这两者各自的优势结合起来:首先，在构造相界面时，通过相体积分数的值确定单元内分界面位置，依靠Level-set函数ϕ的梯度确定界面法向方向；然后采用函数ϕ计算相界面的曲率及其相关量，提高相界面参数的求解精度；该时间步结束时，通过相体积分数确定的界面位置更新流场中的Level-set函数，以克服Level-set方法求解对流输运方程时的质量不守恒问题。

1.2 流动计算

表面张力、壁面浸润性以及相变过程对两相流流动的影响，通过对基本动量方程添加源项来实现。式(4)为流场内动量方程的表达式

等式左边依次为非稳态项以及对流项，右边前3项分别为压力项、黏性项和体积力项。等式右边的-σκδ（ϕ）∇ϕ为基于连续表面力(CSF)模型添加的表面张力源项，以考虑表面张力σ和自由界面的曲率κ对流动的影响。SF为凝固引起的流动阻力。

壁面处控制体内的自由界面法向可通过定义壁面接触角来确定，即通过修正壁面附近的界面法向向量n，令其等于式(5)从而改变界面曲率及表面张力源项。式(5)中nw和tw分别表示壁面的法向和切向单位向量，而θw为定义的壁面接触角

1.3 凝固-融化模型

实验表明，水在过冷状态下结冰与水从常温下逐渐冷却结冰的过程不同，其包含两个过程[25]:①快速结冰过程，当受到某种触发作用（如降温、振动、碰撞等）导致冰核成形，接着冰枝迅速生长并充满水滴内部，水滴相变部分变成疏松的海绵状结构，并升温至冰-水混合平衡温度；②缓慢冻结过程，若冷表面或冷空气继续从水滴中抽走热量，则水滴会逐渐提高固体的比例，这一过程是等温变化，过程的长短由传热速率决定。

根据实验观测的结果，采用一种焓-孔隙度相变模型来模拟水滴内部的凝固-融化相变过程[26]，计算求解的能量方程为

考虑相变后，能量方程中的焓值为显热及潜热之和，H=cpT+ΔH。其中潜热部分与水的凝固/融化相变潜热的关系为

焓-孔隙度相变模型将冰水混合区域视为多孔结构，孔隙率即为液相的体积分数γ，取值区间为[0,1]，其与温度的关系为

假设冰-水混合区内的运动遵守多孔介质内流动的Darcy定律，且符合Carman-Koseny假设，则可推出动量方程源项为[26]

其中，ε为极小值（ε=0.001），以避免分母为零；Cmush为黏糊系数，与多孔介质的形态相关，合理的Cmush常数取值，既可以在液相比例较高时保持流动性，又可以在液相比例降低后抑制流动速率。

从动量源项的表达式可以看出，在液相区，即γ=1，动量源项为0，动量方程与流体N-S方程形式一致；凝固开始变成冰-水混合物后，液相体积分数γ降低，源项在动量方程中的比重逐渐超过瞬态项、对流项和扩散项，直至完全转变为固相后，速度基本为零。

图1　计算区域Fig.1 Droplet impingement case setup

2　计算模型

通过模拟冷水滴撞击不同表面上的动态相变过程，验证本方法的可行性。本文的模拟区域为高4倍水滴直径、宽3倍水滴直径的二维轴旋转对称区域，计算模型及边界条件如图1所示，模型底部为定温无滑移壁面边界，顶部及侧边界均为大气边界条件。对网格进行无关性检验后，决定以沿初始直径分布200个网格为准进行划分，最终计算区域由480000个正方形结构化网格组成，能够准确地捕捉到气液自由界面以及相界面的移动变化。计算中压力速度耦合采用适合非定常流动问题的PISO方法，气液自由界面重构使用Geo-Reconstruct方法。除梯度和压力分别使用Least Squares Cell Based和PRESTO！格式进行离散外，剩余各项均使用二阶迎风格式。为节省计算资源并保证迭代结果的准确性，计算过程采用变时间步长，控制Courant数始终小于0.25，时间离散格式为一阶隐式格式，每个时间步内迭代20次。

本文使用FLUENT 14.5软件进行模拟，计算工质为空气和水，空气相为主相，水为第2相。为了同实验结果进行对比，选择直径3 mm的水滴，以初速度1.4 m·s-1自1倍直径高度下落撞击固体壁面，初始时刻，空气温度为-5℃，水滴温度为相平衡温度。计算了3种典型浸润性表面，其接触角分别为30°、114°和163°，壁面温度为20、-15、-25℃共9种碰撞情形。计算过程中使用的水滴物性参数如表1所示。

表1　水滴物性参数Table 1 Material properties for water

3　结果与讨论

3.1 表面温度对水滴碰撞动力学行为的影响

图2给出了3种典型浸润性表面，在不同温度条件下水滴碰撞后的形态演化过程。在Blake等[16]对纯水结冰所做相关研究的基础上，本文计算中使用的黏糊系数 Cmush如表 2所示，从中可以看出当表面温度越低或浸润性越强，即条件越有利于加快凝固时，则Cmush越大，这与Cmush象征的物理意义相符合。

图2　不同温度及不同浸润性表面上液滴碰撞后的形态变化过程Fig.2 Dynamic behavior of droplet impinging on surfaces of different wettability and temperature

首先对比常温（20℃）条件无相变发生时，不同浸润性表面上水滴的铺展和回缩过程，结果表现出不同浸润性表面上液滴碰撞的典型动态特征。在亲水性表面（θ=30°）上，水滴可达到最大的铺展面积，且回缩过程缓慢，最后水滴收缩至圆盘状并停留在表面上。增大表面接触角后，从水滴的边缘形状可以看出，水滴的铺展润湿行为受阻，不仅最大铺展直径减小，达到最大铺展的时间，即回缩过程开始的时间，也相应提前。在疏水（θ=114°）及超疏水（θ=163°）表面上，水滴同表面之间摩擦损耗以及水滴内部的黏性耗散减少，因此，水滴剩余的动能较大，其回缩过程加快，最终形成完全反弹离开表面。

表3　黏糊系数Cmush的取值Table 2 Value of mushy zone constant Cmush

对比相同浸润性、不同温度条件下水滴的动力学行为发现，水滴前半段的铺展过程都十分相似，且在相同浸润性的表面上，水滴最大铺展形状也十分相似。这说明水滴的撞击热以及液体流动时产生的黏性摩擦热能够一定程度上延缓相变的发生。由于在水滴达到最大铺展直径前后，也即是水滴动能最小的阶段，水滴与表面之间换热面积最大，且黏性及摩擦产生的加热热流最低，故而形成了最有利于结冰发生的条件。从水滴形态对比图（图 2）中可以看出，在常温表面和远低于冰点的表面上，水滴碰撞后半段的回缩过程有明显不同的动态行为。在低于冰点的表面上，由于相变是从水滴底部开始逐渐向上发展，故可以看到水滴回缩时分成了上下两层流动，下层有明显的黏滞感甚至停止了运动，而上层仍保持有液体的流动性。

从亲水性表面（θ=30°）的水滴形态对比图中可以看到，在壁面温度为-25℃时，水滴撞击铺展后几乎没有回缩过程，仅表层薄水膜有轻微振荡，最后水滴以最大铺展的形状冻结到了表面上，这一结果也同文献[14]中实验观察到的现象一致。同样降至-25℃时，水滴碰撞结果由反弹变为停留在表面上。比较疏水表面（θ=114°）上的两个回缩过程发现，常温壁面上，水滴收缩时呈现外缘凸起中心凹陷的形态，来自三相交界线处的表面张力是促进水滴回缩的主要动力；而在-25℃的表面上，由于水滴底部被冻住，表面张力变成了阻止水滴弹离表面的作用力，水滴上层的液体部分仅依靠惯性回缩至中心区域，产生弹起的趋势，随后被拉回表面上。若进一步提高表面的疏水性，在超疏水（θ=163°）表面上，由于最大铺展直径减小，水滴铺展过程进一步缩短，摩擦及黏性造成的损耗更小。水滴在冻结之前保有一定动能，能够产生部分回缩，缩小冻结面积，且最终水滴的上层液体有足够的惯性脱离冻结的底部，实现部分反弹。

此外，从不同浸润性表面上水滴最大铺展发生的时刻可以看出，接触角越大，最大铺展发生的时刻越靠前，表明水滴同表面之间换热的概率也越低。总的来说增强表面疏水性对于延缓或防止水滴冻结有两个方面的优势:一方面是通过阻止水滴润湿表面，减小其铺展面积，降低了散热速率；另一方面是减少水滴碰撞过程中的动能损耗，加快水滴回缩过程，缩短水滴与表面接触的时间。故而在表面温度低于冰点时，超疏水表面更能保持表面的洁净。

3.2 表面浸润性对水滴湿润系数的影响

水滴与壁面之间的接触面积决定了换热量的大小，故而监测水滴碰撞过程中润湿面积的变化能够反映表面浸润性及表面温度对水滴碰撞结冰过程的影响。首先，图3比较了不同浸润性和不同温度下能够达到的最大铺展系数，其定义为水滴湿润面积同初始时刻水滴在表面上的投影面积之比。由于相变过程是从水滴的最大铺展状态前后开始，所以在表面温度相同时，水滴的最大铺展系数越大，水滴冻结的可能性就越高。从图3中可以看出最大铺展系数主要决定于表面浸润性的大小，接触角相同的表面上，最大铺展系数的变化量在 5%以内。超疏水表面上最大铺展系数为亲水性表面上的 1/2，即水滴的最大换热面积可减少 1/2。此外，从图 3中还可看出，在亲水性表面（θ=30°）上，随着壁面温度降低，最大铺展系数减小，表明水滴在完全铺展之前就冻结在了表面上。而在超疏水表面（θ=163°）上，水滴最大铺展系数随着壁面温度降低略有提高，表明在此情况下结冰发生在完全铺展之后，由于水滴底部的相变影响了三相接触线附近表面张力的作用，使得水滴收缩过程有所推迟。

图4和图5定量地比较了各个碰撞表面上水滴湿润面积的动态变化过程。其中，不同过程中的湿润面积以各自的最大湿润面积进行量纲1化，结果表明与文献[14]实验测得的量纲1湿润面积变化趋势一致，最终时刻的湿润面积比相符。从图4中可看到，水滴的铺展过程随着表面浸润性的提高而延长。在θ=30°的亲水性表面上，水滴持续铺展至10 ms才开始回缩，且回缩过程缓慢，而疏水表面和超疏水表面上达到最大铺展面积的时间依次为 6 ms和4 ms。当表面降低至-15℃后，亲水表面上水滴达到最大铺展的时间略提前，且回缩过程消失；同样地，疏水表面上的水滴最终产生部分冻结。但是超疏水表面上水滴的湿润面积同常温下的变化过程基本一致，这说明在-15℃时，超疏水表面能够实现防止水滴碰撞结冰。

图3　不同温度和浸润性壁面上水滴最大铺展系数的对比Fig.3 Comparison of maximum wetting factor with different temperatures and contact angles

图4　-15℃壁面温度下液滴润湿因子的变化Fig.4 Variation of wetting area ratio on different wettability surfaces at -15℃

图5　-25℃壁面温度下液滴润湿因子的变化Fig.5 Variation of wetting area ratio on different wettability surfaces at -25℃

当表面温度进一步降低至-25℃时，如图 5所示，水滴碰撞后的铺展过程与常温壁面下仍然保持一致，并且达到最大铺展的时间也基本不随表面温度的变化而改变。但是降低表面温度后，各个表面上的冻结面积均有所增加，在-25℃时，最终疏水表面上的冻结面积比为0.8左右，而超疏水表面上的冻结面积比大约为0.5。比较两者的回缩过程发现，超疏水表面能够延长回缩的时间，说明提高表面疏水性能够减少水滴同表面间的换热，且一定程度上使水滴保持了更高的动能以克服相变带来的流动阻力。

值得注意的是，疏水表面上数值模拟的水滴收缩过程均比实验结果偏快，这是由于本文研究中使用静态接触角模型表征表面浸润性，忽略了水滴滞后角的影响。实验中使用的疏水表面约有30°的滞后角，而最终模拟的结果表明，在滞后角较小（Δθ=5°）的超疏水表面上与实验结果符合得更好。

3.3 疏水性表面的减冰/防冰机理

同亲水性表面及弱疏水性表面相比，超疏水表面能够进一步减少水滴与壁面之间的接触面积，降低碰撞过程中的黏性耗散，从而防止或者延缓水滴冻结在表面上。但是水滴碰撞冷表面并完全反弹，并不意味着在此过程中水滴内部没有相变发生。分析水滴内的相变云图对探究表面特性对水滴碰撞后动态行为的影响有重要意义。以-15℃的温度下，3种不同浸润性表面上的水滴碰撞相变过程为例，模拟结果表明随着表面疏水性的提高，水滴碰撞后的结果逐渐由完全铺展变成部分反弹，进而是完全反弹，如图6～图8所示。

相变云图能够解释这种变化产生的原因。从图6中可以看到，在亲水表面上，相变发生在水滴达到最大铺展(10 ms)之前。由之前的研究可知，在水滴的边缘区域，由水滴表面张力及界面曲率引起的内部高压是促使水滴回缩的主要动力。但在亲水性表面上，由于接触角很小，在三相接触线附近没有足够的动力克服相变带来的阻力，最后水滴剩余的动能被黏性振荡消耗，直至静止。而在疏水表面上，如图 7所示，虽然相变也从较早的时刻t=5.0327 ms 发生，但是水滴边缘区域提供了较强的表面张力促进液体向中心区域回流。从12.992 ms 和27.835 ms的相变云图中能够看到，水滴底部虽然已经完全变成了固体，红色区域的液体体积分数为 0，但其边缘处被稍微拉离了壁面。若进一步提高表面疏水性，则水滴回缩的动力会更强，从图 8 中7.9915 ms时刻的相变云图可以看到，水滴回缩过程中就可以将底部中间的部分相变区域卷起。在水滴离开壁面后，仍未见有完全固化的区域，表明在超疏水表面上，水滴与壁面接触时间十分短暂，足以令水滴相变过程的时间长度不可忽略。由于同超疏水表面的接触时间很短，水滴底部发生的相变并不完全，加上表面张力的作用较强，最终水滴剩余的动能足以将发生相变的部位裹挟起，共同反弹离开表面。

图6　-15℃亲水表面上水滴碰撞结冰过程Fig.6 Dynamic behavior and phase transition of droplet impinging on hydrophilic surface at -15℃

图7　-15℃疏水表面上水滴碰撞结冰过程Fig.7 Dynamic behavior and phase transition of droplet impinging on hydrophobic surface at -15℃

图8　-15℃超疏水表面上水滴碰撞结冰过程Fig.8 Dynamic behavior and phase transition of droplet impinging on super-hydrophobic surface at -15℃

4　结 论

（1）基于耦合VOF和Level-set的方法，结合焓-孔隙度模型，模拟了冷水滴碰撞不同温度、不同浸润性的表面之后的结冰过程，准确捕捉到水滴碰撞过程中的动态行为，与实验结果符合较好。且本方法能反映水滴内部的相变细节，为研究表面浸润性对水滴碰撞结冰的影响提供依据。

（2）从水滴内部相变的发展过程分析不同浸润性表面上水滴冻结情况，结果表明提高表面疏水性不仅能够减小水滴碰撞后的铺展湿润面积，而且能够缩短水滴同壁面的接触时间，从而降低水滴内部的相变速率。

（3）在水滴碰撞的回缩过程中，疏水或超疏水表面能够产生较大的表面张力克服部分相变引起的黏滞阻力，从而起到减小冻结面积，降低结冰危害程度的效果。在高于某一温度的范围内(-15℃)，超疏水表面能使水滴完全弹离表面。

（4）为提高数值模拟的准确性，下一步工作可以关注在模型中考虑液体运动及相变引起的接触角和物性参数的变化。

符 号 说 明

Cmush——黏糊系数

cp——比热容，kJ·(kg·K)-1

d ——流场内控制体到自由界面的距离

H ——焓，kJ·K-1

k ——热扩散系数，kJ·(K·m·s)-1

L ——潜热，kJ·K-1

n,t ——分别为面积法向、面积切向

P ——压力，Pa

T ——温度，℃

Tsolid,Tliquid——分别为凝固温度、融化温度，℃

t ——时间，s

v ——速度，m·s-1

α ——体积分数

γ ——液相体积分数

θ ——接触角，(°)

κ ——曲率，m-1

μ ——黏度，kg·(m·s)-1

ρ ——密度，kg·m3

σ ——表面张力，N·m-1

下角标

P ——两相流中的主相

S ——两相流中的次相

w ——壁面

References

[1] LAMRAOUI F, FORTIN G, BENOIT R, et al. Atmospheric icing impact on wind turbine production[J]. Cold Regions Science and Technology, 2014, 100: 36-49. DOI: 10.1016/j.coldregions.2013. 12.008.

[2] 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.

[3] MORITA K, OKAMOTO K, AOKI A, et al. Hydrophobic coating study for anti-icing aircraft[R]. SAE Technical Paper, 2011.

[4] ANTONINI C, INNOCENTI M, HORN T, et al. Understanding the effect of superhydrophobic coatings on energy reduction in anti-icing systems[J]. Cold Regions Science and Technology, 2011, 67(1/2): 58-67. DOI: 10.1016/j.coldregions.2011.02.006.

[5] KIMURA S, YAMAGISHI Y, SAKABE A, et al. A new surface coating for prevention of icing on airfoils[R]. SAE Technical Paper, 2007.

[6] MESSINGER B L. Equilibrium temperature of an unheated icing surface as a function of air speed[J]. Journal of the Aeronautical Sciences, 1953, 20(1): 29-42. DOI:10.2514 /8.2520.

[7] BRAKEL T W, CHARPIN J P F, MYERS T G. One-dimensional ice growth due to incoming supercooled droplets impacting on a thin conducting substrate[J]. International Journal of Heat and Mass Transfer, 2007, 50(9/10): 1694-1705. DOI:10.1016/j.ijheatmasstransfer. 2006.10.014.

[8] JUNG S, DORRESTIJN M, RAPS D, et al. Are super-hydrophobic surfaces best for icephobicity?[J]. Langmuir, 2011, 27(6): 3059-3066. DOI:10.1021/la104762g.

[9] FEUILLEBOIS F, LASEK A, CREISMEAS P, et al. Freezing of a subcooled liquid droplet[J]. Journal of Colloid and Interface Science, 1995, 169(1): 90-102. DOI:10.1006 /jcis.1995.1010.

[10] 王皆腾, 刘中良, 勾昱君, 等. 冷表面上水滴冻结过程的研究[J].工程热物理学报, 2007, 28(6): 989-991. DOI:10.3321/j.issn: 0253-231X.2007.06.028.

WANG J T, LIU Z L, GOU Y J, et al. Study of freezing process of water droplet on cold surface[J]. Journal of Engineering Thermophysics, 2007, 28(6): 989-991. DOI:10.3321/j.issn: 0253-231X. 2007.06.028.

[11] CHAUDHARY G, LI R. Freezing of water droplets on solid surfaces: an experimental and numerical study[J]. Experimental Thermal and Fluid Science, 2014, 57: 86-93. DOI:10.1016/j.expthermflusci. 2014.04.007.

[12] HINDMARSH J P, RUSSELL 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. DOI:10.1016/S0017-9310(02)00399-X.

[13] WANG H, HE G, TIAN Q. Effects of nano-fluorocarbon coating on icing[J]. Applied Surface Science, 2012, 258(18): 7219-7224. DOI:10.1016/j.apsusc.2012.04.043.

[14] MISHCHENKO L, HATTON B, BAHADUR V, et al. Design of ice-free nanostructured surfaces based on repulsion of impacting water droplets[J]. ACS Nano, 2010, 4(12): 7699-7707. DOI:10.1021/ nn102557p.

[15] LI H, ROISMAN I V, TROPEA C. Influence of solidification on the impact of supercooled water drops onto cold surfaces[J]. Experiments in Fluids, 2015, 56(6) DOI:10.1007/s00348-015-1999-2.

[16] BLAKE J, THOMPSON D, RAPS D, et al. Simulating the freezing of supercooled water droplets impacting a cooled substrate[J]. AIAA Journal, 2015, 53(7): 1725-1739. DOI:10.2514/1.J053391.

[17] 梁超, 王宏, 朱恂, 等. 液滴撞击不同浸润性壁面动态过程的数值模拟[J]. 化工学报, 2013, 64(8): 2745-2751. DOI: 10.3969/j.issn. 0438-1157.2013.08.006

LIANG C, WANG H, ZHU X, et al. Numerical simulation of droplet impact on surfaces with different wettabilities[J]. CIESC Journal, 2013, 64(8): 2745-2751. DOI: 10.3969/j.issn. 0438-1157.2013.08.006

[18] 郑志伟, 李大树, 仇性启, 等. 液滴碰撞球形凹曲面复合 level set-VOF法的数值分析[J]. 化工学报, 2015, 66(5): 1667-1675. DOI: 10.11949/j.issn.0438-1157.20141116.

ZHENG Z W, LI D S, QIU X Q, et al. Numerical analysis of coupledlevel set-VOF method on droplet impact on spherical concave surface[J]. CIESC Journal, 2015, 66(5): 1667-1675. DOI: 10.11949/j.issn.0438-1157.20141116.

[19] 徐爽, 赵宁, 王春武, 等. 低速下水滴撞击固体表面运动模拟[J].航空动力学报, 2013, 28(3): 695-700. DOI: 10.13224/j.cnki. jasp.2013.03.030.

XU S, ZHAO N, WANG C W, et al. Simulation of droplets impacting on solid surface at low velocity[J]. Journal of Aerospace Power,2013, 28(3): 695-700. DOI: 10.13224/j.cnki.jasp.2013.03.030.

[20] BRIONES A M, ERVIN J S, PUTNAM S A, et al. Micrometer-sized water droplet impingement dynamics and evaporation on a flat dry surface[J]. Langmuir, 2010, 26(16): 13272-13286. DOI:10.1021/ la101557p.

[21] STROTOS G, ALEKSIS G, GAVAISES M, et al. Nondimensionalisation parameters for predicting the cooling effectiveness of droplets impinging on moderate temperature solid surfaces[J]. International Journal of Thermal Sciences, 2011, 50(5): 698-711. DOI:10.1016/j.ijthermalsci.2010. 11.021.

[22] KUMAR A. Solidification of impinging molten metal droplet on a cold substrate[J]. International Journal of Mechanical Engineering and Robotics Research, 2014, 3(2): 86-497.

[23] SUSSMAN M, PUCKETT E G. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows[J]. Journal of Computational Physics, 2000, 162(2): 301-337. DOI: 10.1006/jcph.2000.6537.

[24] 车得福, 李会雄. 多相流及其应用[M]. 西安: 西安交通大学出版社, 2007: 260.

CHE D F, LI H X. Multiphase Flow and Its Applications[M]. Xi’an: Xi’an Jiaotong University Press, 2007: 260.

[25] JUNG S, TIWARI M K, DOAN N V, et al. Mechanism of supercooled droplet freezing on surfaces[J]. Nature Communications, 2012, 3: 615. DOI:10.1038/ncomms1630.

[26] VOLLER V R, PRAKASH C. A fixed grid numerical modelling methodology for convection-diffusion mushy region phase-change problems[J]. International Journal of Heat and Mass Transfer, 1987, 30(8): 1709-1719. DOI:10.1016/0017-9310(87)90317-6.

2015-12-14收到初稿，2016-04-26收到修改稿。

联系人:常士楠。第一作者:冷梦尧（1991—），女，博士研究生。

Received date: 2015-12-14.

中图分类号:TK 01; O 359.1

文献标志码:A

文章编号:0438—1157（2016）07—2784—09

DOI:10.11949/j.issn.0438-1157.20151888

基金项目:国家重点基础研究发展计划项目（2015CB755803）；国家自然科学基金项目（11372026）。

Corresponding author:Prof. CHANG Shinan, sn_chang@buaa.edu.cn supported by the National Basic Research Program of China (2015CB755803) and the National Natural Science Foundation of China (11372026).

Numerical simulation of droplet impinging and freezing on cold surfaces with different wettability

LENG Mengyao, CHANG Shinan, DING Liang
(School of Aeronautics Science and Engineering, Beihang University, Beijing 100191, China)

Abstract:A strategy is presented to simulate the impact and solidification of the water droplet on different substrates. Simulations were performed using a coupled volume-of-fluid and level-set method to tracking the air-water interface and an enthalpy-porosity method to capture the phase transition. Three surface types were investigated: hydrophilic surface (contact angle 30°), hydrophobic surface (contact angle 114°) and superhydrophobic surface (contact angle 163°). The results showed that decreasing the wettability would reduce the contact time and area with the cold surface, and delay the freezing of the droplet. When temperature of the surface is higher than -15℃, the superhydrophobic surface can remain entirely ice-free under the simulated conditions. By comparing the maximum droplet spread, retraction response and time for solidification with the experimental results, the effectiveness and precision of the simulation strategy were demonstrated.

Key words:surface; droplet impact; dynamic behavior; multiphase flow; phase change; numerical simulation