# 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 

>>> 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 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 -periodicity, 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 > 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
>>> 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. thetafaceGradNoMod 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__) + text_to_native_str('.gz'))
>>> testData = CellVariable(mesh=mesh, value=testData)
>>> print(theta.allclose(testData))
1

Last updated on Jun 15, 2022. Created using Sphinx 5.0.1.