Non-intrusive, imaging-based method for shock wave characterization in bubbly gas–liquid fluids

A novel experimental imaging-based method is presented for the non-intrusive determination of shock wave characteristics (i.e. shock wave speed and magnitude, and shock-induced liquid velocity) in a bubbly flow solely from gas bubble velocities. Shock wave speeds are estimated by the relative motion between gas bubbles at two locations by splitting the camera field-of-view using a mirror construction, increasing the dynamic spatial range of the measurement system. Although gas bubbles have in general poor tracing properties of the local fluid velocity, capturing the relative dynamics provides accurate estimates for the shock wave properties. This proposed imaging-based method does not require pressure transducers, the addition of tracer particles, or volumetric reconstruction of the gas bubbles. The shock wave magnitude and shock-induced liquid velocity are computed with a hydrodynamic model, which only requires non-intrusively measured variables as input. Two reference measurements, based on pressure transducers and the liquid velocity field by particle image velocimetry, show that the proposed method provides reliable estimates for the shock wave front speed and the shock-induced liquid velocity within the experimental range of 70 <Us<\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$< U_\textrm{s}<$$\end{document} 400 m/s.  


Introduction
Shock waves may cause severe damage to hydraulic systems by transient pressure variations (Ghidaoui et al. 2005;Schmitt et al. 2006). The dynamics of propagating shock waves through single-phase liquids and multiphase liquids in confined geometries have been studied extensively over recent decades (Noordzij and Van Wijngaarden 1974;Jakeman et al. 1984;Kameda et al. 1998;Tijsseling 2007). Imaging-based methods have been employed to study shock wave characteristics of propagating shock waves through two-phase gas-liquid mixtures (Ando et al. 2011;Frolov et al. 2022) and laser-induced shock waves (Vogel et al. 1996).
In case of laser-induced shock waves, shock waves are emitted spherically, by approximation, when a sufficiently high-intensity focused laser beam evaporates liquid locally. In these highly-controlled measurements, shock wave front locations and shock wave speeds are measured by highspeed camera (Lee et al. 2011;Hayasaka et al. 2016;Horvat et al. 2022;Mur et al. 2022). Emitted shock wave pressures can be estimated reliable by the shock wave speed near the emission center (Vogel et al. 1996), which allows to estimate shock wave pressures solely from camera images.
For shock waves propagating through two-phase gas-liquid mixtures, Campbell et al. (1958) were among the first to use optical devices to determine shock wave propagation speeds. They applied two photoelectric cells at different heights and measured the passage of the shock wave by the change in light transmission through the mixture. However, this method is only applicable to relative large void fractions while the bubble dynamics cannot be studied due to the limited spatial resolution. Ando et al. (2011) used high-speed camera images to validate the theoretically predicted wave speed by superimposing the expected pressure wave front on the images. Frolov et al. (2017) investigated the momentum transfer from the shock wave to a bubbly air-water mixture by tracking the motion of bubbles and a polyethylene thread to quantify the bubble and liquid velocities respectively, while obtaining the shock wave speed from pressure transducers. Recently, Gluzman and Thomas (2022) studied unsteady shock wave propagation in aviation fuel cavitation by high-speed imaging and developed an image processing technique, denoted enhanced gradient shadowgraphy, to enhance the appearance of shock waves in the images. Frolov et al. (2022) used high-speed imaging to study the shock wave front with non-reacting and reacting gas bubbles for the application in pulsed detonation hydro ramjet, and extracted the shock-induced bubble velocity from the images. Shock wave pressures are commonly measured by high-frequency pressure transducers. However, disadvantages of pressure transducers include the limitation to measure only at the walls (not inside the flow domain of interest), they are intrusive, and, in case of high temperatures or restricted areas pressure transducers cannot be used at all.
In this research, we propose an imaging-based method that non-intrusively estimates shock wave pressure and velocity, and shock-induced liquid velocity, for shock waves propagating through an aerated liquid, solely from the observed change in bubble velocity and without the need for pressure transducers or tracer particles. Only a single camera (100 kHz) is required. Although cameras are commonly used to image the response of bubbles (see above), to the best of the authors' knowledge no method has been developed yet that determines shock wave pressures directly from the observed change in bubble velocity during the shock wave passages. Shock waves are generated by dropping a free-falling weight on a cylinder, submerged in an elastic tube, whereby the momentum is transferred to the gas-liquid mixture. The shock wave front speed U s is determined from the relative motion of the bubbles between two separated locations. The relative motion of the gas bubbles is used to indicate the arrival of the shock wave front. Since shock waves induce large transient velocity gradients in the bubbly liquid, bubbles may not be considered ideal candidates at first. Indeed, bubbles with sizes in the range of several millimeters have poor fidelity as flow tracers (Mei 1996). However, their response to sudden changes in ambient pressure is highly consistent. We focus on the relative motions between similar bubbles during the initial interaction between the proximal bubble side and the shock wave front. In a second step, a hydrodynamic model is used to estimate the shock wave pressure change ΔP and shock-induced liquid velocity Δu l from U s (from step 1) and the evolution in slip velocity u slip . Additional two-phase PIV and high-frequency pressure measurements are performed in a controlled experiment of a propagating shock wave through a bubbly air-water mixture in a vertical pipe to validate the proposed method.
The structure of this paper is as follows: The next section describes the hydrodynamic model. The experimental setup and data processing of the bubble motions, the determination of the shock wave speed by a split field-of-view (FOV), the reference two-phase PIV images, and pressure sensor data are described in Sect. 3. Results for the measured shock wavefront speeds, and the shock-induced pressures and liquid velocities, are shown and compared with the reference measurements. The final section summarizes our main findings and conclusions.

Relevant equations
The steady shock wave speed U s in a bubbly liquid with initial density 0 in a deformable tube with initial cross-section A 0 can be derived from the quasi one-dimensional mixture conservation equations with fluid-structure interaction and is given by Ando et al. (2011): where subscript 0 represents the undisturbed state prior to the shock wave, and subscript 1 the state for elevated pressures, E the Young's modulus, h the wall thickness, and a 0 the mean radius as defined by (a outer + a inner )∕2 . For a rigid tube ( E → ∞ ), this reduces to U S = √ (P l,1 − P l,0 )∕( 0 (1 − 0 ∕ 1 )) , while finite values of E yield lower shock wave velocities. The shock-induced liquid velocity Δu l = u 1 − u 0 on the aft side of the shock is given by Ando et al. (2011): Although Eq. (1) provides a direct relation between the steady shock wave speed U s and the shock wave pressure amplitude ΔP s = P l,1 − P l,0 , we do not have information on the mixture density 1 and mid-plane cross-section A 1 after the shock, and thus we cannot directly apply this equation to determine the shock wave speed. For very dilute bubble mixtures, the change in dispersed gas phase volume might possibly be determined by resolving the variation in the diameter of individual bubbles. However, here, we assume that volumetric information is not available. Combining Eqs.
(1-3) yields: The last term on the r.h.s. of Eq. (4) can be neglected for small to moderate elevated pressures, as: a 0 ∕(Eh) ≈ O(10 −9 Pa −1 ) , so that the last term O(1 m/s) ≪ O(U s ) , which yields: Typical values for the fraction in the denominator of Eq. (5) are O(10 3 )∕O(10 7 ) ≪ 1 , which yields the original Joukowsky or water hammer equation: The acoustic wave speed c is commonly used as independent parameter in the Joukowsky equation instead of the shock wave speed U s . For finite-amplitude shock waves, (4) ΔP s = 0 U s Δu l this only holds for single-phase fluids or very dilute bubbly gas-liquid mixtures where U s ≈ c , i.e. for unity Mach number, but this is invalid for higher void fractions where U s > c . The shock wave speed U s is determined accurately from the split-FOV imaging method, and the liquid density can be estimated by 0 = g + (1 − ) l ≈ l for low aeration levels , where g is the gas density and l the liquid density. However, Δu l cannot be measured optically as we purposely avoid the use of seeding particles. Hence, Δu l is determined by the bubble slip velocity ( u slip ) between the liquid and bubble velocity, u slip = u l − u b , and coupled to ΔP s via Eq. (6). The equation of motion for isolated deformable non-spherical bubbles is given by Salibindla et al. (2021): where the six terms on the right-hand side represent addedmass, drag, pressure gradient, history, lift and buoyancy forces, respectively. The added mass coefficient is C A , and ∇P is the pressure gradient around the bubble Salibindla et al. (2021).

Model framework
In the present study, we apply a model to determine the properties of a shock wave propagating through a quiescent aerated liquid, confined in a vertical elastic tube. The following assumptions are made: (1) the flow is one-dimensional, (2) the Basset history force F h is ignored in the model, as the effect of the Basset force is negligible for bubble Reynolds numbers Re b = |u l − u b |D b ∕ b > 50 (Takagi and Matsumoto 1996;Magnaudet and Eames 2000;Salibindla et al. 2021), (3) the gas density is neglected since b ≪ l , (4) the lift force term (based on ∇ × u l = 0 ) is neglected, and (5) the thermodynamic behavior is adiabatic, since ∕ R 2 ≈ O(10 −4 ) ≪ 1 ≪ l g ∕R ≈ O(10 2 ) for bubbles with a typical diameter of 5 mm (Van Wijngaarden 2007). Following Kalra and Zvirin (1981), velocities and gravitational acceleration are taken positive in downward direction. The term on the l.h.s. in Eq. (7) is neglected based on assumption (3). For u l ∕ t = U s u l ∕ y (since y = U s t ) and u = u l , the material derivative Du l ∕Dt = u l ∕ t + u ⋅ u l ∕ y can be written as Du l ∕Dt = U s + u l u l ∕ y . Since U s ≫ u l , the term u l ⋅ u l ∕ y can be neglected, and the material derivative reduces to Du l ∕Dt ≈ du l ∕dt , resulting in du l ∕dt − du b ∕dt = du slip ∕dt . For the response of a bubble to a shock wave, Eq. (7) thus simplifies to: where the change of the bubble slip velocity depends on the shock wave pressure, the drag force and, the buoyancy force. The pressure gradient term is estimated by , and t * from the elapsed time between the arrival of the shock wave t 0 and the moment of maximum absolute bubble velocity (see Fig. 4). We assume that the shock wave profile has a constant slope of dP∕dt = (P l,1 − P l,0 )∕t * . The shock wave pressure ΔP s is derived from Eq. (6). An a priori estimate of Δu l in the range from 0 to u b m/s (because | Δu l |<| u b | ) in steps of 0.01 m/s is used to compute ΔP s , while we solve Eq. 8 numerically by a forward Euler method for t = 0 to t * in time steps ( Δt ) of 10 −8 seconds, since the transient event is of the order of milliseconds. For each Δu l , the maximum bubble velocity u b,model = Δu l − u slip is compared to the measured u b and the estimated Δu l that minimized | u b,model − u slip | . This Δu l is selected and used to compute ΔP s . Furthermore, the added mass coefficient is taken as C A ≈ 1∕2 Batchelor (2000). The bubble diameter D b is updated for each time step Δt for changes in the ambient pressure by the shock wave passage via the polytropic gas relation with the initial average bubble equivalent diameter D b,0 following directly from Eq. (8) for du slip ∕dt = 0 and dP∕dy = 0, and assuming the l ≫ b (Assumption 3): where the initial bubble velocity u b,0 is determined from the images prior to the shock wave passage. The drag coefficient C D is modelled as Turton and Levenspiel (1986) : After the shock wave passage, C D is set to 2.6 (Kalra and Zvirin 1981).
In summary, the framework of the model is as follows: 1. estimate D b,0 from the bubble rising velocity before the impact; 2. an a priori estimate by looping over Δu l from 0 to u b in steps of 0.01 m/s; 3. calculate ΔP s from Δu l , aerated liquid density and the measured shock wave speed U s (from Sect. 3.3); 4. numerically solve the nonlinear differential equation (Eq. 8) from t = 0 to t * with time step Δt = 10 −8 seconds by a forward Euler numerical scheme and update D b for each time step; 5. calculate the bubble velocity u b,model by subtracting u b (the outcome of the numerically solved non-linear differential equation) from the induced liquid velocity Δu l ; 6. determine Δu l and ΔP s by min Fig. 1 Schematic overview of the shock wave passage. An example of experimentally measured liquid and bubble velocities during the shock wave passage is shown in Fig. 11 (left), and the bubble response to the shock wave in Fig. 12 3 Experimental set-up

Facility
Experiments were performed to validate the proposed measurement method and model. Figure 2 shows the experimental set-up, which consists of a vertical transparent acrylic plastic tube (Perspex, E = 2855 MPa) with inner diameter D in of 60 mm and 5.0 mm wall thickness. Shock waves are generated by the impact of a free-falling weight ( m = 3.65 kg), released by an electromagnet, on a submerged piston ( m = 0.246 kg), whereby the momentum of the free-falling weight (the impactor in Fig. 2) is transferred to the piston and subsequently to the liquid. Tests are performed with impact velocities of 0.85 and 1.70 m/s, so that the maximum theoretical shock pressures are within the measurement range of the pressure transducers. Specifically, ΔP s = l u imp U s ≈ 680 kPa < 689 kPa, based on the shock wave speed of 400 m/s in the confined tube. With complete immersion of the piston, effects of entrapped gas between the impactor and the liquid surface are largely avoided. No significant friction occurs between the piston and the inner wall as the D piston = 59 mm < D in , which also allows the gas from the aeration to escape freely around the piston. Bubbles are generated by forcing pressurized air through a porous block (Pentair) at the bottom of the tube. To ensure similar-sized bubbles, only bubbles that emerge near the core of the porous block are let into the vertical tube. The liquid phase consists of filtered tap water and the gas phase is ambient air.
A single high-speed CMOS camera (Phantom VEO 640L), equipped with a Nikon 200 mm lens (f # = 11 ), is aligned with the optical construct of one prism and two mirrors (Thorlabs) to capture bubbles in high spatial and temporal resolution (99 μm/pix, Δt = 10 μ s) for determining the shock wave speed in the top and bottom field-of-view (FOV). The centers of the top and bottom FOVs are separated by 720 mm. The recorded images have a 512 × 56-px format, which is the maximum image size (in pixels) at the frame rate of 100 kfps. The upper part of the image ( 256 × 56 px) displays the top FOV, and the lower part ( 256 × 56 px) the bottom FOV (see Fig. 2). Fluctuations in light intensity are negligible during the transient event ( ≈ 1 ms).
Three high-speed pressure transducers (PCI Piezoelectric PCB102) with resonance frequency ≥ 500 kHz are flushmounted in the Perspex tube wall with an equal spacing of 6 D tube . High-frequency pressure transducers are commonly used to determine the shock wave speed in tubes and serve here as reference measurements. Pressure transducers P1 and P3 correspond to the top and bottom FOV respectively, and P2 to the optical axis of the camera. The signals are sampled by National Instruments (NI) LabVIEW 2018 (version 18.0.1f4) and NI Data Acquisition (DAQ) USB-6212 with a sampling rate of 100 kHz, resulting in a Nyquist frequency of 50 kHz. Once the electromagnet releases the impactor, a TTL signal triggers both the acquisition of the pressure (LabVIEW 2018) and images (DAVIS 8.4) simultaneously.
The aeration level is measured by a differential pressure transducer (Validyne model DP45), which compares the hydrostatic pressure of the aerated column with a reference water-only column of the same height. The sampling frequency of the differential pressure transducer is 100 Hz and is only employed before and after the impacts, as the membrane cannot withstand the larger shock wave amplitudes upon impact. The aeration level is calculated by = ΔP m ∕gL( l − g ) , where ΔP m is the pressure difference over the membrane. The sensor is calibrated in situ using different water heights.

Bubble velocity by imaging
Shadowgraphy is used to capture the motion of bubble images. Planar tracking of individual bubbles is compromised by the possibility of overlapping bubble images in the recorded image. Segregation of overlapping bubbles has been addressed extensively in literature (Lau et al. 2013;Kim et al. 2016;Li et al. 2020). However, for cropped image sizes (to enable high-speed recording) it is often not possible to capture entire bubbles (or only at very low spatial resolutions) which makes volumetric reconstruction of individual bubbles even more problematic. Instead, the displacement of the proximal bubble side is tracked by a 1D-correlation, parallel to the direction of the shock wave, and volumes are not reconstructed. Figure 3 summarizes the image processing steps. First, the raw images are corrected for the background image (without bubbles, Fig. 3b). Second, the bubbles are extracted by applying an intensity threshold, where I binary (i, j) = 1 for I(i, j) ≤ I th and zero otherwise. The threshold level is set to 30 % of the average background intensity scale (Fig. 3 left) for all images. In the binary images, the objects are morphologically filled using the MATLAB (version R2018b) function imfill. Proximal bubble edges are detected by evaluating a 3 × 3 px array around each possible line correlation center. For each image with size ( N y , N x ) there are (N y − 2)(N x − 2) possible line correlation centers that are evaluated. Proximal bubble edges are detected at pixel ( i + 1, j ) when the condition I(i + 1, j) − I(i, j) = 1 is satisfied, together with constraints to prevent (1) that vertical line correlations are applied over the bubbles' edges that are parallel to the incoming shock wave, and (2) that no other bubbles are present in the proximity that interfere with the 1D-correlation, to ensure that only valid line correlations are being processed. In the second step, a Gaussian 2D filter is applied to the background-corrected image (Fig. 3a) to remove high frequency spatial noise, and the intensity gradient operator in vertical direction is calculated to emphasize the edges of the bubbles images using the MATLAB function imgradientxy. 1D-correlations are used to calculate the displacement of the proximal bubble Note that bubble 'sides' are not taken into account due to the sharp local curvature. The shock wave arrives from the top edge between images, where these correlations have a window size of 21 × 1 px and are centered at pixels that are indicated as proximal bubble edges (from step 1). Interpolation using a three-point Gaussian fit is applied to determine the displacement at sub-pixel level (Adrian and Westerweel 2011).

Shock wave front speed by imaging
Once a shock wave front arrives at the proximal side of a bubble, the bubble surface deforms by the change in external pressure and it accelerates in the direction of the shock wave. Note that the differences in scales between travelled distance of the shock wave front ( U s ∕f cam ), and bubble diameter ( D b ) cause a trade-off between high temporal resolutions (but cropped images) or high spatial resolutions (but low frame rates). To circumvent this, we split the field-of-view (FOV) of a single camera into two FOVs, separated by distance ΔY . Each FOV focuses on a local flow region with vertical length L y,im = 256 px ≪ ΔY , so that bubble motions on two separated positions can be studied simultaneously. This results in an increased dynamic spatial range (DSR) = 2L y,im + ΔY ∕Δy b,max , where Δy b,max is the bubbles' maximum displacement (Adrian and Westerweel 2011;Westerweel et al. 2013), with a similar dynamic velocity range (DVR) as if the FOVs were not split, i.e. ΔY = 0 . Therefore, by splitting the FOV by distance ΔY , we increased the capability of the measurement system by DVR × DSR = 2L y,im + ΔY ∕ u , where u represents the rms error in the velocity measurement. Splitting one image into two FOVs is done by an optical alignment of a prism and two mirrors. For finite-amplitude pressure waves, the arrival time of the shock wave is uniquely defined by the time instance at which the bubble velocity intersects the zero velocity threshold (see Fig. 4). The shock wave speed U S is determined using the time difference Δt between the moment when the shock wave affects bubbles in each of the two FOVs and defined by U s = ΔY∕Δt . Increasing ΔY has no effect on the spatial resolution of the images, as the spatial resolution only depends on the size of the FOV. The largest measurement uncertainty is introduced by the arrival time of the shock wavefront. With the simultaneous recording of the two FOVs on the same image sensor, the time difference in the passing of the shock wave is uniquely registered. The magnitude of the displacement is less relevant, as long as the bubble deformation can be measured. Incidentally, for the most severe impacts, precursory waves through the frame (to which the mirrors are attached to) may cause oscillations in the velocity profile. These oscillations have negligible effect on the shock wave speed determination as the velocity fluctuations of the precursory waves are much smaller than the change in velocity by the passage of the shock wave.

Reference measurements by two-phase particle image velocimetry (PIV)
Two-phase particle image velocimetry (PIV) measurements in the gas-liquid flow were performed to measure simultaneously the liquid and gas phase velocities during the shock wave passage. Fluorescent tracer particles (Flu-oStar, 1.1 g/cm 3 , 13 μ m diameter, = 579 nm) closely follow the liquid motion (with a Stokes number St ≪ 1 ). A Nd:YLF laser (Litron LDY300 PIV, 4 mJ/pulse, = 527 nm) produces a laser sheet with thickness of 5 mm at a synchronized frequency of 10 kHz that passes through the center of the cylinder and aligned normal to the camera optical axis (see Fig. 5). To avoid overexposure by reflections, the camera lens are equipped with an optical high-pass filter (Schott OG590 with a 590 nm cut-off wavelength). This filter passes the emitted orange light by the fluorescent tracer particles and blocks the green light emitted by the laser. An LED array opposite the camera is used for shadowgraphy to record the bubbles as dark objects. Two distinct intensity levels are applied, which allows the segmentation of the bubble and liquid velocity following the method of Lindken  Fig. 4 Example of the shock wave speed determination by images with the split FOVs. The mean velocity of the bubbles in the top FOV (black curve) and bottom FOV (orange curve) changes abruptly during the passage of the shock wave front. As expected, initially the shock wave front is detected by the bubbles located in the top FOV, and 3.88 ms later by the bubbles in the bottom FOV separated 0.72 m from the top FOV. Because the two images are recorded by one camera, there are no synchronization issues, i.e. 'perfect synchronization' is achieved. The shock wave speed U s is directly calculated from the time difference. In this example the shock wave speed is 186 m/s ( = 0.72 m / 3.88 ms) Interrogation windows that contain at least one pixel with a lower intensity than the threshold (800 counts, see Fig. 5, right) are labeled as 'bubble', and otherwise labeled as 'fluid' (see Fig. 6b). The average liquid velocity is taken as the mean displacement of the 'fluid' interrogations windows. The spatial and temporal resolutions are 71.6 μm/px and 100 μ s respectively. Second, the bright intensity spots of the tracer particle images are removed by setting these pixels to the average background level and then applying a 2D Gaussian smoothing filter using the MATLAB function imgaussfilt to reduce any intensity jumps that may occur. The filtered image is now comparable to the image with split FOV and the same algorithm from Sect. 3.2 is applied for further processing, where the bubble velocity is determined by 1D-correlations. However, larger objects with more than 1000 connected pixels, such as the markers and larger bubbles with equivalent radius of 2.8 mm (36 pixels), are removed.

Reference measurements by pressure transducers
Three flush-mounted high-frequency pressure transducers (P1, P2 and P3) serve as reference measurements to validate the imaging-based measurements and are separated by a distance ΔY = 0.72 m (see Fig. 2). Since the PCB model 102B18 pressure transducers have a rise time ≤ 1 μ s, resonance frequency > 500 kHz, and a useful range of 689 kPa (100 psi), the shock pressures can be fully resolved. Figure 7 shows a typical example of a recorded pressure profile during the shock wave passage. The oscillations in the pressure profiles are caused by the presence of gas bubbles (Fig. 7, right), while absent in case for (nearly) zero aeration (Fig. 7, left), which agrees with the findings by Ando et al. (2011). Precursory waves through the tube material are visible prior to the shock passage. The dotted black line indicates the constant threshold value (in kPa) that is used to determine the  Fig. 6 Processing of the two-phase PIV images: a raw image, the smaller bright intensity spots are the seeding particle images and the darker regions the bubble shadows; b PIV interrogation windows (24 × 24 pixels, 50 percent overlap) that only contain liquid and seeding particles; c selection of bubble images after applying the threshold; d detected bubbles in the gas phase by applying threshold and filling operators to step (c) by isolating the gas phase, the image can be processed starting at step (d1) of Fig. (3) time of arrival of the shock wave front. Linear interpolation is applied to enhance the temporal resolution. As expected, the passage of the shock wave front is first detected by the top pressure transducer (P1) that causes a sudden increase in pressure magnitude (at t = 0). Consecutively, the shock wave front is detected by the center (P2) and bottom (P3) pressure transducer. The reference shock wave speed U s,ref is computed by ΔY∕Δt = ΔY∕(t top − t bottom ) , where the subscripts 'top' and 'bottom' indicate the arrival time at the top (P1) and bottom (P3) pressure transducers, respectively. Pressure transducer P2 is located in between the top and bottom transducer and provides further validation of the measurement. Figure 8 shows the experimental results for the shock wave front speeds by the non-intrusive optical method U s and the pressure transducers U s,ref . In total, 350 impact measurements were performed, grouped into two sets of 170 and 180 measurements with 0.85 and 1.70 m/s impact velocity, respectively. These sets cover 17 and 18 unique void fraction levels (see Table 1), with ten tests per aeration level.

Shock wave speed
Shock wave speeds are varied by the aeration level of the bubbly liquid, where higher aeration levels result in lower shock wave speeds. The overall mean absolute error (MAE) for the smaller and larger impact velocities are 4.2 m/s and 3.8 m/s for impact velocities 0.85 m/s and 1.70 m/s, respectively, where the MAE is calculated by This corresponds to 2.1 and 1.9 % of the average shock speed of 200 m/s. Even for relatively high aeration levels around 3.5 %, the image processing is still able to determine the time instances of shock wave arrival robustly and to provide accurate shock wave speeds. The variance around the black dotted line increases for higher U s , as this region indicates lower aeration level and fewer bubble images are present. This imaging method is limited by the availability of observed bubbles in both FOVs, and thus cannot be applied to single-phase liquids.

Shock pressures
Shock pressures are estimated from the model and compared with the reference pressure transducers (Fig. 9). The mean absolute error (MAE) is computed by (1∕N) ∑ N i �P calc,i − P exp,i � and is 12.6 kPa for the lower impact velocity (159 measurements) and 28.3 kPa for the  higher impact velocity (173 measurements). Deviations in shock pressures are relatively small for lower pressure magnitudes, while the variance increases for larger pressure magnitudes. PIV measurements are performed to investigate the source of the variances in the shock pressure in more detail. The slip velocity between the bubbles and liquid ( u slip ) forms an important mechanism in the model, and two-phase PIV measurements allow for the simultaneous investigation of the bubble and surrounding liquid velocities upon the arrival of the shock front.

Two-phase PIV measurements
PIV measurements were performed to validate the shockinduced liquid acceleration. A typical recording is shown in the supplementary material, both in a raw and processed format (with u l and u b shown as vectors). Necessary changes to the optical arrangement in the experimental setup are made to perform these planar PIV measurements (Fig. 2), but other components of the system remained identical. In total, 54 impacts were recorded ranging from 0 to 1.0 percent aeration and two impact velocities of 0.85 and 1.70 m/s. Note that the shock-induced liquid velocities are clustered (see Fig. 10). Also, the model predictions show a smaller bias, i.e. ū model ∕ū PIV = 1.062 and 0.916, while P model ∕P ref = 1.095 and 0.828, for 0.85 and 1.70 m/s respectively. Since ΔP l = ΔP l ( l , U s , Δu l ) , i.e. the Joukowsky equation, where U s is accurately measured (see Fig. 8) and the maximum error for l is 1.0 percent (the upper range of the void fraction), most deviations in Fig. 9 are expected to originate from the computed shockinduced liquid velocities. Also, the emitted pressures from collapsing and expanding bubbles may have an effect on the overall pressure signal (see Fig. 7). This is supported by Fig. 10, right, where the model consistently under predicts the shock pressure, as single bubble dynamics models (such as the Rayleigh-Plesset model) are not included. On the other hand, Fig. 9 shows that the largest deviations occur for lower aeration levels, where fewer bubbles are present. Also, Fig. 7 (left) shows that no typical higher frequency pressure variations of smaller bubbles are observed in that specific measurement with = 0.084%. Only for the PIV measurements, the shock speed is derived from the reference pressure transducers. The shock wave speed U s cannot be determined by these larger nonsplit images, as the back-light illumination would overexpose the emitted light by the tracer particles. Since U s can accurately be determined by the reference pressure  Fig. 8 Comparison between the shock wave speeds measured by the imaging-based method (vertical axis) and by the pressure transducers (horizontal axis). The green error bars illustrate ± 2 standard devia-tions and are computed based on the realisations of ten measurements within one void fraction group transducers (Fig. 8), this model input variable is expected to contribute marginally to the model output variance.
The synchronized bubble and liquid velocities during the shock passage are shown Fig. 11, left, where the arrows (a-e) refer to the corresponding images in Fig. 12. Time t = 0 ms (b) corresponds to the image in which the average bubble velocity intersects with the zero velocity threshold. As expected, bubbles accelerate initially faster than the surrounding liquid and decelerate afterwards, and compress in response to the elevated shock pressures. Furthermore, the shock wave front passage through the gas-liquid fluid is observed within one image by the acceleration of the liquid (Fig. 12c), from which the shock wave speed is estimated from the liquid velocity profiles along the vertical location of the FOV in Fig. 11 (right). For Δy = 42 mm and Δt = 100 μ s (consecutive frames), the estimated shock wave speed is 420 m/s, which reasonably approximates the measured shock wave speed of 411.4 m/s by the reference measurement. Fig. 9 Comparison between the pressure magnitude (in kPa) obtained directly from the pressure transducers (horizontal axes) and the computed pressure magnitude (in kPa) by the present model for u impact = 0.85 m/s (left) and 1.70 m/s (right). The 7 error bars (green verti-cal lines) represent ± 2 standard deviations, and are constructed by first sorting the data in ascending order for P exp , then divided into 7 groups of circa 25 measurements each, followed by calculating the mean and standard deviations for each group pressures that follow from the reference pressure transducers (horizontal axis), and the calculated shock pressures by the model

Conclusions and outlook
This article presents a novel experimental imaging-based method to determine shock characteristics (i.e. shock wave speed and magnitude, and shock-induced liquid velocity), non-intrusively and solely from gas bubble velocities. The model is validated by pressure transducers and twophase PIV measurements. We conclude that the proposed method with a split FOV is capable of accurately measuring the shock wave speed U s in the range of 70 < U s < 400 m/s (with MAE of 4.2 and 3.8 m/s). Compared to the PIV measurements, the model is able to estimate the shock-induced liquid velocity within a reasonable margin of error, as validated by PIV for two different impact velocities. Measuring the maximum shock wave pressure of transient shocks remains challenging and incurs the largest relative uncertainties. As the model does not include single bubble dynamics, and therefore pressure oscillations by single bubbles (such as the observed highfrequency perturbations on the 'global' pressure profile in Fig. 7 , right) cannot be computed by the current model. This may explain why the maximum shock wave pressure deviates more than the shock-induced liquid velocity from the reference measurements (even thought they are interdependent by the Joukowsky equation). Furthermore, the liquid velocity is mostly affected by the 'global' pressure wave profile, irrespective of the radiated pressure by oscillating bubbles (which does affect the maximum observed pressure), so that the deviation between observed a b c d e Fig. 11 Left: liquid (blue), bubble (green) and slip (purple) velocity profiles during the shock wave passage. Letters (a) to (e) correspond to the images in Fig. 12. Right: instantaneous velocity profile of the fluid tracer particles. Interestingly, the shock wave front is captured in the center of image t + 0.2 ms at Y = 30 mm. The shock wave front is also indicated in Fig. 12b by the purple dotted line (a) t + 0 ms (b) t + 0.2 ms (c) t + 0.3 ms (d) t + 0.6 ms (e) t + 0.9 ms Fig. 12 Five consecutive two-phase PIV images at the shock wave passage. Bright spots indicate fluid tracer particles, while darker spots represent bubbles. The shock wave front is observed in frame (b) at t = 0.2 ms (indicated by the purple dotted line) and corresponds to t + 0.2 ms in Fig. 11 (right). The liquid accelerates at the top of the image, while being still quiescent at bottom. The inset enlarges bubbles (from the red encircled region) and their response to the shock wave passage corresponds to the schematic representation in Fig. 1