Cantera with internal combustion engine code

[Shah Saud Alam]

Hi,

I need some direction/help on a project that I have been working on. A single-cylinder internal combustion engine (compression-ignition) operates at 1800 rpm and utilizes methane as the only fuel. In order to solve a thermodynamic model of the engine, we need the following three things:

1. GRI3.0 reaction mechanism + Cantera to handle the reaction kinetics.

2. A slider-crank mechanism to predict the instantaneous piston position, velocity, in-cylinder volume, and the combustion surface area.

3. The thermodynamic equations to predict the pressure, temperature, gas molecular weight, gas constant, mixture density, etc.

I have written a MatLab code to handle the above, however, the "combustion" never happens. All that happens is that the gas mixture gets compressed without combustion ever happening.

Am I missing something here? Do I need nested ODE solvers here - the first one handles the engine motion ODEs whereas the second one handles the temperature, pressure, and other thermodynamic ODEs?

I can send you the code to your email separately, as it might violate the group policy if I post it here.

Thanks in advance for helping me out with this.

-ps- I have written codes for constant-volume and constant-pressure combustion kinetics, and they work flawlessly.

[Ingmar Schoegl]

Please have a look at the Python ic_engine.py example. For compression-ignition to work with natural gas, you need unrealistically high compression ratios (the example uses 50!). The updated example that will ship with the upcoming release 2.5 uses a more realistic fuel with compression ratios in line with expectations. Hope this helps,

-Ingmar-

[Shah Saud Alam]

Thanks a lot, Ingmar. I checked the ic_engine.py example that you mentioned. Do you know why Cantera requires unrealistically high compression ratios for compression-ignition simulations? Is this a Cantera-specific or reaction mechanism-specific problem?

Also, do you know when is Cantera 2.5 coming out?

Thank you.

[bryan.a....@gmail.com]

I modeled HCCI combustion of natural gas at 1800 RPM with Cantera and GRI 3.0 just fine, with compression ratios between 16 and 20.  I did not have any difficulty accurately predicting ignition timing.  Combustion duration and efficiency, of course, require a relatively complex model, which I was not running.  For the 16-20 compression ratios I was looking at, elevated pressure and temperature at BDC are necessary.  For example, my HCCI model was pretty happy at 50% residuals, 1070 K, 165 bar at TDC, with ignition following shortly thereafter.  This requires fairly aggressive boost and benefits from internal EGR, limited charge cooling, etc.  I would guess that the aforementioned "unrealistic" compression ratios are without boost and/or elevated intake temperature.

Some affects might have to do with your system setup.  For example, direct injection later in the stroke will require more aggressive TDC conditions for ignition, since there is less time for the charge to "cook".

If your constant volume code shows combustion occurring at conditions you should be achieving around TDC, obvious things to check are:

1) Are you really achieving that condition at TDC?  In my work, heat transfer was significantly overpredicted using literature values; I had to calibrate it downward.  Easy enough to export and examine both pressure and temperature traces for the cycle.

2) Ignition delay.  Obviously, if charge reexpands too much before combustion really starts moving, it will tend to fizzle.

Hope that helps.

On Sunday, August 30, 2020 at 4:34:17 PM UTC-5 saud...@gmail.com wrote:

[Ingmar Schoegl]

Dear Shah,

In a nutshell, your conditions at TDC need to allow for chemistry to be sufficiently fast for ignition to happen, which is mainly a function of temperature (methane ignition is somewhat sluggish, I.e. it needs higher temperatures, which high compression ratios accomplish). Boosting/preheating at BDC like Bryan mentioned will do the same for you. You probably could reduce the compression ratio  in the example shipped with 2.4 (despite the original comment) somewhat which I however haven’t tested; when I updated the existing example for 2.5, my main motivation was to make it run more efficiently. 

At the moment 2.5 is in late alpha and quite stable (simply download the conda cantera-devel version). It shouldn’t take too long before the final release at this point.

-Ingmar-

--


You received this message because you are subscribed to a topic in the Google Groups "Cantera Users' Group" group.


To unsubscribe from this topic, visit https://groups.google.com/d/topic/cantera-users/AXpFPywlE60/unsubscribe.


To unsubscribe from this group and all its topics, send an email to cantera-user...@googlegroups.com.


To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/f79a7477-1d0b-4587-bd15-b64513b3c9e6o%40googlegroups.com.

[Shah Saud Alam]

Hello Bryan,

Sorry for not being able to keep up with the thread. I got busy with the schools and universities reopening, and this issue got pushed to the last. I did try to check the code and the physics (as well as chemistry) involved here, and it seems that I am correct with the methodology. It appears to me that this is either a reaction mechanism issue or something to do with Cantera. I have modeled this in Matlab, however, I have also translated this model into Python, and I will attach the code at the end of the next reply.

-Saud

[Shah Saud Alam]

Hello Ingmar,

Thanks for the input. I also implemented your suggestions and tested the model for increased temperatures. You are correct when you say that the GRI mechanism is sluggish and responds to higher temperatures better.

My advisor used GRI with a Fortran code and Chemkin, and was able to achieve combustion for the given conditions at lower temperatures. Similar conditions with Cantera just give a motoring curve for the pressure and combustion does not happen at all.

To unsubscribe from this group and all its topics, send an email to canter...@googlegroups.com.

[Shah Saud Alam]

Here are the operating conditions:

T0 = 387 K, P0 = 1.47e5 Pa, X = 'CH4:0.0476,H2:0.0290,O2:0.1096,N2:0.8138'

Here is the original code in Matlab:

----------------------------------------------------------------------------------------------------------------------------------

clear; close all; clc;

% set up the variables array

N = 3600;

b = 86/1000;

s = 75/1000;

l = 118/1000;

a = s/2;

Spm = 4.5;

Twall = 400;

ep = 0.65;

sig = 5.67e-8;

var = [N,a,l,b,Spm,Twall,ep,sig];

T = 387;

p = 1.47e5;

% set up the gas object with Cantera

gas = IdealGasMix('gri30.xml');

% set(gas,'T',T,'P',p,'X','H2:1,O2:0.5,N2:1.88')

set(gas,'T',T,'P',p,'X','CH4:0.0476,H2:0.0290,O2:0.1096,N2:0.8138');

% set up the initial conditions array

th  = 4.154;

x   = 0.0938;

As  = 2.827e-2;

V   = 3.802e-4;

rho = density(gas);

Vel = 0;

u   = intEnergy_mass(gas);

T   = temperature(gas);

Mmix = meanMolecularWeight(gas);

R   = gasconstant/Mmix;

p   = pressure(gas);

Y0  = massFractions(gas);

X0  = moleFractions(gas);

IC = [th,x,As,V,rho,Vel,u,T,Mmix,R,p,Y0',X0'];

% dt = 1.85185e-5;

% tend = 0.02463;

dt = 1.1574e-5;

tend = 0.012315;

options = odeset('RelTol',1e-5,...      % relative tolerance

                    'AbsTol',1e-8,...  % absolute tolerance

                    'MaxStep',dt,...    % maximum time step

                    'Stats','On',...    % displays the stats from the iteration

                    'Refine',1);        % produces smoother output

[t,Z] = ode15s(@rxnode,[0 tend],IC,options,gas,var);

for i = 1:11

    figure

    plot(t,Z(:,i),'linewidth',2)

    switch i

        case 1

            str = 'Crank angle [rad]';

        case 2

            str = 'Piston position [m]';

        case 3

            str = 'Surface area [m^2]';

        case 4

            str = 'Volume [m^3]';

        case 5

            str = 'Density [kg/m^3]';

        case 6

            str = 'Velocity [m/s]';

        case 7

            str = 'Internal energy [J/kg]';

        case 8

            str = 'Temperature [K]';

        case 9

            str = 'Mixture weight [gm/mol]';

        case 10

            str = 'Gas constant [J/kg.K]';

        case 11

            str = 'Pressure [bar]';

            plot(t,Z(:,i)/1e5,'linewidth',2)

    end

    ylabel(str)

    xlabel('Time [s]')

    box on, grid on, grid minor

end

figure

hold on

for i = [1,4,6,14]

    plot(t,Z(:,11+i),'linewidth',2)

end

hold off

xlabel('Time [s]')

ylabel('Species mass fractions [-]')

legend('H_2','O_2','H_2O','CH_4','location','best')

box on, grid on, grid minor

----------------------------------------------------------------------------------------------------------------------------------

Here is the function file that handles the ODEs:

function out = rxnode(t,y,gas,vars)

NM = nSpecies(gas);         % number of species in the gas mixture

Wk = molecularWeights(gas); % species molecular weights [gm/mol or kg/kmol]

% unpack the initial conditions array

q   = y(1);     % crank angles [rad]

x   = y(2);     % piston velocity [m]

As  = y(3);     % surface area [m^2]

V   = y(4);     % volume [m^3]

rho = y(5);     % gas mixture density [kg/m^3]

Vel = y(6);     % gas flow velocity [m/s]

u   = y(7);     % internal energy [J/kg]

T   = y(8);     % gas temperature [K]

Mmix = y(9);    % gas mixture mean molecular weight [gm/mol or kg/kmol]

R   = y(10);    % specific gas constant [J/kg.K]

p   = y(11);    % gas pressure [Pa]

Y0  = y(12:12+NM-1); % species mass fractions [-]

X0  = y(12+NM:end);  % species mole fractions [-]

% unpack the variable array

N = vars(1);    % engine speed [rpm]

a = vars(2);    % crank pin radius [m]

l = vars(3);    % connecting rod length [m]

b = vars(4);    % bore [m]

Spm     = vars(5);  % mean piston speed [m/s]

Twall   = var(6);   % engine wall temperature [K]

ep  = var(7);   % emissivity [-]

sig = var(8);   % Stefan-Boltzmann constant [W/m^2.K^4]

set(gas,'T',T,'P',p,'Y',Y0);

% final crank angle

dqdt = 2*pi*N/60;

% final piston position

dxdt = (-a*sin(q)- a^2*sqrt(2)*cos(q)*sin(q)/sqrt(2*l^2-2*a^2*sin(q)^2))*dqdt;

% final surface area

dAsdt = -pi*b*dxdt;

% final volume

dVdt = -pi*b^2/4*dxdt;

% final conservation of mass

drdt = -rho/V*dVdt;

% final conservation of momentum

dVeldt = 0;

% final conservation of species

wdot = netProdRates(gas);   % net production or destruction rates [kmol/m^3.s]

dYdt = wdot.*Wk/rho;

% final conservation of energy

hc = 3.26*(V^-0.06)*((p/1000)^0.8)*(T^-0.4)*(Spm+1.40)^0.8; % Hohenberg convective heat transfer coefficient [W/m^2.K]

dudt = -p/V/rho*dVdt - hc*As/V/rho*(T-Twall) - ep*sig*As/rho/V*(T^4-Twall^4);

% dudt = 0;

% % final temperature

% dTdt = dudt/cv_mass(gas);

dTdt = (-temperature(gas)*gasconstant*(enthalpies_RT(gas)-ones(NM,1))'*wdot/density(gas) ...

    -p/density(gas)/V*dVdt - hc*As/density(gas)/V*(T-Twall) - ep*sig*As/density(gas)/V*(T^4-Twall^4))/cv_mass(gas);

disp(dTdt)

% final mole fractions

dXdt = zeros(NM,1);

for j = 1:NM

    sum1 = 0;

    sum2 = 0;

    for i = 1:NM

        sum1 = sum1 + Y0(i)/Wk(i);

        sum2 = sum2 + dYdt(i)/Wk(i);

    end

    dXdt(j) = (1/Wk(j)*dYdt(j)*sum1 - Y0(j)/Wk(j)*sum2)/sum1^2;

end

% final molecular weight

dMmixdt = 0;

for i = 1:NM

    dMmixdt = dMmixdt + Wk(i)*dXdt(i);

end

% final gas constant

dRdt = -gasconstant/Mmix^2*dMmixdt;

% final ideal gas law

dpdt = rho*R*dTdt + R*T*drdt + rho*T*dRdt;

out = [dqdt;dxdt;dAsdt;dVdt;drdt;dVeldt;dudt;dTdt;dMmixdt;dRdt;dpdt;dYdt;dXdt];

----------------------------------------------------------------------------------------------------------------------------------

My Python code looks like this:

----------------------------------------------------------------------------------------------------------------------------------

import cantera as ct

import matplotlib.pyplot as plt

from scipy.integrate import odeint

import numpy as np

from math import *

# define the function for the ODEs

def rxnode(y,t,gas,var):

    # extract important information from the gas object

    Wk = gas.molecular_weights

    NM = gas.n_total_species

    N = var[0] # engine speed [rpm]

    a = var[1] # crank pin radius [m]

    l = var[2] # connecting rod length [m]

    b = var[3] # bore [m]

    Spm = var[4] # mean piston speed [m/s]

    Twall = var[5] # engine wall temperature [K]

    ep = var[6]    # emissivity [-]

    sig = var[7]   # Stefan-Boltzmann constant [W/m^2.K^4]

    

#    q,x,As,V,rho,Vel,u,T,Mmix,R,p,Y0,X0 = y

    q = y[0]    # crank angles [rad]

    x = y[1]    # piston velocity [m/s]

    As = y[2]   # surface area [m^2]

    V = y[3]    # volume [m^3]

    rho = y[4]  # gas mixture density [kg/m^3]

    Vel = y[5]  # gas flow velocity [m/s]

    u = y[6]    # gas (specific) internal energy [J/kg]

    T = y[7]    # gas temperature [K]

    Mmix = y[8] # gas mixture mean molecular weight [gm/mol or kg/kmol]

    R = y[9]    # specific gas constant [J/kg.K]

    p = y[10]   # gas pressure [Pa]

    Y0 = y[11:11+NM]  # species mass fractions [-]

    X0 = y[11+NM:]  # species mole fractions [-]

    # set gas properties

    gas.TPX = T, p, X0

    ## Start writing ODEs below:

    # final crank angle

    dqdt = 2*pi*N/60

    # final piston position

    dxdt = (-a*sin(q)- a**2*sqrt(2)*cos(q)*sin(q)/sqrt(2*l**2-2*a**2*sin(q)**2))*dqdt

    # final surface area

    dAsdt = -pi*b*dxdt

    # final volume

    dVdt = -pi*b**2/4*dxdt

    # final conservation of mass

    drdt = -rho/V*dVdt

    # final conservation of momentum

    dVeldt = 0

    # final conservation of species

    wdot = gas.net_production_rates

    dYdt = wdot*Wk/rho

    # final conservation of energy

    # internal energy basis

    hc = 3.26*(V**(-0.06))*((p/1000)**0.8)*(T**(-0.4))*(Spm+1.40)**(0.8)

    dQdt = -p/V/rho*dVdt - hc*As/V/rho*(T-Twall) - ep*sig*As/V/rho*(T**4-Twall**4)

    dudt = -p/V/rho*dVdt - hc*As/V/rho*(T-Twall) - ep*sig*As/V/rho*(T**4-Twall**4)

    # temperature basis

    h = gas.standard_enthalpies_RT

    cv = gas.cv_mass

    dTdt = 0

    for i in range(NM):

        dTdt = dTdt - T*ct.gas_constant*(h[i]-1)*wdot[i]/rho

    dTdt = (dTdt + dQdt)/cv

    # final mole fractions

    dXdt = np.zeros((NM))

    for j in range(NM):

        sum1 = np.dot(Y0,1/Wk)

        sum2 = np.dot(dYdt,1/Wk)

        dXdt[j] = (1/Wk[j]*dYdt[j]*sum1 - Y0[j]/Wk[j]*sum2)/sum1**2

    # final molecular weight

    dMmixdt = np.dot(Wk,dXdt)

    # final gas constant

    Mmix = gas.mean_molecular_weight

    dRdt = -ct.gas_constant/Mmix**2*dMmixdt

    # final Ideal Gas Law

    dpdt = rho*R*dTdt + R*T*drdt + rho*T*dRdt

    return np.hstack((dqdt,dxdt,dAsdt,dVdt,drdt,dVeldt,dudt,dTdt,dMmixdt,dRdt,dpdt,dYdt,dXdt))

# write the driver code below

N = 1800 # engine speed [rpm]

b = 86/1000 # engine bore [m]

s = 75/1000 # engine stroke [m]

l = 118/1000 # engine connecting rod length [m]

a = s/2 # engine crank pin radius [m]

Spm = 4.5 # mean piston speed [m/s]

Twall = 400 # engine wall temperature [K]

ep = 0.65 # emissivity of the engine [-]

sig = 5.67e-8 # Stefan-Boltzmann constant [W/m^2.K^4]

var = N, a, l, b, Spm, Twall, ep, sig   # package the var array

T = 600 # initial temperature [K]

p = 1e5 # initial pressure [Pa]

gas = ct.Solution('gri30.xml')

gas.TPX = T, p, 'CH4:0.0476,H2:0.0290,O2:0.1096,N2:0.8138'

# setup the initial conditions array

th = 4.154 # initial crank angle [rad]

x = 0.0938 # initial piston position [m]

As = 2.827e-2 # initial combustion surface area [m^2]

V = 3.802e-4 # initial volume [m^3]

rho = gas.density # initial density [kg/m^3]

Vel = 0 # initial gas flow velocity [m/s]

u = gas.int_energy_mass # initial internal energy [J/kg]

T = gas.T # initial gas temperature [K]

Mmix = gas.mean_molecular_weight # initial gas mean molecular weight [gm/mol or kg/kmol]

R = ct.gas_constant/Mmix # initial gas constant [J/kg.K]

p = gas.P # initial pressure [Pa]

Y0 = gas.Y # initial mass fractions [-]

X0 = gas.X # initial mole fractions [-]

y0 = np.hstack((th, x, As, V, rho, Vel, u, T, Mmix, R, p, Y0, X0))

dt = 1.85185e-5 # time step, corresponds to 1 reading per 0.2 deg of crank angle at 1800 RPM

tend = 0.0243 # end time for simulation [s]

t = np.arange(0,tend,dt)

y = odeint(rxnode,y0,t,args=(gas,var,),full_output=1,rtol=1e-5,atol=1e-8,hmax=dt)

----------------------------------------------------------------------------------------------------------------------------------

I get the following output during the simulation: "ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. warnings.warn(warning_msg, ODEintWarning)"

Since I am new to Python, I am not sure what this means. Perhaps, it's related to the scaling of the step size or something, but I would appreciate your inputs in this regard.

Thank you.

On Monday, 31 August 2020 17:47:33 UTC-5, Ingmar Schoegl wrote:

To unsubscribe from this group and all its topics, send an email to canter...@googlegroups.com.

[bryan.a....@gmail.com]

I'm not sure I follow everything in here, but it seems that you're trying to build a lot that Cantera can provide for you.  When I did this, I just used a reactor with a moving wall (the piston) and a thermal conductivity function.  Things like internal energy and specific heat should be readily addressed by the reactor object.

[Shah Saud Alam]

Hello Bryan,

Yeah, I guess I am using a sledgehammer approach with the model. Would you mind telling me the equivalent functions in Cantera that can be used to handle the odes? Or, if you could point me to specific examples, that would be wonderful. I did check the ic_engine.py, and to be honest, half of the things in it just went over my head because, 1. because my model is not that complex, and 2. the syntax is in Python whereas I prefer Matlab.

Also, if possible, would you mind sharing a MWE of the reactor and wall model that you did earlier? I might be able to compare the two models and build upon the information I can gather.

Thanks a lot for your time, Bryan. Truly, appreciate it.

Best Regards

Saud

[Bryan Callaway]

I don't have a working copy of Matlab any more and it's been a while.  The core in-cylinder code looks like this in Python.  ## are notes that I have added here, as opposed to those in the original code, which is a couple of years old now.

I wrote functions for piston speed, volume, heat transfer, etc:

def piston_speed(t):
  scaler = omega*stroke/2
  a_2=(1/(4*rod_ratio)+1/(16*rod_ratio**3)+15/(512*rod_ratio**5))
  a_4=(1/(64*rod_ratio**3)+3/(256*rod_ratio**5))
  speed = scaler*(math.sin(omega*t)-2*a_2*math.sin(omega*2*t)+4*a_4*math.sin(omega*4*t))
  return  speed ##This uses a third order Fourier function for piston speed as outlined in Hoag, "vehicular engine design".  As I recall, there is an error in Hoag's formula, corrected here.

def piston_position(t):
  sin_alpha=math.sin(math.pi-omega*t)/rod_ratio
  cos_alpha=math.sqrt(1-(sin_alpha)**2)
  len=stroke*rod_ratio/2
  rad = stroke/2
  X=rad+len-len*cos_alpha+rad*math.cos(omega*t)
  return X #Positive velocity compresses. ##This position corresponds to the above velocity.

def cylinder_volume(t):
  return piston_position(t)*area_piston+V_TDC

def heat_transfer(phase):
  P=phase.P
  T=phase.T
  V=cylinder_volume(t)
  Q=C_1*(V**(-.06))*((P/100000)**(0.8))*(T**(-.4))*(RPM*stroke/30+C_2)**(.8)
  return Q

## I believe that this heat transfer model is after Chang, SAE 2004-01-2996.

There are some other functions for things like ringing index, we define some boundary conditions, and we initialize an Excel file into which to output data.  Then we build a reactor net:

#Reservoir into which heat is transferred through "bore ends"; also the "other side" of piston travel.
end_gas = ct.Solution('Air.xml')
end_gas.TP=(T_ends, 100000)
end_reservoir = ct.Reservoir(end_gas)

#Reservoir into which heat is transferred through the cylinder liner.
liner_gas = ct.Solution('Air.xml')
liner_gas.TP=(T_liner, 100000)
liner_reservoir = ct.Reservoir(liner_gas)

Stripping out a lot of work to define initial conditions, set reference conditions for total energy flow, etc.,

charge.TPX=(T_0, P_0, composition_blended)

cylinder = ct.IdealGasReactor(charge)
cylinder.volume = V_BDC

piston = ct.Wall(end_reservoir, cylinder, A=area_piston, velocity=piston_speed) #Moving boundary for compression / expansion.  No heat transfer occurs through this wall.
ends = ct.Wall(end_reservoir, cylinder, A=area_piston*end_area_multiplier) ##This wall exists exclusively to model heat transfer through the head and piston.  Note that the area is effectively divided by two; heat flux is doubled.
liner = ct.Wall(liner_reservoir, cylinder, A=(stroke+piston_clearance)*math.pi*bore/2) ##As does this, for the variable area liner.
net = ct.ReactorNet([cylinder])
net.set_max_time_step(t_final/minimum_steps) ##Sets a maximum time step.  The solver is free to set a smaller one, which is of course important for the very fast portion of combustion.  Looking at your code, which appears to have a fixed time step, it's possible that that's your problem.
t=0

Then we start to use it:

while (t<t_final):
    t_previous=t #this records the temperature at the end of the last time step, or 0 at the beginning of each case.  This lets us know the length of the time step, which is useful for indicated work and so forth.
    t=net.step() #this takes the reactor "net" through one time step of kinetics.  T advances accordingly.
    liner.area=(piston_position(t)+piston_clearance)*math.pi*bore/2 #This calls a function to determine piston position, and sets the area of a reactor wall ("liner") accordingly.  The /2 term: I was using the same wall ("piston") for both changing the reactor volume and the fixed part of the heat transfer area.  So I didn't have to also create a wall for the head, I doubled the heat transfer coefficient and then divided the wall area by two.  It made sense a couple of years ago.

    Q=heat_transfer(charge) #This calls another function which determines the correct heat flux.
    ends.heat_transfer_coeff=Q #This sets the heat transfer coefficient through a reactor boundary "ends" (the piston and head, which have constant area) equal ot the aforedetermined heat flux.
    liner.heat_transfer_coeff=Q #And the same for the reactor boundary "liner", which has a variable area as defined above.

#There are then a lot of rows recording pressure, heat flux, work, etc. from this time step in an Excel sheet for future reference.  At the next time step, the reactor net will see the piston move a certain distance, kinetics proceed a certain amount, heat transfer according to the "Q" variable, and so forth; then everything is recorded and we go again.  Note that this setup allows the solver to determine the time step, although I think there's a maximum time step somewhere in my initialization code.

Everything in here should also be available in Matlab.

Hope that helps.

Bryan

To unsubscribe from this group and all its topics, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/0b21cd87-79e2-4e1a-9c66-302e137c244co%40googlegroups.com.

[Ingmar Schoegl]

I’m happy that you were able to confirm. The discrepancy to Chemkin is odd as the mechanism shouldn’t behave differently (as long as the physics remain the same, and you’re using the same rpm’s). Without a direct comparison it’s difficult to pinpoint, though.

Otherwise, I’m with Bryan: many of the items you need are already part of Cantera.

-Ingmar-

To unsubscribe from this group and all its topics, send an email to cantera-user...@googlegroups.com.


To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/23973a75-1e84-46b7-903b-162bede620ceo%40googlegroups.com.