A UNIFIED FRAMEWORK OF THE SAV-ZEC METHOD FOR AMASS-CONSERVED ALLEN

1. Introduction. When magnetized nanoparticles are dispersed in a nonmagnetic liquid carrier (such as an organic solvent or water), a ferromagnetic fluid, also known as a ferrofluid, is created. This one-of-a-kind substance acts like a colloidal solution on the outside, but when an external magnetic field is applied, its macroscopic behavior shows a clear difference from other ferromagnets, with a very unique high magnetic polarization saturation, while the remanent magnetization inside the substance immediately becomes zero when an external magnetic field is withdrawn. Due to such remarkable characteristics, ferrofluids have been increasingly utilized in industrial technology, and biological and medical clinical fields, e.g., in the treatment of microplastics in sewage, cardiovascular diseases, and even cancer (cf. [6,8,17,18,24,25,26,30,31,33,35,50]). However, there are still many issues that need to be studied in depth in regard to the mathematical modeling and numerical simulation of ferrofluids, and one of the core problems is the multiphase flow interface effect of ferromagnetic fluids, which will be considered in this paper.

The current theoretical framework of the monophase ferromagnetic fluid hydrodynamic model, known as the ferrohydrodynamics (FHD), is relatively complete and contains two main generally accepted hydrodynamic models, that of Rosensweig [35,36] and that of Shliomis [47,48]. For detailed discussions of the well-posedness, regularity, and long-term behavior of the two models, we refer the reader to [2,3,4,29,39,49,51] and the references therein. The extension of the monophase FHD model to the two-phase or multiphase case, utilizing the phase-field approach, is more advantageous and is one of the most widely discussed or applied multiphase FHD models [5,19,27,28,60], because the phase-field approach provides a simple modeling tool that can track complex changes at the free interface through an energetic variational approach without imposing complicated interface jump conditions or using the complicated interface tracking methods. Another advantage of using the phase-field method to construct a multiphase flow model is that the resulting model can satisfy the energy dissipation law, on the basis of which some theoretical validation, such as the well-posedness, can be carried out. This feature, therefore, gives rise to a natural requirement on the designed numerical scheme, which is to establish a numerical scheme that satisfies the energy dissipation law at the discrete level. This not only allows for a flexible treatment in dealing with the stiffness problem embedded in the phasefield model, but it also ensures the reliability of the obtained numerical algorithm. Therefore, the crucial question of how to build a numerical scheme with appropriate temporal and spatial accuracy and easy implementation in practice should be carefully answered.

The goal of this paper is two-fold. First, in the existing two-phase phase-field FHD model, it is noted that the governing equation used to trace the free-motion interface is the Cahn-Hilliard type equation, which formally refers to the spatial fourth-order diffusion equation. The fourth-order Cahn-Hilliard equation is relatively more difficult to solve compared to the other signature equation of the phase-field model, the spatial second-order Allen-Cahn equation. Although the Allen-Cahn equation lacks the volume conservation property that the Cahn-Hilliard equation has, one actually has a considerable number of volume conservation methods that can make it possess this property; see [9,11,20,23,38]. Therefore, in this paper, one of these known conservation techniques, namely the conserved Allen-Cahn equation developed in [38], which not only has volume-conservation property but also follows the energy law, is utilized instead of the Cahn-Hilliard equation to build a new phase-field FHD model, thereby reducing the difficulty of solving the fourth-order interface governing equations. However, the conservative Allen-Cahn equation also introduces an additional numerical difficulty, as it uses a volume-preserving approach by including an extra nonlocal term in the equation, which requires special treatment to avoid solving the nonlocal type equation at each time step when designing numerical schemes.

Second, it is worth noting that despite the use of the conservative Allen-Cahn equation, the resulting two-phase FHD model remains a highly nonlinear and coupled complex system. Hence, the development of numerical algorithms for this model remains challenging, particularly when our goal is to construct an efficient numerical algorithm that not only has the ability to maintain the discrete energy law but also possesses properties that are easy to implement in practice. To this end, it is advantageous to recall some of the existing numerical approaches used to solve the two-phase phase-field FHD model. The pioneering work in [28] employed the B79 Cahn-Hilliard phase-field method to model the two-phase ferrofluid and also came up with an energy-stable numerical scheme which is nonlinear, coupled, and first-order in time. The first linearized type numerical scheme for the two-phase FHD model was obtained in [60], where some auxiliary intermediate variables are introduced and combined with the stabilized technique. However, the resulting scheme is still partially coupled and first-order accurate in time. Hence, the remaining key challenge is how to construct a fully decoupled second-order time-accurate energy-stable scheme, while maintaining a linear and easy-to-implement structure. Moreover, it is worth noting that the numerical simulations in [28] and [60] are both limited to 2D. Therefore, the second main goal of this paper is to design a fully discrete numerical scheme with second-order time accuracy, linearity, mass conservation, unconditional energy stability, and fully decoupled structure based on the continuous finite element method and use it to carry out 3D simulations. However, achieving this goal entails addressing several challenges, including the linear/nonlinear couplings, nonlinearities, nonlocal integration, and ensuring mass conservation and unconditional energy stability.

To obtain the desired numerical scheme, it is crucial to effectively deal with the nonlinear terms, which can be classified into three types. The first type has a symmetric positive definite structure, leading to positive diffusion in the energy law, and thus can be linearized and decoupled using a symmetric implicit-explicit format, while preserving unconditional stability. The second type is the nonlinear potential in the phase-field equation, contributing to the system energy. The third type of nonlinear terms, on the other hand, does not contribute any energy in the energy law. Recently, the energy quadratization SAV method [41,42,43] was proposed in designing linear and energy stable schemes for the phase-field problem, and the ZEC decoupling approach [52,53,54,55,57,58] was invented to handle the zero-energy-contribution nonlinear terms to achieve linear and stable algorithms. By directly combining the SAV method for the second type of nonlinearity with the ZEC method for the third type of nonlinearity, the desired scheme may be obtained. However, the introduction of two scalar variables via two ordinary differential equations (ODEs)-one for SAV and one for ZEC-and coupling them with the original system can indeed add complexity to the PDE system. Moreover, to incorporate the two scalar variables introduced by the SAV and ZEC methods, the unknowns need to be split twice, which further increases the complexity of decoupled implementation; cf. [54,55] for other relatively simple two-phase flow models.

Therefore, in this paper, we propose a unified framework that incorporates the SAV and ZEC methods together to handle both the second and the third kinds of nonlinear terms simultaneously, representing another main contribution of our work. The key to unifying the SAV and ZEC ideas is to incorporate the second and third kinds of nonlinearities together into a special designed ODE for a nonlocal scalar auxiliary variable. With this unified framework, only one scalar variable and one ODE are introduced, so the unknowns only need to be decomposed once in the algorithm implementation. This unified SAV-ZEC combined method framework is also applicable to other two-phase fluid flow systems in designing stable and efficient schemes.

We reformulate the Allen-Cahn equation and magnetic potential equation to eliminate the undesired linear couplings, adopt the second-order pressure correction method [16,34,40], and transform the saddle point system into elliptic equations. Furthermore, we use a special test function in the chemical potential equation to transform the implicit nonlocal integration into an explicit computation. By applying implicit treatments to the nonlocal integration and linear terms and, most importantly, with the aid of the newly introduced scalar variable, we obtain the mass conservation and unconditional energy stability.

The aforementioned numerical techniques, in combination with the finite element method for spatial discretization, allow us to obtain an efficient fully discrete numerical scheme that possesses the properties of full decoupling, second-order accuracy in time, unconditional energy stability, mass conservation, and linearity. It is important to highlight that all the reformulations of the original model, introduction of the new auxiliary variable, and construction of the corresponding ODE are aimed towards the ultimate goal of developing such a numerical scheme. We also demonstrate the well-posedness of the proposed scheme and rigorously prove its unconditional energy stability and mass conservation. The scheme is highly efficient, as it splits the nonlinear coupled saddle point system into a series of independent elliptic problems, and has been verified through various numerical examples, including accuracy tests, energy stability verification, and some 2D/3D benchmark simulations, exhibiting signature "ferrofluid hedgehog" phenomena of two-phase ferrofluid drops.

The rest is organized as follows. In section 2, we develop a two-phase phasefield FHD model using the nonlocal conserved Allen-Cahn dynamics and derive its energy dissipation law. In section 3, we introduce the unified framework of SAV and ZEC approaches to construct a fully discrete finite element numerical scheme and rigorously prove its unconditional energy stability. The section also provides a detailed explanation of the decoupled type of implementation for each step. Section 4 presents various numerical experiments in 2D and 3D to demonstrate the accuracy and efficiency of the proposed numerical scheme. Finally, section 5 offers concluding remarks.

2. Conserved Allen-Cahn type FHD model. We consider a fluid flow system confined in a bounded convex polygon/polyhedron domain Ω ⊂ R d with d = 2 or 3. The well-established monophase Shliomis model for a viscous, homogeneous ferrofluid flow system reads as [47,48]

in which the unknown physical variables are the velocity field u, pressure p, magnetization field m, effective magnetic field h(:= ∇ϕ), and magnetic potential ϕ. Besides, h a is an applied smooth harmonic magnetic field (∇× h a = 0, ∇ • h a = 0), ν is the kinematic fluid viscosity, χ 0 is magnetic susceptibility, µ is permeability of free space, τ is relaxation time constant, β = 1 6νϑ , ϑ is volume fraction of dispersed solid phase, n ∂Ω is the outward normal on the boundary ∂Ω, and the term (m • ∇)h is the socalled Kelvin force. Note that the no flow boundary condition of the fluid prevents the necessity of using boundary conditions for the magnetization equations; see [32] and the references therein.

To extend the monophase model (2.1) to the two-phase case of an immiscible mixture consisting of the ferrofluid and non-ferromagnetic viscous medium, the framework of the phase-field approach requires a labeling variable Φ, which is defined as Φ(t, x) = 1 ferrofluid phase, 0 non-ferromagnetic viscous fluid, with a thin smooth transition layer of width O( ) connecting the two fluid components. Thus, the interface of the mixture can be traced by the level set Γ = {x : Φ(t, x) = 1/2}. Using the conserved Allen-Cahn dynamics [38], the evolution of the phase-field variable follows the following governing equation:

where M > 0 is the mobility parameter, W is the chemical potential, λ accounts as the surface tension parameter, and 2 is the Ginzburg-Landau double-well potential. The nonlocal term in (2.2) is used to eliminate the total variance of mass (or volume).

By coupling the monophase Shiliomis model (2.1) with the conserved Allen-Cahn system (2.2), and assuming that the fluid is incompressible and follows the generalized Fick's law, that is, the mass flux is proportional to the gradient of the chemical potential, we obtain the two-phase FHD model that reads as

Here ν(Φ) = ν w + (ν f -ν w ) 1 1+e -(2Φ-1)/ , ν f and ν w are viscosities for the ferrofluid flow and non-ferromagnetic viscous medium, respectively, χ(Φ) = χ 0 1 1+e -(2Φ-1)/ (one can also use χ(Φ) = Φ 2 χ 0 for simplicity), D(u) = 1 2 (∇u + (∇u) ), and the term Φ∇W is the induced elastic stress by the mixing energy [28,44,45,46,56].

We introduce some function spaces. For two vector functions v, w, we denote the L 2 inner product as (v, w) = Ω v • wdx and L 2 norm w 2 = (w, w). We use H 1 (Ω) to denote the usual Sobolev space and define

Then the system (2.3)-(2.10) admits the following energy dissipation law and has the mass (volume) conservation property as follows.

Theorem 2.1. The system (2.3)-(2.10) possesses the following energy estimate:

11)

where h b = (h a ) t . If the applied magnetic field h a = 0, there holds the energy dissipative law

Moreover, the mass conservation property holds as Ω Φ(t)dx = Ω Φ 0 dx for t ∈ (0, T ].

Proof. Taking the L 2 inner product of (2.3) with W , (2.4) with Φ t , (2.5) with u, (2.7) with µh, (2.8) with µ τ ϕ, respectively, and applying integration by parts, noticing h = ∇ϕ, we have

Taking temporal derivative of (2.8) and taking the L 2 inner product of it with µϕ, we get

Furthermore, taking the L 2 inner product of (2.7) with µ χ0 m, and applying the identities

By using ∇× h = 0 and integration by parts, we derive

Then, by combining (2.13)-(2.19), and using (2.20), we derive

(2.21)

B83

We estimate the terms on the right-hand side by

which completes the proof of estimate (2.11) and also implies the energy dissipative law (2.12) when assuming h a = 0. By taking the L 2 inner product of (2.3) with 1, we obtain d dt Ω Φdx = 0, which implies the mass conservation property. 3. Numerical scheme. In this section, we aim to construct the linear, unconditionally energy-stable, second-order accurate in time, mass conserved, and fully decoupled type numerical algorithms for the system (2.3)-(2.10). In general, there are four difficulties to be overcome in establishing such a numerical format, including (i) how to linearize the nonlinear terms; (ii) how to decouple the couplings among Φ, W, u, p, ϕ, m; (iii) how to discretize the nonlocal term in (2.3); and (iv) how to preserve the mass conservation and energy stability unconditionally at the discrete level.

3.1. Equivalent reformulation. This subsection is the preparation phase. Namely, before we proceed to establish the numerical scheme, we will convert the system (2.3)-(2.10) into an equivalent form in several steps using equation deformation, auxiliary variables, or other means to facilitate the design of numerical algorithms. Considering the complexity of the original system, this effort is worthwhile to obtain the desired type numerical scheme.

3.1.1. Magnetostatic equation. The magnetostatic equation (2.8) poses two difficulties in designing numerical methods. One is that the linear coupling relation between the magnetic field ϕ and the magnetization field m needs to be implicitly discretized at the same time for the energy stability. This makes it very difficult to reach full decoupling. The other difficulty is that the energy estimation of h or (∇ϕ) needs a hybrid test, which complicates the design of the discrete space. One way to reformulate (2.8) to tackle these two numerical issues was proposed in [62] and will be briefly recalled as follows.

For any ψ ∈ H 1 (Ω) ∩ L 2 0 (Ω), by testing 1 τ ψ on (2.8), we obtain

We further take the time derivative of (2.8) and formulate the obtained equation in the weak form to get

By taking the L 2 inner product of the magnetization equation (2.7) with ∇ψ, we derive

By summing up (3.1)-( 3.3), we arrive at a weak formulation of the magnetostatic equation, which reads as

We replace (2.8) with (3.4) and note that the linear coupling of ϕ and m in (3.4) does disappear, and the energy estimate of h (:= ∇ϕ) could be naturally derived in (3.4), thus avoiding the two numerical difficulties described above. See more explanations in the following remark.

Remark 3.1. If the linear coupling of ϕ and m in (2.8) is explicitly decoupled, it would cause numerical instability since the explicit treatment of ∇ • m could not be balanced by any other term when the energy stability is derived for the discrete scheme. And the implicit treatment of ∇•m will present a coupled type scheme which is not the aim of this paper. On the other hand, in (3.4), the linear coupling of ϕ and m disappears, thus providing an opportunity to construct decoupled algorithms.

Besides, to obtain the energy estimate, one needs to take the test function h(:= ∇ϕ) in the magnetization equation (2.7); see the proof of Theorem 2.1. This will require the match of the discrete spaces of m and ϕ in the schemes. For instance, ∇Ψ h ⊂ N h ; here Ψ h and N h are discrete spaces for ϕ and m. This will increase the complexity in the spatial discretizations of the algorithm design. On the other hand, in (3.4), the energy estimate can be obtained by taking ψ = ϕ; see the proof in Theorem 3.1.

Kelvin force.

The Kelvin force term (m • ∇)h in (2.5) also needs some special treatment. When we consider the weak form of (2.5) and take the L 2 inner product of (2.5) with a test function v ∈ H 1 0 (Ω), the weak form of the Kelvin force term becomes µ((m • ∇)h, v). Since h = ∇ϕ, this term involves the second-order derivative of ϕ, which is not feasible for the continuous finite element method; cf. [28].

We overcome this issue by rewriting the Kelvin force as

where we use the fact that ∇× h = 0 and integration by parts; cf. [62]. In (3.5), we note that there are only first-order spatial derivatives. Then the weak form of (2.5) reads as

3.1.3. Nonlinear couplings. Based on the energy law in (2.11), the regularity of the Allen-Cahn system [13], and the monophase Shliomis model [2], we assume

provided smooth initial data and a finite time T . With the aid of (3.4) (to replace (2.8)) and (3.6) (to replace the weak form of (2.5)), the weak form of the system (2.3)-(2.10) is to find (Φ, W, u, p, ϕ, m) satisfying (3.7), such that for all (X, Y, v, q, ψ,

We now study the large number of coupled nonlinear terms present in the system (3.8)-(3.13), which pose significant difficulties in designing the desired type numerical scheme. We can see that there are three types of nonlinear terms as follows.

• The first kind is the symmetric term β(m × ∇ϕ, m × ∇ψ) in (3.12) that builds into the positive diffusion in the energy law, which can be discretized by the symmetric implicit-explicit combination method. • The second kind is the nonlinear potential f (Φ) in (3.9) that builds into the system energy in the energy law, which can be discretized by the linear energy quadratization type approaches, e.g., the so-called SAV method [41,42,43]. • The third kind is the remaining 11 nonlinear terms, and we find that these nonlinear terms contribute zero energy; namely, when treated separately or partially combined, they satisfy the following property:

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 03/31/24 to 129.252.33.201 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy These equalities are obtained in the derivation process of the energy stability if we set X = W , v = u, ψ = µϕ, and n = µ χ0 m in (3.8), (3.10), (3.12), and (3.13), respectively. Thanks to the ZEC approach for the other coupled phase-field type models (see [52,53,54,55,61]), these equations imply that we can extend the ZEC decoupling method to treat these terms. Instead of the direct way of using the SAV method for the second kind of nonlinear terms and the ZEC approach for the third kind of nonlinear terms, respectively, in this paper, we develop a unified framework for the SAV method and ZEC approach to process the two kinds of nonlinear terms together. The motivation for this idea is to facilitate the design of the scheme and improve computational efficiency. Our key strategy to unify the SAV method and ZEC approach is to incorporate the two kinds of nonlinear terms in an appropriate manner to define a special ODE for a nonlocal scalar auxiliary variable, which will be presented as follows.

We denote S(t) := Ω F (Φ)dx + B, where B is a shifting constant such that the radicand is positive, and define a scalar variable r(t) through the following ODE:

(Φ∇W, u)

We notice the weighed parameter 1 2λS multiplying the third kind of nonlinear terms is essential to derive the energy law. From (3.14), it can be seen that the above equation is equivalent to

It is easy to derive that r(t) = Ω F (Φ)dx + B = S(t) after integrating the first equation in (3.16) and applying the initial condition of r(0).

Using the scalar variable r and its ODE, we continue to transform the system (3.8)-(3.13) into another equivalent form: given the initial data (2.10) and r(0), find (Φ, W, u, p, ϕ, m) satisfying the regularity requirements in (3.7) and r ∈ R, such that for all (X, Y, v, q, ψ,

Note that in (3.17)-(3.22), we multiply the second and third kinds of nonlinear terms by r S . It is important to emphasize that this modification does not change the system from a PDE point of view, as r = S. Therefore, the system (3.17)-(3.23) is equivalent to the system (3.8)-(3.13). Meanwhile, since the system (3.17)- (3.23) is the equivalent weak form of the original PDE system, it is clear that it complies with the energy dissipation law. As a result, we do not provide a separate proof of the energy dissipation law for this new system, as it is similar to that of the energy stability of the numerical scheme (see Theorem 3.1). Remark 3.2. If one utilizes the SAV method for the second kind of nonlinearity and extends the ZEC approach to the third kind of nonlinearities, respectively, two scalar variables and two ODEs will be introduced. The weak form (3.8)-(3.13) coupled with the scalar variables and ODEs would result in a more complicated system. Moreover, in the decoupled implementation, as described in section 3.3, the unknowns need to be split twice in terms of the two introduced scalar variables, which also leads to more problems to solve, thereby reducing computational efficiency to some extent; cf. [54,55] for simpler two-phase fluid flows. On the other hand, in our unified framework of the SAV method and ZEC approach, only one scalar variable and one ODE are introduced, and the unknowns need to be split only once in decoupled implementation, which not only alleviates the complexity of the PDE system but also reduces computational costs. Remark 3.3. Although the system (3.17)-(3.23) is equivalent to the original system (2.3)-(2.10) in the weak form, formally, it seems to be more complex. However, it is worth noting that the format of (3.17)-(3.23) is more "algorithm-friendly" than the original system. This formulation allows us to discretize coupled nonlinear terms in a simpler way (to be given in the next subsection), i.e., with a decoupled structure while guaranteeing unconditional energy stability.

3.2. Construction of numerical scheme. In this subsection, we construct the numerical scheme for solving the equivalent system (3.17)- (3.23). Letting N > 0 denote the total number of time steps, we define the uniform time step size as δt = [ T N ] and set t n = nδt. We introduce several conforming finite element spaces for spatial discretization as follows:

The pair of spaces (V h , Q h ) needs to satisfy the inf-sup condition [14]:

∇v for all q ∈ Q h , where the constant β 0 only depends on Ω. Some wellknown inf-sup stable pairs (V h , Q h ) are discussed in [14]. For simplicity, we denote

The numerical scheme of system (3.17)-(3.23) reads as follows: find

Some remarks are in order.

Remark 3.4. We explain the strategy behind developing the above scheme. We discretize the time derivatives by the two-step BDF2 format. The second-order pressure projection method [16,40] is used to decouple the linear coupling of the velocity field u and the pressure p in the fluid momentum equation. The nonlinear coupling term β(m×∇ϕ, m×∇ψ) with the symmetric positive definite structure is discretized by a symmetric implicit-explicit format in (3.30), while all other nonlinear terms multiplied by r are treated in explicit extrapolation, and r is discretized implicitly. These particular discretizations are due to the fact that our goal is to construct a linear, decoupled, and energy-stable scheme. The nonlocal term in (3.25) is discretized in an implicit manner to ensure mass conservation. For ODE (3.23), some subtle combinations of implicit and explicit discretization are applied to achieve unconditional energy stability in order to maintain correlation with the nonlinear coupling term discretized to form the decoupling structure (as shown in subsection 3.3).

The scheme (3.25)-(3.32) may appear to be a coupled version, but in fact, due to the explicit approach used, all direct coupling between variables is eliminated, and instead, all variables are coupled to r. Therefore, with a method that can decouple the r-coupling, the decoupled structure is achieved. In addition, computing the implicit nonlocal integral Ω W n+1 dx needs considerable computational cost. In subsection 3.3, we will propose an effective implementation to address these two issues. Remark 3.5. It can be verified that the final velocity field u n+1 in above scheme (3.25)-(3.32) satisfies the following weakly discrete divergence-free condition:

For simplicity, we denote W n+1 = 1

|Ω| Ω W n+1 dx. The scheme (3.25)-(3.32) holds the energy law unconditionally and the mass conservation property, shown as follows.

Theorem 3.1. The scheme (3.25)-(3.32) is unconditionally energy stable in the sense that

if the imposed magnetic field h a = 0, there holds the energy dissipative law unconditionally,

where

Moreover, the mass conservation property holds as Ω Φ n+1 dx = Ω Φ 0 dx.

Proof. By taking

From (3.29), we derive the orthogonal identity:

Using (3.39) and (3.33), we derive

(3.41)

We rewrite (3.29) as

By taking the L 2 inner product of the above equation with itself, using (3.33) and (3.41), we derive

(3.43)

The combination of (3.40) and (3.43) gives

By combining (3.38) with (3.44), we deduce

By taking ψ = µϕ n+1 in (3.30), we derive

Thus, by combining (3.49), (3.50), and (3.51), we obtain

After dropping several unnecessary positive terms on the left-hand side of (3.52), we derive (3.34). Meanwhile, (3.35) can be derived by simply setting h a = 0 in (3.34).

For the mass conservation property, since it is a cumulative process, we must prove that the first step to calculate Φ 1 is also mass-conserved, which is very easy to show since the first-order version of (3.25) is simply to use the backward Euler d t Φ 1 to replace D t Φ 1 for the time marching, and u = u 0 , Φ * = Φ 0 , namely,

By setting X = 1, we get Ω Φ 1 dx = Ω Φ 0 dx. Then, by taking X = 1 in (3.25), we have (D t Φ n+1 , 1) = 0, which yields 3 Ω Φ n+1 dx = 4 Ω Φ n dx -Ω Φ n-1 dx for n = 1, 2, . . . , N -1. Therefore, we derive Ω Φ n+1 dx = Ω Φ n dx for n = 0, 1, . . . , N -1. The proof is completed.

Decoupled implementation.

In this subsection, we present an efficient implementation method of the proposed scheme (3.25)-(3.32).

The nonlocal term.

We expect to avoid any type of iterations involving nonlocal terms. Hence, we first deal with the nonlocal term Ω W n+1 dx in (3.25). Note that if we take Y = 1 in (3.26), we get

Then, by taking Y = X in (3.26), we also get

Then, by applying (3.53) and (3.54), we can transform (3.25) into the following form: solve Φ n+1 ∈ Y h such that for all X ∈ Y h , there holds

(3.55)

It can be seen that the nonlocal and nonlinear terms in (3.55) are explicitly discretized and they involve only previous time steps, so in fact (3.55) is an elliptic equation with constant coefficients.

The splitting technique.

Then, to obtain the fully decoupled type calculation, we split the unknown variables using the nonlocal variable r as follows:

b are the unknown variables for the split.

Using the split form in (3.56) to replace variables in (3.55), and according to r n+1 , we can decompose the resulting form into two substeps as follows.

• Step 1: find Φ n+1 a ∈ Y h such that for all X ∈ Y h , there holds

Using the split form in (3.56) to replace variables in (3.26), and according to r n+1 , we can decompose the resulting form into two substeps as follows.

•

Using the split form in (3.56) to replace variables in (3.27), and according to r n+1 , we can decompose the resulting form into two substeps as follows.

• Step 5: find ũn+1 a ∈ V h such that for all v ∈ V h , there holds

(3.61)

Using the split form in (3.56) to replace variables in (3.30), and according to r n+1 , we can decompose the resulting form into two substeps as follows.

• Step 7: find ϕ n+1 a ∈ Ψ h such that for all ψ ∈ Ψ h , there holds

Using the split form in (3.56) to replace variables in (3.31), and according to r n+1 , we can decompose the resulting form into two substeps as follows.

• Step 9: find m n+1 a ∈ N h such that for all n ∈ N h , there holds

Using the split form in (3.56) to replace variables in (3.32), we can get a linear algebraic equation for r n+1 that reads as follows.

• Step 11: find r n+1 by

, where for k=a, b, Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. Downloaded 03/31/24 to 129.252.33.201 . Redistribution subject to SIAM license or copyright; see https://epubs.siam.org/terms-privacy

With Steps 1-10 above, we get all variables with subscripts a, b, and also r n+1 from (3.67). Hence, by using the split form given in (3.56), we get the unknown variables Φ n+1 , W n+1 , ũn+1 , ϕ n+1 , and m n+1 . The final unknown variables p n+1 and u n+1 are obtained from Step 12 as follows.

•

Step 12: we update p n+1 from (3.28) and u n+1 from (3.29).

As can be seen from Steps 1-12, the implementation of the scheme (3.25)-(3.32) is completely decoupled. In addition, Steps 1-12 require solving only a few linearly independent elliptic problems.

So far, we have proposed a linear, fully decoupled, second-order in time, massconserved, and unconditionally energy-stable scheme. The final issue is to determine the unique solvability of the equations in (3.57)-(3.67), shown as follows. Proof. The well-posedness of problems (3.57)-(3.66) in Steps 1-10 can be proved by the Lax-Milgram theorem [7], where some of them also use inverse inequality and Korn's inequality [7]. We omit the details here.

We show the unique solvability (3.67) as follows. By taking

Then, by combining (3.71) with (3.72), we get ).

By using the Cauchy-Schwarz inequality, we estimate the last term above as

Therefore, we obtain

Thus, we have 3 -η b = 0, which implies the well-posedness of the linear algebraic equation (3.67).

4. Numerical simulations. In this section, we implement a series of numerical simulations to verify the accuracy and stability of our scheme and show some benchmark simulations of ferrofluids. For spatial discretizations, the first-order (linear) polynomials are used for Y h , Q h , and N h , and second-order (quadratic) polynomials are applied for V h and Ψ h .

We denote e w = w(t n , x)-w n as the approximation error at the recorded moment t n and " " the relation of a ≤ Cb for some constant C. From the chosen finite element spaces, the optimal error orders of the scheme (3.25)- (3.32) To observe the convergence orders, we set δt = 1 2 h, and from the expected optimal error estimates (4.1), there hold

We show the accuracy tests at t = 0.5 and t = 1.0 in Figure 4.1, where the L 2 errors of Φ, u, h, and m all display second-order accuracy, the H 1 error of Φ has first-order accuracy, and the L 2 error of p and the H 1 error of u do not present second-order accuracy, but slightly higher than first order. These convergence results are consistent with the theoretical expectations given in (4.3).

Energy stability.

In this test, we carry out a benchmark coarsening effects simulation (cf. [10,59]) to verify the energy stability of the scheme (3.25)- (3.32). The computational domain is set as Ω = [0, 2π] 2 , and the initial conditions are set as

), We set the time step δt = 1 1000 and plot the profiles of Φ at various times in Figure 4.2(a). It can be seen that the coarsening effect makes the small circle absorbed by the large circle. At around t = 0.5, the small circle disappears completely. We also verify the energy stability of our developed scheme. To get a more intuitive impression of the stability, we compare our scheme with a second-order accurate implicit-explicit scheme, where all nonlinear terms are treated explicitly except for β(m × ∇ϕ n+1 , m × ∇ψ), and all linear terms are treated implicitly. The implicit-explicit type scheme does not guarantee any energy stability; however, it has been widely used for other types of phase-field models due to its ease of implementation (see [21]). Using the implicit-explicit scheme, we plot the total free energy E(Φ n+1 , u n+1 , h n+1 , m n+1 ), defined in Theorem 2.1, in Figure 4.2(b). It can be seen, only when δt ≤ 1 2000 , that the implicit-explicit scheme presents the energy stability, while the energy blows up when δt ≥ 1 1000 . For comparison, using our developed scheme, we plot the discrete energy E n+1 h , defined in Theorem 3.1, in Figure 4.2(b) by varying time steps δt = 1 2 i , i = 2, . . . , 10, where the energy curve indicated by the dashed line is the total energy calculated from the implicit-explicit scheme with δt = 1 2000 that is used as the reference solution. It can be seen that the energy profile calculated using our scheme is stable over all tested time steps and approaches the reference energy as the time step is refined. These numerical results confirm the energy stability of our scheme stated in Theorem 3.1.

The thickness of the diffusive interface is proportional to the parameter , and we also study the interplay between the parameter and mesh size h. We choose = 0.015, M = 10, ν f = 0.1, ν w = 0.05, µ = 10, τ = 0.1, β = 10, χ 0 = 1, λ = 1, and δt = 1 2000 , and we load the initial value in (4.4) but with u 0 = 0, m 0 = 0. The profiles of Φ, computed by h = 1 64 , 1 128 , and 1 256 , at various moments are illustrated in Figure 4.3, which shows that for a fixed , the finer grid could produce a more accurate result. 4.3. 2D ferrofluid hedgehog. When a pool of ferrofluid is subjected to gravity and an external magnetic field h a pointing upward, experiments show that the fluid interface becomes unstable and a regular pattern of peaks and troughs emerges. This phenomenon results from the competition between three forces: gravity, surface tension, and the magnetic field, where gravity and surface tension try to return the surface to a flat state, but the strength of the magnetic field tends to create a vertical surface. In this subsection, we aim to simulate this phenomenon of the two-phase fluid system, i.e., a mixture of a ferrofluid and a nonferromagnetic ambient viscous fluid with different viscosities and almost matching densities, under a nonuniform applied magnetic field; see also [15,27,28]. In [12,28], some analytical results of the interpeak distance based on linear stability analysis with small magnetic susceptibility χ 0 are provided, which is l p = 2π σ ∆ρg , where l p denotes the distance between peaks, σ is the surface tension coefficient, g = |g| is the magnitude of gravity, and ∆ρ is the jump of the density across the interface.

We add the gravity force f g in the fluid momentum equation (2.5) by using the Boussinesq approximation, i.e., f g = (1 + rg 1+e 1-2Φ )g, where r g is a positive constant that depends on the fluid density, and |g| stands for the magnitude of gravity. The computational domain is set as Ω = [0, 1] × [0, 0.2]. The initial shape of the ferromagnetic fluid as a semicircular droplet located in the bottom plane of the computed domain reads as Φ| t=0 = 0.5-0.5 tanh( √

), where x 1 = 0.5, y 1 = -0.01, r 1 = 0.2. All other variables are set as zero, u| t=0 = 0, p| t=0 = 0, m| t=0 = 0. The applied magnetic field h a is generated by a linear combination of dipoles as follows:

where |d| = 1 indicates the direction of the dipole, and x s is the dipole's position. It is easy to calculate that h a is a harmonic field (i.e., ∇× h a = 0, ∇ • h a = 0); cf. [28]. To generate a nonuniform applied magnetic field, we set h a = 5 s=1 α s ∇φ s (x) by placing five dipoles below the container Ω at close range. The positions x s of dipoles are (0.4, -1), (0.45, -1), (0.5, -1), (0.55, -1), and (0.6, -1), and the directions d of the five dipoles are all (0, 1). The intensities α s (t) = In the next few simulations, we study the effect of the gravity magnitude on the number of peaks from a qualitative point of view. First, by setting g = (0, -30000), we plot the obtained profiles of the phase-field variable Φ at various times in Figure 4.4. It can be seen that the droplet slowly becomes flat (from t = 0 to t = 0.4). Starting from t = 0.5, instability begins to appear on the droplet surface as the gradually increasing applied magnetic field exceeds gravity. After t = 0.8, a stable and regular hedgehog pattern containing five peaks is formed. Second, by increasing the gravity magnitude to g = (0, -60000), the computed profiles of Φ at various times are shown in Figure 4.5. It starts with three peaks (t = 0.7), which soon increase to five peaks (t = 0.8), and finally form seven peaks (t = 1). Third, we continue to increase the gravity to g = (0, -90000) and snapshots of Φ at various times are shown in Figure 4.6. Four peaks initially appear (t = 0.8), which quickly become six peaks at t = 0.9s and form the hedgehog pattern of eight peaks after t = 1.3. We further plot the velocity field u, pressure p, magnetization field m, and effective magnetic field h for the third simulation at t = 1.5 in Figure 4.7.

From the above numerical simulations of the benchmark problem of "ferrofluid hedgehog," we conclude that the stronger the magnitude of gravity, the more peaks appear at the ferrofluid interface, which is qualitatively consistent with the theoretical formula given in [12,28].

4.4. 3D ferromagnetic droplet. In this example, we simulate the deformation of a 3D ferrofluid droplet suspended in a viscous medium under a uniformly applied magnetic field; cf. [1,5,19,22,37]. Due to the competition between surface tension, which favors a spherical shape, and the magnetic interfacial force, which creates a shape parallel to the field, the droplet undergoes deformation.

We set the computed domain as Ω = [0.3, 0.7]×[0.3, 0.7]×[0, 1]. A uniform applied magnetic field h a is generated by (4.5) by placing 25 dipoles far below the domain, with h a = 25 s=1 α s ∇φ s (x), where the directions d of all dipoles are (0, 0, 1), the positions x s of dipoles are (-0.5 + 0.5i, -0.5, -15), (-0.5 + 0.5i, 0, -15), (-0.5 + 0.5i, 0.5, -15), (-0.5 + 0.5i, 1, -15), and (-0.5 + 0.5i, 1.5, -15) for i = 0, 1, 2, 3, 4. The intensity is fixed as α s = 1000 for all dipoles. The initial conditions of Φ are set as Φ| t=0 = 0.5 -0.5 tanh (x -x 1 ) 2 + (y -y 1 ) 2 + (z -z 1 ) 2 -r 1 1.2 , . We plot the profiles of the phasefield variable Φ at various times in Figure 4.8. As the magnetic force governs the surface tension, the droplet is elongated in the direction of the applied magnetic field. The simulation is consistent with the results presented in [1,5,19,22,37]. We further plot the velocity field u and the magnetization field m at t = 0.03 in 5. Concluding remarks. In this article, we present a mass-conserved Allen-Cahn type phase-field model of the two-phase ferrofluid flow and construct an efficient numerical algorithm for solving the model. By developing a unified framework of the SAV method and ZEC approach to discretize the nonlinear couplings for linearization and decoupling, and by eliminating linear couplings through the combinations of equations and the projection method, we have constructed a very efficient numerical scheme for the resulting system. The scheme is linear, second-order accurate in time, fully decoupled, mass-conserved, and unconditionally energy stable.

Its implementation is also efficient and requires solving several independent elliptic problems per time step. We prove the scheme's mass conservation, unconditional energy stability, and well-posedness and carry out a number of numerical simulations to verify the effectiveness of the developed model and the scheme's effectiveness and robustness. Furthermore, it is important to note that the proposed unified framework for the SAV-ZEC method is not restricted to the two-phase ferrofluid flow model studied in this article but can also be extended to other coupled-type phase-field models involving fluid flow or other applied fields.