ldeyun/flowMeter
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
75ec21431b
flowMeter/NGToolsPC/Therm.cs

535 lines
24 KiB
C#
Raw Blame History

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()
}
}
Reference in New Issue View Git Blame Copy Permalink