Computational Evaluation of Shock Wave Interaction with a Cylindrical Water Column

: Computational ﬂuid dynamics was employed to predict the early stages of the aerodynamic breakup of a cylindrical water column, due to the impact of a traveling plane shock wave. The unsteady Reynolds-averaged Navier–Stokes approach was used to simulate the mean turbulent ﬂow in a virtual shock tube device. The compressible ﬂow governing equations were solved by means of a ﬁnite volume-based numerical method, where the volume of ﬂuid technique was employed to track the air–water interface on the ﬁxed numerical mesh. The present computational modeling approach for industrial gas dynamics applications was veriﬁed by making a comparison with reference experimental and numerical results for the same ﬂow conﬁguration. The engineering analysis of the shock–column interaction was performed in the shear-stripping regime, where an acceptably accurate prediction of the interface deformation was achieved. Both column ﬂattening and sheet shearing at the column equator were correctly reproduced, along with the water body drift.


Introduction
The aerodynamic breakup process of liquid droplets induced by the interaction with plane shock waves is of crucial importance for many industrial gas dynamics applications. The latter include, just to mention a couple of them, supersonic combustion airbreathing jet engines (scramjets) and rain impact damage to supersonic aircraft [1,2]. The shock wave interaction with a water droplet actually represents the initial stage of the breakup process induced by the high-speed gas stream. Since this process plays an important role in the droplet breakup, a number of studies have been conducted to investigate the shock/droplet interaction for correctly interpreting the mechanism of aerobreakup. In particular, this process has been the subject of several experimental studies employing shock tube devices, wherein a traveling planar shock wave is reproduced, with uniform gaseous flow conditions being established, e.g., [3,4]. In these experiments, the shock front passes over the droplet and causes its deformation and successive breakup.
The physical mechanism of aerobreakup is controlled by the Ohnesorge number (Oh), which compares liquid viscosity to surface tension effects, as well as the Weber number (We), which compares disruptive aerodynamic forces to restorative surface tension ones. However, at low Ohnesorge numbers (Oh < 0.1), the influence of liquid viscosity can be neglected, and the process is governed essentially by the Weber number [5]. The interpretation and identification of breakup mechanisms have been recently reviewed, based on experimental studies using highly resolved images of liquid droplets suddenly exposed to supersonic gas streams [6]. According to this reclassification, at low Weber numbers (10 < We < 10 2 ), Rayleigh-Taylor instability has to be considered the main driving mechanism for the aerobreakup, which is referred to as Rayleigh-Taylor piercing (RTP). Differently, at high Weber numbers (We > 10 3 ), significant shear-induced motions exist, and the breakup process is governed by the instability of the stretched liquid sheets that are formed at the periphery of the deforming droplet. This regime is referred to as shearinduced entrainment (SIE) [7]. Another non-dimensional parameter that is sometimes considered when examining this process is the ratio between liquid and gas densities.
As for other challenging engineering applications, computational fluid dynamics (CFD) can be effectively employed to predict the early stages of the aerobreakup of a water droplet due to the impact of a traveling shock wave. In the context of preliminary analyses, the shock-droplet interaction can be evaluated by considering a water cylinder with a circular cross-section in the high-speed flow behind a planar shock front. The cylindrical geometry of the water body can be efficiently modeled using a two-dimensional computational domain that provides shorter turnaround times compared to three-dimensional simulations. Several numerical studies exist in the literature dealing with two-dimensional shock-column interactions, e.g., [8][9][10], where the sheet-thinning process was simulated and the column deformation and drift were measured. Usually, such works focus on the early stages of aerobreakup and neglect the effects of both surface tension and molecular viscosity. The two-dimensional approach is also justified based on experimental investigations of liquid columns that allow for easier visualization of the wave structures [11,12]. Indeed, in the early stages of the interaction process, the deformation of two-dimensional liquid columns has been found to be very similar to that of three-dimensional spherical droplets [13].
Compressible multi-phase solvers based on a five-equation model [14] are often developed to computationally study shock-water column interactions, where the interface is modeled using volume fractions. The governing equations consist of two continuity equations for the two different phases, the mixture momentum and energy equations, and the volume fraction advection equation, e.g., [8,15]. As far as flow turbulence modeling is concerned, both under-resolved direct numerical simulation (DNS) and large eddy simulation (LES) approaches are commonly followed, e.g., [4].
The main goal of the present work was the computational evaluation of the initial stages in the interaction process between the air shock front and the cylindrical water column following a relatively light approach, which was the unsteady Reynolds-averaged Navier-Stokes (RANS) modeling. The mean turbulent flow in a virtual shock tube device was simulated by supplying the calculation with a suitable turbulence closure model. The interaction process was simulated in the shear-stripping regime that was at relatively high Weber and low Ohnesorge numbers. Numerical calculations were conducted by employing one of the CFD solvers that are commonly and successfully used for building virtual wind tunnels in industrial gas dynamics research, e.g., [16,17]. The compressible flow governing equations were solved by means of a finite volume (FV)-based numerical technique, while the transient tracking of the air-water interface was approximated on the fixed FV mesh using the volume of fluid (VOF) method [18]. The present CFD modeling procedure was validated against reference data that were provided by both high-resolution numerical solutions and experiments.
The rest of this manuscript is organized as follows. After describing the overall computational model in Section 2, the results of the numerical simulations are presented and discussed in Section 3. Finally, in Section 4, some concluding remarks are drawn.

Computational Model
In this section, after introducing the industrial gas dynamics application under study, the details of the CFD model are provided, including the method used for tracking the liquid-gas interface and the CFD solver settings.

Case Study
According to the two-dimensional approach, the simplified geometry of the shock tube device is represented by the rectangular domain sketched in Figure 1, where the planar shock, and thus the post-shock air flow, moves from left to right. The reference coordinates system (x, y) was chosen with the x-axis aligned with the streamwise direction, while the origin corresponded to the leading-edge of the water cylinder in its initial position. Based on previous numerical simulations conducted for the same flow configuration, the spatial domain size was chosen such that −8 < x/d 0 < 16 and |y|/d 0 < 10, where d 0 stands for the initial diameter of the cylindrical water column, which represents the natural reference length. The shock tube was divided into two parts by means of a virtual cross diaphragm, which was initially located at x/d 0 = −2 (upstream of the column). In the following, the gas conditions at the right and left sides of the tube are indicated by the subscripts 1 and 4, respectively. The two tube sections were initialized with the same temperature (T 4 /T 1 = 1) and very different pressure and density levels (p 4 /p 1 = ρ 4 /ρ 1 1), with the ideal gas air being at rest on either side (V 4 = V 1 = 0). Starting from these initial conditions, a planar shock wave develops and travels towards the right side (driven section), while a set of expansion waves propagate towards the left one (driver section). The virtual diaphragm, which corresponds to the contact surface between the two different density regions, travels in the same direction of the shock, though at a lower velocity. Given and knowing the pressure level p 1 on the right side, as well as the compression ratio across the shock β = p 2 /p 1 , corresponding to the prescribed Mach number, the thermodynamic parameters to be imposed on the left side can be analytically determined. In fact, combining the moving shock wave equations with the expansion region relations, it is possible to derive an analytical expression for the ratio p 4 /p 1 , that is: where γ = 1.4 was assumed for the specific heat ratio. Therefore, based on the above estimation, the initial conditions in the driver section were uniquely determined, e.g., [19]. In this study, the shock tube flow was simulated at a Mach number of 1.73, which corresponds to one of the cases recently examined in [20], where the water cylinder was immersed in supersonic air crossflow. The thermodynamic parameters that were imposed are summarized in Table 1, while V S = 593.7 m/s and V 2 = 329.3 m/s represent the constant velocities of the traveling shock wave and post-shock air flow, which were determined from theoretical considerations. The high-speed airstream behind the shock represents the ambient conditions to which the water column is exposed, which are indicated by the subscript 2. For clarity, the water column and post-shock air flow parameters of interest are tabulated in Table 2. The compressible aerodynamics of the interaction process is governed by the flow Mach and Reynolds numbers, which are: and: respectively. Moreover, as discussed in the previous section, the aerobreakup of the water column induced by the impact of the shock wave is controlled by the Ohnesorge and Weber numbers. Based on the liquid water properties and the surface tension σ, which was assumed constant, these dimensionless parameters are defined as: and: respectively. Given the current values of the physical parameters that were involved, these four characteristic numbers took the values provided in Table 3. Owing to the very low Ohnesorge and the relevant Weber numbers, the current breakup configuration belonged to the SIE regime, where the gas was expected to go around the liquid mass producing a surface layer peeling-and-ejection action [7]. Capillary effects could be neglected with respect to the inertial ones in the early stages of the aerodynamic breakup, where shear stripping was the dominant mechanism. Moreover, based on the relevant Reynolds number, viscous effects could be considered negligible. However, viscous effects and surface tension were not ignored in the present work, different from most numerical simulations in the reference literature, e.g., [9,20].

CFD Analysis
The present CFD analysis employed the VOF method for the transient tracking of the liquid-gas interface during the interaction process, which allowed approximating the boundary between the two different phases on the fixed numerical mesh. According to this classical methodology, volume fractions were introduced to represent the volume occupied by each phase in each computational cell. The VOF technique was shown to be more efficient and flexible than other methods when treating the complex interface between two immiscible fluids [18] and is often employed to study the interaction of water droplets with shock waves, e.g., [4,12].
Practically, the compressible Navier-Stokes equations were written for volume averaged fields that were shared by the two different phases. These equations were solved for an effective fluid, whose averaged properties were evaluated according to the volume fractions. For example, the averaged density was determined by: where ρ g represents the variable gas density and α stands for the volume fraction of liquid water (secondary phase). The field variable α was evaluated throughout the computational domain by solving the associated continuity equation, while the volume fraction of the ideal gas air (primary phase) was trivially computed. In any computational cell, depending on the local volume fractions, the variables and properties were representative of either one of the two phases or a mixture of them. Among the advantages of the VOF method, the effect of surface tension forces along the air-water interface can be simply simulated. The present analysis was carried out by employing the continuum surface force (CSF) model proposed in [21], with the momentum equation being supplied with an additional source term due to surface tension, as done in similar studies, e.g., [15].

CFD Solution
The mean turbulent flow in the virtual shock tube device was simulated by solving the unsteady compressible RANS equations. The shear-stress transport (SST) k-ω two-equation eddy-viscosity model was used for turbulence modeling, owing to its proven suitability for complex industrial engineering applications [22]. The flow governing equations, which were not reported here for brevity, can be found, for instance, in [23]. The computational evaluation was performed by using the industrial solver Ansys Fluent 19, which has been successfully employed in analogous works, e.g., [16], as well as in previous CFD studies using RANS models by the same research group [24][25][26].
The present pressure-based solver utilizes the FV approach to discretize the governing equations and approximate the unsteady mean flow solution, where the conservation principles are applied over each computational cell, e.g., [27]. The system of integral balance equations for mass, momentum, and energy was completed by employing two equations of state to model the density variations of the two different fluids. Specifically, along with the ideal gas law for the gas phase, the following Tait's equation of state for the liquid phase was considered, e.g., [28], where p 0 and ρ l0 represent the reference pressure and density levels, while B and κ are constant parameters. Following [12], the pressure-like parameter B and the so-called adiabatic index κ were set to the values of 305 MPa and 6.68, respectively. This way, given the maximum pressure levels that were expected in the simulations, the liquid density can be considered practically constant for the present engineering analysis. The FV mesh was suitably refined in the flow region that was expected to be occupied by the deforming water column, including the near wake. Three different meshes with varying resolution were tested, as reported in Table 4, where h min stands for the minimum mesh size. The total number of FV cells ranged between two and five-hundred thousand, where the highest local resolution resulted in being about 60 and 430 cells per cylinder diameter, respectively. Actually, the computational domain Ω that was used corresponded to half the spatial domain sketched in Figure 1, where a symmetry boundary condition at the centerline (y = 0) was imposed. In fact, according to previous experimental visualizations, the flow can be assumed symmetric at the early stages of the interaction process, e.g., [11]. As for the other boundaries, wall boundary conditions were employed for simulating the solid walls of the virtual shock tube device. Note that the present boundary conditions coincided with those adopted in other similar simulations, e.g., [12]. Finally, the constant time-step of ∆t = 0.1 µs was used for the explicit transient calculation.

Shock Tube Flow
The computational setup described in Section 2.1, without the presence of the water column, practically corresponded to Sod's shock tube problem [29], which admits an exact solution and has become a standard test case for compressible flow solvers. Therefore, while demonstrating the air flow in which the water column was immersed, the initial simulation of the shock tube device also represented a preliminary verification and validation test which was satisfactorily completed, as was done in [8], for example.
The robustness of the present CFD model and the accuracy of the transient numerical solution of the shock tube problem were assessed through comparison with the analytical discontinuous solution [19]. Here, the resolved flow was illustrated by looking at the profiles across the shock tube of some flow field variables at the moment, say t = t 0 , when the incident shock impacted the column. The pressure, density, temperature, and streamwise velocity profiles along the tube are reported in Figure 2, for the three different CFD solutions with varying resolution. Note that, as happens in real flows, theoretical discontinuities were smoothed out by viscous effects.
In Table 4, the results in terms of the traveling shock wave and contact surface velocities are tabulated. By inspection of these results, grid convergence was practically demonstrated, since the matching between numerical and analytical solutions improved with increasing mesh resolution. In Figure 3, a close-up view of the FV mesh in the space region occupied by the water body was shown for the intermediate resolution case, that is CFD II. In the following section, the computational evaluation of the shock wave interaction with the cylindrical water column is presented and discussed for the two cases CFD II and CFD III.

Shock-Column Interaction
The complete setup described in Section 2.1 was employed to simulate the air shock wave interaction with a water droplet. As a consequence of the impact of the shock front, distortion and successive breakup of the water column were observed. In the following discussion, as is usually done in the relevant literature [5,30], the time variable was nondimensionalized as follows: where ε = ρ l /ρ 2 corresponds to the water/air density ratio.

Qualitative Analysis
In this section, the qualitative analysis is presented for the solution with the finest resolution, that is CFD III. At the very early stages of the interaction process, the passage of the shock wave did not produce appreciable cylinder deformation. As observed in [31], there existed a reaction time, which depended on both the drop size and air stream velocity, during which the water column appeared undistorted, with the air flow field resembling the flow past a circular cylinder.
This fact was clearly demonstrated by looking at the contours of the volume fraction of water depicted on the left side of Figure 4, where the time history of the interaction process corresponding to the non-dimensional time interval 0 < t * < 0.03, immediately after the impact, is illustrated. For clarity, the different contour maps were drawn in the reduced spatial domain defined by −1 < x/d 0 < 2 and y/d 0 < 5/3. Note that, due to the finite numerical resolution, the interface between the two immiscible fluids appeared to be slightly diffuse. Apparently, the column showed no change in shape after being struck by the air shock. Moreover, by looking at the mean pressure, temperature, and streamwise velocity contours, both original and reflected moving wave fronts were evident. The impacting wave reflection that was initially regular became a Mach reflection at a certain angle, as observed at the last time instants. The stagnation flow region behind the reflected wave was particularly evident by looking at the velocity contours, while the pressure contours clearly revealed the wave that was transmitted inside the water body. As was illustrated by the temperature contours, the liquid phase was practically maintained at the constant temperature of the driven section.  As the interaction process continued, after the reaction time elapsed and interfacial instabilities appeared, the water column began to be distorted due to the dominant action of pressure forces, while still maintaining its coherence. Specifically, the water cylinder was flattened in the streamwise direction, which corresponded to a decreasing centerline width and an increasing cross-stream breadth. The deformation of the water column is illustrated in Figure 5, where the pressure and streamwise velocity contours are reported for a later time interval in the spatial domain −0.01 < x < 0.05, y < 0.015. For clarity, the volume fraction field is illustrated in the limited range defined by −0.5 < x/d 0 < 4. It is worth stressing that the present qualitative analysis for the different phases of 238 the interaction process in terms of unsteady mean flow is fully consistent with the results 239 provided by more sophisticated high-resolution numerical simulations, e.g. [8,10,20]. By inspection of the volume fraction contours, the drift and distortion of the water column were evident. In fact, the streamwise flattening of the liquid body and the subsequent formation of a peripheral tip were clearly observed for the first temporal phase, say 0.17 < t * < 0.45. Then, the deformation progresses and the incipient breakup of the water column were observed at later time instants, say 0.55 < t * < 0.85. Specifically, as the peripheral tip was drawn downstream by the surrounding flow, due to the dominant action of shear forces, the liquid sheet became thinner and thinner until, finally, the breakup process began. The water column was continuously eroded at the periphery, with the appearance in the near wake of microdroplets stripped from the parent body. Basically, these pictures confirmed the qualitative features of the shear-stripping breakup mechanism, as was expected at high Weber numbers [7]. During the successive stage, the water column would disintegrate into fragments distributed widely in the flow field, which was out of the scope of this study.
It is worth stressing that the present qualitative analysis for the different phases of the interaction process in terms of unsteady mean flow was fully consistent with the results provided by more sophisticated high-resolution numerical simulations, e.g., [8,10,20].

Quantitative Validation
The quantitative validation of the proposed CFD model was performed by examining some geometrical parameters of the deformed water column versus time, while making a comparison with reference numerical data [20], as well as experimental findings [11]. Here, the instantaneous shape of the distorted column was determined as that associated with the isoline at α = 0.5 for the volume fraction of water. In order to demonstrate the effect of increasing the numerical resolution, the results corresponding to solutions CFD II and CFD III are presented.
The time history of the decreasing centerline width of the column w, normalized by the initial diameter d 0 , is shown on the left side of Figure 6. The two different lines that were reported for the reference numerical solution corresponded to α = 0.25 (blue dashed line) and 0.99 (red dashed-dotted line), whereas green dotted and black solid lines were used for CFD II and CFD III, respectively. Apparently, a good agreement was achieved between the present solutions and both references. When making an analogous comparison for the increasing cross-stream breadth b, which is depicted on the right side of the same figure, CFD analysis appeared to provide overestimated results with respect to the experimental findings. However, at the present relatively high Mach number, the reference numerical solution provided similar results. Furthermore, the deformation of the coherent water body was examined in terms of cross-section area A, whose time history is given on the left side of Figure 7. Again, a good agreement with both references was obtained for the decreasing area, which is shown normalized by the initial circular section area. It is worth stressing that the comparison with experimental measures, which were obtained using holographic interferograms, was affected by the inherent uncertainty associated with the definition of the boundary of the deforming water column. As observed in [20], while it was not possible to unambiguously identify the method used in the experiments [11] to define the boundary, the presence of a diffuse interface provided by the numerical methods further worsened the issue.
The time evolution of the leading-edge drift ∆x le is examined on the right side of Figure 7, wherein CFD data appeared overestimated with respect to experimental visualizations, again according to what was found in [20]. However, as observed in [7], the leading-edge drift did not coincide with ∆x cm , that is the center-of-mass (CM) drift. The time-dependent position of the CM is defined by: which inherently involves only the flow region that is actually occupied by the liquid phase, and the CM velocity in the streamwise direction can be analogously evaluated, namely: By inspection of Figure 8, wherein the time histories of the CM drift and velocity are given, the present CFD data appeared to be slightly overestimated with respect to reference numerical results [20], with the positive effect of increasing the numerical resolution being evident.
In fact, while the present resolution sufficed to practically achieve converged evolution for the leading-edge drift, a finer FV grid should be used to better capture the CM drift, according to what was found in [10]. The streamwise acceleration a cm of the CM was obtained by differentiating the above discrete velocity data using second-order finite difference approximations, for both the present CFD and reference numerical solutions. Looking at the time evolution of this variable, which is reported suitably normalized in Figure 9, the agreement between the different solutions was acceptable. Note that a similar result also holds for the aerodynamic drag coefficient. In practice, since the above integral parameters resulted in being quite insensitive to small-scale interface structures, which would not be resolvable by means of the current approach, the overall evolution of the water column drift was well captured. Finally, based on the acceptable results that were obtained in terms of both qualitative analysis and quantitative validation with respect to reference solutions, the proposed CFD model was shown to be able to capture the salient features of this particular gas dynamics application, without being computationally cumbersome. This highlighted the benefits of the present relatively simple and straightforward approach.

Conclusions
The present study was intended as a proof-of-concept, namely the preliminary development and demonstration of an industrial CFD-based prediction tool to simulate the distortion and breakup of water droplets induced by the passage of a shock wave. Following a two-dimensional unsteady RANS approach, CFD analysis of the early stages of the plane shock wave interaction with a cylindrical water column was performed. The computational method was able to reproduce shock tube flow conditions, as well as distortion and incipient breakup features observed in experiments. An acceptably accurate prediction of the interface deformation was achieved. Both stages of the process, which were water cylinder flattening and sheet shearing at the droplet equator, were correctly reproduced. In fact, the present results were in good agreement, both qualitatively and quantitatively, with reference solutions provided by high-resolution numerical studies.
Definitely, there remains a possibility of developing more sophisticated computational models for this particular gas dynamics application, depending on the level of accuracy that is required and the computational cost that is affordable to industrial researchers. Following the pioneering studies in [15,32], for example, one chance would be using wavelet-based adaptive numerical methods [33,34]. This innovative methodology, which exploits the wavelet transform to dynamically and automatically adjust the grid resolution to the local flow conditions, has been recently further developed for wall-bounded supersonic flows [35]. Wavelet-based adaptive unsteady RANS modeling procedures [36][37][38] could be possibly explored in the present context. Moreover, in order to simulate the appearance and dynamics of microdroplets during the breaking-up process, VOF-Lagrangian hybrid methods [39] could be used in the framework of either the present industrial CFD modeling or the wavelet-based approach [40].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

RTP
Rayleigh-Taylor piercing SIE shear-induced entrainment CFD computational fluid dynamics DNS direct numerical simulation LES large eddy simulation RANS Reynolds-averaged Navier-Stokes FV finite volume VOF volume of fluid CSF continuum surface force SST shear stress transport CM center-of-mass