Fueling protocol for h2 into a Type IV cylinder.
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 10 10:51:14 2023
@author: THOMAS-SCBoduch
"""
import math
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
f=open("wenger_init.txt","r")
lines=f.readlines()
for line in lines:
if "dispenser_temperature" in line:
dispenser_temperature=float(line.split("=")[1])
if "dispenser_pressure" in line:
dispenser_pressure_init=float(line.split("=")[1])
if "ambiant_temperature" in line:
ambiant_temperature=float(line.split("=")[1])
if "Tank_volume" in line:
Tank_Volume=float(line.split("=")[1])
if "target_pressure" in line:
target_pressure=float(line.split("=")[1])
if "APRR_target" in line:
APRR_target=float(line.split("=")[1])
if "DPhigh" in line:
DPhigh=float(line.split("=")[1])
if "DPlow" in line:
DPlow=float(line.split("=")[1])
lenght=1.83
radius=0.275
max_APRR_allowed = 10.0 # bars/s
APRR_low=0.3
Base_Volume = 0.833 #m^3
P_Max_Border = dispenser_pressure_init + 40.0 #barg
dt = 1e-1 #s
fueling_time=0
def calculate_pressure_ramp_rate(APRR_target,max_APRR_allowed,Base_Volume,Tank_Volume):
return min(APRR_target* 10.0 / 60.0,max_APRR_allowed * Base_Volume / Tank_Volume * 10.0 / 60.0); # 10.0 / 60.0 MPa/min to Bar/s
def soave_redlich_kwong(pressure, volume, temperature):
#Calculates the Soave-Redlich-Kwong equation to determine pressure based on temperature.
#param gas_constant: Gas constant J/(mol·K).
gas_constant=8.314
#param acentric_factor: Acentric factor no dimension
acentric_factor=0.234
#param critical_temperature: Critical temperature K
critical_temperature=33.19
#param critical_pressure: Critical pressure bar
critical_pressure=12.9
#return: Calculated pressure
Tr = temperature / critical_temperature
alpha = (1 + (0.48 + 1.574*acentric_factor - 0.176*acentric_factor**2) * (1 - np.sqrt(Tr)))**2
a = 0.42748 * (gas_constant**2) * (critical_temperature**2) / critical_pressure
b = 0.08664 * gas_constant * critical_temperature / critical_pressure
return gas_constant * temperature / (volume - b) - alpha * a / (volume**2 + 2 * b * volume - b**2)
def heat_transfer(m_io, h_io, m, U, temperature, density, enthalpy, volume_tank, time_step):
"""
Calculates heat transfer in the storage using the heat conduction equation.
:param m_io: Mass flow rate
:param h_io: Enthalpy
:param m: Gas mass
:param U: Internal energy
:param temperature: Temperature
:param density: Density
:param enthalpy: Enthalpy
:param volume_tank: Tank volume
:param time_step: Time step for calculation
:return: Heat transfer
"""
surface_area = 2 * volume_tank / radius # Assuming the tank has a cylinder shape , diameter= 55 cm so radius= 55/2 = 27,5 cm = 0,275 m
thermal_conductivity = 0.2 # Thermal conductivity for hydrogen
temperature_change = 10 # Assuming temperature change in the storage
# Heat conduction equation: Q = (m_io * h_io - m * enthalpy) + U * density * volume_tank * temperature_change
q_convection = (m_io * h_io - m * enthalpy) + U * density * volume_tank * temperature_change
# Heat transfer by conduction
q_conduction = (thermal_conductivity * surface_area * temperature_change) / volume_tank
# Total heat transfer
q_total = q_convection + q_conduction
return q_total * time_step
def tank_diffusion():
max_iter_time = 50
num_r=50
num_theta=50
alpha = 1.1
delta_r = radius / num_r
delta_theta = 2*np.pi / num_theta
delta_t = (delta_r ** 2)/(4 * alpha)
gamma = (alpha * delta_t) / (delta_r ** 2)
r=np.linspace(0, radius, num_r)
theta=np.linspace(0, 2*np.pi, num_r)
z=np.linspace(0, lenght, num_r)
T=np.empty((max_iter_time, num_r, num_theta))
R,Theta = np.meshgrid(r,theta)
#Boundarie conditions#
T.fill(0)
T[:,num_r-1,] = ambiant_temperature
def FCTS_polar(T):
for t in range(0,max_iter_time-1):
for i in range(1,num_r-1):
r_i = i* delta_r
T[t+1,i,:]= T[t,i,:] + gamma* ( 1 / r_i) * T[t,i+1,:]
return T
def refueling(APRR,Plimit_high,Plimit_low,Pressure,fueling_time,dt):
Pressure_new= APRR*fueling_time+dispenser_pressure_init
Plimit_high_new = min(Plimit_high + APRR * dt/60, Pressure + DPhigh+ DPlow)
Plimit_high_new = Plimit_high + max_APRR_allowed * dt/60
Plimit_low_new= dispenser_pressure_init + dt * APRR_low/60.0 - DPlow
fueling_time=fueling_time+dt
return Pressure_new,fueling_time,Plimit_high_new,Plimit_low_new
#
# Values Selection
# APRR are expected to have same values in com/no com tables -> all no com
# Implement the Wenger protocol rules
# INITIALISATION
#INITIALISATION CYCLE
APRR=calculate_pressure_ramp_rate(APRR_target,max_APRR_allowed,Base_Volume,Tank_Volume)
Plimit_low= dispenser_pressure_init + dt * APRR_low/60.0 - DPlow
Plimit_high= dispenser_pressure_init + DPhigh
print("APRR is",APRR,"barg/s")
sim=refueling(APRR , Plimit_high , Plimit_low, dispenser_pressure_init , 0, dt)
#PLOT PREPARATION
pressure_plot=[sim[0]]
t=[sim[1]]
corridor_down=[sim[3]]
corridor_up=[sim[2]]
#LOOP FOR CALCULATION
while sim[0]
sim=refueling( APRR , sim[2], sim[3], sim[0], sim[1], dt)
pressure_plot.append(sim[0])
t.append(sim[1])
corridor_down.append(sim[3])
corridor_up.append(sim[2])
print("duration of the refueling",t[-1]/60,"min")
print("last corridor down",corridor_down[-1],"MPa")
plt.plot(t,pressure_plot)
plt.plot(t,corridor_down)
plt.plot(t,corridor_up)
tank_diffusion()
Plot essai