[PATCH] gas model: simplify and improve our Z factor calculations
Linus Torvalds
torvalds at linux-foundation.org
Thu Mar 3 14:07:42 PST 2016
From: Linus Torvalds <torvalds at linux-foundation.org>
Date: Thu, 3 Mar 2016 13:57:54 -0800
Subject: [PATCH] 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 <neolit123 at gmail.com>
Signed-off-by: Linus Torvalds <torvalds at linux-foundation.org>
---
Plotting this with google really makes it all look very good all the way
to 500 bar. And the linear mixing for air is so good and matches the air
polynomial we used so well that there is no longer any point in having a
separate function for air even if we have the tables for it.
And all credit for the heavy lifting really obviously goes to Lubomir, I
just wrote the trivial code.
subsurface-core/gas-model.c | 149 +++++++++-----------------------------------
1 file changed, 30 insertions(+), 119 deletions(-)
diff --git a/subsurface-core/gas-model.c b/subsurface-core/gas-model.c
index b90467bbd7c0..420233612cb9 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;
--
2.7.2.334.g7c0da37
More information about the subsurface
mailing list