research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Examples] General cleanup of C++ flamespeed example

Browse source
This commit is contained in:
Ray Speth 2016-04-14 15:21:54 -04:00
parent b50ef03838
commit f2b68d8c34
2 changed files with 45 additions and 71 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

113
samples/cxx/flamespeed/flamespeed.cpp
View file

@ -8,10 +8,10 @@
#include "cantera/oneD/StFlow.h" #include "cantera/oneD/StFlow.h"
#include "cantera/IdealGasMix.h" #include "cantera/IdealGasMix.h"
#include "cantera/transport.h" #include "cantera/transport.h"
#include <fstream>
using namespace Cantera; using namespace Cantera;
using std::cout; using fmt::print;
using std::endl;
int flamespeed(double phi) int flamespeed(double phi)
{ {
@ -20,45 +20,31 @@ int flamespeed(double phi)
doublereal temp = 300.0; // K doublereal temp = 300.0; // K
doublereal pressure = 1.0*OneAtm; //atm doublereal pressure = 1.0*OneAtm; //atm
doublereal uin=0.3; //m/sec doublereal uin = 0.3; //m/sec
gas.setState_TPX(temp, pressure, "CH4:1.0, O2:2.0, N2:7.52");
size_t nsp = gas.nSpecies(); size_t nsp = gas.nSpecies();
vector_fp x(nsp, 0.0);
vector_fp x(nsp); 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);
x[gas.speciesIndex("CH4")] = 1.0;
x[gas.speciesIndex("O2")] = 0.21 / phi / fa_stoic;
x[gas.speciesIndex("N2")] = 0.79 / phi/ fa_stoic;
doublereal C_atoms=1.0; gas.setState_TPX(temp, pressure, x.data());
doublereal H_atoms=4.0; doublereal rho_in = gas.density();
doublereal ax=C_atoms+H_atoms/4.0;
doublereal fa_stoic=1.0/(4.76*ax);
for (size_t 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.data());
doublereal rho_in=gas.density();
vector_fp yin(nsp); vector_fp yin(nsp);
gas.getMassFractions(&yin[0]); gas.getMassFractions(&yin[0]);
try { gas.equilibrate("HP");
gas.equilibrate("HP");
} catch (CanteraError& err) {
std::cout << err.what() << std::endl;
}
vector_fp yout(nsp); vector_fp yout(nsp);
gas.getMassFractions(&yout[0]); gas.getMassFractions(&yout[0]);
doublereal rho_out = gas.density(); doublereal rho_out = gas.density();
doublereal Tad=gas.temperature(); doublereal Tad = gas.temperature();
cout << phi<<' '<<Tad<<endl; print("phi = {}, Tad = {}\n", phi, Tad);
//============= build each domain ======================== //============= build each domain ========================
@ -69,11 +55,11 @@ int flamespeed(double phi)
// create an initial grid // create an initial grid
int nz = 6; int nz = 6;
doublereal lz=0.1; doublereal lz = 0.1;
vector_fp z(nz); vector_fp z(nz);
doublereal dz=lz/((doublereal)(nz-1)); doublereal dz = lz/((doublereal)(nz-1));
for (int iz=0; iz<nz; iz++) { for (int iz = 0; iz < nz; iz++) {
z[iz]=((doublereal)iz)*dz; z[iz] = ((doublereal)iz)*dz;
} }
flow.setupGrid(nz, &z[0]); flow.setupGrid(nz, &z[0]);
@ -112,7 +98,7 @@ int flamespeed(double phi)
vector_fp locs{0.0, 0.3, 0.7, 1.0}; vector_fp locs{0.0, 0.3, 0.7, 1.0};
vector_fp value; vector_fp value;
double uout=inlet.mdot()/rho_out; double uout = inlet.mdot()/rho_out;
value = {uin, uin, uout, uout}; value = {uin, uin, uout, uout};
flame.setInitialGuess("u",locs,value); flame.setInitialGuess("u",locs,value);
value = {temp, temp, Tad, Tad}; value = {temp, temp, Tad, Tad};
@ -129,10 +115,10 @@ int flamespeed(double phi)
flame.showSolution(); flame.showSolution();
int flowdomain=1; int flowdomain = 1;
double ratio=10.0; double ratio = 10.0;
double slope=0.08; double slope = 0.08;
double curve=0.1; double curve = 0.1;
flame.setRefineCriteria(flowdomain,ratio,slope,curve); flame.setRefineCriteria(flowdomain,ratio,slope,curve);
@ -149,58 +135,49 @@ int flamespeed(double phi)
flame.solve(loglevel,refine_grid); flame.solve(loglevel,refine_grid);
double flameSpeed_mix = flame.value(flowdomain,flow.componentIndex("u"),0); double flameSpeed_mix = flame.value(flowdomain,flow.componentIndex("u"),0);
cout << "Flame speed with mixture-averaged transport: " << print("Flame speed with mixture-averaged transport: {} m/s\n",
flame.value(flowdomain,flow.componentIndex("u"),0) << " m/s" << endl; flameSpeed_mix);
// now switch to multicomponent transport // now switch to multicomponent transport
flow.setTransport(*trmulti); flow.setTransport(*trmulti);
flame.solve(loglevel, refine_grid); flame.solve(loglevel, refine_grid);
double flameSpeed_multi = flame.value(flowdomain,flow.componentIndex("u"),0); double flameSpeed_multi = flame.value(flowdomain,flow.componentIndex("u"),0);
cout << "Flame speed with multicomponent transport: " << print("Flame speed with multicomponent transport: {} m/s\n",
flame.value(flowdomain,flow.componentIndex("u"),0) << " m/s" << endl; flameSpeed_multi);
// now enable Soret diffusion // now enable Soret diffusion
flow.enableSoret(true); flow.enableSoret(true);
flame.solve(loglevel, refine_grid); flame.solve(loglevel, refine_grid);
double flameSpeed_full = flame.value(flowdomain,flow.componentIndex("u"),0); double flameSpeed_full = flame.value(flowdomain,flow.componentIndex("u"),0);
cout << "Flame speed with multicomponent transport + Soret: " << print("Flame speed with multicomponent transport + Soret: {} m/s\n",
flame.value(flowdomain,flow.componentIndex("u"),0) << " m/s" << endl; flameSpeed_full);
int np=flow.nPoints();
vector_fp zvec,Tvec,COvec,CO2vec,Uvec; vector_fp zvec,Tvec,COvec,CO2vec,Uvec;
printf("\n%9s\t%8s\t%5s\t%7s\n","z (m)", "T (K)", "U (m/s)", "Y(CO)"); print("\n{:9s}\t{:8s}\t{:5s}\t{:7s}\n",
for (int n=0; n<np; n++) { "z (m)", "T (K)", "U (m/s)", "Y(CO)");
for (size_t n = 0; n < flow.nPoints(); n++) {
Tvec.push_back(flame.value(flowdomain,flow.componentIndex("T"),n)); Tvec.push_back(flame.value(flowdomain,flow.componentIndex("T"),n));
COvec.push_back(flame.value(flowdomain,flow.componentIndex("CO"),n)); COvec.push_back(flame.value(flowdomain,flow.componentIndex("CO"),n));
CO2vec.push_back(flame.value(flowdomain,flow.componentIndex("CO2"),n)); CO2vec.push_back(flame.value(flowdomain,flow.componentIndex("CO2"),n));
Uvec.push_back(flame.value(flowdomain,flow.componentIndex("u"),n)); Uvec.push_back(flame.value(flowdomain,flow.componentIndex("u"),n));
zvec.push_back(flow.grid(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]); print("{: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; print("\nAdiabatic flame temperature from equilibrium is: {}\n", Tad);
cout << "Flame speed for phi="<<phi<<" is "<<Uvec[0]<<" m/s."<<endl; print("Flame speed for phi={} is {} m/s.\n", phi, Uvec[0]);
std::string reportFile = "flamespeed.csv"; std::ofstream outfile("flamespeed.csv", std::ios::trunc);
FILE* FP = fopen(reportFile.c_str(), "w"); outfile << " Grid, Temperature, Uvec, CO, CO2\n";
if (!FP) { for (size_t n = 0; n < flow.nPoints(); n++) {
printf("Failure to open file\n"); print(outfile, " {:11.3e}, {:11.3e}, {:11.3e}, {:11.3e}, {:11.3e}\n",
exit(-1); flow.grid(n), Tvec[n], Uvec[n], COvec[n], CO2vec[n]);
} }
fprintf(FP," Flame speed (mixture-averaged ) = %11.3e m/s\n", flameSpeed_mix);
fprintf(FP," Flame speed (multicomponent ) = %11.3e m/s\n", flameSpeed_multi);
fprintf(FP," Flame speed (multicomponent + Soret) = %11.3e m/s\n", flameSpeed_full);
fprintf(FP," Grid, Temperature, Uvec, CO, CO2\n");
for (int n = 0; n < np; n++) {
fprintf(FP," %11.3e, %11.3e, %11.3e, %11.3e, %11.3e\n",
flow.grid(n), Tvec[n], Uvec[n], COvec[n], CO2vec[n]);
}
fclose(FP);
} catch (CanteraError& err) { } catch (CanteraError& err) {
std::cerr << err.what() << std::endl; std::cerr << err.what() << std::endl;
std::cerr << "program terminating." << endl; std::cerr << "program terminating." << std::endl;
return -1; return -1;
} }
return 0; return 0;
@ -209,7 +186,7 @@ int flamespeed(double phi)
int main() int main()
{ {
double phi; double phi;
cout << "Enter phi: "; print("Enter phi: ");
std::cin >> phi; std::cin >> phi;
return flamespeed(phi); return flamespeed(phi);
} }

3
samples/cxx/flamespeed/flamespeed_blessed_0.csv
View file

@ -1,6 +1,3 @@
Flame speed (mixture-averaged ) = 3.461e-01 m/s
Flame speed (multicomponent ) = 3.516e-01 m/s
Flame speed (multicomponent + Soret) = 3.516e-01 m/s
Grid, Temperature, Uvec, CO, CO2 Grid, Temperature, Uvec, CO, CO2
0.000e+00, 3.000e+02, 3.516e-01, 4.128e-11, 2.311e-12 0.000e+00, 3.000e+02, 3.516e-01, 4.128e-11, 2.311e-12
1.202e-03, 3.000e+02, 3.516e-01, 8.802e-10, 6.402e-11 1.202e-03, 3.000e+02, 3.516e-01, 8.802e-10, 6.402e-11

1 Flame speed (mixture-averaged ) = 3.461e-01 m/s Grid Temperature Uvec CO CO2
Flame speed (mixture-averaged ) = 3.461e-01 m/s
Flame speed (multicomponent ) = 3.516e-01 m/s
Flame speed (multicomponent + Soret) = 3.516e-01 m/s
1 Grid, Temperature, Uvec, CO, CO2 Grid Temperature Uvec CO CO2
2 0.000e+00, 3.000e+02, 3.516e-01, 4.128e-11, 2.311e-12 0.000e+00 3.000e+02 3.516e-01 4.128e-11 2.311e-12
3 1.202e-03, 3.000e+02, 3.516e-01, 8.802e-10, 6.402e-11 1.202e-03 3.000e+02 3.516e-01 8.802e-10 6.402e-11