Modelling of Liquid Hydrogen Boil-Off

Abstract

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 densification data obtained from an LH₂ 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₂, and boil-off from a previously pressurized tank as it is being relieved of vapor.

1. 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₂) [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₂ 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³ cryogenic tank and 0.06% per day for a 20,000 m³ LH₂ 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].

For planned future large-scale LH₂ storage and transport applications involving land-based tanks and sea-borne vessels with capacities from around 50,000 to 160,000 m³ and above, the management of BOG is crucial from both an economic and safety perspective. However, enabling BOG management on such scales represents a significant technical challenge that is yet to be resolved [24,25,26]. Tools that can accurately estimate BOG production are particularly important in the context of industrial-scale LH₂ storage, as the capital expense to design and implement such systems is large, and experimental testing at a relevant pilot-scale can be cost-prohibitive with the results obtained making it difficult to translate to full-scale if their description is entirely empirical [27,28]. A high-fidelity, physically-based BOG model would allow developers to iterate through full-scale vessel layouts completely within the design space, enabling cost and time savings by better informing decisions on the overall development path and strategy [7].

Currently, the space industry is an important user of LH₂, with long duration storage of cryogenic hydrogen being essential in orbital transfer vehicles, the upper stages of rockets, and deep space exploration probes [29,30,31,32,33,34]. Boil-off is particularly important during atmospheric departure, where frictional drag forces result in substantial heat ingress into the containment vessel, thereby increasing vapor temperature and pressure [31]. If the warmer stratified layers enter the rocket propellant feed system, pump cavitation can occur, which could potentially destroy the rocket [35]. Additionally, as this stratified mass is unusable as fuel, it becomes a liability to the payload capacity of the rocket [36]. The impact of boil-off limitations on future long-term Mars missions is described in detail by Schaffer [37], who found that the payload weight of the cryogenic propulsion stages is considerably influenced by the boil-off rate. To reduce the impact of boil-off, designs incorporating multi-layer insulation, minimal tank thermal connection points, and deployable sun shields have been proposed [37]. An accurate and robust BOG model is thus required to model the storage of LH₂ in rockets and spacecraft, and thereby verify the performance of newly designed LH₂ storage systems.

Many computational fluid dynamics (CFD) and thermodynamic boil-off models have been developed for cryogenic liquids such as liquefied natural gas (LNG), including models based on heat ingress and subsequent BOG rates and those accounting for heat transfer between the vapor and liquid phases [38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55]. Migliore et al. [44] developed a non-equilibrium numerical model that allows for the vapor phase to be at a higher temperature than the liquid phase. This was then expanded by Huerta and Vesovic [39] to include more realistic heat transfer mechanisms in the vapor phase. Both models consider heat ingress in the vapor and liquid phases separately and enable heat transfer between the two phases across the interfacial surface area. More recently, a non-equilibrium model was proposed by Wang et al. [46] to simulate LNG during storage and refuelling. However, these models were developed for LNG mixtures in specific scenarios and limited experimental or industrial data were available to verify simulations conducted using these models [56,57]. Furthermore, these models tend to use cubic equations of state for thermodynamic calculations, which are not optimised for the most accurate description of LNG, let alone liquid hydrogen.

We have previously [58,59] presented a boil-off model based on a superheated vapor assumption similar to the non-equilibrium methods presented in the literature. This model is flexible with regards to the fluid type, utilises reference Helmholtz equations of state [60,61] for the thermodynamic property calculations and, most importantly, has physically-based parameters that are anchored by high-quality experimental data for pure fluids measured using a custom BOG apparatus. The model’s performance for LNG was then validated against data measured for binary and multi-component mixtures of light hydrocarbons and N₂. However, to correctly use this model for the description of liquid hydrogen storage and boil-off, it must adequately capture the influence of conversions between the two spin isomers, orthohydrogen. and parahydrogen. In particular, at the cryogenic storage conditions where liquid hydrogen is stored, the equilibrium fraction of parahydrogen is around 99.8%. However, non-equilibrium models need to also describe a vapor phase at a higher temperature, for which the concentration of orthohydrogen may be non-negligible. The temperature dependent equilibrium orthohydrogen–parahydrogen ratio can be calculated using Equations (1) and (2) [62].

$$K=\frac{{{\displaystyle \sum }}_{N=1,3,5,\dots }3\left(2N+1\right)·\exp \left(-\frac{E}{kT}\right)}{{{\displaystyle \sum }}_{N=1,3,5,\dots }3\left(2N+1\right)·\exp \left(-\frac{E}{kT}\right)+{{\displaystyle \sum }}_{N=0,2,4,\dots }\left(2N+1\right)·\exp \left(-\frac{E}{kT}\right)}$$

(1)

$$E=B_{0}N\left(N+1\right)-D_0{N}^2{\left(N+1\right)}^2+H_0{N}^3{\left(N+1\right)}^3$$

(2)

Here, N is the rotational quantum number, k is the Boltzmann constant, E is the energy of the Nth rotational state, B₀ = 7.36 meV is the rotational constant of the hydrogen molecule, and D₀ = 7.69 × 10⁻³ meV and H₀ = 6.45 × 10⁻³ 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% ortho- and 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₂ 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 boil-off losses along the entire LH₂ pathway, from liquefaction to dispensing. Boil-off from the storage tank during the transfer of LH₂ 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₂ fill, with heat fluxes varying from approximately 2.2 W/m² when nearly full, to 1.0 W/m² when empty for a 12.5 m³ 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].

2. 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₂. 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.

Figure 2 and Equations (3) through (13) convey the key concepts and equations underlying the SHV model by reference to a vertically-oriented cylindrical storage tank. The relationships defining heat ingress in each phase, Q_L (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_{relief}$. The boil-off rate, ${\dot{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, ${\dot{n}}_{relief}$, and amount, $n_{relief}$, 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).

$$Q_L=\left(U_L{A}_L\left(T_{boundary}-T_L\right)+Q_{floor}\right)\mathsf{\Delta }t$$

(3)

$$Q_V=\left(U_V{A}_V\left(T_{boundary}-T_V\right)+Q_{roof}\right)\mathsf{\Delta }t$$

(4)

$$Q_{VL}=U_{VL}{A}_{VL}\left(T_V-T_L\right)\mathsf{\Delta }t$$

(5)

$$\frac{d{H}_L}{dt}=Q_L+Q_{VL}-n_{boil}{h}_{boil}$$

(6)

$$\frac{d{H}_V}{dt}=Q_V-Q_{VL}-n_{relief}{h}_V+n_{boil}{h}_{boil}$$

(7)

$${\dot{n}}_{boil}=-\frac{d\left(n_L\right)}{dt}$$

(8)

$$n_{boil}=n_{L,j-1}-n_{L,j}$$

(9)

$${\dot{n}}_{relief}=-\frac{d\left(n_F\right)}{dt}$$

(10)

$$n_{relief}=n_{F,j-1}-n_L-n_V$$

(11)

$$n_F=n_L+n_V$$

(12)

$$L_f=\frac{n_L}{n_F}$$

(13)

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

3. Results and Discussion

3.1. 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, and the BOG calculations of Migliore et al. [44] for a standard LNG storage vertical tank with a 165,000 m³ 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₂ boil-off and tank self-pressurization are compared with corresponding estimations made using BoilFAST to validate its use for liquid hydrogen storage applications.

3.2. LH₂ Refrigeration and Storage System at NASA

BoilFAST was used to simulate LH₂ boil-off experiments conducted in a vacuum-jacketed 125 m³ 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. LH₂ tank dimensions and simulation parameters.

Tank Properties		Fluid Properties
Inner Diameter	2.896 m	Fill Levels	33%, 67%
Insulation Thickness	0.021 m	Ambient Temperature	300 K
Head Height	0.7715 m	UVL	0.6 W/m2/K
Initial Pressure	111 kPa	Refrigeration Rate	Max.: 880 W

. 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:

$$\frac{1}{U}=\frac{A_o}{A_i{h}_i}+\frac{A_o}{A_e}\frac{x}{k_e}+\frac{1}{h_o};A_e=\frac{A_o-A_i}{\ln \left(\frac{A_o}{A_i}\right)}$$

(14)

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⁻¹·K⁻¹; 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⁻²·K⁻¹. In this work, the heat transfer coefficient across the vapor–liquid interface U_VL (0.6 W·m⁻²·K⁻¹) was determined 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 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₂ by cooling below the triple point (13.8033 K [64]) conducted using the IRAS system [78]. The first densification experiment involved 57.5 m³ of LH₂ 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₂ applications in microgravity environments, including the cryogenic storage systems needed for future long-term space missions [90].

3.3. Self Pressurization of Spherical LH₂ 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₂ storage tank subject to three heat fluxes between 0.35 and 3.5 W/m²-, hereafter referred to as Tests 1–3. The experimental setup consisted of an ellipsoidal 4.89 m³ storage tank, modelled here as a sphere that was 83–84% full of LH₂. 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. Summary of the BoilFAST model inputs used to simulate the LH₂ self-pressurization experiments conducted at NASA’s Glenn Research Centre.

Heat Flux (W/m2)	Ambient T. (K)	UL (W/m2/K)	UV (W/m2/K)	UVL (W/m2/K)
0.35	83	0.014	0.014	1.04
2	294	0.0245	0.0245	1.04
3.5	350	0.0351	0.0351	1.04

. 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 vapor–liquid heat transfer coefficient, U_VL = 1.04 W/m²/K, was chosen as this best matched the average vapor temperature behaviour presented by Hasan et al. [53] for the 3.5 W/m² case. In our previous work modelling LNG boil-off [59], we found that a U_VL of 4.0 W/m²/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³ vertical cylinder, while the cases modelled here were conducted in a 4.89 m³ ellipsoidal tank.

The self-pressurization results for Tests 1 to 3 for 0.35 to 3.5 W/m² 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²); 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², 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.

3.4. Boil-Off Losses along the LH₂ Pathway

In the work of Petitpas [27], a 12.5 m³ LH₂ 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₂ was 10.2 m³. 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⁻²·K⁻¹ for both U_V and U_L, a value consistent with the range of U (0.0027–0.0102) W·m⁻²·K⁻¹ that was estimated for similar multi-layer insulation systems [88]. A vapor–liquid heat transfer coefficient, U_VL = 0.15 W/m²/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.

4. 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.