Problem Setup

Larry Aagesen (@laagesen), David Montiel, Sourabh Kadambi (@sourabhkadambi), Sudipta Biswas (@SudiptaBiswas)

After the phase-field model formulation is developed, implemented in code, and verified, it can be set up to solve the scientific/engineering problem of interest. The purpose of this page is to give guidance on some of the important considerations when setting up your code to solve your specific problem, organized in the following sections:

Spatial Dimension (1D vs. 2D vs. 3D)

Boundary conditions

Interface width

Convergence

Impact of orientation

Kinetics and how long to run

Spatial Dimension

In setting up a phase-field simulation, it is important to understand the role of spatial dimensions in the physics governing the problem, and how that might impact the simulation results and its interpretation. This is particularly important if simulations are being performed in reduced dimensions compared to the actual dimension of the problem, which is a common practice employed to reduce the computational costs.

It is likely that the physical behavior of a certain phenomenon might scale differently in 1D, 2D and 3D. This could be related to the physical processes in the bulk regions of the system or at the interfacial regions or both. For instance, the amount of interfacial region relative to the bulk phase region differs in 1D, 2D and 3D, and also governs many physical aspects of material behavior.

Therefore, it is important to determine how reduced dimensionality might affect predictions or comparison of simulation results with reality. In certain situations, on the other hand, it might be appropriate to setup the problem in reduced dimensions of 1D or 2D due to the inherent symmetry or directionality in the problem. Moreover, for simple microstructure geometries, it might be possible to setup the model in cylindrical or spherical coordinates, which can allow the simplicity of 1D computations while also capturing 2D and 3D behavior of the model accurately.

Example

To demonstrate a case where simulation outcomes can significantly differ in different dimensions, we consider the single seed case of the homogeneous nucleation benchmark problem. A simple phase-field model with a single non-conserved order parameter describes an isothermal pure substance with one liquid phase (order parameter = $0$) and one solid phase (order parameter = $1$). The nucleation driving force for the solid phase is $\Delta f=\sqrt2/30$. The free energy of the diffuse interface is $\gamma = 1/3\sqrt2$.

We consider a seed nucleus of $7.5$ unit radius in 2D and 3D. The 2D simulation domain is of size $100\times100$, and the 3D simulation domain is of size $40 \times 40 \times 40$. The nucleus is centered at $x=y=z=0$. Uniform mesh of element size $\Delta x = 0.09765$ is used. Neumann boundary conditions, given by the zero normal-derivative of the order parameter, are applied on all domain boundaries. The figure below shows the starting geometry of the nucleus in 2D (left) and 3D (right).

The time evolution of the order parameter, given by the Allen-Cahn equation, was solved using the MOOSE framework. Details of the numerical method can be found in Wu et al. (2021). The simulations results of nucleus evolution are shown in the figures below: (left) order parameter profile measured radially from the center of domain at different evolution times; (right) radius as a function of evolution time. We see that the seed of same starting radius $r_\circ = 7.5$ units evolves starkly differently in 2D and 3D. While the nucleus grows in 2D, it shrinks and dissolves in 3D.

The role of dimensionality in homogeneous nucleation can be understood from the classical nucleation theory where the solid-liquid interface is modeled as a mathematically sharp interface. The interface is a 1D line in a 2D system and a 2D plane in a 3D system. The free energy of the nucleus particle $\Delta G(r)$ is a balance between the energy cost in forming the interface, and the energy released due to the driving force in forming the bulk particle. For a 2D system, $\Delta G(r) = 2\pi r\gamma - \pi r^2 \Delta f$. For a 3D system, $\Delta G(r) = 4 \pi r^2 \gamma - (4 \pi /3) r^3 \Delta f$. When the rate of change of free energy with respect to particle size is negative, the particle is favored to grow; otherwise it will shrink and dissolve. By setting $d \Delta G / d r = 0$, we obtain the critical radius as $r_c = \gamma/ \Delta f$ in 2D and $r_c = 2\gamma/ \Delta f$ in 3D.

We can now apply the above sharp-interface analysis to our diffuse interface approximation in the phase-field model. For the given model parameters, we obtain $r_c = 5$ units in 2D and $r_c = 10$ units in 3D. In our simulation setup of $r_\circ = 7.5$ units, $r_\circ$ $>$ $r_c$ in 2D, but $r_\circ < r_c$ in 3D. Therefore, the nucleus is favored to grow in 2D but shrinks in 3D as observed in the simulations.

The dependence on dimensionality is further illustrated by considering cases where the initial radius is close to the critical radius: $r_\circ = 0.99 r_c$, $r_\circ = r_c$ and $r_\circ = 1.01 r_c$. The simulation results of radius evolution are shown above for 2D (left) and 3D (right). As expected from the the classical homogeneous nucleation theory, the sub-critical nucleus with $r_\circ = 0.99 r_c$ shrinks and the super-critical nucleus with $r_\circ = 1.01 r_c$ grows.

Since the nucleus in the phase-field model is a diffuse-interface approximation of the classical sharp interface nucleus, $r_\circ = r_c$ is fairly close to an unstable equilibrium. Ideally, the radius would remain constant with time. However, since the system is unstable, small numerical errors accumulate with time, eventually leading to growth or shrinkage of the nucleus. During initial time steps, the interface profile is expected to undergo some changes from the starting profile due to equilibration. While the $tanh$ function is a common choice for the initial condition of the interface, it is an exact solution only for a planar interface, representative of a 1D scenario.

Boundary Conditions

Try different boundary conditions and check their impact

We used the results from Benchmark Problem 1a and Problem 1b to analyze the effect of boundary conditions. In addition, we solved for Cahn-Hilliard dynamics under the same initial condition and simulation parameters as problem 1 but using mixed boundary conditions, i.e., different boundary conditions for each boundary. We compare the results for each case at simulation time t = 1000. The simulations were carried out in the PRISMS-PF framework using a uniform square mesh with $N_x = N_y = 128$ linear elements and a time step, $\Delta t$ = 0.005. The results are shown in the figure below.

As can be seen in the figure above, for periodic boundary boundaries the $\alpha$ - $\beta$ domains are continuous on opposite sides of the system. For No-flux boundaries, the $\alpha$ - $\beta$ interfaces are normal to the boundary. For Dirichlet boundaries (bottom boundary of the right panel), the value of c is fixed along the boundary.

Interface width

Phase-field modeling is a diffuse interface approach, meaning that interfaces are represented by a smooth variation of one or more order parameters across the interface. The width of the interface is a function of the model parameters and the interface width is determined by model parameter choices. In some cases, an analytical expression is available that relates interface width to phase-field model parameters (such as free energy barrier height and gradient energy coefficient, which also impact the interfacial energy). In other cases, no analytical solution is available and the interface width must be approximated or determined numerically based on parameter choices. Such details are specific to the formulation being used.

Once the relationship between model parameters and interface width is understood, an appropriate selection of interface width needs to be made (while maintaining the correct interfacial energy for the system being studied). In some cases, the interface width in the phase-field model can be chosen to match the actual physical width of the interface being studied (such typically sub-nanometer). However, resolving a physically realistic interface width requires a sub-nanometer grid/mesh (grid/mesh convergence is described further in the next section). Using a physically realistic interface width often makes it computationally unfeasible to simulate systems large enough to be statistically representative. Therefore, in many cases an interfacial width that is much larger than the the physical width of an interface is used. In these cases, a careful balance between computational efficiency and model accuracy is needed; the interface width should be chosen to be small enough that the physics of the system are accurately represented, while maintaining adequate computational performance. A useful rule of thumb as a starting point is that the interface width should be at least an order of magnitude smaller than the smallest microstructural feature size of interest. Starting from this guideline, simulations of microstructural evolution with varying interface widths can be run to ensure that the choice of interface width does not affect the simulation results. A small test problem may be useful for testing convergence with respect to interface width; for example a shrinking circular grain embedded in another grain may be used for testing convergence of a grain growth model with respect to interface width, rather starting with large, costly simulations of hundreds of grains.

Convergence studies

Carry out grid/mesh convergence study AFTER you have finalized your interfacial width.

As an example of a mesh convergence study, we can consider Benchmark Problem 1 from the Phase-Field Community Hub. In this problem, which models spinodal decomposition using the Cahn-Hilliard equation, the width of the diffuse interface is 4.47, as defined by the Cahn-Hilliard equation and physical parameters in the problem statement. Given this interface width, we need to ensure there are a sufficient number of grid points (for finite difference schemes) or mesh elements (for finite element or finite volume schemes) through the diffuse interface to adequately resolve the variation of the composition order parameter.

In Problem 1b, a square domain with dimensions $200 \times 200$ and no-flux boundary conditions is considered. An example of a mesh convergence study for this problem using the MOOSE framework phase-field module (a finite element code) is described. The $200 \times 200$ domain is discretized using increasing numbers of elements $N_x$ in the $x$ direction and $N_y$ in the $y$ direction, maintaining $N_x = N_y$ for square elements (using linear Lagrange shape functions). The simulations with varying numbers of elements are run until a significant amount of microstructural evolution has occurred. We need to increase the number of elements until the microstructure at a fixed time no longer changes with further increases in the number of elements; at this point the simulation is converged with respect to the mesh resolution.

The simulation initial conditions and the microstructures at $t = 10,000$ are shown in the figure below.

As the number of elements in each direction is increased from 40 to 80 to 160, changes in the microstructure are observed. However, once the number of elements increases to 320, no further changes are observed in the microstructure. Therefore, the problem is converged with respect to mesh resolution at $N_x = N_y = 160$. For this number of elements, each element has size $\Delta x = \Delta y$ = 200 / 160 = 1.25. Therefore, the number of elements through the diffuse interface width is 4.47 / 1.25 = 3.6. Practical experience in the phase-field community has shown that somewhere between 3 to 5 elements through the interface are usually required to obtain mesh convergence; however, the appropriate resolution is problem-specific, and convergence should be tested for the specific physics and parameters at hand.

Carry out time step convergence studies. Try higher-order schemes, adaptive time stepping, etc. For explicit time integration, know your stability limit.

For codes that use explicit time integration, there is a maximum value of the time step beyond which the solution becomes numerically unstable. This stability limit can be determined by the Courant–Friedrichs–Lewy (CFL) condition and, in general, it strongly depends on the order of the spatial derivatives. Below we show the simulation results of Benchmark Problem 1b using time step values slightly below and slightly above the stability limit. We employed a forward Euler time-integration scheme and spatially discretized the system using $128 \times 128$ first order elements. The simulations were performed in the PRISMS-PF framework. The left panel of the figure below shows a snapshot of the concentration at time $t=33$, which was obtained using a stable time step of $\Delta t=0.015$. If the time step is increased to $\Delta t=0.016$, the time step goes above the stability limit. The right panel shows the concentration at time $t \simeq 32.39$ for $\Delta t=0.016$ and features a numerical instability that appears near the bottom boundary. After only four time steps the amplification of this instability caused the simulation to fail.

Codes that use implicit time stepping schemes may have fewer restrictions with respect to stability as time step size is increased, depending on the problem. However, discretization error still occurs in implicit schemes and increases with the size of the time step taken. Therefore, a convergence study should be carried out to ensure that the size of the time step does not affect the simulation results. As an example, we can again consider Benchmark Problem 1 from the Phase-Field Community Hub. Problem 1b was solved with the phase-field module from the MOOSE framework using $N_x = N_y = 160$. The results are shown below. The microstructure at $t = 1000$ remains the same for $\Delta t = 0.1$, $\Delta t = 0.25$, and $\Delta t = 0.5$. Some differences in the microstructure are observed at $\Delta t = 1$, and further differences become apparent at $\Delta t = 4$. In practice, it is recommended to start with a given time step and gradually decrease it until no changes are apparent with further decreases in time step size.

Adaptive time stepping can be useful to increase the time step size during time periods in the simulation where there are fewer microstructural changes, particularly for codes that use implicit time stepping schemes. However, convergence must still be checked for the parameters of the time stepping scheme being used. An example is the IterationAdaptiveDt time stepping scheme used in the MOOSE framework. This scheme attempts to increase or decrease the time step to keep the solver using a certain number of nonlinear iterations (controlled by the parameter optimal_iterations), within a window or plus or minus the parameter iteration_window. Higher values of optimal_iterations parameter result in higher time steps, and with that comes the risk of discretization error changing simulation results. As shown in the following, for Problem 1b, optimal_iterations values of 6 or 8 gave results consistent with the converged time step of 0.5, but when optimal_iterations was significantly increased to 15, changes in the microstructure resulted.

Impact of Orientation

Oftentimes, the results of the phase field simulations are sensitive to the orientation and alignment of the key microstructural features with the mesh/grid points. To demonstrate this we pick a simple solidification problem with dendritic structure formation. In this case, we utilize the solidification example from the MOOSE-based phase field module. This example problem can quickly demonstrate dendritic structure formation including the formation of secondary dendritic arms in a computationally cost-effective manner. Here, we use 4-fold symmetry of the structure and vary the reference angles to misorient the dendritic arms with respect to the mesh. Dendritic structures corresponding to 0 and 45 degree reference angle is presented below:

It is noteworthy that the shape of the dendrite varies with orientation (as observed by the difference in the dendrite center). For better comparison, we rotate the $45^o$ dendrite to align the primary dendrite arms with the reference orientation dendrite:

This highlights the slight differences between the dendrite shapes, especially the center and the secondary dendritic arms. Additionally, the growth rate of the solid also varies with orientation as observed by the change in solid area fraction over time:

Thus, it is important to evaluate the effect of orientation on the results by misaligning the grids. Furthermore, these effects are influenced by the discretization of the mesh. Hence, it is important to ensure that the mesh is refined enough to properly resolve the interfaces (if necessary, run a mesh convergence study). For more information about strategies to simulate multiple dendrites with varying orientation, please refer to the following works by Biswas et al. [1], Warren et al. [3], Dorr et al. [4], and Pusztai et al. [5].

Kinetics and how long to run

The amount of time the simulation should be run depends on the science or engineering question to be answered, and on the system being studied. For many classic phase-field problems such as grain growth and coarsening, a characteristic feature size of the system increases with time, and the progress of microstructural evolution slows as the characteristic feature size increases. For example, in grain growth, at long times the mean grain diameter $D$ increases with time as $D \propto t^{1/2}$. Another example is coarsening of spherical particles in a matrix; for sufficiently low volume fraction of particles in a matrix where solute is transported by diffusion, the mean particle radius $R$ increases with time as $R \propto t^{1/3}$.

It is also useful to be aware that the system may reach a metastable or stable equilibrium state with respect to system energy, in which case no further microstructural evolution will occur for increasing simulation time. In the grain growth example, if the system evolves to a single grain, stable equilbrium has been reached. If the system evolves to a two-grain configuration with a flat grain boundary between the grains, it has reached a metastable equilibrium; the system could still lower it energy by removing the grain boundary, but in that configuration, there is no kinetic driving force to remove it from the metastable state. To monitor for such possibilities, it is useful for the simulation to periodically output the total free energy of the system; a stable or metastable equilibrium state is indicated by a constant free energy with respect to time.

References

[1] Biswas et al., Solidification and grain formation in alloys: a 2D application of the grand-potential-based phase-field approach, Modelling and Simulation in Materials Science and Engineering, Volume 30, Number 2, 025013 (2022). https://doi.org/10.1088/1361-651X/ac46dc