We develop and analyse the first second-order phase-field model to combine melting and dissolution in multi-component flows. This provides a simple and accurate way to simulate challenging phase-change problems in existing codes. Phase-field models simplify computation by describing separate regions using a smoothed phase field. The phase field eliminates the need for complicated discretizations that track the moving phase boundary. However, standard phase-field models are only first-order accurate. They often incur an error proportional to the thickness of the diffuse interface. We eliminate this dominant error by developing a general framework for asymptotic analysis of diffuse-interface methods in arbitrary geometries. With this framework, we can consistently unify previous second-order phase-field models of melting and dissolution and the volume-penalty method for fluid–solid interaction. We finally validate second-order convergence of our model in two comprehensive benchmark problems using the open-source spectral code Dedalus.

1. Introduction

Many scientific and industrial questions involve fluid flows coupled with phase changes; including sea-ice formation [1], semiconductor crystal manufacture [2], binary alloy solidification [3] and geophysical mantle dynamics [4]. Multi-phase interaction combines the challenges of nonlinear multi-component convection [5] and evolution of phase boundaries [6], creating entirely new effects. Quantifying this complexity demands appropriate mathematical tools.

Moving boundary problems are the standard method to model phase change phenomena. Separate partial differential equations (PDEs) exist in the liquid and solid regions and moving boundary conditions are applied at the interface (figure 1a). A dynamically shifting interface means that the boundary conditions form an essential (often nonlinear) part of the solution [7]. Moving boundaries present many challenges, complicating numerical algorithms [8] and mathematical proofs [9].

Figure 1. (a) A moving boundary formulation in which the domain is partitioned into the fluid Ω⁺, solid Ω⁻ and interface ∂Ω. (b) A phase-field approximation, where the phase ϕ smoothly varies between fluid and solid. The ϕ = 1/2 contour approximates the true interface. We illustrate the asymptotic fluid Ω⁺, solid Ω⁻ and size $O (ε)$ boundary ΔΩ regions of our analysis. (c) The signed-distance coordinate system used to analyse the boundary region ΔΩ. Points on the interface p(s) (blue) are parameterized by arbitrary surface coordinates s. A point off the boundary x can be reached by moving a distance σ in the normal direction $\hat{n} (s)$ from the closest point on the manifold p(s). Coordinate singularities (the corners of the outermost red curves) will occur, but the coordinate system remains well behaved in the interface region (b). The figure is in two dimensions for clarity, though the analysis of §3 is done in three dimensions. (Online version in colour.)

As a possible remedy, it is useful to recall that boundary conditions are a mathematical abstraction; they result from limiting cases of rapid transitions in material properties. There is a long history of reinterpreting discontinuous boundary conditions as smoothed phenomena. Where Gibbs treated capillarity with infinitesimal surfaces [10], van der Waals understood the importance of smoothness at phase boundaries [11]. Where Stefan treated solid–liquid phase boundaries as discontinuous [12], Cahn and Hilliard modelled phase separation as smoothed [3]. Readopting a physics-based viewpoint of boundary conditions allows new possible techniques for addressing complex multi-phase problems. As well as providing a firmer mathematical and physical foundation, smoothed models also simplify numerical implementations by removing the need to track the infinitesimal boundary.

This paper focuses on phase-field models, one of the foremost examples of this smoothed approach. Phase-field models represent distinct phases using a single smoothed phase field ϕ, illustrated in figure 1b [13,14]. The evolution of the phases is then determined by a single set of equations that apply over the entire domain. Many other methods also model phase changes, such as enthalpy methods [15,16], level set methods [17,18], diffuse-domain approaches [19,20] or some immersed-boundary methods [21]. Yet phase-field models stand out for combining several key benefits:

—	They are physically motivated, introduced by Fix [22] and Langer [23] to model free energy near phase boundaries (following from Hohenburg and Halperin’s model C [24]).
—	They generalize canonical models of phase separation, reducing to Allen–Cahn and Hele-Shaw flow (among others) in various asymptotic limits [25,26].
—	They are easily extensible to more general systems, including two-component alloys [27–30], convection [31,32] or multi-phase flows [20,33].
—	They can be made thermodynamically consistent [34–38].
—	They are mathematically rigorous, with well-posedness and convergence results [39,40].
—	They are simple to simulate as they avoid explicit tracking of the interface [22,41–49].

We emphasize this last point. The simulation of moving boundary problems requires specialized algorithms designed to track and apply boundary conditions at the interface. These algorithms can be difficult or impossible to implement in existing codes. For example, spectral methods are popular for their efficiency, but cannot easily handle non-trivial geometries or topologies. Phase-field models (and other diffuse-interface methods) alleviate these difficulties by removing boundary conditions from the problem formulation. By replacing boundary conditions with simple source terms, they can be implemented in general codes for little effort. More general effects can be modelled by changing source terms, as opposed to developing and integrating new algorithms into the codebase. Phase-field models extend the range of phenomena that existing codes can simulate, and accelerate the development of codes to study new scientific problems.

However, phase-field models possess one important drawback for simulation: they must resolve the diffuse interface. For small-scale simulations, this is feasible. But there is a vast disparity between the microscopic scale of the smoothed interface and the macroscopic scale of interest in most problems. This disparity is what makes discontinuous boundary conditions appropriate models in most circumstances. Throughout this paper, we denote this ratio, also known as the Cahn number [50,51], by ε

ε = Cahn number = \frac{microscopic interface width}{bulk system size} .

1.1

Feasible numerical values of ε are orders of magnitude larger than reality. Straightforward analysis implies a commensurate

O (ε)

error with respect to the limiting boundary conditions.

The only way to perform accurate phase-field simulations with achievable values of ε is to accelerate the convergence of the model itself. This can be done through second-order asymptotic analysis in the limit ε → 0. While the first order is sufficient to determine the limiting behaviour, it is the second order that reveals the dominant error of the model. It is then possible to find optimal prescriptions that cancel the dominant error and boost convergence from $O (ε)$ to $O (ε^{2})$ .

This strategy leads to various ‘quantitative’ (i.e. second order) phase-field models, beginning with a correction for arbitrary interface kinetics in pure materials [52,53], and since extended to unequal diffusivities [36,54], multiple components [29], and the combination thereof [55,56]. An introduction to this asymptotic procedure can be found in [57]. Despite much success, progress is difficult. Second-order asymptotic analysis has not yet ascertained a quantitative phase-field model of multi-component convection.¹

This paper presents the first second-order phase-field model of buoyancy-forced convecting binary mixtures. The model, given in §2, builds on first-order models of multi-component convection [31], second-order models of pure melts [59], the diffuse-domain method for Robin boundary conditions [60], and the smooth volume-penalty method for no-slip boundary conditions [61]. We verify second-order convergence in §3 by developing a straightforward asymptotic procedure suitable for general equations and geometries in three dimensions. This procedure allows us to consistently analyse and unify previous second-order phase-field and diffuse-interface methods. We also implement this procedure in the symbolic computing language Mathematica. For brevity, we assume somewhat simplified thermodynamical properties in our model, such as uniform temperature diffusivity, negligible solute within the solid, and Boussinesq buoyancy. Each assumption could be relaxed and analysed using the framework of §3. We finally validate the improved convergence in two comprehensive benchmark problems implemented in the Dedalus numerical code [62] in §4.

2. Models of melting in binary mixtures

(a) Conventional moving boundary formulation

Melting in binary mixtures, such as ice in sea water, is often modelled as a moving boundary problem. We partition the domain into fluid Ω⁺ and solid Ω⁻ regions, pose separate PDEs on each subdomain, and apply boundary conditions at the evolving interface ∂Ω (as in figure 1a).

In the fluid, the temperature T⁺ and dissolved solute concentration C satisfy advection–diffusion equations, and the fluid velocity u and pressure p satisfy incompressible Navier–Stokes equations with Boussinesq buoyancy forcing $g ρ (T^{+}, C) \hat{z}$ ,

\begin{aligned} \partial_{t} T^{+} + u \cdot \nabla T^{+} - κ \nabla^{2} T^{+} & = 0, \\ \partial_{t} C + u \cdot \nabla C - μ \nabla^{2} C & = 0, \\ \partial_{t} u + u \cdot \nabla u - ν \nabla^{2} u + \nabla p + \frac{g ρ (T^{+}, C)}{ρ_{0}} \hat{z} & = 0, \\ \nabla \cdot u & = 0, in Ω^{+}, \end{aligned}

2.1

where κ, μ and ν are the thermal, solutal and momentum diffusivity (each assumed constant), and g, ρ and ρ₀ are the gravitational acceleration, the density of the fluid and a reference density. In the solid only the temperature T⁻ is defined, which follows a diffusive equation,

\partial_{t} T^{-} - κ \nabla^{2} T^{-} = 0, in Ω^{-} .

2.2

We require several boundary conditions at the moving interface [7]. The Gibbs–Thompson relation relates departure of thermosolutal equilibrium (T + mC, where m is the liquidus slope) to a mean-curvature $\bar{K}$ dependent surface energy, and kinetic undercooling proportional to the interfacial normal velocity v. The Stefan condition expresses heat conservation, equating latent heat L release with a discontinuity in temperature flux. The Robin concentration condition ensures total solute conservation. Zero velocity boundary conditions maintain mass conservation,

T + m C - γ \bar{K} + α v = 0, {[κ \hat{n} \cdot \nabla T]}_{-}^{+} + L v = 0, μ \hat{n} \cdot \nabla C^{+} + C^{+} v = 0, u = 0.

2.3

This model involves several idealizations, namely Boussinesq buoyancy forcing, constant diffusivities and densities that are phase, temperature, and concentration invariant, and a linear liquidus relation. We can include more general thermodynamic properties into our framework, but we continue with the current model as it captures many aspects of melting and dissolution in multi-component flows. Kinetic undercooling is only relevant for rapidly solidifying supercooled liquids so we set α = 0. Note that the neglect of density change during melting means the fluid velocity is equal to zero at the moving interface.

(b) Phase-field model

Phase-field models are an alternative approach that is physically motivated and simple to simulate. They represent distinct phases with a smoothed phase field ϕ. This field obeys an Allen–Cahn type equation which forces the phase to ϕ ≈ 1 in the solid and ϕ ≈ 0 in the fluid [31]. The interface is represented implicitly by the level set ϕ = 1/2. The new equations are

\begin{aligned} \partial_{t} T + u \cdot \nabla T - κ \nabla^{2} T & = L \partial_{t} ϕ, \\ \partial_{t} C + u \cdot \nabla C - μ \nabla^{2} C & = \frac{C \partial_{t} ϕ - μ \nabla ϕ \cdot \nabla C}{1 - ϕ + δ}, \\ \partial_{t} u + u \cdot \nabla u - ν \nabla^{2} u + \nabla p + \frac{g ρ (T, C)}{ρ_{0}} \hat{z} & = - \frac{ν}{{(β ε)}^{2}} ϕ u, \\ \nabla \cdot u & = 0, \\ ε \frac{5}{6} \frac{L}{κ} \partial_{t} ϕ - γ \nabla^{2} ϕ & = - \frac{1}{ε^{2}} ϕ (1 - ϕ) (γ (1 - 2 ϕ) + ε (T + m C)) . \end{aligned}

2.4

The new source terms of equation (2.4) can be understood heuristically. At leading order, the phase-field equation develops a tanh-like profile around the interface with thickness ε. Beyond this distance, the phase ϕ tends to its limiting values of zero in the fluid and one in the solid. In the fluid, equation (2.4) reproduces equation (2.1). In the solid, the advective and diffusive solute flux tend to zero (with δ ≪ 1 regularizing the concentration equation for numerical stability), and the velocity is damped by Darcy drag terms for porous media. At next order, thermosolutal forcing perturbs the interface to generate latent heat. For a clear derivation of a similar first-order model, see [31].

Our goal is to surpass this approximate understanding and demonstrate improved convergence of this optimized phase-field model to the original moving boundary formulation. We achieve this with an asymptotic analysis of equation (2.4) as the interface length-scale ε tends to zero. Some approaches focus on a variational derivation of phase-field equations [34–38]; however, we follow the view of [31], which states ‘phase-field equations are only quantitatively meaningful in the sharp-interface limit where they can be ultimately related to experiment’. While variational derivations are useful, it is asymptotic analysis, tested with numerical experiments, that best determines the accuracy of phase-field models. Proving these equations converge to the moving boundary formulation at $O (ε^{2})$ in general geometries is nontrivial, but builds on second-order models of each individual boundary condition; a phase-field model which optimizes the mobility term for zero interface kinetics [59], a concentration equation similar to [31] and the diffuse domain method for Robin boundary conditions [60], and the smooth volume penalty method (which gives β = 1.51044385) [61].

3. Analysis of the phase-field model

In order to understand the phase-field model and demonstrate second-order accuracy, we use a multiple-scales matched-asymptotics framework, which we break into several modular steps:

(i)	Partition the solid Ω⁻, fluid Ω⁺ and size $O (ε)$ boundary ΔΩ regions (figure 1b).
(ii)	Adopt signed distance coordinates in the boundary region (§3a(i), figure 1c).
(iii)	Rescale normal coordinate and operators by ε near the interface (§3a(ii), figure 1b).
(iv)	Expand the variables in an asymptotic power series in ε in each region (§3a(iii)).
(v)	Connect regions with asymptotic matching conditions in the limit ε → 0 (§3a(iv)).
(vi)	Iterate to solve the zeroth-, first- and second-order problems (§3b–d).

This procedure follows our previous analysis of the volume-penalty method [61]. The philosophy of our approach is to determine the evolution of the phase-field model up to and including $O (ε^{2})$ in each region. We show that the variables in the fluid and solid regions, as well as the location of the interface itself, evolve with only $O (ε^{2})$ divergence from the moving boundary formulation. Errors of $O (ε)$ in the temperature and tangential velocity do occur in the boundary region ΔΩ. But this deviation is a necessary consequence of smoothly approximating discontinuous gradients across the interface, and is localized to the boundary region. We thereby derive a second-order accurate phase-field model. We now summarize the key components of this procedure.²

(a) Summary of asymptotic procedure

(i) (Signed-distance coordinate system)

We build a simple orthogonal coordinate system in the boundary region ΔΩ using the signed-distance function from the ϕ = 1/2 level set. The signed distance σ is the minimum distance of a point x to the interface. It follows that the point x must lie in the direction of the unit normal vector $\hat{n}$ from the nearest point on the interface p, which we label with surface coordinates s,

x = p (s) + σ \hat{n} (s) .

3.1

The surface coordinates induce a tangent vector basis t_i. Given orthogonal surface coordinates, we also derive the dual vector basis

\nabla s_{i}

and a unique orthonormal tangent basis

{\hat{t}}_{i}

t_{i} = \frac{\partial p}{\partial s_{i}}, {\hat{t}}_{i} = \frac{t_{i}}{| t_{i} |} and \nabla s_{i} = \frac{{\hat{t}}_{i}}{| t_{i} |} .

3.2

We can then write the surface area measure dA = |t₁||t₂| ds₁ ds₂ and surface gradient

\nabla_{⊥}

\nabla_{⊥} = \nabla s_{1} \frac{\partial}{\partial s_{1}} + \nabla s_{2} \frac{\partial}{\partial s_{2}} = {\hat{t}}_{1} \nabla_{1} + {\hat{t}}_{2} \nabla_{2} .

It is not difficult to show the normal

\hat{n}

is everywhere equal to the gradient of the signed distance. The gradient of the normal is therefore symmetric and diagonalizable. The eigenvectors are the principal directions of curvature (which must align with orthogonal surface coordinates), and the eigenvalues are the principal curvatures κ_i,

\hat{n} = \nabla σ, \nabla_{⊥} \hat{n} = - κ_{1} {\hat{t}}_{1} {\hat{t}}_{1} - κ_{2} {\hat{t}}_{2} {\hat{t}}_{2} = - K .

3.3

We use this orthonormal frame to describe all geometric quantities near the interface. We can then express the gradient using the surface and normal derivatives and the scale tensor J

\nabla = \hat{n} \partial_{σ} + J^{- 1} \cdot \nabla_{⊥}, where J = I - σ K .

3.4

It is straightforward to derive the remaining vector calculus operators from the gradient, which we list in appendix A. Finally, the Cartesian partial time derivative ∂_t can be rewritten in a moving coordinate system using the signed-distance partial time derivative

\partial_{τ}

\partial_{t} = \partial_{τ} - v \partial_{σ} + σ \nabla_{⊥} v \cdot J^{- 1} \cdot \nabla_{⊥} .

3.5

(ii) (Rescaling interfacial coordinates)

We analyse the size ε interfacial region using the rescaled coordinate ξ and derivative $\partial_{ξ}$ ,

σ = ε ξ, \partial_{σ} = \frac{1}{ε} \partial_{ξ} .

3.6

This rescales vector calculus operators (listed in appendix A) through scale factors of the gradient

J = I - ε ξ K, J^{- 1} = \sum_{k = 0}^{\infty} ε^{k} ξ^{k} K^{k} .

3.7

(iii) (Variable expansions with formal power series)

After splitting the domain into the fluid Ω⁺, solid Ω⁻ and interfacial ΔΩ regions, each variable f in each region (fluid f⁺, solid f⁻ and interfacial f) is expressed as a power series in ε,

f^{+} (x, t) = \sum_{k = 0}^{\infty} ε^{k} f_{k}^{+} (x, t), f^{-} (x, t) = \sum_{k = 0}^{\infty} ε^{k} f_{k}^{-} (x, t) and f (ξ, s, τ) = \sum_{k = 0}^{\infty} ε^{k} f_{k} (ξ, s, τ) .

3.8

We do so for the temperature T, solute concentration C, fluid velocity

u = u_{σ} \hat{n} + u_{⊥}

, pressure p, phase field ϕ and interface velocity

v \hat{n}

(which does not depend on ξ). We substitute these series into the hierarchy of equations generated by §3a(ii) to derive a system of equations at each order of ε. Solving each order requires a matching procedure between adjacent regions.

(iv) (Asymptotic matching)

To ensure agreement between different regions, we specify asymptotic matching boundary conditions. This subtle notion requires asymptotic agreement in intermediate zones ξ ∼ ε^−1/2 in the limit that ε → 0. We can then let ξ approach infinity for the inner problem without encountering coordinate singularities (provided $ε ≪ min_{i} | κ_{i}^{- 1} |$ ), and let σ approach zero for the outer problem without entering the interfacial region. That is, for any variable f, we require

lim_{ε \to 0} f (\pm ε^{- 1 / 2} ξ, s, t) \sim lim_{ε \to 0} f^{\pm} (p (s, t) \pm ε^{+ 1 / 2} ξ \hat{n} (s, t), t) .

3.9

Each variable is already expressed as an asymptotic series. Each series term in the outer variables can be further expanded as a Taylor series about the interface at ϕ = 1/2 using the signed distance σ. The matching conditions at each order of ε then simplify to

lim_{ξ \to \pm \infty} f (ξ) = lim_{ξ \to \pm \infty} \sum_{k = 0} ε^{k} f_{k} \sim lim_{ξ \to \pm \infty} \sum_{k = 0} ε^{k} (\sum_{ℓ = 0}^{k} \frac{ξ^{ℓ}}{ℓ!} \partial_{σ}^{ℓ} f_{k - ℓ}^{\pm} |_{σ = 0}) .

3.10

Equipped with our multiple-scales matched-asymptotics procedure, we can now verify second-order convergence of our phase-field model in general smooth geometries.

(b) Zeroth order

At leading order, the phase-field equation reduces to the condition $ϕ_{0}^{\pm} (1 - ϕ_{0}^{\pm}) (1 - 2 ϕ_{0}^{\pm}) = 0$ . We define the solid by $ϕ_{0}^{-} = 1$ , the liquid by $ϕ_{0}^{+} = 0$ , with the boundary region separating them.

Fluid problem—Noting that the phase is zero in the fluid, we recover the desired sharp interface equations for the remaining leading order variables,

\begin{aligned} \partial_{t} T_{0}^{+} + u_{0}^{+} \cdot \nabla T_{0}^{+} - κ \nabla^{2} T_{0}^{+} & = 0, \\ \partial_{t} C_{0}^{+} + u_{0}^{+} \cdot \nabla C_{0}^{+} - μ \nabla^{2} C_{0}^{+} & = 0, \\ \partial_{t} u_{0}^{+} + u_{0}^{+} \cdot \nabla u_{0}^{+} - ν \nabla^{2} u_{0}^{+} + \nabla p_{0}^{+} - B (T_{0}^{+} - N C_{0}^{+}) \hat{z} & = 0, \\ \nabla \cdot u_{0}^{+} & = 0. \end{aligned}

For brevity, we approximate the buoyancy relationship as being linear in temperature and concentration, with proportionality constants B and −NB respectively. This simplification does not affect the convergence results. The solutions may also require external boundary conditions to complete the system.

Solid problem— Within the solid, we reproduce the diffusion equation for the temperature, and zero velocity in the solid. The divergence of the momentum equation reveals a Poisson equation for the pressure. The concentration equation depends sensitively on the decay of the phase to zero, but in the region 1 − ϕ ≪ δ the concentration forcing terms vanish, giving

\partial_{t} T_{0}^{-} - κ \nabla^{2} T_{0}^{-} = 0, \partial_{t} C_{0}^{-} - μ \nabla^{2} C_{0}^{-} = 0, u_{0}^{-} = 0, \nabla^{2} p_{0}^{-} = B \partial_{z} (T_{0}^{-} - N C_{0}^{-}) .

Similarly, internal ice boundary conditions may be necessary.

Boundary problem—The boundary region is defined relative to the interface ϕ = 1/2. To hold true for all ε when expanded into its power series, this implies that

ϕ_{0} (ξ = 0) = \frac{1}{2} and ϕ_{k} (ξ = 0) = 0 for all k > 0.

In this region, the phase-field equation balances the diffusion and reaction terms. The limiting boundary conditions ϕ₀(ξ → +∞) = 0 and ϕ₀(ξ → −∞) = 1 imply a tanh profile for ϕ₀,

\partial_{ξ}^{2} ϕ_{0} = ϕ_{0} (1 - ϕ_{0}) (1 - 2 ϕ_{0}), lim_{ξ \to \pm \infty} ϕ_{0} (ξ) \sim ϕ_{0}^{\pm} |_{σ = 0}, ⟹ ϕ_{0} (ξ) = \frac{1}{2} (1 - \tanh \frac{ξ}{2}) .

Matching the inner heat equation to the outer temperature requires T₀ to be constant in ξ,

\partial_{ξ}^{2} T_{0} = 0, lim_{ξ \to \pm \infty} T_{0} (ξ) \sim T_{0}^{\pm} |_{σ = 0}, ⟹ T_{0} = T_{0}^{\pm} |_{σ = 0} .

The identity

\partial_{ξ} ϕ_{0} = - ϕ_{0} (1 - ϕ_{0})

simplifies the concentration equation. The limiting behaviour of the operator implies exponential growth of the kernel into the solid. This is unphysical, implying that the operand

\partial_{ξ} C_{0}

is zero. Matching to the outer variables requires a constant concentration,

(1 - ϕ_{0}) (\partial_{ξ} + ϕ_{0}) \partial_{ξ} C_{0} = 0, lim_{ξ \to \pm \infty} C_{0} (ξ) \sim C_{0}^{\pm} |_{σ = 0} ⟹ C_{0} = C_{0}^{\pm} |_{σ = 0} .

The divergence constraint and matching to the solid implies zero normal velocity,

\partial_{ξ} u_{σ 0} = 0, lim_{ξ \to \pm \infty} u_{σ 0} (ξ) \sim u_{σ 0}^{\pm} |_{σ = 0}, ⟹ u_{σ 0} = u_{σ 0}^{\pm} |_{σ = 0} = 0.

The tangential momentum equation balances diffusion with damping. The kernel is spanned by an unphysical solution that grows exponentially into the solid, and a physical solution

U (ξ)

that decays exponentially into the solid and is affine in the fluid. Calibrating β [61] allows linear (rather than affine) behaviour of

U (ξ)

in the fluid. Matching to the fluid implies u_⊥0 is zero,

(\partial_{ξ}^{2} - \frac{1}{β^{2}} ϕ_{0}) u_{⊥ 0} = 0, lim_{ξ \to \pm \infty} u_{⊥ 0} (ξ) \sim u_{⊥ 0}^{\pm} |_{σ = 0}, ⟹ u_{⊥ 0} = u_{⊥ 0}^{\pm} |_{σ = 0} = 0.

In summary, the zeroth-order asymptotics shows that the phase field has a tanh profile near the boundary, and that the leading order outer velocities satisfy no-slip boundary conditions. The next order problem reproduces the remaining boundary conditions.

(c) First order

The first-order perturbation of the phase-field equation away from the interface takes the form

(1 - 6 ϕ_{0}^{\pm} + {6 ϕ_{0}^{\pm}}^{2}) ϕ_{1}^{\pm} = (T_{0}^{\pm} + m C_{0}^{\pm}) ϕ_{0}^{\pm} (1 - ϕ_{0}^{\pm}) .

At subsequent orders, the leading order operator

(1 - 6 ϕ_{0}^{\pm} + {6 ϕ_{0}^{\pm}}^{2})

pairs with the highest order

ϕ_{k}^{\pm}

. The forcing terms contain lower-order factors. As the first inhomogeneity is zero, all remaining forcing terms, and therefore all outer phase-field expansions

ϕ_{k \geq 1}^{\pm}

, are zero.

Fluid problem—The remaining fluid equations are linear and homogeneous,

\begin{aligned} \partial_{t} T_{1}^{+} + u_{0}^{+} \cdot \nabla T_{1}^{+} + u_{1}^{+} \cdot \nabla T_{0}^{+} - κ \nabla^{2} T_{1}^{+} & = 0, \\ \partial_{t} C_{1}^{+} + u_{0}^{+} \cdot \nabla C_{1}^{+} + u_{1}^{+} \cdot \nabla C_{0}^{+} - μ \nabla^{2} C_{1}^{+} & = 0, \\ \partial_{t} u_{1}^{+} + u_{0}^{+} \cdot \nabla u_{1}^{+} + u_{1}^{+} \cdot \nabla u_{0}^{+} - ν \nabla^{2} u_{1}^{+} + \nabla p_{1}^{+} - B (T_{1}^{+} - N C_{1}^{+}) \hat{z} & = 0, \\ \nabla \cdot u_{1}^{+} & = 0. \end{aligned}

The external boundary conditions are satisfied by the zeroth-order outer solutions. The first-order corrections therefore have homogeneous external boundaries. They are only perturbed at the melting interface. We show this interfacial perturbation is zero at first order.

Solid problem—The solid equations are the same at first order,

\partial_{t} T_{1}^{-} - κ \nabla^{2} T_{1}^{-} = 0, \partial_{t} C_{1}^{-} - μ \nabla^{2} C_{1}^{-} = 0, u_{1}^{-} = 0, \nabla^{2} p_{1}^{-} = B \partial_{z} (T_{1}^{-} - N C_{1}^{-}) .

Boundary problem—The differential operator of the first-order phase-field equation has a two-dimensional kernel spanned by $\partial_{ξ} ϕ_{0}$ and 6ϕ₀(1 − ϕ₀)(ξ + sinhξ) + sinhξ. The inhomogeneity must be orthogonal to the kernel, implying the Gibbs–Thomson condition,

γ (\partial_{ξ}^{2} - (1 - 6 ϕ_{0} + 6 ϕ_{0}^{2})) ϕ_{1} = - \partial_{ξ} ϕ_{0} (T_{0} + m C_{0} - γ \bar{K}) ⟹ (T_{0}^{\pm} + m C_{0}^{\pm}) |_{σ = 0} - γ \bar{K} = 0.

Matching implies ϕ₁ = 0. Integrating the temperature equation recovers energy conservation,

κ \partial_{ξ}^{2} T_{1} = v_{0} L \partial_{ξ} ϕ_{0}, ⟹ T_{1} = \frac{- v_{0} L}{κ} \int_{ξ}^{\infty} ϕ_{0} d η + \partial_{σ} T_{0}^{+} (0) ξ + T_{1}^{+} (0), {[κ \partial_{σ} T_{0}^{\pm}]}_{σ = 0} + L v_{0} = 0.

Solving the homogeneous concentration equation and matching gives the solute conservation condition,

[(\partial_{ξ} + ϕ_{0}) (μ \partial_{ξ} C_{1} + v_{0} C_{0}) = 0, ⟹ C_{1} = \frac{- v_{0}}{μ} C_{0}^{+} (0) ξ + C_{1}^{+} (0), μ \partial_{σ} C_{0}^{+} |_{σ = 0} + C^{+} |_{σ = 0} v_{0} = 0.

The divergence condition when matched with the solid again implies zero normal velocity,

\partial_{ξ} u_{σ 1} = 0, ⟹ u_{σ 1} = u_{σ 1}^{\pm} |_{σ = 0} = 0.

The tangential velocity is proportional to the physical solution

U (ξ)

. We choose β [61] to ensure linear behaviour in the fluid,

(\partial_{ξ}^{2} - \frac{ϕ_{0}}{β^{2}}) u_{⊥ 1} = 0 ⟹ u_{⊥ 1} (ξ) = \partial_{σ} u_{⊥ 0}^{+} |_{σ = 0} U (ξ), u_{⊥ 1} (ξ \to \infty) \sim \partial_{σ} u_{⊥ 0}^{+} |_{σ = 0} ξ, u_{⊥ 1}^{+} |_{σ = 0} = 0.

We finally use the normal momentum equation to show the pressure is constant at leading order

\partial_{ξ} p_{0} + ν \frac{ϕ_{0}}{β^{2}} u_{σ 1} = 0 ⟹ p_{0} = p_{0}^{\pm} (0) .

The first-order asymptotics have reproduced the Gibbs–Thomson condition with zero interface kinetics and the solute and energy conservation boundary conditions. Hence the phase-field equations will tend to the sharp interface equations in the limit. Calibrating β ensures the first-order outer fluid velocity satisfies homogeneous boundary conditions

u_{σ 1}^{+} |_{σ = 0} = u_{⊥ 1}^{+} |_{σ = 0} = 0

. We now solve the second-order asymptotics to show that the mobility coefficient ε(5/6)(L/κ) ensures homogeneous boundary conditions for the first-order outer temperature and concentration, and that the first-order interfacial velocity error is zero—allowing second-order convergence.

(d) Second order

Fluid problem—At second order, the fluid equations are sourced by the first-order errors,

\begin{aligned} \partial_{t} T_{2}^{+} + u_{0}^{+} \cdot \nabla T_{2}^{+} + u_{2}^{+} \cdot \nabla T_{0}^{+} - κ \nabla^{2} T_{2}^{+} & = - u_{1}^{+} \cdot \nabla T_{1}^{+}, \\ \partial_{t} C_{2}^{+} + u_{0}^{+} \cdot \nabla C_{2}^{+} + u_{2}^{+} \cdot \nabla C_{0}^{+} - μ \nabla^{2} C_{2}^{+} & = - u_{1}^{+} \cdot \nabla C_{1}^{+}, \\ \partial_{t} u_{2}^{+} + u_{0}^{+} \cdot \nabla u_{2}^{+} + u_{2}^{+} \cdot \nabla u_{0}^{+} - ν \nabla^{2} u_{2}^{+} + \nabla p_{2}^{+} - B (T_{2}^{+} - N C_{2}^{+}) \hat{z} & = - u_{1}^{+} \cdot \nabla u_{1}^{+}, \\ \nabla \cdot u_{2}^{+} & = 0. \end{aligned}

Solid problem—The solid velocity is now non-zero from interior pressure and forces,

\begin{aligned} \partial_{t} T_{2}^{-} = κ \nabla^{2} T_{2}^{-}, \partial_{t} C_{2}^{-} = μ \nabla^{2} C_{2}^{-}, \\ \nabla^{2} p_{2}^{-} = B \partial_{z} (T_{2}^{-} - N C_{2}^{-}), \frac{ν}{β^{2}} u_{2}^{-} = B (T_{0}^{-} - N C_{0}^{-}) \hat{z} - \nabla p_{0}^{-} . \end{aligned}

Boundary problem—The phase-field equation now has a more complex inhomogeneous term:

γ (\partial_{ξ}^{2} - (1 - 6 ϕ_{0} + 6 ϕ^{2})) ϕ_{2} = - \partial_{ξ} ϕ_{0} (\frac{5}{6} \frac{L}{κ} v_{0} + (T_{1} + m C_{1}) - ξ γ \bar{K^{2}}) .

We again apply a solvability condition. The odd terms drop out, and integration shows

\int_{- \infty}^{\infty} - \partial_{ξ} ϕ_{0} ϕ_{0} (1 - ϕ_{0}) d ξ = \frac{1}{6}

and

\int_{- \infty}^{\infty} - \partial_{ξ} ϕ_{0} ϕ_{0} (1 - ϕ_{0}) \int_{ξ}^{\infty} ϕ_{0} d η d ξ = \frac{5}{36}

. The interfacial velocity term thus vanishes, giving a constraint between first-order concentration and temperature:

\int_{- \infty}^{\infty} \partial_{ξ} ϕ_{0}^{2} (T_{1}^{\pm} (0) + m C_{1}^{\pm} (0) + \frac{L v_{0}}{κ} (\frac{5}{6} - \int_{ξ}^{\infty} ϕ_{0} d η)) d ξ = 0 ⟹ (T_{1}^{\pm} + m C_{1}^{\pm}) |_{σ = 0} = 0.

An explicit formula for ϕ₂ (using variation of parameters) is non-trivial but it decays exponentially to zero in either direction. The temperature equation is then integrated. If initially

T_{1}^{\pm} |_{σ = 0} = 0

, then matching prevents linear asymptotic behaviour of T₂, implying v₁ = 0,

\begin{aligned} κ \partial_{ξ}^{2} T_{2} & = L v_{1} \partial_{ξ} ϕ_{0} + \partial_{t} T_{0} - κ \nabla_{⊥} \cdot \nabla_{⊥} T_{0} + (κ \bar{K} - v_{0}) \partial_{ξ} T_{1}, \\ κ T_{2} & = (\partial_{t} - κ \nabla_{⊥} \cdot \nabla_{⊥} + (κ \bar{K} - v_{0}) \partial_{σ}) T_{0}^{+} |_{σ = 0} \frac{ξ^{2}}{2} + (κ \bar{K} - v_{0}) \frac{v_{0} L}{κ} \int_{ξ}^{\infty} \int_{η}^{\infty} ϕ_{0} d ζ d η + T_{2}^{\pm} |_{σ = 0} . \end{aligned}

The concentration equation is similarly integrated. The inhomogeneity is constant in ξ, and explicitly solvable. Using previous solutions and

C_{1}^{\pm} |_{σ = 0} = T_{1}^{\pm} |_{σ = 0} = v_{1} = 0

, we find C₂,

\begin{aligned} (\partial_{ξ} + ϕ_{0}) (μ \partial_{ξ} C_{2} + v_{0} C_{1} + v_{1} C_{0}) & = \partial_{t} C_{0} + μ \bar{K} \partial_{ξ} C_{1} - μ \nabla_{⊥} \cdot \nabla_{⊥} C_{0}, \\ μ C_{2} = \frac{v_{0}^{2}}{μ^{2}} C_{0}^{+} |_{σ = 0} \frac{ξ^{2}}{2} + (\partial_{t} + \bar{K} v_{0} & - μ \nabla_{⊥} \cdot \nabla_{⊥}) C_{0}^{\pm} |_{σ = 0} \int_{0}^{ξ} \frac{- \log ϕ}{1 - ϕ} d η + C_{2}^{\pm} |_{σ = 0} . \end{aligned}

The divergence equation implies the second-order normal velocity is now no longer zero,

\partial_{ξ} u_{σ 2} = - \nabla_{⊥} \cdot u_{⊥ 1}, ⟹ u_{σ 2} = u_{σ 2} (- \infty) - \nabla_{⊥} \cdot (\partial_{σ} u_{0}^{+} (0)) \int_{- \infty}^{ξ} U d ζ .

The second-order tangential velocity can then be solved using variation of parameters

\begin{aligned} ν (\partial_{ξ}^{2} - \frac{ϕ_{0}}{β^{2}}) u_{⊥ 2} & = \nabla_{⊥} p_{0} + (ν \bar{K} - v_{0}) \partial_{ξ} u_{⊥ 1} + B (T_{0} - N C_{0}) {\hat{z}}_{⊥} = R, \\ u_{⊥ 2} & = (Q + c) U, where Q \equiv \int_{0}^{ξ} - \frac{\int_{- \infty}^{η} R U d ζ}{U^{2}} d η, \end{aligned}

where c is a constant of integration that cancels the linear part of the solution into the fluid. The normal momentum equation can then be integrated to determine p₁. Matching requires no constant term in the limiting behaviour of the pressure into the fluid,

\begin{aligned} \partial_{ξ} p_{1} & = ν \nabla_{⊥} \cdot (\partial_{σ} u_{⊥ 0}^{+} |_{σ = 0}) (\frac{ϕ_{0}}{β^{2}} (\int_{- \infty}^{ξ} U d η - U^{'} (ξ))) - \frac{ν ϕ_{0}}{β^{2}} u_{σ 2}^{-} |_{σ = 0} - B (T_{0}^{+} - N C_{0}^{+}) |_{σ = 0} {\hat{z}}_{σ}, \\ p_{1} & = ν \nabla_{⊥} \cdot (\partial_{σ} u_{⊥ 0}^{+} |_{σ = 0}) U (ξ) - B (T_{0}^{+} - N C_{0}^{+}) |_{σ = 0} {\hat{z}}_{σ} ξ \\ + \frac{ν}{β^{2}} (\nabla_{⊥} \cdot (\partial_{σ} u_{⊥ 0}^{+} |_{σ = 0}) (\int_{ξ}^{\infty} ϕ_{0} \int_{- \infty}^{ζ} U (η) d η d ζ) + u_{σ 2}^{-} |_{σ = 0} \int_{ξ}^{\infty} ϕ_{0} d η) . \end{aligned}

The second-order asymptotic analysis shows that calibrating the mobility and damping parameters ensures homogeneous boundary conditions of the first-order outer solutions at the interface. Combined with homogeneous external boundary conditions and homogeneous linear evolution equations at first order, this implies that if the outer solutions are initialized correct to $O (ε^{2})$ , then the fields and interfacial velocity will evolve accurate to $O (ε^{2})$ over time. In reality, the chaotic nature of many fluid dynamics problems prevents convergence beyond the Lyapunov timescale of the flow, but at each point in time the system behaves correctly to within second-order accuracy.

4. Numerical validation of the model

We now validate the asymptotic arguments of §3 in two benchmark problems. In each problem, we calculate a numerical reference solution corresponding to the ‘sharp interface’ equations of §2a. We then show the optimal phase-field equations of §2b achieve convergence of $O (ε^{2})$ to these reference solutions. Both the reference and phase-field problems are simulated using the flexible and efficient spectral code Dedalus [62].³ Dedalus is a general-purpose code that takes PDEs in plain text and automatically discretizes them using spectral series (including Fourier series, Chebyshev polynomials and other orthogonal polynomials). Linear terms in the equations are converted into sparse banded matrices and nonlinear terms are treated pseudospectrally. Included implicit–explicit timesteppers integrate the linear terms implicitly, and the nonlinear terms explicitly. Dedalus is written in the Python programming language to aid rapid development, and uses compiled libraries (including BLAS, LAPACK and MPI) to perform efficient parallelized simulations at scale.

(a) Melting and dissolution at a stagnation point

The first benchmark problem examines warm liquid with dissolved solute flowing towards a melting interface at a stagnation point (similar to §5 of [63] and 5.6 of [61]). We compare the sharp interface and phase-field approximations to this problem, as illustrated in figure 2. The symmetries of the system allow us to significantly simplify equation (2.1), revealing a steady travelling wave similarity solution for the moving boundary and phase-field formulations. We do this by transforming to a frame moving leftward at the steady interface melting speed −v. We solve the system as a nonlinear boundary value problem in Dedalus.

Figure 2. On the left, we illustrate stagnation point flow towards a melting interface. The fluid velocity decreases to zero at the interface (black arrows x > 0), and the vertical velocity is proportional to y. The interface velocity (grey arrow at x = 0) uniformly recedes to the left. The temperature (in colour) decreases towards the left, and is invariant in y (as are all quantities excepting the vertical fluid velocity). The concentration (not shown) decreases similarly. We also plot reference and phase-field solutions for temperature T, solute concentration C, normal velocity u and phase field ϕ as a function of x in the second and third columns. Only the temperature exists inside the solid in the reference solution. We negate the normal fluid velocity u for clarity. (Online version in colour.)

Sharp interface model—We solve a diffusion equation for the solid temperature T⁻, advection–diffusion equations for the liquid temperature T⁺ and solute concentration C and a nonlinear third-order equation for the horizontal fluid velocity u,

\begin{aligned} T^{-} (- 1) = - 1, & T^{+} (0) = T^{-} (0), {[\partial_{x} T (0)]}_{-}^{+} = - \frac{L v}{κ}, u (0) = 0, & \partial_{x} u (0) = 0, \\ T^{+} (1) = 1, & T^{\pm} (0) = - m C (0), \partial_{x} C (0) = - \frac{C (0) v}{μ}, C (1) = 0, & \partial_{x} u (1) = - 1. \end{aligned}

4.2

Phase-field model—The phase field instead implicitly models the interfacial boundary conditions through various equation terms, which reduce equation (2.4) to

\begin{aligned} ν \partial_{x}^{2} T & = ((1 - ϕ) u - v) \partial_{x} T + L v \partial_{x} ϕ, \\ μ \partial_{x}^{2} C & = (u - v) \partial_{x} C - \partial_{x} \log (1 - ϕ + δ) (μ \partial_{x} C + v C), \\ ν \partial_{x}^{3} u & = 1 + (u - v) \partial_{x}^{2} u - {(\partial_{x} u)}^{2} + \frac{ν}{{(β ε)}^{2}} ϕ \partial_{x} u, \\ γ \partial_{x}^{2} ϕ & = - ε \frac{5}{6} \frac{L}{κ} v \partial_{x} ϕ + \frac{γ}{ε^{2}} ϕ (1 - ϕ) (1 - 2 ϕ) + \frac{1}{ε} ϕ (1 - ϕ) (T + m C) . \end{aligned}

4.3

We now only require the Dirichlet outer boundary conditions of before to solve the problem

T (- 1) = - D, C (- 1) = 0, ϕ (- 1) = 1, u (- 1) = 0, \partial_{x} u (- 1) = 0

4.4

and

T (1) = 1, C (1) = 1, ϕ (1) = 0, ϕ (0) = \frac{1}{2}, \partial_{x} u (1) = - 1.

4.5

Requiring ϕ(0) = 1/2 is necessary to fix the problem in the moving frame to determine v.

Results—We solve each problem as a nonlinear boundary value problem in Dedalus. We discretize each variable on the solid (−1 < x < 0) and fluid (0 < x < 1) domains using Chebyshev polynomials. Newton–Kantorovich iteration then converges on a solution with a tolerance of 10⁻¹² (see [64, appendix C]). We reproduce example phase field and reference snapshots in figure 2 for ε = 0.05, κ = μ = ν = 1/10 and D = m = L = γ = 1, using 64 grid points for each subdomain. We set δ = 2 × 10⁻⁵ to regularize the solute equation within the solid. The disagreement (though small) is visible by eye for the temperature, concentration and normal velocity.

We then perform a quantitative analysis of convergence in figure 3. We test seven logarithmically spaced values of ε from 10⁻¹ to 10⁻³, for the previous control parameters, and using 256 grid points in each subdomain. We quantify the difference between the reference and phase-field solutions with the L¹ and L^∞ error norms of each variable (u, T⁺, T⁻, C, v) as a function of ε. We find clear $O (ε^{2})$ convergence in the L¹ error norm of each field variable, as well as for the L^∞ error norm of the u, C and v variables. The L¹ convergence demonstrates the quantitative accuracy of the phase-field model. The reason for the apparently restricted $O (ε)$ L^∞ convergence of the temperature fields is the jump in temperature gradient of the reference solution at the interface. The smooth phase-field model cannot follow this kink leading to an $O (ε)$ disagreement that is localized to the boundary. Tangential velocities also suffer from this reduced continuity in general.

(b) Double diffusive melting

In our second benchmark, we examine buoyancy driven flow of warm liquid with dissolved solute underneath a melting solid layer. This problem develops non-trivial geometries from the evolution of the flow. To simulate the reference formulation in Dedalus, we use an evolving coordinate system that maps the fluid and solid regions to a stationary rectangular domain. This transformation allows an efficient spectral discretization using Dedalus, and is used in similar spectral solvers [65]. We repeat this remapping for the phase-field formulation as it concentrates resolution near the ϕ = 1/2 level set. This allows us to efficiently and accurately simulate much smaller ε. By comparing the phase-field simulation to the reference problem as we decrease ε, we demonstrate second-order accuracy of the model.

Sharp differential equations—The full domain exists between 0 < z < 2 and 0 < x < 4. It is partitioned by the interface at height z = h(t, x). Above the interface, we solve the heat equation for the solid temperature field T⁻ with diffusivity κ. Below the interface, we solve incompressible Boussinesq hydrodynamics (using a first-order formulation in terms of velocity u, vorticity q and augmented pressure $p = P + \frac{1}{2} | u |^{2}$ ), with advection and diffusion of temperature T⁺ and solute concentration C,

\begin{aligned} \partial_{t} T^{-} - κ \nabla^{2} T^{-} = 0, \nabla \cdot u = 0, \\ \partial_{t} T^{+} - κ \nabla^{2} T^{+} = - u \cdot \nabla T^{+}, q - \nabla \times u = 0, \\ \partial_{t} C - μ \nabla^{2} C = - u \cdot \nabla C, \partial_{t} u + \nabla p - ν \nabla \times q = q \times u + (T^{+} - N C) {\hat{e}}_{z}, \end{aligned}

4.6

where ν, κ and μ are the momentum, temperature and solute diffusivity. This first-order reformulation is required in the Dedalus solver, and has the added benefit of simplifying the mathematical details of the coordinate transformation. We set the buoyancy parameter B = 1.

Boundary conditions—At the top, we specify conservative temperature boundary conditions. At the bottom, we specify no-slip boundaries with no-flux temperature and solute conditions. The zeroth mode for the vertical velocity is replaced by a choice of pressure gauge,

\partial_{z} T^{-} (t, x, 2) = 0, \partial_{z} T^{+} (t, x, 0) = 0, \partial_{z} C (t, x, 0) = 0, u (t, x, 0) = 0, p (t, x, 0) = 0.

At the interface z = h(t, x), we have a Gibbs–Thomson boundary condition, matching temperature boundary conditions, energy and solute conservation boundary conditions and zero velocity boundary conditions required for mass conservation,

\begin{aligned} T^{+} + m C = - γ \nabla \cdot \hat{n}, & \hat{n} \cdot \nabla T^{+} - \hat{n} \cdot \nabla T^{-} + \frac{L}{κ} \hat{n} \cdot \partial_{t} h {\hat{e}}_{z} = 0, & u_{x} = 0, \\ T^{+} - T^{-} = 0, & \hat{n} \cdot \nabla C + \frac{C}{μ} \hat{n} \cdot \partial_{t} h {\hat{e}}_{z} = 0, & u_{z} = 0. \end{aligned}

4.7

Initial conditions—We initialize the problem with h(0, x) = 1 and zero velocity and pressure, and a decreasing concentration profile with height. We apply a large perturbation to the linearly decreasing temperature field to initiate convection,

[T (0, x, z) = 1 - z + \exp (5^{2} ({(x - 2)}^{2} + {(z - 0.5)}^{2}))

and

[C (0, x, z) = 0.05 + (1 - 0.05) \frac{1}{2} (1 - \tanh (10 (z - 0.5))) .

Phase-field equations—The phase-field equations are a similar reformulation of equation (2.4):

\begin{aligned} \partial_{t} T - κ \nabla^{2} T - L \partial_{t} ϕ = - u \cdot \nabla T, & \partial_{t} C - μ \nabla^{2} C = - u \cdot \nabla C + \frac{\partial_{t} ϕ C - \nabla C \cdot \nabla ϕ}{1 - ϕ + δ}, \\ \nabla \cdot u = 0, & ε \frac{5}{6} \frac{L}{κ} \partial_{t} ϕ - γ \nabla^{2} ϕ = - \frac{1}{ε^{2}} ϕ (1 - ϕ) (γ (1 - 2 ϕ) + ε (T + m C)), \\ q - \nabla \times u = 0, & \partial_{t} u + \nabla p - ν \nabla \times q = q \times u + (T^{+} - N C) {\hat{e}}_{z} - \frac{ν}{{(β ε)}^{2}} ϕ u . \end{aligned}

4.8

We use the optimal damping prescription β = 1.51044385 [61] and choose δ = 10⁻⁴. We specify the same insulating no-slip boundary conditions at z = 0 and z = 2. The boundary conditions at z = h(t, x) require continuity of each field and its derivatives. We initialize using the same initial conditions as for the remapped simulation, plus the initial conditions for the phase field

ϕ (0, x, z) = \frac{1}{2} (1 + \tanh (\frac{1}{2 ε}) (z - 1)) .

Remapped coordinates—To solve the melting problem, we remap our evolving domain in Cartesian space to a fixed rectangular domain with the new coordinates τ, ξ and ζ^±,

τ (t, x, z) = t, ξ (t, x, z) = x, ζ^{+} (t, x, z) = \frac{z}{h (t, x)}, ζ^{-} (t, x, z) = \frac{z - h}{2 - h (t, x)},

for which we plot equally spaced level set contours in figure 4. The differential geometry of the new coordinates is presented in appendix B. We solve the phase-field equations on the remapped domains of the reference problem to concentrate resolution near the phase-field interface and speed up comparison of quantities between simulations.

Model and numerical parameters—We simulate these equations using model parameter values from table 1. We simulate the reference equations using 64 Chebyshev polynomials in the vertical direction and 128 Fourier modes in the periodic horizontal direction. After determining the reference solution, we discretize the phase-field equations on the same evolving domain of the reference simulation for several values of ε. To interpolate the reference geometries into the phase-field simulations between the saved time and grid points, we use third-order interpolating splines. The phase-field simulations discretize the vertical ζ basis using three compound Chebyshev bases $[0, 10 ε] \cup [10 ε, 1 - 10 ε] \cup [1 - 10 ε, 1]$ , with resolutions of 32, 64 and 32 modes, respectively. This greatly reduces the simulation cost. We use a time step size of Δt = 10⁻³, 5 × 10⁻⁴, 2.5 × 10⁻⁴, 2 × 10⁻⁴, 5 × 10⁻⁵ for decreasing choice of ε, and integrate in time using a second-order multistep semi-implicit backwards difference formula (SBDF2).

Table 1. Model parameters used in the second benchmark problem.
Collapse
ν	κ	μ	γ	L	m	N	ε
10⁻²	10⁻²	10⁻²	10⁻²	1 $^{}$	0.2	0	$\sqrt{5} \times 10^{- 3}, 10^{- 2.5}, \sqrt{2} \times 10^{- 2.5}, \sqrt{5} \times 10^{- 2.5}, 10^{- 2}$

Results—A time series of the temperature and concentration fields of the reference solution is given in figure 5, which illustrates a rising buoyant plume from the initial temperature anomaly in the fluid. The warm solute-laden liquid melts the interface in the middle more rapidly than the ambient liquid at the sides, causing a trough to develop.

In figure 6, we plot several error metrics of the phase-field simulations. In the first two columns, we plot the spatial error normalized by the L¹ error for the liquid temperature T⁺, solute concentration C, Cartesian velocity components u_x and u_z, and true pressure $P = p - \frac{1}{2} | u |^{2}$ . We plot these normalized spatial errors for the smoothest (ε = 10⁻²) and sharpest ( $ε = \sqrt{5} \times 10^{- 3}$ ) phase-field simulations at the final time t = 10. These plots reveal a consistent spatial error profile between simulations. To understand the amplitude of the spatial error, we plot the L¹ and L^∞ error norms of each variable as a function of ε in the third and fourth columns of figure 6. We see clear second-order convergence of all variables in L¹ norm. We note that the structure of the boundary layer of the tangential velocity and temperature affects the L^∞ norm. The phase-field model causes a kink in the temperature and tangential velocity near the interface, which leads to an O(ε) error of these variables in the interfacial region. However, this error is localized to the boundary, and does not propagate outward. (This trend is difficult to notice in the temperature plot as the error within the fluid still dominates the $O (ε)$ boundary error for the moderate choices of ε chosen.) We therefore achieve second-order convergence in the fluid and solid regions due to our optimal calibration of the phase-field model parameters.

5. Conclusion

In this paper, we provide a general framework for analysing the convergence of phase-field models. We use this procedure to develop a second-order phase-field model of melting and dissolution in multi-component flows. This is a concrete advancement that showcases second-order accurate approximations of many common boundary conditions; no-slip Dirichlet boundaries, Neumann boundaries, Robin boundaries and Stefan boundaries. We also verify these prescriptions in two thorough benchmark problems with accurate reference solutions. By developing a framework to validate this model, we now possess the machinery required to create second-order accurate extensions to more general thermodynamic properties. We can also consider yet higher-order analyses of this and other models. An automated approach of Richardson sequence extrapolation could also be considered to generate higher-order accurate models, as was done by the authors for the volume-penalty method [61]. To be clear, phase-field models are not necessarily the most appropriate choice for any problem. Remapping was also shown to be an effective strategy for sufficiently simple geometries in §4b. However, our second-order phase-field model is simple to implement, is more accurate than standard phase-field models, and is applicable in much more challenging geometries than remapping approaches.

Data accessibility

The Dedalus code is free and open-source and available at http://dedalus-project.org/. All simulation code, data and analysis scripts and plots are available at https://github.com/ericwhester/phase-field-code. The Mathematica asymptotics script is also available at that link.

Authors' contributions

E.W.H. project conception and design, writing, mathematics (asymptotics and remapping), simulation code, data analysis and visualization; L.-A.C. project design and conception; B.F. project design and conception, phase-field expertise; K.J.B. domain remapping simulation expertise, Dedalus code; G.M.V. mathematics of signed distance coordinates, Dedalus code. All authors provided critical revision of the manuscript, approve the final version of the paper and agree to be held accountable for the work.

Competing interests

We declare that we have no competing interests.

Funding

E.W.H. acknowledges support from The University of Sydney Phillip Hofflin International Research Travel Scholarship, Postgraduate Research Support Scheme and William and Catherine McIlrath Scholarship. L.-A.C. acknowledges support from the European Union’s Horizon 2020 research and innovation programme under Marie Skłodowska-Curie grant agreement no. MIMOP 793450.

Acknowledgements

We acknowledge PRACE for awarding us access to Marconi HPC at CINECA, Italy, and The University of Sydney for access to the Artemis HPC.

Appendix A. Signed-distance coordinates

The signed distance σ of the point x is the minimum distance to the surface. Labelling surface points p and unit normals $\hat{n}$ with orthogonal surface coordinates s, we have

[x = p (s) + σ \hat{n} (s) .

This gives a tangent vector basis t_i, dual basis

\nabla s_{i}

and orthonormal tangent basis

{\hat{t}}_{i}

[t_{i} = \frac{\partial p}{\partial s_{i}}, {\hat{t}}_{i} = \frac{t_{i}}{| t_{i} |} and \nabla s_{i} = \frac{{\hat{t}}_{i}}{| t_{i} |} .

These give the surface gradient

\nabla_{⊥}

, area measure dA = |t₁||t₂| ds₁ ds₂ and surface divergence,

\nabla_{⊥} = \nabla s_{1} \frac{\partial}{\partial s_{1}} + \nabla s_{2} \frac{\partial}{\partial s_{2}} = {\hat{t}}_{1} \nabla_{1} + {\hat{t}}_{2} \nabla_{2} and \nabla_{⊥} \cdot u_{⊥} = \frac{1}{| t_{1} | | t_{2} |} (\frac{\partial}{\partial s_{1}} (| t_{2} | u_{1}) + \frac{\partial}{\partial s_{2}} (| t_{1} | u_{2})) .

The normal

\hat{n}

is everywhere equal to the gradient of the signed distance, and therefore its gradient is symmetric and diagonalizable. The eigenvectors are the principal directions of curvature (which align with orthogonal surface coordinates), and the eigenvalues are the principal curvatures κ_i,

\hat{n} = \nabla σ, \nabla_{⊥} \hat{n} = - κ_{1} {\hat{t}}_{1} {\hat{t}}_{1} - κ_{2} {\hat{t}}_{2} {\hat{t}}_{2} = - K .

We express the gradient using the orthonormal frame

(\hat{n}, {\hat{t}}_{1}, {\hat{t}}_{2})

and the scale tensor J

\nabla = \hat{n} \partial_{σ} + J^{- 1} \cdot \nabla_{⊥}, where J = I - σ K .

The determinant |J| relates mean curvature

\bar{K}

, Gaussian curvature |K | and volume measure dV,

\bar{K} = κ_{1} + κ_{2}, | K | = κ_{1} κ_{2}, | J | = 1 - σ \bar{K} + σ^{2} | K | and d V = | J | d σ d A .

We define the adjugate tensors

\hat{J}

and

\hat{K}

\hat{K} = | K | K^{- 1} = κ_{2} {\hat{t}}_{1} {\hat{t}}_{1} + κ_{1} {\hat{t}}_{2} {\hat{t}}_{2}, \hat{J} = | J | J^{- 1} = I - σ \hat{K} .

We also define the cross product with the unit normal, which admits simple identities,

\begin{aligned} \nabla^{⊥} = \hat{n} \times \nabla_{⊥}, & \nabla^{⊥} \cdot \nabla_{⊥} = \nabla_{⊥} \cdot \nabla^{⊥} = 0, \hat{n} \times (J u_{⊥}) = \hat{J} u^{⊥}, \\ u^{⊥} = \hat{n} \times u_{⊥}, & u^{⊥} \cdot u_{⊥} = u_{⊥} \cdot u^{⊥} = 0, \hat{n} \times u^{⊥} = - u_{⊥} . \end{aligned}

We next find the gradient of the basis vectors, and define Ricci rotation coefficients,

\begin{aligned} \nabla_{⊥} \hat{n} = - κ_{1} {\hat{t}}_{1} {\hat{t}}_{1} - κ_{2} {\hat{t}}_{2} {\hat{t}}_{2}, & R_{i}^{j k} = {\hat{t}}_{j} \cdot (\nabla {\hat{t}}_{i}) \cdot {\hat{t}}_{k}, R_{1}^{12} = {\hat{t}}_{1} \cdot (\nabla {\hat{t}}_{1}) \cdot {\hat{t}}_{2} = ω_{1}, \\ \nabla_{⊥} {\hat{t}}_{i} = κ_{i} {\hat{t}}_{i} \hat{n} + R_{i}^{j k} {\hat{t}}_{j} {\hat{t}}_{k}, & R_{i}^{j k} = - R_{k}^{j i}, R_{2}^{21} = {\hat{t}}_{2} \cdot (\nabla {\hat{t}}_{2}) \cdot {\hat{t}}_{1} = - ω_{2} . \end{aligned}

These relations allow us to calculate all relevant vector calculus operators,

\begin{aligned} \nabla \cdot u & = \frac{\partial_{σ} (| J | u_{σ})}{| J |} + \frac{\nabla_{⊥} \cdot (\hat{J} u_{⊥})}{| J |}, \\ \nabla^{2} f & = \frac{\partial_{σ} (| J | \partial_{σ} f)}{| J |} + \frac{\nabla_{⊥} \cdot (\hat{J} J^{- 1} \nabla_{⊥} f)}{| J |}, \\ \nabla \times u & = - \hat{n} \frac{\nabla_{⊥} \cdot (\hat{J} u^{⊥})}{| J |} + {\hat{J}}^{- 1} (\partial_{σ} (\hat{J} u^{⊥}) - \nabla_{⊥} u_{σ}) \\ - \nabla \times \nabla u = \frac{\hat{n}}{| J |} (- \nabla_{⊥} \cdot (\hat{J} J^{- 1} (\partial_{σ} (J u_{⊥}) - \nabla_{⊥} u_{σ}))) \\ + {\hat{J}}^{- 1} \partial_{σ} (\hat{J} J^{- 1} (\partial_{σ} (J u_{⊥}) - \nabla_{⊥} u_{σ})) + {\hat{J}}^{- 1} \nabla_{⊥} (\frac{\nabla^{⊥} \cdot (J u_{⊥})}{| J |}), \\ \nabla u & = \hat{n} \hat{n} \partial_{σ} u_{σ} + \hat{n} \partial_{σ} u_{⊥} + J^{- 1} (\nabla_{⊥} u_{σ} + K u_{⊥}) \hat{n} + J^{- 1} (\nabla_{⊥} u_{⊥} - K u_{σ}), \\ u \cdot \nabla f & = u_{σ} \partial_{σ} f + u_{⊥} \cdot J^{- 1} \nabla_{⊥} f, \\ u \cdot \nabla u & = \hat{n} (u_{σ} \partial_{σ} u_{σ} + u_{⊥} \cdot J^{- 1} (\nabla_{⊥} u_{σ} + K u_{⊥})) + u_{σ} \partial_{σ} u_{⊥} + u_{⊥} \cdot J^{- 1} (\nabla_{⊥} u_{⊥} - K u_{σ}), \\ \partial_{t} f & = \partial_{τ} f - v \partial_{σ} f + σ \nabla_{⊥} v \cdot J^{- 1} \cdot \nabla_{⊥} f, \\ \partial_{t} u & = (\partial_{τ} u_{σ} - v \partial_{σ} u_{σ} + σ \nabla_{⊥} v \cdot J^{- 1} \cdot (\nabla_{⊥} u_{σ} + K u_{⊥})) \hat{n} \\ + (\partial_{τ} u_{⊥} - v \partial_{σ} u_{⊥} + σ \nabla_{⊥} v \cdot J^{- 1} \cdot (\nabla_{⊥} u_{⊥} - K u_{σ})) . \end{aligned}

Rescaling—Rescaling the normal coordinate by ξ = σ/ε gives the following operators:

\begin{aligned} \nabla & = ε^{- 1} \hat{n} \partial_{ξ} + \sum_{k = 0}^{\infty} ε^{k} ξ^{k} K^{k} \nabla_{⊥}, \\ | J | \nabla \cdot u & = ε^{- 1} \partial_{ξ} u_{σ} - \partial_{ξ} (ξ \bar{K} u_{σ}) + \nabla_{⊥} \cdot u_{⊥} + ε (\partial_{ξ} (ξ^{2} | K | u_{σ}) - ξ \nabla_{⊥} \cdot (\hat{K} u_{⊥})), \\ \nabla^{2} f & = ε^{- 2} \partial_{ξ}^{2} f + ε^{- 1} (- \bar{K} \partial_{ξ} f) + ε^{0} (- ξ \bar{K^{2}} \partial_{ξ} f + \nabla_{⊥} \cdot \nabla_{⊥} f) + O (ε) \\ - \nabla \times \nabla \times u = ε^{- 2} \partial_{ξ}^{2} u_{⊥} + ε^{- 1} (- \bar{K} \partial_{ξ} u_{⊥} - \partial_{ξ} \nabla_{⊥} u_{σ} - \hat{n} \nabla_{⊥} \cdot (\partial_{ξ} u_{⊥})) + O (ε^{0}), \\ u \cdot \nabla f & = ε^{- 1} u_{σ} \partial_{ξ} f + \sum_{k = 0}^{\infty} ε^{k} ξ^{k} K^{k} u_{⊥} \cdot \nabla_{⊥} f, \\ u \cdot \nabla u & = ε^{- 1} (u_{σ} \partial_{ξ} u_{⊥} + \hat{n} u_{σ} \partial_{ξ} u_{σ}) \\ + \sum_{k = 0}^{\infty} ε^{k} ξ^{k} u_{⊥} K^{k} \cdot ((\nabla_{⊥} u_{⊥} - K u_{σ}) + (\nabla_{⊥} u_{σ} + K u_{⊥}) \hat{n}), \\ \partial_{t} & = - ε^{- 1} v \partial_{ξ} + ε^{0} \partial_{τ} + ε^{1} (ξ \nabla_{⊥} v \cdot \sum_{k = 0}^{\infty} ε^{k} ξ^{k} K^{k} \nabla_{⊥}) . \end{aligned}

Appendix B. Differential geometry of remapped coordinates

We solve §4b by remapping to a fixed rectangular domain with coordinates τ, ξ and ζ^±,

τ (t, x, z) = t, ξ (t, x, z) = x, ζ^{+} (t, x, z) = \frac{z}{h (t, x)}, ζ^{-} (t, x, z) = \frac{z - h (t, x)}{2 - h (t, x)}, η (τ, ξ) = h (t, x),

where η is the interface height in the new coordinates. We now develop the differential geometry required to write equations (2.1) and (2.4) in the new coordinates. We give explicit formulae for the ζ⁺ remapping. The ζ⁻ remapping is analogous and straightforward to derive. We use a tangent vector basis derived from the Jacobian

J

of our transformation and its inverse

K

J^{+} \equiv [\begin{matrix} \partial_{t} τ & \partial_{t} ξ & \partial_{t} ζ \\ \partial_{x} τ & \partial_{x} ξ & \partial_{x} ζ \\ \partial_{z} τ & \partial_{z} ξ & \partial_{z} ζ \end{matrix}] = [\begin{matrix} 1 & 0 & - ζ \frac{\partial_{τ} η}{η} \\ 0 & 1 & - ζ \frac{\partial_{ξ} η}{η} \\ 0 & 0 & \frac{1}{η} \end{matrix}] and K^{+} \equiv [\begin{matrix} \partial_{τ} t & \partial_{τ} x & \partial_{τ} z \\ \partial_{ξ} t & \partial_{ξ} x & \partial_{ξ} z \\ \partial_{ζ} t & \partial_{ζ} x & \partial_{ζ} z \end{matrix}] = [\begin{matrix} 1 & 0 & ζ \partial_{τ} η \\ 0 & 1 & ζ \partial_{ξ} η \\ 0 & 0 & η \end{matrix}] .

The tangent vector components are the rows of the spatial component of the inverse Jacobian $e_{i} = {K_{i}}^{j} {\hat{e}}_{j}$ . We use Einstein notation to sum over the spatial indices i = 1, 2. These tangent vectors induce dual vectors $ω^{i} = {J_{j}}^{i} {\hat{e}}^{j}$ which satisfy $ω^{i} \cdot e_{j} = δ_{j}^{i}$ , with components from the column vectors of the spatial component of the Jacobian. We record the length of the vectors using the metric, with co/contravariant components g_ij = e_i · e_j, g^ij = Ωⁱ · Ω^j,

\begin{aligned} e_{1}^{+} = {\hat{e}}_{1} + ζ \partial_{ξ} η {\hat{e}}_{2}, & e_{2}^{+} = η {\hat{e}}_{2}, {ω^{1}}^{+} = {\hat{e}}_{1}, {ω^{2}}^{+} = \frac{{\hat{e}}_{2} - ζ \partial_{ξ} η {\hat{e}}_{1}}{η}, \\ g_{11}^{+} = (1 + ζ^{2} \partial_{ξ} η^{2}), & g_{12}^{+} = ζ η \partial_{ξ} η, {g^{11}}^{+} = 1, {g^{12}}^{+} = - ζ \frac{\partial_{ξ} η}{η}, \\ g_{21}^{+} = ζ η \partial_{ξ} η, g_{22}^{+} = η^{2}, & {g^{21}}^{+} = - ζ \frac{\partial_{ξ} η}{η}, {g^{22}}^{+} = \frac{1}{η^{2}} (1 + ζ^{2} \partial_{ξ} η^{2}) . \end{aligned}

These give ‘Jacobian’ $| | K | |$ determinants of $| | K^{+} | | \equiv \sqrt{| | g_{i j} | |} = η$ . The completely antisymmetric tensor can be transformed from Cartesian coordinates,

E = [i j] {\hat{ω}}^{i} {\hat{ω}}^{j} = [i j] K_{k}^{i} K_{l}^{j} ω^{k} ω^{l} = | | K | | [k l] ω^{k} ω^{l} = \frac{1}{| | K | |} [k l] e_{k} e_{l},

where [ij] is the antisymmetric symbol. We project the gradient to the tangent basis with the covariant derivative

\nabla_{i} = e_{i} \cdot \nabla ⟺ \nabla = ω^{i} \nabla_{i}

. We measure spatial variation of the basis with the connection coefficients

Γ_{i j}^{k} = ω^{k} \cdot \nabla_{j} e_{i} ⟺ e_{j} \cdot \nabla e_{i} = Γ_{i j}^{k} e_{k}

. For the tangent basis, the connection coefficients take a particularly simple and symmetric form

[{Γ_{11}^{1}}^{+} = {Γ_{12}^{1}}^{+} = {Γ_{21}^{1}}^{+} = {Γ_{22}^{1}}^{+} = 0, {Γ_{11}^{2}}^{+} = \frac{ζ \partial_{ξ}^{2} η}{η}, {Γ_{12}^{2}}^{+} = {Γ_{21}^{2}}^{+} = \frac{\partial_{ξ} η}{η}, {Γ_{22}^{2}}^{+} = 0.

We note the connection coefficients for the dual vectors are related to those for the tangent basis by $0 = \nabla_{k} (ω^{i} \cdot e_{j}) = \nabla_{k} (ω^{i}) \cdot e_{j} + ω^{i} \cdot \nabla_{k} (e_{j}) = \nabla_{k} (ω^{i}) \cdot e_{j} + Γ_{j k}^{i}$ . The partial time derivative ∂_t changes as defined in zeroth row of the Jacobian, and the basis vectors also evolve in time,

\partial_{τ} e_{1}^{+} \equiv {Γ_{10}^{i}}^{+} e_{i}^{+} = \frac{ζ}{η} \partial_{τ} \partial_{ξ} η e_{2}^{+}, \partial_{τ} e_{2}^{+} \equiv {Γ_{20}^{i}}^{+} e_{i}^{+} = \frac{\partial_{τ} η}{η} e_{2}^{+} .

These geometric quantities allow us to calculate all the relevant vector calculus quantities,

\begin{aligned} \nabla u & = \nabla_{i} (u^{j} e_{j}) ω^{i} = (\nabla_{i} u^{j} + u^{k} Γ_{k i}^{j}) e_{j} ω^{i}, \\ \nabla \cdot u & = \nabla u : I = (\nabla_{i} u^{j} + u^{k} Γ_{k i}^{j}) e_{j} \cdot ω^{i} = \nabla_{i} u^{i} + u^{k} Γ_{k i}^{i}, \\ \nabla \times u & = E^{⊤} : \nabla u = (E_{i j} ω^{i} ω^{j}) : ω^{k} (\nabla_{k}) (u^{l} e_{l}) = | | K | | [i j] u_{; k}^{l} (ω^{i} \cdot ω^{k}) (ω^{j} \cdot e_{l}) = | | K | | [i j] g^{i k} u_{; k}^{j}, \\ \nabla \times q & = E^{⊤} \cdot \nabla q = (E^{i j} \nabla_{k} q) (e_{j} \cdot ω^{k}) e_{i} = \frac{1}{| | K | |} [i j] q_{, j} e_{i} = \frac{1}{| | K | |} (q_{, 2} e_{1} - q_{, 1} e_{2}), \\ \nabla^{2} f & = \nabla_{i} (\nabla_{j} f ω^{j}) ω^{i} : I = (\nabla_{i} \nabla_{j} f) ω^{j} \cdot ω^{i} - \nabla_{j} f Γ_{k i}^{j} ω^{k} \cdot ω^{i} = g^{j i} \nabla_{i} \nabla_{j} f - g^{k i} \nabla_{j} f Γ_{k i}^{j}, \\ \nabla p & = g^{j i} p_{, j} e_{i}, \\ u \times q & = E \cdot u q = (E_{i j} ω^{i} ω^{j}) \cdot (q u^{k} e_{k}) = | | K | | g^{j k} [i j] u^{i} e_{k} = | | K | | q (g^{2 k} u^{1} - g^{1 k} u^{2}) e_{k}, \\ \nabla \times \nabla \times u & = [i j] \partial_{j} ([k l] \partial_{k} u_{l}) {\hat{e}}_{i} = (δ_{i k} δ_{j l} - δ_{i l} δ_{j k}) \partial_{j} \partial_{k} u_{l} {\hat{e}}_{i} = \nabla (\nabla \cdot u) - \nabla^{2} u, \\ (\nabla \times u) \times u & = [i j] [k l] \partial_{k} u_{l} u_{j} {\hat{e}}_{i} = (\partial_{i} u_{j} - \partial_{j} u_{i}) u_{j} {\hat{e}}_{i} = \frac{1}{2} \nabla (u \cdot u) - (u \cdot \nabla) u, \\ \nabla \cdot (u u) & = (u \cdot \nabla) u + (\nabla \cdot u) u, \\ \partial_{t} & = \partial_{τ} + {J_{0}}^{i} \nabla_{i}, \\ \partial_{t} u & = \partial_{t} u^{i} e_{i} + u^{i} \partial_{t} e_{i} = (\partial_{τ} + J_{0}^{i} \nabla_{i}) u^{j} e_{j} + {J_{0}}^{i \geq 0} Γ_{k i \geq 0}^{j} u^{k} e_{j} . \end{aligned}

The normal vector at the interface is proportional to the ζ dual vector. This gives us the normal gradient and normal velocity at the interface,

\hat{n} = \frac{ω^{2}}{\sqrt{g^{22}}} \hat{n} \cdot \nabla = \frac{g^{2 i}}{\sqrt{g^{22}}} \nabla_{i}, v = \partial_{t} h {\hat{e}}_{z}, \hat{n} \cdot v = \frac{1}{\sqrt{g^{22}}} \frac{\partial_{t} h}{h} .

We finally write the interfacial curvature as $κ = \partial_{ξ}^{2} η / {(1 + \partial_{ξ} η^{2})}^{3 / 2}$ . This completes the relations used to simulate the sharp interface and phase-field equations in remapped geometries.

Footnotes

1 Shortly before submission, we became aware of the work [58]. In it, Subhedar et al. perform second-order analysis of a phase-field model combining melting and advection. Our work is more general as we also account for dissolution, give a more comprehensive analytical treatment, and use more challenging computational benchmarks.

2 We also provide a Mathematica script that automates each step of this analysis at github.com/ericwhester/phase-field-code.

3 The full simulation code, saved data and analysis scripts are freely available at github.com/ericwhester/phase-field-code.

References

Reference

1.
Epstein M, Cheung FB. 1983 Complex freezing-melting interfaces in fluid flow. Annu. Rev. Fluid Mech. 15, 293–319. (doi:10.1146/annurev.fl.15.010183.001453) Crossref, Web of Science, Google Scholar
2.
Glicksman ME, Coriell SR, McFadden GB. 1986 Interaction of flows with the crystal–melt interface. Annu. Rev. Fluid Mech. 18, 307–335. (doi:10.1146/annurev.fl.18.010186.001515) Crossref, Web of Science, Google Scholar
3.
Cahn JW, Hilliard JE. 1958 Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 28, 258–267. (doi:10.1063/1.1744102) Crossref, Web of Science, Google Scholar
4.
Huppert HE. 2002 Geological fluid mechanics. In Perspectives in fluid dynamics: a collective introduction to current research (ed. GK Batchelor). Cambridge, UK: Cambridge University Press. Google Scholar
5.
Turner JS. 1985 Multicomponent convection. Annu. Rev. Fluid Mech. 17, 11–44. (doi:10.1146/annurev.fl.17.010185.000303) Crossref, Web of Science, Google Scholar
6.
Knobloch E, Krechetnikov R. 2015 Problems on time-varying domains: formulation, dynamics, and challenges. Acta Appl. Math. 137, 123–157. (doi:10.1007/s10440-014-9993-x) Crossref, Web of Science, Google Scholar
7.
Worster MG. 2002 Solidification of fluids. In Perspectives in fluid dynamics: a collective introduction to current research (ed. GK Batchelor). Cambridge, UK: Cambridge University Press. Google Scholar
8.
Donea J, Huerta A, Ponthot JP, Rodríguez-Ferran A. 2004 Arbitrary Lagrangian–Eulerian methods. In Encyclopedia of computational mechanics. Atlanta, GA: American Cancer Society. Google Scholar
9.
Hadžić M, Shkoller S. 2017 Well-posedness for the classical Stefan problem and the zero surface tension limit. Arch. Ration. Mech. Anal. 223, 213–264. (doi:10.1007/s00205-016-1041-8) Crossref, Web of Science, Google Scholar
10.
Gibbs JW. 1878 On the equilibrium of heterogeneous substances. Am. J. Sci. s3-16, 441–458. (doi:10.2475/ajs.s3-16.96.441) Crossref, Google Scholar
11.
van der Waals JD. 1979 The thermodynamic theory of capillarity under the hypothesis of a continuous variation of density. J. Stat. Phys. 20, 200–244. (doi:10.1007/BF01011514) Crossref, Web of Science, Google Scholar
12.
Stefan J. 1891 Ueber Die Theorie Der Eisbildung, Insbesondere Über Die Eisbildung Im Polarmeere. Annalen der Physik 278, 269–286. (doi:10.1002/andp.18912780206) Crossref, Google Scholar
13.
Boettinger WJ, Warren JA, Beckermann C, Karma A. 2002 Phase-field simulation of solidification. Annu. Rev. Mater. Res. 32, 163–194. (doi:10.1146/annurev.matsci.32.101901.155803) Crossref, Web of Science, Google Scholar
14.
Steinbach I. 2009 Phase-field models in materials science. Modell. Simul. Mater. Sci. Eng. 17, 073001. (doi:10.1088/0965-0393/17/7/073001) Crossref, Web of Science, Google Scholar
15.
Voller VR, Swaminathan CR, Thomas BG. 1990 Fixed grid techniques for phase change problems: a review. Int. J. Numer. Methods Eng. 30, 875–898. (doi:10.1002/nme.1620300419) Crossref, Web of Science, Google Scholar
16.
Ulvrová M, Labrosse S, Coltice N, Råback P, Tackley PJ. 2012 Numerical modelling of convection interacting with a melting and solidification front: application to the thermal evolution of the basal magma ocean. Phys. Earth Planet. Inter. 206–207, 51–66. (doi:10.1016/j.pepi.2012.06.008) Crossref, Web of Science, Google Scholar
17.
Osher S, Fedkiw RP. 2001 Level set methods: an overview and some recent results. J. Comput. Phys. 169, 463–502. (doi:10.1006/jcph.2000.6636) Crossref, Web of Science, Google Scholar
18.
Chen S, Merriman B, Osher S, Smereka P. 1997 A simple level set method for solving Stefan problems. J. Comput. Phys. 135, 8–29. (doi:10.1006/jcph.1997.5721) Crossref, Web of Science, Google Scholar
19.
Li X, Lowengrub J, Ratz A, Voigt A. 2009 Solving PDEs in complex geometries: a diffuse domain approach. Commun. Math. Sci. 7, 81–107. (doi:10.4310/CMS.2009.v7.n1.a4) Crossref, PubMed, Web of Science, Google Scholar
20.
Aland S, Lowengrub J, Voigt A. 2010 Two-phase flow in complex geometries: a diffuse domain approach. Comput. Model. Eng. Sci.: CMES 57, 77–106. PubMed, Web of Science, Google Scholar
21.
Mac Huang J, Shelley MJ, Stein DB. 2020 A stable and accurate scheme for solving the Stefan problem coupled with natural convection using the immersed boundary smooth extension method. (http://arxiv.org/abs/2006.04736 [physics]) Google Scholar
22.
Fix GJ. 1982 Phase field methods for free boundary problems. Google Scholar
23.
Langer JS. 1986 Models of pattern formation in first-order phase transitions. In Directions in condensed matter physics (eds G Grinstein, G Mazenko), pp. 165–186. Singapore: World Scientific. (doi:10.1142/9789814415309_0005) Google Scholar
24.
Hohenberg PC, Halperin BI. 1977 Theory of dynamic critical phenomena. Rev. Mod. Phys. 49, 435–479. (doi:10.1103/RevModPhys.49.435) Crossref, Web of Science, Google Scholar
25.
Caginalp G. 1989 Stefan and Hele-Shaw type models as asymptotic limits of the phase-field equations. Phys. Rev. A 39, 5887–5896. (doi:10.1103/PhysRevA.39.5887) Crossref, Web of Science, Google Scholar
26.
Caginalp G, Chen X. 1998 Convergence of the phase field model to its sharp interface limits. Eur. J. Appl. Math. 9, 417–445. (doi:10.1017/S0956792598003520) Crossref, Web of Science, Google Scholar
27.
Wheeler AA, Boettinger WJ, McFadden GB. 1992 Phase-field model for isothermal phase transitions in binary alloys. Phys. Rev. A 45, 7424–7439. (doi:10.1103/PhysRevA.45.7424) Crossref, PubMed, Web of Science, Google Scholar
28.
Bi Z, Sekerka RF. 1998 Phase field model of solidification of a binary alloy. Physica A 261, 12. (doi:10.1016/S0378-4371(98)00364-1) Crossref, Web of Science, Google Scholar
29.
Karma A. 2001 Phase-field formulation for quantitative modeling of alloy solidification. Phys. Rev. Lett. 87, 115701. (doi:10.1103/physrevlett.87.115701) Crossref, PubMed, Web of Science, Google Scholar
30.
Ramirez JC, Beckermann C, Karma A, Diepers HJ. 2004 Phase-field modeling of binary alloy solidification with coupled heat and solute diffusion. Phys. Rev. E 69, 051607. (doi:10.1103/PhysRevE.69.051607) Crossref, Web of Science, Google Scholar
31.
Beckermann C, Diepers HJ, Steinbach I, Karma A, Tong X. 1999 Modeling melt convection in phase-field simulations of solidification. J. Comput. Phys. 154, 468–496. (doi:10.1006/jcph.1999.6323) Crossref, Web of Science, Google Scholar
32.
Anderson DM, McFadden GB, Wheeler AA. 2000 A phase-field model of solidification with convection. Physica D 135, 175–194. (doi:10.1016/S0167-2789(99)00109-8) Crossref, Web of Science, Google Scholar
33.
Mokbel D, Abels H, Aland S. 2018 A phase-field model for fluid–structure interaction. J. Comput. Phys. 372, 823–840. (doi:10.1016/j.jcp.2018.06.063) Crossref, Web of Science, Google Scholar
34.
Penrose O, Fife PC. 1990 Thermodynamically consistent models of phase-field type for the kinetic of phase transitions. Physica D 43, 44–62. (doi:10.1016/0167-2789(90)90015-H) Crossref, Web of Science, Google Scholar
35.
Wang SL, Sekerka RF, Wheeler AA, Murray BT, Coriell SR, Braun RJ, McFadden GB. 1993 Thermodynamically-consistent phase-field models for solidification. Physica D 69, 189–200. (doi:10.1016/0167-2789(93)90189-8) Crossref, Web of Science, Google Scholar
36.
McFadden GB, Wheeler AA, Anderson DM. 2000 Thin interface asymptotics for an energy/entropy approach to phase-field models with unequal conductivities. Physica D 144, 154–168. (doi:10.1016/S0167-2789(00)00064-6) Crossref, Web of Science, Google Scholar
37.
Ohno M, Takaki T, Shibuta Y. 2016 Variational formulation and numerical accuracy of a quantitative phase-field model for binary alloy solidification with two-sided diffusion. Phys. Rev. E 93, 012802. (doi:10.1103/PhysRevE.93.012802) Crossref, PubMed, Web of Science, Google Scholar
38.
Bollada P, Jimack P, Mullis A. 2017 Bracket formalism applied to phase field models of alloy solidification. Comput. Mater. Sci. 126, 426–437. (doi:10.1016/j.commatsci.2016.09.036) Crossref, Web of Science, Google Scholar
39.
Caginalp G. 1986 An analysis of a phase field model of a free boundary. Arch. Ration. Mech. Anal. 92, 205–245. (doi:10.1007/BF00254827) Crossref, Web of Science, Google Scholar
40.
Caginalp G, Fife P. 1988 Dynamics of layered interfaces arising from phase boundaries. SIAM J. Appl. Math. 48, 506–518. (doi:10.1137/0148029) Crossref, Web of Science, Google Scholar
41.
Lin JT. 1988 The numerical analysis of a phase field model in moving boundary problems. SIAM J. Numer. Anal. 25, 1015–1031. (doi:10.1137/0725058) Crossref, Web of Science, Google Scholar
42.
Wheeler AA, Murray BT, Schaefer RJ. 1993 Computation of dendrites using a phase field model. Physica D 66, 243–262. (doi:10.1016/0167-2789(93)90242-S) Crossref, Web of Science, Google Scholar
43.
Fabbri M, Voller VR. 1997 The phase-field method in the sharp-interface limit: a comparison between model potentials. J. Comput. Phys. 130, 256–265. (doi:10.1006/jcph.1996.5585) Crossref, Web of Science, Google Scholar
44.
Tong X, Beckermann C, Karma A, Li Q. 2001 Phase-field simulations of dendritic crystal growth in a forced flow. Phys. Rev. E 63, 061601. (doi:10.1103/PhysRevE.63.061601) Crossref, Web of Science, Google Scholar
45.
Lu Y, Beckermann C, Ramirez J. 2005 Three-dimensional phase-field simulations of the effect of convection on free dendritic growth. J. Cryst. Growth 280, 320–334. (doi:10.1016/j.jcrysgro.2005.03.063) Crossref, Web of Science, Google Scholar
46.
Jokisaari AM, Voorhees PW, Guyer JE, Warren J, Heinonen OG. 2017 Benchmark problems for numerical implementations of phase field models. Comput. Mater. Sci. 126, 139–151. (doi:10.1016/j.commatsci.2016.09.022) Crossref, Web of Science, Google Scholar
47.
Favier B, Purseed J, Duchemin L. 2019 Rayleigh–Bénard convection with a melting boundary. J. Fluid Mech. 858, 437–473. (doi:10.1017/jfm.2018.773) Crossref, Web of Science, Google Scholar
48.
Purseed J, Favier B, Duchemin L, Hester EW. 2020 Bistability in Rayleigh–Benard convection with a melting boundary. Phys. Rev. Fluids 5, 023501. (doi:10.1103/PhysRevFluids.5.023501) Crossref, Web of Science, Google Scholar
49.
Couston LA, Hester E, Favier B, Taylor JR, Holland PR, Jenkins A. 2020 Topography generation by melting and freezing in a turbulent shear flow. (http://arxiv.org/abs/2004.09879 [physics]) Google Scholar
50.
Magaletti F, Picano F, Chinappi M, Marino L, Casciola CM. 2013 The sharp-interface limit of the Cahn–Hilliard/Navier–Stokes model for binary fluids. J. Fluid Mech. 714, 95–126. (doi:10.1017/jfm.2012.461) Crossref, Web of Science, Google Scholar
51.
Soligo G, Roccon A, Soldati A. 2019 Breakage, coalescence and size distribution of surfactant-laden droplets in turbulent flow. J. Fluid Mech. 881, 244–282. (doi:10.1017/jfm.2019.772) Crossref, Web of Science, Google Scholar
52.
Karma A, Rappel WJ. 1996 Phase-field method for computationally efficient modeling of solidification with arbitrary interface kinetics. Phys. Rev. E 53, R3017–R3020. (doi:10.1103/PhysRevE.53.R3017) Crossref, Web of Science, Google Scholar
53.
Karma A, Rappel WJ. 1998 Quantitative phase-field modeling of dendritic growth in two and three dimensions. Phys. Rev. E 57, 4323–4349. (doi:10.1103/PhysRevE.57.4323) Crossref, Web of Science, Google Scholar
54.
Almgren RF. 1999 Second-order phase field asymptotics for unequal conductivities. SIAM J. Appl. Math. 59, 2086–2107. (doi:10.1137/S0036139997330027) Crossref, Web of Science, Google Scholar
55.
Garcke H, Stinner B. 2006 Second order phase field asymptotics for multi-component systems. Interfaces Free Boundaries 8, 131–157. (doi:10.4171/IFB/138) Crossref, Web of Science, Google Scholar
56.
Ohno M, Matsuura K. 2009 Quantitative phase-field modeling for dilute alloy solidification involving diffusion in the solid. Phys. Rev. E 79, 031603. (doi:10.1103/PhysRevE.79.031603) Crossref, Web of Science, Google Scholar
57.
Fife PC. 1988 Dynamics of internal layers and diffusive interfaces. Number 53 in CBMS-NSF Regional Conference Series in Applied Mathematics. Philadelphia, PA: Society for Industrial and Applied Mathematics. Google Scholar
58.
Subhedar A, Galenko PK, Varnik F. 2020 Diffuse interface models of solidification with convection: the choice of a finite interface thickness. Eur. Phys. J. Spec. Top. 229, 447–452. (doi:10.1140/epjst/e2019-900099-5) Crossref, Web of Science, Google Scholar
59.
Chen X, Caginalp G, Eck C. 2006 A rapidly converging phase field model. Discrete Continuous Dyn. Syst. 15, 1017–1034. (doi:10.3934/dcds.2006.15.1017) Crossref, Web of Science, Google Scholar
60.
Kockelkoren J, Levine H, Rappel WJ. 2003 Computational approach for modeling intra- and extracellular dynamics. Phys. Rev. E 68, 037702. (doi:10.1103/PhysRevE.68.037702) Crossref, Web of Science, Google Scholar
61.
Hester EW, Vasil GM, Burns KJ. 2019 Improving convergence of volume penalized fluid-solid interactions. (http://arxiv.org/abs/1903.11914 [math]) Google Scholar
62.
Burns KJ, Vasil GM, Oishi JS, Lecoanet D, Brown BP. 2020 Dedalus: a flexible framework for numerical simulations with spectral methods. Phys. Rev. Res. 2, 023068. (doi:10.1103/PhysRevResearch.2.023068) Crossref, Google Scholar
63.
Le Bars M, Worster MG. 2006 Interfacial conditions between a pure fluid and a porous medium: implications for binary alloy solidification. J. Fluid Mech. 550, 149–173. (doi:10.1017/S0022112005007998) Crossref, Web of Science, Google Scholar
64.
Boyd JP. 2001 Chebyshev and Fourier spectral methods, 2nd revised edn. New York, NY: Dover Publications. Google Scholar
65.
Subich CJ, Lamb KG, Stastna M. 2013 Simulation of the Navier–Stokes equations in three dimensions with a spectral collocation method. Int. J. Numer. Methods Fluids 73, 103–129. (doi:10.1002/fld.3788) Crossref, Web of Science, Google Scholar

VIEW FULL TEXT DOWNLOAD PDF

This Issue

Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences cover image

October 2020

Volume 476Issue 2242

Article Information

DOI:https://doi.org/10.1098/rspa.2020.0508
PubMed:33223948
Published by:Royal Society
Online ISSN:1471-2946

History:

Manuscript received29/06/2020
Manuscript accepted23/09/2020
Published online21/10/2020
Published in print28/10/2020

License:

Citations and impact

Keywords

Subjects

Your password must have 8 characters or more and contain 3 of the following:

Improved phase-field models of melting and dissolution in multi-component flows

Abstract