CompressedGasTankModel/Paintball_thermo.py at main · petermaginot/CompressedGasTankModel · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
# cycle_calc.py
# Python 3.14.0 translation of CycleCalc.bas
# Requires: pip install CoolProp pandas

import pandas as pd
from CoolProp.CoolProp import PropsSI
import math

def expansion_calc(target_pressure, U01, n01, V1, U02, n02, V2, fluidname, imposedPhase):
    """
    Perform stepwise gas transfer between two systems until target pressure is reached.
    All units are SI (Pa, J, mol, m³).
    Returns: n1, n2, MolarU1, MolarU2
    """
    step_size = -0.001  # moles transferred from system 1 to system 2 per step
    max_steps = 1000
    moles_transferred = 0.0

    molar_U1, molar_U2 = U01, U02
    n1, n2 = n01, n02
    step_pressure = 0.0
    prev_pressure = 0.0

    for _ in range(max_steps):
        U1i = molar_U1 * n1
        U2i = molar_U2 * n2

        molar_H1 = PropsSI("HMOLAR", "UMOLAR" + imposedPhase, molar_U1, "DMOLAR", n1 / V1, fluidname)
        n1 += step_size
        n2 -= step_size

        molar_U1 = (U1i + molar_H1 * step_size) / n1
        molar_U2 = (U2i - molar_H1 * step_size) / n2
        moles_transferred += step_size

        prev_pressure = step_pressure
        step_pressure = PropsSI("P", "UMOLAR" + imposedPhase, molar_U2, "DMOLAR", n2 / V2, fluidname)

        # stop when target pressure exceeded
        if step_pressure > target_pressure:
            overshoot = -step_size * (step_pressure - target_pressure) / (step_pressure - prev_pressure)
            n1 += overshoot
            n2 -= overshoot
            molar_U1 = (U1i + molar_H1 * (step_size + overshoot)) / n1
            molar_U2 = (U2i - molar_H1 * (step_size + overshoot)) / n2
            moles_transferred += overshoot
            break

    return n1, n2, molar_U1, molar_U2


def test():
    """
    molarH1 = molarU1 + P * molarV1
    molarV1 = vessel_volume/n1
    Mass balance
    dn1/dt = -ndot
    Energy balance:
    dU1/dt = -ndot*molarH1=molarH1*dn1/dt
    dU1/dn1 = Hhat1

    U1f = (U1o+Hhat1*deltan)/(n1o+deltan)
    """

    fluidname = "HEOS::Nitrogen"

    TankTemp = 300.0                       # K
    TankP = 4500 * 101325 / 14.696         # Pa
    VCTemp = 300.0                         # K
    TargetVCPressure = 115 * 101325 / 14.696
    AtmoPressure = 101325.0
    TankMass = 1.0                         # kg
    TankCP = 900                           # J/kg·K
    VCVol = 2.01628984537637E-05           # m³
    TankVol = 0.001261804                  # m³
    imposedPhase = ''                      # change to '|gas' to impose the gas phase, which may speed up calculations


    # Initial molar properties
    MolarU1 = PropsSI("UMOLAR", "T" + imposedPhase, TankTemp, "P", TankP, fluidname)
    MolarU2 = PropsSI("UMOLAR", "T" + imposedPhase, VCTemp, "P", AtmoPressure, fluidname)
    n1 = TankVol * PropsSI("DMOLAR", "T" + imposedPhase, TankTemp, "P", TankP, fluidname)
    n2 = VCVol * PropsSI("DMOLAR", "T" + imposedPhase, VCTemp, "P", AtmoPressure, fluidname)

    data = []

    for i in range(1, 4001):
        if i % 10 == 0:
            print(i)
        n1, n2, MolarU1, MolarU2 = expansion_calc(
            TargetVCPressure, MolarU1, n1, TankVol, MolarU2, n2, VCVol, fluidname, imposedPhase
        )

        GasTemp = PropsSI("T", "UMOLAR" + imposedPhase, MolarU1, "DMOLAR", n1 / TankVol, fluidname)
        GasCv = PropsSI("CVMOLAR", "UMOLAR" + imposedPhase, MolarU1, "DMOLAR", n1 / TankVol, fluidname)

        TankEquilibriumTemp = (n1 * GasTemp * GasCv + TankTemp * TankMass * TankCP) / (
            n1 * GasCv + TankMass * TankCP
        )

        TankTemp = TankEquilibriumTemp
        MolarU1 = PropsSI("UMOLAR", "T" + imposedPhase, TankTemp, "DMOLAR", n1 / TankVol, fluidname)

        # Isentropic expansion of chamber gas
        VCTempMax = PropsSI("T", "UMOLAR" + imposedPhase, MolarU2, "DMOLAR", n2 / VCVol, fluidname)
        S2 = PropsSI("SMOLAR", "UMOLAR" + imposedPhase, MolarU2, "DMOLAR", n2 / VCVol, fluidname)
        MolarU2 = PropsSI("UMOLAR", "SMOLAR" + imposedPhase, S2, "P", AtmoPressure, fluidname)
        n2 = PropsSI("DMOLAR", "SMOLAR" + imposedPhase, S2, "P", AtmoPressure, fluidname) * VCVol

        TankP = PropsSI("P", "T" + imposedPhase, TankTemp, "DMOLAR", n1 / TankVol, fluidname)
        VCTemp = PropsSI("T", "SMOLAR" + imposedPhase, S2, "P", AtmoPressure, fluidname)

        data.append([i, n1, n2, TankTemp, TankP, VCTempMax, VCTemp])

        if TankP < TargetVCPressure:
            print(f"Stopped after {i} shots — tank pressure below target.")
            break

    # Save to CSV
    df = pd.DataFrame(data, columns=[
        "Shot", "n1", "n2", "TankTemp[K]", "TankP[Pa]", "VCTempMax[K]", "VCTemp[K]"
    ])
    df.to_csv("ModelOutput.csv", index=False)
    print("Results saved to ModelOutput.csv")


if __name__ == "__main__":
    test()