mirror of
https://github.com/subsurface/subsurface.git
synced 2025-02-19 22:16:15 +00:00
a8b1c86139
The curve fitting for our gas compressibility was only done in the sane range of 0-500 bar, which is what a scuba cylinder can reasonably be expected to perhaps have. But the planner ends up happily using negative cylinder pressures when you run out of gas, and then the compressibility gives nonsensical results. That's clearly a planner bug, but the nonsensical gas compressibility values made it harder to see what could be wrong. So we just clamp the inpot range to the range we have verified against experimental data. If you try to get compressibility for negative pressures, you get the compressibility for an ideal and imaginary gas. And if you try to get compressibility for pressures over 500 bar, we'll just assume that it's 500 bar. Signed-off-by: Linus Torvalds <torvalds@linux-foundation.org>
92 lines
2.6 KiB
C
92 lines
2.6 KiB
C
// SPDX-License-Identifier: GPL-2.0
|
||
/* gas-model.c */
|
||
/* gas compressibility model */
|
||
#include <stdio.h>
|
||
#include <stdlib.h>
|
||
#include "dive.h"
|
||
|
||
/* "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)
|
||
|
||
/*
|
||
* Z = pV/nRT
|
||
*
|
||
* 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.
|
||
*/
|
||
double gas_compressibility_factor(struct gasmix gas, double bar)
|
||
{
|
||
static const double o2_coefficients[3] = {
|
||
-7.18092073703e-04,
|
||
+2.81852572808e-06,
|
||
-1.50290620492e-09
|
||
};
|
||
static const double n2_coefficients[3] = {
|
||
-2.19260353292e-04,
|
||
+2.92844845532e-06,
|
||
-2.07613482075e-09
|
||
};
|
||
static const double he_coefficients[3] = {
|
||
+4.87320026468e-04,
|
||
-8.83632921053e-08,
|
||
+5.33304543646e-11
|
||
};
|
||
int o2, he;
|
||
double x1, x2, x3;
|
||
double Z;
|
||
|
||
/*
|
||
* The curve fitting range is only [0,500] bar.
|
||
* Anything else is way out of range for cylinder
|
||
* pressures.
|
||
*/
|
||
if (bar < 0) bar = 0;
|
||
if (bar > 500) bar = 500;
|
||
|
||
o2 = get_o2(gas);
|
||
he = get_he(gas);
|
||
|
||
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) * (1000 - 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.
|
||
*
|
||
* The * 0.001 is because we did the linear mixing using the
|
||
* raw permille gas values.
|
||
*/
|
||
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);
|
||
}
|
||
|
||
double gas_density(struct gasmix gas, int pressure)
|
||
{
|
||
int density = gas.he.permille * HE_DENSITY + gas.o2.permille * O2_DENSITY + (1000 - gas.he.permille - gas.o2.permille) * N2_DENSITY;
|
||
|
||
return density * (double) pressure / gas_compressibility_factor(gas, pressure / 1000.0) / SURFACE_PRESSURE / 1000000.0;
|
||
}
|