极地低温复杂环境中的船舶或海洋平台经常面临上层建筑结冰问题，积冰不仅会影响结构物稳性，还会降低起重机、雷达等甲板设备的可靠性，甚至危及作业人员安全[1]．上层建筑结冰现象主要由海浪与船体碰撞引起，海浪飞溅产生的细小喷雾团会在重力及风力作用下在船舶上层建筑上方运动，细小水滴会附着在冷表面而冻结，不断累积就形成积冰[2]，因此单个水滴的冻结机理研究对深入理解整个船舶上层建筑表面结冰过程至关重要．冷表面上单个水滴冻结是一个既包含水滴运动、铺展、回缩及反弹等动力学行为，又包含过冷、传热传质等热力学行为的复杂热物理过程[3]．水滴动力学行为是一个典型的气液两相流运动界面追踪问题，目前已有较多运动界面追踪方法，包括线段法、边界积分法、移动网格方法、质点网格(particle-in-cell，PIC)方法、流体体积(volume of fluid，VOF)方法、水平集(level-set)方法等．其中流体体积方法和水平集方法应用最广泛，二者均可处理具有复杂界面拓扑变化的两相流动问题，前者可保持质量守恒但无法保证界面处物理量的连续性，而后者可准确计算界面曲率及相关物理量却容易导致质量不守恒问题．文献[4]将流体体积方法良好的守恒性和水平集方法处理界面局部尖角的优势相结合，提出了流体体积-水平集耦合方法(CLSVOF)．文献[5]提出了比CLSVOF方法计算更为高效准确的耦合流体体积方法和水平集方法的VOSET方法，并在两相流计算领域得到了广泛应用．文献[6]基于VOSET方法实现了高温金属熔池中自由液面的准确追踪．固体表面粗糙度等因素会引起表面滞后效应，进一步表现为水滴在表面上运动时三相接触线处接触角的动态变化．已有研究者对动态接触角建模进行了有益探索，其中Yokoi模型[7]和Kistler模型[8]应用较为广泛，本研究将通过具体算例讨论两种模型适用性．已有较多静置液滴冻结实验研究表明[9]，水滴在凝固点不一定冻结，水滴在温度低于凝固点时仍能保持液相的现象称为过冷现象，对应水滴为过冷水滴，凝固点与过冷水滴温度的差值即为过冷度．过冷水滴冻结可分为过冷、成核再辉、凝固冻结和冷却四个阶段，其中成核再辉和主冻结两个阶段与水滴冻结直接相关．与水滴主冻结阶段相比，过冷引起的成核再辉过程经历时间极短，多数研究者选择忽略该阶段仅研究主冻结．而本研究弥补了过往研究的不足，考虑了过冷度对水滴主冻结开始时物性参数的影响．水滴主冻结过程同样是一类界面追踪问题，已有分子动力学、光滑粒子法、格子玻尔兹曼法、相场法及焓孔隙度等方法被用来追踪相变过程中的固液界面．其中焓孔隙度方法基于焓变关系引入液相分数，迭代计算每个网格的液相分数即可体现固液界面的推进变化，避免了跟踪固液界面，从而更利于编程实现和节约计算资源．该方法已被广泛应用于水滴冻结、金属熔融、储能材料相变等研究领域，计算流体力学商业软件FLUENT中凝固与融化模型亦采用了此方法处理相变过程．综上所述，水滴运动和冻结过程各自仅涉及两种界面，分别得到广泛研究，而水滴撞击冷壁面后冻结是一个涉及多界面的质量和热量传递复杂物理过程．文献[10]用FLUENT软件模拟了冷水滴撞击不同表面时的动力学及相变行为，分析了不同壁面温度条件下水滴撞击3种典型浸润性表面的结冰情况；文献[3]用耦合焓孔隙度方法和流体体积方法数值模拟了过冷度、韦伯数等因素对水滴冻结过程的影响；文献[11]建立了水滴撞击倾斜壁面后冻结数值模型，分析了倾斜角度对水滴冻结行为的影响．虽然上述研究对水滴撞击冻结物理问题进行了有益探索，但水滴过冷度及表面滞后效应的具体影响和作用机理仍须深入研究．因此，本研究采用VOSET气液界面追踪方法和焓孔隙度相变方法建立模拟单个水滴撞击冷壁面冻结过程的数值模型，该模型可同时考虑过冷度和动态接触角的影响，将本研究数值结果与相关实验结果对比，验证模型准确性，同时进一步分析撞击速度、过冷度、壁面接触角等因素对水滴撞击-冻结过程的影响．1 数值模型及步骤1.1　气液界面的追踪1.1.1　气液两相流动基本控制方程本研究将水滴视为不可压缩、层流、满足纳维-斯托克斯方程的黏性流体，质量及动量守恒方程为∇∙(u)=0；    ∂u/∂t+∇∙(uu)=[1/ρ(φ)]{-∇p+∇∙[μ(φ)∙(∇u)+(∇u)T]+ρ(φ)g-σκ(φ)∇Hε(φ)+SF},式中：ρ(φ)为密度，μ(φ)为黏度，κ(φ)为界面曲率，Hε(φ)为光顺函数，均与符号距离函数φ有关；σ为表面张力系数；p为压力；g为重力加速度；SF为达西源项，当物理场景涉及相变时出现达西源项，其起到联系气液两相流动中的动量方程和固液相变过程中的能量方程的作用．1.1.2　VOSET方法流体体积方法通过引入流体体积分数C来表征区域内气体、液体及气液界面分布情况，C=0代表气相，C=1代表液相，C在区间(0，1)内时则表示气液界面，对于不可压缩流体，C满足以下输运方程∂C/∂t+∇∙(uC)=0．(1)方程(1)可以保证C在计算域内质量守恒，而水平集方法则是将气液界面用符号距离函数φ表示，φ(i,j,t)=d    ((i,j)在气相内) ;0       ((i,j)在界面上) ;-d     ((i,j)在液相内) ,式中d为计算区域内某一点距离界面的距离，其值为0时即表示该点位于界面上．由于函数φ(i,j,t)是一个光滑连续变化函数，因此由其计算而得的界面周边物理量也是准确、连续的．VOSET方法充分利用了流体体积和水平集方法的各自优点，首先用体积分数C确定界面的大概位置，保证质量守恒性；然后用几何重构方法求得符号距离函数φ，保证界面连续性；最后根据新确定的φ值求解界面处物性参数和界面曲率等．一个时间步内的基本求解步骤如下．a. 基于当前时刻的C值求解界面法向量，确定界面的大概位置，采用分段线性(PLIC)方法将气液界面划分为16种类型，并利用界面法向量方向进一步归纳为4种类型．b. 分别计算4种类型界面的网格(i，j)内流量的流入流出情况，其他12种界面类型流量可以根据翻转等操作获得，将界面流量值代入离散控制方程获得当前时刻网格(i，j)的流体体积分数值．c. 再次采用分段线性方法重构形成16种界面，并根据当前时刻的Ci,jn+1或φi,jn+1(首次执行该步骤时使用步骤b中获得的Ci,jn+1)，利用九点差分格式求得新的界面法向量，更新4种典型界面．d. 用文献[5]中描述的几何方法求解界面附近的符号距离函数φi,jn+1，同时结合体积函数Ci,jn+1值判断φ的正负号．e. 重复以上步骤c和d达到3次，便可求得当前时刻准确的符号距离函数φi,jn+1，从而进一步求得准确连续的界面曲率κ(φ)、法向量及界面物性参数(包括密度ρ(φ)和黏度μ(φ)等)．1.2　固液界面的推进与气液界面追踪中常用的流体体积方法类似，焓孔隙度方法引入液相分数α表征区域内固液相分布，核心思路是将固相看作黏度很高的液相，即流动速度和潜热基本为零，固液混合区域为多孔介质．α=0的区域内均为固相，α=1时则全为液相．针对牛顿流体凝固/融化问题，能量方程为∂(ρτ)/∂t+∇∙(ρuτ)=∇∙((λ/c)∇τ)+Sc，(2)式中：τ为温度；λ和c分别为材料的热传导系数和比热；Sc为潜热源项．又有Sc=(L/c)[∂(ρα)/∂t+∇∙(ρuα)]，式中L为液体潜热．上式中包括关于α的非稳态项，即能量方程(2)中包括两个待求变量(温度τ和液相分数α)，二者同时满足以下关系α=0   (ττs) ;(τ-τs)/(τl-τs)  (τs≤τ≤τl) ;1   (ττl) , (3)式中τl和τs为液相及固相临界温度，分别对应0.1 ℃和-0.1 ℃．正式离散求解能量方程前，须要先完成潜热迭代获取当前时刻的液相分数αn+1，具体迭代方式及步骤参见文献[12]．1.3　VOSET方法和焓孔隙度方法耦合1.3.1　各主要变量之间的相互影响及关系本研究模拟场景主要包括水滴动力学过程和热力学过程，前者主要变量为u，后者主要变量为τ．基于能量方程求得的温度τ可进一步获得液相分数α，液相分数α可直接反映区域内固液状态，并通过动量源项SF影响相变过程中两相流动，SF=[(1-α)2/(α3+ε)]Amu，(4)式中：ε为避免分母为0的极小值，取0.001；Am为糊状系数常值，既能在α比例较高时保持液体流动性，又能在α较低时阻碍黏性流动．Am的取值与水滴过冷度有关，当过冷度为15 ℃时，Am=5×104．式(4)中包含速度分量u，离散时须将SF并入主变量离散系数．可见α建立了动量方程源项SF和能量方程源项Sc的联系，并进一步体现u和τ的相互影响．1.3.2　气液固三相的统一显示VOSET方法中的体积分数C能区分气液两相，不考虑固相；而焓孔隙度方法中的液相分数α仅考虑固液两相，不考虑气相．因此，单独使用C和α不能同时显示三相．本研究引入三相分数F，F=(C+α)/2，成功区分气液固三相，见图1和表1．10.13245/j.hust.238458.F001图1三种分数对应云图示意10.13245/j.hust.238458.T001表1三种体积分数的取值参数气相液相固相C011α010F010.52 求解方法水滴以一定初速度撞击水平壁面后会迅速向四周铺展，在未发生破碎的前提下，相同径向处水滴沿周向流动条件相同，即垂直于水平表面的水滴径向剖面沿中轴旋转一周，便可获得整个水滴铺展形态及内部流场、温度分布．本研究采用图2所示的2D轴对称模型模拟水滴铺展及回缩等运动，在保证计算准确的同时，可比3D模型节省更多计算资源，图中V0为水滴初始速度．10.13245/j.hust.238458.F002图2水滴撞击壁面问题物理模型2.1　初始物性条件在特定过冷温度范围内，水能以稳定非晶状态存在，即亚稳态平衡形式[13]，只需要微小的扰动即可触发水的凝固，导致水从亚稳态液相变成稳态固相，热/能量传导、与冷表面的接触、机械撞击等都可能成为扰动源．过冷水滴与冷表面接触后可触发成核及冰晶生长，该阶段经历时间与水滴主冻结阶段相比极短，模拟中可忽略该时间．但必须考虑过冷对物性参数的影响，因为再辉结束后，水滴由纯水状态变成冰水混合状态，其潜热值、热传导系数等均发生变化．水滴主冻结开始前的液相分数可由下式求得αmix=1-(cl/cs)St, St=cs(Δτ/L)，式中：cl和cs分别为液态水和固态冰比热容；St为无量纲斯蒂芬数；Δτ为过冷度．由此，水滴主冻结初始潜热常值及主要物性参数可插值获得：Lmix=Lαmix；ρmix=ρlαmix+ρs(1-αmix)；cmix=clαmix+cs(1-αmix)；λmix=λlαmix+λs(1-αmix)．过冷度引起的密度变化对毫米级水滴直径影响较小，忽略凝固膨胀影响．本研究中空气、纯水、纯冰的潜热、动力黏度、热传导系数等参数详见文献[14]，当Δτ=15 ℃时，主冻结开始时水滴内αmix=81.01%，对应潜热值为270.1 kJ/kg，比纯水滴对应值小，而热传导系数为1 451.1 W/(m·K)，比纯水滴对应值大，这意味着考虑过冷度后的水滴冻结时间更短，也与相关实验值更加贴合．2.2　边界条件由于空气和水滴之间的热传导系数极小，可忽略气液界面的传热，即将气液界面视为绝热，水滴冻结从冷表面开始．如图2所示的场景中，对称轴(Y轴)设置绝热、速度滑移边界条件，底部设置定温、速度无滑移、动态接触角边界条件．接触角作为表面润湿性的定量表征边界条件，对水滴撞击、铺展等动力学行为有较大影响，因此动态接触角模型的选取至关重要，其中Kistler模型通过迭代求解实时动态接触角θd，θd=fHoff(Ca+fHoff-1(θe))，(5)式中：θe为平衡接触角；Ca为无量纲毛细数，Ca=μUcl/σ，μ为材料动力黏度，Ucl为三相接触线速度，σ为表面张力系数；fHoff(χ)为Hoffman函数，    fHoff(χ)=arccos{1-2tanh[5.16[χ/(1+1.31χ0.99)]0.706]}.式(5)中θe根据Ucl符号而取前进接触角θadv或后退接触角θrec，θe=θadv  (Ucl≥0);θrec  (Ucl0).Yokoi模型的形式更加简洁，θd=min[θe+(Ca/ka)1/3,θadv]  (Ucl≥0);max[θe+(Ca/kr)1/3,θrec]  (Ucl0),式中ka和kr为与材料物性相关的经验常数，分别取ka=9.0×10-9，kr=9.0×10-8．图3给出了上述两种模型分别预测的动态接触角曲线和文献中实验及数值结果的对照，由图可见，Kistler模型在较大接触线速度范围内仍较为准确，因此本研究选用该模型作为动态接触角边界条件．10.13245/j.hust.238458.F003图3不同动态接触角模型与相关文献结果对比[3, 7, 8]2.3　计算条件本研究用有限体积法在交错网格中离散控制方程，选用IDEAL算法[15]耦合求解速度与压力场，并用延迟修正技术和MUSCL格式实现对流项的高阶离散．采用交替方向隐式(ADI)方法解代数方程组，时间步长满足收敛条件判断数(CFL)，其中库朗数取0.1．3 结果分析与讨论3.1　网格无关性验证选取五种不同网格尺寸，模拟了10 μL水滴在-16.5 ℃的冷壁面(接触角为90°)上的冻结过程，结果如图4所示．当网格尺寸为35 μm×35 μm时，无量纲中心冻结高度h(固液界面中点高度与水滴冻结完成时顶端高度的比值)演变曲线已趋光顺，且继续减小尺寸变化较小(可忽略)，因此选取35 μm×35 μm作为静置水滴模拟网格尺寸．文献[16]已基10.13245/j.hust.238458.F004图4网格无关性分析于水滴运动算例进行了网格无关性讨论，确定20 μm×20 μm为适宜模拟运动水滴冻结的网格尺寸．3.2　模型准确性验证过冷水滴撞击冷壁面冻结模型的验证主要包括静置水滴冻结过程和水滴撞击后铺展回缩过程．静置水滴冻结过程的典型特征是水滴内部会有一条近似水平的固液分界线从底端向上推进，到达顶端后水滴完成冻结，该过程主要受冷表面温度、接触角及水滴体积影响．本研究选取张旋等[17]的冻结实验数据与本研究数值结果对比，实验中铝板接触角为78°，冷表面温度为-18.4 ℃，水滴初始体积为20 μL，模拟时选用上述条件，各参数与实验条件一致．图5可见本文数值结果与对应时刻文献[17]的实验数据基本符合，证实本研究静置水滴冻结模型较准确．10.13245/j.hust.238458.F005图5冷表面上过冷水滴冻结过程结果对比无量纲铺展因子D/D0(水滴铺展直径与初始直径的比值)可以准确反映水滴撞击水平壁面后的铺展、回缩及反弹等动力学行为．本文验证算例选择与文献[18]一致的实验条件(水滴初始直径为2.84 mm，初始速度为0.7 m/s，过冷度为30 ℃，超疏水表面的静态接触角为160°)，实验中无量纲铺展因子随时间变化及本研究对应模拟结果如图6所示，可见本文模拟结果与相关实验结果基本一致．10.13245/j.hust.238458.F006图6过冷水滴撞击固体表面后铺展因子随时间变化3.3　撞击-冻结过程影响因素分析3.3.1　过冷度过冷度是水滴冻结的驱动力，对水滴主冻结阶段初始物性参数有较大影响．船舶上层建筑结冰环境温度一般不低于-29 ℃，为系统分析过冷度的影响，每隔5 ℃设置一个算例，模拟结果见图7．10.13245/j.hust.238458.F007图7过冷度对水滴撞击及冻结过程的影响图7(a)表明：过冷度增加到一定程度时(如25 ℃)，水滴不再像常温或低过冷度水滴那样经历铺展、回缩后弹离冷壁面，而是经过少许回缩后粘附在冷表面，且过冷度越大最大铺展因子越小．由图7(b)可知：过冷度越大水滴冻结越快，这是因为当过冷度较大时水滴内液相分数较小，其导致水滴内部剩余潜热减少和热传导系数增大，进而加速冻结，上述现象与过冷度越大冻结越快规律相一致．3.3.2　壁面接触角壁面接触角是表面润湿性的定量表征，当表面接触角θsta＜90°时，该表面呈现亲水性，水滴在表面上部分或完全润湿；当表面接触角θsta ≥ 90°时，该表面为疏水性表面，水滴在表面上不润湿．图8(a)表明：表面接触角对水滴铺展及回缩行为有显著影响，较小的壁面接触角导致了更大的最大铺展因子，且亲水性表面上的水滴基本不会弹离表面；同时发现接触角对在超疏水表面上运动的水滴影响较小．图8(b)曲线显示表面接触角越小则冻结越快，一方面当接触角较小时，水滴在表面润湿面积更大，意味着固体和水滴之间导热面积更大；另一方面同样体积的水滴在亲水性表面上平衡时的高度更小，其导热热阻也更小．10.13245/j.hust.238458.F008图8表面接触角对水滴撞击及冻结过程的影响3.3.3　撞击速度初始撞击速度能够使水滴获得初始动能，该动能在水滴撞击后转变为表面能及内部动能，还有一部分动能在铺展及回缩过程中被黏性耗散掉．由图9可见：当其他条件不变时(D0=2.50 mm，Δτ=5 ℃，θsta=135°)，撞击速度对水滴铺展、回缩等动力学行为的影响较为直观；撞击速度越大，则水滴的最大铺展直径也越大；撞击速度较小时，水滴初始动能及惯性不足以克服水滴自身的表面能，因此无法在回缩后弹离表面．由图9还可发现撞击速度对水滴弹离表面的时间基本没有影响．10.13245/j.hust.238458.F009图9速度对水滴撞击及冻结过程的影响4 结论本研究结合VOSET气液界面追踪方法和焓孔隙度相变方法，对考虑过冷度和动态接触角影响的水滴撞击冻结过程进行了数值模拟，模拟结果与相关实验值具有较高一致性，得到以下结论．a. 过冷度能够加快水滴的冻结进程，并阻碍水滴的铺展及回缩过程，较大的过冷度导致较小的水滴最大铺展直径．b. 表面滞后效应导致水滴与固体表面之间的接触角呈现动态变化特点，纳入Kistler动态接触角模型的数值结果与相同条件下的实验结果更加匹配，且接触角越小，水滴冻结越快，最终冻结时间越短．c. 水滴初始撞击速度越大，则水滴在固体表面上的最大铺展直径也越大，也更易弹离表面，当撞击速度达到临界速度时即可发生破碎、飞溅等行为，可见速度越大，水滴越不易在物体表面附着冻结．综上所述，水滴过冷度、表面接触角及撞击速度均会影响单个水滴在结构表面的附着冻结过程，进而影响船舶上层建筑表面多水滴组成的飞沫团撞击冻结过程；当过冷度及接触角较大时，水滴以一定速度撞击结构表面立即冻结，后续水滴撞击后也发生黏附冻结，结构表面易形成霜冰；而当过冷度及表面接触角较小、撞击速度较大时，水滴撞击后粘附未冻结，后续水滴继续撞击到结构表面形成水膜，水膜内伴随冰层生长．本研究得出的相关影响规律可推广应用于船舶上层建筑构件表面大量水滴附着结冰机理研究，进而为提升船舶上层建筑结冰预测方法准确性提供参考．