**Numerical Modelling of Air Pollutant Dispersion in Complex Urban Areas: Investigation of City Parts from Downtowns Hanover and Frankfurt**

#### **Mohamed Salah Idrissi 1, Nabil Ben Salah <sup>2</sup> and Mouldi Chrigui 1,3,\***


Received: 7 May 2019; Accepted: 16 July 2019; Published: 18 July 2019

**Abstract:** Hazardous gas dispersion within a complex urban environment in 1:1 scaled geometry of German cities, Hanover and Frankfurt, is predicted using an advanced turbulence model. The investigation involves a large group of real buildings with a high level of details. For this purpose, Computer Aided Design (CAD) of two configurations are cleaned, then fine grids meshed in. Weather conditions are introduced using power law velocity profiles at inlets boundary. The investigation focused on the effects of release locations and material properties of the contaminants (e.g., densities) on the convection/diffusion of pollutants within complex urban area. Two geometries demonstrating different topologies and boundaries conditions are investigated. Pollutants are introduced into the computational domain through chimney and/or pipe leakages in various locations. Simulations are carried out using Large Eddy Simulation (LES) turbulence model and species transport for the pollutants. The weather conditions are accounted for using a logarithmic velocity profile at inlets. CH4 and CO2 distributions, as well as turbulence quantities and velocity profiles, show important influences on the dispersion behavior of the hazardous gas.

**Keywords:** buildings; urban area; pollution dispersion; turbulence; Large Eddy Simulation (LES)

#### **1. Introduction**

A growing concern about air pollution in urban environments results in stringent legislation for devices certification and energy utilization [1]. Several studies have been interested in the pollutant emissions reduction and the optimization of combustion systems by the use of new burners [2,3] or by capturing target species [4]. To assess the problem in detail, accurate prediction of the transport and dispersion of airborne contaminants in urban environments is needed. Several methods have been used to investigate air pollution [5,6]. Wind tunnels only provide an understanding of the dynamics and scaled data. Such approaches are able to perform parameters sampling points and to deliver well-documented data sets for the evaluation of numerical models [7]. Yet, due to scale limitations, it is a challenge for such models to fully be tested in large scale atmospheric turbulent flow conditions [8]. Numerical approaches have numerous advantages, e.g., high controllability, lower cost, three-dimensional (3D) data output, etc. [9].

Many preceding studies showed that Computational Fluid Dynamics (CFD) simulations, based on the Reynolds averaged Navier-Stokes (RANS) equations, delivered unaccrued results in reproducing the wind-flow pollutant dispersion around buildings [10–14]. Britter and Hanna [11] pointed out that the RANS turbulence model can be used in simple geometry (isolated building) to produce reasonable qualitative results for mean flows. When compared with laboratory or field experiments, this model showed a better performance than that of simple operational models. Blocken et al. [13] tested three cases of tracer gas dispersion in a neutrally stable atmospheric boundary layer and compared their results against wind tunnel full-scale measurements. For all three cases, they concluded that RANS simulations significantly underestimate the cross-wind (lateral) dispersion of the pollutants. Tominaga and Stathopoulos [14] numerically studied the dispersion around an isolated cube. They found that RANS models underpredict horizontal concentrations at the roof of simulated buildings.

To reach more accurate results, several researchers have turned to the Large-Eddy Simulation (LES) approach, which demonstrated its capability and potential in microscale atmospheric dispersion modelling [15–25]. Liu and Barth [15] studied the flow in 3D streets using scalar transport. Their results showed the limitations of eddies developments with sizes larger than 0.5H (H height of the building) in the so-called free surface layer and produce less turbulent kinetic energy (TKE) above the rooftop. Sada and Sato [18] simulated the stack-gas diffusion around a cubical building in order to predict the instantaneous concentration fluctuation of a plume. Their results were compared with data from wind tunnel experiments and showed a good agreement. Liu et al. [25] used the one-equation subgrid-scale (SGS) parameterization and the LES model to study the flow and pollutant transport in 2D urban street canyons. Their results showed that the combined model is an appropriate method for predicting wind field and pollutant dispersion in the crowded urban area.

Recent works [26–30] focused on the study of more complex geometries. Gausseau et al. [26] used LES to investigate the effect of two wind directions on near-field pollutant dispersion from stacks located at a low-rise building in downtown Montreal Canada. Their results showed some similarities with the case of dispersion around an isolated building. Van Hoof and Blocken [28] studied the CO2 dispersion in a semi-enclosed stadium and validated their results using onsite measurements. The validated model was used to evaluate the carbon dioxide concentration in the stadium. Amorim et al. [29] chose two European urban areas (Lisbon and Aveiro, Portugal) to evaluate the effect of trees over the dispersion of carbon monoxide (CO) emitted by road traffic. They found that trees are responsible for an increase of 12% CO concentration in the first configuration (nonaligned wind 45◦). In the second configuration (aligned wind), they observed a decrease of 16% of CO concentration.

In light of presented studies, there are many qualitative, as well as quantitative parameters, which need to be studied when investigating pollutants dispersion in an urban area, especially for complex geometries. Hence, the present study aims at providing more insight into the dispersion process by analyzing the effect of buildings arrangement and the material properties (such as density) on the dispersion of hazardous gases. For this purpose, two different configurations, which exhibit complex geometries of Hanover and Frankfurt cities, Germany, are investigated. The most challenging point is the high density of the surrounding buildings, which makes the carrier phase flow and pollutants concentration field complex, anisotropic, and difficult for prediction. Different scenarios of pollutant injection into the computational domain are tested. Dense buildings are included up to a distance of 600 m in diameter and a high-resolution refined grid with over 40 million cells was used. LES simulations are performed, then results are compared against published data from Leitl, B. and Schatzmann. [30] in order to validate the numerical models [31].

#### **2. Methodology**

#### *2.1. Turbulence Model*

In LES, the flow variables are decomposed into the resolved-scale components and the (SGS) components using a filtering operation [32]. The filtering process effectively filters out eddies whose scales are smaller than the filter width or grid spacing used in the computations. Resulting in filtered, continuity and momentum equations, the incompressible Navier–Stokes equations are given as follows:

$$\frac{\partial \overline{\rho}}{\partial t} + \frac{\partial \overline{\rho} \overline{u\_i}}{\partial x\_i} = 0 \tag{1}$$

$$\frac{\partial}{\partial t}(\overline{\rho}\overline{u}\_{i}) + \frac{\partial}{\partial \mathbf{x}\_{j}}(\overline{\rho}\overline{u}\_{i}\overline{u}\_{j}) = -\frac{\partial \overline{p}}{\partial \mathbf{x}\_{i}} + \overline{\rho}g\_{i} + \frac{\partial}{\partial \mathbf{x}\_{j}} \left[\overline{\rho}\overline{\mathbf{v}}\bigg| \frac{\partial \overline{u}\_{i}}{\partial \mathbf{x}\_{j}} + \frac{\partial \overline{u}\_{j}}{\partial \mathbf{x}\_{i}}\right] - \frac{2}{3}\overline{\rho}\overline{\nu}\frac{\partial \overline{u}\_{k}}{\partial \mathbf{x}\_{k}}\delta\_{ij} - \overline{\rho}\tau^{SGS}\_{ij}\bigg] + \overline{S} \tag{2}$$

where bars and over-bars express mean and filtered quantities. The variables *u*, ρ, *p*, *g*, and δ*ij* denote the velocity, the density, the hydrostatic pressure, the gravitational force, and the kronecker delta, respectively.

The quantity <sup>+</sup><sup>ν</sup> represents the molecular viscosity and *<sup>S</sup>* is the source term. The Smagorinsky SGS model is applied to close the system of equations and determine the SGS stresses τ*SGS ij* using the Equation (3)

$$
\tau\_{ij}^{\text{SGS}} = \overline{u\_i}\overline{u\_j} - \overline{u}\_i\overline{u}\_j = -2\nu\_{\text{SGS}}\overline{S}\_{ij} + \frac{1}{3}\delta\_{ij}\tau\_{kk}^{\text{SGS}}\tag{3}
$$

#### *2.2. Dispersion Model*

Numerous models are accessible to predict the dispersion of airborne material. In this study, pollution dispersion was modelled using the species transport approach by calculating convection-diffusion as described by Equation (4). This equation is solved for N-1 species, where N is the total number of species in the fluid flow. The sum of the mass fractions of all species must be equal to unity, the Nth mass fraction is given from the sum of the N-1 solved mass fractions.

The conservation equation takes the following general formula:

$$\frac{\partial}{\partial t}(\rho \overline{Y}\_i) + \frac{\partial}{\partial \mathbf{x}\_j}(\rho \overline{u}\_j \overline{Y}\_i) = \frac{\partial}{\partial \mathbf{x}\_j}(\rho D\_i \frac{\partial Y\_i}{\mathbf{x}\_j}) + \overline{R}\_i + \overline{S}\_i \tag{4}$$

where ρ is the density, +*Ri* is the net rate of production of species by chemical reaction, +*Si* is the rate of creation by addition from the dispersed phase, *Di* is the diffusion flux of species, *Y* +*<sup>i</sup>* is the mass fraction of specie *<sup>I</sup>*, and <sup>+</sup>*uj* is the flow velocity.

#### *2.3. Numerical Procedure*

The commercial CFD software ANSYS FLUENT 17.0 is employed to simulate the airflow and pollutants dispersion. For all configurations, the dispersion of the pollutant is simulated using LES and dynamic Smagorinsky, SGS, along with the species transport model. The bounded central-differencing scheme is used to discretize the momentum equation. A second order upwind scheme is applied for the pollutant concentration and energy equation. Time integration is bounded to second order implicit and pressure interpolation is second order. For both cases, the time step advancement is set to Δ*t* = 10−3*s*, which produces a Courant number "C" less than one within the computational domain. The presented results are averaged over a period of 220 s, which corresponds to five flow-through times.

#### *2.4. Validation of the Numerical Model*

In order to validate the numerical model, the (CEDVAL B1-1) configuration [30] (Figure 1) is investigated using the same boundaries conditions and physical models. This study was conducted in the BLASIUS wind tunnel at the Meteorological Institute of the University of Hamburg [30]. A modified Standing-Spires and uniform LEGO-roughness was used to generate a neutral atmospheric boundary layer ABL at a scale of 1:200 with *Re* = 37,252. First, the boundary layer flow was validated based on detailed flow measurements and comparison with full-scale data. Then, the building models were placed in the test section. Inside the wind tunnel, a 3 × 7 array of buildings composed of rectangular blocks with the same dimensions H = 0.125 m are used to model the dispersion of carbon dioxide (CO2). The CO2 was emitted by four sources, located at the leeward of building source, with exhaust velocity *we* = 0.033 m/s. Validation is achieved by comparing LES results against the published measurement [30]. The horizontal plane at z = 1.5 m (full-scale) height was used as a measurement plane in which pollutants data were acquired. The same experiment is also simulated

by Longo et al. [31]. Identical computational domain, as well as boundary conditions, are used for the present work. The inlet boundary is set to 1 m upstream of the buildings, in which the ABL profiles are measured within the wind tunnel. The outlet boundary is located at 4 m downstream of the last array of the buildings. The width, depth and height of the domain are 1 m, 5 m, and 1 m, respectively, which correspond to the wind tunnel size.

**Figure 1.** Dimensions of the building models, the source building, and the source [30].

In the simulation, a power law velocity profile is used at the inlet, a zero-pressure gradient is used at the outlet and top, no-slip velocity is used at the ground and walls, and symmetry is used for lateral boundary conditions. All the boundaries and initial conditions are summarized in Table 1. Although the temperature is isothermal (300 k), due to thermal stratification, the buoyancy is not present. Yet, we do have stratification due to density difference between carbon dioxide and air.


**Table 1.** Initial and boundaries conditions used for CEDVAL configuration.

An unsteady, 3D, double precision, pressure-based solver with second order schemes for momentum, turbulence quantities, and concentrations are used. The coupled algorithm is selected for pressure-velocity coupling. A transport equation approach is used for the dispersion study. A grid-independent performance of the numerical simulation results is verified by comparing the results of test cases having different grid numbers. Three meshes with different levels of resolution are investigated. A fine resolved mesh of 10 million elements is used for all computations. Simulations are performed using parallels computation of 40 CPUs. The total flow period is set to 12.0 s, which allows the flow to sweep the domain five times.

Figure 2 shows a comparison between the dimensionless turbulent kinetic energy, dimensionless velocity components, and the dimensionless concentration of CO2 obtained by numerical simulation against the wind tunnel experimental data. This figure demonstrates a good agreement between LES numerical results and measurements. The values of turbulent kinetic energy, velocity components, and CO2 concentration are well-predicted and follow the experimental measurements, demonstrating the ability of LES Smagorinsky, SGS, and species transport models in predicting the pollutants dispersion in atmospheric conditions.

**Figure 2.** Comparison between the numerical and experimental results of the normalized turbulent kinetic energy (**a**,**b**), non-dimensional velocity (**c**,**d**), and dimensionless concentration (**e**,**f**).

#### **3. Computational Domain and Boundaries Conditions**

#### *3.1. Geometry and Boundaries Conditions*

In order to represent real case scenarios, two computational domains are generated, one for each case (Figure 3a,b). In case 1 (Hanover city), the focus is placed on the effect of buildings arrangement, whereas in case 2 (Frankfurt city), the behaviors of two hazardous gases with different densities are studied. At the inlet, flow is perpendicular to the injection plane. The streamwise, spanwise, and vertical coordinates are represented by x, y, and z, respectively. The buildings are modeled based on the existing full-scale data. To limit the number of cells, vegetation is not considered. Figure 3a depicts the computational domain used to simulate the configuration of Hanover city sited at the Latitude and Longitude coordinates of 52◦22 18.3"N and 9◦44 23.3"E, respectively. The dimensions of the computational domain are L = 600 m in length, W = 580 m in width, and H = 200 m in height, inside of which various buildings having different dimensions and showing multiple features are included. In the case of Frankfurt city, shown in Figure 3b, the configuration is located at the coordinates 50◦6 49.19"N and 8◦39 51.13"E. The dimensions of the computational domain are L = 500 m in length, W = 500 m in width, and H = 300 m in height.

**Figure 3.** Geometry, boundaries conditions and release locations: (**a**) Geometry and boundaries conditions used for Hanover city, (**b**) Release locations used for Hanover city, (**c**) Geometry and boundaries conditions used for Frankfurt city, (**d**) Release locations used for Frankfurt city.

In case 1, shown in Figure 3a, results highlight the effect of buildings arrangement. For this purpose, two different scenarios dealing with different release locations are investigated. The first one is an accidental leakage of methane. The hazardous gas (methane) is introduced through an opening located at a confined geometry (x = −27 m, y = 26 m, z = 1 m). The choice of this location is based on its closeness to the buildings and/or distance to the exit. CH4 is injected toward the wall of a building from an inlet having 1 m2, a velocity of 8 m/s, and a temperature of 310 K. The second leakage scenario focuses on the dispersion of methane within an open surface, i.e., along the street and far from buildings (x = 110 m, y = −29 m, z = 1 m). Identical boundary conditions for methane are applied at the location 2. For both leakages, CH4 is injected with 100% mass-fraction. Figure 3b shows the effect of gas density for case 2. Two gases having different densities, i.e., methane and carbon dioxide, are used for the investigation. The choice of these gases is based on the density ratio to the ambient air. CO2 is heavier than ambient air (ρ (CO2) = 1.84 kg/m3), while the methane is lighter than ambient air (ρ (CH4) = 0.71 kg/m3). The CO2 and CH4 gases are injected from a chimney sited at (x = 0 m, y = 0 m, z = 20 m) by applying an inlet velocity of 2 m/s and a temperature equal to 310 K.

For the carrier phase boundary conditions, a power law profile is used to describe the variation of the inlet wind speed as a function of the height given in Equation (5) [33]

$$\frac{\imath \mathcal{U}(z)}{\imath \mathcal{U}\_{ref}} = \left(\frac{z}{z\_{ref}}\right)^a \tag{5}$$

where Uref = 9 m/s is the velocity reference at zref = 10 m and *a* = 0.21.

The turbulence was generated using the inlet boundary condition. For LES simulation, the vortex method was used with a number of vortices Nv = 190 in order to establish a time-dependent inlet profile for the inlet velocity. In this method, a number of vortices Nv are generated and convected randomly at each time step.

With this approach, a perturbation is added on a specified mean velocity profile via a fluctuating vorticity field (i.e., two-dimensional in the plane normal to the streamwise direction). The vortex method is based on the Lagrangian form of the 2D evolution equation of the vorticity and the Biot-Savart law. A particle discretization is used to solve this equation. These particles, or "vortex points',' are convected randomly and carry information about the vorticity field. If *N* is the number of vortex points and A is the area of the inlet section, the amount of vorticity carried by a given particle *i* is represented by the circulation Γ*<sup>i</sup>* and an assumed spatial distribution η:

$$\Gamma\_i(\mathbf{x}, y) = \sqrt[4]{\frac{\pi A k(\mathbf{x}, y)}{3N[2\ln(3) - 3\ln(2)]}} \tag{6}$$

$$\eta(\overrightarrow{\dot{x}\prime}) = \frac{1}{2\pi\sigma^2} \Big( 2e^{-|\mathbf{x}\prime|^2/2\sigma^2} - 1 \Big) 2e^{-|\mathbf{x}\prime|^2/2\sigma^2} \tag{7}$$

where *k* is the turbulence kinetic energy. The parameter σ provides control over the size of a vortex particle. The resulting discretization for the velocity field is given by

$$\overrightarrow{\mu}\left(\overrightarrow{\mathbf{x}}\right) = \frac{1}{2\pi} \sum\_{i=1}^{N} \Gamma\_i \frac{\left(\left(\overrightarrow{\mathbf{x}}\_i - \overrightarrow{\mathbf{x}}\right) \times \overrightarrow{\mathbf{z}}\right) \left(1 - \left|\overrightarrow{\mathbf{x}}^\prime - \overrightarrow{\mathbf{x}}^\prime\right|^2 / 2\sigma^2\right)}{\left|\overrightarrow{\mathbf{x}} - \overrightarrow{\mathbf{x}}^\prime\right|^2} \tag{8}$$

where <sup>→</sup> *z* is the unit vector in the streamwise direction. Originally, the size of the vortex was fixed by an adhoc value of σ. To make the vortex method generally applicable, a local vortex size is specified through a turbulent mixing length hypothesis.

The turbulent kinetic energy, *k*, and the turbulence dissipation rate, ε, used for RANS simulations are determined using the Equations (9) and (10) [26]:

*Iu* is the streamwise turbulence intensity (*Iu* = σ*u*/*U,* with σ*<sup>u</sup>* the standard deviation of *u*)

$$k(z) = \left(I\_{\mu} \, \mathcal{U}\right)^2 \tag{9}$$

$$\kappa(z) = \frac{\left(\mu^\*\right)^3}{\kappa(z+z\_0)}\tag{10}$$

The variable *u\** is the friction velocity, κ = 0.41 represents the von Karman constant and *z*<sup>0</sup> = 1 m (full scale) is the aerodynamic roughness length. All the boundaries and initial conditions are summarized in Table 2.

**Table 2.** Initial and boundaries conditions used for Hanover and Frankfurt city.


#### *3.2. Computational Grid*

Figure 4a exhibits the mesh used for Hanover city simulations. A mesh having 40 million elements is generated by the patch independent algorithm. This algorithm is based on the top-down approach, in which volumes are meshed first and cells are projected to faces and edges. It is mainly used for complex geometries that demonstrate various features having different tolerances. In order to well-capture the boundary layer and well=predict the continuous phase physics, an appropriate refinement near the source of CH4 and close to the ground was performed. The grid quality is improved by increasing the cell number. Figure 4b shows the grid used for the Frankfurt city case. In the building zones, the high-resolution computational grids are composed of prism and tetra elements only. Faraway from the building areas, a hexahedral mesh is applied. For both meshes, the resolution is determined using a grid-sensitivity analysis. It is worth mentioning that grid quality varies between 0.05 and 0.95. The metrics used to quantify the quality is the orthogonality, which relates the angles between adjacent elements. The optimum angle is 90◦ for quadrilateral faced cells and 60◦ for triangular faces elements. The mesh quality plays a significant impact on the accuracy of the numerical prediction, as well as the stability of the simulation.

**Figure 4.** (**a**) Mesh used for Hanover city case 40 million cells, (**b**) mesh used for Frankfurt city case 10 million cells.

#### **4. Results and Discussion**

#### *4.1. Hanover City Results*

Figure 5 shows the characteristics of the flow at the section *z* = 3 m (Figure 5a) and at vertical section *y* = −250 m, (Figure 5b,c). The maximum velocity magnitude, 7 m/s, is at the central location section, which is the narrowest passage along the street. This is explained by the conservation of mass, i.e., if a fluid flow having a constant flow rate passes through smaller section, it will be accelerated. As the flow is incompressible, the transport and dispersion of pollutant species is governed by the wind velocity. Behind buildings, zones of recirculation are noticed. Figure 5b,c depict the resulting vertical profiles of the mean streamwise velocity <*u*> and the turbulence intensity (*I*), obtained at (*y* = −250 m) using LES model. The low power profile is well-reproduced and clearly shows the

frictional effect of the ground. For the turbulence intensity, one notices a low variation from 4% (*z* = 4 m) to 6% (*z* = 200 m), which are smaller values compared to the literature [23] that predicted 12%. Consequently, one concludes that inlet turbulence has a relatively low influence on the dispersion of pollutants. In fact, the turbulence at the release location is mainly governed by the presence of the surrounding buildings, which generate velocity gradients, in turn enhancing the turbulence.

**Figure 5.** Characteristics of approaching flow: (**a**) Velocity at *z* = 3 m, (**b**) mean streamwise velocity <u> at *y* = −250 m, (**c**) turbulent intensity at *y* = −250 m.

Figure 6 shows qualitative results of methane dispersion. Two different injection locations are investigated. The release location significantly affects the dispersion of pollutants. When methane is injected for a time period of 220 s, the burnable methane-air mixture, represented by an iso-surface equal 4% of CH4-mass fraction, reaches the wall of the nearby building (Figure 6b). The burnable mixture shows a well-developed vertical plume. Parts of this mixture are in touch with the ground. In this zone, the buildings interfere with the wind velocity, delivering a flow that is mainly characterized by the urban canopy and loosely dependent on atmospheric wind flow. As a result, the methane

concentration increases in the wakes of the building where wind speed is considerably reduced and the recirculation zone becomes remarkably larger. Figure 6a depicts the concentration on buildings surfaces and demonstrates a critical value of 4% close to the wall of the nearby building. From a practical point of view, if a leakage occurs at this location, the ventilation intakes of these buildings should preferably not be located on this face. Indeed, pollutants can contaminate the indoor air of the nearby building if the windows of the leeward faces are open.

**Figure 6.** Qualitative results showing the dispersion of methane CH4: (**a**,**b**) From location 1 sited at (x = −27 m, y = 26 m, z = 1 m) (**c**,**d**) from location 2 sited at (x = 110 m, y = −29 m, z = 1 m).

Figure 6c,d show a different behavior of methane dispersion due to the new location of the pollutant source, which is placed in an open area in which convective transport prevails. Thus, a change in the stack location results in a very different plume behavior. The pollutant species are subjected to a carrier phase flow having anisotropic convective fluctuations. The turbulence intensity along the y-direction is more important than that in all other directions, making the species iso-surfaces show a narrowly opened cone (Figure 6c).

Figure 6d shows the streamlines starting from the CH4 source (location 2), in which the velocity of the carrier phase is the highest and exceeds 7 m/s. CH4 mass fraction shows different pattern compared to location 1. Pollutants are transported horizontally more than that in the vertical direction. The important carrier phase momentum demonstrates a dominant effect. The horizontal time scale of the mixture momentum is clearly smaller than that in the vertical direction. Despite the importance of buoyancy forces, it appears that the dynamic/momentum of the carrier phase is the parameter that controls the motion and dispersion of CH4.

Figure 7a,b depicts the quantitative result of the CH4 dispersion obtained in location 1-simulation. The hazardous zone of an explosion and/or toxicity (mass fraction between 0.04 and 0.14) is bounded between (29 m < x < −28 m), (26 m < y < 27 m) and (2 m < z < 10 m). At the level z = 5 m, the curve bandwidth is less important in "y-direction" than that in "x-direction". The flow diffusion (dominated by the turbulence) is clearly not isentropic. This can be explained by the vortex effect behind buildings, which transports methane in the positive x-direction. Faraway from the source (z = 10 m), the concentration of methane decreases significantly to reach 0.1. Yet, it is strong enough to cause an explosion. Figure 7d,c shows the quantitative distribution of CH4 dispersion produced in the simulation of location 2-case. The hazardous zone is bounded between (109 m < x < 111 m), (−29 m < y < 24 m) and (0 m < z < 2 m). Moving toward the outlet in the positive y-direction, the concentration decreases considerably and falls below 10% after 5 m. This demonstrates that the wind flow transports rapidly CH4 and prevents a pollutants accumulation.

**Figure 7.** Methane (CH4) mass fraction distribution in x and y directions: (**a**,**b**) For location 1 (at z = 2, 5, 10 m); and (**c**,**d**) for location 2 (at z = 2 m).

#### *4.2. Frankfurt City Results*

Figure 8a shows an instantaneous velocity contour of air at level z = 10 m. The velocity magnitude varies between zero and 9 m/s. Close to the buildings, obvious stagnation zones are registered. Velocity gradients are also observed downstream the buildings. Figure 8b depicts the streamlines obtained at level z = 10 m. The streamlines demonstrate a separation in front of the first raw of buildings and an acceleration of the carrier phase. The form of streamlines changes as we move toward the center of the geometry, and recirculation zones are observed behind buildings. Figure 8c,d display an instantaneous concentration of CH4 and CO2 at the ground, respectively. The CO2 concentration is one order of magnitude larger than that of CH4. This concentration difference demonstrates that the deposition of pollutants is mainly controlled by the density ratio and not by the air-velocity. The maximum concentration of both gases is observed behind the source building. A significant pollutant concentration is observed in the vicinity of computational domain outlet. This is due to the horizontal convective transport of the pollutant toward the outlet due to high wind velocity and negligible diffusion.

**Figure 8.** (**a**) Instantaneous velocity contour at z = 10 m, (**b**) streamline at z = 10 m, (**c**) methane (CH4) concentration at the ground z = 0 m, (**d**) carbon dioxide (CO2) concentration at the ground.

Figure 9 depicts the CO2 and CH4 mass fraction variations along vertical lines defined by (x, y) coordinates. Six different locations are chosen to investigate the effect of the density on the dispersion of pollutants. They are located behind the building source sited at (x = 0 m, y = 0 m). All profiles of CO2 mass fraction show larger concentration than CH4. This is caused by the buoyancy force. The heavier gas is dragged downward (to the ground) by the effect of gravitational force. Carbon dioxide CO2 is heavier than ambient air. As a result, its higher concentration is close to the ground.

**Figure 9.** (**a**) CO2 and CH4 mass fraction distribution at different vertical lines, (**b**) vector showing vortex located behind building source (red roof).

Pollutant transport is also influenced by the vortex located behind the building source (see Figure 9c). This zone is created by the separation of fluid flow from the edges of the buildings. The appearance of this zone strongly affects the transport and diffusion of pollutants. Due to the recirculation of the carrier phase (air), methane is dragged downward, and CH4 is accumulated at the ground.

A maximum concentration of CO2 mass fraction is recorded at (x = 11 m, y = 0), which is the closest location to the wall boundary. In this region, the effect of the non-slip condition reduces the velocity magnitudes, resulting in the formation of stagnating zone in which pollutants are collected. Faraway from the recirculation zone (x = 13 m, y = −17 m) and (x = 13 m, y = 17 m), the concentration of pollutants, CO2 and CH4, is very low and mostly negligible.

Figure 10 shows the effect of the density ratio on the dispersion of CH4 and CO2. The density of pollutant-air mixture increases (CH4)/decreases (CO2) rapidly and reaches the density of ambient air at a height of z = 5 m, far from the injection source (Figure 10a). The density ratio controls the shape and the volume of the pollutant plume. Figure 10b,c show iso-surfaces of CH4 and CO2 concentration. In order to compare the dispersion of pollutant gases, an iso-value of 4% mass fraction is represented. The plume of CH4 is first transported vertically by the inlet flows and then carried by the wind horizontally to cover the roof of nearby buildings. In the case of CO2, the dispersion of the plume is located immediately at the top of the building that include the pollutant source. Compared to CH4, the dispersion is remarkably reduced.

**Figure 10.** (**a**) Density ratio ρ*CO*2/ρ*air*, ρ*CH*4/ρ*air* distribution at the vertical line (x = 0 m, y = 0 m) as a function of the height z, (**b**) iso-surface, representing the plume of CH4, (**c**) iso-surface representing the plume of CO2.

#### **5. Conclusions**

High-resolution CFD simulations of near- and far-field pollutant dispersion in a building group in downtown Hanover and Frankfurt are predicted using LES turbulence model along with the dynamic Smagorinsky SGS. The simulations focus on the concentration of pollutants in the vicinity of the release sources and the spatial distribution of the transported species following different building structures and boundary conditions. The following remarks can be drawn:


(4) The main parameter that controls the dispersion in Frankfurt city configuration is the density ratio of CO2 and CH4 to ambient air. Despite the injection of both gases in the vertical direction, only CO2 shows a significant concentration at the ground.

Future works should focus on the existence of vegetation, as well as the implementation of new sources that account for additional species creation, e.g., CO, NOx, SOx, which describe the transport means in the cities.

**Author Contributions:** Conceptualization, M.S.I. and M.C.; Methodology, M.S.I. and M.C.; Software, N.B.S.; Validation, M.S.I., M.C. and N.B.S.; Formal Analysis, M.S.I.; Investigation, M.S.I.; Resources, N.B.S.; Data Curation, N.B.S.; Writing—Original Draft Preparation, M.S.I.; Writing—Review & Editing, M.S.I. and M.C.; Visualization, M.S.I.; Supervision, M.C.; Project Administration, N.B.S.; Funding Acquisition, N.B.S.

**Funding:** This research received no external funding.

**Acknowledgments:** The numerical simulations reported in this paper were supported by the laboratory of the Research Unit Mechanical Modelling, Energies and Materials National School of engineer of Gabes and the CADFEM-AN company.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Synchronized Multiple Drop Impacts into a Deep Pool**

#### **Manfredo Guilizzoni 1,\*, Maurizio Santini <sup>2</sup> and Stephanie Fest-Santini <sup>3</sup>**


Received: 22 June 2019; Accepted: 25 July 2019; Published: 27 July 2019

**Abstract:** Drop impacts (onto dry or wet surfaces or into deep pools) are important in a wide range of applications, and, consequently, many studies, both experimental and numerical, are available in the literature. However, such works are focused either on statistical analyses of drop populations or on single drops. The literature is heavily lacking in information about the mutual interactions between a few drops during the impact. This work describes a computational fluid dynamics (CFD) study on the impact of two, three, and four synchronized drops into a deep pool. The two-phase finite-volume solver interFoam of the open source CFD package OpenFOAM® was used. After validation with respect to high speed videos, to confirm the performance of the solver in this field, impact conditions and aspects that would have been difficult to obtain and to study in experiments were investigated: namely, the energy conversion during the crater evolution, the effect of varying drop interspace and surface tension, and multiple drop impacts. The results show the very significant effect of these aspects. This implies that an extension of the results of single-drop, distilled-water laboratory experiments to real applications may not be reliable.

**Keywords:** multiple drop impact; computational fluid dynamics (CFD) simulation; volume-of-fluid; OpenFOAM; crater dimensions; vorticity

#### **1. Introduction**

A large number of phenomena, both of natural and technological interest, involve the interaction between one or—much more commonly—many liquid drops and a solid surface or a gas-liquid interface. This interaction may be mechanical, thermal, or chemical, and it usually starts with an impact of the drop onto the surface or interface. The impact velocity may be low or high, normal or oblique, and in most cases, a series of drops impact the surface at the same time or nearly simultaneously. For example, this happens in internal combustion engines, firefighting systems, surface cooling, or painting by sprays, pesticide deposition in agriculture, pollution and microbial transport in air and water, forensic analysis of blood stains, and rain effects on planes, wind turbines, and buildings. This, therefore, demonstrates the importance of a deep understanding of this phenomenon, which is still not fully understood despite its apparent simplicity and more than a century of investigation. In particular, many studies can be found in the literature about:


Detailed introductory information and reviews can be found in [1–5]. On impacts into deep pools, the work by Cole [6] is a specifically valuable reference. On the contrary, studies are quite scarce about the intermediate step, i.e., the outcomes of the combined impact of more than one drop in parallel (with simultaneous or delayed impacts) [7–10] or in a series [11–13].

As it may involve solid, one, or more liquids, and a surrounding gas, drop impact is a multiphase fluid mechanics problem, where inertial, viscous, and capillary forces merge their effects. It may also become a multi-physics problem in the case of heat or mass transfer, or chemical reactions may play a significant role (e.g., for hot or cold drops or surfaces or for drops impacting onto a chemically different liquid, for reactive wetting). Therefore, a rigorous model of the phenomenon can theoretically be created, but its analytical solution is not viable. Consequently, very simplified models and experiments have been the primary method of investigation for a long time. In recent years, the increase in the power of calculus has also made computational fluid dynamics (CFD) simulations feasible. Many remarkable studies were developed (e.g., [14–18]) using different approaches, ranging from the most used Eulerian models in finite volume or finite element frameworks, to Lattice Boltzmann and molecular dynamic simulations. Commercial, open-source, and in-house software packages were used. Specifically, the Volume-of-Fluid [19] approach is one of the most used algorithms, even though other techniques were also developed, e.g., Level-Set [20], markers [21], combined models [22–26], and (as already noted) Lattice Boltzmann [27,28] and molecular dynamics simulations (for nanodrops) [29].

In this study, CFD was used to investigate the outcomes of double and multiple synchronized drop impacts into a deep pool with the same temperature as the drops. The numerical campaign consisted of two steps: the first for validation, with simulations carried out to replicate experimental conditions that had been previously investigated by means of high-speed video acquisitions, and the second, to simulate conditions that would have been difficult to obtain in experiments. Energy conversion during the crater evolution and the effects of varying drop interspace and surface tension and multiple drop impacts were investigated. For the first phase, the agreement between simulations and experiments was satisfactory, particularly for some impact conditions. The results of the second phase evidenced the very significant effect of the investigated aspects. This outcome implies that an extension of the results from single-drop, distilled-water laboratory experiments to real applications may not be reliable.

#### **2. Materials and Methods**

#### *2.1. Experimental Setup*

The experimental set-up and procedures, including uncertainty analysis of the main parameters, were described in detail in a previous paper [10], so they will be only briefly summarized here. A thorough analysis of the experimental data against which the simulation results will be validated is available in the same paper. Drop impact experiments were performed by high-speed visualization with continuous back-light illumination. Millimeter-size deionized water droplets were produced with an on-demand droplet generator [7], which allows the synchronous detaching of two drops of equal size. Falling drops are accelerated by gravity, and during their free fall, they are detected by a light sensor, triggering the high-speed (1000 fps) charge-coupled device (CCD) camera acquisitions. Then, the drops impact the pool, which is contained in a cubic box with a side length of 40 mm. The pool depth is one order of magnitude larger than the drop diameter, to avoid any influence on the studied relevant part of the crater evolution. The drops are pigmented by a solution with well-known rheological behavior [30] and with a concentration of 1% in weight, so that the drop water can be distinguished from the pool water in the high-speed video frames, without significantly affecting the water's density and viscosity. For each recorded sequence, the images before drop impact are post-processed by an algorithm that comprises background subtraction, contrast enhancement, image cleaning, edge detection using the Laplacian of Gaussian (LoG) method, and feature detection to determine the drop section and center of each drop [31]. Using consecutive images coupled with the camera acquisition rate and the pixel-to-meter conversion factor, the diameter and velocity of the impinging droplets can

be evaluated with an accuracy of ±0.03 mm and ±0.02 m/s, respectively. The drop sphericity S is also calculated for each falling drop, defined here as the aspect ratio of the vertical to horizontal axis (thus, S > 1 represents a prolate-shaped and S < 1 represents an oblate-shaped drop) before the instant of impact, and assuming a rotational symmetry around the vertical axis. The sphericity of the falling drops before impact is crucial for obtaining accurate and repeatable experimental results, as sphericity heavily affects crater dimensions and impact-generated vorticity [32–35]. Another critical aspect is the oscillation mode of the drop at the instant of impact [36,37]. The experiments were performed at ambient pressure, with the temperature of the pool equal to 30 ± 1 ◦C and the temperature of the drops equal to 27 ± 1 ◦C. The investigated experimental conditions in terms of drop diameter D, inter-drop distance L0, and impact velocity w are summarized in Table 1, along with the values of the dimensionless groups typically used for drop impact characterization.


**Table 1.** Investigated conditions.

D: drop diameter; w: impact velocity, L0: inter-drop distance; S: drop sphericity; We: Weber number ρDw2/σ; Fr: Froude number w2/(gD), where ρ density, σ surface tension, g gravity.

#### *2.2. Numerical Simulations*

Concerning the numerical setup, simulations were performed using the interFoam solver of the OpenFOAM® open source CFD toolbox [38]. interFoam is an isothermal finite volume (FV) solver based on the volume-of-fluid (VOF) method and implementing a model to include surface tension at the interface [39]. OpenFOAM® was selected due to its free and open source nature and because of its many favorable reviews [40–43] and successful cases of use, which are described in the literature [42–48]. The model implemented in interFoam includes the continuity and momentum equations for a Newtonian and incompressible fluid, whose density and viscosity are calculated on the basis of an indicator function named volume fraction γ. The latter assumes a value of 0 for one phase, 1 for the other, and between 0 and 1 in the interfacial regions, and it is transported by the fluid velocity field w. Volume tracking and interface reconstruction (typically as the isosurface at γ = 0.5), with no explicit interface tracking, is, therefore, performed. The advantage of the VOF method with respect to the two-fluid models is that a single set of equations must be solved. The expressions of all the used equations, a discussion of the meaning of the included terms, and further details about the interFoam solver and models can be found in [47,49,50]. Very promising modified versions of interFoam were also presented in the literature [48,51], but the source code was not made publicly available, so the version included in the official OpenFOAM® distribution was used.

As the investigation is at a constant temperature and the impact velocities are low, both phases were assumed to be incompressible and in laminar motion, in agreement with all the previously cited papers in this field. For each phase, the values of all the relevant thermophysical properties were taken at the respective temperatures measured during the experiments. The initial conditions of the two phases and their positions within the domain were fixed using the setFields OpenFOAM® utility. The implicit Euler scheme (first-order accurate) was used for the time derivative, as it proved to offer better results in comparison with the second-order Crank-Nicholson discretization schemes that were tested during preliminary simulations, selecting different blending factors between 0.5 and 1. The conventional advection term was discretized using Gauss schemes: limited Van Leer for the volume fraction and limited linear for the velocity. For the latter, variations of the limiter parameter were tested, but the best results were obtained by keeping it equal to 1. Finally, the OpenFOAM® specific "interface compression" scheme [52] proved to be the best choice for the discretization of the compression term. The maximum allowed Courant–Friedrichs–Lewy (CFL) number was 0.3. In

fact, the solver is allowed to adapt the time step to keep the CFL number under the desired limit, and the requirement for the interface compression scheme used in interFoam for 3D cases is to limit it under 0.3 for 3D cases [40,52]. Some volume fraction sub-cycles are also performed to further improve the accuracy. Three-dimensional simulations were performed within a domain that is a rectangular cuboid having dimensions 0.015 m × 0.015 m × 0.020 m. This domain represents 1/4 of the real pool system, as the computational effort was reduced by exploiting two symmetry planes: the longitudinal plane passing through the two impact points and the transversal plane orthogonal to the previous one, passing through the midpoint between the impact points [10]. To reduce the domain height, drops' detachment from the drop generator and their fall towards the pool surface were not simulated. The drops were directly initialized as spheres near the free surface of the pool, with an initial velocity corresponding to the experimental one. This approach has a weakness in the fact that drop oscillations after the detachment from the needle are not considered. As already mentioned, drop shape and oscillatory behavior have an important influence on the impact outcomes (as it was confirmed also by the preliminary simulations described in the following section), but in the experiments, the drops' shapes were also nearly spherical, so this approximation seemed to be acceptable for these cases. The used mesh was a purely structured one, with hexahedral cells and grading along the axis connecting the drop centers (with a size ratio of 10 between the last and the first cell) to better capture the thin neck that separates the craters created by neighboring impact drops. In the other directions, the aspect ratio was fixed to 1 for all the cells. For the simulations involving three drops with one centered in the origin of the axes, no grading was used, and the mesh resolution was doubled along both the horizontal axes (x,z).

Domain characteristics, together with the imposed boundary conditions, the initial position of the interfaces, and the mesh details are summarized in Figure 1.

**Figure 1.** Domain dimensions, boundary conditions, and mesh details.

Single drop impacts were also performed for comparison, leaving everything equal to the double impact cases apart from the initial drop position, which was centered on the y-axis.

#### **3. Results and Discussion**

#### *3.1. Validation against Experimental Data*

To validate the numerical simulations, crater evolution and drop-formed vortex structures were compared with those observed by means of high-speed video acquisitions. Details about the latter can be found in [10]. All these simulations were performed using the properties of distilled water for the liquid phase. Simulation results were sliced along the longitudinal plane and saved as images with the same pixel-to-meter ratio of the experimental images, then, the time steps corresponding to the experimental frames were selected and superposed to the latter. Figures 2 and 3 show the results, up to time τ = 8 ms for all three investigated cases and for the entire simulated time for cases A (w = 1.0 m/s) and B (w = 1.4 m/s). This comparison is the most severe in evaluating the agreement, as it allows a direct visualization of the correspondence between the experimental features of the craters and the numerical ones. The agreement is very good for all the investigated cases during the inertial phase of the impact (Figure 2) and for case B during the crater recession (Figure 3, column b). On the contrary, some significant discrepancies appear during this second phase for case A (Figure 3, column a). This will be commented on in the following section, when discussing the crater depth profiles. The results after τ = 8 ms are not shown for case C because, with all the used meshes, the two craters merge, such that the crater depth is still reproduced well, but the crater shape is substantially altered, and the agreement is completely lost. With the VOF approach, air regions coming into contact merge instantaneously because the underlying model is not able to represent retarded coalescence (modified models would be needed, e.g., see [53]). On the other hand, real-world interfaces may resist even with direct contact—at least for a certain time—if the pressure of the contacting regions is not too different (e.g., bouncing bubbles or drops). Thus, in the experiments, the very thin neck between the two crates in case C is conserved, while in the simulations, a much finer mesh than that used for the simulations would have been needed to avoid crater contact and merging (OpenFOAM® is very slow in adaptive remeshing, so this feature could not be exploited).

Different mesh sizes were evaluated to analyze the mesh independence of the solution, which, for fine details, is always an issue when using the VOF technique. Mesh sizes were varied between 176 kcell and 12.35 Mcell, but meshes under 1 Mcell proved unsuitable because the already mentioned neck between the craters also disappeared (resulting in a coalescence of the cavities not observed in the experiments) for cases A and B. Thus, only the finest two meshes (1.4 Mcell and 12.35 Mcell) were used for the final simulations. For the latter, a good overall mesh independence was found, even if some details are still captured differently. No mesh is always the most accurate (as can be seen in Figures 2 and 3 for case B).

The vortices generated by drop impacts—peculiar and different from the single-drop impact case [10]— were visualized in the experiments by dying the drop water. In this way, the dye highlights the drop water, distinguishing it from the pool water, and allowing to follow its motion. In the simulations, this could be accomplished by numerically "dying" the drop water with a passive tracer, by using a multiphase solver (e.g., multiphaseInterFoam) or simply by tracing the flow streamlines [18]. An example of the latter approach is shown in Figure 4, which includes two pictures referring to case A at τ = 18 ms. These two images were obtained by tracing streamlines with different origins, thereby evidencing how the streamline aspect heavily depends on the latter. The complex structure of the simulated flow streamlines is consistent with the experimental observation of vortex structures in the same region [10]. This behaviour is also in agreement with the experimental findings by other authors [54], who showed that after a certain time from the drop impact (for the investigated case, this time should be around 14 ms) the original azimuthal vorticity is progressively tilted and becomes a stream-wise vorticity. Despite the validity of this approach, using it to cover all the domain regions with good resolution would require a very large number of streamlines, such that they would be difficult to follow. Therefore, another visualization technique was preferred. Pictures showing the calculated vorticity (the curl of the velocity vector field) on the longitudinal plane were compared with

the experimental images. Even if the vortex and vorticity are not related in a biunique way (due to irrotational vortices or vortex parts), vorticity offers a picture of the regions where vortices are most likely to occur, with a resolution equal to the number of mesh cells. Figure 5 shows the results for cases A and B. As can be seen, the agreement between simulations and experiments is also good for regions of high vorticity. The discrepancies are due to the differences in the crater evolution between the simulations and experiments. These differences influence the motion of the drop water under the crater itself. In any case, they are not significant enough to undermine the correct identification of the regions reached by the drop water and its swirling motion.

**Figure 2.** Superposition of the numerical crater contours to the experimental high-speed photos, for the three investigate cases A (w = 1.0 m/s), B (w = 1.4 m/s), C (w = 2.0 m/s), time up to 8 ms. In case B, the red line shows the results of the 12.35 Mcell simulation, while the yellow line shows the results of the 1.4 Mcell simulation.

**Figure 3.** Superposition of the numerical crater contours to the experimental high-speed photos, for case A and case B (for the latter: red line, results from the 12.35 Mcell simulation; yellow line, results for the 1.4 Mcell simulation); time up to 18 ms.

**Figure 4.** Simulated flow streamlines in a region where vortices are experimentally observed (case A, w = 1.0 m/s, τ = 18 ms).

**Figure 5.** Comparison between the experimentally detected vortex structures and the vorticity calculated by numerical simulations, for cases A and B, time up to 18 ms.

For case A in Figure 5, the drop fluid is symmetrically distributed below the liquid–air interface (τ = 4 ms). Later (τ = 6 ms), the drop fluid begins to roll up. The continuous deformation of the liquid–air interface generates an outer vortex structure in the form of a vortex filament, which is also visible in the numerically calculated vorticity in the same region. The roll-up is continuous while the expanding crater causes downward convection of the dyed drop fluid. Asymmetries of the outer vortex rings, having their origins in different convection velocities, are also found in the numerical results. When the crater recedes (τ = 9 ms), a central vortex structures is formed. Both vortex structures, the central vortex and the outer vortex, separate from the receding crater and form a unique complex vortex blob that continuous to move downward. This complex structure could not be captured in the simulation, but the region of vorticity is represented very well.

For case B in Figure 5, the generation of the outer vortex structure is superposed, and the dyed liquid spreads around the cavity. The interaction between the two evolving craters causes a local enlargement of the spreading dyed drop liquid, showing a vortical structure in both investigations (the experiment and simulation (τ = 8 ms)). When the cavities recede, apparently disorganized structures separate and move downward. The numerical results clearly show these downward moving structures.

In terms of the characteristic dimensions of the drop craters, the crater's maximum penetration depth zmax, the maximum crater penetration depth in the transverse plane including the impact point zc (i.e., parallel to the yz-plane and at a x-coordinate equal to half of the drop interspace [10]), and the crater penetration depth along the vertical line including the impact point zuip, were measured from both the experiments and the simulations. In this work, these quantities are always presented in a dimensionless form, by dividing them by the drop diameter before impact. Therefore, their symbols are completed by a "dl" subscript in all the following charts. For double drop impacts, these quantities do not always coincide during the crater evolution (differently from what happens for single drop impacts). This can be clearly seen, e.g., in Figures 2 and 3. zmax and zc are accessible both to experiments and simulations, while zuip can be evaluated only by simulation because its measurement requires a cross-section of the craters (along the xy plane) and not just a side view. Figures 6 and 7 show the values of these quantities for case A (w = 1.0 m/s) and case B (w = 1.4 m/s) as a function of time. As for all the charts in this paper, the initial time τ = 0 was fixed at the instant of impact. In the experimental results, the two craters are not perfectly equal, but their difference is very slight (as can be seen for zmax for case B in Figure 7). Therefore, only the values for the left crater are shown in Figure 6. More specifically, Figure 6 reports the results for case A (w = 1.0 m/s), comparing simulations and experiments for the two values of the water–air surface tension set in the simulations.

**Figure 6.** Dimensionless zmax, zc and zuip extracted from the experiments and from the numerical simulations for case A (w = 1.0 m/s): (**a**) the results using σ = 71 mN/m; (**b**) the results using σ = 65 mN/m.

**Figure 7.** Dimensionless zmax, zc and zuip extracted from the experiments and from the numerical simulations for case B (w = 1.4 m/s). volume-of-fluid (VOF) results using σ = 71 mN/m.

As can be seen in Figure 6, the experimental values are better reproduced when σ = 65 mN/m (Figure 6a) instead of σ = 71 mN/m (distilled water value, Figure 6b). This seems reasonable as the dye used for the droplet in the experiments [10] contains surfactants, which may cause—despite the very dilute solution—a reduction in the drop surface tension. For case B, reported in Figure 7, the agreement is already good for simulations with σ = 71 mN/m. This can be explained by the fact that in case B, the impact velocity is higher, so the importance of the capillary actions and surface tension is reduced with respect to case A.

In any case, the overall maximum depth reached by the crater during its evolution—which is one of the most important quantities because it limits the water layer directly reached by the crater—is always reproduced very well by the simulations. The deviation when using σ = 71 mN/m is −6.1% for case A, −1.7% for case B, and −0.14% for case C. The deviation for case A is reduced to −0.68% when using σ = 65 mN/m. In all cases, the numerical simulation underestimates the experimental value, with a discrepancy that decreases with an increasing impact velocity. This is again consistent with the fact that the set surface tension is probably larger than the correct experimental value (thus, the free surface more strongly opposes pulling).

#### *3.2. Energy Conversion During Impact*

After validation, the numerical simulations were used to extract information about some aspects that would be difficult or very time-consuming to study experimentally. Post-processing of the numerical results was first carried out to extract the kinetic and surface energies within the domain for all the simulated time steps. Given the small involved heights, the changes in the potential energy (of the drop liquid and of the pool liquid displaced by the drop) are minor, so the latter will not be described. Kinetic energy was calculated for each cell as:

$$\mathbf{E}\_{\mathbf{k}} = 1/2 \left| \mathbf{w}\_{\text{cell}} \right|^2 \left[ \left( \mathbf{y} \text{ \textquotedblleft} \mathbf{w}\_{\text{water}} + (1 - \mathbf{y}) \text{ \textquotedblright} \right) \right] \tag{1}$$

where wcell is the velocity magnitude in the cell, ρwater and ρair are the water and air densities, and γ is the volume fraction. Summation over all cells was then performed to calculate the total Ek over the whole domain, which, before impact, is obviously equal to the kinetic energy of the drop alone. Surface energy Es was calculated as the product of the air–water interfacial energy σ (taken at its reference value for distilled water at the drop temperature, 71 N/m) by the total area of the air–water interface obtained as the isosurface of the color function γ = 0.5. This value was then decremented by the value of the total interfacial energy of the undisturbed pool surface, so that a direct comparison between the interfacial energy of the drop before the impact and the interfacial energy during crater evolution can be evaluated.

For both kinetic and surface energy, the values extracted from the numerical results at the time step immediately before impact were compared with the theoretical predictions for drops of known diameters and velocities to check their correctness and reliability. The average deviation was 0.78% for kinetic energy and 1.31% for surface energy, which was considered satisfactory. Figure 8 shows the results for both the double drop impacts and the corresponding (i.e., with same impact velocity) single drop impact, after normalization of the values with respect to the kinetic energy and the surface energy immediately before the impact. The time evolution of the two energy contributions appears consistent with the experimentally observed behavior and with each other. When kinetic energy increases, surface energy decreases, and vice versa. Surface energy decreases during the first instants after impact because the drop coalesces with the pool (so a part of the drop and the pool's external surfaces are merged, with a consequent reduction in the total interface area), but the crater is not formed yet. The profiles also show a very regular evolution for impact velocity. The most interesting outcome is that the evolutions for single and double impacts are extremely similar, even if minor differences are obviously present due to the more complex evolution of the crater shape in double drop cases (e.g., causing the mutual intersections between the curves). This suggests that the latter has only a minor influence on energy conversion during the crater's expansion and recession.

**Figure 8.** Time evolution of the normalized kinetic energy (**a**) and normalized surface energy (**b**) in the domain after double drop impact, as a function of time for single (dash-dot lines) and double (continuous lines) drop impacts.

#### *3.3. Sensitivity Analysis on Surface Tension*

As water is easily contaminated, particularly in real applications, its surface energy may be quite different from the distilled water case. Given the observed influence of σ on the results, a sensitivity analysis was carried out on its value. Case B (w = 1.4 m/s) was selected as a reference case, changing only the surface tension value. It can be expected that larger surface tension values cause a more significant influence of the capillary action and thus a higher internal pressure within the original drop and a "tenser" pool interface. As can be observed from Figure 9, the trend of the inertial phase, in which the crater is expanding, is influenced only slightly by the water surface tension. The maximum expansion of the crater is reduced when the surface tension is larger (meaning that the tensing effect dominates with respect to the larger capillary pressure in the drop). The receding phase is heavily changed, and for larger surface tension values, the crater closes much faster. A few sudden changes in the path can be observed. These changes correspond to instants when the capillary waves caused very fast modification in the crater evolution. In fact, once the crater reaches its maximum depth and width, its starts receding, and the surface tension reduces the interface area. This causes the closure of the crater and perturbations of the interface that propagate as waves along the crater interface, resulting in peculiar shapes (e.g., see the time instants between τ = 10 ms and τ = 16 ms in Figure 3) and seemingly irregular trends in the charts. In particular, the changes in the maximum crater depth evidence the instants when different crater regions become deepest.

**Figure 9.** Time evolution of zmax dl as a function of the water surface tension for case A (w = 1.0 m/s, (**a)**) and for case B (w = 1.4 m/s, (**b)**).

#### *3.4. Sensitivity Analysis on Drop Spacing*

The second parameter that was varied is the distance between the drop centers. This distance was changed in the range of 2.58 to 6.58 mm by modifying the reference drop distance (L0 = 4.58 mm) of ±1 and ±2 mm. This resulted in a series of dimensionless drop distances equal to L/D = 1.127 ÷ 2.873, L/L0 =0.345 ÷ 1.437. Figure 10 shows the results in terms of zmax (a) and zuip (b), including also single drop impacts for comparison. Here, too, the sudden changes in the paths are evident. Paths are consistent between the different cases, with a progressive shift of the instants at which the crater's main parameters change their trends due to the combined actions of inertia and capillarity. As expected, when the drop distance is large, the crater tends to behave as in a single drop impact. On the other hand, reduction of the interspace causes the deformation of the craters, which interact with each other first indirectly and then directly merge, affecting more and more the crater dimensions. A peculiar result is that for drops nearly touching each other, the time at which the maximum value of zmax is reached tends to the value for a single drop with double the volume and mass (as the density is not changed), but the maximum value of zmax dl is much larger than that in the corresponding single drop case. Thus, it can be concluded that direct interaction (merging) between the craters enhances crater penetration even more than a simple summation of inertial contributions would predict.

**Figure 10.** Time evolution of zmax dl (**a**) and zuip dl (**b**) as a function of the distance L between the drops. Reference case (L0/D = 2) is case B; the other cases are with the drop interspace decreased and increased by 1 mm and 2 mm. Single drop impact (sdi) includes the same drop diameter and impact velocity; the single drop impact for a drop with the same impact velocity but a volume and mass double of the original one (sdi2V) is also added.

#### *3.5. Three-Dimensional Visualization*

Unlike the experiments, where a single side or top view of the drop impacts can be observed, from the numerical simulations, a fully 3D description of the phenomenon can be obtained. First, the air–water interface can be extracted and visualized. Figure 11 shows a rendering of the crater evolution at different time instants—for case B (w = 1.4 m/s), simulated with 12.35 Mcell. For a more intuitive visualization, the simulated domain was mirrored once to reconstruct half of the real domain. It can be noticed how not only the major features, but also the details (including the neck and the surface and capillary waves) are very well captured.

**Figure 11.** Rendering of some frames showing the crater evolution from the simulation of case B (w = 1.4 m/s), simulated with 12.35 Mcells.

#### *3.6. Impact of Three and Four Synchronized Drops.*

Finally, simulations of multiple drop impacts were tested in all cases for impact velocity w = 1.4 m/s. For these cases, no corresponding experimental tests were performed. Three drops in a single row, aligned along the x-axis, and four drops, two aligned along x and two along z, were tested. Four cases were simulated—one with a drop interspace L equal to L0 for both three- and four-drop cases, a three-drop case with a drop interspace equal to 0.75 L0, and a four-drop case with a drop interspace of 1.5 L0. Figure 12 shows some rendered frames of the crater's evolution for three of these simulated cases. It can be observed how for the three-drop case with L = L0 the craters do not merge, while for the other cases they do, resulting in a different evolution and a maximum crater depth that becomes much larger, as reported in Figure 13. In these cases, too, the interaction between the craters appears to have a significant effect on the maximum depth. Further studies will be needed to verify how multiple drop interactions would affect the crater depth when there is no "free space" left for expansion. On the other hand, each crater finds its way blocked by other interfaces in all directions. In this case, the crater depths might be reduced. Further studies will also be needed to investigate the effect of non-synchronous impacts.

**Figure 12.** Rendering of some frames showing the crater evolution from the simulations of three-drop and four-drop cases (w = 1.4 m/s).

**Figure 13.** Dimensionless zmax as a function of time from the numerical simulations for multiple drop impact cases (three drops, tdi, and four-drops, qdi) having w = 1.4 m/s and different inter-drop spacing. Results for single drop impacts for drops having the same volume (sdi) and two times the volume (sdi2V) are also added for comparison.

#### **4. Conclusions**

Numerical simulations of double and multiple synchronised drop impacts were performed using the interFoam solver of the OpenFOAM® open-source CFD package. Results for double drop impacts were validated against experimental data acquired by high-speed imaging. Crater shape and depth and vortex generation were compared between the simulation and experiments for double drop impacts, and a good agreement was found, particularly in the inertial phase (leading to the maximum crater penetration depth). Simulations were then performed to investigate energy conversion during the crater evolution, the effect of varying surface tension, and drop interspace and multiple drop impacts. The main findings are:


In summary, direct extension of results from the study of single drop impacts to the case of multiple impacts with many unevenly spaced drops and non-simultaneous contact times (much more common in applications), appears not to be straightforward; on the contrary, neglecting mutual interactions can be quite misleading. It seems, therefore, mandatory to perform further studies including both more drops and different drop interactions (delayed impacts, different drop interspaces, etc.) to cast further light on these aspects. This would lead to the possibility of developing new models and correlations for multiple drop impacts, which are more reliable for use in all the applications fields where the phenomenon is of importance.

For such studies, improvements in the VOF approach may also be needed:


In addition, comparison with other methods (e.g., diffuse interface methods) could be interesting.

**Author Contributions:** Conceptualization, M.G., M.S., and S.F.-S.; Formal analysis, M.G., M.S., and S.F.-S.; Funding acquisition, S.F.-S.; Investigation, M.G., M.S., and S.F.-S.; Methodology, M.G.; Software, M.G.; Validation, M.G., M.S., and S.F.-S.; Visualization, M.G.; Writing—original draft, M.G., M.S., and S.F.-S.; Writing—review and editing, M.G., M.S., and S.F.-S.

**Funding:** This research received no external funding. The Research Grant of Stephanie Fest-Santini was financed within "Progetto ITALY® —Azione Giovani in Ricerca" of the University of Bergamo (Italy).

**Acknowledgments:** The contribution of Viola Knisel in performing the numerical simulations is gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
