diff options
| -rw-r--r-- | subsurface-core/gas-model.c | 98 |
1 files changed, 42 insertions, 56 deletions
diff --git a/subsurface-core/gas-model.c b/subsurface-core/gas-model.c index 420233612..81765e003 100644 --- a/subsurface-core/gas-model.c +++ b/subsurface-core/gas-model.c @@ -4,72 +4,58 @@ #include <stdlib.h> #include "dive.h" -/* - * Generic cubic polynomial - */ -static double cubic(double bar, const double coefficient[]) -{ - double x0 = 1.0, - x1 = bar, - x2 = x1*x1, - x3 = x2*x1; - - return x0 * coefficient[0] + - x1 * coefficient[1] + - x2 * coefficient[2] + - x3 * coefficient[3]; -} +/* "Virial minus one" - the virial cubic form without the initial 1.0 */ +#define virial_m1(C, x1, x2, x3) (C[0]*x1+C[1]*x2+C[2]*x3) /* - * Cubic least-square coefficients for O2/N2/He based on data from + * Cubic virial least-square coefficients for O2/N2/He based on data from * * PERRY’S CHEMICAL ENGINEERS’ HANDBOOK SEVENTH EDITION * * with the lookup and curve fitting by Lubomir. * + * The "virial" form of the compression factor polynomial is + * + * Z = 1.0 + C[0]*P + C[1]*P^2 + C[2]*P^3 ... + * + * and these tables do not contain the initial 1.0 term. + * * NOTE! Helium coefficients are a linear mix operation between the * 323K and one for 273K isotherms, to make everything be at 300K. */ -static const double o2_coefficients[4] = { - +1.00117935180448264158, - -0.00074149079841471884, - +0.00000291901111247509, - -0.00000000162092217461 -}; - -static const double n2_coefficients[4] = { - +1.00030344355797817778, - -0.00022528077251905598, - +0.00000295430303276288, - -0.00000000210649996337 -}; - -static const double he_coefficients[4] = { - +1.00000137322788319261, - +0.000488393024887620665, - -0.000000054210166760015, - +0.000000000010908069275 -}; - -static double o2_compressibility_factor(double bar) { return cubic(bar, o2_coefficients); } -static double n2_compressibility_factor(double bar) { return cubic(bar, n2_coefficients); } -static double he_compressibility_factor(double bar) { return cubic(bar, he_coefficients); } - -/* - * We end up using a simple linear mix of the gas-specific functions. - */ double gas_compressibility_factor(struct gasmix *gas, double bar) { - double o2, n2, he; // Z factors - double of, nf, hf; // gas fractions ("partial pressures") - - o2 = o2_compressibility_factor(bar); - n2 = n2_compressibility_factor(bar); - he = he_compressibility_factor(bar); - - of = get_o2(gas) / 1000.0; - hf = get_he(gas) / 1000.0; - nf = 1.0 - of - nf; - - return o2*of + n2*nf + he*hf; + static const double o2_coefficients[3] = { + -0.00071809207370164567, + +0.00000281852572807643, + -0.00000000150290620491 + }; + static const double n2_coefficients[3] = { + -0.00021926035329221337, + +0.00000292844845531647, + -0.00000000207613482075 + }; + static const double he_coefficients[3] = { + +0.00047961098687979363, + -0.00000004077670019935, + +0.00000000000077707035 + }; + double o2, he; + double x1, x2, x3; + double Z; + + o2 = get_o2(gas) / 1000.0; + he = get_he(gas) / 1000.0; + + x1 = bar; x2 = x1*x1; x3 = x2*x1; + + Z = virial_m1(o2_coefficients, x1, x2, x3) * o2 + + virial_m1(he_coefficients, x1, x2, x3) * he + + virial_m1(n2_coefficients, x1, x2, x3) * (1.0 - o2 - he); + + /* + * We add the 1.0 at the very end - the linear mixing of the + * three 1.0 terms is still 1.0 regardless of the gas mix. + */ + return Z + 1.0; } |