775 lines
18 KiB
C++
775 lines
18 KiB
C++
|
||
/*************************************************************************
|
||
* File: therm.cpp
|
||
* Description: Contains thermodynamic functions for the meter object.
|
||
* heat capacity, enthalpy, entropy, sound speed
|
||
* Contains the functions:
|
||
* Therm(), ~Therm(), Run(), coth(), CpiMolar(), Ho(), So(),
|
||
* CprCvrHS(), HS_Mode(), H(), S()
|
||
* Version: ver 1.7 2002.11.17
|
||
* Author: W.B. Peterson
|
||
*Revisions:
|
||
*Copyright (c) 2002 American Gas Association
|
||
|
||
**************************************************************************/
|
||
|
||
#include "aga10.h" #include "detail.h" #include "therm.h" #include <math.h>
|
||
|
||
/**************************************************************************
|
||
|
||
* Function : Therm::Therm()
|
||
* Arguments : void
|
||
* Returns :
|
||
* Purpose : default constructor
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
Therm::Therm(void)
|
||
|
||
{
|
||
// initialize 3 history-sensitive variables dSi = 0.0 ;
|
||
|
||
dTold = 0.0 ; dMrxold = 0.0 ;
|
||
} // Therm::Therm()
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::~Therm()
|
||
* Arguments :
|
||
* Returns : default destructor
|
||
* Purpose : void
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
Therm::~Therm()
|
||
|
||
{
|
||
|
||
}//Therm::~Therm()
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : coth()
|
||
* Arguments : double
|
||
* Returns : double
|
||
* Purpose : calculate hyperbolic cotangent; used in Ho calculations
|
||
* Revisions :
|
||
* Notes : Not a Therm object class member, just a utility for this
|
||
* file. The C++ language has no intrinsic support for
|
||
* hyperbolic cotangent
|
||
**************************************************************************/
|
||
|
||
double coth (double x)
|
||
|
||
{
|
||
return cosh(x)/sinh(x);
|
||
|
||
}// coth()
|
||
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::Run()
|
||
* Arguments : AGA10STRUCT *, Detail *
|
||
* Returns : void
|
||
* Purpose : overall execution control; top level math for SOS and k
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
void Therm::Run(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
|
||
{
|
||
//local variables double c, x, y, z ;
|
||
|
||
//first run basic set of functions within AGA 8 (1994) Detail Method ptD->Run(ptAGA10) ;
|
||
|
||
//find first partial derivative of Z wrt D
|
||
|
||
ptD->dZdD(ptAGA10->dDf) ;
|
||
|
||
//find real gas cv, cp, specific enthalpy and entropy CprCvrHS(ptAGA10, ptD) ;
|
||
|
||
//ratio of real gas specific heats
|
||
|
||
ptAGA10->dk = ptAGA10->dCp / ptAGA10->dCv ;
|
||
|
||
//solve c in three steps, for clarity and ease of debugging x = ptAGA10->dk * RGAS * 1000.0 * ptAGA10->dTf ;
|
||
|
||
y = ptAGA10->dMrx ;
|
||
z = ptAGA10->dZf + ptAGA10->dDf * ptD->ddZdD ;
|
||
|
||
//calculate c, which is SOS^2
|
||
|
||
c = (x / y) * z ;
|
||
|
||
//speed of sound ptAGA10->dSOS = sqrt(c) ;
|
||
|
||
//calculate the real gas isentropic exponent
|
||
|
||
//using expression functionally equivalent to Equation 3.2 ptAGA10->dKappa = (c * ptAGA10->dRhof) / ptAGA10->dPf ;
|
||
|
||
|
||
|
||
|
||
return ;
|
||
|
||
}//Therm::Run()
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::CpiMolar()
|
||
* Arguments : AGA10STRUCT *
|
||
* Returns : double
|
||
* Purpose : Calculate constant pressure ideal gas molar heat capacity
|
||
* in (J/mol-K), applying eqns from Aly, Lee, McFall
|
||
* Notes : For continuity, the original constants and eqn's have been
|
||
* retained. Conversion from thermochemical calories(th) to
|
||
* Joules is applied after the primary calculations are complete.
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
double Therm::CpiMolar(AGA10STRUCT *ptAGA10)
|
||
|
||
{
|
||
double Cp = 0.0 ; double Cpx ;
|
||
|
||
double DT, FT, HT, JT ; double Dx, Fx, Hx, Jx ; double T ;
|
||
|
||
int i ;
|
||
|
||
//to maximize readability of this section, use intermediate variable T T = ptAGA10->dTf ;
|
||
|
||
//calculate heat capacity for each component
|
||
|
||
for (i= 0; i< NUMBEROFCOMPONENTS; i++)
|
||
{
|
||
|
||
//skip species whose concentration is zero if (ptAGA10->adMixture[i] <= 0.0) continue ;
|
||
|
||
//initialize Cp of species to zero
|
||
|
||
Cpx = 0.0 ;
|
||
|
||
// calculate species intermediate terms DT = ThermConstants[i][coefD] / T ;
|
||
|
||
|
||
|
||
|
||
FT = ThermConstants[i][coefF] / T ;
|
||
|
||
HT = ThermConstants[i][coefH] / T ;
|
||
JT = ThermConstants[i][coefJ] / T ;
|
||
|
||
// use intermediate terms to avoid redundant calcs Dx = DT/sinh(DT) ;
|
||
Fx = FT/cosh(FT) ;
|
||
Hx = HT/sinh(HT) ;
|
||
Jx = JT/cosh(JT) ;
|
||
|
||
Cpx += ThermConstants[i][coefB] ;
|
||
|
||
Cpx += ThermConstants[i][coefC] * Dx * Dx ;
|
||
Cpx += ThermConstants[i][coefE] * Fx * Fx ;
|
||
Cpx += ThermConstants[i][coefG] * Hx * Hx ;
|
||
Cpx += ThermConstants[i][coefI] * Jx * Jx ;
|
||
|
||
//use current mole fraction to weight the contribution Cpx *= ptAGA10->adMixture[i];
|
||
|
||
//add this contribution to the sum
|
||
|
||
Cp += Cpx ;
|
||
|
||
}
|
||
|
||
// convert from cal(th)/mol-K to J/mol-K Cp *= CalTH ;
|
||
|
||
return Cp ;
|
||
|
||
} // Therm::CpiMolar()
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::Ho()
|
||
* Arguments : AGA10STRUCT *
|
||
* Returns : double
|
||
* Purpose : Calculate ideal gas specific enthalpy (J/kg)
|
||
* Notes : For continuity, the original constants and eqn's have been
|
||
* retained. Conversion from thermochemical calories(th) to
|
||
* Joules is applied after the primary calculations are complete.
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
double Therm::Ho(AGA10STRUCT *ptAGA10)
|
||
|
||
{
|
||
double H = 0.0 ; double Hx ;
|
||
|
||
double DT, FT, HT, JT ;
|
||
|
||
double cothDT, tanhFT, cothHT, tanhJT ; double T ;
|
||
|
||
int i ;
|
||
|
||
// to maximize readability of this section, use intermediate variable T T = ptAGA10->dTf ;
|
||
|
||
for (i= 0; i< NUMBEROFCOMPONENTS; i++)
|
||
|
||
{
|
||
// skip species whose concentration is zero if (ptAGA10->adMixture[i] <= 0.0) continue ;
|
||
|
||
Hx = 0.0 ;
|
||
|
||
// calculate species intermediate terms DT = ThermConstants[i][coefD] / T ;
|
||
|
||
FT = ThermConstants[i][coefF] / T ; HT = ThermConstants[i][coefH] / T ; JT = ThermConstants[i][coefJ] / T ;
|
||
|
||
cothDT = coth(DT) ; tanhFT = tanh(FT) ; cothHT = coth(HT) ; tanhJT = tanh(JT) ;
|
||
|
||
|
||
|
||
Hx += ThermConstants[i][coefA] ;
|
||
|
||
Hx += ThermConstants[i][coefB] * T ;
|
||
Hx += ThermConstants[i][coefC] * ThermConstants[i][coefD] * cothDT;
|
||
Hx -= ThermConstants[i][coefE] * ThermConstants[i][coefF] * tanhFT;
|
||
Hx += ThermConstants[i][coefG] * ThermConstants[i][coefH] * cothHT;
|
||
|
||
Hx -= ThermConstants[i][coefI] * ThermConstants[i][coefJ] * tanhJT;
|
||
|
||
//use current mole fraction to weight the contribution Hx *= ptAGA10->adMixture[i];
|
||
|
||
//add this contribution to the sum
|
||
|
||
H += Hx ;
|
||
}
|
||
|
||
//convert from cal(th)/g-mol to kJ/kg-mol H *= CalTH ;
|
||
|
||
//convert from kJ/kg-mol to J/kg
|
||
|
||
H /= ptAGA10->dMrx ;
|
||
|
||
// return in J/kg return H * 1.e3; } // Therm::Ho()
|
||
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::So()
|
||
* Arguments : AGA10STRUCT *
|
||
* Returns : double
|
||
* Purpose : ideal gas specific entropy (J/kg-K)
|
||
* Notes : For continuity, the original constants and eqn's have been
|
||
* retained. Conversion from thermochemical calories(th) to
|
||
* Joules is applied after the primary calculations are complete.
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
double Therm::So(AGA10STRUCT *ptAGA10)
|
||
|
||
{
|
||
double S = 0.0 ; double Sx ;
|
||
|
||
double DT, FT, HT, JT ;
|
||
|
||
double cothDT, tanhFT, cothHT, tanhJT ; double sinhDT, coshFT, sinhHT, coshJT ; double T ;
|
||
|
||
int i ;
|
||
|
||
// to improve readability of this section, use intermediate variable T T = ptAGA10->dTf ;
|
||
|
||
for (i= 0; i< NUMBEROFCOMPONENTS; i++)
|
||
|
||
{
|
||
|
||
// skip species whose concentration is zero if (ptAGA10->adMixture[i] <= 0.0) continue ;
|
||
|
||
Sx = 0.0 ;
|
||
|
||
// calculate species intermediate terms DT = ThermConstants[i][coefD] / T ;
|
||
|
||
FT = ThermConstants[i][coefF] / T ; HT = ThermConstants[i][coefH] / T ; JT = ThermConstants[i][coefJ] / T ;
|
||
|
||
cothDT = coth(DT) ; tanhFT = tanh(FT) ; cothHT = coth(HT) ;
|
||
|
||
|
||
|
||
|
||
tanhJT = tanh(JT) ;
|
||
|
||
sinhDT = sinh(DT) ; coshFT = cosh(FT) ; sinhHT = sinh(HT) ; coshJT = cosh(JT) ;
|
||
|
||
Sx += ThermConstants[i][coefK] ;
|
||
|
||
Sx += ThermConstants[i][coefB] * log(T) ;
|
||
Sx += ThermConstants[i][coefC] * (DT * cothDT - log(sinhDT)) ;
|
||
|
||
Sx -= ThermConstants[i][coefE] * (FT * tanhFT - log(coshFT)) ;
|
||
|
||
Sx += ThermConstants[i][coefG] * (HT * cothHT - log(sinhHT)) ;
|
||
Sx -= ThermConstants[i][coefI] * (JT * tanhJT - log(coshJT)) ;
|
||
|
||
//use current mole fraction to weight the contribution Sx *= ptAGA10->adMixture[i];
|
||
|
||
//add this contribution to the sum
|
||
|
||
S += Sx ;
|
||
}
|
||
|
||
//convert cal(th)/mol-K basis to to kJ/kg mol-K S *= CalTH ;
|
||
|
||
//convert from kJ/kg mol-K to kJ/kg-K
|
||
|
||
S /= ptAGA10->dMrx ;
|
||
|
||
// return in J/kg-K return S * 1.e3;
|
||
|
||
} // Therm::So()
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::CprCvrHS()
|
||
* Arguments : AGA10STRUCT *, Detail *
|
||
* Returns : void
|
||
* Purpose : reasonably efficient group calculation of Cp, Cv, H and S
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
void Therm::CprCvrHS(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
|
||
{
|
||
double Cvinc, Cvr, Cpr ; double Hinc ;
|
||
|
||
double Sinc ; double Smixing ; double Cp, Si ; double a, b, x ; int i ;
|
||
|
||
//initialize integrals to zero Cvinc = 0.0 ;
|
||
|
||
Hinc = 0.0 ; Sinc = 0.0 ;
|
||
|
||
//initialize entropy of mixing Smixing = 0.0 ;
|
||
|
||
//find ideal gas Cp
|
||
|
||
Cp = CpiMolar(ptAGA10) ;
|
||
|
||
//find ideal gas enthalpy ptAGA10->dHo = Ho(ptAGA10) ;
|
||
|
||
//find ideal gas entropy Si = So(ptAGA10) ;
|
||
|
||
//calculate ideal gas specific heat capacity at constant pressure in J/kgK ptAGA10->dCpi = (Cp * 1000.0) / ptAGA10->dMrx ;
|
||
|
||
//integrate partial derivatives from D=0 to D=D, applying Gauss-Kronrod quadrature for ( i= 0; i < GK_points; i++)
|
||
|
||
|
||
|
||
|
||
{
|
||
|
||
// set calculation point at + abscissa
|
||
x = ptAGA10->dDf * (1.0 + GK_root[i]) / 2.0 ;
|
||
|
||
//get Z at D ptD->zdetail(x) ; ptD->dZdT(x) ; ptD->d2ZdT2(x) ;
|
||
|
||
//gather contributions at + abscissa; applying weighting factor Hinc += GK_weight[i] * ptD->ddZdT / x ;
|
||
|
||
Cvinc += GK_weight[i] * (2.0 * ptD->ddZdT + ptAGA10->dTf * ptD->dd2ZdT2) / x ; Sinc += GK_weight[i] * (ptD->dZ + ptAGA10 ->dTf * ptD->ddZdT - 1.0) / x ;
|
||
|
||
//set calculation point at - abscissa
|
||
|
||
x = ptAGA10->dDf * (1.0 - GK_root[i]) / 2.0 ;
|
||
//get Z at D ptD->zdetail(x) ;
|
||
//calculate 1st and 2nd partial derivatives of Z wrt T ptD->dZdT(x) ;
|
||
|
||
ptD->d2ZdT2(x) ;
|
||
|
||
//gather contributions at - abscissa; applying weighting factor Hinc += GK_weight[i] * ptD->ddZdT / x ;
|
||
|
||
Cvinc += GK_weight[i] * (2.0 * ptD->ddZdT + ptAGA10->dTf * ptD->dd2ZdT2) / x ; Sinc += GK_weight[i] * (ptD->dZ + ptAGA10 ->dTf * ptD->ddZdT - 1.0) / x ;
|
||
|
||
}
|
||
|
||
//return Z and partial derivatives to full molar density ptD->zdetail(ptAGA10->dDf) ;
|
||
|
||
ptD->dZdT(ptAGA10->dDf) ; ptD->d2ZdT2(ptAGA10->dDf) ;
|
||
|
||
//complete Cv molar
|
||
|
||
Cvr = Cp - RGAS * (1.0 + ptAGA10->dTf * Cvinc * 0.5 * ptAGA10->dDf) ;
|
||
|
||
//intermediate values for Cp, containing 2 partial derivatives a =(ptAGA10->dZf + ptAGA10->dTf * ptD->ddZdT) ;
|
||
|
||
b =(ptAGA10->dZf + ptAGA10->dDf * ptD->ddZdD) ;
|
||
|
||
//calculate dPdT, the partial derivative of P wrt T, at D
|
||
|
||
|
||
|
||
|
||
dPdT = RGAS * ptAGA10->dDf * a ;
|
||
|
||
//calculate dPdD, the partial derivative of P wrt D, at T dPdD = RGAS * ptAGA10->dTf * b ;
|
||
|
||
//equation completing molar Cp, cancelling appropriate terms Cpr = Cvr + RGAS * ((a * a)/b) ;
|
||
|
||
//change from molar to mass basis
|
||
|
||
Cpr /= ptAGA10->dMrx ;
|
||
Cvr /= ptAGA10->dMrx ;
|
||
|
||
// write to the data stucture
|
||
|
||
ptAGA10->dCv = Cvr * 1000.0 ; // convert from joules/kgK to kilojoules/kgK ptAGA10->dCp = Cpr * 1000.0 ;
|
||
|
||
// calculate specific enthalpy
|
||
|
||
ptAGA10->dH = ptAGA10->dHo + 1000.0 * RGAS * ptAGA10->dTf *
|
||
(ptAGA10->dZf - 1.0 - ptAGA10->dTf * Hinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx ;
|
||
|
||
// calculate entropy of mixing
|
||
|
||
for (i= 0; i< NUMBEROFCOMPONENTS; i++)
|
||
{
|
||
if (ptAGA10->adMixture[i]) Smixing += ptAGA10->adMixture[i] * log(ptAGA10->adMixture[i]) ;
|
||
|
||
}
|
||
|
||
Smixing *= RGAS ;
|
||
|
||
// calculate specific entropy
|
||
|
||
ptAGA10->dS = Si – Smixing - 1000.0 * RGAS * (log(ptAGA10->dPf/101325.0) - log(ptAGA10->dZf) + Sinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx ;
|
||
|
||
return ;
|
||
|
||
} // Therm::CprCvrHS()
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::HS_Mode()
|
||
* Arguments : AGA10STRUCT *, Detail *, double, double, bool
|
||
* Returns : void
|
||
* Purpose : Calculates a pressure & temperature from known enthalpy & entropy,
|
||
* with or without prior estimates.This function has a role in the
|
||
* calculation of C*.
|
||
* Solution based on a doubly-nested trial & error algorithm and Newton's
|
||
* method.
|
||
*
|
||
* For illustrative purpose, two approaches are supported by this example.
|
||
* If you are starting without advance knowledge of P & T, set the input parm
|
||
* bGuess to false, thus specifying a conservative search approach.
|
||
* If, however, you have a basis for guessing P & T (plenum conditions of a
|
||
* critical flow nozzle, for example) set P & T via AGA10STRUCT and set
|
||
* bGuess = true. The initial guess allows the search function to be more
|
||
* aggressive and, typically, faster.
|
||
*
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
void Therm::HS_Mode(AGA10STRUCT *ptAGA10, Detail *ptD, double H, double S, bool bGuess)
|
||
|
||
{
|
||
double s0, s1, s2, t0, t1, t2, tmin, tmax ; double h0, h1, h2, p0, p1, p2, px, pmin, pmax ; double delta1, delta2 ;
|
||
|
||
double tolerance = 0.001 ;// convergence tolerance (used for both H and S searches) int i, j ;
|
||
|
||
//s0and h0 are our real gas reference points s0 = S ;
|
||
|
||
h0 = H ;
|
||
|
||
//calling function specifies whether search parameters are supplied thru ptAGA10 or unknown if (bGuess)
|
||
{
|
||
t1 = ptAGA10->dTf ; px = ptAGA10->dPf ; pmax = px * 2.0 ; pmin = px * 0.1 ;
|
||
|
||
|
||
|
||
|
||
tmax = t1 * 1.5 ; tmin = t1 * 0.67 ;
|
||
}
|
||
else // use arbitrary, generic limits
|
||
{
|
||
|
||
t1 = 273.15 ;
|
||
|
||
px = 1013250.0 ; // 10 atmospheres pmax = P_MAX ;
|
||
|
||
pmin = 10000.0 ; // 10 kPa tmax = T_MAX ;
|
||
|
||
tmin = T_MIN ;
|
||
|
||
}
|
||
|
||
// set the temperature differential t2 = t1 + 10.0 ;
|
||
|
||
///////////////////////////////////////////
|
||
|
||
//begin double trial-and-error, searching for T & P
|
||
//run the calculation with initial guesses
|
||
ptD->Run(ptAGA10) ;
|
||
|
||
//h1 is difference between h given and h@Tf, Pf
|
||
|
||
h1 = this->H(ptAGA10, ptD) - h0;
|
||
|
||
//outer loop: search for a t2 which will satisfy constant enthalpy for ( i= 0; i < MAX_NUM_OF_ITERATIONS; i++)
|
||
{
|
||
ptAGA10->dTf = t2 ;
|
||
p1 = px ;// reset one bracket
|
||
p2 = px * 0.1 ;// set other bracket to 0.1x the upper bracket ptAGA10->dPf = p1 ;
|
||
ptD->Run(ptAGA10) ;
|
||
s1 = this->S(ptAGA10, ptD) - s0;
|
||
//inside loop: search for a p2 which will satisfy constant entropy for (j= 0; j < MAX_NUM_OF_ITERATIONS; j++)
|
||
{
|
||
ptAGA10->dPf = p2 ; ptD->Run(ptAGA10) ;
|
||
|
||
s2 = this->S(ptAGA10, ptD) - s0 ;
|
||
//calculate our proportional change
|
||
|
||
|
||
|
||
|
||
delta2 = fabs(s1 - s2) / s0 ; // close enough?
|
||
|
||
if (delta2 < tolerance) break ;
|
||
//revise our estimate to p2 p0 = p2 ;
|
||
|
||
p2 = (p1 * s2 - p2 * s1) / (s2 - s1) ;
|
||
|
||
//check for negative pressure and clamp to pmin for safety if (p2 <= pmin)
|
||
|
||
{
|
||
p2 = pmin ;
|
||
}
|
||
|
||
//check if we've created an unrealistic pressure
|
||
|
||
if (p2 >= pmax ) p2 = pmax ; // swap values
|
||
|
||
p1 = p0 ;
|
||
s1 = s2 ;
|
||
|
||
}
|
||
// check for failure to converge
|
||
if (j >= MAX_NUM_OF_ITERATIONS) ptAGA10->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED ;
|
||
|
||
//calc enthalpy at guessed P & current iter T h2 = this->H(ptAGA10, ptD) - h0 ;
|
||
|
||
//calculate our proportional change
|
||
|
||
delta1 = fabs(h1 - h2) / h0 ;
|
||
|
||
// close enough?
|
||
|
||
if (delta1 < tolerance && i > 0) break ;
|
||
|
||
//revise our estimate to t2 t0 = t2 ;
|
||
|
||
t2 = (t1 * h2 - t2 * h1) / (h2 - h1) ;
|
||
|
||
//check if we've created an unrealistic temperature if (t2 >= tmax ) t2 = tmax ;
|
||
|
||
//revise t2, if necessary
|
||
|
||
if (t2 <= tmin )
|
||
{
|
||
|
||
t2 = t0 + 10.0 ;
|
||
|
||
|
||
|
||
|
||
ptAGA10->dTf = t2 ; ptD->Run(ptAGA10) ;
|
||
h2 = this->H(ptAGA10, ptD) - h0 ;
|
||
}
|
||
|
||
t1 = t0 ;
|
||
|
||
h1 = h2 ;
|
||
}
|
||
// check for failure to converge
|
||
if (i >= MAX_NUM_OF_ITERATIONS) ptAGA10->lStatus = MAX_NUM_OF_ITERATIONS_EXCEEDED ;
|
||
|
||
}//Therm::HS_Mode()
|
||
|
||
|
||
/**************************************************************************
|
||
|
||
* Function : Therm::H()
|
||
* Arguments : AGA10STRUCT *, Detail *
|
||
* Returns : double
|
||
* Purpose : real gas specific enthalpy
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
double Therm::H(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
|
||
{
|
||
double Hinc ; double x ; int i ;
|
||
|
||
//initialize integral Hinc = 0.0 ;
|
||
|
||
//find ideal gas enthalpy ptAGA10->dHo = Ho(ptAGA10) ;
|
||
|
||
//integrate partial derivatives from D=0 to D=D, applying Gauss-Kronrod quadrature for ( i= 0; i < GK_points; i++)
|
||
{
|
||
//calculate 1st and 2nd partial derivatives of Z wrt T
|
||
|
||
x = ptAGA10->dDf * (1.0 + GK_root[i]) / 2.0 ; ptD->zdetail(x) ;
|
||
|
||
ptD->dZdT(x) ;
|
||
|
||
|
||
|
||
|
||
ptD->d2ZdT2(x) ;
|
||
|
||
Hinc += GK_weight[i] * ptD->ddZdT / x ; if (i == 10) break;
|
||
|
||
x = ptAGA10->dDf * (1.0 - GK_root[i]) / 2.0 ; ptD->zdetail(x) ;
|
||
|
||
ptD->dZdT(x) ; ptD->d2ZdT2(x) ;
|
||
|
||
Hinc += GK_weight[i] * ptD->ddZdT / x ;
|
||
|
||
}
|
||
|
||
ptD->zdetail(ptAGA10->dDf) ; ptD->dZdT(ptAGA10->dDf) ; ptD->d2ZdT2(ptAGA10->dDf) ;
|
||
|
||
// calculate specific enthalpy
|
||
|
||
ptAGA10->dH = ptAGA10->dHo + 1000.0 * RGAS * ptAGA10->dTf *
|
||
(ptAGA10->dZf - 1.0 - ptAGA10->dTf * Hinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx ;
|
||
|
||
return(ptAGA10->dH) ; } // Therm::H()
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
/**************************************************************************
|
||
* Function : Therm::S()
|
||
* Arguments : AGA10STRUCT *, Detail *
|
||
* Returns : double
|
||
* Purpose : real gas specific entropy
|
||
* Revisions :
|
||
**************************************************************************/
|
||
|
||
double Therm::S(AGA10STRUCT *ptAGA10, Detail *ptD)
|
||
|
||
{
|
||
double Sinc ; double Smixing ; double x ;
|
||
|
||
int i ;
|
||
|
||
//initialize integral Sinc = 0.0 ;
|
||
|
||
//initialize entropy of mixing Smixing = 0.0 ;
|
||
|
||
//integrate partial derivatives from D=0 to D=D, applying Gauss-Kronrod quadrature for ( i= 0; i < GK_points; i++)
|
||
|
||
{
|
||
//calculate 1st and 2nd partial derivatives of Z wrt T
|
||
|
||
x = ptAGA10->dDf * (1.0 + GK_root[i]) / 2.0 ; ptD->zdetail(x) ;
|
||
|
||
ptD->dZdT(x) ; ptD->d2ZdT2(x) ;
|
||
|
||
Sinc += GK_weight[i] * (ptD->dZ + ptAGA10 ->dTf * ptD->ddZdT - 1.0) / x ; if (i == 10) break;
|
||
|
||
x = ptAGA10->dDf * (1.0 - GK_root[i]) / 2.0 ; ptD->zdetail(x) ;
|
||
|
||
ptD->dZdT(x) ; ptD->d2ZdT2(x) ;
|
||
|
||
Sinc += GK_weight[i] * (ptD->dZ + ptAGA10 ->dTf * ptD->ddZdT - 1.0) / x ;
|
||
|
||
}
|
||
|
||
|
||
|
||
|
||
//reset Z and partial deivatives dZdT and d2ZdT2 ptD->zdetail(ptAGA10->dDf) ; ptD->dZdT(ptAGA10->dDf) ; ptD->d2ZdT2(ptAGA10->dDf) ;
|
||
|
||
//find ideal gas entropy, but only if temperature or composition have changed if (ptAGA10->dTf != dTold || ptAGA10->dMrx != dMrxold)
|
||
|
||
{
|
||
dSi = So(ptAGA10) ; dTold = ptAGA10->dTf ; dMrxold = ptAGA10->dMrx ;
|
||
}
|
||
|
||
//calculate entropy of mixing
|
||
|
||
for (i= 0; i< NUMBEROFCOMPONENTS; i++)
|
||
|
||
{
|
||
if (ptAGA10->adMixture[i]) Smixing += ptAGA10->adMixture[i] * log(ptAGA10->adMixture[i]) ;
|
||
}
|
||
Smixing *= RGAS ;
|
||
|
||
// calculate specific entropy
|
||
|
||
ptAGA10->dS = dSi – Smixing - 1000.0 * RGAS * (log(ptAGA10->dPf/101325.0) - log(ptAGA10->dZf) + Sinc * 0.5 * ptAGA10->dDf) / ptAGA10->dMrx ;
|
||
|
||
return(ptAGA10->dS) ; } // Therm::S()
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|