Modelling of Liquid Hydrogen Boil-Off

: A model has been developed and implemented in the software package BoilFAST that allows for reliable calculations of the self-pressurization and boil-off losses for liquid hydrogen in different tank geometries and thermal insulation systems. The model accounts for the heat transfer from the vapor to the liquid phase, incorporates realistic heat transfer mechanisms, and uses reference equations of state to calculate thermodynamic properties. The model is validated by testing against a variety of scenarios using multiple sets of industrially relevant data for liquid hydrogen (LH2), including self-pressurization and densiﬁcation data obtained from an LH 2 storage tank at NASA’s Kennedy Space Centre. The model exhibits excellent agreement with experimental and industrial data across a range of simulated conditions, including zero boil-off in microgravity environments, self-pressurization of a stored mass of LH 2 , and boil-off from a previously pressurized tank as it is being relieved of vapor.


Introduction
As the world's population grows and nations further develop, global energy demand is set to increase substantially. To minimize the associated increase in greenhouse gas emissions, a transition away from fossil fuels is required. Liquid hydrogen (LH 2 ) [1][2][3] will play a key role in this transition [4] as it is clean-burning and can function as a store of energy generated from zero-or low-carbon sources, such as renewables or natural gas combined with carbon capture and storage. However, at ambient temperatures, hydrogen is a gas and thus occupies substantial volume. To overcome this, hydrogen can be liquefied by cooling to approximately 20 K (−253 • C) at atmospheric pressure [5][6][7]. The resulting temperature difference between the volume of stored liquid and the ambient environment ensures that heat ingress into the LH 2 is inevitable, which can cause it to evaporate [8]. The evaporated vapor is termed boil-off gas (BOG) [9][10][11][12][13] and its generation leads to a pressure increase in the relevant storage tank in a process known as self-pressurization, necessitating venting of the tank into the atmosphere and loss of valuable hydrogen. Depending on the insulation quality and the tank's surface-to-volume ratio, BOG generation can be on the order of 0.4% per day for a 50 m 3 cryogenic tank and 0.06% per day for a 20,000 m 3 LH 2 tank [14]. Additionally, heat ingress causes the vapor temperature to increase faster than that of the liquid due to the vapor's higher thermal diffusivity (α = κ/ρ·c p ), resulting in heat conduction across the vapor-liquid interface, and thus a temperature gradient in the top layer of the liquid phase with the interface at a higher temperature than the bulk liquid: this is known as thermal stratification [15][16][17][18][19][20][21][22][23].
Here, N is the rotational quantum number, k is the Boltzmann constant, E is the energy of the Nth rotational state, B 0 = 7.36 meV is the rotational constant of the hydrogen molecule, and D 0 = 7.69 × 10 −3 meV and H 0 = 6.45 × 10 −3 meV are the rotational energy distortion constants [63]. Unfortunately, most implementations of the reference Helmholtz equations of state for parahydrogen and orthohydrogen developed by Leachman et al. [64] require the user to independently evaluate Equations (1) and (2) to determine the ortho-para ratio at each temperature.
At room temperature at equilibrium, so-called normal hydrogen comprises 75% orthoand 25% para-hydrogen molecules (See Equations (1) and (2)). At absolute zero, all the hydrogen molecules must be in a rotational ground state and thus 100% para-hydrogen [65]. For liquid hydrogen at its normal boiling point (20.3 K), the equilibrium para-hydrogen content is 99.8%. Liquefaction processes use heterogeneous catalysts to ensure this equilibrium ratio is attained to prevent the exothermic ortho-to para-hydrogen conversion from occurring during storage [66]. However, because of its larger thermal diffusivity, the gas phase above the stored liquid hydrogen will generally be warmer than the liquid leading to the reverse, endothermic conversion of para-to ortho-hydrogen [67]. In the absence of suitable heterogeneous catalysts, the conversion is very slow and occurs over the order of a few days. Hence, any dynamic model that assumes the hydrogen vapor is at the equilibrium ortho-para ratio corresponding to its temperature might incorrectly represent the system's dynamic energy balance and, consequently, the boil-off rate. Several proposals based on exploiting the endothermic para-ortho conversion for cooling purposes and/or to reduce boil-off [67][68][69][70] have been considered; these efforts would be facilitated by a BOG model able to account for a vapor that is not at compositional or thermal equilibrium with the liquid.
To our knowledge, the only model developed for BOG studies of LH 2 is that reported by Petitpas [7], who modified a MATLAB (MathWorks, Natick, MA, USA) code previously developed by the National Aeronautics and Space Administration (NASA) to estimate boiloff losses along the entire LH 2 pathway, from liquefaction to dispensing. Boil-off from the storage tank during the transfer of LH 2 from a supply trailer to the storage tank was found to be one of the most important sources of hydrogen loss. While many studies [71][72][73][74][75][76][77][78][79][80][81][82][83][84] have investigated how to minimize boil-off losses and empirically determined that top fill (spray) is the most effective way to reduce transfer losses during loading, the underlying physics remains poorly understood. Petitpas used his model to empirically determine the dependence of heat transfer into the storage tank on the level of LH 2 fill, with heat fluxes varying from approximately 2.2 W/m 2 when nearly full, to 1.0 W/m 2 when empty for a 12.5 m 3 vertical storage tank. It was also shown that appreciable temperature gradients (thermal stratification) may exist within the fluid contained by the vessel.
Here, we present an improved version of the superheated vapor (SHV) model developed previously [58,59] that is capable of simulating liquid hydrogen storage. The model is flexible with regards to the spin isomer composition of the hydrogen in each phase, tank geometry, and associated heat transfer, allowing the simulation of a wide range of scenarios. This more rigorous modelling approach is then used to simulate a range of liquid hydrogen storage experiments conducted at NASA and the Lawrence Livermore National Laboratory [27,53,75,78].

Model Description
The freely available software package BoilFAST (UWA Fluid Sciences and Resources Division, Perth, Australia) [85] contains an implementation of the extended SHV model within a Dynamic Link Library (DLL) that is accessed via a Graphical User Interface (GUI) in a Windows environment. This allows the user to easily set up a simulation, selecting from a range of tank geometries, thermodynamic models, boundary conditions, and compositions, and run the simulation using the SHV model, as well as view or export the results. An image of the BoilFAST GUI is shown in Figure 1. The GUI allows the user to simulate boil-off for many industrially important fluids including LH 2 . To this end, a range of fluid components can be selected in the GUI, including both spin isomers of hydrogen and equilibrium and normal hydrogen. The latter is included to facilitate future comparisons of BOG generation in systems that have been cooled too rapidly and have not equilibrated the spin isomer ratio.

Model Description
The freely available software package BoilFAST (UWA Fluid Sciences and Resources Division, Perth, Australia) [85] contains an implementation of the extended SHV model within a Dynamic Link Library (DLL) that is accessed via a Graphical User Interface (GUI) in a Windows environment. This allows the user to easily set up a simulation, selecting from a range of tank geometries, thermodynamic models, boundary conditions, and compositions, and run the simulation using the SHV model, as well as view or export the results. An image of the BoilFAST GUI is shown in Figure 1. The GUI allows the user to simulate boil-off for many industrially important fluids including LH2. To this end, a range of fluid components can be selected in the GUI, including both spin isomers of hydrogen and equilibrium and normal hydrogen. The latter is included to facilitate future comparisons of BOG generation in systems that have been cooled too rapidly and have not equilibrated the spin isomer ratio.   (3)) and QV (Equation (4)), are calculated as a function of the phase contact area with the tank walls, AL and AV; the overall heat transfer coefficients through the tank walls, UL and UV; and the differences between the ambient temperature, Tboundary, and phase temperatures, TL and TV. This represents the combined effect of all the modes of heat transfer between the tank contents and the environment. The overall heat transfer coefficients in Equations (3) and (4) are estimated by considering the heat transfer by convection at the interior and exterior walls, by conduction through the wall, and radiant heat transfer from the inner tank wall to the fluid (using standard correlations given in Ref [45]). These values are largely dependent on the tank geometry, material properties of the tank walls and insulation, fluid properties, and the temperature differential between the phases and ambient conditions. If details of the tank construction are insufficient, one or both may be adjusted to match the experimental data-in particular, the BOG rate. The final heat transfer process considered is from the vapor to the liquid phase, denoted as QVL, and is evaluated using (Equation   (3)) and Q V (Equation (4)), are calculated as a function of the phase contact area with the tank walls, A L and A V ; the overall heat transfer coefficients through the tank walls, U L and U V ; and the differences between the ambient temperature, T boundary , and phase temperatures, T L and T V . This represents the combined effect of all the modes of heat transfer between the tank contents and the environment. The overall heat transfer coefficients in Equations (3) and (4) are estimated by considering the heat transfer by convection at the interior and exterior walls, by conduction through the wall, and radiant heat transfer from the inner tank wall to the fluid (using standard correlations given in Ref. [45]). These values are largely dependent on the tank geometry, material properties of the tank walls and insulation, fluid properties, and the temperature differential between the phases and ambient conditions. If details of the tank construction are insufficient, one or both may be adjusted to match the experimental data-in particular, the BOG rate. The final heat transfer process considered is from the vapor to the liquid phase, denoted as Q VL , and is evaluated using (Equation (5)), which considers the net heat transfer across the interface by taking into account the evaporation and condensation effects. The heat transfer coefficient, U VL , is dependent on a variety of factors, including the tank geometry, fluid properties, and the orientation of the phase interface; it is this parameter that Perez et al. [58] adjusted to calibrate the results observed for liquid nitrogen in their custom BOG apparatus. This then allowed for the reliable estimation of the boil-off rate and self-pressurization to be made for LNG mixtures studied in the same apparatus. The model also allows the user to specify a fixed heat transfer rate though the tank floor (Q floor ) and/or roof (Q roof ) if required, allowing for a better approximation of the industrial conditions-where tanks are often bottom heated to prevent ground freezing. If these terms are not set to zero, the BoilFAST implementation of the SHV model automatically adjusts its calculation of A L and/or A V to exclude the areas of the floor and roof, respectively. Mass transfer is considered at two locations, the liquid (molar) boil-off to the vapor phase, n boil , and vapor relieved through the vent valve, n relie f . The boil-off rate, . n boil , and quantity, n boil , are determined by Equations (8) and (9), respectively, as the change in liquid quantity between time steps. Similarly, the relief rate, . n relie f , and amount, n relie f , in Equations (10) and (11) are estimated by the change in the total molar quantity of the fluid within the tank between time steps. With the heat and mass transfer defined, the model consolidates these variables into the phase energy balances given by Equations (6) and (7). The right-hand sides of Equations (6) and (7) contain the bulk energy balances for the liquid and vapor phases, respectively; each of which is considered to equal the total enthalpy change of the respective phase over a given time step, which can be calculated using an equation of state (EOS).
The initial conditions of the liquid and vapor in the tank are defined based on the system pressure, assuming a saturated liquid in equilibrium with the vapor phase. The reference EOS models of Leachman et al. [64] for both hydrogen spin isomers and normal hydrogen have been implemented together with Equations (1) and (2) for the equilibrium ortho-para ratio. These are used to calculate the necessary thermodynamic properties: enthalpy, density, and saturated liquid temperature. Using these properties, the amounts of liquid, n L , vapor, n V , total amount of fluid within the system, n F , and the fraction of the system contents in the liquid phase, L f , can be calculated via Equations (12) and (13), respectively. The specified time step interval, ∆t, is used to determine the total amount of heat entering the liquid and vapor phases from the boundary and the amount of heat transferred to the liquid phase from the vapor phase, as given by Equations (3)-(5). . .
Assuming a saturated liquid temperature, T L , the liquid energy balance (Equation (6)) iterates until convergence by applying a root-finding algorithm, as shown in Figure 3. This yields the liquid fraction, L f , which in turn gives the amount of liquid, n L . By solving the liquid energy and mass balances, the volumes of liquid and vapor, liquid (x i ) and boil-off (w i ) composition, and the amount of vapor leaving the liquid (n boil ), Equation (9) during one timestep can be found. To find the vapor temperature, the vapor energy balance is iteratively solved (Equation (7)). During this process, the amount of vapor relieved, n relief , vapor composition (y i ), and the thermodynamic properties of the vapor phase are also found. The vapor temperature is taken to be the average temperature for the vapor phase.
During self-pressurization, the calculation procedure of the system pressure includes the phase energy balances, adjustment of pressure using a root-finding algorithm (as shown in Figure 3), and calculating the overall system mole balance as defined in Equation (12). This ensures that n relief is zero (no venting) during the pressure build-up stage (i.e., n F,j−1 -n L -n V = 0 for all time steps).
vapor composition (yi), and the thermodynamic properties of the vapor phase are also found. The vapor temperature is taken to be the average temperature for the vapor phase.
During self-pressurization, the calculation procedure of the system pressure includes the phase energy balances, adjustment of pressure using a root-finding algorithm (as shown in Figure 3), and calculating the overall system mole balance as defined in Equation (12). This ensures that nrelief is zero (no venting) during the pressure build-up stage (i.e., nF,j−1-nL-nV = 0 for all time steps).

Model Testing and Validation
In our related work [58,59], we tested and validated the SHV model against our experimental data for pure liquid nitrogen [58] and LNG-like binary mixtures of methane and ethane [59]. The model calculations were also compared against literature data reported by the Gas Research Institute [56] for six different LNG weathering tests, and those reported by Harper and Powars [57] for a self-pressurization test from 283 kPa to 1585 kPa in an advanced on-vehicle LNG fuel system. The model has also been compared against the self-pressurization calculations reported by Wang et al. [46] for a horizontal storage tank,  [44] for a standard LNG storage vertical tank with a 165,000 m 3 capacity. Good agreement exists between the SHV model and each of the sources from all of the comparisons listed above. Recently [86], the model's performance was compared against BOG data obtained for ternary mixtures of methane, ethane, and nitrogen using our apparatus [58] and an LNG mixture used as rocket fuel conducted at Kennedy Space Centre (NASA) [87]. However, no quantitative comparisons have yet been made for liquid hydrogen applications. Here, the industrial and experimental data of LH 2 boil-off and tank self-pressurization are compared with corresponding estimations made using BoilFAST to validate its use for liquid hydrogen storage applications.

LH 2 Refrigeration and Storage System at NASA
BoilFAST was used to simulate LH 2 boil-off experiments conducted in a vacuumjacketed 125 m 3 storage tank capacity at NASA's Kennedy Space Centre, equipped with an 80-layer multilayer insulation (MLI) system. This tank incorporates an Integrated Refrigeration and Storage (IRAS) system [75], which maintains the liquid in a zero boil-off state through the use of an internal cryogenic heat exchanger connected to an external Brayton cycle refrigerator. The IRAS geometric properties of the horizontal capsule tank with a 2:1 ellipsoidal head [88] and parameters used by the simulation are outlined in Table 1. At the initial tank pressure of 111 kPa, hydrogen has a saturation temperature of 20.6 K and para-hydrogen content of 99.8%. The experiments conducted involved a 'duty-cycle' mode, where the refrigeration system was turned off for approximately 80 h, then turned on for 92 h to remove the accumulated heat leak. Comparisons between the experimental data (symbols) and simulation (solid curves) are shown in Figure 4 for 33% and 67% liquid fill levels. The overall heat transfer coefficients U L and U V were predicted as the sum of two convection terms to account for the heat transfer between the tank wall and adjacent fluid layer, and a conductive term describing heat ingress through the tank wall: where A o is the outer surface area, A e is the heat transfer area, A i is the inner surface area, x is the insulation layer thickness, k e is the effective thermal conductivity of the multilayer insulation, and h i and h o are the internal and external convection coefficients. Fesmire and Johnson presented comprehensive heat flux and effective thermal conductivity data for 26 different MLI systems across the full relevant pressure range (from high vacuum to ambient pressure) [89]. The closest test series to the actual MLI used in the horizontal IRAS tank is test FP3-A 128 (80 layers, 21 mm thick). The mean effective thermal conductivity of this IRAS MLI system was 0.14 mW·m −1 ·K −1 ; this value was used in Equation (14). The internal and external convection coefficients were predicted using an empirical correlation relating the Nusselt and Rayleigh numbers [45]. Using Equation (14), the average values of U L and U V during the simulation were found to be 0.0066 W·m −2 ·K −1 . In this work, the heat transfer coefficient across the vapor-liquid interface U VL (0.6 W·m −2 ·K −1 ) was determined empirically from a fit to the experimental vapor temperature data; this will be resolved in future work. of this IRAS MLI system was 0.14 mW•m −1• K −1 ; this value was used in Equation (14). The internal and external convection coefficients were predicted using an empirical correla tion relating the Nusselt and Rayleigh numbers [45]. Using Equation (14), the average val ues of UL and UV during the simulation were found to be 0.0066 W•m −2 •K −1 . In this work the heat transfer coefficient across the vapor-liquid interface UVL (0.6 W•m −2 •K −1 ) was de termined empirically from a fit to the experimental vapor temperature data; this will be resolved in future work. Refrigeration On of this IRAS MLI system was 0.14 mW•m K ; this value was used in Equation (14). The internal and external convection coefficients were predicted using an empirical correla tion relating the Nusselt and Rayleigh numbers [45]. Using Equation (14), the average val ues of UL and UV during the simulation were found to be 0.0066 W•m −2 •K −1 . In this work the heat transfer coefficient across the vapor-liquid interface UVL (0.6 W•m −2 •K −1 ) was de termined empirically from a fit to the experimental vapor temperature data; this will be resolved in future work.  Overall, the BoilFAST calculations are in good agreement with the experimental dat for both the 33% (Figure 4a) and 67% (Figure 4b) fill experiments. The model represent the pressurization and subsequent depressurization for both experiments, well onl slightly deviating in the pressurization rate during the initial heat leak and at the activa tion of the IRAS. The temperatures measured closest to the vapor-liquid interface in eac case are compared with the model calculations. In Figure 4a, sensors TT3 and TT4 wer located several centimetres below and above the vapour-liquid interface, respectivel [78], with the average of their readings taken to be the best experimental representatio of the interface's temperature. This estimate of the interface temperature is in good agree ment with the model's estimation of the dynamic liquid temperature during refrigeration noting that the model represents the liquid as a homogenous saturated fluid with a tem perature set by the tank pressure (i.e., a uniform liquid with the interface's temperature For Figure 4b, sensor TT10, located at the 61% fill level [78], was closest to the vapor liquid interface (approximately 13 cm below); the measured temperature was lower tha Overall, the BoilFAST calculations are in good agreement with the experimental data for both the 33% (Figure 4a) and 67% (Figure 4b) fill experiments. The model represents the pressurization and subsequent depressurization for both experiments, well only slightly deviating in the pressurization rate during the initial heat leak and at the activation of the IRAS. The temperatures measured closest to the vapor-liquid interface in each case are compared with the model calculations. In Figure 4a, sensors TT3 and TT4 were located several centimetres below and above the vapour-liquid interface, respectively [78], with the average of their readings taken to be the best experimental representation of the interface's temperature. This estimate of the interface temperature is in good agreement with the model's estimation of the dynamic liquid temperature during refrigeration, noting that the model represents the liquid as a homogenous saturated fluid with a temperature set by the tank pressure (i.e., a uniform liquid with the interface's temperature). For Figure 4b, sensor TT10, located at the 61% fill level [78], was closest to the vapor-liquid interface (approximately 13 cm below); the measured temperature was lower than that estimated because the liquid also thermally stratified as discussed by Notardonato et al. [75]. The liquid temperature estimated by the model is the saturation value corresponding to the tank pressure, which will always be higher than temperatures measured by sensors immersed in the bulk liquid phase (below the interface) due to stratification. The SHV model currently implemented in BoilFAST assumes the liquid phase is saturated and thermally homogenous. Nevertheless, the agreement between the current model and the experimental data at two different filling levels is encouraging, particularly given the application of refrigeration, which could be modelled in BoilFAST by using an ambient temperature of 20 K (below that of the tank contents).
Additionally, BoilFAST was used to simulate experiments investigating the densification of stored LH 2 by cooling below the triple point (13.8033 K [64]) conducted using the IRAS system [78]. The first densification experiment involved 57.5 m 3 of LH 2 at atmospheric pressure (46% full) and was subjected to refrigeration for 14 days total, reaching the triple point after roughly 10 days. BoilFAST was used to simulate the period of IRAS refrigeration until the triple point was reached; a comparison of the model's results with the experimental data is shown in Figure 4c. The liquid temperature estimation closely aligns with the experimental results, except for an irregularity in the data occurring in the range of 125 to 250 h, which was caused by temporary issues with the refrigeration system. The resultant pressure decline in the tank is also well represented, with the same deviation observed in the middle of the experiment. This capability of the model to simulate storage at sub-atmospheric pressures is highly relevant to LH 2 applications in microgravity environments, including the cryogenic storage systems needed for future long-term space missions [90].

Self Pressurization of Spherical LH 2 Storage Tank at NASA
An experiment conducted at NASA's Glenn Research Centre (formerly Lewis Research Centre) by Hasan et al. [53] investigated the self-pressurization and thermal stratification of a LH 2 storage tank subject to three heat fluxes between 0.35 and 3.5 W/m 2 -, hereafter referred to as Tests 1-3. The experimental setup consisted of an ellipsoidal 4.89 m 3 storage tank, modelled here as a sphere that was 83-84% full of LH 2 . Liquid nitrogen and electrical resistance heaters were used to cool and warm, respectively, the encompassing cylindrical cryoshroud chamber and achieve the desired heat flux [53]. The tank is insulated with two blankets of MLI, each having seventeen Mylar layers; however, the thickness of the insulation layer was not stated by Hasan et al. [53]. The heat transfer parameters used when simulating the three cases are given in Table 2. The ambient temperatures were sourced from Hasan et al. [53]. The ambient-vapor and ambient-liquid heat transfer coefficients were simply estimated based on the measured experimental heat fluxes, the prevailing ambient (boundary) temperature, and the initial system temperature response. It was not possible to use existing correlations to estimate these values in this case as there are insufficient details concerning the insulation system used and uncertainly as to how the use of the electrical resistance heaters can be incorporated into such calculations. A vaporliquid heat transfer coefficient, U VL = 1.04 W/m 2 /K, was chosen as this best matched the average vapor temperature behaviour presented by Hasan et al. [53] for the 3.5 W/m 2 case. In our previous work modelling LNG boil-off [59], we found that a U VL of 4.0 W/m 2 /K produced the best match for the experiment. However, this value depends in part on the tank, fluid, and vapor-liquid interface geometry-the LNG experiments were conducted in a 0.0067 m 3 vertical cylinder, while the cases modelled here were conducted in a 4.89 m 3 ellipsoidal tank. The self-pressurization results for Tests 1 to 3 for 0.35 to 3.5 W/m 2 are shown in Figure 5a and are in excellent agreement with the BoilFAST estimates denoted by the solid curves. Some deviation exists with the data's curvature in Test 2 (2 W/m 2 ); however, this potentially reflects that the heat flux at the beginning of the experiment was initially slightly larger than the average value. Stratification in both the liquid and vapor phases is observed for all tests, particularly for the highest heat flux of 3.5 W/m 2 , as shown in Figure 5b. The temperatures recorded by the Resistance Temperature Detector (RTD) sensors, which measure temperature using the electrical resistance of the sensor material at varying heights from the bottom of the tank, clearly depict this phenomenon-where sensors above a height of 140 cm are stated to be in the ullage region (vapor space) [53]. As a result, the sensor at a height of 139 cm is close to the liquid-vapor interface and is well represented by the model's liquid temperature estimation. The tank has an internal height of 183 cm, thus the sensor located at a height of 162 cm is located approximately in the centre of the vapor space. The temperatures recorded by this sensor are in excellent agreement with the vapor temperature calculated by BoilFAST.

Boil-Off Losses along the LH2 Pathway
In the work of Petitpas [27], a 12.5 m 3 LH2 Dewar tank with a 2 m internal diameter at the Lawrence Livermore Cryogenic Hydrogen Test Facility was subjected to ambient heat ingress to measure the change in liquid volume over summer versus winter. For this experiment, the relief pressure was set to 310 kPa and the initial fill volume of LH2 was

Boil-Off Losses along the LH 2 Pathway
In the work of Petitpas [27], a 12.5 m 3 LH 2 Dewar tank with a 2 m internal diameter at the Lawrence Livermore Cryogenic Hydrogen Test Facility was subjected to ambient heat ingress to measure the change in liquid volume over summer versus winter. For this experiment, the relief pressure was set to 310 kPa and the initial fill volume of LH 2 was 10.2 m 3 . Throughout the experiment, the system remained 'static', meaning no mass transfer across the system boundary (such as refuelling) occurred. Petitpas also presented a model to estimate the boil-off losses based on their experimental data. This model ignores the insulation in the energy balance, instead assuming a constant rate of heat ingress into the tank; 50 W was found to yield the closest estimate to the experimental data. Here, instead of using a constant heat ingress rate of 50 W in BoilFAST, this rate was used to estimate the overall heat transfer coefficients between each phase and the ambient environment of 0.0045 W·m −2 ·K −1 for both U V and U L , a value consistent with the range of U (0.0027-0.0102) W·m −2 ·K −1 that was estimated for similar multi-layer insulation systems [88]. A vapor-liquid heat transfer coefficient, U VL = 0.15 W/m 2 /K was chosen as this best matched the vapor temperature data reported by Petitpas [27]. The experimental results are shown in Figure 6, alongside estimates made using BoilFAST and the model of Petitpas [27]. The resulting simulation matched the experimental data well, capturing the variable heat ingress rate evident in the slightly non-linear experimental data.

Summary
A robust, efficient, and flexible simulation tool, BoilFAST, has been developed to es timate the boil-off from liquid hydrogen storage using a non-equilibrium approach. The results presented in this work show that the results from the BoilFAST simulations are in good agreement with the experimental data from multiple large-scale liquid hydrogen storage tests measured across a range of conditions. The SHV model implemented in Boil FAST has the potential to become a valuable tool for the hydrogen industry through the reliable estimation of self-pressurization and boil-off rates for a variety of scenarios and

Summary
A robust, efficient, and flexible simulation tool, BoilFAST, has been developed to estimate the boil-off from liquid hydrogen storage using a non-equilibrium approach. The results presented in this work show that the results from the BoilFAST simulations are in good agreement with the experimental data from multiple large-scale liquid hydrogen storage tests measured across a range of conditions. The SHV model implemented in BoilFAST has the potential to become a valuable tool for the hydrogen industry through the reliable estimation of self-pressurization and boil-off rates for a variety of scenarios and customisable tank geometries. This will enable for better design and simulation of the liquid hydrogen storage infrastructure-in particular, BOG management systems such as compressors-and by informing tank insulation requirements. This work also demonstrates that the existing engineering correlations for manually estimating wall-fluid heat transfer coefficients (U L and U V ) are sufficiently accurate to allow a robust representation of BOG phenomena across a range of tank sizes, geometries, and fluids-including liquid hydrogen, which provided adequate detail regarding the insulation used. This tool will enable the storage of liquid hydrogen to be further optimized in both the energy and space industries. Future improvements to BoilFAST will focus on calculating the vapour-liquid heat transfer coefficients for different tank geometries and insulation systems, describing any thermal stratification of the liquid phase, and being able to simulate hydrogen storage filling and draining operations.