Grand-potential phase field simulations of droplet growth and sedimentation in a two-phase ternary fluid
Werner Verdier Alain Cartalade Mathis Plapp




































Grand-potential phase field simulations of droplet growth and sedimentation 
in a two-phase ternary fluid
    Werner Verdier^1,2, Alain Cartalade^1, Mathis Plapp^2[Corresponding author, mathis.plapp@polytechnique.fr]
    ^1 Université Paris-Saclay, CEA, Service de Thermo-hydraulique

et de Mécanique des Fluides, 91191 Gif-sur-Yvette, France

^2 Laboratoire de Physique de la Matière Condensée, CNRS,
École Polytechnique, 
Institut Polytechnique de Paris, 91120 Palaiseau, France
=======================================================================================================================================================================================================================================================================



§ ABSTRACT 
A methodology is built to model and simulate the dynamics of domain coarsening
of a two-phase ternary liquid with an arbitrary phase diagram. High numerical 
performance is obtained through the use of the phase field-method for interface
capturing, a lattice Boltzmann method numerical scheme for all the model equations, 
and a portable, parallel simulation code running on multiple GPUs. The model is 
benchmarked against an analytic solution for a ternary diffusion couple. 
It also reproduces the well-known power law for droplet coarsening during
Ostwald ripening without fluid flow. Large-scale simulations with flow illustrate 
the effects of momentum transport and buoyancy, as well as droplet coalescence
and sedimentation. 


Keywords: Phase-field models, Phase separation, Multicomponent mixtures, Lattice-Boltzmann method





§ INTRODUCTION

When a homogeneous solid or liquid mixture is rapidly quenched into a thermodynamically 
unstable state, phase separation occurs. After an initial stage of spinodal decomposition,
during which microscopic fluctuations are amplified with time, spatial domains of the new 
equilibrium phases form and are separated by well-defined interfaces with a characteristic thickness. 
At later times, domain coarsening occurs: under the driving force of capillarity, the total interface
area is progressively reduced by the elimination of geometric features with high 
curvature, such as protrusions or small domains.

This phenomenon has been extensively studied because it provides an example of a system 
that never reaches equilibrium, but instead exhibits scaling laws <cit.>. 
Indeed, in a system of infinite size, domain coarsening continues indefinitely 
and leads to patterns that are scale-invariant, that is, the structure and geometry 
of the pattern is statistically invariant in time, and only its overall scale grows 
with time as a power law. The growth exponent depends on the underlying transport 
processes (diffusion, hydrodynamics <cit.>, or elasticity)
and on the nature of the order parameter that describes the domain (scalar, vector or tensor,
conserved or non-conserved) <cit.>.

Besides its theoretical interest, phase separation is also important for a large number
of industrially relevant processes, such as the formation of porous glasses and membranes.
One example that motivated the present work is the conditioning of nuclear waste in glass: the radioactive substances are mixed with a glass matrix, a mixture of glass formers 
that is optimized for long-term resistance to environmental stresses <cit.>.

In order to increase the amount of waste, glass ceramics are being studied as an alternative
to conventional glass. For a simplified ternary system, a liquid-liquid phase separation can occur
when its global composition lies within a miscibility gap. The composition of the first liquid phase is representative of waste,
whereas the composition of the second one is representative of the glass matrix <cit.>.


For the study of generic features of phase separation, such as spinodal decomposition
dynamics and scaling laws in coarsening, the Cahn-Hilliard equation <cit.> and its 
various generalizations have been extensively used. Its appealing feature is that it can be
simply derived from out-of-equilibrium thermodynamics, taking as a starting point
the description of heterogeneous systems by a free-energy functional. The Cahn-Hilliard
equation can describe the complete time evolution of a phase-separating system, from
the initial spinodal decomposition to the late-stage coarsening, if transport is
governed by diffusion. It can be coupled with the Navier-Stokes 
equation <cit.> to model 
coarsening that is mediated by hydrodynamics, as well as the crossover between different 
coarsening regimes <cit.>.

While the Cahn-Hilliard equation is thus an excellent tool to explore generic features
of phase-separating systems, it is rather difficult to adapt it for the description
of specific materials. There are two reasons for this. Firstly, the Cahn-Hilliard
equation is formulated for a binary system, with a concentration as the only dynamic
variable, and a free-energy functional that consists of a double-well potential and a
gradient term. In this formulation, the interface free energy is controlled by the interplay 
of these two terms. It is therefore difficult to modify the bulk thermodynamic properties, 
which are set by the curvature of the free-energy function, while leaving the interface energy 
invariant <cit.>. Secondly, for multi-component
mixtures, there are several independent chemical compositions, and the gradient energy
coefficient becomes a symmetric matrix, with entries that are in principle determined
by the pairwise interaction between the different components. However, since those are
not directly known, the system is underdetermined, that is, for a target interface
energy, there are many possible choices for the gradient coefficients. This choice
has non-trivial consequences for the interface structure: interface adsorption of
certain components can occur, which modifies the interface energy and the macroscopic 
conservation laws for moving interfaces <cit.>.

An alternative way for describing interfaces in multicomponent mixtures has been
developed in the framework of phase-field theory <cit.>. 
The interface is described by a
scalar phase field, which can be seen as a smoothed indicator function (with 1
and 0 corresponding to presence or absence of the phase). The equation of motion
for this field derives from a Ginzburg-Landau free energy functional and is 
coupled to the concentration fields. The latter are then governed by their own 
free energy functions, which do not come into play for the determination of the 
interface properties, and which can therefore be chosen at will. The connection
between the phase-field model and the more traditional free-boundary problem can
be made using the well-established technique of matched asymptotic expansions
<cit.>. 

For the quantitative description of moving interfaces in multi-component mixtures, 
the grand-potential model
(which is a reformulation of the earlier Kim-Kim-Suzuki model <cit.>) 
has been widely used and benchmarked <cit.>. 
It corresponds to a grand-canonical formulation of the
mixture thermodynamics, with the chemical potentials as basic variables, which 
is suitable because two phases that coexist in a mixture have different 
compositions but equal chemical potentials.

Here, we couple a grand-canonical phase-field model for multicomponent mixtures
with a description of fluid flow using the lattice-Boltzmann method 
(LBM) <cit.> to achieve 
a model for domain coarsening that can be used for specific substances with given 
thermodynamics (obtained, for example, from a CALPHAD database <cit.>). This model 
cannot describe the initial stages of spinodal decomposition, because the grand-potential 
formulation requires a monotonic relation between composition and chemical potentials, but it
provides a quantitative description of the late-stage domain coarsening regime. Since, in
real substances, the two phases generally do not have the same density, we include
buoyancy in the model, which makes it possible to describe sedimentation of droplets
of the heavier phase.

We develop an implementation in which all the equations, including the phase-field
model, are solved within the LBM formalism <cit.>. 
Indeed, while LBM was initially developed for fluid flow, it was soon recognized 
that it is a general method for the time integration of partial differential 
equations, and we have previously used this formulation for solidification 
problems <cit.>. The advantage is that the simulation
code is entirely formulated in the LBM framework, and therefore high-performance
algorithms developed for LBM can be used. In particular, LBM is well adapted
for GPU parallelization.

In the present manuscript, we will first describe the grand-potential phase-field
model we use, which is similar to previously published work. Specifically, we
model a two-phase three-component fluid. We will then detail our
numerical implementation, and present a detailed benchmark study on a ternary 
diffusion couple, which shows that non-trivial interface conditions can be accurately
resolved by our formulation. Then, we present simulations of domain coarsening,
in two dimension without fluid flow as a benchmark, and in three dimensions
with fluid flow, in a large system which contains many droplets. These simulations
are carried out with the high-performance computing code 
LBM_Saclay <cit.>. These illustrative simulations 
demonstrate that our code can be used to perform large-scale studies, and thus is a 
suitable tool for the investigation of practical problems.



§ MODEL




 §.§ Thermodynamics

We consider a ternary mixture, and denote by c^A(x⃗,t), c^B(x⃗,t),
and c^C(x⃗,t) the local molar fractions of the A, B, and C components.
Since c^A+c^B+c^C=1, only two concentration fields are independent; we choose
to eliminate c^C by writing c^C=1-c^A-c^B, and write the (Helmholtz) free energy
density of a homogeneous system as f(c^A,c^B). For a specific substance, this
free energy function can be typically obtained from a CALPHAD <cit.> database.
Such a database usually gives the molar Gibbs free energy, but we will suppose that
the molar volume is constant and independent of the composition, so that Gibbs and
Helmholtz free energies are equivalent. The variables that are thermodynamically
conjugate to the compositions are the diffusion potentials, that is, the difference 
of the chemical potentials of A and C and B and C, respectively,

    μ^α = V_a ∂ f/∂ c^α     (α=A,B)

where V_a is the atomic volume (the molar volume divided by Avogadro's number). 
This factor has been included to give the diffusion potentials (which will also be 
called chemical potentials in the following) their usual dimension of energy <cit.>.

For a phase-separating system, the “free energy landscape” f(c^A,c^B) has a 
double-well structure, with two convex regions separated by a concave one; the latter
corresponds to the compositions for which a homogeneous system is thermodynamically
unstable. At equilibrium, the system is in a coexistence of two phases, called
“phase 0” and “phase 1” in the following, which have different compositions,
each one being located in one of the wells of the free energy function.
Minimization of the free energy under the constraint of global mass conservation
yields the conditions for phase coexistence:

    V_a . ∂ f_0∂ c^α|_c^α,_0
    	= V_a . ∂ f_1∂ c^α|_c^α,_1 =
    	μ^α,, 
        f_0(c^A,_0, c^B,_0) - f_1(c^A,_1,
    		c^B,_1) 
            = 1V_a∑_α=A,Bμ^α,
    		(c^α,_0 - c^α,_1) ,

where μ^α,, c_0^α,, and c_1^α, are the chemical
potentials and the compositions of phase 0 and 1 at two-phase equilibrium.
Geometrically, this corresponds to a “common tangent plane” to the free energy landscape.
From this geometric picture, or from a simple counting of degrees of freedom (three equations
for four unknown equilibrium compositions), it is clear that the solution to these
equations is not unique. In the space of compositions, which is usually visualized in
the Gibbs simplex, the various possible equilibria are indicated by tie lines
(see Figure <ref> below).
The global equilibrium corresponds to the unique tie line which contains the
composition inventory (the average composition) of the closed system.

The conditions for phase coexistence can be conveniently reformulated in the 
grand-canonical framework for open system, in which the fundamental variables
are the diffusion potentials, and the relevant thermodynamic potential is the
grand potential,

    ω (μ^A, μ^B) = f(c^A(μ^A, μ^B), c^B(μ^A,μ^B))
    		 - 1V_a∑_α=A,Bμ^α c^α(μ^A, μ^B),

which is the Legendre transform of the free energy. In writing down this formula,
we have supposed that the functions μ^α(c^A,c^B) can be inverted to yield
c^α(μ^A,μ^B). This is guaranteed only if these functions are monotonous,
which corresponds to a convex free energy landscape. In the case of a phase-separating
system, the Legendre transform must be taken separately for each convex region in
composition space (see below for details), which yields two different grand potential
functions for the two phases, ω_0 and ω_1. Phase coexistence is then 
simply given by

    ω_0(μ^A,μ^B)=ω_1(μ^A,μ^B).

This defines a coexistence line in the space of the intensive variables (the two diffusion
potentials), each point of which corresponds to a tie line in composition space.

In general, the calculation of the Legendre transform in Eq. (eq:grandpotential)
cannot be performed analytically, since free energy functions typically involve both polynomials and logarithms of the concentrations. It can always be performed numerically, but in 
the case of interest here
we can exploit the fact that during the late stages of the phase separation 
process, the concentrations will be close to the global equilibrium concentrations.
Therefore, only the vicinity in concentration space of the reference equilibrium
is relevant, and we can perform a second-order Taylor expansion of the free energy
around the equilibrium concentrations for each phase, 
as was already done in previous works <cit.>,

    f_π(c^A,c^B)   =   f_π(c^A,_π,c^B,_π)+∑_αμ^α,/V_a(c^α-c^α,_π)
           + 1/2∑_α,β K^αβ_π(c^α-c^α,_π)(c^β-c^β,_π).

Here and in the remainder of the text, exponents in Greek letters
and the index π identify components A or B, and phase 0 or 1, respectively.
In Eq. (<ref>), we have used Eq. (eq_mudef) in the first order terms, and 
the second order coefficients are given by

    K^αβ_π=.∂^2 f/∂ c^α∂ c^β|_c^A,_π,c^B,_π.

In this quadratic approximation, Eq. (<ref>) yields a linear relation between the
chemical potentials and the compositions,

    μ^α_π(c^A,c^B)=μ^α,+K^α A_π(c^A-c^A,_π)
                                             +K^α B_π(c^B-c^B,_π).

This equation can easily be inverted to obtain c^α(μ^A,μ^B), which, together with Eq. (<ref>), yields an anaytic approximation for the grand potential.

The procedure outlined above is valid in the late stage of phase separation for arbitraty
free energy functions. Since, in the present contribution, we do not intend to model a
particular material, but are rather interested in benchmark calculations, we choose as
an example a simple model system, in which the matrix of second order coefficients is 
diagonal, with equal values for the two phases,

    K_π^αβ = K δ_αβ,

where δ_αβ is the Kronecker symbol. This corresponds to circular parabolic 
free energy wells with equal curvatures for both phases.
Furthermore, we will switch to the
new variables μ̃^α=μ^α-μ^α, and 
ω̃_π=ω_π(μ^A,μ^B)-ω_π(μ^A,,μ^B,), which amounts
to choosing the reference values for the chemical potentials and the grand potential.
In these variables, the grand potentials can be expressed as

    ω̃_π (μ̃^A, μ̃^B)=- 12 K V_a^2∑_α=A,B(μ̃^α)^2- 1V_a∑_α=A,Bμ̃^α c^α,_π.

The phase diagram obtained from this model with c^A,_0=c^B,_0=0.3 and
c^A,_1=c^B,_1=0.4 is displayed in Figure <ref>.





 §.§ Grand potential phase field formalism

In a multi-component Cahn-Hilliard model, the free energy landscape defined above
would be supplemented by square gradient terms in the free energy functional
to capture the free energy cost of inhomogeneities, which then adds non-local terms
(containing Laplace operators) to the diffusion potentials. Instead, we choose to
describe phases and interfaces by an additional scalar phase field
φ(x⃗, t), with values in [0, 1], whose extrema identify one of
the two bulk liquid phase. The smooth variations of φ will locate the
diffuse interfaces.

The model is based on a phenomenological grand potential functional
<cit.> over the volume V of the system, 
function of the phase field φ and of the fields of intensive thermodynamical 
variables μ^A, μ^B,

    Ω[φ, μ^A, μ^B] = 
       ∫_V (ω_ int(φ)+ω_ th(φ,μ^A,μ^B)).

The two contributions describe the interfacial and bulk (thermodynamic) contributions to
the total grand potential, respectively.

The interfacial part,

    ω_ int(φ) = ζ2 |∇φ|^2 + H(φ),

contains a square gradient energy term and the double-well function (φ) =
8φ^2(1-φ)^2. They are parametrized by the
characteristic energy scales H ([E]  ·  L^-3, height of the double-well)
and ζ ([E]  ·  L^-1, gradient energy coefficient). The system
evolves to minimize Ω: the double-well term, favoring a sharp profile
for φ at the interface, will then play against the gradient energy
term, which favors smooth variations of φ.
Assuming a plane interface at equilibrium (where the bulk term ω_ th is
absent), the interface solution is given by the hyperbolic tangent profile

    φ_0(x) = 12( 1 + tanh(2x / W)),

with

    W = √(ζ / H)

the characteristic interface width. Reintroducing this profile in the integral
(<ref>) gives a characteristic value of the surface tension
due to the diffuse interface,

    σ = 23 HW .

This relation also identifies H as the characteristic energy-per-volume
scale of the excess free energy stored in the diffuse interface.
More details about these calculations can be found in introductory
texts about the phase-field method, for example in Refs. <cit.>.

The bulk contribution is an interpolation between the grand potentials of the
two phases,

    ω(φ, μ^A, μ^B) = [1-p(φ)] ω_0(μ^A, μ^B) + p(φ) ω_1(μ^A,μ^B) ,

with p(φ)=3φ^2-2φ^3 an interpolation function odd around φ = 1/2, 
satisfying p({0, 1}) = {0, 1} and p'({0, 1}) = 0, that can be seen as a smoothed
step function <cit.>.

The interface-tracking equation, also known as the Allen-Cahn equation, is obtained
by relating linearly the time evolution of φ and the decrease in Ω,
expressed as the variational derivative,

    ∂_t φ = - M_φδΩδφ.

Here, M_φ is a kinetic coefficient (phase-field mobility). The terms in this
equations that are generated from the interfacial energy stabilize the interface
profile, whereas the bulk term creates a driving force for interface motion if
the grand potentials of the two phases differ.

It is convenient to remove the energy dimensions: define the dimensionless 
chemical potentials μ^α = μ̃^α/ (K V_a), and

    ω_π (μ^A, μ^B)     = ω̃_π (μ̃^A,
    		μ̃^B)/K 
        = - 12∑_α=A,B( (μ^α)^2 +
    		2 μ^α c^α,_π).

We can now write the interface tracking equation (<ref>)
in terms of the fields φ, μ^A and μ^B as

    τ_φ∂_t φ = W^2 ∇^2φ - '(φ) + λ p'(φ) Δω(μ^A,
    		μ^B)

with the phase-field relaxation time τ_φ=1/(M_φ H) and the difference of the 
grand potential densities

    Δω    = ω_0 - ω_1 
        = - ∑_α=A,Bμ^α( c^α,_0 -
    		c^α,_1 ) ,

where the thermodynamical coupling parameter is λ = K / H. We also define the scaled
phase-field mobility M_φ=M_φ H=W^2/τ_φ, which has the dimension
of a diffusion coefficient.



 §.§ Species diffusion and mixed formulation


In the absence of hydrodynamic flow, the redistribution of chemical species 
occurs by diffusion. The concentration of each species is a locally conserved
quantity, hence

    ∂_t c^α(x⃗,t) = - ∇·j⃗^α(x⃗,t).

According to the thermodynamics of irreversible processes, the species currents write

    j⃗^α(x⃗,t)=-∑_β=A,B M^αβ∇μ^β,

where M^αβ are the components of the atomic mobility matrix.
Furthermore, in the grand-canonical setting, where the natural variables are
the chemical potentials, the composition can be obtained from the grand-potential
functional by

    c^α = - V_aδΩδμ^α = - V_a( p(1-φ)
    	∂ω_0∂μ^α + p(φ)
    	∂ω_1∂μ^α) .


In the original grand-potential method for a binary alloy <cit.>, the 
concentration was eliminated from Eq. (eq:conservation) in favor of the chemical 
potential; however, for multi-component alloys, this requires numerous matrix 
inversions <cit.>. Therefore, we prefer
to use a mixed formulation, as in some previous works <cit.>: we keep
both the concentrations and the compositions as dynamic variables.
This allows to solve Eq. (<ref>), using the explicit algorithm of
one's choice to update c^α ; μ^α is then updated by inverting
Eq. (<ref>). Remark that the inversion step has a simple
analytical expression with the quadratic free energies, namely

    μ^α = c^α - p(1 - φ) c^α,_0 - p(φ) c^α,_1 .

Here, we will use a diagonal mobility matrix
and remove the dimension of energy by defining M^α=KM^αα;
M^α then has the dimension of a diffusion coefficient.
If the diffusion coefficient is different in the two phases, the mobility is
linearly interpolated, M^α(φ)=M_0^α(1-φ)+M_1^αφ for both components α=A,B.
Notice that in the bulk phases, Eq. (<ref>) reduces
to Fick's laws with fluxes M^α∇ c^α because of the quadratic
free energies. This will be helpful in preliminary tests since analytical
solutions are available in this case ( sec. <ref>).



 §.§ Equivalent sharp-interface model



The fact that phase field models only track an implicit, diffuse interface
greatly facilitates the numerical treatment of two-phase problems. However,
the behavior of the physical variables other than φ at the phase 
interface is at first glance unspecified. This is in contrast to the
sharp-interface formulations of free-boundary problems, whose well-posedness
depends on the presence of explicit boundary condition at the interface
for values and fluxes of the relevant transport fields.

An important effort in the phase field literature is made to bridge both
formulations by extracting implicit interface conditions from the phase field
equations in the asymptotic limit of small W (“thin-interface limit"),
and to add corrections if necessary to obtain adjustable interface properties. 
This is done using the formalism of matched asymptotic analysis in terms of an asymptotic
parameter ε expressing the scale separation between the interface thickness
and the physically relevant “outer” scales. As this process
is lengthy and has already been presented in much detail in other works,
we will only present its results here. The interested reader can refer to
<cit.> for a step-by-step detail of the calculations ; to
<cit.> and <cit.> for a more detailed
definition of the curvilinear coordinates and subsequent reexpression of
the differential operators ; and to <cit.> and then
<cit.> for a discussion of what can be considered 
a “thin” interface.

The asymptotic analysis of the model without flow was realized
following the formalism of Almgren <cit.>. It was seen
that the case of a ternary system or the presence of a closure relation,
Eq. (<ref>), introduce little novelty in the calculations.
The result is a Gibbs-Thomson relation for the grand potential at the
interface,

    Δω = -δκ - β V

with κ the local curvature of the interface, V its normal
velocity, and the associated coefficients are

    δ    =2/3W/λ
    β    =2/3W/λ M_φ-19/120∑_α=A,BW(c_0^α,-c_1^α,)^2/M^α


The numerical factors are due to integrations across the diffuse interface
profile and depend on the choice of the double-well function and the interpolation
function p(φ). The asymptotics also yields the composition balance at the interface,

    V[c^α]_-^+=-[M^α∂_nμ^α]_-^+

where [ · ]^+_- denotes the jump of a field across the interface,
and ∂_n the gradient projected on the interface normal. 









 §.§ Flow coupling

To include hydrodynamic flow dynamics, we couple our phase-field model to the
incompressible Navier-Stokes equation,

    ∇·u    =0
    ρ_0[∂_tu+∇·(uu)]     =-∇p_h+∇·[ρ_0ν(φ)(∇u+∇u^T)]+F_tot
    ∂_tφ+u·∇φ    =W^2/τ_φ∇^2φ-1/τ_φω_dw^'(φ)+λ/τ_φp^'(φ)Δω(μ^A,μ^B)
    ∂_tc^α+u·∇c^α    =∇·[M^α(φ)∇μ^α]
    μ^α    =c^α-p(1-φ)c_0^α,-p(φ)c_1^α,

where u⃗ is a phase-averaged velocity field <cit.>, 
ρ_0 is the constant density, p_h is the hydrodynamic
pressure enforcing the condition ∇·u=0
(Lagrange multiplier of that condition), ν(φ) is the (phase-dependent) kinematic
viscosity, the component index α is A,B and the total force F_tot 
is the sum of the gravity force F_g=φΔρg
in the Boussinesq approximation (with Δρ the density different between the
phases, phase 1 being denser than phase 0, and g⃗ the constant 
acceleration vector due to gravity) and the surface tension 
force F_σ defined by <cit.>:


    F_σ=3/2σ W[ω_dw^'(φ)/W^2-∇^2φ]∇φ

For the equilibrium solution φ=φ_0, the term inside
the bracket is equivalent to κ|∇φ|
where κ is the curvature defined by κ=-∇·n.
The surface tension force can be expressed by F_s=δ_dσκn,
i.e. the quantity σκn is spread over δ_d=(3/2)W|∇φ|^2
where W is the thickness of the diffuse interface. 

In the phase-field equation Eq. (<ref>), for simplifying
the notations in the next Section, we define the following source
term:


    𝒮_φ(φ,μ^A,μ^B)=-1/τ_φω_dw^'(φ)+λ/τ_φp^'(φ)Δω(μ^A,μ^B)

which contains the contributions of the double-well (first term) and
the thermodynamic imbalance Δω (second term)
responsible for the displacement of the interface. For u=0,
the interface is displaced by diffusion until the thermodynamic equilibrium is
reached, i.e. when the grand-potential densities of each phase are
equal ω_0(μ^A,μ^B)=ω_1(μ^A,μ^B).



§ LATTICE BOLTZMANN METHOD

The numerical resolution of the model is done through
the LBM, which consists of a discretization of the Boltzmann equation
in phase-space. This defines a regular lattice, and the time-explicit
resolution of the equation is assimilable to the distribution functions
undergoing a collision step on the nodes followed by a transport step
along the edges. According to the choice of collision operator (here
the simple BGK operator), of equilibrium functions, and of additional
source and force terms, the moments of the distribution function thus
solved can be shown to be solution of conservative time-evolution
PDEs such as the ones of the present model. This is coherent with
a discretized Chapman-Enskog expansion.

LBM is a very popular method to simulate various conservative PDEs and more particularly the Navier-Stokes equations. For those latter ones, more traditional methods use a prediction-correction algorithm that requires solving a time-consuming Poisson equation. The LBM benefits from the efficiency of the equivalent articifial compressibility algorithm. The collision operator is local and each discrete distribution function follows an identical evolution equation. The method is therefore very simple to implement and is parallel in nature by efficiently exploiting the shared memory. Moreover, coupled with a distributed memory parallelism (e.g. MPI), that algorithm is very efficient. The interested reader can refer to references such as <cit.> for more details on the LBM.

In this work, the standard D2Q9 and D3Q19 lattices are used for all 2D or 3D simulations,
respectively. Those lattices correspond to a discretization of the
velocity-space with 9 or 19 discrete velocities c_k
defined by c_k=(δ x/δ t)e_k
with δ x the space step, δ t the time step, and e_k
the direction vectors. Each set of directions is weighted by a scalar
value w_k. The directions e_k and their weights
are listed in Table <ref> for D2Q9 and Table <ref>
for D3Q19.

We also define a characteristic speed c_s defined by c_s=(1/√(3))δ x/δ t.
The distribution functions and the discrete Lattice
Boltzmann equations (LBE) with the BGK collision operator are detailed
next in Section <ref> for fluid flow, in section
<ref> for phase-field and in Section <ref>
for the composition equations. Many alternative collision operators exist in LBM literature. The most popular operators are the "Two-Relxation-Times" (TRT) and "Multiple-Relaxation-Times" (MRT) which define respectively two (for TRT) or more (for MRT) additional collision rates to tune for improving stability and accuracy. Here, we present a proof of concept of a thermodynamic model that is based on the grand-potential functional and coupled with fluid flow. The ratio of diffusion coefficients are almost 1 and the kinematic viscosity of each phase is identical. The BGK operator is sufficient in this work.






 §.§ LBM for fluid flow


For simulating the fluid flow, the numerical scheme works on the distribution
function f_k(x,t) for which its evolution is given
by the lattice Boltzmann equation:


    f_k(x+c_kδ t,t+δ t)=f_k(x,t)-1/τ_f(φ)+0.5[f_k(x,t)-f_k^eq(x,t)]+δ tℱ_k(x,t)

where the relaxation rate τ_f is related to the kinematic
viscosity by ν(φ)=τ_f(φ)c_s^2δ t and
the source term ℱ_k(x,t) contains the
total force term F_tot of Eq. (<ref>). Several
forcing schemes exist in the LBM literature, here we choose one of
the two most popular methods <cit.>:


    ℱ_k(x,t)     =Γ_k(x,t)(c_k-u)·F_tot
    Γ_k(x,t)     =w_k[1+c_k·u/c_s^2+(c_k·u)^2/2c_s^4-u^2/2c_s^2]

The moment of order zero of Eq. (<ref>) is ∑_kℱ_k=0
and its moment of first-order is equal to ∑_kℱ_kc_k=F_tot.
For recovering the incompressible Navier-Stokes equation with the
Chapman-Enskog expansion, the equilibrium distribution function f_k^eq
has to be defined by <cit.>:


    f_k^eq(x,t)=w_k[p_h+ρ_0c_s^2(c_k·u/c_s^2+(c_k·u)^2/2c_s^4-u^2/2c_s^2)]-δ t/2ℱ_k

where the term ℱ_kδ t/2 has been substracted for
capturing the second-order accuracy when the external force term is
taking into account in Eq. (<ref>). That trick is
equivalent to a change of variables in the distribution function to conserve
an explicit algorithm after the trapezoidal integration of the Boltzmann
equation. Because of that change of variables, 
the term ∑_kℱ_kc_kδ t/2
must be added for updating the velocity for each time-step. After
collision and streaming, the hydrodynamic pressure is obtained by
the moment of order zero of f_k, and the velocity by the moment
of first order:


    p_h    =∑_kf_k
    ρ_0u    =1/c_s^2∑_kf_kc_k+δ t/2F_tot


When a high contrast of density exists between both phases, another
popular method is widely applied for defining the equilibrium distribution
function f_k^eq. In that case, the density becomes an interpolation
of the bulk density of each phase (e.g. ϱ(φ)), and
a dimensionless pressure p^⋆=p_h/ϱ(φ)c_s^2
is introduced in the equilibrium Eq. (<ref>) <cit.>. In that
case, two supplementary forces must be added in F_tot,
the pressure and viscosity forces <cit.>, in
order to recover the incompressible Navier-Stokes. Here, we assume
that the density is identical in both phases, so that the classical
equilibrium Eq. (<ref>) is sufficient for the simulations.

The surface tension force Eq. (<ref>) requires
computing a gradient term and a laplacian term. The gradient is evaluated
by the directional derivatives:


    e_k·∇φ|_x=1/2δ x[φ(x+e_kδ x)-φ(x-e_kδ x)]

where the number of directional derivatives is equal to the number
of moving directions e_k on the lattice i.e. N_pop.
The gradient is obtained by:


    ∇φ|_x=3∑_k=0^N_popw_ke_k(e_k·∇φ|_x).


For the calculation of the ∇^2φ, all
directions of propagation are taken into account by


    (e_k·∇)^2φ|_𝐱=1/δ x^2[φ(x+e_kδ x)-2φ(x)+φ(x-e_kδ x)]

which are used to compute the laplacian:


    ∇^2φ|_x=3∑_k≠0w_k(e_k·∇)^2φ|_x.





 §.§ LBM for phase-field


For the phase-field equation, we introduce a new distribution function
g_k(x,t) evolving with the LBE: 


    g_k(x+c_kδ t,t+δ t)=g_k(x,t)-1/τ_g+0.5[g_k(x,t)-g_k^eq(x,t)]+δ t𝒢_k(x,t)

Eq. (<ref>) is an Advection-Diffusion type Equation
(ADE). The equilibrium distribution function g_k^eq is designed
such as its moments of order zero, one, and two are respectively equal
to φ, uφ and Iφ
where I is the identity tensor
of second order:


    g_k^eq(x,t)=w_kφ[1+c_k·u/c_s^2]-δ t/2𝒢_k

The first term in the right-hand side of Eq. (<ref>) is the
classical term to recover the ADE after the Chapman-Enskog procedure,
and the source term 𝒢_k is simply defined such as its
moment of order zero is equal to 𝒮_φ(φ,μ^A,μ^B):


    𝒢_k=w_k𝒮_φ(φ,μ^A,μ^B)

where 𝒮_φ, defined by Eq. (<ref>).
The interface mobility M_φ=W^2/τ_φ is a constant which is related
to the relaxation rate τ_g by M_φ=τ_gc_s^2δ t.
Finally, after collision and streaming, the phase-field φ
is updated at each time-step by:


    φ=∑_kg_k+δ t/2𝒮_φ





 §.§ LBM for composition equations


For the treatment of the composition equations, we introduce two distribution functions
h_k^α(x,t) (for α=A,B) which evolve
with the LBE:


    h_k^α(x+c_kδ t,t+δ t)=h_k^α(x,t)-1/τ_h^α(φ)+0.5[h_k^α(x,t)-h_k^eq,α(x,t)]

Eqs (<ref>) look like advection-diffusion equations,
but the equilibrium functions h_k^eq,α must be defined
such that its moment of order two is equal to Iμ^α
where μ^α is related to the composition c^α
(moment of order zero) by Eq. (<ref>). In that relation,
c_0^α, and c_1^α, are two scalar input values which
are the two thermodynamic equilibrium compositions of each phase 0
and 1. Because the moments of order 2 differs from the moment of order
0, the equilibrium distribution function must be slightly modified
compared to the classical equilibrium for an ADE by:


    h_k^eq,α(x,t)=
    c^α-(1-w_0)μ^α   for k=0
    
    w_kμ^α+w_kc^αc_k·u/c_s^2   for k≠0

In that equilibrium distribution, the moment of order zero is given
by the first term of the first line and the moment of order one is
given by the second term of the second line. That method is extensively
used for simulating the advective Cahn-Hilliard equation (<cit.>).
The mobility coefficients M^α(φ) are related to the relaxation rate τ_h^α(φ) by M^α(φ)=τ_h^α(φ)c_s^2δ t.
After collision and streaming, the composition is obtained by the
moment of order zero of h_k^α:


    c^α=∑_kh_k^α





§ NUMERICAL SIMULATIONS


The methods of the previous section have been implemented in a C++ code: LBM_Saclay <cit.>. Its main feature is its multi-architecture portability by simple modifications
of compilation options in the makefile. This is possible thanks to the Kokkos library <cit.>. Thus the same code can be run either on multi-CPU partition or on multi-GPU partition of a supercomputer. Time performances of LBM_Saclay have been compared between graphic cards and standard CPU in our previous paper <cit.>. Because of the efficiency of the former ones, the simulations of this work have run on the multi-GPU partitions of the supercomputers Jean-Zay (IDRIS, France) and Topaze (CCRT, France).



 §.§ Ternary diffusion couple


We start by presenting a quantitative validation of the model and its
implicit Gibbs-Thomson condition for a plane interface without
flow. This demonstrates that the phase field model is able
to reproduce the known solution of a sharp-interface problem, the
symmetrical ternary diffusion couple. Phase field studies of this
problem already exist in the literature, and we use this as a test for our model <cit.>. An interface
splits the infinite one-dimensional space in two domains at x_I(t), with
phase 0 on its left and phase 1 on its right. The interface is displaced
at velocity ẋ_I(t) by the interdiffusion of the components. This 
problem can be stated as the one-dimensional free boundary problem


    ∂_tc^α    =M^α∂_xx^2μ^α,     x<x_I(t) or x>x_I(t)
    μ^α    =μ_int^α,     x=x_I(t)
    ẋ_I(c^α|_x_I^--c^α|_x_I^+)     =-M^α(∂_xc^α|_x_I^--∂_xc^α|_x_I^+)     x=x_I(t)
    μ^α    =μ_-∞^α    x=-∞
    μ^α    =μ_+∞^α    x=+∞

with the initial conditions



    x_I(0)=0;     μ^α=μ_-∞^α   (x<0);    μ^α=μ_+∞^α   (x>0).


The problem corresponds to evaluating our phase field model in 1D
with φ=0 or φ=1 and taking the interface conditions
(<ref>), (<ref>) with β=0 and κ=0. 
In the sharp-interface view, the relation between μ^α and c^α
is linear, allowing the problem to be rewritten into the classical
Stefan problem in terms of either c^α or μ^α
only. It is then known <cit.> that a self-similar
solution exists with


    x_I(t)=ξ√(t)

and


    μ^α(x,t)=μ_-∞^α+(μ_int^α-μ_-∞^α)erfc(-x/2√(M^αt))/erfc(-ξ/2√(M^α))    -∞<x<x_I(t)
    μ_+∞^α+(μ_int^α-μ_+∞^α)erfc(x/2√(M^αt))/erfc(ξ/2√(M^α))    x_I(t)<x<+∞


The statement and solution of the problem in terms of the
composition fields is analogous but with a discontinuity at the interface,
in coherence with chemical equilibrium. The couple of interface values
(μ_int^A,μ_int^B) are a solution
of the system, meaning


    ω_0(μ_int^A,μ_int^B)-ω_1(μ_int^A,μ_int^B)=0.


A ternary system has access to a continuous set of such equilibrium
couples (phase diagram tie-lines). The additional dynamical constraints
select one particular equilibrium in this continuum: the coefficient ξ is determined
as the solution of a transcendental equation that couples the thermodynamic
equilibrium and dynamical parameters, see Eq. (22) in <cit.>.





1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 
3cDomain          3cPhase-field parameters          3cTransport parameters1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 11-11 
Symbol     Value     Dim          Symbol     Value     Dim          Symbol     Value     Dim[-L_x,L_x]     [-1,1]     [L]          W     1.2×10^-3     [L]          (c_0^A,,c_0^B,)     (0.3,0.3)     [–][-L_y,L_y]     [-0.002,0.002]     [L]          λ     155.94541910≅155.95     [–]          (c_1^A,,c_1^B,)     (0.4,0.4)     [–]N_x nodes     3000     [–]          M_φ     1.2     [L]^2/[T]          (M^A,M^B)     (1,0.8)     [L]^2/[T]N_y nodes     6     [–]                              (c^A,c^B)_-∞     (0.4,0.175)     [–]δ x     1/1500     [L]                              (c^A,c^B)_+∞     (0.225,0.6)     [–]δ t     1.481481×10^-9     [T]                                        

Parameters for the simulation
of symmetrical diffusion couple. λ^⋆ is the value that
cancels the kinetic coefficient β (Eq. (<ref>)).
In all simulations, that coefficient will be simply noted λ
and approximated by 155.95. The values of (c^A,c^B)_-∞
and (c^A,c^B)_+∞ are respectively equivalent to (μ^A,μ^B)_-∞=(0.1,-0.125)
and (μ^A,μ^B)_+∞=(-0.175,0.2).
The values of collision rates corresponding respectively to M_φ, M^A and M^B are: τ_g=0.6333, τ_h^A=0.6111 and τ_h^B=0.5888.



To reproduce this problem numerically, the phase-field model is solved
on a thin 2D domain (3000 × 6 LBM lattice nodes). The phase
field is initialized with the equilibrium hyperbolic tangent profile.
The composition fields are initialized as step functions. The parameters used are 
referenced in table <ref> and they produce an analytical
solution with ξ=-0.269824. Figure <ref>
compares the interface velocity and the concentration fields simulated with the ones 
expected from the analytical solution. Both are in excellent agreement, numerically
confirming the interface condition (<ref>)–(<ref>)
reconstructed by the phase field model. Note that we have chosen the system's half-length
L and the diffusion time t_D=min_α(L^2/M^α) as
units in these comparisons.






 §.§ Simulations of Ostwald ripening


Next, the phase-field model is exploited to simulate the Ostwald ripening
of a set of droplets (phase 1) inside a continuous matrix (phase 0). 
As the Allen-Cahn model cannot describe the initial regime of phase separation
from homogeneous mixtures, the initial condition starts from already developed droplets,
as detailed below. In this problem, the characteristic
length scale is the average spacing between droplets. However, this
length is not static and increases as the growth proceeds and the
droplet number decreases. In the numerical simulations, the distance between
droplets is bounded by the domain size. For this reason, we once again
set our numerical units of length and time as the half-length of the
domain and the diffusion time, respectively, ℓ=L and 
t_D=min_α(ℓ^2/M^α).



  §.§.§ Initial condition


As already mentioned above, to study the growth kinetics,
we must start from pre-existing droplets. The initial
condition must be carefully constructed to satisfy the condition for
growth, namely that the global composition inventory lies in the miscibility
gap and that the phase fraction be sufficiently
high to measure a statically relevant average. In addition, before
the growth, there will be a transient regime during which the droplets
reach local equilibrium with the surrounding matrix (analogous to
a diffusion couple). We want the droplets to always tend to grow during
this regime, and not to shrink and risk disappearance.

We therefore choose to parametrize the initialization of the composition
c(x,0) and the phase-field φ(x,0)
with three scalar values: the initial global composition of the system
c_g=(c_g^A,c_g^B) and the initial phase
fraction Φ of droplets. The phase-field φ(x,0)
can be easily initialized by creating spherical droplets with randomly
generated positions and radii. For the compositions, we could then
initialize them inside and outside those interfaces using the reference
tie-line’s end-points. However, specifying Φ,
c_g and the tie-line would overspecify the initialization
because of the lever rule. Besides, if the initial composition of
the matrix is at equilibrium c_0^eq, the droplets
can disappear during the transient regime. It is necessary to design
an initial condition on c^A(x,0) and c^B(x,0)
such that the droplets grow from the outset of the time evolution.

When the thermodynamic equilibrium is reached, because of the conservation
rule, the global composition c_g are related to
the equilibrium compositions and the equilibrium phase fraction Φ^eq
by:


    c_g=(1-Φ^eq)c_0^eq+Φ^eqc_1^eq

 That relation can be inverted and yields:


    Φ^eq=|c_0^eq-c_g|/|c_0^eq-c_1^eq|


One way to avoid droplet disappearance is to initialize the
system with a phase fraction of droplets Φ which is lower than
the equilibrium Φ^eq. In that case, we need to initialize
the composition in the matrix and/or in the droplets out of equilibrium
(with an additional degree of freedom due to the ternary case). Here,
the droplet compositions are considered at equilibrium c_1^eq
and we choose to supersaturate the matrix by offsetting its composition
along the tie-line, as


    c_0^ini=c_0^eq-δ(c_1^eq-c_0^eq)

where the coefficient δ is defined by:


    δ    =Φ^eq-Φ/1-Φ

where Φ^eq is given by Eq. (<ref>). When δ>0
(i.e. when Φ<Φ^eq), the matrix composition is offset
inside the miscibility gap. This ensures that the initial transient leading
towards local equilibrium always corresponds to an increase of the global 
phase fraction of droplets Φ. Finally, the initial conditions of compositions write


    c(x,0)=
    (1-δ)c_0^eq+δc_1^eq   if φ(x)<1/2
    c_1^eq   if φ(x)≥1/2

which represents an initial supersaturation in the matrix, whereas
the droplet composition is supposed to be at equilibrium.


The details of the initialization routine we
used are given next. For a specified target phase fraction s,
we initialize the droplet geometry according to the procedure illustrated in
Figure <ref> which results in lists of positions and
volumes for each droplet. Then, to initialize φ(x) on each
lattice node, search for the closest droplet where we define the hyperbolic
tangent profile using its radius and center. Next, the composition fields
are initialized on each lattice node as Eq. (<ref>) for a given composition
inventory c_g, followed by the chemical potential fields
at equilibrium everywhere, μ^α(x)=0. Finally,
for simulations with flow, the flow field is initialized at rest and the pressure
field at zero.






  §.§.§ Geometry measurements


To measure the droplet count N, we use an algorithm for the measure
of the Euler characteristic, to be understood here as the number of
connected regions of φ=1. The algorithm is adapted from reference
<cit.> and it only requires a single pass over the
lattice. The average radius ⟨ R⟩ is measured
using the integral estimate


    ∫4/Wφ(1-φ)d^2x   ≈∑_i2π∫_0^+∞4/Wφ_0(1-φ_0)rdr
       ≈2π W∑_i[R_i/W+𝒪(e^-4R_i/W)]
       ≈2π N⟨ R⟩

where the index i stands for the droplet number. Indeed, the integrand is a function that
has a sharp peak in the interfaces, and the integral approximately gives the total length
of interface present in the system. Assuming the droplets are spherical and sufficiently far
from each other, the integrals can be calculated exactly and a Taylor
expansion gives the desired result, plus an error terms in 𝒪(e^-4R_i/W);
the numerical calculation of this integral is easy to implement 
and very fast (linear in lattice points and no stencil involved).



  §.§.§ Ripening without flow


We first consider the ripening of a population of two-dimensional droplets without
flow to test our initialization procedure and to reproduce the known growth law
for ⟨ R⟩. The domain ranges between [-1,1]×[-1,1] and the boundary conditions are periodic. The parameters of the simulation are listed in Table <ref>.
The initialization produced 2737 droplets. Figure <ref>
shows the phase field and chemical potential fields in a small part of
the domain between [0,0.3]×[0,0.45] containing only a small fraction of the droplet count. The phase 1 domains follow the expected evolution: the smaller
droplets shrink and vanish while the larger ones grow. The chemical
potentials are initialized at zero (equilibrium). The Gibbs-Thomson
condition creates differences between droplets and their surrounding
matrix and neighbours of different size. As the grains grow and rarify,
the chemical potentials gradually homogeneize. Figure <ref>
presents the evolution of the mean grain radius compared to a power
law with exponent 1/3. With this, we confirm the adequacy of the
model to reproduce the growth kinetics.





1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 11-11 
3cDomain (D2Q9)          3cPhase-field parameters          3cTransport parameters1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 11-11 
Symbol     Value     Dim          Symbol     Value     Dim          Symbol     Value     Dim[-L_x,L_x]     [-1,1]     [L]          W     3δ x=1.464×10^-3     [L]          (c_0^A,,c_0^B,)     (0.3,0.3)     [–][-L_y,L_y]     [-1,1]     [L]          λ     155.95     [–]          (c_1^A,,c_1^B,)     (0.4,0.4)     [–]N_x nodes     4096     [–]          M_φ     1.2     [L]^2/[T]          (M_0^AA,M_0^BB)     (1,0.8)     [L]^2/[T]5-7 6-7 7-7 
N_y nodes     4096     [–]          3cDroplet initialization          (M_1^AA,M_1^BB)     (1,0.8)     [L]^2/[T]5-7 6-7 7-7 
δ x     4.882×10^-4     [L]          (c_g^A,c_g^B)     (0.31,0.31)     [–]                    δ t     2.5×10^-8     [T]          Φ     0.08     [–]                                       S_avg±Δ S     (1.02±0.957)×10^-4     [L]^2                    1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 11-11 


Parameters for the simulation
of 2D Ostwald ripening. The initial volumes of the grains are uniformly
sampled in the range [S_avg-Δ S,S_avg+Δ S]. The values of collision rates corresponding respectively to M_φ, M^AA and M^BB are: τ_g=0.8776, τ_h^A=0.8146 and τ_h^B=0.7517.













  §.§.§ Ripening and sedimentation


Next, we look at a three-dimensional ripening process including the effect of buoyancy.
Table <ref> lists the
parameters. The gravity acts along the x-axis and the bounce-back
LBM algorithm is used on the top and bottow walls, enforcing u=0
and ∇φ=∇c=∇μ=0.
The lateral boundaries are periodic.





1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 11-11 
3cDomain (D3Q19)          3cPhase-field parameters          3cDroplet initialization1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 11-11 
Symbol     Value     Dim          Symbol     Value     Dim          Symbol     Value     Dim[-L_x,L_x]     [-16,16]     [L]          W     4δ x=6.25×10^-2     [L]          (c_g^A,c_g^B)     (0.32,0.32)     [–][-L_y,L_y]     [-4,4]     [L]          λ     155.95     [–]          Φ     0.08     [–][-L_z,L_z]     [-4,4]     [L]          M_φ     1.2     [L]^2/[T]          V_avg±Δ V     0.08±0.072     [L]^35-7 6-7 7-7 9-11 10-11 11-11 
N_x nodes     2048     [–]          3cFlow parameters          3cTransport parameters5-7 6-7 7-7 9-11 10-11 11-11 
N_y nodes     512     [–]          ρ_0=Δρ     1     [M]/[L]^3          (c_0^A,,c_0^B,)     (0.3,0.3)     [–]N_z nodes     512     [–]          ν     1     [L]^2/[T]          (c_1^A,,c_1^B,)     (0.4,0.4)     [–]δ x     1.5625×10^-2     [L]          σ     10^-3     [MLT^-2]/[L]          (M_0^AA,M_0^BB)     (1,0.8)     [L]^2/[T]δ t     4.2×10^-5     [T]          Δρ g_x     1     [M]/[LT]^2          (M_1^AA,M_1^BB)     (1,0.8)     [L]^2/[T]1-3 2-3 3-3 5-7 6-7 7-7 9-11 10-11 11-11 


Parameters for the simulation
of 3D Ostwald ripening with a gravity-driven flow. The initial grains
volumes are uniformly sampled in the range [V_avg-Δ V,V_avg+Δ V]. The values of collision rates corresponding respectively to ν, M_φ, M^AA and M^BB are: τ_f=0.6111, τ_g=0.6333, τ_h^A=0.6111 and τ_h^B=0.5888.



Figure <ref> shows successive snapshot pictures of
the droplet geometry. The initial condition (Figure <ref>)
is composed of 2035 droplets. At the beginning, diffusive Ostwald ripening
takes place, and after some time (Figure <ref>), the smallest droplets 
have vanished and the remaining ones start to sediment. The droplets then 
accelerate (Figure <ref>), and some have already reached the bottom
wall. A droplet can be seen hanging from the top wall due to capillarity.
Many droplets lose their spherical shape due to coalescence events. Much
later (Figure <ref>), the top of the domain contains almost
no droplets anymore due to the combined effect of ripening, sedimentation
and coalescence. A large drop of the dense phase forms. This huge mass
continues to sediment, but the surrounding smaller droplets rise due to a recirculation
of the flow field. At the final time of the simulation (Figure <ref>),
only a few droplets are left at the bottom. They are still subject
to evaporation due to the very large drop that formed in the center.
The flow recirculation is clearly seen at the final time step on Figure
<ref>. At earlier times, small droplets may also
be seen rising around the large domain of phase 1 before they quickly
evaporate.





Figure <ref>(a) tracks the number of droplets in the simulation as a
function of time. Around t=4t_D, the very large drop 
forms in the center and its size becomes comparable to the system size.
Hence, this single drop totally dominates the droplet distribution: the remaining droplets
are few in number, and are very small in comparison (as seen in Figure <ref>).
Because of the small system size and since all small droplets
evaporate in favor of the largest drop, the growth dynamics rapidly stops to obey a power law.
The development of the sedimentation dynamics is illustrated in Figure fig:Time-evolution(b),
which displays the maximum fluid velocity in the domain as a function of time. In the
beginning, this maximum velocity increases monotonically: the average droplet size 
increases with time, which leads to higher sedimentation velocities. After the formation
of the large drop at t=4t_D, the velocity saturates and oscillates. In this regime,
the finite system size limits the sedimantation velocity of the central
drop, and the highest velocities occur in the upward recirculation.





§ CONCLUSION AND OUTLOOK

We have constructed and tested a model for the simulation of late-stage coarsening in
phase-separating liquids, taking into account buoyancy effects which lead to droplet
sedimentation. Two-phase coexistence and the motion of interfaces are described by
a grand-canonical phase-field model, which can be adapted to arbitrary thermodynamic
properties as described by a free-energy model. The coupled model is simulated
numerically by the lattice-Boltzmann method, which is used to integrate in time
all equations of the model. This makes the simulation code monolithic despite its
multiphysics nature. Therefore, it is particularly well adapted for high-performance
simulations on GPU architectures.

We have demonstrated the capabilities of our approach
by several test cases. In Section <ref> we started by verifying that this model reconstructs the correct coupling between interface kinetics and thermodynamics of a three-components system when compared with the equivalent sharp-interface problem. Next, simulations of Ostwald ripening have been carried out in two dimensions in a large system, and we have verified that the mean droplet radius follows the expected power law with time. Finally, one large-scale three-dimensional simulation of simultaneous coarsening and droplet sedimentation has been performed. It has demonstrated that our model can be
used to study the interplay between coarsening and fluid flow created by droplet sedimentation.

With respect to a multi-component Cahn-Hilliard model, our approach has at least
two advantages: as already stated, it can be adapted to simulate any desired substance
for which free-energy data are available, without restrictions on the choice of
the surface tension. Furthermore, the equations of the phase-field model are second
order in space, rather than the fourth order Cahn-Hilliard equation, which makes its
numerical integration more efficient. The price to pay is that the initial stages of
phase separation cannot be described by this model, which makes it necessary to construct
an initial state ad hoc by making hypotheses on the initial droplet size distribution.
However, since the memory of the initial state is rapidly lost during the coarsening of
a disordered droplet distribution, this should not be a severe problem.

We have written down a model for a ternary mixture here, but there is no difficulty
in generalizing the approach to mixtures with more than three components, as long
as the necessary thermodynamic and kinetic data (free energies and mobilities)
are available. An interesting question that could be studied in the future is the
interplay between coarsening and sedimentation in the macroscopic redistribution
of components: since the sedimentation of the droplets leads to an accumulation
of the denser phase at the bottom of the system, and the sedimentation velocity
is linked to the droplet size, there is a non-trivial coupling between the two
phenomena.



§ ACKNOWLEDGEMENTS


This work was granted access to the HPC resources of IDRIS (super-computer Jean-Zay, partition V100) and CCRT (Topaze, partition A100). Alain Cartalade wishes to thank the SIVIT project involving Orano and EDF for the financial support.




10

gunton1983phase
JD Gunton, M San Miguel, and PS Sahni.
The dynamics of first-order phase transitions.
In C. Domb and J. L. Lebowitz, editors, Phase Transitions and Critical Phenomena, Vol. 8, pages 267–466, New York, 1983. Academic Press.

siggia_1979
Eric D. Siggia.
Late stages of spinodal decomposition in binary mixtures.
Phys. Rev. A, 20:595–605, Aug 1979.

Bray94
A.J. Bray.
Theory of phase-ordering kinetics.
Advances in Physics, 43(3):357–459, 1994.

Bray02
A. J. Bray.
Theory of phase-ordering kinetics.
Advances in Physics, 51(2):481–587, 2002.

gin2017radionuclides
Stephane Gin, Patrick Jollivet, Magaly Tribet, Sylvain Peuget, and Sophie
  Schuller.
Radionuclides containment in nuclear glasses: an overview.
Radiochimica Acta, 105(11):927–959, 2017.

Schuller_etal_JACS2011
Sophie Schuller, Olivier Pinet, and Bruno Penelon.
Liquid–liquid phase separation process in borosilicate liquids
  enriched in molybdenum and phosphorus oxides.
Journal of the American Ceramic Society, 94(2):447–454, 2011.

Pinet_etal_JNM2019
O. Pinet, J.-F. Hollebecque, I. Hugon, V. Debono, L. Campayo, C. Vallat, and
  V. Lemaitre.
Glass ceramic for the vitrification of high level waste with a high
  molybdenum content.
Journal of Nuclear Materials, 519:121–127, 2019.

Cahn58
J. W. Cahn and J. E. Hilliard.
Free energy of a non-uniform system. 1. interfacial free energy.
J. Chem. Phys., 28:258–267, 1958.

brackbill_kothe_zemach_1992
J.U Brackbill, D.B Kothe, and C Zemach.
A continuum method for modeling surface tension.
Journal of Computational Physics, 100(2):335–354, 1992.

jacqmin_1999
David Jacqmin.
Calculation of two-phase navier–stokes flows using phase-field
  modeling.
Journal of Computational Physics, 155(1):96–127, 1999.

henry_tegze_2018
Hervé Henry and György Tegze.
Self-similarity and coarsening rate of a convecting bicontinuous
  phase separating mixture: Effect of the viscosity contrast.
Phys. Rev. Fluids, 3:074306, Jul 2018.

henry_tegze_2019
Hervé Henry and György Tegze.
Kinetics of coarsening have dramatic effects on the microstructure:
  Self-similarity breakdown induced by viscosity contrast.
Phys. Rev. E, 100:013116, Jul 2019.

semprebon_ciro_kruger_2016
Ciro Semprebon, Timm Krüger, and Halim Kusumaatmaja.
Ternary free-energy lattice boltzmann model with tunable surface
  tensions and contact angles.
Phys. Rev. E, 93:033305, Mar 2016.

rasolofomanana_et_al_2022
M.A. Rasolofomanana, C. Cardon, M. Plapp, T. Philippe, H. Henry, and R. Le
  Tellier.
Diffuse-interface modelling of multicomponent diffusion and phase
  separation in the u-o-zr ternary system.
Computational Materials Science, 214:111650, 2022.

ProvatasElder
N. Provatas and K. Elder.
Phase-field methods in materials science and engineering.
Wiley-VCH, Weinheim, 2010.

Steinbach09
I. Steinbach.
Phase-field models in materials science.
Model. Simul. Mater. Sci. Eng., 17(7):073001, 2009.

PlappHandbook
M. Plapp.
Phase-field models.
In T. Nishinaga, editor, The Handbook of Crystal Growth, 2nd
  edition, Vol. 1B, pages 631–668, Amsterdam, 2015. Elsevier.

karma_rappel_1998
Alain Karma and Wouter-Jan Rappel.
Quantitative phase-field modeling of dendritic growth in two and
  three dimensions.
Phys. Rev. E, 57:4323–4349, Apr 1998.

almgren_1999
Robert F. Almgren.
Second-order phase field asymptotics for unequal conductivities.
SIAM Journal on Applied Mathematics, 59(6):2086–2107, 1999.

echebarria_et_al_2004
Blas Echebarria, Roger Folch, Alain Karma, and Mathis Plapp.
Quantitative phase-field model of alloy solidification.
Phys. Rev. E, 70:061604, Dec 2004.

badillo_2012
Arnoldo Badillo.
Quantitative phase-field modeling for boiling phenomena.
Phys. Rev. E, 86:041603, Oct 2012.

kim_kim_suzuki_1999
Seong Gyoon Kim, Won Tae Kim, and Toshio Suzuki.
Phase-field model for binary alloys.
Phys. Rev. E, 60:7186–7197, Dec 1999.

plapp_2011_2
Mathis Plapp.
Unified derivation of phase-field models for alloy solidification
  from a grand-potential functional.
Phys. Rev. E, 84:031601, Sep 2011.

Choudhury12
A. Choudhury and B. Nestler.
Grand-potential formulation for multicomponent phase transformations
  combined with thin-interface asymptotics of the double-obstacle potential.
Phys. Rev. E, 85:021602, 2012.

Plapp16
M. Plapp.
Phase-field modelling of solidification microstructures.
Journal of the Indian Institute of Science, 96(3):179–198,
  2016.

the_lattice_boltzmann_method
Timm Krüger, Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt, Goncalo
  Silva, and Erlend Magnus Viggen.
The Lattice Boltzmann Method: Principles and Practice.
Springer International Publishing, Cham, 2017.

calphad
Larry Kaufman and Henry L. Bernstein.
Computer calculation of phase diagrams with special reference to
  refractory metals.
New York : Academic Press, 1970.

cartalade2016
Alain Cartalade, Amina Younsi, and Mathis Plapp.
Lattice boltzmann simulations of 3d crystal growth: Numerical schemes
  for a phase-field model with anti-trapping current.
Computers & Mathematics with Applications, 71(9):1784–1798,
  2016.

verdier_kestener_cartalade_2020
Werner Verdier, Pierre Kestener, and Alain Cartalade.
Performance portability of lattice boltzmann methods for two-phase
  flows with phase change.
Computer Methods in Applied Mechanics and Engineering,
  370:113266, 2020.

LangerBegRohu
J. S. Langer.
An introduction to the kinetics of first-order phase transitions.
In C. Godrèche, editor, Solids far from equilibrium, Edition
  Aléa Saclay, pages 297–363, Cambridge, UK, 1991. Cambridge University
  Press.

bayle_2020
R. Bayle, O. Cueto, S. Blonkowski, T. Philippe, H. Henry, and M. Plapp.
Phase-field modeling of the non-congruent crystallization of a
  ternary Ge–Sb–Te alloy for phase-change memory applications.
Journal of Applied Physics, 128(18):185101, November 2020.

bayle_2020_phd
Raphael Bayle.
Simulation des mécanismes de changement de phase dans des
  mémoires PCM avec la méthode multi-champ de phase.
PhD thesis, 2020.
Thèse de doctorat dirigée par Plapp, Mathis, Institut
  polytechnique de Paris 2020.

folch_et_al_1999
R. Folch, J. Casademunt, A. Hernández-Machado, and L. Ramírez-Piscina.
Phase-field model for hele-shaw flows with arbitrary viscosity
  contrast. ii. numerical study.
Phys. Rev. E, 60:1734–1740, Aug 1999.

sun_beckermann_2007
Y. Sun and C. Beckermann.
Sharp interface tracking using the phase-field equation.
Journal of Computational Physics, 220(2):626–653, jan 2007.

He-Shan-Doolen_PRE-Rapid1998
Xiaoyi He, Xiaowen Shan, and Gary D. Doolen.
Discrete boltzmann equation model for nonideal gases.
Phys. Rev. E, 57:R13–R16, Jan 1998.

He-Luo_Incompressible_JSP1997
X. He and L.-S. Luo.
Lattice boltzmann model for the incompressible navier-stokes
  equation.
Journal of Statistical Physics, 88(3/4):pp. 927–944, 1997.

zu_he_2013
Y. Q. Zu and S. He.
Phase-field-based lattice boltzmann model for incompressible binary
  fluid systems with density and viscosity contrasts.
Phys. Rev. E, 87:043301, Apr 2013.

Fakhari_etal_PRE2017
Abbas Fakhari, Travis Mitchell, Christopher Leonardi, and Diogo Bolster.
Improved locality of the phase-field lattice-boltzmann model for
  immiscible fluids at high density ratios.
Phys. Rev. E, 96:053301, Nov 2017.

fakhari_rahimian_2010
Abbas Fakhari and Mohammad H. Rahimian.
Phase-field modeling by the method of lattice boltzmann equations.
Phys. Rev. E, 81:036707, Mar 2010.

LBMsaclaycode
 code.
, 2018.

kokkos
H. Carter Edwards, Christian R. Trott, and Daniel Sunderland.
Kokkos: Enabling manycore performance portability through polymorphic
  memory access patterns.
Journal of Parallel and Distributed Computing, 74(12):3202 –
  3216, 2014.
Domain-Specific Languages and High-Level Frameworks for
  High-Performance Computing.

heulens_blanpain_moelans_2011
J. Heulens, B. Blanpain, and N. Moelans.
Phase-field analysis of a ternary two-phase diffusion couple with
  multiple analytical solutions.
Acta Materialia, 59(10):3946 – 3954, 2011.

lahiri_abinandanan_choudhury_2017
Arka Lahiri, T. A. Abinandanan, and Abhik Choudhury.
Theoretical and numerical study of growth in multi-component alloys.
Metallurgical and Materials Transactions A, 48:4463–4476,
  October 2017.

maugis_et_al_1997
P. Maugis, W.D. Hopfe, J.E. Morral, and J.S. Kirkaldy.
Multiple interface velocity solutions for ternary biphase infinite
  diffusion couples.
Acta Materialia, 45(5):1941 – 1954, 1997.

wiemker_2013
Rafael Wiemker.
Total Euler Characteristic as a Noise Measure to aid Transfer
  Function Design.
In Mario Hlawitschka and Tino Weinkauf, editors, EuroVis - Short
  Papers. The Eurographics Association, 2013.




