A method to quantify the supersonic discharge of airbag cold gas inflators

We discuss a method to quantify the supersonic discharge of airbag cold gas inflators. Since one primary quantity of interest, the mass flow with time, is not directly measurable, a combined experimental and numerical approach was chosen. Shadowgraph and schlieren images visualize the gas dynamic process. Pressure measurements were conducted inside the inflator and downstream of the outlets in the supersonic jet. In this context, a method to measure the pressure of the flow without effects of shock reflection is presented. The temperature inside the inflator was estimated using a fast response heat flux probe and assuming different scenarios for the unknown heat transfer coefficient. Then, a numerical model of the inflator was created. The experimental results served as boundary conditions and some basic sensitivities remaining unknown from the measurements were studied to verify the numerical outcome. The numerical model was verified using experimental results. Finally, the mass flow rate was derived from the numerical model and compared to an analytical model. The method can reconstruct the temporal evolution of the mass flow discharging from the inflator, the pressure and the topology of the flow field within acceptable bounds. Furthermore, the method can deliver inflow data for subsequent airbag inflation studies.


Introduction
While the number of traffic accidents in Germany remained approximately constant at 2.3 million a year during the last 30 years, the number of deaths continuously decreased from 11,300 in 1991 to 2719 in 2020 (Statistisches Bundesamt 2021). The improving safety systems made an important contribution to this development. An essential part of the safety system of automobiles are airbags. There are various airbags in an automobile, e.g., the front airbag, side airbag, knee airbag or curtain airbags. Airbags can also be part of the safety concept of other vehicles like planes, motorcycles or bicycles.
In state-of-the art airbags, the fabric bag itself is inflated by two basic principles: hot or cold gas inflators. Also hybrid gas inflators are commonly used. Hot gas inflators generate the gas by the combustion of a solid propellant. They are typically installed in the steering wheel for example. The alternative are cold gas inflators, where the filling gas is contained in a high-pressure reservoir sealed with a membrane. The inflation process is initiated by an igniter, which weakens the membrane until it ruptures, such that the expanding gas flow inflates the airbag. Cold gas inflators are mostly used in curtain airbags, which are typically located above the b-or c-pillar and protect the heads of the occupants in the event of a side impact.
While airbags are a major driver for increasing safety, in some cases they may be the cause for severe injuries, particularly in so-called out-of-position scenarios, when in the first few milliseconds the passenger hits the airbag before it is fully deployed (Kramer 2013). For front airbags, occupants in out-of-position scenarios may experience severe injuries to the cervical spine due to shear and bending forces to the neck (Hahn et al. 2021). Curtain airbags in out-ofposition scenarios have shown to reduce the acceleration on the head compared to no airbag scenarios, while increasing the forces on the neck (Potula et al. 2012).
During design and development, numerical modeling is used to simulate and improve the inflation process of airbags in detail. Particularly the early stage of airbag deployment is important to study out-of-position cases. A method to simulate airbag inflation is the "uniform pressure" method, which assumes a uniform pressure distribution inside the airbag. This method is not able to predict the correct inflation process, because the jet exiting the gas inflator is not modeled (Mestreau and Loehner 1996). Fluid motion is not part of this method. Other approaches are coupled fluids-structure methods, which are more complex. Meduri et al. (2022) and Rieger (2006) demonstrated for front airbags, that coupled fluid-structure methods capture the inflation process more accurately. However, curtain airbags have a more complex structure than front airbags, consisting of multiple channels and subdivisions of the bag, making it even harder to simulate the inflation process accurately. Potula et al. (2012) used the uniform pressure and a SPH (smoothed particle hydrodynamics) model for the simulation of a curtain airbag. For all scenarios, the airbag was fully inflated when hit by the occupant. The correct deployment of the curtain airbag was not discussed. Fokin and Dessarud (2007) have used a stepwise constant pressure approach to simulate the inflation process of a curtain airbag. The comparison between simulation and experiment showed that the inflation delay times between the compartments/subdivisions of the bag were consistent, while the impact accelerations partially deviated. In addition, some subdivisions did not inflate properly in the simulations. They concluded that one problem of the simulation was a lack of information about the gas inflator exit jet, so that they had to estimate a velocity as input parameter.
Input data for numerical simulations are usually generated by analytical models or tank tests. For the latter, the outflow of the gas inflator is collected in a tank and the pressure rise is measured. The total mass flow rate can then be estimated based on the pressure rise, assuming, e.g., a homogeneous distribution. The results may vary depending on the tank, Rieger (2006), because the flow and the distribution of the physical quantities like pressure and temperature may not be homogeneous. Analytical models are also often too simple and inaccurate to map the complex process adequately. Voges et al. (2007) did a first step to analyze the flow field of cold gas inflators by using PIV (particle image velocimetry) to quantify the supersonic jets of a replica of a gas inflator into a pressurized environment. Later, Klinner et al. (2008) characterized the jets of cold gas inflators with the background oriented schlieren technique (BOS). To our knowledge, no other published works related to the flow field of cold gas inflators exist, while accurate numerical modeling of the inflation process requires proper boundary conditions and capturing of the gas dynamical processes.
Herein we describe a method to quantify the supersonic discharge of cold gas inflators. The focus is to quantify the main parameters of the discharge process, which are the pressure distribution and the mass flow rate. The method employs a dual approach: First, experiments are conducted. Shadowgraph and schlieren images visualize the flow field. The pressure is measured inside the gas inflator and downstream of the outlets in the supersonic jet. Since the flow field features moving shock fronts, a technique to measure the pressure without spurious effects of shock reflections is presented. The temperature of the discharged gas is estimated by using heat flux probes and assuming different scenarios for the unknown heat transfer coefficient. However, the mass flow rate cannot be measured precisely and all measurements are point measurements. Therefore and to further break down and understand the process, numerical simulations support the analysis. Since the geometry of the gas inflator is complex, the numerical model avoids to discretize the entire gas inflator by using the measured data as boundary conditions. The numerical model is verified by comparison with the experiments. Finally, the mass flow rate is derived from the numerical model and compared with a simplified analytical model.

General overview and pressure measurements
All studies herein are based on the cold gas inflators from iSi Automotive. This type contains 13.2 g helium with a pressure of 660 bar and is designed for curtain airbags. Since herein the discharge process itself and the flow field downstream of the gas inflator outlets is the primary interest, the flow into a free environment was studied; the airbag itself was not attached to the gas inflator. A detailed sectional view of the gas inflator is given in Fig. 1b. Figure 1a illustrates the overall experimental setup. The cold gas inflator is mounted in a clamping and covered by a plate, so that ultimately only the so-called "outlet tube" protrudes through the plate. The cover plate of the clamping also serves to establish a well-defined and simple boundary condition (wall) for the numerical simulations.
All gas inflators are off-the-shelf, standard devices delivered by the manufacturer. Most measurements were done sequentially and, moreover, we have acquired mostly five, sometimes ten individual samples of each test point for ensemble averaging. Therefore, it was an important aspect during the design of the test stand to allow efficient and quick exchange of the gas inflators and to use the devices without (or with only minor) adaptions. In sum, more than 200 gas inflators were fired for the studies explained herein.
In the outlet tube, the static and total pressure and heat flux were measured at several positions, see Figure 3. To minimize the need for adaptions of the gas inflators, the pressure and heat flux sensors were mounted in "sensor brackets," which were clamped over the outlet tube from outside. These are shown in Fig. 2b and also in some of the subfigures of Fig. 4. The advantages and disadvantages of this solution will be discussed in the text below.
To measure the static pressure in the outlet tube, p s , four 1.2 mm holes were drilled into the outlet tube for each gas inflator sample. The outside of the taps is then "covered" by the sensor brackets, which have the sensors mounted in a respective counterpart. This setup requires a proper sealing between the outlet tube and the sensor brackets, which is complex due to the curvature of the outlet tube, the small radius and the high pressure. The use of FKM (fluoro-elastomer) 90 shore seals at the passage from tube to sensor brackets (see also Fig. 4 II) and PTFE tape suitable for high pressures on the thread of the pressure sensors proved to seal the pressure taps without leakage. Tightness was evaluated by several tests, which, for the sake of conciseness, will not be presented herein. One disadvantage of this design is that a finite volume exists in front of the sensors, which may lead of spurious amplification and damping, that will be discussed in Sect. 2.3. Moreover, moving shock fronts can enter the volume and falsify the pressure measurements, which will be discussed in Sect. 2.4.
To measure total pressure in the outlet tube, p 0 , a 1.9 mm outer diameter pitot-probe with inner diameter of 1.2 mm and 17.4 mm length was inserted from the top face of the outlet tube into the center axis of the outlet tube. A very similar setup is used to insert the cylindrical heat flux probe (see also Sect. 2.5), which also has a diameter of 1.9 mm and was inserted from the top face of the outlet tube. Downstream of the outlets, the flow field becomes supersonic. To measure the total pressure profile in this supersonic jet, cone-shaped pitot probes with a blunt front face were used, where the pressure sensors are integrated in the probe tip. These probes were attached to a traverse, see Fig. 2a. The traverse arms are wedge-shaped to avoid backlash on the flow. The pitot probes were traversed symmetrically on both sides of the outlets in a distance 2 ⋅ H and 3 ⋅ H , where H is the height of the rectangular outlets. The supersonic jet transversal profile was resolved with four positions with a spacing of 5 7 H , see Fig. 3. Since the flow field is supersonic, the pitot probes measure the total pressure behind their own bow shock, p 02 , which is assumed to be a normal shock. The total pressure in the flow itself, i.e., the undisturbed pressure upstream of the sensor shock, will be denoted p 01 .
For all pressure measurements, fast response pressure transducers of type Kulite HEM-312M with a so-called "B-Screen" were used. The sensor has an outer diameter of 6.35 mm and it is worth to note that smaller sensors are actually available. The reasoning behind the screen and why a rather large sensor is beneficial will be discussed in Sect. 2.4.
An amplifier from Imtron with an integrated 4th order Butterworth filter and a cutoff frequency of 30 kHz amplified the signal. A M2i.4652 Spectrum card sampled the analog signal with a sampling frequency of 50 kHz. Additionally, the results were low pass filtered after A/D conversion according to the f 5% falsification boundary that will be discussed in Sect. 2.3.
Since it is likely that high-frequency pressure fluctuations exist, the setup may feature alias noise in the signals. Therefore, only transient ensemble mean values are accessible. The ensemble means were derived from mostly five consecutive samples (i.e., five different gas inflator ignitions). The measurements of the static pressure p S in the outlet tube are based on ten samples. Since the individual runs may vary, a 95% confidence range will be shown in all figures, which have been derived from the Student's t-distribution.

Gas dynamical process
Many aspects of the measurements that will be discussed in the following subsections are due to the complex instationary processes during the gas inflator start-up and discharge. Therefore, it is worthwhile to illustrate the process schematically, Fig. 4a-g, to explain the gas dynamical phenomena and associated challenges for the measurements.
(a) The discharge process starts with the activation of the igniter, which weakens the membrane until it ruptures. Since the installed flow nozzle is the smallest cross section within the system, it effectively controls the total mass flow rate. Due to the high pressure, a moving shock front is immediately produced, which propagates (faster than the local speed of sound) through the air into the volume of the cold gas inflator. Because of the complex interior geometry, multiple shock reflections are to be expected. The helium, which begins to dis-

Fig. 4
Gas dynamical outflow process of the cold gas inflator (not drawn to scale); thick lines indicate shock structures charge just behind the shock front, and the surrounding air are divided by a contact discontinuity, similar to the Sod shock tube problem, ref. Sod (1978). The shock front, the displaced air and finally the helium all pass through a filter, which basically is a cylindrical component with several holes to filter out any bigger particles, particularly the remnants of the membrane, so that they cannot harm the airbag. Note that the expansion wave is created at the same time as the moving shock front that propagates into the gas reservoir. Since it does not represent a challenge for the measurements, it will not be discussed in detail. (b) Once the shock front and the helium-air interface have passed the filter, they propagate into the outlet tube.
Due to the complex geometry, the shock front is not necessarily perpendicular to the outlet tube center axis and also multiple shocks may be present due to reflections. But in any case, at the very inflow of the pressure tap the shock passing in the outlet tube creates a new moving shock front inside the tap as a result of the sudden pressure rise. This secondary shock front propagates into the cavity of the sensor and reflects either directly on the sensor chip or, if present, on the sensor screen. (c) The primary shock front passes through the outlet tube, finally leaves the outlet tube and propagates into the free environment. The same is true for the interface of helium and air. Since the cold gas inflator has a comparably small volume downstream of the membrane, we can assume that the helium-air interface is still closely following the initial shock and that, once the shock has reached the outlets, all air has been pushed out of all inner volumes at this point. (d) The moving shock front eventually reaches the pitot probes. In contrast to the well-known classic "bow shock," which is a quasi-stationary phenomenon, the moving shock front propagates through air at rest and therefore interacts with the pitot probe face at the impingement area A I . (e) The shock front is partially reflected at the impingement area A I of the probe. Additionally, a secondary moving shock front inside the pressure tap of the pitot probe propagates through the tap to the sensor. The shock reflection at the impingement area A I of the probe face results in a significant pressure amplification, which will be discussed in Sect. 2.4. The rest of the primary shock front propagates further into the environment. Note that this process also happens to the pitot and heat flux probes, when inserted into the outlet tube, because the shock reflection and the associated amplification of the pressure are not due to a supersonic velocity of the flowing gas stream, but it is a result of the propagation of the initial shock wave through a medium at rest.
(f) Since the initial shock front travels faster than the speed of sound it quickly leaves the flow field around the pitot probe. The helium-air interface then reaches the sensor, and thereafter, a supersonic flow field is established with the classic bow shock in front of the pitot probes. (g) The whole start-up process of the gas inflator takes about 0.5 ms. Finally, a flow field with two underexpanded jet cells on either sides of the outlet tube is created. The type of flow topology is also known as a "barrel shock" with the characteristic Mach disk. More information can be found, e.g., in Franquet et al. (2015). In contrast to the start-up process, the state with the two supersonic jets is a more-or-less quasisteady process, where the pressure in the gas inflator and consequently the static and total pressure in the outlet tube is decreasing due to the discharged mass. The total discharge process takes about 100 ms. However, approximately 42% of the gas mass is discharged within the first 10 ms.
Time will herein always be relative to the trigger of ignition or, more specifically, the rising edge of the ignition voltage, which thus represents t = 0 ms.

Frequency response characteristics
The signal response characteristics were studied utilizing the system response model (SRM) by Bergh and Tijdeman (1965). It is applicable also to short tube lengths, as found by Semaan and Scholz (2012). The SRM assumes a transducer volume in front of the chip, that is connected via a tube to the desired measuring point, like schematically shown in Fig. 4d. Generally speaking, for a high-frequency response, a small transducer volume and a short pressure tap length l t with large diameter are desirable. Note that the SRM model is not fully applicable, because the presence of a screen in front of the sensor is not included in the model. Despite this arguable simplification, it was used to estimate the frequency response of the system. Figure 5 illustrates the installation of the pressure sensors in the pitot probes and the position of the screen. The pressure sensors in the sensor brackets to measure the pressures p S and p 0 in the outlet tube are essentially similar; they only differ in terms of the tube length l t . In all cases, the transducer volume has been seen as the combined volume of the cavity in front and behind the screen. While the volume behind the screen is specific to the transducer, the volume in front was reduced as much as technically possible by fabricating a precision thread with stop for the sensors/ transducers in, both, the pitot probes and the sensor brackets.
For the sensors mounted in the sensor brackets to measure p S and p 0 in the outlet tube, the tube length l t must be long enough to connect the sensor with the positions shown in Fig. 3, which is 3.35 mm for the static pressure taps and 17.40 mm for the pitot probe. These tube lengths are as short as technically possible. We chose a diameter of 1.2 mm, resulting in a first resonance frequency f 1.R and a 5% falsification frequency f 5% given in Table 1, as based on the SRM.
However, for the pitot probes to measure p 02 in the supersonic jet, an additional trade-off is required: Pitot probes for supersonic flows must be conical. A small half-cone angle is desirable in terms of detached shocks and slim probes. We have chosen a half-cone angle of 10 • , which is kept always the same. Since the size (i.e., the diameter) of the sensors is also given, the only design choice for the pitot probes is the length of the probe head; different versions are shown in Fig. 5. A short probe head with small l t is beneficial in terms of frequency response, but it increases the area of impingement A I . The probes will be used to measure p 02 in a field with large spatial gradients, where a small tip is desirable to increase the precision and avoid physical "averaging" effects over the area A I . Figure 6 shows the frequency response for different tube lengths l t . The amplitude ratio is the ratio of the amplified or damped signal to the original signal. As expected, for short tube lengths higher frequencies can be resolved. The first resonance and 5% falsification frequency f 5% are also summarized in Table 1. The latter was chosen as the maximum allowed amplification. As based on this study, all signals were low pass filtered after the A/D conversion utilizing the f 5% given in Table 1.
Since the SRM may not be fully applicable, the different pitot probes with lengths of 3 mm, 10 mm and 21 mm were used in a preliminary test. The data will not be shown here for the sake of conciseness, but it was found that the signals of the probes with l t = 3 mm and 10 mm had similar characteristics, while for the 21 mm probe more pronounced resonance and damping effects could be found, so that the results of the SRM could be essentially verified. Furthermore, the 21 mm tube length has led to a considerably longer response time. Therefore, the pitot probe with a length of 10 mm was chosen for the total pressure measurements in the supersonic jet.

Shock reflection free pressure measurement
During the initial start-up of the outflow process, the moving shock front dominates the flow physics. When a shock front reflects on a surface, the surface experiences the so-called reflected pressure p r . As will be shown in this subsection, the pressure reflections must be taken into account, or better be avoided, for accurate measurements. The Rankine-Hugoniot relations deliver an expression for the peak reflected pressure p r , ref. Needham (2018):  where p 1 , 1 and u 1 are the static pressure, density and velocity behind the moving shock front and is the ratio of specific heats. As can be seen, the peak reflected pressure is at least two times higher than the static pressure, like in an elastic collision. In practice, the reflected pressure may be several times higher, because the molecules that first move alongside the shock front eventually reflect from the surface, but are not able to reflect linear-elastically since they are obstructed by incoming gas molecules. This leads to a higher pressure than in a linear-elastic scenario, ref. Karlos and Solomos (2013). The reflected pressure becomes maximal for a frontal reflection, whereas some angle of incidence reduces the effect (Karlos and Solomos 2013). Figure 7 qualitatively illustrates the expected process: A shock front propagates through a stagnant ambient medium and arrives at a certain point P at the arrival time t a . An idealized, non-intrusive, massless sensor at the point P would measure the true total pressure of the fluid behind the shock front p 0,SF . As the shock front moves away from P, the pressure subsides to the ambient pressure at the time t ∞ .
In contrast, with a physical sensor at point P, the shock front reflects on the sensor surface and the reflected pressure p r is measured instead. The peak reflected pressure p r is several times higher than the "true" pressure of the fluid p 0,SF . For an infinite surface (i.e., total shock reflection), the decay of the reflected pressure to the ambient pressure p ∞ is rather slow and takes until time t ∞,inf . For a finite surface, the reflected pressure decays quicker. The sensor at P actually measures the "true" pressure of the fluid after some clearing time t c . The concept of the clearing time was actually developed to define loads on structures after explosions, UFC (2008), but can be qualitatively transferred to the case herein. Accordingly, the clearing time is proportional to: (1) where c r is the speed of sound in the reflected area, D is the diameter of the probe head and A I is the impingement area, see also Karlos and Solomos (2013). With the data from Table 1 for the l t = 10 mm pitot probe, an order of magnitude estimation with Eq. (2) yields a clearing time of 34.5 μ s, if we assume a proportionality constant of the order unity (which is justified by UFC (2008)). However, as will be shown later, the clearing times practically experienced in our setup are two orders of magnitude larger.
The reason for the excessively long clearing times experienced herein is that the sensor chip, i.e., the sensible membrane of the transducer, is surrounded by a constrained volume: As explained in Sect. 2.1, the shock front passing the pressure taps (of the static pressure taps or any of the pitot probes) creates a new, secondary shock front that propagates inside the cylindrical cavity to the sensor. Due to the inner geometry, several reflections can be expected. Once the transducer volume is pressurized to the peak pressure, the molecules are "trapped" inside the cavity and can only be released via the pressure tap itself, resulting in long clearing times.
For practical measurements, the first step is to prevent the secondary shock front to reflect on the sensor chip, which would lead to excessive amplification and can easily damage the sensor. As will be shown, the sensor screens play a major role for this task. At the same time, the clearing time should be as small as possible. Preferably, the clearing time is of similar order of magnitude as the response time of the sensor from Sect. 2.3. We cannot see a way to theoretically or numerically predict the clearing time for the pressure taps and probes. The later approach requires very small cells to discretize the sensor cavity adequately that will lead to very long simulation run times. Therefore, an experimental study was conducted with different sensor diameters and different screen types.
For the static pressure measurements p s in the outlet tube (OT), three different sensor types were used: Two Kulite XTL type pressure transducer with 3.8 mm diameter, one without any screen, one with a so-called "B-Screen," and the HEM type pressure transducer with 6.35 mm diameter with a B-Screen. Figure 8 illustrates the size of the sensors and the screens of the respective sensors relative to the size of the pressure tap. Figure 9a shows different results of static pressure measurements p s in the outlet tube. The arrival time is about 0.5 ms after ignition. For the sensor without a screen, the specific curvature of a shock reflection is recognizable with massive amplification (peak pressures > 400 bar) and slow clearing. Note that in our case, in contrast to the illustration Reflected pressure over time for infinite and finite surfaces, adapted from Karlos and Solomos (2013) in Fig. 7, the reflected pressure will not be cleared to ambient pressure, because the helium-air interface closely follows the shock front. Instead, the reflected pressure subsides to the static pressure of the helium, which is the quantity we are interested in. The same sensor type with a B-Screen reduces the peak reflected pressure by approximately 50% and the clearing time from 4 ms to 2.4 ms. Due to the B-Screen, the transducer volume is split into two cavities, as shown in Fig. 5b. The main reflection now happens on the front face of the screen and not on the sensor chip itself. It seems that a more winding cavity is beneficial.
To follow this idea, the HEM sensor was used, which is much larger than the XTL-type. Therefore, the holes of the screen are further away from the pressure tap compared to the screen of the XTL type, see Fig. 8, which creates an even more "winding" cavity, and it is less likely that shock fronts will enter the cavity between the chip and screen and reflect back and forth. Also the transducer volume is larger, which may be beneficial to diffuse the reflections as much as possible. The HEM type shows no signs of a shock reflection.
A similar study, though with a different sensor for comparison, was done for the pitot probes in the supersonic jet for p 02 . Here, as the only comparison, a XTL type pressure transducers with a so-called "M-Screen" and the HEM type with a B-Screen were used. The M-Screen is basically a very fine mesh, see Hurst et al. (2014), which does not separate the cavity too much. Figure 9b shows the results of the measurements. The pressure reflection phenomena for the measurements with the M-Screen are pronounced significantly. Compared to the measurements with the HEM type, the pressure is approximately 10 times higher and the clearing time about 17 ms. The HEM type still shows no signs of a shock reflection.
The outcome of this comparison is that a larger sensor with a B-screen seems beneficial, although we cannot quantify the relative importance of individual parameters to predict the "optimal" sensor for such a test. The design of a shock reflection free setup is (at least partly) contradictory to the requirements for good frequency response that have been discussed in Sect. 2.3, since a large sensor size, larger transducer volume and a screen are typically   Clearing time with different sensor size and screens for measurements of (a) static pressure p S in the outlet tube and (b) total pressure p 02 in the supersonic jet detrimental for the frequency response, but seem necessary to avoid shock reflections. For an optimal design, the transducer volume would be increased in favor of the clearing time, but at the expense of the response time, until both are approximately equal. In any case, there will be a "delay time," which we define as the time, over which a physical sensor does not read the true pressure of the fluid, due to limited response time (minimum delay time) or due to shock reflection phenomena, or combinations thereof. Table 2 lists the delay times for the various sensors tested. We cannot exactly determine the true delay time for the HEM type sensor. From the curves shown in Fig. 9, we estimate that the it is around 0.5 ms to 1 ms. To put them into context for the overall objective, the table also gives the mass of helium that has left the system during the delay time. This value has been estimated from Eq. (6) that will be introduced in Sect. 4.
For the measurement of p s in the outlet tube, the HEM type sensor with a B-Screen reduces the error from 19.92% to the level of the response time between 2.73 and 5.38%. For the measurement of p 02 in the supersonic jet, the B-Screen and HEM type sensor reduce the error from 61.06 to 2.73-5.38%. All further data shown in this paper have finally been measured using the HEM sensor with the B-screen and the l t and A I as given before.

Gas temperature estimation from heat flux measurements
To estimate the mass flow discharged by the cold gas inflator, additionally temperature measurements have been conducted inside the outlet tube. Classic thermocouples have a response time in the range of milliseconds and thus are too slow to directly measure the gas temperature. A physically close quantity is the heat flux, where sensors are available with response times in the range of microseconds. Equation (3) gives the relation between gas temperature T G and heat flux ̇q . Gas temperature can be calculated, if the heat flux and surface temperature T S are measured and if the heat transfer coefficient is known.
Heat flux sensors effectively measure the surface temperature of the probe tip. Using the theory of one-dimensional heat conduction into a semi-infinite body, heat flux is then derived by (Schultz and Jones 1973): where T S is the surface temperature at the probe tip and t is time. ck are known thermal properties of the sensor, which are the density , specific heat capacity c and thermal conductivity k.
Finally, from Eq. (3) only the heat transfer coefficient remains unknown. To our knowledge, no accurate models to estimate the heat transfer coefficient are known for these conditions. Since the total gas temperature serves as a boundary condition for the numerical simulations, different scenarios were assumed for the heat transfer coefficient and the sensitivities were studied in the numerical simulations, which will be discussed in Sect. 4.
For the measurements, the MCT 19 coaxial thermocouple from Dr. Müller Instruments of type E with a response time of 3 μs was used, see Fig. 10. The advantage of coaxial thermocouples is the low response time, because one thermocouple material is swagged over the other one, resulting in a very thin junction, ref. Olivier (n.d). The conversion of the voltage into the temperature was done using the IEC (1995) tables. The same amplifier and A/D converter were used as for the pressure measurements.

Shadowgraph setup
Schlieren and shadowgraph techniques were used to visualize the supersonic flow field downstream of the outlets.
The flow field is highly turbulent and shock structures, like the barrel shock, are sometimes hard to contrast against the turbulent eddies. Also, the flow field changes continuously, so that averaging of the images cannot be applied. Therefore, both, schlieren and a "focused" shadowgraph technique were tested. Shadowgraphy is sensible to the second spatial derivative of the refraction index and thus density, while schlieren is sensitive to the first spatial derivative (Settles 2001). Since shocks usually feature larger spatial derivatives than turbulence, this method eventually produced better visible shock structures. Furthermore, the sensitivity can be adjusted in a "focused" setup (Settles 2001). Figure 11 illustrates the setup of the "focused" shadowgraph setup. The schlieren setup was mostly the same, but with a knife edge added. A Phantom v611 high-speed camera recorded the process with a sampling rate of 5 kHz and 1280 × 800 pixel resolution. A 100 W halogen lamp illuminated the flow field.

Numerical simulation methodology
Since the experiments only give punctual data points, numerical simulations were used to reconstruct the full flow field. Particularly, the mass flow cannot be measured directly, but must be extracted from accompanying simulations. On the other hand, the inner geometry of the cold gas inflator is very complex and simulation/discretization of the entire geometry comes with many challenges and unknowns. Therefore, for the numerical simulations described herein, the numerical domain was reduced to the region, where experimental data are available and the experimental data were actually used as inlet conditions. Figure 12 shows the computational domain and the mesh. We assume the flow is relatively homogeneous inside the outlet tube; therefore, a quarter-symmetric setup with two symmetry planes was used. The mesh consists mostly of hexahedrons with low aspect ratio. The medium mesh with a nominal cell size of 250 μm , which is equal to 28 cells per H and approximately 40,000 cells inside the outlet tube, has been found reasonable in a grid sensitivity study that will be discussed in Sect. 3.2. In total, the mesh consists of 4.97 million cells.
The density-based solver rhoCentralFoam from the open source CFD package OpenFOAM was used. Detailed information can be found in Greenshields et al. (2010). The central upwind scheme of Kurganov et al. (2001) was used for spatial discretization. Time discretization was done using an implicit second-order backward scheme. The Gauss scheme with second-order linear interpolation was chosen to approximate gradients, divergences and Laplacian terms. For interpolations, the limiter of van Leer (1974) was chosen, which switches from linear to first-order upwind interpolation in areas of strong gradients, e.g., near shocks.
For the density-based solver the pressure is based on the equation of state from Peng and Robinson (1976). The flow field is assumed to be pure helium. Since the flow experiences a broad range of temperatures, seventh-order polynomial equations were used to model the molecular viscosity and the thermal conductivity depending on temperature. The polynomials are based on data from multiple authors in a range from 5.5 to 4000 K. For the molecular viscosity, data from Johnson (1960), Maitland and Smith (1972) and Amdur and Mason (1958) were used. For the thermal conductivity, the data from Vargaftik and Yakush (1977) and Collins et al. (1965) were used.
Turbulence is modeled with the k − -SST DES model of Menter et al. (2003). It is effectively a hybrid RANS/LES approach, since the model reduces the modeled turbulent viscosity in areas, where turbulence can be resolved depending on the grid size and the turbulent length scale. It turned out that the mesh needed to be fine enough to resolve shock structures is also fine enough to resolve the larger scales of turbulence. Despite the arguable simplification with two symmetry planes and thus assuming the resolved turbulence is axis symmetric, the DES model increases the quality of the simulation while adding practically no computational costs.
The measured time varying total pressure and total temperature curves shown in Fig. 13 were used as inlet conditions. In addition, the direction of the flow velocity was assumed to be parallel to the tube axis. For the sake of stability, the oscillations of the measurements were removed and a moving average was actually used. The simulation run time is 4 ms, covering the most dynamic part of the discharge process of the gas inflator. In Fig. 13, two additional inputs, one for the total pressure and one for total temperature, are shown, respectively. These are variations of the boundary conditions that will be discussed at the end of this section.
It should be noted that the inlet plane for the numerical simulations shown in Fig. 12 is not equivalent to the section where total pressure and total temperature have actually been measured, see Fig. 3. Nevertheless, the (low pass filtered) measured values were used as-is, with no further adaptions. As will be shown in Sect. 4 the variations of pressure or temperature between these two stations are effectively neglectable. As will also be shown, it is, however, very important to discretize the "step" visible in both, Figs. 3 and 12, since this step creates a shear layer and ultimately leads to turbulent structures in the upper outlet tube that are important to catch the general flow topology that was observed in the measurements.
Finally, for the outlet (far field) a non-reflective boundary condition (waveTransmissive) was chosen and the walls fulfill the no-slip condition, although no care was taken to actually resolve boundary layers; wall functions are used for all variables.
Note that static pressure in the outlet tube, the velocity flux, density and, thus, the mass flow are then an outcome of the simulation.
To discuss various sensitivities in Sect. 4, different cases and inlet boundary conditions were studied with the numerical simulations: To see the impact of varying the total pressure boundary condition, in one case the input data for total pressure were increased by 4%, which is the estimated uncertainty of the experimental data. The respective input data are added in Fig. 13a.
For the total temperature boundary condition, the specific problem arises that the heat transfer coefficient between the gas and the probe tip is unknown, as has been discussed in Sect. 2.5. Therefore, the true gas temperature is unknown. To cope with this problem, we have assessed three different scenarios: At first, = ∞ W m 2 ⋅K , so that the gas temperature is equivalent to the measured surface temperature of the probe. This is the limiting case, i.e., the maximum temperature that could exist in the flow for negative heat flux values and the minimum temperature for positive heat flux values. The heat flux is first positive and then becomes and stays negative at approximately 1.5 ms. Secondly, we assumed a constant = 2 ⋅ 10 5 W m 2 ⋅K , which is based on an order-of-magnitude estimation for stagnation flows. The respective two curves of temperature over time that were taken as a boundary condition are shown in Fig. 13b. To generate a third limiting case, we also did a simulation with a constant total temperature of 293.15 K at the inlet. The latter represents the (hypothetical) case, that the total temperature rise generated by the moving shock fronts and/or the heat input of the igniter and/or total temperature decrease due to expansion waves subside quickly and have effectively no influence on the flow.
As already mentioned, the "step" in the outlet tube, shown in, both, Figs. 3 and 12, is quite relevant to the flow. To demonstrate this, a case "no step" was run, where the inlet plane for the numerical simulation is the same section where the pressure and temperature data were actually measured, but the step is not existent.
We will refer to our standard case as the "base case," which is the medium mesh, the "Sim. input" curve for total pressure p 0 and the "Sim. Input ( = ∞ W m 2 ⋅K )" for total temperature from Fig. 13. In the base case, the "step" shown in Fig 12 is resolved and the inlet section is slightly upstream of this step.

Grid convergence study
Three different meshes were generated to study the sensitivity of the mesh resolution on the flow topology, pressures and the final mass flow rate. The meshes "coarse," "medium" and "fine" feature a nominal cell size of 500 μm (H/14), 250 μm (H/28) and 125 μm (H/56), respectively. Figure 14 illustrates the resolution of the Mach field at t = 0.9 ms. For the medium and fine mesh, a vortical structure can be seen, which originates from the upper region of the outlet tube and then extends across the supersonic region. This vortical structure is highly dynamic and it varies continuously in shape and position, introducing large-scale, low-frequency fluctuations that will be visible throughout all data to be discussed herein. It also distorts the barrel shock, making it much more dynamic. On the coarse mesh, this vortical structure is not captured adequately and as a consequence the whole flow field and the barrel shock appear as a "classical" barrel shock shape without distortions. The vortical structure can actually be seen in shadowgraph images, as will be discussed in Sect. 4 and was found to be an important phenomenon of the discharge process. Figure 15a shows the sensitivity of the mass flow rate on the mesh resolution. All meshes, even the coarse one, give the same mean mass flow rate. The fluctuations produced by the medium and fine meshes are induced by the vortical structure originating from the outlet tube. Since it is not resolved properly by the coarse mesh, also no fluctuations are resolved.
In terms of total pressure, the biggest difference between the mesh resolutions was found for the total pressure p 02 in the barrel shock at the position z = 0 and x = 2H shown in Fig. 15 (b): The fine mesh actually resolves more fluctuations of the total pressure which are stronger. However, still the medium and fine meshes have the same overall level of the pressure of approximately 16.5 bar between 1 ms and In summary, the coarse mesh does not capture the vortical structure adequately and also gives slightly offset total pressure p 02 in the barrel shock. The medium and fine meshes result in generally similar temporal evolutions, however with different instances of the fluctuations. Considering the high computational demands of the fine mesh solutions, all further results discussed herein are based on the medium mesh.

Sensitivities of the flow field
This section compares the numerical results with the experiments to verify the numerical model. The pressure measurements inside the outlet tube and in the supersonic jet serve as quantitative comparison. Shadowgraph and schlieren images are used to compare the topology of the flow field. Some sensitivities of the flow field will be shown as already indicated in Sect. 3.
Two major characteristics of the discharge process are the moving shock front that propagates through the entire gas inflator and eventually leaves through the outlets and the vortical structure that emerges in the upper region of the outlet tube. Figure 16 gives an overview of the flow topology as seen in shadowgraph images. In Fig. 16a, the moving shock front and its reflection on the lower wall followed by the bright helium jet can be seen. We assume that the brightness is caused by hot particles that emit light. Since the camera is always overexposed in this region, the contrast line seen in the image is not necessarily the "true" interface between the helium and the surrounding air, but the image gives a good impression that the interface follows the shock rather closely. Figure 16b gives an overview over the supersonic jet topology. The Mach disk can be seen and the jet boundary next to the outlets. Also, the vortical structure emerging from the outlet tube is just visible. Often, multiple vortical arms are present. The topology is highly dynamic-shape and size of the barrel shock change continuously, as does the vortical structure. The region of the barrel shock features a high level of smaller scale turbulence, such that the Mach disk and the jet boundaries cannot be exactly traced.
First, we will discuss the different pressures measured in the outlet tube, which is the total pressure p 0 in the center of the outlet tube, which is also used as a boundary condition for the numerical simulation, and the static pressure p S at the tube walls. Both, total and static pressure were measured at the position labeled "measurement plane" in Fig. 3, which is 12.75 mm from the top face of the outlet tube and 12.5 mm downstream of the inlet plane.
As discussed in Sect. 3.1, the total pressure measured in the measurement plane was taken as boundary condition for the "inlet plane" without any adaptions in regards to the "delay time." The total pressure then computed in the "measurement plane" is compared to the experimental data in Fig. 17. As can be seen, the data coincide very well, except that large oscillations exist in the numerical simulation. Recall that the curve labeled "Experiment" is effectively the input data for the numerical simulation and has been filtered additionally to remove the oscillations, i.e., these oscillations are not fed into the simulation, but are a result of the evolving flow field.
In none of the other computational cases, significant deviations were found between the experimental data and the result of the simulations for the total pressure, which is not surprising-it means that the total pressure changes barely over the short distance between the inlet and the measuring position. The same procedure has been conducted for the total temperature boundary conditions. Here, the total temperature prescribed at the inlet plane and the calculated numerically one at the measurement plane were identical and no oscillations occurred. Thus, the total temperature does not change of the short distance either. Hence, the measured data can be used without adaptions as an input.
In contrast to the total pressure, the static pressure p s in the outlet tube is a result of the simulation, and was not prescribed as a boundary condition. Figure 18a shows two different simulation results. The base case reproduces the experimental data fairly well. The numerical simulations consistently show a steeper edge during the ramp-up than the experiments, but this is a direct result of the low pass filtering of the experimental data and the response time of the sensor (see Sect. 2.3). The comparison between the two simulations reveals that the large oscillations, that were already visible for the total pressure, also exist in static pressure, but do not exist for the case with "no step." We will discuss this aspect in more detail below. Nevertheless, both simulations show the same trend for the pressure level, which tends to be slightly lower than what was measured in the experiment.
The pressure measurements are subject to certain repeatability inaccuracies of the cold gas inflators. These are represented in the confidence interval. Since experimental data of the total pressure were used as input parameters with an inaccuracy of approx. 4%, this sensitivity was also investigated. As expected, it was found that an increase of the total pressure by 4% led to an increase of the pressure field by about 4% and thus also the static pressure in the outlet tube increased by about 4% (see Fig. 18b). Thus, there is an appearingly direct linear correlation between the input parameter and the results of the simulation.
To discuss the temporal development of the supersonic jet and the barrel shock, Fig. 19 illustrates the topology of the barrel shock as a sequence of snapshots. The image sequence (a) is the "base case," and (b) the case "no step" and the inlet at the same position as the measurement plane.
Due to the changing total pressure over time as in the outlet tube, the barrel shock increases size, with the strongest expansion at approximately t = 1.5 ms (maximum total pressure, as shown in Fig. 13), and then slowly decreases. The same is true for the case with no step, which is therefore not shown here in full detail.
For the "base case", the barrel shock is strongly distorted, particularly at the outer upper region. In contrast, for the case "no step" the flow topology follows the "classical" barrel-shape. Thus, the distortions of the barrel shock are induced by the vortical structure emerging from the outlet tube as are the pressure fluctuations and this is consequently the main driver for the strong dynamic motion of the process. We note that we also performed a simulation with an inlet plane in the same location as in the "base case," but still without the step, and the results were practically identically to the case "no step." For a more detailed insight, Fig. 20 shows the threedimensional topology of the vortical structure emerging from the outlet tube using the Q-criteria. For the case "base case" with step, a shear layer beginning at the edge of the step and a dominant longitudinal vortex in the top level of the outlet tube are created. From this vortex, three vortical arms featuring vorticity with one dominant arm evolve into the barrel shock and induce oscillations of the barrel shock that lead to the non-classical shape. We refer to the structure consisting of the longitudinal vortex and the vortical arms as "vortical structure." Also, the flow inside the outlet tube is influenced by the vortex, as can be seen from the vector field. It has been found that the vortical structure is highly dynamic, changing shape, strength and location continuously. Also, the quantity of vortical arms is changing, while the dominant one mostly persists. These findings are consistent with the shadowgraph results. For the case "no step" also a longitudinal vortex is created at the top level of the inlet tube. But the vortex is smaller, quasi-stationary and no vortical arms evolve into the barrel shock. Also, since there is no step, there is also no shear layer. Instead, the flow follows smooth streamlines from the inlet into the barrel shock, as can be seen from the vector field in Fig. 20. This is mainly to show that the presence of the step is a prerequisite for the vortical structure to develop in the numerical model, which in turn is the main driver for the unsteadiness of the barrel shock.
Note that we have not yet discussed the influence of the temperature boundary condition from Fig. 13. The reason is that the choice of the temperature profile at the inlet has virtually no significant influence on the flow topology or the moving average pressures. Only the vortical structure is influenced and thus different oscillations are produced while the moving average values remain the same. However, the temperature has a strong influence on the actual mass flow rate, which will be discussed later.
Also note that none of the simulations actually resolved the initial moving shock. Ultimately, this is due to the fact that the starting process of the shock wave originated from the gas reservoir is not resolved in this simplified model. However, the shock propagates with several times the local speed of sound and a numerical simulation would require excessively fine discretization of time and space. Our analyses showed this would massively increase the computational costs. Since the numerical simulation data are still consistent with the experimental data, we conclude that the moving shock is not particularly relevant for the development of the flow field-except that it affects the sensor readings as discussed in Sect. 2.

Total pressure in the supersonic jet
The evolution of the pressure distribution in the underexpanded jet will be discussed, and in this framework, the numerical data are compared to the experiments in this section.
Pitot probes were used to measure the total pressure p 02 downstream of outlets 1 and 2 simultaneously. Figure 21 shows the results for two different positions in the supersonic jet. The probes near the barrel shocks' edge at z = 15∕7 H measure a peak pressure of about 6.5 bar and are exposed to practically ambient pressure at approximately t = 15 ms. This is due to the shrinking size of the barrel shock due to the decreasing pressure in the outlet tube, so that the pitot probes are no longer in the supersonic region after approximately t = 5 ms, until they are completely outside the jet. However, the probes at the centerline ( z = 0 ) measure a peak pressure of about 18.5 bar and a total pressure higher than ambient pressure until approximately t = 85 ms. Furthermore, a slight increase of the pressure can be seen starting at approximately t = 18 ms. From the schlieren images it is observed that this is the time, when the pitot probes are just behind the Mach disk. Hence, the pressure loss due to the Mach disk (at that time and position) is slightly smaller than the prior Fig. 19 Temporal development of the barrel shock for (a) "base case" and (b) case "no step" pressure loss due to the bow shock in front of the pitot probe.
Additionally, the pressure p 02 normalized by the total pressure inside the outlet tube p 0 is presented. It can be seen that p 02 ∕p 0 at the centerline at z = 0 is rather constant for the first 5 to 10 ms, indicating proportionality for the center region of the barrel shock. This is not true for the probe at z = 15∕7 H at the edge of the barrel shock. Furthermore, it can be seen that the moving average total pressure for outlets 1 and 2 coincides very well.
From this we can conclude that the moving average flow through the two outlets is equivalent, which confirms our assumed symmetry boundary condition. Note that this does not indicate physical symmetry of the flow including oscillations and dynamics at every timestep. Furthermore, the total pressure in the barrel shock is more-or-less proportional to the total pressure in the outlet tube. In other words, once the initial flow has ramped up, the process is quasi-stationary for the first few milliseconds.
For some further comparison with the numerical results, the data from outlet 1 will be used in the following. Figure 22 shows the measured total pressure p 02 with the data from the numerical simulations "base case" at position x = 2 H. The conversion from p 01 to p 02 for the numerical data was done assuming a normal shock. A challenge in this context is the unknown and varying shock stand-off distance s. This is important, because total pressure and Mach number feature notable spatial gradients in the supersonic region and therefore the location, where the total pressure   Fig. 21 Comparison of the total pressure p 02 inside the supersonic jet at a distance x = 2H of both outlets, respectively, and normalized p 02 by the total pressure inside the outlet tube p 0 at different z-positions p 01 is extracted and converted into p 02 has some impact on the result. The shock stand-off distance was estimated from schlieren images and Fig. 22 shows some of the images as an example. Within the first 4 ms, the shock stand-off distance varied between approximately 1.2 mm and 1.6 mm. For the position z = 0 and 5/7 H, the shock is normal, but becomes rather oblique for positions at the upper edge. Since the data are not required as input data, but only as a comparison for validation, we have not taken the effort to generate a database with true stand-off distances for every probe position and time. Instead, we converted the data from the numerical simulations with three different "virtual" stand-off distances, s = 0 mm, s = 1 mm and s = 2 mm as examples. Figure 22a shows the pressure p 01 and three different profiles of p 02 assuming the different shock stand-off distances s. The profiles are approximately symmetrical to a line z = 2∕7 H, since the supersonic jet is directed upwards, see Fig. 19. Furthermore, the total pressure in the shear layers is higher than at the edge of the supersonic region. The total pressure p 02 is approximately one third of p 01 because of the losses due to the shock. Also, experimental mean values with confidence intervals are shown.
For a more detailed comparison of the numerical and experimental results, p 02 over time is also shown for the four measurement points at x = 2 H in Fig. 22b. For the positions close to the symmetry plane at z = 0 and z = 5∕7 H, the results for s = 0 mm are at the lower confidence boundary of the experiments and for s = 2 mm at the upper confidence boundary, while s = 1 mm reproduces the mean value of the experiments quite well. This result is consistent, since the shock stand-off distance seen in the schlieren images is mostly closer to 1 mm than to 2 mm, see Fig. 22c Profiles of total pressure p 01 and total pressure behind a normal shock p 02 at t = 2 ms (a) and total pressure p 02 measurements for four positions at x = 2 H with different shock stand-off distances s (b). Schlieren images of the pitot probes in the supersonic region at t ≈ 2 ms (c) the schlieren images were done for the pitot probes of the XTL type sensors, but the cone angle and impingement area of the pitot probes are practically the same. At the edge near the barrel shock, no sensitivity regarding the shock standoff distance can be seen. Since the shock is rather oblique it was to be expected that numerical model overestimates the pressure losses, because oblique shocks create less pressure loses than normal shocks. However, the numerical results for z = 10∕7 H agree satisfactorily with the experimental mean value. The reason might be the small angle of the oblique shock of about 71 • which does not impact the result significantly. In contrast, for z = 15∕7 H the numerically predicted pressure losses overestimate the experiment. Here, the angle is about 50 • and as seen in the schlieren image, sometimes an additional shock reflection is formed due to the interaction of the jet boundary and the detached shock. Because of that and since the position is at the edge of the barrel shock, which is highly dynamic and constantly moving, it is rather not suited for a quantitative comparison. This can be also seen in the confidence interval, since it is drastically increasing after 4 ms. The reason is that after t=4 ms the pitot probes were inside the jet for some experiments and outside for others.
In summary, the numerical model is consistent with the experimental results and, thus, can be used as a means to quantify the mass flow.

Mass flow rate estimation
Before discussing the outcome of the mass flow estimation based on the numerical model, it is worthwhile to introduce a simple analytical model, which can be used as a comparison.
In effect, the mass flow rate is a function of the state inside the gas reservoir and the smallest cross-sectional area of the system, which is the throat of the gas reservoir, ref Fig. 4a. The mass flow rate ṁ is the product of the minimum cross-sectional area A * , which is a known throat diameter of 4.5 mm for this cold gas inflator, the critical speed of sound c * and the critical density * .
The critical density and speed of sound are a function of the state inside the reservoir assuming an isentropic process and therefore can be expressed as a function of the known volume of the reservoir V, total temperature T 0 , the remaining mass of helium m(t) inside the gas reservoir, heat capacity ratio and the specific gas constant R s . This equation represents a homogeneous differential equation and can be solved with an exponential approach: The required integration constant C 1 is the total mass m 0 of 13.2 g helium stored in the gas inflator. Application of (5) and (6) assumes a constant total temperature inside the gas reservoir and the density depending only on the amount of mass left inside the gas reservoir. In practice, a propagating expansion wave will continuously decrease the total temperature inside the gas reservoir, while reflecting back and forth inside it. For a simple Sod shock tube, Sod (1978), with a constant diameter, the reduction of the total temperature can be calculated. However, since the cross-sectional area of the reservoir is much larger than the throat, where the expansion wave is created, the reduction of the total temperature due to the expansion wave is significantly smaller than for the idealized Sod tube. In order to take this effect into account for the Hypersonic Test Facility Delft (HTFD), Schrijer and Bannink (2010) have derived an analytical approach from the wave theory assuming an instantaneous start of the process and thus a centered expansion wave (further information in Liepmann and Roshko (2002)). This approach has also been used for the design of the Hypersonic Ludwieg Tube Braunschweig(HLB), that has a similar design (see in Wu and Radespiel 2013).
In terms of an order of magnitude estimation, we have used the approach of Schrijer and Bannink (2010) to estimate the decrease of the total temperature due to the initial expansion wave. According to it, the initial expansion wave decreases the total temperature from 293.15 K to approximately 290 K. For comparison, in a straight Sod tube the total temperature would decrease to approximately 220 K. For a perfect gas, the expansion wave's front needs approximately 0.32 ms to reach the back end of the reservoir, where it will reflect and propagate to the front end and reflect again. Therefore, multiple reflections are to be expected within 4 ms. Due to the fan like structure of the expansion wave, the curved geometry of the reservoir's ends and wave superimpositions, the reflection process is complex. Therefore, we have conducted axis symmetrical numerical simulations of the discharge of the cold gas inflator reservoir into a free environment to calculate the decrease of the total temperature in the reservoir for the first 4 ms. For a perfect gas, the total temperature inside the gas reservoir reduced to approximately 263 K and for a real gas (Peng-Robinson equation of state) to 277 K in 4 ms; a reduction of 10.3% and 5.5%, respectively. Since the total temperature is inside the root in Eq. (5), the error is also rooted. Thus, the error due to assuming a constant total temperature in Eq. (5) rises up to approximately 3.2%, respectively, 2.5% in 4 ms. We did not take the effort to correct equation (5); instead, the analytical (6) m(t) = C 1 G ⋅ e −Gt and C 1 = m 0 model is used as an upper limit for the mean mass flow rate and the true mass flow will be smaller. This is even more true since the theoretical model does not account for the ramp-up, compressible effects or accumulation of the gas between the reservoir and the outlet tube, since the geometry is not a straight pipe. Pragmatically we have set the peak mass flow (i.e., the starting point t = 0 for Eq. (6) at 0.46 ms.
As indicated before, the choice of the temperature at the inflow boundary has virtually no effect on the topology of the barrel shock and thus flow field and the pressure. It nevertheless has a strong influence on the mass flow rate, the reason being simply that the temperature scales the density of the fluid at the inflow. It also scales the Mach number and hence the velocity, but the effect is smaller since only the root of the temperature is proportional to the Mach number. Since the fluid temperature itself cannot be measured, the two limiting cases were studied, namely the base case, which assumes = ∞ W m 2 ⋅K , such that the surface temperature of the sensor equals the fluid temperature, and the case with a constant temperature of 293.15 K, assuming that the temperature of the flow is not strongly affected by the igniter and the moving shock front that leaves the gas inflator increasing the total temperature and the expansion waves reducing the total temperature in the reservoir. Since these are the two limiting cases, the true temperature and, hence, the true mass flow are somewhere in between as long as the real temperature is not strongly reduced by expansion waves. In that case, the analytical model represents the upper limit.
The resulting mass flow rates are shown and compared with the analytical model (6) in Fig. 23. As expected, a constant total temperature of 293.15 K gives a similar mass flow rate as the analytical model. Due to compressibility effects and the complex geometry of the gas inflator, deviations are to be expected. The "base case" gives a somewhat slower ramp up process. This is because an initial peak with high temperatures is assumed for this case (ref. Fig. 13), such that the density and, thus, mass flow are reduced. Once the influence of the initial temperature peak has relaxed after some milliseconds, both cases converge to a similar mass flow rate.
We note that all other cases, e.g., = 2 ⋅ 10 5 W m 2 ⋅K , resulted in a mass flow rate in between these two limiting cases, as expected, except the case with "no step": If no step is present the outlet tube, the vortical structure and the shear layer at the edge of the outlet tube that have been discussed in Sect. 4.1 are not present. Consequently, the case "no step" gives a larger mass flow rate, because the total pressure loss generated in the system is smaller and particularly the effective cross-sectional area is bigger. This demonstrates the necessity to resolve the subsonic outlet tube with the step in accompanying numerical simulation.

Conclusion
A method to quantify the discharge of a cold gas inflator is discussed. The method is based on experimental point measurements and the measured data is then taken as input data for a numerical simulation to finally quantify the mass flow rate.
It was found that in the initial phase a moving shock front propagates through the cold gas inflator and eventually leaves through the outlets. This shock front is not particularly relevant for the final flow field, but it has large impact on the pressure measurements due to shock reflections for several milliseconds. During this time, the real pressure of the flow cannot be measured. This error was quantified by introducing a "delay time." Relatively large sensors with screens were found to utilize the clearing effect of reflected shock waves, so that the reflected pressure is cleared within the response time of the sensors. With this approach, the delay time was reduced to the level of the response time of the sensor, which is finally a temporal mass flow uncertainty of about 2.7%-5.4%.
The need to reduce pressure reflections contradicts to the demand for high-frequency resolution. In our case, the frequency resolution was deemed less important and was intentionally traded to reduce the pressure reflection impacts. Finally, frequency resolution was in the range 1.4 to 3.7 kHz. Ensemble averages over multiple gas inflator ignitions are used to analyze the flow field.
The temperature of the flow medium was estimated with fast response heat flux probes and assuming different heat transfer coefficients to derive the gas temperature. Since the heat transfer coefficient can finally not be determined, we can only determine a range of possible temperatures. However, the temperature does not influence the general flow topology or the pressure fields, but only the mass flow rate, especially during the first 2 ms. The experimental data are finally used as a boundary condition for numerical simulations to reconstruct the full flow field including all (other) flow variables and also to find the final mass flow. Some core sensitivities were studied. One of the unexpected findings was a vortical structure emerging from the outlet tube that cause low-frequency disturbances to the interior and outside flow. This vortical structure and the resulting disturbances cannot be found using tank tests, but are clearly visible in the measurements.
In conclusion, the numerical model reproduces the experiments satisfactorily and adequately captures the topology of the flow. The method described herein can be used to quantify the discharge of cold gas inflators-at least for similar designs-or for similar applications involving the discharge from a high-pressure vessel. The method is an alternative to more simplifying tank tests to find the mass flow over time, which is the major input parameter that is required for airbag inflation simulations. Our simulation-based method has also the advantage of delivering a complete data set of all needed physical quantities, in contrast to tank tests. Besides a more detailed quantification of the discharge process, the method delivers a suitable simulation strategy to respect the complex topology and strongly temporally and spatially changing physical quantities for future airbag inflation simulations. Thus, the method and the flow field found herein is viewed as an important step toward a more robust design of airbag inflation.