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:

${\dfrac {\partial X_{i}}{\partial t}}=\nabla \cdot \left(D_{i}\nabla X_{i}\right)$

subject to boundary conditions:

$X_{i}(x=0)=1X_{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

$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:

$\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,

${\dfrac {\partial X_{i}}{\partial t}}=\nabla \cdot \left(D\nabla X_{i}\right)$

and the diffusion coefficient is a function of $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