High Resolution 3-D Simulations of Venting in 18650 Lithium-Ion Cells

High temperature gases released through the safety vent of a lithium-ion cell during a thermal runaway event contain flammable components that, if ignited, can increase the risk of thermal runaway propagation to other cells in a multi-cell pack configuration. Computational fluid dynamics (CFD) simulations of flow through detailed geometric models of four vent-activated commercial 18650 lithium-ion cell caps were conducted using two turbulence modeling approaches: Reynolds-averaged Navier-Stokes (RANS) and scale-resolving simulations (SRS). The RANS method was compared with independent experiments of discharge coefficient through the cap across a range of pressure ratios and then used to investigate the ensemble-averaged flow field for the four caps. At high pressure ratios, choked flow occurs either at the current collector plate when flow through the current collector plate is more restrictive or the positive terminal vent holes when flow through the current collector plate is less restrictive. Turbulent mixing occurred within the vent cap assembly, in the jets emerging from the vent holes, and in recirculating zones directly above the vent cap assembly. The global maximum turbulent viscosity ratio ( μ T / μ ) of the MTI, LG MJ1, K2, and LG M36 caps at pressure ratio of P 1 /P 2 = 7 were 4,575, 3,360, 3,855, and 2,993, respectively. SRS and RANS simulations showed that both velocity magnitude and fluctuating velocity magnitude were lower for vent holes which are obstructed by the burst disk. SRS showed high levels of fluctuating velocity in the jets, up to 48.5% of the global maximum velocity. The present CFD models and the resulting insights provide the groundwork for future studies to investigate how jet structure and turbulence levels influence combustion and heat transfer in propagating thermal runaway scenarios.


INTRODUCTION
Lithium-ion batteries (LIB) are a prevailing energy storage device widely used in mobile applications like portable electronic devices, electrical vehicles and aircraft because of their high energy and power density Wang et al., 2019). However, safety remains a challenge for LIBs, especially with increasing power and energy density at the module-and system-level. The fire hazard of LIBs is driven by the risk of thermal runaway, during which a large quantity of heat and gas is produced by a series of exothermal reactions, including: decomposition of the solid electrolyte interphase (SEI) layer, reaction of anode with electrolyte, cathode decomposition, electrolyte decomposition, and others (Kim et al., 2007). The high level of heat and gas generation rates during thermal runaway cause the temperature and pressure inside the cell increase, nearly instantaneously. Common 18650 format LIB cells have a safety vent mechanism built into the positive terminal cap and when the internal pressure of the cell increases to the vent activation pressure, the safety vent opens releasing flammable and toxic gases as well as electrolyte vapors (Nedjalkov et al., 2016). The vented flammable gas may ignite, causing fire and increasing risk of thermal runaway propagation within a multi-cell module or system (Wang et al., 2012). Developing methods, either numerical or experimental, to investigate the flow field during venting can help provide new insights into the combustion and heat transfer during propagating failures.
The following information is needed for development of models which are capable of simulating a full thermal runaway, with effects of combustion and heat transfer resulting from vented gases included: decomposition reaction sequence, kinetic model and parameters for each reaction, heat and gas generation amounts associated with each reaction, and composition of gases emerging from the cell during venting. Significant work has been accomplished in developing the reaction sequence, kinetic model and parameters, and heat generation amounts. The classical paper by Hatchard et al. laid the groundwork for thermal abuse models of Li-ion cells using three decomposition reactions: SEI decomposition, reaction of anode and electrolyte, and cathode decomposition (Hatchard et al., 2001). Since then, many studies have built on the original model framework to account for additional reactions (Kim et al., 2007;Ren et al., 2018;Xu et al., 2018), to consider different cell chemistries (Peng and Jiang, 2016;Dong et al., 2018;Bugryniec et al., 2020), to account for aging (Prada et al., 2012;Abada et al., 2018;Larsson et al., 2018), to investigate different thermal management strategies (Chen et al., 2016;Lopez et al., 2016;Bugryniec et al., 2020), and to account for gas generation from decomposition reactions (Coman et al., 2017;Ostanek et al., 2020).
Relatively few studies have considered the effects of venting in numerical thermal abuse models. Coman et al. (Coman et al., 2016) developed a thermal abuse model of an 18650 cylindrical cell which considered the effects of electrolyte evaporation and venting. The authors reported that time-to-runaway is not captured accurately when endothermic effects such as evaporation of electrolyte and enthalpy carried away from the cell by ejecta are neglected. Coman et al. (Coman et al., 2017) then included gas generation from the SEI decomposition reaction, which contributes to pressure buildup in the cell along with evaporation of electrolyte, and found that time-to-venting was captured more accurately. The evaporation model of Coman et al. (Coman et al., 2017) assumed the electrolyte is in a state of vaporliquid-equilibrium (VLE) prior to venting. Our previous work (Ostanek et al., 2020;Parhizi et al., 2021) built upon the lumped model approach and included a model for non-VLE evaporation after venting using a constant-rate and decaying-rate drying model. For oven tests and other slow heating abuse tests, the cell temperature decreases after vent-activation because of the heat absorbed during non-VLE evaporation. Mao et al. (Mao et al., 2021) considered the reacting flow after venting and developed a lumped model to estimate the flame height for flow emerging from an Nickel Cobalt Manganese (NCM) chemistry 18650 cell. The model assumes that the safety vent is a circular orifice and calculates the vent gas velocity assuming the isentropic flow rate.
The venting of gases emerging from Li-ion cells, and subsequent combustion, has important effect on the thermal runaway propagation within a module or system. Combustion of vented gases increase heat transfer to neighboring cells by impingement and thermal radiation. Walker et al. (Walker et al., 2019) investigated the fraction of the energy released through cell case and the ejecta and vented gases during thermal runaway. The results showed that the energy released through the ejecta and gases was larger than that released by conduction through the cell case for four types of Li-ion caps (Samsung 18650-3Q, LG 18650-MJ1, 3.35 Ah LG 18650, and MOLiCEL 18650-J). Using LG 18650 MJ1 as an example, 78.3% of the total energy released in thermal runaway (74.9 kJ) was transferred by the ejecta and vented gases. Zhong et al. (Zhong et al., 2018) studied the thermal runaway propagation in an 18650 LIB module with 3 × 3 cells. If the cells were displaced with 4 mm distance to each other, the thermal runaway can still propagate to the surrounding cells by the thermal shock created by the flame of the vented gases. Archibald et al. (Archibald et al., 2020) found that the thermal runaway propagation in the 10 × 10 Ah module was faster than that in 10 × 5 Ah module. The larger thermal runaway propagation speed was due to preheating of the unfailed cells by the high temperature vented gas generated from failed cells. The higher chamber temperature produced by the vented gas also decreased the heat loss of cells. Cheng et al. (Cheng et al., 2019) used video analysis to study the thermal runaway procedure of a LIB module with 12 commercial NCM prismatic cells. The thermal runaway propagation was attributed to heat conduction between cells as well as convection and radiation caused by combustion of the vented gas. The flames ignited the wire harness of cells on the other side of the module and initiated thermal runaway in those cells. Mishra et al. (Mishra et al., 2021) investigated the impact of vent gas flow on thermal runaway propagation in a 25-cell module. The primary heat transfer mechanism considered in this study was the flow of hot gases from a failed cell impinging on neighboring cells and nearby surfaces. The cell-to-cell spacing, overhead gap spacing, and location of the vent hole had significant effect on thermal runaway propagation.
Although venting and combustion play an important role in module-level protection, there is a substantial lack of literature on numerical modeling of the external flow phenomena. Existing studies which consider venting and combustion typically simplify the vent geometry to a circular orifice which may not capture the spatial distribution and velocity fields of vented gases within the module. The present work investigates the flow field during venting of an 18650 LIB cell using a three-dimensional (3-D) CFD model based on highly detailed cap geometries. The ensemble-averaged flow field was simulated using Reynoldsaveraged Navier-Stokes (RANS) simulations. Model benchmarking and a mesh sensitivity analysis were conducted for a sharp-edge orifice. From a series of RANS simulations across a range of pressure ratios, sharp-edge orifice equivalent areas were calculated for each battery cap and compared with experimental data. The restriction imposed by each vent cap was quantified using the equivalent sharp-edge orifice flow area. The unsteady flow field was then analyzed using scale-resolving simulations (SRS). Mean and fluctuating velocity distributions obtained from SRS were investigated and the turbulent structures emerging from the vent were visualized.
where C 2 is constant, , and σ ε is the turbulent Prandtl number for turbulence dissipation rate.
Four commercial LIB vent cap assemblies were considered in this paper: MTI, LG MJ1, K2, and LG M36. All four vent caps are used in 18650 format cylindrical Li-ion cells (18 mm diameter, 65 mm length). The MTI vent cap assembly was purchased as an individual component from the manufacturer, while the remaining 3 cells were extracted from live Li-ion cells following the procedure of our previous work . The LG MJ1 cap was extracted from an LG INR18650 MJ1 cell, the K2 cap was extracted from a K2 LFP18650E-1500-03 cell, and the LG M36 cap was extracted from an LG INR18650 M36 cell. The nominal voltage of the LG INR18650 MJ1 cell is 3.6 V, the nominal capacity is 3.5 Ah, and the electrodes are LiNiMnCoO 2 -graphite. The nominal voltage of the K2 LFP18650E-1500-03 cell is 3.2 V, the nominal capacity is 1.5 Ah, and the electrodes are LiFePO 4 -graphite. The nominal voltage of the LG INR18650M36 cell is 3.63 V, the nominal capacity is 3.6 Ah, and the electrodes are LiNiMnCoO 2graphite.
After extracting the vent cap assemblies from the commercial cells, computed tomography (CT) scans were conducted following the procedure of our previous work . From the CT scans, which are publicly available for download (Stoll et al., 2020b;a;2021a;, detailed 3-D geometric models of the caps in a vent activated state were developed for use in CFD simulations. An isometric view of the geometric model for each commercial LIB vent cap assembly is shown in Figure 1. The CFD model was developed using the detailed vent cap geometric models shown in Figure 1. The computational domain, boundary conditions, and mesh used for 3-D simulations of flow through the vent caps are shown in Figure 2. The computational domain was designed to emulate an experimental setup that was developed previously to measure discharge coefficients of commercial LIB caps . The apparatus uses a flat plate with a rounded inlet which holds the LIB cap specimen. The flat plate also serves as a divider between upstream and downstream chambers. The pressure ratio across the cap is set by adjusting the pressure in the two chambers. The computational domain in Figure 2A emulates the experimental apparatus but consists of only the fluid region. In other words, the solid parts (LIB cap and separator plate) were subtracted from the full domain, leaving just the fluid region as shown in Figure 2B. A hemispherical surface having radius of 180 mm was assigned a pressure outlet boundary condition. Impermeable, no-slip, and adiabatic wall boundary conditions were assigned to the surfaces of the LIB cap and the fixture plate. A smaller hemispherical surface having radius of 45 mm was assigned as a pressure inlet boundary condition.
The working fluid was assumed to be single-phase, dry air. The total pressure at the inlet was varied such that P 1 /P 2 ranged from 1.2 to 8. Total temperature at the inlet was T 1 300 K. The outlet total pressure was set to P 2 101.3 kPa. It was assumed the incoming flow was straightened and conditioned, with a turbulence intensity of 0.5% and turbulent length scale of 5 mm.
Two different unstructured meshing approaches were used to discretize the computational domain: tetrahedral elements and polyhedral elements. Tetrahedral meshes were generated using the meshing utility with ANSYS Workbench (ANSYS, 2019c). Polyhedral meshes were obtained by conversion of the tetrahedral meshes using the polyhedral conversion utility within ANSYS Fluent (ANSYS, 2019b). A mesh sensitivity was conducted using coarse, medium, and fine meshes for both tetrahedral and polyhedral meshes. Mesh constraints included a global element size of: 5 mm for the coarse mesh, 3.5 mm for the medium mesh, and 2.5 mm for the fine mesh. All meshes were generated using: global element growth rate of 1.2, edge proximity capture with a minimum element size of 80 μm, curvature capture with a minimum element size of 100 μm, and wall-adjacent inflation layers having five layers with a transition ratio of 0.272. Examples of fine meshes using tetrahedral and polyhedral elements are shown in Figures Figure 2C shows the plane which is normal to the z-axis and Figure 2D shows the plane normal to the x-axis.

SRS CFD Model for LG MJ1 Cap
The scale-resolved simulations (SRS) technique, specifically the stress-blended eddy simulation (SBES) turbulence model, was used to simulate unsteady flow through the LG MJ1 cap. The geometric model and computational domain were modified from the domain used in RANS simulations. Instead of a domain which replicated the experimental setup for discharge coefficient measurement, a new domain was created which includes an 18650 cell surrounded by an air space, as shown in Figure 3A.
The element sizing constraints within the vent cap were the same as that of the fine polyhedral mesh described in section 2.1. The mesh density outside the cell near the vent holes was refined using conical regions placed above the vent holes. The mesh density was adjusted by using three different conical regions with progressively increasing element size, 0.2 mm to 0.5 mm to 0.8 mm, moving from the inner to outer cones, respectively. The mesh density in the three conical regions is shown in Figure 3B. The pressure inlet boundary condition was set in the plane which would be located below the vent cap, as shown in Figure 3A. The outer surfaces of the cylindrical domain of the fluid were assigned pressure outlet boundaries to represent a single cell sitting in a quiescent environment. The wall surfaces of the cap and cell case were assigned no-slip, adiabatic, and impermeable boundary conditions. The time step for the simulation was set to 1E-6 s. The same solver settings in section 2.1 were used: pressure-based coupled algorithm, temperature-dependent polynomial functions for transport properties, and ideal gas equation of state. Advection terms in the momentum equation were discretized using the default setting for the SBES model: bounded, second-order differencing scheme (ANSYS, 2019b).

Sharp-Edge Orifice Benchmarking
The CFD procedures described in section 2.1 were implemented for a benchmark case of compressed air flowing through a sharpedge orifice. The model setup and results of the benchmarking case are shown in the Supplementary Material. A mesh sensitivity analysis was conducted for the benchmarking case, using the same global element sizing and inflation layer settings as those used for the LIB cap simulations. Results of the benchmarking cases agreed with literature data and, thus, the CFD solution procedure was considered suitable for the LIB cap simulations.

RANS Simulations for LIB Caps
A series of RANS CFD simulations were conducted to analyze the ensemble-averaged flow field emerging from LIB vent caps under a range of pressure ratios. First, a mesh sensitivity analysis was conducted for the RANS CFD model. Next, the discharge coefficients for the LIB caps were calculated over a range of pressure ratios. The discharge coefficients were compared to previous experimental data and also used to compare the flow restriction imposed by each vent cap. Finally, details of the simulated flow field, including Mach number and turbulent viscosity ratio distributions, were analyzed and compared for each of the four battery caps.

Mesh Sensitivity of Mach Number Distribution Analysis
A mesh sensitivity analysis was conducted for the LG MJ1 vent cap. The sensitivity analysis was conducted for the coarse, medium and fine meshes for the two element types, tetrahedral and polyhedral. For tetrahedral elements the coarse mesh had ∼4.93E6 elements, the medium mesh had ∼7.69E6 elements, and the fine mesh had ∼12.6E6 elements. The polyhedral mesh sizes were: coarse ∼2.35E6 elements, medium ∼3.14E6 elements, and fine ∼4.14E6 elements.
The velocity magnitude, represented by Mach number, was extracted along two different lines for a flow simulation with a pressure ratio of P 1 /P 2 7: one aligned with the jet and one extending across one of the vent holes in the horizontal (transverse) direction. Along the jet axis, as shown in the inset of Figure 4A, the Mach number distribution showed little sensitivity to the mesh refinement level. In particular, the simulation results for the three polyhedral meshes showed very little change between the three mesh refinements. In Figure 4B, the Mach number distribution also showed little sensitivity to the mesh refinement level. At this pressure ratio, the maximum Mach number was just over 2.0.
In addition to comparing the local Mach number distribution, the mass flow rate through the cap was monitored for the P 1 /P 2 7 flow simulation. The simulated mass flow rate obtained with coarse, medium, and fine meshes using tetrahedral elements had less than 0.4% variation. The simulated mass flow rate obtained with coarse, medium, and fine meshes using polyhedral elements had less than 0.1% variation. The difference in the mass flow rate between the tetrahedral and polyhedral meshes was 0.77% for the coarse mesh, 0.62% for the medium mesh, and 0.97% for the fine mesh.
Polyhedral elements have the following advantages relative to tetrahedral elements: improved mesh quality, faster convergence, and lower overall mesh count (ANSYS, 2019b). Considering these advantages, and the less than 0.97% difference in flow rate with the defined mesh sizing parameters, polyhedral elements with the fine mesh density were used for discharge coefficient comparison and flow field analysis for all four caps. Using the fine polyhedral

Sharp-Edge Orifice Equivalent Area Analysis: Models Compared to Experiments
The flow restriction of each cap was quantified by calculating a sharp-edge orifice equivalent area by a method of least squares fitting of discharge coefficients to those of a sharp-edge orifice. The discharge coefficient is the ratio of actual mass flow rate _ m through the restriction divided by the ideal, isentropic flow through the restriction, _ m ideal , , as shown in Eq. 9.
The minimum flow area required to calculate _ m ideal is difficult to determine for the LIB cap due to the complex geometry of the current collector plate, positive terminal plates, and deformed burst disk . The area appears in Eq. 10, for the ideal mass flow rate through a restriction: where A is the flow area of the restriction, P 1 is the upstream stagnation pressure, P 2 is the downstream stagnation pressure, R is the gas constant for air, T 1 is the upstream stagnation temperature, and c is the ratio of specific heats.
In previous work, the problem was addressed by defining an equivalent flow area, A eq , which is the area that results in the LIB cap discharge coefficients matching that predicted by a semi-empirical model for flow through a sharp-edged orifice via a least squares fitting process . Once the equivalent area is identified, either through CFD or experimentation, it provides a quantitative comparison of the relative flow restriction of different battery caps .
The discharge coefficients of the CFD models, as well as those from previously reported experimental data , of all four commercial caps after the sharp-edge orifice equivalent area least squares fitting process are shown in Figures 5A-D. A semi-empirical model for flow through a sharp-edge orifice is also shown (Jobson, 1955 The equivalent areas obtained from CFD simulations were compared with those obtained from experiments. For the MTI cap, the equivalent area obtained through simulations was A eq 8.06 mm 2 while the experimentally-identified area at 95% confidence was A eq 8.14 ± 0.337 mm 2 . For the LG MJ1 cap, the equivalent area obtained through simulations was A eq 8.99 mm 2 while the experimentally-identified area at 95% confidence was A eq 8.62 ± 0.941 mm 2 . For the K2 cap, the equivalent area obtained through simulations was A eq 7.90 mm 2 while the experimentally-identified area at 95% confidence was A eq 8.50 ± 0.304 mm 2 . For the LG M36 cap, the equivalent area obtained through simulations was A eq 13.7 mm 2 while the experimentally-identified area at 95% confidence was A eq 11.28 ± 0.196 mm 2 . Discrepancy between the measured and simulated equivalent areas may be caused by variation in the position and shape of the deformed burst disk. The equivalent areas from, both CFD simulations and experiments showed that the LG M36 cap had the least restriction since its equivalent area was significantly larger than the other three caps. This larger equivalent area is attributed to the large openings that can allow the burst disk to deflect upwards to a greater extent than the other cap designs, resulting in a larger effective area. The MTI, K2, and LG MJ1 all had similar equivalent areas which is attributed to similar size vent holes in the cap.

Ensemble-Averaged Flow Field Analysis
RANS simulations provide the ensemble-averaged flow field which is equivalent to the time-averaged flow field when the process is stationary (i.e. boundary conditions are not changing). The jet Reynolds number was calculated as: where ρ is air density, μ is molecular viscosity, D is the hydraulic diameter of each individual hole/opening in the positive terminal, and U is the average velocity through each opening. Note that temperature at the jet exit was calculated from isentropic relations and then used in calculating density and molecular viscosity. The average velocity through each opening was calculated using: where _ m is the total mass flow rate through the vent cap, N is the number of openings in the positive terminal, and A is the flow area of a single opening. For a pressure ratio of P 1 /P 2 7, the resulting Reynolds numbers for the MTI, LG MJ1, K2, and LG M36 caps were 1.06E5, 1.12E5, 7.75E4, and 1.27E5, respectively. The large Reynolds numbers indicate fully turbulent jet flow for each cap (Lee and Chu, 2003).
Contours of velocity magnitude, represented as Mach number, on two orthogonal cross-section planes for each cap are shown in Figure 6 for a pressure ratio of P 1 /P 2 7. Mach number was calculated as: In Eq. 13, U is the instantaneous velocity magnitude and c is the speed of sound. The speed of sound for an ideal gas is: where c is the specific heat capacity ratio, R is the specific gas constant, and T is temperature. Figure 6A, B show the Mach number contours for the MTI cap in the x-normal and z-normal views, respectively. The maximum Mach number observed in the planes shown in    Figure 7 shows a representation of the jet trajectories emerging from each cap for a pressure ratio of P 1 /P 2 7 by Mach number isosurfaces of Ma 0.1 and Ma 0.3. The jets emerging from the MTI cap, Figure 7A, and K2 cap, Figure 7C, were relatively symmetric compared to the jets emerging from the LG MJ1, Figure 7B, and LG M36 cap, Figure 7D. The jet emerging in the negative z-direction of the LG MJ1 cap in Figure 7B was partially blocked by the presence of the burst disk. No jet emerged in the negative z-direction for the LG M36 in Figure 7D due to asymmetry of the three vent slots in the positive terminal. Thus, the flow of gases through the safety vents may be asymmetrical depending on multiple aspects of the cap design, including the burst disk and positive terminal.
The Mach number distribution in the x-normal viewing plane for the P 1 /P 2 7 pressure ratio was analyzed further by inserting a streamline which traces the flow originating at the inlet boundary and passes through the hole on the right side of the cap and then   A choke point in compressible flow systems occurs when the pressure ratio across a restriction is large enough to accelerate the flow to Ma 1. When choked flow occurs, the mass flow rate through the restriction is independent of the pressure downstream of the restriction, as shown in Eq. 2. The acceleration of flow through the vent caps plays an important role on the overall flow restriction of the cap. For example, the MTI and K2 caps both experience choked flow at the current collector and increasing the area of the positive terminal will not affect the flow rate. For example, the K2 cap has five holes in the positive terminal with an opening area of 20.6 mm 2 while the MTI cap has four holes with an opening area of 18.5 mm 2 . Despite the 10.7% difference in flow area at the positive terminal cap, the two caps have only a 3% difference in the experimentally-measured equivalent area .
The absolute pressure extracted along the streamlines in Figure 8 was consistent with observed choke plane locations for the P 1 /P 2 7 pressure ratio. For the MTI and K2 caps, a large decrease in pressure was observed across the current collector (around y/D −0.2) where the flow accelerated to Ma 1. The reduction in pressure exceeded the critical pressure ratio, which is a condition required for choked flow to occur. For those two caps, the pressure remains nearly constant between y/D −0.2 to −0.1. For y/D > −0.1, the pressure decreased to the ambient pressure level as the flow expanded out into the surrounding fluid. For the LG MJ1 and LG M36 caps, the pressure decreased monotonically through the cap. Figure 9 shows the turbulent viscosity ratio, μ T /μ, distribution on two orthogonal cross-section planes for each cap at a pressure ratio of P 1 /P 2 7. In-plane streamlines are also overlaid on Figure 9 to illustrate secondary flow patterns. The left column of images in Figure 9 shows the x-normal view and the right column shows the z-normal view.
For the k-ε family of turbulence models (ANSYS, 2019a), turbulent viscosity is calculated as: where ρ is the fluid density, C μ is a model constant, k is turbulent kinetic energy, and ε is the turbulent dissipation rate. Turbulent viscosity is added to the molecular viscosity that appears in the momentum transport equations, effectively increasing the transport of momentum due to the presence of turbulence (ANSYS, 2019a). The turbulent viscosity ratio distribution for the MTI cap, Figures 9A,B, shows a large degree of turbulent mixing in three regions: between the current collector and the positive terminal cap, along the jets, and in the region between the jets directly above the vent cap assembly. Turbulence is generated between the current collector and positive terminal cap by two mechanisms: from wall shear as flow passes through the current collector and from flow recirculation (as indicated in the figure).
Turbulence is generated along the jets from the shear force generated between the high momentum flow emerging from the vent holes and the low momentum surrounding fluid. Additionally, the steep angle of the jets emerging from the MTI cap causes a recirculation region to form directly above the vent cap, leading to increased turbulent mixing. Surrounding fluid is also entrained in the jets, as indicated in Figures 9A,B.
The turbulent viscosity distribution for the LG MJ1 cap, Figures 9C,D, shows a large degree of turbulent mixing in two regions: between the current collector and positive terminal and along the jets. Streamlines in Figures 9C,D show the shallow angle of the jets emerging from the LG MJ1 cap. There was no recirculation zone above the cap and no significant turbulent mixing in that region. Additionally, the current collector of the LG MJ1 was less restrictive (simulated A eq 8.99 mm 2 ) than the MTI current collector (A eq 8.06 mm 2 ). The MTI cap has six small holes around the periphery and small holes in the middle of the current collector from spot-weld tearing during CID-activation . On the other hand, the LG MJ1 current collector has three large slots around the periphery and one large opening in the middle of the current collector from CID-activation . The geometry of the LG MJ1 cap results in less turbulent mixing from shear and recirculation between the current collector and positive terminal. Generally, turbulent viscosity ratio was lower for the LG MJ1 cap relative to the MTI cap. The LG MJ1 cap showed a global maximum of μ T /μ 3,360 and the MTI cap showed a global maximum of μ T /μ 4,757.
The turbulent viscosity ratio distribution for the K2 cap shown Figures 9E,F was similar to that of the MTI cap which is attributed to the similar geometry of the two caps. The global maximum turbulent viscosity ratio for the K2 cap was μ T /μ 3,855. Likewise, the turbulent viscosity ratio distribution for the LG M36 cap, Figures 9G,H, was similar to that of the LG MJ1 cap due to similar vent cap geometry. The global maximum turbulent viscosity ratio for the LG M36 cap was μ T /μ 2,993. Overall the results show that different cell cap geometries have different degrees and spatial distributions of turbulent mixing with gases flowing through them. These flow aspects will influence combustion of vented flammable gases in a simulation of thermal runaway event (Turns, 2011), and thus can have a significant impact on the predicted heat transfer to the cell's surroundings.

Scale-Resolved Simulations of Flow Through LG MJ1 Cap
Scale-resolved simulations (SRS) were conducted on the LG MJ1 cap to analyze the unsteady flow field and dynamic evolution of turbulent structures. The SRS approach is, generally, more expensive but more accurate than RANS turbulence models. One reason that SRS provide improved accuracy is because large scale turbulent eddies, which account for most of the turbulent transport, tend to be anisotropic (ANSYS, 2019a). The anisotropy may be captured in certain RANS models, such as the Reynolds-Stress model, but this model requires additional transport equations and is more computationally expensive than the one-or two-equation models (ANSYS, 2019a). The one-and two-equation RANS models consider Frontiers in Energy Research | www.frontiersin.org December 2021 | Volume 9 | Article 788239 the turbulence to be isotropic and provide a good balance of computation speed and accuracy (ANSYS, 2019a). SRS capture the spectrum of turbulent scales to varying degrees, depending on the model. Direct numerical simulation (DNS) is the most expensive but most accurate method. There is no turbulence modeling in DNS and even the smallest eddies are resolved. DNS is typically reserved for simplified cases and lower Reynolds number flows. Large edge simulation (LES) captures the larger scales but models the smaller scales, operating under the assumption that the smaller scales are more isotropic and more easily modeled (ANSYS, 2019a). Both DNS and LES are computationally expensive, requiring a large number of elements and small time steps. A relatively new class of hybrid RANS-LES turbulence models has been developed, which combines the best features of RANS and LES (ANSYS, 2019a). In regions where the mesh size is fine enough, the LES model is used. In other regions where the mesh size is too coarse for LES, then the RANS model is used.
The SBES hybrid turbulence model was used to simulate flow through LG MJ1 cap to gain additional insight into the flow field beyond what is capable from a purely RANS approach. The SBES model is similar to the detached eddy simulation (DES) model, but uses an improved shielding function to transition between RANS and LES modeled regions (ANSYS, 2019a). However, when using the SBES models, the mesh size has a significant effect on the shielding function which determines the transition between RANS and scale-resolved regions (ANSYS, 2019a). The series of conical refinement regions, described previously in section 2.2, ensured transition from RANS to LES in the region above the vent cap. The final mesh for SRS had 6.33E6 polyhedral elements.
SRS were conducted at a pressure ratio of P 1 /P 2 23.7, which represents the typical pressure difference across the cap at the moment of venting . This differs from the RANS simulations which were conducted at lower pressure ratios, 1.2 ≤ P 1 /P 2 ≤ 8, to maintain consistency with the pressure ratios used in the independent experiments for LIB cap discharge coefficients .
To ensure that the SRS solution was statistically-stationary, the simulation was carried out for 4 ms. The solution was monitored to ensure that mean field variables were independent of time. Between 4 and 7 ms, flow field statistics were recorded to analyze mean and RMS quantities. Figure 10 shows the SRS results. Figure 10A shows an x-normal view and Figure 10B shows a z-normal view of instantaneous velocity magnitude contours taken at a flow time of 7 ms. Note that the Mach disk is visible at the right jet in Figure 10A, and to a lesser extent in Figure 10B. The instantaneous velocity contours show the highly turbulent flow field, as evidenced by the break-up of the jet and mixing of the high momentum jet with the quiescent surrounding air.
The instantaneous velocity magnitude at two points in the jet flow emerging from the holes of the cap were monitored to illustrate the magnitude and frequency of the turbulence in the jet. Point 1 was located at the inner edge of the right jet in the x-normal plane. Point 2 was located at the inner edge of the right jet in the z-normal plane. The instantaneous velocity magnitude at these two points was plotted over the course of 3 ms of flow time in Figures 10C,D, respectively. The mean velocity magnitudes at point 1 and 2 were 230.66 and 196.41 m/s, respectively. The root mean square (RMS) velocity magnitudes at point 1 and 2 were 102.75 and 86.96 m/s, respectively. Turbulence intensity, which is equal to RMS velocity divided FIGURE 10 | SRS simulated instantaneous velocity magnitude for the LG MJ1 cap at P 1 /P 2 23.7 in the x-normal (A) and z-normal (B) planes and instantaneous velocity magnitude plotted for 3 ms of flow time for point 1 (C) and point 2 (D).
Frontiers in Energy Research | www.frontiersin.org December 2021 | Volume 9 | Article 788239 by mean velocity, was 50.4 and 44.3% at points 1 and 2, respectively. The observed turbulence intensity of ∼50% indicates highly turbulent flow (ANSYS, 2019b). The contours of the mean velocity magnitude, represented by Mach number, are shown for the x-normal and z-normal planes in Figures 11A,B, respectively. In the x-normal cross-sectional plane, the effects of the presence of the deformed burst disk was observed as previously in the RANS simulations. The burst disk remained attached and resulted in asymmetric jet structure for the left and right jets in Figure 11A. The maximum velocity in Figure 11A reached Ma 4.03 immediately downstream of the right vent hole. The global maximum velocity was Ma 4.18. The two vent holes shown in the z-normal plane were more symmetric, as shown in Figure 11B.
Contours of normalized RMS velocity (u′/ U max ) in the x-normal and z-normal planes are shown in Figures 11C,D. The highest level of fluctuating velocity observed was in the right jet of the x-normal plane. The value of normalized RMS velocity of the right jet was u′/ U max 0.41. Due to the blockage from the burst disk, the maximum fluctuating velocity of the left jet was u′/ U max 0.29. The global maximum was u′/ U max 0.48. In general, higher levels of fluctuating velocity were observed at the edges of the jet, corresponding to the position of the shear layer between the jet and the quiescent fluid. The shear created between the high momentum and low momentum fluid creates unsteadiness and enhances the mixing between the two fluids.
Isosurfaces of a normalized Q-criterion from the SRS simulation was used to visualize the turbulence developing in the jet region. The normalized Q-criterion is derived based on the second invariant of the velocity gradient tensor and is useful for identifying coherent turbulent structures (ANSYS, 2019a). Figure 11E shows an instantaneous 3-D isosurface taken at a constant normalized Q-criterion value of 0.3. The coherent FIGURE 11 | SRS simulated mean velocity, represented by Mach number, at P 1 /P 2 23.7 for the LG MJ1 cap in the x-normal (A) and z-normal (B) planes, normalized fluctuating velocity magnitude in the x-normal (C) and z-normal planes (D), and a 3-D isosurface taken at normalized Q-criterion 0.3 colored by vorticity magnitude which illustrates coherent turbulence structures in the jet flow (E). The SRS simulation results, conducted at a pressure ratio of P 1 / P 2 23.7, showed the structure of the high velocity and highly turbulent jets emerging from the vent holes. The large pressure ratio accelerated the flow to a global maximum velocity magnitude of Ma 4.18. Unsteadiness introduced by the shear between the high momentum jet and low momentum surrounding fluid resulted in a high value of normalized fluctuating velocity, up to u′/ U max 0.48. The 3-D isosurfaces taken at normalized Q-criterion 0.3 further illustrate the highly turbulent flow by the large number of coherent vortex structures in the jet. The coherent structures grow larger in scale while vorticity decreases as the vortices travel downstream in the jet.

CONCLUSION
In the present work, a RANS CFD model was developed to simulate the steady state venting flow emerging from four LIB caps. The sharp-edge orifice equivalent area for each battery cap obtained from CFD was compared with experiments. The simulated equivalent area agreed well for the MTI, LG MJ1, and K2 caps but underpredicted equivalent area for the LG M36 which was attributed to possible variation in the burst disk position and shape. The steady state velocity field of the flow through the LIB cap obtained from the CFD model showed that at higher pressure ratios, choked flow occurred at varying points for each of the four caps. For the MTI and K2 caps, the choke point occurred near the current collector of the cap which was attributed to low cross-sectional flow area through that component. For the LG MJ1 and LG M36 caps, the choke point was closer to the holes in the positive terminal plate. The choke point location plays an important role in the design of the safety vent since an overly restrictive current collector plate will minimize the benefit from having a positive terminal plate with large flow area. Turbulent viscosity ratio distributions for all four vent caps showed increased turbulent mixing in two regions: between the current collector and positive terminal and along the jets. For the MTI and K2 caps, the steep jet angle caused flow recirculation directly above the vent cap which led to a third region of increased turbulent mixing. The global maximum turbulent viscosity ratio (μ T /μ) of the MTI, LG MJ1, K2, and LG M36 caps at pressure ratio of P 1 /P 2 7 were 4,575, 3,360, 3,855, and 2,993, respectively. Increased turbulent mixing may play an important role in combustion of vented gases (Turns, 2011).
To investigate the turbulent flow field in more detail, scaleresolved simulations (SRS) using the SBES turbulence model were implemented using the LG MJ1 cap. The mean and fluctuating velocity distributions showed that the shear between the high velocity jets and the surrounding quiescent fluid resulted in high levels of fluctuating velocity and, thus, high levels of turbulence intensity. The high levels of fluctuating velocity in the jets agreed with high levels of turbulent viscosity in the jets predicted by RANS simulations. The deformation geometry of the burst disk partially blocked one of the four holes of the positive terminal. The velocity magnitude and the RMS velocity fluctuations for the jet from this partially blocked hole were both dampened. The size and rotation of coherent turbulent structures were also visualized from the SBES simulation results. The coherent structures grew in size as they lost momentum and travelled downstream and the vorticity of those structures also decreased as they traveled downstream. The turbulence intensity and length scale both play an important role in battery safety as both parameters will influence the magnitude of impingement heat transfer from the jets to neighboring surfaces or cells (Kataoka et al., 1987). Thus, the complex turbulence structure of jets emerging from an 18650 Li-ion battery cap, as demonstrated by these SBES simulations, can have a significant impact on model simulations of heat transfer from a Li-ion cell undergoing thermal runaway.
The present work provides a foundation for using highly detailed, 3-D CFD modeling to simulate the venting flow in commercial LIB cells. In the future work, the CFD model will be used to investigate the influence of the flow field on heat transfer to neighboring surfaces. Ultimately, this modeling approach can help improve the design of battery modules to protect against propagating failures driven by venting, combustion, and fire.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication. WL contributed to CFD model development, conducted simulations, performed data reduction, analyzed the results, and prepared the manuscript. VQ contributed to the CFD model development and conducted simulations. KC prepared the LIB cap samples, edited the manuscript, and was responsible for project management. JO analyzed the results, edited the manuscript, and was responsible for funding acquisition and project management. copy-right notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the Naval Surface Warfare Center, Crane Division or the U.S. Government.