View source for 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:

<math>
\dfrac{\partial X_i}{\partial t} = \nabla \cdot \left( D_i \nabla X_i \right)
</math>

subject to boundary conditions:

<math>
X_i ( x=0 ) = 1
X_i ( x=1 ) = 0
</math>

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

<math>
D_i \neq D_i(X_i)
</math>

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:

<math>
\sum_i X_i = 1
</math>

=The Script=

Begin with import statements:

<pre>
from fipy import *
from numpy import *
import Cantera
import matplotlib.pylab as plt
</pre>

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

<pre>
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title('Ar')
</pre>

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, <code>a</code> and <code>ari</code>, which are defined later in the script. These represent the gas phase object and the argon species index, respectively.)

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

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

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

Now we define Fipy equations/grid:

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

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

<math>
\dfrac{\partial X_i}{\partial t} = \nabla \cdot \left( D \nabla X_i \right) 
</math>

and the diffusion coefficient is a function of <math>X_i</math>. 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).

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

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.

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

=Result=

This will result in the following plot:

[[Image:OneVarCanteraVariableD.png]]

=References=

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

{{CanteraFlag}}
[[Category:Fipy]]