examples.phase.impingement.mesh20x20¶

Solve for the impingement of four grains in two dimensions.

In the following examples, we solve the same set of equations as in examples.phase.impingement.mesh40x1 with different initial conditions and a 2D mesh:

>>> from fipy.tools.parser import parse

>>> numberOfElements = parse('--numberOfElements', action = 'store',
...                          type = 'int', default = 400)
>>> numberOfSteps = parse('--numberOfSteps', action = 'store',
...                       type = 'int', default = 10)

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

>>> steps = numberOfSteps
>>> N = int(numerix.sqrt(numberOfElements))
>>> L = 2.5 * N / 100.
>>> dL = L / N
>>> mesh = Grid2D(dx=dL, dy=dL, nx=N, ny=N)

The initial conditions are given by $\phi = 1$ and

$\theta = \begin{cases} \frac{2 \pi}{3} & \text{for $x^2 - y^2 < L / 2$,} \\ \frac{-2 \pi}{3} & \text{for $(x-L)^2 - y^2 < L / 2$,} \\ \frac{-2 \pi}{3}+0.3 & \text{for $x^2 - (y-L)^2 < L / 2$,} \\ \frac{2 \pi}{3} & \text{for $(x-L)^2 - (y-L)^2 < L / 2$.} \end{cases}$

This defines four solid regions with different orientations. Solidification occurs and then boundary wetting occurs where the orientation varies.

The parameters for this example 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 = 10.

and is initialized to liquid everywhere

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

The orientation is initialized to a uniform value to denote the randomly oriented liquid phase

>>> theta = ModularVariable(
...     name='theta',
...     mesh=mesh,
...     value=-numerix.pi + 0.0001,
...     hasOld=1
...     )

Four different solid circular domains are created at each corner of the domain with appropriate orientations

>>> x, y = mesh.cellCenters
>>> for a, b, thetaValue in ((0., 0.,  2. * numerix.pi / 3.),
...                          (L, 0., -2. * numerix.pi / 3.),
...                          (0., L, -2. * numerix.pi / 3. + 0.3),
...                          (L, L,  2. * numerix.pi / 3.)):
...     segment = (x - a)**2 + (y - b)**2 < (L / 2.)**2
...     phase.setValue(1., where=segment)
...     theta.setValue(thetaValue, where=segment)

The phase equation is built in the following way. The source term is linearized in the manner demonstrated in examples.phase.simple (Kobayashi, semi-implicit). Here we use a function to build the equation, so that it can be reused later.

>>> def buildPhaseEquation(phase, theta):
...
...     mPhiVar = phase - 0.5 + temperature * phase * (1 - phase)
...     thetaMag = theta.old.grad.mag
...     implicitSource = mPhiVar * (phase - (mPhiVar < 0))
...     implicitSource += (2 * s + epsilon**2 * thetaMag) * thetaMag
...
...     return TransientTerm(phaseTransientCoeff) == \
...               ExplicitDiffusionTerm(alpha**2) \
...               - ImplicitSourceTerm(implicitSource) \
...               + (mPhiVar > 0) * mPhiVar * phase

>>> phaseEq = buildPhaseEquation(phase, theta)

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. The source term requires the evaluation of the face gradient without the modular operators.

>>> def buildThetaEquation(phase, theta):
...
...     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)
...
...     thetaGradDiff = theta.faceGrad - theta.faceGradNoMod
...     sourceCoeff = (diffusionCoeff * thetaGradDiff).divergence
...
...     return TransientTerm(thetaTransientCoeff * phaseModSq * pFunc) == \
...                DiffusionTerm(diffusionCoeff) \
...                + sourceCoeff

>>> thetaEq = buildThetaEquation(phase, theta)

If the example is run interactively, we create viewers for the phase and orientation variables. Rather than viewing the raw orientation, which is not meaningful in the liquid phase, we weight the orientation by the phase

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

The solution will be tested against data that was created with steps = 10 with a FORTRAN code written by Ryo Kobayashi for phase field modeling. The following code opens the file mesh20x20.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')).flat

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

>>> from builtins import range
>>> for i in range(steps):
...     theta.updateOld()
...     thetaEq.solve(theta, dt=timeStepDuration, solver=GeneralSolver(iterations=2000, tolerance=1e-15))
...     phaseEq.solve(phase, dt=timeStepDuration, solver=GeneralSolver(iterations=2000, tolerance=1e-15))
...     if __name__ == '__main__':
...         phaseViewer.plot()
...         thetaProductViewer.plot()

The solution is compared against Ryo Kobayashi’s test data

>>> print(theta.allclose(testData, rtol=1e-7, atol=1e-7))
1

The following code shows how to restart a simulation from some saved data. First, reset the variables to their original values.

>>> phase.setValue(0)
>>> theta.setValue(-numerix.pi + 0.0001)
>>> x, y = mesh.cellCenters
>>> for a, b, thetaValue in ((0., 0.,  2. * numerix.pi / 3.),
...                          (L, 0., -2. * numerix.pi / 3.),
...                          (0., L, -2. * numerix.pi / 3. + 0.3),
...                          (L, L,  2. * numerix.pi / 3.)):
...     segment = (x - a)**2 + (y - b)**2 < (L / 2.)**2
...     phase.setValue(1., where=segment)
...     theta.setValue(thetaValue, where=segment)

Step through half the time steps.

>>> from builtins import range
>>> for i in range(steps // 2):
...     theta.updateOld()
...     thetaEq.solve(theta, dt=timeStepDuration, solver=GeneralSolver(iterations=2000, tolerance=1e-15))
...     phaseEq.solve(phase, dt=timeStepDuration, solver=GeneralSolver(iterations=2000, tolerance=1e-15))

We confirm that the solution has not yet converged to that given by Ryo Kobayashi’s FORTRAN code:

>>> print(theta.allclose(testData))
0

We save the variables to disk.

>>> (f, filename) = dump.write({'phase' : phase, 'theta' : theta}, extension = '.gz')

and then recall them to test the data pickling mechanism

>>> data = dump.read(filename, f)
>>> newPhase = data['phase']
>>> newTheta = data['theta']
>>> newThetaEq = buildThetaEquation(newPhase, newTheta)
>>> newPhaseEq = buildPhaseEquation(newPhase, newTheta)

and finish the iterations,

>>> from builtins import range
>>> for i in range(steps // 2):
...     newTheta.updateOld()
...     newThetaEq.solve(newTheta, dt=timeStepDuration, solver=GeneralSolver(iterations=2000, tolerance=1e-15))
...     newPhaseEq.solve(newPhase, dt=timeStepDuration, solver=GeneralSolver(iterations=2000, tolerance=1e-15))

The solution is compared against Ryo Kobayashi’s test data

>>> print(newTheta.allclose(testData, rtol=1e-7))
1

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