Why my graph is blocked at 4.2 bars whithout use the valve on my modelisation ?

[tudy62]

Hello, i want to make a modelisation of hydrogen liquid transformation in gas by change of temperature and pressure, augmentation by change of different velocity. I put a valve who open when pressure max is get but in my graph it is look at 4.2 bar (420 000 Pa), and i don't know how to make the change, can you help me ?

Also maybe my equation of hydrogen are not really good.

my code :

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.interpolate import interp1d


# Load Excel data
file_path = r"C.\Book1.xlsx"
df = pd.read_excel(file_path, sheet_name='Sheet1')

time_data = df["temps (s)"].to_numpy()
speed_data = df["vitesse (km/h)"].to_numpy() * (1000 / 3600)  # Convert to m/s

acceleration_data = np.gradient(speed_data, time_data)

num_points = len(time_data) * 10
time_interp = np.linspace(time_data[0], time_data[-1], num_points)
speed_interp_func = interp1d(time_data, speed_data, kind='cubic')
speed_interp = speed_interp_func(time_interp)
accel_interp_func = interp1d(time_data, acceleration_data, kind='cubic')
accel_interp = accel_interp_func(time_interp)

T0 = 20  # Initial temperature (K)
P0 = 4e5  # Initial pressure (Pa)
mass_h2 = 2.0  # Initial LH2 mass (kg)
dt = np.gradient(time_interp)
R = 8.314  # Universal gas constant

# Molar mass of H2 in g/mol
M_H2 = 2.016

# Number of moles of H2
n_H2 = mass_h2 * 1000 / M_H2

C_H2_L = 9.5  # Specific heat capacity of LH2 (kJ/kg·K)
L_H2 = 446  # Latent heat of vaporization for LH2 (kJ/kg)
R_H2 = 4124  # Specific gas constant for H2 (J/kg·K)
M_H2 = 2.016e-3  # Molar mass in kg/mol
V_tank = 0.05  # Tank volume (m^3)

# TEMP_MAX = 35.0  # Max temperature (K)
PRESS_MAX = 5.5e5  # Max pressure (Pa)
RESET_PRESS = 5e5  # Pressure after relief valve opens (Pa)
MIN_PRESSURE = 4e5  # Minimum pressure (Pa)

temperature = [T0]
pressure = [P0]
mass_LH2 = [mass_h2]
mass_GH2 = [0.01]

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

def update(frame):
    acceleration = accel_interp[frame]
    kinetic_energy = 0.5 * mass_h2 * (speed_interp[frame])**2
    energy_transferred = kinetic_energy
    evaporated_mass = energy_transferred / (L_H2 * 1000)  # Convert to J

    new_mass_LH2 = max(0, mass_LH2[-1] - evaporated_mass)
    new_mass_GH2 = max(0.01, mass_GH2[-1] + evaporated_mass)

    dT = (energy_transferred / (new_mass_LH2 * 1000 * C_H2_L)) if new_mass_LH2 > 0 else 0
    new_temp = temperature[-1] + dT

    n_gas = new_mass_GH2 / M_H2
    new_pressure = max(MIN_PRESSURE, (n_gas * R_H2 * new_temp) / V_tank / 100000)

    if new_pressure > PRESS_MAX:
        new_pressure = RESET_PRESS
        new_temp = (new_pressure * V_tank) / (R * n_H2)
        new_mass_GH2 *= 0.9

    temperature.append(new_temp)
    pressure.append(new_pressure * 1.05)  # Slightly increase pressure each step
    mass_LH2.append(new_mass_LH2)
    mass_GH2.append(new_mass_GH2)

    if len(temperature) > len(time_interp):
        temperature.pop(0)
        pressure.pop(0)
        mass_LH2.pop(0)
        mass_GH2.pop(0)

fig1, ax1 = plt.subplots(figsize=(8, 4))
fig2, ax2 = plt.subplots(figsize=(8, 4))

def animate_temperature(frame):
    update(frame)
    ax1.clear()
    ax1.plot(time_interp[:len(temperature)], temperature, color='red', linewidth=2)
    ax1.set_xlabel('Time (s)')
    ax1.set_ylabel('Temperature (K)', color='red')
    ax1.set_title("Temperature Evolution")

def animate_pressure(frame):
    ax2.clear()
    ax2.plot(time_interp[:len(pressure)], pressure, color='blue', linestyle='dashed', linewidth=2)
    ax2.set_xlabel('Time (s)')
    ax2.set_ylabel('Pressure (Pa)', color='blue')
    ax2.set_title("Pressure Evolution")

ani1 = animation.FuncAnimation(fig1, animate_temperature, frames=len(time_interp), interval=50, repeat=False)
ani2 = animation.FuncAnimation(fig2, animate_pressure, frames=len(time_interp), interval=50, repeat=False)

plt.show()

Larz60+ write Apr-15-2025, 09:08 AM:
Please post all code, output and errors (it it's entirety) between their respective tags. Refer to BBCode help topic on how to post. Use the "Preview Post" button to make sure the code is presented as you expect before hitting the "Post Reply/Thread" button. I have added tags this time, please use BBCode on future posts.

[Axel_Erfurt]

Are you sure you posted the right file Book1.xlsx?

[tudy62]

yes but i translated the name of the elements because english is not my native language, maybe the name in tha code have a little difference