Phase Field Model: Friction Force on Interface

Title

Depending of the focus of the paper the title can be different.

• Effective Driving Force Model for the Particle Pinning and Solute Drag Using Phase Field Method : This is a typical title with obscure terms like "effective driving force".
• Modelling of Interface Drag caused by Particle Pinning and Solute Drag in a Unified Model : if we want to emphasis on a generalized model for the interface drag.

Introduction

• Thechnological importance of modelling interface drag phenomena.
• Theories of Particle pinning
• Theories of Solute Drag
• Emphasis on modelling both processes in a single drag model
• Phase Field Models

Controlling motion of interfaces in materials is fundamental aspect of many materials and process designs. As an example in High strength low alloy steels (HSLA) addition of small amount of alloying elements like Nb creates small NbC particles. These particles limit motion of grain boundaries by particle pinning mechanism. In addition, Nb in solute segregates at the grain boundaries and as interface moves it generates solute drag force.

Applying Force on Interface

In polycrystalline materials with dispersion of second phase particles often size of the particles are one to three orders of magnitude smaller than the grains. In these structures resolving of the particles in a discretization grid will be computationally very expensive or impossible. The proposed formulation incorporates Zener pinning force on the grain boundaries as a negative driving force. Depending on the direction of movement of the interface, local energy density shifts to create a negative driving force and decrease speed of the interface.

Shifting of the local energy density of a grain to create artificial force on the interface

The value of Pz can be found from mean field approximations e.g. Zener limiting case or from an independent simulation in smaller scale. When driving pressure for interface movement is smaller than $P_f$ then interface completely stop or mobility coefficient in phase field formulation is set to zero. Systems energy functional can be written as:

$F=\int f_0+\sum _1^p\frac{\kappa }{2}{\left(\mathrm{\nabla }{\eta }_i\right)}^2-\sum _1^{p}6({\eta }_i^2/2-{\eta }_i^3/3)sgn(f_i^{\prime })P_f(v)d^{3}r$

The last term shifts energy of ${\eta }_i$ depending on the sign of $f_i^{\prime }$. The value of $f^{\prime }$ can be interpreted as local driving pressure. Assume an interface with non-zero value of ${\eta }_i$ and ${\eta }_j$. In this formulation if $f_i^{\prime }$ is positive then energy of ${\eta }_i$ is increased by half of the friction force and accordingly energy of ${\eta }_j$ decreased by half of the $P_f$. Note that at interface with equilibrium profile $f_j^{\prime }=-f_i^{\prime }$. Therefore $f^{\prime }$ can be written in following form:

$f_i^{\prime }=\frac{\partial f_0}{\partial {\eta }_i}-\kappa {\mathrm{\nabla }}^2{\eta }_i$

Consequently for interfaces with high curvature the net driving pressure is $f^{\prime }-P_f$. For interfaces with lower curvature a condition occurs that $f^{\prime }<P_f$ therefore the interface stagnates (pinned). In this condition net mobility of the interface is zero. Combination of friction force and variable mobility leads to the evolution equation in the form of

$\frac{\partial {\eta }_i}{\partial t}=-L^{\prime }(\frac{\partial f_0}{\partial {\eta }_i}-{\mathrm{\nabla }}^2{\eta }_i-6{\eta }_i(1-{\eta }_i)sgn(f_i^{\prime })\frac{P_f}{2})$

Where $L^{\prime }=L$ has following simple form:

$L^{\prime }=L|f^{\prime }|>P_f$

$L^{\prime }=0|f^{\prime }|<P_f$

Code Implementation

// gradient square term
del2=delx2*(sumeta2+2*sumeta-24*currenteta);
// driving force part of detadt
detadtM=m*(currenteta*currenteta*currenteta-currenteta)+gamma*(currenteta*sumterm-currenteta*currenteta*currenteta)-kappa*del2;
pzi=3*currenteta*(1-currenteta)*sign(detadtM)*Pz;
// setting mobility to zero if driving force is smaller than pinning force
if (fabs(detadtM)<fabs(pzi)){
    Lf=0;
}
else{
    Lf=L;
}
//actual detadt which is substraction of driving force and pinning components
detadt=-Lf*(detadtM-pzi); 
// one time step increment
eta2[i+jn+kn+pnn]=currenteta+delt*detadt; 

Variable Friction Force

If $P_f$ is function of velocity then one can find velocity of the interface and relate force to the velocity by a formula (e.g. Cahn theory of solute drag):

$P_f=\frac{\alpha v}{1+{\beta }^2{v}^2}$

Finding velocity field at each point is possible using by-part derivation:

$\frac{\partial \eta }{\partial t}=-L{f}^{\prime }=\frac{\partial \eta }{\partial x}\frac{\partial x}{\partial t}$

Therefore generalization to 3D gives:

$v=\frac{-L{f}^{\prime }}{\mathrm{\nabla }\eta }$

Application to Interface Drag Phenomenon

In general one can define dependence of friction force to the velocity an in two particular case one can obtain Zener pinning and solute drag phenomena separately. For particle pinning the friction force is constant (i.e. pinning force) and changes its direction depending on the sign of the velocity. For solute drag however the friction force functionality is more complex depending on the velocity of the interface.

File:Friction Force Velocity.png

Application in Particle Pinning

In this section the friction force methodology is applied on particle pinning. Classical Zener pinning considers inert particles with rarius $P_r$ which are uniformly distributed in the grain structure. The first step is to obtain the friction force (pressure) on the interface. One can obtain the drag force by any simulation tecnique. In here, we use a phase field model with driving force and second phase particle. The geometry consists of a flat interface moving under constant driving pressure within an array of particles in three dimenssion.

File:Nparticle 3D 100.png

File:Nparticle 3D 2800.png

An interface inteacting with an ensumble of particles.

Now that the drag force is known we can insert the value to the Eq. 3 adn simulate microstructures. However, for testing purpose a dome-shape grain was considered first. In this geometry curvature of the the upper side of the grain remains constant and provide constant driving pressure. Volume of a dome-shape grain vs time decreases linearly with different speed depending on the value of the $P_f$:

File:U grain shrinks.png

Volume of a U shape grain is plotted vs time for different values of $P_f$. If $P_f$ is higher than deriving pressure of interface motion interface completely stops.

%-------------------------------------------------------------------------------------------------------------------- \pagebreak

\subsection{Simulation of 2D Grain Growth}

\begin{figure}[h]

\centering
\includegraphics[width=\textwidth]{Grains_Snapshot}

\end{figure} Simulation of grain growth for $P_f=0.01$ until saturation of grain size. Domain size is $500\times500$ grid points.

%--------------------------------------------------------------------------------------------------------------------

File:Grain A Saturates.png

Kinetics of grain growth for different value of $P_f$. Y-axis is plotted in logarithmic scale. Normal grain growth ($P_f=0$) is shown in red. In practical cases volume fraction of second phase particles can vary from 0.001 to 0.1 causing large variation in $P_f$.

File:Nparticle Pf0.1 Pr5 DelG.png

Change in volume of a grain as interface moves through array of particles. Each curves represent a simulation starting with a completely flat interface. Driving force in varried from 0.3 to 0.015.

Application in Solute Drag

• Use a 1D phase field model to get the function
• show the friction force model generates the same velocity-driving force dependent function.
• Microstructural simulation shows a reduced kinetics for grain growth or recrystallization.

Appendix 1. Derivation of Weighting function

In here we prove that the weighting function in the functional equation reproduce the fundamental relationship of interface velocity and driving pressure ($v=M(\mathrm{\Delta }G-Pz)$).

Assuming the interface is moving forward therefore $sgn$ function is positive, the evolution equation is written as the following:

$\frac{\partial \eta }{\partial t}=-L\frac{\delta F}{\delta \eta }=-L[m\frac{\partial f_0}{\partial \eta }-\kappa {\mathrm{\nabla }}^2\eta +g^{\prime }(\eta )P]$ ......................... (A1)

The gradient squared term in curvilinear coordinate system can be rearranged to the gradient of a vector $\mathbf{𝐫}$ normal to iso-$\eta$ profiles. $\mathbf{𝐫}$ is also parallel to the direction of interface velocity. In curvilinear system $\mathrm{\nabla }\eta =\partial \eta \partial r\mathbf{𝐫}$ and ${\mathrm{\nabla }}^2\eta =\mathrm{\nabla }\cdot \mathrm{\nabla }\eta =({\partial }^2\eta /\partial r^2+\mathrm{\nabla }\cdot \mathbf{𝐫}(\partial \eta /\partial r))$. Therefore Eq. (A1) can be written as following:

$\frac{\partial \eta }{\partial t}=-L[m\frac{\partial f_0}{\partial \eta }-\kappa (\frac{{\partial }^2\eta }{\partial r^2}+(\mathrm{\nabla }\cdot \mathbf{𝐫})\frac{\partial \eta }{\partial r})+g^{\prime }(\eta )P]$ ......................... (A2)

The term $\mathrm{\nabla }\cdot \mathbf{𝐫}$ represents curvature of the interface [?].The weighting function $g^{\prime }(\eta )$ is choosen to satisfy two conditions: (i) application of the driving force ($P\ne 0$) does not change the shape of the interface profile, (ii) relationship between the interface velocity and the driving force is linear ($v=MP$ where $M$ is the interface interinsic mobility).

A straight intetrface ($\mathrm{\nabla }\cdot \mathbf{𝐫}=0$) with no driving force remains stationary, i.e. $\partial \eta /\partial t=0$ therfore following equality can be obtained from equation (A2):

$0=\frac{\partial \eta }{\partial t}=-L[m\frac{\partial f_0}{\partial \eta }-\kappa (\frac{{\partial }^2\eta }{\partial r^2}+0)]\to m\frac{\partial f_0}{\partial \eta }=\kappa \frac{{\partial }^2\eta }{\partial r^2}$ ......................... (A3)

We will later show that applying $P$ will not change order parameter profile. Equality of Eq. (A3) and using ($\partial \eta /\partial t=\partial \eta /\partial r\cdot \partial r/\partial t$ )simplifies Eq (A2) to the following equation:

$v=\frac{\partial r}{\partial \eta }=-L[-\kappa \mathrm{\nabla }\cdot \mathbf{𝐫}+g^{\prime }(\eta )P\frac{\partial r}{\partial \eta }]$ ......................... (A4)

where $v$ is velocity of the intetrface. In sharp interface limit

$v=M(\sigma /R-P)$ ............................(A5)

where $1/R$ is the mean curvature of the interface and $\sigma$ is the interfacial energy. As Moelans et. al. showed $L\kappa =\sigma$. By combining Eqs. A4 and A5 $g^{\prime }$ can be obtained:

$g^{\prime }(\eta )=-M/L\frac{\partial \eta }{\partial r}$ ...........................(A6)

If the interface profile doesn't change when an additional driving force is applied then we can use the profile obtained for stationary straight interface. The interface profile can be obtained by solving Eq. A3 for two order parameters. The solution is given in Ref [?]:

$\eta =\frac{1}{2}[1-tanh\left(\sqrt{\frac{m}{2\kappa }}r\right)]$ ...........................(A7)

The derivative $\partial \eta /\partial r$ as a function of $\eta$ is

$\frac{\partial \eta }{\partial r}=-2\sqrt{\frac{m}{2\kappa }}(\eta -{\eta }^2)$ ...........................(A8)

The relationship between phase field model parameters and physical properties of the interface was derived in Ref [1]. Accordingly, the interface mobility is given by

$M=3L\sqrt{\frac{2\kappa }{m}}$ ...........................(A9)

Consequently, by substituting Eqs. A8 and A9 into Eq. A6 the weighting function can be obtained.

$g^{\prime }(\eta )=6(\eta -{\eta }^2)$ ...........................(A10)

The wiething function $g(\eta )$ in the functional equation (Eq. ?) is obtained by integration of $g^{\prime }$.

$g(\eta )=\int g^{\prime }(\eta )d\eta =6(\frac{{\eta }^2}{2}-\frac{{\eta }^3}{3})$ ...........................(A11)

Profile of a moving interface under driving pressure

We show that addition of driving or drag pressure does not change interface profile shape. Profile of a stationary interface is given in Eq. A7. If interface shape remains constant while moving then following equation should satisfy the phase field PDE.

$\eta =\frac{1}{2}[1-tanh(\sqrt{\frac{m}{2\kappa }}r-vt)]$ ..................................(A12)

where $v$ is velocity of the interface and $t$ is time.