KKS Binary Solidification with FiPy

This notebook contains code for binary solidification using the Kim-Kim-Suzuki model [1]
for interfacial energy. This allows easy specification of gamma, but requires constant
chemical potential through the interface. The implementation involves iteratively solving
for chemical composition in pure phases such that the chemical potential constraint is
satisfied [2].

Questions/comments to trevor.keller@nist.gov (Trevor Keller). Cite with the following DOI:

References:

1. Kim, Kim, and Suzuki. "Phase-field model for binary alloys."
   Physical Review E 60:6;7186-7197 (1999).
2. Provatas and Elder. Phase-Field Methods in Materials Science and Engineering,
   Chapter 6, Section 9. Wiley VCH: Weinheim, Germany. 2010.

Model Description (paraphrased after Provatas & Elder)

We are setting out to simulate solidification of two-component alloy with a lenticular
phase diagram, or "binary isomorphous" system, such as Cu-Ni. The free energy curves for
pure phases, $f_S$ and $f_L$, are generated from a CALPHAD database, rather than a regular
solution model.

In addition to this thermodynamic description, we are adopting the KKS treatment of
diffuse interfaces. This simply means that at equilibrium, chemical potential is constant
through the interface, and composition varies to make it so. More concretely, composition
is defined by the phase fraction $\phi$ and two fictitious concentration fields, $C_S$ and
$C_L$, representing composition of the pure phase, as

$$c = h(\phi)C_S + (1-h(\phi))C_L,$$

where the interpolation function $h(\phi)=\phi^3(6\phi^2-15\phi+10)$ takes the values in
solid $h(\phi=1)=1$ and liquid $h(\phi=0)=0$. At equilibrium,

$$\mu = \left.\frac{\partial f_S}{\partial c}\right|_{c=C_S} = \left.\frac{\partial f_L}{\partial c}\right|_{c=C_L}.$$

Taken together, and introducing a double-well function $g(\phi)=\phi^2(1-\phi)^2$, the
thermodynamic and interfacial treatments provide the bulk free energy,

$$f(\phi,c,T) = \omega g(\phi) + h(\phi)f_S(C_S,T) + (1-h(\phi))f_L(C_L,T).$$

Now, assuming nonconserved (Allen-Cahn) dynamics for $\phi$ and conserved (Cahn-Hilliard)
dynamics for $c$, we can write the equations of motion

$$\frac{\partial\phi}{\partial t} = -M_\phi\frac{\delta\mathcal{F}}{\delta\phi} \rightarrow \tau\frac{\partial\phi}{\partial t} = \epsilon_\phi^2\nabla^2\phi -\omega g'(\phi) + h'(\phi)\left[f_L(C_L) - f_S(C_S) - \frac{\partial f_L(C_L)}{\partial c}(C_L - C_S)\right]$$

$$\frac{\partial c}{\partial t} = \nabla\cdot M_c\nabla\frac{\delta\mathcal{F}}{\delta c} = M_c\nabla\cdot Q(\phi)\left[h(\phi)\nabla C_S + (1-h(\phi))\nabla C_L\right]$$

with phase-dependent mobility $Q(\phi)=\frac{1-\phi}{(1+k) - (1-k)\phi}$, partition
coefficient $k=\frac{C_S^e}{C_L^e}$, and time constant $\tau = M_\phi^{-1}$.

The Wrinkle

$C_S$ and $C_L$ are not constants, they are field variables whose values depend on $\phi$
and $c$. Determining their values requires solving for the common tangent, or the coupled
roots

$$f_1(C_S,C_L) = h(\phi)C_S + (1-h(\phi))C_L -c = 0$$

$$f_2(C_S,C_L) = \frac{\partial f_S(C_S)}{\partial c} - \frac{\partial f_L(C_L)}{\partial c} = 0.$$

While these equations can be solved using Newton's Method (cf. Provatas & Elder Appendix
C.3), it's better to invoke a library, which we'll do a little later on. Even with a
highly optimized root solver, determining $(C_S,C_L)$ at every grid point gets
prohibitively expensive. The standard approach is to construct a lookup table for $C_S$
and $C_L$ covering $\phi=[-\delta,1+\delta]$ and $c=[-\delta,1+\delta]$ with a reasonably
high number of points, then interpolating from the LUT at runtime. For best results, use
the interpolated values as initial guesses for a touch-up iteration or two.

Running Locally

It is strongly recommended that you download this software using
git, the distributed version control software:

$ git clone https://github.com/tkphd/KKS-binary-solidification.git

Python

The Python version of this software depends on FiPy
and pycalphad. It is probably best to install these
in a conda environment. Please follow
the directions provided by those packages to configure your system. On Linux, installation
will look something like the following:

$ wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ chmod +x Miniconda3-latest-Linux-x86_64.sh 
$ ./Miniconda3-latest-Linux-x86_64.sh
# Complete the installation process
$ conda create -n kks -c conda-forge python=3 fipy notebook numpy pycalphad scipy
$ conda activate kks
$ cd KKS-binary-solidification
$ jupyter notebook

This should open a web browser in the folder you cloned this repository into. Click on the
file named "FiPy-KKS.ipynb" and tinker to your heart's content.

C++

The C++ version of this code depends on MMSP and
GSL. MMSP is a header-only library, so you need only
download it and set an environmental variable. In Linux, this is simple:

$ git clone https://github.com/mesoscale/mmsp.git
$ echo "export MMSP_PATH=${PWD}/mmsp" >> ~/.bashrc
$ . ~/.bashrc

Please follow the MMSP documentation if you wish to build utilities for
conversion of output between various file formats.

If you do not already have them installed, you will need to install
Make, and headers for
libpng and zlib.
With these dependencies satisfied, you should be able to build and run the program:

$ cd ~/KKS-binary-solidification
$ make
$ ./KKS --help
...
$ ./KKS --example 2 start.dat
System has 13.22% solid, 86.78% liquid, and composition 46.77% B. Equilibrium is 50.00% solid, 50.00% liquid.

Equilibrium Cs=0.54, Cl=0.39. Timestep dt=9.00e-04

At this point, the lookup table for Cs and Cl has been written to consistentC.lut,
and visualizations of the various fields were exported as PNG images. Details are in the
generate function in KKS.cpp. To evolve the system, run

$ ./KKS start.dat 10000 1000

or similar; use ./KKS --help again for details. If you built the utilities, you can
convert all the checkpoint files to PNG images using

$ for f in *.dat; do mmsp2png $f; done