using System; using System.Collections.Generic; using System.Text; namespace NG_Tools { public class Therm { // member data //double dT;// current temperature, in Kelvins //double dP;// current pressure, in Pascals //double dD;// molar density, in moles/dm3 //double dRho;// mass density, in kg/m3 double dPdD;// partial deriv of P wrt D double dPdT;// partial deriv of P wrt T double dSi;// ideal gas specific entropy, kJ/kg.K double dTold;// temperature previously used double dMrxold;// mixture molar mass previously used double[] GK_root = new double[5] { 0.14887433898163121088, 0.43339539412924719080, 0.67940956829902440263, 0.86506336668898451073, 0.97390652851717172008 }; double[] GK_weight = new double[5] { 0.29552422471475286217, 0.26926671930999634918, 0.21908636251598204295, 0.14945134915058059038, 0.066671344308688137179 }; //set the number of points for quadrature int GK_points = 5; //equation constants for ideal gas heat capacity, enthalpy and entropy double[,] ThermConstants = new double[21, 11] {{-29776.4, 7.95454, 43.9417, 1037.09, 1.56373, 813.205, -24.9027, 1019.98,-10.1601, 1070.14,-20.0615}, {-3495.34, 6.95587, 0.272892, 662.738,-0.291318,-680.562, 1.78980, 1740.06, 0.0, 100.0, 4.49823}, { 20.7307, 6.96237, 2.68645, 500.371,-2.56429,-530.443, 3.91921, 500.198, 2.13290, 2197.22, 5.81381}, {-37524.4, 7.98139, 24.3668, 752.320, 3.53990, 272.846, 8.44724, 1020.13,-13.2732, 869.510,-22.4010}, {-56072.1, 8.14319,37.0629,735.402,9.38159,247.190,13.4556, 1454.78,-11.7342, 984.518,-24.0426}, {-13773.1, 7.97183, 6.27078,2572.63,2.05010,1156.72,0.0,100.0,0.0,100.0,-3.24989}, {-10085.4, 7.94680,-0.08380,433.801,2.85539, 843.792,6.31595, 1481.43,-2.88457, 1102.23,-0.51551}, {-5565.60, 6.66789, 2.33458,2584.98,.749019, 559.656,0.0,100.0,0.0,100.0,-7.94821}, {-2753.49, 6.95854, 2.02441,1541.22,.096774, 3674.81,0.0,100.0,0.0,100.0,6.23387}, {-3497.45, 6.96302, 2.40013, 2522.05, 2.21752, 1154.15, 0.0,100.0,0.0,100.0,9.19749}, {-72387.0, 17.8143, 58.2062, 1787.39, 40.7621, 808.645, 0.0,100.0,0.0,100.0,-44.1341}, {-72674.8, 18.6383, 57.4178, 1792.73, 38.6599, 814.151, 0.0,100.0,0.0,100.0,-46.1938}, {-91505.5, 21.3861, 74.3410, 1701.58, 47.0587, 775.899, 0.0,100.0,0.0,100.0,-60.2474}, {-83845.2, 22.5012, 69.5789, 1719.58, 46.2164, 802.174, 0.0,100.0,0.0,100.0,-62.2197}, {-94982.5, 26.6225, 80.3819, 1718.49, 55.6598, 802.069, 0.0,100.0,0.0,100.0,-77.5366}, {-103353.0, 30.4029, 90.6941, 1669.32, 63.2028, 786.001, 0.0,100.0,0.0,100.0,-92.0164}, {-109674.0, 34.0847, 100.253, 1611.55, 69.7675, 768.847, 0.0,100.0,0.0,100.0,-106.149}, {-122599.0, 38.5014, 111.446, 1646.48, 80.5015, 781.588, 0.0,100.0,0.0,100.0,-122.444}, {-133564.0, 42.7143, 122.173, 1654.85, 90.2255, 785.564, 0.0,100.0,0.0,100.0,-138.006}, {0.0,4.9680,0.0,100.0,0.0,100.0,0.0,100.0,0.0,100.0,0.0}, {0.0,4.9680,0.0,100.0,0.0,100.0,0.0,100.0,0.0,100.0,0.0}}; // enumerations for indexing of coefficients //public enum CoefficientList { coefA = 0, coefB, coefC, coefD, coefE, coefF, coefG, coefH, coefI, coefJ, coefK } ; public int coefA = 0; public int coefB = 1; public int coefC = 2; public int coefD = 3; public int coefE = 4; public int coefF = 5; public int coefG = 6; public int coefH = 7; public int coefI = 8; public int coefJ = 9; public int coefK = 10; // conversion constant for thermochemical calories to Joules: 1 cal(IT) = 4.1840 J const double CalTH = 4.1840; public Therm() { // initialize 3 history-sensitive variables dSi = 0.0; dTold = 0.0; dMrxold = 0.0; }//Therm() /************************************************************************** *Function:Math.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 Math.Cosh(x) / Math.Sinh(x); }// Math.Coth() /************************************************************************** *Function:Run() *Arguments:ref AGA10.GasPropsSTRUCT , Detail * *Returns:void *Purpose:overall execution control; top level math for SOS and k *Revisions: **************************************************************************/ public void Run(ref NG_Cal.GasPropsSTRUCT ptAGA10, ref Detail ptD) { //local variables double c, x, y, z; //first run basic set of functions within AGA 8 (1994) Detail Method ptD.Run(ref ptAGA10); //find first partial derivative of Z wrt D ptD.dZdD(ptAGA10.dDf); //find real gas cv, cp, specific enthalpy and entropy CprCvrHS(ref ptAGA10, ref 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 * NG_Cal.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 = Math.Sqrt(c); //calculate the real gas isentropic exponent //using expression functionally equivalent to Equation 3.2 ptAGA10.dKappa = (c * ptAGA10.dRhof) / ptAGA10.dPf; return; }//Run() /************************************************************************** *Function:CpiMolar() *Arguments:ref AGA10.GasPropsSTRUCT *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 CpiMolar(ref NG_Cal.GasPropsSTRUCT 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 < NG_Cal.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 / Math.Sinh(DT); Fx = FT / Math.Cosh(FT); Hx = HT / Math.Sinh(HT); Jx = JT / Math.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; }//CpiMolar() /************************************************************************** *Function:Ho() *Arguments:ref AGA10.GasPropsSTRUCT *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 Ho(ref NG_Cal.GasPropsSTRUCT 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 < NG_Cal.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 = Math.Tanh(FT); cothHT = coth(HT); tanhJT = Math.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.0e3; } // Ho() /************************************************************************** *Function:So() *Arguments:ref AGA10.GasPropsSTRUCT *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 So(ref NG_Cal.GasPropsSTRUCT 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 < NG_Cal.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 = Math.Tanh(FT); cothHT = coth(HT); tanhJT = Math.Tanh(JT); sinhDT = Math.Sinh(DT); coshFT = Math.Cosh(FT); sinhHT = Math.Sinh(HT); coshJT = Math.Cosh(JT); Sx += ThermConstants[i, coefK]; Sx += ThermConstants[i, coefB] * Math.Log(T); Sx += ThermConstants[i, coefC] * (DT * cothDT - Math.Log(sinhDT)); Sx -= ThermConstants[i, coefE] * (FT * tanhFT - Math.Log(coshFT)); Sx += ThermConstants[i, coefG] * (HT * cothHT - Math.Log(sinhHT)); Sx -= ThermConstants[i, coefI] * (JT * tanhJT - Math.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.0e3; }//So() /************************************************************************** *Function:CprCvrHS() *Arguments:ref AGA10.GasPropsSTRUCT , Detail * *Returns:void *Purpose:reasonably efficient group calculation of Cp, Cv, H and S *Revisions: **************************************************************************/ void CprCvrHS(ref NG_Cal.GasPropsSTRUCT ptAGA10, ref 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(ref ptAGA10); //find ideal gas enthalpy ptAGA10.dHo = Ho(ref ptAGA10); //find ideal gas entropy Si = So(ref 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 - NG_Cal.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 = NG_Cal.RGAS * ptAGA10.dDf * a; //calculate dPdD, the partial derivative of P wrt D, at T dPdD = NG_Cal.RGAS * ptAGA10.dTf * b; //equation completing molar Cp, cancelling appropriate terms Cpr = Cvr + NG_Cal.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 * NG_Cal.RGAS * ptAGA10.dTf * (ptAGA10.dZf - 1.0 - ptAGA10.dTf * Hinc * 0.5 * ptAGA10.dDf) / ptAGA10.dMrx; // calculate entropy of mixing for (i = 0; i < NG_Cal.NUMBEROFCOMPONENTS; i++) { if (ptAGA10.adMixture[i] != 0) Smixing += ptAGA10.adMixture[i] * Math.Log(ptAGA10.adMixture[i]); } Smixing *= NG_Cal.RGAS; // calculate specific entropy ptAGA10.dS = Si - Smixing - 1000.0 * NG_Cal.RGAS * (Math.Log(ptAGA10.dPf / 101325.0) - Math.Log(ptAGA10.dZf) + Sinc * 0.5 * ptAGA10.dDf) / ptAGA10.dMrx; return; }//CprCvrHS() /************************************************************************** *Function:HS_Mode() *Arguments:ref AGA10.GasPropsSTRUCT , 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 GasPropsSTRUCT and set *bGuess = true. The initial guess allows the search function to be more *aggressive and, typically, faster. * *Revisions: **************************************************************************/ public void HS_Mode(ref NG_Cal.GasPropsSTRUCT ptAGA10, ref 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 = NG_Cal.P_MAX; pmin = 10000.0; // 10 kPa tmax = NG_Cal.T_MAX; tmin = NG_Cal.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(ref ptAGA10); //h1 is difference between h given and h@Tf, Pf h1 = this.H(ref ptAGA10, ref ptD) - h0; //outer loop: search for a t2 which will satisfy constant enthalpy for (i = 0; i < NG_Cal.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(ref ptAGA10); s1 = this.S(ref ptAGA10, ref ptD) - s0; //inside loop: search for a p2 which will satisfy constant entropy for (j = 0; j < NG_Cal.MAX_NUM_OF_ITERATIONS; j++) { ptAGA10.dPf = p2; ptD.Run(ref ptAGA10); s2 = this.S(ref ptAGA10, ref ptD) - s0; //calculate our proportional change delta2 = Math.Abs(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 >= NG_Cal.MAX_NUM_OF_ITERATIONS) ptAGA10.lStatus = NG_Cal.MAX_NUM_OF_ITERATIONS_EXCEEDED; //calc enthalpy at guessed P & current iter T h2 = this.H(ref ptAGA10, ref ptD) - h0; //calculate our proportional change delta1 = Math.Abs(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(ref ptAGA10); h2 = this.H(ref ptAGA10, ref ptD) - h0; } t1 = t0; h1 = h2; } // check for failure to converge if (i >= NG_Cal.MAX_NUM_OF_ITERATIONS) ptAGA10.lStatus = NG_Cal.MAX_NUM_OF_ITERATIONS_EXCEEDED; }//HS_Mode() /************************************************************************** *Function:H() *Arguments:ref AGA10.GasPropsSTRUCT , Detail * *Returns:double *Purpose:real gas specific enthalpy *Revisions: **************************************************************************/ double H(ref NG_Cal.GasPropsSTRUCT ptAGA10, ref Detail ptD) { double Hinc; double x; int i; //initialize integral Hinc = 0.0; //find ideal gas enthalpy ptAGA10.dHo = Ho(ref 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 * NG_Cal.RGAS * ptAGA10.dTf * (ptAGA10.dZf - 1.0 - ptAGA10.dTf * Hinc * 0.5 * ptAGA10.dDf) / ptAGA10.dMrx; return (ptAGA10.dH); } // H() /************************************************************************** *Function:S() *Arguments:ref AGA10.GasPropsSTRUCT , Detail * *Returns:double *Purpose:real gas specific entropy *Revisions: **************************************************************************/ double S(ref NG_Cal.GasPropsSTRUCT ptAGA10, ref 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(ref ptAGA10); dTold = ptAGA10.dTf; dMrxold = ptAGA10.dMrx; } //calculate entropy of mixing for (i = 0; i < NG_Cal.NUMBEROFCOMPONENTS; i++) { if (ptAGA10.adMixture[i] != 0) Smixing += ptAGA10.adMixture[i] * Math.Log(ptAGA10.adMixture[i]); } Smixing *= NG_Cal.RGAS; // calculate specific entropy ptAGA10.dS = dSi - Smixing - 1000.0 * NG_Cal.RGAS * (Math.Log(ptAGA10.dPf / 101325.0) - Math.Log(ptAGA10.dZf) + Sinc * 0.5 * ptAGA10.dDf) / ptAGA10.dMrx; return (ptAGA10.dS); } // S() } }