A Computational Fluid Dynamics Analysis of Hydrogen Leakage and Nitrogen Purging of a Solid Oxide Fuel Cell Stack

: A computational study of the nitrogen purging of a solid oxide fuel cell stack enclosed in a hot box is presented. The stack operates on ammonia as a fuel, and in the case of a hydrogen leakage, the entire compartment is immediately purged with nitrogen to ensure that there are no regions with high oxygen concentrations. In addition to this, the speed at which a hydrogen leak can be detected is determined. The results are then compared to a case with a relocated nitrogen inlet. A computational ﬂuid dynamics (CFD) model is developed using the Reynolds-averaged Navier–Stokes equations for compressible ﬂow in combination with conservation of energy and species equations in OpenFOAM. The results suggest that for the maximum concentration of oxygen to be below 5%, the hot box should be purged for 35 s, corresponding to 1.1 kg of nitrogen, if the hot box was already heated. If the hot box was at T = 300 K, it should be purged for 95 s, corresponding to 3.0 kg of nitrogen. The purge of the heated hot box results in a heat loss of 18 kW on average. A leak could be detected in 3.2 s during open circuit voltage tests. Changing the location of the outlet does not affect the cold purge, but results in a minimum purge period of 48 s during the hot purge, and the leak could be detected in 2 s. This paper demonstrates how CFD methods can be employed in order to address questions related to hydrogen safety.


Introduction
Fuel cell (FC) technology will become a key technology in order to reach the climate goals set by the European Union [1], the Paris Agreement [2] and the International Maritime Organization IMO [3].Emissions from fuel cells are considered net zero, provided that the fuels are produced using sustainable technologies.Emissions of products, such as NOx and SOx, are considered negligible when using fuel cells [4].Fuel cells commonly use hydrogen as their fuel, but other fuels, carrying hydrogen, such as methanol and liquid ammonia are also being investigated.One of the reasons is the low volumetric energy density of hydrogen compared to traditional fuels.For hydrogen to be utilized as a fuel, it has to be stored at a high pressure (≈700 bar) or in the liquid phase (≈−250 °C), and it will still have a low volumetric energy density compared to other fuels [5].In addition to this, there is still a lack of a distribution infrastructure of hydrogen, compared to e.g.methanol [6] or ammonia that can be used as direct fuel in high temperature fuel cells [7].Ammonia has the advantage of having an existing infrastructure, large production and proven safety procedures.Another benefit of ammonia, compared to hydrocarbons used in FCs, is that the only emission is water.Finally, ammonia can be created via various sustainable ways, as summarized in Figure 1.According to a recently published report by the European Marine Safety Agency, there could be enough green ammonia based on hydrogen generated via electrolysis to power the maritime fleet alone by 2040, although this sector will compete with others for green hydrogen [8].According to the same report, engine developments related to the use of ammonia are ongoing.However, issues related to concerns on nitrogen oxide (NO x ) and nitrous oxide (N 2 O) emissions, as well as the detrimental effects of ammonia slip from engines would need to be addressed [8].A different way of harvesting the energy of ammonia is to employ solid oxide fuel cells (SOFCs) to generate electricity.For the operation of SOFCs using ammonia, ammonia is almost completely cracked at 400 °C [7], meaning the high operational temperature of SOFCs in the range of 600-800 °C is very suitable.The practical efficiency when using ammonia was found to be around 40% in the first 1500 h of operation at a current density of 1 A/cm 2 [7].As the test conditions were not optimised, even higher efficiencies should be achievable.
In previous work studying SOFCs operating on ammonia, preliminary studies assessed the feasibility of the concept and demonstrated that no production of thermal NOx occurs in the anode [9][10][11].The degradation rates of SOFCs operating on both hydrogen and ammonia were investigated by Cinti et al., and they found similar values of 11.7 and 12.5 mV/100 h, respectively, at a temperature of around 745 °C [12].
The high operating temperature raises a question of safety, as hydrogen that invariably results from ammonia cracking has an autoignition temperature of 510 °C [13].The lower flammable level (LFL) of hydrogen is the lowest in the upward direction at 4% in air at standard temperature and pressure [13].The lower explosive level (LEL) is a hydrogen concentration of 10%.The lower oxygen concentration (LOC) is 5% at standard temperature and pressure.The lower limits decrease at higher temperatures [14].
In recent years, the Danish company Alfa Laval in Aalborg has investigated ammoniafueled SOFCs for maritime applications.For the SOFC system, an insulated hot box is used that contains an air-air heat exchanger, a fuel-fuel heat exchanger, an electric heater and the SOFC stack.Figure 2 is a CAD model of the hot box.It can be seen that the nitrogen inlet is at the front, lower right and the outlet at the upper back.To avoid explosive scenarios the hot box is purged with nitrogen.In case of a leak, a sensor is equipped at the outlet, detecting hydrogen concentrations above 0.4%.Several authors have investigated the risk of different hydrogen leaks using computational fluid dynamics (CFD), including hydrogen leaks from hydrogen fuelling stations [15], a detonation after a hydrogen leak in an enclosure [16] and hydrogen leaks in a laboratory [9], vehicles [10] and a ventilation facility [11].
In reference [15], Kikukawa et al. simulated a hydrogen leak from a fueling station nozzle.FLUENT 6.2 software (ANSYS Inc., Canonsburg, PA, USA) was used to set up the CFD model.The Reynolds-averaged Navier-Stokes (RANS) equations, using the standard k-ε as a turbulence model, for a compressible flow were solved, in combination with the conservation of species for an ideal gas mixture.The effects of a barrier were investigated, and it was concluded that a barrier would enhance the diffusion and shorten the travel distance of the hydrogen.
In reference [9], a hydrogen leak from a pipe in a laboratory was simulated.Ansys CFX was used to set up the CFD model.The RANS method was used, applying the Shear Stress Transport (SST) k-ω model due to supersonic flow.The transient terms were discretised using the second-order backward Euler scheme.All of the species, i.e., hydrogen, oxygen and nitrogen, were modelled individually.The ideal gas equation was used as the equation of state.Heat transfer at the walls was also included in this simulation.Prior to the hydrogen leakage simulation, the flow field in the laboratory in a normal state was simulated.These results were then used as the initial state of the leak simulation.The leak was simulated as a time-dependent flow.Experimental tests showed that a hydrogen leak would be detected 2-3 s after its occurrence by a sensor detecting a hydrogen concentration above 0.4%.
In reference [11], CFD in combination with experimental data was used to analyse a hydrogen leak in a semi-closed ventilation facility.The CFD code STARCCM+ was used.The turbulence was modelled using the realisable k-ε model in combination with the conservation of species and energy.The simulation was run at different ventilation flow rates and hydrogen leak pressures.
It is noted that all of the above studies relied on commercial tools that require the purchase of a user license.In contrast, OpenFOAM is a free, open-source CFD software (OpenFOAM v10) developed primarily by OpenCFD Ltd in 2004.It thus has the advantages of having no license cost, an ever-growing user base and complete access to the source code.The possibility of using OpenFOAM for gas dispersion was investigated by Mack and Spruijt [17] and Fiates and Vianna [18].
The former authors validated OpenFOAM-2.1-0 for heavy gas dispersions and compared their results to ANSYS Fluent [17].The model solved the RANS equations in combination with the conservation of species and energy, without considering reactions.The turbulence model used was the standard k-ε model.It was concluded that OpenFOAM was suitable to predict heavy gas dispersions at least on par with commercial codes.
Fiates and Vianna [18] used OpenFOAM for gas jet and gas dispersion modelling.The model solved the RANS equations in combination with the conservation of species and energy, without considering reactions.The convection terms were discretised using central difference and temporal upwind diffusion using the first-order Euler scheme.Pressure-velocity coupling was performed via the PIMPLE algorithm.Similar to this work, the standard k-ε and SST k-ω were investigated as turbulence models.No significant difference was found between the two models, and, importantly, good agreement with experimental data was obtained.
In the present article, OpenFOAM has been employed in order to: • Simulate the purging of the hot box with nitrogen and assess the amount of nitrogen required and how long it takes until the maximum oxygen concentration reaches 5%, given that it is filled with air at atmospheric pressure and initial temperatures of 800 °C and 300 °C, respectively.• Assess the cooling effect of the purge on the hot box.• Assess how long it takes for the hydrogen concentration to exceed 0.4% at the outlet if a leak occurs from the fuel-to-fuel heat exchanger during the OCV test.

•
Compare the above results with a case where the outlet of the hot box has been relocated.The defined leak location is illustrated in Figure 3 and the relocated outlet is illustrated in Figure 4.

Mathematical Model 2.1. Main Assumptions
The list below gives an overview of the main assumptions: • The gases are assumed as ideal.

•
The flow is assumed as incompressible.

•
Buoyancy effects are neglected during the purge simulations.

•
For the modelling of diffusive transport, both the laminar and turbulent Schmidt numbers are equal to one.

Governing Equations
The solvers used in this project are ReactingFoam and RhoReactingBuoyantFoam, which are transient solvers that can handle turbulent, compressible and reacting flows.The governing equations solved are the conservation equations for mass, momentum, species and energy.The conservation equations can be written as stated in Equations ( 1)-( 5) [19,20]. Energy where ρ is the density, t is the time, u is the velocity vector, p is the pressure, g is the gravitational acceleration vector, τ is the stress tensor, µ is the dynamic mixture viscosity, I is the identity matrix, y i is the mass fraction of species i, h is the enthalpy, K is the kinematic energy and α eff is the effective thermal diffusivity.The conservation of species uses the effective viscosity, instead of the effective mass diffusivity, given as ν eff = ν lam + ν turb .This results in the assumption that the laminar and turbulent Schmidt numbers, Sc lam and Sc turb , are equal to one, as stated above.The effect of the turbulent Schmidt number has been reviewed by Gualtieri et al. [21].No universal method to determine the value was concluded, but fitted values for different flow configurations were found and suggest that the value varies locally.
In this work, the conservation of species (Equation ( 4)) is modified to use a specified constant value of Sc turb and Sc lam , meaning D eff is given as: where D eff is the effective mass diffusivity.
During the purge, buoyancy is neglected (g = 0).In OpenFOAM, the pressure gradient and gravity force terms are rearranged as [22]: Here, p rgh = p − ρg • r and r is the position vector.

Turbulence Models
The diameter of the nitrogen purge inlet is 18.5 mm.Subtracting the volume of the components of the HB, the volume of the HB is 1.69 m 3 .The volume flow of nitrogen corresponds to the total volume over a minute.This results in the inlet Reynolds number (Re) given in Equation (8): where V is the volume flow, D is the diameter and ν is the kinematic viscosity of nitrogen.The Reynolds number is above 10,000; this means that the entering flow is fully turbulent.The characteristic length of the hot box is however much larger than that of the inlet, meaning the turbulence is expected to have little impact.During a leak, Re leak,in < 2300, meaning the flow is laminar.A laminar simulation is therefore also compared in this case.
In the literature, the commonly used turbulence models are the SST k − ω, the standard k − ε and the realisable k − ε models [9][10][11]15,17,18]. Based on the literature, this suggests that a RANS model is sufficient and LES is not required.
In this paper three different turbulence models are tested, including the standard k − ε, the SST k − ω and the k − ε − φ − f model.The transport equations for the standard k − ε model are where D k is the effective diffusivity of k, P is the turbulent kinetic energy production rate, D ε is the effective diffusivity of ε and C 1−2 are model coefficients.The standard coefficients of OpenFOAM are used.The SST k − ω model substitutes ω = ε/k into Equation (10) and introduces the blending function F 1 .The transport equation is given as [22]: +S ω (11) where D ω−kω is the effective diffusivity, γ and β are model coefficients, G is the turbulence generation and S ω is a source term.The difference between the transport equation of ε and ω is the term ρCD kω .When The k − ε − φ − f model introduces the transport equation for φ, which is the normalized wall-normal fluctuating velocity scale, and f , the elliptic relaxation factor.The transport equations are, respectively [22,23]: where T and L are the time and length scale used by the k − ε − φ − f model.

Boundary and Initial Conditions
Three scenarios are to be investigated: nitrogen purge of the hot box filled with hot air, nitrogen purge of the hot box filled with ambient air and a fuel leak during OCV testing.The most important boundary and initial conditions will be determined.
For the inlet, a fixed volume flow of 1669 L/min during the purge simulations is applied.This volume flow corresponds to the entire volume of the hot box over 60 s.During the leak simulation, a fixed mass flow of 1.35 kg/s is used.
Fixed temperatures of 300 K and 1073 K are used for the components and at the inlet during the purge simulations and leak simulations, respectively.At the walls and outlet, a zero gradient is used, corresponding to adiabatic conditions (∂φ/∂n = 0).
The molar composition at the inlet is pure nitrogen during the purge.During the purge, it is composed of x H 2 = 0.6 and x N 2 = 0.4.Initially, inside the hot box, the composition is x N 2 = 0.79 and x O 2 = 0.21 during the purge and x N 2 = 0.97 and x O 2 = 0.03 during the leak.
Finally, the inlet values of the turbulence parameters are the standard boundary conditions used in the OpenFOAM tutorials.Turbulent intensities of 0.034 and 0.01 are defined at the inlet during the purge and leak, respectively.Similarly, a mixing length of l m = 0.07D is defined at the inlet during the purge and leak, respectively [24].
At the outlet, the pressure is fixed at 0.995 bar, and elsewhere a zero gradient is used.During the leak simulation, buoyancy is considered, introducing p rgh .In this case, p rgh uses the boundary conditions of p and p is calculated at the boundaries based on p rgh .

Numerical Verification
The PIMPLE algorithm was used, which is a combination of the SIMPLE and PISO algorithm [24].The maximum Courant number (Co = ∆t 2V i ∑ j∈ f aces u i • n i,j A i,j ) used during the purge simulations is Co max = 50 and Co max = 0.5 during the leak simulations.The mean Co is however much lower (in the order of 1 × 10 −2 during the purge).Different values of Co max were tested to ensure the results were independent of the time step.The time derivatives are discretised using first-order Euler scheme.The upwind scheme is used for the discretisation of the turbulent advection terms.The rest of the terms are discretised using second-order accurate schemes.Different schemes were tested for the simulations.
Grid-independent studies were carried out for the purge in the case of a high initial temperature and for the leak.The inspected parameters for the purge were the volume averaged x O 2 in the entire hot box, near the fuel cell stack and near the fuel-fuel heat exchanger, |u| out and p in .The inspected parameters for the leak were the volume-averaged x H 2 concentration in the entire hot box, the mean x H 2 at the outlet and |u| out .The resulting mesh had a size of 7.3 × 10 5 and 6.5 × 10 5 cells for the purge and leak simulations, respectively.

Hot Purge
The simulations for this purge use Sc lam = Sc turb = 1.The maximum molar oxygen concentration is illustrated in Figure 5 for all three turbulence models.The time during the purge simulations was non-dimensionalised by dividing by T purge = 60.This is the initial purge period suggested by Alfa Laval and the period for a volume of nitrogen to be injected corresponding to the volume of the hot box.The maximum x O 2 is below the LOC when the purge is already 60% complete.This corresponds to 1.33 kg of nitrogen being used.The threshold is crossed first for the standard k − ε model and lastly for the SST k − ω model in the period 0.552 < t/T purge < 0.582.
Figure 6 illustrates the distribution of x O 2 in the hot box at different timestamps.At t/T purge = 0.5, the probability of a cell having x O 2 > 0.5 is less than 15.5% for the standard k − ε model, 37.8% for the SST k − ω model and 36.5% for the k − ε − φ − f model.At the end of the purge period and in the steady state, x O 2 is nearly constant.
Figure 7 illustrates the oxygen concentration at different timestamps in the hotbox, using a lower threshold of 5%, meaning only concentrations above 5% are shown.
As x O 2 ,max < 0.5, at t/T purge = 1 and in the steady state, the box appears empty.At t/T purge = 0.5, regions above the inlet and between the FC and HXs still have an x O 2 > 0.5.The visualised distribution is shown for the standard k − ε.
Figure 8 is the normalised change in energy of the hot box on the left, the normalised temperature at the outlet in the middle and the normalised energy gradient on the right.The energy is normalised by dividing by E initial = 520 kJ, and the temperature is normalised by dividing by T initial = 1073 K.The dashed curve for comparison is if the entire volume of nitrogen was heated to T initial , meaning the maximum possible loss.The system loses energy corresponding to 72% of the total mass of the purging nitrogen when heated from T = 300 K to T = 1073 K.This could affect the performance of the FC and HXs, as the gas inside the hot box would cool down the surface of the different components in the hot box.The temperature distribution inside the hot box is illustrated in Figure 9.
The mean temperatures at t/T purge = 0.083 and t/T purge = 1.0 are T/T HB,initial = 0.71 and T/T HB,initial = 0.65, respectively.The temperature scales are fixed to 0.79 < T/T HB,initial < 0.93 and 0.65 < T/T HB,initial < 0.84, respectively, to better visualise the distribution.

Cold Purge
The maximum and mean x O 2 is illustrated in Figure 10 for the standard k − ε and SST k − ω models.The simulation using the standard k − ε model was modified to use Sc turb = 0.9 and Sc lam = 0.76.The mean concentration of oxygen is below 5% when t/T purge > 1.43 and the maximum concentration occurs when t/T purge > 1.58.The standard k − ε achieves x O 2 < 0.05 and x O 2 ,max < 0.05 first and the SST k − ω model achieves this last.The k − ε model using a specified value of Sc (Sc lam = 0.76 and Sc turb = 0.9) is slightly faster at achieving concentrations below 5%; the relative difference between all the simulations is, however, <2%.
Figure 11 illustrates the oxygen concentration at different times in the hot box, using a lower threshold of 5%.When t/T purge = 1.5, there is only a spot near the FC which is above x O 2 > 0.05.It can be seen that the oxygen distribution is very uniform at any given time.

OCV Leak
The leak simulations were performed using Sc lam = 0.3 and Sc turb = 0.8, as well as simulations using the standard k − ε model with Sc lam = Sc turb = 1. Figure 12 shows the mean and maximum concentration of hydrogen at the outlet, and Figure 13 shows the distribution of x H 2 visualised inside the hot box, with a threshold of x H 2 > 0.04.
The molar fraction of hydrogen exceeds x H 2 > 0.004 in the period 2.81 s < t < 3.20 s, meaning it will take 2.81 s at minimum to detect the leak, according to the simulation.In the time period 3.05 s < 3.65 s, the molar fraction exceeds x H 2 > 0.04.The difference between the maximum and mean values for each simulation is very low and the detection time is not affected using one or the other (∆t < 0.01 s).At t > 15.9 s, x The distribution of hydrogen inside the hot box is shown at time t = 2 s, t = 4 s, t = 8 s and t = 16 s in Figure 12.The hydrogen leak starts to disperse and form a cloud at the top of the hot box due to buoyancy.At t = 16 s, the flammable region is the top half of the hot box.For the hydrogen to ignite or explode, however, a sufficient concentration of oxygen is also required.

Alternate Outlet
Figure 14 shows x O 2 ,max and x O 2 for the original outlet and the new outlet, using the SST k − ω turbulence model, when T HB,initial = 1073 K. Figure 15 shows x O 2 ,max and x O 2 for the original outlet and the new outlet, using the SST k − ω turbulence model when T hot box,initial = 300 K.For the new outlet, x O 2 ,max<0.05when t/T purge > 0.795, and this occurs in the original outlet when t/T purge > 0.582.This means the original is 36.6%faster and needs 36.6% less nitrogen.Finally, a leakage during the OCV test was also simulated for the new outlet.Figure 16 shows x H 2 ,out,max and x H 2 ,out for the original outlet and the new outlet, using the standard k − ε turbulence model.The new outlet can detect the leakage after 1.94 < t < 1.95 s.This is 1.2 s faster and corresponds to 62.7% of the original time.

Figure 2 .
Figure 2. CAD drawing of the hot box investigated in this study.

Figure 3 .
Figure 3. Hot box with a hydrogen leak.The left figure is from the front and the right is zoomed in on the back.The hole diameter is 46.7 mm.

Figure 4 .
Figure 4.The hot box with a relocated outlet.The new outlet is marked in green, the original in red and the inlet in blue.The inlet has a diameter of 18.5 mm.

Figure 5 .
Figure 5. Maximum molar oxygen fraction during the purge.The red dotted line represents the LOC.

Figure 6 .
Figure 6.Histogram with x O 2 on the x-axis and probability on the y-axis.Bar width is 0.0005.Top left is at time t/T purge = 0.25, top right is t/T purge = 0.5, bottom left is t/T purge = 1 and bottom right is in the steady state.

Figure 7 .Figure 8 .
Figure 7. Molar fraction of oxygen at different time stamps in the hot box, with a threshold above 0.05, for the simulation using the standard k − ε model.

Figure 9 .
Figure 9.The temperature distribution in the hot box at t/T purge = 0.083 and t/T purge = 1.0.The temperature scales are fixed to 0.79 < T/T HB,initial < 0.93 and 0.65 < T/T HB,initial < 0.84, respectively, to better visualise the distribution.

Figure 10 .
Figure 10.The mean and maximum molar oxygen fraction during the purge.The red dotted line represents the LOC at x O 2 = 5%.The standard k − ε model uses Sc lam = 0.76 and Sc turb = 0.9 instead of Sc = 1.

Figure 11 .
Figure 11.Molar oxygen fraction at different time steps in the hot box, with a threshold above 0.05%, for the simulation using the SST k − ω model.

Figure 12 .
Figure12.The mean and maximum molar hydrogen fraction at the outlet in the time period 0 < t < 16 s.The results are displayed for the three turbulence models, a laminar model and an unmodified k − ε model.The red dotted line is the minimum concentration for the sensor to detect hydrogen and the red line is when the sensor will give sound an alarm.

Figure 13 .
Figure 13.Distribution of hydrogen in the hot box, with a threshold of x H 2 > 0.04.The top left figure is at t = 2 s, top right is at t = 4 s, bottom left is at t = 8 s and bottom right is at t = 16 s.
H 2 > 0.1 for all the models, except the unmodified model and the k − ε − φ − f model.The model in which the outlet molar fraction increases the slowest is the unmodified standard k − ε model and the fastest is the SST k − ω model.The k − ε − φ − f model starts to become unstable when t > 4.4 s, as mass balance is not achieved.

Figure 14 .Figure 15 .
Figure 14.Mean and max molar concentration of oxygen for the new and original outlet vs. normalised time for T HB,initial = 1073 K.