[Kinetics] Added LMRR cantera examples

doc/example-keywords.txt

@@ -18,6 +18,7 @@ fuel cell
 heat transfer
 ignition delay
 internal combustion engine
+jet-stirred reactor
 kinetics
 Matlab
 mixture
@@ -41,6 +42,7 @@ saving output
 sensitivity analysis
 strained flame
 surface chemistry
+shock tube
 thermodynamic cycle
 thermodynamics
 transport

samples/python/kinetics/jet_stirred_reactor.py

@@ -0,0 +1,107 @@
+"""
+Jet-stirred reactor temperature and species profiles
+====================================================
+
+Simulate temperature profiles and species profiles in a jet-stirred reactor across a range of initial temperatures, and observe the impact of incorporating the reduced-pressure linear mixture rule (LMR-R) in such calculations.
+
+Here we will consider a mixture of H2/O2/NH3/Ar (with 10% NH3) at 1.2 atm, and compare results against the experimental measurements of Sabia et al. [1] Two models are compared in this example:
+
+1. A 2023 model of H2 and NH3 chemistry published by Alzueta et al. [2]
+2. An adapted version of this model that has applied the reduced-pressure linear mixture rule (LMR-R) and ab initio third-body efficiencies. [3]
+
+References:
+
+    [1] P. Sabia, M. V. Manna, R. Ragucci, M. de Joannon, Mutual inhibition effect of hydrogen and ammonia in oxidation processes and the role of ammonia as “strong” collider in third-molecular reactions, Int. J. Hydrogen Energy 45 (2020) 32113 -- 32127.
+
+    [2] M. U. Alzueta, I. Salas, H. Hashemi, P. Glarborg, CO-assisted NH3 oxidation, Combust. Flame 257 (2023) 112438.
+
+    [3] P. J. Singal, J. Lee, L. Lei, R. L. Speth, M. P. Burke, Implementation of New Mixture Rules Has a Substantial Impact on Combustion Predictions for H2 and NH3, Proc. Combust. Inst. 41 (2025).
+
+Requires: cantera >= 3.1
+.. tags:: jet-stirred reactor, kinetics
+"""
+
+import numpy as np
+import pandas as pd
+import time as time
+import cantera as ct
+import matplotlib.pyplot as plt
+
+f, ax = plt.subplots(1, 3)
+plt.subplots_adjust(wspace=0.6)
+colours = ["xkcd:grey",'xkcd:purple']
+file = 'example_data/ammonia-CO-H2-Alzueta-2023.yaml'
+models = {'Original': 'baseline', 'LMR-R': 'linear-Burke'}
+
+inputs = {
+    'X': {'H2': 0.03, 'O2': 0.03, 'Ar': 0.846, 'NH3':0.094},
+    'T_range': np.linspace(800,1050,50), # [K]
+    'Tin': 1000, # reactor temperature [K]
+    'P': 1.2, # reactor pressure [atm]
+    'tau': 0.5, # residence time [s]
+    'V': 0.000113, # reactor volume [m3]
+    'K': 2e-5, # 'pressureValveCoefficient'
+    't_max': 50,  # [s]
+    'h': 79.5, # 'heatTransferCoefficient' [W/m2/K]
+    'data': { # experimental data from Sabia et al.
+        'T_range': [807,843,855,870,884,904,925,945,965,995,1018],
+        'deltaT': [0.051,0.051,0.051,0.051,0.101,0.606,1.414,2.626,4.091,6.768,8.586],
+        'X_O2': [3.076,3.053,3.050,3.037,3.024,3.015,2.966,2.924,2.794,2.597,2.261],
+        'X_H2': [3.030,3.038,3.038,3.038,3.030,2.993,2.948,2.829,2.693,2.434,2.126]
+    }
+}
+
+def getStirredReactor(gas,inputs):
+    reactorRadius = (inputs['V']*3/4/np.pi)**(1/3) # [m3]
+    reactorSurfaceArea =4*np.pi*reactorRadius**2 # [m3]
+    fuelAirMixtureTank = ct.Reservoir(gas)
+    exhaust = ct.Reservoir(gas)
+    env = ct.Reservoir(gas)
+    reactor = ct.IdealGasReactor(gas, energy='on', volume=inputs['V'])
+    ct.MassFlowController(upstream=fuelAirMixtureTank,
+                          downstream=reactor,
+                          mdot=reactor.mass/inputs['tau'])
+    ct.Valve(upstream=reactor,
+             downstream=exhaust,
+             K=inputs['K'])
+    ct.Wall(reactor, env, A=reactorSurfaceArea, U=inputs['h'])
+    return reactor
+
+def getTemperatureDependence(gas, inputs):
+    stirredReactor = getStirredReactor(gas,inputs)
+    columnNames = ['pressure'] + [stirredReactor.component_name(item) for item in range(stirredReactor.n_vars)]
+    tempDependence = pd.DataFrame(columns=columnNames)
+    for T in inputs['T_range']:
+        gas.TPX = T, inputs['P']*ct.one_atm, inputs['X']
+        stirredReactor = getStirredReactor(gas,inputs)
+        reactorNetwork = ct.ReactorNet([stirredReactor])
+        t = 0
+        while t < inputs['t_max']:
+            t = reactorNetwork.step()
+        state = np.hstack([stirredReactor.thermo.P,
+                        stirredReactor.mass,
+                        stirredReactor.volume,
+                        stirredReactor.T,
+                        stirredReactor.thermo.X])
+        tempDependence.loc[T] = state
+    return tempDependence
+
+for k,m in enumerate(models):
+    gas = ct.Solution(file, name=models[m])
+    gas.TPX = inputs['Tin'], inputs['P']*ct.one_atm, inputs['X']
+    tempDependence = getTemperatureDependence(gas,inputs)
+    ax[0].plot(tempDependence.index, np.subtract(tempDependence['temperature'],tempDependence.index), color=colours[k],label=m)
+    ax[1].plot(tempDependence.index, tempDependence['O2']*100, color=colours[k])
+    ax[2].plot(tempDependence.index, tempDependence['H2']*100, color=colours[k])
+ax[0].plot(inputs['data']['T_range'],inputs['data']['deltaT'],'o',fillstyle='none',color='k',label="Sabia et al.")
+ax[1].plot(inputs['data']['T_range'],inputs['data']['X_O2'],'o',fillstyle='none',color='k')
+ax[2].plot(inputs['data']['T_range'],inputs['data']['X_H2'],'o',fillstyle='none',color='k')
+ax[0].legend(fontsize=8,frameon=False,loc='upper left')
+ax[0].set_ylabel(r'$\Delta$ T [K]')
+ax[1].set_xlabel(r'Temperature [K]')
+ax[1].set_ylabel(r'O$_2$ mole fraction [%]')
+ax[2].set_ylabel(r'H$_2$ mole fraction [%]')
+ax[0].set_xlim([780,1070])
+ax[1].set_xlim([780,1070])
+ax[2].set_xlim([780,1070])
+plt.show()

samples/python/kinetics/shock_tube.py

@@ -0,0 +1,62 @@
+"""
+Shock-tube species profiles as a function of time
+=================================================
+
+Simulate species profiles for a shock tube as a function of time, and observe the impact of incorporating the reduced-pressure linear mixture rule (LMR-R) in such calculations.
+
+Here we predict the H2O mole fraction time profiles for a mixture of 1163 ppm H2O2/1330 ppm H2O/665 ppm O2/20% CO2/Ar following reflected shock waves (1196 K, 2.127 atm) and compare results against the experimental measurements of Shao et al. [1] Two models are compared in this example:
+
+1. A 2023 model of H2 and NH3 chemistry published by Alzueta et al. [2]
+2. An adapted version of this model that has applied the reduced-pressure linear mixture rule (LMR-R) and ab initio third-body efficiencies. [3]
+
+References:
+
+    [1] J. Shao, R. Choudhary, D. F. Davidson, R. K. Hanson, Shock tube/laser absorption measurement of the rate constant of the reaction: H2O2+CO2 = 2OH+CO2, Proc. Combust. Inst. 39 (2023) 735 -- 743.
+
+    [2] M. U. Alzueta, I. Salas, H. Hashemi, P. Glarborg, CO-assisted NH3 oxidation, Combust. Flame 257 (2023) 112438.
+
+    [3] P. J. Singal, J. Lee, L. Lei, R. L. Speth, M. P. Burke, Implementation of New Mixture Rules Has a Substantial Impact on Combustion Predictions for H2 and NH3, Proc. Combust. Inst. 41 (2025).
+
+Requires: cantera >= 3.1
+.. tags:: shock tube, kinetics
+"""
+
+import cantera as ct
+import matplotlib.pyplot as plt
+import numpy as np
+
+fig, ax = plt.subplots()
+file = 'example_data/ammonia-CO-H2-Alzueta-2023.yaml'
+models = {'Original': 'baseline', 'LMR-R': 'linear-Burke'}
+colours = ["xkcd:grey",'xkcd:purple']
+
+for k,m in enumerate(models):
+    X_H2O2 = 1163e-6
+    X_H2O = 1330e-6
+    X_O2 = 665e-6
+    X_CO2= 0.2*(1-X_H2O2-X_H2O-X_O2)
+    X_Ar = 1-X_CO2
+    gas = ct.Solution(file, name=models[m])
+    gas.TPX = 1196, 2.127*ct.one_atm, {'H2O2':X_H2O2, 'H2O':X_H2O, 'O2':X_O2, 'CO2':X_CO2, 'AR':X_Ar}
+    r = ct.Reactor(contents=gas,energy="on")
+    reactorNetwork = ct.ReactorNet([r])
+    timeHistory = ct.SolutionArray(gas, extra=['t'])
+    estIgnitDelay = 1
+    t = 0
+    counter = 1
+    while t < estIgnitDelay:
+        t = reactorNetwork.step()
+        if counter % 10 == 0:
+            timeHistory.append(r.thermo.state, t=t)
+        counter += 1
+    ax.plot(timeHistory.t*1e6, timeHistory('H2O').X*100, color=colours[k],label=m)
+expData = {
+    't': [12.3,20.3,26.4,39.6,58.5,79.2,96.1,113.8,131.6,145.7,161.2,181.6,195.3,219.9,237.2,248.6,262.4,272.2,280.9],
+    'X_H2O': [1.47E-03,1.59E-03,1.66E-03,1.78E-03,1.98E-03,2.06E-03,2.15E-03,2.22E-03,2.26E-03,2.30E-03,2.39E-03,2.38E-03,2.40E-03,2.42E-03,2.47E-03,2.53E-03,2.51E-03,2.50E-03,2.47E-03]
+}
+ax.plot(expData['t'],np.array(expData['X_H2O'])*100,'o',fillstyle='none',color='k',label='Shao et al.')
+ax.legend(frameon=False, loc='lower right')
+ax.set_ylabel(r'$\rm H_2O$ mole fraction [%]')
+ax.set_xlabel(r'Time [$\mathdefault{\mu s}$]')
+ax.set_xlim([0,300])
+plt.show()

samples/python/onedim/flame_speed.py

@@ -0,0 +1,60 @@
+"""
+Flame speed as a function of equivalence ratio
+==============================================
+
+Simulate flame speeds across a range of equivalence ratios, and observe the impact of incorporating the reduced-pressure linear mixture rule (LMR-R) in such calculations.
+
+Here we will consider a mixture of NH3/air (1 atm, 296 K) and compare results against the experimental measurements of Ronney. [1] Two models are compared in this example:
+
+1. A 2023 model of H2 and NH3 chemistry published by Alzueta et al. [2]
+2. An adapted version of this model that has applied the reduced-pressure linear mixture rule (LMR-R) and ab initio third-body efficiencies. [3]
+
+References:
+
+    [1] P. D. Ronney, Effect of chemistry and transport properties on near-limit flames at microgravity, Combust. Sci. Tech. 59 (1988) 123 -- 141.
+
+    [2] M. U. Alzueta, I. Salas, H. Hashemi, P. Glarborg, CO-assisted NH3 oxidation, Combust. Flame 257 (2023) 112438.
+
+    [3] P. J. Singal, J. Lee, L. Lei, R. L. Speth, M. P. Burke, Implementation of New Mixture Rules Has a Substantial Impact on Combustion Predictions for H2 and NH3, Proc. Combust. Inst. 41 (2025).
+
+Requires: cantera >= 3.1
+.. tags:: flame speed, kinetics
+"""
+
+import cantera as ct
+import matplotlib.pyplot as plt
+import numpy as np
+
+fig, ax = plt.subplots()
+file = 'example_data/ammonia-CO-H2-Alzueta-2023.yaml'
+models = {'Original': 'baseline', 'LMR-R': 'linear-Burke'}
+colours = ["xkcd:grey",'xkcd:purple']
+Tin = 296  # unburned gas temperature [K]
+p=760  # pressure [torr]
+n=16 # number of points to simulate
+phi_list = np.linspace(0.6,2.0,n) # equivalence ratios to simulate across
+for k, m in enumerate(models):
+    vel_list = []
+    gas = ct.Solution(file, name=models[m])
+    for j, phi in enumerate(phi_list):
+        gas.set_equivalence_ratio(phi, 'NH3', {'O2':1, 'N2': 3.76})
+        gas.TP = Tin, (p/760)*ct.one_atm
+        f = ct.FreeFlame(gas, width=0.03)
+        f.set_refine_criteria(ratio=3, slope=0.06, curve=0.10)
+        # f.transport_model = 'multicomponent' # optionally enable
+        # f.soret_enabled = True  # optionally enable
+        f.solve(loglevel=1, auto=True)
+        vel_list.append(f.velocity[0] * 100) # cm/s
+    ax.plot(phi_list, vel_list, color=colours[k],label=m)
+expData = {
+   'X_NH3': [16.3,16.4,17.0,18.0,19.0,20.0,21.9,24.0,26.0,28.5,29.0,30.0,31.0,31.5],
+   'vel': [1.35,1.48,2.30,3.36,4.01,5.88,6.80,8.14,6.73,5.00,4.78,3.3,2.9,3.0]
+}
+X_NH3 = np.divide(expData['X_NH3'],100)
+X_O2 = np.multiply(np.subtract(1,X_NH3), 0.21)
+phi_data = np.divide(np.divide(X_NH3,X_O2),np.divide(4,3))
+ax.plot(phi_data,expData['vel'],'o',fillstyle='none',color='k',label='Ronney')
+ax.legend(frameon=False, loc='upper right')
+ax.set_ylabel(r'Burning velocity [cm $\rm s^{-1}$]')
+ax.set_xlabel(r'Equivalence Ratio')
+plt.show()