Resolving Janaf lookup unboundedness problems

This commit is contained in:
Hrvoje Jasak 2010-10-12 23:05:08 +01:00
parent 7c80ada340
commit 801cfd90e9
6 changed files with 61 additions and 24 deletions

View file

@ -911,7 +911,9 @@ Tolerances
patchToPatchProjectionTol 0.05;
// Thermophysical models
specieThermoTol 1e-4;
specieThermoTol 1e-4;
speciesThermoTJump 20;
speciesThermoMaxIter 100;
// Intersection tolerance
intersectionPlanarTol 0.2;

View file

@ -114,10 +114,12 @@ private:
// Private member functions
//- Check given temperature is within the range of the fitted coeffs
inline void checkT(const scalar T) const;
// Note: bounding T within range. HJ, 12/Oct/2010
inline void checkT(scalar& T) const;
//- Return the coefficients corresponding to the given temperature
inline const coeffArray& coeffs(const scalar T) const;
// Note: bounding T within range. HJ, 12/Oct/2010
inline const coeffArray& coeffs(scalar& T) const;
public:
@ -145,19 +147,23 @@ public:
// Member Functions
//- Heat capacity at constant pressure [J/(kmol K)]
inline scalar cp(const scalar T) const;
// Note: bounding T within range. HJ, 12/Oct/2010
inline scalar cp(scalar T) const;
//- Enthalpy [J/kmol]
inline scalar h(const scalar T) const;
// Note: bounding T within range. HJ, 12/Oct/2010
inline scalar h(scalar T) const;
//- Sensible enthalpy [J/kmol]
inline scalar hs(const scalar T) const;
// Note: bounding T within range. HJ, 12/Oct/2010
inline scalar hs(scalar T) const;
//- Chemical enthalpy [J/kmol]
inline scalar hc() const;
//- Entropy [J/(kmol K)]
inline scalar s(const scalar T) const;
// Note: bounding T within range. HJ, 12/Oct/2010
inline scalar s(scalar T) const;
// Member operators

View file

@ -54,17 +54,21 @@ inline Foam::janafThermo<equationOfState>::janafThermo
template<class equationOfState>
inline void Foam::janafThermo<equationOfState>::checkT(const scalar T) const
inline void Foam::janafThermo<equationOfState>::checkT(scalar& T) const
{
if (T < Tlow_ || T > Thigh_)
{
FatalErrorIn
// Improvements: graceful exit with recovery. HJ, 11/Oct/2010
InfoIn
(
"janafThermo<equationOfState>::checkT(const scalar T) const"
"janafThermo<equationOfState>::checkT(scalar& T) const"
) << "attempt to use janafThermo<equationOfState>"
" out of temperature range "
<< Tlow_ << " -> " << Thigh_ << "; T = " << T
<< abort(FatalError);
<< endl;
// Bracket T to avoid out-of-range error
T = Foam::min(Thigh_, Foam::max(T, Tlow_));
}
}
@ -73,9 +77,11 @@ template<class equationOfState>
inline const typename Foam::janafThermo<equationOfState>::coeffArray&
Foam::janafThermo<equationOfState>::coeffs
(
const scalar T
scalar& T
) const
{
// Note: T will be bounded by checkT in coeffs(T). No longer const
// HJ, 12/Oct/2010
checkT(T);
if (T < Tcommon_)
@ -103,7 +109,7 @@ inline Foam::janafThermo<equationOfState>::janafThermo
Thigh_(jt.Thigh_),
Tcommon_(jt.Tcommon_)
{
for (register label coefLabel=0; coefLabel<nCoeffs_; coefLabel++)
for (register label coefLabel = 0; coefLabel < nCoeffs_; coefLabel++)
{
highCpCoeffs_[coefLabel] = jt.highCpCoeffs_[coefLabel];
lowCpCoeffs_[coefLabel] = jt.lowCpCoeffs_[coefLabel];
@ -116,7 +122,7 @@ inline Foam::janafThermo<equationOfState>::janafThermo
template<class equationOfState>
inline Foam::scalar Foam::janafThermo<equationOfState>::cp
(
const scalar T
scalar T
) const
{
const coeffArray& a = coeffs(T);
@ -127,7 +133,7 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::cp
template<class equationOfState>
inline Foam::scalar Foam::janafThermo<equationOfState>::h
(
const scalar T
scalar T
) const
{
const coeffArray& a = coeffs(T);
@ -142,7 +148,7 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::h
template<class equationOfState>
inline Foam::scalar Foam::janafThermo<equationOfState>::hs
(
const scalar T
scalar T
) const
{
return h(T) - hc();
@ -167,9 +173,12 @@ inline Foam::scalar Foam::janafThermo<equationOfState>::hc() const
template<class equationOfState>
inline Foam::scalar Foam::janafThermo<equationOfState>::s
(
const scalar T
scalar T
) const
{
// Note: T will be bounded by checkT in coeffs(T). No longer const
// HJ, 12/Oct/2010
const coeffArray& a = coeffs(T);
return
this->RR*

View file

@ -41,6 +41,13 @@ const Foam::scalar Foam::specieThermo<thermo>::tol_
);
template<class thermo>
const Foam::scalar Foam::specieThermo<thermo>::TJump_
(
debug::tolerances("speciesThermoTJump", 20)
);
template<class thermo>
const int Foam::specieThermo<thermo>::maxIter_
(

View file

@ -98,13 +98,16 @@ class specieThermo
//- Convergence tolerance of energy -> temperature inversion functions
static const scalar tol_;
//- Max temperature jump of energy -> temperature inversion functions
static const scalar TJump_;
//- Max number of iterations in energy->temperature inversion functions
static const int maxIter_;
// Private member functions
//- return the temperature corresponding to the value of the
//- Return the temperature corresponding to the value of the
// thermodynamic property f, given the function f = F(T) and dF(T)/dT
inline scalar T
(
@ -131,7 +134,7 @@ public:
// Member Functions
// Fundamaental properties
// Fundamental properties
// (These functions must be provided in derived types)
// Heat capacity at constant pressure [J/(kmol K)]

View file

@ -50,23 +50,33 @@ inline Foam::scalar Foam::specieThermo<thermo>::T
scalar Test = T0;
scalar Tnew = T0;
scalar Ttol = T0*tol_;
int iter = 0;
int iter = 0;
do
{
// Limit the temperature jump in a single corrector to TJump_
// HJ, 12/Oct/2010
Test = Tnew;
Tnew = Test - ((this->*F)(Test) - f)/(this->*dFdT)(Test);
Tnew = Test
- Foam::min(((this->*F)(Test) - f)/(this->*dFdT)(Test), TJump_);
if (iter++ > maxIter_)
{
FatalErrorIn
// Improvements: graceful exit with recovery. HJ, 11/Oct/2010
InfoIn
(
"specieThermo<thermo>::T(scalar f, scalar T0, "
"scalar (specieThermo<thermo>::*F)(const scalar) const, "
"scalar (specieThermo<thermo>::*dFdT)(const scalar) const"
") const"
) << "Maximum number of iterations exceeded"
<< abort(FatalError);
) << "Maximum number of iterations exceeded. Rescue by HJ"
<< endl;
// Use value where dFdT is calculated using T0. HJ, 11/Oct/2010
Tnew = f/(this->*dFdT)(T0);
return Tnew;
}
} while (mag(Tnew - Test) > Ttol);