Decoupled, linear, unconditionally energy stable and charge-conservative finite element method for a inductionless magnetohydrodynamic phase-field model

by   Xiaorong Wang, et al.
Chinese Academy of Science

In this paper, we consider the numerical approximation for a diffuse interface model of the two-phase incompressible inductionless magnetohydrodynamics problem. This model consists of Cahn-Hilliard equations, Navier-Stokes equations and Poisson equation. We propose a linear and decoupled finite element method to solve this highly nonlinear and multi-physics system. For the time variable, the discretization is a combination of first-order Euler semi-implicit scheme, several first-order stabilization terms and implicit-explicit treatments for coupling terms. For the space variables, we adopt the finite element discretization, especially, we approximate the current density and electric potential by inf-sup stable face-volume mixed finite element pairs. With these techniques, the scheme only involves a sequence of decoupled linear equations to solve at each time step. We show that the scheme is provably mass-conservative, charge-conservative and unconditionally energy stable. Numerical experiments are performed to illustrate the features, accuracy and efficiency of the proposed scheme.



page 14

page 16

page 17

page 18


An energy stable C^0 finite element scheme for a quasi-incompressible phase-field model of moving contact line with variable density

In this paper, we focus on modeling and simulation of two-phase flow wit...

Analysis of a fully discrete approximation for the classical Keller–Segel model: lower and a priori bounds

This paper is devoted to constructing approximate solutions for the clas...

Broadband Finite-Element Impedance Computation for Parasitic Extraction

Parasitic extraction is a powerful tool in the design process of electro...

Assessment of a non-conservative four-equation multiphase system with phase transition

This work focuses on the formulation of a four-equation model for simula...

A Conservative Finite Element Solver for MHD Kinematics equations: Vector Potential method and Constraint Preconditioning

A new conservative finite element solver for the three-dimensional stead...

A theoretical and experimental investigation of a family of immersed finite element methods

In this article we consider the widely used immersed finite element meth...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In recent years, the phase field method, also called the diffuse interface method, has been widely used to simulate the motion of multiphase incompressible immiscible fluids for both numerical and theoretical research anderson_diffuse-interface_1998; ding_diffuse_2007; liu_phase_2003, such as in fields like material sciences, fracture mechanics and fluid mechanics. Unlike the traditional sharp interface model, the phase field method describes the interface by a balance of molecular forces in a very thin layer rather than a free curve that evolves over time. To indicate different phases, an auxiliary function called phase field or order parameter is introduced rayleigh_theory_1892; waals_thermodynamic_1893. The phase function takes a distinct constant value in each bulk phase, and varies smoothly in the interfacial region. In the diffuse interface theory, the surface motion of the order parameter is driven by a gradient flow. There are two popular gradient flow-type governing equations: Allen-Cahn equations allen_1979 and Cahn-Hilliard equations cahn_free_1958. Comparing with the classical sharp interface model, one principal advantage of phase field method is its ability to capture the interface implicitly and automatically, and allow for topological changes such as self-intersection, pinch-off and splitting. Thus, phase field method has become an attractive technique, for extensive study on the phase field approach, we refer to the recent reviews kim_phase-field_2012; shen_modeling_2011 and the references therein.

Magnetohydrodynamic (MHD) equations are used to describe the behavior of electrically conductive fluids under the influence of magnetic fields and electric currents, which couple the Navier-Stokes equations and Maxwell equations. The theoretical analysis and numerical simulation of incompressible MHD equations is an area of research currently undergoing intense study gerbeau_stabilized_2000; gerbeau_mathematical_2006; moreau_magnetohydrodynamics_1990. However, in most industrial and laboratory cases, the magnetic Reynolds number is small, then the induced magnetic field can be neglected compared with the applied magnetic field. In these cases, Maxwell equations are replaced by Poisson equation for the electric scalar potential, which yields the inductionless magnetohydrodynamic (IMHD) equations badia_block_2014; li_charge-conservative_2019; ni_current_2007. It is well-known that numerical simulation of multiphase flow is a challenging problem in computational fluid dynamics. This problem becomes more difficult for MHD due to the interaction between multiple physical fields. The multiphase MHD flow often is involved in a wide range of scientific and industrial problems, such as astrophysics, aluminum electrolysis, liquid metal magnetic pumps, MHD power generators, fusion reactor blankets abdou_2001; davidson_introduction_2001. Thus developing fast and effective algorithms for incompressible multiphase MHD equations has great theoretical significance and applicable value.

In this paper, we focus on the two-phase incompressible IMHD problem in a general Lipschitz domain, which is frequently performed to simulate the movement in aluminum electrolysis cell gerbeau_mathematical_2006; munger_level_2006 and explore the theoretical study and numerical experiments of liquid lithium-lead cladding abdou_us_2005. Traditionally, most of the existing models for multiphase incompressible IMHD problem are devoted to sharp interface models, for instance, the front-tracking method samulyak_numerical_2007; zhang_direct_2018, the volume-of-fluid method huang_3d_2002; pan_development_2012; takatani_mathematical_2007, the level-set method ki_level_2010; xie_tracking_2007. There are few works have focused on the two-phase incompressible IMHD equations by diffuse-interface method in the literature. In 2014, Ding et al. (Ding2014) presented a two-phase incompressible inductionless magnetohydrodynamics model and simulated the deformation of melt interface in an aluminum electrolytic cell. In 2020, Chen and his collaborators (Chen2020) proposed a linear, decoupled, unconditionally energy stable and second order time-marching schemes to simulate the two-phase incompressible Cahn-Hilliard-IMHD (CHIMHD) conducting flow. Zhang (Zhang2021) derived a diffuse interface model for IMHD fluids systematically and analyzed sharp interface limits for different choice of the mobilities. Mao et al. mao_fully_2021 analyzed the fully discrete finite element approximation of a three-dimensional CHIMHD model and proved the well-posedness of weak solution to the phase field model by using the classical compactness method.

Studying efficient numerical methods for solving incompressible CHIMHD flows is the main focus of this article. The CHIMHD model is a complex system that involves three coupled physical processes: the phase field model (Cahn-Hilliard equation), fluid dynamics (Navier-Stokes equations), and electric field (Possion equation). Due to the highly nonlinear with multi-physics fields coupled in this system, it is relatively complicated and time-consuming to solve it by nonlinear and coupled scheme (Chen2020; mao_fully_2021). Our primary goal is to propose a decoupled and efficient fully discrete mixed finite element method to solve the system. The time discretization is an Euler semi-implicit scheme with some extra first order stabilization terms in Cahn-Hilliard equations and Poisson equation, and implicit-explicit treatments for coupling terms. The subtle time splitting technique allows us to fully decouple system into three processes at each time step. For spacial discretization, we adopt classical inf-sup stable finite element to discretize the phase field and chemical potential, standard velocity-pressure stable finite element to discretize the hydrodynamic unknowns, and stable face-volume finite element pairs to discretize the current density and electric potential. With these techniques, we obtain an efficient fully discrete scheme, which is decoupled, linear while still mass-conservative, charge-conservative and unconditionally energy stable. Numerical experiments are performed to demonstrate the features, accuracy and efficiency of the proposed scheme.

We emphasize that the proposed scheme is charge-conservative and unconditionally energy stable. The charge conservation law, namely, is a physical law in electromagnetics, which plays an important role in keeping the calculation accuracy for the simulation of MHD fluid. It was believed that only when the numerical schemes numerically maintains the physical conservation laws, such as momentum conservation and charge conservation, can we get accurate results for MHD flows at high Hartmann number ni_consistent_2012; ni_current_2007. Due to the rapid changes near the interface, the non-compliance of energy dissipation laws of the scheme may lead to spurious numerical solutions (Hua2011; Johnston2004). Thus, the design of unconditionally energy stable schemes are of great importance for solving phase field models. To the best of our knowledge, the scheme presented in this paper is the first decoupled, linear and unconditionally energy stable scheme for a phase field model of two-phase inductionless MHD flows.

The rest of this paper is organized as follows. In section 2, we present the CHIMHD model, derive a weak formulation and the dissipative energy law. In section 3, we propose a fully discrete charge-conservation finite element method and prove its unconditional energy stability. In section 4, we present some numerical simulations to validate our schemes. Some concluding remarks are given in section 5.

2 The diffuse interface model

In this section, we first introduce the CHIMHD model, then present a weak formulation, and finally show the dissipative energy law in the PDE level.

2.1 Model system

We consider a phase-field model for mixture of two immiscible, incompressible and conducting fluids with the same density in a bounded domain with Lipschitz-continuous boundary in . Following (mao_fully_2021; Zhang2021), the dynamic behavior of the flow is governed by the coupling model of the Cahn-Hilliard equations and the IMHD equations. In this paper, we consider the incompressible CHIMHD system as follows:


where is the density of fluids, denote the velocity and pressure of fluids, represent the current density and electric scalar potential, denote the phase function and chemical potential, and is the applied magnetic field which is assumed to be given.

is the symmetric gradient tensor. In the Cahn-Hilliard equations (

5)-(6), , with is the double-well potential, is the surface tension coefficient, and is the interface thickness, represents the diffusional mobility related to the relaxation time scale. In the IMHD equations (1)-(4), the viscosity of the fluids and the electric conductivity are depending on the phase function . They are Lipschitz-continuous functions of satisfying

where and (i=1,2) are the viscosity and electric conductivity of the fluid on each phase. The system of equations are supplemented with the following initial values and boundary conditions

where is the outer unit normal of . In this paper, we consider two-phase flows with matching density . Without lose of generality, we set in the above system (2.1).

Remark 2.1.

The first term on the right-hand side of (1), is the Lorentz force. The last term on the left-hand side of (1) the continuum surface tension force in the potential form, which denotes the phase introduced force. This term appears differently in literature shen_decoupled_2015; shen_modeling_2011: or It can be shown that these three expressions are equivalent by redefining the pressure .

Remark 2.2.

Though the model considered here is a two-phase model with matched density, one can still employ this model with Boussinesq approximation to model the effect of density difference by a gravitational force in case of small density ratio. The large density ratio case is reserved for future work.

In this paper, we consider the truncated double well potential ,

The original definition of the potential is logarithmic (cahn_free_1958), which guarantees the phase field stays within (-1,1). The Ginzburg-Landau potential is an popular alternative to approximate the one suggested by Cahn and Hilliard. Following shen_decoupled_2015; shen_numerical_2010, we further restrict the growth of the potential to quadratic away from the range [-1,1]. It is proved by caffarelli_bound_1995 that the truncated also ensures the boundness of in the Cahn-Hilliard equation. Therefore, it is a common practice to consider the Cahn-Hilliard equation with the truncated double-well potential shen_phase-field_2010; shen_decoupled_2015. It is clear that the second derivative of is continuous and bounded,

This property is of great help in handling the nonlinear double-well potential by using stabilization method. For the treatment of the general double-well potential, we give some comments in Remark 3.7.

2.2 Weak formulation and energy estimate

Firstly, we introduce some useful notations and Sobolev spaces. Let be the space of square-integrable functions that equipped with the inner product and norm:

Its subspace with mean zero over is written as . Let be the subspaces of with square integrable gradients and square integrable divergences respectively. The equipped norm in is defined by

Their subspaces with vanishing traces and vanishing normal traces on are denoted by . We refer to monk_finite_2003 for their definitions and inner products.

For convenience, we introduce some notations for function spaces

To derive the weak formulation of (2.1), we define a trilinear form:

for any .

Armed with the above notation, a weak formulation of the system (2.1) amounts to find such that


for all .

Now, we are in a position to establish the energy estimate for the CHIMHD system.

Theorem 2.3.

Let be the solution of (7). The following energy dissipation law holds



Taking in (7), we have

Combining the these equalities, we obtain

Integrating both sides over yields the theorem. ∎

The energy law describes the evolution of the total energy caused by energy conversion. Since the induced magnetic field can be neglected and the electric field is considered to be quasi-static, the total energy only consists of the fluid kinetic energy and the Cahn-Hilliard free energy. The dissipation of stems from the friction losses , the Ohmic losses and the diffusion transport term .

3 Decoupled energy stable finite element method

In this section, we propose a decoupled, energy stable, mixed finite element scheme for continuous problem (7).

Let be a shape-regular simplex subdivision of . As usual, we introduce the local mesh size and global mesh size . Here we choose conforming finite element space pairs to discretize velocity and pressure , to approximate current density and electric potential , and to discretize the phase field function and chemical potential . In addition, we assume these spaces satisfy the following inf-sup conditions.

Proposition 3.4 (inf-sup condition).

The finite element pairs , and satisfy the following uniform inf-sup conditions:


where , and only depend on .

With these discrete spaces, the semi-discretization formulation of the system (7) is to find such that


for all .

With similar arguments in Theorem 2.3, one can easily get the energy stability of the semi-discrete scheme, thus the details are omitted here.

Theorem 3.5.

Let be the solution of (7). The energy dissipation law holds


Let be an equidistant partition of the time interval For any time dependent function , the full discrete approximation to will be denoted by . A fully discrete mixed finite element scheme for problem (7) reads as follows:

Given the initial guess and , we compute

by the following three steps.
Step 1: Find such that


for all .
Step 2: Find such that


for all .
Step 3: Find such that


for all .

Before we get into the discussion on the properties of this scheme, several remarks about the scheme are given in order.

Remark 3.6.

To decouple the nonlinear coupled multiphysics system, we introduce two first -order stabilization terms. The first stabilization term in step 1 is to decouple Cahn–Hilliard equations and Navier–Stokes equations shen_numerical_2010. The second stabilization term in step 2 is to decouple Poisson equation and Navier–Stokes equations Zhang2020. These extra two stabilization terms are vital to keep the couping term explicitly while preserving the energy stability, see the proof of Theorem 3.5.

Remark 3.7.

In step 1, we employ the stabilized method shen_numerical_2010 to treat the double-well potential explicitly without suffering from any time step constraint. Note that this stabilizing term introduces an extra consistent error of order , which is the same order as the overall truncation error of the scheme. Note that there are many other efficient methods on constructing energy stable schemes for the Cahn-Hilliard equation, such as, convex-splitting method eyre_unconditionally_nodate, invariant energy quadratization method Yang2017, and scalar auxiliary variable method Shen2019. Here we adopt the stabilized explicit method only for simplicity, and it is easy to extend the scheme to the methods mentioned above.

Remark 3.8.

In this paper, we only focus on the decoupling of multiphysics problem rather the decoupling of all variable in each physical problem. The velocity and pressure in step 3 can be further decoupled by using the first order pressure correction scheme shen_decoupled_2015, and we leave it to the interested readers.

It is clear that the scheme given by (12)-(14) is a decoupled, linear scheme. Next we want to show that the scheme is mass-conservative and charge-conservative.

Proposition 3.9.

Let solve (12)-(14) for any , then the scheme is mass-conservative and charge-conservative, namely,


Letting in the first equation of (12), we have mass conservation . Then, we note that for all , there holds

and . Taking , we obtain . ∎

Now, we are in a position to prove the unconditional energy stability of the decoupled scheme as follows, which is analogous to that of the original problem in Theorem 2.3.

Theorem 3.10.

The decoupled scheme is unconditionally energy stable in the sense that the following energy estimate


holds, where


Letting in (12), we have


Next, setting , we get


Then, taking , and using the identity

it yields


By making the summations of (16)-(19), we obtain


Using Young inequality, we derive the right hand side of the first term of (20) has the following estimate,


For the last two term in (20), using the Taylor expansion,

then we have


Plugging (21) and (22) into (20), and dropping the positive term and , we get the required estimate (15). The proof is thus complete. ∎

4 Numerical Experiments

In this section, we present a series of 2D numerical experiments to illustrate the features of the proposed algorithms. The finite element method is implemented on the finite element software FreeFEM developed by hecht_new_2012. For any integer let be the space of polynomials of degree on element , and denote . We employ the Mini-element arnold_stable_1984 to approximate the velocity and pressure

where , is a bubble function. We choose the lowest-order Raviart-Thomas raviart_mixed_1977 element space

to approximate the current density, the discontinuous and piecewise constant finite element space to approximate the electric potential

The phase field and chemical potential are discretized by first order Lagrange finite element space , where

Example 4.11 (Convergence and accuracy).

The first example is used to verify the convergence rates in both time and space. The computational domain is set as , and the external magnetic field is . The physical parameters are given by with terminal time . The right-hand side functions, initial conditions and Dirichlet boundary conditions are chosen such that the given solutions satisfy the system.

Let the approximation errors at the final time be denoted by

First, we test the temporal convergence orders. The analytic solutions are chosen as

Note that the exact solutions are linear or constant in space, the main error comes from the discretization of the time variable. We fix a mesh size with and test the convergence rate with respect to the time step. Then the errors and orders are displayed in Tables 3-3. From these tables, we observe that the errors of all variable decrease as the mesh is refined, with convergence order of , which accords with our theoretical analysis completely.

0.2 1.7227e-4(—) 1.3384e-3(—) 3.4903e-2(—)
0.1 7.5918e-05(1.18) 5.9037e-4(1.18) 1.6788e-2(1.06)
0.05 3.5535e-05(1.10) 2.7649e-4(1.09) 8.1861e-3(1.04)
0.025 1.7206e-05(1.05) 1.3391e-4(1.05) 4.0362e-3(1.02)
0.0125 8.4693e-06(1.02) 6.5926e-05(1.02) 2.0033e-3(1.01)
0.00625 4.2022e-06(1.01) 3.2712e-05(1.01) 9.9790e-4(1.01)
Table 2: Time convergence rates of the scheme for
0.2 6.9951e-3(—) 3.6822e-2(—) 3.68218e-12
0.1 3.6099e-4(0.95) 1.9026e-2(0.95) 1.90247e-12
0.05 1.8227e-3(0.99) 9.6627e-3(0.98) 9.66331e-13
0.025 9.1439e-4(1.00) 4.8671e-3(0.99) 4.86857e-13
0.0125 4.5781e-4(1.00) 2.4422e-3(0.99) 2.44362e-13
0.00625 2.2904e-4(1.00) 1.2233e-3(1.00) 1.22111e-13
Table 3: Time convergence rates of the scheme for
0.2 1.1298e-1(—) 1.8314e-2(—) 1.1941e-1 (—) 6.1937e-3(—)
0.1 5.1171e-2(1.14) 5.9434e-3(1.62) 7.0687e-2(0.77) 3.2132e-3(0.95)
0.05 2.4156e-2(1.08) 2.2442e-3(1.40)