examples.elphf package¶

Subpackages¶

examples.elphf.diffusion package

Submodules¶

examples.elphf.input module¶

This example adds two more components to examples/elphf/input1DphaseBinary.py one of which is another substitutional species and the other represents electrons and diffuses interstitially.

Parameters from 2004/January/21/elphf0214

We start by defining a 1D mesh

>>> from fipy import PhysicalField as PF

>>> RT = (PF("1 Nav*kB") * PF("298 K"))
>>> molarVolume = PF("1.80000006366754e-05 m**3/mol")
>>> Faraday = PF("1 Nav*e")

>>> L = PF("3 nm")
>>> nx = 1200
>>> dx = L / nx
>>> # nx = 200
>>> # dx = PF("0.01 nm")
>>> ## dx = PF("0.001 nm") * (1.001 - 1/cosh(arange(-10, 10, .01)))
>>> # L = nx * dx
>>> mesh = Grid1D(dx = dx, nx = nx)
>>> # mesh = Grid1D(dx = dx)
>>> # L = mesh.facesRight[0].center[0] - mesh.facesLeft[0].center[0]
>>> # L = mesh.cellCenters[0,-1] - mesh.cellCenters[0,0]

We create the phase field

>>> timeStep = PF("1e-12 s")

>>> phase = CellVariable(mesh = mesh, name = 'xi', value = 1, hasOld = 1)
>>> phase.mobility = PF("1 m**3/J/s") / (molarVolume / (RT * timeStep))
>>> phase.gradientEnergy = PF("3.6e-11 J/m") / (mesh.scale**2 * RT / molarVolume)

>>> def p(xi):
...     return xi**3 * (6 * xi**2 - 15 * xi + 10.)

>>> def g(xi):
...     return (xi * (1 - xi))**2

>>> def pPrime(xi):
...     return 30. * (xi * (1 - xi))**2

>>> def gPrime(xi):
...     return 4 * xi * (1 - xi) * (0.5 - xi)

We create four components

>>> class ComponentVariable(CellVariable):
...     def __init__(self, mesh, value = 0., name = '', standardPotential = 0., barrier = 0., diffusivity = None, valence = 0, equation = None, hasOld = 1):
...         self.standardPotential = standardPotential
...         self.barrier = barrier
...         self.diffusivity = diffusivity
...         self.valence = valence
...         self.equation = equation
...         CellVariable.__init__(self, mesh = mesh, value = value, name = name, hasOld = hasOld)
...
...     def copy(self):
...         return self.__class__(mesh = self.mesh, value = self.value,
...                               name = self.name,
...                               standardPotential = self.standardPotential,
...                               barrier = self.barrier,
...                               diffusivity = self.diffusivity,
...                               valence = self.valence,
...                               equation = self.equation,
...                               hasOld = 0)

the solvent

>>> solvent = ComponentVariable(mesh = mesh, name = 'H2O', value = 1.)
>>> CnStandardPotential = PF("34139.7265625 J/mol") / RT
>>> CnBarrier = PF("3.6e5 J/mol") / RT
>>> CnValence = 0

and two solute species

>>> substitutionals = [
...     ComponentVariable(mesh = mesh, name = 'SO4',
...                       diffusivity = PF("1e-9 m**2/s") / (mesh.scale**2/timeStep),
...                       standardPotential = PF("24276.6640625 J/mol") / RT,
...                       barrier = CnBarrier,
...                       valence = -2,
...                       value = PF("0.000010414586295976 mol/l") * molarVolume),
...     ComponentVariable(mesh = mesh, name = 'Cu',
...                       diffusivity = PF("1e-9 m**2/s") / (mesh.scale**2/timeStep),
...                       standardPotential = PF("-7231.81396484375 J/mol") / RT,
...                       barrier = CnBarrier,
...                       valence = +2,
...                       value = PF("55.5553718417909 mol/l") * molarVolume)]

and one interstitial

>>> interstitials = [
...     ComponentVariable(mesh = mesh, name = 'e-',
...                       diffusivity = PF("1e-9 m**2/s") / (mesh.scale**2/timeStep),
...                       standardPotential = PF("-33225.9453125 J/mol") / RT,
...                       barrier = 0.,
...                       valence = -1,
...                       value = PF("111.110723815414 mol/l") * molarVolume)]

>>> for component in substitutionals:
...     solvent -= component
...     component.standardPotential -= CnStandardPotential
...     component.barrier -= CnBarrier
...     component.valence -= CnValence

Finally, we create the electrostatic potential field

>>> potential = CellVariable(mesh = mesh, name = 'phi', value = 0.)

>>> permittivity = PF("78.49 eps0") / (Faraday**2 * mesh.scale**2 / (RT * molarVolume))

>>> permittivity = 1.
>>> permittivityPrime = 0.

The thermodynamic parameters are chosen to give a solid phase rich in electrons and the solvent and a liquid phase rich in the two substitutional species

>>> solvent.standardPotential = CnStandardPotential
>>> solvent.barrier = CnBarrier
>>> solvent.valence = CnValence

Once again, we start with a sharp phase boundary

>>> x = mesh.cellCenters[0]
>>> phase.setValue(x < L / 2)
>>> interstitials[0].setValue("0.000111111503177394 mol/l" * molarVolume, where=x > L / 2)
>>> substitutionals[0].setValue("0.249944439430068 mol/l" * molarVolume, where=x > L / 2)
>>> substitutionals[1].setValue("0.249999982581341 mol/l" * molarVolume, where=x > L / 2)

We again create the phase equation as in examples.elphf.phase.input1D

>>> mesh.setScale(1)

>>> phase.equation = TransientTerm(coeff = 1/phase.mobility) \
...     == DiffusionTerm(coeff = phase.gradientEnergy) \
...     - (permittivityPrime / 2.) * potential.grad.dot(potential.grad)

We linearize the source term in the same way as in example.phase.simple.input1D.

>>> enthalpy = solvent.standardPotential
>>> barrier = solvent.barrier
>>> for component in substitutionals + interstitials:
...     enthalpy += component * component.standardPotential
...     barrier += component * component.barrier

>>> mXi = -(30 * phase * (1 - phase) * enthalpy +  4 * (0.5 - phase) * barrier)
>>> dmXidXi = (-60 * (0.5 - phase) * enthalpy + 4 * barrier)
>>> S1 = dmXidXi * phase * (1 - phase) + mXi * (1 - 2 * phase)
>>> S0 = mXi * phase * (1 - phase) - phase * S1

>>> phase.equation -= S0 + ImplicitSourceTerm(coeff = S1)

and we create the diffusion equation for the solute as in examples.elphf.diffusion.input1D

>>> for Cj in substitutionals:
...     CkSum = ComponentVariable(mesh = mesh, value = 0.)
...     CkFaceSum = FaceVariable(mesh = mesh, value = 0.)
...     for Ck in [Ck for Ck in substitutionals if Ck is not Cj]:
...         CkSum += Ck
...         CkFaceSum += Ck.harmonicFaceValue
...
...     counterDiffusion = CkSum.faceGrad
...     # phaseTransformation = (pPrime(phase.harmonicFaceValue) * Cj.standardPotential
...     #         + gPrime(phase.harmonicFaceValue) * Cj.barrier) * phase.faceGrad
...     phaseTransformation = (pPrime(phase).harmonicFaceValue * Cj.standardPotential
...             + gPrime(phase).harmonicFaceValue * Cj.barrier) * phase.faceGrad
...     # phaseTransformation = (p(phase).faceGrad * Cj.standardPotential
...     #         + g(phase).faceGrad * Cj.barrier)
...     electromigration = Cj.valence * potential.faceGrad
...     convectionCoeff = counterDiffusion + \
...         solvent.harmonicFaceValue * (phaseTransformation + electromigration)
...     convectionCoeff *= (Cj.diffusivity / (1. - CkFaceSum))
...
...     Cj.equation = (TransientTerm()
...                    == DiffusionTerm(coeff=Cj.diffusivity)
...                    + PowerLawConvectionTerm(coeff=convectionCoeff))

>>> for Cj in interstitials:
...     # phaseTransformation = (pPrime(phase.harmonicFaceValue) * Cj.standardPotential
...     #         + gPrime(phase.harmonicFaceValue) * Cj.barrier) * phase.faceGrad
...     phaseTransformation = (pPrime(phase).harmonicFaceValue * Cj.standardPotential
...             + gPrime(phase).harmonicFaceValue * Cj.barrier) * phase.faceGrad
...     # phaseTransformation = (p(phase).faceGrad * Cj.standardPotential
...     #         + g(phase).faceGrad * Cj.barrier)
...     electromigration = Cj.valence * potential.faceGrad
...     convectionCoeff = Cj.diffusivity * (1 + Cj.harmonicFaceValue) * \
...         (phaseTransformation + electromigration)
...
...     Cj.equation = (TransientTerm()
...                    == DiffusionTerm(coeff=Cj.diffusivity)
...                    + PowerLawConvectionTerm(coeff=convectionCoeff))

And Poisson’s equation

>>> charge = 0.
>>> for Cj in interstitials + substitutionals:
...     charge += Cj * Cj.valence

>>> potential.equation = DiffusionTerm(coeff = permittivity) + charge == 0

If running interactively, we create viewers to display the results

>>> from fipy import input
>>> if __name__ == '__main__':
...     phaseViewer = Viewer(vars=phase, datamin=0, datamax=1)
...     concViewer = Viewer(vars=[solvent] + substitutionals + interstitials, ylog=True)
...     potentialViewer = Viewer(vars = potential)
...     phaseViewer.plot()
...     concViewer.plot()
...     potentialViewer.plot()
...     input("Press a key to continue")

Again, this problem does not have an analytical solution, so after iterating to equilibrium

>>> solver = LinearLUSolver(tolerance = 1e-3)

>>> potential.constrain(0., mesh.facesLeft)

>>> phase.residual = CellVariable(mesh = mesh)
>>> potential.residual = CellVariable(mesh = mesh)
>>> for Cj in substitutionals + interstitials:
...     Cj.residual = CellVariable(mesh = mesh)
>>> residualViewer = Viewer(vars = [phase.residual, potential.residual] + [Cj.residual for Cj in substitutionals + interstitials])

>>> tsv = TSVViewer(vars = [phase, potential] + substitutionals + interstitials)

>>> dt = substitutionals[0].diffusivity * 100
>>> # dt = 1.
>>> elapsed = 0.
>>> maxError = 1e-1
>>> SAFETY = 0.9
>>> ERRCON = 1.89e-4
>>> desiredTimestep = 1.
>>> thisTimeStep = 0.
>>> print("%3s: %20s | %20s | %20s | %20s" % ("i", "elapsed", "this", "next dt", "residual"))
>>> residual = 0.
>>> from builtins import range
>>> from builtins import str
>>> for i in range(500): # iterate
...     if thisTimeStep == 0.:
...         tsv.plot(filename = "%s.tsv" % str(elapsed * timeStep))
...
...     for field in [phase, potential] + substitutionals + interstitials:
...         field.updateOld()
...
...     while True:
...         for j in range(10): # sweep
...             print(i, j, dt * timeStep, residual)
...             # raw_input()
...             residual = 0.
...
...             phase.equation.solve(var = phase, dt = dt)
...             # print phase.name, phase.equation.residual.max()
...             residual = max(phase.equation.residual.max(), residual)
...             phase.residual[:] = phase.equation.residual
...
...             potential.equation.solve(var = potential, dt = dt)
...             # print potential.name, potential.equation.residual.max()
...             residual = max(potential.equation.residual.max(), residual)
...             potential.residual[:] = potential.equation.residual
...
...             for Cj in substitutionals + interstitials:
...                 Cj.equation.solve(var = Cj,
...                                   dt = dt,
...                                   solver = solver)
...                 # print Cj.name, Cj.equation.residual.max()
...                 residual = max(Cj.equation.residual.max(), residual)
...                 Cj.residual[:] = Cj.equation.residual
...
...             # print
...             # phaseViewer.plot()
...             # concViewer.plot()
...             # potentialViewer.plot()
...             # residualViewer.plot()
...
...         residual /= maxError
...         if residual <= 1.:
...             break       # step succeeded
...
...         dt = max(SAFETY * dt * residual**-0.2, 0.1 * dt)
...         if thisTimeStep + dt == thisTimeStep:
...             raise FloatingPointError("step size underflow")
...
...     thisTimeStep += dt
...
...     if residual > ERRCON:
...         dt *= SAFETY * residual**-0.2
...     else:
...         dt *= 5.
...
...     # dt *= (maxError / residual)**0.5
...
...     if thisTimeStep >= desiredTimestep:
...         elapsed += thisTimeStep
...         thisTimeStep = 0.
...     else:
...         dt = min(dt, desiredTimestep - thisTimeStep)
...
...     if __name__ == '__main__':
...         phaseViewer.plot()
...         concViewer.plot()
...         potentialViewer.plot()
...         print("%3d: %20s | %20s | %20s | %g" % (i, str(elapsed * timeStep), str(thisTimeStep * timeStep), str(dt * timeStep), residual))

we confirm that the far-field phases have remained separated

>>> ends = take(phase, (0, -1))
>>> allclose(ends, (1.0, 0.0), rtol = 1e-5, atol = 1e-5)
1

and that the concentration fields has appropriately segregated into into their respective phases

>>> ends = take(interstitials[0], (0, -1))
>>> allclose(ends, (0.4, 0.3), rtol = 3e-3, atol = 3e-3)
1
>>> ends = take(substitutionals[0], (0, -1))
>>> allclose(ends, (0.3, 0.4), rtol = 3e-3, atol = 3e-3)
1
>>> ends = take(substitutionals[1], (0, -1))
>>> allclose(ends, (0.1, 0.2), rtol = 3e-3, atol = 3e-3)
1

examples.elphf.phase module¶

A simple 1D example to test the setup of the phase field equation.

We rearrange Eq. elphf:phase to

$\frac{1}{M_\xi}\frac{\partial \xi}{\partial t} &= \kappa_{\xi}\nabla^2 \xi + \frac{\epsilon'(\xi)}{2}\left(\nabla\phi\right)^2 \\ &\qquad - \left[ p'(\xi) \Delta\mu_n^\circ + g'(\xi) W_n \right] - \sum_{j=2}^{n-1} C_j \left[ p'(\xi) \Delta\mu_{jn}^\circ + g'(\xi) W_{jn} \right] - C_{\text{e}^{-}} \left[ p'(\xi) \Delta\mu_{\text{e}^{-}}^\circ + g'(\xi) W_{\text{e}^{-}} \right]$

The single-component phase field governing equation can be represented as

$\frac{1}{M_\xi} \frac{\partial \xi}{\partial t} = \kappa_\xi \nabla^2 \xi - 2\xi(1-\xi)(1-2\xi) W$

where $\xi$ is the phase field, $t$ is time, $M_\xi$ is the phase field mobility, $\kappa_\xi$ is the phase field gradient energy coefficient, and $W$ is the phase field barrier energy.

We solve the problem on a 1D mesh

>>> from fipy import CellVariable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, Viewer
>>> from fipy.tools import numerix

>>> nx = 400
>>> dx = 0.01
>>> L = nx * dx
>>> mesh = Grid1D(dx = dx, nx = nx)

We create the phase field

>>> phase = CellVariable(mesh = mesh, name = 'xi')
>>> phase.mobility = numerix.inf
>>> phase.gradientEnergy = 0.025

Although we are not interested in them for this problem, we create one field to represent the “solvent” component (1 everywhere)

>>> class ComponentVariable(CellVariable):
...     def copy(self):
...         new = self.__class__(mesh = self.mesh,
...                              name = self.name,
...                              value = self.value)
...         new.standardPotential = self.standardPotential
...         new.barrier = self.barrier
...         return new

>>> solvent = ComponentVariable(mesh = mesh, name = 'Cn', value = 1.)
>>> solvent.standardPotential = 0.
>>> solvent.barrier = 1.

and one field to represent the electrostatic potential (0 everywhere)

>>> potential = CellVariable(mesh = mesh, name = 'phi', value = 0.)
>>> permittivityPrime = 0.

We’ll have no substitutional species and no interstitial species in this first example

>>> substitutionals = []
>>> interstitials = []

>>> for component in substitutionals:
...     solvent -= component

>>> phase.equation = TransientTerm(coeff = 1/phase.mobility) \
...     == DiffusionTerm(coeff = phase.gradientEnergy) \
...     - (permittivityPrime / 2.) \
...        * potential.grad.dot(potential.grad)

>>> enthalpy = solvent.standardPotential
>>> barrier = solvent.barrier
>>> for component in substitutionals + interstitials:
...     enthalpy += component * component.standardPotential
...     barrier += component * component.barrier

We linearize the source term in the same way as in example.phase.simple.input1D.

>>> mXi = -(30 * phase * (1. - phase) * enthalpy \
...         +  4 * (0.5 - phase) * barrier)
>>> dmXidXi = (-60 * (0.5 - phase) * enthalpy + 4 * barrier)
>>> S1 = dmXidXi * phase * (1 - phase) + mXi * (1 - 2 * phase)
>>> S0 = mXi * phase * (1 - phase) - phase * S1

>>> phase.equation -= S0 + ImplicitSourceTerm(coeff = S1)

Note

Adding a Term to an equation formed with == will add to the left-hand side of the equation and subtracting a Term will add to the right-hand side of the equation

We separate the phase field into electrode and electrolyte regimes

>>> phase.setValue(1.)
>>> phase.setValue(0., where=mesh.cellCenters[0] > L / 2)

Even though we are solving the steady-state problem ( $M_\phi = \infty$ ) we still must sweep the solution several times to equilibrate

>>> from builtins import range
>>> for step in range(10):
...     phase.equation.solve(var = phase, dt=1.)

Since we have only a single component $n$ , with $\Delta\mu_n^\circ = 0$ , and the electrostatic potential is uniform, Eq. elphf:phase reduces to

$\frac{1}{M_\xi}\frac{\partial \xi}{\partial t} = \kappa_{\xi}\nabla^2 \xi - g'(\xi) W_n$

which we know from examples.phase.simple.input1D has the analytical solution

$\xi(x) = \frac{1}{2}(1 - \tanh\frac{x - L/2}{2d})$

with an interfacial thickness $d = \sqrt{\kappa_{\xi}/2W_n}$ .

We verify that the correct equilibrium solution is attained

>>> x = mesh.cellCenters[0]

>>> d = numerix.sqrt(phase.gradientEnergy / (2 * solvent.barrier))
>>> analyticalArray = (1. - numerix.tanh((x - L/2.)/(2 * d))) / 2.

>>> phase.allclose(analyticalArray, rtol = 1e-4, atol = 1e-4).value
1

If we are running interactively, we plot the error

>>> if __name__ == '__main__':
...     viewer = Viewer(vars = (phase - \
...         CellVariable(name = "analytical", mesh = mesh,
...                      value = analyticalArray),))
...     viewer.plot()

error in solution to steady-state phase field equation

examples.elphf.phaseDiffusion module¶

This example combines a phase field problem, as given in examples.elphf.phase.input1D, with a binary diffusion problem, such as described in the ternary example examples.elphf.diffusion.input1D, on a 1D mesh

>>> from fipy import CellVariable, FaceVariable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, PowerLawConvectionTerm, Viewer
>>> from fipy.tools import numerix

>>> nx = 400
>>> dx = 0.01
>>> L = nx * dx
>>> mesh = Grid1D(dx = dx, nx = nx)

We create the phase field

>>> phase = CellVariable(mesh=mesh, name='xi', value=1., hasOld=1)
>>> phase.mobility = 1.
>>> phase.gradientEnergy = 0.025

>>> def pPrime(xi):
...     return 30. * (xi * (1 - xi))**2

>>> def gPrime(xi):
...     return 4 * xi * (1 - xi) * (0.5 - xi)

and a dummy electrostatic potential field

>>> potential = CellVariable(mesh = mesh, name = 'phi', value = 0.)
>>> permittivityPrime = 0.

We start with a binary substitutional system

>>> class ComponentVariable(CellVariable):
...     def __init__(self, mesh, value = 0., name = '',
...                  standardPotential = 0., barrier = 0.,
...                  diffusivity = None, valence = 0, equation = None,
...                  hasOld = 1):
...         self.standardPotential = standardPotential
...         self.barrier = barrier
...         self.diffusivity = diffusivity
...         self.valence = valence
...         self.equation = equation
...         CellVariable.__init__(self, mesh = mesh, value = value,
...                               name = name, hasOld = hasOld)
...
...     def copy(self):
...         return self.__class__(mesh = self.mesh,
...                               value = self.value,
...                               name = self.name,
...                               standardPotential =
...                                   self.standardPotential,
...                               barrier = self.barrier,
...                               diffusivity = self.diffusivity,
...                               valence = self.valence,
...                               equation = self.equation,
...                               hasOld = 0)

consisting of the solvent

>>> solvent = ComponentVariable(mesh = mesh, name = 'Cn', value = 1.)

and the solute

>>> substitutionals = [
...     ComponentVariable(mesh = mesh, name = 'C1',
...                       diffusivity = 1., barrier = 0.,
...                       standardPotential = numerix.log(.3/.7) - numerix.log(.7/.3))]
>>> interstitials = []

>>> for component in substitutionals:
...     solvent -= component

The thermodynamic parameters are chosen to give a solid phase rich in the solute and a liquid phase rich in the solvent.

>>> solvent.standardPotential = numerix.log(.7/.3)
>>> solvent.barrier = 1.

We create the phase equation as in examples.elphf.phase.input1D and create the diffusion equations for the different species as in examples.elphf.diffusion.input1D

>>> def makeEquations(phase, substitutionals, interstitials):
...     phase.equation = TransientTerm(coeff = 1/phase.mobility) \
...         == DiffusionTerm(coeff = phase.gradientEnergy) \
...         - (permittivityPrime / 2.) \
...             * potential.grad.dot(potential.grad)
...     enthalpy = solvent.standardPotential
...     barrier = solvent.barrier
...     for component in substitutionals + interstitials:
...         enthalpy += component * component.standardPotential
...         barrier += component * component.barrier
...
...     mXi = -(30 * phase * (1 - phase) * enthalpy
...             +  4 * (0.5 - phase) * barrier)
...     dmXidXi = (-60 * (0.5 - phase) * enthalpy + 4 * barrier)
...     S1 = dmXidXi * phase * (1 - phase) + mXi * (1 - 2 * phase)
...     S0 = mXi * phase * (1 - phase) - phase * S1
...
...     phase.equation -= S0 + ImplicitSourceTerm(coeff = S1)
...
...     for Cj in substitutionals:
...         CkSum = ComponentVariable(mesh = mesh, value = 0.)
...         CkFaceSum = FaceVariable(mesh = mesh, value = 0.)
...         for Ck in [Ck for Ck in substitutionals if Ck is not Cj]:
...             CkSum += Ck
...             CkFaceSum += Ck.harmonicFaceValue
...
...         counterDiffusion = CkSum.faceGrad
...         phaseTransformation = (pPrime(phase.harmonicFaceValue) \
...                 * Cj.standardPotential
...                 + gPrime(phase.harmonicFaceValue) \
...                     * Cj.barrier) * phase.faceGrad
...         electromigration = Cj.valence * potential.faceGrad
...         convectionCoeff = counterDiffusion + \
...             solvent.harmonicFaceValue \
...                 * (phaseTransformation + electromigration)
...         convectionCoeff *= \
...             (Cj.diffusivity / (1. - CkFaceSum))
...
...         Cj.equation = (TransientTerm()
...                        == DiffusionTerm(coeff=Cj.diffusivity)
...                        + PowerLawConvectionTerm(coeff=convectionCoeff))
...
...     for Cj in interstitials:
...         phaseTransformation = (pPrime(phase.harmonicFaceValue) \
...             * Cj.standardPotential \
...             + gPrime(phase.harmonicFaceValue) \
...                 * Cj.barrier) * phase.faceGrad
...         electromigration = Cj.valence * potential.faceGrad
...         convectionCoeff = Cj.diffusivity \
...             * (1 + Cj.harmonicFaceValue) \
...             * (phaseTransformation + electromigration)
...
...         Cj.equation = (TransientTerm()
...                        == DiffusionTerm(coeff=Cj.diffusivity)
...                        + PowerLawConvectionTerm(coeff=convectionCoeff))

>>> makeEquations(phase, substitutionals, interstitials)

We start with a sharp phase boundary

$\xi = \begin{cases} 1& \text{for $x \le L/2$,} \\ 0& \text{for $x > L/2$,} \end{cases}$

or

>>> x = mesh.cellCenters[0]
>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L / 2)

and with a uniform concentration field $C_1 = 0.5$ or

>>> substitutionals[0].setValue(0.5)

If running interactively, we create viewers to display the results

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=([phase, solvent]
...                           + substitutionals + interstitials),
...                     datamin=0, datamax=1)
...     viewer.plot()

This problem does not have an analytical solution, so after iterating to equilibrium

>>> dt = 10000
>>> from builtins import range
>>> for i in range(5):
...     for field in [phase] + substitutionals + interstitials:
...         field.updateOld()
...     phase.equation.solve(var = phase, dt = dt)
...     for field in substitutionals + interstitials:
...         field.equation.solve(var = field,
...                              dt = dt)
...     if __name__ == '__main__':
...         viewer.plot()

phase and two concentration fields in equilibrium

we confirm that the far-field phases have remained separated

>>> numerix.allclose(phase(((0., L),)), (1.0, 0.0), rtol = 1e-5, atol = 1e-5)
1

and that the solute concentration field has appropriately segregated into solute-rich and solute-poor phases.

>>> print(numerix.allclose(substitutionals[0](((0., L),)), (0.7, 0.3), rtol = 2e-3, atol = 2e-3))
1

The same system of equations can model a quaternary substitutional system as easily as a binary. Because it depends on the number of substitutional solute species in question, we recreate the solvent

>>> solvent = ComponentVariable(mesh = mesh, name = 'Cn', value = 1.)

and make three new solute species

>>> substitutionals = [
...     ComponentVariable(mesh = mesh, name = 'C1',
...                       diffusivity = 1., barrier = 0.,
...                       standardPotential = numerix.log(.3/.4) - numerix.log(.1/.2)),
...     ComponentVariable(mesh = mesh, name = 'C2',
...                       diffusivity = 1., barrier = 0.,
...                       standardPotential = numerix.log(.4/.3) - numerix.log(.1/.2)),
...     ComponentVariable(mesh = mesh, name = 'C3',
...                       diffusivity = 1., barrier = 0.,
...                       standardPotential = numerix.log(.2/.1) - numerix.log(.1/.2))]

>>> for component in substitutionals:
...     solvent -= component
>>> solvent.standardPotential = numerix.log(.1/.2)
>>> solvent.barrier = 1.

These thermodynamic parameters are chosen to give a solid phase rich in the solvent and the first substitutional component and a liquid phase rich in the remaining two substitutional species.

Again, if we’re running interactively, we create a viewer

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=([phase, solvent]
...                           + substitutionals + interstitials),
...                     datamin=0, datamax=1)
...     viewer.plot()

We reinitialize the sharp phase boundary

>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L / 2)

and the uniform concentration fields, with the substitutional concentrations $C_1 = C_2 = 0.35$ and $C_3 = 0.15$ .

>>> substitutionals[0].setValue(0.35)
>>> substitutionals[1].setValue(0.35)
>>> substitutionals[2].setValue(0.15)

We make new equations

>>> makeEquations(phase, substitutionals, interstitials)

and again iterate to equilibrium

>>> dt = 10000
>>> from builtins import range
>>> for i in range(5):
...     for field in [phase] + substitutionals + interstitials:
...         field.updateOld()
...     phase.equation.solve(var = phase, dt = dt)
...     for field in substitutionals + interstitials:
...         field.equation.solve(var = field,
...                              dt = dt)
...     if __name__ == '__main__':
...         viewer.plot()

phase and four concentration fields in equilibrium

We confirm that the far-field phases have remained separated

>>> numerix.allclose(phase(((0., L),)), (1.0, 0.0), rtol = 1e-5, atol = 1e-5)
1

and that the concentration fields have appropriately segregated into their respective phases

>>> numerix.allclose(substitutionals[0](((0., L),)), (0.4, 0.3), rtol = 3e-3, atol = 3e-3)
1
>>> numerix.allclose(substitutionals[1](((0., L),)), (0.3, 0.4), rtol = 3e-3, atol = 3e-3)
1
>>> numerix.allclose(substitutionals[2](((0., L),)), (0.1, 0.2), rtol = 3e-3, atol = 3e-3)
1

Finally, we can represent a system that contains both substitutional and interstitial species. We recreate the solvent

>>> solvent = ComponentVariable(mesh = mesh, name = 'Cn', value = 1.)

and two new solute species

>>> substitutionals = [
...     ComponentVariable(mesh = mesh, name = 'C2',
...                       diffusivity = 1., barrier = 0.,
...                       standardPotential = numerix.log(.4/.3) - numerix.log(.4/.6)),
...     ComponentVariable(mesh = mesh, name = 'C3',
...                       diffusivity = 1., barrier = 0.,
...                       standardPotential = numerix.log(.2/.1) - numerix.log(.4/.6))]

and one interstitial

>>> interstitials = [
...     ComponentVariable(mesh = mesh, name = 'C1',
...                       diffusivity = 1., barrier = 0.,
...                       standardPotential = numerix.log(.3/.4) - numerix.log(1.3/1.4))]

>>> for component in substitutionals:
...     solvent -= component
>>> solvent.standardPotential = numerix.log(.4/.6) - numerix.log(1.3/1.4)
>>> solvent.barrier = 1.

These thermodynamic parameters are chosen to give a solid phase rich in interstitials and the solvent and a liquid phase rich in the two substitutional species.

Once again, if we’re running interactively, we create a viewer

>>> if __name__ == '__main__':
...     viewer = Viewer(vars=([phase, solvent]
...                           + substitutionals + interstitials),
...                     datamin=0, datamax=1)
...     viewer.plot()

We reinitialize the sharp phase boundary

>>> phase.setValue(1.)
>>> phase.setValue(0., where=x > L / 2)

and the uniform concentration fields, with the interstitial concentration $C_1 = 0.35$

>>> interstitials[0].setValue(0.35)

and the substitutional concentrations $C_2 = 0.35$ and $C_3 = 0.15$ .

>>> substitutionals[0].setValue(0.35)
>>> substitutionals[1].setValue(0.15)

We make new equations

>>> makeEquations(phase, substitutionals, interstitials)

and again iterate to equilibrium

>>> dt = 10000
>>> from builtins import range
>>> for i in range(5):
...     for field in [phase] + substitutionals + interstitials:
...         field.updateOld()
...     phase.equation.solve(var = phase, dt = dt)
...     for field in substitutionals + interstitials:
...         field.equation.solve(var = field,
...                              dt = dt)
...     if __name__ == '__main__':
...         viewer.plot()

phase and four concentration fields (one like electrons) in equilibrium

We once more confirm that the far-field phases have remained separated

>>> numerix.allclose(phase(((0., L),)), (1.0, 0.0), rtol = 1e-5, atol = 1e-5)
1

and that the concentration fields have appropriately segregated into their respective phases

>>> numerix.allclose(interstitials[0](((0., L),)), (0.4, 0.3), rtol = 3e-3, atol = 3e-3)
1
>>> numerix.allclose(substitutionals[0](((0., L),)), (0.3, 0.4), rtol = 3e-3, atol = 3e-3)
1
>>> numerix.allclose(substitutionals[1](((0., L),)), (0.1, 0.2), rtol = 3e-3, atol = 3e-3)
1

examples.elphf.poisson module¶

A simple 1D example to test the setup of the Poisson equation.

>>> from fipy import CellVariable, Grid1D, DiffusionTerm, Viewer
>>> from fipy.tools import numerix

>>> nx = 200
>>> dx = 0.01
>>> L = nx * dx
>>> mesh = Grid1D(dx = dx, nx = nx)

The dimensionless Poisson equation is

$\nabla\cdot\left(\epsilon\nabla\phi\right) = -\rho = -\sum_{j=1}^n z_j C_j$

where $\phi$ is the electrostatic potential, $\epsilon$ is the permittivity, $\rho$ is the charge density, $C_j$ is the concentration of the $j^\text{th}$ component, and $z_j$ is the valence of the $j^\text{th}$ component.

We will be solving for the electrostatic potential

>>> potential = CellVariable(mesh = mesh, name = 'phi', value = 0.)
>>> permittivity = 1.

We examine a fixed distribution of electrons with $z_{\text{e}^{-}} = -1$ .

>>> class ComponentVariable(CellVariable):
...     def __init__(self, mesh, value = 0., name = '',
...                  standardPotential = 0., barrier = 0.,
...                  diffusivity = None, valence = 0, equation = None):
...         CellVariable.__init__(self, mesh = mesh,
...                               value = value, name = name)
...         self.standardPotential = standardPotential
...         self.barrier = barrier
...         self.diffusivity = diffusivity
...         self.valence = valence
...         self.equation = equation
...
...     def copy(self):
...         return self.__class__(mesh = self.mesh,
...                               value = self.value,
...                               name = self.name,
...                               standardPotential =
...                                   self.standardPotential,
...                               barrier = self.barrier,
...                               diffusivity = self.diffusivity,
...                               valence = self.valence,
...                               equation = self.equation)

Since we’re only interested in a single species, electrons, we could simplify the following, but since we will in general be studying multiple components, we explicitly allow for multiple substitutional species and multiple interstitial species:

>>> interstitials = [
...     ComponentVariable(mesh = mesh, name = 'e-', valence = -1)]
>>> substitutionals = []

Because Poisson’s equation admits an infinite number of potential profiles, we must constrain the solution by fixing the potential at one point:

>>> potential.constrain(0., mesh.facesLeft)

>>> charge = 0.
>>> for Cj in interstitials + substitutionals:
...     charge += Cj * Cj.valence

>>> potential.equation = DiffusionTerm(coeff = permittivity) \
...                      + charge == 0

First, we obtain a uniform charge distribution by setting a uniform concentration of electrons $C_{\text{e}^{-}} = 1$ .

>>> interstitials[0].setValue(1.)

and we solve for the electrostatic potential

>>> potential.equation.solve(var = potential)

This problem has the analytical solution

$\psi(x) = \frac{x^2}{2} - 2x$

We verify that the correct equilibrium is attained

>>> x = mesh.cellCenters[0]
>>> analyticalArray = (x**2)/2 - 2*x

>>> print(potential.allclose(analyticalArray, rtol = 2e-5, atol = 2e-5))
1

If we are running the example interactively, we view the result

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer = Viewer(vars = (charge, potential))
...     viewer.plot()
...     input("Press any key to continue...")

Next, we segregate all of the electrons to right side of the domain

$C_{\text{e}^{-}} = \begin{cases} 0& \text{for $x \le L/2$,} \\ 1& \text{for $x > L/2$.} \end{cases}$

>>> x = mesh.cellCenters[0]
>>> interstitials[0].setValue(0.)
>>> interstitials[0].setValue(1., where=x > L / 2.)

and again solve for the electrostatic potential

>>> potential.equation.solve(var = potential)

which now has the analytical solution

$\psi(x) = \begin{cases} -x& \text{for $x \le L/2$,} \\ \frac{(x-1)^2}{2} - x& \text{for $x > L/2$.} \end{cases}$

We verify that the correct equilibrium is attained

>>> analyticalArray = numerix.where(x < L/2, -x, ((x-1)**2)/2 - x)

>>> potential.allclose(analyticalArray, rtol = 2e-5, atol = 2e-5).value
1

and again view the result

>>> from fipy import input
>>> if __name__ == '__main__':
...     viewer.plot()
...     input("Press any key to continue...")

Finally, we segregate all of the electrons to left side of the domain

$C_{\text{e}^{-}} = \begin{cases} 1& \text{for $x \le L/2$,} \\ 0& \text{for $x > L/2$.} \end{cases}$

>>> interstitials[0].setValue(1.)
>>> interstitials[0].setValue(0., where=x > L / 2.)

and again solve for the electrostatic potential

>>> potential.equation.solve(var = potential)

which has the analytical solution

$\psi(x) = \begin{cases} \frac{x^2}{2} - x& \text{for $x \le L/2$,} \\ -\frac{1}{2}& \text{for $x > L/2$.} \end{cases}$

We again verify that the correct equilibrium is attained

>>> analyticalArray = numerix.where(x < 1, (x**2)/2 - x, -0.5)

>>> potential.allclose(analyticalArray, rtol = 2e-5, atol = 2e-5).value
1

and again view the result

>>> if __name__ == '__main__':
...     viewer.plot()

examples.elphf.test module¶

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