ruoyi-api/ruoyi-ngtools/src/main/java/com/ruoyi/ngCalTools/service/DetailService.java

1388 lines
50 KiB
Java
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

package com.ruoyi.ngCalTools.service;
import com.ruoyi.ngCalTools.model.GasProps;
import com.ruoyi.ngCalTools.utils.GasConstants;
import org.springframework.stereotype.Service;
@Service
public class DetailService {
private static final double RGASKJ = 8.314510e-3;
// 成员变量转换
// 组件数量
private int iNCC;
// 组件ID数组长度为21
private int[] aiCID = new int[21];
// 五个历史变量,用于在重复计算时提高效率
// 上一次计算的混合物ID
double dOldMixID;
// 上一次计算的Pb值
double dOldPb;
// 上一次计算的Tb值
double dOldTb;
// 上一次计算的Pf值
double dOldPf;
// 上一次计算的Tf值
double dOldTf;
// 来自表4第1列的EOS参数
// 长度为58的数组adAn
double[] adAn = new double[58];
// 长度为58的数组adUn
double[] adUn = new double[58];
// 来自表5的特征参数
// 第i个组件的分子量
double[] dMri = new double[21];
// 第i个组件的特征能量参数
double[] dEi = new double[21];
// 第i个组件的尺寸参数 - m^3/kg-mol ^1/3
double[] dKi = new double[21];
// 取向参数
double[] dGi = new double[21];
// 四极矩参数
double[] dQi = new double[21];
// 高温参数
double[] dFi = new double[21];
// 偶极矩参数
double[] dSi = new double[21];
// 关联参数
double[] dWi = new double[21];
// 维里系数能量二元相互作用参数
double[][] dEij = new double[21][21];
// 共形能量的二元相互作用参数
double[][] dUij = new double[21][21];
// 尺寸的二元相互作用参数
double[][] dKij = new double[21][21];
// 取向的二元相互作用参数
double[][] dGij = new double[21][21];
// 表6常量
double[][] adTable6Eij = new double[21][21];
// 表6常量
double[][] adTable6Uij = new double[21][21];
// 表6常量
double[][] adTable6Kij = new double[21][21];
// 表6常量
double[][] adTable6Gij = new double[21][21];
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
// 组件i的摩尔分数数组长度为21
double[] dXi = new double[21];
// 由pdetail()方法计算得到的压力
double dPCalc;
// 当前温度
double dT;
// 当前压力
double dP;
// 在温度T和压力P下的摩尔密度
double dRhoTP;
// 第二维里系数B
double dB;
// 用于计算B的18个系数的数组
double[] adBcoef = new double[18];
// 密度系数的函数数组长度为58
double[] adFn = new double[58];
// 用于3个导数的修正系数数组长度为58
double[] fx = new double[58];
// 混合能量参数
double dU;
// 混合尺寸参数的三次方
double dKp3;
// 混合取向参数
double dW;
// 混合四极矩参数的平方
double dQp2;
// 高温参数
double dF;
// 摩尔密度
double dRho;
// 在braket函数中使用的低密度
double dRhoL;
// 在braket函数中使用的高密度
double dRhoH;
// 在braket函数中使用的低压
double dPRhoL;
// 在braket函数中使用的高压
double dPRhoH;
// 也用于高级流体性质计算的公共变量
// 当前压缩因子
public double dZ;
// Z对T的一阶偏导数
public double ddZdT;
// Z对T的二阶偏导数
public double dd2ZdT2;
// Z对摩尔密度的一阶偏导数
public double ddZdD;
// B对T的一阶偏导数
public double ddBdT;
// B对T的二阶偏导数
public double dd2BdT2;
// 构造函数
public DetailService()
{
//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 < GasConstants.NUMBEROFCOMPONENTS; i++) dXi[i] = 0;
// function table() populates tables of static constants
table();
}
public void run(GasProps gasProps) {
// 实现转换后的逻辑
int i;
// Check for gas composition change
gasProps.bForceUpdate = (gasProps.bForceUpdate || compositionChange(gasProps));
// assign component IDs and values
if (gasProps.bForceUpdate) {
iNCC = -1;
for (i = 0; i < GasConstants.NUMBEROFCOMPONENTS; i++) {
if (gasProps.adMixture[i] > 0.0) {
iNCC = iNCC + 1;
aiCID[iNCC] = i;
dXi[iNCC] = gasProps.adMixture[i];
}
}
iNCC = iNCC + 1;
//calculate composition dependent quantities; ported from original
//FORTRAN functions paramdl() and chardl()
paramdl();
chardl(gasProps);
}
//evaluate T & P dependent parms at base pressure and temperature,
//but only if necessary
if (Math.abs(gasProps.dPb - dOldPb) > GasConstants.P_CHG_TOL || Math.abs(gasProps.dTb - dOldTb) > GasConstants.T_CHG_TOL || gasProps.bForceUpdate) {
dP = gasProps.dPb * 1.0e-6; // AGA 8 uses MPa internally
dT = gasProps.dTb;
//calculate temperature dependent parms
temp();
//determine molar density
ddetail(gasProps);
gasProps.dDb = dRho;
//determine compressibility
gasProps.dZb = zdetail(dRho);
// calculate mass density
dRhoTP = (dP * gasProps.dMrx) / (gasProps.dZb * GasConstants.RGASKJ * dT);
//calculate relative density
relativedensity(gasProps);
//copy density to data structure member
gasProps.dRhob = dRhoTP;
//update history and clear the ForceUpdate flag
dOldTb = gasProps.dTb;
dOldPb = gasProps.dPb;
gasProps.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 = gasProps.dPf * 1.0e-6;
dT = gasProps.dTf;
//check whether to calculate temperature dependent parms
if (Math.abs(gasProps.dTf - dOldTf) > GasConstants.T_CHG_TOL || gasProps.bForceUpdate) {
//if temperature has changed, we must follow through
temp();
//force ForceUpdate flag to true
gasProps.bForceUpdate = true;
}
// check whether to calculate other parms
if (Math.abs(gasProps.dPf - dOldPf) > GasConstants.P_CHG_TOL || gasProps.bForceUpdate) {
//determine molar density
ddetail(gasProps);
gasProps.dDf = dRho;
//determine compressibility
gasProps.dZf = zdetail(dRho);
//calculate mass density
dRhoTP = (dP * gasProps.dMrx) / (gasProps.dZf * GasConstants.RGASKJ * dT);
//copy density to data structure member
gasProps.dRhof = dRhoTP;
//update history
dOldTf = gasProps.dTf;
dOldPf = gasProps.dPf;
}
//calculate legacy factor Fpv
//NOTE: as implemented here, Fpv is not constrained to 14.73 psi and 60F
if (gasProps.dZb > 0.0 && gasProps.dZf > 0.0) {
gasProps.dFpv = Math.sqrt(gasProps.dZb / gasProps.dZf);
} else {
//if either Zb or Zf is zero at this point, we have a serious unexpected problem
gasProps.dFpv = gasProps.dZb = gasProps.dZf = 0.0;
gasProps.lStatus = GasConstants.GENERAL_CALCULATION_FAILURE;
}
//we are now up to date; toggle off the update flag
gasProps.bForceUpdate = false;
// 其他计算逻辑...
}
private boolean compositionChange(GasProps gasProps) {
double dMixID = 0.0;
for (int i = 0; i < GasConstants.NUMBEROFCOMPONENTS; i++) {
dMixID += ((i + 2) * gasProps.adMixture[i]);
}if (dMixID != dOldMixID) {
dOldMixID = dMixID;
return true;
}
return false;
}
public void paramdl() {
int j, k;
// table 5 parameters; declared locally to this function
double[] adTable5Mri;
double[] adTable5Ei;
double[] adTable5Ki;
double[] adTable5Gi;
// 初始化adTable5Mri数组
adTable5Mri = new double[]{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};
// 初始化adTable5Ei数组
adTable5Ei = new double[]{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};
// 初始化adTable5Ki数组
adTable5Ki = new double[]{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};
// 初始化adTable5Gi数组
adTable5Gi = new double[]{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 < GasConstants.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]];
}
}
}
public void chardl(GasProps gasProps) {
//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
gasProps.dMrx = 0.0;
for (j = 0; j < iNCC; j++) {
gasProps.dMrx = gasProps.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;
}
// 其他方法转换...
public void table() {
int j, k;
GasProps gasProps;
// 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 < GasConstants.NUMBEROFCOMPONENTS; j++) {
for (k = j; k < GasConstants.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;
}
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;
}
}
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;
}
public void ddetail(GasProps gasProps) {
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(gasProps);
//check value of "lStatus" returned from subroutine braket
if (gasProps.lStatus == GasConstants.MAX_NUM_OF_ITERATIONS_EXCEEDED || gasProps.lStatus == GasConstants.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
gasProps.lStatus = GasConstants.MAX_NUM_OF_ITERATIONS_EXCEEDED;
dRho = x2;
}// ddetail()
public void braket(GasProps gasProps) {
//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 = GasConstants.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 && gasProps.lStatus != GasConstants.MAX_DENSITY_IN_BRAKET_EXCEEDED) {
// density in braket exceeds maximum allowable density
gasProps.lStatus = GasConstants.MAX_DENSITY_IN_BRAKET_EXCEEDED;
del = 0.01 * (rhomax - rho1) + (dP / (GasConstants.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 and p2>p)
dRhoL = rho1;
dPRhoL = p1;
dRhoH = rho2;
dPRhoH = p2;
gasProps.lStatus = GasConstants.NORMAL;
return;
} else if (p2 > p1) {
if (gasProps.lStatus == GasConstants.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
gasProps.lStatus = GasConstants.NEGATIVE_DENSITY_DERIVATIVE;
dRho = rho1;
return;
}
}
// maximum number of iterations exceeded if we fall through the bottom
gasProps.lStatus = GasConstants.MAX_NUM_OF_ITERATIONS_EXCEEDED;
dRho = rho2;
return;
}// braket()
public void pdetail(double dD) {
dPCalc = zdetail(dD) * dD * GasConstants.RGASKJ * dT;
}// pdetail()
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()
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;
}
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()
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;
}
public void relativedensity(GasProps gasProps) {
double dBX, dZa;
double dMWair = 28.96256;
dBX = -0.12527 + 5.91e-4 * gasProps.dTb - 6.62e-7 * gasProps.dTb * gasProps.dTb;
// calculate compressibility of air
dZa = 1.0 + (dBX * dP) / (GasConstants.RGASKJ * gasProps.dTb);
// calculate ideal gas and real gas relative densities
gasProps.dRD_Ideal = gasProps.dMrx / dMWair;
gasProps.dRD_Real = gasProps.dRD_Ideal * (dZa / gasProps.dZb);
}
}