############################################################### # # ADIABATIC_FLAME - A freely-propagating, premixed flat flame # For cantera V1.7 # ############################################################### #import : from cantera import * from matplotlib.pylab import * import numpy #Functions : ################################################################# # Prepare your run ################################################################# #Mechanisms used for the process matrice_cas = ['Peters','2STEP_mech','1S_C3H8_CM1'] #matrice_cas = ['Spec_30_Reac_300'] n_cas = len(matrice_cas) #Plot symbols type_plot = ['-','--','+'] #Parameter values : #General p = 1e5 # pressure tin = 300.0 # unburned gas temperature phi_min = 0.6 phi_max = 1.8 npoints = 13 phi = zeros(npoints,'d') sl_cas = [] t_cas = [] fuel_species = 'C3H8' #stoich_O2 = gas.n_atoms(fuel_species,'C') + 0.25*gas.n_atoms(fuel_species,'H') air_N2_O2_ratio = 3.76 #Initial grids, chosen to be 0.02cm long : # - Refined grid at inlet and outlet, 6 points in x-direction : initial_grid = 2*array([0.0, 0.001, 0.01, 0.02, 0.029, 0.03],'d')/3 # m # - Uniform grid, 6 points in x-direction (import numpy): #initial_grid = 0.02*array([0.0, 0.2, 0.4, 0.6, 0.8, 1.0],'d') # m # - Uniform grid of 300 points using numpy : #initial_grid = numpy.linspace(0,0.02 , 300) #Tolerance properties tol_ss = [1.0e-5, 1.0e-9] # [rtol atol] for steady-state problem tol_ts = [1.0e-5, 1.0e-9] # [rtol atol] for time stepping loglevel = 1 # amount of diagnostic output (0 # to 5) refine_grid = True # True to enable refinement, False to # disable #Import gas phases with mixture transport model for i, cas in enumerate(matrice_cas): print 'Cas ' + cas +' :\n' #Import gas phases with mixture transport model gas = Solution(cas + '.cti') stoich_O2 = gas.n_atoms(fuel_species,'C') + 0.25*gas.n_atoms(fuel_species,'H') ################# #Stoechiometry : m=gas.n_species ifuel = gas.species_index(fuel_species) io2 = gas.species_index('O2') in2 = gas.species_index('N2') sl_cas.append(zeros(npoints,'d')) t_cas.append(zeros(npoints,'d')) ################################################################# #FIRST SIMULATION : ################################################################# #Stoechiometry : phi[0] = phi_min x = zeros(m,'d') x[ifuel] = phi[0] x[io2] = stoich_O2 x[in2] = stoich_O2*air_N2_O2_ratio gas.TPX = tin, p, x #Create the free laminar premixed flame f = FreeFlame(gas, initial_grid) f.flame.set_steady_tolerances(default=tol_ss) f.flame.set_transient_tolerances(default=tol_ts) f.inlet.X = x f.inlet.T = tin f.transport_model = 'Mix' #First flame: #No energy for starters f.energy_enabled = False #Refinement criteria f.set_refine_criteria(ratio = 10.0, slope = 1, curve = 1) #Max number of times the Jacobian will be used before it must be re-evaluated f.set_max_jac_age(50, 50) #Set time steps whenever Newton convergence fails f.set_time_step(1.0e-6, [2, 5, 10, 20, 80]) #s #Calculation f.solve(loglevel, refine_grid) ################# #Second flame: #Energy equation enabled f.energy_enabled = True #Refinement criteria when energy equation is enabled f.set_refine_criteria(ratio = 3.0, slope = 0.1, curve = 0.2) #Calculation and save of the results f.solve(loglevel, refine_grid) ################# #Third flame and so on ...: #Refinement criteria should be changed ... f.set_refine_criteria(ratio = 2.0, slope = 0.05, curve = 0.05) f.solve(loglevel, refine_grid) sl_cas[i][0] = f.u[0] t_cas[i][0] = f.T[f.flame.n_points-1] #...and saving of the results f.save('c3h8-air_adiabatic-'+str(cas)+'.xml','energy_0', 'solution with the energy equation enabled') ################################################################# #THEN : ################################################################# for j in range(1,npoints): phi[j] = phi_min + (phi_max - phi_min)*j/(npoints - 1) x = zeros(m,'d') x[ifuel] = phi[j] x[io2] = stoich_O2 x[in2] = stoich_O2*air_N2_O2_ratio #Set gas state to that of the unburned gas gas.TPX = tin, p, x #Create the free laminar premixed flame f = FreeFlame(gas, initial_grid) f.transport_model = 'Mix' f.restore('c3h8-air_adiabatic-'+str(cas)+'.xml', 'energy_' + str(j-1)) f.inlet.X = x f.inlet.T = tin f.P = p f.solve(loglevel, refine_grid) sl_cas[i][j] = f.u[0] t_cas[i][j] = f.T[f.flame.n_points-1] #...and saving of the results f.save('c3h8-air_adiabatic-'+str(cas)+'.xml','energy_'+str(j), 'solution with the energy equation enabled') ################################################################# # Save your results if needed ################################################################# #Write the velocity, temperature, density, and mole fractions to a CSV file f.write_csv('c3h8-air_adiabatic-'+str(cas)+'_phi'+str(phi[j])+'.csv', quiet=False) #Write the velocity and temperature for comparison csv_file = 'compa_data_c3h8.csv' with open(csv_file, 'w') as outfile: writer = csv.writer(outfile) for i, cas in enumerate(matrice_cas): writer.writerow([cas]) for j in range(npoints): writer.writerow([phi[j], sl_cas[i][j], t_cas[i][j]]) ################################################################# # Plot your results ################################################################# # create plot fig=figure(1) # create first subplot - adiabatic flame temperature a=fig.add_subplot(111) for i, cas in enumerate(matrice_cas): a.plot(phi,sl_cas[i], type_plot[i], label = cas) hold(True) title(r'Sl vs. $\Phi$') xlabel(r'$\Phi$', fontsize=20) ylabel("Laminar flame speed [m/s]") #a.axis([0.5,1.7,0.0,0.4]) #ax = gca() #ax.set_autoscale_on(False) a.xaxis.set_major_locator(MaxNLocator(13)) # this controls the number of tick marks on the axis a.yaxis.set_major_locator(MaxNLocator(9)) # this controls the number of tick marks on the axis hold(False) legend(bbox_to_anchor=(0.9, 0.9,1,1),loc=8) #fig.text(0.5,0.95,r'Laminar flame speed for $C_{2}H_{2}$ + Air, at P = '+ str(p)+'Pa, Tin = '+str(tin)+'K.',fontsize=22,horizontalalignment='center') grid() subplots_adjust(left=0.08, right=0.96, wspace=0.25) #show() savefig('plot_flamespeed-'+str(tin)+'-'+str(p)+'.png', bbox_inches='tight')