GasFlowMeter/User/AGA10/therm.cpp

/*************************************************************************
* 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()