changes by Karl Meredith to implement adiabatic, freely-propagating flames
This commit is contained in:
Dave Goodwin
parent
71caff1d2b
commit
f18a811cf2
6 changed files with 416 additions and 140 deletions
|
|
@ -90,7 +90,7 @@ namespace Cantera {
|
|||
* Initialize. Base class method does nothing, but may be
|
||||
* overloaded.
|
||||
*/
|
||||
virtual void init(){}
|
||||
virtual void init(){ }
|
||||
|
||||
virtual void setInitialState(doublereal* xlocal = 0){}
|
||||
virtual void setState(int point, const doublereal* state, doublereal* x) {}
|
||||
|
|
@ -186,6 +186,34 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
|
||||
//added by Karl Meredith
|
||||
void setTolerances(doublereal rtol, doublereal atol,int ts=0) {
|
||||
for (int n=0;n<m_nv;n++){
|
||||
if(ts>=0) {
|
||||
m_rtol_ss[n] = rtol;
|
||||
m_atol_ss[n] = atol;
|
||||
}
|
||||
if (ts <= 0) {
|
||||
m_rtol_ts[n] = rtol;
|
||||
m_atol_ts[n] = atol;
|
||||
}
|
||||
}
|
||||
}
|
||||
//added by Karl Meredith
|
||||
void setTolerancesTS(doublereal rtol, doublereal atol) {
|
||||
for (int n=0;n<m_nv;n++){
|
||||
m_rtol_ts[n] = rtol;
|
||||
m_atol_ts[n] = atol;
|
||||
}
|
||||
}
|
||||
//added by Karl Meredith
|
||||
void setTolerancesSS(doublereal rtol, doublereal atol) {
|
||||
for (int n=0;n<m_nv;n++){
|
||||
m_rtol_ss[n] = rtol;
|
||||
m_atol_ss[n] = atol;
|
||||
}
|
||||
}
|
||||
|
||||
/// Relative tolerance of the nth component.
|
||||
doublereal rtol(int n) { return (m_rdt == 0.0 ? m_rtol_ss[n] : m_rtol_ts[n]); }
|
||||
|
||||
|
|
@ -426,7 +454,11 @@ namespace Cantera {
|
|||
// "base class method _finalize called!");
|
||||
//}
|
||||
|
||||
//added by Karl Meredith
|
||||
doublereal m_zfixed;
|
||||
doublereal m_tfixed;
|
||||
|
||||
bool m_adiabatic;
|
||||
protected:
|
||||
|
||||
doublereal m_rdt;
|
||||
|
|
@ -448,6 +480,8 @@ namespace Cantera {
|
|||
string m_id, m_desc;
|
||||
Refiner* m_refiner;
|
||||
|
||||
|
||||
|
||||
private:
|
||||
|
||||
};
|
||||
|
|
|
|||
|
|
@ -33,6 +33,7 @@ namespace Cantera {
|
|||
m_xnew.resize(size(), 0.0);
|
||||
for (int n = 0; n < m_nd; n++) {
|
||||
domain(n)._getInitialSoln(m_x.begin() + start(n));
|
||||
domain(n).m_adiabatic=false;
|
||||
}
|
||||
|
||||
// set some defaults
|
||||
|
|
@ -45,6 +46,21 @@ namespace Cantera {
|
|||
|
||||
}
|
||||
|
||||
|
||||
// added by Karl Meredith
|
||||
void Sim1D::setInitialGuess(string component, vector_fp& locs, vector_fp& vals){
|
||||
|
||||
for (int dom=0;dom<m_nd;dom++){
|
||||
Domain1D& d = domain(dom);
|
||||
int ncomp=d.nComponents();
|
||||
for (int comp=0;comp<ncomp;comp++){
|
||||
if(d.componentName(comp)==component){
|
||||
setProfile(dom,comp,locs,vals);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Set a single value in the solution vector.
|
||||
|
|
@ -396,6 +412,126 @@ namespace Cantera {
|
|||
}
|
||||
|
||||
|
||||
/**
|
||||
* Add node for fixed temperature point of freely propagating flame
|
||||
*/
|
||||
//added by Karl Meredith
|
||||
int Sim1D::setFixedTemperature(doublereal t) {
|
||||
int np = 0;
|
||||
vector_fp znew, xnew;
|
||||
doublereal xmid;
|
||||
doublereal zfixed,interp_factor;
|
||||
doublereal z1,z2,t1,t2;
|
||||
int strt, n, m, i;
|
||||
int m1,m2;
|
||||
vector_int dsize;
|
||||
bool addnewpt=false;
|
||||
|
||||
|
||||
for (n = 0; n < m_nd; n++) {
|
||||
strt = znew.size();
|
||||
Domain1D& d = domain(n);
|
||||
|
||||
int comp = d.nComponents();
|
||||
|
||||
// loop over points in the current grid to determine where new point is needed.
|
||||
int npnow = d.nPoints();
|
||||
int nstart = znew.size();
|
||||
for (m = 0; m < npnow-1; m++) {
|
||||
cout << "T["<<m<<"]="<<value(n,2,m)<<endl;
|
||||
if(value(n,2,m)==t){
|
||||
zfixed=d.grid(m);
|
||||
//set d.zfixed, d.ztemp
|
||||
d.m_zfixed=zfixed;
|
||||
d.m_tfixed=t;
|
||||
cout << "T already fixed at "<<d.grid(m)<<endl;
|
||||
addnewpt=false;
|
||||
break;
|
||||
}
|
||||
else if((value(n,2,m)<t) && (value(n,2,m+1)>t)){
|
||||
cout << "T in between "<<value(n,2,m)<<" and "<<value(n,2,m+1)<<endl;
|
||||
z1=d.grid(m);
|
||||
m1=m;
|
||||
m2=m+1;
|
||||
z2=d.grid(m+1);
|
||||
t1=value(n,2,m);
|
||||
t2=value(n,2,m+1);
|
||||
|
||||
zfixed=(z1-z2)/(t1-t2)*(t-t2)+z2;
|
||||
cout << zfixed<<endl;
|
||||
//set d.zfixed, d.ztemp;
|
||||
d.m_zfixed=zfixed;
|
||||
d.m_tfixed=t;
|
||||
addnewpt=true;
|
||||
break;
|
||||
//copy solution domain and push back values
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
for (m = 0; m < npnow; m++) {
|
||||
// add the current grid point to the new grid
|
||||
znew.push_back(d.grid(m));
|
||||
|
||||
// do the same for the solution at this point
|
||||
for (i = 0; i < comp; i++) {
|
||||
xnew.push_back(value(n, i, m));
|
||||
}
|
||||
if(m==m1&&addnewpt){
|
||||
//add new point at zfixed
|
||||
znew.push_back(zfixed);
|
||||
np++;
|
||||
interp_factor=(zfixed-z2)/(z1-z2);
|
||||
// for each component, linearly interpolate
|
||||
// the solution to this point
|
||||
for (i = 0; i < comp; i++) {
|
||||
xmid = interp_factor*(value(n, i, m) - value(n, i, m+1))+value(n,i,m+1);
|
||||
xnew.push_back(xmid);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
dsize.push_back(znew.size() - nstart);
|
||||
}
|
||||
|
||||
// At this point, the new grid znew and the new solution
|
||||
// vector xnew have been constructed, but the domains
|
||||
// themselves have not yet been modified. Now update each
|
||||
// domain with the new grid.
|
||||
|
||||
int gridstart = 0, gridsize;
|
||||
for (n = 0; n < m_nd; n++) {
|
||||
Domain1D& d = domain(n);
|
||||
// Refiner& r = d.refiner();
|
||||
gridsize = dsize[n]; // d.nPoints() + r.nNewPoints();
|
||||
d.setupGrid(gridsize, znew.begin() + gridstart);
|
||||
gridstart += gridsize;
|
||||
}
|
||||
|
||||
// Replace the current solution vector with the new one
|
||||
m_x.resize(xnew.size());
|
||||
copy(xnew.begin(), xnew.end(), m_x.begin());
|
||||
|
||||
// resize the work array
|
||||
m_xnew.resize(xnew.size());
|
||||
|
||||
copy(xnew.begin(), xnew.end(), m_xnew.begin());
|
||||
|
||||
resize();
|
||||
finalize();
|
||||
return np;
|
||||
}
|
||||
|
||||
//added by Karl Meredith
|
||||
void Sim1D::setAdiabaticFlame(void){
|
||||
int n;
|
||||
for (n = 0; n < m_nd; n++) {
|
||||
Domain1D& d = domain(n);
|
||||
d.m_adiabatic=true;
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Set grid refinement criteria. If dom >= 0, then the settings
|
||||
* apply only to the specified domain. If dom < 0, the settings
|
||||
|
|
|
|||
|
|
@ -1,127 +1,135 @@
|
|||
/**
|
||||
* @file Sim1D.h
|
||||
*/
|
||||
|
||||
#ifndef CT_SIM1D_H
|
||||
#define CT_SIM1D_H
|
||||
|
||||
#include "OneDim.h"
|
||||
#include "../funcs.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* One-dimensional simulations. Class Sim1D extends class OneDim
|
||||
* by storing the solution vector, and by adding a hybrid
|
||||
* Newton/time-stepping solver.
|
||||
*/
|
||||
class Sim1D : public OneDim {
|
||||
|
||||
public:
|
||||
|
||||
/**
|
||||
* Default constructor. This constructor is provided to make
|
||||
* the class default-constructible, but is not meant to be
|
||||
* used in most applications. Use the next constructor
|
||||
* instead.
|
||||
*/
|
||||
Sim1D();
|
||||
|
||||
|
||||
/**
|
||||
* Standard constructor.
|
||||
* @param domains A vector of pointers to the domains to be linked together.
|
||||
* The domain pointers must be entered in left-to-right order --- i.e.,
|
||||
* the pointer to the leftmost domain is domain[0], the pointer to the
|
||||
* domain to its right is domain[1], etc.
|
||||
*/
|
||||
Sim1D(vector<Domain1D*>& domains);
|
||||
|
||||
/// Destructor. Does nothing.
|
||||
virtual ~Sim1D(){}
|
||||
|
||||
|
||||
/**
|
||||
* @name Setting initial values
|
||||
*
|
||||
* These methods are used to set the initial values of
|
||||
* solution components.
|
||||
*/
|
||||
//@{
|
||||
|
||||
/// Set one entry in the solution vector.
|
||||
void setValue(int dom, int comp, int localPoint, doublereal value);
|
||||
|
||||
/// Get one entry in the solution vector.
|
||||
doublereal value(int dom, int comp, int localPoint) const;
|
||||
|
||||
doublereal workValue(int dom, int comp, int localPoint) const;
|
||||
|
||||
/// Specify a profile for one component of one domain.
|
||||
void setProfile(int dom, int comp, const vector_fp& pos,
|
||||
const vector_fp& values);
|
||||
|
||||
/// Set component 'comp' of domain 'dom' to value 'v' at all points.
|
||||
void setFlatProfile(int dom, int comp, doublereal v);
|
||||
|
||||
//@}
|
||||
|
||||
void save(string fname, string id, string desc);
|
||||
|
||||
/// Print to stream s the current solution for all domains.
|
||||
void showSolution(ostream& s);
|
||||
void showSolution();
|
||||
|
||||
const doublereal* solution() { return m_x.begin(); }
|
||||
|
||||
void setTimeStep(doublereal stepsize, int n, integer* tsteps);
|
||||
|
||||
//void setMaxTimeStep(doublereal tmax) { m_maxtimestep = tmax; }
|
||||
|
||||
void solve(int loglevel = 0, bool refine_grid = true);
|
||||
|
||||
void eval(doublereal rdt=-1.0, int count = 1) {
|
||||
OneDim::eval(-1, m_x.begin(), m_xnew.begin(), rdt, count);
|
||||
}
|
||||
|
||||
/// Refine the grid in all domains.
|
||||
int refine(int loglevel=0);
|
||||
|
||||
/// Set the criteria for grid refinement.
|
||||
void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
|
||||
doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
|
||||
void setMaxGridPoints(int dom = -1, int npoints = 300);
|
||||
|
||||
void restore(string fname, string id);
|
||||
void getInitialSoln();
|
||||
|
||||
void setSolution(const doublereal* soln) {
|
||||
copy(soln, soln + m_x.size(), m_x.begin());
|
||||
}
|
||||
|
||||
const doublereal* solution() const { return m_x.begin(); }
|
||||
|
||||
protected:
|
||||
|
||||
vector_fp m_x; // the solution vector
|
||||
vector_fp m_xnew; // a work array used to hold the residual
|
||||
// or the new solution
|
||||
doublereal m_tstep; // timestep
|
||||
vector_int m_steps; // array of number of steps to take before
|
||||
// re-attempting the steady-state solution
|
||||
|
||||
|
||||
private:
|
||||
|
||||
/// Calls method _finalize in each domain.
|
||||
void finalize();
|
||||
|
||||
void newtonSolve(int loglevel);
|
||||
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
/**
|
||||
* @file Sim1D.h
|
||||
*/
|
||||
|
||||
#ifndef CT_SIM1D_H
|
||||
#define CT_SIM1D_H
|
||||
|
||||
#include "OneDim.h"
|
||||
#include "../funcs.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
||||
/**
|
||||
* One-dimensional simulations. Class Sim1D extends class OneDim
|
||||
* by storing the solution vector, and by adding a hybrid
|
||||
* Newton/time-stepping solver.
|
||||
*/
|
||||
class Sim1D : public OneDim {
|
||||
|
||||
public:
|
||||
|
||||
/**
|
||||
* Default constructor. This constructor is provided to make
|
||||
* the class default-constructible, but is not meant to be
|
||||
* used in most applications. Use the next constructor
|
||||
* instead.
|
||||
*/
|
||||
Sim1D();
|
||||
|
||||
|
||||
/**
|
||||
* Standard constructor.
|
||||
* @param domains A vector of pointers to the domains to be linked together.
|
||||
* The domain pointers must be entered in left-to-right order --- i.e.,
|
||||
* the pointer to the leftmost domain is domain[0], the pointer to the
|
||||
* domain to its right is domain[1], etc.
|
||||
*/
|
||||
Sim1D(vector<Domain1D*>& domains);
|
||||
|
||||
/// Destructor. Does nothing.
|
||||
virtual ~Sim1D() {}
|
||||
|
||||
/**
|
||||
* @name Setting initial values
|
||||
*
|
||||
* These methods are used to set the initial values of
|
||||
* solution components.
|
||||
*/
|
||||
//@{
|
||||
|
||||
/// Set initial guess based on equilibrium
|
||||
//added by Karl Meredith
|
||||
void setInitialGuess(string component, vector_fp& locs, vector_fp& vals);
|
||||
|
||||
/// Set one entry in the solution vector.
|
||||
void setValue(int dom, int comp, int localPoint, doublereal value);
|
||||
|
||||
/// Get one entry in the solution vector.
|
||||
doublereal value(int dom, int comp, int localPoint) const;
|
||||
|
||||
doublereal workValue(int dom, int comp, int localPoint) const;
|
||||
|
||||
/// Specify a profile for one component of one domain.
|
||||
void setProfile(int dom, int comp, const vector_fp& pos,
|
||||
const vector_fp& values);
|
||||
|
||||
/// Set component 'comp' of domain 'dom' to value 'v' at all points.
|
||||
void setFlatProfile(int dom, int comp, doublereal v);
|
||||
|
||||
//@}
|
||||
|
||||
void save(string fname, string id, string desc);
|
||||
|
||||
/// Print to stream s the current solution for all domains.
|
||||
void showSolution(ostream& s);
|
||||
void showSolution();
|
||||
|
||||
const doublereal* solution() { return m_x.begin(); }
|
||||
|
||||
void setTimeStep(doublereal stepsize, int n, integer* tsteps);
|
||||
|
||||
//void setMaxTimeStep(doublereal tmax) { m_maxtimestep = tmax; }
|
||||
|
||||
void solve(int loglevel = 0, bool refine_grid = true);
|
||||
|
||||
void eval(doublereal rdt=-1.0, int count = 1) {
|
||||
OneDim::eval(-1, m_x.begin(), m_xnew.begin(), rdt, count);
|
||||
}
|
||||
|
||||
/// Refine the grid in all domains.
|
||||
int refine(int loglevel=0);
|
||||
|
||||
//added by Karl Meredith
|
||||
int setFixedTemperature(doublereal t);
|
||||
//added by Karl Meredith
|
||||
void setAdiabaticFlame(void);
|
||||
|
||||
/// Set the criteria for grid refinement.
|
||||
void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
|
||||
doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
|
||||
void setMaxGridPoints(int dom = -1, int npoints = 300);
|
||||
|
||||
void restore(string fname, string id);
|
||||
void getInitialSoln();
|
||||
|
||||
void setSolution(const doublereal* soln) {
|
||||
copy(soln, soln + m_x.size(), m_x.begin());
|
||||
}
|
||||
|
||||
const doublereal* solution() const { return m_x.begin(); }
|
||||
|
||||
protected:
|
||||
|
||||
vector_fp m_x; // the solution vector
|
||||
vector_fp m_xnew; // a work array used to hold the residual
|
||||
// or the new solution
|
||||
doublereal m_tstep; // timestep
|
||||
vector_int m_steps; // array of number of steps to take before
|
||||
// re-attempting the steady-state solution
|
||||
|
||||
|
||||
private:
|
||||
|
||||
/// Calls method _finalize in each domain.
|
||||
void finalize();
|
||||
|
||||
void newtonSolve(int loglevel);
|
||||
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -241,7 +241,15 @@ namespace Cantera {
|
|||
}
|
||||
|
||||
|
||||
/* void StFlow::init() {
|
||||
|
||||
|
||||
cout << m_do_energy.begin()<< endl;
|
||||
// this->_getInitialSoln();
|
||||
cout << "Initializing StFlow\n";
|
||||
}
|
||||
|
||||
*/
|
||||
void StFlow::setupGrid(int n, const doublereal* z) {
|
||||
resize(n);
|
||||
int j;
|
||||
|
|
@ -466,6 +474,7 @@ namespace Cantera {
|
|||
}
|
||||
rsd[index(4,j)] = 1.0 - sum;
|
||||
diag[index(4,j)] = 0;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -487,9 +496,36 @@ namespace Cantera {
|
|||
// d(\rho u)/dz + 2\rho V = 0
|
||||
//
|
||||
//------------------------------------------------
|
||||
rsd[index(c_offset_U,j)] =
|
||||
-(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
|
||||
-(density(j+1)*V(x,j+1) + density(j)*V(x,j));
|
||||
|
||||
//added by Karl Meredith
|
||||
if(!m_adiabatic){
|
||||
rsd[index(c_offset_U,j)] =
|
||||
-(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
|
||||
-(density(j+1)*V(x,j+1) + density(j)*V(x,j));
|
||||
}
|
||||
else{
|
||||
//we want mdot to propagate outward from fixed T point.
|
||||
if(grid(j)>m_zfixed){
|
||||
rsd[index(c_offset_U,j)] =
|
||||
-(rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
|
||||
-(density(j+1)*V(x,j+1) + density(j)*V(x,j));
|
||||
//algebraic constraint
|
||||
diag[index(c_offset_U, j)] = 0;
|
||||
}
|
||||
else if(grid(j)==m_zfixed){
|
||||
rsd[index(c_offset_U,j)] = 0.001*(T(x,j)-m_tfixed);
|
||||
//algebraic constraint
|
||||
diag[index(c_offset_U, j)] = 0;
|
||||
}
|
||||
else if(grid(j)<m_zfixed){
|
||||
rsd[index(c_offset_U,j)] =
|
||||
-(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
|
||||
-(density(j+1)*V(x,j+1) + density(j)*V(x,j));
|
||||
//algebraic constraint
|
||||
diag[index(c_offset_U, j)] = 0;
|
||||
}
|
||||
}
|
||||
//end of 'added by Karl Meredith'
|
||||
|
||||
|
||||
//------------------------------------------------
|
||||
|
|
@ -750,6 +786,26 @@ namespace Cantera {
|
|||
}
|
||||
|
||||
|
||||
//added by Karl Meredith
|
||||
int StFlow::componentIndex(string name) const {
|
||||
|
||||
|
||||
if(name=="u") {return 0;}
|
||||
else if (name=="V") {return 1;}
|
||||
else if (name=="T") {return 2;}
|
||||
else if (name=="lambda") {return 3;}
|
||||
else {
|
||||
for (int n=4;n<m_nsp+4;n++){
|
||||
if(componentName(n)==name){
|
||||
return n;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return -1;
|
||||
}
|
||||
|
||||
|
||||
void StFlow::restore(const XML_Node& dom, doublereal* soln) {
|
||||
|
||||
vector<string> ignored;
|
||||
|
|
|
|||
|
|
@ -86,7 +86,10 @@ namespace Cantera {
|
|||
thermo_t& phase() { return *m_thermo; }
|
||||
kinetics_t& kinetics() { return *m_kin; }
|
||||
|
||||
/**
|
||||
virtual void init(){
|
||||
}
|
||||
|
||||
/**
|
||||
* Set the thermo manager. Note that the flow equations assume
|
||||
* the ideal gas equation.
|
||||
*/
|
||||
|
|
@ -169,6 +172,9 @@ namespace Cantera {
|
|||
|
||||
|
||||
virtual string componentName(int n) const;
|
||||
|
||||
//added by Karl Meredith
|
||||
int componentIndex(string name) const;
|
||||
|
||||
|
||||
virtual void showSolution(const doublereal* x);
|
||||
|
|
@ -236,6 +242,14 @@ namespace Cantera {
|
|||
void setGas(const doublereal* x,int j);
|
||||
void setGasAtMidpoint(const doublereal* x,int j);
|
||||
|
||||
//Karl Meredith
|
||||
// doublereal density_unprotected(int j) const {
|
||||
// return m_rho[j];
|
||||
// }
|
||||
doublereal density(int j) const {
|
||||
return m_rho[j];
|
||||
}
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
|
|
@ -329,10 +343,6 @@ namespace Cantera {
|
|||
return m_wtm[j]*Y(x,k,j)/m_wt[k];
|
||||
}
|
||||
|
||||
doublereal density(int j) const {
|
||||
return m_rho[j];
|
||||
}
|
||||
|
||||
doublereal flux(int k, int j) const {
|
||||
return m_flux(k, j);
|
||||
}
|
||||
|
|
@ -447,6 +457,7 @@ namespace Cantera {
|
|||
|
||||
doublereal m_efctr;
|
||||
|
||||
|
||||
private:
|
||||
|
||||
vector_fp m_ybar;
|
||||
|
|
|
|||
|
|
@ -171,7 +171,22 @@ namespace Cantera {
|
|||
doublereal *xb, *rb;
|
||||
|
||||
// residual equations for the two local variables
|
||||
r[0] = m_mdot - x[0];
|
||||
|
||||
//added by Karl Meredith
|
||||
if (m_adiabatic)
|
||||
// For the adiabatic case, the mass flow rate is not known. For
|
||||
// this case, set mdot (x[0]) to match rho*u in the flow domain
|
||||
// dgg: I think this formulation will only work if the adiabatic
|
||||
// inlet is on the left, i.e., the flow is left-to-right.
|
||||
|
||||
r[0] = m_flow->density(0)*x[2] - x[0];
|
||||
else
|
||||
// Specified mass flow rate
|
||||
r[0] = m_mdot - x[0];
|
||||
|
||||
// The inlet temperature is always specified. For the adiabatic
|
||||
// case, this is the temperature of the gas far upstream of the
|
||||
// flame.
|
||||
r[1] = m_temp - x[1];
|
||||
|
||||
// both are algebraic constraints
|
||||
|
|
@ -189,11 +204,18 @@ namespace Cantera {
|
|||
// spreading rate. Flow domain sets this to V(0),
|
||||
// so for finite spreading rate subtract m_V0.
|
||||
rb[1] -= m_V0;
|
||||
|
||||
rb[3] += x[0]; // lambda
|
||||
|
||||
//added by Karl Meredith
|
||||
if(m_adiabatic){
|
||||
rb[3]=xb[3]; //zeroo lambda (Avoids last species on last node being singular ???)
|
||||
}
|
||||
else{
|
||||
rb[3] += x[0]; // lambda
|
||||
}
|
||||
for (k = 1; k < m_nsp; k++) {
|
||||
if (m_flow->doSpecies(k)) {
|
||||
rb[4+k] += x[0]*m_yin[k];
|
||||
//writelog("Left "+int2str(k)+" "+fp2str(m_yin[k])+"\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -210,6 +232,7 @@ namespace Cantera {
|
|||
if (m_flow->doSpecies(k)) {
|
||||
// rb[4+k] += x[0]*(-xb[4+k] + m_yin[k]);
|
||||
rb[4+k] += x[0]*(m_yin[k]);
|
||||
//writelog("Right "+int2str(k)+" "+fp2str(m_yin[k])+"\n");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -438,7 +461,15 @@ namespace Cantera {
|
|||
db = diag - nc;
|
||||
|
||||
// zero Lambda
|
||||
rb[0] = xb[3]; // zero Lambda
|
||||
|
||||
//added by Karl Meredith
|
||||
if (m_adiabatic) {
|
||||
rb[0] = xb[0] - xb[0-nc]; //zero U gradient
|
||||
//(This makes it so that U at last node is not undefined)
|
||||
}
|
||||
else{
|
||||
rb[0] = xb[3]; // zero Lambda
|
||||
}
|
||||
rb[2] = xb[2] - xb[2 - nc]; // zero T gradient
|
||||
for (k = 5; k < nc; k++) {
|
||||
rb[k] = xb[k] - xb[k - nc]; // zero mass fraction gradient
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue