examples.phase.impingement.mesh40x1

Solve for the impingement of two grains in one dimension.

In this example we solve a coupled phase and orientation equation on a one dimensional grid. This is another aspect of the model of Warren, Kobayashi, Lobkovsky and Carter [WarrenPolycrystal]

>>> from fipy import *

>>> nx = 40
>>> Lx = 2.5 * nx / 100.
>>> dx = Lx / nx
>>> mesh = Grid1D(dx=dx, nx=nx)

This problem simulates the wet boundary that forms between grains of different orientations. The phase equation is given by

where

and the orientation equation is given by

where

The initial conditions for this problem are set such that for and

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

The parameters for these equations are

>>> timeStepDuration = 0.02
>>> phaseTransientCoeff = 0.1
>>> thetaSmallValue = 1e-6
>>> beta = 1e5
>>> mu = 1e3
>>> thetaTransientCoeff = 0.01
>>> gamma= 1e3
>>> epsilon = 0.008
>>> s = 0.01
>>> alpha = 0.015

The system is held isothermal at

>>> temperature = 1.

and is initially solid everywhere

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

Because theta is an -valued variable (i.e. it maps to the circle) and thus intrinsically has -peridocity, we must use ModularVariable instead of a CellVariable. A ModularVariable confines theta to by adding or subtracting where necessary and by defining a new subtraction operator between two angles.

>>> theta = ModularVariable(
...     name='theta',
...     mesh=mesh,
...     value=1.,
...     hasOld=1
...     )

The left and right halves of the domain are given different orientations.

>>> theta.setValue(0., where=mesh.cellCenters[0] > Lx / 2.)

The phase equation is built in the following way.

>>> mPhiVar = phase - 0.5 + temperature * phase * (1 - phase)

The source term is linearized in the manner demonstrated in examples.phase.simple (Kobayashi, semi-implicit).

>>> thetaMag = theta.old.grad.mag
>>> implicitSource = mPhiVar * (phase - (mPhiVar < 0))
>>> implicitSource += (2 * s + epsilon**2 * thetaMag) * thetaMag

The phase equation is constructed.

>>> phaseEq = TransientTerm(phaseTransientCoeff) \
...   == ExplicitDiffusionTerm(alpha**2) \
...      - ImplicitSourceTerm(implicitSource) \
...      + (mPhiVar > 0) * mPhiVar * phase

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.

>>> phaseMod = phase + ( phase < thetaSmallValue ) * thetaSmallValue
>>> phaseModSq = phaseMod * phaseMod
>>> expo = epsilon * beta * theta.grad.mag
>>> expo = (expo < 100.) * (expo - 100.) + 100.
>>> pFunc = 1. + numerix.exp(-expo) * (mu / epsilon - 1.)

>>> phaseFace = phase.arithmeticFaceValue
>>> phaseSq = phaseFace * phaseFace
>>> gradMag = theta.faceGrad.mag
>>> eps = 1. / gamma / 10.
>>> gradMag += (gradMag < eps) * eps
>>> IGamma = (gradMag > 1. / gamma) * (1 / gradMag - gamma) + gamma
>>> diffusionCoeff = phaseSq * (s * IGamma + epsilon**2)

The source term requires the evaluation of the face gradient without the modular operator. thetagetFaceGradNoMod() evelautes the gradient without modular arithmetic.

>>> thetaGradDiff = theta.faceGrad - theta.faceGradNoMod
>>> sourceCoeff = (diffusionCoeff * thetaGradDiff).divergence

Finally the theta equation can be constructed.

>>> thetaEq = TransientTerm(thetaTransientCoeff * phaseModSq * pFunc) == \
...           DiffusionTerm(diffusionCoeff) \
...           + sourceCoeff

If the example is run interactively, we create viewers for the phase and orientation variables.

>>> if __name__ == '__main__':
...     phaseViewer = Viewer(vars=phase, datamin=0., datamax=1.)
...     thetaProductViewer = Viewer(vars=theta,
...                                 datamin=-pi, datamax=pi)
...     phaseViewer.plot()
...     thetaProductViewer.plot()

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

>>> steps = 10
>>> for i in range(steps):
...     theta.updateOld()
...     thetaEq.solve(theta, dt = timeStepDuration)
...     phaseEq.solve(phase, dt = timeStepDuration)
...     if __name__ == '__main__':
...         phaseViewer.plot()
...         thetaProductViewer.plot()

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

>>> import os
>>> testData = numerix.loadtxt(os.path.splitext(__file__)[0] + '.gz')
>>> testData = CellVariable(mesh=mesh, value=testData)
>>> print theta.allclose(testData)
1