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 here. For this test case you need only one matlab file, pf3d_single_iso.m. You can download it into a directory of your choice 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 files ??? and ??? 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-1.11.0/images/phasefield.tcl.txt -O phasefield.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.  

In the Cahn-Hillard model, the phase field eolves 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$,

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \phi(x) = \frac{1}{2} [ 1 + \tanh(x/\xi) ] }

where Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \xi = 2\sqrt{\varepsilon/U}}
 is the interface width.

Note that Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathbf{x} = (x,y,z)} but the phase field in this solution is independent of Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle y} and Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle z} , i.e. the interface is normal to the Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle x} axis.

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

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \sigma = \int_{-\infty}^{\infty} f(\phi,\nabla\phi) dx = \frac{1}{3} \sqrt{\varepsilon U} }

Algorithm

We represent the field Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \phi(\mathbf{x})} by a three-dimensional array, phi(i,j,k), of size Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_x \times N_y \times N_z} . In Matlab, i goes from 1 to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_x} , j goes from 1 to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_y} , k goes from 1 to Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle N_z} . For simplicity, a uniform grid size Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle 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, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nabla_x \phi} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nabla_y \phi} , Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \nabla_z \phi} , and the second derivative, Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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 Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \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

Planar interface

Matlab

Spherical interface

Matlab

MD++