竖直平板上含活性剂液膜排液过程不稳定现象的数值模拟

2022-06-02李春曦刘成治叶学民

电力科学与工程 2022年5期

关键词：月面活性剂波纹

江蕊，李春曦，刘成治，叶学民

（河北省低碳高效发电技术重点实验室（华北电力大学），河北保定 071003）

0 引言

浸涂是将基板从溶液池中浸渍提拉，经排液干燥，使表面形成一层均匀薄膜的过程[1]。该工艺的作用是保护基板材料、改善表面性能，故而被广泛用于工业领域[2,3]。

在排液过程中，液膜表面的均匀性和稳定性是影响涂层质量的关键因素[4]。深入研究该过程中液膜的流动特性和不稳定现象的发展机理，对提高涂层质量具有重要意义。

目前，关于平板液膜排液不稳定现象的理论研究大多采用线性稳定性分析[5,6]，确定出液膜的临界稳定状态及最危险波数，而很少直接模拟其不稳定现象。虽然稳定性分析能筛选出令液膜失稳的因素，但为验证理论分析结果，还需通过数值模拟再现平板液膜排液中的不稳定现象。

与平板液膜排液具有相似数理模型的线框液膜排液问题已有成功算例可供借鉴。文献[7]对附着于线框表面的液膜建立了二维排液模型，在初始条件中加入余弦波展向扰动，模拟了液膜表面的不稳定现象，发现Marangoni效应和表面黏性会影响液膜流动的稳定性。文献[8]针对线框液膜排液的不稳定现象进行了进一步分析，指出：活性剂分布的不均匀性随排液发展而加剧，Marangoni效应对扰动发展具有重要影响。

上述研究表明：Marangoni效应与液膜表面的不稳定现象紧密相关。虽然平板液膜的排液模型与线框具有相似性；但相比线框排液过程，平板液膜还需满足表面无滑移和无通量等边界条件[9]。此外，对于含活性剂液膜，其排液过程中可能出现的不稳定性现象及其内在机理也尚待研究。

为此，本文拟建立平板表面含活性剂液膜的排液模型，推导表面速度、液膜厚度及活性剂浓度的二维方程组；通过在弯月面处引入随机波初始扰动来模拟液膜表面的不稳定现象，分析液膜流动特性、Marangoni效应对扰动过程的影响以及不稳定现象的发展机理，以期为提高浸涂工艺中液膜的涂覆质量提供理论基础。

1 理论模型建立

自液池中提拉出的竖直平板表面会附着一层液膜，液膜受重力作用排液，过程如图1所示。图1中，平板悬垂静止，液膜顶端厚度为零，底端与液池相连形成弯月面；x、y、z分别为展向、流向及与平板垂直的方向。假设液膜表面含有不溶性活性剂，假设流体为不可压缩黏性流体，并忽略排液中液膜蒸发的影响。排液过程中，平板左右两侧的液膜对称，因而仅需分析其一侧液膜的变化特征。

图1 竖直平板上液膜排液过程Fig. 1 Film draining on a vertical plate

1.1 控制方程

排液过程的控制方程包括连续方程和运动方程：

式中：含“*”的物理量表示有量纲值；“·”表示2个向量的内积； ∇*为哈密顿算子； ∇*2为拉普拉斯算子；V*=(u*,v*,w*)；j=(0,1,0)；μ*为动力黏度；ρ*为密度；p*为压力；t*为时间；g表示重力加速度。

在z*=0处，边界条件为：

在z*=h*处，运动学条件为：

液膜气液界面满足应力平衡边界条件[7]：

式中：“×”表示2个向量相乘的外积；“：”表示张量的并联双点积；n*、t*分别代表液膜表面的单位外法向和切向矢量；T*为液膜应力张量；H*为曲面平均曲率；σ*为表面张力；κs*、μs*分别为液膜表面的膨胀黏度和剪切黏度为液膜的法向曲率；为表面单位张量。

假设表面张力与活性剂浓度满足线性关系[10]：

在液膜自由面z*=h*(x*,y*,t*)处，活性剂浓度的输运方程为[11]：

1.2 方程无量纲化

方程无量纲变换采用式（9），变换后得到的无量纲物理量不带上标“*”。

式中：ε=d*/l*；d*为液膜特征厚度；l*为动态弯月面长度，l*=(d*D*2)1/3，其中D*为弯月面半径；特征速度U*=ρ*gd*2/μ。S为液膜表面黏度；S1为剪切黏度；M为Marangoni数；Pe为Peclet数，表征活性剂扩散程度。

将式（9）代入式（1）～（8）中进行无量纲化。结合润滑理论（ε<<1），忽略数量级为ε2的小量，可得无量纲后x和y方向上表面速度、液膜厚度和活性剂浓度演化方程为：

2 数值模拟

本模型的求解区域设定为Lx=12π，Ly=12π；计算网格划分为100×100；时间步长设定为0.01。所取参数和范围参照了本课题组前期文献[11]，如表1所示。

表1 典型尺度参数和无量纲参数取值或变化范围Tab. 1 Value or variation range of typical scale parameters and dimensionless parameters

采用“FreeFEM++v4.7”求解方程组（10）～（13）。该软件在高精度数值求解液膜流动问题中具有显著优势[12,13]。

计算中边界条件设置如下：

平板顶端y=0处，厚度为0，无液体流动，即：

底端y=Ly处，液膜与液池相接，形成弯月面。根据 Young-Laplace方程[14]，液膜厚度的一阶和二阶导数满足以下关系：

靠近液池区域的液膜沿y方向无活性剂流动，即：

平板展向两侧采用相同的边界条件：

含活性剂液膜排液过程的数学描述是一个与初始条件相关的初值问题，将扰动施加在液膜厚度和浓度的初始条件中[15]。考虑到实验过程中扰动往往随机产生，并未出现确定的主导波数，故本文参考了文献[16]中描述含活性剂液膜水平铺展过程中不稳定现象所采用的扰动形式。本研究中，为模拟竖直平板上排液过程中液膜表面出现的展向扰动波，在液膜厚度初始条件的弯月面附近yd=35处添加随机波数扰动，如式（18）所示：

活性剂浓度的初始条件如式（19）：

式中：Rand[x,y]表示[x,y]闭区间内的随机数；yd是弯月面位置，取yd=35；y1=36.69；y1

式（19）的含义为：液膜主体部分活性剂浓度为0，自接近弯月面处（y≈33）活性剂浓度急剧增加到1，与液池中浓度相同，且在yd=35处引入随机扰动。

在式（19）中，人为施加了一极大的活性剂浓度梯度，其结果与实际工程中快速提拉平板形成的液膜表面活性剂的浓度分布接近。

3 结果与分析

弯月面附近为液膜最薄区域，当其受到微小扰动后极易出现不稳定特征。为探究局部扰动能否随排液持续而维持不稳定性的问题，在初始条件中施加展向扰动后：需要首先分析液膜厚度和流动特性的变化，以此验证平板液膜二维模型的合理性；然后模拟液膜表面的不稳定特征，揭示局部展向扰动的成因及发展过程。

3.1 M=1时液膜流动特征

当M=1、S=10、S1=1时，模拟所得x=15处的液膜厚度和流向表面速度如图2所示，液膜厚度和活性剂浓度的二维演化过程如图3所示。色条颜色反映了活性剂浓度高低。

图2 M=1时主流方向的液膜厚度及表面速度分布Fig. 2 Film thickness and surface velocity distribution in the mainstream direction when M=1

结合图2（a）和图3知，重力作用促使液膜向下排液，上端液膜厚度逐渐减少，弯月面以下厚度则快速增加。随排液进行，动态弯月面位置上移到y=32.5处，液膜厚度变薄、速率减缓；自t=80开始，厚度几乎不变。在图2（b）中，当t=20时排液速度较初始时刻大幅降低，底端弯月面附近的表面速度首先出现负值，这是因初始条件中的大浓度梯度诱发了 Marangoni效应，从而使液膜逆流；t=40时逆流范围有所扩大，逆流速度提高，同时逆流也携带底部的活性剂向上输运，减弱了流向的Marangoni效应，因而后续几个时刻的表面速度呈现出主流速度和逆流速度同时减小的特点。图 2所示排液过程中，液膜厚度变化及逆流特征与本课题组研究一维平板浸涂过程所得的结果[9]4-6一致。

图3 二维液膜表面的厚度和浓度分布Fig. 3 Thickness and concentration distribution of 2D film surface

由于本文重点关注弯月面附近区域（y值位于30～35）液膜表面的不稳定现象，故截取液膜厚度0

由图3可知，在排液过程中，液膜沿流向变化的同时也沿展向变化。在t=0初始时刻，活性剂集中在液膜底部，在弯月面位置处存在着具有随机扰动波幅大小的扰动波纹；t=20时，逆流促使活性剂分布较均匀，扰动波纹消退，弯月面附近沿展向相对均匀；自t=60时开始，活性剂浓度出现局部峰值，展向的浓度分布不再均匀，弯月面处衍生出微弱的起伏波纹；而到排液末期t=100时，活性剂浓度分布的不均匀度更加显著，弯月面处波纹起伏程度也明显加剧，形成多个凸起结构。

3.2 液膜表面的不稳定特征

为进一步考察液膜表面不稳定现象的演化过程，参考光学干涉法测量液膜厚度的实验原理，将液膜厚度在平板上自上而下的变化用渐变色条表示，色条沿流向的疏密程度反映液膜厚度大小，沿展向的弯曲度表示扰动发展特征，结果如图 4所示。

图4 M=1时液膜表面扰动发展过程Fig. 4 Disturbance evolution of the film surface when M=1

由图4知，在t=0时，与初始条件设置的随机波数扰动相对应，液膜底端弯月面处出现了扰动波纹；但在t=20时，扰动波纹消退，液膜表面在展向上呈现相对平整的状态。排液后期自t=60时开始，液膜表面在y≈27处出现微弱波纹，且随排液进行扰动波纹不断发展；到t=100时，条纹紊乱，扰动波纹起伏剧烈，在y=30、x位于20～30的区域出现了2个液膜厚度的高峰。图中液膜表面的条纹特征反映了排液过程中扰动波纹的发展动态。

在排液初期与后期，液膜表面的扰动特征具有明显的差异性，这表现在：t=0排液初期，在厚度初始条件中施加扰动波纹的波幅大小随机，其影响范围仅限于底部弯月面y>35位置，且扰动迅速消退；而在t>60排液后期，重新出现的扰动波纹波幅较小，但却覆盖了y值在15～30的大范围区域；随着对排液进行扰动的波纹波幅的增大，扰动的影响范围不断扩展；到t=100时，存在扰动的区域基本覆盖了y值位于10～35的绝大部分液膜表面。

局部扰动在排液后期液膜极薄的状态下发展更为充分，并随排液逐渐演化为全局扰动。

3.3 不稳定现象分析

3.3.1 Marangoni效应的影响

在S=10、S1=1时，M分别取为1、10和30，分析 Marangoni效应对液膜排液不稳定现象的影响。

图5示出了不同M值时液膜表面速度演化过程。由图 5（a）所示的x=15处的流向速度vs分布可知：vs分布在初始时刻并无差别；但在t=60时，出现大小不同的逆流区；逆流区随M增大而扩大（如M=30时，对应的逆流区上边界接近y=20；液膜逆流区占据近一半的流动域），逆流速度的最大值也随M增大而提高；在t=100排液末期时，Marangoni效应相对t=60时减弱，导致除弯月面区域外，M=1和M=10所对应的逆流区消失，M=30所对应的逆流区也缩小至 15

图5 不同M时液膜表面速度演化过程Fig. 5 Evolution of film surface velocity at different M

图5（b）所示为y=35处液膜表面展向上的速度分布。

（1）在t=0时，展向速度us为0，液膜受重力主导沿着竖直方向流动。

（2）当t=60时，液膜表面衍生出展向速度。M=1情形下，us≈0.03，且沿展向相对均匀。随M增大，us均匀性变差。当M=30时，us在x值位于10～25范围内出现多个峰谷，且速度在正负数值间变化——这表明沿展向上的活性剂浓度分布不同，液膜受到了不同方向的 Marangoni作用。此时，液膜沿流向的逆流现象依旧显著，见图5（a）；而展向上，因 Marangoni效应，液膜表面也呈现出速度不均匀的特征。

（3）在t=100时，us出现震荡特征，不同x位置处的速度数值差异显著。随M增大，us峰值和均值呈增大趋势，液膜呈现如图4所示的波浪起伏特征。

排液后期，液膜表面逐渐发展的扰动波纹并非只出现在底端弯月面附近，而是覆盖了绝大部分流动域（y>20）。当t=100时，液膜表面的扰动现象尤其明显，如图6所示。

图6 t=100时刻液膜表面不稳定现象Fig. 6 Instability of film surface at t=100

另外，当M值较大时，液膜厚度在流向和展向上均发生较大的起伏变化，扰动波纹尺度增大，不规则性增强。特别是在M=30时，波纹起伏最为突出，这表明液膜处于极不稳定的流动状态。

3.3.2 扰动现象的成因

排液过程初始时刻，弯月面附近区域因施加扰动而出现不稳定波纹，液膜受重力作用竖直向下流动，扰动的展向发展受到阻碍，因而扰动发展微弱甚至消退。随排液进行，活性剂逆流使液膜底端（y>33）处的活性剂浓度相对初期减小，而中上部浓度从零逐渐增加（见图7），因而弯月面附近区域（30

图7 M=30液膜活性剂浓度分布Fig. 7 Distribution of film surfactant concentration when M=30

为进一步量化不稳定现象，定义瞬时扰动能量He如式（20）：

式中：hb为未施加扰动时的基态液膜厚度。

某一时刻液膜表面扰动总能量Eh如式（21）：

图8示出了液膜表面扰动能量Eh随时间的演化过程。

图8 液膜表面扰动能量随时间变化Fig. 8 Variation of film disturbing energy with time

（1）t<20排液前期，人为施加的扰动在弯月面附近（30

（2）到t>60排液后期，弯月面h<0.4h0处，受 Marangoni力作用，弯月面区域整体活性剂浓度低，浓度中施加的扰动影响强烈。该区域扰动波纹随排液进行扩展到大部分液膜区域。整个液膜表面扰动能量大幅增加，如Et=100≈20Et=0。

排液后期，弯月面区域液膜厚度小，且存在低浓度的活性剂分布，施加的展向扰动影响突出。由图8知，随M增大，t>60后液膜扰动能量增长的速率加快；尤其是在M=30时，扰动能量增长速率较M=1时大幅增加。这是由于当M增大时，活性剂沿流向的逆流区域扩大；又因展向扰动作用，活性剂浓度梯度诱发的 Marangoni效应也强烈作用于展向方向，于是：扰动波纹的影响不局限于弯月面位置，而是向整个液膜表面发展，从而导致液膜表面不稳定性增强，整体扰动能量大幅增长。