Single Phase Field Model with Isotropic Surface Energy

Single Phase Field Model with Isotropic Surface Energy

Seunghwa Ryu and Wei Cai

This tutorial describes the theory and numerical implementation of a basic 3D phase field model. The model has been implemented in both matlab and in MD++. We describe the algorithm of these codes and how to use them in detail. A brief introduction of the phase field model can be found in Chapter 11 of

Computer Simulations of Dislocations, V. V. Bulatov and W. Cai, Oxford University Press (2006).

---

Download files

The matlab program can be downloaded from this link: pf3d_single_iso.m. For this test case you need only one matlab file. You can download it into a directory of your choice, using the following command

wget http://micro.stanford.edu/mediawiki/images/5/5b/Pf3d_single_iso.m.txt -O pf3d_single_iso.m

and run it within Matlab.

To run the MD++ test cases, first you need to install MD++ on your computer by following the instructions on MD++ Manuals.

Second, copy file Pf3d.tcl to your input file directory of MD++. You can do so by the following commands (assuming you have installed MD++ in ~/Codes/MD++).

 export MDPP=~/Codes/MD++
 mkdir -p ${MDPP}/scripts/work/phasefield
 cd ${MDPP}/scripts/work/phasefield
 wget http://micro.stanford.edu/mediawiki/images/6/6b/Pf3d.tcl.txt -O pf3d.tcl

Compile executable file

This step is only necessary for running MD++ test cases. In this tutorial, we describe how to run these examples on a linux computer (using mc2.stanford.edu as an example).

You can can compile the executable phasefield by typing

 cd ${MDPP}
 make phasefield build=R SYS=mc2

If compilation is successful, this will create binary file sworig_mc2 file in your ${MDPP}/bin directory. To compile MD++ on a different system, you should set variable SYS to a different value.

Theory

In this simple phase field model, the fundamental degree of freedom is a phase field, $\phi (\mathbf {x} )$. The range of the value for $\phi$ is between 0 and 1. For example, we can think of the region where $\phi =1$ as the solid phase and the region where $\phi =0$ as the liquid phase. The region where $\phi$ transitions (rapidly) from 0 to 1 is then the liquid-solid interface.

A free energy functional $F[\phi (\mathbf {x} )]$ is defined as the integral of the local free energy density $f(\phi ,\nabla \phi )$,

$F[\phi (\mathbf {x} )]=\int f(\phi (\mathbf {x} ),\nabla \phi (\mathbf {x} ))d^{3}\mathbf {x}$
$f(\phi ,\nabla \phi )=\varepsilon (\nabla \phi )^{2}+U\phi ^{2}(1-\phi )^{2}$

The system evolves in the direction that reduces the free energy. In the Ginzburg-Landau model, the phase field evolves with the following equation of motion that does not conserve the total amount of the solid phase (or liquid phase),

${\frac {\partial \phi }{\partial t}}=-M{\frac {\delta F}{\delta \phi }}$  (GL)
where ${\frac {\delta F}{\delta \phi }}={\frac {\partial f}{\partial \phi }}-\nabla \left({\frac {\partial f}{\partial \nabla \phi }}\right)$ is the variational derivative.  

For the free energy density function given above, the explicit expression of the variational derivative is,

${\frac {\delta F}{\delta \phi }}=2U\phi (1-\phi )(1-2\phi )-2\varepsilon \nabla ^{2}\phi$

In the Cahn-Hillard model, the phase field evolves with the following equation of motion that conserve the total amount of the solid phase (or liquid phae),

${\frac {\partial \phi }{\partial t}}=\nabla \cdot \left[M\nabla \left({\frac {\delta F}{\delta \phi }}\right)\right]$ (CH)

Note here that we have been treating $f$ as a function with two independent input variables, which take the values of $\phi$ and $\nabla \phi$, respectively.

A stationary, planar interface (between liquid and solid) is described by the following solution that satisfies the boundary conditions: $\phi (x=+\infty )=1$, $\phi (x=-\infty )=0$,

$\phi (x)={\frac {1}{2}}[1+\tanh(x/\xi )]$
where $\xi =2{\sqrt {\varepsilon /U}}$ is the interface width.

Note that $\mathbf {x} =(x,y,z)$ but the phase field in this solution is independent of $y$ and $z$, i.e. the interface is normal to the $x$ axis.

The interface energy (per unit area) of this solution is

$\sigma =\int _{-\infty }^{\infty }f(\phi ,\nabla \phi )dx={\frac {1}{3}}{\sqrt {\varepsilon U}}$

Algorithm

We represent the field $\phi (\mathbf {x} )$ by a three-dimensional array, phi(i,j,k), of size $N_{x}\times N_{y}\times N_{z}$. In Matlab, i goes from 1 to $N_{x}$, j goes from 1 to $N_{y}$, k goes from 1 to $N_{z}$. For simplicity, a uniform grid size $h$ is used in all three directions. By default, periodic boundary conditions (PBC) are applied in all three directions.

At each iteration, the free energy and its variational derivative are computed from the second and first spatial derivatives of the phase field, respectively. These derivatives are computed from centered numerical difference. First, the following indices are introduced to represent the field point neighboring the point of interest in the +x, -x, +y, -y, +z, -z directions, respectively (satisfying PBC).

pxind=[2:N_grid(1),1]; nxind=[N_grid(1),1:N_grid(1)-1];
pyind=[2:N_grid(2),1]; nyind=[N_grid(2),1:N_grid(2)-1];
pzind=[2:N_grid(3),1]; nzind=[N_grid(3),1:N_grid(3)-1];

The first derivatives, $\nabla _{x}\phi$, $\nabla _{y}\phi$, $\nabla _{z}\phi$, and the second derivative, $\nabla ^{2}\phi$, are computed the following manner.

dphidx = (phi(pxind,:,:)-phi(nxind,:,:))/(2*h);
dphidy = (phi(:,pyind,:)-phi(:,nyind,:))/(2*h);
dphidz = (phi(:,:,pzind)-phi(:,:,nzind))/(2*h);
d2phi  = (phi(pxind,:,:)+phi(nxind,:,:)+phi(:,pyind,:)+phi(:,nyind,:) ...
         +phi(:,:,pzind)+phi(:,:,nzind)-6*phi)/(h^2);

Based on these derivatives, the free energy and its variational derivative are computed as,

F = (U * sum(sum(sum( phi.^2.*(1-phi).^2 ))) ...
   + vareps * sum(sum(sum( dphidx.^2 + dphidy.^2 + dphidz.^2)))) * (h^3);
dFdphi = 2*U*phi.*(1-phi).*(1-2*phi) - 2*vareps * d2phi;

The forward-Euler method is used for time integration, i.e.

phi = phi + dt* dphidt;

For Ginzburg-Landau model (equation_type == 0), we use

dphidt = -M*dFdphi;

For Cahn-Hillard model (equation_type == 1), we use

d2dFdphi = (dFdphi(pxind,:,:)+dFdphi(nxind,:,:) ...
           +dFdphi(:,pyind,:)+dFdphi(:,nyind,:) ...
           +dFdphi(:,:,pzind)+dFdphi(:,:,nzind)-6*dFdphi)/(h^2);
dphidt = M*d2dFdphi;

We have also implemented a Ginzburg-Landau model with a constraint on the total volume of each phase (equation == 2), in which we use

dphidt = -M*(dFdphi - mean(mean(mean(dFdphi))));

The average value of dFdphi is subtracted off (similar to a Lagrangian multiplier) so that the simulation preserves the volume integral of $\phi (\mathbf {x} )$ in space. This is equivalent to adjusting the relative chemical potential of the solid/liquid phases through a feedback mechanism in order to keep the total amount of solid (liquid) constant.

Test cases

Random initial configuration

A simple test case is to initialize the phase field values to random numbers (uniformly distribution in [0,1]) and to watch how it evolves with time. This can be accomplished by typing the following commands in Matlab,

clear all
equation_type = 0; dt = 1e-3; pf3d_single_iso

Here equation_type = 0 specifies the use of Ginzburg-Landau equation (no constraint on phases). dt = 1e-3 specifies the time step. You should see the isosurfaces where $\phi (\mathbf {x} )=0.5$, depicting the solid-liquid interface, which evolves with time. Typing pf3d_single_iso again (without clear all) continues the simulation for another Niter = 5000 steps. If the simulation is continued for many iterations, ultimately one phase will overtake the other and there will be no interface visible in the simulation cell.

To perform the simulation using the Cahn-Hillard equation (again starting from random configuration), type the following commands in Matlab,

clear all
equation_type = 1; dt = 1e-4; pf3d_single_iso

Note that a smaller time step is required for the Cahn-Hillard equation. You should see that the solid and liquid regions interpenetrate each other, forming a tortuous microstructure that coarsens very slowly.

We can continue this simulation with the Ginzburg-Landau equation but with the constraint that the total amount of solid phase is conserved. Type the following commands in Matlab,

equation_type = 2; dt = 1e-2; pf3d_single_iso

Notice that here we can afford a much larger time step (without encountering instability) because the starting point of this simulation is a partially relaxed (smooth) phase field from a previous simulation.

If we continue the simulation in this setting for a long time, one possible final outcome is for the system will evolve into two planar phase boundaries.

(MD++ test cases. The code is ready. Just need to prepare the tcl script and write up the description. In MD++, we cannot plot isosurfaces, but can only visualize the phase field by plotting a circle at grid points where the phase field value lies within a window (e.g. [0.45, 0.55]). However, MD++ can save simulation snapshots into .cn files, which can be loaded into matlab and plotted using plot_mdpp_phasefield.m.)

Planar interface

Depending on your initial condition, it may take a long time for the phase field simulation (with equation_type == 2) to reach the state of planar interfaces. If we wish to examine the state of planar interfaces, we can facilitate this by preparing the initial condition to a similar state.

clear all
N_grid = [50 50 50];
phi = rand(N_grid)* 0.6;
for i=round(N_grid(1)/4):round(N_grid(1)*3/4), phi(i,:,:)=phi(i,:,:)+0.3; end
equation_type = 2; dt = 1e-3; Niter = 2000; pf3d_single_iso

After the relaxation has complted, you can plot the field variation normal to the phase boundary by

plot(x,phi(:,1,1))

Spherical interface

We can also initialize the solid phase to be a cube and watch it evolve into a sphere (with equation_type == 2).

clear all
N_grid = [50 50 50];
phi = zeros(N_grid);
phi(12:37,12:37,12:37) = 1.0;
equation_type = 2; dt = 1e-3; Niter = 2000; pf3d_single_iso