Python Forum

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()

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