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 [2]:
from IPython.display import HTML

HTML('''
<a href="{{ site.links.github }}/raw/nist-pages/hackathons/hackathon2/problem2.ipynb"
   download>
<button type="submit">Download Notebook</button>
</a>
''')





In [1]:
from IPython.display import HTML

HTML('''{% include jupyter_benchmark_table.html num="[4]" revision=0 %}''')


Out[1]:
{% include jupyter_benchmark_table.html num="[4]" revision=0 %}

Benchmark Problem 4: Linear Elasticity in 3D

The linear elastic energy of a body is

\begin{equation} E_{\rm elastic}=\frac{1}{2}\int \sigma_{ij}\epsilon_{ij}\,dV=\frac{1}{2}\int C_{ijkl}\epsilon_{ij}\epsilon_{kl}\,dV, \end{equation}

where $\sigma_{ij}=C_{ijkl}\epsilon_{kl}$ is the stress, $C_{ijkl}$ is the elastic tensor, and $\epsilon_{ij}$ is the strain, \begin{equation} \epsilon_{ij}=\frac{1}{2}\left[\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right], \end{equation} with $u_i$ the displacement field. The indices $i,j,k,l$ run from 1 to 3, and $x_i,\,i=1,2,3$, are Cartesian coordinates; we are using Einstein summation convention so repeated indices are summed over.

The elastic tensor obeys symmetries $C_{ijkl}=C_{jikl}=C_{ijlk}=C_{jilk}$ and $C_{ijkl}=C_{klij}$. These symmetries imply that there are only 21 independent entries in the elastic tensor. Usually Voigt notation is used, in which the four indices $ijkl$ are replaced by two indices $IJ$. The mapping for each pair $ij$ (or $kl$) is $11\to1$, $22\to2$, $33\to3$, $23\to4$ (and $32\to4$), $13\to5$ (and $31\to5$), and $12\to6$ (and $21\to6$). The crystal symmetry may further reduce the number of independent entries. In an orthorombic crystal, there are only nine independent entries, and they are (in Voigt notation) $C_{11}, C_{22}, C_{33}, C_{44}, C_{55}, C_{66}, C_{12}, C_{13}$, and $C_{23}$. The tensor $C_{IJ}$ thus has the form

\begin{equation} \left( \begin{matrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0\\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{13} & C_{23} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0\\ 0 & 0 & 0 & 0 & C_{55} & 0\\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{matrix} \right). \end{equation}

For tetragonal symmetry, there are six independent entries, $C_{11}, C_{33}, C_{44}, C_{66}, C_{12}$, and $C_{13}$.

Aluminum silicate, Al$_2$SiO$_5$ is a crystal with orthorombic symmetry and unit cell parameters $a=7.738\;\unicode{x212B}$ , $b=7.857\;\unicode{x212B}$, and $c=5.534\;\unicode{x212B}$. The elastic tensor is given by $C_{11}=233.4$, $C_{22}=289.0$, $C_{33}=380.1$, $C_{44}=99.5$, $C_{55}=87.8$, $C_{66}=112.3$; $C_{11}+C_{22}-2C_{12}=233.4$, $C_{11}+C_{33}-2C_{13}=380.9$, and $C_{22}+C_{33}-2C_{23}=506.3$, all in units of GPa.

(a) (Potato in space) What is the equilibrium shape of a 0.0042~$\mu$m$^3$ volume of Al$_2$SiO$_5$ in free space (stress-free boundaries)? Take the surface energy, $\gamma$, to be equal to 200 mJ/m$^2$. The crystalline axes $a$, $b$, and $c$ are aligned with the $x$, $y$, and $z$-axes of a Cartesian lab coordinate system.

This problem can be cast as a phase field problem, where the phase field $\varphi\in[0,1]$ takes the value of 1 in one phase (the "potato"), and a value of 0 in the other (the surrounding vacuum). Thus, the total free energy can be written \begin{equation} {\mathcal F}=\int \left[f_{\rm elastic}+\frac{\kappa}{2}|\nabla\varphi|^2+h_0f(\varphi)\right]\,dV, \end{equation} with the integral extended over all space, where $f(\varphi)$ is a (dimensionless) double-well function \begin{equation} f(\varphi)=\varphi^2\left[\varphi-1\right]^2, \end{equation} and $h_0$ has the dimension of energy per unit volume. The interface width $W$ between approximately $\varphi=0.1$ and $\varphi=0.9$ in this model is given by $W=2\sqrt{2\kappa/h_0}$, while $\gamma=\sqrt{\kappa h_0/18}$, $\kappa=1.5\gamma W$, and $h_0=12\gamma/W$ (ignoring modification of the phase field order parameter $\varphi$ by the elastic interactions through the interface). Use $\kappa=3\times10^{-9}$~J/m, and $h_0=2.4\times10^{8}$~J/m$^3$.

We use a simple interpolation of the elastic constants,

\begin{equation} C_{ijkl}=h(\varphi)C_{ijkl}^{\rm potato},%+\left[1-h(\varphi)\right]C_{ijkl}^{\rm matrix}, \end{equation}

where $h(\varphi)$ is a smooth interpolation function, \begin{equation} h(\varphi)=\varphi^3\left[6\varphi^2-15\varphi+10\right], \end{equation} that interpolates between $h(\varphi=0)=0$ and $h(\varphi=1)=1$.

Hint: find time-evolution equations for $\varphi$ that monotonically drive the total energy to a minimum while preserving the volume. One way to do this is to set up a Cahn-Hilliard equation for $\varphi$.

(b) (Compressed potato in space) What is the equilibrium shape of a 0.0042 $\mu$m$^3$ volume Al$_2$SiO$_5$ with uniaxial compressive stress of 500~MPa applied in the $z$-direction? Note that the strain must go to zero far from the "potato."

(c) (Potato in a stew) A volume of Al$_2$SiO$_5$ is embedded coherently in a matrix of tetragonal symmetry with unit cell parameters $a_m=b_m=7.6918\;\unicode{x212B}$ and $c_m=5.5674\;\unicode{x212B}$ such that $a$ and $a_m$, $b$ and $b_m$ and $c$ and $c_m$ are pairwise aligned. The elastic tensor of the matrix is given by $C_{11}=C_{22}=269.0$, $C_{12}=177.0$, $C_{13}=146.0$, $C_{33}=480.0$, $C_{44}=124.0$, and $C_{66}=192.0$, all in units of GPa.

This problem can also be cast as a phase field problem, where now the phase field $\varphi\in[0,1]$ takes the value of 1 in one phase (the inclusion), and a value of 0 in the other (the matrix). The total free energy is again \begin{equation} {\mathcal F}=\int \left[f_{\rm elastic}+\frac{\kappa}{2}|\nabla\varphi|^2+h_0f(\varphi)\right]\,dV, \end{equation} where $f(\varphi)$ is given by above.

The interpolation of the elastic constants is now,

\begin{equation} C_{ijkl}=h(\varphi)C_{ijkl}^{\rm inclusion}+\left[1-h(\varphi)\right]C_{ijkl}^{\rm matrix}, \end{equation}

where again $h(\varphi)$ is a smooth interpolation function given above.

In the elastic energy, the relevant strain is the total strain $\epsilon_{ij}$ minus the local misfit strain $\epsilon^0_{ij}$, so

\begin{equation} f_{\rm elastic}=\frac{1}{2}C_{ijkl}\left(\epsilon_{ij}-\epsilon^0_{ij}\right)\left(\epsilon_{kl}-\epsilon^0_{kl}\right). \end{equation}

We will use Vegard's law to determine the local misfit strain, for which we just interpolate the crystallographic misfit strain tensor, $\epsilon^T_{ij}$:

\begin{equation} \epsilon^0_{ij}=h(\varphi)\epsilon^T_{ij}, \end{equation}

where $h(\varphi)$ is the interpolation function, and $\epsilon^T_{ij}$ is \begin{equation} \epsilon^T_{ij}=\left( \begin{matrix} \frac{a-a_m}{a_m} & 0 & 0\\ 0 & \frac{b-b_m}{b_m} & 0\\ 0& 0 & \frac{c-c_c}{c_m} \end{matrix} \right). \end{equation}

Find the equilibrium shape of an isolated inclusion with a volume of 0.0042 $\mu$m$^3$ for $\kappa=3\times10^{-9}$ J/m, and $h_0=2.4\times10^{8}$ J/m$^3$. Use as an initial condition a prolate ellipsoid with $a=c=0.155$ $\mu$m and $b=0.042$ $\mu$m. Because of elastic strain energy and interfacial energy within the system, you will need to increase the intial values of the phase field variable in the matrix and the precipitate (e.g., $\varphi$ in the precipitate $\approx$ 1.05; $\varphi$ in the matrix $\approx 0.05$). Note that the strain must go to zero far from the "potato".

For parts (a) - (c), present a rendering of the final shape as well as the lengths of the "potato" along the principal axes, and also track the total free energy as function of "time" and plot it.