基于SPH方法水滴撞击恒温壁面流动与换热的研究

基于SPH方法水滴撞击恒温壁面流动与换热的研究
严圣林;陈榴;戴韧
【摘要】基于连续表面张力模型,提出了一种液滴和壁面表面张力的处理方法,结合光滑粒子流体动力学方法,建立了模拟大密度差表面张力问题的SPH算法.通过模拟水滴撞击壁面后的运动过程,与实验结果进行对比,验证了算法的有效性.应用该方法,先后模拟了单个低温大水滴,以不同直径和不同速度,撞击不同表面接触特性恒温壁面的运动过程.结果表明:当壁面的接触角在一定范围内,水滴的吸热率随接触角的增大而减小;而水滴的吸热率随着撞击速度和初始直径的增大均增大.
【期刊名称】《能源工程》
【年(卷),期】2018(000)001
【总页数】8页(P1-8)
【关键词】水滴;SPH;表面张力;CSF方法;撞击特性
【作者】严圣林;陈榴;戴韧
【作者单位】上海理工大学能源与动力工程学院,上海200093;上海理工大学能源与动力工程学院,上海200093;上海理工大学能源与动力工程学院,上海200093【正文语种】中文
【中图分类】O359
0 引言
湿蒸汽在叶片流通道中的流动和传热十分复杂，目前对液滴与壁面碰撞后的流动和
传热特性尚缺乏深入认识，不同运动参数的液滴，在不同表面特性的固壁表面的过程传热强化机制,仍有待研究。

液滴撞击固壁流动是一个瞬态的、多相的复杂流动
现象。

环境流体与液滴的密度相差三个数量级以上，环境气相的流动速度与液滴速度相差一个数量级以上。

在撞击过程中，首先液滴会发生塑性变形，液滴与环境气相的界面是不稳定的、可变的。

其次，液滴运动中的受力也是复杂的，除了常规流体微团之间的粘性剪切力，还有液滴与壁面和环境气相微团之间的表面张力作用。

实验展现了液滴撞击加热固壁时的表观运动规律。

施明恒[1]实验研究液滴撞击加
热壁面，发现液滴在铺展阶段的表面热流与液滴最大铺展系数和铺展时间是相关的。

CHANDARA等[2]实验研究丙烷液滴撞击不同温度的加热壁面，发现液滴撞击低
温壁面后铺展但不会反弹，而撞击高温壁面铺展后会反弹。

CASTANET等[3]研究水滴撞击镍制光滑壁面，发现液滴在铺展阶段中水滴的温度变化取决于水滴撞击壁面的法向速度。

但是由于实验技术所限，实验过程中很难量化测量液滴流动参数，难以建立其中的传热与传质规律。

数值模拟是研究液滴微观运动规律的重要技术手段。

自由液面运动模拟的
VOF(volume of fluid, VOF)方法是液滴运动研究的常用方法之一。

STROTOS等[4]数值模拟液滴撞击热壁面，指出液滴的换热量随着铺展面积增加而增加。

高珊
等[5]研究单个液滴撞击热表面的流动和传热过程，指出增加液滴撞击速度可以提
高液滴的吸热量。

但VOF数值算法的稳定性和准确性存在问题，随着时间的推近，界面变得模糊，需要计算量更大的网格自适应技术弥补[6]。

光滑粒子流体动力学[7-8](smoothed particle hydrodynamics，SPH)方法是液
滴运动模拟的另一类常用方法，它是基于描述质点运动的拉格朗日思想的无网格算法。

在流体界面不存在数值扩散，对界面不连续问题和离散相问题，具有很好的自适应性，适合模拟大变形问题,例如水下爆破、微注射成型、涌浪等[9-11]。

使用SPH方法模拟液滴撞击固壁时，首先需要解决液滴与环境流体表面张力以及液滴
与壁面浸润性问题。

目前，SPH方法中添加表面张力方法主要有连续表面张力(continuum surface force，CSF)模型和粒子间作用力(inter-particle interaction force，IIF)模型，模拟真空中液滴变形问题效果都比较好。

但是模拟
液滴撞击固壁时，不仅涉及到液滴粒子间的相互作用力，还需考虑液滴与壁面的浸润性问题。

由于液滴撞击壁面时，液滴粒子在对冲作用下很容易发生飞溅，表面张力难以施加,使得如何处理液滴与固壁间的浸润性成为应用SPH模拟液滴撞击固壁问题的难点。

强洪夫等[12]在含CSF表面张力的SPH方法中引入BRACKBILL等[13]提出的壁面附着力边界条件处理方法，但使用该固壁处理方法模拟液滴撞击固壁时会产生非物理性的飞溅；苏铁熊等[14]的粒子间相互作用力模型的固壁处理方法的作用力参数受到粒子间的初始间距影响而难以确定，给模拟带来了更多的工作量。

基于连续表面张力模型，本文中提出一种液滴和壁面浸润性的新处理方法，建立了一种气相氛围中液滴撞击壁面的SPH计算方法。

该SPH方法处理液滴和壁面浸润性时，计算量小，操作简单。

计算过程中只需要将固壁粒子分为与液滴接触范围内的粒子和未与液滴接触的粒子，定义其色函数分别为不同的值，然后计算表面张力。

本文中应用该方法模拟单个水滴撞击恒温壁面后的运动过程，分析了撞击壁面的浸润性以及水滴的初始直径和速度对撞击过程和换热的影响，可为今后深入研究气流中液滴撞击固壁的传热和冲击力，提供基本的理论分析工具。

1 计算模型
1.1 SPH控制方程
气相环境中大水滴撞击壁面是水和环境气体的两相流动。

由于密度上存在差异，因此传统的SPH方程无法适用于该问题。

OTT等[15]对传统SPH方法中连续性方程进行了改进，用粒子密度代替质量密度，解决了模拟大密度差多相流流动时交界面计算不稳定的问题。

该方法适合本研究大水滴与环境气体之间的密度差较大的情况，
方程为：
(1)
式中：mi为粒子i的质量；uij=ui-uj，ui、uj分别为粒子i、j的速度；
Wij=W(xi-xj,h),xi、xj分别为粒子i、j的空间坐标，h为光滑长度；为核函数对粒子i坐标的空间导数。

动量方程采用Adami多相流的动量方程[16]：
(2)
式中：第一项为压力项，其中
(3)
在同一相中对粒子i和粒子j的压力取平均值，且能保证在不同相的密度不连续交界面处的的连续性；Vi为粒子i的体积，
根据Adami的多相流方程推导出压力、粘性和表面张力引起能量变化，再加上对流传热引起的能量变化，方程如下：
(4)
式中：κ为热传导系数，Ti为粒子i的温度。

1.2 表面张力模型
BRACKBILL等[13]提出CSF表面张力模型，MORRIS[17]根据CSF模型思想，将不同流体属性的粒子标定为不同色值，以此来区分界面处流体的不同属性，并利用色函数得到表面方向及曲率，进而得到表面张力。

定位表面
(5)
式中：cj为色函数，液滴的标值为1，空气的标值为0。

界面法向
(6)
求解曲率
(7)
最后得到单位质量的表面张力公式为
(8)
式中：ρ0、ρ1为交界面处的两相流体的密度，σ为表面张力系数。

1.3 壁面浸润性模型
在液滴的表面张力和壁面的浸润性的共同作用下，液滴在撞击壁面过程中达到最大铺展直径后开始出现回缩和反弹等现象。

以液滴与壁面的接触角的大小来反映液滴与壁面的浸润性。

图1为液滴与壁面接触时示意。

图1中的θ为液滴和壁面的接触角，壁面粒子的属性和液滴粒子的属性相同。

本文中将固体壁面的粒子分为与液滴产生接触和未与液滴接触两部分，如图1所示a点和b点之间的固体壁面的粒子是和液滴产生接触的粒子，其余的壁面粒子是未接触粒子。

对壁面粒子进行色函数定义时，将式(5)中的色函数c分为内外两部分cinterior和coutside。

图1中a点和b点之间的固体壁面粒子的色函数为cinterior，其余粒子的色函数为coutside。

固体壁面为亲水表面(接触角小于90°)时，设定cinterior为1，coutside值为待定；固体壁面为疏水表面(接触角大于90°)时，设定coutside为0，cinterior值为待定。

亲
水表面的coutside值和疏水表面的cinterior值与接触角的关系的确定是通过设置不同亲水表面的coutside和疏水表面的cinterior的值进行计算，待液滴稳定后，测量液滴和壁面的接触角，再线性拟合出亲水表面的coutside值和疏水表面cinterior值与接触角的关系式。

图1 水滴与壁面接触示意
图2为亲水表面(cinterior=1)时，通过对计算值的拟合等得到coutside和接触角的线性关系。

图3为疏水表面(coutside=0)时，通过对计算值的拟合得到cinterior和接触角的线性关系。

图2 水平亲水表面，静态接触角与coutside的关系
图3 水平疏水表面，静态接触角与cinterior的关系
2 计算方法模型验证
应用本文中开发的SPH算法模拟实验结果[18]中的水滴直径2.7 mm，撞击速度0.99 m/s，壁面的接触角为121°的工况作为算法验证。

图4为水滴撞击接触角为120°壁面时SPH计算结果与实验结果[18]的对比。

从图4的对比中可以发现，SPH计算所得到的水滴的铺展过程比实验提前，但两者几乎同时达到最大铺展系数。

SPH计算的最大铺展系数与实验相差1.6%。

图5为水滴以0.99 m/s撞击接触角为121°固体壁面的SPH计算的水滴形态与实验对比(已经将空气粒子隐去),图中的物理时间后的括号内的数值为与其对应的无量纲时间。

从对比图中可以看出，在撞击过程中的18 ms之前，液滴形状存在着一定差异，其他时间数值模拟和实验结果基本吻合。

通过上述对比得出，SPH计算结果与文献实验结果基本相同，但在铺展阶段中SPH方法与实验结果存在着一定差异，因为本文中的固壁处理方法并未考虑固壁表面粗糙度的影响，这也是以后对该方法进行改进的方向。

图4 水滴撞击接触角为121°壁面，SPH计算结果与实验结果[18]对比
图5 水滴以0.99 m/s撞击接触角为121°壁面的SPH计算的水滴形态与实验对比3 计算结果及分析
数值模拟中，水滴的温度为20 ℃，壁面的温度为40 ℃，水滴的密度为1×103 kg/m3，水滴的表面张力为0.0728 N/m，水滴的动力黏度为1×10-3 Pa·s。

结果分析中的无量纲数β、t*和ε中的定义分别如下：
无量纲时间
(9)
铺展系数
(10)
式中：u为水滴的撞击速度；t为水滴撞击壁面后物理时间；d为水滴撞击壁面后的铺展直径；D是液滴撞击初始的直径。

吸热率ε为实际吸热量与最大显热吸热量之比，表示液滴达到最大铺展直径时的吸热量占由于液滴与壁面温差所引起的最大换热量的比重：
(11)
式中：Qsp为水滴达到最大铺展直径的传热量，为方程(4)右端第四项产生的能量变化；m为水滴的质量；cp为水滴的比定压热容；ΔT为Tw-Tl，Tw和Tl分别为壁面和水滴的初始温度。

3.1 固壁材料浸润性对撞击过程的影响
在研究固壁材料浸润性对撞击过程的影响时，水滴的初始撞击速度为1.5 m/s，水滴的初始直径为1 mm，水滴与壁面的接触角分别选取25°、60°、75°、90°、120°和150°。

由此研究从水滴撞击壁面到达最大铺展直径时的流动和换热过程。

图6为从上述计算结果中选取的接触角为25°、75°和120°水滴撞击壁面后的瞬时状态。

图7为水滴撞击不同浸润性壁面后铺展系数随时间的变化关系。

水滴撞击
壁面后，铺展阶段水滴的铺展系数增大，壁面的浸润性对铺展阶段中水滴的形态的影响不是很明显。

随后水滴在表面张力的作用下克服粘性力开始回缩，由于壁面的浸润性不同，造成水滴回缩时的形态存在差别。

当壁面接触角为25°时，水滴边缘不回缩，只有水环向水滴的中心逐渐聚拢；当壁面接触角为75°时，水滴的边缘和水环同时向水滴的中心聚拢，水滴中心逐渐增厚；当壁面接触角为120°时，水滴
回缩幅度比75°壁面上的水滴更大，并且在回缩过程中出现了振荡。

这是因为：水滴与壁面的内聚力随接触角的增大而减小，水滴达到最大铺展系数时，25°壁面与水滴的内聚力较大，水滴无法回缩；而75°壁面和120°壁面的内聚力小于水滴与
空气的表面张力，水滴出现了回缩；120°壁面的内聚力最小，水滴回缩幅度较大，并且在回缩过程中水滴的表面势能和动能相互转换，使得水滴的铺展系数出现了振荡。

图6 直径1 mm的水滴以1.5 m/s速度撞击不同壁面的瞬时状态
图7 水滴撞击不同接触角壁面后铺展系数随时间变化
图8为水滴撞击不同接触角壁面后热流量Φ随时间变化关系。

从图8中发现：在铺展阶段，接触角25°壁面上的水滴的热流量的增长速度最快，接触角75°壁面上的水滴的热流量的增长速度次之，接触角120°壁面上的水滴的热流量的增长速度
最慢；水滴达到最大铺展系数后进入回缩阶段，在该阶段接触角120°壁面上的水
滴的热流量的衰减速度最快，接触角75°壁面上的水滴的热流量的衰减速度次之，接触角25°壁面上的水滴的热流量的衰减速度最慢。

这是因为：在铺展阶段，由于壁面对水滴的内聚力随接触角的减小而增大，使得水滴在小接触角壁面铺展更为迅速，水滴与壁面的接触面积也大于大接触角壁面，因此水滴和小接触角壁面间的热流量也高于大接触角壁面；在回缩阶段，水滴与壁面接触部分温度升高以及接触面
积的减小，导致水滴的热流量减小，由于水滴与壁面的内聚力随接触角的增大而减小，使得水滴在大接触角壁面回缩更为迅速，水滴与壁面的接触面积减小速度也高于小接触角壁面，因此水滴和大接触角壁面间的热流量也低于小接触角壁面。

图8 水滴撞击不同接触角壁面后热流量随时间变化
图9为接触角与最大铺展系数βmax及对应的铺展时间的关系。

图10为接触角与水滴达到最大铺展直径的传热量Qsp及对应的吸热率ε的关系。

从图9和图10
中发现：当接触角在25°～90°的范围内，随着水滴与壁面的接触角增大，水滴的
最大铺展系数和对应的铺展时间、传热量和对应的吸热率均随之减小；当接触角大于90°时，水滴的最大铺展系数和对应的铺展时间、传热量和对应的吸热率基本不变。

这是由于壁面对水滴的内聚力随接触角的增大而减小，内聚力对铺展阶段水滴的铺展过程影响也在减小，致使水滴的最大铺展系数较小和其对应铺展过程时间也较短,也就造成了水滴达到最大铺展直径时的传热量减小，使得水滴达到最大铺展
直径的吸热率降低。

图9 接触角与最大铺展系数及对应无量纲时间的关系
图10 接触角与达到最大直径的传热量及对应的吸热率的关系
3.2 水滴的速度对撞击过程的影响
为了研究水滴的速度对水滴撞击壁面到达铺展最大直径时的流动和换热过程的影响，在数值模拟中，将水滴和壁面的接触角设为75°，水滴的初始直径设为1 mm，水滴的撞击速度分别选取为1、1.25、1.5、1.75和2 m/s。

图11为水滴以不同速度撞击壁面时，铺展系数随无量纲时间的变化。

图12为水
滴以不同速度撞击壁面时，热流量随时间的变化。

图13为水滴以不同速度撞击壁面时，撞击速度与最大铺展系数及对应的无量纲时间的关系。

图14为水滴以不同速度撞击壁面时，撞击速度与达到最大铺展系数时的吸热量和吸热率的关系。

从图11到图14可以看出：随着水滴的撞击速度的增大，在铺展阶段，同一时刻的铺
展系数和热流量同时增大，水滴的最大铺展系数及对应的铺展时间、传热量和吸热率也增大；在回缩阶段，水滴的回缩幅度随着撞击速度的增大而增大，同时水滴与壁面的热流量减小速度也随之增大。

这是因为：水滴撞击速度的增大使得水滴具有更大冲量，促使水滴在壁面具有较高铺展速度，致使水滴的铺展系数增长较快，也就造成水滴与壁面具有更大的接触面积，使得水滴的热流量较高；当水滴达到最大铺展系数时，水滴的最大铺展系数随着撞击速度增大而增大，与其对应的传热量也增大，因此吸热率也得到了提高；水滴速度增大造成水滴达到最大的铺展系数时表面势能也随之增大，致使速度大的水滴在回缩阶段能释放较多的表面势能，从而加速水滴的回缩，也就造成了水滴接触面积迅速减小，使得水滴的热流量迅速减小。

图11 直径为1 mm的水滴以不同速度撞击壁面时，铺展系数随时间的变化
图12 直径为1 mm的水滴以不同速度撞击壁面时，热流量随时间的变化
图13 水滴直径为1 mm，撞击速度与最大铺展系数以及对应的无量纲时间的关系图14 水滴直径为1 mm，撞击速度与达到最大铺展系数时的吸热量及吸热率的关系
3.3 水滴的直径对撞击过程的影响
为了研究水滴的直径对水滴撞击壁面到达铺展最大直径时的流动和换热过程的影响，在数值模拟中，将水滴和壁面的接触角设为75°，把水滴的撞击速度设为1 m/s，水滴的直径分别选取为1、1.25、1.5、1.75和2 mm。

图15为不同直径的水滴撞击壁面时，铺展系数随无量纲时间的变化。

图16为不
同直径的水滴撞击壁面时，热流量随时间的变化。

图17为不同直径的水滴撞击壁面时，初始直径与最大铺展系数以及与其对应的无量纲时间的关系。

图18为不同直径的水滴撞击壁面时，初始直径与达到最大铺展系数时的传热量和吸热率的关系。

从图15到图18中可以看出：随着水滴的初始直径的增大，在铺展阶段同一时刻
的铺展系数和热流量同时增大，水滴的最大铺展系数以及与其对应的铺展时间、传
热量和吸热率也增大；在回缩阶段回缩幅度随着水滴的初始直径的增大而增大，同时水滴与壁面的热流量减小速度也随之增大。

这是因为：水滴撞击速度一定时，初始直径的增大使得其具有较大的动能，致使水滴的铺展系数增长较快，也就造成水滴与壁面具有更大的接触面积，水滴的热流量较高；当水滴达到最大铺展系数时，水滴具有较大表面势能，从而使得水滴的铺展系数较大，与其对应的吸热率也得到提高；较大直径的水滴在回缩阶段能释放较多的表面势能，从而加速水滴的回缩，也就造成水滴接触面积迅速减小，使得水滴的热流量迅速减小。

图15 不同直径的水滴以1 m/s速度撞击壁面时，铺展系数随时间的变化
图16 不同直径的水滴以1 m/s速度撞击壁面时，热流量随时间的变化
图17 水滴撞击速度为1 m/s，初始直径与最大铺展系数以及对应的无量纲时间的关系
图18 水滴撞击速度为1 m/s，初始直径与达到最大铺展系数时的吸热量和吸热率的关系
4 结论
本文基于SPH方法，提出了一种液滴和壁面表面张力的处理方法，即将固壁粒子
分为与液滴接触范围内的粒子和未与液滴接触的粒子，定义其色函数分别为不同的值，并给出色函数值与接触角的关系。

应用该方法模拟水滴撞击石蜡表面并将该方法的计算与实验结果进行对比，发现本方法计算出最大铺展系数与实验仅相差
1.6%，从而验证了本方法的可行性，但是铺展过程和实验结果还存在一定的误差，还需要进一步研究和改进。

应用SPH方法模拟低温水滴撞击恒温壁面的撞击过程
和吸热率，得出以下结论：
(1)当接触角在25°～90°范围内，水滴的最大铺展系数及对应吸热率随着接触角的增大而减小；当接触角大于90°时，水滴的最大铺展系数及对应的吸热率基本不变。

(2)随着水滴的撞击速度的增大，水滴的最大铺展系数及对应无量纲铺展时间均增
大。

水滴速度增大，水滴到达最大铺展直径时，吸热量增大，同时吸热率得到提高。

(3)随着水滴的初始直径的增大，水滴的最大铺展系数及对应无量纲铺展时间均增大。

水滴直径增大时，水滴到达最大铺展直径时，吸热量增大，同时吸热率也得到提高。

参考文献：
[1] 施明恒.液滴碰击加热壁面时的传热分析[J].工程热物理学报，1985，6(4)：68-71.
[2] CHANDARA S, AVEDISIAN C. On the collision of a droplet with a solid surface[J].Proceedings：Mathematical and Physical Sciences，1991，432：13-41.
[3] CASTANET G, LIENART T, LEMOINE F. Dynamics and temperature of droplets impacting onto a heated wall[J].International Journal of Heat and Mass Transfer，2009，52(3)：670-679.
[4] STROTOS G, GAVAISES M, THEODORAKAKOS A, et al.Numerical investigation of the cooling effectiveness of a droplet impinging on a heated surface[J]. International Journal of Heat and Mass Transfer，2008，51(19-20)：4728-4742.
[5] 高珊，曲伟，姚伟.喷雾冷却中液滴冲击壁面的流动和换热[J].工程热物理学报，2007，28(S1)：221-224.
[6] 强洪夫，陈福振，高巍然.基于SPH方法的表面张力修正算法及其应用[J].计算力学学报，2011，S1：37-42.
[7] MONAGHAN J J. An introduction of SPH[J]. Computer Physics Communication，1988，48：89-96.
[8] LIU G R, LIU M B.光滑粒子流体动力学——一种无网格粒子法[M].韩旭译.长沙：
湖南大学出版社，2005.
[9] SWEGLE J W, ATTAWAY S W.On the feasibility of using smoothed particle hydrodynamics for underwater explosion calculation[J]. Computational Mechanics, 1995,17:151-168.
[10] 张响，李倩，石凡，等.基于SPH方法的微注射成型数值模拟[J].化工学报，2012，63(1)：157-162.
[11] MONAGHAN J J.Simulating free surface flow with SPH[J]. Journal of Computational Physics，1992，110：399.
[12] 强洪夫，刘开，陈振福.液滴在气固交界面变形移动问题的光滑粒子流体动力学模拟[J].物理学报，2012，61(20)：204701.
[13] BRACKBILL J U, KOTHE D B, ZEMACH C. A continuum method for modeling surface tension[J].Journal of Computational Physics，1992，100：335-354.
[14] 苏铁熊，马理强，刘谋斌，等.基于光滑粒子动力学方法的液滴冲击固壁面问
题数值模拟[J].物理学报，2013，62(6)：331-337.
[15] OTT F, SCHNETTER E. A modified SPH approach for fluids with large density differences[J].Arxiv Physics/0303112v3 [p-ph],2003. [16] ADAMI S, HU X Y, ADAMS N A. A new surface-tension formulation for multiphase SPH using a reproducing divergence approximation[J]. Journal of Computational Physics，2010，229(13)：5011-5021.
[17] MORRIS J P. Simulating surface tension with smoothed particle hydrodynamics[J].International Journal for Numerical Methods in Fluids，2000，33：333-353.
[18] 李维仲，朱卫英，权生林，等.液滴撞击水平固体表面的可视化实验研究[J].热
科学与技术，2008(2)：155-160.。