examples.phase.impingement package¶

Submodules¶

examples.phase.impingement.mesh20x20 module¶

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

examples.phase.impingement.mesh40x1 module¶

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

$\tau_{\phi} \frac{\partial \phi}{\partial t} = \alpha^2 \nabla^2 \phi + \phi ( 1 - \phi ) m_1 ( \phi , T) - 2 s \phi | \nabla \theta | - \epsilon^2 \phi | \nabla \theta |^2$

where

$m_1(\phi, T) = \phi - \frac{1}{2} - T \phi ( 1 - \phi )$

and the orientation equation is given by

$P(\epsilon | \nabla \theta |) \tau_{\theta} \phi^2 \frac{\partial \theta}{\partial t} = \nabla \cdot \left[ \phi^2 \left( \frac{s}{| \nabla \theta |} + \epsilon^2 \right) \nabla \theta \right]$

where

$P(w) = 1 - \exp{(-\beta w)} + \frac{\mu}{\epsilon} \exp{(-\beta w)}$

The initial conditions for this problem are set such that $\phi = 1$ for $0 \le x \le L_x$ and

$\theta = \begin{cases} 1 & \text{for $0 \le x < L_x / 2$,} \\ 0 & \text{for $L_x / 2 \le x \le L_x$.} \end{cases}$

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

examples.phase.impingement.test module¶

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