## 1 Introduction  
Due to the continuous advancement in computing power and technology, computational fluid dynamics (CFD) has become an established industrial design tool, which helps to reduce design times and improve fluids engineering processes [1]. The present CFD study deals with a particular industrial fluid dynamics problem that has several very important applications, namely, the aerodynamic breakup of liquid droplets under shock wave loading, which is commonly referred to as aerobreakup [2]. The exposure to the high-speed post-shock airflow is responsible for the drift, distortion, and breakup of the original droplet into much smaller fragments [3]. This phenomenon has a crucial role in many fluids engineering applications including, just to mention, raindrop fragmentation in the flow field around supersonic aircrafts. Indeed, in the aerospace industry, a proper aerodynamic design can help greatly reduce the erosion caused by the impingement of small liquid fragments at high relative speeds on the aircraft surfaces [4]. Theoretically, the physics of aerobreakup is controlled by the Ohnesorge number (Oh) of the droplet, which compares the effects of liquid viscosity and surface tension, and the Weber number (We) at freestream conditions, which describes the ratio of disruptive inertial to restorative interfacial forces. Previously, based on the definition of various different breakup modes, the interpretation and identification of the aerobreakup mechanism have been recently reviewed and simplified by [PERSON], who introduced only two principal modes of fragmentation that are based on the underlying physical mechanism of the dominant interfacial instabilities [5]. Basically, at low Ohnesorge numbers (Oh \(\ll 1\)), the process is practically governed by the Weber number. For low Weber numbers (\(10<\text{We}<10^{2}\)), the Rayleigh-Taylor (RT) instability represents the main driving mechanism for the aerobreakup, which is referred to as the RT piercing regime, corresponding to a moderate fragmentation. Differently, at high Weber numbers (We \(>10^{3}\)), the breakup process is governed by the Kelvin-Helmholtz instability of the stretched liquid filaments that are formed at the periphery of the deforming droplet, with the subsequent liquid sheet stripping and thinning, and the extended entrainment of liquid mist. This terminal mechanism, which is referred to as the shear-induced entrainment (SIE) regime, corresponds to a strong fragmentation, where capillary forces and turbulent mixing play an important role [2]. Several experimental studies have been conducted to investigate the relevant physical mechanisms and the main features of the droplet aerobreakup. In practical experiments, a traveling planar shock front is often reproduced through shock tube devices, where uniform airflow conditions are established around the water body [6]. In these experiments, sophisticated visualization techniques are used to analyze the deformation and breakup of droplets under shock wave loading [7]. In particular, the physical mechanisms determining the dynamics of the dispersed particles that are forming as a result of the chaotic fragmentation in the SIE regime have been investigated by different research groups [8; 9]. However, the experimental quantification of the characteristic sizes and velocities of the liquid fragments inside the mist remains very challenging. From a numerical point of view, since the droplet aerobreakup under the SIE regime is characterized by a very broad range of spatial and temporal scales, the accurate prediction of the droplet deformation and fragmentation is particularly demanding in terms of computing power. Due to the high computational cost of fully three-dimensional (3D) simulations, several CFD studies have been conducted in two spatial dimensions, by considering that the shear-induced breakup of a cylindrical liquid column [10; 11; 12].