From a173a3ce79de6b53d7871cba8acfa1f07565d829 Mon Sep 17 00:00:00 2001 From: Linus Torvalds Date: Thu, 3 Mar 2016 13:57:54 -0800 Subject: gas model: simplify and improve our Z factor calculations Lubomir found better compressibility data for the pure gases that we need for scuba, making the air table superfluous: we get good values from just regular linear mixing of the Oxygen, Nitrogen and Helium calculations. Also, rather than using a quintic polynomial, a cubic one does sufficiently well, making for smaller code and fewer coefficients. And judging by the reactions from people on G+ (as well as just looking at how good the fit is with the air data), this is all the right way to do this, and this thus removes the Redlich-Kwong equation. All-credit-goes-to: Lubomir I. Ivanov Signed-off-by: Linus Torvalds Signed-off-by: Dirk Hohndel --- subsurface-core/gas-model.c | 149 +++++++++----------------------------------- 1 file changed, 30 insertions(+), 119 deletions(-) (limited to 'subsurface-core/gas-model.c') diff --git a/subsurface-core/gas-model.c b/subsurface-core/gas-model.c index b90467bbd..420233612 100644 --- a/subsurface-core/gas-model.c +++ b/subsurface-core/gas-model.c @@ -5,159 +5,70 @@ #include "dive.h" /* - * This gives an interative solution of hte Redlich-Kwong equation for the compressibility factor - * according to https://en.wikipedia.org/wiki/Redlich–Kwong_equation_of_state - * in terms of the reduced temperature T/T_crit and pressure p/p_crit. - * - * Iterate this three times for good results in our pressur range. - * + * Generic cubic polynomial */ - -static double redlich_kwong_equation(double t_red, double p_red, double z_init) -{ - return (1.0/(1.0 - 0.08664*p_red/(t_red * z_init)) - - 0.42748/(sqrt(t_red * t_red * t_red) * ((t_red*z_init/p_red + 0.08664)))); -} - -/* - * At high pressures air becomes less compressible, and - * does not follow the ideal gas law any more. - */ -#define STANDARD_TEMPERATURE 293.0 - -static double redlich_kwong_compressibility_factor(struct gasmix *gas, double bar) -{ - /* Critical points according to https://en.wikipedia.org/wiki/Critical_point_(thermodynamics) */ - - double tcn2 = 126.2; - double tco2 = 154.6; - double tche = 5.19; - - double pcn2 = 33.9; - double pco2 = 50.5; - double pche = 2.27; - - double tc, pc; - - tc = (tco2 * get_o2(gas) + tche * get_he(gas) + tcn2 * (1000 - get_o2(gas) - get_he(gas))) / 1000.0; - pc = (pco2 * get_o2(gas) + pche * get_he(gas) + pcn2 * (1000 - get_o2(gas) - get_he(gas))) / 1000.0; - - return (redlich_kwong_equation(STANDARD_TEMPERATURE/tc, bar/pc, - redlich_kwong_equation(STANDARD_TEMPERATURE/tc, bar/pc, - redlich_kwong_equation(STANDARD_TEMPERATURE/tc, bar/pc,1.0)))); -} - -/* - * Generic quintic polynomial - */ -static double quintic(double bar, const double coefficient[]) +static double cubic(double bar, const double coefficient[]) { double x0 = 1.0, x1 = bar, x2 = x1*x1, - x3 = x2*x1, - x4 = x2*x2, - x5 = x2*x3; + x3 = x2*x1; return x0 * coefficient[0] + x1 * coefficient[1] + x2 * coefficient[2] + - x3 * coefficient[3] + - x4 * coefficient[4] + - x5 * coefficient[5]; + x3 * coefficient[3]; } /* - * These are the quintic coefficients by Lubomir I. Ivanov that have - * been optimized for the least-square error to the air - * compressibility factor table (at 300K) taken from Wikipedia: + * Cubic least-square coefficients for O2/N2/He based on data from * - * bar z_factor - * --- ------ - * 1: 0.9999 - * 5: 0.9987 - * 10: 0.9974 - * 20: 0.9950 - * 40: 0.9917 - * 60: 0.9901 - * 80: 0.9903 - * 100: 0.9930 - * 150: 1.0074 - * 200: 1.0326 - * 250: 1.0669 - * 300: 1.1089 - * 400: 1.2073 - * 500: 1.3163 - */ -static const double air_coefficients[6] = { - +1.0002556612420115, - -0.0003115084635183305, - +0.00000227808965401253, - +1.91596422989e-9, - -8.78421542e-12, - +6.77746e-15 -}; - -/* - * Quintic least-square coefficients for O2/N2/He based on tables at + * PERRY’S CHEMICAL ENGINEERS’ HANDBOOK SEVENTH EDITION * - * http://ww.baue.org/library/zfactor_table.php + * with the lookup and curve fitting by Lubomir. * - * converted to bar and also done by Lubomir. + * 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[6] = { - +1.0002231211532653, - -0.0007471497056767194, - +0.00000200444854807816, - +2.91501995188e-9, - -4.48294663e-12, - -6.11529e-15 +static const double o2_coefficients[4] = { + +1.00117935180448264158, + -0.00074149079841471884, + +0.00000291901111247509, + -0.00000000162092217461 }; -static const double n2_coefficients[6] = { - +1.0001898816185364, - -0.00030793319362077315, - +0.00000327557417347714, - -1.93872574476e-9, - -2.7732353e-12, - -2.8921e-16 +static const double n2_coefficients[4] = { + +1.00030344355797817778, + -0.00022528077251905598, + +0.00000295430303276288, + -0.00000000210649996337 }; -static const double he_coefficients[6] = { - +0.9998700261301693, - +0.0005452130351730479, - -2.7853712233619e-7, - +5.5935404211e-10, - -1.35114572e-12, - +2.00476e-15 +static const double he_coefficients[4] = { + +1.00000137322788319261, + +0.000488393024887620665, + -0.000000054210166760015, + +0.000000000010908069275 }; -static double air_compressibility_factor(double bar) { return quintic(bar, air_coefficients); } -static double o2_compressibility_factor(double bar) { return quintic(bar, o2_coefficients); } -static double n2_compressibility_factor(double bar) { return quintic(bar, n2_coefficients); } -static double he_compressibility_factor(double bar) { return quintic(bar, he_coefficients); } +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 specialized functions for known gases, because - * we have special tables for them. - * - * For air we use our known-good table. For other mixes we use a - * linear interpolation of the Z factors of the individual gases. + * 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") - if (gasmix_is_air(gas)) - return air_compressibility_factor(bar); - o2 = o2_compressibility_factor(bar); n2 = n2_compressibility_factor(bar); he = he_compressibility_factor(bar); - of = gas->o2.permille / 1000.0; - hf = gas->he.permille / 1000.0; + of = get_o2(gas) / 1000.0; + hf = get_he(gas) / 1000.0; nf = 1.0 - of - nf; return o2*of + n2*nf + he*hf; -- cgit v1.2.3-70-g09d2