Fipy and Cantera/Diffusion 1D

Overview

This example illustrates how to solve a simple, time-dependent 1D diffusion problem using Fipy with diffusion coefficeints computed by Cantera

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

"""
1D Diffusion:
Cantera Example with Variable Diffusivity

Charles Reid
December 2013

Solve a 1D diffusion problem in which 
two species diffuse into each other.

Diffusion coefficients vary with composition.

Species: 
    Ar
    O2 

BCs: 
    Ar(x=0) = 1
    Ar(x=1) = 0
    O2(x=0) = 0
    O2(x=1) = 1

See also:
    http://www.ctcms.nist.gov/fipy/examples/diffusion/generated/examples.diffusion.mesh1D.html
"""


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


###################################### 
# 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


###################################### 
# 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)

# 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)]

# 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)

# 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()