Starting from Du et. al.'s work in [1,2], the phase-field method or diffuse interface method has been used to simulate the shape transitions of lipid vesicles under various conditions, see [3][4][5][6][7]. The basic idea of the phase-field method is to introduce one or multiple labeling functions (called phase-field variables) to label the internal and external components surrounded by the vesicle membrane. Then, by using the elastic bending energy represented by the phase-field variables to replace the average curvature of the membrane surface, and minimizing the total free energy in a certain metric space, the so-called phase-field elastic bending energy model (EBE, for short) for lipid vesicles has arrived.
The derivation of the classical EBE model normally uses the so-called Allen-Cahn (L 2 -gradient flow) type dynamical system. Because the free energy formally contains a second-order Laplacian term, the governing system consists of one (for single-phase-field EBE model) or multiple coupled (for multi-phase-field EBE model) highly nonlinear fourth-order equation(s). To obtain the conservation of global volume and surface area, it adopts a penalty method, that is, adding two penalty energy potentials to the total energy, so that the volume and surface area can be approximately conserved. However, the penalty method can cause a problem, that is, a larger penalty parameter can better conserve the volume and surface area, but the system becomes stiffer, resulting in a small time step required. Certainly, some numerical techniques can be used to achieve the conservation of volume and surface area through numerical iterations, such as the Lagrangian multiplier method [3,8]. However, it leads to some essential difficulties to develop a numerical scheme with provable unconditional energy stability. Therefore, a natural question arises, that is, how to achieve precise conservation of volume and surface area while deriving the partial differential equation (PDE) system.
Remarkably, the formula for expressing the global volume of the vesicle using the phase-field variable is much simpler than the surface area (see (2.2) and (2.12), and also their definitions in [1][2][3][4][5][6][7]). Therefore, there are indeed more ways to achieve accurate volume conservation than to achieve accurate surface area conservation, which is still an unresolved problem so far. For example, one can reformulate the model using the well-known volumeconserved Cahn-Hilliard dynamics, see [9]. However, it brings a disadvantage that the generated system has two more orders than the Allen-Cahn system, so it is relatively difficult to solve. If the low-order advantage of the Allen-Cahn equation is expected to be retained, inspired by a widely-known conservative Allen-Cahn equation developed by Rubinstein and Sternberg in [10] (see also the application in [11][12][13][14][15][16][17][18][19][20] for different models), we can add a simple nonlocal term to the PDE system, so that the precise conservation of volume can be achieved from the perspective of modeling instead of using advanced numerical techniques. Hence, the first goal of this paper is to use the conservative Allen-Cahn dynamics to reconstruct the multiple-phase-field EBE model thereby conserving the volume in a precise manner.
Next, since both the inside and the outside of the vesicles are fluids, we consider the numerical approximation of the new volume-conserved multi-phase-field vesicle model coupled with the incompressible flow field, i.e., the Naiver-Stokes equations. The focus is to construct an efficient and accurate scheme to solve this highly nonlinear coupling system, that is, to establish a second-order time-accurate, energy-stable, fully-decoupled, and linear scheme. We recall that, regarding algorithm development for the single-phase-field EBE model without coupling the flow field, some available numerical schemes can successfully achieve unconditional energy stability, for example, the Invariant Energy Quadratization (IEQ) method [21], Scalar Auxiliary Variable (SAV) method [22], linear stabilization method [23], nonlinear functional derivative method [9], Exponential Time Differencing method (ETD) [8], etc. Therefore, it may be considered that when dealing with the coupling between the vesicles and the flow field, an effective numerical method can be easily developed by simply combining the above-mentioned numerical methods for the EBE model and the numerical method of solving the Navier-Stokes equation.
However, the fact is quite the opposite. It is not an easy task to develop a scheme that can simultaneously have energy stability, second-order time accuracy, and full decoupling structure. So far, as far as the author knows, among various temporal discretization approach for Navier-Stokes coupled phase-field models (cf. [9,[24][25][26][27][28][29][30][31][32][33]), the only available energy-stable fully-decoupled scheme was developed in [9,23], in which the way to achieve the full decoupling is to add a stabilization term to the explicit advection velocity term. However, the decoupling type scheme in [9,23] is only first-order accurate in time, and it requires more calculations at each time step because the phase field equation needs to be solved with variable coefficients. Moreover, it seems to be very challenging to upgrade the stabilization method to the second-order version. Therefore, how to establish a fully-decoupled numerical scheme with second-order time accuracy is still an unsolved problem. The main difficulty lies in how to discretize the coupled nonlinear terms between the phase-field variable and the flow field (e.g. the advection and stress). At the same time, it is even more unfortunate that these kinds of terms exist in almost all flow-coupled phase-field type models (for example, the phase-field surfactant model [34][35][36], the ternary flow-coupled phasefield model [37,38], the two-phase complex fluid model [31,39,40], etc.). This makes how to establish an ideal fully-decoupled scheme has become a universal difficulty.
Therefore, the second purpose of this paper is to overcome this challenge and develop a novel full decoupling method for the flow-coupled volume-conserved multi-phase-field EBE model. To this end, we note that advection and stress terms satisfy a so-called "zero-energy-contribution" feature. That is, when deducing the energy law, after applying the inner products of some appropriate functions, the results of these two terms will completely cancel out. Thus, using this property, we introduce a nonlocal variable and design an ordinary differential equation (ODE) containing the inner products of the advection and stress with some specific functions. This ODE is trivial at the continuous level because all the terms in it are zero. But after discretization, it can help eliminate all the troublesome nonlinear terms to obtain unconditional energy stability. Meanwhile, the introduction of the nonlocal variable can decompose each discrete equation into multiple sub-equations that can be solved independently, thereby obtaining a fully-decoupled structure.
By combining this novel method with the existing proven effective methods (including the projection method for the Navier-Stokes equations, and the SAV method that linearizes the nonlinear energy potential), we finally arrive at an unconditionally energy stable, linear, fully-decoupled, second-order time-accurate scheme. The implementation of this scheme is very simple. It only needs to solve a few linear independent equations with constant coefficients at each time step, which means that computation is very efficient. We also give a rigorous proof of the solvability and unconditional energy stability of the scheme and further simulate various numerical examples in 2D and 3D to demonstrate the stability and accuracy numerically. In addition to the above-mentioned benefits, another important advantage of the new decoupling method is that it is universally applicable. For instance, first, it can be combined with other linear methods (such as the linear stabilization, IEQ, etc.) to form various types of fully-decoupled energy-stable schemes. Second, it can also be applied to any nonlinear coupling type model to form an effective complete decoupling scheme, as long as the coupling terms follow the "zero-energy-contribution" property.
The rest of the paper is organized as follows. In Section 2, we establish the volume-conserved flow-coupled multi-phase-field EBE model. In Section 3, we develop a numerical scheme to solve the model, propose a detailed process to realize fully-decoupled implementation, and prove the solvability and unconditional energy stability. In Section 4, we perform numerical simulations, including the accuracy/stability tests and other examples in 2D and 3D to show the accuracy and efficiency of the developed scheme. Section 5 provides some concluding remarks.
Now, based on the classical approximate volume conservation multi-phase-field EBE model developed in [1,2,41,42], a new flow-coupled multi-phase-field EBE model with accurate volume conservation is established.
We introduce N scalar phase-field variables φ i (x) = tanh (
, where m is the dimension, d i (x) is the signed distance between a point x and the membrane surface Γ i , positive inside and negative outside, and ϵ ≪ 1 is a transition parameter, which is used to characterize the width of the diffusive interface or transition layer.
Therefore, the phase-field elastic bending energy for the ith vesicle is formulated as follows (cf. [1,2]):
where 2 is the Ginzburg-Landau double-well potential and f
. While using the elastic bending energy to describe vesicles, a fixed surface area for each vesicle is often required. The surface area function A(φ i ) is defined as
2)
The surface area can be given as 3 2 √ 2 A(φ i ) that converges to the surface area quadratically as ϵ → 0. To conserve the surface area of the ith vesicle, a penalized potential is introduced into the total free energy, cf. [1,2], that reads as
where M is a positive penalty parameter, and β i denote the constant related to the surface area, respectively. Throughout in this paper, we set
which is the initial surface area for each type of vesicle.
Regarding the interaction of different type of vesicles i and j with i < j, the adhesion energy between different vesicles is formulated as
where C i j > 0 is the magnitude of adhesion potential.
After coupling with the fluid momentum, the total free energy reads as
where
x is the kinetic energy, u is the fluid velocity, the energy density functional W (φ 1 , . . . , φ N ) (including the elastic bending potential, surface area potential, and adhesion potential) reads as:
and λ is the relative ratio of the kinetic energy to the elastic energy. By adopting the gradient flow approach in L 2 -space, i.e., the Allen-Cahn relaxation dynamics, and assuming the generalized Fick's law, i.e., the mass flux is proportional to the gradient of the chemical potential, we obtain the governing dynamical system that reads as
)
where
and γ i > 0 is the relaxation time scale parameter, the term (u • ∇)φ i is the advection for ith type vesicle, µ i ∇φ i is the induced stress term.
We consider one of the following boundary conditions:
or all variables are periodic, (2.10) where n is the unit outward normal on the boundary ∂Ω . The initial conditions read as
Remark 2.1. Note that the global volume of the vesicle is defined as (see [1,2,43,44]):
(2.12)
Hence, by computing the L 2 -inner product of (2.6) with 1 and using the divergence-free condition (2.9) for the velocity field, we derive d dt
which means the model (2.6) retains the exact volume. It can be seen that the nonlocal term -1
∫ Ω µ i d x added in (2.6) can accurately conserve the total volume of φ. This idea was originally ingeniously proposed in [10] to develop the conservative Allen-Cahn equation.
Remark 2.2. For the sake of completeness, here we also give the classical multi-phase-field EBE model (no flow-coupled case) developed in [1,2]. Similar to the surface area potential, a penalization potential is added to the total free energy to enforce the volume constraint approximately as
where M v ≫ 1 is the penalty parameter, α i = V (φ 0 i ) is the initial volume for the ith type vesicle. Note that when M v is very large, the volume remains approximately the same, but a small time step is required due to the increased stiffness.
The system (2.6)-(2.9) admits the law of energy dissipation, which can be obtained by the following process. By multiplying the inner product of (2.6) with λµ i in L 2 and taking the summation for i = 1, . . . , N , we get
where we use
By multiplying the inner product of (2.7) with -λφ it in L 2 and taking summation for i = 1, . . . , N , we get
Taking the inner product of (2.8) with u in L 2 , and using integration by parts and (2.9), we obtain
(2.17)
Combining the above three equations (2.15)-(2.17) and using the divergence-free condition for u, we obtain the energy dissipation law as
where the two negative terms on the right end specify the energy diffusion rate.
Remark 2.3. We note that when deriving (2.18), the nonlinear integrals related to advection and stress are all canceled out. More precisely, the following two identities hold
where the second one is due to the divergence-free condition (2.9), the boundary condition for u specified in (2.10), and the integration by parts. The two identities mean that these nonlinear terms do not contribute to the total free energy or energy diffusivity, that is, they satisfy the "zero-energy-contribution" property. We will take advantage of this feature to develop the full decoupling type scheme in the next section.
We now construct a fully-decoupled scheme for solving the highly nonlinear, flow-coupled, volume-conserved multi-phase-field EBE model (2.6)-(2.9). Meanwhile, considering the efficiency and accuracy of the algorithm in practice, we also expect the scheme to satisfy linearity, second-order time accuracy, and unconditional energy stability. The detailed process is given as follows.
First, we introduce a nonlocal variable Q(t) and an ODE system related to it that reads as:
u| ∂Ω = 0, or all variables are periodic.
(3.1)
By utilizing the "zero-energy-contribution" property satisfied by the advection and stress terms, we find that the ODE (3.1) is equivalent to
Second, we extract two linear terms from the free energy density W (φ 1 , . . . , φ N ) and define
where 0 < e 1 < 1 2 , and e 2 > 0. Thus the total free energy (2.5) is rewritten as
x can be shown to be bounded from below (see Remark 3.1), thus we define a nonlocal variable U (t) as the square root of it, that reads as
where B is a positive constant such that the radicand is positive. This is the so-called SAV method [11,22,[45][46][47][48] which is an efficient method to linearize the nonlinear energy potential. Using the variable U , we rewrite (2.7) as the following equivalent form:
where the function H i is defined as:
(3.6) Remark 3.1. We roughly outline the reasons why ∫ Ω W (φ 1 , . . . , φ N )d x is bounded from below. We expand the bending potential as
Therefore, we can see that there are several negative terms in W (φ 1 , . . . , φ N ) which are needed to be bounded from below, including the adhesion potential, the quadratic term -
2 |φ| 2 d x, as well as the term
Moreover, for the adhesion potential, some fourth-order polynomial terms will be generated from it after applying the Cauchy-Schwarz inequality (i.e., -C i j (φ 2 i -1)(φ 2 j -1) ≥ - 2 to bound the quadratic terme 2 2 |φ| 2 and those negative fourth-order polynomial terms from the adhesion potential. Moreover, the combination of two terms
x from below after using the Hölder's inequality and the Sobolev inequality ∥φ∥ L 4 ≤ C Ω ∥φ∥ H 2 where C Ω is a constant (see [22]).
Third, by combining (3.1) and (3.5), we rewrite the PDE system (2.6)-(2.9) to:
)
)
with the boundary conditions as
and initial conditions as
(3.15)
Remark 3.2. The new system (3.8)-(3.13) is obtained by combining (3.1), (3.5) and (2.6)-(2.9), where (2.7) is replaced by (3.5). Another modification is that we multiply the coupled nonlinear terms satisfying the "zero-energycontribution" property with the nonlocal variable Q in (3.8) and (3.11). This modification is also reasonable and will not change the equivalence of the system since we (3.13) still implies Q = 1. Thus, the new PDE system (3.8)-(3.13) using the variables (u, p, µ i , φ i , U, Q) is equivalent to the original PDE system (2.6)-(2.9) using the variables (u, p, µ i , φ i ).
The new system (3.8)-(3.13) also holds the energy dissipation law that can be obtained through a similar process obtaining (2.18). Since the discrete energy stability proof process follows the same principle, we show the following detailed process to make it clear.
We multiply the L 2 inner product of (3.8) with λµ i and take the summation for i = 1, . . . , N to derive
By multiplying the inner product of (3.9) with -λφ it in L 2 and taking the summation for i = 1, . . . , N , we get
(3.17)
Multiplying (3.10) with 2λϵU , we get
Taking the inner product of (3.11) with u in L 2 , and using integration by parts and (3.12), we obtain
Multiplying (3.13) with Q, we get
Combining (3.16)- (3.20) and noting that all two terms marked with the same Roman numerals are canceled, we derive the energy law as follows:
where ). In other words, when we try to develop a decoupled and energy stable discrete scheme, we can use different ways to discretize the advection and stress terms to generate a fully-decoupled scheme.
We are now ready to develop a second-order semi-discrete scheme to solve the system (3.8)-(3.13). Given (u,
In the above scheme,
, S 2 , S 3 are positive stabilization parameters, and the boundary conditions are either periodic or the physical boundary conditions as
We explain some details of the scheme (3.23)-(3.29) in the following remarks.
Remark 3.4. The scheme is linear, and it uses implicit and explicit combination method to deal with all nonlinear terms. For the hydrodynamical equations, we use the second-order pressure correction scheme (3.26)-(3.28)-(3.29) (see the overview of projection methods in [49]), of which ũn+1 is the intermediate velocity following the Dirichlet boundary conditions (or periodic) and the final velocity field u n+1 follows the divergence-free condition. To obtain the pressure, we just apply the divergence operator to (3.28) and then obtain the following Poisson equation for p n+1 , i.e.,
with the periodic boundary condition or ∂ n p n+1 | ∂Ω = 0. Once p n+1 is computed from (3.31), we update u n+1 by using (3.28), i.e.,
Remark 3.5. The initialization of the second-order scheme requires all values at t = t 1 , which can be obtained by constructing the first-order scheme based on the backward Euler method. In the above second-order scheme (3.23)-(3.29), as long as we set a = 2, b = 2, c = 0, ψ * = ψ 0 for any variable ψ, the first-order scheme can be easily obtained. Moreover, by using mathematical induction, it is easy to conclude that the following volume conservation property holds:
Remark 3.6. Three additional second-order linear terms are the commonly linear stabilizers (associated with S 1 , S 2 , and S 3 ) used in the SAV type schemes when adopting larger time steps or the system has very high stiffness issue due to the model parameters. Although the SAV method is formally unconditionally energy stable, but exceedingly small time steps are needed to achieve reasonable accuracy. To fix such an inherent deficiency, the SAV approaches are combined with the stabilization technique to construct stabilized-SAV (S-SAV) methods (see also a review for SAV method in [50]). The reasons to use three different stabilizers for the particular EBE model are as follows. Note that the coefficients H i , i = 1, . . . , N contain almost all explicitly processed terms, including the fourth-order terms ∆ 2 φ i . As we all know, the explicit processing method used for higher-order linear terms is unstable, so we have to restore the higher-order terms by using the second-order stabilizer with comparable magnitude, which is the reason for adding the S 3 term. The use of S 1 and S 2 is for similar reasons, and these two terms are used to balance the explicit processed terms f 2 (φ i ) and ∆ f (φ i ) contained in H i , respectively. In Section 4, we provide numerical evidence that these stabilizers are important for maintaining accuracy while enhancing the stability of the numerical schemes with larger time steps, see Fig. 4.2(b).
Now, we discuss the specific approach of implementing the scheme (3.23)-(3.29). Since (3.28)-(3.29) is the standard step of the projection method, we only need to consider the implementation of (3.23)- (3.27).
Note that the scheme (3.23)-(3.27) is not in the decoupled form, on the contrary, all unknowns are still coupled together. This means that a direct method to solve the scheme will bring high computational costs to actual calculations. Therefore, in practice, we implement the scheme through the following steps, which make full use of the nonlocal properties of the auxiliary variables U and Q. The following implementation method can decouple all variables and delete all nonlocal calculations, as shown below.
First, we use the nonlocal variable Q n+1 to split (φ i , µ i , U ) n+1 into a linear combination form that reads as
(3.34)
Then the scheme (3.23)-(3.24) can be rewritten as
According to Q n+1 , the system (3.35) can be split into two subsystems as
and
(3.37)
By taking the L 2 inner product of the first equation in (3.36) and (3.37) with 1, using (3.33), and noting ∇ • u * = 0 and (3.30), we immediately get
The boundary conditions of the φ n+1 i1 and φ n+1 i2 are either periodic or
Second, using the nonlocal variables U n+1 1 and U n+1
, we split the variables (φ i1 , φ i2 , µ i1 , µ i2 ) n+1 into the following form
(3.40)
Replacing (φ i1 , φ i2 , µ i1 , µ i2 ) n+1 in (3.36) and (3.37) using (3.40), and splitting the obtained equations according to U n+1 1 and U n+1
, we get the following four subsystems,
and
(3.44)
The boundary conditions of the above four system are either periodic or
By taking the L 2 inner product of the first equation in the above four subsystems with 1, using (3.33), and noting ∇ • u * = 0 and (3.30), we get
To solve the above four subsystems (3.41)-(3.44), we simply combine the two equations in each subsystem and use (3.46) to get the following four equations
where
i22 . Hence we can easily solve three independent biharmonic equations in (3.47) to get (φ i11 , φ i12 , φ i21 , φ i22 ) n+1 . Third, we rewrite (3.25) into the following form
where
We replace U n+1 and φ n+1 i using the split form given in (3.34) to get
We split the above result according to Q n+1 to derive
(3.51)
We replace (φ i1 , φ i2 ) n+1 using the split form given in (3.40) to get
(3.52)
By applying the simple factorization for each equality in (3.52), we derive
We need verify (3.53) is solvable by showing the denominator 1 -
thus the two denominators in (3.53) are the same). By taking the L 2 inner product of the fourth equation in (3.47) with φ n+1 i22 , and using Fourth, we use the nonlocal variable Q n+1 to split the velocity field ũn+1 as the following form:
By replacing the variables ũn+1 in (3.26), and then splitting the obtained equation according to Q n+1 , we arrive at a system that includes two linear elliptic sub-equations with constant coefficients as follows:
The two split variables ũn+1 1 , ũn+1 2 follow the boundary conditions described in (3.30), i.e., they are either periodic or satisfy:
Fifth, we solve the auxiliary variable Q n+1 . Using the split form for the variables ũn+1 in (3.55) and µ n+1 i in (3.34), one can rewrite the scheme (3.27) as the following form:
where
(3.59)
We need to verify the solvability of (3.58) by showing a 2δt -ϑ 2 ̸ = 0 as follows. By multiplying the L 2 inner product of the second equation in (3.56) with ũn+1 2 , we get
By taking the L 2 inner product of the first equation in (3.37) with λµ n+1 i2 , of the second equation with -λ a 2δt φ n+1 i2 , and combining the two obtained equations, we derive
(3.61)
From the second equation in (3.51), we derive
The following theorem ensures that the developed scheme (3.23)-(3.29) satisfies the energy stability unconditionally.
Theorem 3.1. The following discrete energy dissipation law holds for the scheme (3.23)-(3.29),
where
) .
(3.65)
Proof. We multiply the inner product of (3.26) with 2δt ũn+1 in the L 2 space, we obtain
(3.66)
From (3.28), by taking the L 2 inner product of any variable v with ∇ • v = 0 and v • n| ∂ Ω = 0 (or all variables are periodic), we have
(3.67) Using (3.67), we derive following equality
where we use the following identity
We reformulate the projection step (3.28) as
By taking the square of both sides of the above equation, we get
Hence, by multiplying 2δt 2 /3 of the above equation, we derive
By taking the inner product of (3.28) with 2δtu n+1 in the L 2 space, we have
We combine (3.66), (3.68), (3.72), and (3.73) to obtain
(3.74)
Computing the inner product of (3.23) with 2λδtµ n+1 i in the L 2 space, we have
(3.75)
Computing the L 2 inner product of (3.24) with -λ(3φ
))
(3.76)
By multiplying (3.25) with 2λϵU n+1 and using (3.69), we obtain
(3.77)
By multiplying (3.27) with 2δt Q n+1 and using (3.69), we obtain 1 2
Hence, by combining (3.74)-(3.78) and taking the summation for (3.75)-(3.76) with i = 1, . . . , N , we arrive at
where we use the following identity:
Finally, we obtain (3.64) from (3.79) after dropping the positive terms in { } and dividing both sides by 2. □
In this section, we implement numerical simulations to verify the accuracy and energy stability of the proposed fully-decoupled scheme (3.23)-(3.29) (denoted by DSAV for short). Numerical simulations include accuracy/stability tests, vesicle deformation in 2D and 3D, etc. In all numerical examples, we consider a rectangular computational domain. For directions with periodic boundary conditions, the Fourier-spectral method is used for discretization. For directions with boundary conditions specified in (3.30), the Legendre-Galerkin method is adopted for discretization, where the inf-sup stable pair (P N , P N -2 ) is used for the velocity ( ũ and u) and pressure p, respectively, and P N is used for variables φ i , µ i . (a) Time evolution curves of the ratio of surface area change that is computed using δt = 0.01/2 2 , and (b) the numerical errors in L 2 for all variables computed by using the scheme DSAV with different time steps. (For simplicity, we plot the average of the errors of the two components of the phase-field variables (φ 1 , φ 2 ) and the velocity field u = (u, v). Moreover, for the nonlocal variable Q, we directly compute its error with the exact solution Q(t) = 1.).
We now test the convergence rates of DSAV by refining the time step. We set the 2D domain as (x, y) ∈ Ω = [0, 2π ] 2 and assume the periodic boundary conditions for each direction which is discretized by sufficient fine meshes so the error in spatial directions can be ignored compared with the time discretization errors. We test the two-phase-field model and set the initial condition to be two adjacent circular vesicles with
and set the model parameters to be
Using 129 Fourier modes for each direction and the time step size δt = 0.01/2 2 , we implement the designed scheme DSAV until the steady-state solution is obtained. In Fig. 4.1(a), we plot the evolution curve of the surface area difference ratio (defined as
over time, where the initial and final steady-state profiles of φ 1 -φ 2 are also attached therein. We can see that the surface area difference ratio has been maintained around 1e-3. The two vesicles changed from a state of slight contact at the initial moment to a state of tightly adhering together, where the contact part is a straight line due to the adhesion potential. In Fig. 4.1(b), we plot the L 2 errors for all variables when t = 0.2 computed by using various time steps δt. Since the exact solution is not known (except the auxiliary variable Q and Q(t) = 1 is the exact solution), we treat the numerical solution calculated by DSAV with a small time step of δt = 1e-8 as the approximate exact solution (the error of Q n+1 is computed by using its exact solution of Q(t) = 1). It can be seen that the scheme DSAV presents almost perfect second-order convergence rates for all variables including the auxiliary variable Q.
Furthermore, we use different time steps to test the energy stability of the scheme DSAV. In Fig. 4.2(a), we plot the evolution curve of the total free energy in the original form (2.5) calculated using different time steps. It can be seen that the obtained energy evolution curves always show monotonic decays, which means that DSAV is unconditionally energy stable. Moreover, as the time step gets smaller and smaller, the obtained energy curves overlap with each other, indicating that the obtained solution is getting closer and closer to the accurate solution. We also append a small inset figure to show the comparisons between the original free energy (2.5) and the discrete energy (3.65) computed using δt = 0.01/2 3 and there is almost invisible difference between them. In Fig. 4.2(b), we use different stabilization parameters to verify their effects on improving energy stability by plotting the evolution of the original free energy (2.5). We find that when S 1 = S 2 = 0, the energy decreases only when the time step is very small (δt ≤ 0.01/2 5 ), and when the time step is large, the energy even increases reflecting the instability. We also compare the two energy curves obtained with (S 1 = S 2 = 0, δt = 0.01/2 5 ) and (S 1 = S 2 = 4, δt = 0.01/2 2 ). The coincidence shows that after using the stabilizers, the time step can be increased around 8 times. In Fig. 4.3, we plot the velocity field and the pressure p at various times. In this example, we study the dynamic deformation of vesicles using the single-phase-field model (i.e., N = 1). For this model, the adhesion energy vanishes. We set the computational domain to [0, L] 3 where L = 1.5π , and the initial conditions are given as follows:
where k is the number of vesicles. For k = 1, we set the initial shape of vesicle to be a wide or narrow ellipsoid.
For the former, we set r 11 = 0.75, We use 129 Fourier modes to discretize each of the x and y directions which are assumed to be periodic and use the Legendre polynomials up to the degree of 128 to discretize the z direction which is assumed to satisfy the boundary conditions specified in (3.30).
For k = 1 with the initial shape as a wide ellipsoid, in Fig. 4.4(a), we plot the time evolution of surface area change rate under two surface area parameters M, and attach the snapshots of an isosurface of {φ = 0} at different times. When M = 0, although the volume remains unchanged, the ellipsoid gradually shrinks, and the ratio of surface area change eventually reaches about 16%. When M = 1e5, the final shape of the ellipsoid vesicle becomes the pancake shape, where the middle region is slightly thinner, and the ratio of surface area change is always around 1e-4. In Fig. 4.4(b), the surface area change and profiles of the narrow ellipsoid for M = 0 and M = 1e5 are plotted. When M = 0, the narrow vesicle shrinks to the spherical shape, and when M = 1e5, it becomes a capsule shape. The surface area changes are approximately 25% and 5e-4, respectively.
For k = 2 (two stacked vesicles), as shown in Fig. 4.5(a), we can see that when M = 0, the vesicles fuse and shrink, and the surface area ratio changes up to around 14%; when M = 1e5, the steady-state shape is displayed as a capsule, and the surface area ratio is approximately 1e-4. For k = 4 (four stacked vesicles), as shown in Fig. 4.5(b), we can see that when M = 0, the vesicles fuse and shrink to form a structure with multiple pores, and the surface area ratio changes up to 14%; when M = 1e5, the steady-state shape forms a ring and the surface area ratio is approximately 2e-4.
In this example, we simulate the multi-phase-field model to study the dynamic deformation of vesicles. In each of the following figures, we plot the time evolutions of the ratio of surface area change and attach the topological profiles of φ i at various times using different colors. For the sake of simplicity, we ignore the detailed initial condition formulas of φ i , because they are just set to some simple shapes like circles, spheres, ellipses, or ellipsoids. The computational domain is set to [0, L] m , m = 2, 3 with L = 2π with the periodic boundary conditions, and each direction is discretized using 257 Fourier modes.
The model parameters are set as
In Fig. 4.6(a)-(f) with N = 2, we set the initial conditions of φ i to be multiple slightly touching 2D vesicles stacked in different positions. We can see that, due to the adhesion energy, the vesicles adhere to each other tightly and they finally evolve into different patterns. The ratio of the surface area changes are all kept to be small due to the surface area constraint potential.
In
In this example, we simulate the sedimentation of multiple vesicles driven by the gravity force. The sedimentation of a single vesicle had been studied experimentally in [51], and simulated using the different model in [52]. For simplicity, we approximated the gravity force using the Boussinesq approximation, namely, the momentum equation is replaced by We perform 3D simulations for the multi-phase-field model by setting the computational domain is set to be (x, y, z) ∈ Ω = [0, 2π ] × [0, 2π ] × [0, 4π]. We set periodic boundary conditions along the x and y directions and discretize each of them using 128 Fourier modes. For the z-direction, we use the boundary conditions given in (3.30) and use the Legendre-Galerkin method of Legendre polynomials to the degree of 128 for discretization. In Figs. 4.9-4.11, we simulate the sediment process of two symmetrically placed vesicles (N = 2), two asymmetrically placed vesicles (N = 2), and three asymmetrically placed vesicles (N = 3). We use different colors to plot the interface contour {φ i = 0} at different times of the dynamical process. We find that these vesicles exhibit different dynamic changes of merging, and eventually collide and fuse together due to the adhesion potential, where the vesicles that are asymmetrically descended show a state of tilted angle due to the interactive adhesion.
In this example, we simulate the dynamic motion of multiple vesicles driven by the shear flow imposed on the boundary. The shear flow dynamics for a single vesicle had been studied in [53,54]. We set the computational domain to be (x, y, z) ∈ Ω = [0, 4π] × [0, 2π ] × [0, 2π ]. We set periodic boundary conditions along the x and y directions and discretize each of them using 128 Fourier modes. For the z-direction, we use the boundary conditions given in (3.30) and use the Legendre-Galerkin method of Legendre polynomials to the degree of 128 for discretization. We set the boundary condition of the velocity field u = (u, v, w) as u| z=0,2π = ±u 0 , v| z=0,2π = w| z=0,2π = 0. In Fig. 4.12, we simulate the dynamical process of three vesicles driven by the shear flow on the boundary. We use different colors to plot the interface contour {φ i = 0} at different times of the dynamical process. We find that the In this example, we simulate the dynamic motion of multiple vesicles driven by the Poiseuille flow. We perform 3D simulations for the multi-phase-field model with N = 3, where the computational domain is set to be (x, y, z) ∈ Ω = [0, 4π ] × [0, 2π ] × [0, 2π ]. We set periodic boundary conditions along the x and y directions In Fig. 4.13, we use different colors to plot the interface contour {φ i = 0} at different times of the dynamical process, where we can see that the vesicles form wrinkles due the flow field which had been also observed by experiments in [55,56]. Due to the initial shape of the vesicles and the different angles facing the flow field, different dynamic changes are observed. For example, the surface of the vesicles that were initially ellipsoid show more wrinkles, while the vesicle that is initially spherical were compressed into a flat oblate shape over the time.
In this article, we first modify the classic multi-phase-field EBE model to establish a new volume-conserved model. Then, for the flow-coupled model, we develop a new decoupling method, and combine it with other proven effective numerical methods, to form an efficient numerical scheme. At each time step, one only needs to solve a few linear fully-decoupled equations with constant coefficients. The specialty of the new decoupling method is that the nonlinear coupling terms have been explicitly processed while still ensuring the unconditional energy stability. This method can also be extended to other flow coupling models as long as the coupled nonlinear terms satisfy the "zero-energy-contribution" property. We also conduct numerous numerical tests in 2D and 3D to demonstrate the accuracy and stability of the scheme.