From ab3f5e0a63763f34cf5220ac86ee6730b1813e26 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 30 Jan 2013 22:06:02 +0000 Subject: [PATCH] [Cython] Translated some samples to use the new API --- .../cantera/examples/multiphase/adiabatic.py | 91 +++++++++++++++++++ .../examples/onedim/diffusion_flame.py | 25 ++++- .../examples/thermo/critical_properties.py | 27 ++++++ .../cantera/examples/thermo/isentropic.py | 71 +++++++++++++++ .../cython/cantera/examples/thermo/rankine.py | 75 +++++++++++++++ .../cantera/examples/thermo/sound_speed.py | 52 +++++++++++ 6 files changed, 339 insertions(+), 2 deletions(-) create mode 100644 interfaces/cython/cantera/examples/multiphase/adiabatic.py create mode 100644 interfaces/cython/cantera/examples/thermo/critical_properties.py create mode 100644 interfaces/cython/cantera/examples/thermo/isentropic.py create mode 100644 interfaces/cython/cantera/examples/thermo/rankine.py create mode 100644 interfaces/cython/cantera/examples/thermo/sound_speed.py diff --git a/interfaces/cython/cantera/examples/multiphase/adiabatic.py b/interfaces/cython/cantera/examples/multiphase/adiabatic.py new file mode 100644 index 000000000..2fcc0c161 --- /dev/null +++ b/interfaces/cython/cantera/examples/multiphase/adiabatic.py @@ -0,0 +1,91 @@ +""" +Adiabatic flame temperature and equilibrium composition for a fuel/air mixture +as a function of equivalence ratio, including formation of solid carbon. +""" + +import cantera as ct +import numpy as np +import sys +import csv + +############################################################################## +# Edit these parameters to change the initial temperature, the pressure, and +# the phases in the mixture. + +T = 300.0 +P = 101325.0 + +# phases +gas = ct.Solution('gri30.xml') +carbon = ct.Solution('graphite.xml') + +# the phases that will be included in the calculation, and their initial moles +mix_phases = [(gas, 1.0), (carbon, 0.0)] + +# gaseous fuel species +fuel_species = 'CH4' + +# air composition +air_N2_O2_ratio = 3.76 + +# equivalence ratio range +phi_min = 0.3 +phi_max = 3.5 +npoints = 50 + +############################################################################## + +mix = ct.Mixture(mix_phases) + +# create some arrays to hold the data +phi = np.zeros(npoints) +tad = np.zeros(npoints) +xeq = np.zeros((mix.n_species,npoints)) + +# find fuel, nitrogen, and oxygen indices +ifuel = gas.species_index(fuel_species) +io2 = gas.species_index('O2') +in2 = gas.species_index('N2') + +if gas.n_atoms(fuel_species,'O') > 0 or gas.n_atoms(fuel_species,'N') > 0: + raise "Error: only hydrocarbon fuels are supported." + +stoich_O2 = gas.n_atoms(fuel_species,'C') + 0.25*gas.n_atoms(fuel_species,'H') + +for i in range(npoints): + phi[i] = phi_min + (phi_max - phi_min)*i/(npoints - 1) + X = np.zeros(gas.n_species) + X[ifuel] = phi[i] + X[io2] = stoich_O2 + X[in2] = stoich_O2*air_N2_O2_ratio + + # set the gas state + gas.TPX = T, P, X + + # create a mixture of 1 mole of gas, and 0 moles of solid carbon. + mix = ct.Mixture(mix_phases) + mix.T = T + mix.P = P + + # equilibrate the mixture adiabatically at constant P + mix.equilibrate('HP', solver='gibbs', max_steps=1000) + + tad[i] = mix.T + print('At phi = {:12.4g}, Tad = {:12.4g}'.format(phi[i], tad[i])) + xeq[:,i] = mix.species_moles + +# write output CSV file for importing into Excel +csv_file = 'adiabatic.csv' +with open(csv_file, 'w') as outfile: + writer = csv.writer(outfile) + writer.writerow(['phi','T (K)'] + mix.species_names) + for i in range(npoints): + writer.writerow([phi[i], tad[i]] + list(xeq[:,i])) +print('Output written to {}'.format(csv_file)) + +if '--plot' in sys.argv: + import matplotlib.pyplot as plt + plt.plot(phi, tad) + plt.xlabel('Equivalence ratio') + plt.ylabel('Adiabatic flame temperature [K]') + plt.show() diff --git a/interfaces/cython/cantera/examples/onedim/diffusion_flame.py b/interfaces/cython/cantera/examples/onedim/diffusion_flame.py index 3af0b2306..e402ef87a 100644 --- a/interfaces/cython/cantera/examples/onedim/diffusion_flame.py +++ b/interfaces/cython/cantera/examples/onedim/diffusion_flame.py @@ -6,6 +6,7 @@ import cantera as ct import numpy as np import csv +# Input parameters p = ct.one_atm # pressure tin_f = 300.0 # fuel inlet temperature tin_o = 300.0 # oxidizer inlet temperature @@ -15,6 +16,8 @@ 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 +# Distance between inlets is 2 cm. +# Start with an evenly-spaced 6-point grid. initial_grid = np.linspace(0, 0.02, 6) tol_ss = [1.0e-5, 1.0e-12] # [rtol, atol] for steady-state problem @@ -23,10 +26,17 @@ tol_ts = [5.0e-4, 1.0e-9] # [rtol, atol] for time stepping loglevel = 1 # amount of diagnostic output (0 to 5) refine_grid = 1 # 1 to enable refinement, 0 to disable +# Create the gas object used to evaluate all thermodynamic, kinetic, and +# transport properties. gas = ct.Solution('gri30.xml', 'gri30_mix') gas.TP = gas.T, p +# Create an object representing the counterflow flame configuration, +# which consists of a fuel inlet on the left, the flow in the middle, +# and the oxidizer inlet on the right. f = ct.CounterflowDiffusionFlame(gas, initial_grid) + +# Set the state of the two inlets f.fuel_inlet.mdot = mdot_f f.fuel_inlet.X = comp_f f.fuel_inlet.T = tin_f @@ -35,19 +45,27 @@ f.oxidizer_inlet.mdot = mdot_o f.oxidizer_inlet.X = comp_o f.oxidizer_inlet.T = tin_o +# Set error tolerances f.flame.set_steady_tolerances(default=tol_ss) f.flame.set_transient_tolerances(default=tol_ts) +# construct the initial solution estimate. To do so, it is necessary +# to specify the fuel species. If a fuel mixture is being used, +# specify a representative species here for the purpose of +# constructing an initial guess. f.set_initial_guess(fuel='C2H6') + +# First disable the energy equation and solve the problem without +# refining the grid f.energy_enabled = False f.solve(loglevel, refine_grid=False) +# Now specify grid refinement criteria, turn on the energy equation, +# and solve the problem again. f.energy_enabled = True f.set_refine_criteria(ratio=4, slope=0.2, curve=0.3, prune=0.04) f.solve(loglevel, refine_grid=refine_grid) f.show_solution() - - f.save('c2h6_diffusion.xml') z = f.flame.grid @@ -55,6 +73,7 @@ T = f.T u = f.u V = f.V +# write the velocity, temperature, and mole fractions to a CSV file with open('c2h6_diffusion.csv', 'w') as csvfile: writer = csv.writer(csvfile) writer.writerow(['z (m)', 'u (m/s)', 'V (1/s)', 'T (K)', 'rho (kg/m3)'] + @@ -64,3 +83,5 @@ with open('c2h6_diffusion.csv', 'w') as csvfile: writer.writerow([z[n], u[n], V[n], T[n], gas.density] + list(gas.X)) print('solution saved to c2h6_diffusion.csv') + +f.show_stats(0) diff --git a/interfaces/cython/cantera/examples/thermo/critical_properties.py b/interfaces/cython/cantera/examples/thermo/critical_properties.py new file mode 100644 index 000000000..7fca20b0b --- /dev/null +++ b/interfaces/cython/cantera/examples/thermo/critical_properties.py @@ -0,0 +1,27 @@ +""" +Print the critical state properties for the fluids for which Cantera has +built-in liquid/vapor equations of state. +""" + +import cantera as ct + +fluids = {'water': ct.Water(), + 'nitrogen': ct.Nitrogen(), + 'methane': ct.Methane(), + 'hydrogen': ct.Hydrogen(), + 'oxygen': ct.Oxygen(), + 'carbon dioxide': ct.CarbonDioxide(), + 'heptane': ct.Heptane(), + 'hfc134a': ct.Hfc134a() + } + +print('Critical State Properties') +print('%20s %10s %10s %10s' % ('Fluid','Tc [K]', 'Pc [Pa]', 'Zc')) +for name in fluids: + f = fluids[name] + tc = f.critical_temperature + pc = f.critical_pressure + rc = f.critical_density + mw = f.mean_molecular_weight + zc = pc * mw / (rc * ct.gas_constant * tc) + print('%20s %10.4g %10.4G %10.4G' % (name, tc, pc, zc)) diff --git a/interfaces/cython/cantera/examples/thermo/isentropic.py b/interfaces/cython/cantera/examples/thermo/isentropic.py new file mode 100644 index 000000000..9c532cf56 --- /dev/null +++ b/interfaces/cython/cantera/examples/thermo/isentropic.py @@ -0,0 +1,71 @@ +import cantera as ct +import math +import numpy as np + + +def soundspeed(gas): + """The speed of sound. Assumes an ideal gas.""" + + gamma = gas.cp / gas.cv + return math.sqrt(gamma * ct.gas_constant + * gas.T / gas.mean_molecular_weight) + + +def isentropic(gas=None): + """ + ISENTROPIC isentropic, adiabatic flow example + + In this example, the area ratio vs. Mach number curve is computed. If a + gas object is supplied, it will be used for the calculations, with the + stagnation state given by the input gas state. Otherwise, the calculations + will be done for a 10:1 hydrogen/nitrogen mixture with stagnation T0 = + 1200 K, P0 = 10 atm. + + """ + + if gas is None: + gas = ct.Solution('gri30.xml') + gas.TPX = 1200.0, 10.0*ct.one_atm, 'H2:1,N2:0.1' + + # get the stagnation state parameters + s0 = gas.s + h0 = gas.h + p0 = gas.P + + mdot = 1 # arbitrary + amin = 1.e14 + + data = np.zeros((200,4)) + + # compute values for a range of pressure ratios + for r in range(200): + + p = p0*(r+1)/201.0 + # set the state using (p,s0) + gas.SP = s0, p + + v2 = 2.0*(h0 - gas.h) # h + V^2/2 = h0 + 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[:,0] /= amin + + return data + + +if __name__ == "__main__": + print(isentropic.__doc__) + data = isentropic() + try: + import matplotlib.pyplot as plt + plt.plot(data[:,1], data[:,0]) + plt.ylabel('Area Ratio') + plt.xlabel('Mach Number') + plt.title('Isentropic Flow: Area Ratio vs. Mach Number') + plt.show() + + except ImportError: + print('area ratio, Mach number, temperature, pressure ratio') + print(data) diff --git a/interfaces/cython/cantera/examples/thermo/rankine.py b/interfaces/cython/cantera/examples/thermo/rankine.py new file mode 100644 index 000000000..cafb8d06a --- /dev/null +++ b/interfaces/cython/cantera/examples/thermo/rankine.py @@ -0,0 +1,75 @@ +""" +A Rankine vapor power cycle +""" + +import cantera as ct + +# parameters +eta_pump = 0.6 # pump isentropic efficiency +eta_turbine = 0.8 # turbine isentropic efficiency +p_max = 8.0e5 # maximum pressure + + +def pump(fluid, p_final, eta): + """Adiabatically pump a fluid to pressure p_final, using + a pump with isentropic efficiency eta.""" + h0 = fluid.h + s0 = fluid.s + fluid.SP = s0, p_final + h1s = fluid.h + isentropic_work = h1s - h0 + actual_work = isentropic_work / eta + h1 = h0 + actual_work + fluid.HP = h1, p_final + return actual_work + + +def expand(fluid, p_final, eta): + """Adiabatically expand a fluid to pressure p_final, using + a turbine with isentropic efficiency eta.""" + h0 = fluid.h + s0 = fluid.s + fluid.SP =s0, p_final + h1s = fluid.h + isentropic_work = h0 - h1s + actual_work = isentropic_work * eta + h1 = h0 - actual_work + fluid.HP = h1, p_final + return actual_work + + +def printState(n, fluid): + print('\n***************** State {} ******************'.format(n)) + print(fluid.report()) + + +if __name__ == '__main__': + # create an object representing water + w = ct.Water() + + # start with saturated liquid water at 300 K + w.TX = 300.0, 0.0 + h1 = w.h + p1 = w.P + printState(1, w) + + # pump it adiabatically to p_max + pump_work = pump(w, p_max, eta_pump) + h2 = w.h + printState(2, w) + + # heat it at constant pressure until it reaches the saturated vapor state + # at this pressure + w.PX = p_max, 1.0 + h3 = w.h + heat_added = h3 - h2 + printState(3, w) + + # expand back to p1 + turbine_work = expand(w, p1, eta_turbine) + printState(4, w) + + # efficiency + eff = (turbine_work - pump_work)/heat_added + + print('efficiency = ', eff) diff --git a/interfaces/cython/cantera/examples/thermo/sound_speed.py b/interfaces/cython/cantera/examples/thermo/sound_speed.py new file mode 100644 index 000000000..938d0ba31 --- /dev/null +++ b/interfaces/cython/cantera/examples/thermo/sound_speed.py @@ -0,0 +1,52 @@ +import cantera as ct +import math + + +def equilSoundSpeeds(gas, rtol=1.0e-6, maxiter=5000): + """ + Returns a tuple containing the equilibrium and frozen sound speeds for a + gas with an equilibrium composition. The gas is first set to an + equilibrium state at the temperature and pressure of the gas, since + otherwise the equilibrium sound speed is not defined. + """ + + # set the gas to equilibrium at its current T and P + gas.equilibrate('TP', rtol=rtol, maxiter=maxiter) + + # save properties + s0 = gas.s + p0 = gas.P + r0 = gas.density + + # perturb the pressure + p1 = p0*1.0001 + + # set the gas to a state with the same entropy and composition but + # the perturbed pressure + gas.SP = s0, p1 + + # frozen sound speed + afrozen = math.sqrt((p1 - p0)/(gas.density - r0)) + + # now equilibrate the gas holding S and P constant + gas.equilibrate('SP', rtol=rtol, maxiter=maxiter) + + # equilibrium sound speed + aequil = math.sqrt((p1 - p0)/(gas.density - r0)) + + # compute the frozen sound speed using the ideal gas expression as a check + gamma = gas.cp/gas.cv + afrozen2 = math.sqrt(gamma * ct.gas_constant * gas.T / + gas.mean_molecular_weight) + + return aequil, afrozen, afrozen2 + + +# test program +if __name__ == "__main__": + gas = ct.Solution('gri30.xml') + gas.X = 'CH4:1.00, O2:2.0, N2:7.52' + for n in range(27): + T = 300.0 + 100.0 * n + gas.TP = T, ct.one_atm + print(T, equilSoundSpeeds(gas))