examples.phase package¶

Subpackages¶

Submodules¶

examples.phase.anisotropy module¶

Solve a dendritic solidification problem.

To convert a liquid material to a solid, it must be cooled to a temperature below its melting point (known as “undercooling” or “supercooling”). The rate of solidification is often assumed (and experimentally found) to be proportional to the undercooling. Under the right circumstances, the solidification front can become unstable, leading to dendritic patterns. Warren, Kobayashi, Lobkovsky and Carter [13] have described a phase field model (“Allen-Cahn”, “non-conserved Ginsberg-Landau”, or “model A” of Hohenberg & Halperin) of such a system, including the effects of discrete crystalline orientations (anisotropy).

We start with a regular 2D Cartesian mesh

>>> from fipy import Variable, CellVariable, Grid2D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, Viewer, Matplotlib2DGridViewer
>>> from fipy.tools import numerix
>>> dx = dy = 0.025
>>> if __name__ == '__main__':
...     nx = ny = 500
... else:
...     nx = ny = 20
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

and we’ll take fixed timesteps

>>> dt = 5e-4

We consider the simultaneous evolution of a “phase field” variable $\phi$ (taken to be 0 in the liquid phase and 1 in the solid)

>>> phase = CellVariable(name=r'$\phi$', mesh=mesh, hasOld=True)

and a dimensionless undercooling $\Delta T$ ( $\Delta T = 0$ at the melting point)

>>> dT = CellVariable(name=r'$\Delta T$', mesh=mesh, hasOld=True)

The hasOld flag causes the storage of the value of variable from the previous timestep. This is necessary for solving equations with non-linear coefficients or for coupling between PDEs.

The governing equation for the temperature field is the heat flux equation, with a source due to the latent heat of solidification

$\frac{\partial \Delta T}{\partial t} = D_T \nabla^2 \Delta T + \frac{\partial \phi}{\partial t}$

>>> DT = 2.25
>>> heatEq = (TransientTerm()
...           == DiffusionTerm(DT)
...           + (phase - phase.old) / dt)

The governing equation for the phase field is

$\tau_{\phi} \frac{\partial \phi}{\partial t} = \nabla \cdot \mathsf{D} \nabla \phi + \phi ( 1 - \phi ) m ( \phi , \Delta T)$

where

$m(\phi, \Delta T) = \phi - \frac{1}{2} - \frac{ \kappa_1 }{ \pi } \arctan \left( \kappa_2 \Delta T \right)$

represents a source of anisotropy. The coefficient $\mathsf{D}$ is an anisotropic diffusion tensor in two dimensions

$\mathsf{D} = \alpha^2 \left( 1 + c \beta \right) \left[ \begin{matrix} 1 + c \beta & -c \frac{\partial \beta}{\partial \psi} \\ c \frac{\partial \beta}{\partial \psi} & 1 + c \beta \end{matrix} \right]$

where $\beta = \frac{ 1 - \Phi^2 } { 1 + \Phi^2}$ , $\Phi = \tan \left( \frac{ N } { 2 } \psi \right)$ , $\psi = \theta + \arctan \frac{\partial \phi / \partial y}{\partial \phi / \partial x}$ , $\theta$ is the orientation, and $N$ is the symmetry.

>>> alpha = 0.015
>>> c = 0.02
>>> N = 6.
>>> theta = numerix.pi / 8.
>>> psi = theta + numerix.arctan2(phase.faceGrad[1],
...                               phase.faceGrad[0])
>>> Phi = numerix.tan(N * psi / 2)
>>> PhiSq = Phi**2
>>> beta = (1. - PhiSq) / (1. + PhiSq)
>>> DbetaDpsi = -N * 2 * Phi / (1 + PhiSq)
>>> Ddia = (1.+ c * beta)
>>> Doff = c * DbetaDpsi
>>> I0 = Variable(value=((1, 0), (0, 1)))
>>> I1 = Variable(value=((0, -1), (1, 0)))
>>> D = alpha**2 * (1.+ c * beta) * (Ddia * I0 + Doff * I1)

With these expressions defined, we can construct the phase field equation as

>>> tau = 3e-4
>>> kappa1 = 0.9
>>> kappa2 = 20.
>>> phaseEq = (TransientTerm(tau)
...            == DiffusionTerm(D)
...            + ImplicitSourceTerm((phase - 0.5 - kappa1 / numerix.pi * numerix.arctan(kappa2 * dT))
...                                 * (1 - phase)))

We seed a circular solidified region in the center

>>> radius = dx * 5.
>>> C = (nx * dx / 2, ny * dy / 2)
>>> x, y = mesh.cellCenters
>>> phase.setValue(1., where=((x - C[0])**2 + (y - C[1])**2) < radius**2)

and quench the entire simulation domain below the melting point

>>> dT.setValue(-0.5)

In a real solidification process, dendritic branching is induced by small thermal fluctuations along an otherwise smooth surface, but the granularity of the Mesh is enough “noise” in this case, so we don’t need to explicitly introduce randomness, the way we did in the Cahn-Hilliard problem.

FiPy’s viewers are utilitarian, striving to let the user see something, regardless of their operating system or installed packages, so you won’t be able to simultaneously view two fields “out of the box”, but, because all of Python is accessible and FiPy is object oriented, it is not hard to adapt one of the existing viewers to create a specialized display:

>>> if __name__ == "__main__":
...     try:
...         import pylab
...         class DendriteViewer(Matplotlib2DGridViewer):
...             def __init__(self, phase, dT, title=None, limits={}, **kwlimits):
...                 self.phase = phase
...                 self.contour = None
...                 Matplotlib2DGridViewer.__init__(self, vars=(dT,), title=title,
...                                                 cmap=pylab.cm.hot,
...                                                 limits=limits, **kwlimits)
...
...             def _plot(self):
...                 Matplotlib2DGridViewer._plot(self)
...
...                 if self.contour is not None:
...                     for c in self.contour.collections:
...                         c.remove()
...
...                 mesh = self.phase.mesh
...                 shape = mesh.shape
...                 x, y = mesh.cellCenters
...                 z = self.phase.value
...                 x, y, z = [a.reshape(shape, order='F') for a in (x, y, z)]
...
...                 self.contour = self.axes.contour(x, y, z, (0.5,))
...
...         viewer = DendriteViewer(phase=phase, dT=dT,
...                                 title=r"%s & %s" % (phase.name, dT.name),
...                                 datamin=-0.1, datamax=0.05)
...     except ImportError:
...         viewer = MultiViewer(viewers=(Viewer(vars=phase),
...                                       Viewer(vars=dT,
...                                              datamin=-0.5,
...                                              datamax=0.5)))

and iterate the solution in time, plotting as we go,

>>> if __name__ == '__main__':
...     steps = 10000
... else:
...     steps = 10
>>> from builtins import range
>>> for i in range(steps):
...     phase.updateOld()
...     dT.updateOld()
...     phaseEq.solve(phase, dt=dt)
...     heatEq.solve(dT, dt=dt)
...     if __name__ == "__main__" and (i % 10 == 0):
...         viewer.plot()

phase field and undercooling during solidification of a 6-fold "snowflake" anisotropic seed

The non-uniform temperature results from the release of latent heat at the solidifying interface. The dendrite arms grow fastest where the temperature gradient is steepest.

We note that this FiPy simulation is written in about 50 lines of code (excluding the custom viewer), compared with over 800 lines of (fairly lucid) FORTRAN code used for the figures in [13].

examples.phase.anisotropyOLD module¶

Attention

This example remains only for exact comparison against Ryo Kobayashi’s FORTRAN code. See examples.phase.anisotropy for a better, although not numerically identical implementation.

In this example we solve a coupled phase and temperature equation to model solidification, and eventually dendritic growth, based on the work of Warren, Kobayashi, Lobkovsky and Carter [13].

We start from a circular seed in a 2D mesh:

>>> from fipy import CellVariable, Grid2D, TransientTerm, DiffusionTerm, ExplicitDiffusionTerm, ImplicitSourceTerm, Viewer
>>> from fipy.tools import numerix

>>> numberOfCells = 40
>>> Length = numberOfCells * 2.5 / 100.
>>> nx = numberOfCells
>>> ny = numberOfCells
>>> dx = Length / nx
>>> dy = Length / ny
>>> radius = Length / 4.
>>> seedCenter = (Length / 2., Length / 2.)
>>> initialTemperature = -0.4
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

Dendritic growth will not be observed with this small test system. If you wish to see dendritic growth reset the following parameters such that numberOfCells = 500, steps = 10000, radius = dx * 5. seedCenter = (0. , 0.) and initialTemperature = -0.5.

The governing equation for the phase field is given by:

$\tau_{\phi} \frac{\partial \phi}{\partial t} = \nabla \cdot \left[ D \nabla \phi + A \nabla \xi \right] + \phi ( 1 - \phi ) m ( \phi , T)$

where

$m(\phi, T) = \phi - \frac{1}{2} - \frac{ \kappa_1 }{ \pi } \arctan \left( \kappa_2 T \right).$

The coefficients $D$ and $A$ are given by,

$D = \alpha^2 \left[ 1 + c \beta \right]^2$

and

$A = \alpha^2 c \left[ 1 + c \beta \right] \beta_\psi$

where $\beta = \frac{ 1 - \Phi^2 } { 1 + \Phi^2}$ , $\Phi = \tan \left( \frac{ N } { 2 } \psi \right) `, :math:$ psi = theta + arctan frac{ phi_y } { phi_x }` and $\xi_x = -\phi_y$ and $\xi_y = \phi_x$ .

The governing equation for temperature is given by:

$\frac{\partial T}{\partial t} = D_T \nabla^2 T + \frac{\partial \phi}{\partial t}$

Here the phase and temperature equations are solved with an explicit and implicit technique, respectively.

The parameters for these equations are

>>> timeStepDuration = 5e-5
>>> tau = 3e-4
>>> alpha = 0.015
>>> c = 0.02
>>> N = 4.
>>> kappa1 = 0.9
>>> kappa2 = 20.
>>> tempDiffusionCoeff = 2.25
>>> theta = 0.

The phase variable is 0 for a liquid and 1 for a solid. Here, the ``phase variable is initialized as a liquid,

>>> phase = CellVariable(name='phase field', mesh=mesh, hasOld=1)

The hasOld flag keeps the old value of the variable. This is necessary for a transient solution. In this example we wish to set up an interior region that is solid. The domain is seeded with a circular solidified region with parameters seedCenter and radius representing the center and radius of the seed.

>>> x, y = mesh.cellCenters
>>> phase.setValue(1., where=((x - seedCenter[0])**2
...                           + (y - seedCenter[1])**2) < radius**2)

The temperature field is initialized to a value of -0.4 throughout:

>>> temperature = CellVariable(
...     name='temperature',
...     mesh=mesh,
...     value=initialTemperature,
...     hasOld=1)

The $m(\phi, T)$ variable

is created from the phase and temperature variables.

>>> mVar = phase - 0.5 - kappa1 / numerix.pi * numerix.arctan(kappa2 * temperature)

The following section of code builds up the $A$ and $D$ coefficients.

>>> phaseY = phase.faceGrad.dot((0, 1))
>>> phaseX = phase.faceGrad.dot((1, 0))
>>> psi = theta + numerix.arctan2(phaseY, phaseX)
>>> Phi = numerix.tan(N * psi / 2)
>>> PhiSq = Phi**2
>>> beta = (1. - PhiSq) / (1. + PhiSq)
>>> betaPsi = -N * 2 * Phi / (1 + PhiSq)
>>> A = alpha**2 * c * (1.+ c * beta) * betaPsi
>>> D = alpha**2 * (1.+ c * beta)**2

The $\nabla \xi$ variable (dxi), given by $(\xi_x, \xi_y) = (-\phi_y, \phi_x)$ , is constructed by first obtaining $\nabla \phi$

using getFaceGrad(). The axes are rotated ninety degrees.

>>> dxi = phase.faceGrad.dot(((0, 1), (-1, 0)))
>>> anisotropySource = (A * dxi).divergence

The phase equation can now be constructed.

>>> phaseEq = TransientTerm(tau) == ExplicitDiffusionTerm(D) + \
...     ImplicitSourceTerm(mVar * ((mVar < 0) - phase)) + \
...     ((mVar > 0.) * mVar * phase + anisotropySource)

The temperature equation is built in the following way,

>>> temperatureEq = TransientTerm() == \
...                 DiffusionTerm(tempDiffusionCoeff) + \
...                 (phase - phase.old) / timeStepDuration

If we are running this example interactively, we create viewers for the phase and temperature fields

>>> if __name__ == '__main__':
...     phaseViewer = Viewer(vars=phase)
...     temperatureViewer = Viewer(vars=temperature,
...                                datamin=-0.5, datamax=0.5)
...     phaseViewer.plot()
...     temperatureViewer.plot()

we iterate the solution in time, plotting as we go if running interactively,

>>> steps = 10
>>> from builtins import range
>>> for i in range(steps):
...     phase.updateOld()
...     temperature.updateOld()
...     phaseEq.solve(phase, dt=timeStepDuration)
...     temperatureEq.solve(temperature, dt=timeStepDuration)
...     if i%10 == 0 and __name__ == '__main__':
...         phaseViewer.plot()
...         temperatureViewer.plot()

The solution is compared with test data. The test data was created for steps = 10 with a FORTRAN code written by Ryo Kobayashi for phase field modeling. The following code opens the file anisotropy.gz extracts the data and compares it with the phase variable.

>>> import os
>>> from future.utils import text_to_native_str
>>> testData = numerix.loadtxt(os.path.splitext(__file__)[0] + text_to_native_str('.gz'))
>>> print(phase.allclose(testData))
1

examples.phase.binary module¶

It is straightforward to extend a phase field model to include binary alloys. As in examples.phase.simple, we will examine a 1D problem

>>> from fipy import CellVariable, Variable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, PowerLawConvectionTerm, DefaultAsymmetricSolver, Viewer
>>> from fipy.tools import numerix

>>> nx = 400
>>> dx = 5e-6 # cm
>>> L = nx * dx
>>> mesh = Grid1D(dx=dx, nx=nx)

The Helmholtz free energy functional can be written as the integral [1] [3] [24]

$\mathcal{F}\left(\phi, C, T\right) = \int_\mathcal{V} \left\{ f(\phi, C, T) + \frac{\kappa_\phi}{2}\abs{\nabla\phi}^2 + \frac{\kappa_C}{2}\abs{\nabla C}^2 \right\} dV$

over the volume $\mathcal{V}$ as a function of phase $\phi$ [#phi]_

>>> phase = CellVariable(name="phase", mesh=mesh, hasOld=1)

composition $C$

>>> C = CellVariable(name="composition", mesh=mesh, hasOld=1)

and temperature $T$ [#T]_

>>> T = Variable(name="temperature")

Frequently, the gradient energy term in concentration is ignored and we can derive governing equations

(1)¶ $\frac{\partial\phi}{\partial t} = M_\phi \left( \kappa_\phi \nabla^2 \phi - \frac{\partial f}{\partial \phi} \right)$

for phase and

(2)¶ $\frac{\partial C}{\partial t} = \nabla\cdot\left( M_C \nabla \frac{\partial f}{\partial C} \right)$

for solute.

The free energy density $f(\phi, C, T)$ can be constructed in many different ways. One approach is to construct free energy densities for each of the pure components, as functions of phase, e.g.

$f_A(\phi, T) = p(\phi) f_A^S(T) + \left(1 - p(\phi)\right) f_A^L(T) + \frac{W_A}{2} g(\phi)$

where $f_A^L(T)$ , $f_B^L(T)$ , $f_A^S(T)$ , and $f_B^S(T)$ are the free energy densities of the pure components. There are a variety of choices for the interpolation function $p(\phi)$ and the barrier function $g(\phi)$ ,

such as those shown in examples.phase.simple

>>> def p(phi):
...     return phi**3 * (6 * phi**2 - 15 * phi + 10)

>>> def g(phi):
...     return (phi * (1 - phi))**2

The desired thermodynamic model can then be applied to obtain $f(\phi, C, T)$ , such as for a regular solution,

$f(\phi, C, T) &= (1 - C) f_A(\phi, T) + C f_B(\phi, T) \\ &\qquad + R T \left[ (1 - C) \ln (1 - C) + C \ln C \right] + C (1 - C) \left[ \Omega_S p(\phi) + \Omega_L \left( 1 - p(\phi) \right) \right]$

where

>>> R = 8.314 # J / (mol K)

is the gas constant and $\Omega_S$ and $\Omega_L$ are the regular solution interaction parameters for solid and liquid.

Another approach is useful when the free energy densities $f^L(C, T)$ and $f^S(C,T)$ of the alloy in the solid and liquid phases are known. This might be the case when the two different phases have different thermodynamic models or when one or both is obtained from a Calphad code. In this case, we can construct

$f(\phi, C, T) = p(\phi) f^S(C,T) + \left(1 - p(\phi)\right) f^L(C, T) + \left[ (1-C) \frac{W_A}{2} + C \frac{W_B}{2} \right] g(\phi).$

When the thermodynamic models are the same in both phases, both approaches should yield the same result.

We choose the first approach and make the simplifying assumptions of an ideal solution and that

$f_A^L(T) & = 0 \\ f_A^S(T) - f_A^L(T) &= \frac{L_A\left(T - T_M^A\right)}{T_M^A}$

and likewise for component $B$ .

>>> LA = 2350. # J / cm**3
>>> LB = 1728. # J / cm**3
>>> TmA = 1728. # K
>>> TmB = 1358. # K

>>> enthalpyA = LA * (T - TmA) / TmA
>>> enthalpyB = LB * (T - TmB) / TmB

This relates the difference between the free energy densities of the pure solid and pure liquid phases to the latent heat $L_A$ and the pure component melting point $T_M^A$ , such that

$f_A(\phi, T) = \frac{L_A\left(T - T_M^A\right)}{T_M^A} p(\phi) + \frac{W_A}{2} g(\phi).$

With these assumptions

(3)¶ $\frac{\partial f}{\partial \phi} &= (1-C) \frac{\partial f_A}{\partial \phi} + C \frac{\partial f_B}{\partial \phi} \nonumber \\ &= \left\{ (1-C) \frac{L_A\left(T - T_M^A\right)}{T_M^A} + C \frac{L_B\left(T - T_M^B\right)}{T_M^B} \right\} p'(\phi) + \left\{ (1-C) \frac{W_A}{2} + C \frac{W_B}{2} \right\} g'(\phi)$

and

(4)¶ $\frac{\partial f}{\partial C} &= \left[f_B(\phi, T) + \frac{R T}{V_m} \ln C\right] - \left[f_A(\phi, T) + \frac{R T}{V_m} \ln (1-C) \right] \nonumber \\ &= \left[\mu_B(\phi, C, T) - \mu_A(\phi, C, T) \right] / V_m$

where $\mu_A$ and $\mu_B$ are the classical chemical potentials for the binary species. $p'(\phi)$ and $g'(\phi)$ are the partial derivatives of of $p$ and $g$ with respect to $\phi$

>>> def pPrime(phi):
...     return 30. * g(phi)

>>> def gPrime(phi):
...     return 2. * phi * (1 - phi) * (1 - 2 * phi)

$V_m$ is the molar volume, which we take to be independent of concentration and phase

>>> Vm = 7.42 # cm**3 / mol

On comparison with examples.phase.simple, we can see that the present form of the phase field equation is identical to the one found earlier, with the source now composed of the concentration-weighted average of the source for either pure component. We let the pure component barriers equal the previous value

>>> deltaA = deltaB = 1.5 * dx
>>> sigmaA = 3.7e-5 # J / cm**2
>>> sigmaB = 2.9e-5 # J / cm**2
>>> betaA = 0.33 # cm / (K s)
>>> betaB = 0.39 # cm / (K s)
>>> kappaA = 6 * sigmaA * deltaA # J / cm
>>> kappaB = 6 * sigmaB * deltaB # J / cm
>>> WA = 6 * sigmaA / deltaA # J / cm**3
>>> WB = 6 * sigmaB / deltaB # J / cm**3

and define the averages

>>> W = (1 - C) * WA / 2. + C * WB / 2.
>>> enthalpy = (1 - C) * enthalpyA + C * enthalpyB

We can now linearize the source exactly as before

>>> mPhi = -((1 - 2 * phase) * W + 30 * phase * (1 - phase) * enthalpy)
>>> dmPhidPhi = 2 * W - 30 * (1 - 2 * phase) * enthalpy
>>> S1 = dmPhidPhi * phase * (1 - phase) + mPhi * (1 - 2 * phase)
>>> S0 = mPhi * phase * (1 - phase) - S1 * phase

Using the same gradient energy coefficient and phase field mobility

>>> kappa = (1 - C) * kappaA + C * kappaB
>>> Mphi = TmA * betaA / (6 * LA * deltaA)

we define the phase field equation

>>> phaseEq = (TransientTerm(1/Mphi) == DiffusionTerm(coeff=kappa)
...            + S0 + ImplicitSourceTerm(coeff=S1))

When coding explicitly, it is typical to simply write a function to evaluate the chemical potentials $\mu_A$ and $\mu_B$ and then perform the finite differences necessary to calculate their gradient and divergence, e.g.,:

def deltaChemPot(phase, C, T):
    return ((Vm * (enthalpyB * p(phase) + WA * g(phase)) + R * T * log(1 - C)) -
            (Vm * (enthalpyA * p(phase) + WA * g(phase)) + R * T * log(C)))

for j in range(faces):
    flux[j] = ((Mc[j+.5] + Mc[j-.5]) / 2)           * (deltaChemPot(phase[j+.5], C[j+.5], T)             - deltaChemPot(phase[j-.5], C[j-.5], T)) / dx

for j in range(cells):
    diffusion = (flux[j+.5] - flux[j-.5]) / dx

where we neglect the details of the outer boundaries (j = 0 and j = N) or exactly how to translate j+.5 or j-.5 into an array index, much less the complexities of higher dimensions. FiPy can handle all of these issues automatically, so we could just write:

chemPotA = Vm * (enthalpyA * p(phase) + WA * g(phase)) + R * T * log(C)
chemPotB = Vm * (enthalpyB * p(phase) + WB * g(phase)) + R * T * log(1-C)
flux = Mc * (chemPotB - chemPotA).faceGrad
eq = TransientTerm() == flux.divergence

Although the second syntax would essentially work as written, such an explicit implementation would be very slow. In order to take advantage of FiPy’s implicit solvers, it is necessary to reduce Eq. (2) to the canonical form of Eq. (2), hence we must expand Eq. (4) as

$\frac{\partial f}{\partial C} = \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] p(\phi) + \frac{R T}{V_m} \left[\ln C - \ln (1-C)\right] + \frac{W_B - W_A}{2} g(\phi)$

In either bulk phase, $\nabla p(\phi) = \nabla g(\phi) = 0$ , so we can then reduce Eq. (2) to

(5)¶ $\frac{\partial C}{\partial t} &= \nabla\cdot\left( M_C \nabla \left\{ \frac{R T}{V_m} \left[\ln C - \ln (1-C)\right] \right\} \right) \nonumber \\ &= \nabla\cdot\left[ \frac{M_C R T}{C (1-C) V_m} \nabla C \right]$

and, by comparison with Fick’s second law

$\frac{\partial C}{\partial t} = \nabla\cdot\left[D \nabla C\right],$

we can associate the mobility $M_C$ with the intrinsic diffusivity $D$ by $M_C \equiv D C (1-C) V_m / R T$ and write Eq. (2) as

(6)¶ $\frac{\partial C}{\partial t} &= \nabla\cdot\left( D \nabla C \right) \nonumber \\ &\qquad + \nabla\cdot\left( \frac{D C (1 - C) V_m}{R T} \left\{ \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] \nabla p(\phi) + \frac{W_B - W_A}{2} \nabla g(\phi) \right\} \right).$

The first term is clearly a DiffusionTerm. The second is less obvious, but by factoring out $C$ , we can see that this is a ConvectionTerm with a velocity

$\vec{u}_\phi = \frac{D (1 - C) V_m}{R T} \left\{ \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] \nabla p(\phi) + \frac{W_B - W_A}{2} \nabla g(\phi) \right\}$

due to phase transformation, such that

$\frac{\partial C}{\partial t} = \nabla\cdot\left( D \nabla C \right) + \nabla\cdot\left(C \vec{u}_\phi\right)$

or

>>> Dl = Variable(value=1e-5) # cm**2 / s
>>> Ds = Variable(value=1e-9) # cm**2 / s
>>> D = (Ds - Dl) * phase.arithmeticFaceValue + Dl

>>> phaseTransformationVelocity = (((enthalpyB - enthalpyA) * p(phase).faceGrad
...                                + 0.5 * (WB - WA) * g(phase).faceGrad)
...                                * D * (1. - C).harmonicFaceValue * Vm / (R * T))

>>> diffusionEq = (TransientTerm()
...                == DiffusionTerm(coeff=D)
...                + PowerLawConvectionTerm(coeff=phaseTransformationVelocity))

We initialize the phase field to a step function in the middle of the domain

>>> phase.setValue(1.)
>>> phase.setValue(0., where=mesh.cellCenters[0] > L/2.)

and start with a uniform composition field $C = 1/2$

>>> C.setValue(0.5)

In equilibrium, $\mu_A(0, C_L, T) = \mu_A(1, C_S, T)$ and $\mu_B(0, C_L, T) = \mu_B(1, C_S, T)$ and, for ideal solutions, we can deduce the liquidus and solidus compositions as

$C_L &= \frac{1 - \exp\left(-\frac{L_A\left(T - T_M^A\right)}{T_M^A}\frac{V_m}{R T}\right)} {\exp\left(-\frac{L_B\left(T - T_M^B\right)}{T_M^B}\frac{V_m}{R T}\right) - \exp\left(-\frac{L_A\left(T - T_M^A\right)}{T_M^A}\frac{V_m}{R T}\right)} \\ C_S &= \exp\left(-\frac{L_B\left(T - T_M^B\right)}{T_M^B}\frac{V_m}{R T}\right) C_L$

>>> Cl = ((1. - numerix.exp(-enthalpyA * Vm / (R * T)))
...       / (numerix.exp(-enthalpyB * Vm / (R * T)) - numerix.exp(-enthalpyA * Vm / (R * T))))
>>> Cs = numerix.exp(-enthalpyB * Vm / (R * T)) * Cl

The phase fraction is predicted by the lever rule

>>> Cavg = C.cellVolumeAverage
>>> fraction = (Cl - Cavg) / (Cl - Cs)

For the special case of fraction = Cavg = 0.5, a little bit of algebra reveals that the temperature that leaves the phase fraction unchanged is given by

>>> T.setValue((LA + LB) * TmA * TmB / (LA * TmB + LB * TmA))

In this simple, binary, ideal solution case, we can derive explicit expressions for the solidus and liquidus compositions. In general, this may not be possible or practical. In that event, the root-finding facilities in SciPy can be used.

We’ll need a function to return the two conditions for equilibrium

$0 = \mu_A(1, C_S, T) - \mu_A(0, C_L, T) &= \frac{L_A\left(T - T_M^A\right)}{T_M^A} V_m + R T \ln (1 - C_S) - R T \ln (1 - C_L) \\ 0 = \mu_B(1, C_S, T) - \mu_B(0, C_L, T) &= \frac{L_B\left(T - T_M^B\right)}{T_M^B} V_m + R T \ln C_S - R T \ln C_L$

>>> def equilibrium(C):
...     return [numerix.array(enthalpyA * Vm
...                           + R * T * numerix.log(1 - C[0])
...                           - R * T * numerix.log(1 - C[1])),
...             numerix.array(enthalpyB * Vm
...                           + R * T * numerix.log(C[0])
...                           - R * T * numerix.log(C[1]))]

and we’ll have much better luck if we also supply the Jacobian

$\left[\begin{matrix} \frac{\partial(\mu_A^S - \mu_A^L)}{\partial C_S} & \frac{\partial(\mu_A^S - \mu_A^L)}{\partial C_L} \\ \frac{\partial(\mu_B^S - \mu_B^L)}{\partial C_S} & \frac{\partial(\mu_B^S - \mu_B^L)}{\partial C_L} \end{matrix}\right] = R T\left[\begin{matrix} -\frac{1}{1-C_S} & \frac{1}{1-C_L} \\ \frac{1}{C_S} & -\frac{1}{C_L} \end{matrix}\right]$

>>> def equilibriumJacobian(C):
...     return R * T * numerix.array([[-1. / (1 - C[0]), 1. / (1 - C[1])],
...                                   [ 1. / C[0],      -1. / C[1]]])

>>> try:
...     from scipy.optimize import fsolve 
...     CsRoot, ClRoot = fsolve(func=equilibrium, x0=[0.5, 0.5],
...                             fprime=equilibriumJacobian) 
... except ImportError:
...     ClRoot = CsRoot = 0
...     print("The SciPy library is not available to calculate the solidus and ... liquidus concentrations")

>>> print(Cl.allclose(ClRoot)) 
1
>>> print(Cs.allclose(CsRoot)) 
1

We plot the result against the sharp interface solution

>>> sharp = CellVariable(name="sharp", mesh=mesh)
>>> x = mesh.cellCenters[0]
>>> sharp.setValue(Cs, where=x < L * fraction)
>>> sharp.setValue(Cl, where=x >= L * fraction)

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=(phase, C, sharp),
...                     datamin=0., datamax=1.)
...     viewer.plot()

Because the phase field interface will not move, and because we’ve seen in earlier examples that the diffusion problem is unconditionally stable, we need take only one very large timestep to reach equilibrium

>>> dt = 1.e5

Because the phase field equation is coupled to the composition through enthalpy and W and the diffusion equation is coupled to the phase field through phaseTransformationVelocity, it is necessary sweep this non-linear problem to convergence. We use the “residual” of the equations (a measure of how well they think they have solved the given set of linear equations) as a test for how long to sweep. Because of the ConvectionTerm, the solution matrix for diffusionEq is asymmetric and cannot be solved by the default LinearPCGSolver. Therefore, we use a DefaultAsymmetricSolver for this equation.

We now use the “sweep()” method instead of “solve()” because we require the residual.

>>> solver = DefaultAsymmetricSolver(tolerance=1e-10)

>>> phase.updateOld()
>>> C.updateOld()
>>> phaseRes = 1e+10
>>> diffRes = 1e+10
>>> while phaseRes > 1e-3 or diffRes > 1e-3:
...     phaseRes = phaseEq.sweep(var=phase, dt=dt)
...     diffRes = diffusionEq.sweep(var=C, dt=dt, solver=solver)
>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Stationary phase field. Press <return> to proceed...")

phase and composition fields in equilibrium, compared with phase diagram concentrations

We verify that the bulk phases have shifted to the predicted solidus and liquidus compositions

>>> X = mesh.faceCenters[0]
>>> print(Cs.allclose(C.faceValue[X.value==0], atol=2e-4))
True
>>> print(Cl.allclose(C.faceValue[X.value==L], atol=2e-4))
True

and that the phase fraction remains unchanged

>>> print(fraction.allclose(phase.cellVolumeAverage, atol=2e-4))
1

while conserving mass overall

>>> print(Cavg.allclose(0.5, atol=1e-8))
1

We now quench by ten degrees

>>> T.setValue(T() - 10.) # K

>>> sharp.setValue(Cs, where=x < L * fraction)
>>> sharp.setValue(Cl, where=x >= L * fraction)

Because this lower temperature will induce the phase interface to move (solidify), we will need to take much smaller timesteps (the time scales of diffusion and of phase transformation compete with each other).

The CFL limit requires that no interface should advect more than one grid spacing in a timestep. We can get a rough idea for the maximum timestep we can take by looking at the velocity of convection induced by phase transformation in Eq. (6). If we assume that the phase changes from 1 to 0 in a single grid spacing, that the diffusivity is Dl at the interface, and that the term due to the difference in barrier heights is negligible:

$\vec{u}_\phi &= \frac{D_\phi}{C} \nabla \phi \\ &\approx \frac{Dl \frac{1}{2} V_m}{R T} \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] \frac{1}{\Delta x} \\ &\approx \frac{Dl \frac{1}{2} V_m}{R T} \left(L_B + L_A\right) \frac{T_M^A - T_M^B}{T_M^A + T_M^B} \frac{1}{\Delta x} \\ &\approx \unit{0.28}{\centi\meter\per\second}$

To get a $\text{CFL} = \vec{u}_\phi \Delta t / \Delta x < 1$ , we need a time step of about $\unit{10^{-5}}{\second}$ .

>>> dt = 1.e-5

>>> if __name__ == '__main__':
...     timesteps = 100
... else:
...     timesteps = 10

>>> from builtins import range
>>> for i in range(timesteps):
...     phase.updateOld()
...     C.updateOld()
...     phaseRes = 1e+10
...     diffRes = 1e+10
...     while phaseRes > 1e-3 or diffRes > 1e-3:
...         phaseRes = phaseEq.sweep(var=phase, dt=dt)
...         diffRes = diffusionEq.sweep(var=C, dt=dt, solver=solver)
...     if __name__ == '__main__':
...         viewer.plot()

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Moving phase field. Press <return> to proceed...")

phase and composition fields in during solidification, compared with final phase diagram concentrations

We see that the composition on either side of the interface approach the sharp-interface solidus and liquidus, but it will take a great many more timesteps to reach equilibrium. If we waited sufficiently long, we could again verify the final concentrations and phase fraction against the expected values.

Footnotes

examples.phase.binaryCoupled module¶

Simultaneously solve a phase-field evolution and solute diffusion problem in one-dimension.

It is straightforward to extend a phase field model to include binary alloys. As in examples.phase.simple, we will examine a 1D problem

>>> from fipy import CellVariable, Variable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, LinearLUSolver, Viewer
>>> from fipy.tools import numerix

>>> nx = 400
>>> dx = 5e-6 # cm
>>> L = nx * dx
>>> mesh = Grid1D(dx=dx, nx=nx)

The Helmholtz free energy functional can be written as the integral [1] [3] [24]

$\mathcal{F}\left(\phi, C, T\right) = \int_\mathcal{V} \left\{ f(\phi, C, T) + \frac{\kappa_\phi}{2}\abs{\nabla\phi}^2 + \frac{\kappa_C}{2}\abs{\nabla C}^2 \right\} dV$

over the volume $\mathcal{V}$ as a function of phase $\phi$ [#phi]_

>>> phase = CellVariable(name="phase", mesh=mesh, hasOld=1)

composition $C$

>>> C = CellVariable(name="composition", mesh=mesh, hasOld=1)

and temperature $T$ [#T]_

>>> T = Variable(name="temperature")

Frequently, the gradient energy term in concentration is ignored and we can derive governing equations

(7)¶ $\frac{\partial\phi}{\partial t} = M_\phi \left( \kappa_\phi \nabla^2 \phi - \frac{\partial f}{\partial \phi} \right)$

for phase and

(8)¶ $\frac{\partial C}{\partial t} = \nabla\cdot\left( M_C \nabla \frac{\partial f}{\partial C} \right)$

for solute.

The free energy density $f(\phi, C, T)$ can be constructed in many different ways. One approach is to construct free energy densities for each of the pure components, as functions of phase, e.g.

$f_A(\phi, T) = p(\phi) f_A^S(T) + \left(1 - p(\phi)\right) f_A^L(T) + \frac{W_A}{2} g(\phi)$

where $f_A^L(T)$ , $f_B^L(T)$ , $f_A^S(T)$ , and $f_B^S(T)$ are the free energy densities of the pure components. There are a variety of choices for the interpolation function $p(\phi)$ and the barrier function $g(\phi)$ ,

such as those shown in examples.phase.simple

>>> def p(phi):
...     return phi**3 * (6 * phi**2 - 15 * phi + 10)

>>> def g(phi):
...     return (phi * (1 - phi))**2

The desired thermodynamic model can then be applied to obtain $f(\phi, C, T)$ , such as for a regular solution,

$f(\phi, C, T) &= (1 - C) f_A(\phi, T) + C f_B(\phi, T) \\ &\qquad + R T \left[ (1 - C) \ln (1 - C) + C \ln C \right] + C (1 - C) \left[ \Omega_S p(\phi) + \Omega_L \left( 1 - p(\phi) \right) \right]$

where

>>> R = 8.314 # J / (mol K)

is the gas constant and $\Omega_S$ and $\Omega_L$ are the regular solution interaction parameters for solid and liquid.

Another approach is useful when the free energy densities $f^L(C, T)$ and $f^S(C,T)$ of the alloy in the solid and liquid phases are known. This might be the case when the two different phases have different thermodynamic models or when one or both is obtained from a Calphad code. In this case, we can construct

$f(\phi, C, T) = p(\phi) f^S(C,T) + \left(1 - p(\phi)\right) f^L(C, T) + \left[ (1-C) \frac{W_A}{2} + C \frac{W_B}{2} \right] g(\phi).$

When the thermodynamic models are the same in both phases, both approaches should yield the same result.

We choose the first approach and make the simplifying assumptions of an ideal solution and that

$f_A^L(T) & = 0 \\ f_A^S(T) - f_A^L(T) &= \frac{L_A\left(T - T_M^A\right)}{T_M^A}$

and likewise for component $B$ .

>>> LA = 2350. # J / cm**3
>>> LB = 1728. # J / cm**3
>>> TmA = 1728. # K
>>> TmB = 1358. # K

>>> enthalpyA = LA * (T - TmA) / TmA
>>> enthalpyB = LB * (T - TmB) / TmB

This relates the difference between the free energy densities of the pure solid and pure liquid phases to the latent heat $L_A$ and the pure component melting point $T_M^A$ , such that

$f_A(\phi, T) = \frac{L_A\left(T - T_M^A\right)}{T_M^A} p(\phi) + \frac{W_A}{2} g(\phi).$

With these assumptions

(9)¶ $\frac{\partial f}{\partial \phi} &= (1-C) \frac{\partial f_A}{\partial \phi} + C \frac{\partial f_B}{\partial \phi} \nonumber \\ &= \left\{ (1-C) \frac{L_A\left(T - T_M^A\right)}{T_M^A} + C \frac{L_B\left(T - T_M^B\right)}{T_M^B} \right\} p'(\phi) + \left\{ (1-C) \frac{W_A}{2} + C \frac{W_B}{2} \right\} g'(\phi)$

and

(10)¶ $\frac{\partial f}{\partial C} &= \left[f_B(\phi, T) + \frac{R T}{V_m} \ln C\right] - \left[f_A(\phi, T) + \frac{R T}{V_m} \ln (1-C) \right] \nonumber \\ &= \left[\mu_B(\phi, C, T) - \mu_A(\phi, C, T) \right] / V_m$

where $\mu_A$ and $\mu_B$ are the classical chemical potentials for the binary species. $p'(\phi)$ and $g'(\phi)$ are the partial derivatives of of $p$ and $g$ with respect to $\phi$

>>> def pPrime(phi):
...     return 30. * g(phi)

>>> def gPrime(phi):
...     return 2. * phi * (1 - phi) * (1 - 2 * phi)

$V_m$ is the molar volume, which we take to be independent of concentration and phase

>>> Vm = 7.42 # cm**3 / mol

On comparison with examples.phase.simple, we can see that the present form of the phase field equation is identical to the one found earlier, with the source now composed of the concentration-weighted average of the source for either pure component. We let the pure component barriers equal the previous value

>>> deltaA = deltaB = 1.5 * dx
>>> sigmaA = 3.7e-5 # J / cm**2
>>> sigmaB = 2.9e-5 # J / cm**2
>>> betaA = 0.33 # cm / (K s)
>>> betaB = 0.39 # cm / (K s)
>>> kappaA = 6 * sigmaA * deltaA # J / cm
>>> kappaB = 6 * sigmaB * deltaB # J / cm
>>> WA = 6 * sigmaA / deltaA # J / cm**3
>>> WB = 6 * sigmaB / deltaB # J / cm**3

and define the averages

>>> W = (1 - C) * WA / 2. + C * WB / 2.
>>> enthalpy = (1 - C) * enthalpyA + C * enthalpyB

We can now linearize the source exactly as before

>>> mPhi = -((1 - 2 * phase) * W + 30 * phase * (1 - phase) * enthalpy)
>>> dmPhidPhi = 2 * W - 30 * (1 - 2 * phase) * enthalpy
>>> S1 = dmPhidPhi * phase * (1 - phase) + mPhi * (1 - 2 * phase)
>>> S0 = mPhi * phase * (1 - phase) - S1 * phase

Using the same gradient energy coefficient and phase field mobility

>>> kappa = (1 - C) * kappaA + C * kappaB
>>> Mphi = TmA * betaA / (6 * LA * deltaA)

we define the phase field equation

>>> phaseEq = (TransientTerm(1/Mphi, var=phase) == DiffusionTerm(coeff=kappa, var=phase)
...            + S0 + ImplicitSourceTerm(coeff=S1, var=phase))

When coding explicitly, it is typical to simply write a function to evaluate the chemical potentials $\mu_A$ and $\mu_B$ and then perform the finite differences necessary to calculate their gradient and divergence, e.g.,:

def deltaChemPot(phase, C, T):
    return ((Vm * (enthalpyB * p(phase) + WA * g(phase)) + R * T * log(1 - C)) -
            (Vm * (enthalpyA * p(phase) + WA * g(phase)) + R * T * log(C)))

for j in range(faces):
    flux[j] = ((Mc[j+.5] + Mc[j-.5]) / 2) \
      * (deltaChemPot(phase[j+.5], C[j+.5], T) \
        - deltaChemPot(phase[j-.5], C[j-.5], T)) / dx

for j in range(cells):
    diffusion = (flux[j+.5] - flux[j-.5]) / dx

where we neglect the details of the outer boundaries (j = 0 and j = N) or exactly how to translate j+.5 or j-.5 into an array index, much less the complexities of higher dimensions. FiPy can handle all of these issues automatically, so we could just write:

chemPotA = Vm * (enthalpyA * p(phase) + WA * g(phase)) + R * T * log(C)
chemPotB = Vm * (enthalpyB * p(phase) + WB * g(phase)) + R * T * log(1-C)
flux = Mc * (chemPotB - chemPotA).faceGrad
eq = TransientTerm() == flux.divergence

Although the second syntax would essentially work as written, such an explicit implementation would be very slow. In order to take advantage of FiPy’s implicit solvers, it is necessary to reduce Eq. (2) to the canonical form of Eq. (2), hence we must expand Eq. (4) as

$\frac{\partial f}{\partial C} = \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] p(\phi) + \frac{R T}{V_m} \left[\ln C - \ln (1-C)\right] + \frac{W_B - W_A}{2} g(\phi)$

In either bulk phase, $\nabla p(\phi) = \nabla g(\phi) = 0$ , so we can then reduce Eq. (2) to

(11)¶ $\frac{\partial C}{\partial t} &= \nabla\cdot\left( M_C \nabla \left\{ \frac{R T}{V_m} \left[\ln C - \ln (1-C)\right] \right\} \right) \nonumber \\ &= \nabla\cdot\left[ \frac{M_C R T}{C (1-C) V_m} \nabla C \right]$

and, by comparison with Fick’s second law

$\frac{\partial C}{\partial t} = \nabla\cdot\left[D \nabla C\right],$

we can associate the mobility $M_C$ with the intrinsic diffusivity $D_C$ by $M_C \equiv D_C C (1-C) V_m / R T$ and write Eq. (2) as

(12)¶ $\frac{\partial C}{\partial t} &= \nabla\cdot\left( D_C \nabla C \right) \nonumber \\ &\qquad + \nabla\cdot\left( \frac{D_C C (1 - C) V_m}{R T} \left\{ \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] \nabla p(\phi) + \frac{W_B - W_A}{2} \nabla g(\phi) \right\} \right). \\ &= \nabla\cdot\left( D_C \nabla C \right) \nonumber \\ &\qquad + \nabla\cdot\left( \frac{D_C C (1 - C) V_m}{R T} \left\{ \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] p'(\phi) + \frac{W_B - W_A}{2} g'(\phi) \right\} \nabla \phi \right).$

The first term is clearly a DiffusionTerm in $C$ . The second is a DiffusionTerm in $\phi$ with a diffusion coefficient

$D_{\phi}(C, \phi) = \frac{D_C C (1 - C) V_m}{R T} \left\{ \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] p'(\phi) + \frac{W_B - W_A}{2} g'(\phi) \right\},$

such that

$\frac{\partial C}{\partial t} = \nabla\cdot\left( D_C \nabla C \right) + \nabla\cdot\left(D_\phi \nabla \phi \right)$

or

>>> Dl = Variable(value=1e-5) # cm**2 / s
>>> Ds = Variable(value=1e-9) # cm**2 / s
>>> Dc = (Ds - Dl) * phase.arithmeticFaceValue + Dl

>>> Dphi = ((Dc * C.harmonicFaceValue * (1 - C.harmonicFaceValue) * Vm / (R * T))
...         * ((enthalpyB - enthalpyA) * pPrime(phase.arithmeticFaceValue)
...            + 0.5 * (WB - WA) * gPrime(phase.arithmeticFaceValue)))

>>> diffusionEq = (TransientTerm(var=C)
...                == DiffusionTerm(coeff=Dc, var=C)
...                + DiffusionTerm(coeff=Dphi, var=phase))

>>> eq = phaseEq & diffusionEq

We initialize the phase field to a step function in the middle of the domain

>>> phase.setValue(1.)
>>> phase.setValue(0., where=mesh.cellCenters[0] > L/2.)

and start with a uniform composition field $C = 1/2$

>>> C.setValue(0.5)

In equilibrium, $\mu_A(0, C_L, T) = \mu_A(1, C_S, T)$ and $\mu_B(0, C_L, T) = \mu_B(1, C_S, T)$ and, for ideal solutions, we can deduce the liquidus and solidus compositions as

$C_L &= \frac{1 - \exp\left(-\frac{L_A\left(T - T_M^A\right)}{T_M^A}\frac{V_m}{R T}\right)} {\exp\left(-\frac{L_B\left(T - T_M^B\right)}{T_M^B}\frac{V_m}{R T}\right) - \exp\left(-\frac{L_A\left(T - T_M^A\right)}{T_M^A}\frac{V_m}{R T}\right)} \\ C_S &= \exp\left(-\frac{L_B\left(T - T_M^B\right)}{T_M^B}\frac{V_m}{R T}\right) C_L$

>>> Cl = ((1. - numerix.exp(-enthalpyA * Vm / (R * T)))
...       / (numerix.exp(-enthalpyB * Vm / (R * T)) - numerix.exp(-enthalpyA * Vm / (R * T))))
>>> Cs = numerix.exp(-enthalpyB * Vm / (R * T)) * Cl

The phase fraction is predicted by the lever rule

>>> Cavg = C.cellVolumeAverage
>>> fraction = (Cl - Cavg) / (Cl - Cs)

For the special case of fraction = Cavg = 0.5, a little bit of algebra reveals that the temperature that leaves the phase fraction unchanged is given by

>>> T.setValue((LA + LB) * TmA * TmB / (LA * TmB + LB * TmA))

In this simple, binary, ideal solution case, we can derive explicit expressions for the solidus and liquidus compositions. In general, this may not be possible or practical. In that event, the root-finding facilities in SciPy can be used.

We’ll need a function to return the two conditions for equilibrium

$0 = \mu_A(1, C_S, T) - \mu_A(0, C_L, T) &= \frac{L_A\left(T - T_M^A\right)}{T_M^A} V_m + R T \ln (1 - C_S) - R T \ln (1 - C_L) \\ 0 = \mu_B(1, C_S, T) - \mu_B(0, C_L, T) &= \frac{L_B\left(T - T_M^B\right)}{T_M^B} V_m + R T \ln C_S - R T \ln C_L$

>>> def equilibrium(C):
...     return [numerix.array(enthalpyA * Vm
...                           + R * T * numerix.log(1 - C[0])
...                           - R * T * numerix.log(1 - C[1])),
...             numerix.array(enthalpyB * Vm
...                           + R * T * numerix.log(C[0])
...                           - R * T * numerix.log(C[1]))]

and we’ll have much better luck if we also supply the Jacobian

$\left[\begin{matrix} \frac{\partial(\mu_A^S - \mu_A^L)}{\partial C_S} & \frac{\partial(\mu_A^S - \mu_A^L)}{\partial C_L} \\ \frac{\partial(\mu_B^S - \mu_B^L)}{\partial C_S} & \frac{\partial(\mu_B^S - \mu_B^L)}{\partial C_L} \end{matrix}\right] = R T\left[\begin{matrix} -\frac{1}{1-C_S} & \frac{1}{1-C_L} \\ \frac{1}{C_S} & -\frac{1}{C_L} \end{matrix}\right]$

>>> def equilibriumJacobian(C):
...     return R * T * numerix.array([[-1. / (1 - C[0]), 1. / (1 - C[1])],
...                                   [ 1. / C[0],      -1. / C[1]]])

>>> try:
...     from scipy.optimize import fsolve 
...     CsRoot, ClRoot = fsolve(func=equilibrium, x0=[0.5, 0.5],
...                             fprime=equilibriumJacobian) 
... except ImportError:
...     ClRoot = CsRoot = 0
...     print("The SciPy library is not available to calculate the solidus and ... liquidus concentrations")

>>> print(Cl.allclose(ClRoot)) 
1
>>> print(Cs.allclose(CsRoot)) 
1

We plot the result against the sharp interface solution

>>> sharp = CellVariable(name="sharp", mesh=mesh)
>>> x = mesh.cellCenters[0]
>>> sharp.setValue(Cs, where=x < L * fraction)
>>> sharp.setValue(Cl, where=x >= L * fraction)

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=(phase, C, sharp),
...                     datamin=0., datamax=1.)
...     viewer.plot()

Because the phase field interface will not move, and because we’ve seen in earlier examples that the diffusion problem is unconditionally stable, we need take only one very large timestep to reach equilibrium

>>> dt = 1.e5

Because the phase field equation is coupled to the composition through enthalpy and W and the diffusion equation is coupled to the phase field through phaseTransformationVelocity, it is necessary sweep this non-linear problem to convergence. We use the “residual” of the equations (a measure of how well they think they have solved the given set of linear equations) as a test for how long to sweep.

We now use the “sweep()” method instead of “solve()” because we require the residual.

>>> solver = LinearLUSolver(tolerance=1e-10)

>>> phase.updateOld()
>>> C.updateOld()
>>> res = 1.
>>> initialRes = None
>>> sweep = 0

>>> while res > 1e-4 and sweep < 20:
...     res = eq.sweep(dt=dt, solver=solver)
...     if initialRes is None:
...         initialRes = res
...     res = res / initialRes
...     sweep += 1

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Stationary phase field. Press <return> to proceed...")

documentation/generated/examples/binary/stationary.*

We verify that the bulk phases have shifted to the predicted solidus and liquidus compositions

>>> X = mesh.faceCenters[0]
>>> print(Cs.allclose(C.faceValue[X.value==0], atol=1e-2))
True
>>> print(Cl.allclose(C.faceValue[X.value==L], atol=1e-2))
True

and that the phase fraction remains unchanged

>>> print(fraction.allclose(phase.cellVolumeAverage, atol=2e-4))
1

while conserving mass overall

>>> print(Cavg.allclose(0.5, atol=1e-8))
1

We now quench by ten degrees

>>> T.setValue(T() - 10.) # K

>>> sharp.setValue(Cs, where=x < L * fraction)
>>> sharp.setValue(Cl, where=x >= L * fraction)

Because this lower temperature will induce the phase interface to move (solidify), we will need to take much smaller timesteps (the time scales of diffusion and of phase transformation compete with each other).

The CFL limit requires that no interface should advect more than one grid spacing in a timestep. We can get a rough idea for the maximum timestep we can take by looking at the velocity of convection induced by phase transformation in Eq. (6) (even though there is no explicit convection in the coupled form used for this example, the principle remains the same). If we assume that the phase changes from 1 to 0 in a single grid spacing, that the diffusivity is Dl at the interface, and that the term due to the difference in barrier heights is negligible:

$\vec{u}_\phi &= \frac{D_\phi}{C} \nabla \phi \\ &\approx \frac{Dl \frac{1}{2} V_m}{R T} \left[ \frac{L_B\left(T - T_M^B\right)}{T_M^B} - \frac{L_A\left(T - T_M^A\right)}{T_M^A} \right] \frac{1}{\Delta x} \\ &\approx \frac{Dl \frac{1}{2} V_m}{R T} \left(L_B + L_A\right) \frac{T_M^A - T_M^B}{T_M^A + T_M^B} \frac{1}{\Delta x} \\ &\approx \unit{0.28}{\centi\meter\per\second}$

To get a $\text{CFL} = \vec{u}_\phi \Delta t / \Delta x < 1$ , we need a time step of about $\unit{10^{-5}}{\second}$ .

>>> dt = 1.e-5

>>> if __name__ == '__main__':
...     timesteps = 100
... else:
...     timesteps = 10

>>> from builtins import range
>>> for i in range(timesteps):
...     phase.updateOld()
...     C.updateOld()
...     res = 1e+10
...     sweep = 0
...     while res > 1e-3 and sweep < 20:
...         res = eq.sweep(dt=dt, solver=solver)
...         sweep += 1
...     if __name__ == '__main__':
...         viewer.plot()

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Moving phase field. Press <return> to proceed...")

documentation/generated/examples/binary/moving.*

We see that the composition on either side of the interface approach the sharp-interface solidus and liquidus, but it will take a great many more timesteps to reach equilibrium. If we waited sufficiently long, we could again verify the final concentrations and phase fraction against the expected values.

Footnotes

examples.phase.polyxtal module¶

Solve the dendritic growth of nuclei and subsequent grain impingement.

To convert a liquid material to a solid, it must be cooled to a temperature below its melting point (known as “undercooling” or “supercooling”). The rate of solidification is often assumed (and experimentally found) to be proportional to the undercooling. Under the right circumstances, the solidification front can become unstable, leading to dendritic patterns. Warren, Kobayashi, Lobkovsky and Carter [13] have described a phase field model (“Allen-Cahn”, “non-conserved Ginsberg-Landau”, or “model A” of Hohenberg & Halperin) of such a system, including the effects of discrete crystalline orientations (anisotropy).

We start with a regular 2D Cartesian mesh

>>> from fipy import CellVariable, Variable, ModularVariable, Grid2D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, MatplotlibViewer, Matplotlib2DGridViewer, MultiViewer
>>> from fipy.tools import numerix
>>> dx = dy = 0.025
>>> if __name__ == "__main__":
...     nx = ny = 200
... else:
...     nx = ny = 200
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

and we’ll take fixed timesteps

>>> dt = 5e-4

We consider the simultaneous evolution of a “phase field” variable $\phi$ (taken to be 0 in the liquid phase and 1 in the solid)

>>> phase = CellVariable(name=r'$\phi$', mesh=mesh, hasOld=True)

a dimensionless undercooling $\Delta T$ ( $\Delta T = 0$ at the melting point)

>>> dT = CellVariable(name=r'$\Delta T$', mesh=mesh, hasOld=True)

and an orientation $-\pi < \theta \le \pi$

>>> theta = ModularVariable(name=r'$\theta$', mesh=mesh, hasOld=True)
>>> theta.value = -numerix.pi + 0.0001

The hasOld flag causes the storage of the value of variable from the previous timestep. This is necessary for solving equations with non-linear coefficients or for coupling between PDEs.

The governing equation for the temperature field is the heat flux equation, with a source due to the latent heat of solidification

$\frac{\partial \Delta T}{\partial t} = D_T \nabla^2 \Delta T + \frac{\partial \phi}{\partial t} + c\left(T_0 - T\right)$

>>> DT = 2.25
>>> q = Variable(0.)
>>> T_0 = -0.1
>>> heatEq = (TransientTerm()
...           == DiffusionTerm(DT)
...           + (phase - phase.old) / dt
...           + q * T_0 - ImplicitSourceTerm(q))

The governing equation for the phase field is

$\tau_{\phi} \frac{\partial \phi}{\partial t} = \nabla \cdot \mathsf{D} \nabla \phi + \phi ( 1 - \phi ) m ( \phi , \Delta T) - 2 s \phi | \nabla \theta | - \epsilon^2 \phi | \nabla \theta |^2$

where

$m(\phi, \Delta T) = \phi - \frac{1}{2} - \frac{ \kappa_1 }{ \pi } \arctan \left( \kappa_2 \Delta T \right)$

represents a source of anisotropy. The coefficient $\mathsf{D}$ is an anisotropic diffusion tensor in two dimensions

$\mathsf{D} = \alpha^2 \left( 1 + c \beta \right) \left[ \begin{matrix} 1 + c \beta & -c \frac{\partial \beta}{\partial \psi} \\ c \frac{\partial \beta}{\partial \psi} & 1 + c \beta \end{matrix} \right]$

where $\beta = \frac{ 1 - \Phi^2 } { 1 + \Phi^2}$ , $\Phi = \tan \left( \frac{ N } { 2 } \psi \right)$ , $\psi = \theta + \arctan \frac{\partial \phi / \partial y}{\partial \phi / \partial x}$ , $\theta$ is the orientation, and $N$ is the symmetry.

>>> alpha = 0.015
>>> c = 0.02
>>> N = 4.

>>> psi = theta.arithmeticFaceValue + numerix.arctan2(phase.faceGrad[1],
...                                                   phase.faceGrad[0])
>>> Phi = numerix.tan(N * psi / 2)
>>> PhiSq = Phi**2
>>> beta = (1. - PhiSq) / (1. + PhiSq)
>>> DbetaDpsi = -N * 2 * Phi / (1 + PhiSq)
>>> Ddia = (1.+ c * beta)

>>> Doff = c * DbetaDpsi
>>> I0 = Variable(value=((1, 0), (0, 1)))
>>> I1 = Variable(value=((0, -1), (1, 0)))
>>> D = alpha**2 * Ddia * (Ddia * I0 + Doff * I1)

With these expressions defined, we can construct the phase field equation as

>>> tau_phase = 3e-4
>>> kappa1 = 0.9
>>> kappa2 = 20.
>>> epsilon = 0.008
>>> s = 0.01
>>> thetaMag = theta.grad.mag
>>> phaseEq = (TransientTerm(tau_phase)
...            == DiffusionTerm(D)
...            + ImplicitSourceTerm((phase - 0.5 - kappa1 / numerix.pi * numerix.arctan(kappa2 * dT))
...                                 * (1 - phase)
...                                 - (2 * s + epsilon**2 * thetaMag) * thetaMag))

The governing equation for orientation is given by

$P(\epsilon | \nabla \theta |) \tau_{\theta} \phi^2 \frac{\partial \theta}{\partial t} = \nabla \cdot \left[ \phi^2 \left( \frac{s}{| \nabla \theta |} + \epsilon^2 \right) \nabla \theta \right]$

where

$P(w) = 1 - \exp{(-\beta w)} + \frac{\mu}{\epsilon} \exp{(-\beta w)}$

The theta equation is built in the following way. The details for this equation are fairly involved, see J. A. Warren et al.. The main detail is that a source must be added to correct for the discretization of theta on the circle.

>>> tau_theta = 3e-3
>>> mu = 1e3
>>> gamma = 1e3
>>> thetaSmallValue = 1e-6
>>> phaseMod = phase + ( phase < thetaSmallValue ) * thetaSmallValue
>>> beta_theta = 1e5
>>> expo = epsilon * beta_theta * theta.grad.mag
>>> expo = (expo < 100.) * (expo - 100.) + 100.
>>> Pfunc = 1. + numerix.exp(-expo) * (mu / epsilon - 1.)

>>> gradMagTheta = theta.faceGrad.mag
>>> eps = 1. / gamma / 10.
>>> gradMagTheta += (gradMagTheta < eps) * eps
>>> IGamma = (gradMagTheta > 1. / gamma) * (1 / gradMagTheta - gamma) + gamma
>>> v_theta = phase.arithmeticFaceValue * (s * IGamma + epsilon**2)
>>> D_theta = phase.arithmeticFaceValue**2 * (s * IGamma + epsilon**2)

The source term requires the evaluation of the face gradient without the modular operator. thetafaceGradNoMod evaluates the gradient without modular arithmetic.

>>> thetaEq = (TransientTerm(tau_theta * phaseMod**2 * Pfunc)
...            == DiffusionTerm(D_theta)
...            + (D_theta * (theta.faceGrad - theta.faceGradNoMod)).divergence)

We seed a circular solidified region in the center

>>> x, y = mesh.cellCenters
>>> numSeeds = 10
>>> numerix.random.seed(12345)
>>> for Cx, Cy, orientation in numerix.random.random([numSeeds, 3]):
...     radius = dx * 5.
...     seed = ((x - Cx * nx * dx)**2 + (y - Cy * ny * dy)**2) < radius**2
...     phase[seed] = 1.
...     theta[seed] = numerix.pi * (2 * orientation - 1)

and quench the entire simulation domain below the melting point

>>> dT.setValue(-0.5)

In a real solidification process, dendritic branching is induced by small thermal fluctuations along an otherwise smooth surface, but the granularity of the Mesh is enough “noise” in this case, so we don’t need to explicitly introduce randomness, the way we did in the Cahn-Hilliard problem.

FiPy’s viewers are utilitarian, striving to let the user see something, regardless of their operating system or installed packages, so you the default color scheme of grain orientation won’t be very informative “out of the box”. Because all of Python is accessible and FiPy is object oriented, it is not hard to adapt one of the existing viewers to create a specialized display:

>>> from builtins import zip
>>> if __name__ == "__main__":
...     try:
...         class OrientationViewer(Matplotlib2DGridViewer):
...             def __init__(self, phase, orientation, title=None, limits={}, **kwlimits):
...                 self.phase = phase
...                 Matplotlib2DGridViewer.__init__(self, vars=(orientation,), title=title,
...                                                 limits=limits, colorbar=None, **kwlimits)
...
...                 # make room for non-existent colorbar
...                 # stolen from matplotlib.colorbar.make_axes
...                 # https://github.com/matplotlib/matplotlib/blob
...                 #   /ec1cd2567521c105a451ce15e06de10715f8b54d/lib
...                 #   /matplotlib/colorbar.py#L838
...                 fraction = 0.15
...                 pb = self.axes.get_position(original=True).frozen()
...                 pad = 0.05
...                 x1 = 1.0-fraction
...                 pb1, pbx, pbcb = pb.splitx(x1-pad, x1)
...                 panchor = (1.0, 0.5)
...                 self.axes.set_position(pb1)
...                 self.axes.set_anchor(panchor)
...
...                 # make the gnomon
...                 fig = self.axes.get_figure()
...                 self.gnomon = fig.add_axes([0.85, 0.425, 0.15, 0.15], polar=True)
...                 self.gnomon.set_thetagrids([180, 270, 0, 90],
...                                            [r"$\pm\pi$", r"$-\frac{\pi}{2}$", "$0$", r"$+\frac{\pi}{2}$"],
...                                            frac=1.3)
...                 self.gnomon.set_theta_zero_location("N")
...                 self.gnomon.set_theta_direction(-1)
...                 self.gnomon.set_rgrids([1.], [""])
...                 N = 100
...                 theta = numerix.arange(-numerix.pi, numerix.pi, 2 * numerix.pi / N)
...                 radii = numerix.ones((N,))
...                 bars = self.gnomon.bar(theta, radii, width=2 * numerix.pi / N, bottom=0.0)
...                 colors = self._orientation_and_phase_to_rgb(orientation=numerix.array([theta]), phase=1.)
...                 for c, t, bar in zip(colors[0], theta, bars):
...                     bar.set_facecolor(c)
...                     bar.set_edgecolor(c)
...
...             def _reshape(self, var):
...                 '''return values of var in an 2D array'''
...                 return numerix.reshape(numerix.array(var),
...                                        var.mesh.shape[::-1])[::-1]
...
...             @staticmethod
...             def _orientation_and_phase_to_rgb(orientation, phase):
...                 from matplotlib import colors
...
...                 hsv = numerix.empty(orientation.shape + (3,))
...                 hsv[..., 0] = (orientation / numerix.pi + 1) / 2.
...                 hsv[..., 1] = 1.
...                 hsv[..., 2] = phase
...
...                 return colors.hsv_to_rgb(hsv)
...
...             @property
...             def _data(self):
...                 '''convert phase and orientation to rgb image array
...
...                 orientation (-pi, pi) -> hue (0, 1)
...                 phase (0, 1) -> value (0, 1)
...                 '''
...                 orientation = self._reshape(self.vars[0])
...                 phase = self._reshape(self.phase)
...
...                 return self._orientation_and_phase_to_rgb(orientation, phase)
...
...             def _plot(self):
...                 self.image.set_data(self._data)
...
...         from matplotlib import pyplot
...         pyplot.ion()
...         w, h = pyplot.figaspect(1.)
...         fig = pyplot.figure(figsize=(2*w, h))
...         timer = fig.text(0.1, 0.9, "t = %.3f" % 0, fontsize=18)
...
...         viewer = MultiViewer(viewers=(MatplotlibViewer(vars=dT,
...                                                        cmap=pyplot.cm.hot,
...                                                        datamin=-0.5,
...                                                        datamax=0.5,
...                                                        axes=fig.add_subplot(121)),
...                                       OrientationViewer(phase=phase,
...                                                         orientation=theta,
...                                                         title=theta.name,
...                                                         axes=fig.add_subplot(122))))
...     except ImportError:
...         viewer = MultiViewer(viewers=(Viewer(vars=dT,
...                                              datamin=-0.5,
...                                              datamax=0.5),
...                                       Viewer(vars=phase,
...                                              datamin=0.,
...                                              datamax=1.),
...                                       Viewer(vars=theta,
...                                              datamin=-numerix.pi,
...                                              datamax=numerix.pi)))
...     viewer.plot()

and iterate the solution in time, plotting as we go,

>>> if __name__ == "__main__":
...     total_time = 2.
... else:
...     total_time = dt * 10
>>> elapsed = 0.
>>> save_interval = 0.002
>>> save_at = save_interval

>>> while elapsed < total_time:
...     if elapsed > 0.3:
...         q.value = 100
...     phase.updateOld()
...     dT.updateOld()
...     theta.updateOld()
...     thetaEq.solve(theta, dt=dt)
...     phaseEq.solve(phase, dt=dt)
...     heatEq.solve(dT, dt=dt)
...     elapsed += dt
...     if __name__ == "__main__" and elapsed >= save_at:
...         timer.set_text("t = %.3f" % elapsed)
...         viewer.plot()
...         save_at += save_interval

undercooling and grain orientation during solidification of a collection of anisotropic seeds

The non-uniform temperature results from the release of latent heat at the solidifying interface. The dendrite arms grow fastest where the temperature gradient is steepest.

examples.phase.polyxtalCoupled module¶

Simultaneously solve the dendritic growth of nuclei and subsequent grain impingement.

To convert a liquid material to a solid, it must be cooled to a temperature below its melting point (known as “undercooling” or “supercooling”). The rate of solidification is often assumed (and experimentally found) to be proportional to the undercooling. Under the right circumstances, the solidification front can become unstable, leading to dendritic patterns. Warren, Kobayashi, Lobkovsky and Carter [13] have described a phase field model (“Allen-Cahn”, “non-conserved Ginsberg-Landau”, or “model A” of Hohenberg & Halperin) of such a system, including the effects of discrete crystalline orientations (anisotropy).

We start with a regular 2D Cartesian mesh

>>> from fipy import  CellVariable, Variable, ModularVariable, Grid2D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, PowerLawConvectionTerm, MatplotlibViewer, Matplotlib2DGridViewer, MultiViewer
>>> from fipy.tools import numerix
>>> dx = dy = 0.025
>>> if __name__ == "__main__":
...     nx = ny = 200
... else:
...     nx = ny = 20
>>> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

and we’ll take fixed timesteps

>>> dt = 5e-4

We consider the simultaneous evolution of a “phase field” variable $\phi$ (taken to be 0 in the liquid phase and 1 in the solid)

>>> phase = CellVariable(name=r'$\phi$', mesh=mesh, hasOld=True)

a dimensionless undercooling $\Delta T$ ( $\Delta T = 0$ at the melting point)

>>> dT = CellVariable(name=r'$\Delta T$', mesh=mesh, hasOld=True)

and an orientation $-\pi < \theta \le \pi$

>>> theta = ModularVariable(name=r'$\theta$', mesh=mesh, hasOld=True)
>>> theta.value = -numerix.pi + 0.0001

The hasOld flag causes the storage of the value of variable from the previous timestep. This is necessary for solving equations with non-linear coefficients or for coupling between PDEs.

The governing equation for the temperature field is the heat flux equation, with a source due to the latent heat of solidification

$\frac{\partial \Delta T}{\partial t} = D_T \nabla^2 \Delta T + \frac{\partial \phi}{\partial t} + c\left(T_0 - T\right)$

>>> DT = 2.25
>>> q = Variable(0.)
>>> T_0 = -0.1
>>> heatEq = (TransientTerm(var=dT)
...           == DiffusionTerm(coeff=DT, var=dT)
...           + TransientTerm(var=phase)
...           + q * T_0 - ImplicitSourceTerm(coeff=q, var=dT))

The governing equation for the phase field is

$\tau_{\phi} \frac{\partial \phi}{\partial t} = \nabla \cdot \mathsf{D} \nabla \phi + \phi ( 1 - \phi ) m ( \phi , \Delta T) - 2 s \phi | \nabla \theta | - \epsilon^2 \phi | \nabla \theta |^2$

where

$m(\phi, \Delta T) = \phi - \frac{1}{2} - \frac{ \kappa_1 }{ \pi } \arctan \left( \kappa_2 \Delta T \right)$

represents a source of anisotropy. The coefficient $\mathsf{D}$ is an anisotropic diffusion tensor in two dimensions

$\mathsf{D} = \alpha^2 \left( 1 + c \beta \right) \left[ \begin{matrix} 1 + c \beta & -c \frac{\partial \beta}{\partial \psi} \\ c \frac{\partial \beta}{\partial \psi} & 1 + c \beta \end{matrix} \right]$

where $\beta = \frac{ 1 - \Phi^2 } { 1 + \Phi^2}$ , $\Phi = \tan \left( \frac{ N } { 2 } \psi \right)$ , $\psi = \theta + \arctan \frac{\partial \phi / \partial y}{\partial \phi / \partial x}$ , $\theta$ is the orientation, and $N$ is the symmetry.

>>> alpha = 0.015
>>> c = 0.02
>>> N = 4.

>>> psi = theta.arithmeticFaceValue + numerix.arctan2(phase.faceGrad[1],
...                                                   phase.faceGrad[0])
>>> Phi = numerix.tan(N * psi / 2)
>>> PhiSq = Phi**2
>>> beta = (1. - PhiSq) / (1. + PhiSq)
>>> DbetaDpsi = -N * 2 * Phi / (1 + PhiSq)
>>> Ddia = (1.+ c * beta)

>>> Doff = c * DbetaDpsi
>>> I0 = Variable(value=((1, 0), (0, 1)))
>>> I1 = Variable(value=((0, -1), (1, 0)))
>>> D = alpha**2 * Ddia * (Ddia * I0 + Doff * I1)

With these expressions defined, we can construct the phase field equation as

>>> tau_phase = 3e-4
>>> kappa1 = 0.9
>>> kappa2 = 20.
>>> epsilon = 0.008
>>> s = 0.01
>>> thetaMag = theta.grad.mag
>>> phaseEq = (TransientTerm(coeff=tau_phase, var=phase)
...            == DiffusionTerm(coeff=D, var=phase)
...            + ImplicitSourceTerm(coeff=((phase - 0.5 - kappa1 / numerix.pi * numerix.arctan(kappa2 * dT))
...                                        * (1 - phase)
...                                        - (2 * s + epsilon**2 * thetaMag) * thetaMag),
...                                 var=phase))

The governing equation for orientation is given by

$P(\epsilon | \nabla \theta |) \tau_{\theta} \phi^2 \frac{\partial \theta}{\partial t} = \nabla \cdot \left[ \phi^2 \left( \frac{s}{| \nabla \theta |} + \epsilon^2 \right) \nabla \theta \right]$

where

$P(w) = 1 - \exp{(-\beta w)} + \frac{\mu}{\epsilon} \exp{(-\beta w)}$

The theta equation is built in the following way. The details for this equation are fairly involved, see J. A. Warren et al.. The main detail is that a source must be added to correct for the discretization of theta on the circle.

>>> tau_theta = 3e-3
>>> mu = 1e3
>>> gamma = 1e3
>>> thetaSmallValue = 1e-6
>>> phaseMod = phase + ( phase < thetaSmallValue ) * thetaSmallValue
>>> beta_theta = 1e5
>>> expo = epsilon * beta_theta * theta.grad.mag
>>> expo = (expo < 100.) * (expo - 100.) + 100.
>>> Pfunc = 1. + numerix.exp(-expo) * (mu / epsilon - 1.)

>>> gradMagTheta = theta.faceGrad.mag
>>> eps = 1. / gamma / 10.
>>> gradMagTheta += (gradMagTheta < eps) * eps
>>> IGamma = (gradMagTheta > 1. / gamma) * (1 / gradMagTheta - gamma) + gamma
>>> v_theta = phase.arithmeticFaceValue * (s * IGamma + epsilon**2)
>>> D_theta = phase.arithmeticFaceValue**2 * (s * IGamma + epsilon**2)

The source term requires the evaluation of the face gradient without the modular operator. thetafaceGradNoMod evaluates the gradient without modular arithmetic.

>>> thetaEq = (TransientTerm(coeff=tau_theta * phaseMod**2 * Pfunc, var=theta)
...            == DiffusionTerm(coeff=D_theta, var=theta)
...            + PowerLawConvectionTerm(coeff=v_theta * (theta.faceGrad - theta.faceGradNoMod), var=phase))

We seed a circular solidified region in the center

>>> x, y = mesh.cellCenters
>>> numSeeds = 10
>>> numerix.random.seed(12345)
>>> for Cx, Cy, orientation in numerix.random.random([numSeeds, 3]):
...     radius = dx * 5.
...     seed = ((x - Cx * nx * dx)**2 + (y - Cy * ny * dy)**2) < radius**2
...     phase[seed] = 1.
...     theta[seed] = numerix.pi * (2 * orientation - 1)

and quench the entire simulation domain below the melting point

>>> dT.setValue(-0.5)

In a real solidification process, dendritic branching is induced by small thermal fluctuations along an otherwise smooth surface, but the granularity of the Mesh is enough “noise” in this case, so we don’t need to explicitly introduce randomness, the way we did in the Cahn-Hilliard problem.

FiPy’s viewers are utilitarian, striving to let the user see something, regardless of their operating system or installed packages, so you the default color scheme of grain orientation won’t be very informative “out of the box”. Because all of Python is accessible and FiPy is object oriented, it is not hard to adapt one of the existing viewers to create a specialized display:

>>> from builtins import zip
>>> if __name__ == "__main__":
...     try:
...         class OrientationViewer(Matplotlib2DGridViewer):
...             def __init__(self, phase, orientation, title=None, limits={}, **kwlimits):
...                 self.phase = phase
...                 Matplotlib2DGridViewer.__init__(self, vars=(orientation,), title=title,
...                                                 limits=limits, colorbar=None, **kwlimits)
...
...                 # make room for non-existent colorbar
...                 # stolen from matplotlib.colorbar.make_axes
...                 # https://github.com/matplotlib/matplotlib/blob
...                 #   /ec1cd2567521c105a451ce15e06de10715f8b54d/lib
...                 #   /matplotlib/colorbar.py#L838
...                 fraction = 0.15
...                 pb = self.axes.get_position(original=True).frozen()
...                 pad = 0.05
...                 x1 = 1.0-fraction
...                 pb1, pbx, pbcb = pb.splitx(x1-pad, x1)
...                 panchor = (1.0, 0.5)
...                 self.axes.set_position(pb1)
...                 self.axes.set_anchor(panchor)
...
...                 # make the gnomon
...                 fig = self.axes.get_figure()
...                 self.gnomon = fig.add_axes([0.85, 0.425, 0.15, 0.15], polar=True)
...                 self.gnomon.set_thetagrids([180, 270, 0, 90],
...                                            [r"$\pm\pi$", r"$-\frac{\pi}{2}$", "$0$", r"$+\frac{\pi}{2}$"],
...                                            frac=1.3)
...                 self.gnomon.set_theta_zero_location("N")
...                 self.gnomon.set_theta_direction(-1)
...                 self.gnomon.set_rgrids([1.], [""])
...                 N = 100
...                 theta = numerix.arange(-numerix.pi, numerix.pi, 2 * numerix.pi / N)
...                 radii = numerix.ones((N,))
...                 bars = self.gnomon.bar(theta, radii, width=2 * numerix.pi / N, bottom=0.0)
...                 colors = self._orientation_and_phase_to_rgb(orientation=numerix.array([theta]), phase=1.)
...                 for c, t, bar in zip(colors[0], theta, bars):
...                     bar.set_facecolor(c)
...                     bar.set_edgecolor(c)
...
...             def _reshape(self, var):
...                 '''return values of var in an 2D array'''
...                 return numerix.reshape(numerix.array(var),
...                                        var.mesh.shape[::-1])[::-1]
...
...             @staticmethod
...             def _orientation_and_phase_to_rgb(orientation, phase):
...                 from matplotlib import colors
...
...                 hsv = numerix.empty(orientation.shape + (3,))
...                 hsv[..., 0] = (orientation / numerix.pi + 1) / 2.
...                 hsv[..., 1] = 1.
...                 hsv[..., 2] = phase
...
...                 return colors.hsv_to_rgb(hsv)
...
...             @property
...             def _data(self):
...                 '''convert phase and orientation to rgb image array
...
...                 orientation (-pi, pi) -> hue (0, 1)
...                 phase (0, 1) -> value (0, 1)
...                 '''
...                 orientation = self._reshape(self.vars[0])
...                 phase = self._reshape(self.phase)
...
...                 return self._orientation_and_phase_to_rgb(orientation, phase)
...
...             def _plot(self):
...                 self.image.set_data(self._data)
...
...         from matplotlib import pyplot
...         pyplot.ion()
...         w, h = pyplot.figaspect(1.)
...         fig = pyplot.figure(figsize=(2*w, h))
...         timer = fig.text(0.1, 0.9, "t = %.3f" % 0, fontsize=18)
...
...         viewer = MultiViewer(viewers=(MatplotlibViewer(vars=dT,
...                                                        cmap=pyplot.cm.hot,
...                                                        datamin=-0.5,
...                                                        datamax=0.5,
...                                                        axes=fig.add_subplot(121)),
...                                       OrientationViewer(phase=phase,
...                                                         orientation=theta,
...                                                         title=theta.name,
...                                                         axes=fig.add_subplot(122))))
...     except ImportError:
...         viewer = MultiViewer(viewers=(Viewer(vars=dT,
...                                              datamin=-0.5,
...                                              datamax=0.5),
...                                       Viewer(vars=phase,
...                                              datamin=0.,
...                                              datamax=1.),
...                                       Viewer(vars=theta,
...                                              datamin=-numerix.pi,
...                                              datamax=numerix.pi)))
...     viewer.plot()

and iterate the solution in time, plotting as we go,

>>> eq = thetaEq & phaseEq & heatEq

>>> if __name__ == "__main__":
...     total_time = 2.
... else:
...     total_time = dt * 10
>>> elapsed = 0.
>>> save_interval = 0.002
>>> save_at = save_interval

>>> while elapsed < total_time:
...     if elapsed > 0.3:
...         q.value = 100
...     phase.updateOld()
...     dT.updateOld()
...     theta.updateOld()
...     eq.solve(dt=dt)
...     elapsed += dt
...     if __name__ == "__main__" and elapsed >= save_at:
...         timer.set_text("t = %.3f" % elapsed)
...         viewer.plot()
...         save_at += save_interval

The non-uniform temperature results from the release of latent heat at the solidifying interface. The dendrite arms grow fastest where the temperature gradient is steepest.

examples.phase.quaternary module¶

Solve a phase-field evolution and diffusion of four species in one-dimension.

The same procedure used to construct the two-component phase field diffusion problem in examples.phase.binary can be used to build up a system of multiple components. Once again, we’ll focus on 1D.

>>> from fipy import CellVariable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, PowerLawConvectionTerm, DefaultAsymmetricSolver, Viewer
>>> from fipy.tools import numerix

>>> nx = 400
>>> dx = 0.01
>>> L = nx * dx
>>> mesh = Grid1D(dx = dx, nx = nx)

We consider a free energy density $f(\phi, C_0,\ldots,C_N, T)$ that is a function of phase $\phi$

>>> phase = CellVariable(mesh=mesh, name='phase', value=1., hasOld=1)

interstitial components $C_0 \ldots C_M$

>>> interstitials = [
...     CellVariable(mesh=mesh, name='C0', hasOld=1)
... ]

substitutional components $C_{M+1} \ldots C_{N-1}$

>>> substitutionals = [
...     CellVariable(mesh=mesh, name='C1', hasOld=1),
...     CellVariable(mesh=mesh, name='C2', hasOld=1),
... ]

a “solvent” $C_N$ that is constrained by the concentrations of the other substitutional species, such that $C_N = 1 - \sum_{j=M}^{N-1} C_j$ ,

>>> solvent = 1
>>> for Cj in substitutionals:
...     solvent -= Cj
>>> solvent.name = 'CN'

and temperature $T$

>>> T = 1000

The free energy density of such a system can be written as

$f(\phi, C_0, \ldots, C_N, T) &= \sum_{j=0}^N C_j \left[ \mu^\circ_j(\phi, T) + R T \ln \frac{C_j}{\rho} \right]$

where

>>> R = 8.314 # J / (mol K)

is the gas constant. As in the binary case,

$\mu^\circ_j(\phi, T) = p(\phi) \mu_j^{\circ S}(T) + \left(1 - p(\phi)\right) \mu_j^{\circ L}(T) + \frac{W_j}{2} g(\phi)$

is constructed with the free energies of the pure components in each phase, given the “tilting” function

>>> def p(phi):
...     return phi**3 * (6 * phi**2 - 15 * phi + 10)

and the “double well” function

>>> def g(phi):
...     return (phi * (1 - phi))**2

We consider a very simplified model that has partial molar volumes $\bar{V}_0 = \cdots = \bar{V}_{M} = 0$ for the “interstitials” and $\bar{V}_{M+1} = \cdots = \bar{V}_{N} = 1$ for the “substitutionals”. This approximation has been used in a number of models where density effects are ignored, including the treatment of electrons in electrodeposition processes [25] [26]. Under these constraints

$\frac{\partial f}{\partial \phi} &= \sum_{j=0}^N C_j \frac{\partial f_j}{\partial \phi} \nonumber \\ &= \sum_{j=0}^N C_j \left[ \mu_j^{\circ SL}(T) p'(\phi) + \frac{W_j}{2} g'(\phi) \right] \\ \frac{\partial f}{\partial C_j} &= \left[\mu^\circ_j(\phi, T) + R T \ln \frac{C_j}{\rho} \right] \nonumber \\ &= \mu_j(\phi, C_j , T) \qquad\text{for $ j = 0\ldots M $}$

and

$\frac{\partial f}{\partial C_j} &= \left[\mu^\circ_j(\phi, T) + R T \ln \frac{C_j}{\rho} \right] - \left[\mu^\circ_N(\phi, T) + R T \ln \frac{C_N}{\rho} \right] \nonumber \\ &= \left[\mu_j(\phi, C_j, T) - \mu_N(\phi, C_N, T) \right] \qquad\text{for $ j = M+1\ldots N-1 $}$

where $\mu_j^{\circ SL}(T) \equiv \mu_j^{\circ S}(T) - \mu_j^{\circ L}(T)$ and where $\mu_j$ is the classical chemical potential of component $j$ for the binary species and $\rho = 1 + \sum_{j=0}^{M} C_j$ is the total molar density.

>>> rho = 1.
>>> for Cj in interstitials:
...     rho += Cj

$p'(\phi)$ and $g'(\phi)$ are the partial derivatives of of $p$ and $g$ with respect to $\phi$

>>> def pPrime(phi):
...     return 30. * g(phi)

>>> def gPrime(phi):
...     return 2. * phi * (1 - phi) * (1 - 2 * phi)

We “cook” the standard potentials to give the desired solid and liquid concentrations, with a solid phase rich in interstitials and the solvent and a liquid phase rich in the two substitutional species.

>>> interstitials[0].S = 0.3
>>> interstitials[0].L = 0.4
>>> substitutionals[0].S = 0.4
>>> substitutionals[0].L = 0.3
>>> substitutionals[1].S = 0.2
>>> substitutionals[1].L = 0.1
>>> solvent.S = 1.
>>> solvent.L = 1.
>>> for Cj in substitutionals:
...     solvent.S -= Cj.S
...     solvent.L -= Cj.L

>>> rhoS = rhoL = 1.
>>> for Cj in interstitials:
...     rhoS += Cj.S
...     rhoL += Cj.L

>>> for Cj in interstitials + substitutionals + [solvent]:
...     Cj.standardPotential = R * T * (numerix.log(Cj.L/rhoL)
...                                     - numerix.log(Cj.S/rhoS))

>>> for Cj in interstitials:
...     Cj.diffusivity = 1.
...     Cj.barrier = 0.

>>> for Cj in substitutionals:
...     Cj.diffusivity = 1.
...     Cj.barrier = R * T

>>> solvent.barrier = R * T

We create the phase equation

$\frac{1}{M_\phi}\frac{\partial \phi}{\partial t} = \kappa_\phi \nabla^2 \phi - \sum_{j=0}^N C_j \left[ \mu_j^{\circ SL}(T) p'(\phi) + \frac{W_j}{2} g'(\phi) \right]$

with a semi-implicit source just as in examples.phase.simple and examples.phase.binary

>>> enthalpy = 0.
>>> barrier = 0.
>>> for Cj in interstitials + substitutionals + [solvent]:
...     enthalpy += Cj * Cj.standardPotential
...     barrier += Cj * Cj.barrier

>>> mPhi = -((1 - 2 * phase) * barrier + 30 * phase * (1 - phase) * enthalpy)
>>> dmPhidPhi = 2 * barrier - 30 * (1 - 2 * phase) * enthalpy
>>> S1 = dmPhidPhi * phase * (1 - phase) + mPhi * (1 - 2 * phase)
>>> S0 = mPhi * phase * (1 - phase) - S1 * phase

>>> phase.mobility = 1.
>>> phase.gradientEnergy = 25
>>> phase.equation = TransientTerm(coeff=1/phase.mobility) \
...   == DiffusionTerm(coeff=phase.gradientEnergy) \
...      + S0 + ImplicitSourceTerm(coeff = S1)

We could construct the diffusion equations one-by-one, in the manner of examples.phase.binary, but it is better to take advantage of the full scripting power of the Python language, where we can easily loop over components or even make “factory” functions if we desire. For the interstitial diffusion equations, we arrange in canonical form as before:

$\underbrace{ \frac{\partial C_j}{\partial t} \vphantom{\left\{ \overbrace{ \left[ \mu_j^{\circ SL} \nabla p(\phi) \right] }^{\text{phase transformation}} \right\}} }_{\text{transient}} &= \underbrace{ D_j\nabla^2 C_j \vphantom{\left\{ \overbrace{ \left[ \mu_j^{\circ SL} \nabla p(\phi) \right] }^{\text{phase transformation}} \right\}} }_{\text{diffusion}} \\ & \qquad + \underbrace{ D_j\nabla\cdot \frac{C_j}{1 + \sum_{\substack{k=0\\ k \neq j}}^{M} C_k} \left\{ \overbrace{ \frac{\rho}{R T} \left[ \mu_j^{\circ SL} \nabla p(\phi) + \frac{W_j}{2} \nabla g(\phi) \right] }^{\text{phase transformation}} - \overbrace{ \sum_{\substack{i=0\\ i \neq j}}^{M} \nabla C_i }^{\text{counter diffusion}} \right\} }_{\text{convection}}$

>>> for Cj in interstitials:
...     phaseTransformation = (rho.harmonicFaceValue / (R * T)) \
...       * (Cj.standardPotential * p(phase).faceGrad
...          + 0.5 * Cj.barrier * g(phase).faceGrad)
...
...     CkSum = CellVariable(mesh=mesh, value=0.)
...     for Ck in [Ck for Ck in interstitials if Ck is not Cj]:
...         CkSum += Ck
...
...     counterDiffusion = CkSum.faceGrad
...
...     convectionCoeff = counterDiffusion + phaseTransformation
...     convectionCoeff *= (Cj.diffusivity
...                         / (1. + CkSum.harmonicFaceValue))
...
...     Cj.equation = (TransientTerm()
...                    == DiffusionTerm(coeff=Cj.diffusivity)
...                    + PowerLawConvectionTerm(coeff=convectionCoeff))

The canonical form of the substitutional diffusion equations is

$\underbrace{ \frac{\partial C_j}{\partial t} }_{\text{transient}} &= \underbrace{ D_{j}\nabla^2 C_j \vphantom{\frac{\partial C_j}{\partial t}} }_{\text{diffusion}} \\ & \qquad + \underbrace{ D_{j}\nabla\cdot \frac{C_j}{1 - \sum_{\substack{k=M+1\\ k \neq j}}^{N-1} C_k} \left\{ \overbrace{ \frac{C_N}{R T} \left[ \left(\mu_j^{\circ SL} - \mu_N^{\circ SL}\right) \nabla p(\phi) + \frac{W_j - W_N}{2} \nabla g(\phi) \right] \vphantom{\sum_{\substack{i=M+1\\ i \neq j}}^{N-1} \nabla C_i} }^{\text{phase transformation}} + \overbrace{ \sum_{\substack{i=M+1\\ i \neq j}}^{N-1} \nabla C_i }^{\text{counter diffusion}} \right\} }_{\text{convection}}$

>>> for Cj in substitutionals:
...     phaseTransformation = (solvent.harmonicFaceValue / (R * T)) \
...       * ((Cj.standardPotential - solvent.standardPotential) * p(phase).faceGrad
...          + 0.5 * (Cj.barrier - solvent.barrier) * g(phase).faceGrad)
...
...     CkSum = CellVariable(mesh=mesh, value=0.)
...     for Ck in [Ck for Ck in substitutionals if Ck is not Cj]:
...         CkSum += Ck
...
...     counterDiffusion = CkSum.faceGrad
...
...     convectionCoeff = counterDiffusion + phaseTransformation
...     convectionCoeff *= (Cj.diffusivity
...                         / (1. - CkSum.harmonicFaceValue))
...
...     Cj.equation = (TransientTerm()
...                    == DiffusionTerm(coeff=Cj.diffusivity)
...                    + PowerLawConvectionTerm(coeff=convectionCoeff))

We start with a sharp phase boundary

$\xi = \begin{cases} 1& \text{for $x \le L/2$,} \\ 0& \text{for $x > L/2$,} \end{cases}$

>>> x = mesh.cellCenters[0]
>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L / 2)

and with uniform concentration fields, initially equal to the average of the solidus and liquidus concentrations

>>> for Cj in interstitials + substitutionals:
...     Cj.setValue((Cj.S + Cj.L) / 2.)

If we’re running interactively, we create a viewer

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=([phase]
...                           + interstitials + substitutionals
...                           + [solvent]),
...                     datamin=0, datamax=1)
...     viewer.plot()

and again iterate to equilibrium

>>> solver = DefaultAsymmetricSolver(tolerance=1e-10)

>>> dt = 10000
>>> from builtins import range
>>> for i in range(5):
...     for field in [phase] + substitutionals + interstitials:
...         field.updateOld()
...     phase.equation.solve(var = phase, dt = dt)
...     for field in substitutionals + interstitials:
...         field.equation.solve(var = field,
...                              dt = dt,
...                              solver = solver)
...     if __name__ == '__main__':
...         viewer.plot()

phase and four composition fields in equilibrium

We can confirm that the far-field phases have remained separated

>>> X = mesh.faceCenters[0]
>>> print(numerix.allclose(phase.faceValue[X.value==0], 1.0, rtol = 1e-5, atol = 1e-5))
True
>>> print(numerix.allclose(phase.faceValue[X.value==L], 0.0, rtol = 1e-5, atol = 1e-5))
True

and that the concentration fields have appropriately segregated into their equilibrium values in each phase

>>> equilibrium = True
>>> for Cj in interstitials + substitutionals:
...     equilibrium &= numerix.allclose(Cj.faceValue[X.value==0], Cj.S, rtol = 3e-3, atol = 3e-3).value
...     equilibrium &= numerix.allclose(Cj.faceValue[X.value==L], Cj.L, rtol = 3e-3, atol = 3e-3).value
>>> print(equilibrium)
True

examples.phase.simple module¶

Solve a phase-field (Allen-Cahn) problem in one-dimension.

To run this example from the base FiPy directory, type python examples/phase/simple/input.py at the command line. A viewer object should appear and, after being prompted to step through the different examples, the word finished in the terminal.

This example takes the user through assembling a simple problem with FiPy. It describes a steady 1D phase field problem with no-flux boundary conditions such that,

(13)¶ $\frac{1}{M_\phi}\frac{\partial \phi}{\partial t} = \kappa_\phi \nabla^2\phi - \frac{\partial f}{\partial \phi}$

For solidification problems, the Helmholtz free energy is frequently given by

$f(\phi,T) = \frac{W}{2}g(\phi) + L_v\frac{T-T_M}{T_M}p(\phi)$

where $W$ is the double-well barrier height between phases, $L_v$ is the latent heat, $T$ is the temperature, and $T_M$ is the melting point.

One possible choice for the double-well function is

$g(\phi) = \phi^2(1 - \phi)^2$

and for the interpolation function is

$p(\phi) = \phi^3(6\phi^2 - 15\phi + 10).$

We create a 1D solution mesh

>>> from fipy import CellVariable, Variable, Grid1D, DiffusionTerm, TransientTerm, ImplicitSourceTerm, DummySolver, Viewer
>>> from fipy.tools import numerix

>>> L = 1.
>>> nx = 400
>>> dx = L / nx

>>> mesh = Grid1D(dx = dx, nx = nx)

We create the phase field variable

>>> phase = CellVariable(name = "phase",
...                      mesh = mesh)

and set a step-function initial condition

$\phi = \begin{cases} 1 & \text{for $x \le L/2$} \\ 0 & \text{for $x > L/2$} \end{cases} \quad\text{at $t = 0$}$

>>> x = mesh.cellCenters[0]
>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L/2)

If we are running interactively, we’ll want a viewer to see the results

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer = Viewer(vars = (phase,))
...     viewer.plot()
...     input("Initial condition. Press <return> to proceed...")

We choose the parameter values,

>>> kappa = 0.0025
>>> W = 1.
>>> Lv = 1.
>>> Tm = 1.
>>> T = Tm
>>> enthalpy = Lv * (T - Tm) / Tm

We build the equation by assembling the appropriate terms. Since, with $T = T_M$ we are interested in a steady-state solution, we omit the transient term $(1/M_\phi)\frac{\partial \phi}{\partial t}$ .

The analytical solution for this steady-state phase field problem, in an infinite domain, is

(14)¶ $\phi = \frac{1}{2}\left[1 - \tanh\frac{x-L/2}{2\sqrt{\kappa/W}}\right]$

or

>>> x = mesh.cellCenters[0]
>>> analyticalArray = 0.5*(1 - numerix.tanh((x - L/2)/(2*numerix.sqrt(kappa/W))))

We treat the diffusion term $\kappa_\phi \nabla^2\phi$ implicitly,

Note

“Diffusion” in FiPy is not limited to the movement of atoms, but rather refers to the spontaneous spreading of any quantity (e.g., solute, temperature, or in this case “phase”) by flow “down” the gradient of that quantity.

The source term is

$S = -\frac{\partial f}{\partial \phi} &= -\frac{W}{2}g'(\phi) - L\frac{T-T_M}{T_M}p'(\phi) \\ &= -\left[W\phi(1-\phi)(1-2\phi) + L\frac{T-T_M}{T_M}30\phi^2(1-\phi)^2\right] \\ &= m_\phi\phi(1-\phi)$

where $m_\phi \equiv -[W(1-2\phi) + 30\phi(1-\phi)L\frac{T-T_M}{T_M}]$ .

The simplest approach is to add this source explicitly

>>> mPhi = -((1 - 2 * phase) * W + 30 * phase * (1 - phase) * enthalpy)
>>> S0 = mPhi * phase * (1 - phase)
>>> eq = S0 + DiffusionTerm(coeff=kappa)

After solving this equation

>>> eq.solve(var = phase, solver=DummySolver())

we obtain the surprising result that $\phi$ is zero everywhere.

>>> print(phase.allclose(analyticalArray, rtol = 1e-4, atol = 1e-4))
0
>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Fully explicit source. Press <return> to proceed...")

steady-state phase field zero everywhere

On inspection, we can see that this occurs because, for our step-function initial condition, $m_\phi = 0$ everywhere, hence we are actually only solving the simple implicit diffusion equation $\kappa_\phi \nabla^2\phi = 0$ , which has exactly the uninteresting solution we obtained.

The resolution to this problem is to apply relaxation to obtain the desired answer, i.e., the solution is allowed to relax in time from the initial condition to the desired equilibrium solution. To do so, we reintroduce the transient term from Equation (1)

>>> eq = TransientTerm() == DiffusionTerm(coeff=kappa) + S0

>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L/2)

>>> from builtins import range
>>> for i in range(13):
...     eq.solve(var = phase, dt=1.)
...     if __name__ == '__main__':
...         viewer.plot()

After 13 time steps, the solution has converged to the analytical solution

>>> print(phase.allclose(analyticalArray, rtol = 1e-4, atol = 1e-4))
1
>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Relaxation, explicit. Press <return> to proceed...")

steady-state equilibrium phase field after relaxation

Note

The solution is only found accurate to $\approx 4.3\times 10^{-5}$ because the infinite-domain analytical solution (2) is not an exact representation for the solution in a finite domain of length $L$ .

Setting fixed-value boundary conditions of 1 and 0 would still require the relaxation method with the fully explicit source.

Solution performance can be improved if we exploit the dependence of the source on $\phi$ . By doing so, we can make the source semi-implicit, improving the rate of convergence over the fully explicit approach. The source can only be semi-implicit because we employ sparse linear algebra routines to solve the PDEs, i.e., there is no fully implicit way to represent a term like $\phi^4$ in the linear set of equations $\mathsf{M} \vec{\phi} - \vec{b} = 0$ .

By linearizing a source as $S = S_0 - S_1 \phi$ , we make it more implicit by adding the coefficient $S_1$ to the matrix diagonal. For numerical stability, this linear coefficient must never be negative.

There are an infinite number of choices for this linearization, but many do not converge very well. One choice is that used by Ryo Kobayashi:

>>> S0 = mPhi * phase * (mPhi > 0)
>>> S1 = mPhi * ((mPhi < 0) - phase)
>>> eq = DiffusionTerm(coeff=kappa) + S0 \
...   + ImplicitSourceTerm(coeff = S1)

Note

Because mPhi is a variable field, the quantities (mPhi > 0) and (mPhi < 0) evaluate to variable fields of True and False, instead of single Boolean values.

This expression converges to the same value given by the explicit relaxation approach, but in only 8 sweeps (note that because there is no transient term, these sweeps are not time steps, but rather repeated iterations at the same time step to reach a converged solution).

Note

We use solve() instead of sweep() because we don’t care about the residual. Either function would work, but solve() is a bit faster.

>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L/2)

>>> from builtins import range
>>> for i in range(8):
...     eq.solve(var = phase)
>>> print(phase.allclose(analyticalArray, rtol = 1e-4, atol = 1e-4))
1
>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Kobayashi, semi-implicit. Press <return> to proceed...")

In general, the best convergence is obtained when the linearization gives a good representation of the relationship between the source and the dependent variable. The best practical advice is to perform a Taylor expansion of the source about the previous value of the dependent variable such that $S = S_\text{old} + \left.\frac{\partial S}{\partial \phi}\right\rvert_\text{old} (\phi - \phi_\text{old}) = (S - \frac{\partial S}{\partial\phi} \phi)_\text{old} + \left.\frac{\partial S}{\partial \phi}\right|_\text{old} \phi$ . Now, if our source term is represented by $S = S_0 + S_1 \phi$ , then $S_1 = \left.\frac{\partial S}{\partial \phi}\right|_\text{old}$ and $S_0 = (S - \frac{\partial S}{\partial\phi} \phi)_\text{old} = S_\text{old} - S_1 \phi_\text{old}$ . In this way, the linearized source will be tangent to the curve of the actual source as a function of the dependent variable.

For our source, $S = m_\phi \phi (1 - \phi)$ ,

$\frac{\partial S}{\partial \phi} = \frac{\partial m_\phi}{\partial \phi} \phi (1 - \phi) + m_\phi (1 - 2\phi)$

and

$\frac{\partial m_\phi}{\partial \phi} = 2 W - 30 (1 - 2\phi) L\frac{T-T_M}{T_M},$

or

>>> dmPhidPhi = 2 * W - 30 * (1 - 2 * phase) * enthalpy
>>> S1 = dmPhidPhi * phase * (1 - phase) + mPhi * (1 - 2 * phase)
>>> S0 = mPhi * phase * (1 - phase) - S1 * phase
>>> eq = DiffusionTerm(coeff=kappa) + S0 \
...   + ImplicitSourceTerm(coeff = S1)

Using this scheme, where the coefficient of the implicit source term is tangent to the source, we reach convergence in only 5 sweeps

>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L/2)

>>> from builtins import range
>>> for i in range(5):
...     eq.solve(var = phase)
>>> print(phase.allclose(analyticalArray, rtol = 1e-4, atol = 1e-4))
1
>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Tangent, semi-implicit. Press <return> to proceed...")

Although, for this simple problem, there is no appreciable difference in run-time between the fully explicit source and the optimized semi-implicit source, the benefit of 60% fewer sweeps should be obvious for larger systems and longer iterations.

This example has focused on just the region of the phase field interface in equilibrium. Problems of interest, though, usually involve the dynamics of one phase transforming to another. To that end, let us recast the problem using physical parameters and dimensions. We’ll need a new mesh

>>> nx = 400
>>> dx = 5e-6 # cm
>>> L = nx * dx

>>> mesh = Grid1D(dx = dx, nx = nx)

and thus must redeclare $\phi$ on the new mesh

>>> phase = CellVariable(name="phase",
...                      mesh=mesh,
...                      hasOld=1)
>>> x = mesh.cellCenters[0]
>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L/2)

We choose the parameter values appropriate for nickel, given in [23]

>>> Lv = 2350 # J / cm**3
>>> Tm = 1728. # K
>>> T = Variable(value=Tm)
>>> enthalpy = Lv * (T - Tm) / Tm # J / cm**3

The parameters of the phase field model can be related to the surface energy $\sigma$ and the interfacial thickness $\delta$ by

$\kappa &= 6\sigma\delta \\ W &= \frac{6\sigma}{\delta} \\ M_\phi &= \frac{T_m\beta}{6 L \delta}.$

We take $\delta \approx \Delta x$ .

>>> delta = 1.5 * dx
>>> sigma = 3.7e-5 # J / cm**2
>>> beta = 0.33 # cm / (K s)
>>> kappa = 6 * sigma * delta # J / cm
>>> W = 6 * sigma / delta # J / cm**3
>>> Mphi = Tm * beta / (6. * Lv * delta) # cm**3 / (J s)

>>> if __name__ == '__main__':
...     displacement = L * 0.1
... else:
...     displacement = L * 0.025

>>> analyticalArray = CellVariable(name="tanh", mesh=mesh,
...                                value=0.5 * (1 - numerix.tanh((x - (L / 2. + displacement))
...                                                              / (2 * delta))))

and make a new viewer

>>> if __name__ == '__main__':
...     viewer2 = Viewer(vars = (phase, analyticalArray))
...     viewer2.plot()

Now we can redefine the transient phase field equation, using the optimal form of the source term shown above

>>> mPhi = -((1 - 2 * phase) * W + 30 * phase * (1 - phase) * enthalpy)
>>> dmPhidPhi = 2 * W - 30 * (1 - 2 * phase) * enthalpy
>>> S1 = dmPhidPhi * phase * (1 - phase) + mPhi * (1 - 2 * phase)
>>> S0 = mPhi * phase * (1 - phase) - S1 * phase
>>> eq = TransientTerm(coeff=1/Mphi) == DiffusionTerm(coeff=kappa) \
...                         + S0 + ImplicitSourceTerm(coeff = S1)

In order to separate the effect of forming the phase field interface from the kinetics of moving it, we first equilibrate at the melting point. We now use the sweep() method instead of solve() because we require the residual.

>>> timeStep = 1e-6
>>> from builtins import range
>>> for i in range(10):
...     phase.updateOld()
...     res = 1e+10
...     while res > 1e-5:
...         res = eq.sweep(var=phase, dt=timeStep)
>>> if __name__ == '__main__':
...     viewer2.plot()

and then quench by 1 K

>>> T.setValue(T() - 1)

In order to have a stable numerical solution, the interface must not move more than one grid point per time step, we thus set the timestep according to the grid spacing $\Delta x$ , the linear kinetic coefficient $\beta$ , and the undercooling $\abs{T_m - T}$ Again we use the sweep() method as a replacement for solve().

>>> velocity = beta * abs(Tm - T()) # cm / s
>>> timeStep = .1 * dx / velocity # s
>>> elapsed = 0
>>> while elapsed < displacement / velocity:
...     phase.updateOld()
...     res = 1e+10
...     while res > 1e-5:
...         res = eq.sweep(var=phase, dt=timeStep)
...     elapsed += timeStep
...     if __name__ == '__main__':
...         viewer2.plot()

A hyperbolic tangent is not an exact steady-state solution given the quintic polynomial we chose for the $p$ function, but it gives a reasonable approximation.

>>> print(phase.allclose(analyticalArray, rtol = 5, atol = 2e-3))
1

If we had made another common choice of $p(\phi) = \phi^2(3 - 2\phi)$ , we would have found much better agreement, as that case does give an exact $\tanh$ solution in steady state. If SciPy is available, another way to compare against the expected result is to do a least-squared fit to determine the interface velocity and thickness

>>> try:
...     def tanhResiduals(p, y, x, t):
...         V, d = p
...         return y - 0.5 * (1 - numerix.tanh((x - V * t - L / 2.) / (2*d)))
...     from scipy.optimize import leastsq 
...     x =  mesh.cellCenters[0]
...     (V_fit, d_fit), msg = leastsq(tanhResiduals, [L/2., delta],
...                                   args=(phase.globalValue, x.globalValue, elapsed)) 
... except ImportError:
...     V_fit = d_fit = 0
...     print("The SciPy library is unavailable to fit the interface \
... thickness and velocity")

>>> print(abs(1 - V_fit / velocity) < 4.1e-2) 
True
>>> print(abs(1 - d_fit / delta) < 2e-2) 
True

>>> from fipy import input
>>> if __name__ == '__main__':
...     input("Dimensional, semi-implicit. Press <return> to proceed...")

phase field when solved with physical parameters

examples.phase.symmetry module¶

This example creates four symmetric quadrilateral regions in a box. We start with a CellVariable object that contains the following values:

$\phi(x, y) = x y, 0 \le x \le L, 0 \le y \le L$

We wish to create 4 symmetric regions such that

$\phi(x, y) = \phi(L - x, y) = \phi(L - x, y) = \phi(L - x, L - y), 0 \le x \le L / 2, 0 \le y \le L / 2$

We create a square domain

>>> from fipy import CellVariable, Grid2D, Viewer
>>> from fipy.tools import numerix

>>> N = 20
>>> L = 1.
>>> dx = L / N
>>> dy = L / N

>>> mesh = Grid2D(
...    dx = dx,
...    dy = dy,
...    nx = N,
...    ny = N)

>>> var = CellVariable(name = "test", mesh = mesh)

First set the values as given in the above equation:

>>> x, y = mesh.cellCenters
>>> var.setValue(x * y)

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=var, datamin=0, datamax=L * L / 4.)
...     viewer.plot()

The bottom-left quadrant is mirrored into each of the other three quadrants

>>> q = (x > L / 2.) & (y < L / 2.)
>>> var[q] = var(((L - x)[q],       y[q]))
>>> q = (x < L / 2.) & (y > L / 2.)
>>> var[q] = var((      x[q], (L - y)[q]))
>>> q = (x > L / 2.) & (y > L / 2.)
>>> var[q] = var(((L - x)[q], (L - y)[q]))

>>> if __name__ == '__main__':
...     viewer.plot()

The following code tests the results with a different algorithm:

>>> testResult = numerix.zeros((N // 2, N // 2), 'd')
>>> bottomRight = numerix.zeros((N // 2, N // 2), 'd')
>>> topLeft = numerix.zeros((N // 2, N // 2), 'd')
>>> topRight = numerix.zeros((N // 2, N // 2), 'd')
>>> from builtins import range
>>> for j in range(N // 2):
...     for i in range(N // 2):
...         x = dx * (i + 0.5)
...         y = dx * (j + 0.5)
...         testResult[i, j] = x * y
...         bottomRight[i, j] = var(((L - x,), (y,)))[0]
...         topLeft[i, j] = var(((x,), (L - y,)))[0]
...         topRight[i, j] = var(((L - x,), (L - y,)))[0]
>>> numerix.allclose(testResult, bottomRight, atol = 1e-10)
1
>>> numerix.allclose(testResult, topLeft, atol = 1e-10)
1
>>> numerix.allclose(testResult, topRight, atol = 1e-10)
1

examples.phase.test module¶

Last updated on Jun 27, 2023. Created using Sphinx 6.2.1.