streams.cc

Go to the documentation of this file.

#include "streams.h"
#include <cassert>
#include <iostream>
 
using namespace std;
 
 
streams::streams(shared_ptr<Cantera::Solution> csol,
                 const double _P, 
                 const double _h0, 
                 const double _h1,
                 const vector<double> &_y0, 
                 const vector<double> &_y1) :
    P(_P),
    h0(_h0),
    h1(_h1),
    y0(_y0),
    y1(_y1) {
 
    //----------
 
    gas  = csol->thermo(); 
    nspc = gas->nSpecies();
 
    gCHON.resize(4);
    gCHON[0] =  2.0/gas->atomicWeight(gas->elementIndex("C"));
    gCHON[1] =  0.5/gas->atomicWeight(gas->elementIndex("H"));
    gCHON[2] = -1.0/gas->atomicWeight(gas->elementIndex("O"));
    gCHON[3] = 0.0;
 
    //----------------
 
    setStoicMixf();
 
    getMixtureFraction(&y0[0], true);    // set beta0 and beta1
}
 
 
void streams::getMixingState(const double mixf, vector<double> &ymix,
                             double &hmix, double &Tmix) {
 
    hmix = h1*mixf + h0*(1.0-mixf);
    for(int k=0; k<nspc; k++)
        ymix[k] = y1[k]*mixf + y0[k]*(1.0-mixf);
 
    gas->setMassFractions( &ymix[0] );
    gas->setState_HP(hmix, P, 1.E-10);    // get temperature as Tadiabatic
 
    Tmix = gas->temperature();
}
 
 
void streams::getEquilibrium_HP(const double mixf, vector<double> &yeq,
                                      double &heq, double &Teq) {
 
    //---------- Compute the mixing mass fractions and enthalpy
 
    yeq.resize(nspc);
    for(int k=0; k<nspc; k++)
        yeq[k] = y1[k]*mixf + y0[k]*(1.0-mixf);
    heq = h1*mixf + h0*(1.0-mixf);
 
    gas->setMassFractions(&yeq[0]);
    gas->setState_HP(heq, P, 1.E-10);
 
    gas->equilibrate("HP");
    gas->getMassFractions(&yeq[0]);
 
    Teq = gas->temperature();
 
}
 
 
void streams::getEquilibrium_TP(const double mixf, double Teq, 
                                      vector<double> &yeq, double &heq ) {
 
    //---------- Compute the mixing mass fractions and enthalpy
 
    yeq.resize(nspc);
    for(int k=0; k<nspc; k++)
        yeq[k] = y1[k]*mixf + y0[k]*(1.0-mixf);
 
    gas->setState_TPY(Teq, P, &yeq[0]);
 
    gas->equilibrate("TP");
    gas->getMassFractions(&yeq[0]);
 
    heq = gas->enthalpy_mass();
 
}
 
 
void streams::getProdOfCompleteComb(const double mixf, vector<double> &ypcc,
                                    double &hpcc, double &Tpcc) {
 
    //---------- Compute the mixing mass fractions and enthalpy
 
    ypcc.resize(nspc);
    double d1 = 1.0-mixf;
    for(int k=0; k<nspc; k++) {
        ypcc[k] = y1[k]*mixf + y0[k]*d1;
    }
    hpcc = h1*mixf + h0*d1;
 
    //--------- Set gas and element indicicies
 
    int iC = gas->elementIndex("C");
    int iH = gas->elementIndex("H");
    int iN = gas->elementIndex("N");
    int iCO2 = gas->speciesIndex("CO2");
    int iH2O = gas->speciesIndex("H2O");
    int iN2  = gas->speciesIndex("N2");
    int iO2  = gas->speciesIndex("O2");
 
    //---------- Set ypcc as the mixing mole fractions: Take a basis of one mole:
    // now we are working in moles
    // elemM are moles of each element
    // when stoic:
    // CxHyNz   + (x+y/4)O2  ==>  (z/2)N2 + (x)CO2 + (y/2)H2O
    // otherwise:
    // CxHyNz   + (beta)O2   ==>  (z/2)N2 + (x)CO2 + (y/2)H2O
    // Note this allows fuels with nitrogen and oxygen
 
    gas->setMassFractions( &ypcc[0] );
    gas->getMoleFractions( &ypcc[0] );
 
    double nOnotFromO2;
    double nHnotFromH2O;
    double nCnotFromCO2;
    vector<double> elemM = getElementMoles( &ypcc[0], nOnotFromO2,
            nHnotFromH2O, nCnotFromCO2 );
 
    double x    = elemM[iC];
    double y    = elemM[iH];
    double z    = elemM[iN];
    double beta = elemM[gas->elementIndex("O")] * 0.5;        // moles of O as O2
 
    //--------------------------------------------------------------------------
 
    if(mixf < mixfStoic) {                        // lean: burn all fuel, leftover O2
 
        for(int k=0; k<nspc; k++)
            ypcc[k] = 0.0;
 
        if(iCO2 > 0)                              // CO2 is not in the H2 mechs
            ypcc[iCO2] = x;
        ypcc[iH2O] = y*0.5;
        ypcc[iN2]  = z*0.5;
        ypcc[iO2]  = beta - (x+y/4.0);
 
    }
    else{                                         // rich: burn all O2, leftover fuel
 
        //double eta = beta/(x+y/4.0); // extent of reaction
        double eta = (beta-nOnotFromO2/2)/(x+y/4-nOnotFromO2/2); // extent of reaction
        if(eta > 1.0)
            cout << endl << "eta > 1.0" << endl;
        d1 = 1.0-eta;                            // fraction of fuel unburnt
 
        for(int k=0; k<nspc; k++)
            ypcc[k] *= d1;                       // initialize everything then correct
        if(iCO2 > 0)                             // CO2 is not in the H2 mechs
            ypcc[iCO2] = (x-nCnotFromCO2) + nCnotFromCO2*eta;       // init + formed
        ypcc[iH2O] = (y-nHnotFromH2O)*0.5 + nHnotFromH2O*0.5*eta;   // init + formed
        ypcc[iN2]  = z*0.5;
        ypcc[iO2]  = 0.0;
 
    }
 
    //--------------------------------------------------------------------------
 
    double sum = 0.0;                       // normalize moles
    for(int k=0; k<nspc; k++)
        sum += ypcc[k];
    assert(sum != 0.0);
    for(int k=0; k<nspc; k++)
        ypcc[k] /= sum;
    gas->setMoleFractions( &ypcc[0] );      // set mole fractions
    gas->getMassFractions( &ypcc[0] );      // set ypcc as mass fractions
 
 
    //--------------------------------------------------------------------------
 
    gas->setState_HP(hpcc, P, 1.E-10);    // get temperature as Tadiabatic
    //gas->setState_HP(hpcc, P);    // get temperature as Tadiabatic
    Tpcc = gas->temperature();
 
}
 
 
void streams::setStoicMixf() {
 
    vector<double> x_c(nspc,0.0);
 
    double mc0, mc1, mo0, mo1, mh0, mh1;
 
    vector<double> elemMassFrac0 = getElementMassFracs(&y0[0]);
    vector<double> elemMassFrac1 = getElementMassFracs(&y1[0]);
 
    mc0 = elemMassFrac0[gas->elementIndex("C")]/
        gas->atomicWeight(gas->elementIndex("C"));
    mc1 = elemMassFrac1[gas->elementIndex("C")]/
        gas->atomicWeight(gas->elementIndex("C"));
    mh0 = elemMassFrac0[gas->elementIndex("H")]/
        gas->atomicWeight(gas->elementIndex("H"));
    mh1 = elemMassFrac1[gas->elementIndex("H")]/
        gas->atomicWeight(gas->elementIndex("H"));
    mo0 = elemMassFrac0[gas->elementIndex("O")]/
        gas->atomicWeight(gas->elementIndex("O"));
    mo1 = elemMassFrac1[gas->elementIndex("O")]/
        gas->atomicWeight(gas->elementIndex("O"));
 
    mixfStoic = (2.0*mc0 + 0.5*mh0 - mo0) /
        (mo1-mo0 + 2.0*(mc0-mc1) + 0.5*(mh0-mh1));
 
    cout << endl << "# mixfStoic = m_fuel/(m_fuel+m_air) = " << mixfStoic << endl;
 
}
 
 
vector<double> streams::getElementMassFracs(const double *y) {
 
    vector<double> elemMassFrac(gas->nElements(), 0.0);
    gas->setMassFractions( &y[0] );
    for(int k=0; k<gas->nElements(); k++)
        elemMassFrac[k] = gas->elementalMassFraction(k);
    return elemMassFrac;
}
 
 
vector<double> streams::getElementMoles(const double *x,
                                        double &nOnotFromO2,
                                        double &nHnotFromH2O,
                                        double &nCnotFromCO2) {
 
    nOnotFromO2  = 0.0;
    nHnotFromH2O = 0.0;
    nCnotFromCO2 = 0.0;
 
    //vector<double> atomArr(gas->nElements());
    vector<double> elemM(gas->nElements(), 0.0);
    int iO2  = gas->speciesIndex("O2");
    int iO   = gas->elementIndex("O");
    int iCO2 = gas->speciesIndex("CO2");
    int iC   = gas->elementIndex("C");
    int iH2O = gas->speciesIndex("H2O");
    int iH   = gas->elementIndex("H");
 
    for(int k=0; k<nspc; k++) {
        vector<double> xx(gas->nSpecies(), 0.0);
        xx[k] = 1.0;
        gas->setMoleFractions(&xx[0]);
        for(int m=0; m<gas->nElements(); m++)
            elemM[m] += x[k] * gas->elementalMoleFraction(m);
        if(k != iO2)  nOnotFromO2  += gas->elementalMoleFraction(iO) * x[k];
        if(k != iCO2) nCnotFromCO2 += gas->elementalMoleFraction(iC) * x[k];
        if(k != iH2O) nHnotFromH2O += gas->elementalMoleFraction(iH) * x[k];
    }
    return elemM;
}
 
 
double streams::getMixtureFraction(const double *y, const bool doBeta01) {
 
    vector<double> elemMF;
    double         beta;
 
    elemMF = getElementMassFracs(y);
    beta   = gCHON[0] * elemMF[gas->elementIndex("C")] +
             gCHON[1] * elemMF[gas->elementIndex("H")] +
             gCHON[2] * elemMF[gas->elementIndex("O")] +
             gCHON[3] * elemMF[gas->elementIndex("N")];
 
    if(doBeta01) {
        elemMF = getElementMassFracs(&y0[0]);
        beta0  = gCHON[0] * elemMF[gas->elementIndex("C")] +
                 gCHON[1] * elemMF[gas->elementIndex("H")] +
                 gCHON[2] * elemMF[gas->elementIndex("O")] +
                 gCHON[3] * elemMF[gas->elementIndex("N")];
 
        elemMF = getElementMassFracs(&y1[0]);
        beta1  = gCHON[0] * elemMF[gas->elementIndex("C")] +
                 gCHON[1] * elemMF[gas->elementIndex("H")] +
                 gCHON[2] * elemMF[gas->elementIndex("O")] +
                 gCHON[3] * elemMF[gas->elementIndex("N")];
    }
 
    if( beta1-beta0 == 0 )
        return 0.0;            // to avoid division by zero
    else
        return (beta - beta0)/(beta1 - beta0);
}