A Numerical Validation of 3D Experimental Dam-Break Wave Interaction with a Sharp Obstacle Using DualSPHysics

The presence downstream of a dam of either rigid or erodible obstacles may strongly affect the flood wave propagation, and this complex interaction may lead to further dramatic consequences on people and structures. The open-source Lagrangian-based DualSPHysics solver was used to simulate a three-dimensional dam-break in a closed domain including an oriented obstacle that deflects the flow, thus increasing the complexity of fluid dynamics. By comparing numerical results with experimental data, the effectiveness of the model was evaluated and demonstrated with an extensive sensitivity analysis based on several parameters crucial to the smoothed particle hydrodynamics method, such as the resolution, the boundary conditions, and the properties of the interaction weight function. Charts and summary tables highlight the most suitable conditions for simulating such occurrences in the DualSPHysics framework. The presence of the obstacle, being also an opportunity for observation and study of complex fluid dynamics, opens the way to investigate the fluid interaction with solid objects involved in dam-break events and, possibly, to predict their effect with respect to the relative position between them and the flood and other relevant parameters. Finally, the numerical model presents a good overall agreement.


Introduction
With all their benefits, dams also carry a certain risk of failure that, in most cases, results in catastrophic consequences: from the damaging and/or destruction of buildings and infrastructure to the loss of human lives [1]. The experimental and numerical investigation of dam-break events has always been a central task, and several studies have been carried out to predict the resulting flow dynamics (e.g., [2]).
The presence of obstacles such as buildings, bridge piers, and roads plays an important role in the water depth distribution, velocity propagation, and, more in general, the flow regime. It also implies a greater complexity of the numerical modeling, even in the fixedbed case, due to the interaction of the fluid wave with the impacting structures (e.g., [3,4]). The orientation of these constructions, with respect to the main flow direction, influences the flow peculiarities as well. Over all of these aspects, numerical simulations can provide vital information on the impact force exerted by an approaching dam-break wave [5][6][7][8][9][10][11].
The treatment of the presence of the rigid obstacle was traditionally tackled with simple numerical approaches, such as the shallow water equations (SWEs) [11,12], but further simulations of propagation of dam-break waves in urban areas [13][14][15] have shown that significant errors occur where the flow is strongly three-dimensional (3D). Aureli and co-workers [16] compared the predictions of a two-dimensional (2D) depth-averaged model, a 3D Eulerian two-phase model, and a 3D smoothed particle hydrodynamics (SPH) model in estimating the load exerted by a dam-break wave on a rigid squat structure. The results show that the 2D model based on the shallow-water approach leads to an error in the peak load in the order of 10% compared with the experimental data. The need for accurate predictions on the flow behavior motivated this study with a 3D investigation, thanks also to the increased computing power: the numerical simulation are, hence, executed in a SPH-based framework.
The SPH method belongs to the family of meshfree particle methods (MPMs), which present several advantages over mesh-based methods, such as the natural ability to solve multi-mechanics problems. Its Lagrangian nature allows computing, with relative ease, multi-mechanics problems involving complex interfaces and moving boundaries. By now, the SPH method is well established in the computational fluid dynamics (CFD) discipline, and its effectiveness has been proven in several works, especially when dealing with free-surface flows and large deformations [17][18][19][20][21][22][23][24].
DualSPHysics [25] is an advanced meshless solver with emphasis on free-surface flow modeling and has been shown to be robust and accurate in a wide range of applications: reproducing extreme wave events [26], coastal engineering simulations [27,28], fluid-solid interaction (FSI) [29,30], wave energy converters [31][32][33][34][35], etc. The peculiarities of this solver make it suitable for dam-break simulations, since it is naturally able to handle large deformations and violent impacts that involve solid obstacles [36,37]; it follows that SPH-based numerical models could be as effective as, if not better than, the traditional methods applied to the dam failure records, such as the ones based on 2D SWEs [38,39] or the VOF-or FEM-based [40,41] models.
In this work, small-scale laboratory experiments conducted by Kocaman et al. (2020) [3] are compared to a Lagrangian-based numerical model to examine the problem with a single obstacle. A numerical investigation of a 3D dam-break impacting a sharp obstacle was carried out with the SPH-based solver DualSPHysics. The DualSPHysics code was tested with respect to multiple aspects crucial to the SPH method to prove the effectiveness of the latter in reproducing phenomena usually tackled with grid-based methods. A wide sensitivity analysis based on several parameters, such as the resolution (initial interparticle distance), the interaction domain dimension, and the boundary conditions (BC) treatment, was conducted against the experimental and numerical results depicted in [3]. The effectiveness of the SPH method in reproducing complex free-surface fluid dynamics, shown in this study, can lead the way to the simulation of more complex problems not easily solvable with the traditional finite difference method (FDM), but for which the SPH method is naturally suitable, e.g., debris flow combined with FSI. This paper is structured as follows. In Section 2, the experimental set-up provided by Kocaman et al. (2020) [3] is reported. In Section 3, the SPH formulation for the governing equations implemented in DualSPHysics is presented, along with the boundary conditions treatment. Section 4 exposes the numerical implementation of the model, the crucial parameters taken into account, and the data analysis process. In Section 5, several charts and tables describe the results of the work, divided according to the prevalent fluid dynamics typology, including an introductory FSI analysis. Section 6 summarizes the main results of this research, giving an overall evaluation and detailing on the effects of the analyzed parameters, presenting an outlook of possible future investigations.

Experimental Set-Up
In this section, the reference experimental set-up is presented. The experiments were performed on a rectangular horizontal channel, 1.00 m long, 0.50 m wide, and 0.35 m high (Figure 1a). Bottom and lateral walls were made of thick glass. The upstream water domain at rest was 0.25 m long and h w = 0.15 m high, mimicking the presence of a reservoir. The water domain was confined by a transverse wall s = 0.01 m thick, made up of two lateral fixed parts and a central mobile gate, the motion of which activated the dam-break flow.
The sluice gate, made of acrylic, was placed with its center of gravity aligned with the streamwise plan of symmetry. Waterproofing was achieved by using grease oil over the vertical guides. In order to simulate the sudden removal of the closing gate, a 15 kg mass connected with a steel rope was set to fall down from an elevation of 1 m, generating a fast and complete gate opening. The opening time was about 0.06 s from video recordings, a period smaller than the recommended value for an instantaneous opening [42], 1.25(h w /g) 1/2 = 0.15 s, with g being the gravity acceleration. Before performing any measurement, the upstream reservoir was filled with water. Dry bed conditions were realized at the downstream instead. The rectangular obstacle, 0.15 m long and 0.08 m wide, was located 0.25 m downstream the mobile plate with a rotation of θ = 28.0724 • , to have one of its diagonals overlap the flow direction ( Figure 1a). In order to detect the fluid motion, yellow dye was added to the water before removing the gate. Five ultrasonic sensor (MIC + 35/IU/TC, range 65-300 mm, at location P1, and MIC + 25/IU/TC, range 30-250 mm, at locations P2-P5) sampled the water levels (Figure 1b), and a high-speed 300 fps CCD camera recorded the experiment evolution. Moreover, the numerical investigation [3] produced with a volume of fluid (VOF) and Fractional Area/Volume Obstacle Representation (FAVOR TM ) model, which uses a k − ε method to close Reynolds-Averaged Navier-Stokes (RANS) equations [43], was used as an extra reference in this work.

SPH Formulation
MPMs, in general, refer to the class of meshfree methods that employ a set of finite number of discrete particles to represent the state of a system and record its movement. Each particle can either be directly associated with one discrete physical object or be generated to represent a part of the continuum problem domain. For CFD problems, each particle possesses a set of field variables such as mass, momentum, energy, positions, etc. and other variables (vorticity, etc.) related to the specific problem. The advantages of the MPMs over conventional grid-based numerical methods can be roughly summarized as follows: • The problem domain is discretized with particles without a fixed connectivity, so treatment of large deformation is relatively easier.

•
Discretization of complex geometry is simpler as only an initial discretization is required.

•
It is easy to obtain the features of the entire physical system through tracing the motion of the particles, and, therefore, identifying free surfaces, moving interfaces, and deformable boundaries is no longer a tough task.
Among the MPMs, the smoothed particle hydrodynamics (SPH) method was employed in the present work.

Principles of the SPH Method
The strategy in SPH is to discretize the physical domain (fluid and/or solid objects) into a set of particles, where the physical quantities (position, velocity, density, and pressure) are obtained as an interpolation of the corresponding quantities of the surrounding particles. The contribution of those particles is weighted using a kernel function, with an area of influence that is defined using a characteristic smoothing length. This discretization process is divided into two key steps [44].
The first step is the integral representation or the so-called kernel approximation of field functions, consisting in the integration of an arbitrary function and a smoothing kernel function. The integral representation of a generic spatial function f (r) within an integral volume Ω is given by: where W is the smoothing kernel function or kernel and r is the position vector. In the smoothing function, h is the smoothing length defining the influence area of W ( Figure 2). The second step is the particle approximation. The integral representation is approximated by summing the values of the nearest neighbor particles, which yields the particle approximation of the function at a discrete point (particle). The position vector r b is defined as the position of a particle having a fixed mass m b and a finite volume V b , related by: where ρ b is the density of particle b = 1, . . . , N p , in which N p is the number of particles within the support domain of the particle a. The integral representation can be written in a discrete form for a particle a, being b part of its support domain: where The kernel function, hence, plays a fundamental role in the SPH method. In Dual-SPHysics, the Wendland [45] quintic kernel function can be utilized: and α D,n is a constant depending on the dimension of the problem.

Governing Equations
The governing equations for a fluid domain in the Lagrangian formulation and the relative SPH discretization are presented. The fluid dynamics equations are based on the fundamental physical laws of conservation: • conservation of mass; • conservation of momentum.
These laws yield, respectively, to the continuity equation and the momentum equation where ρ is the density, u is the velocity vector, g is the gravitational acceleration vector, p is the pressure, and Γ is the deviatoric stress tensor. The governing equations in the Weakly Compressible SPH (WCSPH) formulation implemented in DualSPHysics become: The state equation [46]: where γ = 7, ρ 0 = 1000 [kg/m 3 ], and c 0 = c(ρ 0 ), along with Equations (9) and (10), completes the set of governing equations for the fluid domain implemented in DualSPHysics [25]. The artificial viscosity term, Π ab , is added in the momentum equation based on the Neumann-Richtmeyer artificial viscosity, aiming at stabilizing the SPH scheme [47]. A density diffusion term is implemented in DualSPHysics to improve the stability of the scheme by smoothing the density field [48]. This formulation is based on the density diffusion term developed with the name of δ − SPH by Antuono et al. (2012) [49].

Boundary Treatment
Dynamic boundary conditions (DBC) [50] for solid surfaces are standard in Dual-SPHysics [51]. A stationary solid is simply represented by fixed particles with pressure evaluated from the equation of state. Boundaries are easy to set up and computations are relatively stable and efficient, providing robust numerical simulations for complex geometries. However, a small nonphysical gap between the fluid and solid boundaries can form, of the order of the smoothing length h [52], decreasing the accuracy of pressure measured on the boundary, especially when fluid particles approach dry walls.
On the other hand, an advantage of these boundary particles is their computational simplicity, since density, and therefore pressure, can be calculated inside the same loops as fluid particles with a considerable saving of computational time. DBC has been successfully applied to coastal engineering problems, discretizing complex 3D geometries without the need for implementing complex mirroring techniques or semi-analytical wall boundary conditions.
To enhance this formulation, a new boundary method has been developed: the modified dynamic boundary conditions (mDBC) [53]. In this method, the density of solid particles, arranged in the same way as in the original DBC, is obtained from ghost positions within the fluid domain by linear extrapolation. Thanks to an additional boundary interface, located half a particle spacing from the layer of the closest boundary particles to the fluid, a ghost node is mirrored, with respect to this interface, into the fluid domain (similarly to Marrone et al. [54]). In this ghost node, the fluid properties are found through a first-order consistent SPH spatial interpolation [55] over the surrounding fluid particles only.
Once the density and its gradient are computed at the ghost nodes, the density of the boundary particle ρ b is obtained by means of a linear extrapolation with the found values: in this way, the boundary density is presented as part of a fluid continuum and pressure from the equation of state gives smoother and more physical pressure fields, avoiding the non-physical gap between the boundary and the fluid observed when using the DBC.

Numerical Set-Up
In the following paragraphs, the SPH numerical results regarding the water levels at five points, as reported in Figure 1b, are paired with data from ultrasonic sensors and with numerical results [3]. The numerical set-up of DualSPHysics (Figure 3) reproduces the experimental geometry described in Section 2. To reproduce the gate motion, the upward acceleration value of g is assigned in the solver.

Reference Parameters
In order to validate the model built in the DualSPHysics framework, simulations were performed with respect to several parameters: resolution, support domain dimension, and BCs. Let us define the first two quantities: • Resolution: It is the initial interparticle distance, or dp, defined as number of particle n p per separation wall thickness s dp = s n p where n p = {2; 3; 4} and s = 1 cm (12) • Domain multiplicative coefficient: It is the multiplier κ of the smoothing length h, which defines the interaction domain for each particle; the smoothing length depends on dp and κ according to Lastly, the BCs, which corresponds to one of the two treatments introduced in Section 3.3, are used in two ways: 1. The mDBC are applied on every surface in the domain (indicated in the charts with 'mDBC'). 2. The tank walls (wet surfaces at t = 0) are treated with DBC, while the dry walls with mDBC, as visible in Figure 4 (indicated in the charts with 'DBC').
It is necessary to clarify that, to use the mDBC, at least three layers of particles, in the normal direction with respect to the fluid, are requested. The resolution of the simulation, hence, represents the discriminant factor for the choice of the BC treatment, because of the number of particles present in the separation wall, the thinnest element within the solid domain.

Data Analysis
The postprocessing calculation and data analysis are here exposed. To better understand the results, the numerical computation of elevations in DualSPHysics, provided by the postprocessing tool 'MeasureTool', is reported. For a given position (x, y), the mass at several points a, regularly interspersed on the vertical axis z, is computed with: where m b stands for the mass of the particles included in the support domain.
For elevations for which the support domain is completely immersed in the fluid, Equation (14) gives a unitary result, as a percentage of the assigned fluid mass m = max(m a ) ∀a ∈ fluid (reference quantity), while the latter becomes null when the support domain, for certain z values, is not filled with any fluid particle. The free surface of the fluid is calculated when m a assumes half the value of reference quantity m: m a = 0.5m (15) It follows that some divergences might occur and oscillations in the water level measurements, being the SPH particles prone to a certain dispersion, especially in phenomena that involve high velocities and impacts with solid structures, such as the one here investigated. The influence area of the kernel, together with the initial particle spacing, affects the output value of elevations provided by the software.
The performance of the SPH-based model is investigated analyzing the differences in the water elevation time histories. The local error is computed, for a visual purpose, as: where z DSPH t and z EXP t are, respectively, the numerical and experimental water elevations at the t-th time interval.
To evaluate the accuracy of each combination, and to weight each parameter contribution, the index of agreement W * and the linear correlation coefficient R * [56] are reported for every comparison. These quantities are evaluated following the equations: N ∆t being the total number of steps. The overvalue operator is the mean operator, while σ is the standard deviation operator. In Equation (17), first introduced by Willmott et al. [57], the values are bounded in [0, +1], while, in Equation (18), the range is [−1, +1]; in both cases, values close to the upper ceiling correspond to good agreement.

Results
The following charts and tables report a sensitivity analysis in order to appreciate the effectiveness of the SPH approach; the results are divided in three main comparisons to highlight the influence of the aforementioned parameters on various numerical aspects. Each water level chart is matched with the graphic of the local error, to give an immediate visual perception of the distance between the simulations and the experimental references.

Tank Emptying Analysis
The first analysis focuses on the tank emptying following the gate motion, described by the water level at P1 and P2. These points, being aligned with the axis of symmetry of the domain, describe the basic dam-break fluid dynamics, and thus the comparison evaluates the influence of the dimension of the kernel support domain (namely, the κ coefficient). In this case, the mDBC conditions are applied to every surface of the closed tank, and the resolution is: dp = s 3 The DualSPHysics simulations perfectly reproduce the tank emptying, as visible in Figures 5 and 6. In the latter figure, the change in the agreement is evident when the fluid dynamics switches from simple outflow to complex undulations and hydraulic jumps, caused by the wave reflected by the upstream channel-end and sidewalls; the main difference between experimental and numerical results can be observed immediately after the arrival of the reflected wave, approximately after 3 s ( Figure 6). This difference is partly due to the mono-phase simulation performed, while the experimental data are clearly referred to a multi-phase flow motion. The dimension of the support domain (κ) does not influence significantly the emptying of the tank (Figure 5), while a higher value of κ seems to improve the reproduction of the impact with the reflected wave; however, the overall agreement is better with the unitary value of the multiplying coefficient, as highlighted in Table 1.  The registrable error is very limited, and it is interesting to note the small differences in the initial water levels of the two DualSPHysics simulations (Figure 5), as proof of the parameters influence on the measurement of the free surface, also in stationary conditions. The same fluid dynamics, with different value of resolution and with the DBC conditions applied on the wet surfaces, was evaluated, since the current dp does not allow to use the mDBC on every surface. In Figure 7, the elevation chart with dp = s/2 at P2 is presented, since it was noted that at P1 the overall agreement of the simulation with the experimental data is always good. It is clear that a lower resolution costs a divergence between the results, which generally overestimate the elevation; however, in this case, an increment of κ worsens the effectiveness, not only when the fluid dynamics becomes more complex, but on average for all the duration of the simulation. The effect, instead, of the DBC appears to be not relevant, demonstrating what is reported in Section 3.3. To finally conclude this analysis, the influence of resolution is highlighted by the charts in Figure 8, which compares at P2 the simulations with dp = s/3 and dp = s/2; the index of agreement and the linear correlation error are reported in Table 2.

Fluid Propagation Dynamics
In the survey points P3 and P4, being placed laterally with respect to the obstacle, the fluid dynamics cannot be considered as a simple dam-break flow anymore, since the flood is diverted by the oblique surfaces. The fluid propagation is hence studied focusing on the sensitivity of the model with respect to the resolution.
The first 1.5 s at P3 (Figure 9), in which propagation after the impact with the obstacle is observed, both resolution values, dp = s/4 and dp = s/3, return a proper agreement, but, at the stage of the impact with the return wave, higher resolution enhances the performance and reduces the spread and, thus, the linear correlation error (Table 3). It is safe to say that the DualSPHysics model is robust in reproducing the phenomenon, the agreement being very similar, and proper, with different resolutions. At P4 (Figure 10), the complexity of the dynamics, involving several water jumps and lateral jumps, costs higher error and dispersion, but again the simulation can be considered representative of the phenomena. Despite some high peaks of error, the gain of fidelity with a higher resolution appears clear, with obvious costs in terms of computation time.

Impacts on the Downstream Structure
Among the themes mentioned to introduce the dam-break problem, the presence of the obstacle is surely the most interesting: the complex dynamics examined so far is just one of the aspects that such an occurrence encourages to examine in depth. The prediction of the effects on the downstream structures, in fact, could help reduce risky scenarios, deepening the knowledge of the stress conditions the structures are subjected to when involved in flood events. In the following, two different physical quantities, pressure and force, which are missing in the experimental dataset, are studied. Figure 11 shows the time evolution of the pressure on the four lateral surfaces of the obstacle, computed at the midpoint of each one, at a height of 0.03 m (the location of these points is depicted in Figure 4); the data are computed from the case with dp = s/4, κ = 1.2 and mDBC. The pressure swiftly increases at points A and D, placed on the faces impacted by the first wave, which show comparable peak values, with a small shift in time for point D, farther away from the moving gate. Around t = 3.00 s, the impact or the return wave at points A and C can be noticed, as well as in Figure 8, which depicts the water surface close to point A. As the water fills the enclosed domain, the pressure values tend to the hydrostatic one, although the simulated time-span is not long enough to achieve the final stillness condition (about 15 s are needed).
In Figure 12, a sensitivity analysis for the total force exerted by the fluid on the obstacle is presented. Four different configurations are displayed, with two different values for each reference parameter: resolution, support domain dimension, and BCs. The force evolution displays small variation moving among the four series and the main impact events are apparent; it is not immediate to locate the return waves since the force is computed on the whole structure, and each face is interacting with different fluid masses characterized by different velocities. The simulations with higher resolutions present less divergence in the values, as predictable. The peak values are very similar among the different simulations.

Overall Analysis
The last chart ( Figure 13) gives a general overlook on the DualSPHysics model, comparing simulations with the highest resolution and different values of the support domain dimension. At P5, positioned on the axis of symmetry of the tank but behind the obstacle, the dynamics can be considered a hybrid between the two aforementioned ones. The peculiarity of this situation is evidenced in Figure 13, which includes the numerical reference solution from Kocaman et al. [3]. The agreement between the SPH model and this additional reference is meaningful at P5 because of the complexity of the fluid dynamics; it follows that the numerical models give comparable outcome when this phenomenon is reproduced. In fact, the simulation with dp = s/4 and κ = 1, proved to be efficient in tank emptying and pure propagation (Figure 9), appears here, in the second half of the simulation, evidently to underestimate the water level. The same resolution, but with a larger support domain, reproduces very well the mean elevation when the fluid movement tends to quieten. The κ value appears to be fundamental when the flow motion gains significant 3D characteristics, even more than the resolution. With such complex dynamics, high resolution and large support domains are required, because the influence of any parameter cannot be neglected. In Figure 14, a comparison between experimental and numerical simulation at several initial time steps is presented, showing good agreement in the standard dam-break tank emptying and in the propagation after the wave hits the obstacle. The diverted fluid, then impacting the lateral and bottom walls of the enclosed domain, suddenly changes its characteristics, resulting in an series of water jumps; despite being a single-phase simulation, and so not taking into account the air bubbles included after the wave reflection, the numerical model appears to well reproduce the fluid dynamics exhibited by the laboratory experiments. In Table 4, a general overview of the simulations carried out is presented, with the final mean values for index of agreement and linear correlation error. Figure 14. Visual comparison between live experimental frames (left) and DualSPHysics simulation with dp = s/4, κ = 1.0, mDBC (right). White and red lines highlight the fluid mass shape for easier understanding.

Conclusions
The results show that the numerical investigation carried out in the DualSPHysics framework can be considered efficient overall, describing properly the various physical situations that occur: • The tank emptying is very well reproduced with all the attempted values for the parameters ( Figure 5).

•
The FSI is well computed, and the mDBC has demonstrated its suitability for violent impacts of fluid masses with solid dry objects.

•
The agreement in the propagation dynamics is remarkable and increases with the resolution (Figures 9 and 10).

•
The non-negligible multi-phase dynamics observable at P5 is represented with good agreement, especially considering the increasing importance of the air cavities when the dynamics complexity raises, and mainly depends on the resolution and the kernel support domain size κ ( Figure 13).

•
The samples presented in Figures 11 and 12 constitute a valid starting point for further FSI investigation.
The results in Table 4 show that the simulations with the highest values of resolution perform well, but, as exposed, each parameter taken into account influences a particular aspect of the fluid behavior, as shown by the differences in agreement, with the same parameters, at each measurement point. To not increase the computational cost, a balance has to be found in dependence of which aspect is more relevant with respect to the aim of the study. In conclusion, the results indicate that DualSPHysics can be a useful engine to simulate 3D dam-break phenomena, which also involve specific conditions for FSI during flow impacts, and this research provides vital information to initialize a global setup. This tool can be used for modeling fluid-driven debris impacting fixed or moving obstacles, and it clearly overcomes the limitations of mesh-based numerical solvers.
Future analyses will include a deep study of the obstacle behavior and the fluidand debris-induced forces at the time of impact, along with a study about the influence of the position on possible damages to civil structures. The effectiveness of the present framework can be proven with respect to a wider range of terms of comparison, as shown in an introductory way in Section 5.3, using further experimental and numerical references.

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

Abbreviations
The following abbreviations are used in this manuscript: