CFD Model of Refuelling through the Entire HRS Equipment: The Start-Up Phase Simulations

: Refuelling hydrogen-powered cars, buses, trucks, trains, ships, and planes is a technological challenge. The absence of contemporary CFD models of refuelling through the entire hydrogen refuelling station (HRS) equipment is one of the scientiﬁc bottlenecks. Detailed refuelling protocols for more than 10 kg of hydrogen, e


Introduction
Safety is the highest priority for emerging hydrogen systems and infrastructure.Hazards and associated risks of hydrogen-powered vehicles, including cars, buses, trucks, trains, ships, planes, etc. should be at least at the same level or below those for fossil-fuelled vehicles.Hydrogen is compressed for onboard storage to nominal working pressure (NWP) of 35 MPa, e.g., in the case of currently used buses, or 70 MPa, e.g., for cars, to make the driving range comparable to traditional vehicles.During refuelling, hydrogen temperature inside the onboard tanks increases as a result of heat and mass transfer phenomena, including hydrogen compression, the Joule-Thomson effect at the Pressure Control Valve (PCV), and the conversion of kinetic energy into internal energy in onboard tanks.Hydrogen temperature during refuelling should not exceed the regulated threshold of 85 • C [1][2][3].Used in the validation experiment of the National Renewable Energy Laboratory (NREL) Type IV onboard storage tanks are made of carbon fibre-reinforced polymer (CFRP) to bear the pressure load, and hydrogen-tight plastic liner to limit the permeation to the regulated level [2][3][4].The performance and integrity of these tank materials are endangered when they are exposed to a higher-than-regulated temperature threshold of 85 • C during fuelling [1].For these reasons, SAE J2601 [5][6][7], SAE J2579 [4], ISO 15869 [8] (currently withdrawn), UN ECE Regulation GTR#13 [2], and European Regulation R134 [3], limit the temperature inside the tank between −40 • C and 85 • C, the maximum fuelling pressure is limited by 125% of NWP, and the state of charge is restricted by SoC = 100% to ensure tank's safety while refuelling and afterwards.Crossing these thresholds may lead to hydrogen tank failure, including its rupture or unignited or ignited (fire) hydrogen leaks [9,10].Hydrogen pre-cooling in the Heat Exchanger (HE) is used during refuelling to keep the hydrogen temperature in the tank below the 85 • C threshold [11,12].Pre-cooling affects the fuelling station design, performance, reliability and cost [13] and should be avoided when possible.The last requires contemporary and validated modelling tools.
The SAE J2601 [5] is applicable only to light-duty vehicles (LDV) with a maximum hydrogen inventory capacity of 10 kg and NWP = 70 MPa.This document is not valid for refuelling heavy-duty vehicles (HDV), i.e., buses, trucks, coaches, trains, maritime vessels, and aeroplanes, with a larger onboard hydrogen inventory.Unfortunately, SAE J2601-2 [6] provides only high-level safety requirements for HDVs, not refuelling protocols.The same statement is valid for SAE J2601-3 [7] which is developed to provide guidance only for fuelling Hydrogen Powered Industrial Trucks including tractors, forklifts, and pallet jacks.The situation with refuelling protocols becomes critical for the success of the hydrogen economy rollout.The underlying physical phenomena of heat and mass transfer in all components of hydrogen refuelling station (HRS) and their joint effect on refuelling process safety and efficiency are yet to be fully understood and then derived principles to be applied for the development of efficient and cost-effective fuelling protocols.This requires the development and validation of contemporary CFD models with superior qualities compared to a few currently available reduced models.For example, reduced models cannot estimate temperature non-uniformity in onboard tanks, optimise pre-cooling systems, etc. Temperature non-uniformity could lead to a significant increase of localised hydrogen temperature above the calculated by reduced models bulk temperature of 85 • C and thus failure of a composite storage tank during lifetime cycling.
Most of the reduced models of fuelling include only onboard tanks [14 -18] and do not consider other equipment of HRS.To the best authors' knowledge, there is only one reduced model for modelling of refuelling through the entire equipment of HRS, i.e., H2Fills [19].This reduced model is based on the Hydrogen Refuelling Station Dynamic Simulation (HRSDS) software of Kyushu University.There are serious limitations in the use of the H2Fills reduced model compared to the potential of CFD models.These include the absence of modelling capability of the High Pressure (HP) tanks thermal behaviour during blowdown, the inability to assess the non-uniformity of temperature in onboard storage tanks, the impossibility of reproducing experimentally observed temperature drops during leak checks, etc.
Due to the usually long computational time, the CFD hydrogen fuelling studies to mention a few [1,11,12,[20][21][22][23][24][25] were focused on the thermal behaviour of onboard tanks only.In 2018 Bourgeois et al. [26] stressed that it is crucial to take into account the effect of the entire fuelling line while studying the fuelling procedure.In 2021 Kuroki et al. [19] stated that a fuelling model that accounts for the entire fuelling station has not yet been developed and proposed an H2Fills reduced model [19].The filling and emptying process of hydrogen tanks includes various complex and competing phenomena, such as convection, stratification and laminar to turbulent flow transition [13].Contemporary 3D CFD modelling allows us to naturally account for such complex phenomena and develop an in-depth understanding of the heat transfer mechanisms, e.g., using flow and temperature fields reproduction at any specific time at any specific point [13].This advantage differentiates CFD from thermodynamic models, which are essentially limited to operating averaged over a piece of equipment parameters, e.g., the entire tank.This makes thermodynamic models fundamentally unable to predict the distribution of local field variables such as temperature distribution in a single tank or tank assembly.The fundamental difference between 3D CFD and thermodynamic modelling makes CFD particularly useful when the average temperature in a tank reaches 85 • C but localised temperature could be much higher than the 85 • C limit, e.g., conformable tanks and tanks with large L/D ratio.To the best of the authors' knowledge, no CFD model of hydrogen flow through the entire HRS and FCEV onboard storage has been reported before.The presented modelling work was carried out at Ulster University as a part of a doctoral study.
The start-up phase is an important part of the fuelling protocol aiming to detect possible leaks in the system before proceeding to storage tanks refuelling.This study aims to develop a CFD model for simulations of the refuelling process through the entire set of HRS equipment, and compare simulation results against experimental pressure and temperature in onboard tanks during the start-up phase of the refuelling experiment by Kuroki et al. [19] carried out at NREL in the USA.
The Start-Up Phase of Refuelling SAE J2601 [5] states that the start-up phase begins after the user initiates fuelling and ends when the dispenser begins flowing hydrogen to fuel the vehicle (Figure 1).This start-up period can include a connection pulse, determination of vehicle tank capacity category, and a leak check [5].The connection pulse starts when a pressurised nozzle is connected to the receptacle [5].In other words, a pulse of high-pressure hydrogen is sent through the pipes toward the onboard tanks using the PCV.The PCV stops the flow for a few seconds afterwards.This specific stage measures any decrease in pressure in the fuelling line [5] before the main fuelling starts.
The start-up phase is an important part of the fuelling protocol aiming to possible leaks in the system before proceeding to storage tanks refuelling.This stud to develop a CFD model for simulations of the refuelling process through the entir HRS equipment, and compare simulation results against experimental pressu temperature in onboard tanks during the start-up phase of the refuelling experim Kuroki et al. [19] carried out at NREL in the USA.
The Start-Up Phase of Refuelling SAE J2601 [5] states that the start-up phase begins after the user initiates fuell ends when the dispenser begins flowing hydrogen to fuel the vehicle (Figure 1).Th up period can include a connection pulse, determination of vehicle tank capacity ca and a leak check [5].The connection pulse starts when a pressurised nozzle is con to the receptacle [5].In other words, a pulse of high-pressure hydrogen is sent t the pipes toward the onboard tanks using the PCV.The PCV stops the flow fo seconds afterwards.This specific stage measures any decrease in pressure in the f line [5] before the main fuelling starts.From the safety point of view, the initial connection pulse is vital as it determ it is safe to start fuelling or not.Failing to pass this check will terminate the f procedure immediately to prevent potential incidents jeopardising life safe property protection.That is why this study is focused on the start-up phase b important constituent part of the refuelling protocol.

Validation Experiment
The validation of the CFD model is performed in this study against the start-u of the experiment carried out at NREL by Kuroki et al. [19].Figure 2 shows th components used in the experiment [19].The components of HRS include tw pressure (HP) tanks with a 300-litre capacity each working as a hydrogen sou applied in the cascade configuration (one after another).Due to the short duratio connection pulse, only one HP tank is used for the start-up phase.These HP ta connected through 59 m of piping sections with 24 bends (90-degree), PCV, Ma Rate (MFM), HE, breakaway, hose, nozzle, and 5 other valves to the three 36-litre o storage tanks with the internal length to diameter ratio L/D = 3.4.The length, di and material of each piping section, along with the valve flow coefficients are avai the Appendix of the paper [19].The specification of pipes from manifolds to o tanks is not mentioned in the Appendix of the paper [19] but it is assumed to be th as in the H2Fills software demonstration example [19], which replicates the From the safety point of view, the initial connection pulse is vital as it determines if it is safe to start fuelling or not.Failing to pass this check will terminate the fuelling procedure immediately to prevent potential incidents jeopardising life safety and property protection.That is why this study is focused on the start-up phase being an important constituent part of the refuelling protocol.

Validation Experiment
The validation of the CFD model is performed in this study against the start-up phase of the experiment carried out at NREL by Kuroki et al. [19].Figure 2 shows the HRS components used in the experiment [19].The components of HRS include two highpressure (HP) tanks with a 300-litre capacity each working as a hydrogen source and applied in the cascade configuration (one after another).Due to the short duration of the connection pulse, only one HP tank is used for the start-up phase.These HP tanks are connected through 59 m of piping sections with 24 bends (90-degree), PCV, Mass Flow Rate (MFM), HE, breakaway, hose, nozzle, and 5 other valves to the three 36-litre onboard storage tanks with the internal length to diameter ratio L/D = 3.4.The length, diameter, and material of each piping section, along with the valve flow coefficients are available in the Appendix of the paper [19].The specification of pipes from manifolds to onboard tanks is not mentioned in the Appendix of the paper [19] but it is assumed to be the same as in the H2Fills software demonstration example [19], which replicates the NREL refuelling station and is applied in this study.All the valves, except the PCV, are assumed to be fully open during refuelling.
refuelling station and is applied in this study.All the valves, except the PCV, are assumed to be fully open during refuelling.The pressure control during refuelling is carried out by the PCV which regulates the downstream pressure of the PCV at a set Average Pressure Ramp Rate (APRR).Hydrogen temperature is measured in the HP tank (TE1), before and after the PCV (TE2, TE3), immediately after the pre-cooler (TE4), and in the centre of two of three onboard tanks (TE5, TE6) with the measurement accuracy for all of the temperature sensors of ±1.5 K [19].Pressure sensors with ±1 MPa accuracy are located at the exit of the HP tank (PT1), before valve 4 (PT2), and at one of the onboard tank inlets (PT3) [19].The mass flow rate was measured not directly by an MFM but by measuring the change in the total mass of the test tanks [19].Two experiments were carried out [19]: Test No.1 without leak checks during the main fuelling and Test No.2 with initial and two follow-up leak checks during the main fuelling process.The initial conditions and the results of the Test No.2 start-up phase are used to validate the CFD model developed in this study.Experimentally observed localised temperature drop during leak checks at PCV [20] is a subject of our forthcoming study that will be published in due course.
Details of experimental measurements of mass flow rate, pressure and temperature dynamics during the start-up phase are presented in the section "Simulation results and discussion" where they are compared against the performed simulations.

Calculation Domain and Parameters of HRS Components
Figure 3 shows the calculation domain that includes all the HRS components used in the NREL refuelling experiment.The domain boundaries are the internal surfaces of each component of HRS.The modelled geometry uses the component specifications given in [19], except for the HP tanks because only the volume of the HP tanks is given in the experimental paper.Due to this lack of data, the shape of HP tanks is simplified to cylinders.The dimensions of the HP tank are not available but based on the available commercial 300 litres HP tank [27] are assumed to have the length-to-diameter ratio L/D = 3.The total wall thickness of the HP tank, including both CFRP and liner, used in The pressure control during refuelling is carried out by the PCV which regulates the downstream pressure of the PCV at a set Average Pressure Ramp Rate (APRR).Hydrogen temperature is measured in the HP tank (TE1), before and after the PCV (TE2, TE3), immediately after the pre-cooler (TE4), and in the centre of two of three onboard tanks (TE5, TE6) with the measurement accuracy for all of the temperature sensors of ±1.5 K [19].Pressure sensors with ±1 MPa accuracy are located at the exit of the HP tank (PT1), before valve 4 (PT2), and at one of the onboard tank inlets (PT3) [19].The mass flow rate was measured not directly by an MFM but by measuring the change in the total mass of the test tanks [19].Two experiments were carried out [19] Details of experimental measurements of mass flow rate, pressure and temperature dynamics during the start-up phase are presented in the section "Simulation results and discussion" where they are compared against the performed simulations.

Calculation Domain and Parameters of HRS Components
Figure 3 shows the calculation domain that includes all the HRS components used in the NREL refuelling experiment.The domain boundaries are the internal surfaces of each component of HRS.The modelled geometry uses the component specifications given in [19], except for the HP tanks because only the volume of the HP tanks is given in the experimental paper.Due to this lack of data, the shape of HP tanks is simplified to cylinders.The dimensions of the HP tank are not available but based on the available commercial 300 litres HP tank [27] are assumed to have the length-to-diameter ratio L/D = 3.The total wall thickness of the HP tank, including both CFRP and liner, used in simulations is estimated as 33 mm.It is worth mentioning that the onboard tanks are modelled with the same surface area, volume, and L/D as the experimental tanks.The entire computational domain is meshed using 207,252 control volumes.The minimum orthogonal quality of mesh is 0.7 with an average of 0.97.The HP tank has meshed with 62,912 control volumes (CVs) with a minimum maximum sizes of 1 cm and 3 cm, respectively and a growth rate of control volume of 1.2.The onboard tanks are meshed with 37,850 CVs each with minimum and maxi size of 0.05 cm and 2 cm respectively.Figure 4 shows the discretisation of tanks.The valves usually have a complex geometry.In this study, the only available for valves and HE is their flow coefficient.All valves, including the PCV, are modell pipe sections with equivalent internal diameters calculated based on their flow coeffi values [19] and of 0.1 m length.The HE is modelled as a pipe of 1 m in length an internal diameter is calculated based on its flow coefficient equal to 1 [19].The equiv internal diameter of the 5 valves, PCV, HE and MFM is calculated as [28]: The HP tank has meshed with 62,912 control volumes (CVs) with a minimum and maximum sizes of 1 cm and 3 cm, respectively and a growth rate of control volume size of 1.2.The onboard tanks are meshed with 37,850 CVs each with minimum and maximum size of 0.05 cm and 2 cm respectively.Figure 4 shows the discretisation of tanks.The HP tank has meshed with 62,912 control volumes (CVs) with a minimum maximum sizes of 1 cm and 3 cm, respectively and a growth rate of control volume of 1.2.The onboard tanks are meshed with 37,850 CVs each with minimum and maxim size of 0.05 cm and 2 cm respectively.Figure 4 shows the discretisation of tanks.The valves usually have a complex geometry.In this study, the only available for valves and HE is their flow coefficient.All valves, including the PCV, are modelle pipe sections with equivalent internal diameters calculated based on their flow coeffic values [19] and of 0.1 m length.The HE is modelled as a pipe of 1 m in length an internal diameter is calculated based on its flow coefficient equal to 1 [19].The equiva internal diameter of the 5 valves, PCV, HE and MFM is calculated as [28]: The values of flow coefficients,   , were taken from [19] and the discharge coeffic was assumed to be equal to a typical value of Cd = 0.65.The calculated values o The valves usually have a complex geometry.In this study, the only available data for valves and HE is their flow coefficient.All valves, including the PCV, are modelled as pipe sections with equivalent internal diameters calculated based on their flow coefficient values [19] and of 0.1 m length.The HE is modelled as a pipe of 1 m in length and its internal diameter is calculated based on its flow coefficient equal to 1 [19].The equivalent internal diameter of the 5 valves, PCV, HE and MFM is calculated as [28]: The values of flow coefficients, C v , were taken from [19] and the discharge coefficient was assumed to be equal to a typical value of C d = 0.65.The calculated values of an equivalent internal diameter of the PCV, 5 other valves, MFM and HE are shown in Table 1.The increase or decrease of pipe diameter generates minor pressure losses in hydrogen flow.The exact geometry of valves cannot be reproduced in simulations of entire HRS fuelling due to its complexity and computational costs.The transition from one pipe diameter to another was realised in the numerical grid as a gradual change (see Figure 5).
Hydrogen 2023, 4, FOR PEER REVIEW 6 equivalent internal diameter of the PCV, 5 other valves, MFM and HE are shown in Table 1.The increase or decrease of pipe diameter generates minor pressure losses in hydrogen flow.The exact geometry of valves cannot be reproduced in simulations of entire HRS fuelling due to its complexity and computational costs.The transition from one pipe diameter to another was realised in the numerical grid as a gradual change (see Figure 5).All the HRS components between the HP tank and the onboard tanks are treated as cylindrical bodies.Kuroki et al. [19] used in their reduced model calculations of the inner diameters and lengths of components exactly as in the experiment.The applied in the reduced model by Kuroki et al. [19] external diameters of the piping sections, breakaway, hose and nozzle were adjusted to match the real weight of these components in order to have the same thermal mass of components as in the experiment.The reported in [19] internal and external diameters along with the properties of a relevant component, i.e., density, specific heat, and thermal conductivity, are used in this numerical study.Because the external diameter and the thermal mass of the valves, PCV, and MFM are not available in [19], it is assumed that these elements have the same external diameter and thermal properties as their upstream pipes.As stated in the experimental section, the bends are modelled as 90° bends with a radius of 15 mm for 5.1 mm diameter pipes and a bend radius of 24 mm for 7.9 mm diameter pipes (see Figure 6).All the HRS components between the HP tank and the onboard tanks are treated as cylindrical bodies.Kuroki et al. [19] used in their reduced model calculations of the inner diameters and lengths of components exactly as in the experiment.The applied in the reduced model by Kuroki et al. [19] external diameters of the piping sections, breakaway, hose and nozzle were adjusted to match the real weight of these components in order to have the same thermal mass of components as in the experiment.The reported in [19] internal and external diameters along with the properties of a relevant component, i.e., density, specific heat, and thermal conductivity, are used in this numerical study.Because the external diameter and the thermal mass of the valves, PCV, and MFM are not available in [19], it is assumed that these elements have the same external diameter and thermal properties as their upstream pipes.As stated in the experimental section, the bends are modelled as 90 • bends with a radius of 15 mm for 5.1 mm diameter pipes and a bend radius of 24 mm for 7.9 mm diameter pipes (see Figure 6).The increase or decrease of pipe diameter generates minor pressure losses in hydrogen flow.The exact geometry of valves cannot be reproduced in simulations of entire HRS fuelling due to its complexity and computational costs.The transition from one pipe diameter to another was realised in the numerical grid as a gradual change (see Figure 5).All the HRS components between the HP tank and the onboard tanks are treated as cylindrical bodies.Kuroki et al. [19] used in their reduced model calculations of the inner diameters and lengths of components exactly as in the experiment.The applied in the reduced model by Kuroki et al. [19] external diameters of the piping sections, breakaway, hose and nozzle were adjusted to match the real weight of these components in order to have the same thermal mass of components as in the experiment.The reported in [19] internal and external diameters along with the properties of a relevant component, i.e., density, specific heat, and thermal conductivity, are used in this numerical study.Because the external diameter and the thermal mass of the valves, PCV, and MFM are not available in [19], it is assumed that these elements have the same external diameter and thermal properties as their upstream pipes.As stated in the experimental section, the bends are modelled as 90° bends with a radius of 15 mm for 5.1 mm diameter pipes and a bend radius of 24 mm for 7.9 mm diameter pipes (see Figure 6).

Governing Equations and Numerical Details
Governing equations include unsteady conservation equations for mass, momentum, and energy: ∂ρ ∂t Hydrogen 2023, 4 The symbol "overbar" stands for Reynolds averaged parameters and "tilde" for Favre averaged parameters.
The CFD model assumes that only hydrogen substance is present in the geometry and employs the real gas equation of state (EoS) of the National Institute of Standards and Technology (NIST) to account for the thermal properties of hydrogen changing with pressure and temperature.The used CFD engine (ANSYS Fluent) dynamically loads the NIST Thermodynamic and Transport Properties of Refrigerants and Refrigerant Mixtures Database REFPROP v.7.0 [29] into the solver to evaluate transport and thermodynamic hydrogen properties.
Flow turbulence was modelled using the standard k-ε model [30]:

∂(ρε) ∂t
where is the mean rate of strain, The simulations were performed using ANSYS Fluent 2020R2 as a CFD engine [31].The SIMPLE algorithm was applied for pressure-velocity coupling.Convective terms were discretised using the pressure-based implicit solver and first-order upwind numerical scheme.The simulation of 14 s of the start-up phase time takes about 2 h on a 32-core AMD Opteron CPU running at 2.3 GHz.

Initial and Boundary Conditions
The boundaries of all pipelines, valves, and tanks are modelled as non-slip, impermeable walls.The initial conditions, including temperature and pressure inside the tanks, are the same as in the experiment [19].The initial pressure is 88 MPa in HP tanks and pipes upstream of PCV, and the pressure is 5.5 MPa in pipes from PCV down to onboard tanks.The initial temperature of hydrogen in the HP tank and its walls is 17.5 • C (equal to the ambient temperature in the HP tank location of 17.5 • C).The initial HE temperature is 23 • C and the temperature in the rest of the domain is 21 • C both for hydrogen and walls following the experimental conditions [19,32].Heat flux on all pipe surfaces was calculated using the ANSYS Fluent "shell conduction" capability.The shell conduction calculates conjugate heat transfer through walls in both normal and parallel directions to the pipe axis [31].After specifying the material properties, e.g., density, specific heat and thermal conductivity, and wall thickness for each section, Fluent automatically grows specified layers of cells, either prismatic or hexahedral, in the wall surface to simulate 3D heat conduction [31].This model accounts for the specified wall thickness of piping, relevant thermal properties of materials and convective heat transfer on the outer surface of the wall.The conduction heat transfer through the wall is governed by Fourier's law and is computed as: Fluent computes heat flux at the external walls of components as: where the heat transfer coefficient h ext = 7 W/m 2 /K is used in line with the conclusions of the study [19].

Modelling the PCV and the HE
The ANSYS Fluent fixing the values of variables method is used to control the flow through the PCV and temperature in the HE.The method is described in the ANSYS Fluent [31] as "when a variable is fixed in a given cell, the transport equation for that variable is not solved in the cell (and the cell is not included when the residual sum is computed for that variable).The result is a smooth transition between the fixed value of a variable and the values at the neighbouring cells".
Using this method, the mass flow rate is controlled by changing hydrogen velocity in the PCV using the User Defined Function (UDF) capability of Fluent.The UDF is programmed in a way that it computes the mass flow rate in the PCV, The velocity of hydrogen in the PCV is changed dynamically based on the value of dm norm to control the required hydrogen mass flow rate in the PCV.
Details of the cooling system of the HE are not described in the experimental paper [19].For the modelling purposes in this study, the applied equivalent diameter and flow coefficient of the HE were as described in [19].The length of HE was selected as 1 m to expand possible options in the modelling of the HE.For the HE, it is possible to directly control the temperature of cells in the numerical HE.Using the UDF, the temperature in the fluid passing the HE zone is set to match the experimentally measured temperature.

Simulation Results and Discussion
The simulated parameters of the start-up phase were compared against the experimental data [19] to validate the simulations.The raw experimental records were not available and the transients of measured parameters from Test No.2 of the study [19] were digitised and used for comparison with the simulations of this study.
Figure 7 shows the dynamics of the experimental mass flow rate (left) and temperature after the HE (right) during the start-up phase of 14 s duration.These parameters were used as an input to simulate pressure and temperature in onboard hydrogen storage tanks.In the experiment, the mass flow rate is measured by tank mass change.In the simulation, the mass flow rate is calculated in the PCV, and its assessment in the MFM or other locations showed no difference as expected.The UDF can define the mass flow rate in the PCV that closely follows the experimental mass flow rate (see Figure 7, left).This model capability could be useful for the development of refuelling protocols for arbitrary design of HRS, initial parameters of components and ambient conditions.Figure 7 (right) shows the experimental and simulated temperature at the exit from the cooling system (measurement point TE4 in PID, see Figures 2 and 3).The simulated temperature at the exit from the HE reproduces closely the experimental data with an acceptable engineering accuracy of ±0.5 K.
close.This could explain why after the mass flow starts to decrease, the temperature at the HE exit stops the decrease and starts to grow.The integration of mass flow rate over time (Figure 7, left) gives us the total amount of hydrogen mass that entered the onboard tanks.This integration gives us almost 51 g which makes the amount of hydrogen entering each onboard tank during the start-up phase around 17 g, i.e., 9.4% of the initial amount of hydrogen in each tank which is around 180 g.This small amount is sufficient to register an increase in pressure in the onboard tanks.Figure 8 shows the pressure and temperature dynamics in the onboard tanks during the start-up phase of the fuelling.The thermocouples were located in the centre of the onboard tanks in the experiment [19].In the simulations, the pressure and temperature measuring points are located exactly in the same location as in the experiment.The maximum pressure difference between the simulation and the experiment is within 0.1 MPa (Figure 8, left).This is an excellent result bearing in mind that the accuracy of pressure transducers with a maximum pressure of more than 100 MPa is quite high, i.e., 1 MPa [19].As the mass flow rate reached zero at around t = 8 s (Figure 7, left), it is expected that pressure in onboard tanks stays constant (Figure 8, left).In real life, any decrease in pressure after startup is an indication of leaks in the fuelling line, in which case the HRS safety systems prevent the start of the main fuelling stage.The analysis of experimental dynamics of the mass flow rate across the PCV and temperature after the HE allows us to assume that the HE starts to work as the PCV opens to generate the connection pulse and the HE cooling is turned off as the PCV starts to close.This could explain why after the mass flow starts to decrease, the temperature at the HE exit stops the decrease and starts to grow.The integration of mass flow rate over time (Figure 7, left) gives us the total amount of hydrogen mass that entered the onboard tanks.This integration gives us almost 51 g which makes the amount of hydrogen entering each onboard tank during the start-up phase around 17 g, i.e., 9.4% of the initial amount of hydrogen in each tank which is around 180 g.This small amount is sufficient to register an increase in pressure in the onboard tanks.
Figure 8 shows the pressure and temperature dynamics in the onboard tanks during the start-up phase of the fuelling.The thermocouples were located in the centre of the onboard tanks in the experiment [19].In the simulations, the pressure and temperature measuring points are located exactly in the same location as in the experiment.The maximum pressure difference between the simulation and the experiment is within 0.1 MPa (Figure 8, left).This is an excellent result bearing in mind that the accuracy of pressure transducers with a maximum pressure of more than 100 MPa is quite high, i.e., 1 MPa [19].
temperature at the exit from the HE reproduces closely the experimental data with an acceptable engineering accuracy of ±0.5 K.
The analysis of experimental dynamics of the mass flow rate across the PCV and temperature after the HE allows us to assume that the HE starts to work as the PCV opens to generate the connection pulse and the HE cooling is turned off as the PCV starts to close.This could explain why after the mass flow starts to decrease, the temperature at the HE exit stops the decrease and starts to grow.The integration of mass flow rate over time (Figure 7, left) gives us the total amount of hydrogen mass that entered the onboard tanks.This integration gives us almost 51 g which makes the amount of hydrogen entering each onboard tank during the start-up phase around 17 g, i.e., 9.4% of the initial amount of hydrogen in each tank which is around 180 g.This small amount is sufficient to register an increase in pressure in the onboard tanks.Figure 8 shows the pressure and temperature dynamics in the onboard tanks during the start-up phase of the fuelling.The thermocouples were located in the centre of the onboard tanks in the experiment [19].In the simulations, the pressure and temperature measuring points are located exactly in the same location as in the experiment.The maximum pressure difference between the simulation and the experiment is within 0.1 MPa (Figure 8, left).This is an excellent result bearing in mind that the accuracy of pressure transducers with a maximum pressure of more than 100 MPa is quite high, i.e., 1 MPa [19].As the mass flow rate reached zero at around t = 8 s (Figure 7, left), it is expected that pressure in onboard tanks stays constant (Figure 8, left).In real life, any decrease in pressure after startup is an indication of leaks in the fuelling line, in which case the HRS safety systems prevent the start of the main fuelling stage.As the mass flow rate reached zero at around t = 8 s (Figure 7, left), it is expected that pressure in onboard tanks stays constant (Figure 8, left).In real life, any decrease in pressure after startup is an indication of leaks in the fuelling line, in which case the HRS safety systems prevent the start of the main fuelling stage.
The simulations show that the temperature at the centre of the onboard tank (Figure 8, right) reaches its maximum value when the connection pulse is almost finished (see mass flow rate change in time, Figure 7, left).This increase in hydrogen temperature in onboard storage tanks is due to gas compression, the conversion of kinetic energy into internal energy [1], and the lower thermal conductivity of the materials [33].After the termination of the pulse and when the mass flow rate into the tank reaches zero, the temperature inside the onboard tanks decreases due to heat transfer to the colder tank walls, similar to what was observed in other experiments by Kuroki et al. [17].Contrary to this, the experimental temperature transients used for validation Test No.2 of Kuroki et al. [20] show a continuous increase in temperature until the end of the start-up period of 14 s (Figure 8, right).These experimental temperature dynamics are somewhat "unexpected" as hydrogen ingress to the onboard tank, i.e., compression, stops after 8 s of the start-up phase.The experimental temperature continues to increase in the experiment even after 8 s when the mass flow rate into the tanks is zero, and the pressure transient in Figure 8 (left) shows that the pressure inside the tanks remains constant after 8 s from the beginning of the start-up phase.
To understand this unusual behaviour of temperature, let us first assume that the process is simply the result of adiabatic compression of hydrogen in the tank.Using the equation of the adiabatic process, the temperature of hydrogen in onboard tanks can be calculated as: Figure 8 (right) shows by the blue line with diamonds the temperature of hydrogen in the assumption of the adiabatic process (experimental pressure is taken as input, Figure 8, left) in comparison with the experimental and simulated temperature.The maximum temperature is achieved at about 6 s.There is no further temperature increase afterwards in the adiabatic process, and temperature remains constant until the end of the start-up phase.The adiabatic process assumes no heat loss to and through the tank wall, thus the calculated adiabatic temperature in the onboard tank is higher compared to experimental and simulated ones as expected.Now, in an attempt to explain unexpected experimental temperature growth in the onboard tanks after stopping of hydrogen flow, let us assume that the experimental acquisition system uses a smoothing technique similar to the rolling average to eliminate temperature oscillations, or by any other reason.In fact, temperature oscillations at thermocouple locations could be imposed by continuous changing of temperature in this location due to flow turbulence and swirls of the non-uniform temperature inside the tank.The last could "serve" as physical rolling averaging of experimentally measured temperature by an inertial thermocouple.Figure 9 shows a Triple Moving Average (TMA) of the simulated temperature, which is a common practice among experimentalists [34], in two tanks obtained during simulations.The simulation data is reported each 0.5 s and the 10-10-10 TMA method is used for averaging, effectively resulting in averaging over 5 last seconds three times, given the data was recorded every 0.5 s. Figure 9 displays a TMA over the 10 last reported values of the simulated temperature.The deviation of the simulated temperature from the experimental temperature in this case is within the accuracy of the thermocouples used in the experiments, which is +1.5 K [19].The developed CFD model allows us to investigate the underlying fuelling phenomenon in-depth.For example, Figure 10 shows the temperature distribution in the cross-section of one of three onboard tanks.As it was calculated previously, the mass of The developed CFD model allows us to investigate the underlying fuelling phenomenon in-depth.For example, Figure 10 shows the temperature distribution in the cross-section of one of three onboard tanks.As it was calculated previously, the mass of each onboard tank increases during leak check by only 9.4% compared to the initial mass of the tanks.This small amount is not expected to form large temperature non-uniformity inside the tank.Indeed, the simulations demonstrated that during the start-up phase, the maximum difference between the averaged (bulk) temperature inside the onboard tank and the maximum hydrogen temperature is about 1 K (Figure 10).The difference between the minimum and maximum temperatures inside the tank is 2.4 K.This is in line with findings from Ramasamy et al. [25] in which they showed that temperature non-uniformity for tanks with small L/D around 3 is very small.The capability of the CFD model to predict the 3D distribution of temperature can be very useful in cases with tanks with larger L/D and conformable tanks.The developed CFD model allows us to investigate the underlying fuelling phenomenon in-depth.For example, Figure 10 shows the temperature distribution in the cross-section of one of three onboard tanks.As it was calculated previously, the mass of each onboard tank increases during leak check by only 9.4% compared to the initial mass of the tanks.This small amount is not expected to form large temperature non-uniformity inside the tank.Indeed, the simulations demonstrated that during the start-up phase, the maximum difference between the averaged (bulk) temperature inside the onboard tank and the maximum hydrogen temperature is about 1 K (Figure 10).The difference between the minimum and maximum temperatures inside the tank is 2.4 K.This is in line with findings from Ramasamy et al. [25] in which they showed that temperature nonuniformity for tanks with small L/D around 3 is very small.The capability of the CFD model to predict the 3D distribution of temperature can be very useful in cases with tanks with larger L/D and conformable tanks.

Conclusions
The originality of this work is in the development of a CFD model able to reproduce experimentally measured pressure and temperature in the onboard storage tanks during the start-up phase of fuelling through the entire equipment of HRS.The advanced modelling of the performance of the key HRS components, i.e., the PCV and the HE, is carried out using an in-house UDF controlling the mass flow rate in the PCV and temperature at the HE exit.
The significance of the study is in the development of a computationally affordable contemporary CFD tool for the design of fuelling protocols for the whole range of hydrogen-powered vehicles for road, rail, marine, and aeronautical applications.The . Contour plot of temperature in the onboard tank #1 during the start-up phase at t = 6 s.

Conclusions
The originality of this work is in the development of a CFD model able to reproduce experimentally measured pressure and temperature in the onboard storage tanks during the start-up phase of fuelling through the entire equipment of HRS.The advanced modelling of the performance of the key HRS components, i.e., the PCV and the HE, is carried out using an in-house UDF controlling the mass flow rate in the PCV and temperature at the HE exit.
The significance of the study is in the development of a computationally affordable contemporary CFD tool for the design of fuelling protocols for the whole range of hydrogenpowered vehicles for road, rail, marine, and aeronautical applications.The simulation time of 14 s of the start-up phase is about 2 h.Unlike the reduced models which inherently simplify the refuelling process, the CFD model is capable of providing insight into the underlying physics of complex phenomena of heat and mass transfer between hydrogen, components of HRS and the surrounding atmosphere.One of the key advantages of the CFD model over reduced models is the capability to predict temperature non-uniformity in onboard tanks.It was demonstrated that the non-uniformity within onboard tanks during the startup phase is not substantial.However, it is anticipated that the non-uniformity during the main fuelling stage would be significantly elevated.This is essential for the prevention of an onboard tank failure, especially in tanks with large L/D ratios such as those in conformable tanks.
The rigour of this research is defined by the validation of the simulations against the unique experimental data of NREL on pressure and temperature during the start-up phase of the refuelling process through the entire HRS equipment.

Data Availability Statement:
The experimental data presented in this study are publicly available in the cited references.The presented simulation data will be available upon request from readers.

Acknowledgments:
The authors are thankful to Taichi Kuroki from NREL for providing information regarding the experimental details.

Conflicts of Interest:
The authors declare no conflict of interest.
: Test No.1 without leak checks during the main fuelling and Test No.2 with initial and two follow-up leak checks during the main fuelling process.The initial conditions and the results of the Test No.2 start-up phase are used to validate the CFD model developed in this study.Experimentally observed localised temperature drop during leak checks at PCV [20] is a subject of our forthcoming study that will be published in due course.

Hydrogen 2023, 4 589Figure 3 .
Figure 3.The computational domain comprises all elements of HRS in the NREL refu experiment.

Figure 3 .
Figure 3.The computational domain comprises all elements of HRS in the NREL refuelling experiment.

Hydrogen 2023, 4 ,Figure 3 .
Figure 3.The computational domain comprises all elements of HRS in the NREL refue experiment.

Figure 4 .
Figure 4.The discretisation of hydrogen tanks: (top) HP tanks and (bottom) onboard tanks (s are different).

Figure 4 .
Figure 4.The discretisation of hydrogen tanks: (top) HP tanks and (bottom) onboard tanks (scales are different).

Figure 5 .
Figure 5.The transition from one pipe diameter to another is realised as a gradual change.

Figure 6 .
Figure 6.The geometry of modelled bends for (left) 5.1 mm pipes and (right) 7.9 mm pipes.

Figure 5 .
Figure 5.The transition from one pipe diameter to another is realised as a gradual change.

Figure 5 .
Figure 5.The transition from one pipe diameter to another is realised as a gradual change.

Figure 6 .
Figure 6.The geometry of modelled bends for (left) 5.1 mm pipes and (right) 7.9 mm pipes.
m sim , and compares it with the experimental value, .m exp .Based on a comparison of experimental and simulated mass flow rates, the dimensionless parameter dm norm is introduced to calculate the departure of simulated .m sim from the experimental .m exp measured in the onboard tank:

Figure 8 .
Figure 8. (Left): Simulated versus experimental pressure in the onboard tanks.(Right): simulated instantaneous temperature in the tank centre versus experimentally measured temperature in the two onboard tanks, and temperature dynamics in the assumption of the adiabatic process (blue line with diamonds).

Figure 8 .
Figure 8. (Left): Simulated versus experimental pressure in the onboard tanks.(Right): simulated instantaneous temperature in the tank centre versus experimentally measured temperature in the two onboard tanks, and temperature dynamics in the assumption of the adiabatic process (blue line with diamonds).

Figure 8 .
Figure 8. (Left): Simulated versus experimental pressure in the onboard tanks.(Right): simulated instantaneous temperature in the tank centre versus experimentally measured temperature in the two onboard tanks, and temperature dynamics in the assumption of the adiabatic process (blue line with diamonds).

Figure 9 .
Figure 9.Comparison of rolling average simulated temperature against experimentally measured temperature.

Figure 9 .
Figure 9.Comparison of rolling average simulated temperature against experimentally measured temperature.

Figure 10 .
Figure 10.Contour plot of temperature in the onboard tank #1 during the start-up phase at t = 6 s.

Funding:
This research has received funding from the Engineering and Physical Sciences Research Council (EPSRC) of the UK for funding through the EPSRC Centre for Doctoral Training in Sustainable Hydrogen "SusHy" (Grant EP/S023909/1) and the Fuel Cells and Hydrogen 2 Joint Undertaking (now Clean Hydrogen Partnership) under the European Union's Horizon 2020 research and innovation programme through the SH2APED project under grant agreement No. 101007182.Simulations were performed using the Tier 2 High-Performance Computing resources provided by the Northern Ireland High-Performance Computing (NI-HPC) facility, funded by the EPSRC (grant EP/T022175/1, https://www.ni-hpc.ac.uk/Kelvin2/ (accessed on 1 August 2023)).

C
Specific heat [J/kg/K)] S k , S ε User-defined source terms for k [m 2 /s 2 ] and ε [m 3 /s 3 ] C 1ε , C 2ε , C 3ε Standard k-ε model constants t Time [s] C d Discharge coefficient [-] u i , u j , u k Velocity components [m/s]