GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

[Examples] Formatting changes

Browse Source 
Add requirements to all Python examples. Update formatting to resolve
some PEP8 problems.
This commit is contained in:
Bryan W. Weber
2019-12-29 18:31:43 -05:00
committed by Ray Speth
parent d020feb668
commit 6ccefcbcd7
63 changed files with 356 additions and 285 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
interfaces/cython/cantera/examples/kinetics/extract_submechanism.py
View File

@@ -6,11 +6,12 @@ oxidation reactions are extracted from the GRI 3.0 mechanism.
To test the submechanism, a premixed CO/H2 flame is simulated using the original
mechanism and the submechanism, which demonstrates that the submechanism
contains all of the important species and reactions.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
from timeit import default_timer
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
all_species = ct.Species.listFromFile('gri30.yaml')
@@ -22,7 +23,7 @@ for S in all_species:
if 'C' in comp and 'H' in comp:
# Exclude all hydrocarbon species
continue
if 'N' in comp and comp != {'N':2}:
if 'N' in comp and comp != {'N': 2}:
# Exclude all nitrogen compounds except for N2
continue
if 'Ar' in comp:
@@ -54,19 +55,21 @@ gas1 = ct.Solution('gri30.yaml')
gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
species=species, reactions=reactions)
def solve_flame(gas):
gas.TPX = 373, 0.05*ct.one_atm, 'H2:0.4, CO:0.6, O2:1, N2:3.76'
# Create the flame simulation object
sim = ct.CounterflowPremixedFlame(gas=gas, width=0.2)
sim.reactants.mdot = 0.12 # kg/m^2/s
sim.products.mdot = 0.06 # kg/m^2/s
sim.reactants.mdot = 0.12 # kg/m^2/s
sim.products.mdot = 0.06 # kg/m^2/s
sim.set_refine_criteria(ratio=3, slope=0.1, curve=0.2)
sim.solve(0, auto=True)
return sim
t1 = default_timer()
sim1 = solve_flame(gas1)
t2 = default_timer()

8
interfaces/cython/cantera/examples/kinetics/mechanism_reduction.py
View File

@@ -11,6 +11,8 @@ We then create a sequence of reduced mechanisms including only the top reactions
and the associated species, and run the simulations again with these mechanisms
to see whether the reduced mechanisms with a certain number of species are able
to adequately simulate the ignition delay problem.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
@@ -44,8 +46,8 @@ plt.plot(tt, TT, label='K=53, R=325', color='k', lw=3, zorder=100)
R = sorted(zip(Rmax, gas.reactions()), key=lambda x: -x[0])
# Test reduced mechanisms with different numbers of reactions
C = plt.cm.winter(np.linspace(0,1,5))
for i,N in enumerate([40,50,60,70,80]):
C = plt.cm.winter(np.linspace(0, 1, 5))
for i, N in enumerate([40, 50, 60, 70, 80]):
# Get the N most active reactions
reactions = [r[1] for r in R[:N]]
@@ -77,7 +79,7 @@ for i,N in enumerate([40,50,60,70,80]):
tt.append(1000 * t)
TT.append(r.T)
plt.plot(tt,TT, lw=2, color=C[i],
plt.plot(tt, TT, lw=2, color=C[i],
label='K={0}, R={1}'.format(gas2.n_species, N))
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')

4
interfaces/cython/cantera/examples/kinetics/reaction_path.py
View File

@@ -3,8 +3,10 @@ Viewing a reaction path diagram.
This script uses Graphviz to generate an image. You must have Graphviz installed
and the program 'dot' must be on your path for this example to work.
Graphviz can be obtained from http://www.graphviz.org/ or (possibly) installed
Graphviz can be obtained from https://www.graphviz.org/ or (possibly) installed
using your operating system's package manager.
Requires: cantera >= 2.5.0
"""
import os

6
interfaces/cython/cantera/examples/multiphase/adiabatic.py
View File

@@ -1,6 +1,8 @@
"""
Adiabatic flame temperature and equilibrium composition for a fuel/air mixture
as a function of equivalence ratio, including formation of solid carbon.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
@@ -35,7 +37,7 @@ mix = ct.Mixture(mix_phases)
# create some arrays to hold the data
tad = np.zeros(npoints)
xeq = np.zeros((mix.n_species,npoints))
xeq = np.zeros((mix.n_species, npoints))
for i in range(npoints):
# set the gas state
@@ -51,7 +53,7 @@ for i in range(npoints):
tad[i] = mix.T
print('At phi = {0:12.4g}, Tad = {1:12.4g}'.format(phi[i], tad[i]))
xeq[:,i] = mix.species_moles
xeq[:, i] = mix.species_moles
# write output CSV file for importing into Excel
csv_file = 'adiabatic.csv'

8
interfaces/cython/cantera/examples/multiphase/plasma_equilibrium.py
View File

@@ -1,6 +1,8 @@
"""
An equilibrium example with charged species in the gas phase
and multiple condensed phases.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
@@ -11,9 +13,9 @@ import csv
# as incompressible stoichiometric substances. See file KOH.yaml for more
# information.
phases = ct.import_phases('KOH.yaml', ['K_solid', 'K_liquid', 'KOH_a', 'KOH_b',
'KOH_liquid', 'K2O2_solid', 'K2O_solid',
'KO2_solid', 'ice', 'liquid_water',
'KOH_plasma'])
'KOH_liquid', 'K2O2_solid', 'K2O_solid',
'KO2_solid', 'ice', 'liquid_water',
'KOH_plasma'])
# create the Mixture object from the list of phases
mix = ct.Mixture(phases)

7
interfaces/cython/cantera/examples/onedim/adiabatic_flame.py
View File

@@ -1,10 +1,11 @@
"""
A freely-propagating, premixed hydrogen flat flame with multicomponent
transport properties.
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
# Simulation parameters
p = ct.one_atm # pressure [Pa]
@@ -34,10 +35,10 @@ print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.u[0]))
# Solve with multi-component transport properties
f.transport_model = 'Multi'
f.solve(loglevel) # don't use 'auto' on subsequent solves
f.solve(loglevel) # don't use 'auto' on subsequent solves
f.show_solution()
print('multicomponent flamespeed = {0:7f} m/s'.format(f.u[0]))
f.save('h2_adiabatic.xml','multi', 'solution with multicomponent transport')
f.save('h2_adiabatic.xml', 'multi', 'solution with multicomponent transport')
# write the velocity, temperature, density, and mole fractions to a CSV file
f.write_csv('h2_adiabatic.csv', quiet=False)

7
interfaces/cython/cantera/examples/onedim/burner_flame.py
View File

@@ -1,15 +1,16 @@
"""
A burner-stabilized lean premixed hydrogen-oxygen flame at low pressure.
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
p = 0.05 * ct.one_atm
tburner = 373.0
mdot = 0.06
reactants = 'H2:1.5, O2:1, AR:7' # premixed gas composition
width = 0.5 # m
width = 0.5 # m
loglevel = 1 # amount of diagnostic output (0 to 5)
gas = ct.Solution('h2o2.yaml')
@@ -25,7 +26,7 @@ f.solve(loglevel, auto=True)
f.save('h2_burner_flame.xml', 'mix', 'solution with mixture-averaged transport')
f.transport_model = 'Multi'
f.solve(loglevel) # don't use 'auto' on subsequent solves
f.solve(loglevel) # don't use 'auto' on subsequent solves
f.show_solution()
f.save('h2_burner_flame.xml', 'multi', 'solution with multicomponent transport')

6
interfaces/cython/cantera/examples/onedim/burner_ion_flame.py
View File

@@ -1,14 +1,15 @@
"""
A burner-stabilized lean premixed hydrogen-oxygen flame at low pressure.
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
p = ct.one_atm
tburner = 600.0
reactants = 'CH4:1.0, O2:2.0, N2:7.52' # premixed gas composition
width = 0.5 # m
width = 0.5 # m
loglevel = 1 # amount of diagnostic output (0 to 5)
gas = ct.Solution('gri30_ion.yaml')
@@ -26,4 +27,3 @@ f.solve(loglevel=loglevel, stage=2, enable_energy=True)
f.save('CH4_burner_flame.xml', 'mix', 'solution with mixture-averaged transport')
f.write_csv('CH4_burner_flame.csv', quiet=False)

5
interfaces/cython/cantera/examples/onedim/diffusion_flame.py
View File

@@ -1,9 +1,10 @@
"""
An opposed-flow ethane/air diffusion flame
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
# Input parameters
@@ -16,7 +17,7 @@ mdot_f = 0.24 # kg/m^2/s
comp_o = 'O2:0.21, N2:0.78, AR:0.01' # air composition
comp_f = 'C2H6:1' # fuel composition
width = 0.02 # Distance between inlets is 2 cm
width = 0.02 # Distance between inlets is 2 cm
loglevel = 1 # amount of diagnostic output (0 to 5)

22
interfaces/cython/cantera/examples/onedim/diffusion_flame_batch.py
View File

@@ -1,5 +1,3 @@
# -*- coding: utf-8 -*-
# This file is part of Cantera. See License.txt in the top-level directory or
# at https://cantera.org/license.txt for license and copyright information.
@@ -14,14 +12,18 @@ explanation. Also, please don't forget to cite it if you make use of it.
This example can, for example, be used to iterate to a counterflow diffusion flame to an
awkward pressure and strain rate, or to create the basis for a flamelet table.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
import numpy as np
import os
import cantera as ct
import matplotlib.pyplot as plt
class FlameExtinguished(Exception):
class FlameExtinguished(Exception):
pass
@@ -37,7 +39,7 @@ if not os.path.exists(data_directory):
reaction_mechanism = 'h2o2.yaml'
gas = ct.Solution(reaction_mechanism)
width = 18e-3 # 18mm wide
width = 18e-3 # 18mm wide
f = ct.CounterflowDiffusionFlame(gas, width=width)
# Define the operating pressure and boundary conditions
@@ -186,12 +188,10 @@ while np.max(f.T) > temperature_limit_extinction:
# PART 4: PLOT SOME FIGURES
import matplotlib.pyplot as plt
fig1 = plt.figure()
fig2 = plt.figure()
ax1 = fig1.add_subplot(1,1,1)
ax2 = fig2.add_subplot(1,1,1)
ax1 = fig1.add_subplot(1, 1, 1)
ax2 = fig2.add_subplot(1, 1, 1)
p_selected = p_range[::7]
for p in p_selected:
@@ -218,8 +218,8 @@ fig2.savefig(data_directory + 'figure_u_p.png')
fig3 = plt.figure()
fig4 = plt.figure()
ax3 = fig3.add_subplot(1,1,1)
ax4 = fig4.add_subplot(1,1,1)
ax3 = fig3.add_subplot(1, 1, 1)
ax4 = fig4.add_subplot(1, 1, 1)
n_selected = range(1, n, 5)
for n in n_selected:
file_name = 'strain_loop_{0:02d}.xml'.format(n)

15
interfaces/cython/cantera/examples/onedim/diffusion_flame_extinction.py
View File

@@ -1,5 +1,3 @@
# -*- coding: utf-8 -*-
# This file is part of Cantera. See License.txt in the top-level directory or
# at https://cantera.org/license.txt for license and copyright information.
@@ -10,12 +8,16 @@ A hydrogen-oxygen diffusion flame at 1 bar is studied.
The tutorial makes use of the scaling rules derived by Fiala and Sattelmayer
(doi:10.1155/2014/484372). Please refer to this publication for a detailed
explanation. Also, please don't forget to cite it if you make use of it.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
import numpy as np
import os
import cantera as ct
import matplotlib.pyplot as plt
# Create directory for output data files
data_directory = 'diffusion_flame_extinction_data/'
if not os.path.exists(data_directory):
@@ -28,7 +30,7 @@ if not os.path.exists(data_directory):
reaction_mechanism = 'h2o2.yaml'
gas = ct.Solution(reaction_mechanism)
width = 18.e-3 # 18mm wide
width = 18.e-3 # 18mm wide
f = ct.CounterflowDiffusionFlame(gas, width=width)
# Define the operating pressure and boundary conditions
@@ -133,8 +135,8 @@ while True:
break
else:
# Procedure if flame extinguished but abortion criterion is not satisfied
print('Flame extinguished at n = {0}. Restoring n = {1} with alpha = {2}'.format(
n, n_last_burning, alpha[n_last_burning]))
print('Flame extinguished at n = {0}. Restoring n = {1} with '
'alpha = {2}'.format(n, n_last_burning, alpha[n_last_burning]))
# Reduce relative strain rate increase
delta_alpha = delta_alpha / delta_alpha_factor
# Restore last burning solution
@@ -157,7 +159,6 @@ print('Axial strain rate at stoichiometric surface a_stoich={0:.2e} 1/s'.format(
f.strain_rate('stoichiometric', fuel='H2')))
# Plot the maximum temperature over the maximum axial velocity gradient
import matplotlib.pyplot as plt
plt.figure()
plt.semilogx(a_max, T_max)
plt.xlabel(r'$a_{max}$ [1/s]')

4
interfaces/cython/cantera/examples/onedim/flame_fixed_T.py
View File

@@ -1,6 +1,8 @@
"""
A burner-stabilized, premixed methane/air flat flame with multicomponent
transport properties and a specified temperature profile.
Requires: cantera >= 2.5.0
"""
import cantera as ct
@@ -14,7 +16,7 @@ mdot = 0.04 # kg/m^2/s
comp = 'CH4:0.65, O2:1, N2:3.76' # premixed gas composition
# The solution domain is chosen to be 1 cm
width = 0.01 # m
width = 0.01 # m
loglevel = 1 # amount of diagnostic output (0 to 5)
refine_grid = True # 'True' to enable refinement

3
interfaces/cython/cantera/examples/onedim/flamespeed_sensitivity.py
View File

@@ -2,10 +2,11 @@
Sensitivity analysis for a freely-propagating, premixed methane-air
flame. Computes the sensitivity of the laminar flame speed with respect
to each reaction rate constant.
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
# Simulation parameters
p = ct.one_atm # pressure [Pa]

6
interfaces/cython/cantera/examples/onedim/ion_burner_flame.py
View File

@@ -1,14 +1,15 @@
"""
A burner-stabilized premixed methane-air flame with charged species.
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
p = ct.one_atm
tburner = 600.0
reactants = 'CH4:1.0, O2:2.0, N2:7.52' # premixed gas composition
width = 0.5 # m
width = 0.5 # m
loglevel = 1 # amount of diagnostic output (0 to 5)
gas = ct.Solution('gri30_ion.yaml')
@@ -26,4 +27,3 @@ f.solve(loglevel=loglevel, stage=2, enable_energy=True)
f.save('CH4_burner_flame.xml', 'mix', 'solution with mixture-averaged transport')
f.write_csv('CH4_burner_flame.csv', quiet=False)

5
interfaces/cython/cantera/examples/onedim/ion_free_flame.py
View File

@@ -1,9 +1,10 @@
"""
A freely-propagating, premixed methane-air flat flame with charged species.
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
# Simulation parameters
p = ct.one_atm # pressure [Pa]
@@ -12,7 +13,7 @@ reactants = 'CH4:1, O2:2, N2:7.52' # premixed gas composition
width = 0.05 # m
loglevel = 1 # amount of diagnostic output (0 to 8)
# IdealGasMix object used to compute mixture properties, set to the state of the
# Solution object used to compute mixture properties, set to the state of the
# upstream fuel-air mixture
gas = ct.Solution('gri30_ion.yaml')
gas.TPX = Tin, p, reactants

4
interfaces/cython/cantera/examples/onedim/premixed_counterflow_flame.py
View File

@@ -3,6 +3,8 @@ An opposed-flow premixed strained flame
This script simulates a lean hydrogen-oxygen flame stabilized in a strained
flowfield, with an opposed flow consisting of equilibrium products.
Requires: cantera >= 2.5.0
"""
import cantera as ct
@@ -15,7 +17,7 @@ mdot_products = 0.06 # kg/m^2/s
rxnmech = 'h2o2.yaml' # reaction mechanism file
comp = 'H2:1.6, O2:1, AR:7' # premixed gas composition
width = 0.2 # m
width = 0.2 # m
loglevel = 1 # amount of diagnostic output (0 to 5)
# Set up the problem

29
interfaces/cython/cantera/examples/onedim/premixed_counterflow_twin_flame.py
View File

@@ -4,12 +4,15 @@
Simulate two counter-flow jets of reactants shooting into each other. This
simulation differs from the similar premixed_counterflow_flame.py example as the
latter simulates a jet of reactants shooting into products.
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
import sys
# Differentiation function for data that has variable grid spacing Used here to
# compute normal strain-rate
def derivative(x, y):
@@ -23,6 +26,7 @@ def derivative(x, y):
return dydx
def computeStrainRates(oppFlame):
# Compute the derivative of axial velocity to obtain normal strain rate
strainRates = derivative(oppFlame.grid, oppFlame.u)
@@ -38,6 +42,7 @@ def computeStrainRates(oppFlame):
return strainRates, strainRatePoint, K
def computeConsumptionSpeed(oppFlame):
Tb = max(oppFlame.T)
@@ -46,11 +51,12 @@ def computeConsumptionSpeed(oppFlame):
integrand = oppFlame.heat_release_rate/oppFlame.cp
I = np.trapz(integrand, oppFlame.grid)
Sc = I/(Tb - Tu)/rho_u
total_heat_release = np.trapz(integrand, oppFlame.grid)
Sc = total_heat_release/(Tb - Tu)/rho_u
return Sc
# This function is called to run the solver
def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1,
ratio=2, slope=0.3, curve=0.3, prune=0.05):
@@ -73,30 +79,31 @@ def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1,
return np.max(oppFlame.T), K, strainRatePoint
# Select the reaction mechanism
gas = ct.Solution('gri30.yaml')
# Create a CH4/Air premixed mixture with equivalence ratio=0.75, and at room
# temperature and pressure.
gas.set_equivalence_ratio(0.75, 'CH4', {'O2':1.0, 'N2':3.76})
gas.set_equivalence_ratio(0.75, 'CH4', {'O2': 1.0, 'N2': 3.76})
gas.TP = 300, ct.one_atm
# Set the velocity
axial_velocity = 2.0 # in m/s
axial_velocity = 2.0 # in m/s
# Domain half-width of 2.5 cm, meaning the whole domain is 5 cm wide
width = 0.025
# Done with initial conditions
# Compute the mass flux, as this is what the Flame object requires
massFlux = gas.density * axial_velocity # units kg/m2/s
massFlux = gas.density * axial_velocity # units kg/m2/s
# Create the flame object
oppFlame = ct.CounterflowTwinPremixedFlame(gas, width=width)
# Uncomment the following line to use a Multi-component formulation. Default is
# mixture-averaged
#oppFlame.transport_model = 'Multi'
# oppFlame.transport_model = 'Multi'
# Now run the solver. The solver returns the peak temperature, strain rate and
# the point which we ascribe to the characteristic strain rate.
@@ -119,24 +126,24 @@ if '--plot' in sys.argv:
import matplotlib.pyplot as plt
plt.figure(figsize=(8,6), facecolor='white')
plt.figure(figsize=(8, 6), facecolor='white')
# Axial Velocity Plot
plt.subplot(1,2,1)
plt.subplot(1, 2, 1)
plt.plot(oppFlame.grid, oppFlame.u, 'r', lw=2)
plt.xlim(oppFlame.grid[0], oppFlame.grid[-1])
plt.xlabel('Distance (m)')
plt.ylabel('Axial Velocity (m/s)')
# Identify the point where the strain rate is calculated
plt.plot(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint],'gs')
plt.plot(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint], 'gs')
plt.annotate('Strain-Rate point',
xy=(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint]),
xytext=(0.001, 0.1),
arrowprops={'arrowstyle':'->'})
arrowprops={'arrowstyle': '->'})
# Temperature Plot
plt.subplot(1,2,2)
plt.subplot(1, 2, 2)
plt.plot(oppFlame.grid, oppFlame.T, 'b', lw=2)
plt.xlim(oppFlame.grid[0], oppFlame.grid[-1])
plt.xlabel('Distance (m)')

7
interfaces/cython/cantera/examples/onedim/stagnation_flame.py
View File

@@ -13,10 +13,11 @@ resolve the solution. This is important here, since the flamefront moves as
the mass flowrate is increased. Without using 'prune', a large number of grid
points would be concentrated upsteam of the flame, where the flamefront had
been previously. (To see this, try setting prune to zero.)
Requires: cantera >= 2.5.0
"""
import cantera as ct
import numpy as np
import os
# parameter values
@@ -32,7 +33,7 @@ rxnmech = 'h2o2.yaml' # reaction mechanism file
comp = 'H2:1.8, O2:1, AR:7' # premixed gas composition
# The solution domain is chosen to be 20 cm
width = 0.2 # m
width = 0.2 # m
loglevel = 1 # amount of diagnostic output (0 to 5)
@@ -71,7 +72,7 @@ outfile = 'stflame1.xml'
if os.path.exists(outfile):
os.remove(outfile)
for m,md in enumerate(mdot):
for m, md in enumerate(mdot):
sim.inlet.mdot = md
sim.solve(loglevel)
sim.save(outfile, 'mdot{0}'.format(m), 'mdot = {0} kg/m2/s'.format(md))

79
interfaces/cython/cantera/examples/reactors/NonIdealShockTube.py
View File

@@ -28,6 +28,7 @@ import time
import cantera as ct
print('Running Cantera version: ' + ct.__version__)
# Define the ignition delay time (IDT). This function computes the ignition
# delay from the occurrence of the peak concentration for the specified
# species.
@@ -37,8 +38,8 @@ def ignitionDelay(states, species):
# Define the reactor temperature and pressure:
reactorTemperature = 1000 # Kelvin
reactorPressure = 40.0*101325.0 # Pascals
reactorTemperature = 1000 # Kelvin
reactorPressure = 40.0*101325.0 # Pascals
# Define the gas: In this example we will choose a stoichiometric mixture of
# n-dodecane and air as the gas. For a representative kinetic model, we use:
@@ -75,7 +76,7 @@ real_gas.TP = reactorTemperature, reactorPressure
# Define the fuel, oxidizer and set the stoichiometry:
real_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26',
oxidizer={'o2':1.0, 'n2':3.76})
oxidizer={'o2': 1.0, 'n2': 3.76})
# Create a reactor object and add it to a reactor network
# In this example, this will be the only reactor in the network
@@ -90,21 +91,22 @@ t0 = time.time()
estimatedIgnitionDelayTime = 0.005
t = 0
counter = 1;
counter = 1
while(t < estimatedIgnitionDelayTime):
t = reactorNetwork.step()
if (counter%20 == 0):
if counter % 20 == 0:
# We will save only every 20th value. Otherwise, this takes too long
# Note that the species concentrations are mass fractions
timeHistory_RG.append(r.thermo.state, t=t)
counter+=1
counter += 1
# We will use the 'oh' species to compute the ignition delay
tau_RG = ignitionDelay(timeHistory_RG, 'oh')
# Toc
t1 = time.time()
print('Computed Real Gas Ignition Delay: {:.3e} seconds. Took {:3.2f}s to compute'.format(tau_RG, t1-t0))
print("Computed Real Gas Ignition Delay: {:.3e} seconds. "
"Took {:3.2f}s to compute".format(tau_RG, t1-t0))
# Ideal gas IDT calculation
@@ -116,7 +118,7 @@ ideal_gas.TP = reactorTemperature, reactorPressure
# Define the fuel, oxidizer and set the stoichiometry:
ideal_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26',
oxidizer={'o2':1.0, 'n2':3.76})
oxidizer={'o2': 1.0, 'n2': 3.76})
r = ct.Reactor(contents=ideal_gas)
reactorNetwork = ct.ReactorNet([r])
@@ -127,14 +129,14 @@ t0 = time.time()
t = 0
counter = 1;
counter = 1
while(t < estimatedIgnitionDelayTime):
t = reactorNetwork.step()
if (counter%20 == 0):
if counter % 20 == 0:
# We will save only every 20th value. Otherwise, this takes too long
# Note that the species concentrations are mass fractions
timeHistory_IG.append(r.thermo.state, t=t)
counter+=1
counter += 1
# We will use the 'oh' species to compute the ignition delay
tau_IG = ignitionDelay(timeHistory_IG, 'oh')
@@ -142,7 +144,8 @@ tau_IG = ignitionDelay(timeHistory_IG, 'oh')
# Toc
t1 = time.time()
print('Computed Ideal Gas Ignition Delay: {:.3e} seconds. Took {:3.2f}s to compute'.format(tau_IG, t1-t0))
print("Computed Ideal Gas Ignition Delay: {:.3e} seconds. "
"Took {:3.2f}s to compute".format(tau_IG, t1-t0))
print('Ideal gas error: {:2.2f} %'.format(100*(tau_IG-tau_RG)/tau_RG))
# Plot the result
@@ -156,30 +159,30 @@ plt.rcParams['font.family'] = 'serif'
# Figure illustrating the definition of ignition delay time (IDT).
plt.figure()
plt.plot(timeHistory_RG.t, timeHistory_RG('oh').Y,'-o',color='b',markersize=4)
plt.plot(timeHistory_IG.t, timeHistory_IG('oh').Y,'-o',color='r',markersize=4)
plt.plot(timeHistory_RG.t, timeHistory_RG('oh').Y, '-o', color='b', markersize=4)
plt.plot(timeHistory_IG.t, timeHistory_IG('oh').Y, '-o', color='r', markersize=4)
plt.xlabel('Time (s)')
plt.ylabel(r'OH mass fraction, $\mathdefault{Y_{OH}}$')
# Figure formatting:
plt.xlim([0,0.00055])
plt.xlim([0, 0.00055])
ax = plt.gca()
ax.annotate("",xy=(tau_RG,0.005), xytext=(0,0.005),
arrowprops=dict(arrowstyle="<|-|>",color='k',linewidth=2.0),
fontsize=14,)
plt.annotate('Ignition Delay Time (IDT)', xy=(0,0), xytext=(0.00008, 0.00525),
fontsize=16);
ax.annotate("", xy=(tau_RG, 0.005), xytext=(0, 0.005),
arrowprops=dict(arrowstyle="<|-|>", color='k', linewidth=2.0),
fontsize=14)
plt.annotate('Ignition Delay Time (IDT)', xy=(0, 0), xytext=(0.00008, 0.00525),
fontsize=16)
plt.legend(['Real Gas','Ideal Gas'], frameon=False)
plt.legend(['Real Gas', 'Ideal Gas'], frameon=False)
# If you want to save the plot, uncomment this line (and edit as you see fit):
#plt.savefig('IDT_nDodecane_1000K_40atm.pdf',dpi=350,format='pdf')
# plt.savefig('IDT_nDodecane_1000K_40atm.pdf', dpi=350, format='pdf')
# Demonstration of NTC behavior
# Let us use the reactor model to demonstrate the impacts of non-ideal behavior on IDTs in the
# Negative Temperature Coefficient (NTC) region, where observed IDTs, counter to intuition, increase
# with increasing temperature.
# Let us use the reactor model to demonstrate the impacts of non-ideal behavior on IDTs
# in the Negative Temperature Coefficient (NTC) region, where observed IDTs, counter
# to intuition, increase with increasing temperature.
# Make a list of all the temperatures at which we would like to run simulations:
T = np.array([1250, 1225, 1200, 1150, 1100, 1075, 1050, 1025, 1012.5, 1000, 987.5,
@@ -198,7 +201,7 @@ for i, temperature in enumerate(T):
reactorTemperature = temperature
real_gas.TP = reactorTemperature, reactorPressure
real_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26',
oxidizer={'o2':1.0, 'n2':3.76})
oxidizer={'o2': 1.0, 'n2': 3.76})
r = ct.Reactor(contents=real_gas)
reactorNetwork = ct.ReactorNet([r])
@@ -208,17 +211,18 @@ for i, temperature in enumerate(T):
t0 = time.time()
t = 0
counter = 0
counter = 1
while t < estimatedIgnitionDelayTimes[i]:
t = reactorNetwork.step()
if not counter % 20:
if counter % 20 == 0:
timeHistory.append(r.thermo.state, t=t)
counter += 1
tau = ignitionDelay(timeHistory, 'oh')
t1 = time.time()
print('Computed Real Gas Ignition Delay: {:.3e} seconds for T={}K. Took {:3.2f}s to compute'.format(tau, temperature, t1-t0))
print("Computed Real Gas Ignition Delay: {:.3e} seconds for T={}K. "
"Took {:3.2f}s to compute".format(tau, temperature, t1-t0))
ignitionDelays_RG[i] = tau
@@ -230,7 +234,7 @@ for i, temperature in enumerate(T):
reactorTemperature = temperature
ideal_gas.TP = reactorTemperature, reactorPressure
ideal_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26',
oxidizer={'o2':1.0, 'n2':3.76})
oxidizer={'o2': 1.0, 'n2': 3.76})
r = ct.Reactor(contents=ideal_gas)
reactorNetwork = ct.ReactorNet([r])
@@ -240,17 +244,18 @@ for i, temperature in enumerate(T):
t0 = time.time()
t = 0
counter = 0
counter = 1
while t < estimatedIgnitionDelayTimes[i]:
t = reactorNetwork.step()
if not counter % 20:
if counter % 20 == 0:
timeHistory.append(r.thermo.state, t=t)
counter += 1
tau = ignitionDelay(timeHistory, 'oh')
t1 = time.time()
print('Computed Ideal Gas Ignition Delay: {:.3e} seconds for T={}K. Took {:3.2f}s to compute'.format(tau, temperature, t1-t0))
print("Computed Ideal Gas Ignition Delay: {:.3e} seconds for T={}K. "
"Took {:3.2f}s to compute".format(tau, temperature, t1-t0))
ignitionDelays_IG[i] = tau
@@ -263,7 +268,7 @@ ax.plot(1000/T, 1e6*ignitionDelays_IG, '-.', linewidth=2.0, color='r')
ax.set_ylabel(r'Ignition Delay ($\mathdefault{\mu s}$)', fontsize=14)
ax.set_xlabel(r'1000/T (K$^\mathdefault{-1}$)', fontsize=14)
ax.set_xlim([0.8,1.2])
ax.set_xlim([0.8, 1.2])
# Add a second axis on top to plot the temperature for better readability
ax2 = ax.twiny()
@@ -271,12 +276,12 @@ ticks = ax.get_xticks()
ax2.set_xticks(ticks)
ax2.set_xticklabels((1000/ticks).round(1))
ax2.set_xlim(ax.get_xlim())
ax2.set_xlabel('Temperature (K)', fontsize=14);
ax2.set_xlabel('Temperature (K)', fontsize=14)
ax.legend(['Real Gas','Ideal Gas'], frameon=False, loc='upper left')
ax.legend(['Real Gas', 'Ideal Gas'], frameon=False, loc='upper left')
# If you want to save the plot, uncomment this line (and edit as you see fit):
#plt.savefig('NTC_nDodecane_40atm.pdf',dpi=350,format='pdf')
# plt.savefig('NTC_nDodecane_40atm.pdf', dpi=350, format='pdf')
# Show the plots.
plt.show()

6
interfaces/cython/cantera/examples/reactors/combustor.py
View File

@@ -9,6 +9,8 @@ Demonstrates the use of a MassFlowController where the mass flow rate function
depends on variables other than time by capturing these variables from the
enclosing scope. Also shows the use of a PressureController to create a constant
pressure reactor with a fixed volume.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import numpy as np
@@ -43,9 +45,11 @@ exhaust = ct.Reservoir(gas)
# can access variables defined in the calling scope, including state variables
# of the Reactor object (combustor) itself.
def mdot(t):
return combustor.mass / residence_time
inlet_mfc = ct.MassFlowController(inlet, combustor, mdot=mdot)
# A PressureController has a baseline mass flow rate matching the 'master'
@@ -73,7 +77,7 @@ while combustor.T > 500:
Q = - np.sum(states.net_production_rates * states.partial_molar_enthalpies, axis=1)
# Plot results
f, ax1 = plt.subplots(1,1)
f, ax1 = plt.subplots(1, 1)
ax1.plot(states.tres, Q, '.-', color='C0')
ax2 = ax1.twinx()
ax2.plot(states.tres[:-1], states.T[:-1], '.-', color='C1')

3
interfaces/cython/cantera/examples/reactors/custom.py
View File

@@ -8,12 +8,15 @@ Cantera is used for evaluating thermodynamic properties and kinetic rates while
an external ODE solver is used to integrate the resulting equations. In this
case, the SciPy wrapper for VODE is used, which uses the same variable-order BDF
methods as the Sundials CVODES solver used by Cantera.
Requires: cantera >= 2.5.0, scipy >= 0.19, matplotlib >= 2.0
"""
import cantera as ct
import numpy as np
import scipy.integrate
class ReactorOde:
def __init__(self, gas):
# Parameters of the ODE system and auxiliary data are stored in the

8
interfaces/cython/cantera/examples/reactors/fuel_injection.py
View File

@@ -5,6 +5,8 @@ soot precursors.
Demonstrates the use of a user-supplied function for the mass flow rate through
a MassFlowController, and the use of the SolutionArray class to store results
during reactor network integration and use these results to generate plots.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import numpy as np
@@ -25,14 +27,16 @@ gas.equilibrate('TP')
r = ct.IdealGasReactor(gas)
r.volume = 0.001 # 1 liter
# Create an inlet for the fuel, supplied as a Gaussian pulse
def fuel_mdot(t):
"""Create an inlet for the fuel, supplied as a Gaussian pulse"""
total = 3.0e-3 # mass of fuel [kg]
width = 0.5 # width of the pulse [s]
t0 = 2.0 # time of fuel pulse peak [s]
amplitude = total / (width * np.sqrt(2*np.pi))
return amplitude * np.exp(-(t-t0)**2 / (2*width**2))
mfc = ct.MassFlowController(inlet, r, mdot=fuel_mdot)
# Create the reactor network
@@ -73,7 +77,7 @@ labels = {
# Plot the concentrations of some species of interest, including PAH species
# which can be considered as precursors to soot formation.
f, ax = plt.subplots(1,2)
f, ax = plt.subplots(1, 2)
for s in ['o2', 'h2o', 'co2', 'CO', 'h2', 'ch4']:
ax[0].plot(states.t, states(s).X, label=labels.get(s, s))

27
interfaces/cython/cantera/examples/reactors/ic_engine.py
View File

@@ -6,7 +6,7 @@ The simulation uses n-Dodecane as fuel, which is injected close to top dead
center. Note that this example uses numerous simplifying assumptions and
thus serves for illustration purposes only.
Requires: Cantera >= 2.5.0, scipy >= 0.19, matplotlib >= 2.0
Requires: cantera >= 2.5.0, scipy >= 0.19, matplotlib >= 2.0
"""
import cantera as ct
@@ -77,14 +77,17 @@ V_oT = V_H / (epsilon - 1.)
A_piston = .25 * np.pi * d_piston ** 2
stroke = V_H / A_piston
def crank_angle(t):
"""Convert time to crank angle"""
return np.remainder(2 * np.pi * f * t, 4 * np.pi)
def piston_speed(t):
"""Approximate piston speed with sinusoidal velocity profile"""
return - stroke / 2 * 2 * np.pi * f * np.sin(crank_angle(t))
#####################################################################
# Set up Reactor Network
#####################################################################
@@ -150,10 +153,10 @@ cyl.set_advance_limit('temperature', delta_T_max)
#####################################################################
# set up output data arrays
states = ct.SolutionArray(cyl.thermo,
extra=('t', 'ca', 'V', 'm',
'mdot_in', 'mdot_out',
'dWv_dt', 'heat_release_rate'))
states = ct.SolutionArray(
cyl.thermo,
extra=('t', 'ca', 'V', 'm', 'mdot_in', 'mdot_out', 'dWv_dt', 'heat_release_rate'),
)
# simulate with a maximum resolution of 1 deg crank angle
dt = 1. / (360 * f)
@@ -178,6 +181,7 @@ while sim.time < t_stop:
dWv_dt=dWv_dt,
heat_release_rate=heat_release_rate)
#######################################################################
# Plot Results in matplotlib
#######################################################################
@@ -186,17 +190,18 @@ def ca_ticks(t):
"""Helper function converts time to rounded crank angle."""
return np.round(crank_angle(t) * 180 / np.pi, decimals=1)
t = states.t
# pressure and temperature
fig, ax = plt.subplots(nrows=2)
ax[0].plot(t, states.P / 1.e5)
ax[0].set_ylabel('$p$ [bar]')
ax[0].set_xlabel('$\phi$ [deg]')
ax[0].set_xlabel(r'$\phi$ [deg]')
ax[0].set_xticklabels([])
ax[1].plot(t, states.T)
ax[1].set_ylabel('$T$ [K]')
ax[1].set_xlabel('$\phi$ [deg]')
ax[1].set_xlabel(r'$\phi$ [deg]')
ax[1].set_xticklabels(ca_ticks(ax[1].get_xticks()))
plt.show()
@@ -216,12 +221,12 @@ plt.show()
# heat of reaction and expansion work
fig, ax = plt.subplots()
ax.plot(t, 1.e-3 * states.heat_release_rate, label='$\dot{Q}$')
ax.plot(t, 1.e-3 * states.dWv_dt, label='$\dot{W}_v$')
ax.plot(t, 1.e-3 * states.heat_release_rate, label=r'$\dot{Q}$')
ax.plot(t, 1.e-3 * states.dWv_dt, label=r'$\dot{W}_v$')
ax.set_ylim(-1e2, 1e3)
ax.legend(loc=0)
ax.set_ylabel('[kW]')
ax.set_xlabel('$\phi$ [deg]')
ax.set_xlabel(r'$\phi$ [deg]')
ax.set_xticklabels(ca_ticks(ax.get_xticks()))
plt.show()
@@ -233,7 +238,7 @@ ax.plot(t, states('co').X, label='CO')
ax.plot(t, states('c12h26').X * 10, label='n-Dodecane x10')
ax.legend(loc=0)
ax.set_ylabel('$X_i$ [-]')
ax.set_xlabel('$\phi$ [deg]')
ax.set_xlabel(r'$\phi$ [deg]')
ax.set_xticklabels(ca_ticks(ax.get_xticks()))
plt.show()

2
interfaces/cython/cantera/examples/reactors/mix1.py
View File

@@ -14,6 +14,8 @@ contain all species that might be present in any upstream reactor.
Compare this approach for the transient problem to the method used for the
steady-state problem in thermo/mixing.py.
Requires: cantera >= 2.5.0
"""
import cantera as ct

23
interfaces/cython/cantera/examples/reactors/periodic_cstr.py
View File

@@ -16,11 +16,12 @@ example.
Acknowledgments: The idea for this example and an estimate of the conditions
needed to see the oscillations came from Bob Kee, Colorado School of Mines
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
# create the gas mixture
gas = ct.Solution('h2o2.yaml')
@@ -56,8 +57,8 @@ w = ct.Wall(cstr, env, A=1.0, U=0.02)
# Connect the upstream reservoir to the reactor with a mass flow controller
# (constant mdot). Set the mass flow rate to 1.25 sccm.
sccm = 1.25
vdot = sccm * 1.0e-6/60.0 * ((ct.one_atm / gas.P) * ( gas.T / 273.15)) # m^3/s
mdot = gas.density * vdot # kg/s
vdot = sccm * 1.0e-6 / 60.0 * ((ct.one_atm / gas.P) * (gas.T / 273.15)) # m^3/s
mdot = gas.density * vdot # kg/s
mfc = ct.MassFlowController(upstream, cstr, mdot=mdot)
# now create a downstream reservoir to exhaust into.
@@ -73,7 +74,7 @@ network = ct.ReactorNet([cstr])
# now integrate in time
t = 0.0
dt = 0.1
dt = 0.1
states = ct.SolutionArray(gas, extra=['t'])
while t < 300.0:
@@ -83,11 +84,7 @@ while t < 300.0:
if __name__ == '__main__':
print(__doc__)
try:
import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(states.t, states('H2','O2','H2O').Y)
plt.title('Mass Fractions')
plt.show()
except ImportError:
print('Matplotlib not found. Unable to plot results.')
plt.figure(1)
plt.plot(states.t, states('H2', 'O2', 'H2O').Y)
plt.title('Mass Fractions')
plt.show()

5
interfaces/cython/cantera/examples/reactors/pfr.py
View File

@@ -3,10 +3,13 @@
This example solves a plug-flow reactor problem of hydrogen-oxygen combustion.
The PFR is computed by two approaches: The simulation of a Lagrangian fluid
particle, and the simulation of a chain of reactors.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
#######################################################################
# Input Parameters
@@ -137,8 +140,6 @@ for n in range(n_steps):
# Compare Results in matplotlib
#####################################################################
import matplotlib.pyplot as plt
plt.figure()
plt.plot(z1, states1.T, label='Lagrangian Particle')
plt.plot(z2, states2.T, label='Reactor Chain')

29
interfaces/cython/cantera/examples/reactors/piston.py
View File

@@ -16,16 +16,17 @@ The two volumes are connected by an adiabatic free piston. The piston speed is
proportional to the pressure difference between the two chambers.
Note that each side uses a *different* reaction mechanism
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import sys
import cantera as ct
fmt = '%10.3f %10.1f %10.4f %10.4g %10.4g %10.4g %10.4g'
print('%10s %10s %10s %10s %10s %10s %10s' % ('time [s]','T1 [K]','T2 [K]',
'V1 [m^3]', 'V2 [m^3]',
'V1+V2 [m^3]','X(CO)'))
fmt = '{:10.3f} {:10.1f} {:10.4f} {:10.4g} {:10.4g} {:10.4g} {:10.4g}'
print('{:10} {:10} {:10} {:10} {:10} {:10} {:10}'.format(
'time [s]', 'T1 [K]', 'T2 [K]', 'V1 [m^3]', 'V2 [m^3]', 'V1+V2 [m^3]', 'X(CO)'))
gas1 = ct.Solution('h2o2.yaml')
gas1.TPX = 900.0, ct.one_atm, 'H2:2, O2:1, AR:20'
@@ -38,6 +39,7 @@ r1.volume = 0.5
r2 = ct.IdealGasReactor(gas2)
r2.volume = 0.1
# The wall is held fixed until t = 0.1 s, then released to allow the pressure to
# equilibrate.
def v(t):
@@ -46,19 +48,20 @@ def v(t):
else:
return (r1.thermo.P - r2.thermo.P) * 1e-4
w = ct.Wall(r1, r2, velocity=v)
net = ct.ReactorNet([r1, r2])
states1 = ct.SolutionArray(r1.thermo, extra=['t','v'])
states2 = ct.SolutionArray(r2.thermo, extra=['t','v'])
states1 = ct.SolutionArray(r1.thermo, extra=['t', 'v'])
states2 = ct.SolutionArray(r2.thermo, extra=['t', 'v'])
for n in range(200):
time = (n+1)*0.001
net.advance(time)
if n % 4 == 3:
print(fmt % (time, r1.T, r2.T, r1.volume, r2.volume,
r1.volume + r2.volume, r2.thermo['CO'].X[0]))
print(fmt.format(time, r1.T, r2.T, r1.volume, r2.volume,
r1.volume + r2.volume, r2.thermo['CO'].X[0]))
states1.append(r1.thermo.state, t=1000*time, v=r1.volume)
states2.append(r2.thermo.state, t=1000*time, v=r2.volume)
@@ -66,20 +69,20 @@ for n in range(200):
# plot the results if matplotlib is installed.
if '--plot' in sys.argv:
import matplotlib.pyplot as plt
plt.subplot(2,2,1)
plt.subplot(2, 2, 1)
plt.plot(states1.t, states1.T, '-', states2.t, states2.T, 'r-')
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2,2,2)
plt.plot(states1.t, states1.v,'-', states2.t, states2.v, 'r-',
plt.subplot(2, 2, 2)
plt.plot(states1.t, states1.v, '-', states2.t, states2.v, 'r-',
states1.t, states1.v + states2.v, 'g-')
plt.xlabel('Time (ms)')
plt.ylabel('Volume (m3)')
plt.subplot(2,2,3)
plt.subplot(2, 2, 3)
plt.plot(states2.t, states2('CO').X)
plt.xlabel('Time (ms)')
plt.ylabel('CO Mole Fraction (right)')
plt.subplot(2,2,4)
plt.subplot(2, 2, 4)
plt.plot(states1.t, states1('H2').X)
plt.xlabel('Time (ms)')
plt.ylabel('H2 Mole Fraction (left)')

16
interfaces/cython/cantera/examples/reactors/reactor1.py
View File

@@ -1,11 +1,10 @@
"""
Constant-pressure, adiabatic kinetics simulation.
Requires: Cantera >= 2.5.0
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import sys
import numpy as np
import cantera as ct
@@ -24,12 +23,13 @@ dt_max = 1.e-5
t_end = 100 * dt_max
states = ct.SolutionArray(gas, extra=['t'])
print('{:10s} {:10s} {:10s} {:14s}'.format('t [s]','T [K]','P [Pa]','u [J/kg]'))
print('{:10s} {:10s} {:10s} {:14s}'.format(
't [s]', 'T [K]', 'P [Pa]', 'u [J/kg]'))
while sim.time < t_end:
sim.advance(sim.time + dt_max)
states.append(r.thermo.state, t=sim.time*1e3)
print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(sim.time, r.T,
r.thermo.P, r.thermo.u))
print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(
sim.time, r.T, r.thermo.P, r.thermo.u))
# Plot the results if matplotlib is installed.
# See http://matplotlib.org/ to get it.
@@ -41,15 +41,15 @@ if '--plot' in sys.argv[1:]:
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2, 2, 2)
plt.plot(states.t, states.X[:,gas.species_index('OH')])
plt.plot(states.t, states.X[:, gas.species_index('OH')])
plt.xlabel('Time (ms)')
plt.ylabel('OH Mole Fraction')
plt.subplot(2, 2, 3)
plt.plot(states.t, states.X[:,gas.species_index('H')])
plt.plot(states.t, states.X[:, gas.species_index('H')])
plt.xlabel('Time (ms)')
plt.ylabel('H Mole Fraction')
plt.subplot(2, 2, 4)
plt.plot(states.t, states.X[:,gas.species_index('H2')])
plt.plot(states.t, states.X[:, gas.species_index('H2')])
plt.xlabel('Time (ms)')
plt.ylabel('H2 Mole Fraction')
plt.tight_layout()

29
interfaces/cython/cantera/examples/reactors/reactor2.py
View File

@@ -3,29 +3,28 @@ Two reactors connected with a piston, with heat loss to the environment
This script simulates the following situation. A closed cylinder with volume 2
m^3 is divided into two equal parts by a massless piston that moves with speed
proportional to the pressure difference between the two sides. It is
proportional to the pressure difference between the two sides. It is
initially held in place in the middle. One side is filled with 1000 K argon at
20 atm, and the other with a combustible 500 K methane/air mixture at 0.1 atm
(phi = 1.1). At t = 0 the piston is released and begins to move due to the
large pressure difference, compressing and heating the methane/air mixture,
which eventually explodes. At the same time, the argon cools as it expands.
The piston is adiabatic, but some heat is lost through the outer cylinder
walls to the environment.
The piston allows heat transfer between the reactors and some heat is lost
through the outer cylinder walls to the environment.
Note that this simulation, being zero-dimensional, takes no account of shock
wave propagation. It is somewhat artifical, but nevertheless instructive.
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import sys
import os
import csv
import numpy as np
import cantera as ct
#-----------------------------------------------------------------------
# First create each gas needed, and a reactor or reservoir for each one.
#-----------------------------------------------------------------------
# create an argon gas object and set its state
ar = ct.Solution('argon.yaml')
@@ -45,16 +44,14 @@ gas.set_equivalence_ratio(1.1, 'CH4:1.0', 'O2:2, N2:7.52')
# create a reactor for the methane/air side
r2 = ct.IdealGasReactor(gas)
#-----------------------------------------------------------------------------
# Now couple the reactors by defining common walls that may move (a piston) or
# conduct heat
#-----------------------------------------------------------------------------
# add a flexible wall (a piston) between r2 and r1
w = ct.Wall(r2, r1, A=1.0, K=0.5e-4, U=100.0)
# heat loss to the environment. Heat loss always occur through walls, so we
# create a wall separating r1 from the environment, give it a non-zero area,
# create a wall separating r2 from the environment, give it a non-zero area,
# and specify the overall heat transfer coefficient through the wall.
w2 = ct.Wall(r2, env, A=1.0, U=500.0)
@@ -94,21 +91,21 @@ print('Directory: '+os.getcwd())
if '--plot' in sys.argv:
import matplotlib.pyplot as plt
plt.clf()
plt.subplot(2,2,1)
plt.subplot(2, 2, 1)
h = plt.plot(states1.t, states1.T, 'g-', states2.t, states2.T, 'b-')
#plt.legend(['Reactor 1','Reactor 2'],2)
# plt.legend(['Reactor 1','Reactor 2'], 2)
plt.xlabel('Time (s)')
plt.ylabel('Temperature (K)')
plt.subplot(2,2,2)
plt.subplot(2, 2, 2)
plt.plot(states1.t, states1.P / 1e5, 'g-', states2.t, states2.P / 1e5, 'b-')
#plt.legend(['Reactor 1','Reactor 2'],2)
# plt.legend(['Reactor 1','Reactor 2'], 2)
plt.xlabel('Time (s)')
plt.ylabel('Pressure (Bar)')
plt.subplot(2,2,3)
plt.plot(states1.t, states1.V, 'g-', states2.t, states2.V,'b-')
#plt.legend(['Reactor 1','Reactor 2'],2)
plt.subplot(2, 2, 3)
plt.plot(states1.t, states1.V, 'g-', states2.t, states2.V, 'b-')
# plt.legend(['Reactor 1','Reactor 2'], 2)
plt.xlabel('Time (s)')
plt.ylabel('Volume (m$^3$)')

23
interfaces/cython/cantera/examples/reactors/sensitivity1.py
View File

@@ -1,5 +1,7 @@
"""
Constant-pressure, adiabatic kinetics simulation with sensitivity analysis
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import sys
@@ -26,34 +28,34 @@ sim.atol = 1.0e-15
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6
states = ct.SolutionArray(gas, extra=['t','s2','s3'])
states = ct.SolutionArray(gas, extra=['t', 's2', 's3'])
for t in np.arange(0, 2e-3, 5e-6):
sim.advance(t)
s2 = sim.sensitivity('OH', 2) # sensitivity of OH to reaction 2
s3 = sim.sensitivity('OH', 3) # sensitivity of OH to reaction 3
s2 = sim.sensitivity('OH', 2) # sensitivity of OH to reaction 2
s3 = sim.sensitivity('OH', 3) # sensitivity of OH to reaction 3
states.append(r.thermo.state, t=1000*t, s2=s2, s3=s3)
print('%10.3e %10.3f %10.3f %14.6e %10.3f %10.3f' %
(sim.time, r.T, r.thermo.P, r.thermo.u, s2, s3))
print('{:10.3e} {:10.3f} {:10.3f} {:14.6e} {:10.3f} {:10.3f}'.format(
sim.time, r.T, r.thermo.P, r.thermo.u, s2, s3))
# plot the results if matplotlib is installed.
# see http://matplotlib.org/ to get it
if '--plot' in sys.argv:
import matplotlib.pyplot as plt
plt.subplot(2,2,1)
plt.subplot(2, 2, 1)
plt.plot(states.t, states.T)
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2,2,2)
plt.subplot(2, 2, 2)
plt.plot(states.t, states('OH').X)
plt.xlabel('Time (ms)')
plt.ylabel('OH Mole Fraction')
plt.subplot(2,2,3)
plt.subplot(2, 2, 3)
plt.plot(states.t, states('H').X)
plt.xlabel('Time (ms)')
plt.ylabel('H Mole Fraction')
plt.subplot(2,2,4)
plt.subplot(2, 2, 4)
plt.plot(states.t, states('CH4').X)
plt.xlabel('Time (ms)')
plt.ylabel('CH4 Mole Fraction')
@@ -61,7 +63,8 @@ if '--plot' in sys.argv:
plt.figure(2)
plt.plot(states.t, states.s2, '-', states.t, states.s3, '-g')
plt.legend([sim.sensitivity_parameter_name(2),sim.sensitivity_parameter_name(3)],'best')
plt.legend([sim.sensitivity_parameter_name(2), sim.sensitivity_parameter_name(3)],
'best')
plt.xlabel('Time (ms)')
plt.ylabel('OH Sensitivity')
plt.tight_layout()

9
interfaces/cython/cantera/examples/reactors/surf_pfr.py
View File

@@ -2,6 +2,8 @@
This example solves a plug flow reactor problem, where the chemistry is
surface chemistry. The specific problem simulated is the partial oxidation of
methane over a platinum catalyst in a packed bed reactor.
Requires: cantera >= 2.5.0
"""
import csv
@@ -106,10 +108,11 @@ for n in range(NReactors):
upstream.syncState()
sim.reinitialize()
sim.advance_to_steady_state()
dist = n * rlen * 1.0e3 # distance in mm
dist = n * rlen * 1.0e3 # distance in mm
if not n % 10:
print(' {0:10f} {1:10f} {2:10f} {3:10f}'.format(dist, *gas['CH4','H2','CO'].X))
if n % 10 == 0:
print(' {0:10f} {1:10f} {2:10f} {3:10f}'.format(
dist, *gas['CH4', 'H2', 'CO'].X))
# write the gas mole fractions and surface coverages vs. distance
output_data.append(

8
interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py
View File

@@ -1,14 +1,16 @@
"""
CATCOMB -- Catalytic combustion of methane on platinum.
CATCOMB -- Catalytic combustion of methane on platinum.
This script solves a catalytic combustion problem. A stagnation flow is set
up, with a gas inlet 10 cm from a platinum surface at 900 K. The lean,
premixed methane/air mixture enters at ~ 6 cm/s (0.06 kg/m2/s), and burns
premixed methane/air mixture enters at ~6 cm/s (0.06 kg/m2/s), and burns
catalytically on the platinum surface. Gas-phase chemistry is included too,
and has some effect very near the surface.
The catalytic combustion mechanism is from Deutschman et al., 26th
Symp. (Intl.) on Combustion,1996 pp. 1747-1754
Requires: cantera >= 2.5.0
"""
import numpy as np
@@ -31,7 +33,7 @@ comp1 = 'H2:0.05, O2:0.21, N2:0.78, AR:0.01'
comp2 = 'CH4:0.095, O2:0.21, N2:0.78, AR:0.01'
# The inlet/surface separation is 10 cm.
width = 0.1 # m
width = 0.1 # m
loglevel = 1 # amount of diagnostic output (0 to 5)

12
interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py
View File

@@ -4,15 +4,16 @@ A CVD example simulating growth of a diamond film
This example computes the growth rate of a diamond film according to a
simplified version of a particular published growth mechanism (see file
diamond.yaml for details). Only the surface coverage equations are solved here;
the gas composition is fixed. (For an example of coupled gas- phase and
surface, see catalytic_combustion.py.) Atomic hydrogen plays an important
the gas composition is fixed. (For an example of coupled gas-phase and
surface, see catalytic_combustion.py.) Atomic hydrogen plays an important
role in diamond CVD, and this example computes the growth rate and surface
coverages as a function of [H] at the surface for fixed temperature and [CH3].
Requires: cantera >= 2.5.0, pandas >= 0.21.0, matplotlib >= 2.0
"""
import csv
import cantera as ct
import pandas as pd
print('\n****** CVD Diamond Example ******\n')
@@ -54,6 +55,7 @@ with open('diamond.csv', 'w', newline='') as f:
try:
import matplotlib.pyplot as plt
import pandas as pd
data = pd.read_csv('diamond.csv')
plt.figure()
@@ -64,7 +66,7 @@ try:
plt.figure()
for name in data:
if name.startswith('H mole') or name.startswith('Growth'):
if name.startswith(('H mole', 'Growth')):
continue
plt.plot(data['H mole Fraction'], data[name], label=name)
@@ -73,4 +75,4 @@ try:
plt.ylabel('Coverage')
plt.show()
except ImportError:
print("Install matplotlib to plot the outputs")
print("Install matplotlib and pandas to plot the outputs")

30
interfaces/cython/cantera/examples/surface_chemistry/sofc.py
View File

@@ -16,26 +16,26 @@ results.
It is recommended that you read input file sofc.yaml before reading or running
this script!
Requires: cantera >= 2.5.0
"""
import cantera as ct
import math
import csv
import inspect
import os
ct.add_module_directory()
# parameters
T = 1073.15 # T in K
P = ct.one_atm
P = ct.one_atm # One atm in Pa
# gas compositions. Change as desired.
anode_gas_X = 'H2:0.97, H2O:0.03'
cathode_gas_X = 'O2:1.0, H2O:0.001'
# time to integrate coverage eqs. to steady state in
# 'advanceCoverages'. This should be more than enough time.
# 'advance_coverages'. This should be more than enough time.
tss = 50.0
sigma = 2.0 # electrolyte conductivity [Siemens / m]
@@ -54,7 +54,7 @@ def show_coverages(s):
def equil_OCV(gas1, gas2):
return (-ct.gas_constant * gas1.T *
math.log(gas1['O2'].X / gas2['O2'].X) / (4.0*ct.faraday))
math.log(gas1['O2'].X / gas2['O2'].X) / (4.0*ct.faraday))
def NewtonSolver(f, xstart, C=0.0):
@@ -84,14 +84,14 @@ def NewtonSolver(f, xstart, C=0.0):
n += 1
raise Exception('no root!')
#####################################################################
# Anode-side phases
#####################################################################
# import the anode-side bulk phases
gas_a, anode_bulk, oxide_a = ct.import_phases(
'sofc.yaml', ['gas', 'metal', 'oxide_bulk'],
)
'sofc.yaml', ['gas', 'metal', 'oxide_bulk'])
# import the surfaces on the anode side
anode_surf = ct.Interface('sofc.yaml', 'metal_surface', [gas_a])
@@ -140,25 +140,24 @@ def anode_curr(E):
# import the cathode-side bulk phases
gas_c, cathode_bulk, oxide_c = ct.import_phases(
'sofc.yaml', ['gas', 'metal', 'oxide_bulk']
)
'sofc.yaml', ['gas', 'metal', 'oxide_bulk'])
# import the surfaces on the cathode side
cathode_surf = ct.Interface('sofc.yaml', 'metal_surface', [gas_c])
oxide_surf_c = ct.Interface('sofc.yaml', 'oxide_surface', [gas_c, oxide_c])
# import the cathode-side triple phase boundary
tpb_c = ct.Interface(
'sofc.yaml', 'tpb', [cathode_bulk, cathode_surf, oxide_surf_c]
)
tpb_c = ct.Interface('sofc.yaml', 'tpb', [cathode_bulk, cathode_surf, oxide_surf_c])
cathode_surf.name = 'cathode surface'
oxide_surf_c.name = 'cathode-side oxide surface'
def cathode_curr(E):
"""Current to the cathode as a function of cathode
potential relative to electrolyte"""
"""
Current to the cathode as a function of cathode
potential relative to electrolyte
"""
# due to ohmic losses, the cathode-side electrolyte potential is non-zero.
# Therefore, we need to add this potential to E to get the cathode
@@ -176,8 +175,8 @@ def cathode_curr(E):
# being drawn from the cathode (i.e, negative production rate).
return -ct.faraday * w[0] * TPB_length_per_area
# initialization
# initialization
# set the gas compositions, and temperatures of all phases
gas_a.TPX = T, P, anode_gas_X
@@ -200,7 +199,6 @@ for s in [anode_surf, oxide_surf_a, cathode_surf, oxide_surf_c]:
s.advance_coverages(tss)
show_coverages(s)
# find open circuit potentials by solving for the E values that give
# zero current.
Ea0 = NewtonSolver(anode_curr, xstart=-0.51)

12
interfaces/cython/cantera/examples/thermo/isentropic.py
View File

@@ -1,5 +1,7 @@
"""
Isentropic, adiabatic flow example - calculate area ratio vs. Mach number curve
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""
import cantera as ct
@@ -35,7 +37,7 @@ def isentropic(gas=None):
mdot = 1 # arbitrary
amin = 1.e14
data = np.zeros((200,4))
data = np.zeros((200, 4))
# compute values for a range of pressure ratios
for r in range(200):
@@ -48,19 +50,19 @@ def isentropic(gas=None):
v = math.sqrt(v2)
area = mdot/(gas.density*v) # rho*v*A = constant
amin = min(amin, area)
data[r,:] = [area, v/soundspeed(gas), gas.T, p/p0]
data[r, :] = [area, v/soundspeed(gas), gas.T, p/p0]
data[:,0] /= amin
data[:, 0] /= amin
return data
if __name__ == "__main__":
print(isentropic.__doc__)
print(__doc__)
data = isentropic()
try:
import matplotlib.pyplot as plt
plt.plot(data[:,1], data[:,0])
plt.plot(data[:, 1], data[:, 0])
plt.ylabel('Area Ratio')
plt.xlabel('Mach Number')
plt.title('Isentropic Flow: Area Ratio vs. Mach Number')

2
interfaces/cython/cantera/examples/thermo/mixing.py
View File

@@ -6,6 +6,8 @@ is a simpler, steady-state version of the example ``reactors/mix1.py``.
Since the goal is to simulate a continuous flow system, the mixing takes place
at constant enthalpy and pressure.
Requires: cantera >= 2.5.0
"""
import cantera as ct

2
interfaces/cython/cantera/examples/thermo/sound_speed.py
View File

@@ -1,5 +1,7 @@
"""
Compute the "equilibrium" and "frozen" sound speeds for a gas
Requires: cantera >= 2.5.0
"""
import cantera as ct

6
interfaces/cython/cantera/examples/transport/dusty_gas.py
View File

@@ -5,6 +5,8 @@ The Dusty Gas model is a multicomponent transport model for gas transport
through the pores of a stationary porous medium. This example shows how to
create a transport manager that implements the Dusty Gas model and use it to
compute the multicomponent diffusion coefficients.
Requires: cantera >= 2.5.0
"""
import cantera as ct
@@ -14,7 +16,9 @@ import cantera as ct
g = ct.DustyGas('h2o2.yaml')
# set the gas state
g.TPX = 500.0, ct.one_atm, "OH:1, H:2, O2:3, O:1.0E-8, H2:1.0E-8, H2O:1.0E-8, H2O2:1.0E-8, HO2:1.0E-8, AR:1.0E-8"
composition = {"OH": 1, "H": 2, "O2": 3, "O": 1.0E-8, "H2": 1.0E-8, "H2O": 1.0E-8,
"H2O2": 1.0E-8, "HO2": 1.0E-8, "AR": 1.0E-8}
g.TPX = 500.0, ct.one_atm, composition
# set its parameters
g.porosity = 0.2

8
interfaces/cython/cantera/examples/transport/multiprocessing_viscosity.py
View File

@@ -7,6 +7,8 @@ passed between Python processes, it is necessary to set up the computation so
that each process has its own copy of the relevant Cantera objects. One way to
do this is by storing the objects in (module) global variables, which are
initialized once per worker process.
Requires: cantera >= 2.5.0
"""
import multiprocessing
@@ -18,6 +20,7 @@ from time import time
# Global storage for Cantera Solution objects
gases = {}
def init_process(mech):
"""
This function is called once for each process in the Pool. We use it to
@@ -26,6 +29,7 @@ def init_process(mech):
gases[mech] = ct.Solution(mech)
gases[mech].transport_model = 'Multi'
def get_thermal_conductivity(args):
# Pool.imap only permits a single argument, so we pack all of the needed
# arguments into the tuple 'args'
@@ -34,6 +38,7 @@ def get_thermal_conductivity(args):
gas.TPX = T, P, X
return gas.thermal_conductivity
def get_viscosity(args):
# Pool.imap only permits a single argument, so we pack all of the needed
# arguments into the tuple 'args'
@@ -42,6 +47,7 @@ def get_viscosity(args):
gas.TPX = T, P, X
return gas.viscosity
def parallel(mech, predicate, nProcs, nTemps):
"""
Call the function ``predicate`` on ``nProcs`` processors for ``nTemps``
@@ -60,6 +66,7 @@ def parallel(mech, predicate, nProcs, nTemps):
itertools.repeat(X)))
return y
def serial(mech, predicate, nTemps):
P = ct.one_atm
X = 'CH4:1.0, O2:1.0, N2:3.76'
@@ -71,6 +78,7 @@ def serial(mech, predicate, nTemps):
itertools.repeat(X))))
return y
if __name__ == '__main__':
nPoints = 5000
nProcs = 4

26
interfaces/cython/cantera/liquidvapor.py
View File

@@ -4,41 +4,47 @@
from . import PureFluid, _cantera
def Water():
"""
Create a `PureFluid` object using the equation of state for water and the
`WaterTransport` class for viscosity and thermal conductivity."""
`WaterTransport` class for viscosity and thermal conductivity.
"""
class WaterWithTransport(PureFluid, _cantera.Transport):
__slots__ = ()
return WaterWithTransport('liquidvapor.yaml', 'water', transport_model='Water')
def Nitrogen():
"""Create a `PureFluid` object using the equation of state for nitrogen."""
return PureFluid('liquidvapor.yaml','nitrogen')
return PureFluid('liquidvapor.yaml', 'nitrogen')
def Methane():
"""Create a `PureFluid` object using the equation of state for methane."""
return PureFluid('liquidvapor.yaml','methane')
return PureFluid('liquidvapor.yaml', 'methane')
def Hydrogen():
"""Create a `PureFluid` object using the equation of state for hydrogen."""
return PureFluid('liquidvapor.yaml','hydrogen')
return PureFluid('liquidvapor.yaml', 'hydrogen')
def Oxygen():
"""Create a `PureFluid` object using the equation of state for oxygen."""
return PureFluid('liquidvapor.yaml','oxygen')
return PureFluid('liquidvapor.yaml', 'oxygen')
def Hfc134a():
"""Create a `PureFluid` object using the equation of state for HFC-134a."""
return PureFluid('liquidvapor.yaml','hfc134a')
return PureFluid('liquidvapor.yaml', 'hfc134a')
def CarbonDioxide():
"""Create a `PureFluid` object using the equation of state for carbon dioxide."""
return PureFluid('liquidvapor.yaml','carbondioxide')
return PureFluid('liquidvapor.yaml', 'carbondioxide')
def Heptane():
"""Create a `PureFluid` object using the equation of state for heptane."""
return PureFluid('liquidvapor.yaml','heptane')
return PureFluid('liquidvapor.yaml', 'heptane')

8
interfaces/cython/cantera/test/test_onedim.py
View File

@@ -122,7 +122,7 @@ class TestFreeFlame(utilities.CanteraTest):
tol_ts = [1.0e-4, 1.0e-11] # [rtol atol] for time stepping
def create_sim(self, p, Tin, reactants, width=0.05, mech='h2o2.xml'):
# IdealGasMix object used to compute mixture properties
# Solution object used to compute mixture properties
self.gas = ct.Solution(mech)
self.gas.TPX = Tin, p, reactants
@@ -567,7 +567,7 @@ class TestDiffusionFlame(utilities.CanteraTest):
def create_sim(self, p, fuel='H2:1.0, AR:1.0', T_fuel=300, mdot_fuel=0.24,
oxidizer='O2:0.2, AR:0.8', T_ox=300, mdot_ox=0.72, width=0.02):
# IdealGasMix object used to compute mixture properties
# Solution object used to compute mixture properties
self.gas = ct.Solution('h2o2.xml', 'ohmech')
self.gas.TP = T_fuel, p
@@ -952,7 +952,7 @@ class TestIonFreeFlame(utilities.CanteraTest):
Tin = 300
width = 0.03
# IdealGasMix object used to compute mixture properties
# Solution object used to compute mixture properties
self.gas = ct.Solution('ch4_ion.cti')
self.gas.TPX = Tin, p, reactants
self.sim = ct.IonFreeFlame(self.gas, width=width)
@@ -978,7 +978,7 @@ class TestIonBurnerFlame(utilities.CanteraTest):
Tburner = 400
width = 0.01
# IdealGasMix object used to compute mixture properties
# Solution object used to compute mixture properties
self.gas = ct.Solution('ch4_ion.cti')
self.gas.TPX = Tburner, p, reactants
self.sim = ct.IonBurnerFlame(self.gas, width=width)

7
samples/f77/ctlib.f
View File

@@ -26,7 +26,7 @@ c example driver program
c Read in the reaction mechanism. Since this is done differently
c than in Chemkin, this function does not correspond to any CKLIB
c subroutine.
c subroutine.
call newIdealGasMix('gri30.yaml','gri30','')
c get the number of elements, species, and reactions
@@ -39,7 +39,7 @@ c get the number of elements, species, and reactions
c compute the net production rates in cgs units
p = 1.0d6
t = 2500.0d0
call ctwyp(p, t, y, ickwrk, rckwrk, wdot)
do k = 1, kk
write(*,*) k, y(k), wdot(k)
@@ -77,7 +77,7 @@ c temperature, and array of mass fractions.
implicit double precision (a-h,o-z)
double precision y(*), rckwrk(*), wdot(*)
integer ickwrk(*)
c set the state
psi = 0.1*p
call setState_TPY(t, psi, y)
@@ -92,4 +92,3 @@ c convert SI -> cgs
end do
return
end

6
samples/f77/demo.f
View File

@@ -2,7 +2,7 @@ c
c Replace this sample main program with your program
c
c This program uses functions defined in demo_ftnlib.cpp to create
c an ideal gas mixture and print some its properties.
c an ideal gas mixture and print some of its properties.
c
c For a C++ version of this program, see ../cxx/demo.cpp.
c
@@ -27,8 +27,8 @@ c
write(*,*) 'Initial state properties:'
write(*,10) temperature(), pressure(), density(),
$ enthalpy_mole(), entropy_mole(), cp_mole()
c compute the equilibrium state holding the specific
c compute the equilibrium state holding the specific
c enthalpy and pressure constant
call equilibrate('HP')

3
samples/f77/isentropic.f
View File

@@ -58,6 +58,3 @@ c stagnation state properties
$ / meanMolarMass())
return
end

1
samples/f90/demo.f90
View File

@@ -119,4 +119,3 @@ subroutine demo(gas, MAXSP, MAXRXNS)
return
end subroutine demo

12
samples/matlab/equil.m
View File

@@ -4,7 +4,7 @@ function equil(g)
% This example computes the adiabatic flame temperature and
% equilibrium composition for a methane/air mixture as a function of
% equivalence ratio.
help equil;
help equil
if nargin == 1
gas = g;
@@ -19,13 +19,17 @@ ich4 = speciesIndex(gas,'CH4');
io2 = speciesIndex(gas,'O2');
in2 = speciesIndex(gas,'N2');
for i = 1:50
phi(i) = 0.2 + 0.05*i;
nPhis = 50;
phi = linspace(0.2, 2.70, nPhis);
tad(nPhis) = 0;
xeq(nsp,nPhis) = 0;
for i = 1:nPhis
x = zeros(nsp,1);
x(ich4,1) = phi(i);
x(io2,1) = 2.0;
x(in2,1) = 7.52;
set(gas,'Temperature',300.0,'Pressure',101325.0,'MoleFractions', x);
set(gas,'Temperature',300.0,'Pressure',101325.0,'MoleFractions',x);
equilibrate(gas,'HP');
tad(i) = temperature(gas);
xeq(:,i) = moleFractions(gas);

3
samples/matlab/flame.m
View File

@@ -49,7 +49,6 @@ rho0 = density(gas);
equilibrate(gas, 'HP');
teq = temperature(gas);
yeq = massFractions(gas);
rhoeq = density(gas);
z1 = 0.2;
mdot0 = massFlux(left);
@@ -61,8 +60,6 @@ if flametype == 0
else
t1 = temperature(right);
end
zz = gridPoints(flow);
dz = zz(end) - zz(1);
setProfile(f, 2, {'u', 'V'}, [0.0 1.0
mdot0/rho0 -mdot1/rho0
0.0 0.0]);

2
samples/matlab/flame1.m
View File

@@ -3,7 +3,7 @@
% This script simulates a burner-stablized lean hydrogen-oxygen flame
% at low pressure.
help flame1;
help flame1
t0 = cputime; % record the starting time

6
samples/matlab/flame2.m
View File

@@ -1,8 +1,6 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FLAME1 - An axisymmetric stagnation-point non-premixed flame
%
% An axisymmetric stagnation-point non-premixed flame
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This script simulates a stagnation-point ethane-air flame.
t0 = cputime; % record the starting time

6
samples/matlab/ignite.m
View File

@@ -2,7 +2,7 @@ function plotdata = ignite(g)
% IGNITE Zero-dimensional kinetics: adiabatic, constant pressure.
%
% This example solves the same problem as 'reactor1,' but does
% it using on of MATLAB's ODE integrators, rather than using the
% it using one of MATLAB's ODE integrators, rather than using the
% Cantera Reactor class.
%
@@ -14,8 +14,6 @@ else
gas = Solution('gri30.yaml');
end
nsp = nSpecies(gas);
% set the initial conditions
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
@@ -67,7 +65,7 @@ a = 1.0;
function pv = output(s, gas)
times = s.x;
soln = s.y;
[m n] = size(times);
[~,n] = size(times);
pv = zeros(nSpecies(gas) + 4, n);
set(gas,'T',1001.0,'P',oneatm);

1
samples/matlab/ignite_hp.m
View File

@@ -10,7 +10,6 @@ if nargin == 0
end
mw = molecularWeights(gas);
nsp = nSpecies(gas);
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
y0 = [temperature(gas)

1
samples/matlab/ignite_uv.m
View File

@@ -9,7 +9,6 @@ if nargin == 0
end
mw = molecularWeights(gas);
nsp = nSpecies(gas);
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
y0 = [temperature(gas)

2
samples/matlab/lithium_ion_battery.m
View File

@@ -12,7 +12,7 @@
% BinarySolutionTabulatedThermo class):
%
% Reference:
% M. Mayur, S. C. DeCaluwe, B. L. Kee, W. G. Bessler, “Modeling and simulation
% M. Mayur, S. C. DeCaluwe, B. L. Kee, W. G. Bessler, “Modeling and simulation
% of the thermodynamics of lithium-ion battery intercalation materials in the
% open-source software Cantera,” Electrochim. Acta 323, 134797 (2019),
% https://doi.org/10.1016/j.electacta.2019.134797

3
samples/matlab/prandtl1.m
View File

@@ -21,14 +21,11 @@ t = [];
xo2 = [];
io2 = speciesIndex(gas,'O2');
ih2 = speciesIndex(gas,'H2');
ih = speciesIndex(gas,'H');
ih2o = speciesIndex(gas,'H2O');
minT = minTemp(gas);
maxT = maxTemp(gas);
dT = (maxT - minT)/30.0;
atm = oneatm;
t0 = cputime;
for i = 1:31
t(i) = minT + dT*(i-1);

3
samples/matlab/prandtl2.m
View File

@@ -20,14 +20,11 @@ t = [];
xo2 = [];
io2 = speciesIndex(gas,'O2');
ih2 = speciesIndex(gas,'H2');
ih = speciesIndex(gas,'H');
ih2o = speciesIndex(gas,'H2O');
minT = minTemp(gas);
maxT = maxTemp(gas);
dT = (maxT - minT)/30.0;
atm = oneatm;
t0 = cputime;
for i = 1:31
t(i) = minT + dT*(i-1);

9
samples/matlab/rankine.m
View File

@@ -9,8 +9,6 @@ w = Water;
% start with saturated liquid water at t1
set(w,'T',t1,'Liquid',1.0);
h1 = enthalpy_mass(w);
s1 = entropy_mass(w);
p1 = pressure(w);
% pump it to p2
@@ -21,17 +19,14 @@ p2 = pressure(w);
% heat to saturated vapor
set(w,'P',p2,'Vapor',1.0);
h3 = enthalpy_mass(w);
s3 = entropy_mass(w);
heat_added = h3 - h2;
% expand adiabatically back to the initial pressure
work = expand(w, p1, eta_turbine)
h4 = enthalpy_mass(w);
x4 = vaporFraction(w);
work = expand(w, p1, eta_turbine);
% compute the efficiency
efficiency = (work - pump_work)/heat_added
efficiency = (work - pump_work)/heat_added;
function w = pump(fluid, pfinal, eta)

10
samples/matlab/reactor1.m
View File

@@ -18,9 +18,7 @@ else
gas = GRI30('None');
end
nsp = nSpecies(gas);
P = oneatm
P = oneatm;
% set the initial conditions
set(gas,'T',1001.0,'P',P,'X','H2:2,O2:1,N2:4');
@@ -47,10 +45,14 @@ setArea(w, 1.0);
% create a reactor network and insert the reactor:
network = ReactorNet({r});
nSteps = 100;
tim(nSteps) = 0;
temp(nSteps) = 0;
x(nSteps,3) = 0;
t = 0.0;
dt = 1.0e-5;
t0 = cputime;
for n = 1:100
for n = 1:nSteps
t = t + dt;
advance(network, t);
tim(n) = time(network);

6
samples/matlab/reactor2.m
View File

@@ -14,8 +14,6 @@ else
gas = GRI30('None');
end
nsp = nSpecies(gas);
% set the initial conditions
set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
@@ -25,6 +23,10 @@ r = IdealGasReactor(gas);
% create a reactor network and insert the reactor
network = ReactorNet({r});
nSteps = 100;
tim(nSteps) = 0;
temp(nSteps) = 0;
x(nSteps,3) = 0;
t = 0;
dt = 1.0e-5;
t0 = cputime;

14
samples/matlab/surfreactor.m
View File

@@ -19,6 +19,7 @@ surf = importInterface('ptcombust.yaml','Pt_surf', gas);
setTemperature(surf, t);
nsp = nSpecies(gas);
nSurfSp = nSpecies(surf);
% create a reactor, and insert the gas
r = IdealGasReactor(gas);
@@ -50,13 +51,18 @@ setExpansionRateCoeff(w, 1.0);
network = ReactorNet({r});
% setTolerances(network, 1.0e-8, 1.0e-12);
nSteps = 100;
p0 = pressure(r);
names = {'CH4','CO','CO2','H2O'};
x = zeros([nSteps 4]);
tim = zeros(nSteps);
temp = zeros(nSteps);
pres = zeros(nSteps);
cov = zeros([nSteps nSurfSp]);
t = 0;
dt = 0.1;
t0 = cputime;
p0 = pressure(r);
names = {'CH4','CO','CO2','H2O'};
x = zeros([100 4]);
for n = 1:100
for n = 1:nSteps
t = t + dt;
advance(network, t);
tim(n) = t;

2
samples/matlab/tut2.m
View File

@@ -49,7 +49,7 @@ msg = sprintf('time to create gas1b: %f', cputime - t0)
% reaction mechanisms. Under Windows, these files may be located in
% 'C:\Program Files\Common Files\Cantera', or in 'C:\cantera\data',
% depending on how you installed Cantera and the options you
% specified. On a unix/linux/Mac OSX machine, they are usually kept
% specified. On a Unix/Linux/macOS machine, they are usually kept
% in the 'data' subdirectory within the Cantera installation
% directory.

8
samples/matlab/tut4.m
View File

@@ -54,12 +54,12 @@ format short e;
for i = 1:nReactions(g)
if isReversible(g,i)
disp([i, rf(i), rr(i), (rf(i) - rr(i))/rf(i)]);
end
end
end
% You might be wondering how 'equilibrate' works. (Then again, you might
% not, in which case you can go on to the next tutorial now.) Method
% not, in which case you can go on to the next tutorial now.) Method
% 'equilibrate' invokes Cantera's chemical equilibrium solver, which
% uses an element potential method. The element potential method is
% one of a class of equivalent 'nonstoichiometric' methods that all
@@ -75,11 +75,11 @@ end
% does a few other things to generate a good starting guess and to
% produce a reasonably robust algorithm. If you want to know more
% about the details, look at the on-line documented source code of
% Cantera C++ class 'ChemEquil' at http://www.cantera.org.
% Cantera C++ class 'ChemEquil' at https://cantera.org.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
cleanup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% end of tutorial 4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2
samples/matlab/tut5.m
View File

@@ -88,7 +88,7 @@ eqs = reactionEqn(g) % all equations
kc = equil_Kc(g);
for i = 1:nReactions(g)
disp(sprintf('%50s %13.5g', eqs{i}, kc(i)))
fprintf('%50s %13.5g', eqs{i}, kc(i))
end
% 6) Multipliers

4
samples/matlab/tut7.m
View File

@@ -3,7 +3,7 @@
help tut7
% A variety of thermodynamic property methods are provided.
gas = air
gas = Air
set(gas,'T',800,'P',oneatm)
% temperature, pressure, density
@@ -30,4 +30,4 @@ gmass = gibbs_mass(gas)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
cleanup
cleanup