非中心入射流旋转圆盘水跃现象的实验与模拟

任思远，刘春江，2，郭凯，刘辉，项文雨，谢春刚，李雪

（1 天津大学化工学院，天津 300072；2 天津大学化学工程联合国家重点实验室，天津 300072；3 自然资源部天津海水淡化与综合利用研究所，天津 300192）

由于旋转圆盘表面的液膜具有较高的传热和传质效率，入射流旋转圆盘在许多工程中得到广泛应用，如微电子制造、食品加工和航天航空工业[1]。当垂直入射的液体撞击旋转圆盘表面时，液体在离心力、惯性力和黏性力的共同作用下，在圆盘表面沿径向和切向方向流动，此时进口部分会出现液膜较薄的区域，在临界半径处液膜厚度会迅速增加，形成水跃[2]。水跃会影响旋转圆盘表面上的液膜分布，进而影响到圆盘表面传热、传质和反应过程，因此水跃引起研究人员的广泛关注。Rayleigh[3]首先提出了无黏流体在静止固体表面上的水跃模型。此后，Watson[4]研究了入射流在静止表面上的层流和湍流流动，阐明了液体黏度对水跃半径的影响。Bush等[5]在Watson的基础上研究了表面张力对入射流水跃现象的影响，提出的理论模型与实验结果具有更好的一致性。当固体表面开始旋转后，研究工作大部分集中于旋转圆盘表面的液膜厚度及流动形态，对进口区域水跃现象研究较少。Thomas 等[6-7]对旋转圆盘表面的水力学现象进行了分析，预测随着圆盘转速的增加，液膜逐渐变薄，当转速增加到一定的速度时，水跃现象将会消失，并划分了三个不同的流动区域：内部惯性主导区、过渡区和外部旋转主导区。Buyevich等[8]对入射流撞击平面形成的水跃现象进行了数值模拟，他们将进口部分划分为五个区域：入射区域、液膜黏性边界层出现前后的层流区域、水跃区域和亚临界区域。Avedisian等[9]利用高速照相系统研究了不同重力情况下的水跃现象，发现在微重力情况下，水跃直径比在正常重力情况下大。Ozar等[10]利用激光界面反射技术测量了液膜厚度沿圆盘半径的瞬时分布，验证了Thomas等提出的理论。Duchesne 等[11]发现对于相同的表面，当进口高度一定时，Fr数不会随着黏度和表面张力的变化而变化。Mohajer 等[12]的实验验证了Duchesne 的理论，但是发现表面张力的变化会影响Fr数。Wang 等[13]对旋转圆盘表面的进口部分从停滞区、发展的层流边界区、完全湍流区和水跃区四个方面展开研究，分析了旋转对水跃的影响。之后Wang等[14]利用激光多普勒测速仪测定了甘油-水体系在旋转圆盘表面形成的液膜的速度分布，并利用VOF 模型进行了数值模拟，发现水跃通常在低转速、高黏度和入口直径偏小的情况下出现。同时入口直径决定是否有水跃现象的出现，但是对水跃半径的影响很小。

上述对水跃的研究均为进口为圆盘中心的情况下进行的实验与数值模拟，但是目前在许多实际工程中存在非中心入射流的情况，例如在空间站中用于水循环的气相催化脱氨系统（vapor phase catalytic ammonia removal)中的刮膜式旋转圆盘（wiped film rotating disk)[15-16]，入射流中心与圆盘中心位于不同的位置。Niederhaus 和Hall 等[17-19]虽然对刮膜式旋转圆盘的刮刷区域进行了实验研究，但是缺少对入射流冲击旋转圆盘形成的进口区的讨论。入射流进口的偏移使得水跃的形状和水跃的半径发生了变化，液膜的厚度分布极不均匀，进而影响旋转圆盘的传热效率和运转稳定性。因此，本文将对特定条件下的非中心入射流旋转圆盘进口区域的水力学现象进行实验与模拟，研究进口流量和圆盘转速对水跃的影响。

1 实验装置及现象

1.1 实验装置

实验流程如图1 所示，实验装置由进料部分、旋转圆盘部分和图像处理部分组成。进料罐中的液体经过离心泵和流量计后，以一定的流量从进口喷嘴垂直射向旋转的水平放置的圆盘表面，在圆盘表面铺展成液膜，液体以一定的速度离开圆盘表面后，进入液体回收器中，之后流入回收罐中循环使用。高速摄像机（Photron FASTCAM SA4）位于旋转圆盘的上方，具体位置和拍照角度可以根据镜头（Nikon AF 50mm）需要进行调节。实验中照相频率为2000Hz，照片分辨率为1024×1024。实验中的液体为自来水，圆盘直径D为300mm，中心旋转轴直径Dc为50mm，圆盘表面材质为有机玻璃，与水的接触角为70°，进口直径d为3mm，进口位置与圆盘中心的距离e为50mm，进口位置与圆盘表面之间的垂直距离为10mm，水的流量Q的范围是10～30L/h，圆盘的旋转速度N的范围是150～300r/min。为了观察水在转盘表面的流动状态，实验时在水中加入亚甲基蓝固体粉末使自来水呈现深蓝色，亚甲基蓝水溶液浓度为50mg/L，密度为998kg/m3，黏度为1.004mPa·s，表面张力系数为0.072N/m。

图1 实验流程示意图

1.2 实验现象

实验中得到的液膜铺展照片如图2所示，图中进口流量为15L/h，圆盘转速为150r/min，旋转方向为逆时针方向。当圆盘的旋转速度保持稳定后，深蓝色的水由进口喷嘴射向圆盘表面，首先在进口部分聚集，随着时间的增加，水在离心力和黏性力的共同作用下沿着圆盘的旋转方向开始铺展成液膜。实验中的进口流量在10～30L/h 的范围内，在此流量范围下液体并不能完全润湿圆盘表面，会在液膜的最外侧形成水跃，水跃外侧的液膜破裂，以液线流[20]的形式流出圆盘表面。

图2 圆盘表面流动形态

中心入射和非中心入射造成的水跃现象的差异如图3 所示。图中圆盘旋转方向均为逆时针方向，进口流量和旋转速度均保持不变，d为进口直径，h为过渡区液膜厚度，H为水跃高度。本文定义R1为进口旋转方向前方的水跃半径，R为旋转方向后方的水跃半径。中心流入射时，圆盘表面的液膜分布基本呈中心对称，旋转方向前方的水跃半径和旋转方向后方的水跃半径相等。当入射流不在圆盘中心时，水跃会因为离心力和黏性力的作用发生变化，水跃位置发生了迁移，以进口中心与圆盘中心连线的垂直切面的水跃现象为例，旋转方向前方水跃半径明显小于旋转方向后方水跃半径，而且水跃半径受进口流量和圆盘旋转速度影响。

图3 水跃实验照片和示意图

2 数值模拟

2.1 物理模型

计算区域的物理模型如图4所示。研究所关注的区域仅为旋转圆盘的进口区域，而且圆盘上水跃之外的液膜对水跃的影响较小[21]，同时由图2可以看出水跃和铺展的液膜集中分布在圆盘表面靠近中心的进口区域，因此只针对进口部分进行了模型构建与网格划分，选取直径为200mm 的圆盘中包含进口区的1/4扇形区域作为计算域，计算域的高度为5mm，在可以获得进口部分水跃现象的同时减少了网格的数量，节约计算资源，缩短计算时间。

图4 物理模型示意图

2.2 计算模型

本文利用VOF（volume of fluid）多相流模型进行模拟。VOF 模型可以通过求解一组动量方程并跟踪整个计算域中的两种流体的体积分数来模拟空气和水的两相流。在每一个计算单元中，水和空气的体积分数α之和等于1，计算式为式(1)。

如果αq=1，则网格中充满了水；如果αq=0，则网格中充满了空气，如果0＜αq＜1，则表示该网格为气液分界面。VOF 方程中的连续性方程和运动方程为式(2)、式(3)。

式中，ρ为流体的密度，通过体积加权平均得到，如式(4)。

μeff为流体的有效黏度，包括分子黏度μ和湍流黏度μt，计算式为式(5)。

其中分子黏度μ通过体积加权平均得到[式(6)]。

u为速度矢量；t为时间；FV为表面张力，由连续表面CSF模型[22]定义，定义为式(7)、式(8)。

式中，σ为表面张力系数；κ是相界面的曲率；n是表面的法向向量。考虑到液体与旋转圆盘壁面的接触角θw，则壁面计算单元表面的法向为式(9)。

式中，nw和tw为壁面计算单位的法向向量和切向向量。

对于入射流撞击旋转圆盘的三维模型的计算，Benim等[23]经过比较发现Standardk-ω模型能较好地符合实验结果，特别是在进口区域部分模拟效果更好。本文研究内容主要针对进口区域的水跃现象，所以选取Standardk-ω模型进行模拟。Standardk-ω湍流模型[24]的方程见式(10)～式(12)。

式中，a为校正系数；k为湍动能；ω为比湍动能耗散率；Γk和Γω分别为湍动能和比湍动能耗散率的有效扩散系数；Gk和Gω分别为由于平均速度梯度产生的湍动能和比湍动能耗散率；Yk和Yω分别为由湍流引起的湍动能和比湍动能的耗散。

入口边界设置为水相进口，采用速度入口边界条件，出口为压力出口边界条件，操作压力为标准大气压。圆盘表面及旋转轴轴面均为无滑移固体壁面，给定绝对旋转速度，旋转轴为Y轴（圆盘中心轴），圆盘壁面与水的接触角为70°。入口边界的湍动能k和比湍动能耗散率ω计算见式(13)、式(14)。

使用ICEM 软件对物理模型进行网格划分，对靠近壁面的液体流动的区域采用六面体网格，其余空气流动的区域采用四面体网格，对壁面处的网格进行加密，壁面第一层边界层网格厚度为0.01mm，增长率为1.1，层数为10层。网格最大尺寸的范围控制在0.5～1mm 之间，网格数量的范围在486 万～1089万之间，当圆盘转速为300r/min、进口流量为20L/h 时，通过综合比较圆盘表面流体流动形态、旋转方向后方水跃和计算速度等多方面因素，如表1所示，确定网格数量为759万。

表1 网格无关性验证

模拟过程利用商业软件ANSYS Fluent 16.0 双精度浮点数求解器进行瞬态计算，压力速度耦合方式选择PISO 算法，离散格式的选用如下：梯度采用Least Squares Cell Based 方式离散；压力采用PRESTO！格式离散；动量采用二阶迎风格式；体积分数采用Geo-Reconstruct 插值，湍流动能和耗散率方程采用一阶迎风格式。时间步长设置为1×10-5s，计算残差设定为1×10-4。当出口流量稳定后停止计算。

3 计算结果与讨论

3.1 模型验证

将实验得到的照片进行图像锐化处理，再利用ImageJ软件，以圆盘直径为标尺，选取与计算域相同的区域，提取圆盘表面的液膜图像，去除液膜之外的部分，测量液膜的铺展面积。为了减小实验误差，每个操作条件下均随机选取圆盘系统运转稳定后（流体入射圆盘表面后稳定1min 以上）的10 张不同的照片进行处理，将最后得到的数据取平均值作为对应操作条件下的铺展面积。提取计算结果中的气液交界面[25]作为圆盘液膜表面，如图5 所示，从图中可以看到计算结果包含了完整的进口区域的水跃和铺展液膜，并利用与处理实验照片同样的方法统计液膜铺展面积的计算值。图6为不同进口流量下铺展面积实验值与模拟值的对比图，从图中可见计算得到的液膜铺展面积与实验值较为接近，表明只针对进口部分进行建模与数值模拟可以获得较为准确的液膜铺展情况。由图6可以看出相同转速下，随着进口流量的增大，圆盘表面的铺展面积逐渐增大，是因为进口流量的增大使得惯性力的作用增大，促进了液膜的铺展。但是当进口流量相同时，圆盘转速增加，液膜铺展面积减小。这是由于进口位于偏离圆盘中心的位置，靠近圆盘中心的液体和靠近圆盘边缘的液体同时受到离心力的作用，原本有利于液膜铺展的离心力对进口区液体的作用减小，更多的液体随着圆盘的旋转方向运动，从而铺展面积减小。实验值和模拟值的平均相对误差为9%，说明建立的计算模型可以较为准确地再现旋转圆盘进口区域的流体流动现象。

图5 液膜分布图

图6 液膜铺展面积实验值和模拟值对比

3.2 进口流量对水跃的影响

图7表示以圆盘中心为圆心，进口中心到圆心的距离为半径的环面上的液膜分布图，图中横坐标为圆弧线上各点与进口中心之间的相对距离，坐标为负表示位于进口的圆盘旋转方向前方，坐标为正表示位于进口的圆盘旋转方向后方。纵坐标表示液膜表面与圆盘表面之间的距离。不同曲线代表不同的进口流量，圆盘转速均为300r/min。由图7可知，在进口的旋转方向前方，当进口流量较低时，由于惯性作用较小，没有水跃现象；当进口流量大于20L/h 时，惯性作用增大，发生明显的水跃现象，水跃高度变化不大，但明显高于旋转方向后方的水跃高度。这是因为在旋转方向前方，液体受到的离心力和惯性力的作用方向相反，离心力会抑制液膜在旋转方向前方的铺展，在液膜边缘形成水跃时，离心力和表面张力的作用方向相同，而在旋转方向后方离心力有利于液膜的铺展，形成水跃时离心力和惯性力的作用方向相同，所以在旋转方向前方，水跃高度明显高于旋转方向后方的水跃。在进口的旋转方向后方，当流量大于15L/h时有明显的水跃现象，随着进口流量的增大，水跃半径由0逐渐增大到48.77mm，水跃高度由1.17mm 减小至0.45mm后达到稳定，稳定区液膜厚度最小为0.09mm，水跃半径和液膜厚度随流量的变化如图8 和图9 所示。是因为随着进口流量的增大，惯性力的作用增强，表面张力的作用相对减弱，对液体流动的阻碍作用减弱，液体更容易铺展。

图7 液膜分布图

图8 水跃半径随流量变化

图9 液膜高度随流量变化

结合图10 的速度分布曲线可得，当液体由进口入射到旋转圆盘表面后，在惯性力的作用下，液体在圆盘表面由入射中心向四周扩散，由于圆盘表面的黏性力作用，液膜速度略有下降，之后由于离心力的作用液膜的速度逐渐增大，最终不同进口流量的液膜速度达到相同的最大值。在液膜铺展的边缘，由于表面张力的作用，液膜堆积形成水跃，液体的动能转化为位能，速度变小。根据液膜的速度变化可以将旋转圆盘进口区的液膜流动划分为以下4 个区域：惯性力占主导作用的转换区、过渡区、离心力占主导作用的稳定区和水跃区。当没有水跃现象时，旋转圆盘进口区的液膜流动只有转化区和过渡区，液膜较厚，铺展面积较小。当旋转圆盘表面有明显的水跃现象时，流量越大，过渡区过渡越平缓，稳定区长度越长，旋转方向后方水跃半径越大，液膜更加铺展。

图10 液膜速度分布曲线

3.3 圆盘旋转对水跃的影响

图11 表示进口中心环面上的液膜分布图，不同的曲线表示圆盘旋转速度不同，进口流量均为30L/h。由图可知，在进口的旋转方向前方，圆盘旋转速度越大，离心力相对于惯性力的作用越大，水跃半径越小，水跃高度越高，当圆盘转速达到600r/min时，水跃的液膜发生破裂飞溅。在旋转方向后方，圆盘表面液膜的水跃半径随着圆盘转速的增加先增加后减小，如图12 所示，在转速为300r/min 时水跃半径最大，为48.77mm。这是因为当进口流量为30L/h时，此时的流量不足以润湿整个圆盘表面，如图2所示，当液体流动时，存在靠近圆盘中心的液膜边缘和靠外的液膜边缘，两边均受到离心力的作用，外侧受到的离心力有利于液膜向圆盘边缘铺展，内侧受到的离心力不利于液膜向圆盘中心铺展。当圆盘转速较小时，内侧离心力的作用相对于惯性力的作用较小，外侧离心力的作用占主导作用，所以圆盘转速增加，水跃半径增加，液膜更加铺展；当圆盘转速较大时，内侧受到的离心力作用增强，水跃半径随着圆盘转速的增加而减小。图13 表示进口旋转方向后方水跃高度和稳定区液膜厚度随圆盘转速变化的趋势。随着圆盘转速的增加，水跃高度由1.20mm 逐渐减小至0.32mm，稳定区液膜厚度由0.17mm逐渐减小至0.07mm。圆盘转速越大，圆盘对流体的离心作用越大，旋转圆盘上的液膜向出口运动的趋势越强，液膜越薄，水跃高度越低。

图11 液膜分布图

图12 水跃半径随转速变化图

图13 液膜高度随转速变化图

3.4 水跃半径预测

由以上分析可知，非中心入射流旋转圆盘进口区流体流动和水跃现象主要受进口流量和圆盘旋转速度的影响[13]，即主要影响水跃半径的因素为进口雷诺数ReJ和圆盘旋转雷诺数ReD。ReJ和ReD的计算式如式(15)、式(16)所示。

式中，U为进口液体的平均流速；N为圆盘旋转速度；D为圆盘直径。分析进口旋转方向后方的水跃半径与进口流量和圆盘转速的关系，拟合得到水跃半径与ReJ、ReD的经验关联式见式(17)。

对比水跃半径模拟值和预测值，如图14所示，结果表明大部分模拟值和预测值的误差在15%以内，说明经验关联式可以较好地预测水跃半径。

图14 水跃半径模拟值和预测值对比

4 结论

本文通过实验对非中心入射流旋转圆盘进口部分的流体流动进行了观察，并测量了进口部分液膜的润湿面积，建立了旋转圆盘进口部分流体流动过程的CFD模型，得到以下结论。

（1）当圆盘转速一定时，随着进口流量的增大，非中心入射流旋转圆盘进口区的铺展面积逐渐增大。实验值与模拟值的平均相对误差为9%，建立的数值模型可以较为准确地模拟实验现象。

（2）非中心入射流旋转圆盘的进口偏离圆盘中心，导致进口旋转方向前方的水跃和旋转方向后方的水跃现象明显不同。当进口流量大于15L/h时会出现明显的水跃现象，进口旋转方向前方的水跃高度最高可达2.5mm，旋转方向后方的水跃高度稳定在0.45mm。旋转方向后方稳定区液膜厚度最小可达0.07mm。

（3）在进口旋转方向前方，进口流量越大，圆盘旋转速度越小，水跃半径越大。在进口的旋转方向后方，流量越大，水跃半径越大，而旋转速度越大，水跃半径先增大后减小，当圆盘转速为300r/min时水跃半径最大。对进口旋转方向后方的水跃半径与进口流量和圆盘旋转速度的关系进行了拟合并得到关联式，关联式的预测结果与模拟结果的误差在15%以内。

符号说明

D,Dc—— 圆盘直径和中心旋转轴直径，mm

d—— 进口直径，mm

e—— 圆盘表面进口位置中心与圆盘中心的距离，mm

H—— 水跃高度，mm

h—— 过渡区液膜厚度，mm

k—— 湍动能，m2/s2

N—— 圆盘转速，r/min

n—— 法向向量

Q—— 进料流量，L/h

R1,R—— 进口旋转方向前方和后方的水跃半径，mm

Re—— 雷诺数

t—— 时间，s

t—— 切向向量

U—— 进口液体的平均流速，m/s

u—— 速度矢量，m/s

u—— 黏度，Pa·s

v—— 运动黏度，m2/s

α—— 体积分数

κ—— 曲率，1/m

θw—— 接触角，(°)

ρ—— 密度，kg/m3

σ—— 表面张力，N/m

ω—— 比湍动能耗散率，1/s

下角标

D—— 圆盘

J—— 进口

w—— 壁面