Source code
C++ source code for process model
Qmicroturbine.cc
—
C++ source code,
14 kB (15267 bytes)
File contents
/** @file Qmicroturbine.cc
@brief Main file
This file is part of LIBPF
All rights reserved; do not distribute without permission.
@author (C) Copyright 2009-2010 Paolo Greppi 3iP
Repository revision: $Rev: 551 $
File revision date: $Date: 2010-05-15 12:53:03 +0200 (sab, 15 mag 2010) $
*/
static const int verbositySubSystem = 0;
#include "models.h"
#include "stream.h"
#include "graph.h"
#include "usermodels.h"
#include "purecomps.h"
#include <boost/assign/list_of.hpp>
#include <boost/assign/list_inserter.hpp>
void setComponents(void) {
components.addcomp(new purecomps::water);
components.addcomp(new purecomps::N2);
components.addcomp(new purecomps::O2);
components.addcomp(new purecomps::methane);
components.addcomp(new purecomps::CO);
components.addcomp(new purecomps::CO2);
components.addcomp(new purecomps::H2);
} // void setComponents(void)
// /////////////////////////////////////////////////////////////////////////////
// ///////////////////////////////// microturbo ////////////////////////////////
// /////////////////////////////////////////////////////////////////////////////
#define standalone
#ifdef standalone
#define AAA flowsheet<zero_zero>
#else
#define AAA flowsheet<vertex<3,2> >
#endif
class microturbo : public modelBase, public AAA {
private:
static std::string type_;
void maketables(long);
public:
microturbo(const std::string &t, const std::string &d) : persistent_item_interface(-1), persistent_item(-1, t, d), modelBaseInterface(-1), modelBase(t, d) {
maketables(-1);
};
microturbo(long cid) : persistent_item_interface(cid), persistent_item(cid, "", ""), modelBaseInterface(cid), modelBase(cid), AAA(cid) {
maketables(cid);
readtables(cid);
}
Qdouble LHV_CH4; ///< methane Lower Heating Value
Qdouble LHV_H2; ///< hydrogen Lower Heating Value
Qdouble LHV_CO; ///< carbon monoxide Lower Heating Value
Qdouble powerIn; ///< Lower Heating Value input
Qdouble power; ///< Gas Turbine Subsystem power output
Qdouble powerEl; ///< Total net electric output
Qdouble etaEl; ///< Overall electric efficiency
Qdouble beta0; ///< compressor compression ratio (cr) at the nominal point
Qdouble mC0; ///< compressor mass flow at the nominal point (S01)
Qdouble TC0; ///< compressor inlet temperature at the nominal point (S01)
Qdouble PC0; ///< compressor inlet pressure at the nominal point (S01)
Qdouble mrC0; ///< compressor corrected mass flow at the nominal point
Qdouble mT0; ///< turbine mass flow at the nominal point (S05)
Qdouble PT0; ///< turbine inlet pressure at the nominal point (S05)
Qdouble mrT0; ///< turbine corrected mass flow at the nominal point
Qdouble TIT; ///< should be Turbine Inlet Temperature
Qdouble mC; ///< should-be compressor uncorrected mass flow
Qdouble mT; ///< should-be turbine uncorrected mass flow
Qdouble maC; ///< actual compressor uncorrected mass flow
Qdouble maT; ///< actual turbine uncorrected mass flow
Qdouble nr; ///< compressor/turbine reduced frequency
void makeuserassembly(std::list<assignment *>::iterator &p);
void setup(void);
void calculate(int MAXITS = 100, int level = 0);
const std::string &type(void) const {
return type_;
}
int smiter_precalc(void) { return 5; }
bool supports_simultaneous(void) { return true; };
int maxit(void) { return 200; }
}; // microturbo
std::string microturbo::type_("microturbo");
void microturbo::maketables(long cid) {
static const int verbosityLocal = 0;
try {
diagnostic(2, " entered for " << tag());
// Non embedded objects
MAKEQDOUBLE(LHV_CH4, 802618.0, "kJ/kmol", "methane Lower Heating Value");
MAKEQDOUBLE(LHV_H2, 241814.0, "kJ/kmol", "hydrogen Lower Heating Value");
MAKEQDOUBLE(LHV_CO, 282980.0, "kJ/kmol", "carbon monoxide Lower Heating Value");
MAKEQDOUBLE(powerIn, 0.0, "W", "Lower Heating Value input");
MAKEQDOUBLE(power, 0.0, "W", "Gas Turbine Subsystem power output");
MAKEQDOUBLE(powerEl, 0.0, "W", "Total net electric output");
MAKEQDOUBLE(etaEl, 0.0, "", "Overall electric efficiency");
MAKEQDOUBLE(beta0, 3.6, "", "compressor compression ratio at the nominal point");
MAKEQDOUBLE(mC0, 1.66667, "kg/s", "compressor mass flow at the nominal point");
MAKEQDOUBLE(TC0, 296.15, "K", "compressor inlet temperature at the nominal point");
MAKEQDOUBLE(PC0, 101325., "Pa", "compressor inlet pressure at the nominal point");
MAKEQDOUBLE(mrC0, 0.0, "m*s", "compressor corrected mass flow at the nominal point");
MAKEQDOUBLE(mT0, 1.87407, "kg/s", "turbine mass flow at the nominal point");
MAKEQDOUBLE(PT0, 347676., "Pa", "turbine inlet pressure at the nominal point");
MAKEQDOUBLE(mrT0, 0.0, "m*s", "turbine corrected mass flow at the nominal point");
MAKEQDOUBLE(TIT, t0, "K", "should be Turbine Inlet Temperature");
MAKEQDOUBLE(mC, 0.0, "kg/s", "should-be compressor uncorrected mass flow");
MAKEQDOUBLE(mT, 0.0, "kg/s", "should-be turbine uncorrected mass flow");
MAKEQDOUBLE(maC, 0.0, "kg/s", "actual compressor uncorrected mass flow");
MAKEQDOUBLE(maT, 0.0, "kg/s", "actual turbine uncorrected mass flow");
MAKEQDOUBLE(nr, 0.0, "", "compressor and turbine reduced frequency");
if (cid == -1) {
diagnostic(2, "Define unit operations");
makeVertex<mixer>("MIX", "Mixes fuel preheated air and FCS purge", cid);
makeVertex<splitter>("SPLIT", "Compressor outlet splitter", cid, boost::assign::map_list_of<std::string, long>("nStreams", 2));
makeVertex<compressor>("T", "Turbine", cid);
makeVertex<compressor>("C", "Compressor", cid);
std::list<std::string> rHC = boost::assign::list_of("reactionCH4_TC")
("reactionCO_TC")
("reactionH2_TC");
makeVertex<genflashN1>("RX", "Burner", rHC, cid);
makeVertex<hx>("HX", "Turbine recuperator", cid);
diagnostic(2, "Define stream objects and connect");
makeEdge<stream_V>("S01", "Air to GTS", "source", "out", "C", "in", cid);
makeEdge<stream_V>("S02", "NG to GTS burner", "source", "out", "RX", "in", cid);
makeEdge<stream_V>("S03", "Hot air recovery", "source", "out", "MIX", "in", cid);
makeEdge<stream_V>("S04", "RX burner inlet", "MIX", "out", "RX", "in", cid);
makeEdge<stream_V>("S05", "Turbine Inlet", "RX", "out", "T", "in", cid);
makeEdge<stream_V>("S06", "Turbine Outlet", "T", "out", "HX", "hotin", cid);
makeEdge<stream_V>("S07", "Compressed Air", "C", "out", "SPLIT", "in", cid);
makeEdge<stream_V>("S08", "HX cold inlet", "SPLIT", "out2", "HX", "coldin", cid);
makeEdge<stream_V>("S09", "HX cold outlet", "HX", "coldout", "MIX", "in", cid);
makeEdge<stream_V>("S10", "Air to FCS", "SPLIT", "out1", "sink", "in", cid);
makeEdge<stream_V>("S11", "Recuperator exhaust", "HX", "hotout", "sink", "in", cid);
}
} // try
catch (error &e) {
e.append(CURRENT_FUNCTION);
throw;
} // catch
} // microturbo::maketables
void microturbo::setup(void) {
static const int verbosityLocal = 0;
try {
diagnostic(2, " entered for " << tag());
double Pset = 4.0; // max pressure in bar
diagnostic(3, "Calling flowsheet::setup to initialize any embedded flowsheet");
AAA::setup();
diagnostic(3, "Setting input variables");
beta0 = 3.94769306686405;
mC0 = Qdouble(0.306891411, "kg/s");
TC0 = Qdouble(288.15, "K");
PC0 = Qdouble(101325., "Pa");
mT0 = Qdouble(0.309236235, "kg/s");
PT0 = Qdouble(385000., "Pa");
TIT = Qdouble(1144.0, "K");
mrC0 = mC0 * sqrt(TC0 / T0) / PC0;
mrT0 = mT0 * sqrt(TIT / T0) / PT0;
#ifdef standalone
// fresh air
O("S01")->Q("T")->set(15.0 + 273.15, "K");
O("S01")->Q("P")->set(1.0, "atm");
O("S01")->S("flowoption")->set("Nx");
my_cast<stream *>(O("S01"), CURRENT_FUNCTION)->clearcomposition();
O("S01")->O("Tphase")->Q("ndot")->set(38.49159293, "kmol/h");
O("S01")->O("Tphase")->Q("x", "N2")->set(0.792);
O("S01")->O("Tphase")->Q("x", "O2")->set(0.198);
O("S01")->O("Tphase")->Q("x", "H2O")->set(0.01);
// feed natural gas to turbine
O("S02")->Q("T")->set(15.0 + 273.15, "K");
O("S02")->Q("P")->set(Pset, "bar");
O("S02")->S("flowoption")->set("Nx");
my_cast<stream *>(O("S02"), CURRENT_FUNCTION)->clearcomposition();
O("S02")->O("Tphase")->Q("ndot")->set(0.522280634, "kmol/h");
O("S02")->O("Tphase")->Q("x", "CH4")->set(0.99);
O("S02")->O("Tphase")->Q("x", "N2")->set(0.01);
// hot recovered air
O("S03")->Q("T")->set(15.0 + 273.15, "K");
O("S03")->Q("P")->set(1.0, "atm");
O("S03")->S("flowoption")->set("Nx");
my_cast<stream *>(O("S03"), CURRENT_FUNCTION)->clearcomposition();
O("S03")->O("Tphase")->Q("ndot")->set(0.0, "kmol/h");
O("S03")->O("Tphase")->Q("x", "N2")->set(0.792);
O("S03")->O("Tphase")->Q("x", "O2")->set(0.198);
O("S03")->O("Tphase")->Q("x", "H2O")->set(0.01);
#endif // standalone
// Burner
O("RX")->Q("deltaP")->set(100.0, "mbar");
O("RX")->Q("duty")->set(0.0, "W");
O("RX")->S("option")->set("DH");
// total combustion of methane
O("RX")->O("reactions", 0)->Q("z")->set(1.0);
// total combustion of CO
O("RX")->O("reactions", 1)->Q("z")->set(1.0);
// total combustion of H2
O("RX")->O("reactions", 2)->Q("z")->set(1.0);
// Compressor outlet splitter, 1-(stream16/stream15)
O("SPLIT")->Q("outSplit", 0)->set(0.0);
O("SPLIT")->Q("outSplit", 1)->set(1.0);
// Turbine recuperator
O("HX")->Q("deltaPhot")->set(50.0, "mbar");
O("HX")->S("hotoption")->set("D");
O("HX")->Q("deltaPcold")->set(50.0, "mbar");
O("HX")->S("coldoption")->set("D");
// O("HX")->S("option")->set("Thot");
// O("HX")->Q("hotT")->set(324.9144+273.15, "K");
O("HX")->S("option")->set("UA");
O("HX")->Q("U")->set(100.0, "W/(m2*K)");
O("HX")->Q("A")->set(14.2866994561544, "m2");
O("C")->S("option")->set("P");
O("C")->Q("P")->set(Pset, "bar");
O("C")->Q("theta")->set(0.85);
O("C")->Q("etaM")->set(0.99);
O("C")->Q("etaE")->set(1.0);
O("T")->S("option")->set("P");
O("T")->Q("P")->set(1013.25 + 50.0, "mbar");
O("T")->Q("theta")->set(0.719017307023255);
O("T")->Q("etaM")->set(0.99);
O("T")->Q("etaE")->set(1.0);
diagnostic(3, "Initializing cut streams");
// Turbine recuperator (HX) cold outlet, streamB52.17
O("S09")->Q("T")->set(540.2536 + 273.15, "K");
O("S09")->Q("P")->set(Pset, "bar");
O("S09")->S("flowoption")->set("Nx");
my_cast<stream *>(O("S09"), CURRENT_FUNCTION)->clearcomposition();
O("S09")->O("Tphase")->Q("ndot")->set(63.69577 / 2.0, "kmol/h");
O("S09")->O("Tphase")->Q("x", "N2")->set(0.79);
O("S09")->O("Tphase")->Q("x", "O2")->set(0.21);
diagnostic(3, "Defining cut streams");
flowsheetBase::addcut("S09");
diagnostic(3, "Making selected outputs visible in GUI");
powerIn.setOutput();
powerEl.setOutput();
etaEl.setOutput();
TIT.setOutput();
mC.setOutput();
mT.setOutput();
maC.setOutput();
maT.setOutput();
nr.setOutput();
O("RX")->Q("T")->setOutput();
O("RX")->Q("P")->setOutput();
O("RX")->Q("duty")->setOutput();
O("RX")->Q("deltaP")->setOutput();
O("RX")->O("reactions", 0)->Q("conv")->setOutput();
O("RX")->O("reactions", 0)->Q("z")->setOutput();
O("RX")->O("reactions", 1)->Q("conv")->setOutput();
O("RX")->O("reactions", 1)->Q("z")->setOutput();
O("RX")->O("reactions", 2)->Q("conv")->setOutput();
O("RX")->O("reactions", 2)->Q("z")->setOutput();
} // try
catch (error &e) {
e.append(CURRENT_FUNCTION);
throw;
} // catch
} // microturbo::setup
void microturbo::calculate(int MAXITS, int level) {
static const int verbosityLocal = 0;
try {
diagnostic(2, " entered for " << tag());
reset_errors();
AAA::calculate(MAXITS, level);
powerIn = *O("S02:Tphase")->Q("ndotcomps", "CH4") * LHV_CH4 +
*O("S02:Tphase")->Q("ndotcomps", "H2") * LHV_H2 +
*O("S02:Tphase")->Q("ndotcomps", "CO") * LHV_CO;
powerEl = -*O("C")->Q("We") - *O("T")->Q("We");
etaEl = powerEl / powerIn;
} // try
catch (error &e) {
e.append(CURRENT_FUNCTION);
throw;
} // catch
} // microturbo::calculate
void microturbo::makeuserassembly(std::list<assignment *>::iterator &p) {
static const int verbosityLocal = -1;
try {
diagnostic(2, " entered for " << tag());
long i(0);
diagnostic(3, "microturbo flowsheet specific equations");
// turbine matching
nr = sqrt((*O("C")->Q("cr") - One) / (beta0 - One));
mC = mrC0 * nr * *O("S01")->Q("P") / sqrt(*O("S01")->Q("T") / T0);
// mrT0 = mT0*sqrt(TIT/T0)/PT0;
mT = mrT0 * nr * *O("S05")->Q("P") / sqrt(*O("S05")->Q("T") / T0);
maC = *O("S01:Tphase")->Q("mdot");
maT = *O("S05:Tphase")->Q("mdot");
MakeAssignment(*Q("T.P"), *Q("T.P") - (*Q("S11.P") - Qdouble(101325.0, "Pa")), "set discharge pressure");
MakeAssignment(*O("S01:Tphase")->Q("ndot"), mC / *O("S01:Tphase")->Q("AMW"),
"Adjust inlet air to balance GTS");
MakeAssignment(*O("S02:Tphase")->Q("ndot"), *O("S02:Tphase")->Q("ndot") +
(TIT - *O("S05")->Q("T")) * *O("S05:Tphase")->Q("cp") * *O("S05:Tphase")->Q("mdot") / LHV_CH4,
"Adjust NG to GTS burner to fix TIT");
// linear
// Qdouble factor = One - (One - nr)/3.;
// quadratic
Qdouble factor = One - (nr - One) * (nr - One) / 1.5;
// cubic
// Qdouble factor = One + pow(nr-One, 3.)/0.27;
MakeAssignment(*Q("C.theta"), 0.85 * factor, "compressor efficiency");
// MakeAssignment(*Q("T.theta"), 0.719017307023255*pow(factor, 0.1), "turbine efficiency");
/*
// Enable these to match vendor data
MakeAssignment(*Q("HX.A"), *Q("HX.A") * Qdouble(275.+273.15, "K") / *Q("S11.T"),
"Match vendor claimed exhaust temperature");
MakeAssignment(*Q("T.theta"), *Q("T.theta") * sqrt(*Q("S05.T") / Qdouble(1144., "K")),
"Match vendor claimed TIT");
*/
} // try
catch (error &e) {
e.append(CURRENT_FUNCTION);
throw;
} // catch
} // microturbo::makeuserassembly
const char *default_type = "microturbo";
// /////////////////////////////////////////////////////////////////////////////
// //////////////////////////// registrar_mymodels /////////////////////////////
// /////////////////////////////////////////////////////////////////////////////
class registrar_mymodels : public registrar {
private:
static bool done_;
public:
registrar_mymodels(void);
}; // class registrar_mymodels
registrar_mymodels::registrar_mymodels(void) {
static const int verbosityLocal = 0;
diagnostic(2, "Entered");
if (!done_) {
diagnostic(3, "Registering");
model_factory().Register<microturbo>("microturbo");
done_ = true;
}
} // registrar_mymodels::registrar_mymodels
bool registrar_mymodels::done_(false);
static registrar_mymodels r_;
#include "Q.cc"
