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 [13]
>>> 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
where
and the orientation equation is given by
where
The initial conditions for this problem are set such that \(\phi = 1\) for \(0 \le x \le L_x\) 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 \(S^1\)-valued variable (i.e. it maps to the circle) and thus
intrinsically has \(2\pi\)-periodicity,
we must use ModularVariable
instead of a CellVariable
. A
ModularVariable
confines theta
to
\(-\pi < \theta \le \pi\) by adding or subtracting \(2\pi\) 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. theta
faceGradNoMod
evaluates 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=-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
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
>>> from future.utils import text_to_native_str
>>> testData = numerix.loadtxt(os.path.splitext(__file__)[0] + text_to_native_str('.gz'))
>>> testData = CellVariable(mesh=mesh, value=testData)
>>> print(theta.allclose(testData))
1