Diffuse interface methods

The phase-field model for a two-phase system is based on the definition of a phase-field parameter \(C\), which attains values between limiting \(C\) values. Away from the interface, \(C\) attains these limiting values, while there is a smooth transition between these limiting values around the interface. A phase-field method defines an energy functional that is commonly the following Helmholtz (free) energy functional

(46)\[\mathcal{F} = \int_{\Omega} f \, \text{d}V = \int_{\Omega} \frac{\sigma}{\epsilon} \left( \Psi(C) + \frac{\epsilon^2}{2} \left(\nabla C \right)^2\right) \text{d}V,\]

where \(f\) is the Helmholtz energy density, \(\sigma\) denotes the surface tension coefficient, and \(\epsilon\) determines the width of the diffuse interface. The two terms in the energy functional account for the bulk energy density and the interfacial energy density, respectively. The latter is the excess energy due to an inhomogeneous distribution of \(C\) in the interfacial region. The phobic bulk energy density \(\Psi(C)\) is a double-well potential This double-well potential has minima at \(C_{min}\) and \(C_{max}\), and tries to sharpen the interface by separating the phases. On the other hand, the philic gradient term tries to diffuse the interface into the bulk regions. At minimum energy both terms balance, and the interface profile reaches its equilibrium shape. Another way to interpret this minimum energy is through the chemical potential \(\omega\), which is given by the functional derivative of the energy density Eq.46 with respect the phase-field \(C\). The functional derivative is a generalization of the Euler-Lagrange equation, and it can be written as

(47)\[\omega = \frac{\partial f}{\partial C} = \frac{\partial f}{\partial C} - \nabla \cdot \frac{\partial f}{\partial \nabla C}.\]

Once this chemical potential has becomes uniform in space, the energy density has reached its minimum. It is important to note that the original works by Van der Waals [18] and Cahn and Hilliard [4,3,5] corresponds to a one dimensional case (thereby assuming a flat interface). In that case it holds that \(\nabla C=dC/dz\) and \(\Delta C=d^2C/dz^2\), where \(z\) is a spatial location. Later [when? by who?] this was extended to higher dimensions, which usually leads to the following \(\omega_c\) (where \(c\) signifies classic approach)

(48)\[\omega_c = \frac{\partial f}{\partial C} = \frac{\sigma}{\epsilon} \left(\frac{\partial \Psi}{\partial C} - \epsilon^2 \Delta C \right).\]

As intended, when applied to one-dimensional problems, this chemical potential enforces the balance between the sharpening and diffusing terms while conserving (enclosed) mass. At the same time, it is well known that this chemical potential can cause nonphysical (enclosed) mass loss and bulk diffusion when applied to higher-dimensional problems. Therefore, a substantial amount of studies have tried to remedy these nonphysical effects in various ways [2,21,14,9,19,13,24,17,20,23]. Although some of these remedies appear to suppress the nonphysical effects, so far no theoretical cause and prevention of them has been put forward [7]. However, it is possible to attribute the (enclosed) mass loss and bulk diffusion directly to the extension of the original one-dimensional formulation (corresponding to a flat interface) to higher dimensions.

Without loss of generality, it can be assumed that the extension of the energy functional to higher dimensions should consider the gradient of the single phase-field \(C\) normal to the interface. The classic Helmholtz energy functional of Eq.46 is then redefined to

(49)\[\mathcal{F}_b = \int_{\Omega} f_b \, \text{d}V = \int_{\Omega} \frac{\sigma}{\epsilon} \left( \Psi(C) + \frac{\epsilon^2}{2} \left( \nabla_n C \right)^2 \right) \text{d}V,\]

where \(\nabla_n\) is the gradient normal to the interface. The fundamental difference is that here \(\nabla_n C\) (a scalar) is used, compared to \(\nabla C\) (a vector) in the classic Helmholtz energy functional of Eq.46. This subtle modification of the energy functional leads to the following balanced chemical potential

(50)\[\omega_b = \frac{\partial f_b}{\partial C} = \frac{\sigma}{\epsilon} \left( \frac{\partial \Psi}{\partial C} - \epsilon^2 \frac{\partial^2C}{\partial n^2}\right),\]

which is identically zero for an interface (profile) at equilibrium (corresponding to Eq.53).

During the derivation of the classic chemical potential Eq.58, the divergence of \(\nabla C\) leads to an additional curvature term. This can be understood by considering the difference between the two chemical potentials

(51)\[\omega_c - \omega_b = - \sigma \epsilon \left( \Delta C - \frac{\partial^2 C}{\partial n^2} \right) = \sigma \kappa \epsilon \left| \nabla C \right|,\]

where \(\kappa\) is the curvature. This relation also reflects that mass loss and bulk diffusion for the classic chemical potential only vanish if \(\epsilon \to 0\), the so-called sharp interface limit [10]. From a physical perspective, minimization of the redefined Helmholtz energy functional effectuates only the recovery of the equilibrium interface profile. This dynamics corresponds to the original one-dimensional formulation, and does in general not depend on the interface thickness \(\epsilon\). On the other hand, minimization of the classic Helmholtz energy functional effectuates a balance between the recovery of the equilibrium interface profile and an uniform distribution of the right hand side term of Eq.51. However, neither of these goals is accomplished completely, and the overall dynamics strongly depends on the interface thickness \(\epsilon\).

This can be understood better by combining the difference between the chemical potentials Eq.51 with the uncoupled (from the NS equation) CH equation Eq.57. For a curved interface with an equilibrium hyperbolic tangent profile and constant \(M_C\) this gives

(52)\[\frac{\partial C}{\partial t} = M_C \sigma \epsilon \Delta \left( \kappa \left| \nabla C \right| \right),\]

which shows that, in order to make \(\kappa \left| \nabla C \right|\) more uniform in space, \(C\) (and therefore the hyperbolic tangent interface profile) will change. Based on csf, \(\kappa \left| \nabla C \right|\) can be considered as a pressure jump across (or a surface tension force at) the interface in normal direction. From a physical perspective, this interface normal (pressure) term is balanced by surface tension forces parallel to the interface (through the Young-Laplace equation), and must not be diffused. Therefore, as the use of this additional diffusion term is physically problematic, it should be removed. The proposed balanced chemical potential does just that, while respecting the original one-dimensional dynamics. As the multi-dimensional balanced energy functional is essentially identical to the one-dimensional version, application of the energy argument [22] yields an equilibrium for the hyperbolic tangent profile given by Eq.53 without any need for droplet shrinkage. Moreover, by removing the source of the mass loss and bulk diffusion the need for correction terms or special settings has vanished.

Equilibrium interface profile

The equilibrium phase-field \(C\) can be written as a function of the signed distance function \(d(\mathbf{x})\) between the interface and the location \(\mathbf{x}\). In its most general form this results in

(53)\[C(\mathbf{x}) = a \left( b + \tanh \left( \frac{d(\mathbf{x})}{\gamma} \right) \right),\]

where \(a\) and \(b\) are constant values, and \(\gamma\) defines a scaled thickness of the diffuse interface (denoted by \(\epsilon\)). The constant values are determined by the defined limiting \(C\) values, and common settings are:

Table 1 Common settings to define the diffuse interface

option	range of \(C\)	\(a\)	\(b\)	\(\gamma\)
A	\(-1 \le C \le 1\)	1	0	\(\sqrt{2} \epsilon\)
B	\(\phantom{-}0 \le C \le 1\)	1/2	1	\(2 \sqrt{2} \epsilon\)

From a physical perspective both options are identical, however the bulk energy density and some scaling parameters will differ.

There is a strong relation between the Level-Set (LS) and the phase-field function \(C\). The LS is generally defined as a signed distance function, so essentially it is equivalent to the \(d(\mathbf{x})\). By rewriting Eq.53 and application of the inverse hyperbolic tangent function, the following derivation can be obtained

(54)\[\begin{split}\begin{split} \frac{C(\mathbf{x})}{a} - b &= \tanh \left( \frac{d(\mathbf{x})}{\gamma} \right) \notag \\ \frac{1}{2} \ln \left( \frac{1 + \frac{C(\mathbf{x})}{a} - b}{1 - \frac{C(\mathbf{x})}{a} + b} \right) &= \frac{d(\mathbf{x})}{\gamma} \notag \\ d( \mathbf{x} ) = LS( \mathbf{x} ) &= \frac{\gamma}{2} \ln \left( \frac{1 + \frac{C(\mathbf{x})}{a} - b}{1 - \frac{C(\mathbf{x})}{a} + b} \right). \end{split}\end{split}\]

For the two options A and B this results in

Table 2 Relation between phase-field function \(C\) and the Level-Set function

option	range of \(C\)	\(C(\mathbf{x})\)	\(LS(\mathbf{x})\)
A	\(-1 \le C \le 1\)	\(\phantom{\frac{1}{2} + \frac{1}{2}} \tanh \left( \frac{d(\mathbf{x})}{\sqrt{2} \epsilon} \right)\)	\(\frac{\epsilon}{\sqrt{2}} \ln \left( \frac{1+C(\mathbf{x})}{1-C(\mathbf{x})} \right)\)
B	\(\phantom{-}0 \le C \le 1\)	\(\frac{1}{2} + \frac{1}{2} \tanh \left( \frac{d(\mathbf{x})}{2\sqrt{2} \epsilon} \right)\)	\(\sqrt{2} \epsilon \ln \left( \frac{C(\mathbf{x})}{1 - C(\mathbf{x})} \right)\)

Bulk energy density

Depending on the chosen equilibrium interface profile (option A or B), the bulk energy density as defined in energy_functional is defined by

Table 3 Bulk energy density functions

option	range of \(C\)	energy density function \(\Psi\)
A	\(-1 \le C \le 1\)	\((C+1)^2 (1-C)^2\)
B	\(\phantom{-}0 \le C \le 1\)	\(\frac{1}{4} C^2 (1-C)^2\)

Phase-field evolution in time

The evolution of a phase-field \(C(\mathbf{x},t)\) is mathematically given by the following conservation equation

(55)\[\frac{\partial C}{\partial t} + \nabla \cdot \left( \mathbf{u} C \right) = -\nabla_{\mathcal{H}} f(C),\]

where the advection is controlled by the fluid velocity \(\mathbf{u}\) and the right hand side is given by a so-called gradient flow of the energy density \(f(C)\). For \(\mathbf{u}=0\) the uncoupled (from the NS equation) phase-field equation is obtained, where the evolution of \(C\) is only affected by the term on the right hand side. Here \(\nabla_{\mathcal{H}}\) denotes the (Hilbert space) gradient of \(f\) at \(C\), which ensures that \(C\) decreases along the gradient of \(f\) on \(\mathcal{H}\). The precise expression for \(\nabla_{\mathcal{H}}\) can be found from a thermodynamic perspective [15] or by a variational formulation (see appref{ap:gradientFlow}). Two well-known results are the Allen-Cahn (AC) equation (so-called \(L^2\) gradient flow)

(56)\[\frac{\partial C}{\partial t} + \nabla \cdot \left( \mathbf{u} C \right) = -M_A \omega,\]

and the Cahn-Hilliard (CH) equation (so-called \(H^{-1}\) gradient flow)

(57)\[\frac{\partial C}{\partial t} + \nabla \cdot \left( \mathbf{u} C \right) = \nabla \cdot \left( M_C \nabla \omega \right).\]

Here \(M\) is the interface mobility (often set to unity or a function of \(C\)), \(\omega\) the chemical potential, and the subscripts \(A\) and \(C\) denote Allen-Cahn and Cahn-Hilliard, respectively.

Eq.47 reduces to

(58)\[\omega_c = \frac{\sigma}{\epsilon} \left( C^3 - \frac{3}{2}C^2+ \frac{1}{2} C - \epsilon^2 \Delta C \right).\]