In [1]:
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
$('div.prompt').hide();
} else {
$('div.input').show();
$('div.prompt').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Code Toggle"></form>''')
Out[1]:
In [1]:
from IPython.display import HTML
HTML('''
<a href="https://raw.githubusercontent.com/{{ site.links.repo }}/master/benchmarks/benchmark4.ipynb"
download="benchmark4.ipynb">
<button type="submit">Download Notebook</button>
</a>''')
Out[1]:
In [4]:
from IPython.display import HTML
HTML('''{% include jupyter_benchmark_table.html num="[4]" revision=1 %}''')
Out[4]:
See the journal publication entitled ["Phase Field Benchmark Problems for Dendritic Growth and Linear Elasticity"][nolink] for more details about the benchmark problems. Furthermore, read the extended essay for a discussion about the need for benchmark problems.
[nolink]: {{ site.baseurl }}
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]).
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}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.
| 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).
| 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$ |
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.
All benchmark solutions should be run to equilibrium. The following data should be collected for each upload.
All fields from the equilibrium state of the system including $\eta$, $u_x$ and $u_y$.
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,\;\;$ as well as the volume of the precipitate $\;\;\int_V \eta dV$.
The equilibrium field data should be stored in a CSV or JSON file with x, y, eta, u_x and u_y fields. The CSV file should be formatted as a table and the JSON file should be formatted as a list with each list item as a set of key value pairs. The spatial data is not required to be structured, but just a point cloud in space. If the data is not colocated, then values can be dropped from CSV columns (by using two commas together) or from a JSON file by either excluding the key or setting a null value.
The global quantity data should be stored in either a CSV file or JSON file with time, free_energy ,free_energy_grad, free_energy_el and precip_volume fields. The CSV file should be formatted as a table and the JSON file should be formatted as a list with each list item as a set of key value pairs. Values should be collected as frequently as possible, but at least 10 values at a minimum.
Please use the upload form to upload your results.
In [ ]: