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 ,where i=1:N and N is the number of different phases in the system. The range of the value of each phase field is between 0 and 1. We can think of the region 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_i=1} as in the ith phase and the region 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_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
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_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}} }
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_{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 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_k=1-\phi_i} , then 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_{ik}} will be reduced to a formalism similar to the single phase model.
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_{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,
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 \frac{\partial \phi_i}{\partial t} = - M_{ik} \sum_{k\neq i}\frac{\delta F_{ik}}{\delta \phi_i} }
(GL)
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 M_{ik}}
is the mobility of i-k interface, 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 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,
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 \frac{\delta F_{ik}}{\delta \phi_i} = 2 U_{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,
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 \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 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_1} , solid 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_2} , and vapor 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_3} , three phase fields by a 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 1\times3} cell. Each element phi{i}(m,n,p) of the cell is a three-dimensional array, 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} . i goes from 1 to 3, m 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} , n 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} , and p goes from 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. Similar to the previous wiki page, by default, periodic boundary conditions (PBC) are applied in all three directions.
For each phase 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_i} , 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_i} , 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_i} , 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_i} , 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_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 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_{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.