#!/usr/bin/env python """ Constant-pressure, adiabatic kinetics simulation. Author: Antony Misdariis This script permit to run several PSR in a raw over several pressures and temperatures. The timestep is computed dynamically, for each pressure and temperature it is based on the previous conditions AI delay. This assumes that this delay does not vary to much between 2 conditions. The dynamic timestep can be switched off using the parameter ddt Note: if a non-stoichimetric mixture is used, the stop condition must be changed (line 168) """ import sys import os from Cantera import * from Cantera.Reactor import * from Cantera.Func import * import numpy as np import math import matplotlib.pyplot as plt from pylab import * # --------------- Parameters of the simulation ---------------------------------------- npoints = 500 # Number of timesteps to perform each time AI is not detected phi=1 # Equivalence ratio Pmin = 41 # Min pressure to compute Pmax = 41 # Max pressure to compute Pstep = 5 # Step in pressure Tmin = 800 # Min temperature to compute Tmax = 800 # Max temperature to compute Tstep = 20 # Step in temperature isoP = False # If true : constant pressure reactor, else consatnt Volume write = True # Write AI delay in output file (name to modify at line 125) ddt=0 # if =0 dynamic timestep, otherwise ddt=timestep AI0 = 0.0051e0 # only if ddt=0, guess of the first AI time ( used to compute the timestep) # Fuel Definition cti_file='anderlhor.cti' # name of the cti file phase='gas' # phase name in cti # In this case fuel is composed of 3 primary fuel otherwise directly give the fuel molar mass w_isoo=8*12+18 w_nhep=7*12+16 w_tolu=7*12+8 wf=0.428*w_isoo+0.137*w_nhep+0.435*w_tolu # Fuel molar mass a=3.76 # N2/O2 molar ratio in air s=3.36264 # stoichiometric coefficient wN2=28 # N2 molar mass wO2=32 # O2 molar mass yf=(phi/s)/(phi/s+1+a*wN2/wO2) yo=1/(phi/s+1+a*wN2/wO2) yn=1-yf-yo w=1/(yf/wf+yo/wO2+yn/wN2) xo=(w/wO2)*yo xf=(w/wf)*yf xn=(w/wN2)*yn # IFP essence modele x_isoo=0.428*xf x_nhep=0.137*xf x_tolu=0.435*xf y_isoo=(w_isoo/w)*x_isoo y_nhep=(w_nhep/w)*x_nhep y_tolu=(w_tolu/w)*x_tolu yf0=y_isoo+y_nhep+y_tolu # Final composition (mass_fractions) compo="C8H18-1:"+str(y_isoo)+" C7H16-1:"+str(y_nhep)+" toluene:"+str(y_tolu)+" O2:"+str(yo)+" N2:"+str(yn) print "O2 mass fraction :"+str(yo) print "N2 mass fraction :"+str(yn) print "Fuel mass fraction :"+str(yf) print " -> iso-octane mass fraction :"+str(y_isoo) print " -> n-heptane mass fraction :"+str(y_nhep) print " -> toluene mass fraction :"+str(y_tolu) # --------------------------------------------------------------------------------------- # ---------------- Some initialisations ------------------------------------------------- args = sys.argv first = True # --------------------------------------------------------------------------------------- # -------------------------------- Functions ------------------------------------------- def plotPT (P,T,first): if not first: plt.close(); plt.ioff(); tmp=np.loadtxt('xmgrace.txt') clf subplot(2,1,1) title('PSR at P='+str(P)+'bar, T='+str(T)+'K') plt.plot(tmp[:,0],tmp[:,1],marker='o'); xlabel('Time (s)'); ylabel('Temperature (K)'); subplot(2,1,2) plt.plot(tmp[:,0],tmp[:,2],marker='o'); xlabel('Time (s)'); ylabel('Pressure (Pa)'); ion() plt.show(); plt.draw(); def plotX(P,T): tmp=np.loadtxt('mole_fractions.txt') plt.plot(tmp[:,0],tmp[:,1],'k',label='X(C8H18)') plt.plot(tmp[:,0],tmp[:,2],'r',label='X(CO2)') plt.plot(tmp[:,0],tmp[:,3],'g',label='X(H2O)') plt.legend(loc='lower right',prop={'size':12}) plt.title('PSR at P='+str(P)+'bar, T='+str(T)+'K') plt.show(); # --------------------------------------------------------------------------------------- # --------------------------------- Loop on all conditions ------------------------------ for P in range(Pmin,Pmax+1,Pstep): pressure = P*1E5 # Name of the output file if write: output_name="AI_P"+str(P)+"_AND_essence.dat" if not os.path.exists(output_name): print " File :"+output_name+" not found: will be created" fichier = open(output_name, "w") fichier.write( '#Time Temperature\n') fichier.close() AI_prev=AI0 # Boucle sur les temperatures for Tempi in range(Tmin,Tmax+1,Tstep): if ddt == 0: timestep = AI_prev/1E3 else: timestep = ddt # Compute PSR gri3 = importPhase(cti_file, phase) gri3.set(T = Tempi, P = pressure ,Y = compo) r = Reactor(gri3) if isoP: neutre=importPhase('neutre.cti', 'gas') neutre.set(T = Tempi, P = pressure ,Y = 'N2:1') env = Reservoir(neutre) # Define a wall between the reactor and the environment, and # make it flexible, so that the pressure in the reactor is held # at the environment pressure. w = Wall(r,env) w.set(K = 1.0e2) # set expansion parameter. dV/dt = KA(P_1 - P_2) w.set(A = 1.0) sim = ReactorNet([r]) time = 0.0 tim = zeros(npoints,'d') tempe = zeros(npoints,'d') fichier1 = open("xmgrace.txt", "w") fichier2 = open("mole_fractions.txt", "w") yf_end = yf0 tot_av = 0 exact=0 print "Computing auto-ignition delay for T="+str(Tempi)+' K, P='+str(P)+' bar, dt='+str(timestep)+' s' #------------------------------------------------------------------------------------------------- # If your reactor is not stoichiometric the stop condition have to be adjusted # Here the simulation stops when there all the fuel is burned (when its mass fraction # is lower than 1e-6 time its initial value) if phi != 1: print " WARNING: The stop test have to be modified" #------------------------------------------------------------------------------------------------- while yf_end >= 1E-6*yf0: for n in range(npoints): time += timestep sim.advance(time) yf_end=r.massFraction('C8H18-1')+r.massFraction('C7H16-1')+r.massFraction('toluene') if tot_av >= exact: exact +=10 print ' - Progression: '+str(int(tot_av))+'%' tim[n] = time tempe[n] = r.temperature() fichier1.write( '%10.3e %10.6f %10.6f \n' % (sim.time(),r.temperature(),r.pressure())) fichier2.write( '%10.3e %10.6f %10.6f %10.6f \n' % (sim.time(),r.moleFraction('C8H18-1'),r.moleFraction('CO2'),r.moleFraction('H2O'))) fichier1.close() fichier2.close() #Calcul de la derive de la temperature Temperature=np.array(tempe) Time=np.array(tim) dTdt = np.diff(Temperature) dTdt_max=dTdt.max() dTdt_max_index=np.where(dTdt==dTdt_max)[0] t_AI=Time[dTdt_max_index] if write: fichier = open(output_name, "a") fichier.write( '{0[0]:.5} {1} \n'.format(t_AI,Tempi) ) fichier.close() if len(args) > 1 and args[1] == '-plotPT': plotPT(P,Tempi,first) if len(args) > 1 and args[1] == '-plotX': plotX(P,Tempi) first = False print "Temps d'autoallumage:"+str(float(t_AI))+' s' AI_prev=float(t_AI)