Shockwaves in spillways with the particle finite element method

Noname manuscript No. (will be inserted by the editor) Shockwaves in spillways with the particle finite element method. Fernando Salazar ·Javier San-Mauro · Miguel ´ Angel Celigueta ·Eugenio O˜nate Received: date / Accepted: date Abstract Changes in direction and cross-section in supercritical hydraulic channels generate shockwaves which result in increase in flow depth with reF. Salazar ·J. San-Mauro ·M. ´ A. Celigueta ·E. O˜nate Centre International de M`etodes Num`erics en Enginyeria (CIMNE). Universitat Polit`ecnica de Catalunya (UPC), Campus Norte UPC, Gran Capit´an s/n. 08034. Barcelona, Spain Tel.: +34-934-071-495 F. Salazar E-mail: [email protected]c.edu J. San Mauro E-mail: [email protected]c.edu M.´ A. Celigueta E-mail: [email protected]c.edu E. O˜nate E-mail: [email protected]c.edu 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 2 Fernando Salazar et al. gard to that for uniform regime. These disturbances are propagated downstream and need to be considered in the design of the chute walls. In dam spillways, where flow rates are often high, this phenomenon can have significant implications for the cost and complexity of the solution. It has been traditionally analysed by means of reduced-scale experimental tests, as it has a clear three-dimensional character and therefore cannot be approached with two-dimensional numerical models. In this work, the ability of the particle finite element method (PFEM) to reproduce this phenomenon is analysed. PFEM has been successfully applied in previous works to problems involving high irregularities in free-surface. First, simple test cases available in the technical bibliography were selected to be reproduced with PFEM. Subsequently, the method was applied in two spillways of real dams. The results show that PFEM is capable of capturing the shockwave fronts generated both in the contractions and in the expansions that occur behind the spillway piers. This suggests that the method may be useful as a complement to laboratory test campaigns for the design and hydraulic analysis of dam spillways with complex geometries. Keywords Particle finite element method ·Spillways ·Shockwaves 1 Introduction Shockwaves are frequently generated in dam spillways by uniform flow disturbances due to the presence of piers, curves or changes in the cross-section [34]. As a consequence, local maxima in flow depth are produced whose magnitude 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 3 can be much greater than the incoming uniform depth. This has important practical implications in the dimensioning of the chute walls, which must be designed with greater height to adequately convey the flow, with the obvious consequences in terms of magnitude and cost. The design of the spillway must tend to the constant cross-section and slope with the aim of approaching uniform flow and thus avoiding this phenomenon. However, these geometric particularities are frequently unavoidable due to the topographical and/or geological characteristics of the dam and spillway site. The hydrodynamic study of spillways has historically been carried out using small-scale experimental tests. The computational methods available during the period of greatest intensity in dam design and construction did not allow these complex problems to be addressed, since they are not treatable with two-dimensional models [22]. However, the numerical tools currently available can be used for studying this type of 3-D problems with sufficient detail and with an assumable computational cost. Numerical models are often used for the design and analysis of spillways in combination with physical model tests [31], [30]. This allows the numerical model to be validated for the geometry considered, so that the versatility of the latter can subsequently be used to extend the catalogue of situations analysed in alternatives studies. In addition, numerical models allow a more detailed analysis of aspects that are difficult to measure in the laboratory, such as stream lines or pressure in areas at risk of cavitation [32], [36]. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 4 Fernando Salazar et al. Numerical models have been primarily used to evaluate the discharge capacity of spillways [36], [26] and to optimize the approach conditions [31], thanks to the ease with which flow patterns can be analysed. The numerical schemes employed include eulerian [2], [26], and lagrangian approaches such as smoothed particle hydrodynamics (SPH) [29]. However, the analysis of shockwaves is much less frequent. Some examples with quasi-horizontal bottom have been addressed with depth-integrated approaches, which offered accurate results (e.g. [37], [33]). However, spillway chutes often feature steep slopes, which limits the applicability of 2D simplifications. The particle finite element method (PFEM) is a lagrangian approach for the numerical modelling of fluid dynamics problems among many other applications [1]. The fact that the PFEM does not need a background mesh, and thus the user does not need to foresee what parts of the domain will be occupied by the fluid in later stages of the simulation, makes it especially appropriate to consider phenomena with strong changes in the free surface both in time and space [9], [8]. The method has been employed to face a variety of problems in different fields of engineering, such as seakeeping [12], landslide-generated waves [25], [27], industrial forming processes [18], ground excavation [3], fluid-structure interaction [13], among others [15],[14], [7]. PFEM has also been previously applied to solving 3-dimensional free-surface flows, particularly in hydraulic structures [10]. It has also been employed to estimate air flow demand in bottom outlets [28], [20]. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 5 This paper explores the possibilities of using PFEM in full-scale hydrodynamic analysis of spillways with geometric irregularities that generate shockwaves. Firstly, it has been applied to reference cases available in the classical literature including curves and channel contractions. Then, two examples of application to real dam spillways are presented, with different characteristics: the first features a mild slope and an abrupt convergence, reducing the useful width of the channel by more than 50 % in a length similar to the initial width. In the second case, the convergence is smoother, although it allows for a detailed analysis of the negative waves generated behind the piers. The article is completed with a brief reminder of the basis of the PFEM implementation used in this work and a discussion on the possibilities and limitations, as well as ideas for future work. 2 Numerical model The PFEM is a documented method in the literature [19], well suited for continua with large deformations and separation. In the PFEM, the domain is modelled using an Updated Lagrangian formulation [38]. That is, all variables are assumed to be known in the current configuration at time t. The new set of variables in the domain are sought for in the next or updated configuration at time t+∆t. The finite element method (FEM) is used to solve the equations of continuum mechanics. Hence a mesh discretising the domain must be generated in order to solve the governing equations in the standard FEM fashion. When the mesh gets too distorted due to the large deformations it must be 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 6 Fernando Salazar et al. rebuilt to avoid inverted elements. The re-meshing step can be used to improve the mesh quality as well. The PFEM is therefore a compound method consisting in a combination of the Finite Element Method, an Updated Lagrangian formulation and a re-meshing algorithm. The equations to be solved for a single phase fluid (like water) are the Navier-Stokes equations for incompressible fluids: Momentum conservation ρDui Dt =−∂ ∂xi p+µ∂ ∂xj∂ui ∂xj+ρfi(1) for i, j =x, y, z Mass conservation ∂ui ∂xi = 0 (2) for i=x, y, z with u= ¯u(3) for the wall nodes (considered as fixed fluid nodes) and p= 0 (4) for the free surface fluid nodes. In the equations above, the unknowns are p(pressure) and u(velocity). ρ and µare the fluid density and dynamic viscosity, respectively, uiare the velocities along the ith global (cartesian) axis, fiare the volumetric accelerations 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 7 (only gravity acceleration is used in this work), and ¯uis the prescribed velocity. Note that fixing the three components of the velocity to zero is equivalent to assuming a no-slip boundary condition. This type of boundary condition was used in every fluid-wall interphase. The Lagrangian formulation implies that the nodes move with the fluid velocity. In particular, the nodes located at the free surface follow the evolution of the boundary in a very accurate, natural way. According with the PFEM technique [19], equations 1 and 2 are discretised with a standard FEM mesh and then solved by means of an implicit, Fractional Step method. The elements used are simplex, 4-noded tetrahedra, with linear interpolation of both the pressure and the velocity. When the finite elements get very distorted, the mesh is re-generated, but the nodes and their information are conserved. Adaptive mesh refinement techniques can be used to improve the solution in zones where large deformations of the fluid occur. The details of the algorithm and our implementation was described in previous publications [13], [16], [17]. In this paper, only the basic steps of the algorithm are succinctly described, together with some enhancements specifically implemented for the present application. 2.1 Basic steps of the PFEM In the PFEM, the mesh nodes in the fluid and solid domains are treated as particles that contain all the information of geometry, material and mechanical properties of the underlying subdomains. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 8 Fernando Salazar et al. A typical solution with the PFEM involves the following steps. 1. The starting point at each time step is the cloud of points Cin the fluid and the walls. The walls domain simply contains fluid nodes which have an externally imposed velocity and can be moved as a rigid body. For instance, nCdenotes the cloud at time t=nt(Fig. 1). 2. The domain is discretised with a finite element mesh nMusing the particles as the mesh nodes. We use an efficient mesh generation scheme based on the Delaunay tesselation [9]. 3. The free surface is detected by means of the Alpha Shape Method [6], a geometrical criterion which removes the elements coming from the Delaunay Tesselation that are not considered as part of the fluid domain. The Alpha Shape Method compares the diameter of the circumsphere of each element with a reference, maximum allowed diameter. Bigger circumspheres correspond to either big elements (connecting distant nodes) or flat, distorted elements which can only be present in areas where the density of nodes drops considerably (in areas with a homogeneous density of nodes the Delaunay Tesselation generates high quality tetrahedra). Note that the free surface is not the border between two fluids, but the external boundary of the domain considered fluid. 4. The Lagrangian equations of motion for the overall continuum are solved using the standard FEM. The state variables in the next (updated) configuration for nt+∆t are computed: velocities, pressure, strain rate and viscous stresses. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 9 5. The mesh nodes are moved to a new position n+1Cwhere n+ 1 denotes the time nt+∆t, in terms of the time increment size. 6. Go back to step 1 and repeat the solution for the next time step to obtain a new n+1C. Fig. 1: Sequence of steps to update a “cloud” of particles (nodes) representing a domain containing two fluids and a solid boundary from ntto n+1t. 2.2 Inlets Fluid inlets are necessary when modelling spillways or channels as an upstream condition. Being the PFEM a fully Lagrangian method, the inlet condition 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 16 Fernando Salazar et al. Where: θis the angle formed by the chute wall and the direction of the incoming flow, β0is the shockwave angle and F0is the Froude number before the contraction. Reinauer and Hager [24] proposed expressions to calculate the maximum flow depth on the walls and axis, as well as the location of the first maximum: X1= (2 + 0.126α)θ+ (1.4−0.015α) (9) h1=h01 + 1 √2S02 (10) h2=h01 + √2S02 (11) h3=h0w−1+ 1.86S0.5 0−0.2α0.6(12) w=be b0 (13) S0=θF0(14) Where X1is the location of h1,h1..3is the value of the three first maxima, alternatively located on the wall and the axis. αis the channel slope. θis the angle formed by the wall with the initial flow direction. h0,F0and b0are the flow depth, the Froude number and the channel width before the contraction, respectively. beis the channel width downstream of the contraction. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 17 Table 4: Parameters for the contractions analysed Case C1 Case C2 Channel slope θ(rad) 0.186 0.169 Angle of the chute wall R3.25 m 4.25 m Channel width at inlet/outlet b0(m) 0.5/0.3 0.5/0.15 Flow depth at inlet h00.05 m 0.05 m Froude number at inlet F06 6 Inflow velocity v04.2 m/s 4.2 m/s Mesh size 1.2 cm 1.2 cm Number of elements 1.29 1061.11 106 Two of the settings available in the literature have been selected to be modelled with the PFEM to assess its capability to reproduce shockwaves in these configurations, with different conditions. The parameters defining both tests are shown in Table 4. The maximum mesh size in both cases was set to 12 mm. This value ensured having at least 3 elements in depth in the whole domain. Also, only half of the original geometry was considered, taking advantage of the symmetry of the design. 3.2.1 Contraction case C1 Figure 6 shows a general view of the numerical result, where the flow depth maxima can be identified. They are alternatively located on the wall and axis, as expected. The theoretical values of X1,h1,h2and h3have been computed with Equations 9, 10, 11 and 12. A direct comparison with the numerical results 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 18 Fernando Salazar et al. X1=0.532 m h1=0.160 m h2=0.333 m h3=0.178 m Fig. 6: Contraction case C1. The expected empirical values for X1,h1,h2and h3were placed over a view of the numerical result. cannot be made, since the intrinsic properties of PFEM result in irregular shape of the free surface, with some elements separating from the main fluid body. Therefore, planes –horizontal for h1,h2yh3and vertical for X1– have been drawn and superimposed on the numerical results for comparison. It can be seen that all results lie within the correspondent location. The plan view of the numerical result (Figure 7) allows approximating the resulting angle β1formed by the shockwave front, which in this case presents some curvature. The estimated direction (18◦) is 5% lower than the empirical estimate computed with Equation 8 (19.2◦). 3.2.2 Contraction case C2 In this case, the results of the PFEM were compared to the experimental measurements published by Reinauer [21]. For that purpose, the profiles of the 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 19 β1=18º 0.4 0.3 0.2 0.1 0 Flow depth (m) Fig. 7: Contraction C1. Numerical results (plan view). The shockwave front features a small curvature. Its direction can be estimated as forming 18◦with the main flow direction. This is 5% lower than the empirical estimate computed with Equation 8 (19.2◦). free surface along the chute wall and axis were extracted and superimposed over the experimental ones in Figure 8. The results coincide qualitatively in terms of the alternation of maxima between the axis and the wall, including the local minimum in the wall immediately downstream of the contraction. The main discrepancy is observed in the second maximum on the axis, which moves downstream in the numerical model. As in the previous case, the angle of the shockwave front has been approximated and compared to the empirical value given by Equation 8. The result is shown in Figure 9. The difference in β1is coherent with the location of the flow depth maxima in the numerical model, slightly downstream as compared to the empirical estimate (Figure 8). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 20 Fernando Salazar et al. 0.4 0024 0.4 0024 dimensions in meters Fig. 8: Contraction case C2. Side view of the numerical results as compared with the flow depth measured in the experiment (red lines) [21]. Side view. Top: Chute wall. Bottom: Channel axis β1=17º 0.4 0.3 0.2 0.1 0 Flow depth (m) Fig. 9: Contraction case C2. Numerical results (plan view). The direction of the shockwave front can be estimated as forming 17◦with the main flow direction. As in the previous case, the empirical estimate given by Equation 8 (18.2◦) is 5% higher. 4 Real scale tests In this section, two cases are presented in which the PFEM is applied to analyse the hydraulic behaviour of spillways with different degree of convergence. Both share geometric features resulting in the formation of shockwaves. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 21 4.1 Case study #1 4.1.1 Model description The first spillway analysed includes a chute with a mild slope (11.1 %). It features four spans 13 m wide regulated by gates supported by three piers: the central pier is 3.5 m wide and has a rectangular profile at the downstream end, and the other two are 2 m wide, with downstream hydrodynamic profile and shorter length. The total width of the initial section is 59.5 m and after a transition along 55 m it joins the discharge channel, which has a constant useful width of 21 m. In this case, the side channels feature sloping bottom to direct the flow towards the axis. In addition, they are slightly elevated with respect to the central channel, which in turn is aligned with that of the downstream chute. Given that the focus of this study is on the shockwaves generated by the change in cross-section, the energy dissipation structure was not analysed. Figure 10 shows the geometry and the main dimensions of the model. A discharge rate of 1,800 m3/s has been considered. This value corresponds to one of the experimental tests for which data are available (the actual flow in the laboratory was obtained by applying hydraulic similarity to the aforementioned 1,800 m3/s). The gates remain fully open while evacuating this flow, therefore the weir operates in free flow regime. This situation has been reproduced in the numerical model by applying an inlet boundary condition corresponding to the critical regime (flow depth = 4.98m; inflow velocity 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 22 Fernando Salazar et al. 21 m 59.5 m 55 m 350 m Fig. 10: Geometry of the numerical model for Case study #1 = 6.95m/s), as described in Section 2.2. A constant mesh size of 0.5 m has been selected in accordance with the expected flow depth, resulting in about three million elements once the flow is stabilised. These magnitudes constitute a compromise between a sufficient level of detail to reproduce the flow depth variations with a file size that allows for efficient post-process, and with moderate computation time. 4.1.2 Expected behaviour In addition to the convergence of the whole section, which goes from 59 to 21 m as mentioned above, the central channel also has a convergence so that its useful width goes from 26 m at the inlet (two bays) to 21 m of the chute. As there is also a step between the bottom of this area and that of the outer spans, additional shock waves are expected to occur. Both convergences are indicated in Figure 11. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 23 Main convergence Secondary convergence Fig. 11: Left: Convergence of the chute in Case study #1. Right: Expected shockwaves It is well known that contractions of this type with supercritical flows generate shockwaves propagating downstream, similar to those described in Section 3.2. The angle β1formed by the wave front can be obtained with Equation 8. According to this theoretical formulation, two pairs of wave fronts should appear, corresponding to the two width changes shown in Figure 11, which will propagate downstream. 4.1.3 Numerical results The numerical modelling performed with PFEM has been capable to reproduce the formation of shockwaves observed in the physical model and described by Reinauer and Hager [24], both in the walls and in the spillway axis, due to the confluence of the shock waves (Figure 12). The complex geometry of this spillway prevents local maxima from being estimated using the formulas mentioned above. As an example, h2in this case results from the combination of the increase in flow depth produced by the 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 24 Fernando Salazar et al. h1 h2 h3 Fig. 12: Local maxima in the chute alternates between the walls and the axis. convergence and that due to the expansion behind the central pier, known as roostertail [34], [35]. The numerical results have been compared to those obtained in the laboratory. Figures 13 and 14 include different views of the experiments and the corresponding numerical results. In the latter, the magnitude of the velocity component perpendicular (along the Yaxis in the reference system used) to the main flow direction (Xdirection) is plotted for better visualisation of the shockwaves. This and other variables, which require complex instrumentation to be measured in laboratory, can be easily analysed in the numerical models. The results presented for this case study show PFEM’s ability to qualitatively capture the phenomenon. For a quantitative evaluation, they have been compared to the theoretical values of the angle formed by the wave front, and to the flow depth measured in the laboratory with the considered flow rate. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 25 Fig. 13: Comparison between numerical and experimental results for Case Study #1. View from upstream. Fig. 14: Comparison between numerical and experimental results for Case Study #1. View from downstream. The theoretical direction of the wave front β1can be obtained with Equation 8 by entering the value of the incoming Froude number and the angle formed with the main direction of the channel. With the features of this case study, a β1value of 27.2◦is obtained for the secondary convergence. In Fig1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 32 Fernando Salazar et al. it is important to emphasize that the angle of the shockwave front fits with the theoretical estimate (Figure 15). The Case Study #2 has major limitations regarding quantitative evaluation of the results, as there are no experimental data to compare with. Furthermore, the geometry at the end of the piers does not correspond to those for which empirical formulas are available (rectangular, semi-elliptical). In spite of that, the results fit qualitatively with the observed behaviour, and the local depth maxima due to the roostertails behind the piles are observed, with magnitudes within the expected ranges taking into account the specific geometry. 6 Summary and conclusions The possibilities of PFEM to model shock waves in hydraulic structures with supercritical regime have been shown in different configurations: laboratory tests corresponding to curves and convergence have been reproduced, and subsequently 2 spillways of real dams featuring different geometries have been modelled. The results show that PFEM is a useful tool for numerical modelling of this type of phenomena, offering a good approximation with an assumable computational cost, particularly for the first shock wave. The main limitation of the presented approach is the lack of consideration of the interaction with air, which causes increase in flow depth in a supercritical regime [4]. In previous works on hydraulic assessment of bottom outlets, interaction between water and air was considered and meaningful results were obtained [28], [20]. However, consideration of air entrainment in spillways im1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 33 plies using a large domain of analysis to model the surrounding air, together with a dense mesh in the interface between both fluids, to account for their interaction. This results in excessive computational cost with the current version of the code. The authors are currently working in optimizing this PFEM implementation to consider aeration. The mentioned works on bottom outlets and other precedents [12] are the basis of the new developments. The same problem arises in many laboratory tests, where aeration cannot be adequately modelled due to scale effects [11]. In spite of that, laboratory tests have been extensively used in the design and evaluation of dam spillways. Therefore, the model presented can be useful as a complement to experimental tests. The application of PFEM in these settings could reduce costs, enhance the analysis by considering flow patterns and stream lines and, above all, enable a wider spectrum of possible solutions to be considered, thanks to the greater flexibility of numerical models for the generation of geometries that can be laborious to build in the laboratory. Acknowledgements The research was supported by the Spanish Ministry of Economy and Competitiveness (Ministerio de Econom´ıa y Competitividad, MINECO) through the project CALA (RTC-2016-4581-5). 7 Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest. 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 34 Fernando Salazar et al. References 1. Aubry, R., Idelsohn, S., O˜nate, E.: Particle finite element method in fluid-mechanics including thermal convection-diffusion. Computers & Structures 83(17–18), 1459–1475 (2005) 2. Bayon, A., Valero, D., Garc´ıa-Bartual, R., L´opez-Jim´enez, P.A., et al.: Performance assessment of openfoam and flow-3d in the numerical modeling of a low reynolds number hydraulic jump. Environmental modelling & software 80, 322–335 (2016) 3. Carbonell, J.M., O˜nate, E., Su´arez, B.: Modeling of ground excavation with the particle finite-element method. Journal of Engineering Mechanics 136(4), 455–463 (2009) 4. Chanson, H.: Stepped spillway flows and air entrainment. Canadian journal of civil engineering 20(3), 422–435 (1993) 5. Chow, V.: Open-Channel Hydraulics. McGraw-Hill, USA (1959) 6. Edelsbrunner, H., M¨ucke, E.P.: Three-dimensional alpha shapes. ACM Transactions on Graphics 13(1), 43–72 (1994) 7. Idelsohn, S., Mier-Torrecilla, M., O˜nate, E.: Multi-fluid flows with the particle finite element method. Computer Methods in Applied Mechanics and Engineering 198(33), 2750–2767 (2009) 8. Idelsohn, S., O˜nate, E., Del Pin, F.: A lagrangian meshless finite element method applied to fluid-structure interaction problems. Computers and Structures 81(8), 655–671 (2003) 9. Idelsohn, S., O˜nate, E., Pin, F.D.: The particle finite element method: a powerful tool to solve incompressible flows with free–surfaces and breaking waves. International Journal for Numerical Methods in Engineering 61(7), 964–989 (2004) 10. Larese, A., Rossi, R., O˜nate, E., Idelsohn, S.R.: Validation of the particle finite element method (PFEM) for simulation of free surface flows. Engineering Computations 25(4), 385–425 (2008) 11. Lian, J., Qi, C., Liu, F., Gou, W., Pan, S., Ouyang, Q.: Air entrainment and air demand in the spillway tunnel at the jinping-i dam. Applied Sciences 7(9), 930 (2017) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 35 12. Nadukandi, P., Servan-Camas, B., Becker, P.A., Garcia-Espinosa, J.: Seakeeping with the semi-lagrangian particle finite element method. Computational Particle Mechanics 4(3), 321–329 (2017) 13. O˜nate, E., Idelsohn, S.R., Celigueta, M.A., Rossi, R.: Advances in the particle finite element method for the analysis of fluid–multibody interaction and bed erosion in free surface flows. Computer Methods in Applied Mechanics and Engineering 197(19), 1777–1800 (2008) 14. O˜nate, E., Rossi, R., Idelsohn, S.R., Butler, K.M.: Melting and spread of polymers in fire with the particle finite element method. International Journal for Numerical Methods in Engineering 81(8), 1046–1072 (2010) 15. O˜nate, E., Celigueta, M.A., Idelsohn, S.R.: Modeling bed erosion in free surface flows by the particle finite element method. Acta Geotechnica 1(4), 237–252 (2006) 16. O˜nate, E., Celigueta, M.A., Idelsohn, S.R., Salazar, F., Su´arez, B.: Possibilities of the particle finite element method for fluid–soil–structure interaction problems. Computational Mechanics 48(3), 307–318 (2011) 17. O˜nate, E., Franci, A., Carbonell, J.M.: Lagrangian formulation for finite element analysis of quasi-incompressible fluids with reduced mass losses. International Journal for Numerical Methods in Fluids 74(10), 699–731 (2014) 18. O˜nate, E., Franci, A., Carbonell, J.M.: A particle finite element method for analysis of industrial forming processes. Computational Mechanics 54(1), 85–107 (2014) 19. O˜nate, E., Idelsohn, S.: The particle finite element method. an overview. International Journal of Computational Methods 1, 267–307 (2004) 20. Pozo, D., Salazar, F., Toledo, M.: Modeling the hydraulic performance of the aeration system in dam bottom outlets using the particle finite element method. Revista Internacional de M´etodos Num´ericos para C´alculo y Dise˜no en Ingenier´ıa 30(1), 51–59 (2014). [in Spanish] 21. Reinauer, R.: Kanalkontraktionen bei schiessendem abfluss und stosswellenreduktion mit diffraktoren. Ph.D. thesis, ETH Zurich (1995) 22. Reinauer, R., Hager, W.H.: Supercritical flow behind chute piers. Journal of Hydraulic Engineering 120(11), 1292–1308 (1994) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 36 Fernando Salazar et al. 23. Reinauer, R., Hager, W.H.: Supercritical bend flow. Journal of hydraulic engineering 123(3), 208–218 (1997) 24. Reinauer, R., Hager, W.H.: Supercritical flow in chute contraction. Journal of hydraulic Engineering 124(1), 55–64 (1998) 25. Salazar, F., Iraz´abal, J., Larese, A., O˜nate, E.: Numerical modelling of landslidegenerated waves with the particle finite element method (PFEM) and a non-Newtonian flow model. International Journal for Numerical and Analytical Methods in Geomechanics (2015) 26. Salazar, F., Mor´an, R., Rossi, R., O˜nate, E.: Analysis of the discharge capacity of radialgated spillways using cfd and ann–oliana dam case study. Journal of Hydraulic Research 51(3), 244–252 (2013) 27. Salazar, F., O˜nate, E., Mor´an, R.: Numerical modelling of landslides in reservoirs via the particle finite element method (PFEM). Revista Internacional de M´etodos Num´ericos para C´alculo y Dise˜no en Ingenier´ıa 28(2), 112–123 (2012). [in Spanish] 28. Salazar, F., San-Mauro, J., Celigueta, M.´ A., O˜nate, E.: Air demand estimation in bottom outlets with the particle finite element method. Computational Particle Mechanics 4(3), 345–356 (2017). DOI 10.1007/s40571-016-0117-4. URL https://doi.org/10.1007/s40571-016-0117-4 29. Saunders, K., Prakash, M., Cleary, P.W., Cordell, M.: Application of smoothed particle hydrodynamics for modelling gated spillway flows. Applied Mathematical Modelling 38(17-18), 4308–4322 (2014) 30. Teng, P., Yang, J.: Modeling and prototype testing of flows over flip-bucket aerators. Journal of Hydraulic Engineering 144(12), 04018,069 (2018) 31. USBR: El vado dam – service spillway modification – physical model study. Tech. Rep. HL-2017-02 (2017) 32. Valero, D., Bung, D.B., Crookston, B.M.: Energy dissipation of a type iii basin under design and adverse conditions for stepped and smooth spillways. Journal of Hydraulic Engineering 144(7), 04018,036 (2018) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 Shockwaves in spillways with the particle finite element method. 37 33. Valiani, A., Caleffi, V.: Brief analysis of shallow water equations suitability to numerically simulate supercritical flow in sharp bends. Journal of Hydraulic Engineering 131(10), 912–916 (2005) 34. Vischer, D., Hager, W.H., Cischer, D.: Dam hydraulics, vol. 2. Wiley Chichester, UK (1998) 35. Xue, H., Diao, M., Ma, Q., Sun, H.: Hydraulic characteristics and reduction measure for rooster tails behind spillway piers. Arabian Journal for Science and Engineering 43(10), 5597–5604 (2018) 36. Yang, J., Andreasson, P., Teng, P., Xie, Q.: The past and present of discharge capacity modeling for spillways—a swedish perspective. Fluids 4(1), 10 (2019) 37. Zienkiewicz, O., Ortiz, P.: A split-characteristic based finite element model for the shallow water equations. International Journal for Numerical Methods in Fluids 20(8-9), 1061–1080 (1995) 38. Zienkiewicz, O.C., Taylor, R.L.: The finite element method for solid and structural mechanics. Butterworth-Heinemann (2005) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65