mirror of
https://github.com/subsurface/subsurface.git
synced 2024-11-28 13:10:19 +00:00
2a966ac2a9
Replace a macro calculating a degree-three polynomial by an inline function. Moreover, calculate the powers 1, 2 and 3 of the pressure inside the function. The compiler will be smart enough to optimize this to the same code. The only important thing is to write "x*x*x*coeff" instead of "coeff*x*x*x". The compiler can't optimize the latter because ... wonderful floating point semantics. Signed-off-by: Berthold Stoeger <bstoeger@mail.tuwien.ac.at>
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 */
|
||
static double virial_m1(const double coeff[], double x)
|
||
{
|
||
return x*coeff[0] + x*x*coeff[1] + x*x*x*coeff[2];
|
||
}
|
||
|
||
/*
|
||
* 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 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);
|
||
|
||
Z = virial_m1(o2_coefficients, bar) * o2 +
|
||
virial_m1(he_coefficients, bar) * he +
|
||
virial_m1(n2_coefficients, bar) * (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;
|
||
}
|