Use real gas compressibility in planner

Modify formluas for gas use to take into account the
compressibility correction for real gases. This introduces
also the inverse formula to compute the pressure for a given
amount of gas.

Signed-off-by: Robert C. Helling <helling@atdotde.de>
This commit is contained in:
Robert C. Helling 2017-01-12 21:19:40 +01:00 committed by Dirk Hohndel
parent fedadc65db
commit 7725842383
3 changed files with 19 additions and 4 deletions

View file

@ -136,6 +136,7 @@ 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);
extern double isothermal_pressure(struct gasmix *gas, double p1, int volume1, int volume2);
static inline int get_o2(const struct gasmix *mix)

View file

@ -62,3 +62,13 @@ double gas_compressibility_factor(struct gasmix *gas, double bar)
*/
return Z * 0.001 + 1.0;
}
/* Compute the new pressure when compressing (expanding) volome v1 at pressure p1 bar to volume v2
* taking into account the compressebility (to first order) */
double isothermal_pressure(struct gasmix *gas, double p1, int volume1, int volume2)
{
double p_ideal = p1 * volume1 / volume2 / gas_compressibility_factor(gas, p1);
return p_ideal * gas_compressibility_factor(gas, p_ideal);
}

View file

@ -246,7 +246,7 @@ static void update_cylinder_pressure(struct dive *d, int old_depth, int new_dept
if (in_deco)
cyl->deco_gas_used.mliter += gas_used.mliter;
if (cyl->type.size.mliter) {
delta_p.mbar = gas_used.mliter * 1000.0 / cyl->type.size.mliter;
delta_p.mbar = gas_used.mliter * 1000.0 / cyl->type.size.mliter * gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0);
cyl->end.mbar -= delta_p.mbar;
}
}
@ -830,8 +830,11 @@ static void add_plan_to_notes(struct diveplan *diveplan, struct dive *dive, bool
volume = get_volume_units(cyl->gas_used.mliter, NULL, &unit);
deco_volume = get_volume_units(cyl->deco_gas_used.mliter, NULL, &unit);
if (cyl->type.size.mliter) {
deco_pressure = get_pressure_units(1000.0 * cyl->deco_gas_used.mliter / cyl->type.size.mliter, &pressure_unit);
pressure = get_pressure_units(1000.0 * cyl->gas_used.mliter / cyl->type.size.mliter, &pressure_unit);
int remaining_gas = (double)cyl->end.mbar * cyl->type.size.mliter / 1000.0 / gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0);
double deco_pressure_bar = isothermal_pressure(&cyl->gasmix, 1.0, remaining_gas + cyl->deco_gas_used.mliter, cyl->type.size.mliter)
- cyl->end.mbar / 1000.0;
deco_pressure = get_pressure_units(1000.0 * deco_pressure_bar, &pressure_unit);
pressure = get_pressure_units(cyl->start.mbar - cyl->end.mbar, &pressure_unit);
/* Warn if the plan uses more gas than is available in a cylinder
* This only works if we have working pressure for the cylinder
* 10bar is a made up number - but it seemed silly to pretend you could breathe cylinder down to 0 */
@ -840,7 +843,8 @@ static void add_plan_to_notes(struct diveplan *diveplan, struct dive *dive, bool
translate("gettextFromC", "Warning:"),
translate("gettextFromC", "this is more gas than available in the specified cylinder!"));
else
if ((float) cyl->end.mbar * cyl->type.size.mliter / 1000.0 < (float) cyl->deco_gas_used.mliter)
if ((float) cyl->end.mbar * cyl->type.size.mliter / 1000.0 / gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0)
< (float) cyl->deco_gas_used.mliter)
snprintf(warning, sizeof(warning), " &mdash; <span style='color: red;'>%s </span> %s",
translate("gettextFromC", "Warning:"),
translate("gettextFromC", "not enough reserve for gas sharing on ascent!"));