From d671ca745c81436bad62f04ddc6215224bd4bcbf Mon Sep 17 00:00:00 2001
From: Dave Goodwin <dgg@caltech.edu>
Date: Mon, 15 Nov 2004 03:17:15 +0000
Subject: [PATCH] initial import

---
 Cantera/cxx/demos/flamespeed.cpp | 214 +++++++++++++++++++++++++++++++
 1 file changed, 214 insertions(+)
 create mode 100644 Cantera/cxx/demos/flamespeed.cpp

diff --git a/Cantera/cxx/demos/flamespeed.cpp b/Cantera/cxx/demos/flamespeed.cpp
new file mode 100644
index 000000000..d99682a90
--- /dev/null
+++ b/Cantera/cxx/demos/flamespeed.cpp
@@ -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;
+    }
+}