DOI: 10.3303/CET24111072 Paper Received: 24 January 2024; Revised: 29 May 2024; Accepted: 26 July 2024 Please cite this article as: Schiaroli A., Baldassarri A., Ustolin F., 2024, Analysis of Accidental Hydrogen Releases in the Glass Manufacturing Industry, Chemical Engineering Transactions, 111, 427-432 DOI:10.3303/CET24111072 CHEMICAL ENGINEERING TRANSACTIONS VOL. 111, 2024 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Valerio Cozzani, Bruno Fabiano, Genserik Reniers Copyright © 2024, AIDIC Servizi S.r.l. ISBN 979-12-81206-11-3; ISSN 2283-9216 Analysis of Accidental Hydrogen Releases in the Glass Manufacturing Industry Alice Schiarolia,b,*, Alberto Baldassarria,b, Federico Ustolina aNorwegian University of Science and Technology – NTNU, Department of Mechanical and Industrial Engineering, Richard Birkelands vei 2, 7034, Ttrondheim, Norway bAlma Mater Studiorum – Università di Bologna, Dipartimento di Ingegneria Civile, Chimica, Ambientale e dei Materiali , via Terracini 28, 40131, Bologna, Italy alice.schiaroli@unibo.it Glass is one of the most ubiquitous materials in the world. Due to the extremely high temperatures required in the melting process, the glass industry is considered a hard-to-abate sector and poses major challenges to meet the net-zero target in the next decades. Since the highest share of emissions from glass production stems from the combustion of natural gas, its replacement with hydrogen is considered a promising solution to reduce the sector's environmental impact. This is the aim of the H2GLASS project, launched by the European Union at the beginning of 2023. In this context, addressing hydrogen-safety-related aspects is a top priority. In fact, due to hydrogen’s peculiar flammability properties (e.g., wide flammability range, low ignition energy), its utilisation in furnaces may pose significant risks. In this study, a computational fluid dynamic (CFD) model is developed in Ansys Fluent to investigate hydrogen diffusion in enclosures following an accidental release. A grid and time-step sensitivity analyses are carried out to identify the best setup. The model is then validated against experimental data. The results of this work can be used as the starting point to build a CFD model suitable for studying hydrogen releases in large domains, such as glass manufacturing facilities, where obstacles and mechanical ventilation systems are present. 1. Introduction In recent years, hydrogen has gained ever-increasing attention as an alternative to fossil fuels in the framework of the energy transition. Indeed, it has the potential to be employed in the decarbonisation of sectors where other technologies, such as electrification, might prove challenging to implement. One of these sectors is the glass industry, which is extremely energy-consuming due to the high temperatures, around 1500 °C (Jatzwauk, 2021), required in the melting and fining processes. In 2007, the average specific energy consumption in the EU glass industry was estimated to be 7.8 GJ per tonne of saleable glass (Schmitz et al., 2011), compared to 15.7 GJ per tonne of crude steel (Pardo et al., 2012) and 3.7 GJ per tonne of clinker (Ecofys, 2009). Currently, in traditional glass furnaces heat is provided through direct flame, fuelled by the combustion of natural gas in specific burners, resulting in CO2 and other greenhouse gas (GHG) emissions. In 2022, the glass sector produced 22 Mt of CO2,eq in the EU (CINEA, 2022), accounting for 2.7% of the total emissions from industry (Eurostat, 2023), and 95 Mt worldwide (CINEA, 2022). Several decarbonisation options exist to mitigate the environmental impact of the glass sector (Furszyfer Del Rio et al., 2022), but the use of hydrogen as a fuel appears to be the most promising one. For instance, Zier et al. (2023) highlighted that transitioning to hydrogen combustion is the only way the German glass industry has to comply with the loose 2 °C temperature increase target established by the Paris Agreement. Nevertheless, the application of hydrogen technologies in hard-to- abate sectors is challenging. In the attempt to speed up the process, the EU has funded and launched, at the beginning of 2023, the H2GLASS project – advancing Hydrogen (H2) technologies and smart production systems TO decarbonise the GLass and Aluminium SectorS –, whose primary goal is to increase the technology readiness level of hydrogen combustion in glass furnaces (H2GLASS, 2023). The project also places particular emphasis on safety aspects. In fact, the use of hydrogen, if not properly handled, is known to pose significant risks. An accidental release within a confined environment can quickly result in a dangerous build-up of hydrogen 427 with the consequent formation of an extremely explosive atmosphere. Considering that hydrogen has a wide flammability range, 4-75% vol (Yang et al., 2021), a very low minimum ignition energy, 0.018 mJ (Yang et al., 2021), and a strong tendency to detonation (Yang et al., 2021), it is fundamental to carefully assess and evaluate its dispersion in enclosures. For this type of investigation, computational fluid dynamic (CFD) models are valuable and powerful tools. Zhang et al. (2021) and Patel et al. (2023) used CFD models developed in Ansys Fluent to analyse the effect of the leakage position, direction, and vent openings arrangement, finding that lateral vents close to the ceiling and roof vents above the leakage source are the most effective combination to reduce hydrogen concentration. Furthermore, Hussein et al. (2020) simulated in Ansys Fluent a release from the tank of a hydrogen fuel cell vehicle (HFCV) in an underground car park with natural ventilation, to evaluate the effect of the leakage diameter on the dispersion of hydrogen. Li et al. (2021), analysed a similar release scenario using the CFD code GASFLOW-MPI, simulating also immediate and delayed ignition of the hydrogen-air mixture into a tunnel. The effect of ventilation has been studied by Houf et al. (2013), who simulated in FUEGO and FLACS the release from a fuel cell forklift in an industrial warehouse, finding a limited effect in reducing the hydrogen concentration. The same release was simulated in FLACS by Lucas et al. (2021) in a much larger environment and with obstructions, highlighting a modest influence of passive ventilation on reducing blast overpressure. For laminar flows, the effect of ventilation was assessed solely through empirical formulas (Bauwens and Dorofeev, 2014). At present, the development and validation of a CFD model capable of predicting the behaviour of a laminar hydrogen leak within a confined environment are still missing. This work is the starting point of the construction of this model. 2. Methodology In this section, the methodology used for the model setup is presented. Section 2.1 provides the numerical details of the model, while section 2.2 describes the case study used for its validation. 2.1 CFD setup The model was developed in Ansys Fluent (2023), in which the finite volume method is used to solve the governing equations of fluid mechanics. For every fluid flow, the Navier-Stokes equations are solved to ensure mass (Eq. (1)) and momentum (Eq. (2)) conservation. ∂ρ ∂t + ∇̅ ∙ (ρu̅) = 0 (1) ∂ ∂t (ρu̅) + ∇̅ ∙ (ρu̅u̅) = - ∇̅P + ∇̅ ∙ T̿ + ρg̅ (2) Where ρ is the density of the fluid, u̅ is the velocity vector, P is the pressure, and g̅ is the gravitational acceleration. T̿ is the stress tensor, expressed as follows (Eq. (3)): T̿ = μ [(∇̅u̅ + ∇̅u̅T) - 2 3 ∇̅⋅ u̅ I]̿ (3) Where μ is the molecular viscosity of the fluid and I ̿is the unit tensor. In addition, when there is more than one chemical species in the system, the mass conservation equation is introduced for a single component, as shown in Eq. (4): ∂ ∂t (ρYi) + ∇̅ ⋅ (ρu̅Yi) = - ∇̅ ⋅ Ji̅ + Ri (4) Where Yi is the mass fraction of the species i, Ri is its net rate of production by chemical reaction (not present in this case) and Ji̅ is the diffusion flux, calculated through the Fick’s law (Eq. (5)). Ji̅ = - ρDi,m∇̅Yi - Di,T ∇̅T T (5) Where Di,m is the mass diffusion coefficient for the species i in the mixture, Di,T is the thermal diffusion coefficient, and T is the temperature of the fluid. Lastly, when heat transfer has to be modelled, the energy conservation (Eq. (6)) equation is introduced. ∂ ∂t (ρE) + ∇̅⋅ [u̅(ρE+P)] = ∇̅ ⋅ (k∇̅T - ∑ hiJi̅ i + T̿ ⋅ u̅) (6) 428 Where E is the total energy of the fluid, k is the thermal conductivity and hi is the sensible enthalpy, defined as the sum of the enthalpies of each component in the mixture, weighted by their mass fraction. As for the numerical schemes and methods employed in the simulation, the SIMPLE (Semi-Implicit Method for Pressure-Linked Equation) algorithm was selected for the pressure-velocity coupling. For the spatial discretisation, the Second Order Upwind scheme was used for all the variables of interest, to improve the accuracy of the solution. Finally, the transient formulation was set to Bounded Second Order Implicit, which ensures high stability and robustness even for long time steps. This formulation requires higher computational power but leads to very precise results. However, for the sake of comparison, a simulation was also conducted using the SIMPLEC (Semi-Implicit Method for Pressure-Linked Equation-Consistent) algorithm, which improves the SIMPLE algorithm to achieve faster convergence, and the PRESTO! (PREssure STaggering Option) spatial discretisation scheme for pressure, which is considered to give a more physically accurate representation of the pressure field, compared to the Secon Order Upwind scheme. 2.2 Case study The experimental study carried out by De Stefano et al. (2019) was chosen for the validation because it exhibited some of the characteristics of a potential hydrogen release in a glass production facility (e.g., confined environment and low-velocity releases).. The experimental domain was a closed unventilated chamber with dimensions 470 mm × 330 mm × 200 mm and a total volume of 31 litres. A symmetry condition could be exploited to split the computational domain along the y-axis, reducing the computational effort of the simulations. The release hole, from which hydrogen was discharged vertically downward, was located on the ceiling of the box and its diameter was 4 mm. Furthermore, four rods of sensors were positioned into the chamber, but due to the symmetry condition only two of them were included in the simulation. The geometry shown in Figure 1a was employed to reproduce the chamber. Figure 1: (a) Geometrical domain, with hydrogen inlet and sensors; (b) medium mesh considered in the simulation. The xz plane is the symmetry plane. The experiments were conducted, as reported by the authors, at atmospheric pressure and room temperature, which was considered to be 20 °C. This value was used in the simulations both as the chamber’s initial temperature and as a boundary condition for the hydrogen inflow. The release flow rate was 0.1 Nm3/h, and the release duration was 53 seconds. It is worth mentioning that with these data, the calculated total released volume is lower than that declared in the study. As stated by the authors, this discrepancy was due to an error in the mass flow rate control instrument. To compensate for this gap, a higher flow rate, calculated so that the total volume released matched the value of 1.67 ×10 -3 Nm3, was adopted in this study for the first 10 seconds.. The resulting mass flow rate, defined by means of a User-Defined Function (UDF), was halved due to the symmetry condition, leading to the following expression (Eq. (7)): m ̇ = { 2.139×10 -6 kg/s, if t ≤ 10 s 1.249×10 -6 kg/s, if 10 s < t ≤ 53 s (7) Two different diffusion coefficients, that took into account the effect of temperature, were evaluated: one determined experimentally by Scott and Cox (1960), which is 0.763 cm2/s, and the other obtained from Winkelmann's empirical equation, which is 0.679 cm2/s. The effect of the equation of state was assessed as well, testing both the ideal gas and the Peng-Robinson equation of state, to see if a potential pressure increase at the hydrogen inlet could lead to a deviation from ideality. (a) (b) 429 In addition, to ensure that the results were not dependent on the mesh or on the time step employed in the simulations, a sensitivity analysis was performed, testing three different degrees of grid refinement and three different time step lengths. The meshes created for this validation were of the structured type, consisting exclusively of hexahedral elements. The cylindrical area around the inlet was discretised using the C-grid technique. The total number of elements was 33,600 for the coarse mesh, 120,800 for the medium mesh (see Figure 1b), and 299,840 for the fine mesh. The time step lengths were 0.2 s, 0.1 s, and 0.05 s. 3. Results Figure 2a displays the hydrogen concentration profiles obtained with three different levels of grid refinement. The curves represent sensors B4 and B1, selected in order to capture the upper and the lower concentration limits, since they were positioned at the greatest and lowest heights, respectively. Figure 2: (a) Grid sensitivity analysis; (b) results of the simulations compared to experimental data. It can be observed that for sensor B1, the three curves are perfectly overlapped, while for sensor B4, the coarse mesh produces a slightly lower concentration profile compared to those produced by the medium and fine meshes. Therefore, the results of the sensitivity analysis indicate that the medium mesh ensures grid independence. As for the time step, since for both sensors all three curves were almost identical, it can be concluded that the largest of the tested values (0.2 seconds) ensures time step independence. Figure 2b shows the comparison between the concentration profiles obtained in the simulation and the experimental data. The graph includes data from all the sensors in rod B, while for rod C only C4 and C3 are represented. For the sake of brevity, sensors C1 and C2 were excluded because, except for sensor C4, experimental observations showed perfect horizontal homogenisation (i.e., sensors at the same height produced the same concentration curves over time). The curve for sensor C3 confirms that the simulation is able to capture diffusion in the horizontal direction. Overall, it can be observed that the simulation curves reproduce the experimental data with excellent agreement, although they slightly underestimate the concentration detected by some sensors. For instance, for sensor B3, the numerically predicted peak concentration is 6.38%, while in the experiment it was 6.97%, resulting in an error of 8.5%, which is still within an acceptable range. Instead, for sensors B4 and C4, the errors are as low as 0.7% and 1.8%, respectively. For sensors B2 and B3, the peak value corresponds to the final concentration value, when complete homogenisation of the hydrogen-air mixture inside the chamber is achieved. In general, the average concentration value, which is 5.46% in the simulation, is lower than that observed in the experiment, which was about 6%. However, the trend over time of the simulated curves is very similar to that of the experimental ones. Indeed, the concentration gradient (defined as the difference between the highest and the lowest concentration values) 20 seconds after the end of the release is 7.89%, while in the experiment it was 7.86%, resulting in a 0.4% error. Finally, it can be observed that the curves obtained from the simulation take a slightly longer time to converge to the final value, compared to the experimental ones. In fact, the experimental homogenisation time was reported to be 183 seconds, at which the concentration gradient was 0.3%. In the simulation, this concentration difference was barely reached at the final time of the test, which was 220 seconds. It is worth noticing that the above-mentioned data were obtained using the ideal gas equation of state and a diffusion coefficient of 0.763 cm2/s, as explained in Section 2. However, a diffusion coefficient of 0.679 cm2/s, determined through Winkelmann’s empirical equation, was tested as well. For the higher sensors, the differences were negligible, but they became more significant moving downwards, though always remaining (a) (b) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0 25 50 75 100 125 150 175 200 225 H 2 m o la r fr a c ti o n Time (s) B4 C B4 M B4 F B1 C B1 M B1 F 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0 25 50 75 100 125 150 175 200 225 H 2 m o la r fr a c ti o n Time (s) B4 exp C4 exp B3 exp B2 exp B1 exp B4 sim C4 sim B3 sim B2 sim B1 sim C3 sim 430 small. Additionally, with a lower diffusion coefficient, the curves did not even reach full convergence at the end of the simulation. The main figures and errors obtained in these simulations are reported in Table 1. Table 1: Main figures of the simulations, with different diffusion coefficients (D); simulation 1: D = 0.763 cm2/s; simulation 2: D = 0.679 cm2/s. Experiments Simulation 1 Error sim. 1 Simulation 2 Error sim. 2 B4 max molar fraction (%) 13.75 13.66 0.7% 13.87 0.9% C4 max molar fraction (%) 13.02 12.79 1.8% 12.89 1.0% B3 max molar fraction (%) 6.97 6.38 8.5% 6.21 10.9% B2 max molar fraction (%) 6.00 5.35 10.8% 5.13 14.5% B1 max molar fraction (%) 5.90 5.28 10.5% 4.97 15.8% Gradient at 73 s (%) 7.86 7.89 0.4% 8.61 9.54% Homogenisation time (s) 183 220 20.2% The effect of the equation of state was also assessed, conducting a simulation with the Peng-Robinson real gas equation of state to see if a potential pressure increase at the hydrogen inlet could lead to a deviation from ideality, but identical profiles were obtained. Also changing the pressure-velocity coupling algorithm to SIMPLEC and the spatial discretisation scheme for pressure to PRESTO! gave very similar results. 4. Discussion The simulations conducted for the validation study are in excellent agreement with experimental data, resulting in concentration profiles over time that are slightly lower in value, but with a very similar trend. It should be noticed that the final concentration, once complete homogenisation is achieved, is lower than the experimental one. The final average concentration in the simulations perfectly corresponds to the ratio between the moles of hydrogen introduced into the chamber and the total moles of the air-hydrogen gas mixture. Therefore, an error in the mass balance committed by the solver can be ruled out. Any discrepancy between volumetric fractions reported in the paper and molar fractions calculated in the simulation can also be ruled out at room temperature and atmospheric pressure. The discrepancy could be due to a different quantity of hydrogen actually released, compared to what was declared in the study used for the validation, a smaller amount of air present in the chamber before the release, or imperfect sealing of the chamber, resulting in air leakage. This aspect could have been clarified by knowing the transient pressure profile inside the chamber and comparing it with the expected linear increase. Another explanation could be the reduction of the domain due to the non-negligible volume occupied by the sensor rods. Since they were not included in the simulation geometry, this may have caused a lower value of the final hydrogen concentration. However, it is highly improbable that the rods cause a volume reduction sufficient to cause a difference in the concentration from about 6% to 5.46%. In fact, a set of experimental tests was conducted with obstacles placed inside the chamber and the final concentration reported in the paper was the same as the case without them, even though they occupied 10% of the enclosure volume, as stated by the authors. These aspects support the hypothesis of imperfect sealing and consequent air leakage during the tests, which also explains why the discrepancy between the experimental data and the simulated curves becomes noticeable only in the later stage of the simulations. Furthermore, given the longer homogenisation time, it can also be stated that the model moderately underestimates the diffusion phenomenon. This may be due to the fact that the temperature at which the experiment was conducted was slightly higher than what was assumed. Further tests could be carried out with a higher diffusion coefficient even if, from a numerical perspective, this parameter is significant only for very low release velocities. However, the overall differences are relatively small. Thus, it can be stated that the CFD model developed in this work is capable of reproducing with great accuracy experimental data of hydrogen concentration following a laminar release in a confined environment without ventilation. This study represents the first step in building a model for hydrogen releases capable of replicating a wide range of operating conditions. The long-term goal of this work is to apply this model in the risk analysis of an accidental hydrogen release within a glass manufacturing plant. In such a scenario, the release is likely to originate from a hole in a low-pressure pipeline from which the fuel is generally supplied to the burners. In this application, a low- velocity hydrogen flow with a constant flow rate can be expected. The next step for this work is to validate experiments conducted in larger domains, in the presence of obstacles and mechanical ventilation, in order to represent all the features of a potential release inside a glass manufacturing facility. Furthermore, the model can be used to explore different ventilation rates with the aim of suggesting the optimal configuration to prevent the formation of hydrogen explosive atmospheres. This information contributes to deepening the knowledge of 431 hydrogen safety and supports the development of adequate safety measures, promoting the deployment of hydrogen technologies in the glass industry. 5. Conclusions At present, there is scientific evidence that replacing natural gas with hydrogen in glass furnaces is one of the most effective decarbonisation options. A review of existing literature on this topic revealed that very few studies have conducted simulations in large industrial environments, and most of these involved high-pressure releases from HFCVs. In this work, a CFD model was developed and validated to simulate a release scenario typical of the glass industry. The simulation results demonstrated an excellent correspondence with the experimental data, despite a small underestimation (10%) of hydrogen diffusion. Therefore, the developed model is capable of accurately reproducing the transient concentration profile of hydrogen in a confined space following a laminar release and can be employed to conduct consequence analyses (e.g., determine the mass of hydrogen within the flammability zone and calculate the overpressure in case of ignition). The next steps in developing this model will involve a larger domain, the presence of obstacles, higher release velocities, and mechanical ventilation. References Ansys Fluent, 2023, Fluid Simulation Software accessed 01.11.2023. Bauwens C., Dorofeev S., 2014, CFD modeling and consequence analysis of an accidental hydrogen release in a large scale facility, International Journal of Hydrogen Energy, 39.35, 20447–20454. CINEA, 2022, How LIFE is reducing emissions from glass production, European Climate, Infrastructure and Environment Executive Agency accessed 15.09.2023. Ecofys, 2009, Methodology for the free allocation of emission allowances in the EU ETS post 2012 accessed 15.09.2023. Eurostat, 2023, Air emissions accounts for greenhouse gases by NACE accessed 15.09.2023. Furszyfer Del Rio D. D., Sovacool B. K., Foley A. M., Griffiths S., Bazilian M., Kim J., Rooney D., 2022, Decarbonizing the glass industry: A critical and systematic review of developments, sociotechnical systems and policy options, Renewable and Sustainable Energy Reviews, 155, 111885. H2GLASS, 2023, Decarbonising the glass industry with hydrogen technologies accessed 15.09.2023. Houf W., Evans G., Ekoto I., Merilo E., Groethe M., 2013, Hydrogen fuelcell forklift vehicle releases in enclosed spaces, International Journal of Hydrogen Energy, 38.19, 8179–8189. Hussein H., Brennan S., Molkov V., 2020, Dispersion of hydrogen release in a naturally ventilated covered car park, International Journal of Hydrogen Energy, 45.43, 23882–23897. Jatzwauk C., 2021, Design and Operation of Glass Furnaces, Chapter in: Richet P., Conradt R., Takada A. and Dyon J. (Ed.), Encyclopedia of Glass Science, Technology, History, and Culture, Vol II, Wiley, Hoboken, NJ, USA, 1147–1164. Li Y., Xiao J., Zhang H., Breitung W., Travis J., Kuznetsov M., Jordan T., 2021, Numerical analysis of hydrogen release, dispersion and combustion in a tunnel with fuel cell vehicles using all-speed CFD code GASFLOW- MPI, International Journal of Hydrogen Energy, 46.23, 12474–12486. Lucas M., Skjold T., Hisken H., 2020, Computational fluid dynamics simulations of hydrogen releases and vented deflagrations in large enclosures, Journal of Loss Prevention in the Process Industries, 63, 103999. Pardo N., Moya J.A., Vatopoulos K., 2012, Prospective Scenarios on Energy Efficiency and CO2 Emissions in the EU Iron & Steel Industry, Publications Office of the European Union. Patel P., Baalisampang T., Arzaghi E., Garaniya V., Abbassi R., Salehi F., 2023, Computational analysis of the hydrogen dispersion in semiconfined spaces, Process Safety and Environmental Protection, 176, 475–488. Schmitz A., Kamiński J., Maria Scalet B., Soria A., 2011, Energy consumption and CO2 emissions of the European glass industry, Energy Policy, 39.1, 142–155. Scott D. S., Cox K. E., 1960, Temperature dependence of the binary diffusion coefficient of gases, The Canadian Journal of Chemical Engineering, 38.6, 201–205. Yang F., Wang T., Deng X., Dang J., Huang Z., Hu S., Li Y., Ouyang M., 2021, Review on hydrogen safety issues: Incident statistics, hydrogen diffusion, and detonation process, International Journal of Hydrogen Energy, 46.61, 31467–31488. Zhang X., Wang Q., Hou X., Li Y., Miao Y., Li K., Zhang L., 2021, Numerical Analysis of the Hydrogen Dispersion Behavior in Different Directions in a Naturally Ventilated Space, Applied Sciences, 11.2, 615. Zier M., Pflugradt N., Stenzel P., Kotzur L., Stolten D., 2023, Industrial decarbonization pathways: The example of the German glass industry, Energy Conversion and Management, X.17, 100336. 432 57schiaroli.pdf Analysis of Accidental Hydrogen Releases in the Glass Manufacturing Industry