so I am new to both python and cantera but I've found my way around and am able to get my simulation to run how I need it to. The only trouble I am having is being able to save the data that cantera computes as a csv file. I know I need to use the write_csv command, and i did that and was able to create the csv, but the data is not populating and I'm just a little stumped.
import pandas as pd
import numpy as np
import time
import cantera as ct
import matplotlib.pyplot as plt
import matplotlib as mpl
import csv
mpl.style.use('classic')
mpl.rcParams['lines.linewidth'] = 2
mpl.rcParams['patch.force_edgecolor'] = True
mpl.rcParams['patch.facecolor'] = 'k'
mpl.rcParams['legend.edgecolor'] = 'inherit'
mpl.rcParams['legend.fancybox'] = True
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['xtick.labelsize'] = 14
plt.rcParams['ytick.labelsize'] = 14
mpl.rcParams['legend.fontsize'] = 'large'
plt.rcParams['figure.autolayout'] = True
gas = ct.Solution('mech_2020_ammonia_fixed.yaml')
# Specify initial conditions ######################
reactorTemperature = 900 # Kelvin
reactorPressure = 24115350 # in Pa.
equiv_ratio =0.714
###################################################
gas.TP = reactorTemperature, reactorPressure
#set(gas,'T',913.15,'P',reactorPressure,'Y','nh3:0.03, h2o:0.95', 'o2:0.21, n2:0.79');
fuel={"NH3":0.8, "C2H5OH":0.0, "H2O":0.1}
oxidizer={"O2":0.21, "N2":0.79}
gas.set_equivalence_ratio(equiv_ratio, fuel, oxidizer, basis='mass')
phi_mass = gas.equivalence_ratio(fuel, oxidizer, basis='mass')
phi2_mole = gas.equivalence_ratio(fuel, oxidizer, basis='mole')
rho0 = gas.density # initial density of the flow
gas()
# In this example, the result is the same as above
print("phi_mass = {:1.3f}".format(phi_mass))
print("phi_mole = {:1.3f}".format(phi2_mole))
print("mass fraction of ammonia = {:1.3f}".format(gas["NH3"].Y[0]))
print("mass fraction of oxygen = {:1.3f}".format(gas["O2"].Y[0]))
print("mass fraction of nitrogen= {:1.3f}".format(gas["N2"].Y[0]))
print("mass fraction of water = {:1.3f}".format(gas["H2O"].Y[0]))
print("mass fraction of Ethanol = {:1.3f}".format(gas["C2H5OH"].Y[0]))
# Reactor parameters
residenceTime = 300 # s
reactorVolume = 2*(1e-5) # m3
# Instrument parameters
# This is the "conductance" of the pressure valve and will determine its efficiency in
# holding the reactor pressure to the desired conditions.
pressureValveCoefficient =0.0001
# This parameter will allow you to decide if the valve's conductance is acceptable. If there
# is a pressure rise in the reactor beyond this tolerance, you will get a warning
maxPressureRiseAllowed = 1000
fuelAirMixtureTank = ct.Reservoir(gas)
exhaust = ct.Reservoir(gas)
stirredReactor = ct.IdealGasReactor(gas, energy='off', volume=reactorVolume, name='isothermal_reactor')
massFlowController = ct.MassFlowController(upstream=fuelAirMixtureTank,
downstream=stirredReactor,
mdot=stirredReactor.mass/residenceTime)
mdot=stirredReactor.mass/residenceTime
## Verify Mass, flow rate, and volume
print("Volume(m3)",reactorVolume)
print("Mass Flow Rate (kg/s)",mdot)
print("Mass in Reactor (kg)",stirredReactor.mass)
print("Residence Time (s)",residenceTime)
## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## vv
pressureRegulator = ct.Valve(upstream=stirredReactor,
downstream=exhaust,
K=pressureValveCoefficient)
reactorNetwork = ct.ReactorNet([stirredReactor])
# Now compile a list of all variables for which we will store data
columnNames = [stirredReactor.component_name(item) for item in range(stirredReactor.n_vars)]
columnNames = ['pressure'] + columnNames
# Use the above list to create a DataFrame
timeHistory = pd.DataFrame(columns=columnNames)
# Start the stopwatch
maxSimulationTime = 1500 # seconds
tic = time.time()
# Set simulation start time to zero
t = 0
counter = 1
while t < maxSimulationTime:
t = reactorNetwork.step()
# We will store only every 10th value. Remember, we have 1200+ species, so there will be
# 1200 columns for us to work with
if(counter%10 == 0):
# Extract the state of the reactor
state = np.hstack([stirredReactor.thermo.P, stirredReactor.mass,
stirredReactor.volume, stirredReactor.T, stirredReactor.thermo.X])
# Update the dataframe
timeHistory.loc[t] = state
counter += 1
# Stop the stopwatch
toc = time.time()
print('Simulation Took {:3.2f}s to compute, with {} steps'.format(toc-tic, counter))
# We now check to see if the pressure rise during the simulation, a.k.a the pressure valve
# was okay
pressureDifferential = timeHistory['pressure'].max()-timeHistory['pressure'].min()
if(abs(pressureDifferential/reactorPressure) > maxPressureRiseAllowed):
print("WARNING: Non-trivial pressure rise in the reactor. Adjust K value in valve")
print("reactorpressure",reactorPressure)
plt.figure(1)
plt.semilogx(timeHistory.index, timeHistory['N2'],'-o',color='C2', label='N2')
plt.semilogx(timeHistory.index, timeHistory['NH3'],'-o',color='C0', label='NH3')
plt.semilogx(timeHistory.index, timeHistory['O2'],'-o',color='C6', label='O2')
plt.semilogx(timeHistory.index, timeHistory['C2H5OH'],'-o',color='C3', label='C2H5OH')
plt.semilogx(timeHistory.index, timeHistory['H2O'],'-o',color='C1', label='H2O')
plt.xlabel('Time (s)')
plt.ylabel(r'Molar Concentration');
plt.grid(b=None)
plt.legend(fontsize='large',loc = 'best')
plt.tight_layout()
plt.xlim(0, 10**4)
plt.savefig('C:/ammonia_Trials/Temp_vary/Runs/concentration_6.3k.png',orientation='landscape', bbox_inches='tight')
plt.show()
plt.figure(2)
plt.semilogx(timeHistory.index, timeHistory['H2O'],'-o',color='C3', label='H$_{2}$o')
plt.xlabel('Time (s)')
plt.ylabel(r'Mass Fraction for H$_{2}$O');
plt.grid(b=None)
plt.tight_layout()
plt.xlim(0, 10**4)
plt.savefig('C:/ammonia_Trials/Temp_vary/Runs/water_concentraction_6.3H.png', orientation='landscape', bbox_inches='tight')
plt.show()
plt.figure(3)
plt.semilogx(timeHistory.index, timeHistory['temperature'],'-o',color='C3', label='Temperature')
plt.xlabel('Time (s)')
plt.ylabel(r'Temperature (K)');
plt.grid(b=None)
plt.tight_layout()
plt.xlim(0, 10**4)
plt.savefig('C:/ammonia_Trials/Temp_vary/Runs/temperature_6.3H.png',orientation='landscape', bbox_inches='tight')
plt.show()
plt.figure(4)
plt.semilogx(timeHistory.index, timeHistory['pressure'],'-o',color='C8', label='Pressure')
plt.ylim(0,30000000)
plt.xlabel('Time (s)')
plt.ylabel(r'Pressure (Pa)');
plt.grid(b=None)
plt.tight_layout()
plt.xlim(0, 10**4)
plt.savefig('C:/ammonia_trials/Temp_vary/Runs/pressure_6.3H.png',orientation='landscape', bbox_inches='tight')
plt.show()
plt.figure(5)
plt.semilogx(timeHistory.index, timeHistory['mass']*1000,'-o',color='C8', label='Mass')
plt.xlabel('Time (s)')
plt.ylabel(r'Total Mass Fraction in Reactor (g)');
plt.grid(b=None)
plt.tight_layout()
plt.xlim(0, 10**4)
plt.savefig('C:/ammonia_trials/Temp_vary/Runs/mass_6.3H.png', orientation='landscape', bbox_inches='tight')
plt.show()
plt.figure(6)
plt.semilogx(timeHistory.index, timeHistory['CO'],'-o',color='C6', label='CO')
plt.semilogx(timeHistory.index, timeHistory['CO2'],'-o',color='C5', label='CO$_{2}$')
plt.xlabel('Time (s)')
plt.ylabel(r'Species Mass Fraction');
plt.grid(b=None)
plt.tight_layout()
plt.legend(fontsize='large',loc = 'best')
plt.xlim(0, 10**4)
plt.savefig('C:/ammonia_trials/Temp_vary/Runs/co_co2_6.3H.png', orientation='landscape', bbox_inches='tight')
plt.show()
plt.figure(7)
plt.semilogx(timeHistory.index, timeHistory['CH2O'],'-o',color='C6', label='CH2O')
plt.xlabel('Time (s)')
plt.ylabel(r'Species Mass Fraction for CH$_{2}$O');
plt.grid(b=None)
plt.tight_layout()
plt.xlim(0, 10**4)
plt.savefig('C:/ammonia_trials/temp_vary/runs/formaldehyde_6.3H.png', orientation='landscape', bbox_inches='tight')
plt.show()
#writes all mass fractions with respect to time in the CSV file below
#open the file in write mode
columnNames = [stirredReactor.component_name(item) for item in range(stirredReactor.n_vars)]
columnNames = ['mf'] + columnNames
timeHistory = pd.DataFrame(columns=columnNames)
tic = time.time()
timeHistory.write_csv('7_3a_repeat.csv', cols=('t','T','Y','P'))