Fluid structure interactions modelling in Vented lean deflagrations

The standard 20ft ISO containers are studied both experimentally and numerically with model obstacles to ascertain the peak overpressures generated in case of an accidental fast deflagrations. Apart from overpressure its often important to know the maximum deflections of the enclosures to examine the structural integrity. The container walls are not rigid; they not only contribute through acoustic response to the pressure waves but also through structural resonance response to the generated overpressures. Hence to capture these effects its necessary to do the fluid structure interaction (FSI) or in simple terms the coupled Computational Fluid Dynamics (CFD) and Finite element analysis (FE) simulations. Alternatively, the structural response can be coupled to CFD simulation in either one-way or two-way interactions to reduce the associated computation efforts. Numerical simulations have been conducted with a pseudo two-way coupled CFD & FE, to aid our understanding of the combination of various factor contributing for the generation of overpressures using an in-house solver developed based on open source CFD code OpenFOAM named as HyFOAM. The CFD solver solves the compressible Navier-stokes equations along with the coupled fluid structure interaction model equations. The turbulent flame deflagrations are modelled using the flame area combustion model, the necessary modification are incorporated to governing equations to include the dominant flame instabilities present during the vented deflagration process. For lean hydrogen mixtures the Lewis number effects are important hence suitable modelling improvements are added to the combustion model. The numerical results are first validated against the recent experiments conducted at Gexcon, Norway as part of the HySEA project supported by the Fuel Cells and Hydrogen 2 Joint Undertaking (FCH 2 JU) under the Horizon 2020 Framework Programme for Research and Innovation. Numerical predictions of CFD and FE coupled, are assessed against the experimental results to study the contributing factors affecting the generated overpressure with in the ISO containers in view to overall improve the numerical predictions.


Introduction
Hydrogen is a clean and green energy carrier for the future considering currents environmental and health problems associated with the pollutants emitted by the fossil fuels usage. Hydrogen gas handling and storage units typically deals with either pressurizing or depressuring the gas and this equipment's are generally housed in a container like enclosures. The standard shipping ISO containers are also being considered for housing standalone portable power generation units using the latest fuel cell technologies. From process safety point of view knowing the hazards associated with these kind of container applications is of very crucial importance for overall design and operations of hydrogen installations.
Venting the gas during an accidental explosion in an enclosure is one of the simple mitigation techniques to reduce the generated overpressures inside the enclosure and retains the structural integrity of the process equipment with in it. The main design parameters are vent size and the location of the vents. To understand the factors that effects vented explosion process, many experimental and numerical studies have been performed in the last few decades Mccann (1985), Copper (1986), Molkov (1993), Bauwens (2008). The venting of lean hydrogen mixtures is studied to further access the influence of the hydrogen gas properties in the venting process Bauwens (2008), Molkov (1999), Bauwens (2011), Bauwens (2014). The Lewis number effects are important in the lean hydrogen-air combustion processes Abdel (1985), Chakraborty (2011). Hence suitable modelling correction should be done to take into account the Lewis number effects into turbulent flame deflagration modelling. Flame instabilities are also found to have dominating effects on the venting process of the lean hydrogen mixtures, Bauwens (2008), Molkov (1999). The Darrieus-Landau and thermodiffusive instabilities also affects the flame propagation in lean hydrogen mixtures. While venting of hot gases occurs through vent areas, Helmholtz oscillations are generated within the vessel due to expulsion of bulk of the gas and the resulting flame accelerations give rise to Rayleigh-Taylor instabilities of the flame front. These instabilities are studied in numerical modelling either using an algebraic expression or solving a transport equation Bauwens (2008), Molkov (1999). The various factor mentioned above present in the turbulent lean hydrogen deflagrations but they influence are dominant at different stages of the flame development. The Darrieus-Landau and thermodiffusive instabilities are found to be dominant in the initial stage of flame propagation after the start of the ignition process. The Helmholtz oscillations and Rayleigh-Taylor instabilities tend to contribute in the later stage venting process, when the hot gases start to exit from the vent area. The hydrogen deflagration studies in ISO containers are very scarce. Sommersel (2017) reported explosion experiments in inhomogeneous hydrogen air clouds in a standard 20′ ISO container. Test parameter variations included nozzle configuration, jet direction, and reservoir back pressure, time of ignition after release and degree of obstacles. The paper presents the experimental setup, resulting pressure records and high speed videos. The explosion pressures from the experiments without obstacles were in the range of 0.4-7 kPa. In the experiments with obstacles the gas exploded more violently producing pressures in order of 100 kPa. The CFD and FE coupled solution are more frequent in wind energy technologies development studies Lin (2016). The study of the steel structures in natural fires was attempted by employing the CFD and FE coupled numerical methods by Michał (2017). In their work coupling is done between CFD and FE by dedicated scripts, which utilize developed methods and compute the heat transfer between gas and solid phase. The emphasis is put on the proper calculation of temperature field inside the structural members and, in particular, the non-uniform temperature distribution inside the section. Coming to the Vented explosion applications, the CFD and FE coupled numerical methods still rare Salaün (2016) as major focus of the vented explosion studies is been in quantifying the peak overpressures with in the enclosures treating the enclosed as a rigid body.

Numerical Modelling
The open source Computational Fluid Dynamics (CFD) code OpenFOAM libraries are been used to develop the HyFOAM solver for vented lean hydrogen deflagration processes. The flow governing Navier-Stokes equations are solved in Large Eddy Simulation (LES) approach with a collocated, finite-volume method. The pressure-velocity coupling in the momentum equation is solved in a fully compressible Pressure-Implicit Split Operator (PISO) solution method. The Diffusion terms are discretized using a second-order accurate central differencing scheme and the advective terms approximated using a second-order accurate limited linear scheme. The transient term was discretized using a fully implicit, second-order accurate three-time-level method, OpenFOAM (2015). A one equation eddy viscosity model is used for evaluating the subgrid scale (SGS) turbulence, Furby (1997). The main difficulty in LES is the proper treatment of the flame front or the reaction zone, since the characteristic scales for turbulent combustion are in the SGS, for which SGS reaction rate models are required. The Flame Surface Wrinkling Model developed by Weller (1998) is adopted for simulating the turbulent deflagrations. The set of governing equations are solved sequentially with iteration over the explicit coupling terms to obtain convergence. The segregate approach results in a Courant number restriction, Weller (1998). Courant number of 0.1 was used in the present numerical simulations. The following section deals with the description of the Flame Surface Wrinkling Model completing the governing equation for the reacting flows.

Combustion model
The flamelet concept simplifies the turbulent combustion treatment by separating the combustion modelling from the analysis of the turbulent flow field by assuming that reaction takes place in relatively thin layers that separate regions of unburned and fully burned gases. The laminar flamelet approach is used with conditional filtering to create a set of transport equations representing the complex combustion process, Weller (1998), Tabor (2004). The unburnt zone volume fraction is denoted as regress variable ( b ), taking values b = 1 in fresh gases and b = 0 in fully burnt gas region. The transport equation for the resolved part of regress variable (b) is given as: where,  is subgrid flame wrinkling, can be regarded as the turbulent to laminar flame speed ratio and is formally related to the flame surface density by || b     ,  is the density, L S is laminar flame speed and sgs  is the subgrid turbulent diffusion coefficient. Symbols ( ̅ ) and (̃) represent the filtered and the density weighted filtering operations respectively. The subscripts u indicates conditioning on the unburned gases region. The resolved unburned gas volume fraction b is related to b through u bb   . The closure for the sub-grid wrinkling is provided by a balanced transport equation, where, s U is the surface filtered local instantaneous velocity of the flame, which is modelled as The direction of flame propagation is || .
The terms G and ( 1) R  in equation 2 are sub-grid turbulence generation and removal rate, with G and R as rate coefficients requiring modelling. The modelling of these terms is based on flame-speed correlation of [14] are given below where, is the Kolmogorov time scale, û is the sub grid turbulence intensity and is the turbulent Reynolds number. The modelling of the terms in equation (5) for lean turbulent premixed combustion takes into account the Lewis number (Le) effects in the turbulent flame speed correlation The algebraic reaction rate closure, MFSD proposed in Muppala (2005) is adopted in the present study. This model has been successfully applied to both pure and mixed fuels, under varying Lewis number conditions Chakraborty (2011), Burke (2015), Muppala (2009), in both RANS and LES contexts. The MFSD model predictions for the turbulent flame speed (ST) for equivalence ratio between and inclusive of 0.4 and 0.8 along with Goulier's (2016) correlation is compared with the available experimental measured values of Kitagawa (2008). These results are shown in figure 1. The Muppala models (MFSD) analytical predictions of ST trends are consistent for the turbulence levels quantitatively to the experimental results in consideration.  Goulier (2016), Muppala (2005) with the experimental data of Kitagawa (2008).
The Darrieus-Landau and thermodiffusive instabilities also affects the flame propagation in lean mixtures. These instabilities are modelled considering the analytical expression proposed by Bauwens (2011). The flame wrinkling due to the Darreius-landau instabilities is modelled as algebraic expression based on Bauwens (2011) as, where, is cutoff wavelength of unstable scales and ∝ 1 is a coefficient to account for uncertainty in 1 . The unstrained laminar flame speed ( ,0 L S ) for lean hydrogen-air mixture is adopted based on the numerical study carried out by Verhelst (2003), for evaluating L S at a given equivalence ratio ( 1/   ) and reference condition, expressed as power law function of elevated temperature and pressure, where in cm/s , P is pressure in bar and unburnt gas temperature in K. The above correlation is valid for the equivalence ratios ( ) between 0.33 and 0.47 (lean mixtures), pressures range of 1bar ≤ P ≤ 8.5 bar and temperature range of 300 K ≤ T ≤ 800 K, with reference state Tu0 = 300 K. The flame wrinkling factor is equation (1) is evaluated as Thus, equations (1)-(8) describe the combustion model for the lean hydrogen mixtures in HyFOAM solver.

Experiments
The experiments are carried out at the Gexcon test site at Bergen, the details can be found in . The experimental rig consists of a standard 20-ft. ISO container with typical dimensions of 20' x 8' x 8'.6", fitted with a steel frame, shown in Figure 1. The container walls are around 2 mm thick. The purpose of the steel frame is to support obstacles, such as the basket holding the bottles, as well as the pressure probes inside the container. The frame as shown in Figure 1  The parameters that are measured during the tests are the deflagration overpressures and deflection of the container wall. The deflection are measured at the middle of the side walls of the containers (D1 and D2) using a laser displacement probe. A polythene sheet for used to cover the vent to retain the homogenous mixture at the start of the each experiment, which eventually got blown away during the explosion process.

CFD and FE
Finite Element method/Analysis (FE) is numerical method to solve partial differential equations. This method is well suited for the computational structural stress analysis hence became synonymous with the structural analysis. Computational fluid dynamics problems can be solved using Finite Difference method, Finite Element method and Finite Volume Method (FV), but the FV method is the most preferred method and Hence FV method employed in most of commercial and general purpose CFD solvers. Since the fluid dynamics and the structural stress analysis differ in their formulation and discretization methods, an interface modelling is required that will provide physical interaction between the two computational domain (fluid and solid regions). A fully coupled CFD and FE , will have two-way interaction i.e. the influence of the fluid forces on the solid structure and the displacements of the solid boundary on the fluid flow are solved simultaneously in their respective numerical methods, with a suitable interface program to inter transfer the required parameter values. In a one-way interaction, FE are solved independently, the tabulated boundary values from FE used in CFD, vis versa. In the present study the pseudo two way approach is being used, wherein the structural displacements against the overpressures from the experiments are used in CFD simulations. The structural displacements are applied to the container boundaries and scaled according to the generated peak overpressures values during the explosion process.   The initial rise of the overpressure excites the container wall, which sets the wall to vibrate, later oscillating in response mode. The maximum deflection of container wall is occurring during this resonance phase. The resonance oscillation of the container walls is also seems to contributing to the second overpressure peak in the pressure trace curve. This second overpressure peak is much more distinct and evident in the Figure 7 & 8 involving obstacles compared to the empty container cases (Test 1& 2).

Numerical setup
The standard 20-ft ISO container details used in the numerical simulation are shown in Figure  9. To validate the CFD and FE interactions, the scenarios experimental studied is considered for the numerical modelling, configuration of no obstacles, back wall ignition, steel frame inside the container and doors fully open with 15 % hydrogen concentration by volume, shown in Figure 9 (b). The ignition of the hydrogen mixture was initiated by creating a spherical hot patch at the center of the back end wall at the mid height of the container with products composition and temperature, mimicking the electric inductive spark used in the experiments. An hybrid hexagon-tetrahedral computational mesh was generated for the container geometry using the 'SnappyHexMesh' utility in OpenFOAM. The mesh distribution in the computation domain is shown in Figure 10. The volume enclosing the container, 30.0 m × 15.0 m × 35 m was also meshed to capture the venting of the burnt gas, the external explosions and to reduce the effect of boundary conditions on the numerical results. A non-uniform cell size of 0.5 cm was used in the ignition region, a 3 cm cell size inside the chamber and in the area immediately outside the chamber to resolve the external explosion. The total cells in computational mesh are approximately between 3.5 ~ 5.5 million for the three cases considered in the present study.

Figure 10. Computational domain and the mesh distribution in vertical cut plane
The boundary conditions applied to the container walls is moving wall, non-slip, adiabatic walls for the sides and roof walls except the ground, which is fixed with zero displacement. The

Results and discussion
The numerical predictions of overpressures for Case-1, are plotted in Figure 13 along with two experiment tests. These results are generated considering the CFD and FE pseudo two way coupling. Some of the pressure probes have not properly functioned during the experiments, but still whatever available information recorded at the plot are is plotted in the graphs. On each plot the two symmetric probe values ( P1&P2, P3&P4, P5&P6, P7&P8) of the same tests are been averaged to reduce the experimental trail variations. The pressure probe plots in 13 (a)-(d) are inside the container, 13 (e)-(f) are outside the container. The Figure 13 and Figure 14 are numerical predictions for Test-1, with and without the CFD & FE coupling. The numerical simulation results are written with a frequency of 10 ms. The coupling of CFD results with the FE (structural response) has resulted in a significant improvement in the overpressure trends. Figure 13(a) and 14(a), the location for the peak Figure 13. Case-1, pressure trace curves along with Test-1 and Test-2 experimental results, with CFD and FE in pseudo two-way coupling. Figure 14. Case-1, pressure trace curves along with Test-1 and Test-2 experimental results, without the CFD and FE pseudo two-way coupling.
Overpressure in the container for back wall ignition, showed a clear change in the pressure trace curves. The simulation results with CFD and FE coupled, reproduced the experimental trends of positive overpressures within the container. On the other hand, the numerical results without the CFD and FE coupling showed a region of strong negative pressure after the initial pressure peak. The pressure trace curves at the pressure probes outside the container are very similar in both scenarios of with and without CFD & FE coupling. Without the structural response modelling the container wall as rigid resulted in strong Helmholtz oscillations after the initial pressure peak (figure 14 (a)-(f)). The structural response in the present study are obtained from the experiments. But in absence of the experimental results these structural response can be obtained by using well established non-linear dynamics FE software's packages, such as LSDYNA, IMPETUS and so on. The present study shows the viability of the simple pseudo two-way CFD and FE coupling in improving the numerical prediction for vented lean deflagration in flexible structures.

Conclusions
The 20-ft ISO containers are being considered for developing self-contained portable power generation units using the latest fuel cell technologies. The possible scenarios of lean hydrogenair deflagrations in these containers are being studied numerically in the present safety study. The Experimental data from full scale container tests carried out by Gexcon, Bergen, are used for model validation. The vented lean deflagration requires the modification to flame speed correlation and turbulent flame speed along with the models for flame instabilities in the combustion model equations to simulate deflagrations in lean hydrogen-air mixtures. The first overpressure peak is well predicted and with discernable time lags in occurrences of pressure peaks in the later parts of the pressure trace in numerical results, wherein the structural response trends to influence the generated overpressures. Therefore, to further improve the overpressure trends, the CFD simulations are coupled to FE in a pseudo two-way interactions, where in the experimental structural responses are used to dynamically evaluate boundary wall displacements based on the peak overpressure inside the container. The final numerical results of the coupled CFD and FE are very promising in predicting the experimental trends within the experimental uncertainties. In future work, this same CFD and FE coupled solution methods will be used for evaluating overpressures for deflagration in containers having obstacles.