A Finite Element Method for Electrowetting on Dielectric
aa r X i v : . [ phy s i c s . c o m p - ph ] J un A FINITE ELEMENT METHOD FOR ELECTROWETTING ONDIELECTRIC
QUAN ZHAO ∗ AND
WEIQING REN † Abstract.
We consider the problem of electrowetting on dielectric (EWoD). The system involvesthe dynamics of a conducting droplet, which is immersed in another dielectric fluid, on a dielectricsubstrate under an applied voltage. The fluid dynamics is modeled by the two-phase incompressibleNavier-Stokes equations with the standard interface conditions, the Navier slip condition on thesubstrate and a contact angle condition which relates the dynamic contact angle and the contactline velocity, as well as the kinematic condition for the evolution of the interface. The electric forceacting on the fluid interface is modeled by the Maxwell’s equations in the domain occupied by thedielectric fluid and the dielectric substrate. We develop a numerical method for the model based onits weak form. This method combines the finite element method for the Navier-Stokes equations on afixed bulk mesh with a parametric finite element method for the dynamics of the fluid interface, andthe boundary integral method for the electric force along the fluid interface. Numerical examplesare presented to demonstrate the accuracy and convergence of the numerical method, the effect ofvarious physical parameters on the interface profile and other interesting phenomena such as thetransportation of droplet driven by applied non-uniform electric potential difference.
Key words.
Electrowetting, moving contact lines, two-phase flows, the finite element method
AMS subject classifications.
1. Introduction.
Since the pioneer work of Lippmann [23] on electro-capillarity,it has been found that applied electric fields have a great effect on the wetting behaviorof small charged droplets. This phenomenon is referred to as electrowetting and hasreceived much attention in recent years [27, 44, 7]. In the device of electro-wetting ondielectric (EWoD), a dielectric film is placed on the substrate to separate the dropletand the electrode to avoid electrolytic decomposition [4] (see Fig. 1.1 for the set-upof EWoD). EWoD has found many applications in various fields, such as adjustablelenses [20], electronic displays [18], lab-on-a-chip devices [31, 8], suppressing coffeestrain effects [15], etc.The static problem of EWoD has been extensively studied in recent years, forexample, in Refs. [19, 6, 38, 5, 28, 25, 37, 16, 11, 12, 13] and many others. Thesework has revealed the structure of the static interface profile. It was found that theelectric force does not contribute to the force balance at the contact line, thereforethe local static contact angle θ Y still satisfies the Young-Dupr´e equation(1.1) γ cos θ Y = γ − γ , where γ , γ and γ are the surface tension coefficients of the fluid-fluid and fluid-solidinterfaces. On the other hand, the divergent electric force incurs a large curvatureand causes a significant deformation of the fluid interface in a small neighborhood ofthe contact line. The contact angle of the interface outside this small region, calledthe apparent contact angle and denoted by θ B , is well characterized by the Lippmannequation [32, 6, 13](1.2) cos θ B = cos θ Y + ǫφ γ d , ∗ Department of Mathematics, National University of Singapore, Singapore 119076([email protected]). † Corresponding author. Department of Mathematics, National University of Singapore, Singa-pore, 119076 ([email protected]). 1
Q. Zhao and W. Ren where φ is the applied voltage, ǫ and d are the permittivity and thickness of thedielectric substrate, respectively. Except in the extreme case of contact angle satura-tion, this equation also matches experimental results quite well for different types ofdroplets and insulators, and wide range of φ and d (see [39, 40, 26] for example).In this work, we consider the dynamical problem of EWoD. Because of its impor-tance in industrial applications, a lot of efforts have been devoted to this problem andsome numerical methods have been proposed in recent years. These include LatticeBoltzmann methods [9, 22, 36], molecular dynamics simulations [14, 21], the level setmethod [41, 17], the phase-field approach [24, 16, 30], and others [8, 29, 10], etc. Here,we develop a finite element method for EWoD based on our earlier work on movingcontact lines [43].The model we use in this work for EWoD is based on the contact line modeldeveloped by Ren et al. [33, 35, 34]. It contains the incompressible Navier-Stokesequations for the two-phase fluid dynamics, the Navier slip condition on the substrateand a dynamic contact angle condition at the contact line. On the fluid interface,besides the viscous stress and the curvature force, the electric force also contributesto the force balance thus the interface conditions. We assume that the electric chargingtime is negligible compared to the time scale of the fluid motion, therefore model theelectrostatic potential using the Maxwell’s equation [13].Based on the previous work for two-phase fluid dynamics [3] and moving contactlines [43], we develop a finite element method for the EWoD model. The methodcouples the finite element method for the incompressible Navier-Stokes equations anda semi-implicit parametric finite element method for the evolution of the fluid inter-face. We use unfitted mesh such that the discretization of the moving fluid interface isdecoupled from the fixed bulk mesh. Besides, the electric force on the fluid interfaceis computed by using the boundary integral method. The numerical method obeys asimilar energy law as the continuum model when the electric field is absent, which isa desired property for the numerical method.This paper is organzied as follows. In section 2, we present the EWoD model andits dimensionless form. In section 3, we develop the numerical method based on aweak form of the continuum model, and present the full discretized scheme for theNavier-Stokes equations, the dynamics of the fluid interface, and the electric force onthe fluid interface. In section 4, we present numerical examples to demonstrate theconvergence and accuracy of the numerical method, the effect of the various physicalparameters on the interface profiles, as well as the wetting dynamics driven by non-uniform electric potential. Finally, we draw the conclusion in section 5.
2. The mathematical model.
We consider the EWoD problem in the 2d spaceas shown in Fig. 1.1(a). In this setup, a conductive liquid droplet (shaded in blue),is immersed in another dielectric fluid such as vapor or oil, and deposited on a thindielectric substrate/insulator (shaded in orange), below which we place an electrode(shaded in gray). Between the electrode and the droplet, we apply a voltage difference,which generates electric fields and influences the wetting behavior of the chargeddroplet.The corresponding mathematical setup is shown in Fig. 1.1(b). We consider thesystem in a bounded domain and use the Cartesian coordinates, where the fluid-insulator interface is on the x -axis. For simplicity, we assume the periodic structurealong x -direction for the system, and moreover, another electrode is placed on the topΓ . The domain is composed of three regions, the conducting droplet, the dielectricfluid, and the dielectric substrate, which are labeled as Ω , Ω , and Ω , respectively. Finite Element Method for Electrowetting on Dielectric Dielectricfluid Voltage (a) (b)
ElectrodeConductive fluidDielectric substrate
Fig. 1.1 . (a) An illustration of the EWoD, where a conductive liquid droplet (shaded in blue) isdeposited on a dielectric-coated electrode (shaded in gray). (b) Geometric setup of the EWoD systemwith moving contact lines in two-phase flows in a bounded domain where Ω ∪ Ω = [ − L x , L x ] × [0 , L y ] and Ω = [ − L x , L x ] × [ − d, . The fluid interface is denoted by the open curve Γ, which intersects with the flat insula-tor at two contact points, labeled as x l and x r . Γ and Γ represent the correspondingfluid-insulator interfaces, and Γ : y = − d is the insulator-electrode interface.Some relevant parameters are defined as follows: ρ i and µ i are the density andviscosity of fluid i ( i = 1 , ǫ i is permittivity of the dielectric medium inΩ i ( i = 2 , γ denotes the surface tension of the fluid interface Γ, and γ i denotes thesurface tension of the fluid-solid interface Γ i ( i = 1 , n the unit normal vector on Γ pointing to fluid 2, and n w = (0 , − T and t w = (1 , T the unit normal and tangent vector on Γ ∪ Γ respectively. The mathematical modelis an extension of the contact line model proposed by Ren et al. [33, 35, 34] to includethe contribution of electric field in EWoD. It consists of the two-phase incompressibleNavier-Stokes equations for the fluid dynamics and the Maxwell’s equation for theelectrostatic potential.First we consider the two-phase fluid dynamics in the domain Ω = Ω ( t ) ∪ Ω ( t ).Let u ( x , t ) : Ω × [0 , T ] → R be the fluid velocity and p ( x , t ) : Ω × [0 , T ] → R bethe pressure. The dynamic of the system is governed by the standard incompressibleNavier-Stokes equations in Ω i ( t )( i = 1 , ρ i ( ∂ t u + u · ∇ u ) = −∇ p + ∇ · τ d , (2.1a) ∇ · u = 0 , (2.1b)where τ d = 2 µ i D ( u ) is the viscous stress with D ( u ) = ( ∇ u + ( ∇ u ) T ).On the fluid interface Γ( t ), we have the following conditions[ u ] = , (2.2a) [ p I − τ d ] · n = (cid:16) γκ + ǫ |∇ Φ | (cid:17) n , (2.2b) v n = u | Γ( t ) · n , (2.2c)where [ · ] denotes the jump from fluid 1 to fluid 2, I ∈ R × is the identity matrix, κ is the curvature of the fluid interface Γ, Φ is the electrostatic potential, and v n denotes the normal velocity of the fluid interface. Equation (2.2a) states the fluidvelocity is continuous across the interface. Equation (2.2b) states that the tangential Q. Zhao and W. Ren stress is continuous across the interface and the normal stress jump is balanced bythe curvature force γκ n together with the electric force ǫ |∇ Φ | n . Equation (2.2c) isthe kinetic condition for the interface, where the fluid interface evolves according tothe local fluid velocity.At the lower solid wall Γ i ( t ) ( i = 1 , u · n w = 0 , (2.3a) t w · τ d · n w = − β i u s , (2.3b)where β i ( i = 1 ,
2) is the friction coefficient of fluid i at the wall, and u s = u · t w isthe slip velocity of the fluids.At the contact points x l and x r , the dynamic contact angles, denoted by θ ld and θ rd respectively, depend on the contact line velocity [33, 35], γ (cos θ ld − cos θ Y ) = β ∗ ˙ x l , (2.4a) γ (cos θ rd − cos θ Y ) = − β ∗ ˙ x r , (2.4b)where θ Y is the equilibrium contact angle satisfying the Young’s relation (1.1), and β ∗ is the friction coefficient of the fluid interface at the contact line. The contactline velocity is determined by the slip velocity of the fluids: ˙ x l,r = u s | x = x l,r . Thecontact angle condition (2.4) is a force balance at the contact point: the Young stressresulted from the deviation of the dynamic contact angle from the equilibrium angleis balanced by the friction force. Note that the electric force does not contribute tothe force balance at the contact line.Furthermore, we use the no-slip boundary condition on the upper wall Γ and theperiodic boundary conditions along Γ . The applied voltageinduces electrostatic fields in the fluid region Ω and the dielectric substrate Ω , wherethe electric potential, denoted by Φ, satisfies the Laplace equation,(2.5) ∇ Φ = ∂ Φ ∂x + ∂ Φ ∂y = 0 , x ∈ Ω , Ω . On the interface Γ ∪ Γ and at the top/bottom electrode Γ ∪ Γ , the electric potentialsatisfies the Dirichlet boundary conditionsΦ = φ, on Γ ∪ Γ , (2.6a) Φ = 0 , on Γ ∪ Γ , (2.6b)where φ > n w · [ ǫ ∇ Φ] = 0 , where [ · ] denotes the jump from Ω to Ω , ǫ = ǫ in Ω and ǫ = ǫ in Ω . Furthermore,the periodic boundary condition is prescribed along Γ ∪ Γ . Next we present the above governing equationsand the boundary/interface conditions in their dimensionless form. By choosing L Finite Element Method for Electrowetting on Dielectric U as the characteristic length and velocity respectively, we rescale the physicalquantities asˆ ρ i = ρ i ρ , ˆ µ i = µ i µ , ˆ β i = β i β , ˆ β ∗ = β ∗ µ , ˆ γ i = γ i γ , ˆ ǫ i = ǫ i ǫ , ˆ x = x L , ˆ t = U tL , ˆ u = u U , ˆ p = pρ U , ˆ κ = Lκ, ˆΦ = Φ φ .
We define the Reynolds number Re , the Capillary number Ca , the slip length l s , theWeber number W e and the parameter η as(2.8) Re = ρ U Lµ , Ca = µ Uγ , l s = µ β L , W e = Re · Ca, η = ǫ φ γL , where η measures the relative strength of the electric force compared to the surfacetension at the fluid interface.For ease of presentation, we will drop the hats on the dimensionless parame-ters and variables. In the dimensionless form, the governing equations for the fluiddynamics in Ω i ( i = 1 ,
2) read ρ i ( ∂ t u + u · ∇ u ) + ∇ · T = 0 , (2.9a) ∇ · u = 0 , (2.9b)where T = p I − Re τ d is the stress tensor. These equations are supplemented withthe following boundary/interface conditions:(i) The interface conditions on Γ( t ): (cid:2) u (cid:3) = , (2.10a) W e (cid:2) T (cid:3) · n = (cid:16) κ + ǫ η |∇ Φ | (cid:17) n , (2.10b) κ = ∂ ss X · n , (2.10c) v n = u | Γ( t ) · n , (2.10d)where X denotes the fluid interface, s is the arc length parameter with s = 0being at the left contact point.(ii) The condition for the dynamic contact angles:1 Ca (cos θ ld − cos θ Y ) = β ∗ ˙ x l ( t ) , (2.11a) 1 Ca (cos θ rd − cos θ Y ) = − β ∗ ˙ x r ( t ) . (2.11b)(iii) The boundary conditions on Γ ( t ) ∪ Γ ( t ):(2.12) u · n w = 0 , l s t w · τ d · n w = − β i u s . (iv) The no-slip condition on Γ (2.13) u = . (v) Periodic boundary conditions on Γ :(2.14) u | x = − L x = u | x = L x , T | x = − L x = T | x = L x . Q. Zhao and W. Ren
The governing equations for the electric field potential Φ( x, t ) in Ω i ( i = 2 ,
3) read(2.15) ∇ Φ = 0 , subject to the following boundary conditions:(i) The Dirichlet Boundary conditions on Γ, Γ , Γ and Γ :Φ = 1 , on Γ( t ) ∪ Γ ( t ) , (2.16a) Φ = 0 , on Γ , Γ . (2.16b)(ii) The interface condition on Γ ( t ):(2.17) [Φ] = 0 , n w · [ ǫ ∇ Φ] = 0 . (iii) Periodic boundary conditions on Γ ∪ Γ :(2.18) Φ | x = − L x = Φ | x = L x , ∂ Φ ∂ n (cid:12)(cid:12)(cid:12)(cid:12) x = − L x = − ∂ Φ ∂ n (cid:12)(cid:12)(cid:12)(cid:12) x = L x , where n is the unit outward normal to Ω on Γ .Equations (2.9) and (2.15) together with the boundary/interface conditions (2.10)-(2.14) and (2.16)-(2.18) form the complete model for the problem of EWoD. For thisdynamical system, we define the total energy as W ( t ) = X i =1 , Z Ω i ( t ) ρ i | u | d L − cos θ Y W e | Γ ( t ) | + 1 W e | Γ( t ) |− ηW e X i =2 , Z Ω i ǫ i |∇ Φ | d L , (2.19)where | Γ ( t ) | and | Γ( t ) | are the total arc length of Γ ( t ) and Γ( t ), respectively. Thefour terms represents the kinetic energy of the fluids, the interfacial energy at the solidwall, the interfacial energy of the fluid interface and the electrical energy, respectively.The dynamical system obeys the energy dissipation lawdd t W ( t ) = − X i =1 , Re Z Ω i µ i k D ( u ) k F d L − X i =1 , Re · l s Z Γ i β i | u s | d s − β ∗ Re ( ˙ x l + ˙ x r ) , (2.20)where the three terms represent the viscous dissipation in the bulk fluid with k · k being the Frobenius norm, the dissipation at the solid wall due to the friction and thedissipation at the contact points, respectively. Details for the derivation of (2.20) areprovided in the Appendix A.
3. The numerical method.
The numerical method consists of a finite elementmethod for the fluid dynamics, a parametric finite element method for the dynamicsof the fluid interface, and the boundary integral method for the electric field. Thenumerical method is based on a weak formulation for the EWoD model, which we willpresent first.
Finite Element Method for Electrowetting on Dielectric To propose the weak formulation for equations (2.9)-(2.14), we define the function space for the pressure and the velocity, respectively as P := (cid:26) ζ ∈ L (Ω) : Z Ω ζ d L = 0 (cid:27) , (3.1a) U := n ω = ( ω , ω ) T ∈ (cid:0) H (Ω) (cid:1) : ω · n w = 0 on Γ ∪ Γ , ω = on Γ , (3.1b) ω ( − L x , y ) = ω ( L x , y ) , ∀ ≤ y ≤ L y o , with the L -inner product on Ω = Ω ( t ) ∪ Ω ( t )(3.2) ( u, v ) := X i =1 , Z Ω i ( t ) uv d L , ∀ u, v ∈ L (Ω) . We parameterize the fluid interface as X ( α, t ) = ( X ( α, t ) , Y ( α, t )) T , where α ∈ I =[0 , α = 0 , K := L ( I ) = (cid:26) ψ : I → R , and Z I | ψ ( α ) | | ∂ α X | d α < + ∞ (cid:27) , (3.3a) X := n g = ( g , g ) T ∈ ( H ( I )) : g | α =0 , = 0 o , (3.3b)equipped with the L -inner product on I (3.4) (cid:0) u, v (cid:1) Γ := Z I u ( α ) v ( α ) | ∂ α X | d α, ∀ u, v ∈ L ( I ) . Taking the inner product of (2.9a) with ω ∈ U then using the boundary/interfaceconditions in (2.10)-(2.14) and ∇ · u = 0, we obtain [3, 43] (cid:16) ρ [ ∂ t u + ( u · ∇ ) u ] , ω (cid:17) = 12 (cid:20) dd t ( ρ u , ω ) + ( ρ ∂ t u , ω ) (cid:21) + 12 (cid:16) ρ, [( u · ∇ ) u ] · ω − [( u · ∇ ) ω ] · u (cid:17) , (3.5)where ρ = ρ χ Ω1( t ) + ρ χ Ω2( t ) with χ being the characteristic function. With thespecial treatment of the inertia term in (3.5), the numerical scheme enjoys the discretestability for the fluid kinetic energy in the absence of electric field. This will bediscussed later in section 3.4. For the viscous term, integrating by parts then applyingthe boundary/interface conditions yields (cid:16) ∇ · T , ω (cid:17) = − (cid:16) p, ∇ · ω (cid:17) + 2 Re (cid:16) µD ( u ) , D ( ω ) (cid:17) − (cid:16) [ T ] · n , ω (cid:17) Γ + (cid:16) T · n w , ω (cid:17) Γ ∪ Γ = − (cid:16) p, ∇ · ω (cid:17) + 2 Re (cid:16) µD ( u ) , D ( ω ) (cid:17) − W e (cid:16) κ + ǫ η |∇ Φ | , n · ω (cid:17) Γ + 1 Re · l s (cid:16) β u s , ω s (cid:17) Γ ∪ Γ , (3.6) Q. Zhao and W. Ren where µ = µ χ Ω1( t ) + µ χ Ω2( t ) , β = β χ Γ1( t ) + β χ Γ2( t ) , and ω s = ω · t w .Equation (2.10c) for the curvature can be rewritten as κ n = ∂ ss X . Taking theinner product of it with g = ( g , g ) T ∈ X on Γ( t ) then applying integration by partsyields 0 = (cid:16) κ, n · g (cid:17) Γ + (cid:16) ∂ s X , ∂ s g (cid:17) Γ − ( ∂ s X · g ) (cid:12)(cid:12)(cid:12) α =1 α =0 = (cid:16) κ, n · g (cid:17) Γ + (cid:16) ∂ s X , ∂ s g (cid:17) Γ − ( g ∂ s X ) (cid:12)(cid:12)(cid:12) α =1 α =0 = (cid:16) κ, n · g (cid:17) Γ + (cid:16) ∂ s X , ∂ s g (cid:17) Γ + β ∗ Ca [ ˙ x l g (0) + ˙ x r g (1)] − cos θ Y [ g (1) − g (0)] , (3.7)where we have used the fact that ∂ s X | α =0 = cos θ ld , ∂ s X | α =1 = cos θ rd , and thedynamic contact angle conditions in (2.11).Combining these results, we obtain the weak formulation for the dynamic system(2.9)-(2.14) as follows: Given the initial fluid velocity u and the interface X ( α ), findthe fluid velocity u ( · , t ) ∈ U , the pressure p ( · , t ) ∈ P , the fluid interface X ( · , t ) ∈ X ,and the curvature κ ( · , t ) ∈ K such that12 (cid:20) dd t (cid:16) ρ u , ω (cid:17) + (cid:16) ρ ∂ t u , ω (cid:17) + (cid:16) ρ ( u · ∇ ) u , ω (cid:17) − (cid:16) ρ ( u · ∇ ) ω , u (cid:17)(cid:21) + 2 Re (cid:16) µD ( u ) , D ( ω ) (cid:17) − (cid:16) p, ∇ · ω (cid:17) − W e (cid:16) κ n , ω (cid:17) Γ − ǫ ηW e (cid:16) |∇ Φ | n , ω (cid:17) Γ + 1 Re · l s (cid:16) βu s , ω s (cid:17) Γ ∪ Γ = 0 , ∀ ω ∈ U , (3.8a) (cid:16) ∇ · u , ζ (cid:17) = 0 , ∀ ζ ∈ P , (3.8b) (cid:16) n · ∂ t X , ψ (cid:17) Γ − (cid:16) u · n , ψ (cid:17) Γ = 0 , ∀ ψ ∈ K , (3.8c) (cid:16) κ, n · g (cid:17) Γ + (cid:16) ∂ s X , ∂ s g (cid:17) Γ + β ∗ Ca h ˙ x l g (0) + ˙ x r g (1) i − cos θ Y [ g (1) − g (0)] = 0 , ∀ g ∈ X . (3.8d)Eq. (3.8a) is obtained from Eq. (3.5) and Eq. (3.6). Eq. (3.8b) results from theincompressibility condition (2.9b). Eq. (3.8c) is obtained from the kinetic condition(2.10d). Eq. (3.8d) is obtained from (3.7).The above system (3.8a)-(3.8d) is an extension of the weak formulation introducedin Ref. [3] for two-phase flows and Ref. [43] for moving contact lines. In the currentproblem for EWoD, we have the additional term for the electric force in (3.8a). Theelectric force is obtained from solving (2.15)-(2.18). We solve equations (3.8a)-(3.8d) for u and p on the fluid domain Ω and for X and κ on the reference domain I using thefinite element method on unfitted meshes. To that end, we uniformly partition thetime domain as [0 , T ] = ∪ Mm =1 [ t m − , t m ] with t m = mτ , τ = T /M , and the referencedomain as I = ∪ J Γ j =1 I j with I j = [ α j − , α j ], α j = jh and h = 1 /J Γ . We approximate Finite Element Method for Electrowetting on Dielectric K and X by the finite element spaces K h := n ψ ∈ C ( I ) : ψ | I j ∈ P ( I j ) , ∀ j = 1 , , · · · , J Γ o , (3.9a) X h := n g = ( g , g ) T ∈ ( C ( I )) : g | I j ∈ ( P ( I j )) , ∀ j = 1 , , · · · , J Γ ;(3.9b) g | α =0 , = 0 o , where P ( I j ) denotes the space of the polynomials of degree at most 1 over I j . Denoteby Γ m := X m ( · ) ∈ X h the numerical approximation to the fluid interface Γ( t ) at t = t m . Then Γ m (0 ≤ m ≤ M ) are polygonal curves consisting of ordered linesegments. The unit tangential vector t m and normal vector n m of Γ m are constantvectors on each interval I j with possible jumps at the nodes α j , and they can be easilycomputed as(3.10) t mj := t m | I j = h mj | h mj | , n mj := n m | I j = ( t mj ) ⊥ , ≤ j ≤ J Γ where h mj := X m ( α j ) − X m ( α j − ), and ( · ) ⊥ denotes the counter-clockwise rotationby 90 degrees.Let T h = { ¯ o j } Nj =1 be a triangulation of the bounded domain Ω, and S hk := (cid:8) ϕ ∈ C ( ¯Ω) : ϕ | o j ∈ P k ( o j ) , ∀ j = 1 , , · · · , N (cid:9) , (3.11a) S h := { ϕ ∈ L (Ω) : ϕ | o j ∈ P ( o j ) , ∀ j = 1 , , · · · , N } , (3.11b)where k ∈ N + , and P k ( o j ) denotes the space of polynomials of degree k on o j . Let U h and P h denote the finite element spaces for the numerical solutions for the velocityand pressure, respectively. In this work, we choose them as(3.12) (cid:16) U h , P h (cid:17) = (cid:16)(cid:0) S h (cid:1) ∩ U , (cid:0) S h + S h (cid:1) ∩ P (cid:17) , which satisfies the inf-sup stability condition [1, 43](3.13) inf ϕ ∈ P h sup = ω ∈ U h ( ϕ, ∇ · ω ) k ϕ k k ω k ≥ c > , where k·k and k·k denote the L and H -norm on Ω respectively, and c is a constant.In the simulation, the partition of the reference domain I for the interface andthe triangulation T h of the fluid domain Ω are both fixed in time. As a result, thetriangulation of Ω is decoupled from the discretization of the interface Γ m , thus maynot fit the interface, i.e. the line segments comprising of Γ m are in general not theboundaries of the elements in T h . At t = t m , the interface Γ m divides Ω into two sub-domains, Ω m and Ω m . Correspondingly, we split T h into three disjoint subsets: T m being the set of elements in Ω m , T m being the set of elements in Ω m , and T mf beingthe set of elements that intersect with the interface (see Fig. 3.1 for an illustration).The splitting can be easily done by using a recursive algorithm as follows:(1) First form T mf by locating all elements that intersect with Γ m .(2) Locate one element o j ∗ in Ω m , for example, the one that contains the pointwith coordinate ( ( x ml + x mr ) , T m = { o j ∗ } .(3) Check all neighbours of the elements in T m . If a neighbor is not in T mf , thenadd it to T m .0 Q. Zhao and W. Ren -1 -0.5 0 0.5 100.20.40.60.81 (a) (b)
Fig. 3.1 . (a) Illustration of the discretization of the fluid domain Ω and the fluid interface Γ m (red curve) at t = t m . The fluid interface divides Ω into two sub-domains: Ω m being the regionenclosed by Γ m and the lower wall, Ω m being the region outside. (b) Intersection of the interface Γ m with the bulk mesh. The elements of intersection form T mf . The last step is repeated until no element can be added to T m . This gives the set T m . The rest of the elements not belonging to T m ∪ T mf form the set T m .Denote by ρ m , µ m and β m the numerical approximations of the density ρ ( · , t ), theviscosity µ ( · , t ) and the friction coefficient β ( · , t ) at t = t m , respectively. We define ρ m and µ m as ρ m | o j := ρ , if o j ∈ T m ,ρ , if o j ∈ T m , ( ρ + ρ ) , if o j ∈ T mf , µ m | o j := µ , if o j ∈ T m ,µ , if o j ∈ T m , ( µ + µ ) , if o j ∈ T mf , where 1 ≤ j ≤ N , 0 ≤ m ≤ M . Similarly, we denote by Γ m and Γ m the boundary ofΩ m and Ω m at the lower wall respectively, and define β m at the lower wall as(3.14) β m | ∂o j := β , if ∂o j ⊂ Γ m ,β , if ∂o j ⊂ Γ m , ( β + β ) , if x ml ∈ ∂o j or x mr ∈ ∂o j , where ∂o j is the boundary of o j at the lower wall, and x ml = X m | α =0 and x mr = X m | α =1 are the two contact line points.The finite element method is given as follows. First we partition the time domain[0 , T ], the reference domain I for the interface and the fluid domain Ω as describedabove. Let Γ := X ( · ) ∈ X h and u ∈ U h be the initial interface and velocity field,respectively. For m ≥
0, we compute u m +1 ∈ U h , p m +1 ∈ P h , X m +1 ∈ X h , and Finite Element Method for Electrowetting on Dielectric κ m +1 ∈ K h by solving the linear system12 h(cid:16) ρ m u m +1 − ρ m − u m τ , ω h (cid:17) + (cid:16) ρ m − u m +1 − u m τ , ω h (cid:17) + (cid:16) ρ m ( u m · ∇ ) u m +1 , ω h (cid:17) − (cid:16) ρ m ( u m · ∇ ) ω h , u m +1 (cid:17)i − (cid:16) p m +1 , ∇ · ω h (cid:17) + 2 Re (cid:16) µ m D ( u m +1 ) , D ( ω h ) (cid:17) − W e (cid:16) κ m +1 n m , ω h (cid:17) Γ m − ǫ ηW e (cid:16) |∇ Φ | n m , ω h (cid:17) Γ m + 1 Re · l s (cid:16) β m u m +1 s , ω hs (cid:17) Γ m ∪ Γ m = 0 , ∀ ω h ∈ U h , (3.15a) (cid:16) ∇ · u m +1 , ζ h (cid:17) = 0 , ∀ ζ h ∈ P h , (3.15b) (cid:16) X m +1 − X m τ · n m , ψ h (cid:17) h Γ m − (cid:16) u m +1 · n m , ψ h (cid:17) Γ m = 0 , ∀ ψ h ∈ K h , (3.15c) (cid:16) κ m +1 n m , g h (cid:17) h Γ m + (cid:16) ∂ s X m +1 , ∂ s g h (cid:17) h Γ m − cos θ Y (cid:2) g h (1) − g h (0) (cid:3) + β ∗ Caτ (cid:2)(cid:0) x m +1 r − x mr (cid:1) g h (1) + (cid:0) x m +1 l − x ml (cid:1) g h (0) (cid:3) = 0 , ∀ g h ∈ X h . (3.15d)Here, g h = ( g h , g h ) T , ω hs = ω h · t w , u m +1 s = u m +1 · t w , and x ml = X m | α =0 , x mr = X m | α =1 . The electric force ǫ η |∇ Φ | in (3.15a) is computed by solving equations(2.15)-(2.18) in the domain Ω m ∪ Ω ; its computation will be presented in section 3.3.In (3.15d), the derivative ∂ s X m +1 is taken with respect to the arc length of Γ m : ∂ s X m +1 = | ∂ α X m | ∂ α X m +1 , and similarly for ∂ s g h . At the first time step, we set ρ − = ρ .In the numerical method, the inner product ( · , · ) Γ( t m ) is approximated by using ei-ther the trapezoidal rule, denoted by ( · , · ) h Γ m , or the Simpson’s rule, denoted by ( · , · ) Γ m .Since we are using unfitted mesh, the interface Γ m intersects with the bulk mesh. De-note by (cid:8) α ′ j (cid:9) J ′ Γ j =1 the set, in ascending order, of both the α -values of the intersectingpoints and the mesh points of the reference interval I , α j = j/J Γ ( j = 0 , · · · , J Γ ).Then the inner products involving interface and bulk quantities are approximated bythe Simpson’s rule on the mesh (cid:8) α ′ j (cid:9) J ′ Γ j =1 ,(3.16)( u, v ) Γ m = 16 J ′ Γ X j =1 (cid:12)(cid:12)(cid:12) X m ( α ′ j ) − X m ( α ′ j − ) (cid:12)(cid:12)(cid:12)h ( u · v )( α ′ + j − ) + 4( u · v )( α ′ j − ) + ( u · v )( α ′− j ) i , and the inner products involving only quantities on the interface are simply approxi-mated by the trapezoidal rule on the mesh { α j } J Γ j =1 , (cid:0) u, v (cid:1) h Γ m = 12 J Γ X j =1 (cid:12)(cid:12)(cid:12) X m ( α j ) − X m ( α j − ) (cid:12)(cid:12)(cid:12)h(cid:0) u · v (cid:1) ( α + j − ) + (cid:0) u · v (cid:1) ( α − j ) i , (3.17)where α ′ j − = (cid:0) α ′ j − + α ′ j (cid:1) , u ( α ′± j ) are the one-sided limits of u at α ′ j , and similarlyfor u ( α ± j ).2 Q. Zhao and W. Ren (b)(a)
Fig. 3.2 . The boundary integral method is applied to the domain D = Ω m (upper panel) andthe domain D = Ω (lower panel) separately. The electrostatic potential Φ and its directionalderivative Ψ are evaluated at the discrete points indicated by ” ◦ ” and ” × ” respectively along theboundaries. The electric force on the fluid inter-face is computed by solving Eqs. (2.15)-(2.18) using the boundary integral method.Without loss of generality, we assume that Φ satisfies the Laplace equation on a do-main D with boundary Σ. Denote the directional derivative of Φ in the outwardnormal direction of the boundary by Ψ = n · ∇ Φ, where n is the outward unit normalvector on Σ. For any point p lies on the smooth part of Σ, we have the boundaryintegral equation(3.18) 12 Φ( p ) = Z Σ (cid:20) Φ( q ) ∂G ( p , q ) ∂ n ( q ) − Ψ( q ) G ( p , q ) (cid:21) d s ( q ) , where G ( p , q ) = π ln | p − q | is the Green function for the Laplace equation in R ,and ∂G ( p , q ) ∂ n ( q ) = n · ∇ q G ( p , q ).Eq. (3.18) can be used to compute Φ and/or its derivative Ψ on the boundaryΓ. To that end, we partition Σ then approximate it by a collection of line segments:Σ ≃ ∪ Mj =1 Σ ( j ) . We further approximate Φ( p ) and Ψ( p ) by constants on each linesegment: Ψ ≃ Φ j and Ψ ≃ Ψ j on Σ ( j ) , 1 ≤ j ≤ M . Denote by p j the mid-point ofthe line segment Σ ( j ) . It lies on a smooth part of the approximation boundary Σ ( j ) .Then we can apply Eq. (3.18) to p = p j and obtain(3.19) 12 Φ j = M X k =1 h A jk Φ k − B jk Ψ k i , j = 1 , · · · , M, where(3.20) A jk = Z Σ ( k ) ∂G ( p j , q ) ∂ n ( q ) d s ( q ) , B jk = Z Σ ( k ) G ( p j , q )d s ( q ) . For the current problem, we apply the boundary integral method to the domain D = Ω m and the domain D = Ω separately. A discretization of the boundaryof the two domains is shown in Fig. 3.2. Eq. (3.19) is applied at each discretepoint. These equations, together with the prescribed Dirichlet boundary conditionsΦ | Γ m = Φ | Γ = 1 and Φ | Γ = Φ | Γ = 0, the periodic boundary conditions Φ | Γ = Φ | Γ ′ ,Ψ | Γ = − Ψ | Γ ′ , Φ | Γ = Φ | Γ ′ and Ψ | Γ = − Ψ | Γ ′ , as well as the interface conditions Finite Element Method for Electrowetting on Dielectric = 0 and [ ǫ Ψ] = 0 on Γ m form a system of linear algebraic equations for Φ andΨ at the discrete points. After the linear system is solved, Ψ is used to compute theelectric force: |∇ Φ | = Ψ on Γ m . For the numerical method (3.15a)-(3.15d), we can show that it yields a unique solution. Furthermore, in the special casewhen the electric force is absent, the numerical method is unconditionally energystable. We make the following assumptions on the interface Γ m , ∀ ≤ m ≤ M ,i) the interface Γ m does not intersect with itself;ii) the parameterization is such that | ∂ α X m | >
0, andiii) the first and last segments of Γ m are not parallel to the x -axis.These assumptions particularly imply that the mesh points on Γ m do not merge, andthe dynamic contact angle is not 0 or π . Theorem 3.1 (
Existence and Uniqueness ). Let the finite element spaces ( U h , P h ) satisfy the inf-sup stability condition (3.13) and the fluid interface Γ m satisfy the ab-voe assumptions i)–iii). Then the numerical method (3.15a) - (3.15d) admits a uniquesolution. Theorem 3.2 (
Unconditional Energy Stability ). Assume the electric force iszero. Let (cid:0) u m +1 , p m +1 , X m +1 , κ m +1 (cid:1) be the solution to the numerical scheme (3.15a) - (3.15d) . Then the following stability bound holds ˜ W ( ρ m , u m +1 , X m +1 ) + 2 τRe k√ µ m D ( u m +1 ) k + τRe · l s (cid:16) β m u m +1 s , u m +1 s (cid:17) Γ m ∪ Γ m + β ∗ Re · τ h(cid:0) x m +1 r − x mr (cid:1) + (cid:0) x m +1 l − x ml (cid:1) i ≤ ˜ W ( ρ m − , u m , X m ) , (3.21) where ˜ W ( ρ, u , X ) = ( ρ u , u ) − cos θ Y W e | Γ | + W e | Γ | is the total energy of the system.Moreover, for m ≥ , we have ˜ W ( ρ m − , u m , X m ) + 2 τRe m − X k =0 k p µ k D ( u k +1 ) k + τRe · l s m − X k =0 (cid:0) β k u k +1 s , u k +1 s (cid:1) Γ k ∪ Γ k + β ∗ Re · τ m − X k =0 h(cid:0) x k +1 r − x kr (cid:1) + (cid:0) x k +1 l − x kl (cid:1) i ≤ ˜ W ( ρ , u , X ) . (3.22)The three summation terms in (3.22) correspond to the energy dissipation dueto the viscous stress in the bulk, the friction force on the wall and the contact linefriction, respectively.The above theorems are extensions of the previous work by Barrett et al. [3] toproblems with moving contact lines. In another related work [43], similar results wereobtained for moving contact lines but on fitted meshes. There the energy stabilityonly holds locally in time due to the required interpolation of the velocity and densitybetween the fitted meshes at each time step. Here the global energy stability of thediscrete system is established on the unfitted mesh. The proof of the theorems issimilar to that in Ref. [43] and is provided in the Appendix B and C.The discrete scheme (3.15c)-(3.15d) introduces an implicit tangential velocity forthe mesh points along the fluid interface such that they tend to be uniformly dis-tributed according to the arc length in long time [2, 42]. In our numerical experimentspresented below, no re-meshing for the fluid interface is needed during the simulation.4 Q. Zhao and W. Ren -3 -3 -3 -2 -1 s -3 -2 h (c) (d)(b)(a) Fig. 4.1 . Upper panel: the electric force along the interface (a semi-circle) computed usingdifferent mesh sizes for ǫ = ǫ = 1 (left) and ǫ = 0 . , ǫ = 1 (right). Lower left: the log-logplot of the electric force versus the arc length computed using h = 1 / . Lower right: the relativeerrors defined in (4.3) versus the mesh sizes h .
4. Numerical results.
In this section, we present some numerical results forEWoD obtained using the proposed numerical method. We first test the accuracyand convergence of the boundary integral method for the computation of the electricforce on a given fluid interface in section 4.1. Then we test the convergence of thenumerical method (3.15a)-(3.15d) using the example of a spreading droplet in section4.2. Subsequently, we examine the effect of the different parameters in the model onthe equilibrium interface profiles in section 4.3. Finally, in section 4.4 the numericalmethod is applied to the spreading and migration dynamics of a droplet on varioussubstrates.Unless otherwise stated, we will choose ρ = 1, ρ = 0 . µ = 1, µ = 0 . β = 1, β = 0 . β ∗ = 0 . Re = 10, Ca = 0 . l s = 0 . u = in the numerical examples. The computational domain occupied by the fluidsis Ω = [ − , × [0 , The electrostatic potential Φsatisfies the Maxwell’s equations in the domain Ω ∪ Ω . This is outside the domainΩ , which has a wedge-like geometry with an open angle θ at the contact points x l,r .The open angle is equal to the contact angle. In this geometry, the electric force |∇ Φ | in the vicinity of the contact point behaves as [6](4.1) |∇ Φ | ∼ O ((∆ s ) ν − ) , as ∆ s → + , where ∆ s is the distance to the contact point, and ν ∈ ( ,
1) is related to the contactangle and satisfies(4.2) ǫ tan [ ν ( π − θ )] + ǫ tan ( νπ ) = 0 . In particular, we have ν = π π − θ when ǫ = ǫ = 1. Finite Element Method for Electrowetting on Dielectric t -0.4-0.35-0.3-0.25 t -0.4-0.38-0.36-0.340 0.5 1 1.5 2 t t (a)(c) (b)(d) Fig. 4.2 . The evolution of the (left) contact point (upper panels) and the contact angle (lowerpanels) computed using different mesh sizes and time steps. In the coarse mesh, the mesh size is h = h = 1 / for the interface and h Ω = h = 1 / in the bulk, and the time step is τ = 0 . . Leftpanels: η = 0 ; Right panels: η = 0 . . To assess the performance of the boundary integral method for the current prob-lem, we compute the electric force along a given interface. The interface is the semi-circle of radius r = 0 . , d = 0 .
2. The numerical results are shown in Fig. 4.1. In the two upper panels,we plot the electric force along the interface. The different curves are obtained usingdifferent mesh sizes for ǫ = ǫ = 1 (left) and ǫ = 0 . , ǫ = 1 (right). It is evidentthat the electric force exhibits a singular behavior at the contact points: the force atthe contact point (more precisely, the force at half grid point away from the contactpoint) keeps increasing as the mesh is refined.To further examine the singular structure, we depict the log-log plot of the electricforce against the arc-length in the lower-left panel of the figure. It can be seen that theelectric force behaves as |∇ Φ | ∼ O ( s p ) as the contact line is approached. This is ingood agreement with the the theoretical prediction of Eqs. (4.1)-(4.2), which is shownas the straight lines with slope 2( ν − ≈ − .
667 for ǫ = 1 and 2( ν − ≈ − . ǫ = 0 .
25. These lines fit the numerical results very well.To examine the convergence of the numerical solution as the mesh is refined, wecompute the relative error defined as(4.3) e ( h ) = R I (cid:12)(cid:12)(cid:12) |∇ Φ h | − |∇ Φ h | (cid:12)(cid:12)(cid:12) d α R I |∇ Φ h | d α , where ∇ Φ h denotes the numerical solution obtained with mesh size h . The error fordifferent mesh sizes is shown in the lower-right panel of Fig. 4.1. As expected, theerror decreases as the mesh is refined. In this simulation, we used the uniform meshalong the interface. In practice, one may use local mesh refinement near the contactpoint to obtain more accurate solutions. We assess the perfor-mance of the numerical method (3.15a)-(3.15d) for the contact line dynamics by car-6
Q. Zhao and W. Ren rying out simulations under different mesh sizes. We consider the spreading dynamicsof a droplet. Initially the fluid interface is given by a semi-circle with center (0 , r = 0 .
4. The equilibrium contact angle of the droplet is θ Y = 2 π/
3. Thethickness of the dielectric substrate is d = 0 .
2, and the permittivities are ǫ = 1 and ǫ = 1.We use uniform mesh (see Fig. 3.1), and denote the mesh size in the bulk by h Ω = 1 /J Ω and the mesh size on the reference domain of the interface Γ by h = 1 /J Γ .In the boundary integral method, all the boundaries and interfaces are discretizedinto line segments; for example, the fluid-solid interface Γ ∪ Γ is also discretized into J Γ line segments. Different values of J Ω and J Γ will be used in the mesh refinementfor the convergence test. For simplicity, we shall fix the discretization of the outerboundaries.The evolution of the (left) contact point is shown in Fig. 4.2 for η = 0 (upper-leftpanel) and η = 0 . Table 4.1
Relative change of the droplet area at t = 4 . h, h Ω are the mesh size in the discretizationof the interface and Ω , respectively, and τ is the time step. In the coarse mesh, h = h = 1 / , h Ω = h = 1 / and τ = 0 . . ( h, h Ω , τ ) η = 0 η = 0 . η = 0 . | ∆ A | ( t = 4) order | ∆ A | ( t = 4) order | ∆ A | ( t = 4) order( h , h , τ ) 2.16E-3 - 3.38E-4 - 3.11E-4 -( h , h , τ ) 6.28E-4 1.78 7.51E-5 2.17 1.50E-4 1.05( h , h , τ ) 1.70E-4 1.88 1.77E-5 2.09 4.82E-5 1.64( h , h , τ ) 4.33E-5 1.97 4.14E-6 2.10 1.36E-5 1.82 Table 4.2
Convergence of the contact angle to the equilibrium angle at t = 4 . h, h Ω are the mesh size inthe discretization of the interface and Ω , respectively, and τ is the time step. In the coarse mesh, h = h = 1 / , h Ω = h = 1 / , τ = 0 . . ( h, h Ω , τ ) η = 0 η = 0 . η = 0 . | ∆ θ | ( t = 4) order | ∆ θ | ( t = 4) order | ∆ θ | ( t = 4) order( h , h , τ ) 7.12E-2 - 2.00E-1 - 3.79E-1 -( h , h , τ ) 3.54E-2 1.01 1.40E-1 0.51 3.11E-1 0.29( h , h , τ ) 1.77E-2 1.00 9.74E-2 0.52 2.27E-1 0.45( h , h , τ ) 8.90E-3 0.99 6.61E-2 0.55 1.53E-1 0.57Next we examine the conservation of area for the droplet. The finite element space P h for the pressure given in (3.12) contains piecewise constant functions, which ensuresthe local mass conservation particularly over each element. Besides, by choosing ζ = χ Ω1( t ) in (3.8b) and ψ = 1 in (3.8c), we can establish the mass conservation law Finite Element Method for Electrowetting on Dielectric t | Ω ( t ) | = (cid:16) n · ∂ t X , (cid:17) Γ( t ) = (cid:16) u · n , (cid:17) Γ( t ) = Z Ω ( t ) ∇ · u d L = (cid:16) ∇ · u , χ Ω1( t ) (cid:17) = 0 . (4.4)Therefore in this and the following simulations, we further enrich the pressure space P h with the basis function χ Ω m , the characteristic function over the domain occupiedby the droplet. This helps preserving the area of the droplet [3]. In addition, thepressure jump across the interface can be captured with this function. In practicalnumerical implementations, the new contribution of the single basis function to (3.15a)and (3.15b) can be written in terms of the integrals over Γ m as(4.5) (cid:16) ∇ · ω h , χ Ω m (cid:17) = Z Ω m ∇ · ω h d L = (cid:0) ω h , n m (cid:1) h Γ m , ∀ ω h ∈ U h . In Table. 4.1, we present the relative area change | ∆ A | of the droplet at t = 4.At this time the droplet has evolved close to the steady state. Here, ∆ A is defined as(4.6) ∆ A h ( t m ) = | Ω m | − | Ω || Ω | , where | Ω m | is the area of the droplet at time t m . From the table, we observe that thedroplet area is well-preserved. As the mesh is refined, the area change | ∆ A | convergesto zero with order close to 2.After the droplet reaches the steady state, the theoretical value of the contactangle should converge to the equilibrium angle θ Y = 2 π/
3. In table 4.2, we present∆ θ , the deviation of the contact angle from this equilibrium angle at t = 4, obtainedwith different mesh sizes. We observe that in all three cases for the different valuesof η , the contact angle indeed converges to the equilibrium angle. When η = 0, i.e.in the absence of electric field, the convergence order is close to 1; however, the orderis reduced to about 0.5 for η = 0 . , . h = X h and κ h the numerical solution for the interface and its curvatureat the steady state, respectively. Then from (3.15d), we have(4.7) (cid:16) κ h n h , g h (cid:17) h Γ h + (cid:16) ∂ s X h , ∂ s g h (cid:17) h Γ h − cos θ Y (cid:2) g h (1) − g h (0) (cid:3) = 0 , ∀ g h ∈ X h . By choosing g h = ( g h , g h ) T := ( φ ( α ) , T in (4.7), where φ ( α ) is the piecewiselinear function taking the value 1 at α = 0 and 0 at all other nodes, we obtain(4.8) 12 κ h (0) n h , (cid:12)(cid:12) X h ( α ) − X h ( α ) (cid:12)(cid:12) − cos θ h + cos θ Y = 0 , where n h , is the first component of n h , and θ h is the contact angle of X h . This yields(4.9) (cid:12)(cid:12) cos θ h − cos θ Y (cid:12)(cid:12) = 12 (cid:12)(cid:12) κ h (0) sin θ h (cid:12)(cid:12) ∆ s, where ∆ s = (cid:12)(cid:12) X h ( α ) − X h ( α ) (cid:12)(cid:12) . On the other hand, the interface condition im-plies that the curvature of the fluid interface has the same singular structure as the8 Q. Zhao and W. Ren -0.6 -0.4 -0.2 0 0.2 0.400.20.40.60.8 0 0.2 0.4 0.6 0.8 1-0.6-0.4-0.200.20.4 (a) (b)
Fig. 4.3 . Left panel: Equilibrium profiles of the interface. Right panel: The value of cos θ app versus η/d (discrete points with dash lines) and the prediction of Lippmann equation (straight line). electric field at the contact point. Therefore, we have κ h (0) ∼ O (cid:0) (∆ s ) ν − (cid:1) , andconsequently(4.10) (cid:12)(cid:12) θ h − θ Y (cid:12)(cid:12) ≤ C (∆ s ) ν − ≤ C h ν − , where C and C are constants independent of the mesh size h . When η = 0, we have ν = 1 and the curvature is constant along the interface at the steady state, thereforethe convergence order of the contact angle is 1. On the other hand, in the presence ofthe electric field, we have 1 / < ν <
1, and the convergence order is lowered. In thecurrent example, we have ν = 3 / (cid:12)(cid:12) θ h − θ Y (cid:12)(cid:12) ∼ O ( h / ), which isconsistent with the numerical results. In this example, we investigate the influ-ence of the various model parameters, such as η , ǫ and d , on the equilibrium profile ofthe fluid interface. The fluid interface of the droplet is initially given by a semi-circlewith centre (0 ,
0) and radius r = 0 .
4. The equilibrium contact angle of the dropletis θ Y = 2 π/
3. The fluid domain Ω = [ − , × [0 ,
1] is discretized into 4848 triangleswith 2534 vertices, and the fluid interface is discretized into J Γ = 512 elements. Thetime step is τ = 2 × − .The numerical results are shown in Fig. 4.3. The left panel shows the equilibriumprofiles of the interface for different values of η , ǫ and d . Comparing with the interfaceprofile when η = 0 (i.e. in the absence of electric field), we see that the electric forceflattens the interface and make the droplet spread.A more quantitative assessment of the effect of the electric force is shown in theright panel, where we plot cos θ app against η/d for different values of ǫ and d . Theapparent contact angle θ app is computed by fitting the interface by a circular arc usingthe apex of the interface and the given area of the droplet. From the numerical results,we observe that the ratio η/d plays the dominant role here; more specifically, cos θ app increases linearly with η/d with the slope close to 1. In contrast, the permittivityof the fluid outside the droplet, ǫ , has little effect on the interface profile. This isin good agreement with the Lippmann equation (1.2). In terms of the dimensionlessparameters, Eq. (1.2) reads(4.11) cos θ B = cos θ Y + ηd . The contact angle θ B computed using this equation is also shown in the figure, andgood agreement with the numerical results can be observed. The discrepancy might Finite Element Method for Electrowetting on Dielectric -1 -0.5 0 0.5 100.51 (a) -1 -0.5 0 0.5 100.51 (b) -1 -0.5 0 0.5 100.51 (c) -1 -0.5 0 0.5 100.51 (d) Fig. 4.4 . Snapshots of the interface and the velocity field on a hydrophilic dielectric substratewith θ Y = π/ . Parameters are η = 0 . , ǫ = ǫ = 1 and d = 0 . . (a) t = 0 ; (b) t = 0 . , max x ∈ Ω | u | = 1 . (c) t = 0 . , max x ∈ Ω | u | = 0 . ; (d) t = 1 . , max x ∈ Ω | u | = 0 . . be due to the finite system size, the effect of the boundary conditions, or the finitevalue of d and η . After all, the Lippmann equation is an asymptotic result which isderived in the limit of d → η → -1 -0.5 0 0.5 100.51 (a) -1 -0.5 0 0.5 100.51 (b) -1 -0.5 0 0.5 100.51 (c) -1 -0.5 0 0.5 100.51 (d) Fig. 4.5 . Snapshots of the interface and the velocity field on a hydrophobic dielectric substratewith θ Y = 2 π/ . Parameters are η = 0 . , ǫ = ǫ = 1 and d = 0 . . (a) t = 0 ; (b) t = 0 . , max x ∈ Ω | u | = 0 . ; (c) t = 0 . , max x ∈ Ω | u | = 0 . ; (d) t = 1 . , max x ∈ Ω | u | = 0 . . We investigate the detailed dynamics of a droplet on dielec-tric substrates driven by the electric force. We consider a hydrophilic case with thecontact angle θ Y = π/ θ Y = 2 π/ Q. Zhao and W. Ren t -0.4-0.3-0.2-0.10 0 0.5 1 1.5 t (b)(a) Fig. 4.6 . The energy loss ∆ W ( t ) = W ( t ) − W (0) (left panel) and the kinetic energy W k ( t ) = R Ω ρ | u | d L (right panel) as functions of time. Here the discrete energy W ( t ) is computed usingthe dimensionless from of (A.2) . the droplet to spread on the substrate.In Fig. 4.6, we plot the loss of the total energy ∆ W := W ( t ) − W (0) (left panel)and the kinetic energy of the fluids W k ( t ) = R Ω 12 ρ | u | d L (right panel) as functionsof time. The total energy, as given in Eq. (2.19), consists of the kinetic energy, theinterfacial energies and the electrostatic energy. We observe that the total energy ofthe discrete system decays in time, a desired property of the numerical method. -1 -0.5 0 0.5 100.51 (a) -1 -0.5 0 0.5 100.51 (b) -1 -0.5 0 0.5 100.51 (c) -1 -0.5 0 0.5 100.51 (d) -1 -0.5 0 0.5 100.51 (e) -1 -0.5 0 0.5 100.51 (f) Fig. 4.7 . Snapshots of the droplet migrating on a dielectric substrate. The electrostatic potentialon the lower boundary of the substrate is prescribed as Φ | y = − d = ( x + 1) . Other parameters are η = 0 . , ǫ = ǫ = 1 , d = 0 . , and θ Y = 2 π/ . (a) t = 0 ; (b) t = 0 . ; max x ∈ Ω | u | = 0 . ; (c) t = 0 . ; max x ∈ Ω | u | = 0 . ; (d) t = 1 . , max x ∈ Ω | u | = 0 . ; (e) t = 1 . , max x ∈ Ω | u | = 0 . ; (f) t = 2 . , max x ∈ Ω | u | = 0 . . In the last example, we simulate the migration of a droplet on a substrate withnon-uniform electrostatic potential prescribed on the boundary. The non-uniformpotential mimics an array of electrodes placed below the substrate that is used totransport the droplet in experiments. The electrostatic potential on the lower bound-ary of the substrate is given as Φ | y = − d = ( x + 1). Other parameters are chosen as η = 0 . ǫ = ǫ = 1, d = 0 .
2, and θ Y = 2 π/
3. The initial interface of the dropletis given by a semi-ellipse as ( x − . . + y . = 1 , y ≥
0. Numerical results for theinterface profile and the velocity field are shown in Fig. 4.7. We observe that thedroplet first de-wets the substrate to form a near-circular shape with contact angleclose to θ Y = 2 π/
3. In this stage, the dynamics is mainly driven by the surface ten-
Finite Element Method for Electrowetting on Dielectric
5. Conclusions.
In this work, we presented a hydrodynamic model for elec-trowetting on dielectric based on the earlier work on moving contact lines, and devel-oped an efficient numerical method for the model. The numerical method combines asemi-implicit parametric finite element method for the dynamics of the fluid interfaceand the finite element method for the Navier-Stokes equations as well as the boundaryintegral method for the electric field. We proved that the numerical method admits aunique solution. In the case without the electric field, we showed that the numericalmethod obeys a similar energy law as the continuum model.In the numerical experiments, we assessed the accuracy and convergence of thenumerical method and investigated the effect of the different physical parameters onthe interface dynamics and its equilibrium profile. Numerical results for the equilib-rium profile of the interface agree well with the predictions of the Lippmann equation.The numerical solution for the electric force exhibits a singular structure near thecontact point that is consistent with theoretical results. This singularity incurs largecurvature of the interface near the contact point, which deteriorates the convergenceorder of the numerical solution, particularly the convergence of the contact angle tothe equilibrium angle.In this work, we focused on simulations in two dimensions. The numerical methodcan be readily extended to three-dimensional problems. This will be left to our futurework.
Acknowledgement.
The work was partially supported by Singapore MOE RSBgrant, Singapore MOE AcRF grants (No. R-146-000-267-114, No. R-146-000-285-114)and NSFC grant (No. 11871365).
Appendix A. Energy law for the continuum model.
The total energy forthe EWoD model (2.1)-(2.7) in its original dimension form reads(A.1) W ( t ) = X i =1 , Z Ω i ρ i | u | d L − γ cos θ Y | Γ ( t ) | + γ | Γ( t ) |− X i =2 , Z Ω i ǫ i |∇ Φ | d L . Integrating by parts and using the electrostatic potential equation in (2.5) as wellas the boundary conditions for Φ, we can transform the electrical energy in (A.1)into expressions only involving line integrals over Γ and Γ . This gives the followingalternative form of the total energy W ( t ) = X i =1 , Z Ω i ρ i | u | d L − γ cos θ Y | Γ ( t ) | + γ | Γ( t ) | + φ (cid:18)Z Γ ǫ ( n · ∇ Φ) d s + Z Γ ǫ ( n w · ∇ Φ) d s (cid:19) , (A.2)which can be used to compute the discrete electrical energy in view of the boundaryintegral method in (3.19). The dynamic system obeys the energy dissipation law(A.3) dd t W ( t ) = − X i =1 , Z Ω i µ i k D ( u ) k F d L − X i =1 , Z Γ i β i | u s | d s − β ∗ ( ˙ x l + ˙ x r ) ≤ . Q. Zhao and W. Ren
We show the proof of (A.3). The dissipation of the fluid kinetic energy isdd t X i =1 , Z Ω i ( t ) ρ i | u | d L = X i =1 , Z Ω i ρ i u · ( ∂ t u + u · ∇ u ) d L = X i =1 , Z Ω i u · ( −∇ p + ∇ · τ d ) d L = − X i =1 , Z Ω i ∇ u : τ d d L − X i =1 , Z Γ i u · τ d · n w d s + Z Γ u · [ p I − τ d ] · n d s = − X i =1 , Z Ω i µ i k D ( u ) k F d L − X i =1 , Z Γ i β i | u s | d s + Z Γ ( γκ + ǫ |∇ Φ | )( u · n ) d s, (A.4)where we have used the divergence free condition in (2.1b), the interface conditionsin (2.2a)-(2.2b), the boundary condition (2.3), the no-slip boundary condition on theupper wall Γ and the periodic boundary conditions along Γ .The time derivative of the interfacial energies isdd t (cid:16) − γ cos θ Y | Γ ( t ) | + γ | Γ( t ) | (cid:17) = − γ cos θ Y ( ˙ x r − ˙ x l ) − Z Γ γκ v n d s + γ (cid:0) ˙ x r cos θ rd − ˙ x l cos θ ld (cid:1) = − Z Γ γκ ( u · n ) d s − β ∗ (cid:0) ˙ x r + ˙ x l (cid:1) , (A.5)where we have used (2.2c) and (2.4).Denote by ˆ n the outward unit normal vector on the boundary of Ω and Ω .Differentiating the electrical energy yieldsdd t − X i =2 , Z Ω i ǫ i |∇ Φ | d L = − X i =2 , Z Ω i ǫ i ∇ Φ · ∇ ( ∂ t Φ) d L − X i =2 , Z ∂ Ω i ǫ i |∇ Φ | ( u · ˆ n ) d s = − X i =2 , Z Ω i ǫ i ∇ · ( ∂ t Φ ∇ Φ) d L + 12 Z Γ ǫ |∇ Φ | ( u · n ) d s = − X i =2 , Z ∂ Ω i ǫ i ∂ t Φ ( ∇ Φ · ˆ n ) d s + 12 Z Γ ǫ |∇ Φ | ( u · n ) d s = − Z Γ ǫ |∇ Φ | ( u · n ) d s, (A.6)where the first equality results from the Reynolds transport theorem, the secondequality is obtained from (2.5), the third equality comes from the divergence theorem,and for the last equality, we have used the boundary conditions in (2.6)-(2.7), theperiodic boundary condition along Γ ∪ Γ as well as the fact that(A.7) ∂ t Φ + u · ∇ Φ = 0 , on Γ ∪ Γ ∪ Γ ∪ Γ . Finite Element Method for Electrowetting on Dielectric ρ U L ). Appendix B. Proof of Theorem 3.1.
It suffices to show that the correspond-ing homogeneous system has only zero solution. By noting that electric force ǫ η |∇ Φ | in (3.15a) is explicitly evaluated on Γ m , thus we can consider solving the followinghomogeneous system for (cid:0) u h , p h , X h , κ h (cid:1) ∈ (cid:0) U h , P h , X h , K h (cid:1) such that12 h(cid:16) ( ρ m + ρ m − ) u h τ , ω h (cid:17) + (cid:16) ρ m ( u m · ∇ ) u h , ω h (cid:17) − (cid:16) ρ m ( u m · ∇ ) ω h , u h (cid:17)i − (cid:16) p h , ∇ · ω h (cid:17) + 2 Re (cid:16) µ m D ( u h ) , D ( ω h ) (cid:17) − W e (cid:16) κ h n m , ω h (cid:17) Γ m + 1 Re · l s (cid:16) β m u hs , ω hs (cid:17) Γ m ∪ Γ m = 0 , ∀ ω h ∈ U h , (B.1a) (cid:16) ∇ · u h , ζ h (cid:17) = 0 , ∀ ζ h ∈ P h , (B.1b) (cid:16) X h τ · n m , ψ h (cid:17) h Γ m − (cid:16) u h · n m , ψ h (cid:17) Γ m = 0 , ∀ ψ h ∈ K h , (B.1c) (cid:16) κ h n m , g h (cid:17) h Γ m + (cid:16) ∂ s X h , ∂ s g h (cid:17) h Γ m + β ∗ Caτ h x hr g h (1) + x hl g h (0) i = 0 , ∀ g h ∈ X h , (B.1d)where X h = ( X h , Y h ) T , u hs = u h · t w , and x hl := X h | α =0 and x hr := X h | α =1 .Taking ω h = u h , ζ h = p h , ψ h = W e κ h and g h = W e X h , then combining theseequations yields (cid:16) ρ m + ρ m − u h , u h (cid:17) + 2 τRe (cid:16) µ m D ( u h ) , D ( u h ) (cid:17) + τRe · l s (cid:16) β m u hs , u hs (cid:17) Γ m ∪ Γ m + 1 W e (cid:16) ∂ s X h , ∂ s X h (cid:17) h Γ m + β ∗ Re · τ [( x hr ) + ( x hl ) ] = 0 . (B.2)By Korn’s inequality, we have(B.3) k u h k ≤ C h (cid:16) ( ρ m + ρ m − ) u h , u h (cid:17) + 2 τRe (cid:16) µ m D ( u h ) , D ( u h ) (cid:17)i ≤ , thus we immediately obtain u h = . We also have X h = by noting x hr = x hl = 0.Substituting X h = into Eq. (B.1d), we obtain(B.4) (cid:16) κ h n m , g h (cid:17) h Γ m = 0 , ∀ g h ∈ X h . Denote n mj = ( n mj, , n mj, ) T , j = 1 , , · · · , J Γ . Choosing g h in (B.4) as(B.5) g h (cid:12)(cid:12) α j = (cid:2) h mj +1 + h mj ) (cid:3) ⊥ κ h ( α j ) , ≤ j ≤ J Γ − , (cid:0) n m , κ h ( α j ) , (cid:1) T , j = 0 , (cid:16) n mJ Γ , κ h ( α j ) , (cid:17) T , j = J Γ , Q. Zhao and W. Ren and by noting the norm in (3.17), we obtain0 = ( κ h ( α )) | h m | ( n m , ) + ( κ h ( α J Γ )) | h mJ Γ | ( n mJ Γ , ) + 12 J Γ − X j =1 (cid:0) κ h ( α j ) (cid:1) | h mj + h mj +1 | , (B.6)which implies κ h ( α j ) = 0 , ∀ ≤ j ≤ J Γ from the assmuptions i)–iii). We thensubstitute u h = and κ h = 0 into (B.1a) and obtain(B.7) (cid:0) p h , ∇ · ω h (cid:1) = 0 , ∀ ω h ∈ U h . Using the stability condition in (3.13), we consequently obtain p h = 0. This showsthat the homogeneous linear system (B.1a) - (B.1d) has only the zero solution. Thus,the numerical scheme (3.15a)-(3.15d) admits a unique solution. Appendix C. Proof of Theorem 3.2.
Setting ω h = u m +1 , ζ h = p m +1 , ψ h = W e κ m +1 and g h = W e · τ ( X m +1 − X m ) in Eqs. (3.15a)-(3.15d), noting η = 0and then combining these equations yield12 τ h(cid:16) ρ m u m +1 − ρ m − u m , u m +1 (cid:17) + (cid:16) ρ m − (cid:0) u m +1 − u m (cid:1) , u m +1 (cid:17)i + 2 Re (cid:16) µ m D ( u m +1 ) , D ( u m +1 ) (cid:17) + 1 Re · l s (cid:0) β m u m +1 s , u m +1 s (cid:1) Γ m ∪ Γ m + 1 W e · τ (cid:16) ∂ s X m +1 , ∂ s ( X m +1 − X m ) (cid:17) h Γ m − cos θ Y W e · τ (cid:2) ( x m +1 r − x m +1 l ) − ( x mr − x ml ) (cid:3) + β ∗ Re "(cid:18) x m +1 r − x mr τ (cid:19) + (cid:18) x m +1 l − x ml τ (cid:19) = 0 . (C.1)It is easy to see that the following inequalities hold: (cid:16) ρ m u m +1 − ρ m − u m , u m +1 (cid:17) + (cid:16) ρ m − (cid:0) u m +1 − u m (cid:1) , u m +1 (cid:17) ≥ (cid:16) ρ m u m +1 , u m +1 (cid:17) − (cid:16) ρ m − u m , u m (cid:17) . (C.2) (cid:16) ∂ s X m +1 , ∂ s ( X m +1 − X m ) (cid:17) h Γ m ≥ | Γ m +1 | − | Γ m | . (C.3)Using (C.2) and (C.3) in (C.1) and noting x m +1 r − x m +1 l = | Γ m +11 | , x mr − x ml = | Γ m | , weimmediately obtain Eq. (3.21). Replacing m by k in Eq. (3.21) and by summarisingup for k from 0 to m −
1, we obtain the global discrete energy dissipation law (3.22).
REFERENCES[1]
M. Agnese and R. N¨urnberg , Fitted finite element discretization of two-phase Stokes flow ,Int J Numer Methods Fluids., 82 (2016), pp. 709–729.[2]
J. W. Barrett, H. Garcke, and R. N¨urnberg , A parametric finite element method forfourth order geometric evolution equations , J. Comput. Phys., 222 (2007), pp. 441–467.[3]
J. W. Barrett, H. Garcke, and R. N¨urnberg , A stable parametric finite element discretiza-tion of two-phase Navier–Stokes flow , J. Sci. Comput., 63 (2015), pp. 78–117. Finite Element Method for Electrowetting on Dielectric [4] B. Berge , Electrocapillarit´e et mouillage de films isolants par l’eau , Comptes Rendus deL’Academie des Sciences Paris, Serie, II, 317 (1993), pp. 157–163.[5]
M. Bienia, M. Vallade, C. Quilliet, and F. Mugele , Electrical-field-induced curvatureincrease on a drop of conducting liquid , Europhys. Lett., 74 (2006), 103.[6]
J. Buehrle, S. Herminghaus, and F. Mugele , Interface profiles near three-phase contactlines in electric fields , Phys. Rev. Lett., 91 (2003), 086101.[7]
L. Chen and E. Bonaccurso , Electrowetting - from statics to dynamics , Adv. Colloid InterfaceSci., 210 (2014), pp. 2–12.[8]
S. K. Cho, H. Moon, and C.-J. Kim , Creating, transporting, cutting, and merging liquiddroplets by electrowetting-based actuation for digital microfluidic circuits , J. Microelec-tromechanical. Syst., 12 (2003), pp. 70–80.[9]
L. Clime, D. Brassard, and T. Veres , Numerical modeling of electrowetting processes indigital microfluidic devices , Comput & Fluids, 39 (2010), pp. 1510–1515.[10]
L. Corson, N. Mottram, B. Duffy, S. Wilson, C. Tsakonas, and C. Brown , Dynamicresponse of a thin sessile drop of conductive liquid to an abruptly applied or removedelectric field , Phys. Rev. E, 94 (2016), 043112.[11]
L. T. Corson, C. Tsakonas, B. R. Duffy, N. J. Mottram, I. C. Sage, C. V. Brown, andS. K. Wilson , Deformation of a nearly hemispherical conducting drop due to an electricfield: Theory and experiment , Phys. Fluids, 26 (2014), 122106.[12]
D. Crowdy , Exact solutions for the static dewetting of two-dimensional charged conductingdroplets on a substrate , Phys. Fluids, 27 (2015), 061705.[13]
H. Cui and W. Ren , Interface profile near the contact line in electro-wetting on dielectric ,SIAM J. Appl. Math., 80 (2020), pp. 402–421.[14]
C. D. Daub, D. Bratko, K. Leung, and A. Luzar , Electrowetting at the nanoscale , J. Phys.Chem. C, 111 (2007), pp. 505–509.[15]
H. B. Eral, D. M. Augustine, M. H. Duits, and F. Mugele , Suppressing the coffee staineffect: how to control colloidal self-assembly in evaporating drops using electrowetting , SoftMatter, 7 (2011), pp. 4954–4958.[16]
M. Fontelos, G. Gr¨un, U. Kindelan, and F. Klingbeil , Numerical simulation of static anddynamic electrowetting , J. Adhes. Sci. Technol., 26 (2012), pp. 1805–1824.[17]
Y. Guan, B. Li, and L. Xing , Numerical investigation of electrowetting-based droplet splittingin closed digital microfluidic system: Dynamics, mode, and satellite droplet , Phys. Fluids,30 (2018), 112001.[18]
R. A. Hayes and B. J. Feenstra , Video-speed electronic paper based on electrowetting , Nature,425 (2003), 383.[19]
H. K. Kang , How electrostatic fields change contact angle in electrowetting , Langmuir, 18(2002), pp. 10318–10322.[20]
S. Kuiper and B. Hendriks , Variable-focus liquid lens for miniature cameras , Appl. Phys.Lett., 85 (2004), pp. 1128–1130.[21]
A. Kutana and K. Giapis , Atomistic simulations of electrowetting in carbon nanotubes , Nano.Lett., 6 (2006), pp. 656–661.[22]
H. Li and H. Fang , Lattice Boltzmann simulation of electrowetting , Euro. Phys. J. SpecialTopics, 171 (2009), pp. 129–133.[23]
G. Lippmann , Relations entre les ph´enom`enes ´electriques et capillaires , PhD thesis, Gauthier-Villars Paris, France:, 1875.[24]
H.-W. Lu, K. Glasner, A. Bertozzi, and C.-J. Kim , A diffuse-interface model for electrowet-ting drops in a Hele-Shaw cell , J. Fluid. Mech., 590 (2007), pp. 411–435.[25]
J. Monnier, P. Witomski, P. Chow-Wing-Bom, and C. Scheid , Numerical modeling ofelectrowetting by a shape inverse approach , SIAM J. Appl. Math., 69 (2009), pp. 1477–1500.[26]
H. Moon, S. K. Cho, R. L. Garrell, and C.-J. C. Kim , Low voltage electrowetting-on-dielectric , J. Appl. Phys., 92 (2002), pp. 4080–4087.[27]
F. Mugele and J.-C. Baret , Electrowetting: from basics to applications , J. Phys. Condens.Matter, 17 (2005), R705.[28]
F. Mugele and J. Buehrle , Equilibrium drop surface profiles in electric fields , J. Phys.Condens. Matter, 19 (2007), 375112.[29]
M. M. Nahar, G. S. Bindiganavane, J. Nikapitiya, and H. Moon , Numerical modeling of3D electrowetting droplet actuation and cooling of a hotspot , in Proceedings of the 2015COMSOL Conference, Boston, MA, USA, 2015, pp. 7–9.[30]
R. H. Nochetto, A. J. Salgado, and S. W. Walker , A diffuse interface model for elec-trowetting with moving contact lines , Math. Mod. Meth. Appl. Sci., 24 (2014), pp. 67–111.[31]
M. G. Pollack, A. D. Shenderov, and R. B. Fair , Electrowetting-based actuation of droplets Q. Zhao and W. Ren for integrated microfluidics , Appl. Phys. Lett., 2 (2002), pp. 96–101.[32]
C. Quilliet and B. Berge , Electrowetting: a recent outbreak , Curr. Opin. Colloid InterfaceSci., 6 (2001), pp. 34–39.[33]
W. Ren and W. E , Boundary conditions for the moving contact line problem , Phys. Fluids,19 (2007), 022101.[34]
W. Ren and W. E , Derivation of continuum models for the moving contact line problem basedon thermodynamic principles , Commun. Math. Sci., 9 (2011), pp. 597–606.[35]
W. Ren, D. Hu, and W. E , Continuum models for the contact line problem , Phys. Fluids, 22(2010), 102103.[36] ´E. Ruiz-Guti´errez and R. Ledesma-Aguilar , Lattice-Boltzmann simulations of electrowet-ting phenomena , Langmuir, 35 (2019), pp. 4849–4859.[37]
C. Scheid and P. Witomski , A proof of the invariance of the contact angle in electrowetting ,Math. Comput. Model., 49 (2009), pp. 647–665.[38]
B. Shapiro, H. Moon, R. L. Garrell, and C.-J. C. Kim , Equilibrium behavior of sessile dropsunder surface tension, applied external fields, and material variations , J. Appl. Phys., 93(2003), pp. 5794–5811.[39]
M. Vallet, M. Vallade, and B. Berge , Limiting phenomena for the spreading of water onpolymer films by electrowetting , Eur. Phys. J. B, 11 (1999), pp. 583–591.[40]
H. Verheijen and M. Prins , Reversible electrowetting and trapping of charge: model andexperiments , Langmuir, 15 (1999), pp. 6616–6620.[41]
S. W. Walker and B. Shapiro , Modeling the fluid dynamics of electrowetting on dielectric(EWOD) , J. Microelectromech. Syst., 15 (2006), pp. 986–1000.[42]
Q. Zhao, W. Jiang, and W. Bao , An energy-stable parametric finite element method forsimulating solid-state dewetting , (2020), arXiv:2003.01677 [math.NA].[43]
Q. Zhao and W. Ren , An energy-stable finite element method for the simula-tion for moving contact lines in two-phase flows , J. Comput. Phys., 417 (2020),doi: .[44]