倒置平面基底上液膜流动的演化特性研究

邵明玉,张忠伟,马驰骋,荆 栋,华 珍

(山东理工大学 交通与车辆工程学院,山东 淄博 255000)

薄膜流动现象普遍存在于化学工程、生物工程、机械润滑等领域,作为一类典型的流体动力学问题,蕴含着丰富的非线性行为,引起了广大研究者的关注[1-3]。当薄膜在平面、圆柱面或者沿着细丝流动时,含运动接触线的薄膜流动存在多种有趣的形态,比如稳定行波演化,或指突不稳定。影响薄膜流动特性的因素很多,如壁面拓扑结构[4]、静电场力[5]、热毛细力[6]和滑移边界效应[7]等,而且在薄膜流动过程中往往多种因素共同作用。

从流体动力学角度来讲,薄膜流动属于低速开式的交界面流动,自由面处的非线性动力学行为分析和数值模拟一直是薄膜流动研究的重点之一。相对于理论分析和实验验证,数值模拟在薄膜下落问题研究中更加普遍,而且数值方法众多。胡军等[8]采用拉格朗日有限元法,开展了均匀加热下落薄膜的直接数值模拟；Lin等利用有限差分方法离散求解空间,用alternating direction implicit(ADI)方法进行时间迭代,虽然精度很高,但算法施加相对繁琐[9-10]；潘晓军等利用傅里叶伪谱法与Runge-Kutta法相结合开展了水平基底上薄膜流体流动的数值模拟[11-12]；Gao等基于有限体积法研究了薄膜流动[13]；Gu等采用计算流体力学方法研究了沿斜板下落薄膜的流动规律,并分析了两相流工况下的薄膜流动形态[14]。

近年来,更加新颖的仿真技术手段应用到流体动力学仿真中。李杰和Raach等采用OpenFoam软件实现了薄膜流动的二维数值模拟,获得了大量的数值结果[15-16]；Wehinger利用STAR-CCM+模拟软件开展了液膜下降流动的数值模拟,通过局部网格加密获得了比较精确的仿真结果[17]；Hu使用有限单元求解器Dolfin开展了非牛顿流体的三维薄膜流动仿真,并分析了流动的稳定性[18]。从文献检索中可以发现,尽管有限单元法早在1984年便开始用于下降薄膜的数值模拟,但是在流体动力学仿真中,有限差分方法和有限体积法的使用更加广泛。随着薄膜流动问题的深入研究,寻求一种更为简单有效的数值方法将会促进该类问题的研究。

本文在长波近似的理论框架下,推导了倒置平面基底上的薄膜流动润滑方程,然后针对该方程建立有限元模型,使用开源有限元软件Freefem++仿真薄膜二维流动和三维流动中的运动接触线演化过程,基于流动形貌的演化过程,系统地分析了倒置平面基底上的液膜流动演化规律。

1 薄膜运动方程

1.1 基本方程

3种沿平板下落的薄膜流动现象如图1所示,分别为沿斜板流动、竖直流动和倒置基底流动。这3种流动都属于重力驱动流动,前两种现象已经得到了非常广泛的研究,在前人工作的基础上,本文重点研究倒置平面基底上的薄膜流动问题。研究中假设流体为黏性不可压缩的牛顿流体,因此，满足连续方程和Navier-Stokes方程。

图1 平面薄膜流动Fig.1 Falling films on planes

(1)

(2)

(3)

(4)

图1中所示的薄膜与基底的接触面采用无滑移边界,即

u=v=w=0。

(5)

在自由界面处满足动力方程和张力平衡方程，

(6)

(7)

1.2 润滑理论

基于薄膜流动特征,采用润滑理论推导其流动方程。由于薄膜厚度相比其他两个方向上尺寸较小,因此，引入小参数ε=h0/L。其中：h0为薄膜平均厚度；L为其他方向的特征长度。无量纲坐标可以表示为

(8)

令U表示x方向的速度尺度,速度项的无量纲表达式为

(9)

通过特征长度和特征速度,可以确定特征时间,同时将体积力势函数和压强定义为

(10)

(11)

将无量纲参数表达式(8)～(11)带入连续方程和Navier-Stokes可以得到

(12)

(13)

(14)

(15)

其中，Re表示雷诺数。为简化表达，下文省略无量纲符号。本文主要研究重力驱动的薄膜流动,体积力只考虑重力(ρgsinα,0,-ρgcosα),略去式(13)～式(15)中包含ε及其高阶小量的部分，可以得到

(16)

(17)

(18)

对式(18)在z方向上积分,代入边界条件式(7)可以得到压强函数为

p=ρgcosα(h-z)-σκ+p0。

(19)

其中，平均曲率为

κ≈2h。

(20)

得到压强的表达式(19)后,对式(16)和式(17)在z方向上积分两次,利用边界条件可以得到x方向和y方向的速度表达式，

(21)

(22)

动力边界条件(6)也可以用积分形式表示，

(23)

通过平均速度进一步简化式(21)和式(22)可得

(24)

(25)

将式(24)和式(25)代入式(23)，可得倒置平面基底上薄膜流动的润滑方程，

(26)

根据表面张力和重力项的平衡条件[19],液膜厚度的演化控制方程可以进一步化简为

(27)

其中，Dθ=(3Ca)1/3cotθ,Ca=μU/σ代表毛细管数,式(27)虽然比式(26)简单,但仍然是强非线性、高阶偏微分方程,采用差分方法求解该方程依然会比较复杂。

2 有限元方程

根据有限单元法和Freefem++编译环境[20],求解微分方程时，需要首先建立微分方程的变分形式。直接建立方程(27)的变分形式比较困难,为了得到方程(27)的微分形式,引入一个新的参数φ，

φ=2h。

(28)

式(27)可以改写为

(29)

式(27)的求解转换为式(28)和式(29)的求解,虽然方程数目增加,但是求解效率更高。通过分部积分,推导出式(28)和式(29)的变分形式为

(30)

(31)

其中：H和ψ分别为h和φ的试函数；S表示计算域。

初始条件和计算域如图2所示,在Freefem++中对y=0和y=Ly处施加周期边界条件,初始条件为光滑函数连接的两部分薄膜厚度，

(32)

其中：b为前驱膜的厚度；x0为y的函数，

x0=xf-A0cos(2πy/λ)。

(33)

其中：A0和λ分别为初始扰动的幅值和波长；xf表示连接处位置。

3 倒置基底薄膜流动二维模拟

在Freefem++中无法直接进行二维仿真,但是可以通过设置一个较小的y方向尺寸近似开展二维仿真,忽略y方向上函数值的变化，得到的结果与二维模拟十分接近,完全能满足计算精度。有限元模型和边界条件如下图2所示。使用前驱膜假设,假定固体表面完全润湿,前驱膜厚度为b,液膜从左侧计算域流进,在x=0处厚度设置为1,液膜初始形貌根据式(32)给定,仿真模拟中xf为15。

图2 有限元网格和边界条件Fig.2 Finite element grids and boundary conditions

倒置平面基底上的液膜流动形貌由重力和表面张力的共同作用决定。倾斜角度越大,作为驱动力的重力作用越强,相较于表面张力的作用更大,从而更易引起液气交界面的非线性演化。文中Dθ取值为-0.61,-1.56和-1.93, 分别代表倾斜角度122°, 148°和153°。 首先， 计算Dθ=-0.61时的液膜流动铺展过程,有限元模型和不同时刻液膜的铺展演化形貌如图3所示,可以看出,液膜厚度在y方向上的变化完全可以忽略,说明采用这种近似方式开展液膜流动的二维仿真是可行的。

图3 二维问题的近似求解域和数值演化结果Fig.3 Approximate calculation domain for two dimensional problems and time evolution results

只保留图3中y=0.5处的数值结果,即得到如图4中所示的二维演化形状。从图4中可以看出,当Dθ=-0.61,倾斜角度较小时,液膜流动呈现出一种稳定的行波态,经过一段时间的演化,毛细脊的高度不再发生改变,稳定地向前推进。可以推断出,当倒置角度较小时,液气接触表面的铺展是运动稳定的,其流动稳定性与竖直基底上的液膜流动稳定性一致。当Dθ=-1.56,液膜流动规律变得十分复杂,t=10～30时间内,前进波呈现出一种强烈的衰减形式,随着时间增长,液膜呈现出一种强衰减形态和孤立波并存的形式,靠近接触线前缘逐渐演化出一系列孤立波,而靠近左侧边界处依然存在强衰减的前进波,如图4B中t=40～60时所示。在液膜铺展过程中,液膜厚的区域迁移快,而液膜薄的区域迁移慢,因此，在移动接触前缘区域能够观察到两个相邻的波峰融合成一个大的孤立波形态。而当Dθ=-1.93时,液膜运动的不稳定性更强,在较短的时间内就发展成多个孤立波存在的形式,如图4C所示。通过二维的薄膜流动仿真模拟可以确定,θ越大,液膜的流动越不稳定,主要原因是θ越大，重力驱动作用越强,导致表面张力无法约束液气接触面的稳定铺展。

图4 不同倒置角度基底薄膜流动二维演化Fig.4 Two dimensional flow simulation of the falling film on an inverted plane for different angles

4 倒置基底薄膜流动三维模拟

在液膜流动过程中,液气两相界面往往会发生失稳以指状形貌流动铺展,二维仿真可以观察到两相界面的起伏，但是无法研究其黏性指进行为。因此,进一步针对方程(27)开展三维有限元数值模拟。三维演化模拟结果如图5～图7所示。计算域为(0≤x≤100,0≤y≤40)。为了考察液膜流动的三维演化特性,施加初始条件为

x0=xf-Aicos(2πy/λi)。

(34)

其中：i=1，2，…，50；Ai和λi分别表示第i阶的扰动幅值和扰动波长。

图5～图7分别给出了Dθ为-0.61,-1.56和-1.93这3种工况下, 倒置平面基底上薄膜流动的三维演化过程。 对比图4和图5～7可以看出, 液膜的三维流动和二维模拟具有密切的联系性, 二维流动反映了三维流动的一些基本特性, 如移动接触区域的迁移速度、 流动形貌等,通过测量发现， 图4A和图5中毛细脊的移动速度和高度是完全一致的, 图4B和图4C中二维演化形貌也分别与图6和图7的演化形貌在x方向上是相对应的。

图5 Dθ=-0.61工况下不同时刻下落液膜的三维演化模拟Fig.5 Three dimensional simulation of time evolution of the falling film for case of Dθ=-0.61

图6 Dθ=-1.56工况下不同时刻下落液膜的三维演化模拟Fig.6 Three dimensional simulation of time evolution of the falling film for case Dθ=-1.56

图7 Dθ=-1.93工况下不同时刻下落液膜的三维演化模拟Fig.7 Three dimensional simulation of time evolution of the falling film for case Dθ=-1.93

从图5可以看出,当Dθ=-0.61时,液膜以固定的前进速度、稳定的行波形式向下游移动,并没有发生黏性指进现象，初始扰动Aicos(2πy/λi)

完全消失, 说明这种流动形式是稳定的。如图6所示，当Dθ=-1.56时,除了能得到x方向的演化形态,也可以得到y方向的演化形态,在t=20和t=40时,自由液面以强衰减形态为主,而经过长时间的演化,在t=60时出现了大量的孤立波,与图4B中的结果是一一对应的。经过一段时间演化,液气两相接触面在y方向上厚度出现了明显变化,即发生了黏性指进现象,说明当倾斜角度较大时,自由液面的运动是不稳定的。当倾斜角进一步增大,即Dθ=-1.93时,液膜流动的不稳定性更强,如图7所示，可以看出,经过较短时间的演化,在t=30时发生了明显的黏性指进现象,强衰减的流动模态基本消失。从图6和图7可以发现,孤立波是指形结构的重要几何特征,随着时间的增加,孤立波会成为脱落的液滴,导致计算程序的不收敛,这是由于有限单元法对连续性的要求所决定的。

5 结语

本文通过理论分析并结合数值模拟,主要研究了倒置平面基底上含运动接触线的薄膜流动问题。

1) 基于润滑近似理论,建立了倒置平面基底上薄膜流动的润滑方程,采用有限元方法基本原理,引入新参数将薄膜流动演化方程转变为两个方程,推导了薄膜流动强非线性高阶偏微分方程的变分形式。通过开源有限元软件Freefem++对薄膜流动运动接触线的演化开展了二维数值模拟和三维数值模拟,它们具有良好的对应性。

2) 数值结果表明，对于倒置基底上的薄膜流动,倒置角度越大,重力的作用越强,自由液面流动的不稳定性增强。倾斜角度较小时,液膜以稳定的行波态运动,当倾斜角度变大,可以观察到强衰减形态和孤立波。对比三维铺展形貌和二维形貌可以发现，孤立波是黏性指进的重要几何特征,因此，可以说出现孤立波意味着液气接触面有黏性指进现象的发生。

本文采用的有限元分析步骤对于分析该类问题具有一定的推广性,由于薄膜流动的运动微分方程都具有类似的形式,因此，容易将本文的工作扩展到其他类似问题中。