Model Formulation

Nana Ofori-Opuku, Jim Warren, Pierre-Clement Simon

• Review the literature
  • Start with a careful lit review and see what others have done.
  • Take the time to understand the derivations from the literature.
  • Consider good examples from the literature (Can we list these on PFHub?)

General Considerations on the formulation of phase field models

Phase field models are, quite generally, extensions of classical non-equilibrium thermodynamics. There are quite a few treatments of this in the literature. Here we will follow formulations similar to those developed by Sekerka and Bi in "Interfaces for the Twenty-First Century." In that spirit, we will start with a very general formulation, including hydrodynamics, and then simplify the problem. For those interested in starting with a simple model, you can skip over much of the initial formulation and jump to the section on the further reduction of the problem.

In thermodynamics one usually starts from consideration of the internal energy $E$, where $E=E(S,M_i,V)$, and $S$ is the entropy, $V$ is the volume and $M_i$ is the mass of species $i$. If we decide we don't want to consider elasticity we can divide through by the volume and compute the bulk energy density $e=e(s,\rho_i)$, with $\rho_i=M_i/V$, and similarly $s=S/V$. It is then usually postulated that in addition to the thermodynamic variables there are any number of phase fields $\phi_j$ and the internal energy density is written $$e=e(s,\rho_i,\phi_j)$$.

Here we are going to consider a single phase field $\phi$ that distinguishes between two phases (could be liquid-solid, could be solid-gas, etc.) and postulate that the total entropy of a system of volume $V$ can be written as $$S = \int s^{NC} dV,$$ where $$s^{NC}=s-\frac{\epsilon^2}{2}|\nabla\phi|^2-\sum_i \frac{\alpha_i^2}{2}|\nabla\rho_i|^2.$$ Here the superscript $NC$ is meant to indicate that this is a non-classical definition of the entropy density $s$. We have chosen the simplest possible forms of the non-classical extensions (that is, simple square-gradient penalties). The entropy density is itself a function, of course, $$s=s(e,\rho_i,\phi).$$

To close the deal we now need to work our way through the various continuum-versions of the laws of nature.

Evolution equations

In order to write down evolution equations we need to resort to the know laws of physics. For continuum theories, we really have only a few ideas to fall back on

• Conservation of mass
• Conservation of energy (First Law of Thermodynamics)
• Conservation of momentum and maybe angular momentum if absolutely necessary
• The second law of thermodynamics (entropy increases, free energy decreases...)

A complete discussion of how to use the above rules in this context is outside the scope of this best-practice guide, but we offer a highly abbreviated discussion of the ''flavor,'' following the ideas of irreversible thermodynamics (we largely are following works like deGroot and Mazur, although the continuum mechanics community may be more comfortable with Noll, Coleman, and Truesdale), as well as the aforementioned work by Sekerka and Bi.

Mass

The first idea is that we need to specify a flux of particles $${\bf J}_i=\rho_i \left({\bf v}_i-{\bf v}\right),$$ where ${\bf v}_i$ is the velocity of the particles of species $i$. In general we want to work in a center-of-mass frame $${\bf v}=\frac{\sum_i \rho_i {\bf v}_i}{\sum_i \rho_i}.$$ Notice that this implies (by construction) that $\sum_i {\bf J_i}=0$.

The law of conservation of mass can be simply written as $$\frac{\partial \rho_i}{\partial t}+\nabla\cdot(\rho_i {\bf v}_i)=0.$$ Then we use the definition of ${\bf J}_i$ to rewrite that as $$\frac{\partial \rho_i}{\partial t}+ \nabla\cdot\left({\bf J_i}+\rho_i{\bf v}\right)=0.$$ Note that if we sum this equation over $i$, we get the continuity equation $$\frac{\partial \rho}{\partial t}+ \nabla\cdot\left(\rho{\bf v}\right)=0.$$

At this point is is common to introduce the notation $$\frac{D}{Dt}=\frac{\partial}{\partial t}+{\bf v}\cdot\nabla,$$ the "Lagrangian derivative" or "material derivative". The we can write the mass conservation equation as $$\frac{D \rho_i}{Dt}+ \nabla\cdot{\bf J_i}+\rho_i\nabla\cdot{\bf v}=0.$$

If the system in incompressible (a common assumption, although we'd have to abandon that if we want to consider gasses), then $\nabla\cdot{\bf v}=0$ and we can write $$\frac{D \rho_i}{Dt}+ \nabla\cdot{\bf J_i}=0.$$

Momentum

The law of conservation of momentum (Newton's second Law) for a system without body forces (like gravity) can be written $$\frac{\partial\rho {\bf v}}{\partial t}+\nabla\cdot\left( \rho {\bf v}\otimes{\bf v}-\sigma\right)=0.$$ Where $\sigma$ is the stress tensor. Using continuity and incompressibility as assumptions we get $$\rho\frac{D{\bf v}}{Dt}=\nabla\cdot\sigma.$$

Energy

Our final conservation expression is for energy. Just like the equation for momentum and mass, we start with a continuity form, and assume a flux of internal energy ${\bf J_e}$, and distinguish between the total energy density $e_T$, the internal energy density $e$ and the kinetic energy density $e_k=e_T-e= \frac{1}{2}\rho v^2$ and write $$\frac{\partial e_T}{\partial t}+\nabla\cdot\left(e_T{\bf v}+{\bf J_e}\right)=0.$$ We can then use the momentum equation as well as the mass continuity equation (along with incompressibility) to arrive at $$\frac{D e}{D t}+\nabla\cdot{\bf J_e}=\nabla{\bf v}:\sigma.$$

First Law of Thermodynamics

We now write down the laws of thermodynamics for the densities. $$e=Ts-p+\sum\rho_i\mu_i,$$ where $p$ is the hydrostatic pressure, $$p=-\frac{1}{3}{\mathrm Tr}\sigma,$$ and $\mu_i$ and $T$ are defined as derivatives of $e$ through the expression $$de=Tds+\sum\mu_i d\rho_i+\frac{\partial e}{\partial\phi}d\phi.$$ Of course, given where we're headed, let's rewrite this as $$ds=\frac{1}{T}de-\sum\frac{\mu_i}{T}d\rho_i+\frac{\partial s}{\partial \phi}d\phi.$$

Second Law

We have chosen the entropy to "extend" non-classically, as it is the star actor in the second law of thermodynamics, so we can write the entropy production (which must be positive) as $$S^{\mathrm prod}=\frac{dS}{dt}+\int {\bf J_s}\cdot{\bf dA}\ge0,$$ where ${\bf J_s}$ is the flux of entropy. We can compute (for an incompressible system) $$S^{\mathrm prod}=\int \left[\sum -\left(\frac{\mu_i}{T}\right)^{NC}\frac{D\rho_i}{Dt}+S_\phi\frac{D\phi}{Dt}+\frac{1}{T}\frac{De}{Dt}\right]dV+\int \left[{\bf J_s}-\sum \alpha_i^2\nabla\rho_i\frac{D\rho_i}{Dt}-\epsilon^2\nabla\phi\frac{D\phi}{Dt}\right]\cdot{\bf dA},$$ where we have introduced $$-\left(\frac{\mu_i}{T}\right)^{NC}=\frac{\partial s}{\partial \rho_i}+\alpha_i\nabla^2\rho_i,$$ and $$S_\phi=\frac{\partial s}{\partial\phi}-\epsilon^2\nabla^2\phi.$$ At this point we can insert the equations of motion we got from our conservation laws and do a lot of algebra, which we leave as an amusing exercise for the reader. We eventually find that $$S^{\mathrm prod}=\int \left[-\sum\nabla\left(\frac{\mu_i}{T}\right)^{NC}\cdot{\bf J_i}+S_\phi\frac{D\phi}{Dt}+\nabla \frac{1}{T}\cdot{\bf J_e}+\frac{1}{T}Y:\nabla {\bf v}\right]dV,$$ and $${\bf J_s}=s^{NC}{\bf v} + \frac{1}{T}{\bf J_e}-\sum \left(\frac{\mu_i}{T}\right)^{NC}{\bf J_i}+\sum \alpha_i^2\frac{D\rho_i}{Dt}\nabla\rho_i+\epsilon^2\frac{D\phi}{Dt}\nabla\phi.$$ It is worth noting that this form for $J_s$ eliminates the explicit surface terms from the entropy production. We have also introduced the tensor $Y$ which is rather a complicated beast (it came from all those integrations by parts which we have skipped in this presentation), and has the form $$\frac{Y}{T}=\frac{\sigma}{T}+\sum \alpha_i^2\nabla\rho_i\otimes\nabla\rho_i+\epsilon^2\nabla\phi\otimes\nabla\phi+\left(\frac{p}{T}-\sum \alpha_i^2\rho_i\nabla^2\rho_i-\sum\frac{\alpha_i^2}{2}|\nabla\rho_i|^2-\frac{\epsilon^2}{2}|\nabla\phi|^2\right){\bf I},$$ where $\bf I$ is the identity tensor.

Constitutive Equations and Evolution Equations

After all this work we can now find constitutive equations. Inspection of the final expression for $S^{\mathrm prod}$, suggests a that we can guarantee it be positive definite (which is required from the second law) by assuming $$J_i= -M_i\nabla\left(\frac{\mu_i}{T}\right)^{NC},$$ $$J_e=M_T\nabla\frac{1}{T},$$ $$Y_{kl}=\eta\left(\frac{\partial v_k}{\partial x_l}+\frac{\partial v_l}{\partial x_k}\right),$$ $$\frac{D\phi}{Dt}=M_\phi S_\phi,$$ where we have introduced mobilities and the viscosity $\eta$. It should be emphasized that we could have assumed all sorts of "cross-couplings" between the various fluxes here, but instead we have made the simplest assumptions. Inserting these into our equations of motion, we (at last!) arrive at general evolution equations $$\frac{D e}{D t}=M_T\nabla^2\frac{1}{T}+\nabla{\bf v}:\sigma,$$ $$\frac{D \rho_i}{D t}=\nabla\cdot M_i\nabla\left(-\frac{\partial s}{\partial \rho_i}-\alpha_i\nabla^2\rho_i\right)$$ $$\frac{D\phi}{Dt}=M_\phi\left(\frac{\partial s}{\partial\phi}+\epsilon^2\nabla^2\phi\right).$$ $$\frac{\sigma}{T}=\frac{\eta}{T}\left(\frac{\partial v_k}{\partial x_l}+\frac{\partial v_l}{\partial x_k}\right)-\sum \alpha_i^2\nabla\rho_i\otimes\nabla\rho_i+\epsilon^2\nabla\phi\otimes\nabla\phi-\left(\frac{p}{T}-\sum \alpha_i^2\rho_i\nabla^2\rho_i-\sum\frac{\alpha_i^2}{2}|\nabla\rho_i|^2-\frac{\epsilon^2}{2}|\nabla\phi|^2\right){\bf I},$$ with $$\rho\frac{D{\bf v}}{Dt}=\nabla\cdot\sigma.$$

Further reduction of the problem

We have presented the derivation of equations of evolution for an incompressible, two-phase, multicomponent system. Most phase field treatments ignore the velocity terms, but for pedagogical purposes we have retained all of the terms that arise from flow (for an incompressible system). These terms should not just be tossed away without careful consideration of the system of interest! Nonetheless, now that we have taken care to show how these terms arise, and also, for careful readers, how to extend this approach to multiple phases and additional gradient corrections, we will proceed with some more simplifications for a less complex system. For those who are interested in solid state systems that can creep, the work of Mishin, Warren, Sekerka, and Boettinger (2013) extends this framework. Here we eliminate the ${\bf v}$ equations by fiat, assuming that only diffusion controls the evolution of the system, which is often reasonable in microgravity situations. We can also go further, and consider an isothermal system. Then we have $$\frac{\partial \rho_i}{\partial t}=\nabla\cdot M_i\nabla\left(-\frac{\partial s}{\partial \rho_i}-\alpha^2_i\nabla^2\rho_i\right),$$ $$\frac{\partial\phi}{\partial t}=M_\phi\left(\frac{\partial s}{\partial\phi}+\epsilon^2\nabla^2\phi\right).$$

This system of equations should be adequate to describe a diffusion-controlled multicomponent, isothermal, two-phase system. One could, of course, retain thermal diffusion, add more phases, retain the velocity terms and more!

Clearly define the equations, parameters, initial conditions, and boundary conditions

Now that we have specified that we wish to consider, a diffusion-controlled multicomponent, isothermal, two-phase system. We still need to decide on the specifics of the system, and, in particular, the thermodynamic state functions that detail how the energy or entropy vary. All of the variables need to be clearly defined, and ideally should be either parameters that can be traced to thermodynamic variables or other physical constraints. As noted, we wish to consider an isothermal system, and so we can use a Legendre transformation to replace $s$ in the above expressions with the Helmholtz free energy $f=e-Ts$. $$\frac{\partial \rho_i}{\partial t}=\nabla\cdot M_i\nabla\left(\frac{1}{T}\frac{\partial f}{\partial \rho_i}-\alpha^2_i\nabla^2\rho_i\right)$$ $$\frac{\partial\phi}{\partial t}=M_\phi\left(-\frac{1}{T}\frac{\partial f}{\partial\phi}+\epsilon^2\nabla^2\phi\right).$$

Additionally, we will now only consider only a two component systems, so we only have species $\rho_1$ and $\rho_2$. Note that the total density $\rho=\rho_1+\rho_2$ is a constant (for an incompressible system). We could formulate the problem in terms of these variables or we can move to a concentration $c_1=\rho_1/(\rho_1+\rho_2)$. Clearly $c_1+c_2=1$, so we only have one independent variable. In this case we can pick $c_1$ or $c_2$ as the independent variable and label it $c$. As the reader can see, these types of exercises are highly non-trivial, require substantial care, and need to be done with rigor to understand the host of approximations and assumptions that can be made (at a minimum here we have ignored elasticity, different molar volumes for the species, and other issues are external fields, and there will be more to come!). We move to a concentraion picture and redefine the mobilities, $\alpha$ and $\epsilon$ in such a way that the following expressions hold: $$\frac{\partial c}{\partial t}=\nabla\cdot M_i\nabla\left(\frac{\partial f}{\partial c}-\alpha^2\nabla^2c\right)$$ $$\frac{\partial\phi}{\partial t}=M_\phi\left(\epsilon^2\nabla^2\phi-\frac{\partial f}{\partial\phi}+\right).$$

#State function

Dimensional Analysis

Clearly define and pose your physical problem before solving

Once one has decided on the thermodynamic picture. It is best to get a bit more specific in posing your physical problem. Now lets say we want to consider a binary solution where there are just two phases (say liquid and solid). We could have considered three components $n_1, n_2, n_3$, or five phases $\phi_1 ...\phi_5$, but again, lets keep it simple. For a binary we have $n_1$ and $n_2$. We could formulate the problem in terms of these variable or we can move to concentrations by assuming a fixed molar volume for all of the components ($v_m=V/N$). In that case, we have that $n_1+n_2=N/V=1/v_m$, and then $c_1+c_2=1$, where $c_1=n_1*v_m$. In this case we can pick $c_1$ or $c_2$ as the independent variable and label it $c$. As the reader can see, these types of exercises are highly non-trivial, require substantial care, and need to be done with rigor to understand the host of approximations and assumptions that can be made (at a minimum here we have ignored elasticity, different molar volumes for the species, and other issues are external fields, and there will be more to come!).

We're now in a position to write down something that looks like a phase field model. We have a variable $c$ and for two phases we only need a single phase field variable $\phi$ that, in the usual tradition takes on the value $\phi=1$ in the solid, and $\phi=0$ in the liquid. $$w=w(c,\phi,\nabla c,\nabla \phi)$$, and we can go one step further and write $$w=e(s,c,\phi) +\frac{\epsilon^2}{2}\left|\nabla c\right|^2+\frac{\alpha^2}{2}\left|\nabla \phi\right|^2.$$

We have now made a few more major assumptions. We have postulated that the energy density can be written as a ''bulk'' piece that interpolates between two phases $\left(e(c,\phi)\right)$ and a pair of gradient terms that have an energy cost for creating interfaces. We could have made other assumptions, especially around the form of the gradient terms form, were we could have introduced anisotropy into the picture by making the cost of interfaces different in different directions. We choose to not do that here. We now are almost ready to write down our equations of evolution!

• Complete dimensional analysis and understanding relevant scales; state the reference frame
• Plot your free energy before implementing it
• Consider the impact of interpolation functions, barrier functions, etc.
• Consider every boundary term and make sure you aren’t eliminating important terms. Check the fluxes (add benchmark problem or discussion to Poisson problem?)
• (Show examples on PFHub, possibly benchmark problems?)
• Be mindful of your assumptions and approximations Having clearly defined the mathematical framework underlying the problem, it becomes important, prior to progressing further, to contemplate and thoroughly grasp the various assumptions and approximations inherent within the formulation. Often, the exacting mathematical underpinnings from which the model has arisen may have presented challenges of being intractable in their complexity, or alternatively, they might have been of an impractical nature or even lacking contact with empirical observations.

The considerations surrounding assumptions and approximations are distributed throughout various places throughout the model's formulation. Foremost, are the descriptors of thermodynamic energy density, which provide the basis for our understanding of energy relationships, as well as the delineations of boundary conditions and the specifications of initial conditions. Notably, the latter two warrant attention since, in the absence of these conditions, the partial differential equations that determine the temporal evolution of our problem could only hypothetically yield an infinite set of solutions.

• If possible, run a small test problem with and without approximations and to quantify their impact
• Consider how approximations change in different dimensions and symmetries
• If possible, avoid assuming symmetry to reduce the complexity of your problem

Carefully consider your length and time scales

a) Identify the Smallest Feature Size of Interest in Your Problem: The smallest feature size refers to the smallest spatial scale that needs to be accurately resolved in your simulation. This could be, for example, the width of an interface between two distinct materials or the size of a fine structure. Properly resolving the smallest features is essential for capturing detailed interactions and behaviors in your system.

Additional Item 1: Generally the smallest feature is the interface width between the features that need to be resolved. The interface width often represents the smallest detail that demands attention. Neglecting this width may lead to inaccurate results and incomplete understanding of the system's behavior.

b) Identify the Smallest Time Scale of Interest in Your Problem: The smallest time scale corresponds to the fastest processes occurring in your system. It's important to accurately capture these time scales to ensure that rapid events are not overlooked or misrepresented.

Additional Item 2: The time steps you consider shall be small enough to allow capturing the kinetics of the process you want to model. Choosing appropriate time steps is crucial to accurately capture dynamic processes. Smaller time steps are necessary to capture fast kinetics and prevent underestimation or distortion of time-dependent phenomena.

c) Identify the Smallest Length Scales of Interest for All Physics and Consider How They Interact: This involves analyzing the smallest spatial scales relevant to each physical process in your simulation. Understanding how these different length scales interact is essential for capturing multi-scale effects accurately.

Additional Item 3: Depending on the numerical method used for solving the equations, you may need 5-10 mesh elements to resolve the interface. The mesh resolution, i.e., the number of mesh elements used to discretize the domain, plays a crucial role in capturing interfaces accurately. Higher mesh resolution is required to resolve fine interfaces properly.

1. The Sample Size You Consider Will Be Relevant to Your Feature Size: The size of the simulation domain should be appropriately chosen based on the scale of the features you are interested in. Having a domain that is much larger than necessary can lead to unnecessary computational costs, while having a domain that is too small might exclude important interactions.

2. It Is Recommended to Formulate Your Problem in a Non-Dimensional Form If Possible: Non-dimensionalization involves scaling variables in your equations to remove physical units. This can help simplify the equations, reduce the number of parameters, and make the problem more amenable to analysis and comparison.

3. If You Consider an Initial Condition Far from the Expected Steady State Solution, You May Need to Select a Significantly Smaller Time Step: When the initial condition differs significantly from the expected steady state, smaller time steps may be necessary to accurately capture the transient behavior and prevent numerical instability.

4. Be Mindful of Initial Conditions You Select: Poorly chosen initial conditions, especially for interfaces or abrupt changes, can lead to numerical artifacts, excessive computational demands, and inaccurate results. It's important to select physically meaningful initial conditions that smoothly transition into the desired system state.

   Additional Item 4: They may also result in numerical artifacts that lead to artificial stabilization of unstable phases. Inaccurate initial conditions can stabilize phases that would otherwise be unstable, leading to unrealistic outcomes in your simulation.

5. Using Automatic Time Stepping Integrator Algorithms Can Benefit PF Simulations: Automatic time stepping algorithms adjust the time step size during simulation based on the system's dynamics. This can improve simulation efficiency by using larger time steps during stable periods and smaller steps during dynamic events.

In summary, these considerations highlight the importance of careful parameter selection, proper resolution of length and time scales, and accurate initial conditions when performing numerical simulations. Addressing these aspects helps ensure the reliability and validity of your simulation results in various scientific and engineering applications.

• Non-dimensionalization
  • It gets MUCH harder with coupled physics, so be careful