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 
>>> from fipy import CellVariable, ModularVariable, Grid1D, TransientTerm, DiffusionTerm, ExplicitDiffusionTerm, ImplicitSourceTerm, GeneralSolver, Viewer >>> from fipy.tools import numerix
>>> 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
and the orientation equation is given by
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. ... )
is an -valued variable (i.e. it maps to the circle) and thus
intrinsically has -periodicity,
we must use
ModularVariable instead of a
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 > Lx / 2.)
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
phase equation is constructed.
>>> phaseEq = TransientTerm(phaseTransientCoeff) \ ... == ExplicitDiffusionTerm(alpha**2) \ ... - ImplicitSourceTerm(implicitSource) \ ... + (mPhiVar > 0) * mPhiVar * phase
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
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.
evaluates the gradient without modular arithmetic.
>>> thetaGradDiff = theta.faceGrad - theta.faceGradNoMod >>> sourceCoeff = (diffusionCoeff * thetaGradDiff).divergence
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=-numerix.pi, datamax=numerix.pi) ... phaseViewer.plot() ... thetaProductViewer.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): ... 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
steps = 10 with a FORTRAN code written by Ryo Kobayashi for
phase field modeling. The following code opens the file
extracts the data and compares it with the
>>> import os >>> from future.utils import text_to_native_str >>> testData = numerix.loadtxt(os.path.splitext(__file__) + text_to_native_str('.gz')) >>> testData = CellVariable(mesh=mesh, value=testData) >>> print(theta.allclose(testData)) 1