外加电场下奇粘性对薄膜不稳定性影响的数值模拟

Numerical Simulation of the Effect of Odd-Viscosity on Liquid Film Instability with an External Electric Field

摘要: 本文研究了在外加电场作用下，包含奇粘性的流体薄液膜沿倾斜基底下落时的不稳定效应。其中，奇粘性是由于流体打破常规的时间反衍对称性而产生的。利用线方法数值模拟了Bao和Jian的非线性演化方程。研究结果表明，当奇粘性增大时，薄膜的振幅减小，不稳定性减弱。这说明奇数粘性对下落的薄液膜具有稳定性作用，进一步证实了Bao和Jian的研究结论。此外，对于给定的奇粘性，薄膜的振幅随着电参数的增大而增大，这说明电参数增强了薄膜不稳定性。

Abstract: In this paper, the instability effect of a thin liquid film containing odd viscosity falling along an in-clined substrate under an applied electric field is studied. Among them, the odd viscosity is pro-duced by the fluid breaking the conventional time inverse diffraction symmetry. Nonlinear evolu-tion equation of Bao and Jian was numerically simulated by using lines method. The results show that when the odd viscosity increases, the amplitude of the film decreases and the instability weakens. This shows that the odd viscosity has a stabilizing effect on the falling thin liquid film, further confirming the research conclusion of Bao and Jian. In addition, for a given odd-viscosity, the amplitude of liquid film increases with the increase of the electrical parameter, indicating that the electrical parameter enhances the liquid film instability.

1. 引言

薄膜流动不稳定性问题是流体力学中典型的界面不稳定问题。在天体物理学 [1]、海洋学 [2] 和惯性约束聚变 [3] 等方面广泛存在，近年来人们对薄膜流动不稳定的研究日益关注。在涂层工艺中 [4]，这种不稳定性会导致涂层不均匀或液滴滴落。因此控制和抑制薄膜失稳对许多工程应用是不可或缺的，目前已有几种方法来实现这一目标，如外部磁场的存在 [5]、高频机械振荡 [6] 和利用几何特征 [7]。此外，利用外加电场在薄膜界面处产生额外应力去保持应力平衡也是一种被普遍认可的影响薄膜流动稳定性方法。Gonzalez等人 [8] 研究了法向电场作用下薄膜沿斜板的流动，得到了自由表面非线性演化方程。Tseluiko等人 [9] 研究了电场对下落薄膜的影响，发现电场对薄膜稳定性具有失稳作用。2018年，Tudball [10] 考虑了地形影响的情况下研究了法向电场影响下薄膜沿斜面流动的稳定性。

一般的薄膜不稳定性的研究忽略了奇粘性的影响，只考虑了均匀粘度的流体。Avron [11] 突破性地发现奇粘性是由于经典流体打破常规的时间反衍对称性而产生的。为了得到更准确的稳定性分析结果，我们不能忽略奇粘性的影响。2020年，Zhao和Jian [12] 研究了奇粘性在磁场下对薄膜不稳定性的影响，发现奇粘性和磁场稳定了流动。此外，Zhao和Jian [13] 还研究了奇粘性对沿斜壁流动粘弹性薄膜不稳定性的影响，发现奇粘性使流动稳定，而粘弹性使流动失稳。

Bao和Jian [14] 研究了奇粘性在电场下对沿斜面流动的薄膜不稳定的影响，通过线性和弱非线性稳定性分析得出奇粘性具有稳定作用的结论。上述稳定性分析中，仅仅通过线性和弱非线性进行稳定性分析，目前还未利用数值模拟去加强结论。在本文中，我们将在外加电场的情况下，利用线方法数值模拟Bao和Jian [14] 的非线性演化方程，研究奇粘度对薄膜失稳的影响。

2. 数学模型

2.1. 数学模型的描述

图形示意图参照Bao和Jian [14] 中的图1，我们考虑密度为ρ的不可压缩粘性液膜沿倾斜角为β的平行板的流动行为。建立二维笛卡尔坐标系(x, y)来描述流动过程，其中x和y分别沿平行板方向和垂直于平行板方向。h(x, t)表示液膜瞬时位置，未扰动高度为h₀。下板电势为φ₁ = φ_c (φ_c是常数)。上板接地，即φ₂ = 0。液体–介质界面的表面张力系数为γ，重力加速度为g。区域1表示理想导体的流体层，这表明液体薄膜界面处和液体内部的电势均为φ_c。由于电导率很小，我们忽略了感应磁场，因此宜采用麦克斯韦方程的静电近似法。区域2表示介电常数为ε₀的理想电介质。对于该双层流体，我们可以考虑为水–空气系统。外加法向电场的会产生额外的麦克斯韦应力，从而影响薄膜的稳定性。因此在时间反演对称性被打破的情况下，液体薄膜下落问题应同时考虑偶应力张量τ^e和奇应力张量τ^o。此时应力张量可表示为 $\tau ={\tau }^e+{\tau }^o$ [11]，具体形式为 [11] [15]

${\tau }_{ij}^e=-p{\delta }_{ij}+{\eta }^e\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$, (1)

${\tau }_{ij}^o=-{\eta }^o\left({\delta }_{i1}{\delta }_{j1}-{\delta }_{i2}{\delta }_{j2}\right)\left(\frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1}\right)+{\eta }^o\left({\delta }_{i1}{\delta }_{j2}+{\delta }_{i2}{\delta }_{j1}\right)\left(\frac{\partial u_1}{\partial x_1}-\frac{\partial u_2}{\partial x_2}\right)$,(2)

其中下标i, j = 1, 2，η^e和η^o分别为偶数和奇数粘度系数。

设(x, y)坐标系下的速度场为u = (u, v)，区域1的控制方程为具有奇粘性影响的不可压缩Navier-Stokes方程

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$, (3)

$\rho \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\frac{\partial p}{\partial x}+{\eta }^e\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)-{\eta }^o\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+\rho g\sin \beta$, (4)

$\rho \left(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)=-\frac{\partial p}{\partial y}+{\eta }^e\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+{\eta }^o\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)-\rho g\cos \beta$, (5)

其中p表示压强。对于区域2，有如下电场E₂与电势φ₂的关系式

$E_2=-\nabla {\phi }_2$, (6)

并且电势φ₂满足高斯定律

${\phi }_{2xx}+{\phi }_{2yy}=0$. (7)

在倾斜底板y = 0处施加无滑移、无渗透边界条件为

$u=v=0$. (8)

自由面 $y=h\left(x,t\right)$ 处的运动学条件为

$v=h_t+u{h}_x$.(9)

由于自由面上 $y=h\left(x,t\right)$ 的电势是连续的，上板是接地的，有条件

${{\phi }_2|}_{y=h}={{\phi }_1|}_{y=h}={\phi }_c$, (10a)

${{\phi }_2|}_{y=d}=0$. (10b)

考虑其切向应力平衡得到

$2\left(v_y-u_x\right)h_x+\left(1-h_x^2\right)\left(v_x+u_y\right)+\eta \left[\left(1-h_x^2\right)\left(u_x-v_y\right)+2h_x\left(v_x+u_y\right)\right]=0$ 在 $y=h\left(x,t\right)$,(11)

式中η为奇粘度系数η^o与偶粘度系η^e之比，即η = η^o/η^e。此外，得到法向应力平衡

$\begin{array}{l}{p}_g-p+\frac{{\eta }^e}{1+h_x^2}\left\{2v_y-2h_x\left(v_x+u_y\right)+2h_x^2{u}_x+\eta \left[\left(1-h_x^2\right)\left(v_x+u_y\right)-2h_x\left(u_x-v_y\right)\right]\right\}\\ =\frac{\gamma h_{xx}}{{\left(1+h_x^2\right)}^{3/2}}+\frac{{\epsilon }_0}{1+h_x^2}{\left(h_x{\phi }_{2x}-{\phi }_{2y}\right)}^2-\frac{{\epsilon }_0}{2}\left({\phi }_{2x}^2+{\phi }_{2y}^2\right)\end{array}$ (12)

其中p_g为大气压强(区域2的压强)。

为使方程无量纲化，引入以下由上标星号表示的无量纲变量

$x^{\ast }=\frac{x}{h_0}$, $y^{\ast }=\frac{y}{h_0}$, $t^{\ast }=\frac{t}{h_0/u_0}$, $h^{\ast }=\frac{h}{h_0}$, $d^{\ast }=\frac{d}{h_0}$, $u^{\ast }=\frac{u}{u_0}$, $v^{\ast }=\frac{v}{u_0}$, $p^{\ast }=\frac{p}{{\eta }^e{u}_0/h_0}$,

${\phi }^{\ast }=\frac{\phi }{{\phi }_c}$, $u_0=\rho g{h}_0^2\sin \beta /2{\eta }^e$. (13)

为方便起见，在之后的表达式中去掉星号。将无量纲变量代入控制方程和边界条件得到方程参见Bao和Jian [14]。其他无量纲参数雷诺数Re、毛细管数C (测量粘滞力与毛细管力的比值)和电参数E (测量电与引力的比值)，形式分别为

$\mathrm{Re}=\frac{\rho u_0{h}_0}{{\eta }^e}$, $C=\frac{{\eta }^e{u}_0}{\gamma }=\frac{{\mathrm{Re}}^{2/3}{\sin }^{1/3}\beta }{2^{1/3}K}$, $E=\frac{{\phi }_c^2{\epsilon }_0}{2{\eta }^e{u}_0{h}_0}$, $K=\frac{\gamma {\rho }^{1/3}}{g^{1/3}{\eta }^e{}^{4/3}}$. (14)

这里Kapitza数 [10] 是表面张力与惯性力的无量纲之比。对于给定的材料，K在给定的温度下是常数，对于室温下的水，Kapitza数约为3364.5。

2.2. 长波演化方程

首先假设薄膜界面的变形波长λ远大于未扰动厚度h₀，即h₀/λ = δ，其中δ是一个小参数，它通常被称为薄膜参数。在液体层(区域1)中存在有明显差异的长度尺度，引入以下变量 [10]

$\xi =\delta x$, $\tau =\delta t$, $y=y$, $v=\delta w$.

为了在产生的模型中确保表面的影响，重新调整毛细管数C，即 $C={\delta }^2{C}^{\prime }$。将这些新的重新调整的变量代入无量纲方程和边界条件得到含δ的方程，之后通过下列渐近展开可以推导主阶方程和一阶方程

$u=u_0+\delta u_1+\cdots$, $w=w_0+\delta w_1+\cdots$, $p=p_0+\delta p_1+\cdots$, ${\phi }_2={\phi }_{20}+\delta {\phi }_{21}+\cdots$.

忽略Ο(δ)项和更高阶项，直接将展开后的变量代入模型的无量纲方程，推导主阶方程并求解得到主阶解

$u_0=-y^2+2hy$, (15a)

$w_0=-y^2{h}_{\xi }$,(15b)

$p_0=-2\left(\eta +\cot \beta \right)\left(y-h\right)+{\bar{p}}_g-\frac{h_{\xi \xi }}{C^{\prime }}+\frac{E}{{\left(h-d\right)}^2}$.(15c)

同理得到一阶方程并求解得到一阶解

$\begin{array}{c}{u}_1=\mathrm{Re}\left(\frac{1}{3}{y}^3{h}_{\tau }+\frac{1}{6}{y}^{4}h{h}_{\xi }\right)+\cot \beta h_{\xi }{y}^2-\frac{h_{\xi \xi \xi }}{2C^{\prime }}{y}^2+\frac{E{h}_{\xi }}{{\left(h-d\right)}^3}-\mathrm{Re}\left(h^2{h}_{\tau }+\frac{2}{3}{h}^4{h}_{\xi }\right)y\\ -2\cot \beta h{h}_{\xi }y+\frac{h_{\xi \xi \xi }}{C^{\prime }}hy-\frac{2Eh{h}_{\xi }}{{\left(h-d\right)}^3}y-4\eta h{h}_{\xi }y\end{array}$ (16)

流向流量 $q\left(\xi ,\tau \right)$ 是通过如下关系式定义的

$q\left(\xi ,\tau \right)={\displaystyle {\int }_0^{h\left(\xi ,\tau \right)}u\left(\xi ,y,\tau \right)\text{d}y}$, (17)

其中 $u\left(\xi ,y,\tau \right)=u_0\left(\xi ,y,\tau \right)+\delta u_1\left(\xi ,y,\tau \right)+{\rm O}\left({\delta }^2\right)$。运动学边界条件的另一种形式如下

$h_{\tau }+q_{\xi }=0$, (18)

将式(17)代入式(18)后利用变换 $\delta \left({\partial }_{\tau },{\partial }_{\xi }\right)=\left({\partial }_t,{\partial }_x\right)$ 得到薄膜非线性演化方程

$h_t+\frac{\partial }{\partial x}\left[\frac{2}{3}{h}^3+\frac{8}{15}{h}^6{h}_x\text{Re}-\frac{2}{3}\left(3\eta +\cot \beta +\frac{E}{{\left(h-d\right)}^3}\right)h^3{h}_x+\frac{1}{3C}{h}^3{h}_{xxx}\right]=0$. (19)

式(19)是Bao和Jian [14] 推导出的结果。当忽略电场(E = 0)、奇粘性(η = 0)和地形(β = 0)的影响，式(19)与Tseluiko等人 [16] 得到的式(2.2)一致。当忽略地形(β = 0)和电场(E = 0)的影响时，其结果与Tudball [10] 的方程一致。

3. 数值模拟

由于演化方程难以直接求解，因此采用线方法 [17] 来模拟液膜的时空演化。演化方程可以改写为

$h_t+q_x=0$, (20a)

其中

$q=\frac{2}{3}{h}^3+\frac{8}{15}{h}^6{h}_x\mathrm{Re}-\frac{2}{3}\left(3\eta +\cot \beta +\frac{E}{{\left(h-d\right)}^3}\right)h^3{h}_x+\frac{1}{3C}{h}^3{h}_{xxx}$, (20b)

边界条件可以表示为：

$h\left(L,t\right)=0.1$, $h_x\left(L,t\right)=0$, $h_x\left(0,t\right)=0$, $h_{xxx}\left(0,t\right)=0$,(21)

其中L为域的无量纲长度。初始条件可以表示为：

$h\left(x,0\right)=1+0.05\cos \left(kx\right)$, (22)

构造一个空间网格 $x_i=\left(i-1\right)\Delta x$， $i=1,2,\cdots ,N$，其中Δx = L/N。为了满足质量守恒条件，应在i ± 1/2处的离散化q，因此式(20a)离散形式为

$h_{it}+\frac{q_{i+\frac{1}{2}}-q_{i-\frac{1}{2}}}{\Delta x}=0$, (23a)

当i在3和 $N-\text{2}$ 之间时，

$\begin{array}{c}{q}_{i+\frac{1}{2}}=\frac{2}{3}\frac{h_{i+1}^3+h_i^3}{2}+\frac{8}{15}\mathrm{Re}\frac{h_{i+1}^6+h_i^6}{2}\frac{h_{i+1}-h_i}{\Delta x}-\frac{2}{3}\left(3\eta +\cot \beta +\frac{E}{{\left(h-d\right)}^3}\right)\frac{h_{i+1}^3+h_i^3}{2}\frac{h_{i+1}-h_i}{\Delta x}\\ +\frac{1}{3C}\frac{h_{i+1}^3+h_i^3}{2}\frac{h_{i+2}-3h_{i+1}+3h_i-h_{i-1}}{\Delta x^3}\end{array}$ (23b)

$\begin{array}{c}{q}_{i-\frac{1}{2}}=\frac{2}{3}\frac{h_i^3+h_{i-1}^3}{2}+\frac{8}{15}\mathrm{Re}\frac{h_i^6+h_{i-1}^6}{2}\frac{h_i-h_{i-1}}{\Delta x}-\frac{2}{3}\left(3\eta +\cot \beta +\frac{E}{{\left(h-d\right)}^3}\right)\frac{h_i^3+h_{i-1}^3}{2}\frac{h_i-h_{i-1}}{\Delta x}\\ +\frac{1}{3C}\frac{h_i^3+h_{i-1}^3}{2}\frac{h_{i+1}-3h_i+3h_{i-1}-h_{i-2}}{\Delta x^3}\end{array}$ (23c)

由于i = N + 1对应x = L，根据边界条件(21)，该端点的厚度值为0.1，故 $h_N{}_{+\text{1}}=0.\text{1}$。用域外的点(鬼点)来表示边界附近的导数得到 $h_N{}_{+\text{2}}=h_N$，对称性条件得到h₀ = h₂，h₋₁ = h₃，因此可以用标准微分方程求解器求解。在实际工程问题中，首先要确定参数范围。我们根据文献 [14] 给出以下参考范围：ρ ~ 1000 kg m⁻³，μ ~ 0.9 × 10⁻³ Pa∙s，h₀ ~ 3.6 × 10⁻³ m，β ~ π/4，K ~ 3364.5，ε₀ ~ 10⁻¹¹ faradm⁻¹，g ~ 9.8 m∙s⁻²。

图1表示了奇粘性系数η对自由面时空演化的影响。如图显示了t = 5到t = 80时间间隔为5的情况。对于一个固定的η，自由面振幅随着时间t的推移而增加，这表明了薄膜的不稳定。观察图1(a)~(d)，波的振幅随着η的增加而减小，这说明η具有稳定薄膜的作用。此外，随着奇数粘度的增加，振幅随着时间的变化速率变慢，这证明奇数粘度降低了振幅的增长率使得液膜变得更加稳定。

Figure 1. Film thicknesses at various times (Re = 1, β = π/4, d = 3, E = 1) (a) η = 0, (b) η = 0.5, (c) η = 1, (d) η = 1.5

图1. 薄膜厚度在不同时期(Re = 1，β = π/4，d = 3，E = 1) (a) η = 0，(b) η = 0.5，(c) η = 1，(d) η = 1.5

图2描述了在发展时间为t = 80的情况下，不同奇粘性和电参数下薄膜厚度与水平距离的非线性演化。如图，振幅随着奇粘性系数η的增大而减小，这验证了η具有抑制不稳定性的作用。对比图2(a)~(d)，我们观察到随着电参数E的增加薄膜的振幅增加，这说明E增强了薄膜不稳定性。在不同的电参数E的情况下，奇粘性都表现出抑制薄液膜不稳定性的作用。

为了更好的体现奇粘性的影响，图3显示了不同奇粘性系数η下自由表面的时空演化特征。我们清楚地观察到，当η增加时，薄膜的振幅减小，表明η对薄膜的不稳定性有稳定作用。

Figure 2. Film thicknesses at different odd viscosity coefficients η of 0, 0.5, 1, 1.5 (Re = 1, β = π/4, d = 3, t = 80): (a) E = 0, (b) E = 0.5, (c) E = 1, (d) E = 5

图2. 不同奇粘性η为0，0.5，1，1.5 (Re = 1, β = π/4, d = 3, t = 80)：(a) E = 0，(b) E = 0.5，(c) E = 1，(d) E = 5时的薄膜厚度

(a)(b)(c)(d)

Figure 3. Three-dimensional models of film thickness evolution with t and x (Re = 1, β = π/4, d = 3, E = 5): (a) η = 0, (b) η = 0.5, (c) η = 1, (d) η = 1.5

图3. 薄膜厚度随t和x变化的三维模型(Re = 1，β = π/4，d = 3，E = 5)：(a) η = 0，(b) η = 0.5，(c) η = 1，(d) η = 1.5

4. 结论

本研究系统地研究了在外加电场作用下，奇粘度对倾斜基底下薄膜不稳定性的影响。利用线方法数值模拟了Bao和Jian [14] 的非线性演化方程，成功地刻画了自由表面流动的演化过程，并证实了增加奇粘性η可以提高薄膜体系的稳定性。在其他参数固定的情况下，对于给定的奇粘性，薄膜的振幅随着电参数E的增加而增加，这说明E增强了薄膜不稳定性。

基金项目

本研究由国家自然科学基金地区项目(12262026)、内蒙古自然基金面上项目(2021MS01007)、内蒙古自治区党委组织部：2019年“草原英才”滚动支持项目(12000-12102408)和内蒙古自治区高等学校创新团队发展计划支持(NMGIRT2323)资助。

NOTES

^*通讯作者Email: [e-mail]