############################################################### # # Autoignition of a methane air mixture at stoichiometry # and atmospheric pressure, for different # initial temperature # ############################################################### import sys,os import numpy as np from cantera import * ################################################################# #Mechanisms used for the process matrice_cas = ['gri30'] n_cas = len(matrice_cas) #Specify the number of time steps and the time step size nt = 20000 dt = 1.e-7 #s ################# #Enter general parameters phi = 1.5 P = 101325.0 # Pascal cond = 'UV' # Const P or Cons V fuel_species = 'CH4' stoich_O2 = 2 air_N2_O2_ratio = 3.76 Tmin = 3.5 # Kelvin Tmax = 10. # Kelvin npoints = 26 Ti = np.zeros(npoints,'d') Ti2 = np.zeros(npoints,'d') tim = np.zeros(nt,'d') temp_cas = np.zeros(nt,'d') dtemp_cas = np.zeros(nt-1,'d') T_check = np.zeros(nt-1,'d') Autoignition_cas = [] FinalTemp_cas = [] Y_CO = np.zeros([npoints,nt],'d') Y_OH = np.zeros([npoints,nt],'d') ################# #Loop over case for i, cas in enumerate(matrice_cas): print 'Cas ' + cas +' :\n' gas = Solution(cas + '.cti') # find index m = gas.n_species ifuel = gas.species_index(fuel_species) io2 = gas.species_index('O2') in2 = gas.species_index('N2') ico = gas.species_index('CO') ioh = gas.species_index('OH') # create some arrays to hold the data FinalTemp_cas.append(np.zeros(npoints,'d')) Autoignition_cas.append(np.zeros(npoints,'d')) #Loop over initial conditions for j in range(npoints): Ti2[j] = Tmin + (Tmax - Tmin)*j/(npoints - 1) Ti[j] = 10000/Ti2[j] x = np.zeros(m,'d') x[ifuel] = phi x[io2] = stoich_O2 x[in2] = stoich_O2*air_N2_O2_ratio gas.TPX = Ti[j], P, x # Now create a reactor network consisting of the single batch reactor r = IdealGasReactor(gas) sim = ReactorNet([r]) #Run the simulation # Initial simulation time time = 0.0 if (Ti[j]<1480.0): nt = 20000 dt = 1.e-6 if (Ti[j]<1205.0): nt = 20000 dt = 1.e-5 if (Ti[j]<1055.0): nt = 20000 dt = 1.e-4 if (Ti[j]<950.0): nt = 50000 dt = 5.e-4 #Loop for nt time steps of dt seconds. for n in range(nt): time += dt sim.advance(time) tim[n] = time temp_cas[n] = r.T #for sp in range(m): # Y_CO[j][n] = r.Y[ico] # Y_OH[j][n] = r.Y[ioh] # print '%10.3e %10.3f' % (tim[n], r.T) # r.pressure(), r.intEnergy_mass()) ################################################################# # Catch the autoignition timings ################################################################# # Get autoignition timing Dtmax = [0,0.0] for n in range(nt-1): dtemp_cas[n] = (temp_cas[n+1]-temp_cas[n])/dt T_check[n] = (temp_cas[n]+temp_cas[n + 1])/2.0 if ((dtemp_cas[n]>Dtmax[1]) and (T_check[n]>(Ti[j] + 3.0*Ti[j]/100))): Dtmax[0] = n Dtmax[1] = dtemp_cas[n] # Local print Autoignition = (tim[Dtmax[0]]+tim[Dtmax[0] + 1])/2. print 'For ' +str(Ti[j]) +' ' + str(Ti2[j]) +', Autoig time = ' + str(Autoignition), 'Autoig Temp = ' + str((temp_cas[Dtmax[0]]+temp_cas[Dtmax[0] + 1])/2.0) # Posterity Autoignition_cas[i][j] = Autoignition*1000 #ms FinalTemp_cas[i][j] = temp_cas[nt-1] csv_file = 'T_'+str(Ti2[j])+'.csv' with open(csv_file, 'w') as outfile: writer = csv.writer(outfile) #writer.writerow(['time','Temperature','CO', 'OH']) writer.writerow(['time','Temperature']) time = 0.0 for jj in range(nt): #writer.writerow([time, temp_cas[jj], Y_CO[j][jj], Y_OH[j][jj]]) writer.writerow([time, temp_cas[jj]]) time = time + dt #if (Ti2[j] == 0.45): ##if ((Ti2[j] == 0.45) or (Ti2[j] == 0.6) or (Ti2[j] == 0.75) or ( Ti2[j] == 1.1) or (Ti2[j] == 1.25)): # element = 'C' # print "Je recommence l analyse pour rxn" # gas.TPX = Ti[j], P, x # r = IdealGasReactor(gas) # sim = ReactorNet([r]) # # Initial simulation time # time = 0.0 # for n in range(nt): # time += dt # sim.advance(time) # if (n==Dtmax[0]): # print "Time ?", time # print "Temp ?", r.T, gas.T # dutemps = r.T # diagram = ReactionPathDiagram(gas, element) # diagram.title = 'Reaction path diagram following {0}'.format(element) # diagram.label_threshold = 0.05 # diagram.threshold = 0.05 # diagram.scale = -1 # print(diagram.threshold) # diagram.show_details = False # print(diagram.show_details) # dot_file = 'Spec28_phi'+str(phi)+'_Tini'+str(Ti2[j])+'_Tig'+str(dutemps)+'.dot' # diagram.write_dot(dot_file) # print(diagram.get_data()) # # print("Wrote graphviz input file to '{0}'.".format(os.path.join(os.getcwd(), dot_file))) # if (n==(Dtmax[0]-1)): # print "Time ?", time # print "Temp ?", r.T, gas.T # dutemps = r.T # diagram = ReactionPathDiagram(gas, element) # diagram.title = 'Reaction path diagram following {0}'.format(element) # diagram.label_threshold = 0.05 # diagram.threshold = 0.05 # diagram.scale = -1 # print(diagram.threshold) # diagram.show_details = False # print(diagram.show_details) # dot_file = 'Spec28_phi'+str(phi)+'_Tini'+str(Ti2[j])+'_Tig'+str(dutemps)+'.dot' # diagram.write_dot(dot_file) # print(diagram.get_data()) # # print("Wrote graphviz input file to '{0}'.".format(os.path.join(os.getcwd(), dot_file))) ################################################################# # Save results ################################################################# # write output CSV file for importing into Excel csv_file = 'Phi'+ str(phi) +'_P'+str(P)+'_Trange_'+str(cond)+'.csv' with open(csv_file, 'w') as outfile: writer = csv.writer(outfile,delimiter=' ',quotechar='|', quoting=csv.QUOTE_MINIMAL) writer.writerow(['###############']) writer.writerow(['#IGOR']) writer.writerow(['#WAVES','1/T','Tig','Tend','Tini']) writer.writerow(['#BEGIN']) for i, cas in enumerate(matrice_cas): for j in range(npoints): writer.writerow([Ti2[j], Autoignition_cas[i][j], FinalTemp_cas[i][j], Ti[j]]) writer.writerow(['#END']) print 'output written to '+csv_file