Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                

A thermodynamically consistent phase-field lattice Boltzmann method for two-phase electrohydrodynamic flows

Fang Xiong Lei Wang leiwang@cug.edu.cn, wangleir1989@126.com Jiangxu Huang Kang Luo School of Mathematics and Physics, China University of Geosciences, Wuhan 430074, China School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China School of Energy Science and Engineering, Harbin Institute of Technology, Harbin 150001, China
Abstract

In this work, we aim to develop a phase-field based lattice Boltzmann (LB) method for simulating two-phase electrohydrodynamics (EHD) flows, which allows for different properties (densities, viscosities, conductivity and permittivity) of each phase while maintaining thermodynamic consistency. To this end, we first present a theoretical analysis on the two-phase EHD flows by using the Onsager’s variational principle, which is an extension of Rayleigh’s principle of least energy dissipation and, naturally, guarantees thermodynamic consistency. It shows that the governing equations of the model include the hydrodynamic equations, Cahn-Hilliard equation coupled with additional electrical effect, and the full Poisson-Nernst-Planck electrokinetic equations. After that, a coupled lattice Boltzmann (LB) scheme is constructed for simulating two-phase EHD flows. In particular, in order to handle two-phase EHD flows with a relatively larger electric permittivity ratio, we also introduce a delicately designed discrete forcing term into the LB equation for electrostatic field. Moreover, some numerical examples including two-phase EHD flows in planar layers and charge diffusion of a Gaussian bell are simulated with the developed LB method. It is shown that our numerical scheme shares a second-order convergence rate in space in predicting electric potential and charge density. Finally, we used the current model to simulate the deformation of a droplet under an electric field and the dynamics of droplet detachment in reversed electrowetting. Our numerical results align well with the theoretic solutions, and the available experimental/numerical data, demonstrating that the proposed method is feasible for simulating two-phase EHD flows.

keywords:
Electrohydrodynamics , Phase-field method , Onsager’s variational principle , Lattice Boltzmann method , Electrowetting
journal: Springer

1 Introduction

Electrohydrodynamics (EHD) is an interdisciplinary science that deals with the interaction of fluids with electric fields [1]. As one of the fundamental subfields of EHD, a deep understanding of the physics and dynamics of multiphase electrohydrodynamic flow is crucial because of its relevance in a broad range of microfluidics and industrial processes, such as ink-jet printing [2], electrospraying [3], manipulation of micro-drops by continuous electrowetting [4], oil extraction from oil-water emulsions [5], EHD pumps [6], to name a few. In this setting, in order to advance the design and development of the aforementioned technologies and explore new applications, it is crucial to gain a comprehensive understanding of the flow physics related to the multiphase EHD.

Generally, when an electric field is applied to a two-fluid system, a net force is generated at the fluid-fluid interface due to the imbalance of electrical properties between the fluids [7]. For perfectly dielectric or perfectly conducting fluids, the electric surface force solely acts in the direction normal to the fluid interface, which is further balanced by capillary traction due to surface tension [8, 9]. In such ideal situations, the fluids remain motionless at steady state, and the resultant phenomenon is usually called electrohydrostatics [9, 10, 11]. However, as there are no perfectly insulating or perfectly conducting fluids in the real world [12], researchers have long been interested in investigating fluids with finite electrical conductivity. In the pioneering work of Taylor [13], it was observed that the weakly conducting of the fluids allows electrical charges to accumulate at the droplet surface, resulting in a tangential electric stress to be generated, which in turn requires a tangential viscous stress to balance it, leading to an EHD flow formed in the system. The EHD theory proposed in Taylor’s classical study is known as the leaky dielectric model [1, 9], which is now widely employed by researchers in theoretical and numerical investigations of droplet deformation under a uniform electric field [14, 15, 16, 17, 18, 19]. It shows that Taylor’s leaky dielectric model is capable of predicting droplet deformations for small values of capillary number. However, there are some quantitative discrepancies between the leaky dielectric theory and experiment when the droplet undergoes large deformation [12, 20]. To unravel the cause of the discrepancies, in a noteworthy study, Feng et al. [12] investigated the EHD behavior of a drop at finite electric Reynolds number, which defined as the ratio of charge convection to Ohmic conduction. It is found that the nonlinear surface charge convection, which explicitly appears in the leaky-dielectric model but is neglected in the original study, may account for the discrepancies between Taylor’s theoretical prediction and experiment. On the basis of Feng et al.’s work [12], the influence of surface charge convection has garnered significant attention from some researchers in the EHD community [21, 22, 23, 24], and it is now well accepted that incorporating the mechanism of surface charge convection in the corresponding EHD model is more aligned with the actual physical circumstances.

Apart from the behavior of fluid interfaces under the action of electric fields, understanding droplet dynamics on a solid electrode is also of prime importance due to its wide applications in engineering [25, 26], especially in modern-day microfluidic systems [27, 28]. The first work to investigate the effect of an electric field on wetting behavior is attributed to Lippman [29], who found that the presence of a voltage tends to enhance the wettability of the substrate. This phenomenon is then referred to as electrowetting and has received much attention in the last few decades [30, 31, 32, 33, 34, 35, 36]. Roughly speaking, the existing methods for simulation of electrowetting can be classified into two major categories. The first category is the effective contact angle approach [31, 32, 33]. In this approach, the contribution of the electrical field on the contact line is established by the wetting boundary condition characterized by the Young-Lippmann equation. However, since this approach does not consider the governing equations for the electrostatic field, it cannot reproduce the distributions of the electric potential and the free charges in the system. To this end, some researchers have investigated electrowetting using a more realistic mathematical model in which the coupled hydrodynamics-electrostatics equations are included in the system [34, 35, 36]. Following this route, we have noticed that the two fluids used in this approach are typically assumed to be perfect dielectric fluids or follow the leaky dielectric model . Nevertheless, as mentioned in the preceding paragraph, these two models are established by using some assumptions, and thus, they are only applicable to corresponding specific situations. Of late, some researchers have adopted a more general model called the Poisson-Boltzmann equation to investigate electrowetting [37, 38, 39]. While the Poisson-Boltzmann equation being widely employed in studying ionic distribution in the electric double layer, one point that needs to be mentioned is that the Poisson-Boltzmann equation is derived from the Nernst–Planck equation under some assumptions [49]. If these assumptions are not satisfied, such as the cases of the surface being heterogeneously charged [40] or the overlapping electric double layer [41, 42], using the Poisson-Boltzmann equation may lead to incorrect results, and one has to adopt the Nernst–Planck equation in such cases. Many researchers have verified this statement in the single-phase flows [40, 41, 42]. Likewise, we believe that using the Nernst–Planck model to study the electrowetting multiphase problem should also be more physically plausible than the Poisson-Boltzmann equation.

The literature survey above indicates that understanding two-phase flows in an electric field requires full consideration of EHD. Unfortunately, since the evolution of the phase interface is inherently nonlinear and strongly coupled with the fluid flow and the electrostatic field, and simultaneously, various variables (such as velocity, density, and electrical properties) vary significantly or even sharply near the phase interface. These characteristics make investigating the interfacial phenomenon for two-phase EHD flow change challenging. Thus, it is desirable to develop numerical methods for such complex problems. Rooted in kinetic theory [43, 44], the lattice Boltzmann (LB) method has developed into a powerful and attractive method for simulating fluid flows and complex physical processes. In the past decades, the phase-field-based LB method [47, 48] and the pseudopotential LB method [45, 46], two popular methods for simulating multiphase flow in the LB community, have been widely adopted to simulate two-phase EHD flows. However, these existing LB methods are almost constructed based on the above-mentioned simplified models, such as the perfect dielectric model [50, 51], the leaky dielectric model [51, 17], and the Poisson-Boltzmann model [37, 39], which all share some assumptions to some extent, as mentioned above. Apart from the drawbacks mentioned above, it is also noted that the governing equations for two-phase EHD flow in most previous works are directly established by coupling the hydrodynamic and simplified Maxwell’s equations with the interface capturing model. Thus, these LB methods may not be fully consistent with the laws of thermodynamics, which plays a vital role in multiphase flow modeling.

Framed in this general background, the current work tends to propose a thermodynamically consistent phase-field LB method for two-phase EHD flows. To this end, this work first revisits the governing equations for two-phase EHD flows by using the Onsager’s variational principle [52, 53, 54]. It will be seen later that the thermodynamcially consistent physical model for two-phase EHD flows consist of the Navier-Stokes equations, the Cahn-Hilliard equation coupled with electrical effects, and the full Poisson-Nernst-Planck electrokinetic equations. After that, a general LB method is proposed for simulating two-phase EHD flows based on the resultant physical model. In particular, to simulate the two-phase EHD flows with relatively larger ratios of electric permittivity, which is usually entournted in the electrowetting phenomenon, an improved LB method for electric potential equation is developed. The remainder of the present paper is organized as follows. In Sect. 2, we will derive the governing equations for two-phase EHD flows with Onsager’s variational principle. In Sect. 3 and Sect. 4, the corresponding LB method and the boundary conditions are presented. In Sect. 5, we give the iterative algorithm for the current LB method. In Sect. 6, we conduct some numerical simulations to validate the developed LB method. In Sect. 7, the proposed LB method is used to simulate two practical problems. At last, a brief conclusion is drawn in Sect. 8.

2 The Mathematical model

In this section, we intend to develop a thermodynamically consistent phase-field model for two-phase EHD flows, in which the electrical and the Cahn-Hilliard equations are coupled to the hydrodynamic equations through extra stress that mimics surface tension and electrical force. In what follows, we first present the generalized governing equations for the two-phase EHD flows. Then, some undetermined terms that appeared in the macroscopic equations are established using Onsager’s variational principle to ensure thermodynamic consistency.

2.1 The generalized governing equations for the two-phase EHD flows

Assuming the two-fluids are immiscible and incompressible, and neglecting the influence of the gravitational force, the general formation of the mass and momentum conservation equations in the dimensional form can be expressed as [55]

𝐮=0,𝐮0\nabla\cdot\mathbf{u}=0,∇ ⋅ bold_u = 0 , (1)
ρ(𝐮t+𝐮𝐮)=𝚷+𝐅,𝜌𝐮𝑡𝐮𝐮𝚷𝐅\rho\left({\frac{{\partial{\bf{u}}}}{{\partial t}}+{\bf{u}}\cdot\nabla{\bf{u}}% }\right)=\nabla{\bf{\Pi}}+{\bf{F}},italic_ρ ( divide start_ARG ∂ bold_u end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ bold_u ) = ∇ bold_Π + bold_F , (2)

where t𝑡titalic_t, ρ𝜌\rhoitalic_ρ and 𝒖𝒖\bm{u}bold_italic_u are the time, mass density and velocity, respectively. 𝚷𝚷{\bf{\Pi}}bold_Π is a stress tensor that relates to the hydrodynamic pressure p𝑝pitalic_p and the viscous stress tensor μ[𝐮+(𝐮)T]𝜇delimited-[]𝐮superscript𝐮T\mu\left[{\nabla{\bf{u}}+{{\left({\nabla{\bf{u}}}\right)}^{\rm T}}}\right]italic_μ [ ∇ bold_u + ( ∇ bold_u ) start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT ], which can be expressed as

𝚷=p𝐈+12μ[𝐮+(𝐮)T],𝚷𝑝𝐈12𝜇delimited-[]𝐮superscript𝐮T{\bf{\Pi}}=-p{\bf{I}}+\frac{1}{2}\mu\left[{\nabla{\bf{u}}+{{\left({\nabla{\bf{% u}}}\right)}^{\rm T}}}\right],bold_Π = - italic_p bold_I + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_μ [ ∇ bold_u + ( ∇ bold_u ) start_POSTSUPERSCRIPT roman_T end_POSTSUPERSCRIPT ] , (3)

with μ𝜇\muitalic_μ being the fluid dynamics viscosity. The last term in Eq. (2), i.e., 𝐅𝐅\bf{F}bold_F, stands for the external force exerted on the fluids which contains the electrical and interfacial forces and its specific expression will be determined later.

Apart from the above hydrodynamic equations, considering the essential electrical laws in an EHD flow is another problem that must be specified. Since the current density due to the motion of charge carriers is too small, the influence of the magnetic field can be ignored [55]. Consequently, the electric field intensity 𝐄𝐄{\bf{E}}bold_E is irrotational, and Gauss’s law is restated as [1, 8]

(ε𝐄)=q,×𝐄=0,𝐄=φ,formulae-sequence𝜀𝐄𝑞formulae-sequence𝐄0𝐄𝜑\nabla\cdot\left({\varepsilon{\bf{E}}}\right)=q,\;\;\;\;\;\;\;\;\;\nabla\times% {\bf{E}}=0,\;\;\;\;\;\;\;\;\;{\bf{E}}=-\nabla\varphi,∇ ⋅ ( italic_ε bold_E ) = italic_q , ∇ × bold_E = 0 , bold_E = - ∇ italic_φ , (4)
qt+𝐮q+𝐉D=0,𝑞𝑡𝐮𝑞subscript𝐉𝐷0\frac{{\partial q}}{{\partial t}}+{\bf{u}}\cdot\nabla q+\nabla\cdot{{\bf{J}}_{% D}}=0,divide start_ARG ∂ italic_q end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ italic_q + ∇ ⋅ bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT = 0 , (5)

where ε𝜀\varepsilonitalic_ε, q𝑞qitalic_q and φ𝜑\varphiitalic_φ denote the permittivity, charge density and potential, respectively. Noting that the last equation in Eq. (5) is the so-called Nernst–Planck equation with 𝐉Dsubscript𝐉𝐷{{\bf{J}}_{D}}bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT being the current density due to conduction, which will be derived in what follows.

Since two-phase EHD flows are typical binary fluid flows, accurately capturing the interface behavior has a significant influence on numerical performance. The current work utilizes the phase-field method [47, 48], a popular diffuse interface approach for interface tracking, to model the interfacial dynamics. In this setting, the evolution of the order parameter ϕitalic-ϕ\phiitalic_ϕ (a parameter introduced to distinguish different phases) is described by the convective Cahn-Hilliard equation, which can be expressed as

ϕt+𝐮ϕ+𝐉ϕ=0,italic-ϕ𝑡𝐮italic-ϕsubscript𝐉italic-ϕ0\frac{{\partial\phi}}{{\partial t}}+\mathbf{u}\cdot\nabla\phi+\nabla\cdot{\rm{% }}{\mathbf{J}_{\phi}=0},divide start_ARG ∂ italic_ϕ end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ italic_ϕ + ∇ ⋅ bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = 0 , (6)

where 𝐉ϕsubscript𝐉italic-ϕ\mathbf{J}_{\phi}bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT is the undetermined phase flux density.

2.2 Onsager’s variational principle

As shown in Sec. 2.1, one can find that there exist some unknown terms (i.e., 𝐅𝐅\bf{F}bold_F, 𝐉Dsubscript𝐉𝐷{{\bf{J}}_{D}}bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT and 𝐉ϕsubscript𝐉italic-ϕ{{\bf{J}}_{\phi}}bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT) in the generalized governing equations, which must be established before developing numerical approach. To maintain thermodynamic consistency, we will adopt Onsager’s variational principle to determine the governing equations expressed in terms of the above unknowns. To present it more clearly, we first introduce the basic principle of Onsager’s variational principle in an isothermal system. We then use it to determine 𝐅𝐅\bf{F}bold_F, 𝐉Dsubscript𝐉𝐷{{\bf{J}}_{D}}bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT and 𝐉ϕsubscript𝐉italic-ϕ{{\bf{J}}_{\phi}}bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT appeared in the governing equations for two-phase EHD flows.

Onsager’s variational principle is a variational approach based on the reciprocal relations in linear irreversible thermodynamics [56]. For an open system, the principle states that the evolution equations can be obtained by maximizing the Onsager-Machlup action [52, 53, 54],

𝒪=S˙+S˙ΦS(α˙,α˙),𝒪˙𝑆superscript˙𝑆subscriptΦ𝑆˙𝛼˙𝛼\mathcal{O}=\dot{S}+{{\dot{S}}^{*}}-{\Phi_{S}}\left({\dot{\alpha},\dot{\alpha}% }\right),caligraphic_O = over˙ start_ARG italic_S end_ARG + over˙ start_ARG italic_S end_ARG start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT - roman_Φ start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT ( over˙ start_ARG italic_α end_ARG , over˙ start_ARG italic_α end_ARG ) , (7)

where the simple ” \cdot” represents the time derivatives, S˙˙𝑆\dot{S}over˙ start_ARG italic_S end_ARG is the rate of change of the entropy, S˙superscript˙𝑆{{\dot{S}}^{*}}over˙ start_ARG italic_S end_ARG start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT is the rate of entropy given by the system to the environment, ΦS(α˙,α˙)subscriptΦ𝑆˙𝛼˙𝛼{\Phi_{S}}\left({\dot{\alpha},\dot{\alpha}}\right)roman_Φ start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT ( over˙ start_ARG italic_α end_ARG , over˙ start_ARG italic_α end_ARG ) is the dissipation function defined as the half rate of the entropy production with α=(α1,α2,)𝛼subscript𝛼1subscript𝛼2\alpha=\left({{\alpha_{1}},{\alpha_{2}},\ldots}\right)italic_α = ( italic_α start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_α start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … ) being the set of variables that characterizes the non-equilibrium state of a system. In an isothermal system, it is known that the rate of entropy can be expressed as [57]

S˙=Q˙T=U˙T,superscript˙𝑆˙𝑄𝑇˙𝑈𝑇{{\dot{S}}^{*}}=-\frac{{\dot{Q}}}{T}=-\frac{{\dot{U}}}{T},over˙ start_ARG italic_S end_ARG start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT = - divide start_ARG over˙ start_ARG italic_Q end_ARG end_ARG start_ARG italic_T end_ARG = - divide start_ARG over˙ start_ARG italic_U end_ARG end_ARG start_ARG italic_T end_ARG , (8)

where T𝑇Titalic_T is the temperature, Q˙˙𝑄{\dot{Q}}over˙ start_ARG italic_Q end_ARG and U˙˙𝑈{\dot{U}}over˙ start_ARG italic_U end_ARG are the rate of heat transfer from the environment to the system, and the rate of change of the system energy, and they are equal according to the first thermodynamic law. Note that the rate of change in Helmholtz free-energy and free-energy dissipation function (defined as half the rate of free energy dissipation) at a constant temperature are respectively written as [52, 53, 54],

˙=T(S+S˙),ΦF(α˙,α˙)=TΦS(α˙,α˙).formulae-sequence˙𝑇𝑆superscript˙𝑆subscriptΦ𝐹˙𝛼˙𝛼𝑇subscriptΦ𝑆˙𝛼˙𝛼\mathcal{\dot{F}}=-T\left({S+{{\dot{S}}^{*}}}\right),\;\;\;\;\;\;\;\;\;{\Phi_{% F}}\left({\dot{\alpha},\dot{\alpha}}\right)=T{\Phi_{S}}\left({\dot{\alpha},% \dot{\alpha}}\right).over˙ start_ARG caligraphic_F end_ARG = - italic_T ( italic_S + over˙ start_ARG italic_S end_ARG start_POSTSUPERSCRIPT ∗ end_POSTSUPERSCRIPT ) , roman_Φ start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT ( over˙ start_ARG italic_α end_ARG , over˙ start_ARG italic_α end_ARG ) = italic_T roman_Φ start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT ( over˙ start_ARG italic_α end_ARG , over˙ start_ARG italic_α end_ARG ) . (9)

Submitting Eq. (9) into Eq. (7), one can observe that the maximization of the Onsager-Machlup action is equivalent to the minimization of the Rayleighian [52, 53, 54], i.e.,

=˙(α˙,α˙)+ΦF(α˙,α˙).˙˙𝛼˙𝛼subscriptΦ𝐹˙𝛼˙𝛼\mathcal{R}=\mathcal{\dot{F}}\left({\dot{\alpha},\dot{\alpha}}\right)+{\Phi_{F% }}\left({\dot{\alpha},\dot{\alpha}}\right).caligraphic_R = over˙ start_ARG caligraphic_F end_ARG ( over˙ start_ARG italic_α end_ARG , over˙ start_ARG italic_α end_ARG ) + roman_Φ start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT ( over˙ start_ARG italic_α end_ARG , over˙ start_ARG italic_α end_ARG ) . (10)

Because ˙˙\mathcal{\dot{F}}over˙ start_ARG caligraphic_F end_ARG is linear and ΦFsubscriptΦ𝐹{\Phi_{F}}roman_Φ start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT is quadratic in the rates α˙˙𝛼\dot{\alpha}over˙ start_ARG italic_α end_ARG, the above Rayleighian can be restated as

=iαiα˙i+12i,jζijα˙iα˙j,subscript𝑖subscript𝛼𝑖subscript˙𝛼𝑖12subscript𝑖𝑗subscript𝜁𝑖𝑗subscript˙𝛼𝑖subscript˙𝛼𝑗\mathcal{R}=\sum\limits_{i}{\frac{{\partial\mathcal{F}}}{{\partial{\alpha_{i}}% }}{{\dot{\alpha}}_{i}}}+\frac{1}{2}\sum\limits_{i,j}{{\zeta_{ij}}}{{\dot{% \alpha}}_{i}}{{\dot{\alpha}}_{j}},caligraphic_R = ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT divide start_ARG ∂ caligraphic_F end_ARG start_ARG ∂ italic_α start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG over˙ start_ARG italic_α end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∑ start_POSTSUBSCRIPT italic_i , italic_j end_POSTSUBSCRIPT italic_ζ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT over˙ start_ARG italic_α end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT over˙ start_ARG italic_α end_ARG start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT , (11)

where ζijsubscript𝜁𝑖𝑗{{\zeta_{ij}}}italic_ζ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT is the friction coefficient and it satisfies the reciprocal relation ζij=ζjisubscript𝜁𝑖𝑗subscript𝜁𝑗𝑖{\zeta_{ij}}={\zeta_{ji}}italic_ζ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT = italic_ζ start_POSTSUBSCRIPT italic_j italic_i end_POSTSUBSCRIPT. In the end, the kinetic equation is obtained by minimizing \mathcal{R}caligraphic_R with regard to the rates,

jζijα˙j=αi.subscript𝑗subscript𝜁𝑖𝑗subscript˙𝛼𝑗subscript𝛼𝑖\sum\limits_{j}{{\zeta_{ij}}{{\dot{\alpha}}_{j}}}=-\frac{{\partial\mathcal{F}}% }{{\partial{\alpha_{i}}}}.∑ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT italic_ζ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT over˙ start_ARG italic_α end_ARG start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT = - divide start_ARG ∂ caligraphic_F end_ARG start_ARG ∂ italic_α start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG . (12)

Following the above variational approach, it turns out that if we compute the total Helmholtz free-energy \mathcal{F}caligraphic_F and the dissipation function ΦFsubscriptΦ𝐹{\Phi_{F}}roman_Φ start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT for a two-phase EHD flow, then application of Onsager’s variational principle could yield the system’s time evolution equations. For completeness, the wetting boundary condition that accounts for the contact angle between the phase interface and the solid substrate is also incorporated into the following analysis. In such a case, the Landau free-energy functional in the phase field modelling can be written as [58, 8]

ϕ=Ω(ψ0+κ2|ϕ|2)𝑑Ω+Ωψs𝑑A,subscriptitalic-ϕsubscriptΩsubscript𝜓0𝜅2superscriptitalic-ϕ2differential-dΩsubscriptΩsubscript𝜓𝑠differential-d𝐴{{\mathcal{F}}_{\phi}}=\int_{\Omega}{\left({{\psi_{0}}+\frac{\kappa}{2}{{\left% |{\nabla\phi}\right|}^{2}}}\right)}d\Omega+\int_{\partial\Omega}{{\psi_{s}}}dA,caligraphic_F start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + divide start_ARG italic_κ end_ARG start_ARG 2 end_ARG | ∇ italic_ϕ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_d italic_A , (13)

where ΩΩ\Omegaroman_Ω and ΩΩ{\partial\Omega}∂ roman_Ω are the material volume and its solid boundary, ψ0subscript𝜓0\psi_{0}italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is the bulk free-energy density which takes the following double-well form ψ0=βϕ2(ϕ1)2subscript𝜓0𝛽superscriptitalic-ϕ2superscriptitalic-ϕ12{\psi_{0}}=\beta{\phi^{2}}{\left({\phi-1}\right)^{2}}italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = italic_β italic_ϕ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ϕ - 1 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [59], and this form suggests that the order parameter ϕitalic-ϕ\phiitalic_ϕ equal 0 and 1 represent the light and heavy fluids, respectively. ψssubscript𝜓𝑠{{\psi_{s}}}italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT is the wall free-energy per unit area at the solid boundary, and 0.5κ|ϕ|20.5𝜅superscriptitalic-ϕ20.5\kappa{\left|{\nabla\phi}\right|^{2}}0.5 italic_κ | ∇ italic_ϕ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT is excess interfacial energy at the fluid-fluid interface. β𝛽\betaitalic_β and κ𝜅\kappaitalic_κ are two coefficients that related to the interface surface tension γ𝛾\gammaitalic_γ and the interface width W𝑊Witalic_W, and they satisfy β=12γ/W,κ=3γW/2formulae-sequence𝛽12𝛾/𝑊𝜅3𝛾𝑊/2\beta={{12\gamma}\mathord{\left/{\vphantom{{12\gamma}W}}\right.\kern-1.2pt}W},% \kappa={{3\gamma W}\mathord{\left/{\vphantom{{3\sigma W}2}}\right.\kern-1.2pt}2}italic_β = 12 italic_γ start_ID / end_ID italic_W , italic_κ = 3 italic_γ italic_W start_ID / end_ID 2 [59, 48]. In addition to phase-field free energy, the presence of the electric field also provides a contribution to the total free-energy, which can be expressed as [8]

𝐄=12Ω|𝐃|2ε(ϕ)𝑑Ω,subscript𝐄12subscriptΩsuperscript𝐃2𝜀italic-ϕdifferential-dΩ{{\mathcal{F}}_{\bf{E}}}=\frac{1}{2}\int_{\Omega}{\frac{{{{\left|{\bf{D}}% \right|}^{2}}}}{\varepsilon(\phi)}}d\Omega,caligraphic_F start_POSTSUBSCRIPT bold_E end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_D | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_ε ( italic_ϕ ) end_ARG italic_d roman_Ω , (14)

where 𝐃=ε(ϕ)𝐄𝐃𝜀italic-ϕ𝐄{\bf{D}}=\varepsilon(\phi){\bf{E}}bold_D = italic_ε ( italic_ϕ ) bold_E is the electric displacement field. Also, the kinetic energy induced by the fluid inertia is given by [55]

𝐮=12Ωρ|𝐮|2𝑑Ω.subscript𝐮12subscriptΩ𝜌superscript𝐮2differential-dΩ{{\mathcal{F}}_{\bf{u}}}=\frac{1}{2}\int_{\Omega}{\rho{{\left|{\bf{u}}\right|}% ^{2}}}d\Omega.caligraphic_F start_POSTSUBSCRIPT bold_u end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_ρ | bold_u | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d roman_Ω . (15)

Additionally, the charge diffusive motion adds to the total free energy a contribution, which can be given by [58]

Fq=λ2Ωq2𝑑Ω,subscript𝐹𝑞𝜆2subscriptΩsuperscript𝑞2differential-dΩ{F_{q}}=\frac{\lambda}{2}\int_{\Omega}{{q^{2}}}d\Omega,italic_F start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT = divide start_ARG italic_λ end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_q start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d roman_Ω , (16)

in which λ=α/σ𝜆𝛼/𝜎\lambda={\alpha\mathord{\left/{\vphantom{\alpha\sigma}}\right.\kern-1.2pt}\sigma}italic_λ = italic_α start_ID / end_ID italic_σ with α𝛼\alphaitalic_α and σ𝜎\sigmaitalic_σ being the charge diffusion coefficient and electric conductivity, respectively. Framed in the above energy analysis, we can obtain the total Helmholtz free-energy as

=ϕ+𝐄+𝐮+q=Ω(ψ0+κ2|ϕ|2)𝑑Ω+Ωψs𝑑A+12Ω|𝐃|2ε(ϕ)𝑑Ω+12Ωρ|𝐮|2𝑑Ω+λ2Ωq2𝑑Ω.subscriptitalic-ϕsubscript𝐄subscript𝐮subscript𝑞subscriptΩsubscript𝜓0𝜅2superscriptitalic-ϕ2differential-dΩsubscriptΩsubscript𝜓𝑠differential-d𝐴12subscriptΩsuperscript𝐃2𝜀italic-ϕdifferential-dΩ12subscriptΩ𝜌superscript𝐮2differential-dΩ𝜆2subscriptΩsuperscript𝑞2differential-dΩ\begin{split}{\mathcal{F}}&={{\mathcal{F}}_{\phi}}+{{\mathcal{F}}_{\bf{E}}}+{{% \mathcal{F}}_{\bf{u}}}+{{\mathcal{F}}_{q}}\\ &=\int_{\Omega}{\left({{\psi_{0}}+\frac{\kappa}{2}{{\left|{\nabla\phi}\right|}% ^{2}}}\right)}d\Omega+\int_{\partial\Omega}{{\psi_{s}}}dA+\frac{1}{2}\int_{% \Omega}{\frac{{{{\left|{\bf{D}}\right|}^{2}}}}{\varepsilon(\phi)}}d\Omega+% \frac{1}{2}\int_{\Omega}{\rho{{\left|{\bf{u}}\right|}^{2}}}d\Omega+\frac{% \lambda}{2}\int_{\Omega}{{q^{2}}}d\Omega.\end{split}start_ROW start_CELL caligraphic_F end_CELL start_CELL = caligraphic_F start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT + caligraphic_F start_POSTSUBSCRIPT bold_E end_POSTSUBSCRIPT + caligraphic_F start_POSTSUBSCRIPT bold_u end_POSTSUBSCRIPT + caligraphic_F start_POSTSUBSCRIPT italic_q end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + divide start_ARG italic_κ end_ARG start_ARG 2 end_ARG | ∇ italic_ϕ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_d italic_A + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_D | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_ε ( italic_ϕ ) end_ARG italic_d roman_Ω + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_ρ | bold_u | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d roman_Ω + divide start_ARG italic_λ end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_q start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d roman_Ω . end_CELL end_ROW (17)

Taking the variational operator to the total Helmholtz free-energy, we have

δ=Ω(κ2ϕ+ψ0ϕ)δϕ𝑑Ω+Ωκ(ϕ)δϕ𝑑Ω+Ωψsϕδϕ𝑑AΩ𝜀(ϕ)2ε2(ϕ)|𝐃|2δϕ𝑑Ω=Ω[κ2ϕ+ψ0ϕ𝜀(ϕ)2ε2(ϕ)|𝐃|2]δϕ𝑑Ω+Ω(κ𝐧𝐰ϕ+ψsϕ)δϕ𝑑A,𝛿subscriptΩ𝜅superscript2italic-ϕsubscript𝜓0italic-ϕ𝛿italic-ϕdifferential-dΩsubscriptΩ𝜅italic-ϕ𝛿italic-ϕdifferential-dΩsubscriptΩsubscript𝜓𝑠italic-ϕ𝛿italic-ϕdifferential-d𝐴subscriptΩsuperscript𝜀italic-ϕ2superscript𝜀2italic-ϕsuperscript𝐃2𝛿italic-ϕdifferential-dΩsubscriptΩdelimited-[]𝜅superscript2italic-ϕsubscript𝜓0italic-ϕsuperscript𝜀italic-ϕ2superscript𝜀2italic-ϕsuperscript𝐃2𝛿italic-ϕdifferential-dΩsubscriptΩ𝜅subscript𝐧𝐰italic-ϕsubscript𝜓𝑠italic-ϕ𝛿italic-ϕdifferential-d𝐴\begin{split}\delta\mathcal{F}&=\int_{\Omega}{(-\kappa}{\nabla^{2}}\phi+\frac{% {\partial{\psi_{0}}}}{{\partial\phi}})\delta\phi d\Omega+\int_{\Omega}\kappa% \nabla\cdot(\nabla\phi)\delta\phi d\Omega+\int_{\partial\Omega}{\frac{{% \partial{{\psi_{s}}}}}{{\partial\phi}}}\delta\phi dA-\int_{\Omega}{\frac{{% \mathop{\varepsilon}^{\prime}(\phi)}}{{2{\varepsilon^{2}}(\phi)}}}|\mathbf{D}{% |^{2}}\delta\phi d\Omega\\ &=\int_{\Omega}{[-\kappa}{\nabla^{2}}\phi+\frac{{\partial{\psi_{0}}}}{{% \partial\phi}}-\frac{{\mathop{\varepsilon}^{\prime}(\phi)}}{{2{\varepsilon^{2}% }(\phi)}}|\mathbf{D}{|^{2}}]\delta\phi d\Omega+\int_{\partial\Omega}{(-\kappa% \mathbf{{n_{w}}}\cdot\nabla\phi+\frac{{\partial{\psi_{s}}}}{{\partial\phi}}})% \delta\phi dA,\end{split}start_ROW start_CELL italic_δ caligraphic_F end_CELL start_CELL = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( - italic_κ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ϕ + divide start_ARG ∂ italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_ϕ end_ARG ) italic_δ italic_ϕ italic_d roman_Ω + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_κ ∇ ⋅ ( ∇ italic_ϕ ) italic_δ italic_ϕ italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT divide start_ARG ∂ italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_ϕ end_ARG italic_δ italic_ϕ italic_d italic_A - ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG italic_ε start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_ϕ ) end_ARG start_ARG 2 italic_ε start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ϕ ) end_ARG | bold_D | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_δ italic_ϕ italic_d roman_Ω end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT [ - italic_κ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ϕ + divide start_ARG ∂ italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_ϕ end_ARG - divide start_ARG italic_ε start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_ϕ ) end_ARG start_ARG 2 italic_ε start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ϕ ) end_ARG | bold_D | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ] italic_δ italic_ϕ italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT ( - italic_κ bold_n start_POSTSUBSCRIPT bold_w end_POSTSUBSCRIPT ⋅ ∇ italic_ϕ + divide start_ARG ∂ italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_ϕ end_ARG ) italic_δ italic_ϕ italic_d italic_A , end_CELL end_ROW (18)

where the Gauss integral theorem has been adopted, and 𝐧wsubscript𝐧𝑤{{\bf{n}}_{w}}bold_n start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT is the inward unit normal vector of the solid surface. Minimizing \mathcal{F}caligraphic_F with respect to the order parameter ϕitalic-ϕ\phiitalic_ϕ yields the following equilibrium conditions,

κ2ϕ+ψ0ϕ𝜀(ϕ)2ε2(ϕ)|𝑫|2=const,inΩ,formulae-sequence𝜅superscript2italic-ϕsubscript𝜓0italic-ϕsuperscript𝜀italic-ϕ2superscript𝜀2italic-ϕsuperscript𝑫2Planck-constant-over-2-piconstinΩ-\kappa{\nabla^{2}}\phi+\frac{{\partial{\rm{}}{\psi_{0}}}}{{\partial\phi}}-% \frac{{\mathop{\varepsilon}^{\prime}(\phi)}}{{2{\varepsilon^{2}}(\phi)}}|\bm{D% }{|^{2}}=\hbar\equiv{\rm{const}},\;\;\;\;{\rm{in}}\;\;\Omega,- italic_κ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ϕ + divide start_ARG ∂ italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_ϕ end_ARG - divide start_ARG italic_ε start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_ϕ ) end_ARG start_ARG 2 italic_ε start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ϕ ) end_ARG | bold_italic_D | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = roman_ℏ ≡ roman_const , roman_in roman_Ω , (19)
κ𝐧wϕ+ψsϕ=¯λ=0,onΩ,formulae-sequence𝜅subscript𝐧𝑤italic-ϕsubscript𝜓𝑠italic-ϕ¯𝜆0onΩ-\kappa{{\bf{n}}_{w}}\cdot\nabla\phi+\frac{{\partial{\psi_{s}}}}{{\partial\phi% }}=\mathchar 22\relax\mkern-10.0mu\lambda=0,\;\;\;\;{\rm{on}}\;\;\partial\Omega,- italic_κ bold_n start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ⋅ ∇ italic_ϕ + divide start_ARG ∂ italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_ϕ end_ARG = ¯ italic_λ = 0 , roman_on ∂ roman_Ω , (20)

where Planck-constant-over-2-pi\hbarroman_ℏ is the chemical potential in the bulk region, and Eq. (20) is the so-called wetting boundary condition. Based on the form of Eq. (18), the rate of change of the total free-energy can be given by

˙=Ω[ϕt]𝑑Ω+Ω[¯λϕt]𝑑A+Ω(𝐄𝐃t)𝑑Ω+Ω(𝒖𝒖t)𝑑Ω+λΩ(qqt)𝑑Ω.˙subscriptΩdelimited-[]Planck-constant-over-2-piitalic-ϕ𝑡differential-dΩsubscriptΩdelimited-[]¯𝜆italic-ϕ𝑡differential-d𝐴subscriptΩ𝐄𝐃𝑡differential-dΩsubscriptΩ𝒖𝒖𝑡differential-dΩ𝜆subscriptΩ𝑞𝑞𝑡differential-dΩ\dot{\mathcal{F}}=\int_{\Omega}{[\hbar\frac{{\partial\phi}}{{\partial t}}]}d% \Omega+\int_{\partial\Omega}[{\mathchar 22\relax\mkern-10.0mu\lambda}\frac{{% \partial\phi}}{{\partial t}}]dA+\int_{\Omega}{(\mathbf{E}\cdot\frac{{\partial% \mathbf{D}}}{{\partial t}})d\Omega+\int_{\Omega}{(\bm{u}\cdot\frac{{\partial% \bm{u}}}{{\partial t}})d\Omega+\lambda\int_{\Omega}{(q\frac{{\partial q}}{{% \partial t}})d\Omega}}}.over˙ start_ARG caligraphic_F end_ARG = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT [ roman_ℏ divide start_ARG ∂ italic_ϕ end_ARG start_ARG ∂ italic_t end_ARG ] italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT [ ¯ italic_λ divide start_ARG ∂ italic_ϕ end_ARG start_ARG ∂ italic_t end_ARG ] italic_d italic_A + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( bold_E ⋅ divide start_ARG ∂ bold_D end_ARG start_ARG ∂ italic_t end_ARG ) italic_d roman_Ω + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( bold_italic_u ⋅ divide start_ARG ∂ bold_italic_u end_ARG start_ARG ∂ italic_t end_ARG ) italic_d roman_Ω + italic_λ ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_q divide start_ARG ∂ italic_q end_ARG start_ARG ∂ italic_t end_ARG ) italic_d roman_Ω . (21)

With Eq. (5) and Eq. (6), we note that

Ω[(ϕt+𝐮ϕ)]𝑑Ω=Ω𝐉ϕ𝑑Ω=Ω𝐉ϕdΩ,subscriptΩdelimited-[]Planck-constant-over-2-piitalic-ϕ𝑡𝐮italic-ϕdifferential-dΩsubscriptΩPlanck-constant-over-2-pisubscript𝐉italic-ϕdifferential-dΩsubscriptΩPlanck-constant-over-2-pisubscript𝐉italic-ϕ𝑑Ω\int_{\Omega}{\left[{\hbar\left({\frac{{\partial\phi}}{{\partial t}}+{\bf{u}}% \cdot\nabla\phi}\right)}\right]}d\Omega=-\int_{\Omega}{\hbar\nabla\cdot{{\bf{J% }}_{\phi}}}d\Omega=\int_{\Omega}{\nabla\hbar\cdot{{\bf{J}}_{\phi}}}d\Omega,∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT [ roman_ℏ ( divide start_ARG ∂ italic_ϕ end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ italic_ϕ ) ] italic_d roman_Ω = - ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT roman_ℏ ∇ ⋅ bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT italic_d roman_Ω = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ∇ roman_ℏ ⋅ bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT italic_d roman_Ω , (22)
Ω[q(qt+𝐮q)]𝑑Ω=Ωq𝐉D𝑑Ω=Ωq𝐉DdΩ,subscriptΩdelimited-[]𝑞𝑞𝑡𝐮𝑞differential-dΩsubscriptΩ𝑞subscript𝐉𝐷differential-dΩsubscriptΩ𝑞subscript𝐉𝐷𝑑Ω\int_{\Omega}{\left[{q\left({\frac{{\partial q}}{{\partial t}}+{\bf{u}}\cdot% \nabla q}\right)}\right]}d\Omega=-\int_{\Omega}{q\nabla\cdot{{\bf{J}}_{D}}}d% \Omega=\int_{\Omega}{\nabla q\cdot{{\bf{J}}_{D}}}d\Omega,∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT [ italic_q ( divide start_ARG ∂ italic_q end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ italic_q ) ] italic_d roman_Ω = - ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_q ∇ ⋅ bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT italic_d roman_Ω = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ∇ italic_q ⋅ bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT italic_d roman_Ω , (23)

where the impermeability conditions 𝐧w𝐉ϕ=0subscript𝐧𝑤subscript𝐉italic-ϕ0{{\bf{n}}_{w}}\cdot{{\bf{J}}_{\phi}}=0bold_n start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ⋅ bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = 0, 𝐧w𝐉D=0subscript𝐧𝑤subscript𝐉𝐷0{{\bf{n}}_{w}}\cdot{{\bf{J}}_{D}}=0bold_n start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ⋅ bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT = 0 at the solid surface have been utilized. Submitting Eq. (22) into Eq. (21), we obtain

˙=Ω[𝐉ϕ𝐮ϕ]𝑑Ω+Ω[¯λ(ϕ𝐮τ(ϕ)τ)]𝑑A+Ω𝐄(q𝐮𝐉𝐃)𝑑Ω+Ω(𝒖𝐅𝚷:𝐮)dΩ+Ω[𝐧𝐰(𝚷𝐮τ)]dA+λΩ[q𝐉𝐃q𝒖q]dΩ,\begin{split}\dot{\mathcal{F}}&=\int_{\Omega}{[\nabla\hbar\cdot{\mathbf{J}_{% \phi}}-\hbar\mathbf{u}\cdot\nabla\phi]}d\Omega+\int_{\partial\Omega}[{% \mathchar 22\relax\mkern-10.0mu\lambda}(\mathop{\phi}\limits^{\cdot}-\mathbf{u% }_{\tau}\cdot{(\nabla\phi)}_{\tau})]dA+\int_{\Omega}{\mathbf{E}\cdot(-q\mathbf% {u}-\mathbf{{J_{D}}})d\Omega}\\ &\qquad+\int_{\Omega}{(\bm{u}\cdot\mathbf{F}-\mathbf{\Pi}:\nabla\mathbf{u})d% \Omega+\int_{\Omega}[\mathbf{n_{w}}\cdot({\mathbf{\Pi}}\mathbf{u}_{\tau})]dA+% \lambda\int_{\Omega}{[\nabla q\cdot{\rm{}}\mathbf{{J_{D}}}-q\bm{u}\cdot\nabla q% ]d\Omega}},\end{split}start_ROW start_CELL over˙ start_ARG caligraphic_F end_ARG end_CELL start_CELL = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT [ ∇ roman_ℏ ⋅ bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT - roman_ℏ bold_u ⋅ ∇ italic_ϕ ] italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT [ ¯ italic_λ ( italic_ϕ start_POSTSUPERSCRIPT ⋅ end_POSTSUPERSCRIPT - bold_u start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ⋅ ( ∇ italic_ϕ ) start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ) ] italic_d italic_A + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT bold_E ⋅ ( - italic_q bold_u - bold_J start_POSTSUBSCRIPT bold_D end_POSTSUBSCRIPT ) italic_d roman_Ω end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( bold_italic_u ⋅ bold_F - bold_Π : ∇ bold_u ) italic_d roman_Ω + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT [ bold_n start_POSTSUBSCRIPT bold_w end_POSTSUBSCRIPT ⋅ ( bold_Π bold_u start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT ) ] italic_d italic_A + italic_λ ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT [ ∇ italic_q ⋅ bold_J start_POSTSUBSCRIPT bold_D end_POSTSUBSCRIPT - italic_q bold_italic_u ⋅ ∇ italic_q ] italic_d roman_Ω , end_CELL end_ROW (24)

where τ𝜏{\bf{\tau}}italic_τ represents the direction tangent to the solid surface.

Assuming that the no-slip boundary holds at the solid surface and when the system is away from its equilibrium state, the energy dissipation in the system arises from the diffusive currents 𝐉ϕsubscript𝐉italic-ϕ{{{\bf{J}}_{\phi}}}bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT, 𝐉Dsubscript𝐉𝐷{{{\bf{J}}_{D}}}bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT in the bulk region and the material time derivative of ϕitalic-ϕ\phiitalic_ϕ, defined as ϕ˙=ϕ/t+𝐮ϕ˙italic-ϕitalic-ϕ/𝑡𝐮italic-ϕ\dot{\phi}={{\partial\phi}\mathord{\left/{\vphantom{{\partial\phi}{\partial t}% }}\right.\kern-1.2pt}{\partial t}}+{\bf{u}}\cdot\nabla\phiover˙ start_ARG italic_ϕ end_ARG = ∂ italic_ϕ start_ID / end_ID ∂ italic_t + bold_u ⋅ ∇ italic_ϕ, at the solid surface. Consequently, the dissipation function ΦF(α˙i,α˙)subscriptΦ𝐹subscript˙𝛼𝑖˙𝛼{\Phi_{F}}\left({{{\dot{\alpha}}_{i}},\dot{\alpha}}\right)roman_Φ start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT ( over˙ start_ARG italic_α end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , over˙ start_ARG italic_α end_ARG ) for the two-phase EHD flows can be given as follows [52, 58, 53]

ΦF=Ω|𝚷|22μ𝑑Ω+Ω|𝐉ϕ|22M𝑑Ω+Ω|𝐉D|22σ𝑑Ω+Ωϕ˙22Γ𝑑A,subscriptΦ𝐹subscriptΩsuperscript𝚷22𝜇differential-dΩsubscriptΩsuperscriptsubscript𝐉italic-ϕ22𝑀differential-dΩsubscriptΩsuperscriptsubscript𝐉𝐷22𝜎differential-dΩsubscriptΩsuperscript˙italic-ϕ22Γdifferential-d𝐴{\Phi_{F}}=\int_{\Omega}{\frac{{{{\left|{\bf{\Pi}}\right|}^{2}}}}{{2\mu}}}d% \Omega+\int_{\Omega}{\frac{{{{\left|{{{\bf{J}}_{\phi}}}\right|}^{2}}}}{{2M}}}d% \Omega+\int_{\Omega}{\frac{{{{\left|{{{\bf{J}}_{D}}}\right|}^{2}}}}{{2\sigma}}% }d\Omega+\int_{\partial\Omega}{\frac{{{{\dot{\phi}}^{2}}}}{{2\Gamma}}}dA,roman_Φ start_POSTSUBSCRIPT italic_F end_POSTSUBSCRIPT = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_Π | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ end_ARG italic_d roman_Ω + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_M end_ARG italic_d roman_Ω + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_σ end_ARG italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT divide start_ARG over˙ start_ARG italic_ϕ end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 roman_Γ end_ARG italic_d italic_A , (25)

where M𝑀Mitalic_M is the mobility and ΓΓ\Gammaroman_Γ is a phenomenological parameter. This relation states that the shear viscosity in the bulk region (first term), composition diffusion in the bulk region (second and third terms) and composition relaxation at the solid surface (fourth term) all contribute to the changes in the energy dissipation. Apparently, by submitting Eq. (25) and Eq. (24) into Eq. (10), we can immediately obtain the expression for the Rayleighian as

=Ω(ψ0+κ2|ϕ|2)𝑑Ω+Ωψs𝑑Ω+12Ω|𝐃|2ε𝑑Ω+12Ωρ|𝐮|2𝑑Ω+λ2Ωq2𝑑Ω+Ω|𝚷|22μ𝑑Ω+Ω|𝐉ϕ|22M𝑑Ω+Ω|𝐉D|22σ𝑑Ω+Ωϕ˙22Γ𝑑A.subscriptΩsubscript𝜓0𝜅2superscriptitalic-ϕ2differential-dΩsubscriptΩsubscript𝜓𝑠differential-dΩ12subscriptΩsuperscript𝐃2𝜀differential-dΩ12subscriptΩ𝜌superscript𝐮2differential-dΩ𝜆2subscriptΩsuperscript𝑞2differential-dΩsubscriptΩsuperscript𝚷22𝜇differential-dΩsubscriptΩsuperscriptsubscript𝐉italic-ϕ22𝑀differential-dΩsubscriptΩsuperscriptsubscript𝐉𝐷22𝜎differential-dΩsubscriptΩsuperscript˙italic-ϕ22Γdifferential-d𝐴\begin{split}\mathcal{R}&=\int_{\Omega}{\left({{\psi_{0}}+\frac{\kappa}{2}{{% \left|{\nabla\phi}\right|}^{2}}}\right)}d\Omega+\int_{\partial\Omega}{{\psi_{s% }}}d\Omega+\frac{1}{2}\int_{\Omega}{\frac{{{{\left|{\bf{D}}\right|}^{2}}}}{% \varepsilon}}d\Omega+\frac{1}{2}\int_{\Omega}{\rho{{\left|{\bf{u}}\right|}^{2}% }}d\Omega+\frac{\lambda}{2}\int_{\Omega}{{q^{2}}}d\Omega\\ &+\int_{\Omega}{\frac{{{{\left|{\bf{\Pi}}\right|}^{2}}}}{{2\mu}}}d\Omega+\int_% {\Omega}{\frac{{{{\left|{{{\bf{J}}_{\phi}}}\right|}^{2}}}}{{2M}}}d\Omega+\int_% {\Omega}{\frac{{{{\left|{{{\bf{J}}_{D}}}\right|}^{2}}}}{{2\sigma}}}d\Omega+% \int_{\partial\Omega}{\frac{{{{\dot{\phi}}^{2}}}}{{2\Gamma}}}dA.\end{split}start_ROW start_CELL caligraphic_R end_CELL start_CELL = ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + divide start_ARG italic_κ end_ARG start_ARG 2 end_ARG | ∇ italic_ϕ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_d roman_Ω + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_D | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_ε end_ARG italic_d roman_Ω + divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_ρ | bold_u | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d roman_Ω + divide start_ARG italic_λ end_ARG start_ARG 2 end_ARG ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_q start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d roman_Ω end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_Π | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_μ end_ARG italic_d roman_Ω + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_M end_ARG italic_d roman_Ω + ∫ start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT divide start_ARG | bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_σ end_ARG italic_d roman_Ω + ∫ start_POSTSUBSCRIPT ∂ roman_Ω end_POSTSUBSCRIPT divide start_ARG over˙ start_ARG italic_ϕ end_ARG start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 roman_Γ end_ARG italic_d italic_A . end_CELL end_ROW (26)

Thereafter, minimizing \mathcal{R}caligraphic_R with respect to the rates 𝐮𝐮\bf{u}bold_u, 𝐉ϕsubscript𝐉italic-ϕ{{{\bf{J}}_{\phi}}}bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT and 𝐉Dsubscript𝐉𝐷{{\bf{J}}_{D}}bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT under the incompressible condition given by Eq. (1), we obtain

𝐅=ϕ+q𝐄12(ασq2),𝐅Planck-constant-over-2-piitalic-ϕ𝑞𝐄12𝛼𝜎superscript𝑞2{\bf{F}}=\hbar\nabla\phi+q{\bf{E}}-\frac{1}{2}\nabla\left({\frac{\alpha}{% \sigma}{q^{2}}}\right),bold_F = roman_ℏ ∇ italic_ϕ + italic_q bold_E - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ∇ ( divide start_ARG italic_α end_ARG start_ARG italic_σ end_ARG italic_q start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) , (27)
𝐉ϕ=M,𝐉D=σ(𝐄λq).formulae-sequencesubscript𝐉italic-ϕ𝑀Planck-constant-over-2-pisubscript𝐉𝐷𝜎𝐄𝜆𝑞{{\bf{J}}_{\phi}}=-M\nabla\hbar,\;\;\;\;\;{{\bf{J}}_{D}}=\sigma\left({{\bf{E}}% -\lambda\nabla q}\right).bold_J start_POSTSUBSCRIPT italic_ϕ end_POSTSUBSCRIPT = - italic_M ∇ roman_ℏ , bold_J start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT = italic_σ ( bold_E - italic_λ ∇ italic_q ) . (28)

By submitting Eq. (27) and Eq. (28) into the generalized equations described in Sect. 2.1, we obtain the system of system equations for the two-phase EHD flows maintaining thermodynamic consistency,

𝐮=0,𝐮0\nabla\cdot{\bf{u}}=0,∇ ⋅ bold_u = 0 , (29)
ρ(𝐮t+𝐮𝐮)=p^+μ2𝐮+ϕ+q𝐄,𝜌𝐮𝑡𝐮𝐮^𝑝𝜇superscript2𝐮Planck-constant-over-2-piitalic-ϕ𝑞𝐄\rho\left({\frac{{\partial{\bf{u}}}}{{\partial t}}+{\bf{u}}\cdot\nabla{\bf{u}}% }\right)=-\nabla\hat{p}+\mu{\nabla^{2}}{\bf{u}}+\hbar\nabla\phi+q{\bf{E}},italic_ρ ( divide start_ARG ∂ bold_u end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ bold_u ) = - ∇ over^ start_ARG italic_p end_ARG + italic_μ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_u + roman_ℏ ∇ italic_ϕ + italic_q bold_E , (30)
ϕt+𝐮ϕ=M,=4βϕ(ϕ12)(ϕ1)κ2ϕ12ε(ϕ)|𝐄|2,formulae-sequenceitalic-ϕ𝑡𝐮italic-ϕ𝑀Planck-constant-over-2-piPlanck-constant-over-2-pi4𝛽italic-ϕitalic-ϕ12italic-ϕ1𝜅superscript2italic-ϕ12superscript𝜀italic-ϕsuperscript𝐄2\frac{{\partial\phi}}{{\partial t}}+{\bf{u}}\cdot\nabla\phi=\nabla\cdot M% \nabla\hbar,\;\;\;\;\;\hbar=4\beta\phi\left({\phi-\frac{1}{2}}\right)\left({% \phi-1}\right)-\kappa{\nabla^{2}}\phi-\frac{1}{2}\varepsilon^{\prime}\left(% \phi\right){\left|{\bf{E}}\right|^{2}},divide start_ARG ∂ italic_ϕ end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ italic_ϕ = ∇ ⋅ italic_M ∇ roman_ℏ , roman_ℏ = 4 italic_β italic_ϕ ( italic_ϕ - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ) ( italic_ϕ - 1 ) - italic_κ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ϕ - divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_ε start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_ϕ ) | bold_E | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (31)
(ε𝐄)=q,𝐄=φ,formulae-sequence𝜀𝐄𝑞𝐄𝜑\nabla\cdot\left({\varepsilon{\bf{E}}}\right)=q,\;\;\;\;\;{\rm{}}{\bf{E}}=-% \nabla\varphi,∇ ⋅ ( italic_ε bold_E ) = italic_q , bold_E = - ∇ italic_φ , (32)
qt+𝐮q=[αq+σφ].𝑞𝑡𝐮𝑞delimited-[]𝛼𝑞𝜎𝜑\frac{{\partial q}}{{\partial t}}+{\bf{u}}\cdot\nabla q=\nabla\cdot\left[{% \alpha\nabla q+\sigma\nabla\varphi}\right].divide start_ARG ∂ italic_q end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ italic_q = ∇ ⋅ [ italic_α ∇ italic_q + italic_σ ∇ italic_φ ] . (33)

To conclude this section, certain remarks on the current physical model are summarized below:

Remark I: Based on the above governing equations, it is noted that the third term in the chemical potential 0.5ε(ϕ)|𝐄|20.5superscript𝜀italic-ϕsuperscript𝐄20.5\varepsilon^{\prime}\left(\phi\right){\left|{\bf{E}}\right|^{2}}0.5 italic_ε start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_ϕ ) | bold_E | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT is associated with the presence of electric field, which suggests that when one adopts phase-field based method to model two-phase EHD flows, the conventional chemical potential (i.e., 4βϕ(ϕ0.5)(ϕ1)κ2ϕ4𝛽italic-ϕitalic-ϕ0.5italic-ϕ1𝜅superscript2italic-ϕ4\beta\phi\left({\phi-0.5}\right)\left({\phi-1}\right)-\kappa{\nabla^{2}}\phi4 italic_β italic_ϕ ( italic_ϕ - 0.5 ) ( italic_ϕ - 1 ) - italic_κ ∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ϕ ) used in traditional phase-field method is not applicable theoretically. Similar term was also obtained by Eck et al. [58] and Yang et al. [61], who used the phase-field-based numerical methods to study the problem of electrowetting on dielectric.

Remark II: Reduction consistency is a crucial property for phase-field modelling [60]. This means that for a two-phase system, if one fluid component is absent, then the governing equations for the two-phase system should accurately reduce to those of the corresponding single-phase system. With reduction consistency in mind, the permittivity of ε𝜀\varepsilonitalic_ε in a two-phase EHD flow must adhere to the following constraint:

ε(ϕ)|=ϕ=00,ε(ϕ)|=ϕ=10.\varepsilon^{\prime}\left(\phi\right)\left|{{}_{\phi=0}}\right.=0,\;\;\;\;\;\;% \varepsilon^{\prime}\left(\phi\right)\left|{{}_{\phi=1}}\right.=0.italic_ε start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_ϕ ) | start_FLOATSUBSCRIPT italic_ϕ = 0 end_FLOATSUBSCRIPT = 0 , italic_ε start_POSTSUPERSCRIPT ′ end_POSTSUPERSCRIPT ( italic_ϕ ) | start_FLOATSUBSCRIPT italic_ϕ = 1 end_FLOATSUBSCRIPT = 0 . (34)

In such a case, accounting for the unequal permittivity of the two fluids, the Hermite interpolation is employed to handle the permittivity across the fluid-fluid interface,

ε(ϕ)=εv+ϕ2(2ϕ3)(εvεl),𝜀italic-ϕsubscript𝜀𝑣superscriptitalic-ϕ22italic-ϕ3subscript𝜀𝑣subscript𝜀𝑙\varepsilon\left(\phi\right)={\varepsilon_{v}}+{\phi^{2}}\left({2\phi-3}\right% )\left({{\varepsilon_{v}}-{\varepsilon_{l}}}\right),italic_ε ( italic_ϕ ) = italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT + italic_ϕ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( 2 italic_ϕ - 3 ) ( italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT ) , (35)

where the subscripts v𝑣vitalic_v and l𝑙litalic_l denote the vapour and liquid phases, respectively. The other physical variables including density ρ𝜌\rhoitalic_ρ, dynamic viscosity μ𝜇\muitalic_μ, and conductivity σ𝜎\sigmaitalic_σ are all determined using a simple linear function of the order parameter

ρ(ϕ)=ϕ(ρlρv)+ρv,μ(ϕ)=ϕ(μlμv)+μv,σ(ϕ)=ϕ(σlσv)+σv.formulae-sequence𝜌italic-ϕitalic-ϕsubscript𝜌𝑙subscript𝜌𝑣subscript𝜌𝑣formulae-sequence𝜇italic-ϕitalic-ϕsubscript𝜇𝑙subscript𝜇𝑣subscript𝜇𝑣𝜎italic-ϕitalic-ϕsubscript𝜎𝑙subscript𝜎𝑣subscript𝜎𝑣\rho\left(\phi\right)=\phi\left({{\rho_{l}}-{\rho_{v}}}\right)+{\rho_{v}},\;\;% \;\;\mu\left(\phi\right)=\phi\left({{\mu_{l}}-{\mu_{v}}}\right)+{\mu_{v}},\;\;% \;\;\sigma\left(\phi\right)=\phi\left({{\sigma_{l}}-{\sigma_{v}}}\right)+{% \sigma_{v}}.italic_ρ ( italic_ϕ ) = italic_ϕ ( italic_ρ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ) + italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT , italic_μ ( italic_ϕ ) = italic_ϕ ( italic_μ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ) + italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT , italic_σ ( italic_ϕ ) = italic_ϕ ( italic_σ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ) + italic_σ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT . (36)

3 Lattice Boltzmann method

Unlike classical computational fluid dynamics approaches based on the discretization of the macroscopic equations, the LB method is based on the lattice Boltzmann equation, which describes the dynamic evolution process of the distribution function in the discrete velocity space [43, 44]. For two-dimensional cases, the two-dimensional nine-velocity (D2Q9) discrete lattice velocity is commonly adopted [44], which is given by

𝐜i={c(0,0),i=0,c(cos[(i1)π/2],sin[(i1)π/2]),i=1,2,3,4,2c(cos[(2i9)π/4],sin[(2i9)π/4]),i=5,6,7,8,subscript𝐜𝑖cases𝑐00𝑖0formulae-sequence𝑐𝑖1𝜋/2𝑖1𝜋/2𝑖1234formulae-sequence2𝑐2𝑖9𝜋/42𝑖9𝜋/4𝑖5678{{\bf{c}}_{i}}=\left\{\begin{array}[]{l}c\left({0,0}\right),\;\;\;\;\;\;\;\;\;% \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;% \;\;\;\;\;\;\;\;i=0,\\ c\left({\cos\left[{{{\left({i-1}\right)\pi}\mathord{\left/{\vphantom{{\left({i% -1}\right)\pi}2}}\right.\kern-1.2pt}2}}\right],\sin\left[{{{\left({i-1}\right)% \pi}\mathord{\left/{\vphantom{{\left({i-1}\right)\pi}2}}\right.\kern-1.2pt}2}}% \right]}\right),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i=1,2,3,4,\\ \sqrt{2}c\left({\cos\left[{{{\left({2i-9}\right)\pi}\mathord{\left/{\vphantom{% {\left({2i-9}\right)\pi}4}}\right.\kern-1.2pt}4}}\right],\sin\left[{{{\left({2% i-9}\right)\pi}\mathord{\left/{\vphantom{{\left({2i-9}\right)\pi}4}}\right.% \kern-1.2pt}4}}\right]}\right),\;\;\;\;\;\;i=5,6,7,8,\end{array}\right.bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = { start_ARRAY start_ROW start_CELL italic_c ( 0 , 0 ) , italic_i = 0 , end_CELL end_ROW start_ROW start_CELL italic_c ( roman_cos [ ( italic_i - 1 ) italic_π start_ID / end_ID 2 ] , roman_sin [ ( italic_i - 1 ) italic_π start_ID / end_ID 2 ] ) , italic_i = 1 , 2 , 3 , 4 , end_CELL end_ROW start_ROW start_CELL square-root start_ARG 2 end_ARG italic_c ( roman_cos [ ( 2 italic_i - 9 ) italic_π start_ID / end_ID 4 ] , roman_sin [ ( 2 italic_i - 9 ) italic_π start_ID / end_ID 4 ] ) , italic_i = 5 , 6 , 7 , 8 , end_CELL end_ROW end_ARRAY (37)

where i𝑖iitalic_i is the lattice direction, c=Δx/Δt𝑐Δ𝑥/Δ𝑡c={{\Delta x}\mathord{\left/{\vphantom{{\Delta x}{\Delta t}}}\right.\kern-1.2% pt}{\Delta t}}italic_c = roman_Δ italic_x start_ID / end_ID roman_Δ italic_t is the lattice speed with ΔxΔ𝑥{\Delta x}roman_Δ italic_x and ΔtΔ𝑡{\Delta t}roman_Δ italic_t being the lattice spacing and time step. In addition, the sound speed cssubscript𝑐𝑠c_{s}italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT in D2Q9 model is expressed as cs=c/3subscript𝑐𝑠𝑐/3{c_{s}}={c\mathord{\left/{\vphantom{c{\sqrt{3}}}}\right.\kern-1.2pt}{\sqrt{3}}}italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = italic_c start_ID / end_ID square-root start_ARG 3 end_ARG, and the weight coefficients in different lattice directions are given by ω0=4/9,ω14=1/9,ω58=1/36formulae-sequencesubscript𝜔04/9formulae-sequencesubscript𝜔141/9subscript𝜔581/36{\omega_{0}}={4\mathord{\left/{\vphantom{49}}\right.\kern-1.2pt}9},\;\;{\omega% _{1-4}}={1\mathord{\left/{\vphantom{19}}\right.\kern-1.2pt}9},\;\;{\omega_{5-8% }}={1\mathord{\left/{\vphantom{1{36}}}\right.\kern-1.2pt}{36}}italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 4 start_ID / end_ID 9 , italic_ω start_POSTSUBSCRIPT 1 - 4 end_POSTSUBSCRIPT = 1 start_ID / end_ID 9 , italic_ω start_POSTSUBSCRIPT 5 - 8 end_POSTSUBSCRIPT = 1 start_ID / end_ID 36.

3.1 Lattice Boltzmann method for hydrodynamic equations

The LB equation for incompressible fluid flow can be written as [48]

𝑓i(𝐱+𝐜iΔt,t+Δt)=𝑓i(𝐱,t)1𝜏f[𝑓i(𝐱,t)𝑓ieq(𝐱,t)]+Δt(112𝜏f)R¯i(𝐱,t),subscript𝑓𝑖𝐱subscript𝐜𝑖Δ𝑡𝑡Δ𝑡subscript𝑓𝑖𝐱𝑡1subscript𝜏𝑓delimited-[]subscript𝑓𝑖𝐱𝑡superscriptsubscript𝑓𝑖𝑒𝑞𝐱𝑡Δ𝑡112subscript𝜏𝑓subscript¯𝑅𝑖𝐱𝑡\mathop{f}\nolimits_{i}({\bf{x}}+{{\bf{c}}_{i}}\Delta t,t+\Delta t)=\mathop{f}% \nolimits_{i}({\bf{x}},t)-\displaystyle\frac{1}{{\mathop{\tau}\nolimits_{f}}}[% \mathop{f}\nolimits_{i}({\bf{x}},t)-\mathop{f}\nolimits_{i}^{eq}({\bf{x}},t)]+% \Delta t(1-\frac{1}{{2\mathop{\tau}\nolimits_{f}}})\mathop{\bar{R}}\nolimits_{% i}({\bf{x}},t),italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x + bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_Δ italic_t , italic_t + roman_Δ italic_t ) = italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT end_ARG [ italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) ] + roman_Δ italic_t ( 1 - divide start_ARG 1 end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT end_ARG ) start_BIGOP over¯ start_ARG italic_R end_ARG end_BIGOP start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) , (38)

where 𝑓i(𝐱,t)subscript𝑓𝑖𝐱𝑡\mathop{f}\nolimits_{i}({\bf{x}},t)italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) denotes the distribution function at position 𝐱𝐱{\bf{x}}bold_x and time t𝑡titalic_t, ΔtΔ𝑡\Delta troman_Δ italic_t is the time step, R¯i(𝐱,t)subscript¯𝑅𝑖𝐱𝑡\mathop{\bar{R}}\nolimits_{i}({\bf{x}},t)start_BIGOP over¯ start_ARG italic_R end_ARG end_BIGOP start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) is the discrete forcing term defined as

R¯i(x,t)=𝜔i[𝐮ρ+𝐜i𝐅𝑐s2+𝐮ρ:(𝐜i𝐜i𝑐s2𝐈)𝑐s2].subscript¯𝑅𝑖𝑥𝑡subscript𝜔𝑖delimited-[]𝐮𝜌subscript𝐜𝑖𝐅superscriptsubscript𝑐𝑠2:𝐮𝜌subscript𝐜𝑖subscript𝐜𝑖superscriptsubscript𝑐𝑠2𝐈superscriptsubscript𝑐𝑠2\mathop{\bar{R}}\nolimits_{i}(x,t)=\mathop{\omega}\nolimits_{i}[\mathbf{u}% \cdot\nabla\rho+\frac{{\mathbf{c}_{i}\cdot\mathbf{F}}}{{\mathop{c}\nolimits_{s% }^{2}}}+\frac{{\mathbf{u}\nabla\rho:(\mathbf{c}_{i}\mathbf{c}_{i}-\mathop{c}% \nolimits_{s}^{2}\mathbf{I})}}{{\mathop{c}\nolimits_{s}^{2}}}].start_BIGOP over¯ start_ARG italic_R end_ARG end_BIGOP start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_x , italic_t ) = italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT [ bold_u ⋅ ∇ italic_ρ + divide start_ARG bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⋅ bold_F end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG bold_u ∇ italic_ρ : ( bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_I ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ] . (39)

𝑓ieq(𝐱,t)superscriptsubscript𝑓𝑖𝑒𝑞𝐱𝑡\mathop{f}\nolimits_{i}^{eq}({\bf{x}},t)italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) is the local equilibrium distribution function given by

𝑓ieq(𝐱,t)={(𝜔01)p𝑐s2+ρ𝑠i(𝐮),i=0,𝜔ip𝑐s2+ρ𝑠i(𝐮),i0,superscriptsubscript𝑓𝑖𝑒𝑞𝐱𝑡casessubscript𝜔01𝑝superscriptsubscript𝑐𝑠2𝜌subscript𝑠𝑖𝐮𝑖0missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionsubscript𝜔𝑖𝑝superscriptsubscript𝑐𝑠2𝜌subscript𝑠𝑖𝐮𝑖0missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpression\mathop{f}\nolimits_{i}^{eq}({\bf{x}},t)=\left\{{\begin{array}[]{*{20}{c}}{(% \mathop{\omega}\nolimits_{0}-1)\displaystyle\frac{p}{{\mathop{c}\nolimits_{s}^% {2}}}+\rho\mathop{s}\nolimits_{i}(\mathbf{u})},&{i=0},\\ {\mathop{\omega}\nolimits_{i}\displaystyle\frac{p}{{\mathop{c}\nolimits_{s}^{2% }}}+\rho\mathop{s}\nolimits_{i}(\mathbf{u})},&{i\neq 0},\end{array}}\right.italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) = { start_ARRAY start_ROW start_CELL ( italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - 1 ) divide start_ARG italic_p end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_ρ italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_u ) , end_CELL start_CELL italic_i = 0 , end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT divide start_ARG italic_p end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + italic_ρ italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_u ) , end_CELL start_CELL italic_i ≠ 0 , end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW end_ARRAY (40)

with

𝑠i(𝐮)=ωi[𝐜i𝐮𝑐s2+(𝐜i𝐮)22𝑐s4𝐮𝐮2cs2].subscript𝑠𝑖𝐮subscript𝜔𝑖delimited-[]subscript𝐜𝑖𝐮superscriptsubscript𝑐𝑠2superscriptsubscript𝐜𝑖𝐮22superscriptsubscript𝑐𝑠4𝐮𝐮superscriptsubscript2𝑐𝑠2\mathop{s}\nolimits_{i}(\mathbf{u})=\omega_{i}[\frac{{\mathbf{c}_{i}\cdot% \mathbf{u}}}{{\mathop{c}\nolimits_{s}^{2}}}+\displaystyle\frac{{\mathop{(% \mathbf{c}_{i}\cdot\mathbf{u})}\nolimits^{2}}}{{2\mathop{c}\nolimits_{s}^{4}}}% -\frac{{\mathbf{u}\cdot\mathbf{u}}}{{\mathop{2c}\nolimits_{s}^{2}}}].italic_s start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_u ) = italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT [ divide start_ARG bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⋅ bold_u end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG start_BIGOP ( bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⋅ bold_u ) end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG - divide start_ARG bold_u ⋅ bold_u end_ARG start_ARG start_BIGOP 2 italic_c end_BIGOP start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ] . (41)

In addition, the relaxation time τfsubscript𝜏𝑓{\tau_{f}}italic_τ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT is related to the density and dynamic viscosity and is defined as

μ=ρcs2(τf12).𝜇𝜌superscriptsubscript𝑐𝑠2subscript𝜏𝑓12\mu=\rho c_{s}^{2}\left({{\tau_{f}}-\displaystyle\frac{1}{2}}\right).italic_μ = italic_ρ italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_τ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ) . (42)

With the distribution function, the hydrodynamic pressure p𝑝pitalic_p and velocity 𝐮𝐮\bf{u}bold_u are calculated by

ρ𝐮=i𝐜ifi+12Δt𝐅,𝜌𝐮subscript𝑖subscript𝐜𝑖subscript𝑓𝑖12Δ𝑡𝐅\rho\mathbf{u}=\sum\limits_{i}{\mathbf{c}_{i}{f_{i}}}+\displaystyle\frac{1}{2}% \Delta t\mathbf{F},italic_ρ bold_u = ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Δ italic_t bold_F , (43)
p=𝑐s21𝜔0[i0𝑓i+Δt2𝐮ρ+ρ𝑠0(𝐮)].𝑝superscriptsubscript𝑐𝑠21subscript𝜔0delimited-[]subscript𝑖0subscript𝑓𝑖Δ𝑡2𝐮𝜌𝜌subscript𝑠0𝐮p=\frac{{\mathop{c}\nolimits_{s}^{2}}}{{1-\mathop{\omega}\nolimits_{0}}}[\sum% \limits_{i\neq 0}{\mathop{f}\nolimits_{i}}+\displaystyle\frac{{\Delta t}}{2}% \mathbf{u}\cdot\nabla\rho+\rho\mathop{s}\nolimits_{0}(\mathbf{u})].italic_p = divide start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 1 - italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG [ ∑ start_POSTSUBSCRIPT italic_i ≠ 0 end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + divide start_ARG roman_Δ italic_t end_ARG start_ARG 2 end_ARG bold_u ⋅ ∇ italic_ρ + italic_ρ italic_s start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( bold_u ) ] . (44)

3.2 Lattice Boltzmann equation for Cahn-Hilliard equation

To reduce the spurious velocities arising from the force imbalance at the discrete level, the present work adopts a well-balanced LB method [62] to solve the Cahn-Hilliard equation, and the corresponding evolution equation reads

𝑔i(𝐱+𝐜iΔt,t+Δt)=𝑔i(𝐱,t)1𝜏g[𝑔i(𝐱,t)𝑔ieq(𝐱,t)]+Δt𝐺i(𝐱,t)+12Δt2t𝐺i(𝐱,t),subscript𝑔𝑖𝐱subscript𝐜𝑖Δ𝑡𝑡Δ𝑡subscript𝑔𝑖𝐱𝑡1subscript𝜏𝑔delimited-[]subscript𝑔𝑖𝐱𝑡superscriptsubscript𝑔𝑖𝑒𝑞𝐱𝑡Δ𝑡subscript𝐺𝑖𝐱𝑡12superscriptΔ𝑡2subscript𝑡subscript𝐺𝑖𝐱𝑡\mathop{g}\nolimits_{i}(\mathbf{x}+\mathbf{c}_{i}\Delta t,t+\Delta t)=\mathop{% g}\nolimits_{i}(\mathbf{x},t)-\displaystyle\frac{1}{{{}_{\mathop{\tau}% \nolimits_{g}}}}[\mathop{g}\nolimits_{i}(\mathbf{x},t)-\mathop{g}\nolimits_{i}% ^{eq}(\mathbf{x},t)]+\Delta t\mathop{G}\nolimits_{i}(\mathbf{x},t)+\frac{1}{2}% \mathop{\Delta t}\nolimits^{2}\mathop{\partial}\nolimits_{t}\mathop{G}% \nolimits_{i}(\mathbf{x},t),italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x + bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_Δ italic_t , italic_t + roman_Δ italic_t ) = italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - divide start_ARG 1 end_ARG start_ARG start_FLOATSUBSCRIPT italic_τ start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT end_FLOATSUBSCRIPT end_ARG [ italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) ] + roman_Δ italic_t italic_G start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) + divide start_ARG 1 end_ARG start_ARG 2 end_ARG start_BIGOP roman_Δ italic_t end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_G start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) , (45)

with 𝑔ieq(𝐱,t)superscriptsubscript𝑔𝑖𝑒𝑞𝐱𝑡\mathop{g}\nolimits_{i}^{eq}(\mathbf{x},t)italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) being the equilibrium distribution function given by

gieq(𝐱,t)={ϕ(1ω0),i=0,ωi,i0.superscriptsubscript𝑔𝑖𝑒𝑞𝐱𝑡casesitalic-ϕ1subscript𝜔0Planck-constant-over-2-pi𝑖0subscript𝜔𝑖Planck-constant-over-2-pi𝑖0g_{i}^{eq}\left({{\bf{x}},t}\right)=\left\{\begin{array}[]{l}\phi-\left({1-{% \omega_{0}}}\right)\hbar,\;\;i=0,\\ {\omega_{i}}\hbar,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i\neq 0.\end{array}\right.italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) = { start_ARRAY start_ROW start_CELL italic_ϕ - ( 1 - italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) roman_ℏ , italic_i = 0 , end_CELL end_ROW start_ROW start_CELL italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_ℏ , italic_i ≠ 0 . end_CELL end_ROW end_ARRAY (46)

The discrete forcing term 𝐺i(𝐱,t)subscript𝐺𝑖𝐱𝑡\mathop{G}\nolimits_{i}(\mathbf{x},t)italic_G start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) is expressed as

𝐺i(𝐱,t)=𝜔i(𝐮ϕ)[1+𝐈:(𝐜i𝐜i𝑐s2𝐈)2𝑐s2].subscript𝐺𝑖𝐱𝑡subscript𝜔𝑖𝐮italic-ϕdelimited-[]1:𝐈subscript𝐜𝑖subscript𝐜𝑖superscriptsubscript𝑐𝑠2𝐈2superscriptsubscript𝑐𝑠2\mathop{G}\nolimits_{i}(\mathbf{x},t)=\mathop{\omega}\nolimits_{i}(\mathbf{u}% \cdot\nabla\phi)[-1+\frac{{\mathbf{I}:(\mathbf{c}_{i}\mathbf{c}_{i}-\mathop{c}% \nolimits_{s}^{2}\mathbf{I})}}{{2\mathop{c}\nolimits_{s}^{2}}}].italic_G start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) = italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_u ⋅ ∇ italic_ϕ ) [ - 1 + divide start_ARG bold_I : ( bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT bold_I ) end_ARG start_ARG 2 italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ] . (47)

In addition, the mobility M𝑀Mitalic_M is expressed as a function of relaxation time τgsubscript𝜏𝑔\tau_{g}italic_τ start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT through

M=𝑐s2(𝜏g12)Δt,𝑀superscriptsubscript𝑐𝑠2subscript𝜏𝑔12Δ𝑡M=\mathop{c}\nolimits_{s}^{2}(\mathop{\tau}\nolimits_{g}-\frac{1}{2})\Delta t,italic_M = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_τ start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ) roman_Δ italic_t , (48)

and the order parameter in this model is calculated by

ϕ=i𝑔i.italic-ϕsubscript𝑖subscript𝑔𝑖\phi=\sum\limits_{i}{\mathop{g}\nolimits_{i}}.italic_ϕ = ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT . (49)

Note that there some derivatives appeared in the evolution equation, which need to be numerically calculated in the implementation. In such a case, the backward difference scheme is adopted to compute the time derivative,

t(Θ)=Θ(𝐱,t)Θ(𝐱,tΔt)Δt.subscript𝑡ΘΘ𝐱𝑡Θ𝐱𝑡Δ𝑡Δ𝑡\mathop{\partial}\nolimits_{t}(\Theta)=\frac{{\Theta(\mathbf{x},t)-\Theta(% \mathbf{x},t-\Delta t)}}{{\Delta t}}.∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( roman_Θ ) = divide start_ARG roman_Θ ( bold_x , italic_t ) - roman_Θ ( bold_x , italic_t - roman_Δ italic_t ) end_ARG start_ARG roman_Δ italic_t end_ARG . (50)

Moreover, to maintain the numerical accuracy of our LB method, we use conventional second-order isotropic discretization schemes to calculate the gradient and Laplace operator terms [63],

Θ(𝐱)=i0𝜔i𝐜iΘ(𝐱+𝐜iΔt)𝑐s2Δt,Θ𝐱subscript𝑖0subscript𝜔𝑖subscript𝐜𝑖Θ𝐱subscript𝐜𝑖Δ𝑡superscriptsubscript𝑐𝑠2Δ𝑡\nabla\Theta(\mathbf{x})=\sum\limits_{i\neq 0}{\displaystyle\frac{{\mathop{% \omega}\nolimits_{i}\mathbf{c}_{i}\Theta(\mathbf{x}+\mathbf{c}_{i}\Delta t)}}{% {\mathop{c}\nolimits_{s}^{2}\Delta t}}},∇ roman_Θ ( bold_x ) = ∑ start_POSTSUBSCRIPT italic_i ≠ 0 end_POSTSUBSCRIPT divide start_ARG italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_Θ ( bold_x + bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_Δ italic_t ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_t end_ARG , (51)
2Θ(𝐱)=i02ωi[Θ(𝐱+𝐜iΔt)Θ(𝐱)]𝑐s2Δt2,superscript2Θ𝐱subscript𝑖0subscript2𝜔𝑖delimited-[]Θ𝐱subscript𝐜𝑖Δ𝑡Θ𝐱superscriptsubscript𝑐𝑠2superscriptΔ𝑡2\mathop{\nabla}\nolimits^{2}\Theta(\mathbf{x})=\sum\limits_{i\neq 0}{% \displaystyle\frac{{\mathop{2\omega}\nolimits_{i}[\Theta(\mathbf{x}+\mathbf{c}% _{i}\Delta t)-\Theta(\mathbf{x})]}}{{\mathop{c}\nolimits_{s}^{2}\mathop{\Delta t% }\nolimits^{2}}}},∇ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Θ ( bold_x ) = ∑ start_POSTSUBSCRIPT italic_i ≠ 0 end_POSTSUBSCRIPT divide start_ARG start_BIGOP 2 italic_ω end_BIGOP start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT [ roman_Θ ( bold_x + bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_Δ italic_t ) - roman_Θ ( bold_x ) ] end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_BIGOP roman_Δ italic_t end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (52)

with ΘΘ\Thetaroman_Θ being any physical variable.

3.3 Lattice Boltzmann equation for electric potential equation

It is noted that the electric potential is actually governed by the Poisson equation, i.e,

ε(ϕ)φ+q=0.𝜀italic-ϕ𝜑𝑞0\nabla\cdot\varepsilon\left(\phi\right)\nabla\varphi+q=0.∇ ⋅ italic_ε ( italic_ϕ ) ∇ italic_φ + italic_q = 0 . (53)

To solve it using the LB method, the model proposed by Chai et al. [64] is typically adopted in the LB community. In such a case, the relaxation time for electric potential is defined as

ε(ϕ)=χcs2(τφ0.5)Δt,𝜀italic-ϕ𝜒superscriptsubscript𝑐𝑠2subscript𝜏𝜑0.5Δ𝑡\varepsilon\left(\phi\right)=\chi c_{s}^{2}\left({{\tau_{\varphi}}-0.5}\right)% \Delta t,italic_ε ( italic_ϕ ) = italic_χ italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_τ start_POSTSUBSCRIPT italic_φ end_POSTSUBSCRIPT - 0.5 ) roman_Δ italic_t , (54)

in which χ𝜒\chiitalic_χ is an artificial parameter that adjusts the relaxation time. Frankly, this scheme performs well when the diffusion coefficient is uniform or when there are only slight differences in diffusion. However, in EHD flows, the permittivity ratio is always large [65, 67], indicating that the ratio of relaxation times between two fluids is also substantial. Thus, the conventional scheme may share slow and inefficient solution procedures [68]. Moreover, it is worth noting that while an artificial parameter could adjust the relaxation time, using inappropriate values may result in unstable outcomes [69]. In this setting, we aim to develop an improved method for solving the electric potential equation with a large permittivity ratio. With the definition of the variable permittivity, the electric potential equation can be rewritten as

εvφ+ε^(ϕ)φ+q=0,subscript𝜀𝑣𝜑^𝜀italic-ϕ𝜑𝑞0\nabla\cdot{\varepsilon_{v}}\nabla\varphi+\nabla\cdot\hat{\varepsilon}\left(% \phi\right)\nabla\varphi+q=0,∇ ⋅ italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ∇ italic_φ + ∇ ⋅ over^ start_ARG italic_ε end_ARG ( italic_ϕ ) ∇ italic_φ + italic_q = 0 , (55)

where ε^(ϕ)=ϕ2(2ϕ3)(εvεl)^𝜀italic-ϕsuperscriptitalic-ϕ22italic-ϕ3subscript𝜀𝑣subscript𝜀𝑙\hat{\varepsilon}\left(\phi\right)={\phi^{2}}\left({2\phi-3}\right)\left({{% \varepsilon_{v}}-{\varepsilon_{l}}}\right)over^ start_ARG italic_ε end_ARG ( italic_ϕ ) = italic_ϕ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( 2 italic_ϕ - 3 ) ( italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT - italic_ε start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT ). Based on the above recast equation, the LB equation for electric potential can be expressed as

hi(𝐱+𝐜iΔt,t+Δt)=hi(𝐱,t)1τh[hi(𝐱,t)hieq(𝐱,t)]ε^(ϕ)𝑐s2𝜔i𝐜iφτh+Δtϖiq,subscript𝑖𝐱subscript𝐜𝑖Δ𝑡𝑡Δ𝑡subscript𝑖𝐱𝑡1subscript𝜏delimited-[]subscript𝑖𝐱𝑡superscriptsubscript𝑖𝑒𝑞𝐱𝑡^𝜀italic-ϕsuperscriptsubscript𝑐𝑠2subscript𝜔𝑖subscript𝐜𝑖𝜑subscript𝜏Δ𝑡subscriptitalic-ϖ𝑖𝑞{h_{i}}(\mathbf{x}+{{\bf{c}}_{i}}\Delta t,t+\Delta t)={h_{i}}(\mathbf{x},t)-% \frac{1}{{{\tau_{h}}}}[{h_{i}}(\mathbf{x},t)-h_{i}^{eq}(\mathbf{x},t)]-\frac{% \hat{\varepsilon}\left(\phi\right)}{{\mathop{c}\nolimits_{s}^{2}}}\frac{{% \mathop{\omega}\nolimits_{i}\mathbf{c}_{i}\nabla\varphi}}{{\tau_{h}}}{\rm{+}}% \Delta t{\varpi_{i}}q,italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x + bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_Δ italic_t , italic_t + roman_Δ italic_t ) = italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG [ italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) ] - divide start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG divide start_ARG italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∇ italic_φ end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG + roman_Δ italic_t italic_ϖ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_q , (56)

where hi(𝐱,t)subscript𝑖𝐱𝑡{h_{i}}(\mathbf{x},t)italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) is the distribution function and hieq(𝐱,t)superscriptsubscript𝑖𝑒𝑞𝐱𝑡h_{i}^{eq}(\mathbf{x},t)italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) is the equilibrium distribution function given by

hieq(𝐱,t)={(𝜔01)φ(𝐱,t),i=0,𝜔iφ(𝐱,t),i=1,,8,superscriptsubscript𝑖𝑒𝑞𝐱𝑡casessubscript𝜔01𝜑𝐱𝑡𝑖0missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionsubscript𝜔𝑖𝜑𝐱𝑡𝑖18missing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionmissing-subexpressionh_{i}^{eq}(\mathbf{x},t)=\left\{{\begin{array}[]{*{20}{c}}{(\mathop{\omega}% \nolimits_{0}-1)\varphi(\mathbf{x},t)},&{i=0},\\ {\mathop{\omega}\nolimits_{i}\varphi(\mathbf{x},t)},&{i=1,\cdots,8},\end{array% }}\right.italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) = { start_ARRAY start_ROW start_CELL ( italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT - 1 ) italic_φ ( bold_x , italic_t ) , end_CELL start_CELL italic_i = 0 , end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW start_ROW start_CELL italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_φ ( bold_x , italic_t ) , end_CELL start_CELL italic_i = 1 , ⋯ , 8 , end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL start_CELL end_CELL end_ROW end_ARRAY (57)

ϖisubscriptitalic-ϖ𝑖{\varpi_{i}}italic_ϖ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT represents the diffusional-weights given by ϖ0=0,ϖ18=1/8formulae-sequencesubscriptitalic-ϖ00subscriptitalic-ϖ181/8{\varpi_{0}}=0,{\varpi_{1-8}}={1\mathord{\left/{\vphantom{18}}\right.\kern-1.2% pt}8}italic_ϖ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 , italic_ϖ start_POSTSUBSCRIPT 1 - 8 end_POSTSUBSCRIPT = 1 start_ID / end_ID 8. In addition, the electric potential in the current model can be calculated by

φ=i011𝜔0hi.𝜑subscript𝑖011subscript𝜔0subscript𝑖\varphi=\sum\limits_{i\neq 0}{\frac{1}{{1-\mathop{\omega}\nolimits_{0}}}}{h_{i% }}.italic_φ = ∑ start_POSTSUBSCRIPT italic_i ≠ 0 end_POSTSUBSCRIPT divide start_ARG 1 end_ARG start_ARG 1 - italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT . (58)

Through the Chapman-Enskog analysis shown in the Appendix A, the relaxation time is defined as 𝜏h=0.5+𝜀v/𝑐s2Δtsubscript𝜏0.5subscript𝜀𝑣/superscriptsubscript𝑐𝑠2Δ𝑡\mathop{\tau}\nolimits_{h}=0.5+{{\mathop{\varepsilon}\nolimits_{v}}\mathord{% \left/{\vphantom{{\mathop{\varepsilon}\nolimits_{v}}{\mathop{c}\nolimits_{s}^{% 2}}}}\right.\kern-1.2pt}{\mathop{c}\nolimits_{s}^{2}}}\Delta titalic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT = 0.5 + italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT start_ID / end_ID italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_t. In addition, due to the mesoscopic feature of the LB method, the potential gradient that appeared in the discrete forcing term can be calculated locally by using the following equation

φ=i𝐜iiε^(ϕ)+𝑐s2𝜏hΔt.𝜑subscript𝑖subscript𝐜𝑖subscript𝑖^𝜀italic-ϕsuperscriptsubscript𝑐𝑠2subscript𝜏Δ𝑡\nabla\varphi=-\frac{{\sum\limits_{i}{\mathbf{c}_{i}}\mathop{h}\nolimits_{i}}}% {{\hat{\varepsilon}\left(\phi\right)+\mathop{c}\nolimits_{s}^{2}\mathop{\tau}% \nolimits_{h}\Delta t}}.∇ italic_φ = - divide start_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) + italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t end_ARG . (59)

3.4 Lattice Boltzmann equation for Nernst-Planck equation

Unlike single-phase flows, in a two-phase EHD system, charges are confined around the interface in a thin transition region. Capturing the charge density in such an area with the present diffusive approach requires converting the Nernst-Planck equation into the following equivalent form [65],

qt+𝐮q=(αq)+R.𝑞𝑡𝐮𝑞𝛼𝑞𝑅\frac{{\partial q}}{{\partial t}}+\mathbf{u}\cdot\nabla q=\nabla\cdot(\alpha% \nabla q)+R.divide start_ARG ∂ italic_q end_ARG start_ARG ∂ italic_t end_ARG + bold_u ⋅ ∇ italic_q = ∇ ⋅ ( italic_α ∇ italic_q ) + italic_R . (60)

Here, the forcing term R𝑅Ritalic_R is defined as

R=σqε+σεε𝐄σ𝐄,𝑅𝜎𝑞𝜀𝜎𝜀𝜀𝐄𝜎𝐄R=-\frac{{\sigma q}}{\varepsilon}+\frac{\sigma}{\varepsilon}\nabla\varepsilon% \cdot\mathbf{E}-\nabla\sigma\cdot\mathbf{E},italic_R = - divide start_ARG italic_σ italic_q end_ARG start_ARG italic_ε end_ARG + divide start_ARG italic_σ end_ARG start_ARG italic_ε end_ARG ∇ italic_ε ⋅ bold_E - ∇ italic_σ ⋅ bold_E , (61)

in which (ε𝐄)=q𝜀𝐄𝑞\nabla\cdot(\varepsilon\mathbf{E})=q∇ ⋅ ( italic_ε bold_E ) = italic_q has been used. Inspired by previous works [70, 71], the evolution equation for charge density can thus be expressed as [66]

li(𝐱+𝐜iΔt,t+Δt)=li(𝐱,t)1τl[li(𝐱,t)lieq(𝐱,t)]+Δt𝑆i(𝐱,t)+Δt𝑇i(𝐱,t),subscript𝑙𝑖𝐱subscript𝐜𝑖Δ𝑡𝑡Δ𝑡subscript𝑙𝑖𝐱𝑡1subscript𝜏𝑙delimited-[]subscript𝑙𝑖𝐱𝑡superscriptsubscript𝑙𝑖𝑒𝑞𝐱𝑡Δ𝑡subscript𝑆𝑖𝐱𝑡Δ𝑡subscript𝑇𝑖𝐱𝑡{l_{i}}(\mathbf{x}+{{\bf{c}}_{i}}\Delta t,t+\Delta t)={l_{i}}(\mathbf{x},t)-% \frac{1}{{{\tau_{l}}}}[{l_{i}}(\mathbf{x},t)-l_{i}^{eq}(\mathbf{x},t)]+\Delta t% \mathop{S}\nolimits_{i}(\mathbf{x},t)+\Delta t\mathop{T}\nolimits_{i}(\mathbf{% x},t),italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x + bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT roman_Δ italic_t , italic_t + roman_Δ italic_t ) = italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG [ italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) - italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) ] + roman_Δ italic_t italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) + roman_Δ italic_t italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) , (62)

where Sisubscript𝑆𝑖S_{i}italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and Tisubscript𝑇𝑖T_{i}italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are two discrete forcing terms and they are given by

𝑆i(𝐱,t)=(112τl)𝜔iR,subscript𝑆𝑖𝐱𝑡112subscript𝜏𝑙subscript𝜔𝑖𝑅\mathop{S}\nolimits_{i}(\mathbf{x},t)=(1-\frac{1}{{2{\tau_{l}}}})\mathop{% \omega}\nolimits_{i}R,italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) = ( 1 - divide start_ARG 1 end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG ) italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_R , (63)
𝑇i(𝐱,t)=(112τl)𝜔i𝐜it(q𝐮)𝑐s2.subscript𝑇𝑖𝐱𝑡112subscript𝜏𝑙subscript𝜔𝑖subscript𝐜𝑖subscript𝑡𝑞𝐮superscriptsubscript𝑐𝑠2\mathop{T}\nolimits_{i}(\mathbf{x},t)=(1-\frac{1}{{2{\tau_{l}}}})\frac{{% \mathop{\omega}\nolimits_{i}{{\bf{c}}_{i}}\cdot\mathop{\partial}\nolimits_{t}(% q{\bf{u}})}}{{\mathop{c}\nolimits_{s}^{2}}}.italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( bold_x , italic_t ) = ( 1 - divide start_ARG 1 end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG ) divide start_ARG italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⋅ ∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_q bold_u ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG . (64)

In addition, the local equilibrium function is defined as

lieq(𝐱,t)=𝜔iq(1+𝐜i𝐮𝑐s2).superscriptsubscript𝑙𝑖𝑒𝑞𝐱𝑡subscript𝜔𝑖𝑞1subscript𝐜𝑖𝐮superscriptsubscript𝑐𝑠2l_{i}^{eq}(\mathbf{x},t)=\mathop{\omega}\nolimits_{i}q(1+\frac{{{{\bf{c}}_{i}}% \cdot{\bf{u}}}}{{\mathop{c}\nolimits_{s}^{2}}}).italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ( bold_x , italic_t ) = italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_q ( 1 + divide start_ARG bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⋅ bold_u end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) . (65)

The charge density is determined by

q=ili+12ΔtR,𝑞subscript𝑖subscript𝑙𝑖12Δ𝑡𝑅q=\sum\limits_{i}{{l_{i}}}+\frac{1}{2}\Delta tR,italic_q = ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Δ italic_t italic_R , (66)

and the relaxation time τlsubscript𝜏𝑙\tau_{l}italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT is evaluated by 𝜏l=0.5+α/𝑐s2Δtsubscript𝜏𝑙0.5𝛼/superscriptsubscript𝑐𝑠2Δ𝑡\mathop{\tau}\nolimits_{l}=0.5+{\alpha\mathord{\left/{\vphantom{D{\mathop{c}% \nolimits_{s}^{2}}}}\right.\kern-1.2pt}{\mathop{c}\nolimits_{s}^{2}}}\Delta titalic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT = 0.5 + italic_α start_ID / end_ID italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_t.

4 Wetting boundary condition

When two-phase fluids encounter a solid substrate, the substrate’s wettability significantly influences fluid interface dynamics. Therefore, it is imperative to establish a robust wetting boundary condition that considers the contact angle between the phase interface and the solid surface. This section will present the details of implementing the wetting boundary condition with the current LB method.

As shown in Eq. (20), the wetting boundary condition can be defined once the wall free-energy ψssubscript𝜓𝑠{\psi_{s}}italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT is determined. In this work, we are using the cubic wall free-energy approach [72], which only takes into account the interaction at the three-phase junction and ignores the interactions between the solid and bulk phases, and its express can be defined as

ψs=b12ϕ2b13ϕ3,subscript𝜓𝑠subscript𝑏12superscriptitalic-ϕ2subscript𝑏13superscriptitalic-ϕ3{\psi_{s}}=\frac{{{b_{1}}}}{2}{\phi^{2}}-\frac{{{b_{1}}}}{3}{\phi^{3}},italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT = divide start_ARG italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_ϕ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - divide start_ARG italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 3 end_ARG italic_ϕ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT , (67)

where b1subscript𝑏1b_{1}italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT is a model parameter. Additionally, following the route provided in Ref.[73], it is noted that ψssubscript𝜓𝑠{\psi_{s}}italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT in the bulk region must meet with

dψsdϕ=±2κψ0.𝑑subscript𝜓𝑠𝑑italic-ϕplus-or-minus2𝜅subscript𝜓0\frac{{d{\psi_{s}}}}{{d\phi}}=\pm\sqrt{2\kappa{\psi_{0}}}.divide start_ARG italic_d italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_ARG start_ARG italic_d italic_ϕ end_ARG = ± square-root start_ARG 2 italic_κ italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG . (68)

Combing Eq. (67) and Eq. (68), one can find that ψssubscript𝜓𝑠{\psi_{s}}italic_ψ start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT shares two stable solutions, i.e., ψs1=0subscript𝜓𝑠10{\psi_{s1}}=0italic_ψ start_POSTSUBSCRIPT italic_s 1 end_POSTSUBSCRIPT = 0 and ψs2=1subscript𝜓𝑠21{\psi_{s2}}=1italic_ψ start_POSTSUBSCRIPT italic_s 2 end_POSTSUBSCRIPT = 1. Thus, the surface tensions for gas-solid and liquid-solid can be given by

γsg=b12ψs12b13ψs13+0ψs12κψ0𝑑ϕ=0,subscript𝛾𝑠𝑔subscript𝑏12superscriptsubscript𝜓𝑠12subscript𝑏13superscriptsubscript𝜓𝑠13superscriptsubscript0subscript𝜓𝑠12𝜅subscript𝜓0differential-ditalic-ϕ0{\gamma_{sg}}=\frac{{{b_{1}}}}{2}\psi_{s1}^{2}-\frac{{{b_{1}}}}{3}\psi_{s1}^{3% }+\int_{0}^{{\psi_{s1}}}{\sqrt{2\kappa{\psi_{0}}}}d\phi=0,italic_γ start_POSTSUBSCRIPT italic_s italic_g end_POSTSUBSCRIPT = divide start_ARG italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_ψ start_POSTSUBSCRIPT italic_s 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - divide start_ARG italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 3 end_ARG italic_ψ start_POSTSUBSCRIPT italic_s 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ψ start_POSTSUBSCRIPT italic_s 1 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT square-root start_ARG 2 italic_κ italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG italic_d italic_ϕ = 0 , (69)
γsl=b16ψs22b13ψs23+1ψs22κψ0𝑑ϕ=b16.subscript𝛾𝑠𝑙subscript𝑏16superscriptsubscript𝜓𝑠22subscript𝑏13superscriptsubscript𝜓𝑠23superscriptsubscript1subscript𝜓𝑠22𝜅subscript𝜓0differential-ditalic-ϕsubscript𝑏16{\gamma_{sl}}=\frac{{{b_{1}}}}{6}\psi_{s2}^{2}-\frac{{{b_{1}}}}{3}\psi_{s2}^{3% }+\int_{1}^{{\psi_{s2}}}{\sqrt{2\kappa{\psi_{0}}}}d\phi=\frac{{{b_{1}}}}{6}.italic_γ start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT = divide start_ARG italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 6 end_ARG italic_ψ start_POSTSUBSCRIPT italic_s 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - divide start_ARG italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 3 end_ARG italic_ψ start_POSTSUBSCRIPT italic_s 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + ∫ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_ψ start_POSTSUBSCRIPT italic_s 2 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT square-root start_ARG 2 italic_κ italic_ψ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG italic_d italic_ϕ = divide start_ARG italic_b start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 6 end_ARG . (70)

It is known that the local static contact angle θYsubscript𝜃𝑌{\theta_{Y}}italic_θ start_POSTSUBSCRIPT italic_Y end_POSTSUBSCRIPT on a chemically homogeneous wall satisfies the Young’s equation

cosθY=γsgγslγ.subscript𝜃𝑌subscript𝛾𝑠𝑔subscript𝛾𝑠𝑙𝛾\cos{\theta_{Y}}=\frac{{{\gamma_{sg}}-{\gamma_{sl}}}}{\gamma}.roman_cos italic_θ start_POSTSUBSCRIPT italic_Y end_POSTSUBSCRIPT = divide start_ARG italic_γ start_POSTSUBSCRIPT italic_s italic_g end_POSTSUBSCRIPT - italic_γ start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT end_ARG start_ARG italic_γ end_ARG . (71)

Based on Eqs. (20), (69) and (70), we can ultimately obtain the wetting boundary condition as

𝐧wϕ=2βκcosθY(ϕϕ2).subscript𝐧𝑤italic-ϕ2𝛽𝜅subscript𝜃𝑌italic-ϕsuperscriptitalic-ϕ2\mathbf{n}_{w}\cdot\nabla\phi=-\sqrt{\frac{{2\beta}}{\kappa}}\cos{\theta_{Y}}(% \phi-\mathop{\phi}\nolimits^{2}).bold_n start_POSTSUBSCRIPT italic_w end_POSTSUBSCRIPT ⋅ ∇ italic_ϕ = - square-root start_ARG divide start_ARG 2 italic_β end_ARG start_ARG italic_κ end_ARG end_ARG roman_cos italic_θ start_POSTSUBSCRIPT italic_Y end_POSTSUBSCRIPT ( italic_ϕ - italic_ϕ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) . (72)

It is worth noting that the wetting boundary condition mentioned above is established without considering the influence of the electric field. Since Lippmann’s pioneering work on electro-capillarity [29], it has been discovered that the divergent electric force leads to a significant deformation of the fluid interface and causes a large curvature in the vicinity of the contact line. In such situations, the apparent contact angle θBsubscript𝜃𝐵{\theta_{B}}italic_θ start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT can be described by the Lippmann equation [29],

cosθB=cosθY+εsφ22γd,subscript𝜃𝐵subscript𝜃𝑌subscript𝜀𝑠superscript𝜑22𝛾𝑑\cos{\theta_{B}}=\cos{\theta_{Y}}+\frac{{{\varepsilon_{s}}{\varphi^{2}}}}{{2% \gamma d}},roman_cos italic_θ start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT = roman_cos italic_θ start_POSTSUBSCRIPT italic_Y end_POSTSUBSCRIPT + divide start_ARG italic_ε start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_φ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_γ italic_d end_ARG , (73)

where φ𝜑\varphiitalic_φ is the applied voltage, εssubscript𝜀𝑠\varepsilon_{s}italic_ε start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT and d𝑑ditalic_d are the permittivity and thickness of the dielectric substrate, respectively. In the literature, the second term on the right-hand side of Eq. (73) is also referred to as the electrowetting number. This number is used to measure the relative strength of the electric force in contrast to the surface tension at the fluid surface, and it can be rewritten as

εsφ2/[2γd]=η/d,subscript𝜀𝑠superscript𝜑2/delimited-[]2𝛾𝑑𝜂/𝑑{{{\varepsilon_{s}}{\varphi^{2}}}\mathord{\left/{\vphantom{{{\varepsilon_{s}}{% \varphi^{2}}}{\left[{2\gamma d}\right]}}}\right.\kern-1.2pt}{\left[{2\gamma d}% \right]}}={\eta\mathord{\left/{\vphantom{\eta d}}\right.\kern-1.2pt}d},italic_ε start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_φ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_ID / end_ID [ 2 italic_γ italic_d ] = italic_η start_ID / end_ID italic_d , (74)

in which η=0.5εsφ2/γ𝜂0.5subscript𝜀𝑠superscript𝜑2/𝛾\eta={{0.5{\varepsilon_{s}}{\varphi^{2}}}\mathord{\left/{\vphantom{{0.5{% \varepsilon_{s}}{\varphi^{2}}}\gamma}}\right.\kern-1.2pt}\gamma}italic_η = 0.5 italic_ε start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_φ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_ID / end_ID italic_γ. Comparing Eq. (71) and Eq. (73), it is found that the electrowetting boundary condition can be established by using an effective surface tension for liquid-solid, which is given by

γ^sl=γslεsφ22γd.subscript^𝛾𝑠𝑙subscript𝛾𝑠𝑙subscript𝜀𝑠superscript𝜑22𝛾𝑑{{\hat{\gamma}}_{sl}}={\gamma_{sl}}-\frac{{{\varepsilon_{s}}{\varphi^{2}}}}{{2% \gamma d}}.over^ start_ARG italic_γ end_ARG start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT = italic_γ start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT - divide start_ARG italic_ε start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_φ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_γ italic_d end_ARG . (75)

In summary, the electrowetting boundary condition is still implemented with Eq. (72) expect that the previous surface tension for liquid-solid γslsubscript𝛾𝑠𝑙{\gamma_{sl}}italic_γ start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT is replaced by γ^slsubscript^𝛾𝑠𝑙{{\hat{\gamma}}_{sl}}over^ start_ARG italic_γ end_ARG start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT.

5 Implementation scheme of the proposed LB method

Due to the nonlinear coupling characteristic of governing equations, we use an iterative scheme outlined in Algorithm 5 to implement the current LB method. The key algorithm of this approach involves a cyclic sequence of substeps, where each cycle corresponds to one time step. In each iteration, the macroscopic variables in the system are updated by solving the LB equation for every physical field. It is worth noting that because the electric potential equation is governed by a steady elliptic equation, an inner loop is required to obtain the convergent solution for the electric potential at each time step.

  Algorithm 1 Implementation scheme of the proposed LB method

 

  ##\## 1. At the beginning (n=0𝑛0n=0italic_n = 0), the fluid variables (ϕitalic-ϕ\phiitalic_ϕ, ρ𝜌\rhoitalic_ρ, 𝐮𝐮{\bf{u}}bold_u, φ𝜑\varphiitalic_φ, q𝑞qitalic_q) and the distribution functions (gi0,hi0,li0,fi0superscriptsubscript𝑔𝑖0superscriptsubscript𝑖0superscriptsubscript𝑙𝑖0superscriptsubscript𝑓𝑖0g_{i}^{0},h_{i}^{0},l_{i}^{0},f_{i}^{0}italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT , italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT , italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT , italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT) at each grid point 𝐱isubscript𝐱𝑖\mathbf{x}_{i}bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are initialized
  ##\## 2. Solve the Chan-Hilliard equation Eq. (31) using Eq. (45)
  for all 𝐱isubscript𝐱𝑖\mathbf{x}_{i}bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT do
     𝑔in+1𝑔in1𝜏g[𝑔in𝑔ieq,n]+Δt𝐺in+12Δt2t𝐺insuperscriptsubscript𝑔𝑖n1superscriptsubscript𝑔𝑖𝑛1subscript𝜏𝑔delimited-[]superscriptsubscript𝑔𝑖𝑛superscriptsubscript𝑔𝑖𝑒𝑞𝑛Δ𝑡superscriptsubscript𝐺𝑖𝑛12superscriptΔ𝑡2subscript𝑡superscriptsubscript𝐺𝑖𝑛\mathop{g}\nolimits_{i}^{{\rm{n+1}}}\leftarrow\mathop{g}\nolimits_{i}^{n}-% \displaystyle\frac{1}{{\mathop{\tau}\nolimits_{g}}}[\mathop{g}\nolimits_{i}^{n% }-\mathop{g}\nolimits_{i}^{eq,\;n}]+\Delta t\mathop{G}\nolimits_{i}^{n}+\frac{% 1}{2}\mathop{\Delta t}\nolimits^{2}\mathop{\partial}\nolimits_{t}\mathop{G}% \nolimits_{i}^{n}italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_n + 1 end_POSTSUPERSCRIPT ← italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_g end_POSTSUBSCRIPT end_ARG [ italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT ] + roman_Δ italic_t italic_G start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG start_BIGOP roman_Δ italic_t end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_G start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT,
     ##\## where gieq,nsuperscriptsubscript𝑔𝑖𝑒𝑞𝑛g_{i}^{eq,\;n}italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT and Ginsuperscriptsubscript𝐺𝑖𝑛G_{i}^{n}italic_G start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT are calculated with Eq. (46) and Eq. (47),
     ϕn+1igin+1superscriptitalic-ϕ𝑛1subscript𝑖superscriptsubscript𝑔𝑖𝑛1\phi^{n+1}\leftarrow\sum\limits_{i}{g_{i}^{n+1}}italic_ϕ start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT ← ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_g start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT,
     ρn+1(ϕ)ϕn+1(ρlρv)+ρvsuperscript𝜌𝑛1italic-ϕsuperscriptitalic-ϕ𝑛1subscript𝜌𝑙subscript𝜌𝑣subscript𝜌𝑣\rho^{n+1}\left(\phi\right)\leftarrow\phi^{n+1}\left({{\rho_{l}}-{\rho_{v}}}% \right)+{\rho_{v}}italic_ρ start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT ( italic_ϕ ) ← italic_ϕ start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT ( italic_ρ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ) + italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT.
  end for##\## 3. Solve the electric potential equation Eq. (55) using Eq. (56)
  for all 𝐱isubscript𝐱𝑖\mathbf{x}_{i}bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT do
     in+1in1𝜏h[inieq,n]ε^n(ϕ)𝑐s2𝜔i𝐜iϕn𝜏h+Δtϖiqnsuperscriptsubscript𝑖n1superscriptsubscript𝑖𝑛1subscript𝜏delimited-[]superscriptsubscript𝑖𝑛superscriptsubscript𝑖𝑒𝑞𝑛superscript^𝜀𝑛italic-ϕsuperscriptsubscript𝑐𝑠2subscript𝜔𝑖subscript𝐜𝑖superscriptitalic-ϕ𝑛subscript𝜏Δ𝑡subscriptitalic-ϖ𝑖superscript𝑞𝑛\mathop{h}\nolimits_{i}^{{\rm{n+1}}}\leftarrow\mathop{h}\nolimits_{i}^{n}-% \displaystyle\frac{1}{{\mathop{\tau}\nolimits_{h}}}[\mathop{h}\nolimits_{i}^{n% }-\mathop{h}\nolimits_{i}^{eq,\;n}]-\frac{{\hat{\varepsilon}^{n}(\phi)}}{{% \mathop{c}\nolimits_{s}^{2}}}\frac{{\mathop{\omega}\nolimits_{i}\mathbf{c}_{i}% \nabla\phi^{n}}}{{\mathop{\tau}\nolimits_{h}}}+\Delta t\mathop{\varpi}% \nolimits_{i}q^{n}italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_n + 1 end_POSTSUPERSCRIPT ← italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG [ italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT ] - divide start_ARG over^ start_ARG italic_ε end_ARG start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ( italic_ϕ ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG divide start_ARG italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∇ italic_ϕ start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG + roman_Δ italic_t italic_ϖ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_q start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT,
     ##\## where hieq,nsuperscriptsubscript𝑖𝑒𝑞𝑛h_{i}^{eq,\;n}italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT and φ𝜑\varphiitalic_φ are calculated utilizing Eq. (57) and Eq. (58),
     φn+1i011𝜔0hin+1superscript𝜑𝑛1subscript𝑖011subscript𝜔0superscriptsubscript𝑖𝑛1\varphi^{n+1}\leftarrow\sum\limits_{i\neq 0}{\displaystyle\frac{1}{{1-\mathop{% \omega}\nolimits_{0}}}}{h_{i}^{n+1}}italic_φ start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT ← ∑ start_POSTSUBSCRIPT italic_i ≠ 0 end_POSTSUBSCRIPT divide start_ARG 1 end_ARG start_ARG 1 - italic_ω start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT,
     repeat until i[𝜑in+1𝜑in]2i[𝜑in]2<105subscript𝑖superscriptdelimited-[]superscriptsubscript𝜑𝑖𝑛1superscriptsubscript𝜑𝑖𝑛2subscript𝑖superscriptdelimited-[]superscriptsubscript𝜑𝑖𝑛2superscript105\sqrt{\displaystyle\frac{{\sum\limits_{i}{\mathop{[\mathop{\varphi}\nolimits_{% i}^{n+1}-\mathop{\varphi}\nolimits_{i}^{n}]}\nolimits^{2}}}}{{\sum\limits_{i}{% \mathop{[\mathop{\varphi}\nolimits_{i}^{n}]}\nolimits^{2}}}}}<\mathop{10}% \nolimits^{-5}square-root start_ARG divide start_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_BIGOP [ italic_φ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT - italic_φ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ] end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_BIGOP [ italic_φ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ] end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG < 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT.
  end for##\## 4. Solve the Nernst-Planck equation Eq. (60) using Eq. (62
  for all 𝐱isubscript𝐱𝑖\mathbf{x}_{i}bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT do
     𝑙in+1𝑙in1𝜏l[𝑙in𝑙ieq,n]+Δt𝑆in+Δt𝑇insuperscriptsubscript𝑙𝑖n1superscriptsubscript𝑙𝑖𝑛1subscript𝜏𝑙delimited-[]superscriptsubscript𝑙𝑖𝑛superscriptsubscript𝑙𝑖𝑒𝑞𝑛Δ𝑡superscriptsubscript𝑆𝑖𝑛Δ𝑡superscriptsubscript𝑇𝑖𝑛\mathop{l}\nolimits_{i}^{{\rm{n+1}}}\leftarrow\mathop{l}\nolimits_{i}^{n}-% \displaystyle\frac{1}{{\mathop{\tau}\nolimits_{l}}}[\mathop{l}\nolimits_{i}^{n% }-\mathop{l}\nolimits_{i}^{eq,\;n}]+\Delta t\mathop{S}\nolimits_{i}^{n}+\Delta t% \mathop{T}\nolimits_{i}^{n}italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_n + 1 end_POSTSUPERSCRIPT ← italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG [ italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT ] + roman_Δ italic_t italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT + roman_Δ italic_t italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT,
     ##\## where Sinsuperscriptsubscript𝑆𝑖𝑛S_{i}^{n}italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT, Tinsuperscriptsubscript𝑇𝑖𝑛T_{i}^{n}italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT and lieq,nsuperscriptsubscript𝑙𝑖𝑒𝑞𝑛l_{i}^{eq,\;n}italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT are calculated using Eqs. (63), (64) and (65), respectively,
     qn+1ilin+1+12ΔtRnsuperscript𝑞𝑛1subscript𝑖superscriptsubscript𝑙𝑖𝑛112Δ𝑡superscript𝑅𝑛q^{n+1}\leftarrow\sum\limits_{i}{{l_{i}^{n+1}}}+\displaystyle\frac{1}{2}\Delta tR% ^{n}italic_q start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT ← ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Δ italic_t italic_R start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT.
  end for##\## 5. Solve the hydrodynamic equations Eq. (29) and Eq. (30) using Eq. (43
  for all 𝐱isubscript𝐱𝑖\mathbf{x}_{i}bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and t𝑡titalic_t do
     𝑓in+1𝑓in1𝜏f[𝑓in𝑓ieq,n]+Δt(112𝜏f)R¯insuperscriptsubscript𝑓𝑖n1superscriptsubscript𝑓𝑖𝑛1subscript𝜏𝑓delimited-[]superscriptsubscript𝑓𝑖𝑛superscriptsubscript𝑓𝑖𝑒𝑞𝑛Δ𝑡112subscript𝜏𝑓superscriptsubscript¯𝑅𝑖𝑛\mathop{f}\nolimits_{i}^{{\rm{n+1}}}\leftarrow\mathop{f}\nolimits_{i}^{n}-% \displaystyle\frac{1}{{\mathop{\tau}\nolimits_{f}}}[\mathop{f}\nolimits_{i}^{n% }-\mathop{f}\nolimits_{i}^{eq,\;n}]+\Delta t(1-\frac{1}{{2\mathop{\tau}% \nolimits_{f}}})\mathop{\bar{R}}\nolimits_{i}^{n}italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_n + 1 end_POSTSUPERSCRIPT ← italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT end_ARG [ italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT - italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT ] + roman_Δ italic_t ( 1 - divide start_ARG 1 end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT end_ARG ) start_BIGOP over¯ start_ARG italic_R end_ARG end_BIGOP start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT,
     ##\## where R¯insuperscriptsubscript¯𝑅𝑖𝑛{\bar{R}}_{i}^{n}over¯ start_ARG italic_R end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT and fieq,nsuperscriptsubscript𝑓𝑖𝑒𝑞𝑛f_{i}^{eq,\;n}italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q , italic_n end_POSTSUPERSCRIPT are calculated using Eq. (39) and Eq. (40),
     ρ𝐮n+1i𝐜ifin+1+12Δt𝐅n𝜌superscript𝐮𝑛1subscript𝑖subscript𝐜𝑖superscriptsubscript𝑓𝑖𝑛112Δ𝑡superscript𝐅𝑛\rho\mathbf{u}^{n+1}\leftarrow\sum\limits_{i}{\mathbf{c}_{i}{f_{i}^{n+1}}}+% \displaystyle\frac{1}{2}\Delta t\mathbf{F}^{n}italic_ρ bold_u start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT ← ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n + 1 end_POSTSUPERSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG roman_Δ italic_t bold_F start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT.
  end for##\## 6. Advance the time step (n+1n𝑛1𝑛n+1\leftarrow nitalic_n + 1 ← italic_n) and repeat step 2 to 5 until the end time or steady state is reached

 

6 Numerical Validation

In this section, we will validate the proposed thermodynamically consistent phase-field LB method by simulating the EHD flows with two superimposed planar fluids [74], the charge relaxation of a Gaussian bell [74] and the equilibrium interface profile of electrowetting on dielectric. Our primary focus will be on its performance in predicting the electric potential, charge density, and contact angle of a charged droplet.

6.1 EHD flows with two superimposed planar fluids

Refer to caption
Figure 1: Schematic diagram of EHD flows with two superimposed planar fluids.

We first validate the proposed model by simulating the EHD flows with two superimposed planar fluids [74]. The configuration is shown in Fig. 1, in which a square enclosure with width being L𝐿Litalic_L is filled with two conducting fluids, and the height of these two fluids are both set to be L/2𝐿2L/2italic_L / 2. Note that the electric properties of these two fluids are homogeneous, and the potential at the upper wall and bottom wall are given by φLsubscript𝜑𝐿{\varphi_{L}}italic_φ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT and φHsubscript𝜑𝐻{\varphi_{H}}italic_φ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT (φH>φLsubscript𝜑𝐻subscript𝜑𝐿{\varphi_{H}}>{\varphi_{L}}italic_φ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT > italic_φ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT), respectively. According to the work of López-Herrera et al. [74], the analytical solution for the electric potential in each medium can be given by

𝜑1exact(y)=2y+R^1+R^,𝜑2exact(y)=2y+11+R^R^,formulae-sequencesuperscriptsubscript𝜑1𝑒𝑥𝑎𝑐𝑡𝑦2𝑦^𝑅1^𝑅superscriptsubscript𝜑2𝑒𝑥𝑎𝑐𝑡𝑦2𝑦11^𝑅^𝑅\mathop{\varphi}\nolimits_{1}^{exact}(y)=\frac{{-2y+\hat{R}}}{{1+\hat{R}}},\;% \;\;\;\mathop{\varphi}\nolimits_{2}^{exact}(y)=\frac{{-2y+1}}{{1+\hat{R}}}\hat% {R},italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_x italic_a italic_c italic_t end_POSTSUPERSCRIPT ( italic_y ) = divide start_ARG - 2 italic_y + over^ start_ARG italic_R end_ARG end_ARG start_ARG 1 + over^ start_ARG italic_R end_ARG end_ARG , italic_φ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_x italic_a italic_c italic_t end_POSTSUPERSCRIPT ( italic_y ) = divide start_ARG - 2 italic_y + 1 end_ARG start_ARG 1 + over^ start_ARG italic_R end_ARG end_ARG over^ start_ARG italic_R end_ARG , (76)

where R^^𝑅\hat{R}over^ start_ARG italic_R end_ARG denotes the ratio of conductivities between two fluids i.e., R^=σl/σv^𝑅subscript𝜎𝑙/subscript𝜎𝑣{{\hat{R}={\sigma_{l}}}\mathord{\left/{\vphantom{{R={\sigma_{l}}}{{\sigma_{v}}% }}}\right.\kern-1.2pt}{{\sigma_{v}}}}over^ start_ARG italic_R end_ARG = italic_σ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_σ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT. Not that the electric properties used in our simulations are the same as those adopted by López-Herrera et al. [74]. Fig. 2 gives the profiles of φ𝜑\varphiitalic_φ along the y𝑦yitalic_y coordinate, in which the mesh is set to be 100×100100100100\times 100100 × 100. It is shown that our numerical data give reasonably accurate results. To evaluate the convergence rate of the proposed LB model, the relative errors under different grid resolutions are further calculated, which is defined as

Eφ=i(𝜑numerical(𝐱i)𝜑analytical(𝐱i))2i𝜑analytical(𝐱i)2,subscript𝐸𝜑subscript𝑖superscriptsubscript𝜑𝑛𝑢𝑚𝑒𝑟𝑖𝑐𝑎𝑙subscript𝐱𝑖subscript𝜑𝑎𝑛𝑎𝑙𝑦𝑡𝑖𝑐𝑎𝑙subscript𝐱𝑖2subscript𝑖superscriptsubscript𝜑𝑎𝑛𝑎𝑙𝑦𝑡𝑖𝑐𝑎𝑙subscript𝐱𝑖2E_{\varphi}=\sqrt{\frac{{\sum\limits_{i}{\mathop{(\mathop{\varphi}\nolimits_{% numerical}(\mathbf{x}_{i})-\mathop{\varphi}\nolimits_{analytical}(\mathbf{x}_{% i}))}\nolimits^{2}}}}{{\sum\limits_{i}{\mathop{\mathop{\varphi}\nolimits_{% analytical}(\mathbf{x}_{i})}\nolimits^{2}}}}},italic_E start_POSTSUBSCRIPT italic_φ end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_BIGOP ( italic_φ start_POSTSUBSCRIPT italic_n italic_u italic_m italic_e italic_r italic_i italic_c italic_a italic_l end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) - italic_φ start_POSTSUBSCRIPT italic_a italic_n italic_a italic_l italic_y italic_t italic_i italic_c italic_a italic_l end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) ) end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_BIGOP italic_φ start_POSTSUBSCRIPT italic_a italic_n italic_a italic_l italic_y italic_t italic_i italic_c italic_a italic_l end_POSTSUBSCRIPT ( bold_x start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG , (77)

where the subscripts ”analytical” and ”numerical” represent the analytical solution and numerical result, respectively, and the summation is over the entire domain. As seen from Fig. 2, the current LB scheme shares a second order accuracy in space. The above numerical results demonstrate the feasibility of the present treatment for electrostatic field, and verifies that the electric potential equation can be solved properly by the current proposed model.

Refer to caption
Refer to caption
Figure 2: (a) Comparison of the potential between analytical solutions given by Eq. (76) and numerical results for the EHD flows with two superimposed planar fluids; (b) Relative error of the electric potential versus the mesh size, in which the ratios of the droplet’s conductivity and permittivity to surrounding fluids are set to be 3.0 and 2.0, respectively.

6.2 Charge diffusion of a Gaussian bell

Refer to caption
Refer to caption
Figure 3: (a) Comparison of the charge density at diffrtent time between the analytical solution given by Eq. (79) and the numerical results for charge diffusion of a Gaussian bell; (b) relative error of the charge density versus the mesh size, in which the simulated parameters are chosen as σ=1.0𝜎1.0\sigma=1.0italic_σ = 1.0, ε=2.0𝜀2.0\varepsilon=2.0italic_ε = 2.0, a=0.05𝑎0.05a=0.05italic_a = 0.05 and L=1.0𝐿1.0L=1.0italic_L = 1.0.

In order to further verify the proposed model, another test problem discussed is the time diffusion of a charge density distribution [74]. In this problem, a square enclosure with width of L𝐿Litalic_L is occupied by a single-phase fluid, and the electric potential at the four boundaries is set to be 0.0. The initial shape of the charge density obeys a Gaussian bell defined as

q(x,y,t=0)=exp(ι22a2)a2π,𝑞𝑥𝑦𝑡0superscript𝜄22superscript𝑎2𝑎2𝜋q\left({x,y,t=0}\right)=\frac{{\exp\left({\frac{{-{\iota^{2}}}}{{2{a^{2}}}}}% \right)}}{{a\sqrt{2\pi}}},italic_q ( italic_x , italic_y , italic_t = 0 ) = divide start_ARG roman_exp ( divide start_ARG - italic_ι start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_a start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) end_ARG start_ARG italic_a square-root start_ARG 2 italic_π end_ARG end_ARG , (78)

where ι2=x2+y2superscript𝜄2superscript𝑥2superscript𝑦2{\iota^{2}}={x^{2}}+{y^{2}}italic_ι start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT = italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_y start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT and a𝑎aitalic_a is a free parameter determining the width and height of the bell. If the domain boundaries are sufficiently far from the concentrated charge bump, i.e., aLmuch-less-than𝑎𝐿a\ll Litalic_a ≪ italic_L, the analytical solution for the charge density can be expressed as

q(x,y,t)=exp(ι22a2+σtε)a2π.𝑞𝑥𝑦𝑡superscript𝜄22superscript𝑎2𝜎𝑡𝜀𝑎2𝜋q\left({x,y,t}\right)=\frac{{\exp\left({\frac{{-{\iota^{2}}}}{{2{a^{2}}}}+% \frac{{\sigma t}}{\varepsilon}}\right)}}{{a\sqrt{2\pi}}}.italic_q ( italic_x , italic_y , italic_t ) = divide start_ARG roman_exp ( divide start_ARG - italic_ι start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 italic_a start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG + divide start_ARG italic_σ italic_t end_ARG start_ARG italic_ε end_ARG ) end_ARG start_ARG italic_a square-root start_ARG 2 italic_π end_ARG end_ARG . (79)

We now perform some simulations under different times with the conductivity, permittivity and free parameter being σ=1.0𝜎1.0\sigma=1.0italic_σ = 1.0, ε=2.0𝜀2.0\varepsilon=2.0italic_ε = 2.0, and a=0.05𝑎0.05a=0.05italic_a = 0.05, respectively. In addition, the width of the square enclosure is set to be L=1.0𝐿1.0L=1.0italic_L = 1.0, which is large enough to derive the analytical solution. Fig. 3 presents the comparison between the numerical results and analytical solution for the charge density along the horizontal center line, where the lattice size is 200×200200200200\times 200200 × 200. Good consistence can be observed, which verifies the applicability of the current model for Nernst-Planck equation. Further, the convergence rate of the present model is investigated as well. To this end, the relative error defined by Eq. (76) is calculated with grid resolution varying from 50 to 400 at t=2𝑡2t=2italic_t = 2. As shown in Fig. 3, one can clearly seen that our model shares a second-order convergence rate, which is consistent with the theoretic analysis.

6.3 Equilibrium interface profiles

The above two simulations are conducted without considering the influence of the wetting behaviour. In this example, we intend to test the model’s capability in predicting the equilibrium interface profiles and compare them with the analytical solutions predicted by the Lippmann equation. In the simulation, the lattice size of the computation domain is set to be 200×100200100200\times 100200 × 100, which is fine enough to give the grid independence results. Initially, a conductive droplet with a radius of r=25𝑟25r=25italic_r = 25 is located on the bottom wall, and the area surrounding the droplet is filled with an insulating gas. The order profile function is defined as

ϕ=0.5+0.5tanh{2.0×[r(x200)2+(y1)2]W},italic-ϕ0.50.5tanh2.0delimited-[]rsuperscript𝑥2002superscript𝑦12𝑊\phi{\rm{=0}}{\rm{.5+0}}{\rm{.5tanh\{}}\frac{{{\rm{2}}{\rm{.0}}\times[{\rm{r-}% }\sqrt{\mathop{(x-200)}\nolimits^{2}+\mathop{(y-1)}\nolimits^{2}}]}}{W}\},italic_ϕ = 0.5 + 0.5 roman_tanh { divide start_ARG 2.0 × [ roman_r - square-root start_ARG start_BIGOP ( italic_x - 200 ) end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + start_BIGOP ( italic_y - 1 ) end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ] end_ARG start_ARG italic_W end_ARG } , (80)

in which the interface thickness W𝑊Witalic_W is fixed at four lattice units in our simulations. The density ratio (ρl/ρvsubscript𝜌𝑙/subscript𝜌𝑣{{{\rho_{l}}}\mathord{\left/{\vphantom{{{\rho_{l}}}{{\rho_{v}}}}}\right.\kern-% 1.2pt}{{\rho_{v}}}}italic_ρ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT), the kinetic viscosity ratio (μl/μvsubscript𝜇𝑙/subscript𝜇𝑣{{{\mu_{l}}}\mathord{\left/{\vphantom{{{\mu_{l}}}{{\mu_{v}}}}}\right.\kern-1.2% pt}{{\mu_{v}}}}italic_μ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT), and the permittivity ratio are set to be 100, 100, and 81, respectively. The top wall is kept at a low electric potential φLsubscript𝜑𝐿\varphi_{L}italic_φ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT, and the bottom wall is held at a high electric potential φHsubscript𝜑𝐻\varphi_{H}italic_φ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT. In addition, apart from the periodic boundary condition used in the horizontal direction and the no-slip bounce-back boundary condition imposed at the rigid walls [44], the wetting boundary condition described by Eq. (72) should also be adopted for the bottom wall. Moreover, to examine the effect of an electric field on wetting behaviour, it is essential to form a steady droplet pattern before applying the electric field. Thus, in our simulation, we first run the code without an electric field until it reaches the prescribed Young’s contact angle. Following this, the droplet will be relaxed by incorporating the effect of an electric field and ultimately stabilize into an equilibrium pattern.

Refer to caption
Refer to caption
Figure 4: (a) Equilibrium profile of the droplet interface for different values of the dimensionless parameters η𝜂\etaitalic_η and d𝑑ditalic_d; (b) Variation of cosθB𝑐𝑜𝑠subscript𝜃𝐵cos{\theta}_{B}italic_c italic_o italic_s italic_θ start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT versus η/d𝜂𝑑\eta/ditalic_η / italic_d and the comparison with predictions of the Young-Lippmann equation.

In Fig. 4, the graph illustrates the change in the equilibrium interface profile of the droplet for different η𝜂\etaitalic_η and d𝑑ditalic_d with εv=1.0subscript𝜀𝑣1.0{{\varepsilon_{v}}}=1.0italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 1.0. It is evident that when a voltage is applied between the conductive droplet and the solid surface coated with a dielectric material, the change of interfacial energy leads to a decrease in the apparent contact angle compared to the case in the absence of the electric field. To quantify the numerical results, we present the variation of cosθBsubscript𝜃𝐵\cos{\theta_{B}}roman_cos italic_θ start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT for different d𝑑ditalic_d and εvsubscript𝜀𝑣{{\varepsilon_{v}}}italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT. It is observed that cosθBsubscript𝜃𝐵\cos{\theta_{B}}roman_cos italic_θ start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT increases linearly with a slope close to 1, and the permittivity of the fluid outside the droplet has an insignificant influence on the interface profile. The current results align closely with the analytical solutions provided by the Lippmann equation (Fig. 4).

7 Applications

In this section, the proposed method is used to simulate two practical problems: droplet deformation in an electric field and droplet detachment in reversed electrowetting. The first problem demonstrates our model’s ability to capture the fluid-fluid interface dynamics under an electric field. On the other hand, the second problem shows that the proposed method can replicate the motion of the contact line of a charged droplet. The simulated results will be compared against theoretical, experimental, and previous numerical data to validate their accuracy.

7.1 Droplet deformation in an electric field

Refer to caption
Figure 5: The configuration of a droplet suspended in another leaky dielectric fluid under a uniform electric field. The subscripts l𝑙litalic_l and v𝑣vitalic_v represents the internal and external fluids, respectively.

The deformation of a droplet in the presence of an externally applied electric field is a critical scientific concern in the EHD community. It has been found that when a uniform electric field is applied, the droplet can deform either a prolate shape (aligned with the electric field) or an oblate shape (perpendicular to the electric field), which is a result of the combined effects of conductivity and permittivity ratios. To show that the current model is also capable of handling such a complex problem, some simulations are carried out and our numerical results are compared with the theoretical solutions and the available numerical data.

The configuration is depicted in Fig. 5. Initially, the droplet suspended in another leaky dielectric fluid is a circle droplet with a radius of r𝑟ritalic_r, and it is placed in the middle of a square cavity. To induce the droplet deformation, a uniform electric field 𝐄𝐄\bf{E}bold_E is applied in the vertical direction, and the electric potentials at the top wall and the bottom wall are set to be φHsubscript𝜑𝐻{\varphi_{H}}italic_φ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT and φLsubscript𝜑𝐿{\varphi_{L}}italic_φ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT (φH>φLsubscript𝜑𝐻subscript𝜑𝐿{\varphi_{H}}>{\varphi_{L}}italic_φ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT > italic_φ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT), respectively. In this setting, the droplet deformation in the presence of an electric field can be characterized utilizing a parameter D^^𝐷\hat{D}over^ start_ARG italic_D end_ARG given by [13]

D^=L^H^L^+H^,^𝐷^𝐿^𝐻^𝐿^𝐻\hat{D}=\frac{{\hat{L}-\hat{H}}}{{\hat{L}+\hat{H}}},over^ start_ARG italic_D end_ARG = divide start_ARG over^ start_ARG italic_L end_ARG - over^ start_ARG italic_H end_ARG end_ARG start_ARG over^ start_ARG italic_L end_ARG + over^ start_ARG italic_H end_ARG end_ARG , (81)

in which L^^𝐿\hat{L}over^ start_ARG italic_L end_ARG denotes the length of the droplet in the direction parallel to the electric field, H^^𝐻\hat{H}over^ start_ARG italic_H end_ARG is the length perpendicular to the electric field. Based on this equation, it is clear that the droplet deforms into a prolate (oblate) shape when the deformation factor D^^𝐷\hat{D}over^ start_ARG italic_D end_ARG is larger (smaller) than zero. Assuming that the two fluids are extremely viscous and conducting, Taylor’s small deformation theory predicts the deformation factor D^^𝐷\hat{D}over^ start_ARG italic_D end_ARG at steady state as [13]

D^=916CaE(2+R^)2[R^2+12S^+35(R^S^)2+3B^1+B^],^𝐷916subscript𝐶𝑎𝐸superscript2^𝑅2delimited-[]superscript^𝑅212^𝑆35^𝑅^𝑆23^𝐵1^𝐵\hat{D}=\frac{9}{{16}}\frac{{\mathop{Ca}\nolimits_{E}}}{{\mathop{(2+\hat{R})}% \nolimits^{2}}}[\mathop{\hat{R}}\nolimits^{2}+1-2\hat{S}+\frac{3}{5}(\hat{R}-% \hat{S})\frac{{2+3\hat{B}}}{{1+\hat{B}}}],over^ start_ARG italic_D end_ARG = divide start_ARG 9 end_ARG start_ARG 16 end_ARG divide start_ARG start_BIGOP italic_C italic_a end_BIGOP start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_ARG start_ARG start_BIGOP ( 2 + over^ start_ARG italic_R end_ARG ) end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ start_BIGOP over^ start_ARG italic_R end_ARG end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 1 - 2 over^ start_ARG italic_S end_ARG + divide start_ARG 3 end_ARG start_ARG 5 end_ARG ( over^ start_ARG italic_R end_ARG - over^ start_ARG italic_S end_ARG ) divide start_ARG 2 + 3 over^ start_ARG italic_B end_ARG end_ARG start_ARG 1 + over^ start_ARG italic_B end_ARG end_ARG ] , (82)

where R^^𝑅\hat{R}over^ start_ARG italic_R end_ARG, S^^𝑆\hat{S}over^ start_ARG italic_S end_ARG and B^^𝐵\hat{B}over^ start_ARG italic_B end_ARG denotes the ratios of the droplet’s conductivity, permittivity, and viscosity to the surrounding fluid, respectively. The electric capillary number CaE=𝜀v𝐸02r/γsubscript𝐶𝑎𝐸subscript𝜀𝑣superscriptsubscript𝐸02𝑟/𝛾\mathop{Ca}\nolimits_{E}={{\mathop{\varepsilon}\nolimits_{v}\mathop{E}% \nolimits_{0}^{2}r}\mathord{\left/{\vphantom{{\mathop{\varepsilon}\nolimits_{v% }\mathop{E}\nolimits_{0}^{2}r}\gamma}}\right.\kern-1.2pt}\gamma}start_BIGOP italic_C italic_a end_BIGOP start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT = italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT italic_E start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_r start_ID / end_ID italic_γ is the ratio of electric stress (uesubscript𝑢𝑒u_{e}italic_u start_POSTSUBSCRIPT italic_e end_POSTSUBSCRIPT) to capillary stress (ucsubscript𝑢𝑐u_{c}italic_u start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT). It should be noted that Taylor’s small deformation theory is constructed by ignoring the surface charge convection, which suggests that the multiphase interface is charged instantaneously. However, the above assumption cannot hold well when the charge relaxation time is comparable to the flow and capillary time scales. In such a case, Feng considered the influence of the surface charge convection and proposed an improved formation for the deformation factor D^^𝐷\hat{D}over^ start_ARG italic_D end_ARG [20]

D^=CaE3(1+2R^)2[R^2+R^+13S^].^𝐷subscript𝐶𝑎𝐸3superscript12^𝑅2delimited-[]superscript^𝑅2^𝑅13^𝑆\hat{D}=\frac{{\mathop{Ca}\nolimits_{E}}}{{3\mathop{(1+2\hat{R})}\nolimits^{2}% }}[\mathop{\hat{R}}\nolimits^{2}+\hat{R}+1-3\hat{S}].over^ start_ARG italic_D end_ARG = divide start_ARG start_BIGOP italic_C italic_a end_BIGOP start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT end_ARG start_ARG 3 start_BIGOP ( 1 + 2 over^ start_ARG italic_R end_ARG ) end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG [ start_BIGOP over^ start_ARG italic_R end_ARG end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + over^ start_ARG italic_R end_ARG + 1 - 3 over^ start_ARG italic_S end_ARG ] . (83)

Note that when the surface charge convection is incorporated, apart from the above mentioned electric capillary number CaE𝐶subscript𝑎𝐸Ca_{E}italic_C italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT, we also need to consider another three dimensionless numbers [24] , i.e.,

Re=𝑟2𝜌v𝜀v𝐸02𝜇v2,ReE=𝜀v2𝐸02𝜇v𝜎v,α=ω𝑘BT𝜇v𝜀v𝑟2𝐸02,formulae-sequenceResuperscript𝑟2subscript𝜌𝑣subscript𝜀𝑣superscriptsubscript𝐸02superscriptsubscript𝜇𝑣2formulae-sequencesubscriptRe𝐸superscriptsubscript𝜀𝑣2superscriptsubscript𝐸02subscript𝜇𝑣subscript𝜎𝑣𝛼𝜔subscript𝑘𝐵𝑇subscript𝜇𝑣subscript𝜀𝑣superscript𝑟2superscriptsubscript𝐸02{\mathop{\rm Re}\nolimits}=\frac{{\mathop{r}\nolimits^{2}\mathop{\rho}% \nolimits_{v}\mathop{\varepsilon}\nolimits_{v}\mathop{E}\nolimits_{0}^{2}}}{{% \mathop{\mu}\nolimits_{v}^{2}}},\qquad\mathop{{\mathop{\rm Re}\nolimits}}% \nolimits_{E}=\frac{{\mathop{\varepsilon}\nolimits_{v}^{2}\mathop{E}\nolimits_% {0}^{2}}}{{\mathop{\mu}\nolimits_{v}\mathop{\sigma}\nolimits_{v}}},\qquad% \alpha=\frac{{\omega\mathop{k}\nolimits_{B}T\mathop{\mu}\nolimits_{v}}}{{% \mathop{\varepsilon}\nolimits_{v}\mathop{r}\nolimits^{2}\mathop{E}\nolimits_{0% }^{2}}},roman_Re = divide start_ARG italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT italic_E start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , roman_Re start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT = divide start_ARG italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_E start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT italic_σ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT end_ARG , italic_α = divide start_ARG italic_ω italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT end_ARG start_ARG italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT italic_r start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_E start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (84)

where Re𝑅𝑒Reitalic_R italic_e is the Reynolds number defined as the ratio of electric force to viscous force, ReE𝑅subscript𝑒𝐸Re_{E}italic_R italic_e start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT is the electric Reynolds number used to quantify the surface charge convection, α𝛼\alphaitalic_α is the charge diffusion coefficient with kBsubscript𝑘𝐵k_{B}italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT, ω𝜔\omegaitalic_ω and T𝑇Titalic_T being the Boltzmann constant, charge mobility and fluid temperature, respectively.

In the simulation, the square width L𝐿Litalic_L is set to be L=8r𝐿8𝑟L=8ritalic_L = 8 italic_r, and the lattice size of the computation domain is fixed at 200×200200200200\times 200200 × 200. The periodic boundary condition is adopted in the horizontal direction, while on the top and bottom plates, the boundary conditions are realized using the half-way bounce-back scheme [44]. The density ratio is set to be 2.0 in order to avoid remarkable droplet migration driven by buoyancy [12]. The other parameters used in the phase-field simulation are chosen as γ=0.001𝛾0.001\gamma=0.001italic_γ = 0.001 (surface tension), W=5.0𝑊5.0W=5.0italic_W = 5.0 (interface width), M=0.1𝑀0.1M=0.1italic_M = 0.1 (mobility). Additionally, in order to compare the current numerical results with the theoretical and numerical solutions obtained from the leaky dielectric model, the value of ReE𝑅subscript𝑒𝐸Re_{E}italic_R italic_e start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT should be sufficiently small, i.e., ReE1much-less-thansubscriptRe𝐸1{{\mathop{\rm Re}\nolimits}_{E}}\ll 1roman_Re start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT ≪ 1, such that the ohmic conduction mechanism is dominant over the charge convection mechanism in the system. Also, according to the experiment, the Reynolds number is in the order of O(1)O1{\rm O}\left(1\right)roman_O ( 1 ), and the charge diffusion coefficient is about 104superscript104{10^{-4}}10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT. With the above analysis, unless otherwise stated, the following simulations in this subsection are all performed at ReE=103𝑅subscript𝑒𝐸superscript103Re_{E}=10^{-3}italic_R italic_e start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT, Re=1.0𝑅𝑒1.0Re=1.0italic_R italic_e = 1.0, α=104𝛼superscript104\alpha=10^{-4}italic_α = 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT and CaE=0.2𝐶subscript𝑎𝐸0.2Ca_{E}=0.2italic_C italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT = 0.2.

Refer to caption
Refer to caption
Refer to caption
Figure 6: Stream lines (left half), velocity vectors (right half) and the profile of the droplet (red lines) for three representative steady-state droplets: (a) prolate deformation obtained at (R^,S^)=(5,0.5)^𝑅^𝑆50.5(\hat{R},\hat{S})=(5,0.5)( over^ start_ARG italic_R end_ARG , over^ start_ARG italic_S end_ARG ) = ( 5 , 0.5 ); (b) slight prolate deformation at (R^,S^)=(5,5)^𝑅^𝑆55(\hat{R},\hat{S})=(5,5)( over^ start_ARG italic_R end_ARG , over^ start_ARG italic_S end_ARG ) = ( 5 , 5 ); (c) oblate deformation at (R^,S^)=(5,15)^𝑅^𝑆515(\hat{R},\hat{S})=(5,15)( over^ start_ARG italic_R end_ARG , over^ start_ARG italic_S end_ARG ) = ( 5 , 15 ), , in which the simulated parameters are chosen as CaE=0.2𝐶subscript𝑎𝐸0.2Ca_{E}=0.2italic_C italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT = 0.2, 𝜌l/𝜌v=2.0subscript𝜌𝑙/subscript𝜌𝑣2.0{{\mathop{\rho}\nolimits_{l}}\mathord{\left/{\vphantom{{\mathop{\rho}\nolimits% _{l}}{\mathop{\rho}\nolimits_{v}}}}\right.\kern-1.2pt}{\mathop{\rho}\nolimits_% {v}}}=2.0italic_ρ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 2.0, 𝜇l/𝜇v=1.0subscript𝜇𝑙/subscript𝜇𝑣1.0{{\mathop{\mu}\nolimits_{l}}\mathord{\left/{\vphantom{{\mathop{\mu}\nolimits_{% l}}{\mathop{\mu}\nolimits_{v}}}}\right.\kern-1.2pt}{\mathop{\mu}\nolimits_{v}}% }=1.0italic_μ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 1.0.
Refer to caption
Refer to caption
Refer to caption
Figure 7: Variation in electric field lines (left half), electric potential lines (right half), charge distribution (colored parts) and the electric field force (red arrows) for three different values of permittivity ratios: (a) S^=0.5^𝑆0.5\hat{S}=0.5over^ start_ARG italic_S end_ARG = 0.5; (b) S^=5.0^𝑆5.0\hat{S}=5.0over^ start_ARG italic_S end_ARG = 5.0; (c) S^=15.0^𝑆15.0\hat{S}=15.0over^ start_ARG italic_S end_ARG = 15.0, in which the simulated parameters are chosen as CaE=0.2𝐶subscript𝑎𝐸0.2Ca_{E}=0.2italic_C italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT = 0.2, 𝜌l/𝜌v=2.0subscript𝜌𝑙/subscript𝜌𝑣2.0{{\mathop{\rho}\nolimits_{l}}\mathord{\left/{\vphantom{{\mathop{\rho}\nolimits% _{l}}{\mathop{\rho}\nolimits_{v}}}}\right.\kern-1.2pt}{\mathop{\rho}\nolimits_% {v}}}=2.0italic_ρ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 2.0, 𝜇l/𝜇v=1.0subscript𝜇𝑙/subscript𝜇𝑣1.0{{\mathop{\mu}\nolimits_{l}}\mathord{\left/{\vphantom{{\mathop{\mu}\nolimits_{% l}}{\mathop{\mu}\nolimits_{v}}}}\right.\kern-1.2pt}{\mathop{\mu}\nolimits_{v}}% }=1.0italic_μ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 1.0, R^=5.0^𝑅5.0\hat{R}=5.0over^ start_ARG italic_R end_ARG = 5.0.

Fig. 6 presents the streamlines and velocity vectors for different permittivity ratios S^^𝑆\hat{S}over^ start_ARG italic_S end_ARG with conductivity ratio R^^𝑅\hat{R}over^ start_ARG italic_R end_ARG being fixed at 5.0. As shown in Fig. 6, it is seen that for the case of (R^,S^)=(5,0.5)^𝑅^𝑆50.5(\hat{R},\hat{S})=(5,0.5)( over^ start_ARG italic_R end_ARG , over^ start_ARG italic_S end_ARG ) = ( 5 , 0.5 ), the droplet deforms into a prolate shape with four symmetrical and counter-rotating vortices formed inside and outside of the droplet, and the vector diagram in the first quadrant of the droplet displays that the circulation direction is from the equator to the poles along the droplet surface. When (R^,S^)=(5,5)^𝑅^𝑆55(\hat{R},\hat{S})=(5,5)( over^ start_ARG italic_R end_ARG , over^ start_ARG italic_S end_ARG ) = ( 5 , 5 ), the deformation of the droplet is insignificant, and flow patterns are similar to those observed at (R^,S^)=(5,0.5)^𝑅^𝑆50.5(\hat{R},\hat{S})=(5,0.5)( over^ start_ARG italic_R end_ARG , over^ start_ARG italic_S end_ARG ) = ( 5 , 0.5 ) (see Fig. 6). However, when the ratio of the electric permittivities (S^=15.0^𝑆15.0\hat{S}=15.0over^ start_ARG italic_S end_ARG = 15.0) is larger than the ratio of the electrical conductivities (R^=5^𝑅5\hat{R}=5over^ start_ARG italic_R end_ARG = 5), an oblate droplet with reversed flow direction is obtained. To have a better understanding on the influence of the electric field, Fig. 7 illustrates the distributions of the electric potential, charge density, electric field lines, as well as the electric force direction for different permittivity ratios. It is evident that due to the discontinuity in the permittivity at the interface, the electric potential and electric field lines are distorted in this region. Specifically, a fluid with high permittivity polarizes more in response to an applied electric field, resulting in denser potential lines inside the droplet for a relatively smaller permittivity ratio. Additionally, it is noted that for all cases, the maximum charge density appears at the poles while it is negligible at the equator. In particular, we also found that the position of the positive and negative charges is related to the droplet’s shape, which directly leads to a difference in the direction of the corresponding electric field force.

Table 1: Quantitative comparison of the results from the present work with theoretical solutions and existing data.
R^^𝑅\hat{R}over^ start_ARG italic_R end_ARG S^^𝑆\hat{S}over^ start_ARG italic_S end_ARG CaE𝐶subscript𝑎𝐸Ca_{E}italic_C italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT Deformation factor (D^^𝐷\hat{D}over^ start_ARG italic_D end_ARG)
Present Eq.(82) Eq.(83) Ref. [24]
5 5 0.2 0.03503 0.03670 0.02960 0.03150
5 60 0.2 -0.26304 -0.40520 -0.27590 -0.27510
1 2 0.2 -0.04102 -0.04380 -0.05000 -0.05390
1.75 3.5 0.1 -0.01968 -0.02230 -0.02070 -0.02000
3.25 3.5 0.1 0.00848 0.00850 0.00800 0.00900
4.75 3.5 0.1 0.02236 0.02280 0.01800 0.02090

To quantify the results, Table. 1 presents the deformation factor of the droplet at various parameters, and it can be seen that our numerical results fit well with these existing data [24] . Moreover, in order to provide a more convincing evaluation, we also conduct additional simulations using various electric capillary numbers, permittivity ratios, and conductivity ratios. We then plot the results in Fig. 8 for comparison. It is evident from this figure that the numerical results from our current model align well with both existing theoretical and numerical benchmark data. These results demonstrate the viability of the proposed approach for droplet deformation in an electric field and indicate that the present model can accurately simulate two-phase EHD flows.

Refer to caption
Refer to caption
Figure 8: (a) Effect of electric capillary number CaE𝐶subscript𝑎𝐸Ca_{E}italic_C italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT on droplet deformation D^^𝐷\hat{D}over^ start_ARG italic_D end_ARG at different conductivity ratios R^^𝑅\hat{R}over^ start_ARG italic_R end_ARG with S^=3.5^𝑆3.5\hat{S}=3.5over^ start_ARG italic_S end_ARG = 3.5; (b) effect of different permittivity ratios on deformation factor with CaE=0.2𝐶subscript𝑎𝐸0.2Ca_{E}=0.2italic_C italic_a start_POSTSUBSCRIPT italic_E end_POSTSUBSCRIPT = 0.2 and R^=5.0^𝑅5.0\hat{R}=5.0over^ start_ARG italic_R end_ARG = 5.0.

7.2 Droplet detachment in reversed electrowetting

Refer to caption
Figure 9: Schematic diagram of reversed electrowetting, in which a nonconducting droplet is immersed in another conducting fluid.

In classical electrowetting experiments, the Lippmann equation is often used to describe the wetting behaviour of a conductive droplet. According to this equation, the contact angle of the droplet decreases as the applied voltage increases. In this setting, when it comes to removing droplets from a solid surface, a critical process in applications such as enhanced oil recovery [79] and surface cleaning [75] , the only way to do so is by suddenly releasing the applied voltage [76]. However, due to the strict critical conditions for droplet detachment, this strategy would fail in many conditions, such as the applied voltage being too low or the process not being fast enough to restore sufficient energy [77] . To this end, recently, some researchers manipulated the non-conductive droplets by using a ”reversed” electrowetting effect [78] . The setup of this new approach is depicted in Fig. 9, where a droplet is initially placed on the bottom wall and submerged in water. The droplet’s working liquid is a non-conductive fluid, similar to silicone oil. When a potential difference is applied between the bottom wall and the surrounding water, the apparent contact angle of the droplet tends to increase, leading to the contact radius decrease. In this manner, the droplet is expected to separate from the substrate as long as the applied voltage is sufficiently large. Owing to the wettability of the substrate in this setup being increased in applied voltage, the phenomenon is referred to as reversed electrowetting.

We now turn to simulate the detachment of a non-conductivity droplet in reversed electrowetting. Our simulations are performed on a computational domain of 200×100200100200\times 100200 × 100, and a semicircular droplet with a radius of 25 is initialized in contact with the bottom wall. The density ratio, dynamic viscosity ratio and the electric permittivity ratio are chosen as ρl/ρv=1.069subscript𝜌𝑙/subscript𝜌𝑣1.069{{{{\rho_{l}}}\mathord{\left/{\vphantom{{{\rho_{l}}}\rho}}\right.\kern-1.2pt}% \rho}_{v}}=1.069italic_ρ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_ρ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 1.069, μl/μv=0.097subscript𝜇𝑙/subscript𝜇𝑣0.097{{{{\mu_{l}}}\mathord{\left/{\vphantom{{{\mu_{l}}}\mu}}\right.\kern-1.2pt}\mu}% _{v}}=0.097italic_μ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 0.097, and εl/εv=32.0subscript𝜀𝑙/subscript𝜀𝑣32.0{{{{\varepsilon_{l}}}\mathord{\left/{\vphantom{{{\varepsilon_{l}}}\varepsilon}% }\right.\kern-1.2pt}\varepsilon}_{v}}=32.0italic_ε start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 32.0, which approach those of realistic water-oil system adopted in the experiment [79]. The initial droplet contact angle in the absence of voltage is set to be 40superscript40{40^{\circ}}40 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT. Different electric potentials are imposed between the top (φHsubscript𝜑𝐻{\varphi_{H}}italic_φ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT) and bottom walls (φLsubscript𝜑𝐿{\varphi_{L}}italic_φ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT), where φH>φLsubscript𝜑𝐻subscript𝜑𝐿{\varphi_{H}}>{\varphi_{L}}italic_φ start_POSTSUBSCRIPT italic_H end_POSTSUBSCRIPT > italic_φ start_POSTSUBSCRIPT italic_L end_POSTSUBSCRIPT. The boundary conditions of the mesoscopic particles, as well as the other unmentioned parameters, are the same as those adopted in Sec. 4. Moreover, to realize the reversed electrowetting phenomenon, the effective surface tension for liquid-solid should modified as γ^sl=γsl+η/dsubscript^𝛾𝑠𝑙subscript𝛾𝑠𝑙𝜂/𝑑{{\hat{\gamma}}_{sl}}={\gamma_{sl}}+{\eta\mathord{\left/{\vphantom{\eta d}}% \right.\kern-1.2pt}d}over^ start_ARG italic_γ end_ARG start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT = italic_γ start_POSTSUBSCRIPT italic_s italic_l end_POSTSUBSCRIPT + italic_η start_ID / end_ID italic_d.

Refer to caption
Refer to caption
Figure 10: Comparison of shapes during the process of the reversed electrowetting drives a droplet detachment from the substrate: (a) experimental results [79]; (b) the present numerical results (the blue arrows represent the direction of velocity.
Refer to caption
Refer to caption
Figure 11: Evolution of the normalized contact radius on the substrate over time. The present numerical results (black lines) are compared with the experimental results (blue circles) and the available data (red squares): (a) 𝜇l/𝜇v=0.097subscript𝜇𝑙/subscript𝜇𝑣0.097{{\mathop{\mu}\nolimits_{l}}\mathord{\left/{\vphantom{{\mathop{\mu}\nolimits_{% l}}{\mathop{\mu}\nolimits_{v}}}}\right.\kern-1.2pt}{\mathop{\mu}\nolimits_{v}}% }=0.097italic_μ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 0.097; (b) 𝜇l/𝜇v=0.022subscript𝜇𝑙/subscript𝜇𝑣0.022{{\mathop{\mu}\nolimits_{l}}\mathord{\left/{\vphantom{{\mathop{\mu}\nolimits_{% l}}{\mathop{\mu}\nolimits_{v}}}}\right.\kern-1.2pt}{\mathop{\mu}\nolimits_{v}}% }=0.022italic_μ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT start_ID / end_ID italic_μ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT = 0.022.

The comparison of the current simulation with the experiment conducted by Weng et al. [79] is shown in Fig. 10, which illustrates the variation of the droplet shape and contact radius. As seen from this figure, our simulated results are comparable to those of Weng et al.’s experiment [79] . In particular, we find that in the reversed electrowetting, the droplet tends to contract under the action of the electric field. Owing to the storing surface energy of the droplet being large enough, leading to detachment from the solid surface. Fig. 10 also depicts the corresponding velocity vectors near the droplet, and it can be seen that there are two counter-rotating vortexes at the front of the droplet interface. For a quantitative analysis, we recorded the time evolution of the contact radius and compared it with experimental studies and numerical data computed by COMSOL. As shown in Fig. 11, our results are consistent with previous observations, indicating that our model is a suitable candidate for simulating two-phase EHD flows.

8 Conclusion

In this paper, we first use Onsager’s variational principle to develop a thermodynamically consistent phase-field model for two-phase EHD flows. This model allows a two-phase incompressible fluid to have different electrohydrodynamic properties for each phase, such as densities, viscosities, permittivities, and conductivities. In particular, the deduced model incorporates the influence of the surface charge convection, which is often overlooked in previous studies. This feature enables the model to be applied to problems where bulk charge conduction and convection are relevant, such as the characterization of the cone-to-jet transition region in EHD cone-jet electrosprays.

After obtaining the thermodynamically consistent phase-field model, the lattice Boltzmann method is utilized to simulate the two-phase EHD flows. In this approach, it adopts four distribution functions for solving the Cahn-Hilliard equation, the electric potential equation, the Nernst-Planck equation, and the hydrodynamic equations. Additionally, we have developed an improved LB method for the electric potential equation to account for two-phase EHD flows with a larger permittivity ratio. The performance of the LB method is evaluated by simulating EHD flows with two superimposed planar fluids and the charge diffusion of a Gaussian bell. The numerical results indicate that the developed LB method has a second-order convergence rate in predicting electric potential and charge density. Moreover, we consider the equilibrium interface profiles of a conductivity droplet and found that the resulting apparent contact agrees well with those derived by the Lippmann equation. To further test the current LB method’s capability in practical problems, we also conduct some simulations for two typical two-phase EHD problems, including droplet deformation under an electric field and droplet detachment in reversed electrowetting. Based on the comparison between our numerical results and the analytical solutions, as well as experimental/existing numerical data, it is found that the present method enables us to derive comparative results for predicting droplet deformation and contact radius. In conclusion, the simulated results demonstrate that the current method is feasible for simulating two-phase EHD flows.

CRediT authorship contribution statement

Fang Xiong: Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. Lei Wang: Writing – review & editing, Writing – original draft, Supervision, Conceptualization. Jiangxu Huang: Supervision, Software, Methodology, Resources. Kang Luo: Supervision, Software, Methodology, Funding acquisition.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgements

This work is financially supported by the National Natural Science Foundation of China (Grant No. 51906051).

Appendix A Chapman-Enskog analysis on electric potential equation

To show that the electric potential equation Eq. (32) can be recovered with no deviation term using the proposed LB equation Eq. (56), the Chapman-Enskog analysis is performed in this appendix. Base on the definitions of discrete lattice velocity 𝐜isubscript𝐜𝑖\mathbf{c}_{i}bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and the equilibrium distribution function hieqsuperscriptsubscript𝑖𝑒𝑞h_{i}^{eq}italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT given by Eq. (37) and Eq. (57), we can easily obtain

ihieq=0,i𝐜ihieq=0,i𝐜i𝐜ihieq=cs2φ.formulae-sequencesubscript𝑖superscriptsubscript𝑖𝑒𝑞0formulae-sequencesubscript𝑖subscript𝐜𝑖superscriptsubscript𝑖𝑒𝑞0subscript𝑖subscript𝐜𝑖subscript𝐜𝑖superscriptsubscript𝑖𝑒𝑞superscriptsubscript𝑐𝑠2𝜑\sum\limits_{i}{h_{i}^{eq}}=0,\qquad\sum\limits_{i}{\mathbf{c}_{i}h_{i}^{eq}}=% 0,\qquad\sum\limits_{i}{\mathbf{c}_{i}\mathbf{c}_{i}h_{i}^{eq}}=c_{s}^{2}\varphi.∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT = 0 , ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT = 0 , ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_φ . (85)

In addition, the distribution function, the derivatives of time and space, as well as the source term can be expanded in consecutive scales of ξ𝜉\xiitalic_ξ:

hi=hi(0)+ξhi(1)+𝜉2hi(2)+,t=ξt1+𝜉2t2,=ξ1,q=𝜉2𝑞(2),\begin{gathered}h_{i}=h_{i}^{(0)}+\xi h_{i}^{(1)}+\mathop{\xi}\nolimits^{2}h_{% i}^{(2)}+\cdots,\\ \frac{\partial}{{\partial t}}=\xi\frac{\partial}{{\partial t_{1}}}+\mathop{\xi% }\nolimits^{2}\frac{\partial}{{\partial t_{2}}},\qquad\nabla=\xi\mathop{\nabla% }\nolimits_{1},\qquad q=\mathop{\xi}\nolimits^{2}\mathop{q}\nolimits^{(2)},% \end{gathered}start_ROW start_CELL italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT + italic_ξ italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT + ⋯ , end_CELL end_ROW start_ROW start_CELL divide start_ARG ∂ end_ARG start_ARG ∂ italic_t end_ARG = italic_ξ divide start_ARG ∂ end_ARG start_ARG ∂ italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG + italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG ∂ end_ARG start_ARG ∂ italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG , ∇ = italic_ξ ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_q = italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT , end_CELL end_ROW (86)

where ξ𝜉\xiitalic_ξ is the expansion parameter. Then, applying the Taylor expansion to Eq. (56) at time t𝑡titalic_t and space 𝐱𝐱\mathbf{x}bold_x, we get

ΔtDihi+Δt22Di2hi=1τh[hihieq]ε^(ϕ)𝑐s2𝜔i𝐜i1φτh+Δtϖiq,Δ𝑡subscript𝐷𝑖subscript𝑖superscriptΔ𝑡22superscriptsubscript𝐷𝑖2subscript𝑖1subscript𝜏delimited-[]subscript𝑖superscriptsubscript𝑖𝑒𝑞^𝜀italic-ϕsuperscriptsubscript𝑐𝑠2subscript𝜔𝑖subscript𝐜𝑖subscript1𝜑subscript𝜏Δ𝑡subscriptitalic-ϖ𝑖𝑞\Delta tD_{i}h_{i}+\frac{{\mathop{\Delta t}\nolimits^{2}}}{2}D_{i}^{2}h_{i}=-% \frac{1}{{\tau_{h}}}[h_{i}-h_{i}^{eq}]-\frac{\hat{\varepsilon}\left(\phi\right% )}{{\mathop{c}\nolimits_{s}^{2}}}\frac{{\mathop{\omega}\nolimits_{i}\mathbf{c}% _{i}{\nabla_{1}}\varphi}}{{\tau_{h}}}{\rm{+}}\Delta t{\varpi_{i}}q,roman_Δ italic_t italic_D start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + divide start_ARG start_BIGOP roman_Δ italic_t end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 end_ARG italic_D start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG [ italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ] - divide start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG divide start_ARG italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG + roman_Δ italic_t italic_ϖ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_q , (87)

where 𝐷i=ξ𝐷1i+𝜉2𝑡2subscript𝐷𝑖𝜉subscript𝐷1𝑖superscript𝜉2subscriptsubscript𝑡2\mathop{D}\nolimits_{i}=\xi\mathop{D}\nolimits_{1i}+\mathop{\xi}\nolimits^{2}% \mathop{\partial}\nolimits_{\mathop{t}\nolimits_{2}}italic_D start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_ξ italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT + italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT with 𝐷1i=𝑡1+𝐜i1subscript𝐷1𝑖subscriptsubscript𝑡1subscript𝐜𝑖subscript1\mathop{D}\nolimits_{1i}=\mathop{\partial}\nolimits_{\mathop{t}\nolimits_{1}}+% \mathbf{c}_{i}\cdot\mathop{\nabla}\nolimits_{1}italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT = ∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT + bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ⋅ ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. Substituting Eq. (86) into Eq. (87) yields the following equations in the successive order of

O(𝜉0):i(0)=ieq,:𝑂superscript𝜉0superscriptsubscript𝑖0superscriptsubscript𝑖𝑒𝑞\displaystyle O(\mathop{\xi}\nolimits^{0}):\mathop{h}\nolimits_{i}^{(0)}=% \mathop{h}\nolimits_{i}^{eq},italic_O ( italic_ξ start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT ) : italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT = italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT , (88)
O(𝜉1):𝐷1ii(0)=1𝜏hΔti(1)ε^(ϕ)𝑐s2Δt𝜔i𝐜i1φτh,:𝑂superscript𝜉1subscript𝐷1𝑖superscriptsubscript𝑖01subscript𝜏Δ𝑡superscriptsubscript𝑖1^𝜀italic-ϕsuperscriptsubscript𝑐𝑠2Δ𝑡subscript𝜔𝑖subscript𝐜𝑖subscript1𝜑subscript𝜏\displaystyle O(\mathop{\xi}\nolimits^{1}):\mathop{D}\nolimits_{1i}\mathop{h}% \nolimits_{i}^{(0)}=-\frac{1}{{\mathop{\tau}\nolimits_{h}\Delta t}}\mathop{h}% \nolimits_{i}^{(1)}-\frac{\hat{\varepsilon}\left(\phi\right)}{{\mathop{c}% \nolimits_{s}^{2}\Delta t}}\frac{{\mathop{\omega}\nolimits_{i}\mathbf{c}_{i}% \nabla_{1}\varphi}}{{\tau_{h}}},italic_O ( italic_ξ start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT ) : italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT = - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t end_ARG italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT - divide start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_t end_ARG divide start_ARG italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG , (89)
O(𝜉2):𝑡2i(0)+𝐷1ii(1)+Δt2𝐷1i2i(0)=1𝜏hΔti(2)+ϖi𝑞(2).:𝑂superscript𝜉2subscriptsubscript𝑡2superscriptsubscript𝑖0subscript𝐷1𝑖superscriptsubscript𝑖1Δ𝑡2superscriptsubscript𝐷1𝑖2superscriptsubscript𝑖01subscript𝜏Δ𝑡superscriptsubscript𝑖2subscriptitalic-ϖ𝑖superscript𝑞2\displaystyle O(\mathop{\xi}\nolimits^{2}):\mathop{\partial}\nolimits_{\mathop% {t}\nolimits_{2}}\mathop{h}\nolimits_{i}^{(0)}+\mathop{D}\nolimits_{1i}\mathop% {h}\nolimits_{i}^{(1)}+\frac{{\Delta t}}{2}\mathop{\mathop{D}\nolimits_{1i}}% \nolimits^{2}\mathop{h}\nolimits_{i}^{(0)}=-\frac{1}{{\mathop{\tau}\nolimits_{% h}\Delta t}}\mathop{h}\nolimits_{i}^{(2)}+{\varpi_{i}}\mathop{q}\nolimits^{(2)}.italic_O ( italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) : ∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT + italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + divide start_ARG roman_Δ italic_t end_ARG start_ARG 2 end_ARG start_BIGOP italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT = - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t end_ARG italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT + italic_ϖ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_q start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT . (90)

Combining Eq. (85), Eq. (86) and Eq. (88), we get

ii(k)=0(k>1).subscript𝑖superscriptsubscript𝑖𝑘0𝑘1\sum\limits_{i}{\mathop{h}\nolimits_{i}^{(k)}}=0\quad(k>1).∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_k ) end_POSTSUPERSCRIPT = 0 ( italic_k > 1 ) . (91)

Moreover, applying Eq. (89) to Eq. (90), one can obtain

𝑡2i(0)+𝐷1i(112𝜏h)i(1)𝐷1iε^(ϕ)2cs2𝜔i𝐜i1φ𝜏h=1𝜏hΔti(2)+ϖi𝑞(2).subscriptsubscript𝑡2superscriptsubscript𝑖0subscript𝐷1𝑖112subscript𝜏superscriptsubscript𝑖1subscript𝐷1𝑖^𝜀italic-ϕsuperscriptsubscript2𝑐𝑠2subscript𝜔𝑖subscript𝐜𝑖subscript1𝜑subscript𝜏1subscript𝜏Δ𝑡superscriptsubscript𝑖2subscriptitalic-ϖ𝑖superscript𝑞2\mathop{\partial}\nolimits_{\mathop{t}\nolimits_{2}}\mathop{h}\nolimits_{i}^{(% 0)}+\mathop{D}\nolimits_{1i}(1-\frac{1}{{2\mathop{\tau}\nolimits_{h}}})\mathop% {h}\nolimits_{i}^{(1)}-\mathop{D}\nolimits_{1i}\frac{\hat{\varepsilon}\left(% \phi\right)}{{\mathop{2c}\nolimits_{s}^{2}}}\frac{{\mathop{\omega}\nolimits_{i% }\mathbf{c}_{i}\mathop{\nabla}\nolimits_{1}\varphi}}{{\mathop{\tau}\nolimits_{% h}}}=-\frac{1}{{\mathop{\tau}\nolimits_{h}\Delta t}}\mathop{h}\nolimits_{i}^{(% 2)}+{\varpi_{i}}\mathop{q}\nolimits^{(2)}.∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT + italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT ( 1 - divide start_ARG 1 end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG ) italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT - italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT divide start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) end_ARG start_ARG start_BIGOP 2 italic_c end_BIGOP start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG divide start_ARG italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG = - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t end_ARG italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT + italic_ϖ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_q start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT . (92)

After a summation of Eq. (92) over i𝑖iitalic_i, we can get the scale equation as

(112𝜏h)1i𝐜ii(1)1ε^(ϕ)2𝜏h1φ=𝑞(2).112subscript𝜏subscript1subscript𝑖subscript𝐜𝑖superscriptsubscript𝑖1subscript1^𝜀italic-ϕ2subscript𝜏subscript1𝜑superscript𝑞2(1-\frac{1}{{2\mathop{\tau}\nolimits_{h}}})\mathop{\nabla}\nolimits_{1}\sum% \limits_{i}{\mathbf{c}_{i}}\mathop{h}\nolimits_{i}^{(1)}-\mathop{\nabla}% \nolimits_{1}\frac{\hat{\varepsilon}\left(\phi\right)}{{2\mathop{\tau}% \nolimits_{h}}}\mathop{\nabla}\nolimits_{1}\varphi=\mathop{q}\nolimits^{(2)}.( 1 - divide start_ARG 1 end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG ) ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT - ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT divide start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ = italic_q start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT . (93)

On the basis of Eq. (89), we have

i𝐜ii(1)=𝜏hΔt[ε^(ϕ)𝑐s2Δt𝜏hi𝜔i𝐜i𝐜i1φ1i𝐜i𝐜ii(0)𝑡1i𝐜ii(0)]=ε^(ϕ)1φ𝜏hΔt𝑐s21φ.subscript𝑖subscript𝐜𝑖superscriptsubscript𝑖1subscript𝜏Δ𝑡delimited-[]^𝜀italic-ϕsuperscriptsubscript𝑐𝑠2Δ𝑡subscript𝜏subscript𝑖subscript𝜔𝑖subscript𝐜𝑖subscript𝐜𝑖subscript1𝜑subscript1subscript𝑖subscript𝐜𝑖subscript𝐜𝑖superscriptsubscript𝑖0subscriptsubscript𝑡1subscript𝑖subscript𝐜𝑖superscriptsubscript𝑖0^𝜀italic-ϕsubscript1𝜑subscript𝜏Δ𝑡superscriptsubscript𝑐𝑠2subscript1𝜑\begin{split}\sum\limits_{i}{\mathbf{c}_{i}}\mathop{h}\nolimits_{i}^{(1)}&=% \mathop{\tau}\nolimits_{h}\Delta t[-\frac{\hat{\varepsilon}\left(\phi\right)}{% {\mathop{c}\nolimits_{s}^{2}\Delta t\mathop{\tau}\nolimits_{h}}}\sum\limits_{i% }{\mathop{\omega}\nolimits_{i}\mathbf{c}_{i}\mathbf{c}_{i}{\mathop{\nabla}% \nolimits_{1}}\varphi}-\mathop{\nabla}\nolimits_{1}\sum\limits_{i}{\mathbf{c}_% {i}\mathbf{c}_{i}\mathop{h}\nolimits_{i}^{(0)}}-\mathop{\partial}\nolimits_{% \mathop{t}\nolimits_{1}}\sum\limits_{i}{\mathbf{c}_{i}\mathop{h}\nolimits_{i}^% {(0)}}]\\ &=-{\hat{\varepsilon}\left(\phi\right)}{\mathop{\nabla}\nolimits_{1}}\varphi-% \mathop{\tau}\nolimits_{h}\Delta t\mathop{c}\nolimits_{s}^{2}\mathop{\nabla}% \nolimits_{1}\varphi.\end{split}start_ROW start_CELL ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT end_CELL start_CELL = italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t [ - divide start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_t italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT end_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ - ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT - ∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT ] end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL = - over^ start_ARG italic_ε end_ARG ( italic_ϕ ) ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ - italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ . end_CELL end_ROW (94)

With the aid of the above equation, Eq. (93) can be rewritten as

1[𝑐s2(𝜏h12)Δt1φ]=𝑞(2)1(ε^(ϕ)1φ).\mathop{\nabla}\nolimits_{1}[\mathop{c}\nolimits_{s}^{2}(\mathop{\tau}% \nolimits_{h}-\frac{1}{2})\Delta t\mathop{\nabla}\nolimits_{1}\varphi]=-% \mathop{q}\nolimits^{(2)}-\mathop{\nabla}\nolimits_{1}\mathop{(\hat{% \varepsilon}\left(\phi\right)\mathop{\nabla}\nolimits_{1}\varphi}).∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ) roman_Δ italic_t ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ ] = - italic_q start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT - ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_BIGOP ( over^ start_ARG italic_ε end_ARG ( italic_ϕ ) ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_φ end_BIGOP ) . (95)

Multiplying 𝜉2superscript𝜉2\mathop{\xi}\nolimits^{2}italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT on both sides of Eq. (95), we can derive the electric potential equation Eq. (55) with εvsubscript𝜀𝑣\varepsilon_{v}italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT being

𝜀v=cs2(𝜏h12)Δt.subscript𝜀𝑣superscriptsubscriptabsent𝑐𝑠2subscript𝜏12Δ𝑡\mathop{\varepsilon}\nolimits_{v}\mathop{=c}\nolimits_{s}^{2}(\mathop{\tau}% \nolimits_{h}-\frac{1}{2})\Delta t.italic_ε start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT start_BIGOP = italic_c end_BIGOP start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ) roman_Δ italic_t . (96)

Particularly, according to Eq. (89), we obtain a local scheme for computing φ𝜑\nabla\varphi∇ italic_φ, which is expressed as

φ=i𝐜ii(1)ε^(ϕ)+𝑐s2𝜏hΔt=i𝐜iiε^(ϕ)+𝑐s2𝜏hΔt.𝜑subscript𝑖subscript𝐜𝑖superscriptsubscript𝑖1^𝜀italic-ϕsuperscriptsubscript𝑐𝑠2subscript𝜏Δ𝑡subscript𝑖subscript𝐜𝑖subscript𝑖^𝜀italic-ϕsuperscriptsubscript𝑐𝑠2subscript𝜏Δ𝑡\nabla\varphi=-\frac{{\sum\limits_{i}{\mathbf{c}_{i}}\mathop{h}\nolimits_{i}^{% (1)}}}{{\hat{\varepsilon}\left(\phi\right)+\mathop{c}\nolimits_{s}^{2}\mathop{% \tau}\nolimits_{h}\Delta t}}=-\frac{{\sum\limits_{i}{\mathbf{c}_{i}}\mathop{h}% \nolimits_{i}}}{{\hat{\varepsilon}\left(\phi\right)+\mathop{c}\nolimits_{s}^{2% }\mathop{\tau}\nolimits_{h}\Delta t}}.∇ italic_φ = - divide start_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT end_ARG start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) + italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t end_ARG = - divide start_ARG ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_h start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_ARG start_ARG over^ start_ARG italic_ε end_ARG ( italic_ϕ ) + italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_τ start_POSTSUBSCRIPT italic_h end_POSTSUBSCRIPT roman_Δ italic_t end_ARG . (97)

Appendix B Chapman-Enskog analysis on Nernst-Planck equation

Similar to the process derivation, the equilibrium distribution function for Nernst-Planck equation satisfy

ilieq=q,i𝐜ilieq=q𝐮,i𝐜i𝐜ilieq=cs2q.formulae-sequencesubscript𝑖superscriptsubscript𝑙𝑖𝑒𝑞𝑞formulae-sequencesubscript𝑖subscript𝐜𝑖superscriptsubscript𝑙𝑖𝑒𝑞𝑞𝐮subscript𝑖subscript𝐜𝑖subscript𝐜𝑖superscriptsubscript𝑙𝑖𝑒𝑞superscriptsubscript𝑐𝑠2𝑞\sum\limits_{i}{l_{i}^{eq}}=q,\qquad\sum\limits_{i}{\mathbf{c}_{i}l_{i}^{eq}}=% q\mathbf{u},\qquad\sum\limits_{i}{\mathbf{c}_{i}\mathbf{c}_{i}l_{i}^{eq}}=c_{s% }^{2}q.∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT = italic_q , ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT = italic_q bold_u , ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_c start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT = italic_c start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q . (98)

Applying the Taylor series expansion to Eq. (62), we have

Δt𝐷i𝑙i+Δt22𝐷i2𝑙i=1𝜏l(𝑙i𝑙ieq)+Δt𝑆i+Δt𝑇i.Δ𝑡subscript𝐷𝑖subscript𝑙𝑖superscriptΔ𝑡22superscriptsubscript𝐷𝑖2subscript𝑙𝑖1subscript𝜏𝑙subscript𝑙𝑖superscriptsubscript𝑙𝑖𝑒𝑞Δ𝑡subscript𝑆𝑖Δ𝑡subscript𝑇𝑖\Delta t\mathop{D}\nolimits_{i}\mathop{l}\nolimits_{i}+\frac{{\mathop{\Delta t% }\nolimits^{2}}}{2}\mathop{\mathop{D}\nolimits_{i}}\nolimits^{2}\mathop{l}% \nolimits_{i}=-\frac{1}{{\mathop{\tau}\nolimits_{l}}}(\mathop{l}\nolimits_{i}-% \mathop{l}\nolimits_{i}^{eq})+\Delta t\mathop{S}\nolimits_{i}+\Delta t\mathop{% T}\nolimits_{i}.roman_Δ italic_t italic_D start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + divide start_ARG start_BIGOP roman_Δ italic_t end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 end_ARG start_BIGOP italic_D start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = - divide start_ARG 1 end_ARG start_ARG italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG ( italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT ) + roman_Δ italic_t italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + roman_Δ italic_t italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT . (99)

By introducing the following expansions

𝑙i=𝑙i(0)+ξ𝑙i(1)+𝜉2𝑙i(2)+,subscript𝑙𝑖superscriptsubscript𝑙𝑖0𝜉superscriptsubscript𝑙𝑖1superscript𝜉2superscriptsubscript𝑙𝑖2\displaystyle\mathop{l}\nolimits_{i}=\mathop{l}\nolimits_{i}^{(0)}+\xi\mathop{% l}\nolimits_{i}^{(1)}+\mathop{\xi}\nolimits^{2}\mathop{l}\nolimits_{i}^{(2)}+\cdots,italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT + italic_ξ italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT + ⋯ , (100)
t=ε𝑡1+𝜉2𝑡2,=ξ1,formulae-sequencesubscript𝑡𝜀subscriptsubscript𝑡1superscript𝜉2subscriptsubscript𝑡2𝜉subscript1\displaystyle\mathop{\partial}\nolimits_{t}=\varepsilon\mathop{\partial}% \nolimits_{\mathop{t}\nolimits_{1}}+\mathop{\xi}\nolimits^{2}\mathop{\partial}% \nolimits_{\mathop{t}\nolimits_{2}},\qquad\nabla=\xi\mathop{\nabla}\nolimits_{% 1},∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = italic_ε ∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT + italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT , ∇ = italic_ξ ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , (101)
S=ξ𝑆(1),T=ξ𝑇(1)+𝜉2𝑇(2),formulae-sequence𝑆𝜉superscript𝑆1𝑇𝜉superscript𝑇1superscript𝜉2superscript𝑇2\displaystyle S=\xi\mathop{S}\nolimits^{(1)},\qquad T=\xi\mathop{T}\nolimits^{% (1)}+\mathop{\xi}\nolimits^{2}\mathop{T}\nolimits^{(2)},italic_S = italic_ξ italic_S start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT , italic_T = italic_ξ italic_T start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_T start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT , (102)

and substituting Eqs. (100)-(102) into Eq. (99) yields

O(𝜉0):𝑙i=𝑙ieq,:𝑂superscript𝜉0subscript𝑙𝑖superscriptsubscript𝑙𝑖𝑒𝑞\displaystyle O(\mathop{\xi}\nolimits^{0}):\mathop{l}\nolimits_{i}=\mathop{l}% \nolimits_{i}^{eq},italic_O ( italic_ξ start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT ) : italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_e italic_q end_POSTSUPERSCRIPT , (103)
O(𝜉1):𝐷1i𝑙i(0)=1Δt𝜏l𝑙i(1)+𝑆i(1)+𝑇i(1),:𝑂superscript𝜉1subscript𝐷1𝑖superscriptsubscript𝑙𝑖01Δ𝑡subscript𝜏𝑙superscriptsubscript𝑙𝑖1superscriptsubscript𝑆𝑖1superscriptsubscript𝑇𝑖1\displaystyle O(\mathop{\xi}\nolimits^{1}):\mathop{D}\nolimits_{1i}\mathop{l}% \nolimits_{i}^{(0)}=-\frac{1}{{\Delta t\mathop{\tau}\nolimits_{l}}}\mathop{l}% \nolimits_{i}^{(1)}+\mathop{S}\nolimits_{i}^{(1)}+\mathop{T}\nolimits_{i}^{(1)},italic_O ( italic_ξ start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT ) : italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT = - divide start_ARG 1 end_ARG start_ARG roman_Δ italic_t italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT , (104)
O(𝜉2):𝑡2𝑙i(0)+𝐷1i𝑙i(1)+Δt2𝐷1i2𝑙i(0)=1Δt𝜏l𝑙i(2)+𝑇i(2).:𝑂superscript𝜉2subscriptsubscript𝑡2superscriptsubscript𝑙𝑖0subscript𝐷1𝑖superscriptsubscript𝑙𝑖1Δ𝑡2superscriptsubscript𝐷1𝑖2superscriptsubscript𝑙𝑖01Δ𝑡subscript𝜏𝑙superscriptsubscript𝑙𝑖2superscriptsubscript𝑇𝑖2\displaystyle O(\mathop{\xi}\nolimits^{2}):\mathop{\partial}\nolimits_{\mathop% {t}\nolimits_{2}}\mathop{l}\nolimits_{i}^{(0)}+\mathop{D}\nolimits_{1i}\mathop% {l}\nolimits_{i}^{(1)}+\frac{{\Delta t}}{2}\mathop{\mathop{D}\nolimits_{1i}}% \nolimits^{2}\mathop{l}\nolimits_{i}^{(0)}=-\frac{1}{{\Delta t\mathop{\tau}% \nolimits_{l}}}\mathop{l}\nolimits_{i}^{(2)}+\mathop{T}\nolimits_{i}^{(2)}.italic_O ( italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) : ∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT + italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + divide start_ARG roman_Δ italic_t end_ARG start_ARG 2 end_ARG start_BIGOP italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT end_BIGOP start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT = - divide start_ARG 1 end_ARG start_ARG roman_Δ italic_t italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT + italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT . (105)

With Eq. (66), Eq. (98) and Eq. (100), one can obtain

i𝑙i(1)=Δt2R,subscript𝑖superscriptsubscript𝑙𝑖1Δ𝑡2𝑅\displaystyle\sum\limits_{i}{\mathop{l}\nolimits_{i}^{(1)}=-\frac{{\Delta t}}{% 2}}R,∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT = - divide start_ARG roman_Δ italic_t end_ARG start_ARG 2 end_ARG italic_R , (106)
i𝑙i(k)=0,k>1.formulae-sequencesubscript𝑖superscriptsubscript𝑙𝑖𝑘0𝑘1\displaystyle\sum\limits_{i}{\mathop{l}\nolimits_{i}^{(k)}=0,\quad k>1}.∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( italic_k ) end_POSTSUPERSCRIPT = 0 , italic_k > 1 . (107)

where the forcing term R𝑅Ritalic_R is given by Eq. (61) . Combining Eq. (104) and Eq. (105), we get

𝑡2𝑙i(0)+(112𝜏l)𝐷1i𝑙i(1)+Δt2𝐷1i(𝑆i(1)+𝑇i(1))=1Δt𝜏l𝑙i(2)+𝑇i(2).subscriptsubscript𝑡2superscriptsubscript𝑙𝑖0112subscript𝜏𝑙subscript𝐷1𝑖superscriptsubscript𝑙𝑖1Δ𝑡2subscript𝐷1𝑖superscriptsubscript𝑆𝑖1superscriptsubscript𝑇𝑖11Δ𝑡subscript𝜏𝑙superscriptsubscript𝑙𝑖2superscriptsubscript𝑇𝑖2\mathop{\partial}\nolimits_{\mathop{t}\nolimits_{2}}\mathop{l}\nolimits_{i}^{(% 0)}+(1-\frac{1}{{2\mathop{\tau}\nolimits_{l}}})\mathop{D}\nolimits_{1i}\mathop% {l}\nolimits_{i}^{(1)}+\frac{{\Delta t}}{2}\mathop{D}\nolimits_{1i}(\mathop{S}% \nolimits_{i}^{(1)}+\mathop{T}\nolimits_{i}^{(1)})=-\frac{1}{{\Delta t\mathop{% \tau}\nolimits_{l}}}\mathop{l}\nolimits_{i}^{(2)}+\mathop{T}\nolimits_{i}^{(2)}.∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 0 ) end_POSTSUPERSCRIPT + ( 1 - divide start_ARG 1 end_ARG start_ARG 2 italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG ) italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + divide start_ARG roman_Δ italic_t end_ARG start_ARG 2 end_ARG italic_D start_POSTSUBSCRIPT 1 italic_i end_POSTSUBSCRIPT ( italic_S start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT + italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ) = - divide start_ARG 1 end_ARG start_ARG roman_Δ italic_t italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT end_ARG italic_l start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT + italic_T start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT . (108)

Then, summing Eq. (104) and Eq. (108) over i𝑖iitalic_i, we have

𝑡1q+1(q𝐮)=R,subscriptsubscript𝑡1𝑞subscript1𝑞𝐮𝑅\displaystyle\mathop{\partial}\nolimits_{\mathop{t}\nolimits_{1}}q+\mathop{% \nabla}\nolimits_{1}\cdot(q\mathbf{u})=R,∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_q + ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ⋅ ( italic_q bold_u ) = italic_R , (109)
𝑡2q+1[(12𝜏l)Δt𝑐S21q]=0.subscriptsubscript𝑡2𝑞subscript1delimited-[]12subscript𝜏𝑙Δ𝑡superscriptsubscript𝑐𝑆2subscript1𝑞0\displaystyle\mathop{\partial}\nolimits_{\mathop{t}\nolimits_{2}}q+\mathop{% \nabla}\nolimits_{1}[(\frac{1}{2}-\mathop{\tau}\nolimits_{l})\Delta t\mathop{c% }\nolimits_{S}^{2}\mathop{\nabla}\nolimits_{1}q]=0.∂ start_POSTSUBSCRIPT italic_t start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT italic_q + ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT [ ( divide start_ARG 1 end_ARG start_ARG 2 end_ARG - italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT ) roman_Δ italic_t italic_c start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∇ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_q ] = 0 . (110)

Eventually, taking Eq. (109) ×ξabsent𝜉\times\xi× italic_ξ + Eq. (110) ×𝜉2superscript𝜉2\times\mathop{\xi}\nolimits^{2}× italic_ξ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT supplemented with the imcompressibility condition 𝐮=0𝐮0\nabla\cdot\mathbf{u}=0∇ ⋅ bold_u = 0, the Nernst-Planck equation can be derived with

α=(𝜏l12)𝑐S2Δt.𝛼subscript𝜏𝑙12superscriptsubscript𝑐𝑆2Δ𝑡\alpha=(\mathop{\tau}\nolimits_{l}-\frac{1}{2})\mathop{c}\nolimits_{S}^{2}% \Delta t.italic_α = ( italic_τ start_POSTSUBSCRIPT italic_l end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ) italic_c start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_t . (111)

References

  • [1] D.A. Saville, Electrohydrodynamics: the Taylor-Melcher leaky dielectric model, Annu. Rev. Fluid Mech. 29 (1) (1997) 27-64.
  • [2] P.V. Raje, N.C. Murmu, A review on electrohydrodynamic-inkjet printing technology, Int. J. Emerg. Technol. Adv. Eng 4(5) (2014) 174-183.
  • [3] W. Deng, A. Gomez, Electrospray cooling for microelectronics, Int. J. Heat Mass Transf. 54 (11-12) (2011) 2270-2275.
  • [4] D. Bhattacharya, S. Chakraborty, A. Karmakar, S. Chattopadhyay, Understanding the voltage-induced electrowetting and microfluidic droplet movement phenomena on a Teflon-on-flexible (TOF) substrate, Phys. Fluids 36(3) (2024).
  • [5] C.R. Vigo, W.D. Ristenpart, Aggregation and coalescence of oil droplets in water via electrohydrodynamic flows, Langmuir 26 (13) (2010) 10703-10707.
  • [6] F.C. Lai, EHD gas pumping–A concise review of recent development, J. Electrost. 106 (2020) 103469.
  • [7] D.T. Papageorgiou, Film flows in the presence of electric fields, Annu. Rev. Fluid Mech. 51 (2019) 155-187.
  • [8] L.D. Landau, E.M. Lifshitz, Electrodynamics of Continuous Media, second ed., Pergamon, Oxford, 1975.
  • [9] J.R. Melcher, G.I. Taylor, Electrohydrodynamics: a review of the role of interfacial shear stress, Annu. Rev. Fluid Mech. 1 (1) (1969) 111–146.
  • [10] C.T. O’KONSKI, H.C. Thacher, The distortion of aerosol droplets by an electric field. J. Phys. Chem. 57 (9) (1953) 955-958.
  • [11] O.A. Basaran, L.E. Scriven, Axisymmetric shapes and stability of charged drops in an external electric field. Phys. Fluids A 1 (5) (1989) 799-809.
  • [12] J.Q. Feng, T.C. Scott, A computational analysis of electrohydrodynamics of a leaky dielectric drop in an electric field, J. Fluid Mech. 311 (1996) 289-326.
  • [13] G.I. Taylor, Disintegration of water drops in an electric field. Proc. R. Soc. London, Ser. A 280 (1382) (1964) 383-397.
  • [14] R.V. Craste, O.K. Matar. Electrically induced pattern formation in thin leaky dielectric films. Phys. Fluids 17 (3) (2005).
  • [15] Y. Mori, Y.N. Young, From electrodiffusion theory to the electrohydrodynamics of leaky dielectrics through the weak electrolyte limit, J. Fluid Mech. 855 (2018) 67-130.
  • [16] G. Supeene, C.R. Koch, S. Bhattacharjee. Deformation of a droplet in an electric field: Nonlinear transient response in perfect and leaky dielectric media, J. Colloid Interface Sci. 318 (2) (2008) 463-476.
  • [17] J. Zhang, D.Y. Kwok, A 2D lattice Boltzmann study on electrohydrodynamic drop deformation with the leaky dielectric theory, J. Comput. Phys. 206 (1) (2005) 150-161.
  • [18] W. Zhang, J. Wang, W. Hu, H. Liu, B. Li, K. Yu, Enhancement of electric field on bubble dispersion characteristics in leaky-dielectric liquid medium, Int. J. Multiph. Flow 142 (2021) 103743.
  • [19] J. Zhao, Z. Wang, Y. Gu, E. Sauret, Electrohydrodynamic viscous fingering of leaky dielectric fluids in a channel, Phys. Fluids 35 (3) (2023).
  • [20] J.Q. Feng, Electrohydrodynamic behaviour of a drop subjected to a steady uniform electric field at finite electric Reynolds number, Proc. R. Soc. London, Ser. A 455 (1986) (1999) 2245-2269.
  • [21] E. Yariv, Y. Almog, The effect of surface-charge convection on the settling velocity of spherical drops in a uniform electric field, J. Fluid Mech. 797 (2016) 536-548.
  • [22] R. Sengupta, L.M. Walker, A.S. Khair, The role of surface charge convection in the electrohydrodynamics and breakup of prolate drops, J. Fluid Mech. 833 (2017) 29-53.
  • [23] S. Mandal, A. Bandopadhyay, S. Chakraborty, Effect of surface charge convection and shape deformation on the dielectrophoretic motion of a liquid drop, Phys. Rev. E 93 (4) (2016) 043127.
  • [24] K. Luo, J. Wu, H.L. Yi, H.P. Tan, Numerical analysis of two-phase electrohydrodynamic flows in the presence of surface charge convection. Phys. Fluids 32 (12) (2020).
  • [25] R.A. Hayes, B.J. Feenstra, Video-speed electronic paper based on electrowetting, Nature 425 (6956) (2003) 383-385.
  • [26] H.B. Eral, D.M. Augustine, M.H. Duits, F. Mugele, Suppressing the coffee stain effect: how to control colloidal self-assembly in evaporating drops using electrowetting, Soft Matter 7 (10) (2011) 4954–4958.
  • [27] P. Paik, V.K. Pamula, M.G. Pollack, R.B. Fair, Electrowetting-based droplet mixers for microfluidic systems, Lab Chip 3(1) (2003) 28-33.
  • [28] K.L. Wang, T.B. Jones, Electrowetting dynamics of microfluidic actuation, Langmuir 21 (9) (2005) 4211-4217.
  • [29] G. Lippmann, Relations entre les phénomènes électriques et capillaires, Ann. Chim. Phys. 5 (11) (1875) 494-549.
  • [30] G. McHale, C.V. Brown, M.I. Newton, G.G. Wells, N. Sampara, Dielectrowetting driven spreading of droplets, Phys. Rev. Lett. 107 (18) (2011) 186101.
  • [31] J.J. Huang, C. Shu, J.J. Feng, Y.H. Chew, A phase-field-based hybrid lattice-Boltzmann finite-volume method and its application to simulate droplet motion under electrowetting control. J. Adhes. Sci. Technol. 26 (12-17) (2012) 1825-1851.
  • [32] K. Mohseni, A. Arzpeyma, A. Dolatabadi. Behaviour of a moving droplet under electrowetting actuation: Numerical simulation, Can. J. Chem. Eng. 84 (1) (2006) 17-21.
  • [33] X. Huo, L. Li, Y. Yang, X. Liu, Q. Yu, Q. Wang, The dynamics of directional transport of a droplet in programmable electrowetting channel, Phys. Fluids 35 (3) (2023).
  • [34] Q. Zhao, W. Ren. A finite element method for electrowetting on dielectric, J. Comput. Phys. 429 (2021) 109998.
  • [35] J.S. Hong, S.H. Ko, K.H. Kang, I.S. Kang, A numerical investigation on AC electrowetting of a droplet, Microfluid. Nanofluid. 5 (2008) 263-271.
  • [36] D.S. Pillai, K.C. Sahu, R. Narayanan, Electrowetting of a leaky dielectric droplet under a time-periodic electric field, Phys. Rev. Fluids 6 (7) (2021) 073701.
  • [37] H. Aminfar, M. Mohammadpourfard, Lattice Boltzmann method for electrowetting modeling and simulation, Comput. Meth. Appl. Mech. Eng. 198 (47-48) (2009) 3852-3868.
  • [38] S. Goel, D.S. Pillai, Electrokinetic Thin-Film Model for Electrowetting: The Role of Bulk Charges, Langmuir, 39 (37) (2023) 13076-13089.
  • [39] X. Xu, F. Wang, Z. Qin, B. Wen, Electrowetting lattice Boltzmann method for micro-and nano-droplet manipulations, Phys. Rev. E 107 (4) (2023) 045305.
  • [40] L.M. Fu, J.Y. Lin, R.J. Yang, Analysis of electroosmotic flow with step change in zeta potential, J. Colloid Interface Sci. 258 (2) (2003) 266–275.
  • [41] H.M. Park, J.S. Lee, T.W. Kim, Comparison of the Nernst-Planck model and the Poisson-Boltzmann model for electroosmotic flows in microchannels, J. Colloid Interface Sci. 315 (2) (2007) 731–739.
  • [42] S. Bhattacharyya, A.K. Nayak, Electroosmotic flow in micro/nanochannels with surface potential heterogeneity: An analysis through the Nernst-Planck model with convection effect, Colloid Surf. A-Physicochem. Eng. Asp. 339 (1–3) (2009) 167–177.
  • [43] C.K. Aidun, J.R. Clausen, Lattice Boltzmann method for complex flows, Annu. Rev. Fluid Mech. 42 (2010) 439–472.
  • [44] T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E.M. Viggen, The lattice Boltzmann method, Springer International Publishing 10 (978-3) (2017) 4-15.
  • [45] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (3) (1993) 1815.
  • [46] L. Wang, J. Huang, K. He. Thermal lattice Boltzmann model for liquid-vapor phase change, Phys. Rev. E 106 (5) (2022) 055308.
  • [47] X. He, S. Chen, R. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability, J. Comput. Phys. 152 (2) (1999) 642-663.
  • [48] H. Liang, B.C. Shi, Z.L. Guo, Z. Chai, Phase-field-based multiple-relaxation-time lattice Boltzmann model for incompressible multiphase flows, Phys. Rev. E 89 (5) (2014) 053320.
  • [49] X. Yang, B. Shi, Z. Chai, Z. Guo, A coupled lattice Boltzmann method to solve Nernst–Planck model for simulating electro-osmotic flows, J. Sci. Comput. 61 (2014) 222-238.
  • [50] R. Singh, S.S. Bahga, A. Gupta, Electrohydrodynamics in leaky dielectric fluids using lattice Boltzmann method, Eur. J. Mech. B-Fluids 74 (2019) 167-179.
  • [51] X. Liu, Z. Chai, B. Shi. A phase-field-based lattice Boltzmann modeling of two-phase electro-hydrodynamic flows, Phys. Fluids 31 (9) (2019).
  • [52] T. Qian, X.P. Wang, P. Sheng, A variational approach to moving contact line hydrodynamics, J. Fluid Mech. 564 (2006) 333-360.
  • [53] X. Xu, U. Thiele, T. Qian, A variational approach to thin film hydrodynamics of binary mixtures. J. Phys.-Condes. Matter 27 (8) (2015) 085005.
  • [54] X. Xu, T. Qian, Hydrodynamic boundary conditions derived from Onsager’s variational principle, Procedia IUTAM 20 (2017) 144-151.
  • [55] A. Castellanos, Ed. Electrohydrodynamics. Vol. 380. Springer, 2014.
  • [56] L. Onsager, Reciprocal relations in irreversible processes. I., Phys. Rev. 37 (1931) 405.
  • [57] S.J. Blundell, K.M. Blundell, Concepts in thermal physics, Oup Oxford 2010.
  • [58] C. Eck, M. Fontelos, G. Grün, F. Klingbeil, On a phase-field model for electrowetting, Interface Free Bound. 11 (2) (2009) 259-290.
  • [59] J. Shen, Modeling and numerical approximation of two-phase incompressible flows by a phase-field approach, in Multiscale Modeling and Analysis for Materials Simulation, Lect. NotesSer. Inst. Math. Sci. Natl. Univ. Singap. 22, World Scientific, Hackensack, NJ 2012 147-195.
  • [60] S. Dong, Multiphase flows of N immiscible incompressible fluids: a reduction-consistent and thermodynamically-consistent formulation and associated algorithm, J. Comput. Phys. 361 (2018) 1-49.
  • [61] J. Yang, I.C. Christov, S. Dong, Phase Field Modeling and Numerical Algorithm for Two-Phase Dielectric Fluid Flows. arXiv preprint arXiv, (2023) 2305.09020.
  • [62] L. Ju, P. Liu, B. Yan, J. Bao, S. Sun, and Z. Guo, A well-balanced lattice Boltzmann model for binary fluids based on the incompressible phase-field theory, arXiv: 2311.10827 [Commun. Comput. Phys. (to be published)].
  • [63] S.P. Thampi, S. Ansumali, R. Adhikari, S. Succi, Isotropic discrete Laplacian operators from lattice hydrodynamics. J. Comput. Phys. 234 (2013) 1-7.
  • [64] Z. Chai, B. Shi, A novel lattice Boltzmann model for the Poisson equation, Appl. Math. Model. 32 (10) (2008) 2050-2058.
  • [65] G. Tomar, D. Gerlach, G. Biswas, N. Alleborn, A. Sharma, F. Durst, S.W.J. Welach, A. Delgado, Two-phase electrohydrodynamic simulations using a volume-of-fluid approach, J. Comput. Phys. 227 (2) (2017) 1267-1285.
  • [66] Z. Chai, B. Shi, Z. Guo, A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection–diffusion equations, J. Sci. Comput. 69 (2016) 355-390.
  • [67] M. Zagnoni, C.N. Baroud, J.M. Cooper. Electrically initiated upstream coalescence cascade of droplets in a microfluidic flow, Phys. Rev. E 80 (4) (2009) 046303.
  • [68] J. Perko, R.A. Patel, Single-relaxation-time lattice Boltzmann scheme for advection-diffusion problems with large diffusion-coefficient heterogeneities and high-advection transport, Phys. Rev. E 89 (5) (2014) 053309.
  • [69] X. Xiang, Z. Wang, B. Shi, Modified lattice Boltzmann scheme for nonlinear convection diffusion equations. Commun. Nonlinear Sci. Numer. Simul. 17 (6) (2012) 2415-2425.
  • [70] L. Wang, Z. Wei, T. Li, Z. Chai, B. Shi, A lattice Boltzmann modelling of electrohydrodynamic conduction phenomenon in dielectric liquids, Appl. Math. Model. 95 (2021) 361-378.
  • [71] Z. Chai, T.S. Zhao, Nonequilibrium scheme for computing the flux of the convection-diffusion equation in the framework of the lattice Boltzmann method, Phys. Rev. E 90 (1) (2014) 013305.
  • [72] V. Khatavkar, P. Anderson, H. Meijer, Capillary spreading of a droplet in the partially wetting regime using a diffuse-interface model, J. Fluid Mech. 572 (2007) 367-387.
  • [73] H. Liang, H. Liu, Z. Chai, B. Shi, Lattice Boltzmann method for contact-line motion of binary fluids with high density ratio, Phys. Rev. E 99 (6) (2019) 063306.
  • [74] J.M. López-Herrera, S. Popinet, M.A. Herrada, A charge-conservative approach for simulating electrohydrodynamic two-phase flows using volume-of-fluid, Journal of Computational Physics 230(5) (2011) 1939-1955.
  • [75] Y. Zhao, S.K. Chung, U.C. Yi, S.K. Cho, Droplet manipulation and microparticle sampling on perforated microfilter membranes, J. Micromech. Microeng. 18 (2) (2008) 025030.
  • [76] S.J. Lee, J. Hong, K.H. Kang, I.S. Kang, Electrowetting-induced droplet detachment from hydrophobic surfaces, Langmuir 30 (7) (2014) 1805-1811.
  • [77] Q. Vo, T. Tran, Critical conditions for jumping droplets, Phys. Rev. Lett. 123 (2) (2019) 024502.
  • [78] Q. Wang, M. Xu., C. Wang, J. Gu, N. Hu, J. Lyu, W. Yao, Actuation of a nonconductive droplet in an aqueous fluid by reversed electrowetting effect, Langmuir 36 (28) (2020) 8152-8164.
  • [79] W. Weng, Q. Wang, J. Gu, J. Li, C. Wang, W. Yao, The dynamics of droplet detachment in reversed electrowetting (REW), Colloid Surf. A-Physicochem. Eng. Asp. 616 (2021) 126303.