Physical model of non-adiabatic blowdown of cryo-compressed hydrogen storage tanks

(cid:1) Physical model of blowdown dynamics for cryo-compressed hydrogen storage tanks. (cid:1) The model accounts for heat transfer at both the tank and the discharge pipe walls. (cid:1) The model employs the high-accuracy Helmholtz energy formulations. (cid:1) Validation against tests with pressure 0.6 e 20 MPa and temperature 80 e 310 K.

accurately reproducing the temperature distribution in time and space for the cryogenic hydrogen jet experiment.© 2023 The Authors.Published by Elsevier Ltd on behalf of Hydrogen Energy Publications LLC.This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/).

Introduction
The fast-growing market of hydrogen technologies requires competitive techniques to store and transport large quantities of this energy carrier.The cryo-compressed hydrogen (CcH2) storage is being investigated as capable to optimise the gravimetric and volumetric capacities against the energy required for the compression and cooling down of the hydrogen gas in comparison to commercially used compressed gaseous (CGH2) and liquid hydrogen (LH 2 ) respectively [1,2].Furthermore, storage of cryo-compressed hydrogen is not affected by the boil-off as LH 2 storage.The studies [3e5] investigated the application of storage systems refuelled with hydrogen at 80 K and pressures up to 35 MPa for light duty vehicles.For example, Lawrence Livermore National Laboratory assessed the performances of generation 2 prototype of cryogenic pressure vessel at about 34.5 MPa installed on a Toyota Prius [4].The CcH2 stationary storage solution at pressure of 30 MPa is being evaluated for the refuelling of heavy duty vehicles by "H2 Mobility", operating more than 90 refuelling stations [6].
The inherently safer design of CcH2 storage systems and refuelling infrastructure requires an understanding of potential incident consequences.In case of a release through the Thermally Activated Pressure Relief Device (TPRD) or other relief device installed on a storage system, the hydrogen blowdown dynamics, including transient mass transfer, will be affected by the heat transfer in the system.This, in turn, would influence hydrogen parameters at the release nozzle, and, consequently, the hazard distances of unignited and ignited jets.During blowdown of pressurised hydrogen system, temperature in a storage tank decreases due to the gas expansion.This process competes with the tendency of the gas temperature to increase due to the heat transfer through the tank wall from the surrounding atmosphere to hydrogen.In 2007 Schefer et al. [7] highlighted the importance of heat transfer in the storage tank blowdown dynamics.This effect should be even more pronounced in the case of CcH2 storage tanks with a damaged insulation.
Several experimental and analytical studies have assessed the effect of heat transfer during filling procedure of highpressure hydrogen storage tanks at initial ambient temperature.In 2007, Monde et al. [8] proposed a theoretical model to assess the increase of temperature during fuelling of a hydrogen storage tank to pressure of 35 MPa.Their model is based on the coupling of energy balance with unsteady onedimensional (1D) heat conduction at the tank wall, with the assumption of a constant internal heat transfer coefficient during the fuelling procedure.Then, this theoretical model was used in 2008 by Woodfield et al. [9] to simulate experiments on filling of hydrogen storage tanks with volume in the range 39e205 L up to pressure of 35 and 70 MPa with initial ambient temperature.Comparison against experiments showed an overprediction by the theoretical model of the measured gas temperature.In 2012, Monde et al. [10] applied the theoretical model [8] to simulate the temperature rise during filling of hydrogen storage tanks with volume 31e40 L and pressures up to 35 and 70 MPa.The comparison between calculations and experiments highlighted the drawbacks in assuming a constant internal heat transfer coefficient, as this had to be varied among tests in the range 35e200 W/m 2 /K to reproduce experiments.The choice of a constant heat transfer coefficient may prevent applicability of the model to arbitrary initial conditions of the filling or blowdown procedure, e.g. at cryogenic temperatures, but also hinder the correct representation of the fast-changing conditions inside the tank.In 2019, Molkov et al. [11] developed a physical model to reproduce the thermal behaviour of onboard hydrogen tanks during fuelling.The model accounts for the convective heat transfer between hydrogen, tank wall and the atmosphere using Nusselt number correlations and original hypothesis of entrainment into the jet.The non-ideal behaviour of hydrogen gas under high-pressure was taken into account by using the Abel-Noble Equation of State (EoS).The physical model was validated against experiments on fuelling of hydrogen storage tanks with volume in the range 29e74 L and pressures up to 77 MPa.The model accuracy was within temperature deviations measured during experiments.
Fewer studies have been conducted on the modelling of non-adiabatic blowdown of hydrogen storage tanks.In 2021, Molkov et al. [12] developed a physical model accounting for heat transfer through the wall of high pressure hydrogen tanks in an engulfing fire while releasing hydrogen through a TPRD.This model employed the under-expanded jet theory developed earlier by Molkov et al. [13] to calculate parameters at the real and notional nozzle exits.The model was validated against experimental data on the blowdown of 19 L, 70 MPa Type IV tank with 1 mm TPRD orifice filled in by helium, and the destructive fire test with 36 L, 70 MPa Type IV hydrogen tank.While the Abel-Noble EoS was proved to represent well conditions and blowdown dynamics of high-pressure hydrogen storage tanks at initial ambient temperature [12], it may have limited applicability to cryogenic hydrogen gas.In this case, the high-accuracy Helmholtz energy formulations, e.g. by Leachman et al. [14], are generally employed for the EoS, as implemented by the National Institute of Standards and Technology (NIST), which will be hereby indicated as NIST EoS.In the previous works of authors [15,16] it was shown that the deviation in calculated conditions at the release nozzle and dispersion between Abel-Noble EoS and NIST EoS is negligible for CcH2 with storage pressure 0.2e0.6MPa ab.Beyond this range the deviation in calculated hydrogen density becomes significant and grows with the increase of pressure.Thus, NIST EoS shall be used for calculations.Furthermore, the heat transfer through the wall of a release pipe connecting the storage system to the nozzle may affect the cryogenic flow characteristics.The numerical study carried out by Cirrone et al., in 2022 [17] demonstrated that the effect of heat transfer in the release pipe exposed to ambient air on hydrogen flow parameters shall be taken into account to accurately reproduce the characteristics of resulting hydrogen jet fires.In 2021, Venetsanos et al. [18] developed a simplified 1D transient model to account for the discharge line effects, i.e. pressure losses and heat transfer, during blowdown of ambient and cryogenic hydrogen storages.The model predictions were compared against experimental data on blowdown of hydrogen storage with initial temperature equal to 80 K and 300 K, and initial pressure equal to 20 MPa.However, the experimentally measured time history of temperature and pressure in the storage tank were used as an input to the model, preventing validation of the combined non-adiabatic blowdown and discharge line modelling, and its application for arbitrary initial conditions of the hydrogen storage.
No validated models are available to accurately represent the blowdown dynamics of CcH2 tanks.The present study proposes a new physical model expanding the work [12].In the proposed formulation, the non-ideal behaviour of CcH2 is taken into account by implementing the EoS with highaccuracy Helmholtz energy formulations [14] and properties from CoolProp open source database [19].The non-adiabatic blowdown model accounts for the conductive heat transfer through the storage tank wall by solving an unsteady 1D heat transfer equation.The model includes the convective heat transfer at the internal and external tank wall interfaces, i.e. hydrogen/wall and wall/external atmosphere, respectively, in either natural or forced convection regimes.Heat transfer coefficients are calculated through Nusselt correlations.The physical model considers the effect of heat transfer through the release pipe wall exposed to ambient air to accurately evaluate the hydrogen temperature, pressure and mass flow rate at the release nozzle.The model performance is assessed through comparison with experimental measurements of temperature and pressure during blowdown of hydrogen storage tanks at initial ambient and cryogenic (80 K) temperature.Sixteen experimental tests performed within PRESLHY project on "Pre-normative Research for Safe use of Liquid Hydrogen" [20,21], with initial storage pressure in the range 0.6e20.0MPa ab and release diameter in the range 0.5e4.0mm, were used for the model validation.

Validation experiments
Tests performed on the DISCHA facility by Pro-Science within the PRESLHY project [20,21] were used here for the validation of developed non-adiabatic blowdown model.The tank was made of stainless-steel and had volume of 2.81 L. The cylindrical tank had internal diameter D int ¼ 160 mm and internal height of 140 mm.Wall thickness was 30 mm at the top and bottom of the tank, whereas it was 20 mm at the vertical wall.
The tank was exposed to ambient air for the ambient temperature release tests (Fig. 1, left).For the cryogenic release tests, the tank was placed in a cooling box and immerged in a liquid nitrogen (LN 2 ) bath with temperature equal to 77 K (Fig. 1, centre).Fig. 1 (right) shows a scheme of the experimental facility and equipment for the cryogenic tests.The hydrogen release line (overall system) was located at height of 30 mm from the tank bottom.The release line included a tubular connection (kept as short as possible) between the tank and the release valve, which were immerged into the LN 2 bath.This was followed by the discharge pipe, exposed to ambient air for both ambient and cryogenic temperature tests, ending into a circular nozzle.This pipe had length equal to 55 mm, internal and external diameters equal to 10 and 12 mm respectively.Four nozzles with circular apertures of diameter equal to 0.5, 1.0, 2.0 and 4.0 mm were used in the experiments.
Several ports on the top of the tank allowed the connection to pressure and temperature sensors.A static pressure sensor in the filling line was used to measure the pressure inside the tank during the release test.Two sets of thermocouples were installed inside the tank at different heights to capture the temperature behaviour during the tests.The first set included three closed standard type K thermocouples with diameter 0.33 mm and a sensitive tip covered by thin stainless-steel shell (indicated as TC1, TC2, TC3).TC1 was placed at a height of 30 mm from the bottom of the tank, whereas TC2 and TC3 at heights of 70 mm and 110 mm respectively.The second set included three open thermocouples where the stainless-steel shell of the sensitive tip was removed and with diameter 0.25 mm (indicated as TC1o, TC2o and TC3o).Both sets were installed in comparable positions inside the vessel.Two further closed thermocouples of 1 mm diameter with stainless-steel shell cover were placed in the discharge line: one sensor was placed into the hydrogen flow after the release valve (TC4) and another sensor was inserted into the stainless-steel material at the nozzle (TC nz ).More details on the experimental set-up and equipment are available in [20,21].
Sixteen tests out of the available forty-four experiments were selected to maximize the validation domain of the model by choosing the tests with minimum and maximum initial Fig. 1 e The DISCHA facility for ambient (left) and cryogenic (centre) temperature tests including the LN 2 bath cooling box and additional equipment; scheme of the DISCHA facility for cryogenic temperature tests (right) [21].
storage pressure for each release nozzle diameter and storage initial temperature.The selected sixteen tests cover initial storage pressure P s ¼ 0.6e20 MPa ab, initial storage temperature T s ¼ 80e310 K, release nozzle diameter d n ¼ 0.5e4.0mm.Table 1 shows the initial storage conditions and nozzle diameter for each validation test.The initial storage temperature in Table 1, which will be used as initial condition for calculations, is given by the average of the three closed thermocouples readings (TC1-TC3) prior to the blowdown start.

Physical model description
The present physical model advances the non-adiabatic blowdown model accounting for heat transfer through the wall of high pressure hydrogen storage tanks developed in [12,22] to extend its applicability to CcH2 releases and to account for the heat transfer through the discharge line.The non-adiabatic blowdown model is described below by highlighting the novelties of the methodology presented in this work compared to models [12,22].
Non-ideal behaviour of hydrogen is expected at high pressures and cryogenic temperatures.The Abel-Noble EoS was proved to work well for high-pressure hydrogen storage tanks at initial ambient temperature [12,22].The aim of this research is to extend the applicability of the non-adiabatic blowdown model to cryogenic hydrogen storages and validate it against available tests with pressure up to 20 MPa.The developed model accounts for the non-ideal behaviour of hydrogen through the EoS and fluid properties based on highaccuracy Helmholtz energy formulations [14], instead of the Abel-Noble EoS, by employing the opensource CoolProp Cþþ library [19].This library allows to calculate hydrogen properties by knowing two variables of its thermodynamic state.
The first law of thermodynamics is used to assess the change of storage conditions during blowdown.Fig. 2 shows the scheme of a tank wall and the parameters affecting the heat transfer through the wall.The rate of the hydrogen internal energy, U, change in the tank is calculated from the rate of heat, Q, transfer to/from hydrogen through the tank wall and the rate of enthalpy exiting the tank by hydrogen outflow, h out [12],: In Eq. ( 1), the rate of heat transfer by convection at the internal wall is calculated as [12]: where A int is the internal area of the storage tank, T wðintÞ is the temperature of the tank wall in contact with hydrogen, T 1 is the temperature of hydrogen in the tank and k int is the heat transfer coefficient at the internal tank wall as per Eq. ( 3).The following sections will show the equations and procedures required for the calculations of the parameters in Eqs. ( 1) and ( 2).The tank and release system considered in the study are schematically shown in Fig. 3.The under-expanded jet theory [13] allows to calculate hydrogen parameters in the storage tank (location "1" in Fig. 3), and flow conditions at the exit (location "3") of the real nozzle ("2e3") and at the exit (location "4") of the notional nozzle ("3e4") during the tank blowdown.However, in the validation tests hydrogen is released from the vessel into atmosphere through a release system and discharge pipe, as shown in Fig. 1 (right).The connection between the tank and the valve in Fig. 1 (right) was kept as short as possible and all these components were immerged into the LN 2 bath for the cryogenic tests.The heat exchange area for this release line is significantly smaller than the surface of the storage tank exposed to the same LN 2 bath temperature.Thus,  it is considered that the heat exchange taking place through the release line immerged in the LN 2 bath can be neglected.Section Heat transfer through the discharge line wall discusses and justifies this assumption through comparison of experimental measurements of temperature (Fig. 4).The heat transfer between the metal of the valve under LN 2 temperature and metal of the discharge pipe under atmospheric temperature was neglected due to comparatively small area of this conjunction.After the valve, hydrogen flows through the discharge pipe ("1e2") exposed to ambient air, and thus subject to strong heat transfer that cannot be neglected especially for cryogenic tests.Thus, the under-expanded jet theory [13] cannot be applied in a straight forward way and must be expanded to account for the heat transfer through the discharge pipe and non-ideal gas behaviour by the NIST EoS based on high-accuracy Helmholtz energy formulations [14].More details are provided in Sections Heat transfer through the discharge line wall and Under-expanded jet theory with inclusion of NIST EoS respectively.

Convective heat transfer coefficient at wall boundaries
The convective heat transfer inside the tank and within the discharge pipe is calculated according to the convection regime: natural, forced or combined.The regime is defined by the ratio of the Grashof to Reynolds number to determine the corresponding Nusselt number, Nu, following the methodology [12].The convective heat transfer at the internal tank wall is then calculated as [11]: The full set of equations to characterise the convective heat transfer at internal tank wall surface can be found in [11].l g is provided by CoolProp database for the given conditions (Section Calculation procedure and assumptions for details).

Heat conduction through the tank wall
The model solves the unsteady heat conduction equation through a tank wall exposed on one side to hydrogen at temperature T 1 and on another side to the surrounding air at ambient temperature T ext , as shown in Fig. 2. The 1D heat conduction equation is applied [23], and the boundary conditions at internal and external surfaces of the tank are defined as in [12] where the full set of equations is given.

Heat transfer through the discharge line wall
CoolProp library [19] implementing the NIST EoS [14] allows to calculate the thermodynamic state and properties of hydrogen by knowing two parameters of a single-phase fluid.Thus, no further considerations or expressions are needed to determine the flow properties.For each time t during calculations, P 1 and T 1 are used to calculate r 1 , h 1 , s 1 .If the discharge pipe between storage tank and nozzle is not insulated, heat transfer through the pipe wall may strongly affect the flow parameters at the real nozzle exit.For such release at conditions close to the saturation line it is important to determine if the flow is single-phase (gas) or two-phase (liquid and gas).The developed model takes into account the heat transfer through the release pipe wall (see Fig. 3).Due to the presence of a nozzle of smaller diameter at the pipe end downstream, it is assumed that P 2 ¼ P 1 .The first and the second law of thermodynamics are combined to assess the effect of heat transfer on the fluid properties: The heat transfer through the discharge pipe wall is calculated at each time step t as: Here k int;pipe is the convective heat transfer coefficient calculated for either forced, combined, or natural convection [12]; A int;pipe is the internal surface of the pipe; and T w;pipeðintÞ is the temperature at the pipe wall surface interfacing the  hydrogen, which is calculated by the 1D heat conduction equation.For the present case, it will also be considered the assumption of T w;pipeðintÞ as constant in time, without solving the 1D heat conduction equation through the pipe wall, to simplify the model and speed up calculations.The effect of this assumption on results and calculation time will be assessed in Section Effect of heat transfer at the discharge pipe wall and assumptions for calculation T 1 is the temperature of hydrogen in the storage tank, as it is assumed to be equal to that of the flow at the entrance of the pipe exposed at the outer surface to air at ambient temperature.The validity of this assumption for the investigated cryocompressed releases is supported by the experimental evidence.Fig. 4 compares the history of the experimental measurements at the three thermocouples inside the tank (TC1, TC2 and TC3) and their average against measurement at the thermocouple (TC4) located downstream the release line and valve placed inside the cooling box and thus LN 2 bath, prior to the entrance of the 55 mm pipe exposed to ambient air.It is possible to observe that, as expected, temperature at TC4 follows the trend of TC1 during the initial stage of the release (up to about 50 s), as this thermocouple is located closer the bottom of the tank.Afterwards, temperature at TC4 gets closer to the average temperature inside the storage tank as consequence of a stronger mixing inside the vessel.Overall temperature at TC4 is well within the experimental variation of temperature inside the tank, confirming the correctness of assuming temperature at the entrance of the pipe equal to that in the storage tank for the use of the model in the studied scenario.Finally, it is possible to use the energy conservation equation to retrieve the thermodynamic state h 2 at the end of the pipe, prior to enter the nozzle section [24]: with The transient values of hydrogen mass flow rate _ m 3 and hydrogen density at the pipe exit from the previous time step are used to solve Eq. (7).Value of h 2 are calculated from Eq. ( 6).Afterwards, this parameter is used as input to CoolProp database along with pressure P 2 ¼ P 1 to estimate the remaining parameters of the hydrogen flow at the end of pipe prior to enter the nozzle section: T 2 ; s 2 and r 2 .

Under-expanded jet theory with inclusion of NIST EoS
The under-expanded jet theory in [13] is expanded to take into account the heat transfer through the discharge pipe and nonideal gas behaviour through the NIST EoS [14].Hydrogen flows through the pipe in conditions of heat transfer from the surroundings.Then, due to short length of the nozzle at the end of pipe, we assume that hydrogen undergoes isentropic expansion in the short real nozzle, i.e. s 2 ¼ s 3 .The flow is choked at the real nozzle exit, i.e. the velocity is equal to local speed of sound.The energy conservation equation is employed to calculate conditions at the real nozzle exit: The equation is solved using an iterative algorithm.Temperature decreases gradually by DT along the isentropic transformation from the conditions at location 2 to 3, i.e. s n 3 ðT; PÞ ¼ s 2 ðT 2 ;P 2 Þ.For each iteration n, the enthalpy, h n 3 ðT n 3 ;s 3 ¼ s 2 Þ, and the speed of sound, u n 3 ðT n 3 ;s 3 ¼ s 2 Þ, at the real nozzle exit are determined by the CoolProp database by using the condition s 3 ¼ s 2 and T n 3 , i.e. temperature at the iteration n.At each iteration n, the energy conservation equation is solved as follows, until a solution is found: The algorithm stops when the equation of energy conservation is satisfied with a given tolerance, D. Details for the choice of tolerance D and DT are given in Section Calculation procedure and assumptions.Temperature T end 3 , which provides the resolution of the equation for energy conservation, is the jet temperature at the real nozzle exit.Knowing T 3 and s 3 , it is possible to determine from the CoolProp database all other properties of interest, i.e. v 3 , P 3 and r 3 .The mass flow rate is calculated as In the case of the discharge pipe absence or if the pipe is properly vacuum insulated, the model can be amended to evaluate direct expansion from the storage to the nozzle and skip steps described in Section Heat transfer through the discharge line wall.
The notional nozzle concept is proved to be efficient to save computational time for numerical simulations of underexpanded jet dispersion by orders of magnitude.The expansion of the flow in the notional nozzle from location "3" (real nozzle exit) to ambient pressure (P 4 ¼ P amb Þ at location "4" (notional nozzle exit) assumes, as always in our underexpanded jet theory, the conservation of energy and speed of sound at the notional nozzle exit: The equation is solved using an iterative algorithm: temperature T 4 is varied iteration after iteration by a given DT.For each iteration n, the speed of sound v n 4 ðT n 4 ; P 4 Þ is calculated from CoolProp database.The parameters are substituted in Eq. ( 10) and the process is repeated until the balance is satisfied within a given tolerance, e.g.0.1 kJ/kg.

Calculation procedure and assumptions
The first law of thermodynamics differentiated in time can be used to calculate the specific internal energy with advancement of time t þ Dt from parameters calculated at the time step t: where : m tþDt Parameters needed for determination of k t int , i.e., b; m g , l g , c p;g , are provided by CoolProp database for T t 1 and P t 1 .This allows to calculate u tþDt

1
. The density of hydrogen in the tank at the time t þ Dt is calculated as: Now u tþDt  .Sections Heat conduction through the tank wall and Under-expanded jet theory with inclusion of NIST EoS describes how to calculate the temperature at the internal tank wall, T w ðintÞ , and mass flow rate of hydrogen, _ m t 3 .The initial conditions in the storage tank at time t ¼ 0 s are provided in Table 1.Only for the first time step, time ¼ Dt, parameters at the release nozzle are initialised by applying the under-expanded jet theory with NIST EoS and without inclusion of heat transfer in the pipe.The calculated hydrogen mass flow rate and velocity are used as input to calculate the heat transfer rate through the pipe wall and the associated hydrogen mass flow rate in the iterative process.Afterwards, the calculations proceed as shown in Fig. 5.The entire calculation algorithm is implemented in MATLAB [25], with inclusion of open source Cþþ library CoolProp, implementing NIST EoS [14] and transport properties for hydrogen [19].
Calculation of the heat transfer through the tank wall requires the knowledge of the heat transfer coefficient at external tank wall.This parameter is assumed to be constant and equal to 6 W/m 2 /K for air at ambient temperature [22].In the cryogenic tests, the stainless-steel tank is immersed in a LN 2 bath and the corresponding external heat transfer coefficient is considered to be equal to 120 W/m 2 /K [26].The same authors reported a value of 245 W/m 2 /K for nucleate boiling regime.The assessment of the effect of variation of this value for Test 22c did not demonstrate any significant difference in results.The wall thickness was considered to be uniform throughout the tank and equal to 30 mm, as per the tank walls with largest surface exposed to external conditions.An assessment for Test 22w was performed to compare the effect of the wall thickness when changed from 30 mm to 20 mm (thickness of the vertical walls).The resulting dynamics of pressure and temperature in the storage tank did not show to be affected by the change in wall thickness.It is considered that since the largest gradient of temperature is towards the tank internal surface, the change in wall thickness between the two comparable values of 20 and 30 mm is not enough to affect the results.The tank wall is made of stainless steel 1.4571.The properties for the material for the ambient temperature tests are: density r w ¼ 8000 kg/m 3 [27], specific heat c p;w ¼ 500 J/kg/K [28] and thermal conductivity l w ¼ 16.3 W/m/K [29].The stainless-steel specific heat and thermal conductivity decrease with the temperature.For the cryogenic test the material properties at 80 K are: c p;w ¼ 200 J/ kg/K [30,31] and l w ¼ 9.0 W/m/K [30].The experiments were performed in sequence, so, in some of the cryogenic temperature tests, the discharge pipe wall exposed to the atmosphere was at an initial temperature below ambient, being cooled down from the previous tests and by conduction from the valve connection located in the cooling box with LN 2 .To facilitate the model validation, the pipe wall temperature at time t ¼ 0 s, when solving Eq. ( 5), is taken from the experiment.A discharge coefficient, C d , is applied to account for friction and minor losses in the system.The solution of the non-adiabatic blowdown transient, Eqs.11e13, uses a time step in the range Dt ¼ 0.01e0.05s depending on the test initial conditions and expected blowdown time.For the tests with larger nozzle diameter (d n ¼ 2e4 mm) expecting a blowdown time below 10 s, time step Dt ¼ 0.01 s is applied.To assess if this value provides a converged solution, a sensitivity study was performed for Test 22w (P 1 ¼ 20.03 MPa ab, T 1 ¼ 300.9 K, d n ¼ 4 mm) by reducing Dt from 0.01 s to 0.001 s.The results for a C d ¼ 0.8 did not show any appreciable difference, whereas the calculation time increased from 9.3 min to almost 100 min on a quad-core machine.For the tests with smaller diameter (d n ¼ 0.5e1.0mm) expecting a blowdown time in the range 30e200 s depending on the nozzle diameter and initial pressure, a time step Dt ¼ 0.05 s is applied.The convergence study was performed by reducing Dt from 0.05 s to 0.02 s for Test 8c (P 1 ¼ 20.1 MPa ab, T 1 ¼ 84.8 K, d n ¼ 1 mm).Results for a C d ¼ 0.7 did not present any significant difference.Thus, a solution was considered as converged for Dt ¼ 0.05 s (reducing the calculation time by a factor of 1.6).A further parameter to be considered is the number of nodes z in the spatial discretization of the tank wall (Dx).It should be reminded that the time step Dt shall be chosen to fulfil the condition [11]:   5 e Calculation algorithm of the non-adiabatic blowdown model including heat transfer through the pipe wall and NIST EoS.Note: * -P lim is a limit pressure applied to stop the algorithm calculations.In the present case P lim is taken as 0.3% more than P amb .0.75 mm respectively.The analysis was performed for Test 22w (P 1 ¼ 20.03 MPa ab, T 1 ¼ 300.9 K, d n ¼ 4 mm).The resulting pressure and temperature dynamics inside the tank did not show any relevant difference, as well as the calculation time.On the other hand, if Dt ¼ 0.0001 s is applied, the calculation time would double when changing z from 10 to 40.Thus, z ¼ 20 was chosen as a compromise between the accuracy and calculation time to ensure its validity also for longer releases with smaller nozzle diameter.In the present case, for a wall thickness L t ¼ 30 mm, z ¼ 20 (Dx ¼ 1.5 mm), the condition of Eq. ( 14) Dt≪ 0.276 s is respected, as the chosen Dt varies in the range 0.01e0.05s.
A final assessment has been performed on the choice of tolerance D and DT on the solution of the under-expanded jet theory equations in Section Under-expanded jet theory with inclusion of NIST EoS.A parametric analysis has been performed for Test 22w (P 1 ¼ 20.03 MPa ab, T 1 ¼ 300.9 K, d n ¼ 4 mm) by changing D from 1000 J/kg to 10 J/kg and DT from 0.01 K to 0.001 K.The variation did not lead to significant difference in the results, despite an increase in calculation time from 9.3 min to approximately 80 min.However, to adapt the model to also lower pressures (about 0.6 MPa ab), it was considered to have a good compromise between solution accuracy and calculation time for D ¼ 100 J/kg and DT ¼ 0.001 K.It should be highlighted that the optimum combination of discretization parameters depends on the specific problem, thus a convergence analysis should be performed when applying the model to different test conditions.
Table 2 summarizes the parameters required as input to the model.The column "Input value" reports an example of calculation for the cryogenic Test 8c.

Results and discussion
The developed non-adiabatic model for CcH2 provides as output the dynamics of the following quantities during the blowdown of hydrogen from a storage tank: temperature, pressure and density in the tank; mass flow rate; temperature, pressure, density and velocity of hydrogen at the real and notional nozzle exits.The calculations of temperature and pressure dynamics inside the storage tank by the physical model are validated against experimental data in the following sections.

Comparison of three blowdown models
Fig. 6 shows the comparison between the experimental overpressure (left) and temperature (right), where TCave is averaged over the reading of three "closed" thermocouples (TC1, TC2, TC3), transients and three blowdown models: adiabatic model with use of Abel-Noble EoS (Model AD-AN) [13]; previously published non-adiabatic model accounting for heat transfer through the tank wall only and using Abel-Noble EoS (Model HT-AN) [12]; and developed in this study the nonadiabatic blowdown model taking into account heat transfer through the tank and discharge pipe walls and with implementation of NIST EoS (Model HT-NIST).
The overpressure dynamics does not show significant difference between models HT-AN and HT-NIST, signalling that the choice of EoS does not affect greatly the modelled overpressure dynamics for the selected Test 8w with initial storage temperature close to the ambient temperature.Both models closely reproduce experimental pressure transient for C d ¼ 0.7.The adiabatic model AD-AN demonstrates somewhat faster decrease of overpressure.The temperature dynamics is closely reproduced by both models with heat transfer HT-AN and HT-NIST and is not reproduced at all by the adiabatic AD-AN model.The developed model HT-NIST practically overlaps with the experimental average temperature curve in the first 10 s of the release, corresponding to the time of steepest variation of parameters inside the tank, whereas the model HT-AN results in a slightly lower temperature curve during this period.This is considered to be mainly associated with the use of constant parameters in model HT-AN (e.g., g, b, etc.) during the steep change of overpressure inside the tank, rather than with the use of NIST EoS and inclusion of heat transfer through the release pipe wall for this test starting from an initial storage temperature close to the ambient temperature.Model HT-AN converges to results of model HT-NIST after 10 s, and both models show experimentally observed increase of temperature due to heat transfer through the tank wall.The adiabatic blowdown model AD-AN demonstrates continuously decreasing temperature inside the tank down to about 80 K at 30 s, i.e. drastic difference to the models with heat transfer showing close to the test temperature of about 260 K at the same time 30 s.
The calculation time is different for three models.The model AD-AN requires under 1 min for calculations.The calculation time for the model HT-AN is almost 20 min, and it increases by more than an order of magnitude for model HT-NIST to 350 min.This is a consequence of the less efficient iterative algorithms and calculations needed for the solution of the under-expanded jet theory equations with NIST EoS.Despite the advantages in the calculation time, the adiabatic blowdown model demonstrates unacceptable performance in assessment of temperature dynamics within the tank during blowdown.On the other hand, model HT-AN [12] provides a good accuracy of both the pressure and temperature dynamics in the tank, whereas requiring a significantly lower calculation time compared to model HT-NIST.Thus, it is suggested that model HT-AN [12] can be efficiently used for calculation of non-adiabatic blowdown for hydrogen storage tanks with initial temperature close to ambient.However, the Abel-Noble EoS is not applicable to high pressure cryogenic hydrogen storage, when there is a need to use the developed in this study model HT-NIST.Thus, the model HT-NIST will be validated against blowdown experiments at both ambient and cryogenic temperatures in the remaining parts of the paper, and it will be denoted as "the model".

Effect of heat transfer at the discharge pipe wall and assumptions for calculation
The solution of the 1D heat conduction equation for the discharge pipe wall requires a low time step to satisfy the condition of Eq. ( 14) for the wall thickness of 1 mm used in the validation experiments.Let us consider a cryogenic test and 20 nodes across the discharge pipe wall, similarly to the discretization parameter employed for the 30 mm thick tank wall.The resulting Dt is 0.00025 s, which is 40 and 200 times lower than the time steps suggested for the analysis depending on the test selected for calculations, i.e.Dt ¼ 0.01 s and 0.05 s respectively.The calculation time would increase by the same order of magnitude.Given the thickness of the discharge pipe wall of only 1 mm, a lower number of nodes can be used to maintain a spatial discretization (Dx) comparable to that of the 30 mm thick tank wall.An assessment is performed with 5 nodes across the pipe wall, which requires a time step lower than 0.004 s, for Test 8c (P 1 ¼ 20.1 MPa ab, T 1 ¼ 84.8 K, d n ¼ 1.0 mm) and C d ¼ 0.7 (curve "Model, Tw ¼ transient" in Fig. 7).In this case, temperature at both the internal and external pipe wall node boundaries with hydrogen and ambient air respectively, rapidly decreases to that of hydrogen in the storage equal to about 60 K within 6 s, to then follow the same trend beyond this time (Fig. 7d).Fig. 7 compares the overpressure and temperature dynamics inside the storage tank for two cases with the assumption of constant temperature at the discharge pipe wall in solving Eq. ( 5), respectively equal to either 222.9K (experimental measurement of temperature at the nozzle wall at time 0 s) or 300.0 K (extreme case with pipe temperature equal to ambient).The higher is the heat flux through the discharge pipe wall, the slower is the overpressure and temperature decrease inside the storage tank.Maximum relative variation between mass flow rate between the cases with transient or constant temperature T w ¼ 222.9 K of the pipe wall is about 5% up to 15 s.Beyond this time, the mass flow rate is higher for increasing heat flux at the pipe wall.The variation in calculated mass flow rates increases with time, in particular as the calculation approaches the end of the storage tank blowdown.Nevertheless, a significant effect is observed on the calculation time due to the different requirements in time step.The case with transient temperature at the wall (Dt ¼ 0.004 s) requires approximately 46 h for the calculation of the entire storage tank blowdown dynamics, whereas the case with constant temperature at the wall (Dt ¼ 0.05 s) requires about 4.5 h.The relative difference in calculated temperature is within 13% at the real and notional nozzles until the jets are still under-expanded (up to about 55 s).On the other hand, a maximum relative difference for hydrogen velocity is below 8% at the real nozzle and below 6% at the notional nozzle.The relative difference in notional nozzle diameter increases from 0.1% at the start of blowdown at time 0.05 s to 11% at the end of the process at time 55 s.Overall, the difference in calculations is considered to be acceptable within engineering accuracy for the studied case, including pipe parameters.Given the significant saving in computational time and acceptable accuracy, the assumption of constant temperature at the pipe wall equal to experimental measurement at the time 0 s will be accepted in the following calculations.Nevertheless, for different release geometries and scenarios a transient solution of temperature across the discharge pipe wall, which shows exact reproduction of the pressure dynamics, shall be considered to obtain accurate estimations.

Validation of the model
The discharge coefficient, C d , is applied in calculations to account for friction and minor losses in the piping system and real nozzle compared to the ideal case of no losses with C d ¼ 1.For each of the simulated tests, different discharge coefficients are applied to find the optimum characteristic for this experiment.Fig. 8 shows an example of calculations for Test 4c.For this test the optimum value is C d ¼ 0.8.
Figs. 9 and 10 compare pressure and temperature dynamics for eight tests at cryogenic temperature.The calculated temperature is compared with experimental temperature, TCave, averaged arithmetically over the readings of three type K "closed" thermocouples (i.e., with a stainless-steel sensitive tip) located at different heights in the tank (from bottom to top: TC1, TC2, TC3).The developed model reproduces well the experimental pressure and temperature dynamics for all these tests.The arithmetical averaging of the three experimental temperature readings, TCave, instead of mass-averaging may be a cause of the generally slight deviation of TCave from the calculated temperature in the storage tank.Tests 16c and 22c with the largest diameter (d n ¼ 4.0 mm) result in somewhat lower temperature compared to the experimental one (see Fig. 10).This may be associated with the inertia of the "closed" thermocouples, which may not well capture the steeper variation of temperature associated with these faster blowdowns.A deeper analysis of this issue will be carried out in Section The effect of thermocouple inertia.Tests 16c and 22c show a larger deviation among the records of the three thermocouples inside the tank, signalling a larger non-uniformity of temperature distribution within the tank.Tests with lower initial storage pressure (about 0.6 MPa ab) present a certain level of noise when approaching the ambient pressure, e.g.Tests 9c and 16c after about 5 s and 1.5 s respectively, whereas as expected calculations tend to zero.
The optimum discharge coefficients for the whole set of tests are found to be in the range C d ¼ 0.6e0.8.The maximum calculated hydrogen mass flow rate is at the start of blowdown (defines the time step Dt) and varies greatly from 0.1 to 6.6 g/s when increasing diameter from 0.

The effect of thermocouple inertia
Predictions of pressure and temperature dynamics in the storage tank are seen to agree well with experimental measurements.Few exceptions are given by the releases with larger diameter, e.g.Test 16w (P 1 ¼ 0.59 MPa ab, T 1 ¼ 296.0 K, d n ¼ 4.0 mm).This deviation is deemed to be caused by the thermocouple inertia being comparable with the blowdown duration.Indeed, Test 16w is one of the tests showing the largest difference in temperature measurements by the "closed" and "open" type of thermocouples, both used in this experiment conversely to other tests."Open" thermocouples had the stainless-steel tip removed, decreasing the thermocouple inertia but as well reducing measurements accuracy for cryogenic temperatures [20].
Fig. 13 shows the comparison between the predicted temperature in the storage tank and experimental measurements by the two thermocouple types.The "open" thermocouples measurements give better agreement with the model calculations due to reduced sensors inertia.This is important for such a short blowdown duration.However, these sensors may lose accuracy for reduced temperatures, including cryogenic, and thus "closed" thermocouples were used in the experiments and thus in the model validation process.

The use of notional nozzle exit parameters in CFD simulations
In the framework of CFD modelling of under-expanded jets, the notional nozzle concept helps to use a larger size of inflow  The VSM uses a constant size of the volumetric release area, but allows simulation of variable in reality inflow parameters.First of all, the dependences of the flow parameters at the notional nozzle exit, i.e. hydrogen mass flow rate _ m, notional nozzle velocity y 4 , enthalpy h 4 , turbulent kinetic energy k TKE;4 and dissipation rate ε 4 , are obtained and approximated as a function of the blowdown time.Then, the volumetric sources are defined to reflect the influx of all mentioned above transported variables in the form of equalised rate of their generation, in the defined by the VSM volume of the domain, in the conservation equations.The inflow of parameters to be accounted for in the conservation equations of mass, momentum in the jet direction (x-axis), energy, hydrogen mass fraction, turbulent kinetic energy and its dissipation rate are: x À momentum inflow _ m$y where DV is the volume designated by the VSM for application of the source term and S are volumetric source terms for each solved conservation equation.It should be noted that the source terms for mass conservation and for hydrogen mass fraction are the same in this problem formulation as the only mass coming through the nozzle is hydrogen mass.Accordingly, the source terms for the governing equations are: The applicability limits and the validation of the VSM are reported in [13] using an example of under-expanded hydrogen jet experiments performed by HSE [32], particularly Run 7 with hydrogen release from a storage at quasisteady state pressure of P 1 ¼ 10.0 MPa and temperature T 1 ¼ 14 C through a real nozzle of 3 mm diameter and estimated mass flow rate _ m ¼ 0:045 kg/s.Simulations were carried out for different sizes of the volumetric source part of the domain.The measured hydrogen concentrations and simulation results shown in Fig. 14 are in good agreement with the experimental data provided that the ratio of the release volumetric source size to the notional nozzle diameter is up to 4. This means that at the start of CFD simulations the volumetric release source size can be chosen equal to the size of the notional nozzle and stay constant during the blowdown process even though the notional nozzle size decreases with time.The simulation results will be pretty accurate until notional nozzle exit diameter decreases by 4 times below the size of the chosen release volume.The described above VSM was used for simulations of non-adiabatic release from the cryo-compressed hydrogen storage tank.ANSYS Fluent version 2020R2 served as a platform to perform numerical simulations.The results of CFD simulations on release and dispersion of a hydrogen jet generated in Test 22c are compared to experimental data in an inter-comparison study among partners of the PRESLHY project [33].Here only the results of the CFD modelling performed by Ulster University on the temperature distribution in the jet are reported.The goal is to demonstrate the capability of the developed in this study non-adiabatic blowdown model, in the conjunction with the VSM model [13], to properly simulate the release characteristics of transient cryogenic hydrogen jets and thus hazard distances.
Fig. 15 shows the comparison between the CFD simulations and experimental measurements of temperature at different locations within the hydrogen jet.The thermocouples were located at different places relatively to the jet axis, at different heights from the ground and axial distances from the nozzle.At sensors closer to the release point (TC8 and TC5, respectively at axial distances of 200 mm and 250 mm), simulations predict somewhat colder temperatures than measured in the experiment, possibly due to the sensors inertia.However, a good agreement is achieved for the sensors located further downstream the nozzle (TC6 and TC7, respectively at axial distances of 750 mm and 1750 mm), despite the steeper temperature decrease during 0.5 s after the start of the release due to difference between the instantaneous and inertial opening of the release valve in simulations and test correspondingly.[32] and simulated hydrogen volumetric fraction for various sizes of the volume source in the VSM [13].

Fig. 2 e
Fig. 2 e Scheme of a tank wall and parameters used in the conjugate heat transfer calculations.

Fig. 4 e
Fig. 4 e Experimental hydrogen temperature in the tank (thermocouples TC1, TC2, TC3 and their average TC ave ) and after the release valve (TC4) placed inside the LN 2 bath for Test 8c (P 1 ¼ 20.1 MPa ab, T 1 ¼ 84.8 K, d n ¼ 1.0 mm).

Fig. 3 e
Fig. 3 e Schematic of the model: 1 e storage tank with tabular connection to valve under LN2 temperature; 1e2 e discharge pipe under atmospheric temperature; 2 e end of pipe prior to the nozzle; 2e3 e real nozzle; 3 e real nozzle exit; 3e4 e notional nozzle; 4 -notional nozzle exit.

1 and r tþDt 1 can 1 ;
be used as input to CoolProp database to determine the thermodynamic state of hydrogen at the time þDt: T tþDt

Dt≪ 1 2 r
w c p;w l w Dx 2 : (14) Time step Dt ¼ 0.001 s was chosen for a parametric study comparing z ¼ 10, 20 and 40, corresponding to Dx ¼ 3.0, 1.5 and

Fig.
Fig.5e Calculation algorithm of the non-adiabatic blowdown model including heat transfer through the pipe wall and NIST EoS.Note: * -P lim is a limit pressure applied to stop the algorithm calculations.In the present case P lim is taken as 0.3% more than P amb .

Fig. 6 e
Fig. 6 e Experimental pressure and temperature dynamics for Test 8w (P 1 ¼ 20.19 MPa ab, T 1 ¼ 307.7 K, d n ¼ 1.0 mm, C d ¼ 0.7) against calculated by three different blowdown models.

Fig. 7 e
Fig. 7 e Effect of assumptions for calculation of heat transfer at the discharge pipe wall on the pressure (a) and temperature (b) dynamics in the storage tank; hydrogen mass flow rate (c); comparison of temperature in the storage tank, and at the internal and external discharge pipe walls (d) for Test 8c (P 1 ¼ 20.1 MPa ab, T 1 ¼ 84.8 K, d n ¼ 1.0 mm) and C d ¼ 0.7.
5 to 4.0 mm for an initial storage pressure of 0.6 MPa ab.For tests with initial storage pressure of about 20 MPa ab, the hydrogen mass flow rate varies in the range 4.3e241.7 g/s.The calculation time varies depending on the test initial pressure and nozzle diameter.The shortest calculation time is about 15 min and is achieved for Test 16c (P 1 ¼ 0.61 MPa ab, T 1 ¼ 80.2 K, d n ¼ 4.0 mm), being the test with shorter blowdown time due to the lower storage pressure and larger release diameter.Opposite, Test 4c (P 1 ¼ 20.12 MPa ab, T 1 ¼ 80.2 K, d n ¼ 0.5 mm) records the longest calculation time of about 19 h on a quad-core laptop.For the same release diameter but lower initial pressure Test 1c (P 1 ¼ 0.59 MPa ab, T 1 ¼ 83.7 K, d n ¼ 0.5 mm), calculation requires approximately 6 h.Figs.11 and 12 show the comparison of calculations against experiments with initial storage temperature equal to ambient.The comparison confirms the accurate predictive capability of the developed non-adiabatic blowdown model.Similar observations as per cryogenic tests can be made for Tests 16w and 22w with shortest blowdown time in Fig. 12 (see Section The effect of thermocouple inertia for details).The optimum discharge coefficient for all the set of ambient temperature tests is in the same range C d ¼ 0.6e0.8 as for cryogenic temperatures.The calculated hydrogen mass flow rate varies in the range 0.049e2.8g/s for initial storage pressure of 0.6 MPa ab, whereas it varies in the range 1.9e107.0g/s for pressure of 20 MPa ab.

Fig. 8 e
Fig. 8 e Effect of the discharge coefficient, C d , on the blowdown dynamics for Test 4c (P 1 ¼ 20.12 MPa ab, T 1 ¼ 80.2 K, d n ¼ 0.5 mm).

Fig. 9 e 6 Fig. 10 eFig. 11 eFig. 12 e
Fig. 9 e Validation of the developed non-adiabatic blowdown model against experiments at initial cryogenic temperature and nozzle diameters d n ¼ 0.5 mm and 1.0 mm.

Fig. 15 e
Fig. 15 e Top: temperature dynamics in the hydrogen cryogenic jet for Test 22c (P 1 ¼ 20.27 MPa ab, T 1 ¼ 80.2 K, d n ¼ 4.0 mm), CFD simulation results (dashed lines) versus experimental measurements (solid lines).Bottom: location of thermocouples relative to the jet axis (dash-dotted line on plane z ¼ 0) at different heights and axial distances (in mm).

Fig. 14 e
Fig.14e Comparison of measured[32] and simulated hydrogen volumetric fraction for various sizes of the volume source in the VSM[13].

Table 1 e
Parameters (nozzle diameter, pressure and temperature) of the sixteen validation tests.

Table 2 e
Input parameters for the developed model including heat transfer through the pipe and NIST EoS.