Single Phase Field Model with Anisotropic Surface Energy

Single Phase Field Model with Anisotropic Surface Energy

Seunghwa Ryu and Wei Cai

This tutorial describes the theory and numerical implementation of a 3D phase field model with anisotropic surface energy. 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_aniiso.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_aniiso.m.txt -O pf3d_single_iso.m

and run it within Matlab.

To run the MD++ test cases, first you need to install and compile MD++ on your computer. Follow the instructions on [Single_Phase_Field_Model_with_Isotropic_Surface_Energy | the previous wiki page].

Theory

Similar to [Single_Phase_Field_Model_with_Isotropic_Surface_Energy | the previous wiki page], the fundamental degree of freedome in this model is a phase 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})} . The only difference is that the gradient term in the the free energy density function now depends on the local orientation of the gradient.

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 F[\phi(\mathbf{x})] = \int f(\phi(\mathbf{x}), \nabla\phi(\mathbf{x})) d^3 \mathbf{x} }

$f(\phi ,\nabla \phi )={\frac {1}{2}}(\epsilon (\mathbf {n} ))^{2}(\nabla \phi )^{2}+U\phi ^{2}(1-\phi )^{2}$

Notice that we have replaced the constant $\varepsilon$ in the previous formulation by ${\frac {1}{2}}(\epsilon (\mathbf {n} ))^{2}$, where $\mathbf {n}$ is a unit vector specifying the local orientation of $\nabla \phi (\mathbf {x} )$.

(the following are copied from pf3d_single_iso.m wiki page)

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

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 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}) = 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