剪切作用与均匀加热电场共同驱动下倾斜平板奇粘度流体薄膜的线性稳定性分析

Linear Stability Analysis of Odd-Viscous Fluid Films on an Inclined Plane Driven by Combined Shear Action and Uniformly Heated Electric Field

摘要: 倾斜平板液膜流动的稳定性调控是多相流与界面流体力学的核心研究方向，广泛应用于微纳涂层、工业换热等工程场景。本研究聚焦剪切–均匀热毛细–电场–奇粘度四场耦合的倾斜平板液膜流动，通过构建耦合数学模型、长波近似简化、线性与弱非线性稳定性分析及有限差分数值模拟，结合自由表面时空演化观测，系统探究多场协同调控机制。研究构建含奇粘度非对称应力张量的自由表面演化方程，明确马兰戈尼数Ma、电参数E及剪切应力τ的增大显著加剧流动失稳，奇粘度μ₁通过修正应力张量实现稳定化；确定热毛细力与电场麦克斯韦应力相互抵消的临界条件，流动达到动态平衡；发现奇粘度与剪切应力的耦合效应具有波数依赖性，波数k < P₁时剪切增强稳定性，k > P₁时剪切加剧失稳，且μ₁ = 1时可放大τ对自由表面演化的调控差异。通过弱非线性分析划分四类失稳区域，量化各参数对临界振幅ζa的影响规律。本研究建立的四场耦合模型与揭示的调控机制，为液膜流动精准控制提供理论支撑。

Abstract: The stability control of liquid film flow on an inclined plate is a core research direction in multiphase flow and interfacial fluid mechanics, which is widely applied in engineering scenarios such as micro-nano coating and industrial heat transfer. This study focuses on liquid film flow on an inclined plate coupled with four fields: shear, uniform thermocapillary, electric field, and odd viscosity. Through the establishment of a coupled mathematical model, long-wave approximation simplification, linear and weakly nonlinear stability analysis, finite difference numerical simulation, combined with the observation of space-time evolution of the free surface, the synergistic regulation mechanism of multiple fields is systematically explored. The study constructs a free surface evolution equation containing an asymmetric stress tensor of odd viscosity, and clarifies that the increase of Marangoni number Ma, electric parameter E, and shear stress τ significantly intensifies flow instability, while odd viscosity μ₁ achieves stabilization by modifying the stress tensor; the critical condition is determined, at which point the thermocapillary force and the Maxwell stress of the electric field cancel each other out, and the flow reaches a dynamic equilibrium; it is found that the coupling effect between odd viscosity and shear stress has a wave number dependence: when the wave number k is less than the wave number for P₁, shear enhances stability; when k is greater than the wave number for P₁, shear intensifies instability. In addition, when μ₁ = 1, it can significantly amplify the regulation difference of τ On evolution of the free surface. Four types of instability regions are divided through weakly nonlinear analysis, and the influence of each parameter on critical amplitude ζa is quantified. The four-field coupling model and the revealed regulation mechanism established in this study provide theoretical support for the precise control of liquid film flow.

1. 引言

倾斜平板液膜流动是多相流与界面流体力学的经典研究方向，广泛应用于微纳涂层、工业换热、薄膜制备等工程场景，液膜流动稳定性直接决定工艺实施效率与产品质量，因而成为流体力学领域长期关注的研究热点[1] [2]。剪切作用、热毛细效应、外加电场等多物理场耦合下的液膜失稳机制，因涉及复杂的界面动力学行为与动量输运过程，是研究的一大重点[3] [4]。早期围绕常规粘性液膜的剪切耦合失稳研究已形成系统理论体系，而具有时间反演对称性破缺特性的奇粘度流体被引入液膜流动研究后，为液膜稳定性调控开辟了全新路径，相关研究也逐步成为界面流体力学与非牛顿流体力学的交叉研究前沿[5] [6]。

针对剪切作用下的液膜流动稳定性，国内外学者已开展大量基础性研究。Samanta [1]针对剪切施加下的下降液膜，深入分析其最优扰动增长规律与临界雷诺数特性，明确了剪切力对液膜流动稳定性的核心调控机制；Hossain等人[3]结合基底起伏、热毛细效应等复杂工程工况，探究了剪切耦合下液膜的水动力失稳特征。随着非牛顿流体研究的不断深入，奇粘度效应被引入剪切耦合液膜体系，Desai等人[7] [8]将奇粘度与剪切、热毛细效应、基底滑移相结合，研究发现奇粘度可通过修正动量输运过程显著改变液膜失稳阈值，证实了奇粘度在剪切耦合液膜失稳调控中的关键作用。

奇粘度是一类打破经典流体本构关系的特殊粘性系数，由Avron [5]首次从理论层面阐明其物理本质，与传统偶粘度的能量耗散特性不同，奇粘度源于流体微观或介观尺度粒子的手性运动与旋转特性，核心通过修正应力张量结构改变流体的动量输运规律[6]。经典粘性流体的应力张量满足对称特性，而奇粘度的引入会使应力张量产生非对称分量，形成奇应力张量与偶应力张量耦合的本构形式，在不可压缩流体体系中具体表现为对速度梯度项的非对称修正，进而改变流场压力梯度分布与界面演化规律[9] [10]。这种非对称修正的本质，是奇粘度为流体微团提供了定向横向约束力，该力垂直于剪切、热毛细力的作用方向，可有效阻碍微团的拉伸或偏移变形，削弱局部受力不均，这也是奇粘度实现宏观稳定化的核心物理根源[5] [8]。在液膜流动体系中，奇粘度并非单纯的粘性耗散项，而是通过调控界面剪切应力分布、抑制扰动波增长幅值实现对液膜流动的稳定化作用[7]，这一特性已在电场、热毛细效应等单物理场耦合的液膜体系中得到实证验证，也为宏观液膜流动的非接触式精准调控提供了新途径[6] [11] [12]。

在热场与奇粘度耦合的液膜失稳研究方面，Mukhopadhyay等人[13]针对线性温度分布下的倾斜平板液膜，探究了奇粘度对热毛细失稳与波系演化的影响，发现奇粘度可有效抑制热毛细效应诱导的界面扰动增长；Jia & Jian [14]揭示了奇粘度对加热倾斜基底上液膜瑞利–泰勒失稳的调控机制，丰富了热场耦合下奇粘度流体液膜的失稳理论体系。在外加电场与奇粘度耦合的液膜研究方面，Papageorgiou [15]系统综述了电场作用下液膜流动的研究进展，明确了麦克斯韦应力对液膜界面形态与流动稳定性的调控规律；Bao & Jian [9]首次将奇粘度与外电场耦合引入液膜流动体系，发现电场力与奇粘度的协同作用可实现对液膜失稳的精准调控，为多物理场耦合下的液膜调控提供了新思路。

本研究聚焦剪切作用、外加电场、热毛细效应与奇粘度效应耦合下的倾斜平板液膜流动，通过构建多物理场耦合数学模型、细化演化方程求解逻辑，采用线性小扰动法与多尺度方法，分别开展液膜流动的线性稳定性与弱非线性稳定性分析，揭示多场耦合下液膜的失稳规律与非线性失稳区域特征；通过有限差分法开展数值模拟研究，验证理论分析结果的可靠性，系统揭示多场协同对液膜流动稳定性的调控机制，为微纳涂层、换热系统等工程领域的流动控制与性能优化提供理论支撑与方法参考。

2. 数学模型

如图1所示，本文旨在研究剪切作用与外加电场共同作用下，倾斜均匀加热平板上具有奇粘度效应的流体薄膜的自由流动特性。

Figure 1. Schematic diagram of the flow model for odd-viscous fluid films on an inclined uniformly heated plane under shear action and external electric field

图1. 剪切作用与外加电场下，奇粘度流体薄膜在倾斜均匀加热平板上的流动模型示意图

该模型中，基底为刚性无渗透倾斜角为β的平板，平板长度为L，g为重力加速度，坐标系选取笛卡尔坐标系$\left(x,z\right)$ ，其中x轴沿倾斜平板方向延伸，z轴垂直于基底向上，u和w分别表示对应的速度分量。流体薄膜的自由表面由$z=h\left(x,t\right)$ 描述，该表面受到切向外部剪切应力${\tau }_a$ 作用，当${\tau }_a>0$ 时，剪切应力沿下游方向施加；当${\tau }_a<0$ 时，剪切应力沿上游方向作用。流体的物理性质参数设定如下：密度为ρ、表面张力为σ、热扩散率为$\kappa ={\kappa }_{th}/\left(\rho c_p\right)$ (其中${\kappa }_{th}$ 为流体的热导率，$c_p$ 为流体在恒定压力下的比热容)。基底维持恒定温度$T=T_w$ ，液膜被均匀加热后与周围空气形成恒定温度差$\Delta T=\left(T_a-T_{\infty }\right)>0$ ($T_{\infty }$ 为环境空气温度)；本研究假设流体无挥发性，且忽略浮力与蒸发效应的影响。此外，流动过程中施加外加电场，其中区域1 (流体层)视为理想导体，因此流体内部及薄膜界面处的电势$\phi$ 保持恒定值${\phi }_c$ ；采用静电近似处理电场问题(忽略感应磁场)，区域2 (空气层)视为理想电介质，其介电常数为${\epsilon }_0$ ，上平板接地且电势为0，该双层结构可模拟水–空气系统。法向电场的引入会产生额外的麦克斯韦应力，进而对流体流动稳定性产生影响。

对于均匀理想的牛顿流体，应力张量${\tau }_{ij}$ 和应变率张量$u_{kl}$ 满足关系${\tau }_{ij}={\mu }_{ijkl}{u}_{kl}$ ，其中${\mu }_{ijkl}$ 是粘性张量，$u_i$ 表示流体速度分量，$u_{kl}=1/2\left(\partial u_l/\partial x_k+\partial u_k/\partial x_l\right)$ 。当时间反演对称性成立时，应力张量呈对称形式；而当时间反演对称性被破坏时，粘性张量可分解为对称的偶粘度分量与无耗散的反对称奇粘度分量[5]。本研究中，流体薄膜的时间反演对称性被破坏，柯西应力张量${\tau }_{ij}$ 分解为偶部${\tau }_{ij}^e$ 与奇部${\tau }_{ij}^o$ ，即${\tau }_{ij}={\tau }_{ij}^e+{\tau }_{ij}^o$ [5]-[10]。其中，偶部满足${\tau }_{ij}^e=-p{\delta }_{ij}+{\mu }^e\left(u_{ij}+u_{ji}\right)$ ，$i,j=1,2$ ，${\mu }^e$ 为标准偶粘度系数，${\delta }_{ij}$ 为克罗内克符号，p为流体压力；奇部满足${\tau }_{ij}^o=-2{\mu }^o\left({\delta }_i{}_1{\delta }_j{}_1-{\delta }_i{}_2{\delta }_j{}_2\right)u_{12}+{\mu }^e\left({\delta }_i{}_1{\delta }_j{}_2-{\delta }_i{}_2{\delta }_j{}_1\right)\left(u_{11}-u_{22}\right)$ ，$i,j=1,2$ ，${\mu }^o$ 为奇粘度系数。其中奇部应力张量的具体分量为：${\tau }_{xx}^o=-{\mu }^o\left(u_z+u_x\right),{\tau }_{xz}^o={\tau }_{zx}^o={\mu }^o\left(u_x-u_z\right),{\tau }_{zz}^o={\mu }^o\left(u_z+u_x\right)$ 。后续统一采用$u_1=u$ ，$u_2=w$ ，$x_1=x$ 和$x_2=z$ 的符号表示。

2.1. 控制方程

基于流体力学基本原理建立如下的控制方程

$u_x+w_z=0$ ， (1)

$\rho \left(u_t+u{u}_x+w{u}_z\right)=-p_x+\rho g\text{sin}\beta +{\mu }^e\left(u_{xx}+u_{zz}\right)-{\mu }^o\left(w_{xx}+w_{zz}\right)$ ， (2)

$\rho \left(w_t+u{w}_x+w{w}_z\right)=-p_z-\rho g\text{cos}\beta +{\mu }^e\left(w_{xx}+w_{zz}\right)+{\mu }^o\left(u_{xx}+u_{zz}\right)$ ， (3)

$T_t+u{T}_x+w{T}_z=\kappa \left(T_{xx}+T_{zz}\right)$ . (4)

对于区域2 (空气层)，电场满足高斯定律，其控制方程为：

${\phi }_{xx}+{\phi }_{zz}=0$ ， (5)

式中，表面张力σ与温度的关系满足$\sigma \left(T\right)={\sigma }_{\infty }-{\alpha }_T\left(T-T_{\infty }\right)$ 为参考温度下的表面张力，${\alpha }_T$ 为表面张力温度系数)。

2.2. 边界条件

在流体与基底的边界处(即$z=0$ 处)，由于基底表面光滑，我们考虑无滑移边界条件以及无渗透边界条件，即$u=0$ ，$w=0$ 。此外假设基底温度保持在恒定温度$T=T_a$ ，电势满足$\phi ={\phi }_a$ 。

在自由表面$z=h\left(x,t\right)$ 处，切向满足应力平衡$\overrightarrow{n}\cdot \tau \cdot \overrightarrow{t}={\nabla }_s\sigma \cdot \overrightarrow{t}+{\tau }_a\cdot \overrightarrow{t}$ ，法向满足应力平衡$p_{\infty }+\overrightarrow{n}\cdot \tau \cdot \overrightarrow{n}=-\sigma \nabla \cdot \overrightarrow{n}+\overrightarrow{n}\cdot M\cdot \overrightarrow{n}$ 。其中$M={\epsilon }_{0}EE-{\epsilon }_0/2{\left|E\right|}^2\delta$ ，为空气的Maxwell应力张量，δ是克罗内克张量。还须符合运动学条件、牛顿冷却定律以及电势的连续性，具体边界条件如下：

$\begin{array}{l}{\mu }^e\left[\left(u_z+w_x\right)\left(1-h_x^2\right)-2\left(u_x-w_z\right)h_x\right]+{\mu }^o\left[\left(u_x-w_z\right)\left(1-h_x^2\right)+2\left(u_z+w_x\right)h_x\right]\\ =\text{}\left({\sigma }_x+h_x{\sigma }_z\right){\left(1+h_x^2\right)}^{1/2}+{\tau }_a{\left(1+h_x^2\right)}^{1/2}\end{array}$ ， (6)

$\begin{array}{l}{p}_a-p+\left\{2{\mu }^e\left[u_x{h}_x^2-\left(u_z+w_x\right)h_x+w_z\right]+{\mu }^o\left[\left(u_z+w_x\right)\left(1-h_x^2\right)-2\left(u_x-w_z\right)h_x\right]\right\}{\left(1+h_x^2\right)}^{-1}\\ =\sigma \left(T\right)h_{xx}{\left(1+h_x^2\right)}^{-3/2}+{\epsilon }_0{\left(h_x{\phi }_x-{\phi }_z\right)}^2{\left(1+h_x^2\right)}^{-1}-\frac{{\epsilon }_0}{2}\left({\phi }_x^2+{\phi }_z^2\right)\end{array}$ ， (7)

$w=h_t+u{h}_x$ ， (8)

${\kappa }_{th}\left(T_z-h_x{T}_x\right){\left(1+h_x^2\right)}^{-1/2}+{\kappa }_g\left(T-T_{\infty }\right)=0$ ， (9)

${\phi |}_{\text{z}=h}={\phi }_c$ ， (10)

由于倾斜顶板接地，在$z=d$ 处电势$\phi =0$ 。上述定义中${\kappa }_g$ 为液膜与环境空气的热传递系数，$p_a$ 为大气压力。

2.3. 无量纲

为消除物理量单位对分析结果的影响，选取特征尺度对各物理量进行无量纲化处理：

$U_0=\frac{H^{2}g\text{sin}\beta }{\nu }$ ，$\nu =\frac{{\mu }^e}{\rho }$ ，$\left(x,z,h,\delta \right)=H\left(x^{*},z^{*},h^{*},{\delta }^{*}\right)$ ，$\left(u,w\right)=U_0\left(u^{*},w^{*}\right)$ ，$t=\frac{H}{U_0}{t}^{*}$ ，

$p-p_{\infty }=\frac{U_0}{H}{\mu }^e{p}^{*}$ ，$T=T_{\infty }+T^{*}\text{Δ}T$ ，$\phi ={\phi }^{*}{\phi }_c$ ， (11)

通过使用无量纲变量，方程(1)~(5)化为如下形式

$u_x+w_z=0$ ， (12)

$\text{Re}\left(u_t+u{u}_x+w{u}_z\right)=-p_x+1+u_{xx}+u_{zz}-{\mu }_1\left(w_{xx}+w_{zz}\right)$ ， (13)

$\text{Re}\left(w_t+u{w}_x+w{w}_z\right)=-p_z-\text{cos}\beta +w_{xx}+w_{zz}+{\mu }_1\left(u_{xx}+u_{zz}\right)$ ， (14)

$Pe\left(T_t+u{T}_x+w{T}_z\right)=T_{xx}+T_{zz}$ ， (15)

区域2 (空气层)的电场控制方程无量纲形式为：

${\phi }_{xx}+{\phi }_{zz}=0$ ， (16)

相应地，无量纲边界条件如下：

在$z=0$ 平面，满足$u=0$ ，$w=0$ ，$T=1$ ， (17)

在$z=h\left(x,t\right)$ 处，满足无量纲后的边界条件有电势$\phi =1$ ，

$\begin{array}{l}\left[\left(u_z+w_x\right)\left(1-h_x^2\right)-2\left(u_x-w_z\right)h_x\right]+{\mu }_1\left[\left(u_x-w_z\right)\left(1-h_x^2\right)+2\left(u_z+w_x\right)h_x\right]\\ =-Ma\left(T_x+h_x{T}_z\right){\left(1+h_x^2\right)}^{1/2}+\tau {\left(1+h_x^2\right)}^{1/2}\end{array}$ ， (18)

$\begin{array}{l}p=\left\{2\left[u_x{h}_x^2-\left(u_z+w_x\right)h_x+w_z\right]+{\mu }_1\left[\left(u_z+w_x\right)\left(1-h_x^2\right)-2\left(u_x-w_z\right)h_x\right]\right\}{\left(1+h_x^2\right)}^{-1}\\ -\left(We-MaT\right)h_{xx}{\left(1+h_x^2\right)}^{-3/2}-2E{\left(h_x{\phi }_x-{\phi }_z\right)}^2{\left(1+h_x^2\right)}^{-1}+E\left({\phi }_x^2+{\phi }_{\text{z}}^2\right)\end{array}$ ， (19)

$w=h_t+u{h}_x$ ， (20)

$\left(T_z-h_x{T}_x\right){\left(1+h_x^2\right)}^{-1/2}+BiT=0$ ， (21)

在$z=d$ 处平面电势$\phi =0$ 。其中， (22)

$Ma=\frac{\text{Δ}\sigma }{U_0{\mu }^e}$ ，$We=\frac{{\sigma }_{\infty }}{U_0{\mu }^e}=\frac{K}{{\text{Re}}^{2/3}{\text{sin}}^{1/3}\beta }$ ，$Bi=\frac{{\kappa }_g}{{\kappa }_{th}}H$ ，$\text{Re}=\frac{H{U}_0}{\nu }$ ，

$\text{Pr}=\frac{\nu }{\kappa }$ ，${\mu }_1=\frac{{\mu }^o}{{\mu }^e}$ ，$E=\frac{{\phi }_c^2{\epsilon }_0}{2{\mu }^e{U}_{0}H}$ ，$\tau =\frac{{\tau }_{s}H}{{\mu }^e{U}_0}$ ， (23)

即马兰戈尼数Ma (用于表征表面张力梯度引起的热毛细应力与重力诱导粘性应力的比值)、韦伯数We (用于表征自由表面的Blob压力与重力诱导压力的比值) (Kapitza数K在20℃时约为3000 [9])、自由表面Biot数Bi (用于表征表面对流换热与内部导热的相对强度)、雷诺数Re (用于表征惯性力与粘性力的比值)、普朗特数Pr (用于表征动量扩散与热扩散的相对速率)、Péclet数Pe = RePr (用于表征对流换热与导热的相对强度)、无量纲奇粘度系数μ₁、电参数E、无量纲外加剪切应力τ。值得注意的是，电参数E和马兰戈尼数Ma均为0时，上述无量纲控制方程与边界条件退化为Desai等人[8]推导的结果。

3. 长波近似

基于润滑理论核心假设，结合长波流体流动特征，假设扰动波长远大于液膜特征厚度，引入小参数ε (表征长波近似量级，$\epsilon \ll 1$ )，对空间与时间偏导数进行量级缩放：

$\left({\partial }_t,{\partial }_x\right)\to \epsilon \left({\partial }_{\tau },{\partial }_X\right)$ ，${\partial }_{xx}\to {\epsilon }^2{\partial }_{XX}$ ，$w\to \epsilon v$ 。 (24)

将上述量级缩放关系代入无量纲控制方程(12)~(16)，得到长波近似下的简化控制方程组：

$u_X+v_z=0$ ， (25)

$\epsilon \text{Re}\left(u_{\tau }+u{u}_X+v{u}_z\right)=-\epsilon p_X+1+{\epsilon }^2{u}_{XX}+u_{zz}-{\mu }_1\epsilon \left({\epsilon }^2{v}_{XX}+v_{zz}\right)$ ， (26)

${\epsilon }^2\text{Re}\left(v_{\tau }+u{v}_X+v{v}_z\right)=-p_z-\text{cos}\beta +{\epsilon }^3{v}_{XX}+\epsilon v_{zz}+{\mu }_1\left({\epsilon }^2{u}_{XX}+u_{zz}\right)$ ， (27)

$\epsilon Pe\left(T_{\tau }+u{T}_X+v{T}_z\right)={\epsilon }^2{T}_{XX}+T_{zz}$ 。 (28)

区域二满足

${\epsilon }^2{\phi }_{XX}+{\phi }_{zz}=0$ 。 (29)

相应地，无量纲边界条件(17)~(22)经长波近似后转化为如下形式：

在$z=0$ 平面，有$u=0$ ，$v=0$ ，$T=1$ 。 (30)

在$z=h\left(X,\tau \right)$ 处，有电势$\phi =1$ ，

$\begin{array}{l}\left[\left(u_z+{\epsilon }^2{v}_X\right)\left(1-{\epsilon }^2{h}_X^2\right)-2{\epsilon }^2\left(u_X-v_z\right)h_X\right]+{\mu }_1\left[\epsilon \left(u_X-v_z\right)\left(1-{\epsilon }^2{h}_X^2\right)+2\epsilon \left(u_z+{\epsilon }^2{v}_X\right)h_X\right]\\ =-Ma\epsilon \left(T_X+h_X{T}_z\right){\left(1+{\epsilon }^2{h}_X^2\right)}^{1/2}+\tau {\left(1+{\epsilon }^2{h}_X^2\right)}^{1/2}\end{array}$ ， (31)

$\begin{array}{l}p=\left\{2\epsilon \left[{\epsilon }^2{u}_X{h}_X^2-\left(u_z+{\epsilon }^2{v}_X\right)h_X+v_z\right]+{\mu }_1\left[\left(u_z+{\epsilon }^2{v}_X\right)\left(1-{\epsilon }^2{h}_X^2\right)-2{\epsilon }^2\left(u_X-v_z\right)h_X\right]\right\}{\left(1+{\epsilon }^2{h}_X^2\right)}^{-1}\\ -\left(We-MaT\right){\epsilon }^2{h}_{XX}{\left(1+{\epsilon }^2{h}_X^2\right)}^{-3/2}-2E{\left({\epsilon }^2{h}_X{\phi }_X-{\phi }_y\right)}^2{\left(1+{\epsilon }^2{h}_X^2\right)}^{-1}+E\left({\epsilon }^2{\phi }_X^2+{\phi }_y^2\right)\end{array}$ ， (32)

$v=h_{\tau }+u{h}_X$ ， (33)

$\left(T_z-{\epsilon }^2{h}_X{T}_X\right){\left(1+{\epsilon }^2{h}_X^2\right)}^{-1/2}+BiT=0$ 。 (34)

在$z=d$ 处平面电势$\phi =0$ 。 (35)

在式(24)~(34)中，忽略高阶小项$O\left({\epsilon }^2\right)$ 。为保证方程推导的一致性与精度，本研究对关键参数的量级做出如下假设：雷诺数Re、马兰戈尼数Ma、无量纲奇粘度系数${\mu }_1$ 、无量纲剪切应力τ、电参数E及Biot数Bi均为$O\left(1\right)$ 量级；普朗特数Pr为$O\left(\epsilon \right)$ 量级(遵循Desai等人[7]的研究范式)；韦伯数We为$O\left({\epsilon }^{-}{}^2\right)$ 量级，定义无量纲表面张力参数$W={\epsilon }^{2}We$ ，$W\to O\left(1\right)$ ，以适配长波近似下的量级平衡。

通过联立式(25)与(31)，忽略O(ε²)，求解得到长波近似下的压力场分布

$p=-\text{cot}\beta \left(z-h\right)+{\mu }_1{u}_z-W{h}_{XX}-\frac{E}{{\left(h-d\right)}^2}+O\left(\epsilon \right)$ ， (36)

将压力场表达式(35)代入式(25)，结合连续性方程(24)，进一步化简得到长波近似下的流向动量输运方程：

$\epsilon Re\left(u_{\tau }+u{u}_X+v{u}_z\right)=1+u_{zz}+\epsilon W{h}_{XXX}-\epsilon \text{cot}\beta h_X-\frac{\epsilon E{h}_X}{{\left(h-d\right)}^3}$ 。 (37)

4. 构建自由表面方程

本尼方程模型是长波近似下处理重力作用下不可压缩流体薄膜流动的经典方法，适用于长波主导的流动稳定性分析。本研究采用该模型，假设各物理量沿小参数$\epsilon$ 展开为一阶渐近形式$\left(u,v,p,T,\phi \right)=\left(u_0,v_0,p_0,T_0,{\phi }_0\right)+\epsilon \left(u_1,v_1,p_1,T_1,{\phi }_1\right)$ ，其中$\left(u_0,v_0,p_0,T_0,{\phi }_0\right)$ 为前导阶解，$\left(u_1,v_1,p_1,T_1,{\phi }_1\right)$ 为一阶修正解。

将上述渐近展开式代入长波近似下的控制方程组(24)~(28)及边界条件(29)~(34)，通过分离不同量级项，建立前导阶方程组与对应边界条件，求解得到前导阶解：

$T_0=1-\frac{Biz}{1+Bih}$ ，$u_0=\left(h+\tau \right)z-\frac{1}{2}{z}^2$ ，$v_0=-\frac{1}{2}{z}^2{h}_X$ ，$p_0=\tau {\mu }_1-\left(\text{cot}\beta +{\mu }_1\right)\left(z-h\right)-W{h}_{XX}-\frac{E}{{\left(h-d\right)}^2}$ 。(38)

将前导阶解代入一阶修正方程组，结合边界条件约束求解，得到一阶修正解中的流向速度分量：

$\begin{array}{c}{u}_1=Ma{h}_X\frac{Bi}{{\left(1+Bih\right)}^2}z-2{\mu }_1\left(h+\tau \right)h_{X}z+\left(\text{cot}\beta h_X-W{h}_{XXX}+\frac{E}{{\left(h-d\right)}^3}{h}_X\right)\left(\frac{1}{2}{z}^2-hz\right)\\ +\text{Re}\left(h+\tau \right)\left(\frac{z^4}{24}+\frac{z{h}^3}{3}-\frac{z^{3}h}{6}\right)h_X\end{array}$ ， (39)

根据流体连续性原理，液膜厚度演化满足$h_{\tau }+q_x=0$ ，其中$q={\displaystyle {\int }_0^{h}u\text{d}z}$ 为截面体积流量。将u = u₀ + u₁代入流量表达式，结合厚度演化方程，经量级整理$\left({\partial }_t,{\partial }_x\right)\to \epsilon \left({\partial }_{\tau },{\partial }_X\right)$ ，最终得到液膜自由表面的演化方程：

$\begin{array}{l}{h}_t+\left[\frac{1}{3}{h}^3+\frac{1}{2}{h}^2\tau +Ma{h}_x\frac{Bi}{2{\left(1+Bih\right)}^2}{h}^2\right[\\ {-{\mu }_1\left(h+\tau \right)h^2{h}_x-\frac{1}{3}{h}^3\left(\text{cot}\beta h_x-We{h}_{xxx}+\frac{E}{{\left(h-d\right)}^3}{h}_x\right)+\frac{2}{15}Re\left(h+\tau \right)h^5{h}_x]}_x=0,\end{array}$ (40)

其中各项分别代表了平均流量引起的对流项、外加剪切应力的驱动项、表面温度梯度引起的热毛细应力、奇粘性的影响、重力项、表面张力效应、外加电场的调控作用与流体惯性效应。

5. 线性稳定性分析

采用线性小扰动方法对自由表面演化方程开展稳定性分析。在之前的工作中，我们得到了自由表面h的演化方程(39)，整理为如下形式

$h_t+{\left(c_1+c_2{h}_x+c_3{h}_{xxx}\right)}_x=0$ ， (41)

其中

$c_1=\frac{1}{3}{h}^3+\frac{1}{2}{h}^2\tau$ ，$c_2=\frac{MaBi{h}^2}{2{\left(1+Bih\right)}^2}-{\mu }_1\left(h+\tau \right)h^2-\frac{1}{3}\text{cot}\beta h^3-\frac{1}{3}\frac{E}{{\left(h-d\right)}^3}{h}^3+\frac{2}{15}Re\left(h+\tau \right)h^5$ ，$c_3=\frac{1}{3}{h}^{3}We$ 。(42)

假设液膜厚度的基本解为恒定值1，在基本解基础上施加小振幅扰动$h\left(x,t\right)=1+\widehat{h}\left(x,t\right)+c.c.$ ，将液膜厚度扰动$\widehat{h}$ 带入方程(40)，并对方程做泰勒展开，忽略四阶项$O\left({\widehat{h}}^4\right)$ 得到

$\begin{array}{l}{\widehat{h}}_t+{c^{\prime }}_1{\widehat{h}}_x+c_2{\widehat{h}}_{xx}+c_3{\widehat{h}}_{xxxx}+{c^{″}}_1\widehat{h}{\widehat{h}}_x+{c^{\prime }}_2\widehat{h}{\widehat{h}}_{xx}+{c^{\prime }}_2{\widehat{h}}_x^2+{c^{\prime }}_3\widehat{h}{\widehat{h}}_{xxxx}+{c^{\prime }}_3{\widehat{h}}_x{\widehat{h}}_{xxx}\\ +\frac{{c^{‴}}_1}{2}{\widehat{h}}^2{\widehat{h}}_x+\frac{{c^{″}}_2}{2}{\widehat{h}}^2{\widehat{h}}_{xx}+{c^{″}}_2\widehat{h}{\widehat{h}}_x^2+\frac{{c^{″}}_3}{2}{\widehat{h}}^2{\widehat{h}}_{xxxx}+{c^{″}}_3\widehat{h}{\widehat{h}}_x{\widehat{h}}_{xxx}+O\left({\widehat{h}}^4\right)=0\end{array}$ 。 (43)

我们假设$\widehat{h}\left(x,t\right)=\widehat{h}\text{exp}\left(ikx-i\omega t\right)+c.c.$其中$\widehat{h}$ 为扰动振幅(满足$\left|\widehat{h}\right|\ll h_0$ )，k为无量纲波数，忽略扰动项的高阶小量($O\left({\widehat{h}}^2\right)$ 及以上)，对方程进行线性化展开，整理得到色散关系：

$D\left(\omega ,k\right)=\left(1+\tau \right)k+\frac{BiMai{k}^2}{2{\left(1+Bi\right)}^2}-\frac{1}{3}\left(\text{cot}\beta +W{h}_{XXX}{k}^2+\frac{E}{{\left(1-d\right)}^3}\right)i{k}^2+\left(\frac{2}{15}Re-{\mu }_1\right)\left(1+\tau \right)i{k}^2$ 。 (44)

由式(44)，扰动增长率ω_i的符号可以判定流动的稳定性，ω_i = 0时流动临界稳定，得到临界雷诺数Re_c

${\text{Re}}_c=\frac{5}{2\left(1+\tau \right)}\left(3{\mu }_1\left(1+\tau \right)+\text{cot}\beta +k^{2}W\text{e}-Ma\frac{3Bi}{2{\left(1+Bi\right)}^2}+\frac{E}{{\left(1-d\right)}^3}\right)$ 。 (45)

6. 弱非线性稳定性分析

为精准捕捉液膜流动失稳后的非线性演化特征，采用多尺度方法对自由表面演化方程开展弱非线性稳定性分析。将液膜厚度扰动$\widehat{h}$ 沿小参数$\phi$ (表征扰动幅值量级，$\phi \ll 1$ )展开为幂级数形式：

$\widehat{h}\left(x,x_1,\cdots ;t,t_1,t_2\cdots \right)=\phi {\widehat{h}}_1+{\phi }^2{\widehat{h}}_2+{\phi }^3{\widehat{h}}_3+\cdots$ ， (46)

其中t和x是快速变化的尺度，x₁，t₁和t₂是缓慢变化的尺度，各尺度间相互独立。对偏导数进行如下展开：

${\partial }_t={\partial }_t+\phi {\partial }_{t_1}+{\phi }^2{\partial }_{t_2}+\cdots$ ，${\partial }_x={\partial }_x+\phi {\partial }_{x_1}+\cdots$ 。 (47)

将厚度扰动展开式(45)与偏导数展开式(46)代入线性化后的自由表面演化方程(42)，按小参数$\phi$ 的不同量级分离方程，得到各阶扰动分量的控制方程及对应的边界条件。首先考虑一阶项$O\left(\phi \right)$ ，整理得到${\widehat{h}}_1$ 的控制方程$D_0^{11}{\widehat{h}}_1=0$ ($D_0^{11}$ 为一阶偏导算子)，假设其解为单色行波形式${\widehat{h}}_1=\widehat{h}\text{exp}\left[i\left(kx-{\omega }_{r}t\right)\right]+c.c.$，其中，k为波数，${\omega }_r$ 为线性扰动角频率，c. c.表示复共轭项。

进一步考虑二阶项$O\left({\phi }^2\right)$ ，忽略高阶小项后得到${\widehat{h}}_2$ 的非齐次控制方程，假设其解为二次谐波形式：

${\widehat{h}}_2=F_2{\widehat{h}}^2\text{exp}\left[2i\left(kx-{\omega }_{r}t\right)\right]+c.c.$， (48)

将该解代入二阶控制方程，结合一阶解的边界条件约束，求解得到系数F₂的表达式为

$F_2=\frac{2k\left({c^{\prime }}_2-{c^{\prime }}_3{k}^2\right)-i{c^{\prime }}_1{k}^2}{4k\left(4c_3{k}^2-c_2\right)}$ ， (49)

式中，$c_1$ 、$c_2$ 、$c_3$ 由式(42)给出，${c^{\prime }}_1$ 、${c^{\prime }}_2$ 、${c^{\prime }}_3$ 分别为其对液膜厚度h的一阶导数。

针对$O\left({\phi }^3\right)$ ，整理得到${\widehat{h}}_3$ 的控制方程$D_0^{11}{\widehat{h}}_3+D_1^{11}{\widehat{h}}_2+D_2^{11}{\widehat{h}}_1=-{\text{Ν}}_3^1$ ($D_1^{11}$ 、$D_2^{11}$ 为高阶偏导算子，${\text{N}}_3^1$ 为三阶非线性项)。根据三阶方程的可解性条件，得到扰动振幅的标准Ginzburg-Landau方程：

${\widehat{h}}_{t_2}+iG{\widehat{h}}_{x_1}+J_1{\widehat{h}}_{x_1{x}_1}-{\phi }^{-2}{\omega }_i\widehat{h}+\left(J_2+i{J}_4\right){\left|\widehat{h}\right|}^2\widehat{h}=0$ ， (50)

其中，G为色散修正系数，$J_1=c_2-6c_3{k}^2$ 为扩散系数，J₄为相位修正系数，J₂为朗道系数(其表达式为$J_2=i{c}_1{F}_{2}k+1/2i{c^{″}}_{1}k-{c^{\prime }}_2{F}_2{k}^2-1/2{c^{″}}_2{k}^2+7{c^{\prime }}_3{F}_2{k}^4+1/2{c^{″}}_3{k}^4$ )，${\omega }_i$ 为线性扰动增长率。

忽略方程中的扩散效应与相位修正项，假设扰动振幅解的形式为$\widehat{h}=a{\text{e}}^{-ib\left(t_2\right)t_2}$ (a为振幅幅值，$b\left(t_2\right)$ 为相位因子)，代入Ginzburg-Landau方程后化简得到振幅演化方程(J_2,r为朗道系数J₂的实部)：

$\frac{\partial a}{\partial t_2}=\left({\phi }^{-2}{\omega }_i-J_{2,r}{a}^2\right)a$ 。 (51)

当扰动达到平衡状态时，$\partial a/\partial t_2=0$ ，解得反映非线性失稳的强度的平衡振幅$\zeta a={\left({\omega }_i/J_{2,r}\right)}^{1/2}$ 。

7. 数值研究

结合流体连续性原理与自由表面演化方程，建立液膜厚度h与截面体积流量q的耦合控制方程组：

$h_t+q_x=0$ ， (52)

$\begin{array}{c}q=\frac{1}{3}{h}^3+\frac{1}{2}{h}^2\tau +Ma{h}_x\frac{Bi}{2{\left(1+Bih\right)}^2}{h}^2-{\mu }_1\left(h+\tau \right)h^2{h}_x\\ -\frac{1}{3}{h}^3\left(\text{cot}\beta h_x-We{h}_{xxx}+\frac{E}{{\left(h-d\right)}^3}{h}_x\right)+\frac{2}{15}Re\left(h+\tau \right)h^5{h}_x.\end{array}$ (53)

结合倾斜加热平板液膜流动的物理特性，设定如下边界约束条件，以排除边界扰动对流动稳定性的干扰：在出口边界处$\left(x=L\right)$ ：液膜厚度$h\left(\underset{\_}{L},t\right)=0.1$ ，自由表面斜率满足$h_x\left(L,t\right)=0$ ，保证出口处液膜形态稳定；入口边界处$\left(x=0\right)$ ：自由表面斜率$h_x\left(0,t\right)=0$ ，三阶导数$h_{xxx}\left(0,t\right)=0$ ，避免入口处附加扰动诱发液膜失稳。

为模拟小扰动作用下的液膜流动演化，设定初始液膜厚度表达式为：

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

其中，无量纲基态厚度为1，扰动振幅为0.05，满足小扰动假设($\left|\widehat{h}\right|\ll 1$ )，k为无量纲扰动波数。

构建均匀空间网格$x_i=\left(i-1\right)\Delta x,i=1,2,\cdots ,N$ ，其中$\Delta x=L/N$ 为网格步长，N为网格节点数。为保证高波数扰动的捕捉精度，选取网格节点数$N\ge 50$ 。为满足流动过程中的质量守恒条件，将截面体积流量q在网格半节点$i=\pm 1/2$ 处离散，液膜厚度的时间导数采用中心差分格式：

$h_{it}=\frac{q_{i+\frac{1}{2}}-q_{i-\frac{1}{2}}}{\Delta x}$ 。 (55)

对内部网格节点($3\le i\le N-2$ )，采用中心差分格式离散各阶导数项，得到半节点流量的表达式：

$\begin{array}{c}{q}_{i+\frac{1}{2}}=\frac{1}{6}\left(h_i^3+h_{i+1}^3\right)+\frac{\tau }{4}\left(h_i^2+h_{i+1}^2\right)+\frac{MaBi\left(h_i^2+h_{i+1}^2\right)}{{\left(2+Bi\left(h_i+h_{i+1}\right)\right)}^2}\frac{h_{i+1}-h_i}{\Delta x}-\frac{{\mu }_1}{2}\left(h_i^3+h_{i+1}^3\right)\frac{h_{i+1}-h_i}{\Delta x}-\frac{{\mu }_1}{2}\tau \left(h_i^2+h_{i+1}^2\right)\frac{h_{i+1}-h_i}{\Delta x}\\ -\frac{1}{6}\left(h_i^3+h_{i+1}^3\right)\left(\cot \beta \frac{h_{i+1}-h_i}{\Delta x}-We\frac{h_{i+2}-3h_{i+1}+3h_i-h_{i-1}}{\Delta x^3}\right)-\frac{2}{3}\frac{E\left(h_i^3+h_{i+1}^3\right)}{{\left(h_i+h_{i+1}-2d\right)}^3}\frac{h_{i+1}-h_i}{\Delta x}\\ +\frac{1}{15}\mathrm{Re}\left(h_i^6+h_{i+1}^6\right)\frac{h_{i+1}-h_i}{\Delta x}+\frac{\tau }{15}\mathrm{Re}\left(h_i^5+h_{i+1}^5\right)\frac{h_{i+1}-h_i}{\Delta x}\end{array}$ (56)

$\begin{array}{c}{q}_{i-\frac{1}{2}}=\frac{1}{6}\left(h_{i-1}^3+h_i^3\right)+\frac{\tau }{4}\left(h_{i-1}^2+h_i^2\right)+\frac{MaBi\left(h_{i-1}^2+h_i^2\right)}{{\left(2+Bi\left(h_i+h_{i-1}\right)\right)}^2}\frac{h_i-h_{i-1}}{\Delta x}-\frac{{\mu }_1}{2}\left(h_{i-1}^3+h_i^3\right)\frac{h_i-h_{i-1}}{\Delta x}\\ -\frac{{\mu }_1}{2}\tau \left(h_{i-1}^2+h_i^2\right)\frac{h_i-h_{i-1}}{\Delta x}-\frac{1}{6}\left(h_{i-1}^3+h_i^3\right)\left(\cot \beta \frac{h_i-h_{i-1}}{\Delta x}-We\frac{h_{i+1}-3h_i+3h_{i-1}-h_{i-2}}{\Delta x^3}\right)-\frac{2}{3}\frac{E\left(h_{i-1}^3+h_i^3\right)}{{\left(h_i+h_{i+1}-2d\right)}^3}\frac{h_i-h_{i-1}}{\Delta x}\\ +\frac{1}{15}\mathrm{Re}\left(h_{i-1}^6+h_i^6\right)\frac{h_i-h_{i-1}}{\Delta x}+\frac{\tau }{15}\mathrm{Re}\left(h_{i-1}^5+h_i^5\right)\frac{h_i-h_{i-1}}{\Delta x}.\end{array}$ (57)

为满足边界条件的约束，引入域外鬼点补充边界附近的导数计算：出口边界($x=L$ )：根据给定厚度条件，取$h_{N+1}=0.1$ ；基于平坦性假设，取$h_N{}_{+2}=h_N$ ；入口边界($x=0$ )：基于流动对称性条件，取$h_0=h_2$ ，$h_{-1}=h_3$ 。

方程组(54)~(56)可以用标准微分方程求解器求解。

8. 讨论

本研究已推导得到本模型下流体薄膜的线性稳定性弱非线性稳定性解析解。后续分析在参数β = π/2、d = 2、Bi = 1的前提下开展。首先，忽略奇粘度作用(μ₁ = 0)，探究剪切应力τ、马兰戈尼数Ma、电参数E对流动线性稳定性的调控作用，相关结果由图2呈现；其次，固定Ma = 1且忽略外加电场作用(E = 0)，分析奇粘度μ₁及其与剪切应力τ的耦合效应对线性稳定性的影响机制，图3进行了详细展示；接着，基于上述固定参数及多物理场耦合条件，开展弱非线性稳定性分析，结合失稳区域的划分规律探讨非线性作用特征，图4、图5分别呈现了相关分析结果，相关参数设置与图2、图3保持一致；最后，采用数值求解方法，绘制液膜自由表面的时空演化图与时间演化曲线，图6~图8直观呈现了数值模拟结果。

8.1. 线性稳定性分析

图2直观呈现了忽略奇粘度作用(μ₁ = 0)时，流动线性稳定性受参数调控效果，参数设置与前文设定一致。图2(a)~(c)为Re-k平面内流动中性稳定性曲线(ω_i = 0)的参数依赖特性，曲线之上区域对应ω_i < 0，扰动能量经粘性耗散或外力约束发生衰减，流动维持稳定状态；曲线之下区域对应ω_i > 0，扰动能量不断累积叠加，流动逐步趋向失稳。图2(d)~(f)反映了扰动增长率ω_i随波数k的变化规律。由图2(a)与图2(d)可观察到，外加剪切应力τ的增大破坏了流体的粘性平衡状态，显著加剧流动不稳定性，使失稳进程加快，最大扰动增长率同步升高；与此同时，小波数区间(k < 0.015)内流动始终处于失稳状态，且不受τ变化的影响，区间内扰动源持续存在，扰动能量难以通过粘性耗散过程消除，最终形成全局失稳的特征。由图2(b)与图2(e)可知，电参数E增大时，流体内部会诱导产生附加应力场，使得不稳定区域范围扩大，最大扰动增长率提升，进而加速流动失稳；当E ≤ 0.5时，小雷诺数下流体惯性力可忽略不计，电场作用无法突破粘性耗散的约束，流动保持全局稳定。由图2(c)与图2(f)可得出，马兰戈尼数Ma的增大会诱导沿界面的温度梯度驱动流形成，促使不稳定区域扩展、最大扰动增长率提高，推动流动失稳；由式(45)，满足3MaBi(1 − d)³ = 2E(1 + Bi)²这一临界条件时，临界雷诺数Re_c中热毛细作用与外加电场效应相互抵消，流动稳定性达到动态平衡。

图3(a)与图3(b)呈现了Ma = 1且E = 0时，μ₁对线性稳定性的影响。随着μ₁的增大，其对扰动的“阻尼效应”不断增强，能够快速吸收扰动能量并将其转化为热能耗散，促进流动稳定，最大扰动增长率随之降低；小雷诺数下，粘性耗散主导流动，液膜保持全局稳定。图3(c)进一步呈现了μ₁ = 1时，外加剪切应力τ变化的影响，与图2(a)相比，奇粘度抑制剪切应力诱导的扰动增长，扩大了稳定的范围。且小雷诺数下，流动依旧保持全局稳定。图3(c)中，剪切应力对稳定性调控呈现“反转”，不同波数对应扰动尺度下，流场主导机制发生切换：波数k < P₁时为长波扰动，对流输运作用占据主导，奇粘度的粘性耗散高效吸收扰动能量，剪切应力的增大促进流体混合、削弱局部扰动集中，协同提升液膜稳定性；波数k > P₁时为短波扰动，流场界面效应凸显，表面张力与奇应力的耦合产生共振效应，剪切应力引发的扰动频率与流体本征频率匹配，叠加后的扰动能量突破奇粘度的抑制，剪切应力增大为短波扰动补充能量，最终加剧液膜失稳。

Figure 2. On (Re, k) plane, neutral stability curves for (a) different shear stresses τ when the electric parameter E = 0 and Marangoni number Ma = 1; (b) different electric parameters E when Ma = 1 and τ = 0; (c) different Marangoni numbers Ma when E = 0 and τ = 0. When Reynolds number Re = 30, disturbance growth rate ω_i vs wave number k for (d) different τ when E = 0 and Ma = 1; (e) different E when Ma = 1 and τ = 0; (f) different Ma when E = 0 and τ = 0

图2. (Re, k)平面上，(a) 电参数E = 0，马兰戈尼数Ma = 1时，不同剪切应力τ下；(b) Ma = 1时，τ = 0时不同E下；(c) E = 0，τ = 0时，不同Ma下的中性稳定性曲线。雷诺数Re = 30时，随波数k变化，(d) E = 0，Ma = 1时，不同τ下；(e) Ma = 1时，τ = 0时，不同E下；(f) E = 0，τ = 0时不同Ma下的扰动增长率ω_i曲线

Figure 3. When the applied shear stress τ = 0, on (Re, k) plane, (a) neutral stability curves with different μ₁. When Reynolds number Re = 30, with the variation of wave number k, (b) disturbance growth rate ω_i curves with different μ₁. When μ₁ = 1, on (Re, k) plane, (c) neutral stability curves with different τ

图3. 外加剪切应力τ = 0时，(Re, k)平面上，(a) 不同μ₁下的中性稳定性曲线。雷诺数Re = 30时，随波数k变化，(b) 不同μ₁的扰动增长率ω_i曲线。μ₁ = 1时，(Re, k)平面上，(c) 不同τ下的中性稳定性曲线

8.2. 弱非线性稳定性分析

图4呈现了外加剪切应力τ、外加电场电参数E和马兰戈尼数Ma对流动弱非线性稳定性的影响。图4(a)~(c)为Re-k平面上的弱非线性不稳定分析相图(ω_i = 0, J_2,r = 0)，划分为四个特征区域：超临界稳定区域I (ω_i > 0, J_2,r > 0)、亚临界不稳定区域II (ω_i < 0, J_2,r < 0)、无条件稳定区域III (ω_i < 0, J_2,r > 0)及爆炸状态区域IV (ω_i > 0, J_2,r < 0)。图4(d)~(f)反映了临界振幅ζa随波数k的变化。随着τ、E及Ma的增大，中性稳定性曲线与J_2,r曲线均向上偏移，导致Ⅱ、Ⅲ区面积缩小，Ⅰ、Ⅳ区面积扩大，流动趋于失稳，与线性稳定性分析相比更能长迅速，流动快速崩溃；Ⅱ、Ⅲ区为过渡区域，扰动增长速率随参数动态变化。非线性效应下参数的变化使流动更易趋近于“稳定极限”或“失稳极限”。图4(d)~(f)的临界振幅ζa曲线中，曲线左侧对应区域Ⅳ，上侧对应区域Ⅰ；随着τ、E及Ma的增大，区域Ⅳ的临界振幅ζa逐渐减小，区域Ⅰ的临界振幅ζa增大，流动的稳定区域扩展，更大的扰动振幅下仍能维持稳定；失稳区域的收缩则表明流动对扰动的耐受能力提升。

Figure 4. On (Re, k) plane, weakly nonlinear instability phase diagrams with (a) different τ when E = 0 and Ma = 1; (b) different E when Ma = 1 and τ = 0; (c) different Ma when E = 0 and τ = 0. When Re = 30, critical amplitude ζa vs k for (d) different τ when E = 0 and Ma = 1; (e) different E when Ma = 1 and τ = 0; (f) different Ma when E = 0 and τ = 0

图4. (Re, k)平面上，(a) E = 0, Ma = 1时，不同τ下；(b) Ma = 1, τ = 0时，不同E下；(c) E = 0, τ = 0时，不同Ma下的弱非线性不稳定分析相图。雷诺数Re = 30时，随波数k变化，(d) E = 0, Ma = 1时，不同τ下；(e) Ma = 1, τ = 0时，不同E下；(f) E = 0, τ = 0时，不同Ma下的临界振幅ζa曲线

图5(a)与图5(b)呈现了奇粘度μ₁对弱非线性不稳定性的调控。随着μ₁的增大，中性稳定性曲线与J_2,r曲线均向下偏移，导致Ⅱ、Ⅲ区面积扩大，Ⅰ、Ⅳ区面积缩小，临界振幅ζa同步降低，流动趋向稳定。从微观而言，奇粘度增大使横向约束力增强，可有效限制流体微团的侧向偏移与局部变形，减少表面扰动的幅值，临界振幅ζa的减小。图5(c)呈现了μ₁ = 1时，外加剪切应力τ变化对弱非线性不稳定性的影响，与图4(a)相比，相图曲线均右移，Ⅱ、Ⅲ区面积扩大，表明奇粘度弱化了剪切应力对非线性失稳的促进作用。

Figure 5. When τ = 0, on (Re, k) plane, (a) weakly nonlinear instability phase diagrams with different μ₁. When Re = 30, (b) critical amplitude ζa vs k for different μ₁. When μ₁ = 1, on (Re, k) plane, (c) weakly nonlinear instability phase diagrams with different τ

图5. 外加剪切应力τ = 0时，(Re, k)平面上，(a) 不同μ₁的弱非线性不稳定分析相图。雷诺数Re = 30时，随k变化，(b)不同μ₁的临界振幅ζa曲线。μ₁ = 1时，(Re, k)平面上，(c) 不同τ的弱非线性不稳定分析相图。

8.3. 自由表面时空演化分析

图6的演化曲线表明：随着雷诺数Re、马兰戈尼数Ma及电参数E的增大，液膜自由表面的波动幅度持续增强，波峰与波谷之间的差值显著扩大，波动频率同步升高，流动失稳程度不断加剧。Re增大使惯性力相对粘性力增强，削弱粘性力对扰动的抑制作用，Ma增大所诱导的热毛细驱动流以及E增大所引发的麦克斯韦附加应力均会增加扰动能量的输入，促使扰动持续放大并沿液膜流动方向传播。

Figure 6. When τ = 0 and μ₁ = 0, Free surface space-time evolution when (a) Re = 15, E = 0 and Ma = 0; (b) Re = 30, E = 2 and Ma = 0; (c) Re = 30, E = 0 and Ma = 10. Free surface temporal evolution when (d) Re = 15, E = 0 and Ma = 0; (e) Re = 30, E = 2 and Ma = 0; (f) Re = 30, E = 0 and Ma = 10

图6. τ = 0，μ₁ = 0时，(a) Re = 15，E = 0，Ma = 0时；(b) Re = 30，E = 2, Ma = 0时；(c) Re = 30，E = 0，Ma = 10时的自由表面时空演化。(d) Re = 15，E = 0，Ma = 0时；(e) Re = 30，E = 2，Ma = 0时；(f) Re = 30，E = 0，Ma = 10时的自由表面时间演化

图7呈现了剪切应力τ变化时，液膜自由表面的演化规律。当τ由0转为负值时，流体剪切方向与液膜流动方向相反，对扰动的“抑制性剪切效应”占据主导地位，波的频率与幅度减小，震荡强度显著弱化，液膜趋于稳定，下游区域出现的波峰翻折现象也随之消失；当τ由0逐渐增大时，上述现象呈现反向变化，剪切方向与液膜流动方向一致，扰动能量随剪切作用持续累积，波幅快速扩大、频率显著提升，初始时刻扰动向上游传播的趋势明显增强，下游区域的波峰翻折程度同步加剧，液膜流动失稳程度不断升高。

Figure 7. When Re = 30, E = 0, Ma = 0 and μ₁ = 0, Free surface space-time evolution when (a) τ = −0.5; (b)τ = 0; (c) τ = 1. Free surface temporal evolution when (d) τ = −0.5; (e) τ = 0; (f) τ = 1

图7. Re = 30，E = 0，Ma = 0，μ₁ = 0时，(a) τ = −0.5时；(b) τ = 0时；(c) τ = 1时的自由表面时空演化。(d) τ = −0.5时；(e) τ = 0时；(f) τ = 1时的自由表面时间演化

图8呈现了μ₁ = 1时，剪切应力τ对液膜自由表面演化的影响。奇粘度通过修正应力张量的非对称分量，进一步强化了剪切应力对扰动能量的“放大与抑制”效应，使τ对液膜自由表面演化的调控作用更加突出。具体而言，当τ减小时，奇粘度的“阻尼效应”与负向剪切的扰动抑制作用相互叠加，液膜波动的震荡幅度较图7进一步缩小，液膜流动稳定性得到显著提升；当τ增大时，尽管奇粘度可对扰动起到一定的抑制作用，但正向剪切所输入的扰动能量占据主导地位，液膜波动的震荡幅度扩大，且τ = 1时，初始时刻扰动向上游传播的趋势更为显著，下游区域的波峰翻折程度也同步加剧。这与Desai等人[7]提出的“奇粘度通过动量输运修正调控液膜失稳”的研究范式形成呼应，完善了奇粘度与剪切应力的耦合调控理论。

Figure 8. When Re = 30, E = 0, Ma = 0 and μ₁ = 2, free surface space-time evolution when (a) τ = −0.5; (b) τ = 0; (c) τ = 1. Free surface temporal evolution when (d) τ = −0.5; (e) τ = 0; (f) τ = 1

图8. Re = 30，E = 0，Ma = 0，μ₁ = 2时，(a) τ = −0.5时；(b) τ = 0时；(c) τ = 1时的自由表面时空演化。(d) τ = −0.5时；(e) τ = 0时；(f) τ = 1时的自由表面时间演化

9. 结论

本研究针对剪切–均匀热毛细–电场–奇粘度四场耦合下的倾斜平板液膜流动稳定性展开研究，通过理论建模、稳定性分析与数值模拟，发现满足$3MaBi{\left(1-d\right)}^3=2E{\left(1+Bi\right)}^2$ 时，热毛细诱导的界面张力梯度与电场麦克斯韦应力相互抵消。奇粘度μ₁通过修正应力张量的非对称分量增强扰动的能量耗散，拓宽稳定区间：小雷诺数下粘性耗散主导流动，液膜全局稳定。μ₁ = 1、Ma = 1且忽略电场时，剪切应力τ对稳定性的影响呈现波数依赖性：波数k小于交点P₁对应波数时，奇粘度耗散主导，τ增大促进流体混合、削弱扰动集中，液膜稳定性提升：k大于该波数时，剪切扰动频率与流体本征频率易共振，突破奇粘度抑制，τ增大加剧失稳。弱非线性稳定性分析中，Ma、E、τ增大时，Ⅱ、Ⅲ区面积扩大，Ⅰ、Ⅳ区面积扩大，失稳趋势增强；μ₁增大则推动区域演化反向变化，临界振幅ζa减小，流动趋向稳定。液膜自由表面演化显示，Ma、E、τ增大时，液膜波动幅度与频率提升，下游波峰翻折加剧：τ由0转为负值时，波频与幅度减小、震荡弱化、翻折消失，液膜稳定。μ₁ = 1显著放大τ的调控差异：τ减小时，阻尼与负向剪切抑制作用叠加，震荡幅度缩小：τ增大时，正向剪切能量输入占优，震荡幅度扩大，且τ = 1时上游扰动传播与下游翻折更显著。

基金项目

本研究由国家自然科学基金地区项目(12262026、12562027)、内蒙古自治区自然科学基金(2021MS01007、2025MS01004)以及内蒙古自治区高等学校创新团队项目(NMGIRT2323)资助。

NOTES

^*通讯作者。