Hi,
After spending some hours with this and mathematica, I realized that this is in fact a waste of time since it gets the difference of the factor from 1 significantly wrong in the parameter range we are interested in. But the Redlich-Kwong-equation (of whose existence I learned from wikipedia) does a much better job numerically. So here is a patch that implements it. This equation is a cubic equation for Z, so I decided to solve it iteratively (three iterations seem to be good enough). |
From 602aa415d7662b21f4c7405e1c9e7d32868ce36f Mon Sep 17 00:00:00 2001 From: "Robert C. Helling" <[email protected]> Date: Thu, 25 Feb 2016 23:20:23 +0100 Subject: [PATCH] Compute the air compressibility factor To: [email protected]
This computes the compressibility of the gas depending on the mixture. As it turns out, the van der Waals equation gives results that are numerically not really useful, so we use the Redlich Kwong equation which is, according to Wikipedia, much more accurate, which can be confirmed given the empirical values for air. As opposed to the previous approach with a look-up table, this takes into account the actual gasmix and temperature (if known), but for air at 300K gives slightly less accurate results. Signed-off-by: Robert C. Helling <[email protected]> --- qt-models/cylindermodel.cpp | 3 +- subsurface-core/dive.c | 77 ++++++++++++++++----------------------- subsurface-core/dive.h | 6 +-- subsurface-core/divelist.c | 2 +- subsurface-core/libdivecomputer.c | 3 +- subsurface-core/profile.c | 4 +- subsurface-core/statistics.c | 2 +- 7 files changed, 43 insertions(+), 54 deletions(-) diff --git a/qt-models/cylindermodel.cpp b/qt-models/cylindermodel.cpp index 1a587ce..0e46465 100644 --- a/qt-models/cylindermodel.cpp +++ b/qt-models/cylindermodel.cpp @@ -36,7 +36,8 @@ static QString get_cylinder_string(cylinder_t *cyl) // liters: if we don't have a working pressure, we cannot // convert the cylinder size to cuft. if (wp.mbar && prefs.units.volume == units::CUFT) { - double real_value = ml_to_cuft(gas_volume(cyl, wp)); + temperature_t faketemp = {.mkelvin = 0}; + double real_value = ml_to_cuft(gas_volume(cyl, wp, faketemp)); value = ml_to_cuft(ml) * bar_to_atm(wp.mbar / 1000.0); decimals = (value > 20.0) ? 0 : (value > 2.0) ? 1 : 2; unit = QString("(%1)%2").arg(real_value, 0, 'f', 0).arg(CylindersModel::tr("cuft")); diff --git a/subsurface-core/dive.c b/subsurface-core/dive.c index 8c69e9a..fafb1db 100644 --- a/subsurface-core/dive.c +++ b/subsurface-core/dive.c @@ -847,63 +847,50 @@ static void update_min_max_temperatures(struct dive *dive, temperature_t tempera } /* + * 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. + * + */ + +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. - * - * NOTE! The compressibility doesn't just depend on the - * gas, but on temperature too. However, this currently - * just follows the 300K curve for air, and ignores the - * gasmix. And the temperature we don't really even have - * a way to try to take into account. */ -#define ARRAY_SIZE(array) (sizeof(array)/sizeof(array[0])) -double gas_compressibility_factor(struct gasmix *gas, double bar) -{ - static const struct z_factor { - double bar, z_factor; - } air_table[] = { - { 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 } - }; - const struct z_factor *n; - int i; +double gas_compressibility_factor(struct gasmix *gas, double bar, double kelvin) +{ + /* Critical points according to https://en.wikipedia.org/wiki/Critical_point_(thermodynamics) */ - for (i = 0; i < ARRAY_SIZE(air_table); i++) { - const struct z_factor *n = air_table+i; - double frac; + double tcn2 = 126.2; + double tco2 = 154.6; + double tche = 5.19; - if (n->bar < bar) - continue; - if (!i) - return n->z_factor; + double pcn2 = 33.9; + double pco2 = 50.5; + double pche = 2.27; - /* How far from this one? */ - frac = (n->bar - bar) / (n->bar - n[-1].bar); + double tc, pc; - /* Silly linear interpolation */ - return frac*n[-1].z_factor + (1.0-frac)*n->z_factor; - } + 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; - /* Over 500 bar? We make shit up */ - return air_table[ARRAY_SIZE(air_table)-1].z_factor; + return (redlich_kwong_equation(kelvin/tc, bar/pc, + redlich_kwong_equation(kelvin/tc, bar/pc, + redlich_kwong_equation(kelvin/tc, bar/pc,1.0)))); } -int gas_volume(cylinder_t *cyl, pressure_t p) +int gas_volume(cylinder_t *cyl, pressure_t p, temperature_t airtemp) { double bar = p.mbar / 1000.0; - double z_factor = gas_compressibility_factor(&cyl->gasmix, bar); + double z_factor = gas_compressibility_factor(&cyl->gasmix, bar, airtemp.mkelvin ? airtemp.mkelvin / 1000.0 : 300.0); return cyl->type.size.mliter * bar_to_atm(bar) / z_factor; } diff --git a/subsurface-core/dive.h b/subsurface-core/dive.h index 0f5c434..3ad2108 100644 --- a/subsurface-core/dive.h +++ b/subsurface-core/dive.h @@ -127,9 +127,9 @@ extern double get_vertical_speed_units(unsigned int mms, int *frac, const char * extern unsigned int units_to_depth(double depth); extern int units_to_sac(double volume); -/* Volume in mliter of a cylinder at pressure 'p' */ -extern int gas_volume(cylinder_t *cyl, pressure_t p); -extern double gas_compressibility_factor(struct gasmix *gas, double bar); +/* Volume in mliter of a cylinder at pressure 'p' at temperature kelvin*/ +extern int gas_volume(cylinder_t *cyl, pressure_t p, temperature_t airtemp); +extern double gas_compressibility_factor(struct gasmix *gas, double bar, double kelvin); static inline int get_o2(const struct gasmix *mix) diff --git a/subsurface-core/divelist.c b/subsurface-core/divelist.c index 4008678..6546a05 100644 --- a/subsurface-core/divelist.c +++ b/subsurface-core/divelist.c @@ -274,7 +274,7 @@ static double calculate_airuse(struct dive *dive) if (!end.mbar || start.mbar <= end.mbar) continue; - airuse += gas_volume(cyl, start) - gas_volume(cyl, end); + airuse += gas_volume(cyl, start, dive->airtemp) - gas_volume(cyl, end, dive->airtemp); } return airuse / 1000.0; } diff --git a/subsurface-core/libdivecomputer.c b/subsurface-core/libdivecomputer.c index 3990884..337a20c 100644 --- a/subsurface-core/libdivecomputer.c +++ b/subsurface-core/libdivecomputer.c @@ -136,7 +136,8 @@ static int parse_gasmixes(device_data_t *devdata, struct dive *dive, dc_parser_t dive->cylinder[i].type.workingpressure.mbar *= 206.843 / 206.7; char name_buffer[9]; int rounded_size = ml_to_cuft(gas_volume(&dive->cylinder[i], - dive->cylinder[i].type.workingpressure)); + dive->cylinder[i].type.workingpressure, + dive->airtemp)); rounded_size = (int)((rounded_size + 5) / 10) * 10; switch (dive->cylinder[i].type.workingpressure.mbar) { case 206843: diff --git a/subsurface-core/profile.c b/subsurface-core/profile.c index 8d5cebd..be3e2e8 100644 --- a/subsurface-core/profile.c +++ b/subsurface-core/profile.c @@ -190,7 +190,7 @@ static int get_local_sac(struct plot_data *entry1, struct plot_data *entry2, str cyl = dive->cylinder + index; - airuse = gas_volume(cyl, a) - gas_volume(cyl, b); + airuse = gas_volume(cyl, a, dive->airtemp) - gas_volume(cyl, b, dive->airtemp); /* milliliters per minute */ return airuse / atm * 60 / duration; @@ -736,7 +736,7 @@ static int sac_between(struct dive *dive, struct plot_data *first, struct plot_d a.mbar = GET_PRESSURE(first); b.mbar = GET_PRESSURE(last); cyl = dive->cylinder + first->cylinderindex; - airuse = gas_volume(cyl, a) - gas_volume(cyl, b); + airuse = gas_volume(cyl, a, dive->airtemp) - gas_volume(cyl, b, dive->airtemp); if (airuse <= 0) return 0; diff --git a/subsurface-core/statistics.c b/subsurface-core/statistics.c index 373a6a0..b3855eb 100644 --- a/subsurface-core/statistics.c +++ b/subsurface-core/statistics.c @@ -363,7 +363,7 @@ void get_gas_used(struct dive *dive, volume_t gases[MAX_CYLINDERS]) start = cyl->start.mbar ? cyl->start : cyl->sample_start; end = cyl->end.mbar ? cyl->end : cyl->sample_end; if (end.mbar && start.mbar > end.mbar) - gases[idx].mliter = gas_volume(cyl, start) - gas_volume(cyl, end); + gases[idx].mliter = gas_volume(cyl, start, dive->airtemp) - gas_volume(cyl, end, dive->airtemp); } } -- 2.5.4 (Apple Git-61)
From what I learned, this is all black magic when it comes too mixtures of gases (like air or trimix) but the standard approach is to go via critical values and for those compute weighted averages. There is a whole industry for these real gas calculations as those are quite important for the oil/gas industry. Best Robert |
signature.asc
Description: Message signed with OpenPGP using GPGMail
_______________________________________________ subsurface mailing list [email protected] http://lists.subsurface-divelog.org/cgi-bin/mailman/listinfo/subsurface
