Efficient numerical solution of discrete multi-component Cahn–Hilliard systems

Computers and Mathematics with Applications 67 (2014) 106–121

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Efficient numerical solution of discrete multi-component Cahn–Hilliard systems P. Boyanova a,∗ , M. Neytcheva b a

Institute for Information and Communication Technology - BAS, Acad. G. Bonchev Str., bl. 25A, 1113 Sofia, Bulgaria

b

Uppsala University, Department of Information Technology, Box 337, 751 05 Uppsala, Sweden

article

info

Article history: Received 30 April 2013 Received in revised form 20 September 2013 Accepted 26 October 2013

abstract In this work we develop preconditioners for the iterative solution of the large scale algebraic systems, arising in finite element discretizations of microstructures with an arbitrary number of components, described by the diffusive interface model. The suggested numerical techniques are applied to the study of ternary fluid flow processes. © 2013 Elsevier Ltd. All rights reserved.

Keywords: Cahn–Hilliard equation Multi-component systems Phase-field model Preconditioning

1. Introduction Microstructure evolution in complex multiphase–multicomponent systems has been and is of significant importance for numerous industrial processes as well as for academic research, technology improvements and development of new materials and devices. Phase transitions and morphological changes in microstructures take place in a variety of processes, such as spinodal decomposition, grain and crystal growth, solute drag, directional solidification, diffusion-controlled processes, reaction pathways controlling the structural evolution of complex material mixtures and many more (see, for instance, [1,2]). Multiphase, and especially, three-phase flows are present in many important application fields. For example, the behaviour of the gas–oil–water system in porous media is the major problem dealt with in petroleum engineering and, in particular, in developing efficient recovery technologies. The challenging task is to accurately solve the three-phase flow problem on micro (pore) scale in order to allow for better macro-scale models with less heuristics as in the present ones, based on Darcy laws. Today, a deeper understanding of the dynamics of those complex and coupled phenomena is to a large extent gained by computer simulations of the underlying processes, based on adequate mathematical models. In order to treat multicomponent problems, various techniques have been developed, falling into two major classes—interface tracking (see, for instance, [3,4]) and interface capturing (see, e.g., [5–7]). Among the most utilized tools for studying microstructure evolution is the phase-field (PF) model, which is used also in this work. This paper is organized as follows. In Section 2 we describe briefly the basic principles of the phase-field model and some of its forms, which are most used in practice. Section 3 is devoted to phase-field models, based on the Cahn–Hilliard equation, together with some constitutive relations to tune the equations to particular types of physical processes. In Section 4, we

∗

Corresponding author. E-mail address: [email protected] (P. Boyanova).

0898-1221/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.camwa.2013.10.013

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

107

describe the fully discretized Cahn–Hilliard equation with a constitutive relation of polynomial type, the arising nonlinear and linearized algebraic systems to be solved and their properties. We also describe and analyse the numerical solution method and the preconditioning techniques. Section 5 contains numerical illustrations for ternary fluid flow problems. Conclusions are found in Section 6. 2. The phase-field model—general concepts The development of the PF model has started with the works by Cahn and Hilliard [8,9] and since then the model has been developed further to be able to properly treat a broad variety of multiphase–multicomponent systems. We describe first the general idea of the PF model, following [10–12]. The phase-field model is based on thermodynamic principles and on the assumption that interfaces between phases or components in a microstructure are diffused, and each phase component can be represented by some smoothly varying function, referred to as the order parameter, mass fraction or the concentration. The basic idea is to let the order parameter obtain a distinct value in each bulk phase/component of the mixture and then, sharply but smoothly change from one to another within the interfacial regions. In this way the need to explicitly track the interfaces is avoided and, thus, also the difficult task to impose and handle proper boundary conditions on those interfaces. Creation, merging and destruction of interfaces, which are major obstacles in the sharp interface models, are implicitly built in the PF model. Further, PF incorporates in a natural way curvature-driven physics, allows for an easy inclusion of other phenomena in the model, such as surface tension and dynamic wetting, in this way combining the large scale simulation of the hydrodynamic interaction with the constrain of micro scale properties (cf. [13]). It also captures the behaviour of the components that occur far away from the interfaces. As already mentioned, the important variable to determine in PF is the order parameter, denoted here by c = c (x, t ), which is often interpreted as the ‘concentration’ of a certain component at certain point in space and time. The model is based on minimization of the so-called free energy functional E, which involves c. The classical form of E (cf., e.g., [9]) is E (c , ∇ c ) =



f0 (c ) + υ|∇ c |2 dΩ ,

 Ω



where f0 (c ) is the homogeneous free energy density, υ is referred to as the gradient energy coefficient and Ω is the volume of the material system under consideration. The term υ|∇ c |2 describes the gradient energy, which acts as a penalty for sharply varying concentration. Once defining E, a variational method is used to determine evolutionary equations that steer the system towards the minimum of the energy functional. One particular such evolution pathway is the Cahn–Hilliard (C–H) equation,

  δE ∂c = ∇ · L(c )∇ , ∂t δc

(1)

applied to each component of the multiphase system. Here L(c ) is the so-called mobility, which may or may not depend on the concentration. Note that the relation (1) is a (componentwise) conservation law with respect to c. Another also often used evolutionary relation is the so-called Allen–Cahn equation

∂c δE = −L(c ) , ∂t δc which, however, does not conserve c. In this work we deal only with the C–H equation. Remark 2.1. In (1), δ E /δ c is the functional derivative of the functional E. For completeness, we recall that the functional derivative of a functional F of the form F (y) =

 Ω

G(x, y, ∇ y)dΩ =

 Ω

G(x1 , . . . , xn , y, yx1 , . . . , yxn )dΩ ,

where yxi = ∂ y/∂ xi , using the Euler–Lagrange equation, is found to be of the form n δF ∂G  ∂ ∂G ∂G ∂G = − = −∇ · . δy ∂y ∂ xi ∂ yxi ∂y ∂∇ y i=1

In its most general form, the PF model may describe the evolution of the microstructures of mixtures, containing multiple components, for instance, water, air, metals, and at the same time, each component may co-exist in several phases, say, liquid and solid. The model allows also for phase changes during the time evolution (see [10] for more detailed explanations). Here we consider the case, where we have multiple immiscible components, however no phase changes are included in the model. Thus, in the sequel, the words ‘component’ and ‘phase’ are interchangeable.

108

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

Computational burden The basic theoretical assumption, made in PF is that the interfaces are diffusive at nanoscale. Experimental systems can handle at most 10 nanometres (10−8 m) for the interface. The sizes of the bulk phases are usually of the order of micrometres or larger. The dynamics of the physical processes to be captured, occurs at the length scale of the interfaces, which, although being smooth, are rather sharp. As a rule of thumb, to resolve the interfaces numerically, 5–10 discretization points have to be placed within interface regions. This indicates that numerical simulations of interface dynamics pose very high demands on computer resources, both in terms of memory and computing power. As relatively less is happening in the bulk regions, at a first glance, locally refined meshes seem to be the right strategy towards reducing those demands. However, the interfaces do evolve in time and change their position in space rather dynamically, some of them might disappear or new may appear. This means that locally refined meshes, while perhaps reducing the number of degrees of freedom in the bulk regions per time step, may entail quite high overhead costs to handle the meshes throughout the time evolution, in particular, in a parallel environment. Remark 2.2. For problems where the solution undergoes rapid changes in some part of the domain, such as in layer and interfacial regions, the discretization must be done properly so that a good resolution is obtained. The latter implies that a certain accuracy of the discrete solution has to be guaranteed by choosing a sufficient number of discretization points to be placed in the layer regions. Such regions have to be accounted for in singularly perturbed convection–diffusion and parabolic problems, wave propagation and also in diffuse interface models. How to estimate the number of grid points that are necessary to achieve a certain required accuracy depends on the type of problem and on the discretization method. More details on this topic can be found, for instance, in [14] and the references therein, and in [15], where wave propagation problems and construction of layer-adapted (Bakhvalov–Shishkin) meshes for singularly perturbed differential equations are discussed, correspondingly. Another way to lower the computational demands in numerical simulations is to take the width of the interfaces to be bigger than that assumed in the model. However, thick interfaces may introduce unphysical effects, for instance, jumps in the chemical potential (defined in Section 3) at the interfaces [10]. Therefore, to still allow for thicker interfaces, certain adjustments in the model have been made, such as separating the contribution of gradient energy in the bulk from that at the interface, introducing the so-called ‘anti-trapping current’ to cancel the jump in the chemical potential, and [16,17] include the so-called ‘composition gradient energy’ instead of |∇ c |2 . 3. Cahn–Hilliard models for multi-component systems In this work we consider multiphase systems with n components, where n ≥ 3. For two-phase systems (n = 2) much research and algorithm development is already available (see, for instance, [18–27]). Various generalizations of the PF model to any number of components have also been done, see, for instance, [28,29,12,11,27,30–32]. We point out that the transition from binary to systems with three and more components is not straightforward and deserves special attention. Assume that we deal with a mixture of n immiscible components (phases), which occupy an isolated region Ω ∈ Rd , d = 1, 2, 3. In the sequel, n denotes the outward unit normal vector to Ω and e = [1, 1, 1 · · · , 1]T . Denote also by ci = ci (x, t ), x ∈ Ω , i = 1, 2, . . . , n the concentrations (the order parameters) of each of the phases. We denote c = [c1 , c2 , . . . , cn ]T and ∇ c = [∇ c1 , ∇ c2 , . . . , ∇ cn ]T . The variables ci (x, t ) can be seen as mole fractions and, thus, these must obey the relations n 

ci (x, t ) = 1,

x ∈ Ω , t ≥ 0,

and

0 ≤ ci (x, t ) ≤ 1,

x ∈ Ω , t ≥ 0,

(2)

i=1

which means that the only physically relevant values, which the concentration is assumed to attain, lie on the so-called Gibbs simplex

 G≡

n

ci ∈ R :

n 

 ci = 1, 0 ≤ ci ≤ 1 .

(3)

i=1

From (2) it follows that only n−1 of the concentrations are independent. Thus, one concentration can always be computed n−1 as cn = 1 − i=1 ci . In addition to the concentration ci , the so-called chemical potential ηi (x, t ) is another basic physical quantity, defined for each component. The vector η = [η1 , η2 , . . . , ηn ]T of all chemical potentials is taken to be the functional derivative of E, i.e., η = δδEc . Using the chemical potential, the Cahn–Hilliard equation can be written as

∂c = ∇ · (L(c)∇η) , ∂t

x ∈ Ω , t > 0,

(4)

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

109

coupled with no mass flux boundary conditions (L∇η)i · n = 0, for all i, x ∈ ∂ Ω , t > 0. Note that (4) is satisfied componentwise. The mobility matrix L(c) is symmetric and positive semidefinite. In order for (2) to be satisfied and to ensure mass conservation, L should also satisfy the additional constraint L(c)e = 0, i.e., L should have a one-dimensional kernel. The positive semidefiniteness ensures that the total free energy functional is decreasing in time. More details are found in, e.g., [11]. 3.1. The free energy functional and choices of the free energy density The general form of the free energy functional E (c, ∇ c) is E (c, ∇ c) =

 Ω

f (c, ∇ c) dΩ ,

(5)

where f (c) is the so-called free energy density, assumed to be nonnegative and smooth enough. The function f (c) can be chosen differently, depending on the target problem. It can be decomposed as f (c) = f0 (c) + f1 (c, ∇ c) + f2 (c), where f0 (c) is referred to as the homogeneous free energy density, f1 (c, ∇ c)—as the gradient free energy density and f2 (c) accounts for other effects, included in the model. For example, when considering also mechanical effects, f2 (c) = f2 (c, E (u)), where E is the linearized stress tensor and u(c) are the displacements with respect to some reference point in space at a given time. The sum f0 + f1 could be seen as obtained from a truncated Taylor expansion of a function f (c, ∇ c) around a homogeneous point (c, 0) or generalizations of that. Some typical choices of f1 (c, ∇ c) are n  1

2 i =1 n  1 i,j=1

2

υ|∇ ci |2

(6)

υij ∇ ci · ∇ cj

(7)

n  1 i,j=1,i
υij (ci ∇ cj − cj ∇ ci )

(8)

 υij (x)(ci ∇ cj − cj ∇ ci ).

(9)

2

As pointed out in [10], the gradient energy includes terms from the Taylor expansion of f (c , ∇ c , ∇ 2 c , . . .) around the point f0 = f (c , 0, 0, . . .). The particular form of f1 is problem dependent and incorporates symmetries/asymmetries, isotropic/anisotropic properties of the material, etc., see also [33,30], and various extensions to handle for instance nonisothermal multicomponent solidification [34] and elastic behaviour [35]. In general, Υ ≡ [υij ]ni,j=1 is a constant symmetric fourth order tensor of gradient energy coefficients and υij are d × d matrices with positive entries. In the case (9),  υij are convex functions that are homogeneous of degree two and allow for modelling anisotropy surface energy.  of the   n The homogeneous part f0 is assumed to have exactly n local minima on the hypersurface S ≡ ci ∈ Rn : i=1 ci = 1 . Relation (7) illustrates one of the substantial differences between binary and multiphase systems, namely, the possible presence of terms with coupled composition gradients. Remark 3.1. The matrix Υ has to be symmetric to ensure that free energy is invariant to reflections. In [10] it is argued that Υ must be also positive definite in order to ensure well-posedness of the problem and to avoid the appearance of unphysical surface energy. The argument is that if Υ has a negative eigenvalue, there will be a coupling of gradients, for which an increasingly sharp interface would lower the free energy of the system, producing a physically impossible surface energy. The set of evolutionary equations would become ill-posed as solving a diffusion equation with a negative diffusivity is unconditionally unstable. It turns out that for n ≥ 3 it is highly relevant to choose f0 (c) in such a way that the model is algebraically and dynamically consistent with the corresponding diphasic model, obtained when one of the phases is not present. Indeed, in many multicomponent problems, the interfaces most often separate two phases only and triple (or more) junctions may occur accidentally and for a relatively short time. An illustrative example is that of a gas bubble, raising in a stratified two-fluid system, where the problem becomes triphasic only when the bubble crosses the liquid–liquid interface, but during the rest of the time it is described by the diphasic model. Up to the knowledge of the authors, the notion of algebraic consistency is first introduced in [12], in the context of ternary systems with no phase changes. The requirements for algebraic consistency are the following:

110

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

(i) when the ith component of the system with n components is not present, the total free energy E (n) becomes exactly equal to E (n−1) of a system, containing the other n − 1 components; (ii) when the ith component is not present at the initial time, it does not appear during the time evolution. The dynamic consistency, also discussed in [12], means that the reduced (n − 1)-phasic solutions should remain stable with respect to perturbations of initial data. The latter is of particular importance for numerical simulations and the presence of computation-related errors. In this work we do not include any additional effects in the model, i.e., f2 ≡ 0. We also do not consider gradient free energy with couplings between the concentration components, thus Υ is diagonal. Various choices of the free energy density are possible and correspond to different constitutive assumptions. Next, we describe two such choices that have been much analysed and used in the related literature. 3.1.1. Logarithmic free energy A particular choice of the free energy density, referred to as ‘regular’, is used in several works on multiphase C–H models. The description below is based on [11,36,37,29,19]. The regular free energy density is of the form f0 (c) = θ

n 

ci ln ci −

i=1

1 2

cT Ac.

(10)

Here, θ is a constant, corresponding to the absolute temperature and A is a constant symmetric n × n matrix, which is not negative definite, thus, it has at least one non-negative eigenvalue. The requirement for the non-negative definiteness of A follows from physical considerations, namely f0 (c) should have more that one local minimum. Depending on the particular multicomponent system, there exists a critical temperature θc , such that for θ > θc (θ < θc ) f0 (c) is convex (non-convex). This form of f0 is often used when simulating spinodal decomposition. Note, that the logarithmic free energy density for n ≥ 3 is consistent with that for n = 2. The behaviour of the C–H models with logarithmic homogeneous free energy density has been tested numerically in several works, [36,29,38,19] and others. The attention has been focused primarily on the quality of the solution and not on the numerical solution methods used. The arising nonlinear algebraic systems have been solved by a nonlinear Gauss–Seidel method [20], shown to be stable for small time steps (1t ≈ h2 ). Operator splitting techniques for the nonlinear term have been applied (cf. [19]), and found to be faster than the Gauss–Seidel method. In [38] the nonlinear term is treated explicitly, which is shown to require also very small time steps (1t < Ch4 ). Although the numerical solution methods we develop are general applicable in the case of logarithmic form of the homogeneous free energy density, due to the above mentioned drawbacks it is not used in this paper. 3.1.2. Polynomial free energy In other works, the homogeneous free energy density is chosen to be of polynomial type. The evolution of the system is  driven again by the minimization of E under the constraint of mass conservation for each component, Ω ci dΩ = const. For two-phase models f0 is a double-well potential, most often chosen as f0 (c ) = c 2 (1 − c )2

or f0 (c ) = (1 − c 2 )2

(11)

with minima at c = 0 and c = 1 in the case (11) (left) and at c = ±1 in the case (11) (right), see for instance [13,2,39,40] and the references therein. For example, in [13,40,39], the double well potential is of the form 14 (1 − c 2 )2 . Observe, that for the diphasic model it suffices to compute only one of the concentrations c and the other one is recovered as 1 − c. For systems with n ≥ 3 phases, various polynomials are considered as generalizations of (11). One example is f0 (c) =

n 

αij ci2 cj2 +

i ,j = 1

n 

mi (T )ci2 (3 − 2ci ),

i =1

where αij > 0 and mi (T ) are linear in the temperature T (cf. [33]). In [32,41] some rather specific choices for ternary systems have been considered, such as f0 |i =

1 4

ci2 (1 − ci )2 −

1 2

ci

3 

cj .

j =1

A detailed general explanation of the latter form can be found in [33,30]. In some of the works, a scaling, directly related to the width of the diffuse interface, is introduced, f =

1

ϵ

f0 (c) + ϵf1 (c, ∇ c).

In [12,42,43], the free energy functional is chosen to be a polynomial, scaled as above, depending explicitly on the interface width ϵ but also on the surface tension σ. For ternary systems, the corresponding form of the free energy functional

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

111

is as follows: E (c; ϵ, σ) =

 

12

Ω

ϵ

f0 (c) +

3 13

2 i=1 4

 ϵυi |∇ ci |

2

dΩ ,

(12)

where υ1 , υ2 , υ3 are some constant coefficients to be determined as functions of the given in advance surface tensions σ12 , σ23 , σ13 . The formulation of (12) is based on a model for binary systems with the free energy functional (2) Eσ,ϵ

 σ 2 3 2 2 12 c (1 − c ) + σϵ|∇ c | dΩ . ϵ 4

  = Ω

(13)

As pointed out in [44], this form of the free energy functional turns out to be very suitable for computer simulations, since it allows to choose ϵ larger than its theoretical value without modifying the capillary properties at the interfaces. In [12] it is shown that the algebraic consistency conditions (i) and (ii) hold for the PF model with E as in (12) if and only if the following conditions for the coefficients υi are satisfied:

υi = σij +σik −σjk ,

for all i, j, k = 1, 2, 3.

(14)

The coefficients −υi are known as the spreading coefficients of phase i at the interface between the phases j and k (see, for instance, [12,44] and the references therein). If υi < 0, the spreading is total and if υi > 0, the spreading is partial. There is evidence, indicating that the spreading coefficients are very important for three-face flows, as in porous media (cf. [45,46]). If total spreading conditions are acting along a certain interface, the interfacial energy of the three-phase fluid system can be decreased by having a film of oil between the gas phase and the water phase, and thus, oil spreads spontaneously between the gas and water, while in case of partial spreading the above effect is not observed. The triphasic model, analysed in [12,44], is based on the free energy form (12) and could be seen as simulating three times two-phase models, due to the enforced strict compatibility of the triphasic model with the diphasic one, as in (13). One choice of the homogeneous free energy density, which leads to well-posed triphasic models, consistent with the corresponding (1) (2) reduced diphasic models and also handling both partial and total spreading, is of the form f0 (c) = f0 (c) + f0 (c), where (1)

f0 (c) = σ12 c12 c22 +σ13 c12 c32 +σ23 c22 c32 + c1 c2 c3 (υ1 c1 + υ2 c2 + υ3 c3 ),

(15)

(2)

f0 (c) = 3λυ1 υ2 υ3 (1)

with some non-negative real parameter λ. Using the relation (14), the function f0

in (15) is transformed into

 1 υ1 c12 (c2 + c3 )2 + υ2 c22 (c1 + c3 )2 + υ2 c32 (c3 + c2 )2 . 2 The following conditions are shown to be necessary for having an algebraically and dynamically consistent well-posed problem: (1)

f0 (c) =

υ1 υ2 + υ2 υ3 + υ1 υ3 > 0, υi + υj > 0, i ̸= j.

(16)

Note that, the conditions (16) do not define υi in a unique way. A more detailed discussion can be found, for instance, in [12]. As commented in [12], λ = 0 suffices to model partial spreading while for problems with total spreading, λ must be strictly positive and large enough. Further, as in [30,12], we assume that Li are constant and consider the so-called isotropic form of the equations by determining υi in such a way, that the products Li υi are constant, i.e., there holds L = Li υi , i = 1, 2, 3. This corresponds to introducing scaled chemical potentials ηi /υi . Following [12], the ternary C–H system takes the form

 ∂c  i = ∇ · (Li ∇ηi )   ∂t     4υ  1 ∂ f0 ∂ f0 3ϵυi  ηi = T − − 1ci  ϵ j̸=i υj ∂ ci ∂ cj 4

(17)

for i = 1, 2, 3 where υT is the harmonic average, defined via the equality υ3 = υ1 + υ1 + υ1 . T 1 2 3 4. Solution methods and preconditioning techniques 4.1. Formulation of the three-phase model The framework for numerical simulations of multi-phase problems, presented in this work, is based on the structure and properties of the discrete C–H problem and is a generalization of that for diphasic problems, developed in [40,39]. We detail

112

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

it for the case of the algebraically and dynamically consistent triphasic model (17), (15), as also studied in [43]. In addition, in Section 4.4 we outline the extension to problems with more than three phases. (1) (2) Consider the problem (17), where the homogeneous free energy density f0 (c) = f0 (c) + f0 (c) has the form (15). Assuming that the initial data is admissible, that is c(x, 0) ∈ G, then c(x, t ) ∈ G and we can recover the concentration of the third phase as c3 = 1 − c1 − c2 at any time and space instance. The chemical potential η3 can be computed via the relation  3 i=1 ηi /υi = 0, see e.g. [12]. Thus, the main computational burden falls on solving a reduced problem for c1 , c2 , η1 and η2 . In order to take into account the incompressible hydrodynamics of the mixture, the C–H equations are coupled to the Navier–Stokes (N–S) equations, see e.g. [21,13,47,48]. Then, assuming a zero velocity jump between two phases, a convection term is added in the evolution equations for the concentrations and the PF problem takes the following form

 ∂c  i + (u · ∇)ci = L1ηi   ∂t     ∂ f0 3ϵ 4υ  1 ∂ f0  ηi = T − − 1ci ,  ϵυi j̸=i υj ∂ ci ∂ cj 4 for i = 1, 2. Here u is the velocity vector. To make Eqs. (18) dimensionless, we define u t uc t x u′ = , t′ = = , x′ = , lc uc tc lc

(18)

ηi′ =

ηi , ϵ

where lc is the characteristic length and uc is the characteristic velocity. Substituting these variables into Eqs. (18) and dropping the primes, we obtain

 1  ∂ ci  + (u · ∇)ci = 1ηi  ∂t Pe  3  ηi = Fi (c) − Cn2 1ci

(19)

4 where Pe = uc lc ϵ/L, Cn = ϵ/lc are the Peclet and the Cahn number, respectively. In (19), Fi (c) denotes the nonlinear term Fi (c) ≡

4υT 

υi



1



υj

j̸=i

∂ f0 ∂ f0 − ∂ ci ∂ cj



,

i = 1, 2.

(20)

For the particular choice of λ = 0 in (15), which is often used in numerical simulations, the term Fi (c) reduces to Fi (c) = 12(g (ci ) − ξ ci cj ck ) = 12(g (ci ) − ξ ci cj (1 − ci − cj )) = Fi (c1 , c2 ),

(21)

6υj υk

where g (r ) = r (1 − r )(1 − 2r ) and ξ = υ υ +υ υ +υ υ , i = 1, 2, j, k = 1, 2, 3, i ̸= j ̸= k. i j i k j k 4.2. Discrete formulation To discretize (19) in space, we use the finite element method (FEM). In this work we use standard conforming piecewise linear basis functions for the two unknown vector functions c and η. The time discretization can be done in various ways, for instance using the θ -method. For brevity we formulate the time discretization applying the backward Euler method. Consider the FEM approximations for the concentration and the chemical potential at a time step tk = t0 + k1t in the N N k k form ci (x, tk ) = m=1 ci,m χm (x), ηi (x, tk ) = m=1 ηi,m χm (x), where χm (x), m = 1, . . . , N are the FEM basis functions and

ck = [(ck1 )T , (ck2 )T ]T and ηk = [(ηk1 )T , (ηk2 )T ]T , cki = {cik,m }Nm=1 , ηki = {ηik,m }Nm=1 are the corresponding unknown FEM vectors. The fully discretized system (19) (after some rearrangements) takes the form

Mηk − F (ck ) − γ K ck = 0,

(22)

1t ωK ηk + Mck + 1t W ck − Mck−1 = 0,

, ω = Pe1 . The matrices    M K M= , K= ,

where γ =

3 Cn2 4



M

K

 W =



W W

consist of mass matrices M, stiffness matrices K and convection matrices W for each phase component. Note, that M , K , W admit such a block-diagonal form for any choice of gradient free energy function with no couplings between the components of the concentration. The nonlinear term F (ck ) = (F1T (ck ), F2T (ck ))T is a vector with elements

Fi,m (c ) = k



 Ω

Fi

N  l =1

c1k,l

χl (x),

N  l =1

 c2k,l

χl (x) χm (x)dx,

m = 1, . . . , N , i = 1, 2.

(23)

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

113

4.3. Nonlinear solution and preconditioning techniques Next, we use the experience from [40,39] and generalize the numerical solution techniques, which have been shown to be most numerically and computationally efficient. Namely, at each time step, we solve the system in (22) using a quasiNewton method. The explicit form of the exact Jacobian, A, of (22) reads as



M A= 1t ωK

 −γ K − J (F (c)) . M + 1t W

(24)

We now follow the framework from [40,39] and approximate A in two steps: (S1) First, approximate A by A0 ,





−γ K

M 1t ωK

A0 =

M

.

(25)

(S2) Second, permute A0 symmetrically to a block-diagonal form

 P T A0 P =



M 1t ωK



A0

=

A0



−γ K M M 1t ωK

 . −γ K  M

The block A0 admits exactly the form, analysed in [40,39]. Therefore, A0 is further approximated by a matrix in the form:

  A0 =

M 1t ω K

−γ K M + 2 βK

 (26)

where β = 1t ωγ . Finally,  A0 is replaced by its exact factorization

  A0 =

M 1t ωK

0

M+



βK

I 0

 1 −γ M − K , M −1 ( M + β K )

(27)

that is utilized when solving systems with A and leads to optimal order computational complexity. Based on the analysis in [40,39], the validity of step (S2) is straightforward. Since the spectrum of ( A0 )−1 A0 is contained  in the interval [0.5, 1], A0 is an optimal preconditioner for A0 . Although this result is already derived in [49], see also [40], for completeness of the presentation we include it in an Appendix. We next justify (S1), namely, we show that for small enough time steps, A0 is a good approximation of A. Consider the generalized eigenvalue problem

Aq = λA0 q.

(28)

We transform it to

(A − A0 )q = µA0 q,

(29)

where µ = λ − 1. Due to the form

 A − A0 =

0 0

 −J , 1t W

half of the eigenvalues of (29) are equal to zero and correspond to eigenvectors q = (q1 , 0)T . Therefore (28) has the same number of eigenvalues that are equal to 1. From (29) we see also that

−J q2 = µ(Mq1 − γ K q2 ) 1t W q2 = µ(1t ωK q1 + Mq2 ),

(30)

and if we express q1 from the first equation in (30) and substitute it in the second one, after a series of transformations, we find that the other half of the eigenvalues µ are solutions to the system

Qq2 = µq2 , Q = 1t (I + 1t γ ωM −1 KM −1 K )−1 (M −1 W + ωM−1 KM−1 J ).

114

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121 1

1

 = M 2 QM− 2 . Then Consider the similarity transformation Q  = 1t (I + 1t γ ωK 2 )−1 (W  + ωK J ), Q 1

1

1

1

1

1

 = M− 2 KM− 2 is a symmetric and positive definite matrix, W  = M− 2 W M− 2 and J  = M− 2 JM− 2 . where K Then there holds that ∥ ≤ 1t ∥W ∥ + ∥Q

1



2

1t ω ∥, ∥J γ

(31)

where we utilize the inequality |a|/(1 + a2 b2 ) ≤ 1/(2b) for b > 0. ∥. Consider in detail the entries of J (F (c)). We first analyse the term ∥J

 ∂ F1,1 (c)  ∂ c1,1  ..   .   ∂ F1,N (c)   ∂ c1,1 J (F (c)) = J =   ∂ F2,1 (c)   ∂c 1 ,1   ..   .  ∂ F (c) 2,N ∂ c1,1

∂ F1,1 (c) ∂ c1,N .. . ∂ F1,N (c) ∂ c1,N ∂ F2,1 (c) ∂ c1,N .. . ∂ F2,N (c) ∂ c1,N

··· ..

.

··· ··· ..

.

···

∂ F1,1 (c) ∂ c2,1 .. . ∂ F1,N (c) ∂ c2,1 ∂ F2,1 (c) ∂ c2,1 .. . ∂ F2,N (c) ∂ c2,1

··· ..

.

··· ··· ..

.

···

∂ F1,1 (c)  ∂ c2,N   ..   .  ∂ F1,N (c)   ∂ c2,N  . ∂ F2,1 (c)   ∂ c2,N    ..   . ∂ F (c)  2,N

∂ c2,N

The particular expression for the entries can be further worked out. Denote

ρij (c) =

∂ Fi (c) ∂ cj

(32)

and then



ρ11 (c)χ1 χ1

 Ω  ..   .   ρ (c)χ χ  11 N 1  J (F (c)) =  Ω   ρ21 (c)χ1 χ1  Ω  ..   .  Ω

ρ21 (c)χN χ1

 =

 ··· ..

Ω

. 

··· Ω ··· ..

Ω

. 

··· Ω

ρ11 (c)χ1 χN .. . ρ11 (c)χN χN ρ21 (c)χ1 χN .. . ρ21 (c)χN χN

 Ω

 Ω Ω

 Ω

ρ12 (c)χ1 χ1 .. .

 ··· ..

Ω

.

ρ12 (c)χN χ1

···

ρ22 (c)χ1 χ1

···

.. . ρ22 (c)χN χ1

..

.

···

ρ12 (c)χ1 χN



      ρ12 (c)χN χN    Ω   ρ22 (c)χ1 χN   Ω  ..   .   ρ22 (c)χN χN .. .

Ω



J11 J21

J12 . J22

(33)

After the congruence transformation of J we obtain

 = J

1

1

M − 2 J12 M − 2

1

1

M − 2 J22 M − 2

M − 2 J11 M − 2 M − 2 J21 M − 2

1

1

1

1

 .

From (33) it follows, that the blocks Jij , i, j = 1, 2, resemble mass matrices with coefficients that multiply their elements. The following estimate holds true

  1  − 12 M Jij M − 2  ≤ max |ρij (c)|.

(34)

c∈G

Due to the form of Fi (c), the functions ρij (c) are bounded for any c ∈ G. For the case (21),

ρij (c) = 12ξ (ci2 + 2ci cj − ci ),

ρii (c) = 12(6ci2 − 6ci + 1) + 12ξ (cj2 + 2ci cj − cj ).

For example, for ξ = 3.79, as was used in some of the numerical experiments, and since c1 + c2 ≤ 1, we obtain the bound |ρij (c)| ≤ 12. The estimate (34) is independent of the problem and the discretization parameter. We see, thus, that

∥ = O(1). ∥J

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

115

Recall next that γ = (3/4)Cn2 and ω = 1/Pe. In order to insure sufficient resolution of the discretization mesh, Cn = rh and Pei h = ζ have to hold for some r ≥ 2 and ζ < 1. Then, by choosing a time step 1t, that is small enough, compared to ∥ in (31) can be made arbitrary close to zero and thus all eigenvalues λ of (28) will cluster around one. h, the norm ∥Q We collect the results in a proposition. Proposition 1. Consider the matrices A and A0 as defined in (24) and (25), correspondingly. Then, for the eigenvalues {λi }2N i=1 of the generalized eigenvalue problem

Aq = λA0 q the following holds:

λi =



1 1+δ

i = 1, . . . , N i = N + 1, . . . , 2N ,

(35)

where the scalar δ can be made arbitrarily close to zero by choosing the time step 1t small enough with respect to the space discretization parameter h, namely, if 1t < ξh, ξ ≪ 1. Here N is the number of degrees of freedom for the concentration. Remark 4.1. Note that by definition β = 1t ωγ . Since γ < 1, it holds that β < 1γt ω . Thus, the choice of 1t to assure small values of δ guarantee that β < 1 is also small. 4.4. Generalization to problems with arbitrary number of components The approach, outlined in Section 4.3, can be straightforwardly generalized to multi-phase problems with arbitrary number of components n, when the gradient free energy is of the form (6), or (7) with a diagonal Υ , and the mobility matrix L is diagonal. For such problems, when n ≥ 3, the discrete multi-component system has the structure (22) and all three matrices M , K , W are (n − 1) × (n − 1) block-diagonal matrices, where each diagonal block is of dimension, equal to the sum of the number of degrees of freedom for ci and ηi . In order to approximate the exact Jacobian of a discrete (n − 1)-component system as in step (S1) by a matrix of the form (25), the Jacobian of the nonlinear term F (c) has to be bounded. That is true for any choice of a polynomial free energy, as well as for logarithmic free energy, regularized as in e.g. [11]. Step (S2) is then applicable to each of the diagonal blocks. Remark 4.2. In the case of more general forms of the free energy density and/or the mobility matrix, the off-diagonal blocks K in (24) might no longer be diagonal. In that case, in order to allow for the use of the presented preconditioning approach, an additional approximation step needs to be performed. Such studies, however, are a topic of further research and are not discussed in the present paper. 4.5. Solution algorithm We now outline the method for the solution of multi-component problems. At each time step k = 1, 2 . . . , we compute approximations of ck = ((ck1 )T , . . . , (ckn−1 )T )T and ηk = ((ηk1 )T , . . . , (ηkn−1 )T )T using a quasi-Newton method. The algorithm reads as follows (see also [39]): Given c0 = c0 , η0 = η0 , For k = 1, . . . , T /δ t Set ck,0 = ck−1 , ηk,0 = ηk−1 , For s = 0, 1, 2 · · · do until convergence For i = 1, . . . , n k,s k,s psi = −M ηi + Fi (ck,s ) + γ K ci k,s k,s k,s s qi = −1t ωK ηi − Mci − 1tW ci − Mcki −1 s s s zi = ρ pi + qi Solve Z gsi = zsi wsi = Mgsi − ρ psi Solve Z δsc ,i = wsi δsη,i = ρ1 (gsi − δsc ,i ) Update csi +1 = csi + ςs δsc ,i , ηs+1 = ηs + ςs δsη,i

End End End

√ β K , and ςs ∈ (0, 1] is a damping parameter. In our experiments, due to the fast convergence of the nonlinear solver, we use ςs = 1. Here ρ =



1t ω ,Z γ

=M+

116

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

Fig. 1. Initial condition for the three phases.

√

In the proposed algorithm, only solutions with matrices of the form M + β K , matrix–vector multiplications and vector updates need to be performed. The matrix Z is symmetric and positive definite and can be solved using freely available, highly optimized program toolboxes, for sequential, as well as parallel execution, that lead to optimal computational complexity. The other operations are BLAS primitives. Moreover, the computations for the different phases can be decoupled, providing an additional level of parallelism, where communication occurs only for transmitting local contributions when evaluating the terms Fi (ck,s ). Remark 4.3. The ideas, presented in this paper, are applicable for the more general θ -method for any value θ ∈ (0, 1], as well as for adaptive time stepping algorithms and numerical schemes that consider splitting of the nonlinear term over different time steps. 5. Numerical illustrations We test now the proposed solution method numerically using a three-component C–H model with a polynomial free energy, namely, the problem of spreading of a liquid lens between two stratified phases in two dimensions. This is a benchmark problem from [12], and here we adopt the same problem setting. The computational domain is taken to be Ω = [−0.3, 0.3] × [−0.15, 0.15] and the initial conditions are

 1 + tanh min(|x| − 0.1, x2 ) , c1 (x, 0) = 2 ϵ    2 1 1 − tanh max(−|x| + 0.1, x2 ) , c2 (x, 0) = 2 ϵ c3 (x, 0) = 1 − c1 (x, 0) − c2 (x, 0), 1





2

where x = (x1 , x2 ) ∈ Ω , see Fig. 1. In our experiments the interface width is taken as ϵ = 10−2 and the mobility is L0 = 10−4 . As already mentioned, the surface tensions σ12 , σ23 , σ13 determine the character of the process. Here, we examine both the case of partial and total spreading. Problem 1 (Partial Spreading). Consider the problem (17)–(15) with σ12 = 1 , σ23 = 0.8 , σ13 = 1.4. In this case all spreading coefficients are positive and in (15) we take λ = 0. As discussed in [12], the corresponding model is well-posed. Problem 2 (Total Spreading). Consider the problem (17)–(15) with σ12 = 1 , σ23 = 1 , σ13 = 3. In this case one of the spreading coefficients becomes negative and in order to ensure stability of the model, λ has to be large enough, see [12]. In the presented experiments λ = 7. The dynamics of these two flow problems are driven by minimization of the free energy. Starting from a given initial condition, the system seeks to attain an equilibrium state, determined by the model parameters. In the case of partial spreading, the equilibrium state for the lens in the limit ϵ → 0 can be derived analytically. The contact angles at triple junctions, see Fig. 2, are shown to depend on the three surface tensions by the Young relations, namely sin α1

σ23

=

sin α2

σ13

=

sin α3

σ12

.

For Problem 1, the evolution of the interface is shown in Fig. 3. In the case of total spreading, the third phase tends to spread in between the other two. The computed evolution of the interfaces for Problem 2 is illustrated in Fig. 4. All numerical experiments in this article are performed in Matlab. The aim of the conducted numerical experiments is to study the convergence of the nonlinear algorithm, presented in Section 4.5. The C–H equations are discretized using linear conforming finite elements on a regular triangular mesh. We consider a sequence of mesh refinements, where the coarsest triangulation corresponds to a characteristic size h = 0.3/127 and N = 33 153 mesh points, and the finest triangulation corresponds to h = 0.3/512 and N = 525 825. Note, that we

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

117

Fig. 2. Contact angles at the triple junction at equilibrium.

(a) t = 0.23.

(b) t = 5.86. Fig. 3. Partial spreading: evolution of the interface position.

solve a system for two of the concentration components only and corresponding chemical potentials, thus the size of the discrete system is 4N. We use stable time stepping schemes for the discretization of the C–H system. For Problem 1 it suffices to apply the backward Euler method in time and a fully implicit representation of the nonlinear term. As noted in [43], in the case of total spreading dynamics, the existence and convergence when using a fully implicit discretization of the nonlinear term needs further consideration. There, a semi-implicit splitting for the nonlinear term is proposed, that leads to a stable time stepping scheme, and decrease of energy, existence and convergence of the solution are shown for any time step. The splitting is based ∂f on an approximation (dFi ) of the term ∂ c0 in (20), dFi

(c , c m

m+1

)=

υi 4

i

[

+1 cm i

+

cm i

+1 +1 2 m 2 ][(cm + cm ) + (cm j + ck ) ] + j k

υk

υj 4

+1 2 m+1 2 m [(cm ) + (cm + ckm+1 + cm j ) ][ci i + ck ] j

m+1 +1 2 2 m + cjm+1 + cm ) + (cm [(cm k ) ][ci i + cj ] k   1 m+1 m 2 1 m m+1 2 m+1 m+1 m+1 2 m m m 2 + λ(ci + ci ) (cj ck ) + (cj ck ) + (cj ck ) + (cj ck ) ,

+

4

2

2

where cm and cm+1 are the concentration vectors for two consecutive time steps. The framework, presented in Section 4.3 applies successfully for the resulting discrete problem. As stated in Section 4.3, in order for the constructed approximate Jacobian to be a good preconditioner to the exact one, and thus, to lead to fast convergence of the quasi-Newton method, the time step has to be chosen small enough, compared to the spatial discretization step. In order to quantify this constraint for the case of Problems 1 and 2, we perform experiments for 1t equal to h/4 and h/10. Some results for the obtained iteration counts for Problems 1 and 2 are shown in Tables 1 and 2, correspondingly. The iteration counts are averaged over 50 time steps. At each time step, convergence is achieved when the norm of the update δs is smaller than 10−6 . The corresponding number of nonlinear iterations is presented in the table cells as the first number in the pair It nonlin /It lin . The second number represents the average linear iterations count for solving one subsystem with a matrix Z . Those matrices are symmetric positive definite, and for small values of the coefficient β (see Remark 4.1), exhibit properties, similar to those of a mass matrix alone. We have observed in numerical experiments, see also [39], that in many cases it suffices to apply the unpreconditioned conjugate gradient (CG) method to solve the corresponding subsystems. Moreover, to a great extent, the accuracy of the inner linear solver does not affect the convergence properties of the outer nonlinear solver. Table 1 shows a typical behaviour of the convergence of the solution process. There, we show the iterations to solve systems with the matrix Z via the unpreconditioned CG method with a stopping tolerance 10−3 and 10−1 . As we see from the results, an inner stopping tolerance of 10−1 does not lead to a significant change in the behaviour of the nonlinear solver. Experiments, not included in this paper, show that the use of an inner tolerance larger than 10−1 leads to an increase of the number of nonlinear iterations and, therefore, is not recommended. Based on the above, we choose to use 10−1 for the rest of the experiments.

118

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121 Table 1 Partial spreading: number of nonlinear (outer)/linear (inner) iterations.

1t

System size 132 612

527 364

2 103 300

19/18 18/14

20/30 19/23

20/49 20/38

18/8 18/6

19/12 19/10

20/21 20/16

20/3 19/2

20/3 20/3

CG, tol. 10−3 h/4 h/10 CG, tol. 10−1 h/4 h/10

AGMG, tol. 10−1 h/4 h/10

19/2 18/2

(a) t = 0.0012.

(b) t = 0.045.

(c) t = 0.28.

(d) t = 3.94. Fig. 4. Total spreading: evolution of the interface position.

Table 2 Total spreading: number of nonlinear (outer)/linear (inner) iterations.

1t

System size 132 612

527 364

2 103 300

45/8 27/6

30/12 27/10

29/22 27/17

31/2 28/2

30/3 27/3

CG, tol. 10−1 h/4 h/10

AGMG, tol. 10−1 h/4 h/10

61/2 28/2

Using the unpreconditioned CG method as an inner solver does not guarantee an optimal overall complexity, since the number of linear iterations may grow with the system size. However, the increase is not pronounced and unpreconditioned CG becomes a viable option, provided the simplicity of its implementation in serial and parallel, and its minimal computational and communication requirements. The last group of results, presented in Table 1, are obtained when systems with Z are solved using the aggregation-based algebraic multigrid AGMG method. The solution involves the use of flexible conjugate gradient iterations method and one Gauss–Seidel pre- and post-smoothing step in the AMG preconditioner at each level, see [50–53] for more details. When an optimal linear solver is applied to the inner subsystems, the overall solution procedure exhibits optimal computational complexity. As expected, in this case both inner and outer iterations stabilize independently of the number of degrees of freedom.

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

119

For Problem 2, when total spreading is present in the process, the solution turns out to be more computationally challenging. As seen in Table 2, the number of nonlinear iterations, needed for the quasi-Newton method to converge, is generally larger. However, they stabilize and are fewer when the space and time steps provide better resolution. Remark 5.1. As can be seen from the algorithm in Section 4.5, the proposed solution method for three-face models is an upgrade of that for two-phase problems, proposed earlier in [40,39]. Here we only study the convergence of the nonlinear and linear iterations. Detailed performance studies in terms of overall execution times can be found in [39], for both the sequential and the parallel implementation of the solution algorithm. We note that for the considered model problems, when constant time stepping is used, the number of Newton iterations are larger in the beginning of the solution process, and then gradually decreases with time. For Problem 1, the decrease of the nonlinear iterations is monotone. When solving Problem 2, convergence rates have the same general behaviour, except for a short period of time around t = 0.045 when we observe an increase of the number of nonlinear iterations. That is a sign that the time discretization parameters are not small enough to resolve the dynamics of the process, in that case, the appearance of a second interface, perpendicular to the boundaries. A slow down in the nonlinear convergence can be used to detect situations, when the dynamics of the solution undergo a change. In such situations one can use an adaptive time-stepping to better capture the behaviour of the solution. The quasi-Newton method, presented in this article, allows for straightforward inclusion of such techniques. For instance, a smaller time step is used to obtain the solution in Fig. 4(b). 6. Conclusions In this work we present an efficient solution method for multi-component problems, modelled by the Cahn–Hilliard equation. The technique is based on a quasi-Newton method, where the exact Jacobian is replaced by a high-quality approximation. The numerical solution method is a generalization of the two-phase method from [40,39]. The algorithm consists of matrix–vector operations and solutions of systems with a symmetric positive definite matrix, where off-theshelf software modules are directly applicable. The efficiency of the method is demonstrated by numerical experiments for a three-component Cahn–Hilliard model. The framework is applicable to various important classes of multiphase models and for an arbitrary number of components. Acknowledgements The work of the authors has been partly funded by the Swedish Research Council (VR) via grant VR 2008-5072 ‘Finite element preconditioners for algebraic problems as arising in modelling of multiphase microstructures’. The first author has been also partly funded by the Bulgarian NSF Grant DMU 03/62. The support is hereby gratefully acknowledged. Appendix Consider the matrices

−α 2 F T γ 2H

 A=

H β 2F



 −α 2 F T , γ 2 H + αβγ (F + F T )

 and

B=

H β 2F

where H is spd, the symmetric part of F , 12 (F + F T ) is positive semidefinite and α, β and γ are some real constants. Then, the eigenvalues of the generalized eigenvalue problem

λB

  v1 v2

  =A

v1 v2

(36)

satisfy the relations:



λ ∈ [0.5, 1], λ = 1 for v2 in the null space of F + F T .

Proof. We note first, that



0 A=B− 0



0 . αβγ (F + F T )

Then, problem (36) is equivalent to

(λ − 1)B

  v1 v2

 =

0 0

0 −αβγ (F + F T )

 

v1 , v2

120

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

i.e.,

(λ − 1)(Hv1 − α 2 F T v2 ) = 0

(37)

(λ − 1)(β 2 F v1 + (γ 2 H + αβγ (F + F T ))v2 ) = −αβγ (F + F T )v2 .

Clearly, if v2 is in the null space of F + F T , then λ = 1. Next, we consider λ ̸= 1 and express v1 = α 2 H −1 F T v2 from the first equality in (37). Substituting it in the second equality, we obtain

(λ − 1)(β 2 α 2 FH −1 F T + γ 2 H + αβγ (F + F T ))v2 = −αβγ (F + F T )v2 . We combine the similar terms, and after straight-forward transformations obtain



1

λ

 − 1 (β 2 α 2 FH −1 F T + γ 2 H )v2 = αβγ (F + F T )v2 .

Since H is spd, we transform the latter relation as



1

λ



− 1 (I +  F F T ) v2 = ( F + F T ) v2 , αβ

1

1

1

where  F = γ H − 2 FH − 2 and v2 = H 2 v2 . Hence, 1

λ

v2 ( F + F T ) v2 T

−1=

v2 v2 + ( F T v2 )T ( F T v2 ) T

.

We now utilize the assumption that F + F T is positive semidefinite together with the Cauchy–Schwarz inequality, and obtain: T T T T 0 ≤ v2 ( F + F T ) v2 = ( F T v2 )T v2 + v2 ( F T v2 ) = 2v2 ( F T v2 ) ≤ v2 v2 + ( F T v2 )T ( F T v2 ).

The latter shows that 0 ≤ λ1 − 1 ≤ 1, or 1 ≤ λ1 ≤ 2, that is

1 2

≤ λ ≤ 1.



References [1] I. Singer-Loginova, H.M. Singer, The phase field technique for modeling multiphase materials, Reports on Progress in Physics 71 (2008) 106501. [2] A. Novick-Cohen, The Cahn–Hilliard equation, in: C. Dafermos, M. Pokorny (Eds.), Handbook of Differential Equations. IV Evolutionary Partial Differential Equations, Elsevier, 2008, pp. 201–228. [3] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow, Annual Review of Fluid Mechanics 31 (1999) 567–603. [4] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi, W. Tauber, J. Han, S. Nas, Y.-J. Jan, A front tracking method for the computations of multiphase flow, Journal of Computational Physics 169 (2001) 708–759. [5] J.A. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, Cambridge University Press, 2000. [6] S. Osher, R.P. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, in: Applied Mathematical Sciences, Springer, New York, 2003. [7] S. Osher, J.A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations, Journal of Computational Physics 79 (1988) 12–49. [8] J.W. Cahn, On spinodal decomposition, Acta Metallurgica 9 (1961) 795–801. [9] J.W. Cahn, J. Hilliard, Free energy of a nonuniform system, I. Interfacial free energy, Journal of Chemical Physics 28 (1958) 258–267. [10] D.A. Cogswell, A phase-field study of ternary multiphase microstructures, Ph.D. Thesis for obtaining the degree of Doctor of Philosophy in Material Science and Engineering, MIT, 2010. http://www.dancogswell.com/wp-content/uploads/2010/02/Cogswell-Thesis.pdf. [11] C.M. Elliott, S. Luckhaus, A generalized diffusion equation for phase separation of a multi-component mixture with interfacial free energy, Preprint Series 887, IMA, University of Minnesota, 514 Vincent Hall, 206 Church Street S.E., Minneapolis, Minnesota 55455. [12] F. Boyer, C. Lapuerta, Study of a three component Cahn–Hilliard flow model, ESAIM: Mathematical Modelling and Numerical Analysis 40 (2006) 653–687. [13] M. Do-Quang, G. Amberg, The splash of a ball hitting a liquid surface: numerical simulation of the influence of wetting, Physics of Fluids 21 (2009) 022102. [14] B. Gustafsson, High Order Difference Methods for Time Dependent PDE, in: Springer Series in Computational Mathematics, vol. 38, Springer-Verlag, Berlin, 2008. [15] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Convection–Diffusion–Reaction and Flow Problems, second ed., in: Springer Series in Computational Mathematics, vol. 24, Springer-Verlag, Berlin, 2008. [16] N. Moelans, B. Blanpain, P. Wollants, Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems, Physical Review B 78 (2008) 024113. [17] A. Karma, W.J. Rappel, Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics, Physical Review E 53 (1996) 3017. [18] J.W. Barrett, J.F. Blowey, Finite element approximation of an Allen–Cahn/Cahn–Hilliard system, IMA Journal of Numerical Analysis 22 (2002) 11–71. [19] J.W. Barrett, J.F. Blowey, H. Garcke, On fully practical finite element approximations of degenerate Cahn–Hilliard systems, ESAIM: Mathematical Modelling and Numerical Analysis 35 (2001) 713–748. [20] J.F. Blowey, C.M. Elliott, The Cahn–Hilliard gradient theory for phase separation with non-smooth free energy part 1: mathematical analysis, European Journal of Applied Mathematics 2 (1991) 233–280. [21] F. Boyer, A theoretical and numerical model for the study of incompressible mixture flows, Computers & Fluids 31 (2002) 41–68. [22] Q. Du, R.A. Nicolaides, Numerical analysis of a continuum model of phase transition, SIAM Journal on Numerical Analysis 28 (1991) 1310–1322. [23] C.M. Elliott, The Cahn–Hilliard model for the kinetics of phase separation, in: J.F. Rodrigues (Ed.), Mathematical Models for Phase Change Problems, in: International Series of Numerical Mathematics, vol. 88, 1989, pp. 35–73. [24] C.M. Elliott, A.M. Stuart, The global dynamics of discrete semilinear parabolic equations, SIAM Journal on Numerical Analysis 30 (1993) 1622–1663. [25] X. Feng, Fully discrete finite element approximations of the Navier–Stokes–Cahn–Hilliard diffuse interface model for two-phase fluid flows, SIAM Journal on Numerical Analysis 44 (2006) 1049–1072.

P. Boyanova, M. Neytcheva / Computers and Mathematics with Applications 67 (2014) 106–121

121

[26] X. Feng, Y. He, C. Lei, Analysis of finite element approximations of a phase field model for two-phase fluids, Mathematics of Computation 76 (2007) 539–571. [27] H. Garcke, B. Stinner, Second order phase field asymptotics for multi-component systems, Interfaces and Free Boundaries 8 (2006) 131–157. [28] J.W. Barrett, J.F. Blowey, An improved error bound for a finite element approximation of a model for phase separation of a multi-component alloy, IMA Journal of Numerical Analysis 19 (1999) 147–168. [29] J.F. Blowey, M.I.M. Copetti, C.M. Elliott, The numerical analysis of a model for phase separation of a multi-component alloy, IMA Journal of Numerical Analysis 16 (1996) 111–139. [30] H. Garcke, B. Nestler, B. Stoth, A multi phase field concept: numerical simulations of moving phase boundaries and multiple junctions, SIAM Journal on Applied Mathematics 60 (2000) 295–315. [31] J. Kim, J. Lowengrub, Phase field modeling and simulation of three-phase flows, Interfaces and Free Boundaries 7 (4) (2005) 435–466. [32] J. Kim, K. Kang, J. Lowengrub, Conservative multigrid methods for ternary Cahn–Hilliard systems, Communications in Mathematical Sciences 2 (2004) 53–77. [33] H. Garcke, B. Nestler, B. Stoth, On anisotropic order parameter models for multi-phase systems and their sharp interface limits, Physica D 115 (1998) 87–108. [34] A. Miranville, G. Schimperna, Nonisothermal phase separation based on a microforce balance, Discrete and Continuous Dynamical Systems Series B 5 (3) (2005) 753–768. [35] G. Garcke, On Cahn–Hilliard systems with elasticity, Proceedings of the Royal Society of Edinburgh Section A 133 (2003) 307–331. [36] J.F. Blowey, C.M. Elliott, The Cahn–Hilliard gradient theory for phase separation with non-smooth free energy part 2: numerical analysis, European Journal of Applied Mathematics 3 (1992) 147–179. [37] J.W. Barrett, J.F. Blowey, An error bound for the finite element approximation of a model for phase separation of a multi-component alloy, IMA Journal of Numerical Analysis 16 (1996) 257–287. [38] M.I.M. Copetti, Numerical experiments of phase separation in ternary mixtures, Mathematics and Computers in Simulation 52 (2000) 41–51. [39] O. Axelsson, P. Boyanova, M. Kronbichler, M. Neytcheva, X. Wu, Numerical and computational efficiency of solvers for two-phase problems, Computers and Mathematics with Applications 65 (2013) 301–314. [40] P. Boyanova, M. Do-Quang, M. Neytcheva, Efficient preconditioners for large scale binary Cahn–Hilliard models, Computational Methods in Applied Mathematics 12 (2012) 1609–9389. http://dx.doi.org/10.2478/cmam-2012-0001. ISSN (Online). [41] J. Kim, A numerical method for the ternary Cahn–Hilliard system with a degenerate mobility, Applied Numerical Mathematics 59 (2009) 1029–1042. [42] F. Boyer, S. Minjeaud, Fully discrete approximation of a three component Cahn–Hilliard model, in: Proceedings of Algoritmy, 2009, pp. 143–152. [43] F. Boyer, S. Minjeaud, Numerical schemes for a three component Cahn–Hilliard model, ESAIM: Mathematical Modelling and Numerical Analysis 45 (2011) 697–738. [44] F. Boyer, C. Lapuerta, S. Minjeaud, B. Piar, M. Quintard, Cahn–Hilliard/Navier–Stokes model for the simulation of three-phase flows, Transport in Porous Media 82 (2010) 463–483. [45] A.A. Keller, M. Chen, Effect of spreading coefficient on three-phase relative permeability of NAPL, Water Resources Research 39 (2003) 1288. [46] V. Mani, K.K. Mohanty, Effect of the spreading coefficient on three-phase flow in porous media, Journal of Colloid and Interface Science 187 (1997) 45–56. [47] D. Jacqmin, Calculation of two-phase Navier–Stokes flows using phase-field modeling, Journal of Computational Physics 155 (1999) 96–127. [48] J. Kim, K. Kang, J. Lowengrub, Conservative multigrid methods for Cahn–Hilliard fluids, Journal of Computational Physics 193 (2004) 511–543. [49] O. Axelsson, M. Neytcheva, Operator splittings for solving nonlinear, coupled multiphysics problems with an application to the numerical solution of an interface problem, TR 2011-009, Institute for Information Technology, Uppsala University, 2011. [50] Y. Notay, AGMG software and documentation, http://homepages.ulb.ac.be/~ynotay/AGMG. [51] Y. Notay, An aggregation-based algebraic multigrid method, Electronic Transactions on Numerical Analysis 37 (2010) 123–146. [52] A. Napov, Y. Notay, An algebraic multigrid method with guaranteed convergence rate, SIAM Journal on Scientific Computing 34 (2012) A1079–A1109. [53] Y. Notay, Aggregation-based algebraic multigrid for convection–diffusion equations, SIAM Journal on Scientific Computing 34 (2012) A2288–A2316.