using System; using System.Collections.Generic; using System.Text; namespace NG_Tools { /************************************************************************* * File:detail.h * Description:Header file for the 'Detail' class *See 'detail.cpp' for the implementation. * Version:ver 1.72002.11.17 * Author:W.B. Peterson *Revisions: *Copyright (c) 2002 American Gas Association **************************************************************************/ public class Detail { // member data int iNCC;// number of components int[] aiCID = new int[21];// component IDs //five history variables are used to improve efficiency during repeated calculations double dOldMixID; // mixture ID from previous calc double dOldPb; // Pb from previous calc double dOldTb; // Tb from previous calc double dOldPf; // Pf from previous calc double dOldTf; // Tf from previous calc //EOS parameters from table 4, column 1 double[] adAn = new double[58]; double[] adUn = new double[58]; // characterization parameters from table 5 double[] dMri = new double[21];// molecular weight of ith component double[] dEi = new double[21]; // characteristic energy parameter for ith component double[] dKi = new double[21]; // size parameter for ith component - m^3/kg-mol ^1/3 double[] dGi = new double[21]; // orientation parameter double[] dQi = new double[21]; // quadrupole parameter double[] dFi = new double[21]; // high temperature parameter double[] dSi = new double[21];// dipole parameter double[] dWi = new double[21];// association parameter double[,] dEij = new double[21, 21]; // virial coefficient energy binary interaction parm double[,] dUij = new double[21, 21]; // binary interaction parameter for conformal energy double[,] dKij = new double[21, 21]; // binary interaction parameter for size double[,] dGij = new double[21, 21];// binary interaction parameter for orientation double[,] adTable6Eij = new double[21, 21]; // Table 6 constants double[,] adTable6Uij = new double[21, 21]; // Table 6 constants double[,] adTable6Kij = new double[21, 21]; // Table 6 constants double[,] adTable6Gij = new double[21, 21]; // Table 6 constants double[] adTable5Qi = new double[21]; // table 5 constants double[] adTable5Fi = new double[21]; // table 5 constants double[] adTable5Si = new double[21]; // table 5 constants double[] adTable5Wi = new double[21]; // table 5 constants double[] dXi = new double[21];// mole fraction of component i double dPCalc;// pressure calculated by pdetail() double dT;// current temperature double dP;// current pressure double dRhoTP;// molar density at T & P double dB;// 2nd virial coefficient, B double[] adBcoef = new double[18];// 18 coefficients to calculate B double[] adFn = new double[58];// function for coefficients of density double[] fx = new double[58];// modified coefficients used for 3 derivs double dU;// mixture energy parameter double dKp3;// mixture size parameter ^3 double dW;// mixture orientation parameter double dQp2;// mixture quadrupole parameter ^2 double dF;// high temperature parameter double dRho;// molar density double dRhoL;// low density used in braket function double dRhoH;// high density used in braket function double dPRhoL; // low pressure used in braket function double dPRhoH; // high pressure used in braket function //public variables also used for advanced fluid property calculations public double dZ;// current compressibility public double ddZdT;// first partial derivative of Z wrt T public double dd2ZdT2;// second partial derivative of Z wrt T public double ddZdD;// first partial derivative of Z wrt molar density public double ddBdT;// first partial derivative of B wrt T public double dd2BdT2;// second partial derivative of B wrt T /************************************************************************** *Function:Detail() *Arguments:void *Returns: *Purpose:default constructor; includes initialization of *history-sensitive variables & data tables 4 and 6 *Revisions: **************************************************************************/ public Detail() { //initialize history-sensitive variables dOldMixID = 0.0; // mixture ID from previous calc dOldPb = 0.0; // base pressure from previous calc dOldTb = 0.0; // base temperature from previous calc dOldPf = 0.0; // flowing pressure from previous calc dOldTf = 0.0; // flowing temperature from previous calc //initialize gas component array used within this class for (int i = 0; i < NG_Cal.NUMBEROFCOMPONENTS; i++) dXi[i] = 0; // function table() populates tables of static constants table(); }// Detail() /************************************************************************** *Function:compositionchange() *Arguments:GasPropsSTRUCT * *Returns:void *Purpose:Compares new composition to old by creating a semi-unique *numerical ID. It is possible but very unlikely that 2 *sequential & different compositions will produce the same ID *Revisions: **************************************************************************/ public bool compositionchange(ref NG_Cal.GasPropsSTRUCT ptAGA10) { double dMixID = 0.0; int i; // generate the numerical ID for the composition for (i = 0; i < NG_Cal.NUMBEROFCOMPONENTS; i++) dMixID += ((i + 2) * ptAGA10.adMixture[i]); // update the history variable, if different from previous if (dMixID != dOldMixID) { dOldMixID = dMixID; return true; } else { return false; } }// compositionchange() /************************************************************************** *Function:Run() *Arguments:GasPropsSTRUCT * *Returns:void *Purpose:public method to coordinate and run the full calc sequence *Revisions: **************************************************************************/ public void Run(ref NG_Cal.GasPropsSTRUCT ptAGA10) { int i; // Check for gas composition change ptAGA10.bForceUpdate = (ptAGA10.bForceUpdate || compositionchange(ref ptAGA10)); // assign component IDs and values if (ptAGA10.bForceUpdate) { iNCC = -1; for (i = 0; i < NG_Cal.NUMBEROFCOMPONENTS; i++) { if (ptAGA10.adMixture[i] > 0.0) { iNCC = iNCC + 1; aiCID[iNCC] = i; dXi[iNCC] = ptAGA10.adMixture[i]; } } iNCC = iNCC + 1; //calculate composition dependent quantities; ported from original //FORTRAN functions paramdl() and chardl() paramdl(); chardl(ref ptAGA10); } //evaluate T & P dependent parms at base pressure and temperature, //but only if necessary if ((Math.Abs(ptAGA10.dPb - dOldPb) > NG_Cal.P_CHG_TOL) || (Math.Abs(ptAGA10.dTb - dOldTb) > NG_Cal.T_CHG_TOL) || (ptAGA10.bForceUpdate)) { dP = ptAGA10.dPb * 1.0e-6; // AGA 8 uses MPa internally dT = ptAGA10.dTb; //calculate temperature dependent parms temp(); //determine molar density ddetail(ref ptAGA10); ptAGA10.dDb = dRho; //determine compressibility ptAGA10.dZb = zdetail(dRho); // calculate mass density dRhoTP = (dP * ptAGA10.dMrx) / (ptAGA10.dZb * NG_Cal.RGASKJ * dT); //calculate relative density relativedensity(ref ptAGA10); //copy density to data structure member ptAGA10.dRhob = dRhoTP; //update history and clear the ForceUpdate flag dOldTb = ptAGA10.dTb; dOldPb = ptAGA10.dPb; ptAGA10.bForceUpdate = true; } //repeat the process using flowing conditions //begin by loading P & T from data structure //AGA 8 uses MPa internally; converted from Pa here dP = ptAGA10.dPf * 1.0e-6; dT = ptAGA10.dTf; //check whether to calculate temperature dependent parms if ((Math.Abs(ptAGA10.dTf - dOldTf) > NG_Cal.T_CHG_TOL) || (ptAGA10.bForceUpdate)) { //if temperature has changed, we must follow through temp(); //force ForceUpdate flag to true ptAGA10.bForceUpdate = true; } // check whether to calculate other parms if ((Math.Abs(ptAGA10.dPf - dOldPf) > NG_Cal.P_CHG_TOL) || (ptAGA10.bForceUpdate)) { //determine molar density ddetail(ref ptAGA10); ptAGA10.dDf = dRho; //determine compressibility ptAGA10.dZf = zdetail(dRho); //calculate mass density dRhoTP = (dP * ptAGA10.dMrx) / (ptAGA10.dZf * NG_Cal.RGASKJ * dT); //copy density to data structure member ptAGA10.dRhof = dRhoTP; //update history dOldTf = ptAGA10.dTf; dOldPf = ptAGA10.dPf; } //calculate legacy factor Fpv //NOTE: as implemented here, Fpv is not constrained to 14.73 psi and 60F if ((ptAGA10.dZb > 0.0) && (ptAGA10.dZf > 0.0)) { ptAGA10.dFpv = Math.Sqrt(ptAGA10.dZb / ptAGA10.dZf); } else //if either Zb or Zf is zero at this point, we have a serious unexpected problem { ptAGA10.dFpv = ptAGA10.dZb = ptAGA10.dZf = 0.0; ptAGA10.lStatus = NG_Cal.GENERAL_CALCULATION_FAILURE; } //we are now up to date; toggle off the update flag ptAGA10.bForceUpdate = false; }// Run() /************************************************************************** *Function:table() *Arguments:void *Returns:void *Purpose:builds tables of constants *Revisions: **************************************************************************/ //Tables 4 and 6 are filled only during object initialization. // //component ID's, mapped to each species supported by AGA Report#8 //1- methane8- hydrogen15- n-hexane //2- nitrogen9- carbon monoxide16- n-heptane //3- carbon dioxide10- oxygen17- n-octane //4- ethane11- i-butane18- n-nonane //5- propane12- n-butane19- n-decane //6- water13- i-pentane20- helium //7- hydrogen sulfide14- n-pentane21- argon public void table() { int j, k; // 58 constants from table 4 - column A(n) adAn[0] = 0.153832600;adAn[1] = 1.341953000; adAn[2] = -2.998583000; adAn[3] = -0.048312280; adAn[4] = 0.375796500; adAn[5] = -1.589575000; adAn[6] = -0.053588470; adAn[7] = 0.886594630; adAn[8] = -0.710237040; adAn[9] = -1.471722000; adAn[10] = 1.321850350; adAn[11] = -0.786659250; adAn[12] = 2.29129E-09; adAn[13] = 0.157672400; adAn[14] = -0.436386400; adAn[15] = -0.044081590; adAn[16] = -0.003433888; adAn[17] = 0.032059050; adAn[18] = 0.024873550; adAn[19] = 0.073322790; adAn[20] = -0.001600573; adAn[21] = 0.642470600; adAn[22] = -0.416260100; adAn[23] = -0.066899570; adAn[24] = 0.279179500; adAn[25] = -0.696605100; adAn[26] = -0.002860589; adAn[27] = -0.008098836; adAn[28] = 3.150547000; adAn[29] = 0.007224479; adAn[30] = -0.705752900; adAn[31] = 0.534979200; adAn[32] = -0.079314910; adAn[33] = -1.418465000; adAn[34] = -5.99905E-17; adAn[35] = 0.105840200; adAn[36] = 0.034317290; adAn[37] = -0.007022847; adAn[38] = 0.024955870; adAn[39] = 0.042968180; adAn[40] = 0.746545300; adAn[41] = -0.291961300; adAn[42] = 7.294616000; adAn[43] = -9.936757000; adAn[44] = -0.005399808; adAn[45] = -0.243256700; adAn[46] = 0.049870160; adAn[47] = 0.003733797; adAn[48] = 1.874951000; adAn[49] = 0.002168144; adAn[50] = -0.658716400; adAn[51] = 0.000205518; adAn[52] = 0.009776195; adAn[53] = -0.020487080; adAn[54] = 0.015573220; adAn[55] = 0.006862415; adAn[56] = -0.001226752; adAn[57] = 0.002850908; // 58 constants from table 4 - column Un adUn[0] = 0.0; adUn[1] = 0.5; adUn[2] = 1.0; adUn[3] = 3.5; adUn[4] = -0.5; adUn[5] = 4.5; adUn[6] = 0.5; adUn[7] = 7.5; adUn[8] = 9.5; adUn[9] = 6.0; adUn[10] = 12.0; adUn[11] = 12.5; adUn[12] = -6.0; adUn[13] = 2.0; adUn[14] = 3.0; adUn[15] = 2.0; adUn[16] = 2.0; adUn[17] = 11.0; adUn[18] = -0.5; adUn[19] = 0.5; adUn[20] = 0.0; adUn[21] = 4.0; adUn[22] = 6.0; adUn[23] = 21.0; adUn[24] = 23.0; adUn[25] = 22.0; adUn[26] = -1.0; adUn[27] = -0.5; adUn[28] = 7.0; adUn[29] = -1.0; adUn[30] = 6.0; adUn[31] = 4.0; adUn[32] = 1.0; adUn[33] = 9.0; adUn[34] = -13.0; adUn[35] = 21.0; adUn[36] = 8.0; adUn[37] = -0.5; adUn[38] = 0.0; adUn[39] = 2.0; adUn[40] = 7.0; adUn[41] = 9.0; adUn[42] = 22.0; adUn[43] = 23.0; adUn[44] = 1.0; adUn[45] = 9.0; adUn[46] = 3.0; adUn[47] = 8.0; adUn[48] = 23.0; adUn[49] = 1.5; adUn[50] = 5.0; adUn[51] = -0.5; adUn[52] = 4.0; adUn[53] = 7.0; adUn[54] = 3.0; adUn[55] = 0.0; adUn[56] = 1.0; adUn[57] = 0.0; //Most of the tables are filled with 1.0 or 0.0 //It is up to us to set non-zero values for (j = 0; j < NG_Cal.NUMBEROFCOMPONENTS; j++) { for (k = j; k < NG_Cal.NUMBEROFCOMPONENTS; k++) { adTable6Eij[j, k] = 1.0; adTable6Uij[j, k] = 1.0; adTable6Kij[j, k] = 1.0; adTable6Gij[j, k] = 1.0; } } //Lnsert the 132 items of non-zero and non-1.0 data //This looks more cumbersome than it is, considering table 6 has 1764 members adTable6Eij[0, 1] = 0.971640; adTable6Eij[0, 2] = 0.960644; adTable6Eij[0, 4] = 0.994635; adTable6Eij[0, 5] = 0.708218; adTable6Eij[0, 6] = 0.931484; adTable6Eij[0, 7] = 1.170520; adTable6Eij[0, 8] = 0.990126; adTable6Eij[0, 10] = 1.019530; adTable6Eij[0, 11] = 0.989844; adTable6Eij[0, 12] = 1.002350; adTable6Eij[0, 13] = 0.999268; adTable6Eij[0, 14] = 1.107274; adTable6Eij[0, 15] = 0.880880; adTable6Eij[0, 16] = 0.880973; adTable6Eij[0, 17] = 0.881067; adTable6Eij[0, 18] = 0.881161; adTable6Eij[1, 2] = 1.022740; adTable6Eij[1, 3] = 0.970120; adTable6Eij[1, 4] = 0.945939; adTable6Eij[1, 5] = 0.746954; adTable6Eij[1, 6] = 0.902271; adTable6Eij[1, 7] = 1.086320; adTable6Eij[1, 8] = 1.005710; adTable6Eij[1, 9] = 1.021000; adTable6Eij[1, 10] = 0.946914; adTable6Eij[1, 11] = 0.973384; adTable6Eij[1, 12] = 0.959340; adTable6Eij[1, 13] = 0.945520; adTable6Eij[2, 3] = 0.925053; adTable6Eij[2, 4] = 0.960237; adTable6Eij[2, 5] = 0.849408; adTable6Eij[2, 6] = 0.955052; adTable6Eij[2, 7] = 1.281790; adTable6Eij[2, 8] = 1.500000; adTable6Eij[2, 10] = 0.906849; adTable6Eij[2, 11] = 0.897362; adTable6Eij[2, 12] = 0.726255; adTable6Eij[2, 13] = 0.859764; adTable6Eij[2, 14] = 0.855134; adTable6Eij[2, 15] = 0.831229; adTable6Eij[2, 16] = 0.808310; adTable6Eij[2, 17] = 0.786323; adTable6Eij[2, 18] = 0.765171; adTable6Eij[3, 4] = 1.022560; adTable6Eij[3, 5] = 0.693168; adTable6Eij[3, 6] = 0.946871; adTable6Eij[3, 7] = 1.164460; adTable6Eij[3, 11] = 1.013060; adTable6Eij[3, 13] = 1.005320; adTable6Eij[4, 7] = 1.034787; adTable6Eij[4, 11] = 1.004900; adTable6Eij[6, 14] = 1.008692; adTable6Eij[6, 15] = 1.010126; adTable6Eij[6, 16] = 1.011501; adTable6Eij[6, 17] = 1.012821; adTable6Eij[6, 18] = 1.014089; adTable6Eij[7, 8] = 1.100000; adTable6Eij[7, 10] = 1.300000; adTable6Eij[7, 11] = 1.300000; adTable6Uij[0, 1] = 0.886106; adTable6Uij[0, 2] = 0.963827; adTable6Uij[0, 4] = 0.990877; adTable6Uij[0, 6] = 0.736833; adTable6Uij[0, 7] = 1.156390; adTable6Uij[0, 11] = 0.992291; adTable6Uij[0, 13] = 1.003670; adTable6Uij[0, 14] = 1.302576; adTable6Uij[0, 15] = 1.191904; adTable6Uij[0, 16] = 1.205769; adTable6Uij[0, 17] = 1.219634; adTable6Uij[0, 18] = 1.233498; adTable6Uij[1, 2] = 0.835058; adTable6Uij[1, 3] = 0.816431; adTable6Uij[1, 4] = 0.915502; adTable6Uij[1, 6] = 0.993476; adTable6Uij[1, 7] = 0.408838; adTable6Uij[1, 11] = 0.993556; adTable6Uij[2, 3] = 0.969870; adTable6Uij[2, 6] = 1.045290; adTable6Uij[2, 8] = 0.900000; adTable6Uij[2, 14] = 1.066638; adTable6Uij[2, 15] = 1.077634; adTable6Uij[2, 16] = 1.088178; adTable6Uij[2, 17] = 1.098291; adTable6Uij[2, 18] = 1.108021; adTable6Uij[3, 4] = 1.065173; adTable6Uij[3, 6] = 0.971926; adTable6Uij[3, 7] = 1.616660; adTable6Uij[3, 10] = 1.250000; adTable6Uij[3, 11] = 1.250000; adTable6Uij[3, 12] = 1.250000; adTable6Uij[3, 13] = 1.250000; adTable6Uij[6, 14] = 1.028973; adTable6Uij[6, 15] = 1.033754; adTable6Uij[6, 16] = 1.038338; adTable6Uij[6, 17] = 1.042735; adTable6Uij[6, 18] = 1.046966; adTable6Kij[0, 1] = 1.003630; adTable6Kij[0, 2] = 0.995933; adTable6Kij[0, 4] = 1.007619; adTable6Kij[0, 6] = 1.000080; adTable6Kij[0, 7] = 1.023260; adTable6Kij[0, 11] = 0.997596; adTable6Kij[0, 13] = 1.002529; adTable6Kij[0, 14] = 0.982962; adTable6Kij[0, 15] = 0.983565; adTable6Kij[0, 16] = 0.982707; adTable6Kij[0, 17] = 0.981849; adTable6Kij[0, 18] = 0.980991; adTable6Kij[1, 2] = 0.982361; adTable6Kij[1, 3] = 1.007960; adTable6Kij[1, 6] = 0.942596; adTable6Kij[1, 7] = 1.032270; adTable6Kij[2, 3] = 1.008510; adTable6Kij[2, 6] = 1.007790; adTable6Kij[2, 14] = 0.910183; adTable6Kij[2, 15] = 0.895362; adTable6Kij[2, 16] = 0.881152; adTable6Kij[2, 17] = 0.867520; adTable6Kij[2, 18] = 0.854406; adTable6Kij[3, 4] = 0.986893; adTable6Kij[3, 6] = 0.999969; adTable6Kij[3, 7] = 1.020340; adTable6Kij[6, 14] = 0.968130; adTable6Kij[6, 15] = 0.962870; adTable6Kij[6, 16] = 0.957828; adTable6Kij[6, 17] = 0.952441; adTable6Kij[6, 18] = 0.948338; adTable6Gij[0, 2] = 0.807653; adTable6Gij[0, 7] = 1.957310; adTable6Gij[1, 2] = 0.982746; adTable6Gij[2, 3] = 0.370296; adTable6Gij[2, 5] = 1.673090; }// table() /************************************************************************** *Function:paramdl() *Arguments:void *Returns:void *Purpose:sets up characterization & binary interaction parameters *Revisions: **************************************************************************/ public void paramdl() { int j, k; // table 5 parameters; declared locally to this function double[] adTable5Mri = new double[NG_Cal.NUMBEROFCOMPONENTS] { 16.0430, 28.0135, 44.0100, 30.0700, 44.0970, 18.0153, 34.0820, 2.0159, 28.0100, 31.9988, 58.1230, 58.1230, 72.1500, 72.1500, 86.1770, 100.2040, 114.2310, 128.2580, 142.2850, 4.0026, 39.9480 }; double[] adTable5Ei = new double[NG_Cal.NUMBEROFCOMPONENTS] { 151.318300, 99.737780, 241.960600, 244.166700, 298.118300, 514.015600, 296.355000, 26.957940, 105.534800, 122.766700, 324.068900, 337.638900, 365.599900, 370.682300, 402.636293, 427.722630, 450.325022, 470.840891, 489.558373, 2.610111, 119.629900 }; double[] adTable5Ki = new double[NG_Cal.NUMBEROFCOMPONENTS] { 0.4619255, 0.4479153, 0.4557489, 0.5279209, 0.5837490, 0.3825868, 0.4618263, 0.3514916, 0.4533894, 0.4186954, 0.6406937, 0.6341423, 0.6738577, 0.6798307, 0.7175118, 0.7525189, 0.7849550, 0.8152731, 0.8437826, 0.3589888, 0.4216551 }; double[] adTable5Gi = new double[NG_Cal.NUMBEROFCOMPONENTS] { 0.000000, 0.027815, 0.189065, 0.079300, 0.141239, 0.332500, 0.088500, 0.034369, 0.038953, 0.021000, 0.256692, 0.281835, 0.332267, 0.366911, 0.289731, 0.337542, 0.383381, 0.427354, 0.469659, 0.000000, 0.000000 }; //most of the table 5 parameters are zero for (j = 0; j < NG_Cal.NUMBEROFCOMPONENTS; j++) { adTable5Qi[j] = 0.0; adTable5Fi[j] = 0.0; adTable5Si[j] = 0.0; adTable5Wi[j] = 0.0; } //a small number of exceptions adTable5Qi[2] = 0.690000; adTable5Qi[5] = 1.067750; adTable5Qi[6] = 0.633276; adTable5Fi[7] = 1.0000; adTable5Si[5] = 1.5822; adTable5Si[6] = 0.3900; adTable5Wi[5] = 1.0000; // setup characterization parameters for non-zero components for (j = iNCC - 1; j >= 0; j--) { dMri[j] = adTable5Mri[aiCID[j]]; dKi[j] = adTable5Ki[aiCID[j]]; } for (j = 0; j < iNCC; j++) { dGi[j] = adTable5Gi[aiCID[j]]; dEi[j] = adTable5Ei[aiCID[j]]; } for (j = 0; j < iNCC; j++) { dQi[j] = adTable5Qi[aiCID[j]]; dFi[j] = 0.0; if (aiCID[j] == 7) dFi[j] = adTable5Fi[7]; dSi[j] = adTable5Si[aiCID[j]]; dWi[j] = adTable5Wi[aiCID[j]]; } // Binary interaction parameters for arrays: eij, kij, wij, uij for (j = 0; j < iNCC; j++) { for (k = j; k < iNCC; k++) { dUij[j, k] = adTable6Uij[aiCID[j], aiCID[k]]; dKij[j, k] = adTable6Kij[aiCID[j], aiCID[k]]; dEij[j, k] = adTable6Eij[aiCID[j], aiCID[k]]; dGij[j, k] = adTable6Gij[aiCID[j], aiCID[k]]; } } }// paramdl() /************************************************************************** *Function:chardl() *Arguments:GasPropsSTRUCT * *Returns:void *Purpose:computes composition-dependent quantities *Revisions: **************************************************************************/ public void chardl(ref NG_Cal.GasPropsSTRUCT ptAGA10) { //variables local to function int i, j; double tmfrac, k5p0, k2p5, u5p0, u2p5, q1p0; double Xij, Eij, Gij, e0p5, e2p0, e3p0, e3p5, e4p5, e6p0; double e7p5, e9p5, e12p0, e12p5; double e11p0, s3; //normalize mole fractions and calculate molar mass tmfrac = 0.0; for (j = 0; j < iNCC; j++) { tmfrac = tmfrac + dXi[j]; } for (j = 0; j < iNCC; j++) { dXi[j] = dXi[j] / tmfrac; } // reset virial coefficients for (j = 0; j < 18; j++) { adBcoef[j] = 0.0; } // initialize a key subset of the local variables k5p0 = 0.0; k2p5 = 0.0; u5p0 = 0.0; u2p5 = 0.0; dW = 0.0; q1p0 = 0.0; dF = 0.0; // calculate gas molecular weight ptAGA10.dMrx = 0.0; for (j = 0; j < iNCC; j++) { ptAGA10.dMrx = ptAGA10.dMrx + dXi[j] * dMri[j]; } // calculate the composition-dependent quantities, applying a nested loop for (i = 0; i < iNCC; i++) { k2p5 = k2p5 + dXi[i] * dKi[i] * dKi[i] * Math.Sqrt(dKi[i]); u2p5 = u2p5 + dXi[i] * dEi[i] * dEi[i] * Math.Sqrt(dEi[i]); dW = dW + dXi[i] * dGi[i]; q1p0 = q1p0 + dXi[i] * dQi[i]; dF = dF + dXi[i] * dXi[i] * dFi[i]; for (j = i; j < iNCC; j++) { if (i != j) Xij = 2.0 * dXi[i] * dXi[j]; else Xij = dXi[i] * dXi[j]; // proceed while skipping interaction terms which equal 1.0 if (dKij[i, j] != 1.0) k5p0 += Xij * (Math.Pow(dKij[i, j], 5.0) - 1.0) * Math.Pow((Math.Pow(dKi[i], 5.0) * Math.Pow(dKi[j], 5.0)), 0.5); if (dUij[i, j] != 1.0) u5p0 += Xij * (Math.Pow(dUij[i, j], 5.0) - 1.0) * Math.Pow((Math.Pow(dEi[i], 5.0) * Math.Pow(dEi[j], 5.0)), 0.5); if (dGij[i, j] != 1.0) dW += Xij * (dGij[i, j] - 1.0) * ((dGi[i] + dGi[j]) / 2.0); // calculate terms required for second virial coefficient, B Eij = dEij[i, j] * Math.Sqrt(dEi[i] * dEi[j]); Gij = dGij[i, j] * (dGi[i] + dGi[j]) / 2.0; e0p5 = Math.Sqrt(Eij); e2p0 = Eij * Eij; e3p0 = Eij * e2p0; e3p5 = e3p0 * e0p5; e4p5 = Eij * e3p5; e6p0 = e3p0 * e3p0; e11p0 = e4p5 * e4p5 * e2p0; e7p5 = e4p5 * Eij * e2p0; e9p5 = e7p5 * e2p0; e12p0 = e11p0 * Eij; e12p5 = e12p0 * e0p5; s3 = Xij * Math.Pow((Math.Pow(dKi[i], 3.0) * Math.Pow(dKi[j], 3)), 0.5); adBcoef[0] = adBcoef[0] + s3; adBcoef[1] = adBcoef[1] + s3 * e0p5; adBcoef[2] = adBcoef[2] + s3 * Eij; adBcoef[3] = adBcoef[3] + s3 * e3p5; adBcoef[4] = adBcoef[4] + s3 * Gij / e0p5; adBcoef[5] = adBcoef[5] + s3 * Gij * e4p5; adBcoef[6] = adBcoef[6] + s3 * dQi[i] * dQi[j] * e0p5; adBcoef[7] = adBcoef[7] + s3 * dSi[i] * dSi[j] * e7p5; adBcoef[8] = adBcoef[8] + s3 * dSi[i] * dSi[j] * e9p5; adBcoef[9] = adBcoef[9] + s3 * dWi[i] * dWi[j] * e6p0; adBcoef[10] = adBcoef[10] + s3 * dWi[i] * dWi[j] * e12p0; adBcoef[11] = adBcoef[11] + s3 * dWi[i] * dWi[j] * e12p5; adBcoef[12] = adBcoef[12] + s3 * dFi[i] * dFi[j] / e6p0; adBcoef[13] = adBcoef[13] + s3 * e2p0; adBcoef[14] = adBcoef[14] + s3 * e3p0; adBcoef[15] = adBcoef[15] + s3 * dQi[i] * dQi[j] * e2p0; adBcoef[16] = adBcoef[16] + s3 * e2p0; adBcoef[17] = adBcoef[17] + s3 * e11p0; } } //grab the first 18 constants from table 4, completing Bnij for (i = 0; i < 18; i++) adBcoef[i] *= adAn[i]; //final products of chardl are mixture size parameter K, energy parameter U, //and quadrupole parameter Q dKp3 = Math.Pow((k5p0 + k2p5 * k2p5), 0.6); dU = Math.Pow((u5p0 + u2p5 * u2p5), 0.2); dQp2 = q1p0 * q1p0; }// chardl() /************************************************************************** *Function:bvir() *Arguments:void *Returns:void *Purpose:computes 2nd virial coefficient & partial derivs thereof *Revisions: **************************************************************************/ public void bvir() { //variables local to function double t0p5, t2p0, t3p0, t3p5, t4p5, t6p0, t11p0; double t7p5, t9p5, t12p0, t12p5; double t1p5, t4p0; double[] Bx = new double[18]; int i; //reset B and partial devivatives to 0.0 dB = ddBdT = dd2BdT2 = 0.0; //pre-calculate Math .Powers of T t0p5 = Math.Sqrt(dT); t2p0 = dT * dT; t3p0 = dT * t2p0; t3p5 = t3p0 * t0p5; t4p5 = dT * t3p5; t6p0 = t3p0 * t3p0; t11p0 = t4p5 * t4p5 * t2p0; t7p5 = t6p0 * dT * t0p5; t9p5 = t7p5 * t2p0; t12p0 = t9p5 * t0p5 * t2p0; t12p5 = t12p0 * t0p5; t1p5 = dT * t0p5; t4p0 = t2p0 * t2p0; //coefficients for B Bx[0] = adBcoef[0]; Bx[1] = adBcoef[1] / t0p5; Bx[2] = adBcoef[2] / dT; Bx[3] = adBcoef[3] / t3p5; Bx[4] = adBcoef[4] * t0p5; Bx[5] = adBcoef[5] / t4p5; Bx[6] = adBcoef[6] / t0p5; Bx[7] = adBcoef[7] / t7p5; Bx[8] = adBcoef[8] / t9p5; Bx[9] = adBcoef[9] / t6p0; Bx[10] = adBcoef[10] / t12p0; Bx[11] = adBcoef[11] / t12p5; Bx[12] = adBcoef[12] * t6p0; Bx[13] = adBcoef[13] / t2p0; Bx[14] = adBcoef[14] / t3p0; Bx[15] = adBcoef[15] / t2p0; Bx[16] = adBcoef[16] / t2p0; Bx[17] = adBcoef[17] / t11p0; //sum up the pieces for second virial coefficient, B for (i = 0; i < 18; i++) { dB += Bx[i]; } //calculate terms for first derivative of B, wrt T for (i = 0; i < 18; i++) { if (adUn[i] != 0) Bx[i] *= adUn[i]; } //sum up the pieces of first derivative of B //note div by dT; changes exponent of T for (i = 0; i < 18; i++) { if (adUn[i] != 0) ddBdT += Bx[i] / dT; } //sign change here ddBdT = -ddBdT; //calculate terms for second derivative of B, wrt T for (i = 0; i < 18; i++) { if (adUn[i] != 0 && adUn[i] != -1.0) Bx[i] *= (adUn[i] + 1.0); } //sum up the pieces of second derivative of B //note division by dT, thereby changing the exponent of T //loop will ignore Bx[0] which is = 0.0 for (i = 0; i < 18; i++) { if (adUn[i] != 0 && adUn[i] != -1.0) dd2BdT2 += Bx[i] / t2p0; } }// bvir() /************************************************************************** *Function:temp() *Arguments:void *Returns:void *Purpose:computes temperature-dependent quantities *Revisions: **************************************************************************/ public void temp() { //Note: this function was ported from the AGA Report No.8 FORTRAN listing, //retaining as much of the original content as possible //variables local to function double tr0p5, tr1p5, tr2p0, tr3p0, tr4p0, tr5p0, tr6p0; double tr7p0, tr8p0, tr9p0, tr11p0, tr13p0, tr21p0; double tr22p0, tr23p0, tr; /*calculate second virial coefficient B*/ bvir(); //calculate adFn(12) through adFn(57) //adFn(0)-adFn(11) do not contribute to csm terms tr = dT / (dU); tr0p5 = Math.Sqrt(tr); tr1p5 = tr * tr0p5; tr2p0 = tr * tr; tr3p0 = tr * tr2p0; tr4p0 = tr * tr3p0; tr5p0 = tr * tr4p0; tr6p0 = tr * tr5p0; tr7p0 = tr * tr6p0; tr8p0 = tr * tr7p0; tr9p0 = tr * tr8p0; tr11p0 = tr6p0 * tr5p0; tr13p0 = tr6p0 * tr7p0; tr21p0 = tr9p0 * tr9p0 * tr3p0; tr22p0 = tr * tr21p0; tr23p0 = tr * tr22p0; adFn[12] = adAn[12] * dF * tr6p0; adFn[13] = adAn[13] / tr2p0; adFn[14] = adAn[14] / tr3p0; adFn[15] = adAn[15] * dQp2 / tr2p0; adFn[16] = adAn[16] / tr2p0; adFn[17] = adAn[17] / tr11p0; adFn[18] = adAn[18] * tr0p5; adFn[19] = adAn[19] / tr0p5; adFn[20] = adAn[20]; adFn[21] = adAn[21] / tr4p0; adFn[22] = adAn[22] / tr6p0; adFn[23] = adAn[23] / tr21p0; adFn[24] = adAn[24] * dW / tr23p0; adFn[25] = adAn[25] * dQp2 / tr22p0; adFn[26] = adAn[26] * dF * tr; adFn[27] = adAn[27] * dQp2 * tr0p5; adFn[28] = adAn[28] * dW / tr7p0; adFn[29] = adAn[29] * dF * tr; adFn[30] = adAn[30] / tr6p0; adFn[31] = adAn[31] * dW / tr4p0; adFn[32] = adAn[32] * dW / tr; adFn[33] = adAn[33] * dW / tr9p0; adFn[34] = adAn[34] * dF * tr13p0; adFn[35] = adAn[35] / tr21p0; adFn[36] = adAn[36] * dQp2 / tr8p0; adFn[37] = adAn[37] * tr0p5; adFn[38] = adAn[38]; adFn[39] = adAn[39] / tr2p0; adFn[40] = adAn[40] / tr7p0; adFn[41] = adAn[41] * dQp2 / tr9p0; adFn[42] = adAn[42] / tr22p0; adFn[43] = adAn[43] / tr23p0; adFn[44] = adAn[44] / tr; adFn[45] = adAn[45] / tr9p0; adFn[46] = adAn[46] * dQp2 / tr3p0; adFn[47] = adAn[47] / tr8p0; adFn[48] = adAn[48] * dQp2 / tr23p0; adFn[49] = adAn[49] / tr1p5; adFn[50] = adAn[50] * dW / tr5p0; adFn[51] = adAn[51] * dQp2 * tr0p5; adFn[52] = adAn[52] / tr4p0; adFn[53] = adAn[53] * dW / tr7p0; adFn[54] = adAn[54] / tr3p0; adFn[55] = adAn[55] * dW; adFn[56] = adAn[56] / tr; adFn[57] = adAn[57] * dQp2; }// temp() /************************************************************************** *Function:ddetail() *Arguments:GasPropsSTRUCT * *Returns:void *Purpose:calculates density *Revisions: **************************************************************************/ //Note: this function was ported from the AGA Report No.8 FORTRAN listing, //retaining as much of the original content as possible public void ddetail(ref NG_Cal.GasPropsSTRUCT ptAGA10) { int imax, i; double epsp, epsr, epsmin; double x1, x2, x3, y1, y2, y3; double delx, delprv, delmin, delbis, xnumer, xdenom, sgndel; double y2my3, y3my1, y1my2, boundn; //initialize convergence tolerances imax = 150; epsp = 1.0e-6; epsr = 1.0e-6; epsmin = 1.0e-7; dRho = 0.0; //call subroutine braket to bracket density solution braket(ref ptAGA10); //check value of "lStatus" returned from subroutine braket if (ptAGA10.lStatus == NG_Cal.MAX_NUM_OF_ITERATIONS_EXCEEDED || ptAGA10.lStatus == NG_Cal.NEGATIVE_DENSITY_DERIVATIVE) { return; } //set up to start Brent's method //x is the independent variable, y the dependent variable //delx is the current iteration change in x //delprv is the previous iteration change in x x1 = dRhoL; x2 = dRhoH; y1 = dPRhoL - dP; y2 = dPRhoH - dP; delx = x1 - x2; delprv = delx; //solution is bracketed between x1 and x2 //a third point x3 is introduced for quadratic interpolation x3 = x1; y3 = y1; for (i = 0; i < imax; i++) { //y3 must be opposite in sign from y2 so solution between x2,x3 if (y2 * y3 > 0.0) { x3 = x1; y3 = y1; delx = x1 - x2; delprv = delx; } //y2 must be value of y closest to y=0.0, then x2new=x2old+delx if (Math.Abs(y3) < Math.Abs(y2)) { x1 = x2; x2 = x3; x3 = x1; y1 = y2; y2 = y3; y3 = y1; } //delmin is minimum allowed step size for unconverged iteration delmin = epsmin * Math.Abs(x2); //if procedure is not converging or if delprv is less than delmin //use bisection instead //delbis = 0.5d0*(x3 - x2) is the bisection delx delbis = 0.5 * (x3 - x2); // tests to select numerical method for current iteration if (Math.Abs(delprv) < delmin || Math.Abs(y1) < Math.Abs(y2)) { // use bisection delx = delbis; delprv = delbis; } else { if (x3 != x1) { // use inverse quadratic interpolation y2my3 = y2 - y3; y3my1 = y3 - y1; y1my2 = y1 - y2; xdenom = -(y1my2) * (y2my3) * (y3my1); xnumer = x1 * y2 * y3 * (y2my3) + x2 * y3 * y1 * (y3my1) + x3 * y1 * y2 * (y1my2) - x2 * xdenom; } else { // use inverse linear interpolation xnumer = (x2 - x1) * y2; xdenom = y1 - y2; } // before calculating delx check delx=xnumer/xdenom is not out of bounds if (2.0 * Math.Abs(xnumer) < Math.Abs(delprv * xdenom)) { // procedure converging, use interpolation delprv = delx; delx = xnumer / xdenom; } else { // procedure diverging, use bisection delx = delbis; delprv = delbis; } } // check for convergence if ((Math.Abs(y2) < epsp * dP) && (Math.Abs(delx) < epsr * Math.Abs(x2))) { dRho = x2 + delx; return; } //when unconverged, abs(delx) must be greater than delmin //minimum allowed magnitude of change in x2 is 1.0000009*delmin //sgndel, the sign of change in x2 is sign of delbis if (Math.Abs(delx) < delmin) { sgndel = delbis / Math.Abs(delbis); delx = 1.0000009 * sgndel * delmin; delprv = delx; } //final check to insure that new x2 is in range of old x2 and x3 //boundn is negative if new x2 is in range of old x2 and x3 boundn = delx * (x2 + delx - x3); if (boundn > 0.0) { // procedure stepping out of bounds, use bisection delx = delbis; delprv = delbis; } //relable variables for next iteration //x1new = x2old, y1new=y2old x1 = x2; y1 = y2; // next iteration values for x2, y2 x2 = x2 + delx; pdetail(x2); y2 = dPCalc - dP; } // ddetail: maximum number of iterations exceeded ptAGA10.lStatus = NG_Cal.MAX_NUM_OF_ITERATIONS_EXCEEDED; dRho = x2; }// ddetail() /************************************************************************** *Function:braket() *Arguments:GasPropsSTRUCT * *Returns:void *Purpose:brackets density solution *Revisions: **************************************************************************/ //Note: this function was ported from the AGA Report No.8 FORTRAN listing, //retaining as much of the original content as possible public void braket(ref NG_Cal.GasPropsSTRUCT ptAGA10) { //variables local to function int imax, it; double del, rhomax, videal; double rho1, rho2, p1, p2; //initialize imax = 200; rho1 = 0.0; p1 = 0.0; rhomax = 1.0 / dKp3; if (dT > 1.2593 * dU) rhomax = 20.0 * rhomax; videal = NG_Cal.RGASKJ * dT / dP; if (Math.Abs(dB) < (0.167 * videal)) { rho2 = 0.95 / (videal + dB); } else { rho2 = 1.15 / videal; } del = rho2 / 20.0; // start iterative density search loop for (it = 0; it < imax; it++) { if (rho2 > rhomax && ptAGA10.lStatus != NG_Cal.MAX_DENSITY_IN_BRAKET_EXCEEDED) { // density in braket exceeds maximum allowable density ptAGA10.lStatus = NG_Cal.MAX_DENSITY_IN_BRAKET_EXCEEDED; del = 0.01 * (rhomax - rho1) + (dP / (NG_Cal.RGASKJ * dT)) / 20.0; rho2 = rho1 + del; continue; } //calculate pressure p2 at density rho2 pdetail(rho2); p2 = dPCalc; //test value of p2 relative to p and relative to p1 if (p2 > dP) { //the density root is bracketed (p1

p) dRhoL = rho1; dPRhoL = p1; dRhoH = rho2; dPRhoH = p2; ptAGA10.lStatus = NG_Cal.NORMAL; return; } else if (p2 > p1) { if (ptAGA10.lStatus == NG_Cal.MAX_DENSITY_IN_BRAKET_EXCEEDED) del *= 2.0; rho1 = rho2; p1 = p2; rho2 = rho1 + del; continue; } else { //lStatus= NEGATIVE_DENSITY_DERIVATIVEindicates that //pressure has a negative density derivative, since p2 is less than //some previous pressure ptAGA10.lStatus = NG_Cal.NEGATIVE_DENSITY_DERIVATIVE; dRho = rho1; return; } } // maximum number of iterations exceeded if we fall through the bottom ptAGA10.lStatus = NG_Cal.MAX_NUM_OF_ITERATIONS_EXCEEDED; dRho = rho2; return; }// braket() /************************************************************************** *Function:pdetail() *Arguments:double *Returns:void *Purpose:calculates pressure, given D and T. Calls zdetail() *Revisions: **************************************************************************/ public void pdetail(double dD) { dPCalc = zdetail(dD) * dD * NG_Cal.RGASKJ * dT; }// pdetail() /************************************************************************** *Function:zdetail() *Arguments:double *Returns:void *Purpose:calculates compressibility *Revisions: **************************************************************************/ public double zdetail(double d) { // variables local to function double D1, D2, D3, D4, D5, D6, D7, D8, D9, exp1, exp2, exp3, exp4; // Math .Powers of reduced density D1 = dKp3 * d; D2 = D1 * D1; D3 = D2 * D1; D4 = D3 * D1; D5 = D4 * D1; D6 = D5 * D1; D7 = D6 * D1; D8 = D7 * D1; D9 = D8 * D1; exp1 = Math.Exp(-D1); exp2 = Math.Exp(-D2); exp3 = Math.Exp(-D3); exp4 = Math.Exp(-D4); // the following expression for Z was adopted from FORTRAN example in AGA8 dZ = 1.0 + dB * d + adFn[12] * D1 * (exp3 - 1.0 - 3.0 * D3 * exp3) + (adFn[13] + adFn[14] + adFn[15]) * D1 * (exp2 - 1.0 - 2.0 * D2 * exp2) + (adFn[16] + adFn[17]) * D1 * (exp4 - 1.0 - 4.0 * D4 * exp4) + (adFn[18] + adFn[19]) * D2 * 2.0 + (adFn[20] + adFn[21] + adFn[22]) * D2 * (2.0 - 2.0 * D2) * exp2 + (adFn[23] + adFn[24] + adFn[25]) * D2 * (2.0 - 4.0 * D4) * exp4 + adFn[26] * D2 * (2.0 - 4.0 * D4) * exp4 + adFn[27] * D3 * 3.0 + (adFn[28] + adFn[29]) * D3 * (3.0 - D1) * exp1 + (adFn[30] + adFn[31]) * D3 * (3.0 - 2.0 * D2) * exp2 + (adFn[32] + adFn[33]) * D3 * (3.0 - 3.0 * D3) * exp3 + (adFn[34] + adFn[35] + adFn[36]) * D3 * (3.0 - 4.0 * D4) * exp4 + (adFn[37] + adFn[38]) * D4 * 4.0 + (adFn[39] + adFn[40] + adFn[41]) * D4 * (4.0 - 2.0 * D2) * exp2 + (adFn[42] + adFn[43]) * D4 * (4.0 - 4.0 * D4) * exp4 + adFn[44] * D5 * 5.0 + (adFn[45] + adFn[46]) * D5 * (5.0 - 2.0 * D2) * exp2 + (adFn[47] + adFn[48]) * D5 * (5.0 - 4.0 * D4) * exp4 + adFn[49] * D6 * 6.0 + adFn[50] * D6 * (6.0 - 2.0 * D2) * exp2 + adFn[51] * D7 * 7.0 + adFn[52] * D7 * (7.0 - 2.0 * D2) * exp2 + adFn[53] * D8 * (8.0 - D1) * exp1 + (adFn[54] + adFn[55]) * D8 * (8.0 - 2.0 * D2) * exp2 + (adFn[56] + adFn[57]) * D9 * (9.0 - 2.0 * D2) * exp2; return dZ; }// zdetail() /************************************************************************** *Function:dZdT() *Arguments:double *Returns:double *Purpose:calculates the first partial derivative of Z wrt T *Revisions: **************************************************************************/ public double dZdT(double d) { //variables local to function double tmp; int i; double D1, D2, D3, D4, D5, D6, D7, D8, D9, exp1, exp2, exp3, exp4; //set up Math .Powers of reduced density D1 = dKp3 * d; D2 = D1 * D1; D3 = D2 * D1; D4 = D3 * D1; D5 = D4 * D1; D6 = D5 * D1; D7 = D6 * D1; D8 = D7 * D1; D9 = D8 * D1; exp1 = Math.Exp(-D1); exp2 = Math.Exp(-D2); exp3 = Math.Exp(-D3); exp4 = Math.Exp(-D4); // create terms uC*T^-(un+1) from coefficients we've already computed (An[n]) for (i = 12; i < 58; i++) { if (adUn[i] != 0 && adFn[i] != 0) { fx[i] = (adFn[i] * adUn[i] * D1) / dT; } else { fx[i] = 0.0; } } //initial part of equation ddZdT = d * ddBdT; //n=13 evaluates to zero except for hydrogen, for whom fn = 1 if (dF != 0) ddZdT += fx[12] - (fx[12] * (1.0 - 3.0 * D3) * exp3); tmp = (1.0 - 2.0 * D2) * exp2; ddZdT += (fx[13] - (fx[13] * tmp)); ddZdT += fx[14] - (fx[14] * tmp); ddZdT += fx[15] - (fx[15] * tmp); tmp = (1.0 - 4.0 * D4) * exp4; ddZdT += fx[16] - (fx[16] * tmp); ddZdT += fx[17] - (fx[17] * tmp); ddZdT = ddZdT - (fx[18] + fx[19]) * D1 * 2.0 - (fx[21] + fx[22]) * D1 * (2.0 - 2.0 * D2) * exp2 - (fx[23] + fx[24] + fx[25]) * D1 * (2.0 - 4.0 * D4) * exp4 - fx[26] * D1 * (2.0 - 4.0 * D4) * exp4 - fx[27] * D2 * 3.0 - (fx[28] + fx[29]) * D2 * (3.0 - D1) * exp1 - (fx[30] + fx[31]) * D2 * (3.0 - 2.0 * D2) * exp2 - (fx[32] + fx[33]) * D2 * (3.0 - 3.0 * D3) * exp3 - (fx[34] + fx[35] + fx[36]) * D2 * (3.0 - 4.0 * D4) * exp4 - fx[37] * D3 * 4.0 - (fx[39] + fx[40] + fx[41]) * D3 * (4.0 - 2.0 * D2) * exp2 - (fx[42] + fx[43]) * D3 * (4.0 - 4.0 * D4) * exp4 - fx[44] * D4 * 5.0 - (fx[45] + fx[46]) * D4 * (5.0 - 2.0 * D2) * exp2 - (fx[47] + fx[48]) * D4 * (5.0 - 4.0 * D4) * exp4 - fx[49] * D5 * 6.0 - fx[50] * D5 * (6.0 - 2.0 * D2) * exp2 - fx[51] * D6 * 7.0 - fx[52] * D6 * (7.0 - 2.0 * D2) * exp2 - fx[53] * D7 * (8.0 - D1) * exp1 - fx[54] * D7 * (8.0 - 2.0 * D2) * exp2 - fx[56] * D8 * (9.0 - 2.0 * D2) * exp2; return ddZdT; }// dDdT() /************************************************************************** *Function:d2ZdT2() *Arguments:double *Returns:double *Purpose:calculates the second partial derivative of Z wrt T *Revisions: **************************************************************************/ public double d2ZdT2(double d) { //variables local to function double tmp; int i; double D1, D2, D3, D4, D5, D6, D7, D8, D9, exp1, exp2, exp3, exp4; //set up Math .Powers of reduced density D1 = dKp3 * d; D2 = D1 * D1; D3 = D2 * D1; D4 = D3 * D1; D5 = D4 * D1; D6 = D5 * D1; D7 = D6 * D1; D8 = D7 * D1; D9 = D8 * D1; exp1 = Math.Exp(-D1); exp2 = Math.Exp(-D2); exp3 = Math.Exp(-D3); exp4 = Math.Exp(-D4); // create terms uC*T^-(un+1) from coefficients we've already computed (An[n]) for (i = 12; i < 58; i++) { if (adUn[i] != 0 && adFn[i] != 0) { fx[i] = (adFn[i] * D1 * adUn[i] * (adUn[i] + 1.0)) / (dT * dT); } else { fx[i] = 0.0; } } //initial part of equation dd2ZdT2 = d * dd2BdT2; //n=13 evaluates to zero except for hydrogen, for whom fn = 1 if (dF != 0) dd2ZdT2 += fx[12] - (fx[12] * (1.0 - 3.0 * D3) * exp3); tmp = (1.0 - 2.0 * D2) * exp2; dd2ZdT2 += -fx[13] + (fx[13] * tmp); dd2ZdT2 += -fx[14] + (fx[14] * tmp); dd2ZdT2 += -fx[15] + (fx[15] * tmp); tmp = (1.0 - 4.0 * D4) * exp4; dd2ZdT2 += -fx[16] + (fx[16] * tmp); dd2ZdT2 += -fx[17] + (fx[17] * tmp); dd2ZdT2 = dd2ZdT2 + (fx[18] + fx[19]) * D1 * 2.0 + (fx[21] + fx[22]) * D1 * (2.0 - 2.0 * D2) * exp2 + (fx[23] + fx[24] + fx[25]) * D1 * (2.0 - 4.0 * D4) * exp4 + fx[26] * D1 * (2.0 - 4.0 * D4) * exp4 + fx[27] * D2 * 3.0 + (fx[28] + fx[29]) * D2 * (3.0 - D1) * exp1 + (fx[30] + fx[31]) * D2 * (3.0 - 2.0 * D2) * exp2 + (fx[32] + fx[33]) * D2 * (3.0 - 3.0 * D3) * exp3 + (fx[34] + fx[35] + fx[36]) * D2 * (3.0 - 4.0 * D4) * exp4 + fx[37] * D3 * 4.0 + (fx[39] + fx[40] + fx[41]) * D3 * (4.0 - 2.0 * D2) * exp2 + (fx[42] + fx[43]) * D3 * (4.0 - 4.0 * D4) * exp4 + fx[44] * D4 * 5.0 + (fx[45] + fx[46]) * D4 * (5.0 - 2.0 * D2) * exp2 + (fx[47] + fx[48]) * D4 * (5.0 - 4.0 * D4) * exp4 + fx[49] * D5 * 6.0 + fx[50] * D5 * (6.0 - 2.0 * D2) * exp2 + fx[51] * D6 * 7.0 + fx[52] * D6 * (7.0 - 2.0 * D2) * exp2 + fx[53] * D7 * (8.0 - D1) * exp1 + fx[54] * D7 * (8.0 - 2.0 * D2) * exp2 + fx[56] * D8 * (9.0 - 2.0 * D2) * exp2; return dd2ZdT2; }// d2ZdT2() /************************************************************************** *Function:dZdD() *Arguments:double *Returns:double *Purpose:calculates the first partial derivative of Z wrt D *Revisions: **************************************************************************/ //For efficiency and continuity with AGA 8 code example, each term //is evaluated individually rather than through looping through tables. //Temporary storage is used to hold portions of complex equations and //to facilitate debugging. Additional speed optimization is possible. public double dZdD(double d) { double temp, temp1, temp2, temp3; int i; double D1, D2, D3, D4, D5, D6, D7, D8, D9, exp1, exp2, exp3, exp4; // set up Math .Powers of reduced density D1 = dKp3 * d; D2 = D1 * D1; D3 = D2 * D1; D4 = D3 * D1; D5 = D4 * D1; D6 = D5 * D1; D7 = D6 * D1; D8 = D7 * D1; D9 = D8 * D1; exp1 = Math.Exp(-D1); exp2 = Math.Exp(-D2); exp3 = Math.Exp(-D3); exp4 = Math.Exp(-D4); //create terms uC*T^-(un+1) from coefficients we've already computed (An[n]) for (i = 12; i < 58; i++) { fx[i] = adFn[i]; } //initial part of equation ddZdD = dB / dKp3; //evaluate all remaining terms, simplifying where possible //n=13 evaluates to zero except for hydrogen, for whom fn = 1 if (dF != 0) { temp1 = -9.0 * D3 * exp3; temp2 = (1.0 - 3.0 * D3) * exp3; temp3 = -temp2 * 3.0 * D6; temp = temp1 + temp2 + temp3; ddZdD += -fx[12] + fx[12] * temp; } //n = 14..16 temp1 = -4.0 * D2 * exp2; temp2 = (1.0 - 2.0 * D2) * exp2; temp3 = -temp2 * 2.0 * D2; temp = temp1 + temp2 + temp3; ddZdD += -fx[13] + fx[13] * temp; ddZdD += -fx[14] + fx[14] * temp; ddZdD += -fx[15] + fx[15] * temp; // n =17..18 temp1 = -16.0 * D4 * exp4; temp2 = (1.0 - 4.0 * D4) * exp4; temp3 = -temp2 * 4.0 * D4; temp = temp1 + temp2 + temp3; ddZdD += -fx[16] + fx[16] * temp; ddZdD += -fx[17] + fx[17] * temp; // n = 19..20 temp = 4.0 * D1; ddZdD += fx[18] * temp; ddZdD += fx[19] * temp; // n =21..23 temp1 = -4.0 * D3 * exp2; temp2 = (2.0 - 2.0 * D2) * 2.0 * D1 * exp2; temp3 = -temp2 * D2; temp = temp1 + temp2 + temp3; ddZdD += fx[20] * temp; ddZdD += fx[21] * temp; ddZdD += fx[22] * temp; // n =24..27 temp1 = -16.0 * D5 * exp4; temp2 = (2.0 - 4.0 * D4) * 2.0 * D1 * exp4; temp3 = -temp2 * 2.0 * D4; temp = temp1 + temp2 + temp3; ddZdD += fx[23] * temp; ddZdD += fx[24] * temp; ddZdD += fx[25] * temp; ddZdD += fx[26] * temp; // n =28 temp = 9.0 * D2; ddZdD += fx[27] * temp; // n =29..30 temp = -D3 * exp1 + (3.0 - D1) * 3.0 * D2 * exp1; temp -= (3.0 - D1) * D3 * exp1; ddZdD += fx[28] * temp; ddZdD += fx[29] * temp; // n =31..32 temp1 = -4.0 * D4 * exp2; temp2 = (3.0 - 2.0 * D2) * 3.0 * D2 * exp2; temp3 = -(3.0 - 2.0 * D2) * 2.0 * D4 * exp2; temp = temp1 + temp2 + temp3; ddZdD += fx[30] * temp; ddZdD += fx[31] * temp; // n =33..34 temp1 = -9.0 * D5 * exp3; temp2 = (3.0 - 3.0 * D3) * 3.0 * D2 * exp3; temp3 = -(3.0 - 3.0 * D3) * 3.0 * D5 * exp3; temp = temp1 + temp2 + temp3; ddZdD += fx[32] * temp; ddZdD += fx[33] * temp; // n =35..37 temp1 = -16.0 * D6 * exp4; temp2 = (3.0 - 4.0 * D4) * 3.0 * D2 * exp4; temp3 = -(3.0 - 4.0 * D4) * D6 * 4.0 * exp4; temp = temp1 + temp2 + temp3; ddZdD += fx[34] * temp; ddZdD += fx[35] * temp; ddZdD += fx[36] * temp; //n = 38..39 temp = 16.0 * D3; ddZdD += fx[37] * temp; ddZdD += fx[38] * temp; //n = 40..42 temp1 = -4.0 * D5 * exp2; temp2 = (4.0 - 2.0 * D2) * 4.0 * D3 * exp2; temp3 = -(4.0 - 2.0 * D2) * 2.0 * D5 * exp2; temp = temp1 + temp2 + temp3; ddZdD += fx[39] * temp; ddZdD += fx[40] * temp; ddZdD += fx[41] * temp; // n =43..44 temp = -16.0 * D7 * exp4 + (4.0 - 4.0 * D4) * 4.0 * D3 * exp4; temp -= (4.0 - 4.0 * D4) * D7 * 4.0 * exp4; ddZdD += fx[42] * temp; ddZdD += fx[43] * temp; // n =45 temp = 25.0 * D4; ddZdD += fx[44] * temp; // n =46..47 temp = -4.0 * D6 * exp2 + (5.0 - 2.0 * D2) * 5.0 * D4 * exp2; temp -= (5.0 - 2.0 * D2) * D6 * 2.0 * exp2; ddZdD += fx[45] * temp; ddZdD += fx[46] * temp; // n =48..49 temp = -16.0 * D8 * exp4 + (5.0 - 4.0 * D4) * 5.0 * D4 * exp4; temp -= (5.0 - 4.0 * D4) * D8 * 4.0 * exp4; ddZdD += fx[47] * temp; ddZdD += fx[48] * temp; // n =50 temp = 36.0 * D5; ddZdD += fx[49] * temp; // n =51 temp = -4.0 * D7 * exp2 + (6.0 - 2.0 * D2) * 6.0 * D5 * exp2; temp -= (6.0 - 2.0 * D2) * D7 * 2.0 * exp2; ddZdD += fx[50] * temp; // n =52 temp = 49.0 * D6; ddZdD += fx[51] * temp; // n =53 temp = -4.0 * D8 * exp2 + (7.0 - 2.0 * D2) * 7.0 * D6 * exp2; temp -= (7.0 - 2.0 * D2) * D8 * 2.0 * exp2; ddZdD += fx[52] * temp; // n =54 temp = -1.0 * D8 * exp1 + (8.0 - D1) * 8.0 * D7 * exp1; temp -= (8.0 - D1) * D8 * exp1; ddZdD += fx[53] * temp; // n =55..56 temp = -4.0 * D1 * D8 * exp2 + (8.0 - 2.0 * D2) * 8.0 * D7 * exp2; temp -= (8.0 - 2.0 * D2) * D8 * 2.0 * D1 * exp2; ddZdD += fx[54] * temp; ddZdD += fx[55] * temp; // n =57..58 temp = -4.0 * D2 * D8 * exp2 + (9.0 - 2.0 * D2) * 9.0 * D8 * exp2; temp -= (9.0 - 2.0 * D2) * D2 * D8 * 2.0 * exp2; ddZdD += fx[56] * temp; ddZdD += fx[57] * temp; ddZdD *= dKp3; return ddZdD; }// dZdD() /************************************************************************** *Function:relativedensity() *Arguments:GasPropsSTRUCT * *Returns:void *Purpose:calculates relative density via methods listed in AGA 8 *Revisions: **************************************************************************/ public void relativedensity(ref NG_Cal.GasPropsSTRUCT ptAGA10) { double dBX, dZa; const double dMWair = 28.96256; // calculate second virial coefficient for air dBX = -0.12527 + 5.91e-4 * ptAGA10.dTb - 6.62e-7 * ptAGA10.dTb * ptAGA10.dTb; // calculate compressibility of air dZa = 1.0 + (dBX * dP) / (NG_Cal.RGASKJ * ptAGA10.dTb); // calculate ideal gas and real gas relative densities ptAGA10.dRD_Ideal = ptAGA10.dMrx / dMWair; ptAGA10.dRD_Real = ptAGA10.dRD_Ideal * (dZa / ptAGA10.dZb); }// relativedensity() } }