Benchmark Problem 4: Elastic Precipitate

See the journal publication entitled "Phase Field Benchmark Problems for Dendritic Growth and Linear Elasticity" for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems.

Overview

Precipitates are a key microstructural feature impacting the strength of alloys [1], and they are often elastically stressed, which affects their shape and their microstructure evolution during service. Elasticity has long been incorporated into phase field models: indeed, Cahn's seminal paper on spinodal decomposition [2] incorporates elastic strains due to composition fluctuations. Eshelby presents an analytical solution for the elastic field of a single coherent, elastically stressed precipitate in an infinite matrix [3], but the generalized problem of multiple interacting precipitates in a matrix with arbitrary crystal structure, lattice parameter misfit and elastic stiffnesses can only be solved numerically. Sharp-interface approaches provide insight into equilibrium elastic shapes and coarsening under the influence of elastic stress [4, 5, 6, 7, 8], but these approaches have difficulty simulating precipitate splitting or merging. Early phase field formulations studying elastically stressed precipitates demonstrate the power of the method (e.g., Refs. [9, 10, 11]), and present-day studies have expanded to 3D simulations (e.g., Refs. [12, 13, 14, 15]) and formulations that include plasticity (e.g., Refs. [16, 17, 18]).

Model Formulation

In this formulation, one phenomenological order parameter, $\eta$, is evolved, which has a value of 0 in the matrix and a value of 1 in the precipitate for an unstressed system with planar interfaces. This choice makes interpolation of materials properties between phases straightforward. The free energy of the system, $\mathcal{F}$, includes contributions from interfacial and elastic energy and is expressed as

\begin{equation} \mathcal{F}=\int_{V}\left(f_{\text{bulk}}\left(\eta\right) + \frac{\kappa}{2}|\nabla \eta|^{2} + f_{\text{el}}\left(\eta\right) \right)dV, \end{equation}

where $f_{\text{bulk}}$ is the bulk free energy density, $\kappa$ is the gradient energy coefficient, and $f_{\text{el}}$ is the local elastic free energy density. The $f_{\text{bulk}}$ term is a symmetric double-well with minima of zero, such that its contribution is only to the interfacial energy. As discussed in Ref. [1], we choose $f_{\text{bulk}}$ to have a 10th-order polynomial form,

\begin{equation} f_{\text{bulk}}=w\sum_{j=0}^{10}a_j\eta^j, \end{equation}

which makes the energy wells of the matrix and precipitate phases deep and narrow. This prevents the actual value of $\eta$ in each phase from shifting significantly from its equilibrium value due to the presence of a curved interface or elastic strain. The height of the energy barrier is controlled by $w$. The $f_{\text{bulk}}$ coefficients are given in the following table, and ensure that $f_{\text{bulk}}\left(0\right)=f_{\text{bulk}}\left(1\right)=f_{\text{bulk}}'\left(0\right)=f_{\text{bulk}}'\left(1\right)=0$ and that the energy curve remains concave down between the two energy wells.

Parameter	Value
$a_0$	0
$a_1$	0
$a_2$	8.072789087
$a_3$	-81.24549382
$a_4$	408.0297321
$a_5$	-1244.129167
$a_6$	2444.046270
$a_7$	-3120.635139
$a_8$	2506.663551
$a_9$	-1151.003178
$a_{10}$	230.2006355

The large number of significant digits are necessary to ensure that the first derivative of $f_{\text{bulk}}$ is zero at $\eta=0$ and $\eta=1$.

The elastic energy density is given as [1]

\begin{equation} f_{\text{el}}\left(\eta\right)=\frac{1}{2}\sigma_{ij} \epsilon_{ij}^{\text{el}}, \end{equation}

where $\sigma_{ij}=C_{ijkl}\left(\eta\right) \epsilon_{ij}^{\text{el}}$ is the elastic stress, $\epsilon_{ij}^{\text{el}}$ is the elastic strain, and $C_{ijkl}\left(\eta\right)$ is the elastic stiffness tensor such that the system is mechanically stable (the Einstein summation convention is used). To incorporate the dependence of the elastic stiffness on the phase, the stiffness is interpolated smoothly from one phase to the other across the diffuse interface,

\begin{equation} C_{ijkl}\left(\eta\right)= C_{ijkl}^{\text{matrix}}\left[1-h\left(\eta\right)\right]+C_{ijkl}^{\text{precip}} \, h\left(\eta\right), \end{equation}

where $C_{ijkl}^{\text{matrix}}$ and $C_{ijkl}^{\text{precip}}$ are the stiffness tensors of the matrix and precipitate phases, respectively, and $h\left(\eta\right)=\eta^3\left(6\eta^2-15\eta+10\right)$ is a smooth interpolation function that ensures that $h\left( 0 \right ) = h'\left( 0 \right)=h'\left( 1 \right)=0$ and $h\left(1\right) = 1$ [2].

Because the lattice parameters of the two phases are different, the elastic strain differs from the total strain, $\epsilon_{ij}^{\text{total}}$, as [3]

\begin{equation} \epsilon_{ij}^{\text{el}}=\epsilon_{ij}^{\text{total}}-\epsilon_{ij}^{0}\left(\eta\right), \end{equation}

where $\epsilon_{ij}^{0}$ is the local stress-free strain. It is calculated as

\begin{equation} \epsilon_{ij}^0\left(\eta\right)= \epsilon_{ij}^T \, h\left(\eta\right), \end{equation}

where $\epsilon_{ij}^{T}$ is the crystallographic misfit strain tensor between the matrix and precipitate phases defined with respect to the matrix. Finally, the total strain is related to the displacements, $u_i$, as [3]

\begin{equation} \epsilon_{ij}^{\text{total}}=\frac{1}{2}\left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right]. \end{equation}

In this problem, precipitate shapes must evolve to their equilibrium shape while remaining as small particles embedded in a much larger matrix. To do so, we employ the Cahn-Hilliard equation to perform fictive time evolution [2, 1], which conserves the total integral of $\eta$ within the simulation. The evolution of $\eta$ is given as

\begin{equation} \frac{\partial\eta}{\partial t}=\nabla\cdot\left[M\nabla\left\{ \frac{\delta \mathcal{F}}{\delta\eta}\right\} \right], \end{equation}

where $M$ is the mobility and the chemical potential is

\begin{equation} \mu \equiv \frac{\delta \mathcal{F}}{\delta\eta}=\frac{\partial f_{\text{chem}}}{\partial\eta}+\frac{\partial f_{\text{elastic}}}{\partial\eta}-\kappa\nabla^{2}\eta. \end{equation}

We have flexibility in choosing $M$, as we are only interested in the final state of the system. Furthermore, we assume that the relaxation dynamics for elasticity are much faster than for the diffusion of $\eta$, as is generally the case for phase field models. As such, we solve the time-independent equation for mechanical equilibrium at each time step,

\begin{equation} \nabla\cdot\sigma_{ij} = 0. \end{equation}

Parameterization and simulation conditions

This problem is solved in two dimensions to reduce computational costs, but note that we do not utilize symmetry to further reduce the problem size. The matrix and precipitate phases have cubic symmetry, such that three independent elastic stiffnesses exist for each phase: $C_{1111}$, $C_{1122}$, and $C_{1212}$ [1], and we take $C_{ijkl}^{\text{precip}}=1.1C_{ijkl}^{\text{matrix}}$. In addition, the precipitate misfit strain takes the form $\epsilon^T_{11}=\epsilon^T_{22} > 0$, $\epsilon^T_{12}=0$. Because this benchmark problem relies on the balance between interfacial and elastic energy, we use dimensional units of attojoules and nanometers. The diffuse interface width is chosen as 5 nm for $0.05 < \eta < 0.95$ and the interfacial energy is chosen as 50 aJ/nm$^2$ (equivalent to 50 mJ/m$^2$). The model parameters are given in the following table.

Parameter values for all variants

Quantity	Symbol	Value
Gradient energy coefficient	$\kappa$	0.29 aJ/nm
Well height	$w$	0.1 aJ/nm$^3$
Mobility	$M$	5
Misfit strain	$\epsilon^T_{11}$=$\epsilon^T_{22}$	0.5 %
Elastic stiffness matrix	$C^{\text{matrix}}_{1111}$	250 aJ/nm$^3$
Elastic stiffness matrix	$C^{\text{matrix}}_{1122}$	150 aJ/nm$^3$
Elastic stiffness matrix	$C^{\text{matrix}}_{1212}$	100 aJ/nm$^3$

Note that 1 aJ/nm$^3$ is equivalent to 1 GPa.

We utilize both circular and elliptical initial precipitate shapes for a given initial precipitate area [2, 3]; all initial precipitates have a diffuse interface width of 5 nm. To have an equal area for an ellipse as a circle with radius $r$, we choose ellipse axes as $a_{[10]}=r/0.9$ and $a_{[01]}=0.9r$. Simulations are performed for two initial precipitate sizes: a smaller one with an area of $20^2 \pi \textrm{ nm}^2$ and a larger one with an area of $75^2 \pi \textrm{ nm}^2$. The center of each precipitate is embedded in the center of a square computational domain, which is given the coordinate (0,0). The computational domain is $(400 \textrm{ nm})^2$ for the smaller precipitates and $(1500 \textrm{ nm})^2$ for the larger precipitates to allow long-range elastic fields to decay. No-stress boundary conditions are applied for the displacements, and no-flux boundary conditions are applied for $\eta$. Because our implementation is based on solving for displacements rather than strain, we specify $u_{[10]}=0$ at the top, middle, and bottom of the $y=0$ axis ($e.g.,$ in the [01] direction) and $u_{[01]}=0$ at the top, middle, and bottom of the $x=0$ axis ($e.g.,$ in the [10] direction) to remove the nullspace in the solution. Simulations are run until equilibrium is achieved.

The presence of elastic strain energy or a curved interface will increase the final value of $\eta$ in both the matrix and precipitate phases from the equilibrium value given by the common tangent of $f_{\text{bulk}}$. In addition, the precipitate may change size during the course of the energy relaxation because of the shifting balance between the $f_{\text{bulk}}$ and $f_{\text{el}}$ energy contributions. Because the precipitate volume within the computational domain is much smaller than that of the matrix, a precipitate may shrink entirely away in the process achieving the equilibrium value of $\eta$ in the matrix. To avoid this, the initial value of $\eta$ in the matrix should be set slightly greater than zero. For the simulations with the small particles, we set $\eta^{\text{matrix}}_0=0.0065$, while for the large particles, $\eta^{\text{matrix}}_0=0.005$. In addition, we set $\eta^{\text{precip}}_0=1$ for all simulations.

Overall, there are 8 different parameter variations for this problem. These are labeled (a) through (h).

Parameter values for (a) through (h)

Quantity	Symbol	Value (a)	Value (b)	Value (c)	Value (d)	Value (e)	Value (f)	Value (g)	Value (h)
Radius	$r$	20 $\textrm{nm}$	75 $\textrm{nm}$	20 $\textrm{nm}$	75 $\textrm{nm}$	20 $\textrm{nm}$	75 $\textrm{nm}$	20 $\textrm{nm}$	75 $\textrm{nm}$
Eclipse axes (10)	$a_{[10]}$	$r$	$r$	$r$	$r$	$r$ / 0.9	$r$ / 0.9	$r$ / 0.9	$r$ / 0.9
Eclipse axes (01)	$a_{[01]}$	$r$	$r$	$r$	$r$	0.9 $r$	0.9 $r$	0.9 $r$	0.9 $r$
Precipitate area		20$^2 \pi \textrm{ nm}^2$	75$^2 \pi \textrm{ nm}^2$	20$^2 \pi \textrm{ nm}^2$	75$^2 \pi \textrm{ nm}^2$	20$^2 \pi \textrm{ nm}^2$	75$^2 \pi \textrm{ nm}^2$	20$^2 \pi \textrm{ nm}^2$	75$^2 \pi \textrm{ nm}^2$
Domain size		$\text{(400 nm)}^2$	$\text{(1500 nm)}^2$	$\text{(400 nm)}^2$	$\text{(1500 nm)}^2$	$\text{(400 nm)}^2$	$\text{(1500 nm)}^2$	$\text{(400 nm)}^2$	$\text{(1500 nm)}^2$
Order Parameter (matrix)	$\eta^{\text{matrix}}_0$	0.0065	0.005	0.0065	0.005	0.0065	0.005	0.0065	0.005
Order Parameter (precip)	$\eta^{\text{precip}}_0$	1	1	1	1	1	1	1	1
Elastic stiffness precip	$C^{\text{precip}}_{1111}$	$\text{250 aJ/nm}^3$	$\text{250 aJ/nm}^3$	$\text{275 aJ/nm}^3$	$\text{275 aJ/nm}^3$	$\text{250 aJ/nm}^3$	$\text{250 aJ/nm}^3$	$\text{275 aJ/nm}^3$	$\text{275 aJ/nm}^3$
Elastic stiffness precip	$C^{\text{precip}}_{1122}$	$\text{150 aJ/nm}^3$	$\text{150 aJ/nm}^3$	$\text{165 aJ/nm}^3$	$\text{165 aJ/nm}^3$	$\text{150 aJ/nm}^3$	$\text{150 aJ/nm}^3$	$\text{165 aJ/nm}^3$	$\text{165 aJ/nm}^3$
Elastic stiffness precip	$C^{\text{precip}}_{1212}$	$\text{100 aJ/nm}^3$	$\text{100 aJ/nm}^3$	$\text{110 aJ/nm}^3$	$\text{110 aJ/nm}^3$	$\text{100 aJ/nm}^3$	$\text{100 aJ/nm}^3$	$\text{110 aJ/nm}^3$	$\text{110 aJ/nm}^3$

Example Result at Equilibrium

The final morphologies of the precipitates for problems (a), (c), (e) and (g). The dark pink curve is for variant (a) and (e) and the light pink is for variant (c) and (g). The results indicate that the shape of the initial precipitate does not influence the final shape of the precipitate for the smaller precipitate variant.

Submission Guidelines

All benchmark solutions should be run to equilibrium. The following data should be collected for each upload.

• Global quantities as the simulation evolves including

  • the total free energy, $\;\;\mathcal{F}\;\;$

  • the interfacial free energy, $\;\;\mathcal{F}_{\text{grad}}=\int_V \frac{\kappa}{2} |\nabla \eta|^2 \; dV\;\;$

  • the elastic free energy, $\;\;\mathcal{F}_{\text{el}} = \int_V f_{\text{el}} \; dV\;\;$

  • the area of the precipitate $\;\;\int_V \eta dV\;\;$ and

  • the precipitate lengths $a_{10}$, $a_{01}$ and $a_d$ measured from the center of the drop to the $\eta=0.5$ contour in the $x$ ([10]), $y$ ([01]) and diagonal directions, respectively. The angle used for the diagonal direction is given by $\theta_d$ such that $\tan\theta_d=a_{01}/a_{10}$.

• The $\eta=0.5$ level set contour position at equilibrium or the latest time step.

Evolving Data Format

The evolving data should be stored in a CSV file with columns labeled as time, a_01, a_10, a_d,elastic_free_energy,gradient_free_energy, precipitate_area and total_free_energy. The CSV file should be formatted as a table and have the following form (note that the column ordering is inconsequential),

a_01,a_10,a_d,elastic_free_energy,gradient_free_energy,precipitate_area,time,total_free_energy
19.97429316515008,19.974293165149973,20.140688631434397,6.185957746168657,4.510418048831537,1264.0,0.1,17.72178588199252
19.86315763877582,19.863157638775874,20.098536436029132,6.054959230125162,2.9620862374085544,1264.0,1.1,17.207555522195793
19.906346454363486,19.906346454363465,20.134576060150618,6.021500927987024,2.736582255973086,1264.0,2.1,17.204113636263912
...

The data should be collected frequently during the simulation, but greater than 20 data points at a minimum (more than 1000 data points is unnecessary and won't improve resolution). The data should be named all_data in the "Short name of data" box located in the "Data Files" section of the upload form. The 2D radio button should be checked, the entry in the "Name of the x-axis column" box should be time and the entry in the "Name of y-axis column" should be total_free_energy. Only one "Data Files" section upload is required for the evolving data.

Equilibrium Data Format

The equilibrium data should be stored in either a CSV file with columns labeled as x and y. The CSV file should be formatted as a table and have the following form (note that the column ordering is inconsequential),

x,y
-9.5,-18.615967564796016
-8.5,-18.84790477132558
-7.5,-19.030286708741155
-6.5,-19.158255175095714

The contour data should be in a sequence that enables an ordered traversal of the contour line. The data should be named contour in the "Short name of data" box located in the "Data Files" section of the upload form. The 2D radio button should be checked, the entry in the "Name of the x-axis column" box should be x and the entry in the "Name of y-axis column" should be y. Only one "Data Files" section upload is required for the evolving data.

Please use the upload form to upload your results.