Fipy and Cantera/Diffusion 1D

Overview

This example illustrates how to solve a simple, time-dependent 1D diffusion problem using Fipy with diffusion coefficients computed by Cantera. This is a very simple problem. The point is not to demonstrate earth-shaking complexity, the point is illustrating how to make these two packages talk to each other.

The equation being solved is:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \dfrac{\partial X_i}{\partial t} = \nabla \cdot \left( D_i \nabla X_i \right) }

subject to boundary conditions:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X_i ( x=0 ) = 1 X_i ( x=1 ) = 0 }

The mixture being used is argon in oxygen. Although the diffusion coefficient is a very weak function of composition, and therefore the assumption

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle D_i \neq D_i(X_i) }

would normally apply, this example illustrates the sweep solution technique (in which the diffusion coefficient is a function of the composition) anyway.

Note that because this is a binary mixture, we only need to solve for a single component (argon), because:

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \sum_i X_i = 1 }

The Script

Begin with import statements:

from fipy import *
from numpy import *
import Cantera
import matplotlib.pylab as plt

Next, create a plot, which we'll use to visualize our results:

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title('Ar')

Computing the diffusion coefficient using Cantera proceeds as follows:

• Pass in a Fipy CellVariable object
• Create a CellVariable to contain diffusion coefficient values
• Looping over each cell:
  • Translate the value of the Fipy CellVariable at a given cell into a thermochemical state for a Cantera phase object (we want to use a single Cantera phase object, so we don't have to keep loading our CTI/XML mechanism file into memory)
  • Compute the value of the diffusion coefficient
  • Populate the diffusion CellVariable
• Return the diffusion CellVariable

We'll define this in a function, which we can then pass when specifying the diffusion coefficient for the Fipy equation. (Note that I'm going to use two global variables here, a and ari, which are defined later in the script. These represent the gas phase object and the argon species index, respectively.)

###################################### 
# Diffusion coefficient function
def get_diffusion_coeff( phi ):

    D = CellVariable(name="DiffusionCoeff",mesh=mesh)

    for i in range(phi.shape[0]):
        var = phi.value[i]
        X0 = "AR:%0.8f, O2:%0.8f"%(var,1-var)
        a.set(X=X0)

        # mixture-averaged diffusion coefficients
        Dcoeff = a.mixDiffCoeffs()[ari]

        D.put(i,Dcoeff)

    return D

We also want to define a single Cantera phase, which we'll use in the get_diffusion_coeff() function:

###################################### 
# Gas object 
#
# This object will be used to evaluate all thermodynamic, kinetic,
# and transport properties
#
rxnmech    =  'h2o2.cti'       # reaction mechanism file
mix        =  'ohmech'         # gas mixture model
comp       =  'O2:1.0, AR:1.0' # gas composition
a = Cantera.IdealGasMix(rxnmech, mix)
a.set(T = 298.15, P = Cantera.OneAtm, X = comp)

arlab = 'AR'
o2lab = 'O2'
ari = a.speciesIndex(arlab)
o2i = a.speciesIndex(o2lab)

Now we define Fipy equations/grid:

# First, define the mesh
nx = 100
dx = 0.01
L = nx * dx
xs = linspace(0,L,nx)
mesh = Grid1D(nx=nx,dx=dx)

# Define boundary values 
arLeft  = 1.
arRight = 0.

# Define timestep-related parameters
# Courant stability constraint for timestep
D0 = 1.0e-5
safetyFactor = 0.9
timestepDuration = safetyFactor * dx**2 / (2 * D0)

# Simulation duration
steps = 500
print "Nsteps =",steps
t = timestepDuration * steps


# Now define your sweep variables
arvar = [CellVariable(name=arlab, mesh=mesh, value=arRight, hasOld=1)]

At this point, we can define the Fipy equation that we're solving, which is, as given above,

Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle \dfrac{\partial X_i}{\partial t} = \nabla \cdot \left( D \nabla X_i \right) }

and the diffusion coefficient is a function of Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "https://en.wikipedia.org/api/rest_v1/":): {\displaystyle X_i} . This is subject to the boundary conditions set above (1 on the left side of the domain, 0 on the right side of the domain).

# Define the equation being solved
eq = TransientTerm() == DiffusionTerm(coeff=get_diffusion_coeff(arvar[0]))

# Set the boundary conditions
arvar[0].constrain(arLeft, mesh.facesLeft)
arvar[0].constrain(arRight, mesh.facesRight)

Now we timestep through the solutions, and plot the solution every 50 timesteps. At each timestep, we are performing 5 sweeps (that is, we're performing 5 iterations per timestep). Again, our diffusion coefficient is not a strong function of composition, because Ar and O2 are so similar in shape/properties, so even 5 sweeps is probably unnecessary. Again, this is purely for illustrative purposes.

# Timestep through solutions
arvar[0].setValue(arRight)
for step in range(steps):
    # only move forward in time once per time step
    arvar[0].updateOld()
    
    # but "sweep" many times per time step
    for sweep in range(5):
        res = eq.sweep(var=arvar[0],
                       dt=timestepDuration)

    if step%50==0:
        ax.plot(xs,arvar[0].value)
        ax.hold(True)
        plt.show()
        plt.draw()

Result

This will result in the following plot:

References

Fipy documentation covering 1D diffusion problems: http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html