Towards an Acoustic Simulation of a Water Drop Impacting in a Water Pool

The sound which is produced when a water drop impacts into a water pool is a prominent example for acoustics produced by multiphase flow. In this work the feasibility of numerical methods for simulating this challenging test case is evaluated. First the multiphase flow needs to produce the correct physical mechanisms, e.g. the bubble entrapment. For this an in-house block-structured finite-volume solver with the volume-of-fluid method is used. For the curvature computation a standard finite difference method within the continuum surface force model is employed, including some necessary improvements. A high resolution in space and time is essential and therefore the method is parallelized by domain decomposition. The acoustic part is simulated with the linearized Euler equations which are valid in each phase but need to be adapted in the interface region. The results are compared with numerical and experimental data. It is shown, that the methods are suitable for simple test cases. A coupled drop impact test case corresponds with equivalent experiments until the drop detachment. The acoustic pressure shows a significant rise in the vicinity of the bubble detachment within both phases. However, an oscillation of the cavity bottom can not be observed in the multiphase neither in the acoustic outputs of the airborne signal.


Introduction
Since the beginning of the past century people are interested in the sound producing mechanisms caused by a drop impacting into a pool (Mallock 1918). Over the past decades various researches studied the flow regimes during the drop impact experimentally or numerically (Oguz and Prosperetti 1990;Morton et al. 2000;Zhang et al. 2012). Pumphrey and Walton (1988) concluded by their studies, that the "regular bubble entrainment" only occurs under certain conditions, which were later characterized by the Weber number 1 3 and the Froude number where g is the acceleration due to gravity, U I is the impact velocity, D is the drop diameter and is the surface tension. Oguz and Prosperetti (1990) found that the boundaries of the regular entrainment region, shown as shaded area in Fig. 1, could be approximated by the following power law relation: in which the upper boundary is given by A = 48.3 , = 0.247 and the lower boundary by A = 41.3 , = 0.179 . In this region a bubble detaches from the bottom of the cavity.
The formation of a bubble in the liquid was found to coincide with the start of a sound pulse (Franz 1959). Although no experimental work has proven this theory, it was assumed by Leighton (2012) that the airborne sound results from the underwater sound field propagating through the water-air interface. However, progress was recently made by Phillips et al. (2018) in determination of the source of the characteristic so-called "pling" sound. An oscillation of the bottom of the cavity, induced by the detached air bubble, was identified as the source.
To reproduce the sound of a drop impact, simulations of multiphase flows and acoustics are necessary. Finite-volume solvers for multiphase flows can be divided into interface-tracking and interface-capturing methods (Cano-Lozano et al. 2015). Front Tracking techniques, level-set methods, marker and cell (MAC) techniques or the volume-of-fluid (VOF) method capture the position of the interface sharp or diffuse. Its mass conservation and the good handling of topology changes made the VOF method to one of the most popular methods.
(1) We = U 2 I D (2) Fr = U 2 I gD , (3) We = A ⋅ Fr Fig. 1 Froude-Weber diagram with shaded regular entrainment region based on Oguz and Prosperetti (1990); square represents the test case by Phillips et al. (2018); circle represents the test case by Morton et al. (2000) 0 In comparison to multiphase flows the research of computational aeroacoustics (CAA) is in an early stage. Mainly, this is based on the lack of computational power when combining acoustics and flow. The small sizes of the acoustic variables or the high frequencies to be resolved are two examples of difficulties which are less critical in other CFD applications (Hardin and Pope 1994). When using direct noise computations (DNC), the compressible Navier-Stokes equations (NSE) are solved for the aerodynamic and acoustic field simultaneously. This approach is inefficient in low Mach number flows due to high spatial and temporal resolution requirements. A more suitable approach is the acoustic/ viscous splitting technique by Hardin and Pope (1994). Density, velocity and pressure are considered as a composition of an incompressible background flow with a superimposed acoustic perturbation, which is understood as Expansion about Incompressible flow (Shen and Sorensen 1999).
Various methods solving the acoustic and multiphase part are existent, but few works handle a coupled simulation. For example, the sound generation and propagation in two phases were simulated by Tajiri et al. (2010) with the help of a finite difference lattice Boltzmann method.
In this paper we are simulating the acoustics of a water drop impacting into a water pool by using a finite-volume framework solving the incompressible NSE extended by the VOF method for the multiphase flow and the Linearized Euler Equations (LEE) with a highresolution scheme for the acoustics. The outcome of our methods will be compared with the work of Morton et al. (2000) and Phillips et al. (2018).

Numerical Methods
The following methods are implemented in our in-house flow solver FASTEST, which uses the finite-volume method on block-structured meshes. Spatial parallelization is carried out with MPI (Message Passing Interface). At each block boundary one additional layer of control volumes (CVs) is embedded. By the use of this layer, data can be exchanged throughout a simulation. To avoid a loss of efficiency, the number of CVs per processor are distributed equally. The drawback of this parallelization strategy is that numerical methods which need more CVs than number of layers available have to be adapted at block boundaries.
The multiphase flow is modeled using the incompressible NSE with the one-fluid formulation (Prosperetti and Tryggvason 2009) as with the density , the viscosity , the velocity vector u i , the pressure p, the time t and the Cartesian coordinates x i . A main difference to the single phase set of equations is the last term in Eq. (5) that describes the surface tension.

Multiphase Flow
To capture the incompressible and immiscible multiphase flow the VOF methodology is employed. The scalar represents the local volume fraction in a CV, so that Every CV which contains a value of 0 < < 1 is defined as an interface CV. Within the NSE the volume fraction distinguishes the material properties ( = { , }) as a result of A transport equation for advecting the volume fraction is added into the framework: Employing the Crank-Nicolson method for the temporal discretization (Ubbink and Issa 1999) results in the equation where index n denotes the different time levels, V P is the CV volume, Δt is the time step size, n b is the number of neighbors and S i,f is the face area. The assumption that the variation of the velocity field does not influence the volume fraction transport during one time step is reasonable, when the time step size is small enough (Ubbink and Issa 1999). Subsequently u n+1 i,f = u n i,f and Eq. (9) can be reduced and rearranged for the unknown new value of the volume fraction as To guarantee a bounded solution while maintaining the sharpness of the interface, highresolution (HR) schemes are necessary to calculate the face values of the volume fractions f . Here the M-CICSAM (Waclawczyk 2007) consisting of the compressive Hyper-C-Scheme (Leonard 1991) and the high-order FROMM-Scheme (Fromm 1968) is employed. The blending of these two schemes depends on the angle between the normal vector of the interface and the line connecting the donor CV center with the acceptor CV center. Comparisons between different HR-schemes can be found in Ubbink and Issa (1999) and Waclawczyk and Koronowicz (2008).

Surface Tension
A constant surface tension , the curvature , the unit normal vector to the multiphase interface n i and the absolute value of the smoothed volume fraction gradient |∇̃| assemble the (6) (x i , t) = 0 fluid 1, 1 fluid 2.
surface tension term. Brackbill et al. (1992) introduced the continuum surface method (CSF) which rewrites the surface tension force into an equivalent volume force. The volume fraction gradient becomes zero outside the diffuse interface, since surface tension only acts on the interface. To obtain the curvature and the normal vector at each interface position P the first and second spatial derivative of the volume fraction are needed: The gradients, exemplary shown for two dimensions in Eq. (11) are computed with a finite differencing scheme in the direction of the indices. To improve the accuracy of the curvature values in each cell, some additional features are applied. A better estimation of the first derivatives is achieved by filtering the volume fraction with a smooth least-squares polynomial (LSP) filter of degree two (Palacio et al. 2002): The mean value of all spatial directions x 1 = x , x 2 = y and x 3 = z , derived by the sum of the coefficients b k times the corresponding concentrations divided by the sum of the coefficients, leads to a smoother distribution of the concentration. Using the central differencing scheme on the first derivatives to obtain the second derivatives, once again applying a filter on the first derivatives does not improve the outcome. Missing values at block boundaries are extrapolated with known concentrations. For further improvement of the curvature a 5-point Chebyshev filter (Palacio et al. 2002) in the manner of Eq. (12) and a volume fraction dependent weighting (Denner and van Wachem 2014) are applied on the curvature itself. The coefficients for the different filters are listed in Table 1.
As stated by Renardy and Renardy (2002) the curvature is expected to be more accurate when the volume fraction is not close to zero or unity. Therefore, the weighting is defined as with index Q denoting all the direct neighboring CVs and the weighting factor with the exponent as In Eq. (14) it can be seen, that the factor reaches zero if the volume fraction tends zero or unity and becomes maximal when the volume fraction equals 0.5. In contrast to Denner (2013) we choose = 1 , because for dynamic problems an exponent greater one lead to excessive damping.

Acoustics
The generation and propagation of the acoustic waves are computed with the Linearized Euler Equations within an acoustic/viscous splitting technique (Hardin and Pope 1994). Density, velocity and pressure are split into an incompressible part and an acoustic perturbation: Incompressible flow and acoustic quantities are marked with inc and ac , respectively. Replacing the variables in the compressible NSE with this ansatz and assuming that the acoustic quantities are small in comparison to the flow quantities, that the flow within each phase is incompressible and that the higher order acoustic terms are neglected lead to: Equations (18) and (19) represent the evolution of the acoustic density and velocity.
Assuming an isentropic flow (Entropy s = const. ), a general expression for the speed of sound is based on the equation of state p = p( , s) and derived by the partial derivative with respect to the density: The speed of sound or in this context the velocity of the acoustic waves is only dependent on the "stiffness" of the medium, e.g. the ability to respond with compression in present of increasing pressure (Leveque 2002). Finally the acoustic pressure (21) can be derived by combining the general equation for the speed of sound (20), the equation for the acoustic density (18) and the time derivative of the pressure decomposition (17): Further details about the LEE and its derivation can be found in Ewert and Schröder (2003), Hardin and Pope (1994), Kornhaas et al. (2015) and Shen and Sorensen (1999).
Due to the small energy scales in the acoustics, a one-way coupling is employed. Hence, only the flow influences the acoustic computations with a source term, which is represented by the temporal derivative of the pressure on the right hand side of (21). Applying this technique to multiphase systems with surface tension an adaption has to be made. When the pressure changes over the interface and additionally the interface moves, the source term locally becomes non-zero from one time step to the next. The resulting acoustic sources are unphysical and therefore have to be suppressed in the interface region (Friedrich and Schäfer 2018) according to: A second mechanism is employed, in which the acoustic source term computation is corrected with a reference pressure, see Eq. (23). Numerically caused pressure changes in flow regions containing air influence the values in the whole domain. When an acoustic pulse travels from air into water, its amplitude gets increased twofold. Depending on if a CV is located within the air or the water it is corrected with a corresponding reference value, which is taken from a point away from the area of interest: For determination of the acoustics the LEE are transformed into a local coordinate system with the coordinate which is normal to the CV face. A one-dimensional problem is solved at each CV face: In (24) the Jacobi matrix A is defined as and the variable vector U as In multiphase flows the different densities and speeds of sound define the locally changing impedance which is referred as layered media (Leveque 2002). A high-resolution scheme consisting of Godunov's method and the Lax-Wendroff method combined with an Osher limiter ( = 1 ) is applied for the acoustic wave propagation (24) in this layered media. At block boundaries the high-resolution scheme reduces to Godunov's method caused by one ghost cell layer.
As stated in Leveque (2002) the methods used for the acoustic computation are typically stable only with an acoustic CFL number less than 0.5. Due the high values for the speed of sound, this would lead to an unnecessary small time step for the coupled CV in air, simulation. To avoid this problem, N acoustic time steps are performed during one flow time step: The acoustic CFL number in Eq. (27), in which Δt is the time step size and Δh the CV size, is set in advance of a simulation and then N acoustic time steps are performed accordingly after each converged flow iteration. This subcycling technique assumes that the acoustic sources are constant within the corresponding acoustic time steps (Kornhaas et al. 2015).

Preliminary Test Case
In previous work, it was shown that the acoustic transmission and reflection at a stationary multiphase interface is in good agreement with analytical values (Staab et al. 2015). In order to check if the methods described in the previous sections are suitable for the drop impact test case, a preliminary three dimensional test case is set up. According to Phillips et al. (2018) the main driver of the drop impact acoustics is the bubble oscillation after the detachment. Therefore, an ellipse-shaped air bubble, defined by the major axis of 0.923 × 10 −3 m and the minor axis of 0.675 × 10 −3 m , is initialized under a water surface at (0.0025, 0.0015, 0.0025) with a gap of 0.139 × 10 −3 m . The grid consists of 128 × 128 × 128 CVs and the material parameters for air and water are listed in Table 2. The speed of sound is scaled due to the comparable small domain of 5 × 10 −3 m × 5 × 10 −3 m × 5 × 10 −3 m . Surface tension will force the ellipse to return in a circle shape which leads to an oscillation of the bubble and consequently to the surface as it is shown in Fig. 2.
To validate the methods from Sects. 2.1, 2.2 and 2.3, the acoustic pressure above the surface and the surface motion are compared with respect to their frequency and phase (27) CFL ac = c Δt ac Δh , with Δt ac = Δt∕N. shift. The airborne signal taken at a monitoring point (0.003, 0.004, 0.003) and the surface movement right above the bubble are shown in Fig. 3.
First of all, it should be mentioned that the acoustic signal has an oscillatory course. From this it can be concluded that the manipulation of the source term does not suppress the expected behavior. Not taking into account the first peak of each signal, an average frequency of the following four peaks for the airborne regime can be found at 664 Hz, which is less than 2% diverging from the average surface oscillation frequency with 677 Hz. In addition to the frequency, the shift in time of the two signals is used for validation. Between the first upward directed peak of the surface at 1.38 × 10 −3 s and the first positive peak of the acoustic pressure at 2.52 × 10 −3 s , there is a difference of Δt 1 signals = 1.14 × 10 −3 s . Between the first downward directed peak of the surface at 2.13 × 10 −3 s and the first negative peak of the acoustic pressure at 3.46 × 10 −3 s , the difference in time is Δt 2 signals = 1.33 × 10 −3 s . The distance of the first CV above the surface, in which the acoustic source is not suppressed, to the monitoring point is Δd = 1.52 × 10 −3 m . With a speed of sound of 1 m∕s , the phase shift between the surface movement and the acoustic signal is in plausible range so that Δd ≈ c ⋅ Δt 1,2 signals . The preliminary test case provides a promising result and therefore the next step is to apply the methods to the more complex drop impact test case.

Remarks on Resolution
The outcome of the work by Phillips et al. (2018) anticipates some settings and boundary conditions. Others are found by the need of the correct flow behavior. The task is to find an optimal balance between computational time and physical correctness. Morton et al. (2000) stated in his work, that a certain resolution is necessary to capture all the effects from the experiments. To achieve a sufficient resolution within an acceptable amount of computational time, the grid gets adapted in an r-adaptive manner. First the simulation is performed on a uniform grid to identify the area of interest. Afterwards the CVs are concentrated in this area while keeping the overall number of CVs constant. Inside the refined area the grid has a Cartesian structure. To investigate the influence of the grid, the cavity depth over time is computed on three different grid levels using a 2.9 × 10 −3 m diameter drop impacting with U I = 1.55 m∕s ( We = 84 , Fr = 94 , circle in Fig. 1). An example of the coarsest mesh can be seen in Fig. 4 and the grid level details can be found in Table 3. The results of the mesh study in Fig. 5 are scaled in length by the initial drop diameter D and in time by D∕U I .

Spatial Resolution
No bubble entrapment occurs on the coarse grid, in comparison to the higher grid levels. Since the detachment of a bubble is crucial for the acoustic part, the coarse grid will not be included for the coupled simulation. However, the path of the three grid levels are very similar until = 4 where the coarse grid regime stays at a lower cavity depth. The medium and fine resolutions diverge in the interval between = 6 and = 8 . During this period the bottom of the fine grid resolution deflects slightly to the surface.  Morton et al. (2000) and our simulations on three grid levels  Our results of the medium and fine resolutions are compared with the numerical and the corresponding experimental results of Morton et al. (2000). In Fig. 5 the simulation and experimental data by Morton show different courses over time. First there is an offset, while from = 4 to = 6 both lie on top of each other. After that point the simulation data of Morton stays constant until the end, where the cavity depth gets lower. The experimental data in contrast shows some sort of wave from = 6 on with a rise of depth at the end. It can be seen that the results from our simulation on the medium and the fine grid level also behave wave like from = 6 on with a rise at the end. At the beginning our simulations show the same offset as the simulation data of Morton et al. However, from = 5 to = 8 our fine grid level regime agrees within a 2% range to the experimental data of Morton. The results in Sect. 4 are obtained with the fine grid setup.
Finally, the importance of the curvature computation on the fine grid is exposed in Fig. 6. The drop loses its shape during the impact into the surface and some fluid emerges at the top of the drop. With the described improvements to the standard model from Sect. 2.2, the surface tension stays stable as the shape of the drop shown in Fig. 7.

Temporal Resolution
The maximum temporal resolution is either given by the expected acoustic signal or the capillary time-step constraint of the multiphase flow.
From theory the acoustics are produced by an oscillating bubble which introduces oscillations in the bottom of the cavity. The oscillation frequency was found to be in range of the natural oscillation frequency f for a submerged gas bubble derived by (Minnaert 1933): With bubble radius r, specific heat ratio for the gas inside the bubble , static pressure on the exterior surface of the bubble P 0 and the density of the surrounding liquid the frequency is a function of the bubble radius. As well in the experiments as in the simulations the radius can vary. A diameter range of 0.7 × 10 −3 m to 0.9 × 10 −3 m leads to a frequency range of approximately 9500-7250 Hz. To ensure that the multiphase simulation is able to produce this oscillation, a maximum time step size is determined following the Nyquist-Shannon sampling theorem (Shannon 1948), which states that the sampling frequency should be twice as high as the signal frequency. In order to not only ensure a distinguish association of a frequency but rather a clear resolution, the sampling rate is being increased tenfold. These considerations concerning the acoustics lead to a maximum time step size of Δt max = 1 × 10 −6 s. Within multiphase flows with surface tension, the time step size needs to resolve the propagation of capillary waves for a stable simulation (Brackbill et al. 1992). In Denner and van Wachem (2015) it is demonstrated, that the constraint needs to be applied irrespective of an explicit or implicit surface tension implementation. According to Brackbill et al. (1992) the capillary time step constraint leads to a maximum time step size of With the material properties of air and water as listed in Table 2 and the minimum cell size of the fine grid level in Table 3, a maximum time step size of Δt max = 1.3986 × 10 −5 s results.
Comparing both criteria the chosen maximum time step size follows the acoustic considerations. If the time step size is in conflict with the CFL number, which is limited at 0.3, an adaptive time step mechanism reduces the time step size accordingly. Phillips et al. (2018) deliver snapshots of a high-speed video at key stages of a 4.0 mm drop impacting with U I = 1.29 m∕s ( We = 90 , Fr = 42 , square in Fig. 1). In Figs. 8 and 10 the snapshots of the experiment are compared with the equivalent results from our simulation in Figs. 9 and 11.

Results and Discussion
The simulation is carried out with the fine grid resolution from Sect. 3.1, a maximum time step size of Δt = 1 × 10 −6 s and the material properties of Table 2  velocity the moments are shifted in time by 3.3 × 10 −3 s in the first three frames in Fig. 9 and by 4.7 × 10 −3 s in the last frame of Fig. 9 and in every frame in Fig. 10. Every basic feature stated by Franz (1959) during a drop impact can be seen in the simulation: cavity creation ( Fig. 9 at 7.23 × 10 −3 s and 12.6 × 10 −3 s ), begin of recoiling due to surface tension ( Fig. 9 at 17.56 × 10 −3 s and 18.38 × 10 −3 s ), entrapment of a small air bubble ( Fig. 11 at 19.51 × 10 −3 s , 19.84 × 10 −3 s and 20.31 × 10 −3 s).
However, the main key feature of the acoustic source, which is the oscillation of the cavity bottom can not be observed explicitly in the simulation. For a better understanding of the separation process in the simulation, more frames between 19.84 × 10 −3 s and 20.31 × 10 −3 s are shown in Fig. 12. In the first row, the bubble detaches from the cavity bottom while in the second and third row the bubble oscillation starts and the bottom moves upwards. Although the bubble does some oscillatory movement from the last frame in the first row to the last frame of the third row, the cavity bottom seems to move monotonous. There are different possible reasons for the lack of visible oscillation, e.g. insufficient mesh resolution, insufficient output resolution or unfavorable choice of drop parameters.
The multiphase flow does not show the desired detail in the physical behavior, nevertheless the acoustic part will be analyzed. The acoustic pressure path from the experiments of Phillips et al. (2018) is close to zero until bubble detachment and from that point on a damped oscillation of goes through the domain. The values of the amplitudes of the acoustic pressure differ between the fluids about 10 1 to 10 2 with smaller values in air.
In Figs. 13 and 15 the acoustic pressure in air and water at monitoring points (0.0173, 0.0208, 0.0173) and (0.0173, 0.0065, 0.0173) are shown, respectively. From the scales of the figures, it can be seen, that the amplitude difference between the phases is in the  expected range. At the beginning a deflection in the underwater signal occurs, which is caused by the initialization of the simulation and not by the impact of the drop. Other deflections until the drop detachment represent small changes of the time step size, for example in the underwater signal at 4.3 × 10 −3 s or 17.8 × 10 −3 s. To ensure that time step size changes do not influence the acoustic computation during the bubble detachment, it is held constant from 17.78 × 10 −3 s to 22.72 × 10 −3 s.
Within both signals a significant rise can be associated to the time right before the first frame of the first row in Fig. 12. Between the first and the second frame in the first row, which represent the pinch off, the maximum values are reached in the acoustic signals. That the pinch off collides with the maximums of the acoustic signals corresponds to the previous mentioned experimental work. In the airborne signal a wave pattern with a small amplitude occurs from 20.25 × 10 −3 s to 20.6 × 10 −3 s right after the pinch off wave ( 19.9 × 10 −3 s to 20.25 × 10 −3 s ), see Fig. 14. Since there is no visible oscillation of the cavity surface, the small wave pattern can not be associated with absolute certainty to the motion of the cavity bottom. In contrast to this the underwater signal shows an indefinable deflection. A close-up of the acoustic pressure in water during the detachment process is displayed in Fig. 16.
From the theory of the experiments there should be an oscillation in the cavity bottom, thus in the airborne signal. Both phenomena can not be observed within the simulations in the expected manner. To improve our simulations an even finer resolution, different surface tension models or different approaches for the acoustic source computation are part of the ongoing research.

Conclusion
We presented a finite-volume framework for simulating a coupled multiphase acoustic problem. The methods have been described and a preliminary test case achieved promising results. For the complex drop impact test case our simulation did not show the same results as the experimental work of Phillips et al. (2018). Although our simulation reproduced the main features of a drop impact, e.g. building a cavity due to the impact, recoil of the cavity due to surface tension and an entrapped air bubble, the crucial detail was missing. An oscillation of the bottom of the cavity as the main driver of the airborne acoustics was not present in our simulation. However, a significant rise of the acoustic signal in both phases with the start of the bubble detachment process was present, which corresponds to experimental results.
The next step is to simulate this test case on an even higher grid level to resolve the detachment process more precisely. Further investigations on how the surface tension model influences the detachment process as well as uncertainties in the bubble shape or the dimensionless quantities could lead to better results.