Changeset 7725842 in subsurface


Ignore:
Timestamp:
Jan 16, 2017, 3:17:40 AM (2 months ago)
Author:
Dirk Hohndel <dirk@…>
Branches:
master
Children:
a8e8d56
Parents:
fedadc65
git-author:
Robert C. Helling <helling@…> (01/12/17 12:19:40)
git-committer:
Dirk Hohndel <dirk@…> (01/16/17 03:17:40)
Message:

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@…>

Location:
core
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • core/dive.h

    r4359974 r7725842  
    137137extern int gas_volume(cylinder_t *cyl, pressure_t p);
    138138extern double gas_compressibility_factor(struct gasmix *gas, double bar);
     139extern double isothermal_pressure(struct gasmix *gas, double p1, int volume1, int volume2);
    139140
    140141
  • core/gas-model.c

    r7be962b r7725842  
    6363        return Z * 0.001 + 1.0;
    6464}
     65
     66/* Compute the new pressure when compressing (expanding) volome v1 at pressure p1 bar to volume v2
     67 * taking into account the compressebility (to first order) */
     68
     69double isothermal_pressure(struct gasmix *gas, double p1, int volume1, int volume2)
     70{
     71        double p_ideal = p1 * volume1 / volume2 / gas_compressibility_factor(gas, p1);
     72
     73        return p_ideal * gas_compressibility_factor(gas, p_ideal);
     74}
  • core/planner.c

    rbb4bf63 r7725842  
    247247                cyl->deco_gas_used.mliter += gas_used.mliter;
    248248        if (cyl->type.size.mliter) {
    249                 delta_p.mbar = gas_used.mliter * 1000.0 / cyl->type.size.mliter;
     249                delta_p.mbar = gas_used.mliter * 1000.0 / cyl->type.size.mliter * gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0);
    250250                cyl->end.mbar -= delta_p.mbar;
    251251        }
     
    831831                deco_volume = get_volume_units(cyl->deco_gas_used.mliter, NULL, &unit);
    832832                if (cyl->type.size.mliter) {
    833                         deco_pressure = get_pressure_units(1000.0 * cyl->deco_gas_used.mliter / cyl->type.size.mliter, &pressure_unit);
    834                         pressure = get_pressure_units(1000.0 * cyl->gas_used.mliter / cyl->type.size.mliter, &pressure_unit);
     833                        int remaining_gas = (double)cyl->end.mbar * cyl->type.size.mliter / 1000.0 / gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0);
     834                        double deco_pressure_bar = isothermal_pressure(&cyl->gasmix, 1.0, remaining_gas + cyl->deco_gas_used.mliter, cyl->type.size.mliter)
     835                                        - cyl->end.mbar / 1000.0;
     836                        deco_pressure = get_pressure_units(1000.0 * deco_pressure_bar, &pressure_unit);
     837                        pressure = get_pressure_units(cyl->start.mbar - cyl->end.mbar, &pressure_unit);
    835838                        /* Warn if the plan uses more gas than is available in a cylinder
    836839                         * This only works if we have working pressure for the cylinder
     
    841844                                        translate("gettextFromC", "this is more gas than available in the specified cylinder!"));
    842845                        else
    843                                 if ((float) cyl->end.mbar * cyl->type.size.mliter / 1000.0 < (float) cyl->deco_gas_used.mliter)
     846                                if ((float) cyl->end.mbar * cyl->type.size.mliter / 1000.0 / gas_compressibility_factor(&cyl->gasmix, cyl->end.mbar / 1000.0)
     847                                    < (float) cyl->deco_gas_used.mliter)
    844848                                        snprintf(warning, sizeof(warning), " &mdash; <span style='color: red;'>%s </span> %s",
    845849                                                translate("gettextFromC", "Warning:"),
Note: See TracChangeset for help on using the changeset viewer.