initial import
This commit is contained in:
Dave Goodwin
parent
488f84970e
commit
d671ca745c
1 changed files with 214 additions and 0 deletions
214
Cantera/cxx/demos/flamespeed.cpp
Normal file
214
Cantera/cxx/demos/flamespeed.cpp
Normal file
|
|
@ -0,0 +1,214 @@
|
|||
#ifdef WIN32
|
||||
#pragma warning(disable:4786)
|
||||
/*Warning messages that are identified as Warning C4786: are created
|
||||
when MSVC generates extremely long names that it uses for debugging
|
||||
purposes. The long names are generated by the template expansion
|
||||
process and the warning messages normally can be ignored. Since these
|
||||
warnings tend to hide more interesting warning/error messages, you may
|
||||
wish to suppress the warning*/
|
||||
#endif
|
||||
|
||||
#include <cantera/Cantera.h>
|
||||
#include <cantera/IdealGasMix.h>
|
||||
#include <cantera/transport.h>
|
||||
#include <cantera/equilibrium.h>
|
||||
#include <cantera/onedim.h>
|
||||
|
||||
// include additional header files here if needed
|
||||
#ifdef WIN32
|
||||
#include <conio.h>
|
||||
#endif
|
||||
|
||||
int main(int argc, char** argv) {
|
||||
try {
|
||||
int i;
|
||||
// your code goes here
|
||||
IdealGasMix gas("gri30.cti","gri30_mix");
|
||||
//IdealGasMix gas("temp_phi_1.0_sp.cti","gas");
|
||||
|
||||
doublereal temp = 300.0; // K
|
||||
doublereal pressure = 1.0*OneAtm; //atm
|
||||
doublereal uin=0.3; //m/sec
|
||||
|
||||
gas.setState_TPX(temp, pressure, "CH4:1.0, O2:2.0, N2:7.52");
|
||||
int nsp = gas.nSpecies();
|
||||
|
||||
vector_fp x;
|
||||
x.resize(nsp);
|
||||
|
||||
doublereal phi=1.1;
|
||||
cout << "Enter phi: ";
|
||||
cin >> phi;
|
||||
cout << endl;
|
||||
|
||||
doublereal C_atoms=1.0;
|
||||
doublereal H_atoms=4.0;
|
||||
doublereal ax=C_atoms+H_atoms/4.0;
|
||||
doublereal fa_stoic=1.0/(4.76*ax);
|
||||
for(int k=0;k<nsp;k++){
|
||||
if(k==gas.speciesIndex("CH4")){ x[k]=1.0; }
|
||||
else if(k==gas.speciesIndex("O2")){ x[k]=0.21/phi/fa_stoic; }
|
||||
else if(k==gas.speciesIndex("N2")){ x[k]=0.79/phi/fa_stoic; }
|
||||
else{ x[k]=0.0;
|
||||
}
|
||||
}
|
||||
gas.setState_TPX(temp,pressure,x.begin());
|
||||
doublereal rho_in=gas.density();
|
||||
|
||||
//vector<doublereal> yin;
|
||||
double *yin=new double[nsp];
|
||||
gas.getMassFractions(yin);
|
||||
|
||||
try {
|
||||
equilibrate(gas,HP);
|
||||
}
|
||||
catch (CanteraError) {
|
||||
showErrors(cout);
|
||||
}
|
||||
double *yout=new double[nsp];
|
||||
gas.getMassFractions(yout);
|
||||
doublereal rho_out = gas.density();
|
||||
doublereal Tad=gas.temperature();
|
||||
cout << phi<<' '<<Tad<<endl;
|
||||
|
||||
double Tin=temp;
|
||||
double Tout=Tad;
|
||||
double breakpt=0.2;
|
||||
//============= build each domain ========================
|
||||
|
||||
|
||||
//-------- step 1: create the stagnation flow -------------
|
||||
|
||||
AxiStagnFlow flow(&gas);
|
||||
|
||||
// create an initial grid
|
||||
//doublereal z[] = {0.0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009,0.0095, 0.01};
|
||||
int nz=5;
|
||||
doublereal lz=0.02;
|
||||
doublereal *z=new double[nz+1];
|
||||
doublereal dz=lz/((doublereal)(nz-1));
|
||||
for(int iz=0;iz<nz;iz++){
|
||||
z[iz]=((doublereal)iz)*dz;
|
||||
}
|
||||
//add one node onto end of domain to help with zero gradient at outlet
|
||||
z[nz]=lz*1.05;
|
||||
nz++;
|
||||
|
||||
flow.setupGrid(nz, z);
|
||||
|
||||
// specify the objects to use to compute kinetic rates and
|
||||
// transport properties
|
||||
|
||||
Transport* tr = newTransportMgr("Mix", &gas);
|
||||
flow.setTransport(*tr);
|
||||
flow.setKinetics(gas);
|
||||
flow.setPressure(pressure);
|
||||
|
||||
//------- step 2: create the inlet -----------------------
|
||||
|
||||
Inlet1D inlet;
|
||||
|
||||
inlet.setMoleFractions(x.begin());
|
||||
doublereal mdot=uin*rho_in;
|
||||
inlet.setMdot(mdot);
|
||||
inlet.setTemperature(temp);
|
||||
|
||||
|
||||
//------- step 3: create the outlet ---------------------
|
||||
|
||||
Outlet1D outlet;
|
||||
|
||||
//=================== create the container and insert the domains =====
|
||||
|
||||
vector<Domain1D*> domains;
|
||||
domains.push_back(&inlet);
|
||||
domains.push_back(&flow);
|
||||
domains.push_back(&outlet);
|
||||
|
||||
OneDim flamesim(domains);
|
||||
|
||||
Sim1D flame(domains);
|
||||
|
||||
//----------- Supply initial guess----------------------
|
||||
|
||||
vector_fp locs;
|
||||
vector_fp value;
|
||||
|
||||
locs.resize(3);
|
||||
value.resize(3);
|
||||
|
||||
//ramp values from inlet to adiabatic flame conditions
|
||||
// over 70% of domain and then level off at equilibrium
|
||||
double z1=0.7;
|
||||
|
||||
double uout;
|
||||
uout=inlet.mdot()/rho_out;
|
||||
uin=inlet.mdot()/rho_in;
|
||||
locs[0]=0.0; locs[1]=z1; locs[2]=1.0;
|
||||
value[0]=uin; value[1]=uout; value[2]=uout;
|
||||
flame.setInitialGuess("u",locs,value);
|
||||
|
||||
value[0]=temp; value[1]=Tad; value[2]=Tad;
|
||||
flame.setInitialGuess("T",locs,value);
|
||||
|
||||
for(i=0;i<nsp;i++){
|
||||
value[0]=yin[i]; value[1]=yout[i]; value[2]=yout[i];
|
||||
flame.setInitialGuess(gas.speciesName(i),locs,value);
|
||||
}
|
||||
|
||||
inlet.setMoleFractions(x.begin());
|
||||
inlet.setMdot(mdot);
|
||||
inlet.setTemperature(temp);
|
||||
|
||||
flame.showSolution();
|
||||
|
||||
int flowdomain=1;
|
||||
double ratio=10.0;
|
||||
double slope=0.2;
|
||||
double curve=0.02;
|
||||
double prune=-0.00005;
|
||||
|
||||
flame.setRefineCriteria(flowdomain,ratio,slope,curve,prune);
|
||||
|
||||
int loglevel=1;
|
||||
bool refine_grid = true;
|
||||
|
||||
/* Solve species*/
|
||||
flow.fixTemperature();
|
||||
refine_grid=false;
|
||||
flame.solve(loglevel,refine_grid);
|
||||
|
||||
/* Solve freely propagating flame*/
|
||||
/* Linearally interpolate to find location where this temperature would exist
|
||||
temperature at this location will then be fixed for remainder of calculation.*/
|
||||
|
||||
flow.solveEnergyEqn();
|
||||
refine_grid=true;
|
||||
flame.setFixedTemperature(900.0);
|
||||
flame.setAdiabaticFlame();
|
||||
flame.solve(loglevel=1,refine_grid);
|
||||
|
||||
int np=flow.nPoints();
|
||||
vector<doublereal> zvec,Tvec,COvec,CO2vec,Uvec;
|
||||
|
||||
printf("\n%9s\t%8s\t%5s\t%7s\n","z (m)", "T (K)", "U (m/s)", "Y(CO)");
|
||||
for(int n=0;n<np;n++){
|
||||
Tvec.push_back(flame.value(flowdomain,flow.componentIndex("T"),n));
|
||||
COvec.push_back(flame.value(flowdomain,flow.componentIndex("CO"),n));
|
||||
CO2vec.push_back(flame.value(flowdomain,flow.componentIndex("CO2"),n));
|
||||
Uvec.push_back(flame.value(flowdomain,flow.componentIndex("u"),n));
|
||||
zvec.push_back(flow.grid(n));
|
||||
printf("%9.6f\t%8.3f\t%5.3f\t%7.5f\n",flow.grid(n),Tvec[n],Uvec[n],COvec[n]);
|
||||
}
|
||||
|
||||
cout << endl<<"Adiabatic flame temperature from equilibrium is: "<<Tad<<endl;
|
||||
cout << "Flame speed for phi="<<phi<<" is "<<Uvec[0]<<" m/s."<<endl;
|
||||
|
||||
return 0;
|
||||
}
|
||||
catch (CanteraError) {
|
||||
showErrors(cerr);
|
||||
cerr << "program terminating." << endl;
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
Loading…
Add table
Reference in a new issue