Sequential Modeling of Liquid Hydrogen Leakage, Dispersion and Post-Ignition Thermal Hazards in Semi-Open Railway Stations

1 Introduction

Hydrogen-powered transportation is an important pathway for deep decarbonization in the transport sector and has attracted increasing attention in rail transit systems. Liquid hydrogen (LH2), with its high volumetric energy density, high storage and transport efficiency, and rapid refueling capability, is considered a promising energy-storage form for future medium- and long-distance hydrogen fuel-cell locomotives. It can greatly improve driving range within limited onboard space and can serve as a clean-energy alternative for non-electrified railway lines [1,2,3]. However, during external refueling, parking, or maintenance, a rupture-induced leakage in the LH2 storage and supply system may rapidly release ultra-low-temperature LH2 onto the ground outside the station, forming a cryogenic liquid pool. The pool then evaporates continuously under ground heat transfer, air convection, and radiation, forming a low-temperature H2 cloud. This process involves coupled physical phenomena, including pool spreading, phase-change evaporation, buoyancy-driven migration, turbulent mixing, and flammable-cloud evolution. As a result, accident consequences show strong transient behavior and spatial non-uniformity [4,5,6,7].

In addition to cryogenic hazards, H2 has a wide flammability range, low ignition energy, and high diffusivity, which leads to a high fire and explosion risk after mixing with air [8]. In semi-open spaces such as stations, the early H2 cloud is mainly controlled by inertia and ambient wind and disperses downwind. As its density decreases, buoyancy gradually becomes dominant, and local retention zones may form beneath the roof, around the locomotive, and in the leeward regions of structures. When the concentration enters the flammable range and an ignition source is present, rapid combustion or even deflagration may occur. Station roofs, parked locomotives, support columns, and other obstacles further alter the flow field, making hydrogen-cloud migration paths and hazardous-zone distributions more complex [9,10].

Current studies on hydrogen and LH2 leakage, dispersion, and combustion mainly use experimental methods and numerical simulations, covering both open and confined-space scenarios. In experiments, high-pressure hydrogen release tests have been used to analyze the effects of leak diameter, pressure, and ambient wind on dispersion behavior [11,12,13]. For example, Niu et al. [14] found that concentration decay in hydrogen jets shows a degree of self-similarity, whereas turbulence induced by ambient wind and obstacles can significantly change far-field concentration distributions and cause local accumulation. However, due to safety and scale limitations, full-scale LH2 leakage experiments remain limited, and equivalent studies often rely on substitute gases or scaled experiments [15]. In numerical simulations, CFD and commercial software such as PHAST have been widely applied to open-space and hydrogen-refueling-station scenarios to simulate hydrogen leakage, dispersion, and combustion [16,17,18,19]. Zhang et al. [20] used CFD coupled with a k-ε turbulence model and a combustion model to study hydrogen leakage and combustion in fuel-cell bus maintenance workshops, showing that leakage location, leak diameter, and ventilation wind speed significantly affect flammable-cloud accumulation, jet-fire behavior, heat-release duration, and workshop safety design. These studies usually couple multiphase flow, turbulent transport, and combustion models to predict hydrogen-cloud concentration fields, temperature fields, and the evolution of flammable regions. For example, in LH2 refueling station scenarios, barrier-structure parameters can markedly influence downwind expansion and vertical migration of hydrogen clouds, and proper design can effectively reduce the flammable-cloud range and cryogenic impact zone [11,21,22]. However, these studies mainly focus on open or semi-open industrial facilities, whose boundary conditions differ clearly from those of station spaces.

For confined spaces, existing studies mainly focus on tunnels, underground garages, and enclosed industrial buildings. The results show that spatial scale, ventilation mode, and obstacle distribution significantly affect upper-layer accumulation, recirculation structures, and combustion intensity of hydrogen clouds [23,24,25,26]. Xu et al. [27] conducted experimental and numerical studies on hydrogen leakage and dispersion in large-aspect-ratio spaces and identified three typical dispersion stages, including buoyancy-dominated rise, horizontal spreading, and vertical filling, showing that leakage rate and direction significantly affect dispersion time and concentration distribution and providing a basis for warning-threshold and emergency-strategy design. Li et al. [28] analyzed leakage accidents of hydrogen-powered vehicles in parking garages and showed that leakage amount, ventilation-opening arrangement, exhaust-fan operation, and partition-wall layout jointly affect hydrogen dispersion and retention risks. Their findings provide a reference for hydrogen-safety ventilation design and risk control in parking garages. Chen et al. [29] analyzed high-pressure hydrogen leakage and dispersion from hydrogen fuel-cell trains in tunnels and found that leakage direction, hole diameter, and train speed strongly affect hydrogen accumulation risk. When the train speed is not lower than 40 km/h, the piston-wind effect can effectively suppress hydrogen accumulation between the train roof and tunnel ceiling. Overall, recirculating airflow, roof confinement, and obstacle-induced flow in confined or semi-confined spaces can enhance local turbulent mixing, change the spatial distribution of flammable clouds, and promote flame acceleration. However, rail transit stations usually have semi-open boundaries, large-span roofs, support columns, and parked locomotives. Their ventilation organization, buoyancy-driven migration paths, and hydrogen-cloud retention mechanisms differ from those of tunnels, underground garages, and enclosed buildings. Existing studies are therefore insufficient to reveal the coupling mechanisms among liquid-pool evaporation, hydrogen-cloud migration, obstacle disturbance, and flame propagation after LH2 leakage in station scenarios.

Existing studies on LH2 leakage have mainly considered open environments, hydrogen refueling stations, and tunnels, whereas semi-open large spaces such as railway stations have received relatively limited attention. In addition, dispersion and combustion are often investigated separately, and the combined effects of liquid-pool evaporation, low-temperature H2 cloud dispersion, station structures, flammable-cloud distribution, and post-ignition flame development have not been extensively examined within an integrated station scenario. In railway stations, roofs, locomotives, and support columns may influence the upward and lateral transport of released H2 and contribute to spatial variations in local hydrogen accumulation. Therefore, this study considers a representative LH2-release scenario in a semi-open railway station under prescribed release, ventilation, and geometric conditions. The analysis focuses on the temporal and spatial evolution of the H2 cloud from the release region to the upper station space, the local influence of the locomotive and support columns on hydrogen distribution, and the development of dispersion-related and post-ignition thermal-hazard regions. The results provide a scenario-specific evaluation of LH2 leakage and ignition behavior under the selected station configuration and operating conditions, rather than a generalized parametric description of LH2 release characteristics.

2 Numerical Model and Methodology

2.1 Physical Properties and Theoretical Model

During operation, parking, or refueling of a liquid-hydrogen fuel-cell locomotive, an integrity failure of cryogenic storage tanks, connecting pipelines, valves, or safety accessories may cause LH2 to be rapidly released onto the ground and form a cryogenic liquid pool. The catastrophic rupture considered in this study refers to abnormal loss of integrity in the hydrogen storage system, causing a large amount of LH2 to be released in an uncontrolled manner. It is not a controllable release such as a small-hole leakage or safety-valve discharge. Although such accidents have a relatively low probability, they are typical high-consequence events. After the cryogenic liquid is released, it undergoes liquid-pool expansion, intense evaporation, H2-air mixing, flammable-cloud formation, and ignition combustion, which may cause cryogenic injury, fire, and local high-temperature impact. In this study, a 6 kg LH2 release is selected as the representative accident condition. This release amount is determined with reference to reported hydrogen inventory requirements for hydrogen-powered vehicles; for example, approximately 5.6 kg of LH2 can support a hydrogen-powered passenger vehicle over a driving range of about 350 miles. Considering that the onboard hydrogen storage capacity of a fuel-cell locomotive is expected to be larger than that of a passenger vehicle, the 6 kg release scenario is used as a conservative representative case for evaluating leakage, dispersion, and combustion hazards of the locomotive LH2 storage system in a semi-open station space [30].

The physical properties of LH2, H2, and air are the basis for source-term calculation and CFD simulation. In the source-term integral model, LH2 is treated as a cryogenic liquid medium. Its saturation temperature, liquid density, latent heat of vaporization, constant-pressure specific heat, and thermal conductivity are used to describe pool formation and heat-induced evaporation after leakage. The phase-change temperature of LH2 is set to 20.3 K to determine the inlet temperature after evaporation into the gas phase. In the CFD model, gaseous H2 and air are treated as miscible gas species. The density of the H2-air mixture is calculated using the ideal-gas equation of state and updated with local temperature and species mass fraction. Transport properties, including viscosity, specific heat capacity, thermal conductivity, and diffusion coefficient, are assigned according to material properties or temperature-dependent parameters. This treatment can reflect the warming, buoyancy variation, and mixing dispersion of low-temperature H2 after it enters the environment. The properties of H2 are listed in Table 1 [31].

This study couples an integral model with CFD simulation to describe the full-process evolution after LH2 leakage. First, the source-term integral model is used to calculate liquid-pool formation, expansion, and evaporation after LH2 is released onto the ground, yielding the time-varying liquid-pool radius and evaporation mass flow rate. The transient source terms from the integral model are then converted into a user-defined function (UDF) and loaded into the CFD model as the H2 inlet boundary condition. The inlet area is determined by the liquid-pool radius, and the inlet velocity is determined jointly by the evaporation mass flow rate and inlet area. Thus, the mass input in the CFD model can reflect the dynamic variations in liquid-pool expansion and evaporation intensity. In the dispersion stage, the Realizable k-ε turbulence model and species transport model are used to solve the flow, dispersion, and dilution of the H2-air mixture in the station space. The analysis focuses on the effects of low-temperature buoyancy, natural wind, and obstacles such as the locomotive, roof, and support columns on hydrogen-cloud migration and local retention.

Fig. 1 presents the modeling framework for evaporation, dispersion, and combustion following LH2 leakage. The transient H2 concentration, temperature, and velocity fields obtained from the dispersion simulation are transferred to the combustion model as the initial conditions for the ignition case. For the combustion simulation, a detailed H2 chemical reaction mechanism is imported through Chemkin, while a turbulent combustion model and a radiation heat-transfer model are employed to predict flame propagation and the evolution of the temperature field after ignition of the flammable H2-air cloud. Through the coupling of the source-term integral model, CFD dispersion model, and combustion model, the process of LH2 leakage, liquid-pool evaporation, H2 dispersion, ignition, and subsequent thermal effects is numerically represented, providing a modeling basis for hazardous-area identification, monitoring-point arrangement, and emergency protection analysis in railway-station scenarios. A single time coordinate, t, measured from the onset of LH2 leakage, is used throughout the study. At t = 3.0 s, the analysis is divided into two scenario branches. In the non-ignited dispersion branch, the CFD calculation continues without ignition to characterize the subsequent migration and dilution of the H2 cloud. In the ignition/combustion branch, ignition is imposed at t = 3.0 s at the location (−0.7, −3.6, 0.574) m, where the local H2-air mixture satisfies the flammability criterion. The instantaneous H2 concentration, temperature, and velocity fields at t = 3.0 s are used as the initial flow-field conditions for the combustion calculation. The same leakage-based time coordinate is retained after ignition; the time is not reset for the combustion calculation.

Figure 1: LH₂ leakage and dispersion process.

2.2 LH₂ Leakage and Liquid-Pool Source-Term Model

At present, source strength is usually calculated based on the liquid-pool size and evaporation rate after leakage, but delayed phase change and related phenomena are not fully considered. In CFD simulations, the atmospheric boundary layer is usually implemented through specific boundary conditions to reproduce dynamic liquid-pool expansion and subsequent evaporation when coupled with an integral model [32]. The leakage and dispersion modules are used to analyze and calculate material release, evaporation, and dispersion, covering various leakage scenarios, including the formation and expansion of an LH2 liquid pool.

When LH2 is released through a rupture opening, the leakage mass flow rate can be expressed as: m•leak=CdAo2ρl[Pt−Pa+ρlg(ht−ho)](1) where m • l e a k is the leakage mass flow rate, kg/s; ρl is the liquid density; Cd is the discharge coefficient; Ao is the leak-orifice area, m2; Pt is the internal pressure of the storage tank, Pa; Pa is the ambient pressure, Pa; g is the gravitational acceleration, m/s2; and ht − ho is the height difference between the liquid level and the leak orifice, m.

After leaking onto the ground, LH2 forms a cryogenic liquid pool. Assuming that the liquid pool spreads approximately in a circular shape on the ground, the mass conservation equation of the liquid pool can be expressed as: dMpdt=m•leak−m•evap(2) Ap(t)=πR2(t)(3) Mp(t)=ρlAp(t)δp(4) R(t)=[Mp(t)πρδp]12(5) where Mp is the mass of LH2 in the liquid pool, kg; m • e v a p is the evaporation mass flow rate of the liquid pool, kg/s; Ap is the liquid-pool area, m2; R is the liquid-pool radius, m; and δp is the equivalent thickness of the liquid pool.

The evaporation of the liquid pool is jointly determined by heat conduction from the ground, convective heat transfer from the air, and environmental radiation. The evaporation mass flow rate is expressed as: m•evap(t)=Ap(t)Lv[qg″(t)+hc(T∞−Tb)+εpσ(Ts4−Tb4)](6) where the transient conductive heat flux from the ground can be written as: qg″(t)=kg(Tg,0−Tb)παg(t−t0+Δt)(7) where Lv is the latent heat of vaporization of LH2, J/kg; q″g is the conductive heat flux transferred from the ground to the liquid pool, W/m2; hc is the convective heat-transfer coefficient, W/(m2·K); T∞ is the ambient air temperature, K; Tβ is the boiling point of LH2, K, which is taken as 20.3 K in this study; εp is the surface emissivity of the liquid pool; σ is the Stefan-Boltzmann constant, W/(m2·K4); Ts is the temperature of the surrounding walls or the environmental radiation temperature, K; kg is the thermal conductivity of the ground, W/(m·K); Tg,0 is the initial ground temperature, K; αg is the thermal diffusivity of the ground, m2/s; t0 is the time when the liquid pool begins to form, s; and Δt is a small time correction used to avoid the singularity of the initial heat flux, s.

In the CFD model, the hydrogen generated by liquid-pool evaporation is treated as a time-varying velocity inlet, and the inlet velocity is expressed as: uin(t)=m•evap(t)ρH2(Tin)Ap(t)(8) where uin is the normal inlet velocity of hydrogen, m/s; ρH2(Tin) is the hydrogen density at the inlet temperature, kg/m3; and Tin is the inlet hydrogen temperature, K.

2.3 Governing Equations and Dispersion Model for a Semi-Open Space with Obstacles

Computational fluid dynamics was used to simulate the dispersion behavior of H2 in the station area. The gas-phase flow is governed by the conservation equations of mass, momentum, and energy, while the mixing and diffusion of H2 and air are described by the species transport equation. Because the hydrogen cloud generated after leakage is affected simultaneously by buoyancy, natural wind, and station structures, its dispersion process involves turbulent mixing and local recirculation. Therefore, the Realizable k-ε model was adopted for turbulence closure in this study. Compared with the standard k-ε model, this model is more suitable for predicting complex shear flows, obstacle-induced flow fields, and recirculation regions, and can describe the migration, dilution, and local retention of H2 clouds in semi-open station spaces [33]. The governing equation is as follows:

Continuity Equation: ∂ρ∂t+∇·(ρv)=Sm(9) where ρ is the density of the gas mixture, kg/m3; u is the velocity vector, m/s; and Sm is the mass source term. When the evaporation source term is imposed through a velocity inlet, Sm = 0 inside the computational domain.

The momentum equation is: ∂(ρv)∂t+∇·(ρvv)=−∇p+∇·τeff+ρg(10) τeff=μeff∇v+(∇v)T−23(∇·v)I(11) where p is the static pressure, Pa; τeff is the effective stress tensor, Pa; μeff is the effective dynamic viscosity, Pa·s; I is the unit tensor; and the superscript T denotes matrix transpose.

The energy equation is: ∂(ρE)∂t+∇·v(ρE+p)=∇keff∇T−∑ihiJi+τeff·v+Sh+Srad(12) where E is the total energy of the gas mixture, J/kg; keff is the effective thermal conductivity, W/(m·K); T is temperature; hi is the specific enthalpy of species i, J/kg; Ji is the diffusion flux of species i, kg/(m2·s); Sh is the volumetric heat source generated by chemical reactions or other processes. During the combustion stage, Sh = ωH2ΔHc, W/m3; and Srad is the radiation heat-transfer source term, W/m3.

The species transport equation is: ∂(ρYi)∂t+∇·(ρvYi)=−∇·Ji+Ri+Si(13) Ji=−ρDeff,i∇Yi(14) Deff,i=Di+μtρSct(15) where Yi is the mass fraction of species i; Ri is the species production rate caused by chemical reactions, kg/(m3·s), and Ri = 0 in the non-reacting dispersion stage; Si is the user-defined species source term, kg/(m3·s). In this study, the hydrogen release was introduced through a transient velocity inlet boundary condition determined by the LH2 pool evaporation model using a user-defined function (UDF). Therefore, no additional volumetric species source was applied inside the computational domain, and Si was set to zero during the dispersion simulation; Deff,i is the effective diffusion coefficient, m2/s; Di is the molecular diffusion coefficient, m2/s; and Sct is the turbulent Schmidt number.

The density of the gas mixture is calculated using the ideal-gas equation of state: ρ=pRuT∑i(YiMi)(16) where ρ is the density of the gas mixture, kg/m3; p is the static pressure, Pa; Ru is the universal gas contant, J/(mol·K); Mi is the malor mass of species i.

The transport equations for turbulent kinetic energy k and turbulent dissipation rate ε are expressed as: ∂(ρk)∂t+∂∂xj(ρkuj)=∂∂xj(μ+μtσk)∂k∂xj+Gk+Gb−ρε+YM+Sk(17) ∂(ρε)∂t+∂∂xj(ρεuj)=∂∂xj(μ+μtσε)∂ε∂xj+ρC1Sε−ρC2ε2k+vε+C1εεkC3εGb+Sε(18) where Gk is the turbulent kinetic energy generated by the mean velocity gradient; YM denotes the contribution of fluctuating dilatation in compressible turbulence to the total dissipation rate; C1, C2 andC1ε are constants, where C1ε = 1.44 and C2 = 1.9 are empirical constants of the standard k-ε turbulence model; and μt is the turbulent viscosity.

2.4 Chemical Reaction Model

During ignition and combustion of the hydrogen cloud, the combustion process was simulated using the partially premixed combustion model in ANSYS Fluent. For the station ignition case, the H2 concentration, temperature, and velocity fields obtained from the dispersion simulation at t = 3.0 s are used as the initial flow-field conditions for the combustion calculation, and the leakage-based time coordinate is retained without resetting after ignition. The partially premixed model describes the flame evolution based on the mixture fraction and reaction progress variable, while the chemical reaction process is represented using a detailed H2-air kinetic mechanism consisting of 10 species and 19 elementary reactions. This approach enables the prediction of flame propagation and thermal hazard distribution following ignition of the hydrogen-air cloud.

The overall hydrogen combustion reaction is:

In the partially premixed combustion model, the flame propagation rate is controlled by the turbulent flame speed ST, which is expressed as: ST=SL[1+Ct(u′SL)n](20) u′=2k3(21) where ST is the turbulent flame speed, m/s; SL is the laminar flame speed, m/s; u′ is the turbulent fluctuation velocity, m/s; k is the turbulent kinetic energy, m2/s2; and Ct and n are model constants. The laminar flame speed SL is determined by the local mixture composition, temperature, and pressure.

The heat released by combustion is treated as the chemical-reaction heat source term in the energy equation: Schem=ω•H2ΔHc(22) where Schem is the chemical-reaction heat source term, W/m3; ωH2 is the hydrogen consumption rate due to reaction, kg/(m3·s); and ΔHc is the heat of combustion of hydrogen, J/kg.

The imported H2 kinetic mechanism is used to determine the composition-dependent combustion properties required by the partially premixed formulation, including the laminar flame characteristics used in the turbulent flame-speed closure.

2.5 Solver Selection

The numerical simulations of LH2 evaporation and dispersion were performed using ANSYS Fluent. A three-dimensional transient pressure-based solver was adopted to capture the temporal evolution of hydrogen dispersion after LH2 leakage. The species transport model was applied to describe the mixing and diffusion process between hydrogen and air, while the turbulence effect induced by buoyancy, natural ventilation, and structural obstacles was modeled using the standard k-ε turbulence model [34]. The pressure–velocity coupling was achieved using the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm [35]. The momentum, energy, species transport, and turbulence equations were discretized using the second-order upwind scheme, while pressure interpolation was performed using the second-order scheme. The transient terms were discretized using the second-order implicit scheme. The time-step size was determined based on the characteristic flow velocity, mesh resolution, and numerical stability. The Courant number and residual convergence were monitored during the calculation, and a time-step sensitivity analysis was conducted to ensure the reliability of transient results.

2.6 Radiation Model

Fluent provides several radiation models, including DTRM, P-1, Rosseland, S2S, and DO, which can account for wall heating or cooling caused by radiation. Among them, the DTRM model has a relatively simple form, and increasing the number of rays can improve calculation accuracy. The P-1 model is suitable for combustion problems with large optical thickness. The Rosseland model has high computational efficiency and low memory demand. The S2S model is mainly used for surface radiative heat transfer in enclosed spaces without participating media. The DO model has a broad range of applicability and can handle radiation problems under different optical thicknesses [36]. Therefore, the DO radiation model is used in this study to predict the influence of radiative heat transfer.

3 Model Development

3.1 LH₂ Release Model

To systematically evaluate the influence of source-term variation after LH2 leakage on the H2 dispersion scenario, this study uses an integral model to simulate liquid-pool expansion and subsequent evaporation after LH2 is released onto the ground. The liquid-pool radius R(t) and mass flow rate m(t) calculated by the integral model are then converted into transient inlet boundary conditions in the CFD model. The time coordinate shown in Fig. 2 is used as the time reference for the dispersion stage. It should be noted that the liquid-pool radius and mass flow rate in Fig. 2 are both close to 0 during t = 0–2.8 s, indicating that no effective liquid-pool source term has formed during this stage. From t ≈ 2.8 s, the mass flow rate and liquid-pool radius begin to increase significantly, indicating that LH2 enters the stage of ground-pool formation and evaporative release.

Figure 2: Time variation of evaporation rate and radius of the H₂ pool source term (a) The radius changes over time; (b) Evaporation rate varies with time.

After leakage, LH2 is released onto the ground and forms a circular cryogenic liquid pool whose radius varies with time. Bund constraints are not considered in the calculation. Because the temperature of LH2 is much lower than the ambient temperature, continuous heat exchange occurs between the liquid pool and the ground, air, and surrounding solid walls. Liquid hydrogen therefore evaporates continuously and forms an H2 cloud. As shown in Fig. 2, the LH2 source-term evolution has clear unsteady characteristics. During t ≈ 2.8–3.7 s, the liquid pool expands rapidly and the mass flow rate increases quickly. At t ≈ 3.7–3.8 s, the mass flow rate reaches its peak, approximately 3.6 kg/s, indicating that the mass of H2 released into the gas phase per unit time is the largest. The mass flow rate then decreases, while the liquid-pool radius continues to increase and reaches its maximum value of about 1.2 m at t ≈ 4.2 s. This indicates that the maximum liquid-pool expansion lags behind the peak mass flow rate.

After t ≈ 4.2 s, liquid supply weakens and evaporation consumption gradually becomes dominant. The liquid-pool radius begins to shrink, and the evaporation mass flow rate decays rapidly. By t ≈ 7.5–7.6 s, the evaporation mass flow rate has decreased to a negligible level, indicating that the active evaporation source has become very weak. A small residual liquid pool may persist for a short period; by approximately t ≈ 8.5 s, the residual pool nearly disappears and the inlet velocity approaches zero, marking the complete cessation of the evaporation source. Thereafter, no new H2 is introduced into the CFD domain, while the existing H2 cloud continues to migrate, dilute, and disperse under the effects of buoyancy, natural wind, and disturbance from station obstacles. Accordingly, the main source-term evolution can be divided into four stages: t = 0–2.8 s, no effective source term; t ≈ 2.8–3.8 s, rapid release and intense evaporation; t ≈ 3.8–4.2 s, decreasing mass flow rate with continued liquid-pool expansion; and t ≈ 4.2–7.6 s, liquid-pool shrinkage and rapid source-term decay, followed by a short residual-pool disappearance period up to approximately 8.5 s.

In the CFD model, H2 generated by liquid-pool evaporation enters the computational domain through a velocity-inlet boundary. Because Fluent cannot easily describe a time-varying liquid-pool boundary directly, a user-defined function (UDF) is used to define the source term obtained from the integral model as: v(t)=1ρH2A(t)dM(t)dt(23) where ρH2 is the density of gaseous H2 under the inlet condition, kg/m3; and A(t) is the gaseous H2 inlet area, m2. During evaporation, the inlet area is related to the liquid-pool radius R(t) calculated by the integral model, A(t) = πR(t)2.

Thus, the CFD inlet boundary reflects both the change in liquid-pool area and the change in source strength. Specifically, the inlet source term starts at t ≈ 2.8 s, reaches the maximum mass flow rate at t ≈ 3.7–3.8 s, corresponds to the maximum liquid-pool area at t ≈ 4.2 s, and decays to nearly 0 at t ≈ 7.5–7.6 s. This treatment avoids the error caused by assuming a constant liquid-pool radius and more reasonably describes the early strong source release and later source-term decay after LH2 leakage.

3.2 SValidation of the LH₂ Leakage Integral Model

To evaluate the predictive capability of the source-term integral model for LH2 release, pool evaporation, and subsequent gaseous hydrogen cloud evolution, publicly reported NASA large-scale LH2 release experimental data were used for comparison, as shown in Fig. 3 [37]. In the experiment, gas samples were collected using sampling bottles, and hydrogen concentrations were also recorded by hydrogen sensors. Test 6 was selected for validation because, under its wind-field conditions, the hydrogen cloud passed through a relatively large number of monitoring locations, providing sufficient concentration-distribution data for assessing the predicted flammable-cloud envelope. The validation focused on the 4 vol% H2 concentration envelope. The experimental results show that, at 20.94 s after LH2 release, the 4 vol% flammable hydrogen cloud reached a height of approximately 20–25 m and a lateral spreading distance of approximately 37–40 m. The integral model predicted a corresponding cloud height of approximately 20–23 m and a lateral spreading distance of approximately 35–38 m. The predicted flammable-cloud height and lateral extent are generally consistent with the experimental observations, indicating that the integral model can capture the main evolution characteristics of the gaseous hydrogen cloud after LH2 evaporation. The remaining discrepancies may be attributed to wind-speed fluctuations during the experiment, simplified environmental boundary conditions in the model, and the limited spatial resolution of the monitoring points. Overall, the comparison supports the use of the source-term integral model for predicting LH2 pool evaporation and the subsequent formation of the hydrogen cloud. Therefore, the transient liquid-pool radius and evaporation mass flow rate obtained from the integral model can be used as CFD inlet boundary conditions for the subsequent LH2 leakage and dispersion simulations in the station.

Figure 3: H₂ leakage experiment and model verification (a) H₂ concentration contour at 20.94 s in the NASA Test 6 experiment; (b) H₂ concentration contour obtained from the integral-model simulation.

3.3 Station Geometry Model

The layout of the station and its surrounding area can strongly influence hydrogen-cloud dispersion. A simplified geometric model is established based on the main structural features of the station, as shown in Fig. 4. The four sides are open. Natural wind enters from the left side and exits from the other three sides. The inlet wind speed was set to 2 m/s to represent a weak natural-ventilation condition in the semi-open station, consistent with the baseline ventilation condition commonly adopted in locomotive-related hydrogen dispersion simulations [38]. Nine support columns with dimensions of 1 m × 1 m × 6 m are arranged in the station. The locomotive length is 30 m. The center of the LH2 liquid pool is located at the coordinate origin, and the locomotive head is 50 m from the air inlet. The H2 dispersion simulation uses the same source conditions described above. The computational domain is shown in Fig. 4a, and the boundary conditions are listed in Table 2.

Figure 4: Station structure (a,b).

3.4 Basic Verification of the Species Transport Model

To verify the basic capability of the species transport model, the gas leakage and dispersion experiment conducted by Liu et al. [40] in a gas-cabin pipeline is used as a reference. It should be noted that this case is not dynamically or geometrically similar to the full-scale LH2 leakage and H2 dispersion scenario in the semi-open railway station. Instead, it is used only to examine whether the adopted species transport and turbulence models can capture the basic processes of gas advection, accumulation, and dilution under controlled release conditions. The experiment used non-toxic substitute gases (N2 for fuel gas and CO2 for air) in a rectangular pipe gallery with a cross-section of 85 mm × 170 mm and a length of 10 m. Gas was vertically released from an inlet 25 mm above the ground, with a mass flow rate of 1.587 × 10−4 kg/s. In the numerical simulation, multiple monitoring points are arranged. The same species transport and turbulence models as those used in the H2 dispersion simulation are adopted. The experimental and simulated results are compared to evaluate the feasibility of the CFD method for predicting gas dispersion.

Fig. 5 compares the experimental and simulated N2 volume fractions at the monitoring point over time. Both the experimental and simulated values show an overall increasing trend with leakage time, indicating that the model can capture gas accumulation and dispersion in confined spaces after release. The simulated values are lower than the experimental values at the early leakage stage, indicating a certain delay in the model response to early local concentration growth. As dispersion develops, the simulated curve gradually approaches the trend of the experimental data. The difference between them may be related to experimental measurement uncertainty, simplification of inlet release conditions, and local turbulent fluctuations. Overall, this comparison indicates that the adopted species transport and turbulence models can reproduce the basic trend of gas accumulation and dilution in a controlled duct experiment. However, because the N2/CO2 case differs from the present LH2 station scenario in release rate, length scale, temperature, and buoyancy behavior, it is used only as a basic verification of gas transport rather than as a direct validation of the full station-scale LH2 dispersion process.

Figure 5: Comparison of experimental and simulated N₂ volume fractions.

3.5 Flame Combustion Model after LH₂ Evaporation and Dispersion

The H2-air flammable cloud formed during LH2 evaporation and dispersion may burn when exposed to an open flame, high-temperature surface, electric spark, or other ignition source. Because H2 has a wide flammability range and a low minimum ignition energy, local flammable clouds formed during the dispersion stage have high ignition sensitivity. Therefore, after obtaining the H2 concentration and temperature fields in the station space, a combustion calculation model is established to analyze flame propagation and temperature-field evolution after ignition. The combustion calculation inherits the H2-air mixed-cloud distribution obtained from the dispersion stage and sets a local ignition source in the flammable region to trigger combustion. The ignition source is treated as an equivalent high-temperature region. The combustion simulation uses a pressure-based transient solver. The governing equations include mass, momentum, energy, and species transport equations. Turbulence closure uses the Realizable k-ε model. The H2 combustion reaction mechanism is imported through Chemkin and contains 10 species and 19 elementary reactions. The contribution of radiation heat transfer to the flame temperature field and the thermal impact on surrounding walls is also considered. The time step is determined according to flame propagation velocity, mesh size, and convergence stability, and time-step sensitivity checks are performed to ensure the reliability of transient combustion results. To verify the predictive ability of the combustion model for H2-air flame propagation, the H2-air combustion experiment performed by Liu et al. [41] is used for comparison. The experimental device is a cylindrical reaction vessel with an inner diameter of 0.4 m and a height of 2.3 m. The validation calculation establishes a two-dimensional axisymmetric model according to the experimental geometry and initial conditions, and a local ignition source is placed at the same ignition position as in the experiment. In the validation cases, the initial temperature is 373 K and the H2 volume fraction is 8%–12%. The main comparison is the flame-front position over time under different hydrogen concentrations. A small time step is used to capture the rapid propagation of the flame front, and the experimental results are used to evaluate the model prediction of the H2-air flame propagation trend. For the station scenario, the ignition condition is determined from the dispersion-stage concentration field, and a location is treated as ignitable when the local H2 volume fraction is within 4–74 vol%. In the ignition/combustion branch, the t = 3.0 s dispersion field is selected as the ignition initial field, and ignition is imposed at (−0.7, −3.6, 0.5) m. The same time coordinate measured from the onset of LH2 leakage is retained throughout the combustion calculation and is not reset after ignition.

Fig. 6 compares the experimental and simulated flame-front positions under different H2 volume fractions. The results show that when the H2 volume fraction is 10% and 12%, the simulation can reproduce the trend of flame-front propagation over time. When the H2 volume fraction is 6%, there is a certain deviation between the simulated and experimental results, but the overall increasing trend of the flame-front position remains basically consistent. The deviation may arise from simplification of the ignition process, the two-dimensional axisymmetric model approximation, uncertainty in the turbulent combustion model, and experimental measurement errors. Overall, the validation results show that the detailed H2 reaction mechanism and turbulent combustion model used here can reasonably describe flame propagation after ignition of an H2-air mixture and can be used for engineering-scale analysis of combustion consequences in the station scenario.

Figure 6: Comparison of experimental and simulated flame-front positions.

3.6 Grid-Independence Verification

To ensure the accuracy of the numerical results while maintaining reasonable computational efficiency, five mesh models containing 4.8 × 105, 5.0 × 105, 5.8 × 105, 6.4 × 105, and 7.1 × 105 cells were generated. Considering the influence of the support columns on the local flow field and hydrogen dispersion, local mesh refinement was applied in the vicinity of the columns to improve the resolution of the flow-field variations around these obstacles, as illustrated in Fig. 7a. Grid sensitivity was evaluated for both the dispersion and combustion stages. For the dispersion stage, the time-varying H2 volume fraction at the monitoring point (0, 0, 6) was selected as the evaluation parameter. For the combustion stage, the temperature distribution along a 30-m-long monitoring line oriented in the Y direction at a height of 1.6 m was compared under different mesh resolutions. As shown in Fig. 7b, differences in the transient H2 volume fraction at (0, 0, 6) decrease progressively with increasing mesh density. When the number of cells increases from 5.8 × 105 to 6.4 × 105 and 7.1 × 105, only minor differences are observed among the concentration curves. A similar tendency is observed for the combustion calculation, as shown in Fig. 7c, where the temperature distributions along the monitoring line obtained using the three finer meshes are in close agreement. Further mesh refinement therefore produces only limited changes in both the predicted H2 dispersion and temperature distribution. Considering the balance between numerical accuracy and computational cost, the mesh containing 5.8 × 105 cells was selected for the subsequent dispersion and combustion simulations.

Figure 7: Grid-independence verification (a) Section mesh at the column sectionv; (b) H₂ volume fraction at the monitoring point; (c) Temperature distribution along the Y-direction monitoring line.

4 Results and Discussion

4.1 LH₂ Evaporation and Dispersion in the Station

Fig. 8 shows the H2 distribution formed in the vertical mid-plane of the station and above the liquid pool after LH2 evaporation. After being released onto the ground, LH2 first undergoes liquid outflow, ground spreading, and heat-induced evaporation. A clear gas-phase H2 evaporation source term forms after approximately 2.8 s. Low-temperature H2 then forms a high-concentration region above the liquid pool and disperses upward and downwind under the combined effects of buoyancy and lateral wind. During t = 7.6–8.0 s, the already weak evaporation source continues to decay and the residual pool boundary shrinks. At about t = 8.5 s, the residual liquid pool nearly disappears and the inlet velocity approaches zero, marking the complete cessation of the evaporation source. Afterward, no new H2 is introduced into the computational domain. The existing H2 cloud continues to migrate and dilute under buoyancy, natural wind, and station-structure disturbances. This process shows that the early hazard after LH2 leakage is not instantaneous. Instead, it is a dynamic process jointly controlled by liquid-pool spreading, evaporation enhancement, and source-term decay. Fig. 8 shows the H2 distribution formed in the vertical mid-plane of the station and above the liquid pool after LH2 evaporation. After being released onto the ground, LH2 first undergoes liquid outflow, ground spreading, and heat-induced evaporation. A clear gas-phase H2 evaporation source term forms after approximately 2.8 s. Low-temperature H2 then forms a high-concentration region above the liquid pool and disperses upward and downwind under the combined effects of buoyancy and lateral wind. During t = 7.6–8.0 s, the already weak evaporation source continues to decay and the residual pool boundary shrinks. At about t = 8.5 s, the residual liquid pool nearly disappears and the inlet velocity approaches zero, marking the complete cessation of the evaporation source. Afterward, no new H2 is introduced into the computational domain. The existing H2 cloud continues to migrate and dilute under buoyancy, natural wind, and station-structure disturbances. This process shows that the early hazard after LH2 leakage is not instantaneous. Instead, it is a dynamic process jointly controlled by liquid-pool spreading, evaporation enhancement, and source-term decay.

Figure 8: H₂ volume-fraction distribution in the vertical mid-plane and above the pool (a–f).

Fig. 9 shows the time evolution of the front position and length of the H2 dispersion cloud in the station. The results indicate that at 4 s after LH2 leakage, the H2 cloud front has advanced 2.81 m compared with that at 3 s, and the cloud length reaches 10.06 m. The changes in front displacement and cloud length at this stage may be related to wind-speed conditions and the station computational-domain layout. At t = 5 s, the H2 cloud front is located at 9.25 m, and the cloud length is 16.48 m, increasing by approximately 2 m compared with the previous 2 s. This indicates a relatively high dispersion rate at this time. By 9 s, the cloud front has advanced to 22.42 m, and the cloud length has increased to 28.92 m. As dispersion continues, the H2 cloud gradually separates from the ground and further expands outward in the upper region of the station. To analyze H2 dispersion more intuitively, the velocity distribution of the 4% H2 volume-fraction iso-surface at 9 s after leakage is further presented.

Figure 9: Time history of cloud-front position and cloud length.

Fig. 10 shows the velocity distribution of the 4% H2 volume-fraction iso-surface in the station at 9 s after leakage, including front, side, and top views. The results show that the H2 cloud has almost surrounded the locomotive, and the H2 concentration around the locomotive is relatively high. Under buoyancy, most H2 migrates toward the station roof and disperses along the roof. The maximum cloud velocity at this time is 3.11 m/s, while the cloud velocity in some regions near the roof is relatively low. The top view indicates that the dispersion footprint of the H2 cloud is large and has covered the area around the locomotive, which may affect locomotive safety. As dispersion continues, the H2 cloud is expected to continue migrating under the wind field and further expand along the station roof. This t = 9 s result belongs to the non-ignited dispersion branch and is used only to characterize the spatial extent and migration of the flammable H2 cloud. It is not used as the initial field for the combustion calculation. At this time, the locomotive is already located within the lower-flammability envelope of the H2 cloud, indicating a substantial dispersion hazard if ignition does not occur earlier.

Figure 10: Velocity diagram of the station iso-surface at 9 s (gas concentration: 4%) (a–c).

Fig. 11 shows the H2 concentration distribution on the Z = 6 m plane at five time points in the late dispersion stage. After the leakage release stops, the flammable H2 cloud continues to disperse around the station and reaches a certain limiting dispersion range between the source position and the roof. Its overall position remains relatively stable. As the dispersion range of the H2 cloud continuously expands, the internal concentration gradually decreases until the cloud completely dissipates. The dissipation process lasts approximately 200 s. Therefore, the H2 concentration distributions at 40 s, 80 s, 120 s, 160 s, and 200 s are selected for comparison, and the influence of support columns on H2 cloud dispersion is analyzed. At 40 s, the H2 cloud is further diluted compared with the earlier stage, and the concentration decreases significantly, with a maximum volume fraction of 0.056. At 80 s, the cloud volume continues to increase, and air entrainment is enhanced. Due to obstruction by support columns, the dispersion speed in the downwind direction around the columns decreases, and the maximum volume fraction drops to 0.0184. At 120 s, the cloud area continues to expand outward, but the increase in volume is smaller than in the earlier stage. At 160 s and 200 s, the dispersion area of the H2 cloud changes only slightly, indicating that it gradually enters a passive dispersion stage. Due to vortices on the leeward side of support columns, the local H2 concentration decreases rapidly and the flow velocity is low, indicating that support columns have a certain influence on H2 dispersion.

Figure 11: Schematic diagram of H₂ cloud concentration variation in the late dispersion stage (a–e).

Fig. 12 shows the monitoring-point arrangement and the variation in H2 percentage. C7 is located directly above the LH2 liquid pool at (0, 0, 6). Its H2 concentration reaches 0.9 at 4 s and then gradually decreases to 0 at around 31 s. C8 is close to the inlet and is 10 m from the LH2 liquid pool. Its monitored concentration remains 0, indicating that H2 mainly disperses along the station roof and downwind. C9 is located at (0, 10, 6), and its H2 concentration reaches a peak value of 0.405 at 5.6 s, indicating local accumulation in this region during the early dispersion stage. C10 and C11 are located 10 m in the positive and negative X directions, respectively, with a Y coordinate of 0 and a height of 6 m. The results show that the H2 concentration at the monitoring point in the negative X direction is higher than that in the positive X direction, whereas the positive X direction detects H2 concentration changes earlier. This may be related to the blocking effect of the locomotive on dispersion in the negative X direction.

Figure 12: Monitoring-point locations and results (a,b).

4.2 H₂ Combustion Distribution in the Station

As shown in Fig. 13a, a monitoring line in the Y direction is placed beside the locomotive to evaluate flame combustion characteristics around the locomotive and personnel safety risk during the fire. The monitoring line is located at X = 0, with a total length of 30 m and a height of 1.6 m. Fig. 13b shows the temperature measured along the Y line over time. The results show that at 4 s, the flame has reached the front edge of the locomotive, and the local temperature can reach 979 K. At Y = 8.4 m, the temperature is about 300 K, indicating that the region within Y = 8.4 m is affected by high temperature at this time. The temperature peak appears at Y = −2.55 m, and this distribution may be related to the station spatial structure and wind-speed conditions. At 5 s, the average flame temperature decreases, and the maximum temperature is 849.6 K, which is 128.9 K lower than the peak temperature at 4 s. Nevertheless, it may still threaten personnel safety and equipment integrity. At 10 s, local high-temperature regions still exist near the Y line. By 15 s and 20 s, the temperature fluctuation along the Y line has weakened significantly and gradually returns to near ambient temperature. These results provide a reference for personnel safety assessment and protective-measure development during the early stage of a fire.

Figure 13: Y-line location and temperature variation (a,b).

A single monitoring dataset is insufficient for comprehensive assessment of accident hazards. Therefore, this study further analyzes the temperature variation over time at different heights in the direction perpendicular to the locomotive. After a fire occurs, personnel may evacuate from the area around or near the middle of the locomotive. To evaluate the thermal hazard in this region, a Z-direction monitoring line with a height of 6 m is arranged at X = 0 m and Y = 10 m, using the coordinate origin as the reference. Fig. 14a shows the position of the Z line in the station physical model. The temperature distribution at different heights along the Z line during the fire can provide a basis for personnel safety-risk assessment. Fig. 14b shows the temperature distribution along the Z line at different time periods. Unlike the Y-line monitoring results, the Z-line temperature at 4 s is close to ambient temperature, indicating that the flame has not yet propagated to this position. At t = 5 s, the temperature begins to rise at heights of 2 m and above from the ground, and the top temperature reaches the maximum value of 1882.3 K. This indicates that the flame propagates around the locomotive with a certain inclination angle. The temperature at a height of 1.6 m is about 300 K, indicating a weak thermal impact at this height. At 8 s, the flame-temperature peak is about 183 K higher than those at 6 s and 7 s, and the peak position changes. At 15 s and 20 s, the temperature in this region has decreased to near ambient temperature, indicating that most of the flame has moved away from this position and that the thermal hazard has decreased significantly. These results provide a reference for personnel evacuation and emergency response during fire accidents.

Figure 14: Variation in Z-line values and temperature at different times (a,b).

Fig. 15 shows the time variation of turbulent kinetic energy and temperature at a specific point on the station roof. The monitoring point is located at (0, 10, 6), and it is used to analyze the transient variations in local temperature and turbulent kinetic energy during flame combustion. The results show that at t = 5 s, the turbulent kinetic energy at this point reaches a peak of 2.14 J/kg, and the corresponding temperature is 1882.3 K, consistent with the Z-line monitoring results. Subsequently, both turbulent kinetic energy and temperature gradually decrease with time. The minimum turbulent kinetic energy is 0.0156 J/kg, corresponding to a temperature of 310.9 K. Because the station space is relatively open, flame propagation is weakly constrained except by the locomotive, roof, and support columns. The flammable H2 cloud mainly undergoes jet combustion along the roof. Therefore, flame propagation is less confined, and local turbulent kinetic energy and temperature decay rapidly. At 40 s, the turbulent kinetic energy is 0.03 J/kg and the temperature is 346.8 K, indicating that this point is still affected by heat to some extent, possibly due to radiative heat transfer and residual heat diffusion.

Figure 15: Relationship between turbulent kinetic energy and temperature over time.

Fig. 16 shows the peak flame temperature at different times during the station accident, which is used to analyze the maximum temperature that may be reached in the station area during fire development. The high temperature generated by the fire poses a significant threat to personnel safety and surrounding facilities. Therefore, determining the peak temperature at different times is important for risk assessment. The results show that the highest temperature in the station during the entire combustion process is 2193.5 K, appearing at t = 3.5 s at the location (−0.7, −3.6, 0.574). This indicates that the region beneath the locomotive is one of the main high-temperature hazard zones. During the early fire stage of t = 3–7 s, the peak flame temperature remains above 2000 K, indicating that the area around the locomotive remains at a high temperature. At t = 8 s, the highest flame temperature in the station decreases to 1859.87 K. At about t = 40 s, the maximum temperature around the locomotive decreases to 894 K, and the temperature in most regions is close to 300 K.

Figure 16: Peak flame temperature at different time periods.

Fig. 17 shows the temperature variation at the monitoring points over time. To analyze the thermal impact in the locomotive-tail region, four temperature monitoring points are arranged at the locomotive tail: D1(−1, 25, 0.35), D2(−1, 25, 1.9), D3(−1, 25, 3.9), and D4(−1, 25, 6). The monitoring results show that the locomotive-tail temperature increases with height. The temperature near the upper region is higher, whereas the near-ground temperature is close to ambient temperature. Among the monitoring points, D4 reaches the highest temperature of 604.7 K at 9.4 s. Overall, the thermal impact on the locomotive-tail region is relatively weak.

Figure 17: Monitoring-point data.

5 Conclusions

Based on an integral model and CFD numerical simulation, this study develops a full-process model of LH2 leakage, evaporation, dispersion, and combustion for station scenarios and systematically analyzes hydrogen dispersion and combustion behavior. The dispersion and combustion models are validated using experimental data. The results show that the developed model can reasonably reflect H2 cloud dispersion and flame-propagation characteristics. The main conclusions are as follows:

• (1)The source-term calculation shows a strongly transient LH2 evaporation process. The evaporation mass flow rate reaches a maximum of approximately 3.6 kg/s at 3.7–3.8 s, whereas the liquid-pool radius reaches its maximum value of approximately 1.2 m at 4.2 s. The evaporation mass flow rate decreases to a negligible level by approximately 7.5–7.6 s, while the residual liquid pool nearly disappears by approximately 8.5 s. Nevertheless, the previously formed H2 cloud continues to migrate and dilute for approximately 200 s, indicating that the dispersion hazard persists substantially longer than the active evaporation period.
• (2)Under the representative 6 kg LH2 release and 2 m/s wind condition, the flammable H2 cloud expands rapidly around the locomotive and toward the upper station region. In the non-ignited dispersion branch, at t = 9 s the 4 vol% H2 cloud front reaches 22.42 m and the cloud length reaches 28.92 m, with the flammable envelope extending around the locomotive. These results identify the source region, locomotive-adjacent region, and upper station space as important locations for hydrogen monitoring and emergency control.
• (3)The station geometry affects the migration and dilution of the H2 cloud. The maximum H2 volume fraction decreases from 0.056 at 40 s to 0.0184 at 80 s during the late dispersion stage. The roof promotes upper-level migration, while the locomotive and support columns modify the local flow and dilution patterns. Because these results are obtained for the present representative geometry and boundary conditions, the identified local retention characteristics should not be generalized to all station configurations.
• (4)Following ignition at t = 3.0 s in the ignition/combustion branch, the thermal hazard is concentrated mainly around the locomotive and upper station region during the early stage. The maximum flame temperature reaches 2193.5 K at t = 3.5 s and remains above 2000 K during t = 3–7 s. At about t = 40 s, the maximum temperature around the locomotive decreases to approximately 894 K, while the temperature in most regions approaches ambient conditions. These results indicate that early emergency response should prioritize locomotive-adjacent and upper station regions.