Benchmark Problem 2: Ostwald Ripening¶

• Free energy and dynamics
  • Figure 1: Free energy density surface
• Parameter values
• Domain geometries and boundary conditions
  • (a) Square periodic
  • (b) Square no-flux
  • (c) T-shape
  • (d) Sphere
• Initial conditions
  • Initial conditions for (a), (b) and (c)
    • Figure 2: initial $c$ for (a), (b) and (c)
    • Figure 3: initial $\eta_1$ for (a), (b) and (c)
  • Initial conditions for (d)
    • Figure 4: initial $c$ for (d)
    • Figure 5: initial $\eta_1$ for (d)
• Results

See the journal publication entitled "Benchmark problems for numerical implementations of phase field models" for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems.

Free energy and dynamics¶

The atomic fraction of solute is specified by the conserved variable $c$, while the phase is indicated by a structural order parameter, $\eta$. The structural order parameter is non-conserved and is a phenomenological phase descriptor, such that the $\alpha$ phase is indicated by $\eta=0$, while the $\beta$ phase is indicated by $\eta=1$. If multiple energetically equivalent orientation variants exist (for example, due to crystallographic symmetry considerations or ordered and disordered phases), the model may include $p$ number of structural order parameters, $\eta_{p}$, with one for each orientation variant. We include a nontrivial number of order parameters by setting $p=4$, a value commonly used in superalloy models; this will stress the numerical solver while not making the problem intractable.

In this benchmark problem, the free energy of the system is based on a more general formulation of that presented in Ref. ZHU. Specifically, the free energy here includes the concentration gradient energy as outlined in WBM1 and is expressed as

$$ F=\int_{V}\left(f_{chem}\left(c,\eta_{1},...\eta_{p}\right)+\frac{\kappa_{c}}{2}|\nabla c|^{2}+\sum_{i=1}^{p}\frac{\kappa_{\eta}}{2}|\nabla\eta_{i}|^{2}\right)dV $$

where $\kappa_{c}$ and $\kappa_{\eta}$ are the gradient energy coefficients for $c$ and $\eta_{i}$, respectively. While the model in Ref. ZHU follows the Kim-Kim-Suzuki (KKS) formulation for interfacial energy, we use the Wheeler-Boettinger-McFadden (WBM) formulation for simplicity. In the KKS model, the interface is treated as an equilibrium mixture of two phases with fixed compositions such that an arbitrary diffuse interface width may be specified for a given interfacial energy. In the WBM model, interfacial energy and interfacial width are linked with the concentration, such that very high resolution across the interface may be required to incorporate accurate interfacial energies.

The formulation for $f_{chem}$ in Ref. ZHU is adapted for our benchmark problem as

\begin{equation} f_{chem}\left(c,\eta_{1},...\eta_{p}\right)= f^{\alpha}\left(c\right)\left[1-h\left(\eta_{1}, ...\eta_{p}\right)\right]+f^{\beta}\left(c\right)h\left(\eta_{1}, ...\eta_{p}\right)+wg\left(\eta_{1}, ...\eta_{p}\right), \end{equation}

where $f^{\alpha}$ and $f^{\beta}$ are the chemical free energy densities of the $\alpha$ and $\beta$ phases, respectively, $h\left(\eta_{1},...\eta_{p}\right)$ is an interpolation function, and $g\left(\eta_{1},...\eta_{p}\right)$ is a double-well function. The function $h$ increases monotonically between $h(0)=0$ and $h(1)=1$, while the function $g$ has minima at $g(0)=0$ and $g(1)=0$. The height of the double well barrier is controlled by $w$. We choose the simple formulation

\begin{equation} f^{\alpha}\left(c\right)=\varrho^{2}\left(c-c_{\alpha}\right)^{2} \end{equation}

\begin{equation} f^{\beta}\left(c\right)=\varrho^{2}\left(c_{\beta}-c\right)^{2} \end{equation}

\begin{equation} h\left(\eta_{1},...\eta_{p}\right)=\sum_{i=1}^{p}\eta_{i}^{3}\left(6\eta_{i}^{2}-15\eta_{i}+10\right) \end{equation}

\begin{equation} g\left(\eta_{1},...\eta_{p}\right)=\sum_{i=1}^{p}\left[\eta_{i}^{2}\left(1-\eta_{i}\right)^{2}\right]+\alpha\sum_{i=1}^{p}\sum_{j\neq i}^{p}\eta_{i}^{2}\eta_{j}^{2}, \end{equation}

where $f^{\alpha}$ and $f^{\beta}$ have minima at $c_{\alpha}$ and $c_{\beta}$, $\varrho^{2}$ controls the curvature of the free energies, and $\alpha$ controls the energy penalty incurred by the overlap of multiple non-zero $\eta_{i}$ values at the same point. Because the energy values of the minima are the same (Fig. 1), $c_{\alpha}$ and $c_{\beta}$ correspond exactly with the equilibrium atomic fractions of the $\alpha$ and $\beta$ phases.

The time evolution of $c$ is again governed by the Cahn-Hilliard equation, see CAHN and ELLIOT,

\begin{equation} \frac{\partial c}{\partial t}=\nabla\cdot\Bigg\{M\nabla\left(\frac{\partial f_{chem}}{\partial c}-\kappa_{c}\nabla^{2}c\right)\Bigg\}. \end{equation}

The Allen-Cahn equation, see ALLEN, which is based on gradient flow, governs the evolution of $\eta_{i}$,

\begin{equation} \frac{\partial\eta_{i}}{\partial t}=-L\left[\frac{\delta F}{\delta \eta_{i}}\right]=-L\left(\frac{\partial f_{chem}}{\partial\eta_{i}}-\kappa_{\eta}\nabla^{2}\eta_{i}\right), \end{equation}

where $L$ is the kinetic coefficient of $\eta_{i}$. We choose $M=5$ and $L=5$ so that the transformation is diffusion-controlled and the kinetic coefficients and gradient energy coefficients are isotropic. In addition, we again choose $c_{\alpha}=0.3$ and $c_{\beta}=0.7$, and further specify $\kappa_{c}=\kappa_{\eta}=3$, $\varrho=\sqrt[]{2}$, $w=1$, and $\alpha=5$. For these values, the diffuse interface between $0.1<\eta<0.9$ has a width of of 4.2 units.

Figure 1: Free energy density surface¶

Parameter values¶

Domain geometries and boundary conditions¶

Several boundary conditions, initial conditions and computational domain geometries are used to challenge different aspects of the numerical solver implementation. We test four combinations that are increasingly difficult to solve: two with square computational domains, see (a) and (b), with side lengths of 200 units, one with a computational domain in the shape of a "T," see (c), with a total height of 120 units, a total width of 100 units, and horizontal and vertical section widths of 20 units, and one in which the computational domain is the surface of a sphere with a radius of r = 100 units, see (d). While most codes readily handle rectilinear domains, a spherical domain may pose problems, such as having the solution restricted to a two-dimensional curved surface. The coordinate systems and origins are given in Fig. 2. Periodic boundary conditions are applied to one square domain, see (a), while no-flux boundaries are applied to the other square domain, see (b), and the "T"-shaped domain, see (c). Periodic boundary conditions are commonly used with rectangular or rectangular prism domains to simulate an infinite material, while no-flux boundary conditions may be used to simulate an isolated piece of material or a mirror plane. As the computational domain is compact for the spherical surface, no boundary conditions are specified for it. Note that the same initial conditions are used for the square computational domains with no-flux, see (b), and periodic boundary conditions, see (a), such that when periodic boundary conditions are applied, there is a discontinuity in the initial condition at the domain boundaries.

(a) Square periodic¶

A 2D square domain with $L_x = L_y = 200$ and periodic boundary conditions.

(b) Square no-flux¶

A 2D square domain with $L_x = L_y = 200$ and no flux boundary conditions.

(c) T-shape¶

A T-shaped region with zero flux boundary conditions and with dimensions, $a=b=100$ and $c=d=20$.

(d) Sphere¶

The domain is the surface of a sphere with radius 100.

Initial conditions¶

The initial conditions are chosen such that the average value of $c$ over the computational domain is approximately $0.5$.

Initial conditions for (a), (b) and (c)¶

The initial value for $c$ on the square and "T" computational domains is specified by $$ c\left(x,y\right) = c_{0}+\epsilon\left[\cos\left(0.105x\right)\cos\left(0.11y\right)+\left[\cos\left(0.13x\right)\cos\left(0.087y\right)\right]^{2}\right.\nonumber \\ \left.+\cos\left(0.025x-0.15y\right)\cos\left(0.07x-0.02y\right)\right], $$ where $c_{0}=0.5$ and $\epsilon=0.05$. The intial value for $\eta_i$ on the square and "T" domain are given by,

\begin{align*} \eta_{i}\left(x,y\right) = & \epsilon_{\eta}\left\{ \cos\left(\left(0.01i\right)x-4\right)\cos\left(\left(0.007+0.01i\right)y\right)\right.\nonumber \\ & +\cos\left(\left(0.11+0.01i\right)x\right)\cos\left(\left(0.11+0.01i\right)y\right)\nonumber \\ & +\psi\left[\cos\left(\left(0.046+0.001i\right)x+\left(0.0405+0.001i\right)y\right)\right.\nonumber \\ & \left.\left.\cos\left(\left(0.031+0.001i\right)x-\left(0.004+0.001i\right)y\right)\right]^{2}\right\} ^{2} \end{align*}

where $\epsilon_{\eta}$=0.1 and $\psi$=1.5.

Figure 2: initial $c$ for (a), (b) and (c)¶

Figure 3: initial $\eta_1$ for (a), (b) and (c)¶

Initial conditions for (d)¶

The initial value of $c$ for the spherical computational domain is specified by $$ c\left(\theta,\phi\right) = c_{0}+\epsilon_{sphere}\left[\cos\left(8\theta\right)\cos\left(15\phi\right)+\left(\cos\left(12\theta\right)\cos\left(10\phi\right)\right)^{2}\right.\nonumber \\ +\left.\cos\left(2.5\theta-1.5\phi\right)\cos\left(7\theta-2\phi\right)\right], $$ where $\epsilon_{\text{sphere}}=0.05$ and $c_0=0.5$, and $\theta$ and $\phi$ are the polar and azimuthal angles, respectively, in a spherical coordinate system. These angles are translated into a Cartesian system as $\theta=\cos^{-1}\left(z/r\right)$ and $\phi=\tan^{-1}\left(y/x\right)$ dependent upon angle.

The intial value for $\eta_i$ for the spherical domain is specified by

\begin{eqnarray} \eta_{i}\left(\theta,\phi\right) & = & \epsilon_{\eta}^{sphere} \left\{ \cos\left(i\theta-4\right)\cos\left(\left(0.7+i\right)\phi\right)\right.\nonumber \\ & & +\cos\left(\left(11+i\right)\theta\right)\cos\left(\left(11+i\right)\phi\right)\nonumber \\ & & +\psi\left[\cos\left(\left(4.6+0.1i\right)\theta+\left(4.05+0.1i\right)\phi\right)\right.\nonumber \\ & & \left.\left.\cos\left(\left(3.1+0.1i\right)\theta-\left(0.4+0.1i\right)\phi\right)\right]^{2} \right\} ^{2} \end{eqnarray}

where $\epsilon_{\eta}^{sphere}=0.1$ and $i=1,\ldots,4$ enumerates the order parameters corresponding to the different phase variants.

Figure 4: initial $c$ for (d)¶

Figure 5: initial $\eta_1$ for (d)¶

Results¶

Results from this benchmark problem are displayed on the simulation result page for different codes.

Feedback¶

Feedback on this benchmark problem is appreciated. If you have questions, comments, or seek clarification, please contact the CHiMaD phase field community through the Gitter chat channel or by email. If you found an error, please file an issue on GitHub.