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