In a species separation phase field model, a single concentration is used to model a two phase system. The evolution of the concentration is predicted by taking derivatives of a double-well free energy. This free energy can be represented with a simple thermodynamic model such that $$$$f_0(c) = k_b T (c \ln c + (1 - c) \ln(1 - c)) + Wc(1-c)$$$$ However, the natural logs in this expression can significantly complicate the solution, typically requiring small time steps to converge.

An alternative to the thermodynamic free energy is to use polynomial free energy, which simplify the numerical solution but decrease the physical accuracy. We have implemented polynomial free energies of three different orders in the phase field module that can be used to solve this problem.

Polynomial Equations

We have implemented three polynomial free energies of different orders, as shown in the figure.

Plot of the three polynomial free energies, with the thermodynamic free energy shown for reference.

The simplest model is the fourth order free energy, where $$$$f^4_0(c) = 2.0^4 W (c - c_{eq})^2) (1 - c - c_{eq})^2,$$$$ Where $$$W$$$ is the height of the energy barrier and $$$c_{eq}$$$ is the equilibrium concentration.

The sixth and eighth order polynomials are defined by the equation $$$$f^N_0(c) = 2.0^N W \sum_{n=1}^N H^N_n c^n$$$$ where $$$$\begin{eqnarray} H^6_1 &=& \frac{3}{2} c_{eq} (c_{eq} - 1) \\ H^6_2 &=& - \frac{9}{2} c_{eq}^2 + \frac{9}{2} c_{eq} + \frac{3}{4} \\ H^6_3 &=& 6 c_{eq}^2 - 6 c_{eq} - \frac{7}{2} \\ H^6_4 &=& - 3.0 c_{eq}^2 + 3.0 c_{eq} + \frac{27}{4} \\ H^6_5 &=& -6\\ H^6_6 &=& 2 \end{eqnarray}$$$$


$$$$\begin{eqnarray} H^8_1 &=& \frac{3}{4} c_{eq}(c_{eq} - 1) \\ H^8_2 &=& - \frac{15}{4}c_{eq}^2 + \frac{5}{4} c_{eq} + \frac{3}{8} \\ H^8_3 &=& 10 c_{eq} - 10 c_{eq} - \frac{11}{4} \\ H^8_4 &=& - 15 c_{eq}^2 + 15 c_{eq} + \frac{75}{8} \\ H^8_5 &=& 12 c_{eq}^2 - 12 c_{eq} - 18 \\ H^8_6 &=& -4 c_{eq}^2 + 4 c_{eq} + 20 \\ H^8_7 &=& - 12 \\ H^8_8 &=& 3 \end{eqnarray}$$$$

Quantitative Model Development

To create a quantitative model using the polynomial free energies, values for the mobility $$$M$$$, the barrier height $$$W$$$, $$$\kappa$$$ and $$$c_{eq}$$$ must be determined in terms of measurable quantities. The equilibrium concentration $$$c_{eq}$$$ defines the locations of the wells in the free energy according to $$$$\begin{eqnarray} c_{eq} = e^{-\frac{E_f}{k_b T}} \end{eqnarray}$$$$ To determine the other parameters, we follow the derivation from Moelans et al 2008. To simplify the derivation, we consider an interface along the $$$y$$$-direction, such that the interfacial width is defined as $$$$\begin{eqnarray} l_i &=& \frac{1}{|d x/d c|_{x=0}} \\ &=& \sqrt{\frac{\kappa_N}{2\ \mathrm{max}({f^N_{loc}})}}. \end{eqnarray}$$$$ The polynomial free energies have been developed such that $$$\mathrm{max}(f^N_{loc}) = W$$$, thus $$$$\begin{eqnarray} l_i = \sqrt{\frac{\kappa}{2 W}}. \end{eqnarray}$$$$ Therefore, $$$\kappa_N$$$ can be expressed as a function of the barrier height and the interfacial width according to \begin{eqnarray} \kappa = 2 W l_i^2 .\label{eq:kappa} \end{eqnarray}

To define the barrier height $$$W$$$, we first determine the surface energy. The surface energy is calculated according to $$$$\begin{eqnarray} \sigma = \int^{+\infty}_{-\infty} \left( f^N_{loc} + \frac{\kappa}{2} \left( \frac{\partial c}{\partial x} \right)^2 \right) dx.\label{eq:gbenergy1} \end{eqnarray}$$$$ According to the principles of variational calculus, the functions that extremize the function have zero first derivatives, i.e. $$$$\begin{eqnarray} \frac{\partial f_{loc}}{\partial c} - \kappa \frac{\partial^2 c}{\partial x^2} = 0, \end{eqnarray}$$$$ or, the integrated equation $$$$\begin{eqnarray} f_{loc} - \frac{\kappa}{2} \left( \frac{\partial c}{\partial x} \right)^2 = 0. \end{eqnarray}$$$$ By rearranging the equation, we obtain $$$$\begin{eqnarray} \frac{\partial c}{\partial x} = \sqrt{ \frac{2 f^N_{loc}}{\kappa}}. \label{eq:dcdx} \end{eqnarray}$$$$ Combining these equations gives $$$$\begin{eqnarray} \sigma &=& \int^{+\infty}_{-\infty} \left( f^N_{loc} + \frac{\kappa}{2} \frac{2 f^N_{loc}}{\kappa} \right) dx \\ &=& 2 \int^{+\infty}_{-\infty} f_{loc}(c) \ dx. \end{eqnarray}$$$$ After a change of variables and substitution, $$$$\begin{eqnarray} \sigma &=& 2 \int^{1}_{0} f_{loc}(c) \frac{d x}{d c} dc \\ &=& \sqrt{2 \kappa} \int^{1}_{0} \sqrt{f_{loc}(c)}\ dc. \label{eq:gbenergy2} \end{eqnarray}$$$$ To define the value of the integral $$$\int^{1}_0 \sqrt{f_{loc}(c)} \ dc$$$, we express it as $$$$\begin{eqnarray} \int^{1}_{0} \sqrt{f_{loc}(c)}\ dc &=& \sqrt{W} \ 2^{\frac{N}{2}} \int^{1}_{0} \sqrt{\sum_{n=0}^N H^N_n c^{\frac{n}{2}}} dc \\ &=& \sqrt{W} \ K_N \label{eq:int} \end{eqnarray}$$$$ where $$$K_N = 2^{\frac{N}{2}} \int^{1}_{0} \sqrt{\sum_{n=0}^N H^N_n c^{\frac{n}{2}}} dc$$$ is a function of $$$c_{eq}$$$ that varies for each polynomial order $$$N$$$. If we assume that $$$c_{eq} \ll 1$$$, $$$K_N$$$ becomes a constant value which we can calculate for each polynomial order, according to $$$$\begin{eqnarray} K_4 &=& \frac{2}{3} \\ K_6 &=& \frac{3}{16} \sqrt{3} - \frac{9}{64} \sqrt{2} \ln \left(\frac{\sqrt{3} - \sqrt{2}}{\sqrt{3} + \sqrt{2}} \right) \\ K_8 &=& 0.83551 \end{eqnarray}$$$$ with $$$K_8$$$ evaluated numerically. Thus, the surface energy is equal to $$$$\begin{eqnarray} \sigma = K_N \sqrt{2 \kappa W}. \label{eq:surfenergy} \end{eqnarray}$$$$ We can then solve for $$$W$$$ to obtain $$$$\begin{eqnarray} W = \frac{\sigma^2}{2 K_N^2 \kappa_N}. \end{eqnarray}$$$$ By combining these equations, we obtain $$$$\begin{eqnarray} \kappa = \frac{\sigma l_i}{ K_N}. \end{eqnarray}$$$$ This expression is then substituted to give $$$$\begin{eqnarray} W = \frac{\sigma}{2 K_N l_i}. \end{eqnarray}$$$$

The form of the mobility must be derived to produce a quantitative physical behavior. In order for the phase field model to accurately describe vacancy diffusion, the Cahn-Hilliard equation must be equivalent to Fickian diffusion when $$$c \ll 1$$$, i.e. $$$$\begin{eqnarray} \nabla \cdot M \nabla \frac{\partial F}{\partial c} = \nabla \cdot D \nabla c, \end{eqnarray}$$$$ where $$$M$$$ is the mobility for the polynomial free energy and $$$D$$$ is the diffusivity. Thus, for $$$c \ll 1$$$ $$$$\begin{eqnarray} M = D \frac{\nabla c}{\nabla {\partial f^N_{loc}}/{\partial c}} \label{eq:mob} \end{eqnarray}$$$$ For each of the polynomial free energies, the value of $$$\nabla \partial f_{loc}^N/\partial c$$$ for $$$c \ll 1$$$ can be written in the form $$$$\begin{eqnarray} \nabla \frac{\partial f_{loc}^N}{\partial c} = 2^{N+1} H^N_2 W \nabla c \label{eq:mob2} \end{eqnarray}$$$$ where $$$H^N_2$$$ is a polynomial coefficient defined above. Thus, the value for the mobility is determined by combining the equations to obtain $$$$\begin{eqnarray} M &=& \frac{D}{2^{N+1} C^N_2 W_N} \\ &=& \frac{K_N}{ 2^{N+1} C^N_2} \frac{D l_i}{ \sigma}\label{eq:Mi} \end{eqnarray}$$$$

To summarize, the quantitative expressions for the four model parameters required by the polynomial free energies can be obtained using the equations above. We restate them here: $$$$\begin{eqnarray} c_{eq} &=& e^{-\frac{E_f}{k_b T}} \\ \kappa &=& \frac{\sigma l_i}{K_N} \\ W &=& \frac{\sigma}{2 K_N l_i} \\ M &=& \frac{K_N}{ 2^{N+1} H^N_2} \frac{D l_i}{ \sigma} \end{eqnarray}$$$$

Numerical Implementation in the Phase Field Module

The polynomial free energies have been implemented in the PolynomialFreeEnergy material in the MOOSE phase field module. The free energies are all contained in the same material, where the order is an input parameter. To minimize required code and coding errors, the free energy derivatives are obtained using automatic differentiation via the ExpressionBuilder tools. The free energy is used by the parsed kernels, SplitCHParsed and CahnHilliard.

The material properties required for the model have been implemented in the PFParamsPolyFreeEnergy, which also takes the polynomial order as an input.