{'fraction_non_alphanumeric': 0.029832459539548665, 'fraction_numerical': 0.010485525416001824, 'mean_word_length': 4.312897366030881, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compressible air flow through a collapsing liquid cavity

We present a multiscale approach to simulate the impact of a solid object on a liquid surface: upon impact a thin liquid sheet is thrown upwards all around the rim of the impactor while in its wake a large surface cavity forms. Under the influence of hydrostatic pressure the cavity immediately starts to collapse and eventually closes in a single point from which a thin, needle-like jet is ejected. Existing numerical treatments of liquid impact either consider the surrounding air as an incompressible fluid or neglect air effects altogether. In contrast, our approach couples a boundary-integral method for the liquid with a Roe scheme for the gas domain and is thus able to handle the fully \emph{compressible} gas stream that is pushed out of the collapsing impact cavity. Taking into account air compressibility is crucial, since, as we show in this work, the impact crater collapses so violently that the air flow through the cavity neck attains supersonic velocities already at cavity diameters larger than 1 mm. Our computational results are validated through corresponding experimental data.

The sequence of events as a circular disc of 2 cm radius impacts a water surface at 1 m/s: (a) Immediately after impact a crown splash is thrown up into the air and an impact cavity forms below the free surface. The moving disc draws air into the cavity. (b) Hydrostatic pressure pushes the cavity together and the direction of air flow reverses. (c) Eventually the cavity closes in a single point. (d) After closure two violent jets are ejected up-and downwards from the closure location. The blue and red lines represent the free surface and the disc, respectively, from the present numerical simulation with air in (a)-(c). For the jetting in (d) the single-phase simulation from [15] is shown.
In this work we report a seemingly simple and harmless situation which nevertheless requires to model the dynamics of the gas phase in a fully compressible way: a circular disc of 2 cm radius which impacts on a liquid surface with a speed of 1 m/s. Figure 1 illustrates the sequence of events during the disc impact extracted from high-speed video images [16,18] and compared to the results of our simulations. Upon impact, first a thin liquid splash is thrown up all around the circumference of the penetrating disc. In the wake of the impactor a large cavity is created which subsequently starts to collapse due to the hydrostatic pressure of the surrounding liquid. When the cavity closes about half-way down its length two very fast and thin jets are observed shooting up-and down from the closure point [15,[19][20][21].
In the beginning of the process, obviously, air is drawn into the cavity by the moving disc. At a later stage, however, this inward flow is counteracted by the shrinking of the cavity volume itself and the direction of air flow is not a priori clear. We will show that in the competition between cavity expansion (just above the disc) and shrinking (around the neck), eventually the shrinking becomes dominant. Accordingly, the air flow through the neck reverses and air is pushed out of the cavity. The collapse is so violent that the air stream can attain supersonic speeds. To handle this situation, we implement an axisymmetric boundary-integral method (BIM) to simulate the motion of the liquid surface which is two-way coupled with a fully compressible Roe solver for the highly unsteady gas flow.
In Section 2 we will introduce briefly our BIM and Roe implementations and then describe in some detail the two-way coupling between the BIM for the liquid and the Roe scheme for the gas domain. Sections 3.1 and 3.2 show the main physical results extracted from our simulations and Section 3.3 briefly summarizes the experimental validation carried out in [18]. Section 4 concludes the article.

Numerical methods
In the impact process viscous effects are negligible as can be seen by estimating the Reynolds numbers for gas and liquid as Re g,l = V 0 R 0 /ν g,l . Here, V 0 = 1 m/s is the impact speed, R 0 = 2 cm is the disc radius, and νg = 1.46 · 10 −5 m 2 /s and ν l = 1.12·10 −6 m 2 /s are the dynamic viscosities of gas and liquid at 15 • C, respectively. Both Reg and Re l are larger than 10 3 and viscosity is thus negligible. Furthermore, we showed in earlier works [16] that only a negligible amount of vorticity is present in the system. We can thus assume the flow to be inviscid and irrotational as is required for the boundary-integral method and the applicability of the inviscid Euler equations.
Our simulation is split in two stages. In the first, incompressible stage, a twofluid boundary integral method is used to simulate the first part of the impact and cavity collapse where gas velocities are of the order of the disc velocity. In the second, compressible stage, we use a single-fluid BIM for the liquid coupled to a Roe scheme to solve the compressible Euler equations in the gas domain [22].
The transition between both stages is fixed at the moment that the air flow through the cavity neck reverses and the gas axial velocity at the neck is zero, i.e. u neck = 0. 1 At the transition moment the liquid cavity has an elongated shape with a neck located roughly at the middle as can be seen in Fig. 1 (b). Furthermore, the gas flow is to a very good approximation one-dimensional directed along the symmetry axis of the cavity (as can be verified by flow profiles obtained from the two-fluid BIM shown in Section 3.1). This makes the situation reminiscent of gas flow through a convergingdiverging de-Laval nozzle frequently encountered in aerodynamics [23]. In this spirit, we remove the inner potential fluid during the compressible stage and replace it with a compressible gas described by the one-dimensional Euler equations. To integrate Euler's equations in time we use the well-known scheme due to Roe [24]. This description is valid along most of the cavity. Above the initial free surface as well as in a small zone above the disc the air dynamics can be neglected as will be described in Section 2.3. The pronounced difference to the standard nozzle situation, however, is that in our case the "nozzle" geometry is determined by a liquid interface which is changing rapidly in time and thus creates a highly unsteady gas flow.
The two-way coupling between the gas and the liquid domains is accomplished via (i) the interfacial shape and its instantaneous velocity which is provided by the BIM and serves as an input into the gas solver and (ii) the gas pressure which is obtained from the solution of the Euler equations and serves as a boundary condition for the BIM. Above the location of the initial water level the surface pressure of the BIM remains atmospheric. Our combined BIM/Euler method has the advantage that we retain the computational efficiency of the BIM which requires only the modeling of the boundary and not of the entire liquid domain in contrast to, e.g. compressible air/water treatments based on the Volume-of-Fluid method for sloshing tanks [25,26].
All quantities are non-dimensionalized with the disc radius R 0 = 2 cm and the impact velocity V 0 = 1 m/s. As a third quantity for non-dimensionalization we use the density of water ρ l = 998.23 kg/m 3 for the equations concerning the liquid domain and the density of air (at rest under atmospheric pressure) ρg = 1.2 kg/m 3 for the gas equations. Since viscosity is neglected we have four dimensionless parameters which are the Froude number, the Weber number, the Euler number for the liquid, and the Euler number for the gas defined as: with the atmospheric pressure pa = 101.3 kPa, the acceleration of gravity g = 9.81 m/s 2 , and the air/water surface tension σ = 72.8 mN/m. We use cylindrical coordinates r and z with z = 0 at the height of the initial free surface as is appropriate for our axisymmetric setup. In all computations the disc is assumed to have zero thickness. The boundary-integral method for a single fluid (liquid) is briefly sketched in 2.1.1 and its extension to two fluids (liquid and air) is given in 2.1.2. Some important aspects of our specific BIM implementation are described in 2.1.3. The Roe solver is, again briefly, presented in Section 2.2. In Section 2.3 we describe the coupling between the gas and liquid domains during the compressible stage.

Boundary integral method for a single fluid
If liquid flow is inviscid and irrotational as is the case in our setup [15,16] the flow field v can be described as the gradient of a scalar potential φ v = ∇φ (5) which satisfies Laplace's equation throughout the liquid domain From Eq. (6) the boundary-integral equation can be derived using Green's identities [27]: with S denoting the boundary, n the normal vector pointing out of the liquid domain, and φn = n·∇φ the normal derivative of the potential along the boundary. Equation (7) expresses the potential φ at an arbitrary point r inside the domain (then β = 4π) or on its edge (then β = 2π) merely in terms of quantities that are defined on the surface of the liquid domain. This has the major advantage that it is sufficient to evolve the surface quantities φ and φn effectively reducing the computational problem by one spatial dimension.
If on every point along the interface either φ or φn is known, Eq. (7) can be solved for the missing quantity such that afterwards both φ and φn are available along the entire boundary. With this the interfacial velocity can be computed as v = ∂φ ∂s t + φnn (8) with ∂/∂s denoting the tangential derivative and t the tangential vector along the surface. The numerical procedure for this is sufficiently well described in the literature [27][28][29] and the details of the present implementation will be given in Section 2.1.3. The boundary conditions for solving Eq. (7) are provided as follows: on the disc the liquid is required to follow the disc's motion meaning that φn = −1 in nondimensional coordinates. On the air/liquid interface the potential φ is specified by integrating Bernoulli's equation: with D/Dt = ∂/∂t+v·∇ denoting the material derivative, C = ∇·n the local curvature of the interface, and p = pg/pa the dimensionless gas pressure. Here we assume that the pressure at a point far away from the symmetry axis where the flow is quiescent is atmospheric, i.e. 1 in non-dimensional coordinates.
In total, the boundary-integral simulation for a single fluid contains three substeps to advance from time step j to j + 1 [27][28][29]. First, with the velocity v (j) we integrate Bernoulli's equation (9) to compute the potential φ (j+1) f along the free surface. The second step updates the position of the free surface by integrating and moves the disc downwards with its prescribed velocity. Finally, we solve the boundary-integral equation (7) to obtain the liquid potential over the disc φ n,f . Then the process repeats.

Boundary-integral method for two fluids
The BIM can be extended to describe two immiscible fluids with a moving interface separating both phases. Our approach closely follows that of [30,31] and is thus only briefly sketched here. The liquid and gas phase both satisfy the boundary-integral equation (7) within their respective domains. At the interface between the two fluids tangential stresses vanish since the fluids are inviscid, i.e., the fluids can slip freely along the boundary. To ensure continuity of the interface, however, the normal velocities of both phases must exactly balance with φn,g being the normal derivative of the gas potential and φn that of the liquid as defined above. The minus sign is due to the normal vector pointing always out of the respective domains. Furthermore, the pressure jump across the interface is given by the Laplace pressure. These two conditions are sufficient to derive the two-fluid version of the BIM [30,31].

Specific implementation of the boundary-integral method
In our axisymmetric situation the surface integrals in Eq. (7) can be reduced to onedimensional line integrals in the (r, z)-plane after analytical integration over the azimuthal angle. Numerical integration is carried out using 8-point Gaussian quadrature with the weak logarithmic singularities removed analytically as in [29]. Between the nodes the interface shape and potentials are interpolated using cubic splines with the node-to-node distance serving as the spline parameter s. Natural boundary conditions (i.e. a vanishing second derivative with respect to s) are used for all splines, except for the potential at the connection between the disc's edge and the free surface as described below. Once the known integrals in Eq. (7) are evaluated we transform Eq. (7) into a matrix equation which is solved by LU decomposition [28,29]. Time-stepping for the integration of Eqs. (9) and (10) is carried out by an iterative Crank-Nicholson procedure during the incompressible stage. In the compressible stage, we use a simpler forward-Euler scheme for ease of coupling between the BIM and the Roe solver in the two respective domains. The time step in the incompressible stage is determined by the condition that neighboring nodes may not collide even if their velocities were directed exactly towards each other. This leads to: with i running over all nodes, v i the free surface velocity at node i, and d i the distance to the neighboring node. This quantity is multiplied with a safety factor f which is in most simulations chosen to be 0.1. In the compressible stage the time-step is determined by the stability condition of the Roe solver as described in the next section.
To ensure a continuous boundary of the liquid domain the last node of the free surface remains fixed at the disc's edge. This connection point requires some special consideration. First, the surface is not smooth and thus the prefactor β in Eq. (7) must be modified to become β = 2α for the liquid and β = 2π − 2α for the gas equations. Here, α is the angle connecting the horizontal disc and the tangent to the free surface at the connection point through the liquid domain (i.e. α > π in our situation). Second, the connection point is at the same time the last node, F, of the free interface and the first node, D, of the disc. To solve Eq. (7) we consider it part of the disc, i.e. we impose φn = −1 and obtain the corresponding value for φ. Together with the pinning at the disc's edge, the position r and the liquid velocity v at D are thus completely determined. We then need to ensure that F and D remain identical. For this we first copy the spatial coordinates of D on F. To ensure further that also the velocity v is identical on both nodes we project v on the tangential and normal vectors of the free surface at F which determines the values of φn and ∂φ/∂s. The latter is imposed as a boundary condition on the spline function for φ. This procedure ensures that the velocity v of the connection point is identical when seen from the disc or from the free surface, even though the respective normal and tangent vectors are discontinuous.
In our simulations the free surface extends from the edge of the disc at r = 1 out to 100 where any motion is negligible and the surface is cut off. Note that BIMs do not require a closed liquid domain since the portion of S at infinity gives no contribution to Eq. (7) provided that φ goes to zero there. We use an adaptive mesh to ensure that the sensitive areas such as the crown splash or the cavity neck are properly resolved without wasting computation time by placing a large amount of nodes on unimportant parts. During the incompressible stage the local node distance d is inversely proportional to the local curvature C with a proportionality constant between 0.05 and 0.005. We impose a maximum distance (dmax = 10) and minimum distance (d min = 0.01) and construct our regridding algorithm such that large gradients in the node density, which might cause instabilities, are avoided. When coupling with the Euler solver during the compressible stage, the BI mesh corresponds to the grid cells used for the Roe scheme as will be described in Section 2.3. Note that then we also allow for node distances smaller than d min .
Boundary-integral methods are known to be vulnerable to instabilities (see e.g. [2]) due to the lack of a naturally damping viscosity which prevents small numerical disturbances from building-up over time. To handle such instabilities we use a smoothing algorithm as follows: At every nth (usually 2 ≤ n ≤ 10) time step the free surface nodes are redistributed such that new nodes fall exactly half-way between old nodes as illustrated in Fig. 2 (a). This periodic regridding procedure efficiently ensures the stability of the numerical scheme [2].
A specific detail of our physical problem is the thin sheet of liquid thrown up around the rim of the impacting disc (see Fig. 1 (a)). In reality, this sheet will quickly develop non-axisymmetric instabilities leading to the formation of individual droplets earning it the title of a "crown splash" as seen, e.g., in the famous pictures of [1]. As here we are not interested in the details of this splash and our axisymmetric code is not able to handle the droplet formation in any case, we cut off the splash as soon the distance between the two sides of the liquid sheet at a given point falls below the local node distance as illustrated in Fig. 2 (b). The actual surface surgery is similar to the one used in [15] to handle the pinch-off of the cavity prior to jet formation.

Roe method for a compressible gas
In the compressible stage the inner gas is described by the one-dimensional Euler equations for conservation of mass, momentum, and energy: All quantities are dimensionless and S is the cross-sectional area of the cavity, ρ the gas density, u the velocity, and p the gas pressure. The total energy E per unit mass is defined as and the total enthalpy H, again per unit mass, is The fact that we use pa as a reference pressure for p leads to the appearance of Eug in Eqs. (13)- (15). Note the two source terms p∂S/∂z and −p∂S/∂t on the right-hand side which account for a cavity radius which is changing in space and time.
Integration of Eqs. (13)-(15) is carried out using a Roe scheme [23,24] whose implementation is fairly standard and thus omitted here.
The time step during the compressible stage is restricted by the Courant-Friedrich-Lewy (CFL) stability condition for the Roe solver. We fix a constant CFL number C * < 1 (usually C * = 0.5) and determine the time step by: with the cell size ∆z, the local speed of sound in each cell c i , and the index i running over all cells.
One detail that is of interest are the boundary conditions at the upper and lower end of the cavity. In contrast to the standard problem of flow through a nozzle which possesses an inlet on one side and an outlet on the other side, in our case gas leaves the shrinking cavity on both sides. Since both outflows are subsonic we can prescribe one flow quantity at each boundary [23]. To compute the flux through the lower (upper) face of the first (last) computational cell we add a boundary cell at each end whose state values are calculated as follows.
At the upper exit we impose that the pressure in the boundary cell be atmospheric. In order to compute the density and velocity of the boundary cell, let ρ N , u N , and p N be the components of the state vector, and c N = γp N /ρ N the local speed of sound in the last computational cell. Assuming that the discharge proceeds isentropically and that the I + characteristic transporting information out of the computational domain is conserved, we have for the values of the upper boundary (ub) cell: p ub = 1.
At the bottom exit we impose the velocity ue with which the gas leaves the computational (Euler) domain. Since velocities between this point and the disc are small (see Section 2.3) conservation of mass allows us to calculate this velocity from the rate of change in cavity volume between the disc and the boundary cell, i.e. ue · π · r 2 e = S b φndS b with the surface S b including the disc itself. Similarly as at the upper exit, the state of the lower boundary (lb) cell can then be computed from the state of the first computational cell (ρ 1 , u 1 , p 1 ) using conservation of the I − characteristic: These boundary conditions cannot completely prevent reflection of waves at the upper and lower ends which causes some small oscillations in the state variables ρ, u, and p as can be seen for example in Fig. 4 (b). The oscillations however remain sufficiently small so that the averaged evolution of the gas dynamics can still reliably be extracted.

Coupling between boundary-integral and Roe method
During the compressible stage of the simulation, the gas domain is split into three separate zones as illustrated in Fig. 3 (a). The gas pressure which is required for coupling to the BIM is obtained in a different way for each zone. First, in the "atmospheric zone" the pressure is taken to be atmospheric since the gas dynamics is negligible. Next, in the "Euler zone" the pressure is provided by the solution of the Euler equations (13)- (15). Finally, in the "bubble zone" the gas dynamics is again neglected. The pressure, however, cannot be taken to be atmospheric since the zone is not open to the atmosphere. Due to the low gas velocities in that zone (of order of the disc speed) the pressure is taken equal to the pressure at the bottom end of the Euler zone. The upper end of the Euler zone can remain at a fixed vertical position (between 0 and 1 disc radii below the initial free surface). The bottom end is fixed to be at the maximum radial extension of the cavity below the neck as depicted in Fig. 3 (a). Note that the Euler zone cannot be extended all the way down to the disc since the free surface departs almost horizontally from the disc's edge. This leads to large gradients in the cavity radius on the right-hand side of Eq. (14) and can thus cause numerical instabilities. Since the bottom end of the Euler zone is moving downward in time, its length needs to be extended during the simulation. Our algorithm adds a new Roe cell whenever the distance between the desired location (maximum radial extension of the cavity) and the actual position of the last Roe cell becomes larger than the size of a single cell.
For the results presented here, the initial number of cells in the Roe solver is 600 and grows dynamically by extension of the Euler zone. To ensure sufficient resolution for the high air speeds even during the last moments of the simulation, we double the number of cells by splitting each cell in half whenever the minimum cavity radius r neck falls below a certain value. This is done twice: at r neck = 0.2 and r neck = 0.05. Due to extension and node doubling the number of Roe cells at the end of the simulation is  Fig. 3 (a) Illustration of the three zones in which the gas domain is split during the compressible stage: in the "atmospheric zone" on top, the pressure is always atmospheric and the gas dynamics are neglected. In the "Euler zone" we use the Roe scheme to calculate the fully compressible gas dynamics. In the "bubble zone" again the pressure is constant and given by the pressure at the end of the Euler zone. (b) Schematic illustration of the alignment of Euler cells and BI nodes: the BI nodes are always placed at the height of the center of the Euler cells. In the neck area a BI node is placed in every 2nd Roe cell, while away from the neck the spacing is larger to save computation time.
somewhat above 3000. The height of the Roe cells is always constant with ∆z = 0.0012. We find that the number of Roe cells is not crucial for the total computation time which is mainly determined by the BIM calculations. It is crucial to properly align the positions of the BI nodes with the Euler cells. For this we place a BI node always exactly in the center of every nth Euler cell as illustrated in Fig. 3 (b). Usually n = 2 in a refined zone around the cavity neck and n=5 outside this zone. The cross-sectional area S as well as its spatial and temporal derivatives for each cell are calculated at the height of the cell center from the splines interpolating the cavity surface in the BIM. To avoid numerical instabilities of the BIM we use the periodic regridding described above which now makes the BI nodes "jump" between Euler cells. This implies that there must always be at least one Roe cell without a BI node between two cells which contain a node, i.e. n ≥ 2.
The two-way coupling between the gas and the liquid domain is accomplished as follows. At each time step j we first do a BI step to advance the shape and velocity potential of the free surface from j to j + 1 using the gas pressure of step j. The BI step further provides the derivatives ∂S/∂z = 2πrdr/dz and ∂S/∂t = 2πrṙ. This calculation is followed by a Roe step using the new cavity shape j + 1 to obtain the new pressure at j + 1 and so forth.
We performed an extensive set of simulations to verify that the results presented in the next section are numerically robust when changing any of the above mentioned simulation parameters.

Justification of the 1D compressible scheme
In Fig. 4 (a) we show the gas velocity through the cavity neck as a function of the shrinking cavity neck r neck for the incompressible two-fluid BIM. We use r neck instead of time to allow for an easier comparison with experiments in Section 3.3 and [18]. Already at r neck = 0.05 (corresponding to 1 mm) an incompressible gas would surpass  Fig. 4 (a) The gas velocity at the neck as obtained from an incompressible two-fluid BI simulation. The velocity easily surpasses the speed of sound demonstrating the need for our more sophisticated compressible approach. (b) The velocity obtained from the multiscale simulation (red line) agrees so well with the velocity from the two-fluid BIM (blue line) that both are hardly distinguishable for low velocities where compressibility is negligible. The start of the red curve marks the transition from the incompressible to the compressible stage. The small oscillations in the red curve are due to wave reflections at the ends of the Euler zone. (c) The gas flow field obtained from the two-fluid BIM at the moment when the velocity at the neck reverses and the simulation passes from the incompressible to compressible stage. Except for a rather narrow zone around the neck the flow is to a good approximation one-dimensional.
the speed of sound. This clearly demonstrates the need for a CFD method which takes the fully compressible gas dynamics into account if one wants to study the ejected air stream close to cavity collapse.
A full two-way coupling is required as both gas and liquid flows occur on similar time scales: We first estimate the typical time scale for the gas flow Tgas as the length of the cavity L ≈ 10 cm divided by the speed with which a perturbation travels, c = 330 m/s, to obtain Tgas = 0.3 ms. A typical time scale for the variation of the cavity radius can be derived by considering the time that it takes the cavity to collapse from a neck radius of 4 mm 2 down to zero which is Tcav ≈ 1 ms and thus of the same order as Tgas.
For some time after switching from the incompressible to the compressible stage the gas velocities are still moderate and compressibility effects should be negligible. We can thus expect that in the beginning of the compressible stage the gas velocity obtained from our multiscale approach should be similar to the two-fluid BIM in Fig. 4 (a). That this is indeed the case is demonstrated in Fig. 4 (b) which gives us a first indication that the coupling between the Roe solver and the BIM works correctly. It further gives good evidence that the assumption of one-dimensional gas flow in the compressible stage is justified.
To verify the 1D assumption more explicitly, Fig. 4 (c) shows the flow field obtained from the two-fluid BIM at the moment of flow reversal. Except for a small region around the neck where by definition the vertical velocity is zero, our assumption is well justified.

Structure of the compressible gas flow
The intricate structure of the gas flow in the Euler zone is illustrated by the velocity profile in Fig. 5 for various instants of time. Here we normalize velocities with the speed of sound to obtain the Mach number Ma=u/c. At early times as in Fig. 5 (a) one appreciates that at the bottom end of the Euler zone air is entrained by the downward moving disc at velocities of the order of the disc speed. This downflux is however overcompensated by the shrinking of the cavity around the neck so that the total flux is directed upwards as can be seen by the velocity maximum at z ≈ −1.9. Above the maximum the cavity widens again and the velocity decays towards the upper end of the cavity. The consequence of the competition between cavity expansion at the bottom and cavity collapse around the neck is the creation of a stagnation point with Ma = 0 as is clearly visible in Fig. 5 (a).
As the collapse progresses gas is pushed through the rapidly diminishing cavity neck at ever higher and higher speeds leading to a sharp velocity peak at the neck as illustrated in Fig. 5 (b). At the top and bottom boundaries of the Euler zone the velocities remain almost unaltered as compared to Fig. 5 (a).
Finally, as the flow speed increases even further a shock wave develops upstream of the neck as shown by the magnification in Fig. 6. Thanks to the shock-capturing ability of the employed Roe scheme our method is able to handle the shock formation quite well. Figure 7 shows the Mach number at the cavity neck. For our standard configuration of a 2 cm disc impacting at 1 m/s the flow becomes sonic at a neck radius of 0.025 (corresponding to 0.5 mm). Here we also add data for a higher impact speed of 2 m/s which qualitatively shows the same behavior but where sonic flow is attained already at a neck radius of 0.06 (1.2 mm). Once sonic velocities at the neck are reached our numerical scheme becomes unstable which is why do not present any data beyond this point. Note that due to our unsteady situation supersonic speeds are attained even earlier at locations above the neck, cf. Fig. 6.
We now turn to study the pressure distribution in the Euler zone which is illustrated in Fig. 8. In the early stages of the process (Fig. 8 (a)) the pressure remains virtually atmospheric with only a very slight dip around the cavity neck. At a later time, however, the pressure at the neck diminishes substantially as can be seen in Fig. 8 (b). In a steady state situation one would expect the neck pressure to reach a minimum value of with γ = 1.4 the isentropic exponent, as Ma neck becomes unity. As shown in Fig. 9 (a) our situation -although highly unsteady -exhibits a similar behavior with p neck ≈ 0.6 at the final instant before the simulation destabilizes. Notably, below the neck the pressure is very uniform all the way down to the end of the Euler zone ( Fig. 9 (a)) which allows us to define a single pressure value for the "bubble" between the neck and the disc. Figure 9 (b) demonstrates that in our unique situation of air being pushed through a rapidly shrinking nozzle sonic speeds can be attained with a bubble pressure which is merely 2% higher than the surrounding atmosphere.   Fig. 7 The Mach number at the neck for two different impact speeds: 1 m/s (red solid curve) and 2 m/s (green dashed curve). The flow becomes sonic when neck has shrunk to 0.025 or 0.06 disc radii, respectively.  Fig. 8 The pressure profile in the Euler zone at the same instants as in Fig. 5. The low pressure at the neck is caused by the high gas speeds in that region.

Experimental validation
Here we describe briefly the validation of our numerical results with the experimental data of [18] which is achieved in three different ways: (i) We use smoke particles to directly measure the air speed as it is pushed out of the collapsing cavity. We can reliably measure air speeds up to 10 m/s and very good agreement with the numerical data is found as presented in Fig. 2 of [18].
(ii) In order to confirm the validity of our predictions also for higher gas velocities we compare the numerical cavity shape close to pinch-off with experimental images. We find that the experimental shape is not smoothly curved but exhibits a rather pronounced "kink" at its neck. This effect -which is due to the low pressure induced by the high gas speeds -is not present in single-fluid simulations, but can be reproduced rather well by the inclusion of air effects as shown in Fig. 3 of [18].
(iii) Finally, we show that the cavity neck in the experiment exhibits a significant upward motion prior to final collapse as the fast air stream pushes the surface minimum upwards. We find very good agreement for this motion between experiment and our multiscale simulations as shown in Fig. 4 of [18].
The quantitatively consistent observation of the above air effects in our compressible simulations and corresponding experiments gives us strong confidence in the reliability of our numerical scheme.

Conclusions
We presented a multiscale model to simulate the impact of a solid object onto a liquid surface. Our focus was on the fast stream of air that is pushed upwards as the impact cavity collapses due to hydrostatic pressure. We showed that in our case of a 2 cm disc impacting at 1 m/s the air attains supersonic velocities and thus requires the use of a fully compressible computational method -in contrast to existing treatments such as Volume-of-Fluid or Level-Set methods [4][5][6][7]17] in which the air flow was considered incompressible.
In our simulations the impact process is split in an incompressible and a compressible stage. During the incompressible stage which covers the first part of the process, air is entrained into the cavity at relatively low speeds (compared to the speed of sound). This allows us to use a two-fluid boundary-integral method for the gas and the liquid domain. The compressible stage starts as the air flow reverses and air is pushed out through the cavity neck. In this stage we couple two different methods: a Roe scheme to solve the one-dimensional Euler equations in the gas domain and a singlefluid boundary-integral method in the liquid domain. The two domains are connected via the pressure at the free surface.
The predominantly one-dimensional character of the air stream and the shape of the impact cavity make our system reminiscent to the common problem of air flow through a converging-diverging nozzle in aerodynamics. There is, however, an important and fundamental difference: since in our case the confining cavity is a liquid, the "nozzle" shape is rapidly evolving in time.
For low gas velocities we find very good agreement between our simulations and direct experimental measurements [18]. The shape of the cavity as well as a final upward motion of the cavity neck [18] give further strong evidence that our multiscale numerical method faithfully reflects reality.