ODT
Loading...
Searching...
No Matches
streams.cc
Go to the documentation of this file.
1
6#include "streams.h"
7#include "domain.h"
8#include "yaml-cpp/yaml.h"
9#include <cassert>
10
11using namespace std;
12
14
18void streams::init(domain *p_domn, const vector<double> &gammas) {
19
20 domn = p_domn;
21
22 nspc = domn->gas->nSpecies();
23
24 T0 = domn->io->streamProps["T0"].as<double>();
25 T1 = domn->io->streamProps["T1"].as<double>();
26
27 y0 = vector<double>(nspc, 0.0);
28 y1 = vector<double>(nspc, 0.0);
29
30 D0 = vector<double>(nspc, 0.0);
31 D1 = vector<double>(nspc, 0.0);
32
33 hsp0 = vector<double>(nspc, 0.0);
34 hsp1 = vector<double>(nspc, 0.0);
35
36 bool Lmole = domn->io->streamProps["moleOrMass"].as<string>() == "MOLE" ? true : false;
37
38 gCHON = gammas;
39
40 //--------------- stream properties: T,h,y.
41
42 double sum = 0.0; // for normalizing
43 vector<string> y0Labels; // species names in list
44 vector<double> y0short; // value corresponding to y0Labels
45 YAML::Node yy = domn->io->streamProps["comp0"];
46 for(YAML::const_iterator it = yy.begin(); it!=yy.end(); it++) {
47 y0Labels.push_back(it->first.as<string>());
48 y0short.push_back(it->second.as<double>());
49 sum += it->second.as<double>();
50 }
51 for(int i=0; i<y0Labels.size(); i++)
52 y0[domn->gas->speciesIndex(y0Labels[i])] = y0short[i]/sum;
53
54 sum = 0.0;
55 vector<string> y1Labels; // species names in list
56 vector<double> y1short; // value corresponding to y0Labels
57 yy = domn->io->streamProps["comp1"];
58 for(YAML::const_iterator it = yy.begin(); it!=yy.end(); it++) {
59 y1Labels.push_back(it->first.as<string>());
60 y1short.push_back(it->second.as<double>());
61 sum += it->second.as<double>();
62 }
63 for(int i=0; i<y1Labels.size(); i++)
64 y1[domn->gas->speciesIndex(y1Labels[i])] = y1short[i] / sum;
65
66 if(Lmole) {
67 domn->gas->setMoleFractions( &y0[0] );
68 domn->gas->getMassFractions( &y0[0] );
69 domn->gas->setMoleFractions( &y1[0] );
70 domn->gas->getMassFractions( &y1[0] );
71 }
72
73 domn->gas->setState_TPY( T0, domn->pram->pres, &y0[0] );
74 h0 = domn->gas->enthalpy_mass();
75 rho0 = domn->gas->density();
76 M0 = domn->gas->meanMolecularWeight();
77 domn->tran->getMixDiffCoeffs(&D0[0]);
78 domn->gas->getEnthalpy_RT(&hsp0[0]); // non-dimensional enthalpy
79 for (int k=0; k<nspc; k++)
80 hsp0[k] = hsp0[k] * T0 * GasConstant / domn->gas->molecularWeight(k); // J/kg
81
82 domn->gas->setState_TPY( T1, domn->pram->pres, &y1[0] );
83 h1 = domn->gas->enthalpy_mass();
84 rho1 = domn->gas->density();
85 M1 = domn->gas->meanMolecularWeight();
86 domn->tran->getMixDiffCoeffs(&D1[0]);
87 domn->gas->getEnthalpy_RT(&hsp1[0]); // non-dimensional enthalpy
88 for (int k=0; k<nspc; k++)
89 hsp1[k] = hsp1[k] * T1 * GasConstant / domn->gas->molecularWeight(k); // J/kg
90
91 //---------------- stoich mixf
92
93 setStoicMixf();
94
95 getMixtureFraction(&y0[0], true); // set beta0 and beta1
96
97}
99
106void streams::getMixingState(const double mixf, vector<double> &ymix,
107 double &hmix, double &Tmix){
108
109 hmix = h1*mixf + h0*(1.0-mixf);
110 for(int k=0; k<nspc; k++)
111 ymix[k] = y1[k]*mixf + y0[k]*(1.0-mixf);
112
113 domn->gas->setMassFractions( &ymix[0] );
114 domn->gas->setState_HP(hmix, domn->pram->pres, 1.E-10); // get temperature as Tadiabatic
115
116 Tmix = domn->gas->temperature();
117}
118
120
128void streams::getProdOfCompleteComb(const double mixf, vector<double> &ypcc,
129 double &hpcc, double &Tpcc){
130
131 //---------- Compute the mixing mass fractions and enthalpy
132
133 ypcc.resize(nspc);
134 double d1 = 1.0-mixf;
135 for(int k=0; k<nspc; k++) {
136 ypcc[k] = y1[k]*mixf + y0[k]*d1;
137 }
138 hpcc = h1*mixf + h0*d1;
139
140 //--------- Set gas and element indicicies
141
142 int iC = domn->gas->elementIndex("C");
143 int iH = domn->gas->elementIndex("H");
144 //int iO = domn->gas->elementIndex("O"); // !!!!! currently unusedvariable
145 int iN = domn->gas->elementIndex("N");
146 int iCO2 = domn->gas->speciesIndex("CO2");
147 int iH2O = domn->gas->speciesIndex("H2O");
148 int iN2 = domn->gas->speciesIndex("N2");
149 int iO2 = domn->gas->speciesIndex("O2");
150
151 //---------- Set ypcc as the mixing mole fractions: Take a basis of one mole:
152 // now we are working in moles
153 // elemM are moles of each element
154 // when stoic:
155 // CxHyNz + (x+y/4)O2 ==> (z/2)N2 + (x)CO2 + (y/2)H2O
156 // otherwise:
157 // CxHyNz + (beta)O2 ==> (z/2)N2 + (x)CO2 + (y/2)H2O
158 // Note this allows fuels with nitrogen and oxygen
159
160 domn->gas->setMassFractions( &ypcc[0] );
161 domn->gas->getMoleFractions( &ypcc[0] );
162
163 double nOnotFromO2 = 0.0;
164 double nHnotFromH2O = 0.0;
165 double nCnotFromCO2 = 0.0;
166 vector<double> elemM = getElementMoles( &ypcc[0], nOnotFromO2,
167 nHnotFromH2O, nCnotFromCO2 );
168
169 double x = elemM[iC];
170 double y = elemM[iH];
171 double z = elemM[iN];
172 double beta = elemM[domn->gas->elementIndex("O")] * 0.5; // moles of O as O2
173
174 //--------------------------------------------------------------------------
175
176 if(mixf < mixfStoic) { // lean: burn all fuel, leftover O2
177
178 for(int k=0; k<nspc; k++)
179 ypcc[k] = 0.0;
180
181 if(iCO2 > 0) // CO2 is not in the H2 mechs
182 ypcc[iCO2] = x;
183 ypcc[iH2O] = y*0.5;
184 ypcc[iN2] = z*0.5;
185 ypcc[iO2] = beta - (x+y/4.0);
186
187 }
188 else{ // rich: burn all O2, leftover fuel
189
190 //double eta = beta/(x+y/4.0); // extent of reaction
191 double eta = (beta-nOnotFromO2/2)/(x+y/4-nOnotFromO2/2); // extent of reaction
192 if(eta > 1.0)
193 *domn->io->ostrm << endl << "eta > 1.0" << endl;
194 d1 = 1.0-eta; // fraction of fuel unburnt
195
196 for(int k=0; k<nspc; k++)
197 ypcc[k] *= d1; // initialize everything then correct
198 if(iCO2 > 0) // CO2 is not in the H2 mechs
199 ypcc[iCO2] = (x-nCnotFromCO2) + nCnotFromCO2*eta; // init + formed
200 ypcc[iH2O] = (y-nHnotFromH2O)*0.5 + nHnotFromH2O*0.5*eta; // init + formed
201 ypcc[iN2] = z*0.5;
202 ypcc[iO2] = 0.0;
203
204 }
205
206 //--------------------------------------------------------------------------
207
208 double sum = 0.0; // normalize moles
209 for(int k=0; k<nspc; k++)
210 sum += ypcc[k];
211 assert(sum != 0.0);
212 for(int k=0; k<nspc; k++)
213 ypcc[k] /= sum;
214 domn->gas->setMoleFractions( &ypcc[0] ); // set mole fractions
215 domn->gas->getMassFractions( &ypcc[0] ); // set ypcc as mass fractions
216
217
218 //--------------------------------------------------------------------------
219
220 domn->gas->setState_HP(hpcc, domn->pram->pres, 1.E-10); // get temperature as Tadiabatic
221 //domn->gas->setState_HP(hpcc, domn->pram->pres); // get temperature as Tadiabatic
222 Tpcc = domn->gas->temperature();
223
224}
225
227
229void streams::setStoicMixf() {
230
231 vector<double> x_c(nspc,0.0);
232
233 double mc0, mc1, mo0, mo1, mh0, mh1;
234
235 vector<double> elemMassFrac0 = setElementMassFracs(&y0[0]);
236 vector<double> elemMassFrac1 = setElementMassFracs(&y1[0]);
237
238 mc0 = elemMassFrac0[domn->gas->elementIndex("C")]/
239 domn->gas->atomicWeight(domn->gas->elementIndex("C"));
240 mc1 = elemMassFrac1[domn->gas->elementIndex("C")]/
241 domn->gas->atomicWeight(domn->gas->elementIndex("C"));
242 mh0 = elemMassFrac0[domn->gas->elementIndex("H")]/
243 domn->gas->atomicWeight(domn->gas->elementIndex("H"));
244 mh1 = elemMassFrac1[domn->gas->elementIndex("H")]/
245 domn->gas->atomicWeight(domn->gas->elementIndex("H"));
246 mo0 = elemMassFrac0[domn->gas->elementIndex("O")]/
247 domn->gas->atomicWeight(domn->gas->elementIndex("O"));
248 mo1 = elemMassFrac1[domn->gas->elementIndex("O")]/
249 domn->gas->atomicWeight(domn->gas->elementIndex("O"));
250
251 mixfStoic = (2.0*mc0 + 0.5*mh0 - mo0) /
252 (mo1-mo0 + 2.0*(mc0-mc1) + 0.5*(mh0-mh1));
253
254 *domn->io->ostrm << endl << "# mixfStoic = m_fuel/(m_fuel+m_air) = " << mixfStoic << endl;
255
256}
257
259
264vector<double> streams::setElementMassFracs(const double *y) {
265
266
267 vector<double> atomArr(domn->gas->nElements());
268 vector<double> elemMassFrac(domn->gas->nElements(), 0.0);
269 double sum = 0.0;
270
271 domn->gas->setMassFractions( &y[0] );
272
273 for(int k=0; k<nspc; k++) {
274 sum=0.0;
275 domn->gas->getAtoms(k, &atomArr[0]); // [nelements] in sp k
276 for(int m=0; m<(int)atomArr.size(); m++)
277 sum += atomArr[m] * domn->gas->atomicWeight(m);
278 for(int m=0; m<(int)atomArr.size(); m++)
279 elemMassFrac[m] += y[k] * atomArr[m]/sum * domn->gas->atomicWeight(m);
280 // is * mass frac of elem in sp
281 }
282
283 return elemMassFrac;
284
285}
286
288
293vector<double> streams::setElementMoleFracs(const double *y) {
294
295
296 vector<double> atomArr(domn->gas->nElements());
297 vector<double> elemMoleFrac(domn->gas->nElements(), 0.0);
298
299 domn->gas->setMassFractions( &y[0] );
300
301 for(int k=0; k<nspc; k++) {
302 domn->gas->getAtoms(k, &atomArr[0]); // [nelements] in sp k
303 for(int m=0; m<(int)atomArr.size(); m++)
304 elemMoleFrac[m] += domn->gas->moleFraction(k) * atomArr[m];
305 }
306 double sum = 0.0;
307 for(int m=0; m<(int)atomArr.size(); m++)
308 sum += elemMoleFrac[m];
309 assert(sum != 0.0);
310 for(int m=0; m<(int)atomArr.size(); m++)
311 elemMoleFrac[m] /= sum;
312
313 return elemMoleFrac;
314
315}
316
318
326vector<double> streams::getElementMoles(const double *x,
327 double &nOnotFromO2,
328 double &nHnotFromH2O,
329 double &nCnotFromCO2) {
330
331
332 vector<double> atomArr(domn->gas->nElements());
333 vector<double> elemM(domn->gas->nElements(), 0.0);
334 int iO2 = domn->gas->speciesIndex("O2");
335 int iO = domn->gas->elementIndex("O");
336 int iCO2 = domn->gas->speciesIndex("CO2");
337 int iC = domn->gas->elementIndex("C");
338 int iH2O = domn->gas->speciesIndex("H2O");
339 int iH = domn->gas->elementIndex("H");
340
341 for(int k=0; k<nspc; k++) {
342 domn->gas->getAtoms(k, &atomArr[0]); // [nelements] in sp k
343 for(int m=0; m<(int)atomArr.size(); m++)
344 elemM[m] += x[k] * atomArr[m];
345 if(k != iO2) nOnotFromO2 += atomArr[iO] * x[k];
346 if(k != iCO2) nCnotFromCO2 += atomArr[iC] * x[k];
347 if(k != iH2O) nHnotFromH2O += atomArr[iH] * x[k];
348 }
349 return elemM;
350
351}
352
354
363double streams::getMixtureFraction(const double *y, const bool doBeta01) {
364
365 vector<double> elemMF;
366 double beta;
367
368 elemMF = setElementMassFracs(y);
369 beta = gCHON[0] * elemMF[domn->gas->elementIndex("C")] +
370 gCHON[1] * elemMF[domn->gas->elementIndex("H")] +
371 gCHON[2] * elemMF[domn->gas->elementIndex("O")] +
372 gCHON[3] * elemMF[domn->gas->elementIndex("N")];
373
374 if(doBeta01) {
375 elemMF = setElementMassFracs(&y0[0]);
376 beta0 = gCHON[0] * elemMF[domn->gas->elementIndex("C")] +
377 gCHON[1] * elemMF[domn->gas->elementIndex("H")] +
378 gCHON[2] * elemMF[domn->gas->elementIndex("O")] +
379 gCHON[3] * elemMF[domn->gas->elementIndex("N")];
380
381 elemMF = setElementMassFracs(&y1[0]);
382 beta1 = gCHON[0] * elemMF[domn->gas->elementIndex("C")] +
383 gCHON[1] * elemMF[domn->gas->elementIndex("H")] +
384 gCHON[2] * elemMF[domn->gas->elementIndex("O")] +
385 gCHON[3] * elemMF[domn->gas->elementIndex("N")];
386 }
387
388 if( beta1-beta0 == 0 )
389 return 0.0; // to avoid division by zero
390 else
391 return (beta - beta0)/(beta1 - beta0);
392
393}
394
domain::gas
IdealGasPhase * gas
pointer to cantera thermochemistry object (reaction rates, Cp, etc.)
Definition domain.h:69
domain::tran
Transport * tran
pointer to cantera transport object (viscosity, diffusivity, etc.)
Definition domain.h:70
domain::io
inputoutput * io
pointer to input/output object
Definition domain.h:72
domain::pram
param * pram
pointer to the parameters object
Definition domain.h:73
inputoutput::streamProps
YAML::Node streamProps
yaml sub node
Definition inputoutput.h:41
inputoutput::ostrm
ostream * ostrm
Runtime: points to cout or to a file.
Definition inputoutput.h:33
param::pres
double pres
initial pressure (Pa)
Definition param.h:40
streams::y1
vector< double > y1
stream mixf=1 composition vector
Definition streams.h:39
streams::init
void init(domain *p_domn, const vector< double > &gammas)
Definition streams.cc:18
streams::T0
double T0
stream mixf=0 temperature
Definition streams.h:34
streams::setStoicMixf
void setStoicMixf()
Definition streams.cc:229
streams::hsp1
vector< double > hsp1
stream mixf=1 species enthalpies
Definition streams.h:47
streams::beta1
double beta1
mixf = (beta-beta0) / (beta1-beta0)
Definition streams.h:56
streams::rho1
double rho1
stream mixf=1 density
Definition streams.h:43
streams::domn
domain * domn
pointer to domain object
Definition streams.h:32
streams::getElementMoles
vector< double > getElementMoles(const double *x, double &nOnotFromO2, double &nHnotFromH2O, double &nCnotFromCO2)
Definition streams.cc:326
streams::beta0
double beta0
Definition streams.h:54
streams::h1
double h1
stream mixf=1 enthalpy
Definition streams.h:37
streams::M1
double M1
stream mixf=1 mean molecular weight
Definition streams.h:41
streams::T1
double T1
stream mixf=1 temperature
Definition streams.h:35
streams::nspc
int nspc
number of species in gas mechanism
Definition streams.h:51
streams::D0
vector< double > D0
stream mixf=0 diffusivities
Definition streams.h:44
streams::hsp0
vector< double > hsp0
stream mixf=0 species enthalpies
Definition streams.h:46
streams::mixfStoic
double mixfStoic
stoichiometric mixture fraction
Definition streams.h:49
streams::setElementMassFracs
vector< double > setElementMassFracs(const double *y)
Definition streams.cc:264
streams::getMixingState
void getMixingState(const double mixf, vector< double > &ymix, double &hmix, double &Tmix)
Definition streams.cc:106
streams::M0
double M0
stream mixf=0 mean molecular weight
Definition streams.h:40
streams::h0
double h0
stream mixf=0 enthalpy
Definition streams.h:36
streams::setElementMoleFracs
vector< double > setElementMoleFracs(const double *y)
Definition streams.cc:293
streams::getMixtureFraction
double getMixtureFraction(const double *y, const bool doBeta01=false)
Definition streams.cc:363
streams::D1
vector< double > D1
stream mixf=1 diffusivities
Definition streams.h:45
streams::getProdOfCompleteComb
void getProdOfCompleteComb(const double mixf, vector< double > &ypcc, double &hpcc, double &Tpcc)
Definition streams.cc:128
streams::rho0
double rho0
stream mixf=0 density
Definition streams.h:42
streams::gCHON
vector< double > gCHON
gammas, as in beta = sum_i (y_i*gamma_i)
Definition streams.h:58
streams::y0
vector< double > y0
stream mixf=0 composition vector
Definition streams.h:38
domain.h
Header file for class domain.
streams.h
Header file for class streams.
Generated by doxygen 1.10.0