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/problem1.ipynb"
download>
<button type="submit">Download Notebook</button>
</a>
''')
Out[2]:
In [1]:
from IPython.display import HTML
HTML('''{% include jupyter_benchmark_table.html num="[3]" revision=0 %}''')
Out[1]:
The solidification of an undercooled liquid can be described by a phase field model. The free energy density is given by \begin{equation} {\mathcal F}=\int \left[\frac{1}{2} W^2(\hat n)|\nabla \varphi|^2+f(\varphi,u)\right]\,dV. \end{equation} Here, $u$ is a dimensionless temperature field, \begin{equation} u=\frac{T({\mathbf r},t)-T_m}{L \; / \; c_p}, \end{equation} where $T({\mathbf r},t)$ is the space- and time-dependent temperature, $T_m$ the melting temperature, $L$ the latent heat, $c_p$ the specific heat at constant pressure, and $\varphi$ is a phase field variable which describes the phase of the material. The free energy density $f(\varphi,u)$ is a double-well potential, \begin{equation} f(\varphi,u)=-\frac{1}{2}\varphi^2+\frac{1}{4}\varphi^4 +\lambda u\varphi\left[1-\frac{2}{3}\varphi^2+\frac{1}{5}\varphi^4\right], \end{equation} constructed to have minima at $\varphi=\pm1$, corresponding to solid ($\varphi=1$) and liquid phases ($\varphi=-1$), respectively. The variable $\lambda$ is a dimensionless coupling constant, and the last term on the right-hand side tilts the free energy double-well to favor solid or liquid, depending on if the temperature is above or below the melting temperature $T_m$.
The term $W(\hat n)$ controls the interface thickness and contains the crystalline anisotropy or symmetry. In two dimensions, it takes the form $ W(\hat n)=W_0a(\hat n)$, where the unit normal vector $\hat n=\frac{\nabla\varphi}{|\nabla\varphi|}$. We will use a simple form for $a(\hat n)$ to reflect in-plane symmetry, $a(\hat n)=1+\epsilon_m\cos(m(\theta-\theta_0))$, where $m$ is a non-negative integer and $\theta$ is the in-plane azimuthal angle, $\tan\theta=n_y\;/\;n_x$; $\theta_0$ is an offset azimuthal angle that specifies the orientation of the crystalline axis relative to the $x$-axis in the laboratory coordinate system.
The time evolution of the fields $u$ and $\varphi$ are given by \begin{eqnarray} \frac{\partial u}{\partial t}&=&D\nabla^2u+\frac{1}{2}\frac{\partial \varphi}{\partial t}\\ \tau(\hat n)\frac{\partial\varphi}{\partial t} & = & - \frac{\delta{\mathcal F}}{\delta \varphi({\mathbf r},t)}, \end{eqnarray} where $D$ is a thermal diffusion constant and $\tau(\hat n)=\tau_0\left[a(\hat n)\right]^2$. We will take \begin{equation}\lambda=\frac{D\tau_0}{0.6267W_0^2} \end{equation} for obscure technical reasons (it makes kinetic contributions to the temperature on the interface vanish). Because of the dependence of $W$ on $\varphi$, the functional derivative becomes a little complicated:
$$ \begin{align} \tau(\hat n)\frac{\partial\varphi}{\partial t} & = & \left[\varphi-\lambda u\left(1-\varphi^2\right)\right]\left[1-\varphi^2\right]+\nabla\cdot\left[ W(\hat n)^2\nabla\varphi\right]\nonumber\\ &&+\frac{\partial}{\partial x}\left[|\nabla\varphi|^2W(\hat n)\frac{\partial W(\hat n)}{\partial \left(\frac{\partial\varphi}{\partial x}\right)}\right] +\frac{\partial}{\partial y}\left[|\nabla\varphi|^2W(\hat n)\frac{\partial W(\hat n)}{\partial \left(\frac{\partial\varphi}{\partial y}\right)}\right]. \end{align} $$We will use the following values: $W_0=1$, $\tau_0=1$, and $D=10$. The degree of undercooling, $\Delta$, is given by \begin{equation}\Delta=\frac{T_m-T_0}{L/c_p}, \end{equation} where $T_0$ is the initial (uniform) temperature of the system $[T({\mathbf r},t=0)=T_0]$, and also specifies the boundary condition $u({\mathbf r},t)\to-\Delta$ as $|{\mathbf r}|\to\infty$.
Use as initial condition that $\varphi({\mathbf r},t=0)=1$ for $|{\mathbf r}|\leq1$ and $\varphi({\mathbf r},t=0)=-1$ for $|{\mathbf r}|>1$, such that the sphere of solid is centered in the computational domain. This problem involves many different length scales which may stress your meshing and solver. You should be careful to use a mesh that can resolve the interface (e.g., 5 mesh points across the interface) while still being able to include regions where the temperature $T\approx T_0$. The solution will eventually reach a steady state where the tip velocity of the solidification front is constant if the computational domain is large enough. We will embed the system in a bounding box $\Omega$ of size $L\times L$ with the boundaries aligned with the $x$- and $y$-axes, and with boundary $\partial \Omega$. We use Dirichlet boundary conditions $T=T_0$ so that $u=-\Delta$ on the boundary $\partial \Omega$.
Plot the tip velocity as a function of time for (a), (b), and (c) below. Find the velocity both by directly measuring from your solution, and by calculating it via \begin{equation} v_n=D\hat n\cdot\left[\nabla u|_S-\nabla u|_L\right], \end{equation} where $\nabla u|_S$ and $\nabla u|_L$ are the gradients of $u$ evaluated on the solid and liquid side of the interface, respectively.
(a) Plot the tip velocity of the solidification front as a function of time for $m=4$, $\Delta=0.05$, $\epsilon_4=0.025$, $\theta_0=0$, $L=20$, and also $L=100$.
(b) Plot the tip velocity of the solidification front as a function of time for $m=4$, $\Delta=0.05$, $\epsilon_4=0.025$, $\theta_0=\pi/4$, $L=20$, and $L=100$.
(c) Plot the tip velocity of the solidification front as a function of time for $m=6$, $\Delta=0.05$, $\epsilon_6=0.025$, $\theta_0=0$, and $L=100$.