Multi Phase Field Model with Isotropic Interface Energy

Multi Phase Field Model with Isotropic Interface Energy

Seunghwa Ryu, Yanming Wang and Wei Cai

This tutorial describes the theory and numerical implementation of a 3D multi phase phase field model with isotropic interface 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 multi iso.m. For the 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/b/b7/Pf3d_multi_iso.m.txt -O Pf3d_multi_iso.m

and run it within Matlab.

Theory

The fundamental degrees of freedom in multi phase field model are a series of phase fields $\phi _{i}(x)$,where i=1:N and N is the number of different phases in the system. The range of the value of each phase field $\phi _{i}$ is between 0 and 1. We can think of the region where $\phi _{i}=1$ as in the ith phase and the region where $\phi _{i}=0$ as not in the ith phase. Thereby, Similar to the previous wiki page but considering the interfaces between different phases, the free energy functional and the local free energy density are given by

$F[\phi _{1}(\mathbf {x} ),...,\phi _{N}(\mathbf {x} )]=\int \sum _{k\neq {i}}{f_{ik}(\phi _{i}(\mathbf {x} ),\nabla \phi _{i}(\mathbf {x} ),\phi _{k}(\mathbf {x} ),\nabla \phi _{k}(\mathbf {x} ))d^{3}\mathbf {x} }$
$f_{ik}(\phi _{i},\nabla \phi _{i},\phi _{k},\nabla \phi _{k})=\varepsilon _{ik}(\phi _{i}\nabla \phi _{k}-\phi _{k}\nabla \phi _{i})^{2}+U_{ik}\phi _{i}^{2}\phi _{k}^{2}$

Intuitively, if we let $\phi _{k}=1-\phi _{i}$, then $f_{ik}$ will be reduced to a formalism similar to the single phase model.

$f_{ik}(\phi _{i},\nabla \phi _{i},\phi _{k},\nabla \phi _{k})=\varepsilon _{ik}|\nabla \phi _{i}|^{2}+U_{ik}\phi _{i}^{2}(1-\phi _{i})^{2}$

In the Ginzburg-Landau model, The evolution of the ith phase without constraint can be expressed with the following equation of motion,

${\frac {\partial \phi _{i}}{\partial t}}=-M_{ik}\sum _{k\neq i}{\frac {\delta F_{ik}}{\delta \phi _{i}}}$  (GL)

where $M_{ik}$ is the mobility of i-k interface, and $F_{ik}=\int f_{ik}(\phi _{i},\nabla \phi _{i},\phi _{k},\nabla \phi _{k})d^{3}\mathbf {x}$ 

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

${\frac {\delta F_{ik}}{\delta \phi _{i}}}=2U_{ik}\phi _{i}\phi _{k}(\phi _{i}-\phi _{k})-2\varepsilon _{ik}(\phi _{k}\nabla \phi _{i}^{2}-\phi _{i}\nabla \phi _{k}^{2})$

In the Cahn-Hillard model, the system evolves with the following equation of motion that conserve the total amount of each phase,

 ${\frac {\partial \phi _{i}}{\partial t}}=\nabla \cdot \left[M_{ik}\nabla \left({\frac {\delta F_{ik}}{\delta \phi _{i}}}\right)\right]$ (CH)

Algorithm

In matlab, We present liquid $\phi _{1}$, solid $\phi _{2}$, and vapor $\phi _{3}$, three phase fields by a $1\times 3$ cell. Each element phi{i}(m,n,p) of the cell is a three-dimensional array, of size $N_{x}\times N_{y}\times N_{z}$. i goes from 1 to 3, m goes from 1 to $N_{x}$, n goes from 1 to $N_{y}$, and p goes from $\mathbb {N} _{z}$. For simplicity, a uniform grid size $h$ is used in all three directions. Similar to the previous wiki page, by default, periodic boundary conditions (PBC) are applied in all three directions.

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

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

Based on these derivatives, the free energy for each i-k pair and the total free energy are computed as,

for i=1:N_phases,
       for k=i+1:N_phases,
           grad_term = (phi{i}.*dphidx{k} - phi{k}.*dphidx{i}).^2 ...
                     + (phi{i}.*dphidy{k} - phi{k}.*dphidy{i}).^2 ...
                     + (phi{i}.*dphidz{k} - phi{k}.*dphidz{i}).^2 ;
           F_pair(i,k) = vareps(i,k) * sum(sum(sum( grad_term ))) ...
                       + U(i,k) * sum(sum(sum( phi{i}.^2.*phi{k}.^2 )));
           F = F + F_pair(i,k);
       end
end

And the variational derivative of each $F_{ik}$ is computed as,

 for i=1:N_phases,
       for k=1:N_phases,
           if k ~= i
               dFikdphii = 2*U(i,k)*phi{i}.*phi{k}.*(phi{k}-phi{i}) ...
                         + 4*m(i,k)*phi{i}.*phi{k} ...
                         - 2*vareps(i,k)*(phi{k}.*d2phi{i}-phi{i}.*d2phi{k});
               dphidt0{i} = dphidt0{i} - M(i,k) * dFikdphii;
           end
       end
 end

The asymmetric potential can be applied by assigning the parameter m a nonzero value in the above equation.

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

 dphidt{i} = dphidt0{i};

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

d2dphidt0{i} = (dphidt0{i}(pxind,:,:)+dphidt0{i}(nxind,:,:) ...
                          +dphidt0{i}(:,pyind,:)+dphidt0{i}(:,nyind,:) ...
                          +dphidt0{i}(:,:,pzind)+dphidt0{i}(:,:,nzind) ...
                          -6*dphidt0{i})/(h^2);
           dphidt{i} = d2dphidt0{i};

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{i} = dphidt0{i} - mean(mean(mean(dphidt0{i})));

In addition, to calculate the evolution of the whole system, a for loop from 1:N_phases is used to count the changes in all phases.

Matlab test cases

Spherical interface

We can initialize the solid phase to be a sphere,

clear all;
solid_radius = 0.3; equation_type = 0; pf3d_multi_iso

Cylinderical interface

We can also initialize the solid phase to be a cylinder,

clear all;
solid_shape = [0 1 1]; solid_radius = 0.25; equation_type = 0; pf3d_multi_iso

MD++ test cases

In the previous wiki page, we've shown how to download and compile related files in MD++. Here we'd like to show how to run the test cases in MD++.

The following command will give you the single-aniso test case,

bin/phasefield_mc2 scripts/work/phasefield/pf3d.tcl 0 

If you comment out the 28th line of the tcl file, it will give you the single-iso test case.

Running with the command shown below will give you the multi-iso test case.

bin/phasefield_mc2 scripts/work/phasefield/pf3d.tcl 1

The results of MD++ are stored in a series of .cn files, but you can use the matlab file Plot_mdpp_phasefield.m to read and replot the data from MD++.