Abstract
This paper deals with the development of a 2D weakly compressible SPH model to simulate wave pressures acting on vertical and slotted coastal structures. Attention is devoted to investigate the diffusive term in the continuity equation in order to smooth out the high-frequency numerical noise in the pressure field. A hybrid formulation based on two literature diffusive models is proposed. The interaction between regular waves with vertical and perforated breakwaters is analyzed in time, space, and frequency domain. Numerical results are compared with laboratory experiments and other diffusive SPH formulations, varying the magnitude of the adopted diffusive term. On the basis of an error analysis, the results show that the hybrid formulation gives a better agreement with the experimental data for the majority of the investigated cases. Moreover, SPH simulations highlight nonlinear trends of dynamic pressures in correspondence with geometrical singularities, such as the holes of slotted walls, due to strong pressure drops induced by the flow motion.
1. Introduction
Smoothed particle hydrodynamics (SPH), introduced in 1977 for simulating astrophysical problems [1, 2], is arguably the most popular mesh-free method. Along the years, SPH was used to study a wide range of complex hydrodynamic phenomena, including flows impacting on structures. SPH technique was firstly developed to model weakly compressible free-surface flows [3] and successively extended, for instance, to unsteady flows (see, e.g., [4, 5]), flows impacting on rectangular and cylindrical bodies [6], advective-diffusion phenomena induced by pollutants in water [7], free-surface open-channel flows [8], processes of hard rock breaking [9], and jet propagating into water domains [10]. In coastal engineering field, SPH simulations were limited to study the interaction between waves and vertical or rubble mound breakwaters (see, e.g., [11–13]). In any case, little attention has been paid to the evaluation of wave pressures for design purposes.
In this paper, coastal structures widely adopted in harbors, such as vertical and perforated-wall caisson breakwaters, are analyzed in terms of pressure loads. The main purpose of the perforated ones is to minimize the wave reflection and guarantee safe navigation during sea storms. The performances of perforated breakwaters are usually analyzed by a structural and a hydraulic point of view [14]. In particular, pressure distributions for vertical and perforated breakwaters have been mainly assessed in literature using empirical and approximated formulas (see, e.g., [15, 16]), while the use of advanced numerical tools has been rarely adopted to support their design. For the optimization of these coastal structures, the correct evaluation of the flow impacts in the fluid-structure interaction is fundamental.
In this context, standard weakly compressible SPH model is employed for the numerical simulations. It is well known that standard SPH leads to the generation of high-frequency numerical noise in the pressure field that can compromise the numerical solution, avoiding correct wave loads analysis. Some authors have recently introduced a diffusive term in the continuity equation to smooth the numerical noise occurring in the pressure field [17–20]. Even if the formulations of Molteni and Colagrossi [17], Ferrari et al. [19], and Groenenboom and Cartwright [20] show an appropriate smoothing level of spurious high-frequency pressures along solid bodies, they give a decay of potential energy for a long simulation time [21]. Conversely, the diffusive model of Antuono et al. [18] is able to preserve the hydrostatic pressure component, but it has some limitations in evaluating the dynamic pressures when flow impacts on a solid boundary occurs. In the case of a breakwater subjected to several wave trains, the flow field is firstly characterized by a slow motion induced by nonbreaking waves along the flume (from the wavemaker to the breakwater) and, successively, by a faster motion due to the wave impact in which the magnitude of the pressure peaks is dependent on the type of incoming wave attack [16]. In order to obtain an improved model for the interaction between waves and breakwater walls involving different flow dynamics, a hybrid diffusive formulation based on the coupling of models by Molteni and Colagrossi [17] and Antuono et al. [18] through a linear transition between them is introduced here.
The proposed hybrid approach is firstly applied to simulate the interaction between regular waves and a vertical breakwater and, successively, for two types of perforated breakwaters. SPH simulations of dynamic pressures are compared with experimental data and the diffusive models by Molteni and Colagrossi [17] and Antuono et al. [18].
2. SPH Model
2.1. Governing Equations
In the case of a viscous, weakly compressible and barotropic flow, the adopted field equations are the Navier-Stokes equations and a state equation that relates the evolution of the pressure field with the density field: where , , , and are, respectively, the position of a generic material point, its velocity, pressure, and density; represents the mass force acting on the fluid, the initial density at the free surface, the initial sound speed, and the viscous stress tensor.
The weakly compressible assumption in SPH implies a very low sound speed with respect to the actual value. In order to reduce the simulation time, the initial sound speed is evaluated through the following constrain: where is the water depth. The characteristic velocity of (2) is set to be equal to in the case of water wave studies, where is the wave celerity. However, this assumption does not influence the fluid dynamics and it results in density fluctuations lower than 1% [22].
The discretized SPH scheme of the system (1) is characterized in the momentum equation by an artificial pressure to reproduce the fluid viscosity and furnish a good numerical stability [22]. A numerical filter called XSPH is added to correct the updating of particle motion [3]. As a consequence, the discrete formulation of the system (1) is where the subindexes indicate the quantities associated with the generic th or th particles, while the symbol denotes the gradient taken with respect to the coordinates of particle . The vector represents the spatial coordinate of the th and th particle. In particular, , (th particle volume) = , where is the th particle mass (constant during the flow evolution). The symbol represents the smoothing or kernel function depending on the value of = , where is the smoothing length. The artificial viscosity coefficient is linked to the dynamic viscosity and it allows simulating inviscid flows. For hydraulic processes, its magnitude assumes values ranging from 0.01 to 0.05 [3]. The quantity is expressed by
The influence of the XSPH numerical filter on the velocity field is quantified by the parameter , ranging from 0 to 1. The kernel function adopted in the SPH scheme is Gaussian with a cutoff radius set here at , assuming that , where is the initial interparticle distance. The diffusive term added in the continuity equation to stabilize the numerical solution will be further analyzed in Section 2.3. The equations expressed by the system (3) preserve the global mass and both the linear and angular momenta.
2.2. Computational Strategies
A stable integration scheme can be adopted for the solution of system (3). In this context, the evaluation of some diffusive terms can be computationally demanding; consequently, the choice of an appropriate integration scheme can reduce the run time of the simulation. Considering the differential relationships of the governing system in the following form: the generic quantity can be split as , in which contains only the diffusive term. Operating this manipulation, the system is solved with a 4th order Runge-Kutta integration scheme through a frozen diffusive approach [23]:
The scheme (6) differs from the original one by taking the diffusive term to be evaluated once within each time step and the numerical errors are assumed to be small [21].
Coupled with the use of a ramp function for the first wave period, the wavemaker is updated at each subtime step of system (6) to give a gradual movement and avoid shock phenomena during the generation phase. The time evolution of the wavemaker displacement, , reads as where is the velocity of the wavemaker.
To model the solid boundaries, fixed ghost particles [6] are implemented. In the fixed ghost particles framework, the ghost particles are fixed in the frame of reference of the solid and are created only once at the beginning of the simulation with a regular distribution. The fixed ghost particles cover a body region with size equal to the width of the kernel support radius. To compute the quantities attributed to each ghost particle an interpolation point is associated with it, obtained by mirroring the position of the fixed ghost particle into the fluid domain. The main advantage of using the fixed ghost particles instead of the standard ghost ones (see, e.g., [24]) is that their distribution is always uniform and does not depend on the fluid particle positions, allowing for a simple modeling of the solid profiles.
At each time step, the physical quantities of the interpolation point are evaluated through a MLS interpolation of the fluid particle values [25]. For example, the pressure of ghost particles, , simulating a solid boundary is calculated as where is the distance between the ghost particle (or the interpolation point inside the fluid) and the body profile and it proves to be dependent on . The kernel is referring to the one with MLS correction. In the simulations, the solid body is approximated by regular equispaced grid of boundary particles with a prescribed distance equal to .
2.3. Diffusive Formulations
Standard weakly compressible SPH formulations have the drawback of generating spurious numerical oscillations in the pressure field. In recent years, different authors have proposed several diffusive corrections to stabilize the solution and attain more reliable results [17–20]. As noticeable in (3), these formulations are obtained through the introduction of a diffusive term in the continuity equation. Groenenboom and Cartwright [20] derived a diffusive formulation by inspecting the time-discretized version of the momentum equations for an inviscid fluid without external forces and analyzing the backward Euler finite difference formulation: where gives the magnitude of the diffusive term in the continuity equation. As already noticed by the authors, this formulation is similar to that proposed by Molteni and Colagrossi [17], if the linearized state equation is assumed. Specifically, a direct connection between the two formulations can be found when the acoustic time step restriction referring to the weak compressibility assumption is adopted: in which , where is the velocity of the particles of the system.
For the considered fluid problems characterized by a wave-structure interaction, the above constraint represents the most restrictive condition ensuring the temporal stability with respect to other conditions linked to the particle acceleration, the artificial viscosity, and the diffusive term. In this way, it is possible to write a relationship between the parameter of Groenenboom and Cartwright model [20] and the parameter of Molteni and Colagrossi model [17]: The parameter is thus a function of the integration scheme since the CFL number is involved in its evaluation. Moreover, in this formulation it is not necessary to impose any other constraints on the time step length related to the diffusive term, since it is intrinsically satisfied by the time step being explicitly present in this diffusive scheme. Another formulation was proposed by Ferrari et al. [19] that presents a slightly different form of the Molteni and Colagrossi [17] formulation but has the same mathematical structure.
In the mentioned formulations, the diffusive term is approximated by the Morris formula that represents the Laplacian in SPH schemes: where is a scalar quantity. To investigate the behavior of this operator, Antuono et al. [18] studied the convergence of the Laplacian in SPH context and they noticed that this last formula is singular near the free surface. In particular, they found that, for kernels in the form , the following expression holds: in which . From (13) the Morris formula (12) converges to the Laplacian of only if and . This condition is verified inside the fluid domain but not near the free surface [26], and diverges as introducing errors in the numerical solution.
To find an expression that always converges in the fluid domain, Antuono et al. [18] proposed a modified diffusive term to better approximate (13). The improvement of the model comes from the introduction of renormalized density gradients, inside the diffusive term, in order to assure the convergence over the fluid domain and preserve the conservation of mass. This formulation, called -SPH, reads as in which
The quantity represents the renormalized density gradient defined as
The range of variation of parameter has been successively found by Antuono et al. [21] through a linear stability analysis and, in case of a Gaussian kernel, it results in
In the typical range of adopted in hydraulic applications, is smaller than about 0.22.
2.4. Hybrid Formulation
Antuono et al. [21] presented a detailed analysis of the diffusive terms, putting particular attention to Morris-like models (12) and to -SPH formulation (14). In particular, because of the inaccuracy near the free surface, the Morris-like models show an unphysical upwards displacement of the fluid particles near the free surface. This inaccuracy is no longer encountered in the -SPH scheme, as a result of the presence of the renormalized density gradients. Another issue related to the former models is the loss of potential energy, also noticed for purely hydrostatic simulations. Conversely, the latter model preserves the hydrostatic solution. A different situation is encountered when high velocities and rapid changes in the fluid domain occur. In these cases, the Antuono et al. [18] model shows the presence of unphysical traveling sound waves after the impact with a solid boundary.
Here, a small modification to the diffusive term proposed by Antuono et al. [18] is introduced by considering a tuning parameter that activates or deactivates the renormalized density gradients as a function of the position of the particles inside the fluid domain. Considering (15), the expression for becomes in which . In the case that , the Molteni and Colagrossi [17] formulation is recovered, while the Antuono et al. [18] formulation is recovered when . For values in the range , a transition (or hybrid zone) between the two formulations is obtained. This last condition is introduced in order to give a gradual transition between the two models, instead of a step function. A variation law for the -parameter is defined for the transition zone, which in this case is considered as simple linear variation, though smoother function can be taken into account. Referring to (13), the introduction of this parameter acts on the second term of the right hand side of the equation; hence, the transition from the involved two models can be regarded as a gradual introduction of a higher approximation of the diffusive term. The use of the hybrid diffusive term implies that in the fluid domain away from the walls of the breakwater the -SPH formulation is applied, while in the area close to the structure the Molteni and Colagrossi term is enforced. In this way, the 4th order diffusive formulation conserves the fluid properties in the majority of the fluid domain and, in presence of fluid impacts, the 2nd order diffusive term smooths out the spurious oscillations, avoiding the generation of spurious shock waves.
The diffusive terms of Molteni and Colagrossi [17], Antuono et al. [18], and the hybrid formulation act in an area fixed heuristically to the length of the adopted support kernel (). The same size is considered for the case of the vertical wall and for the case of the perforated breakwaters. In the latter case, the transition between the two models is introduced also inside the chamber of the caisson where three vertical walls are considered in the analysis of the wave pressures, as displayed in Figure 1.
3. SPH Numerical Simulations
In this section, the capabilities of the proposed hybrid formulation are compared to the diffusive models proposed by Molteni and Colagrossi [17] and Antuono et al. [18], against laboratory experimental data [27, 28]. SPH simulations of wave pressures are performed for the cases of a vertical breakwater and two types of perforated breakwaters subjected to regular waves propagating along a plane channel in intermediate water conditions. In order to optimize the numerical results, the SPH parameters and have been set equal to 0.01 and 0.25, respectively. For the vertical breakwater and the perforated ones, the spatial resolutions are and 115, respectively.
The evaluation of the dynamic pressures acting on the body profiles of the breakwaters is performed considering an interpolation through a kernel function, , corrected with a first order moving least square interpolator, over the fluid particles on a support area with radius equal to . The adopted numerical pressure gauges are located along the vertical walls, adopting the initial spatial discretization .
3.1. Vertical Breakwater
The first test case deals with the interaction between a regular wave train and a simple vertical breakwater. The incident waves are 3rd order Stokes waves, described by a wave height and a wave period . The value of is fixed to 0.54 m. Numerical results deduced from the adopted three diffusive SPH schemes are compared with the experimental data obtained by Kirca and Kabdasli [27] in terms of dynamic pressures, , acting on the vertical body profile and, particularly, for some specific points located below and above the SWL. In the following simulations, the same order of magnitude of the diffusive term () is imposed. Figure 2 shows the time variations of dynamic pressures, respectively, at , 0.93, and 1.07, where corresponds to the SWL. Since nonbreaking waves are taken into account, the resulting shape is characterized by the occurrence of standing wave pressures [16] with temporal scale equal to the wave period. They are defined by a double positive peak and a greater negative peak. For the considered points of the body profile, a good agreement between all diffusive models and the experimental data can be observed in terms of amplitude and wave shifts of . Since the reference points are located close to the SWL, the diffusive model of Molteni and Colagrossi [17] shows a slight progressive decreasing trend of due to decay of potential energy and a progressive underestimation of positive peaks and overestimation of negative ones, as noticeable in Figure 2(a). It is worth noting that the reference experimental data of dynamic pressures exhibit spiky values of positive peaks and nonphysical negative values acting above the SWL (see Figure 2(c)).
(a)
(b)
(c)
In any case, an evident difference in the application of diffusive SPH models refers to the magnitude of high-frequency spurious oscillations. This feature is here exploited through a spectral analysis, allowing the dynamic pressure signals in frequency domain to be analyzed. For the same three locations mentioned above, the dynamic pressure power spectra, , are highlighted in Figure 3, where the vertical black lines are inserted to separate the real energy contents linked to the real wave pressure and the spurious oscillations occurring at a frequency greater than the different wave harmonics. It can be observed that the occurrence of significant high-frequency noise appear in a range greater than the real wave harmonics (in the considered waves, range between 3 and 4, where is here the frequency and is the peak frequency). The greater spurious energetic contents are associated with the model of Molteni and Colagrossi [17]. For the spectrum of referring to the point above the SWL (Figure 3(c)), high-frequency energetic levels are distributed almost uniformly for values greater that the real wave harmonics.
(a)
(b)
(c)
The level of nonphysical high-frequency noise appearing in the oscillatory feature of is analyzed by evaluating the ratio between the zero-order moment of spurious high-frequency noise in the power spectrum of dynamic pressure, , and the zero-order moment of the power spectrum of dynamic pressure, , as follows: where is the number of wave harmonics occurring in the dynamic pressure wave spectrum and and are the minimum and the maximum wave frequencies.
Figure 4 describes the variation of the ratio along the vertical wall (from the bottom to the maximum wave run-up, ). The magnitude of high-frequency wave pressures furnishes the capabilities of the considered diffusive terms to smooth out the numerical noise and recover a well filtered oscillatory signal. For points located far from the SWL () the diffusive model by Molteni and Colagrossi [17] induces lowest values of resulting in a smoothest effect of the dynamic pressure, while the model by Antuono et al. [18] shows highest values of spurious numerical noise. The vertical range corresponds to measuring points located in the underwater part of the breakwater, that is, below the maximum run-down, . The hybrid approach furnishes values of in between the other formulations. Across the SWL () limited by the maximum run-up and run-down, the hybrid formulation and the model by Antuono et al. [18] furnish lower values of high-frequency wave pressures than the model by Molteni and Colagrossi [17]. This last model shows highest energy associated with the numerical noise due to the lack of the renormalized density gradients which conversely act to regularize the free surface and, consequently, the spatial distribution of fluid particles near the gauge point. However, we remark that the above spectral analysis gives an assessment of the diffusive models only in terms of the oscillatory pressure component.
3.2. Perforated Breakwaters
Concerning the modelling of perforated breakwaters, two different geometries are considered [28]. Both caissons present a top cover plate over the free surface and one internal chamber. The front wall of the first breakwater is perforated by two rectangular holes, while the back wall is solid. It is characterized by a porosity = 20%, a width of chamber = 0.2 m, a vertical base = 0.2 m, and a height of top cover above SWL, = 0.08 m. The second one is perforated by three holes, adopting = 20%, = 0.15 m, = 0.2 m, and = 0.16 m. In a plane channel 0.4 m deep, second-order incident waves are characterized by = 0.08 m and = 1.2 s for the breakwater number 1 and by = 0.1 m and = 1 s for the breakwater number 2.
In Figures 5 and 6 are shown the spatial distributions of dynamic pressures at the front wall and the internal walls of the chamber, occurring when the maximum pressure induced by the wave crest within a regular wave train appears in correspondence with the SWL. The geometrical characteristics of the breakwaters (, , and ) and the water depth, , are illustrated in Figures 5 and 6. Positive dynamic pressures are directed in the external side of the walls and negative ones in their internal side. We note that the pressure peaks at the three reference walls appear at different time instants. Their phase shift is dependent on the wave celerity and the width of chamber. SPH results obtained with the hybrid formulation and the models by Molteni and Colagrossi [17] and Antuono et al. [18] using = 0.1 are compared with the experimental data by Chen et al. [28], showing an overall agreement on the front wall and the right internal wall. The magnitude of along the depth is greater for the front wall and lower for the internal walls due the wave dissipation generated by the presence of perforations, as also observed by Takahashi [16]. A general underestimation of the dynamic pressures is observed when the formulation of Molteni and Colagrossi [17] is applied. Since the hydrostatic solution is not preserved in the above mentioned model, unphysical negative values of appear near the lower horizontal parts of the caisson and the maximum dynamic pressure at the SWL is underestimated with respect to the laboratory data. Conversely, the use of the present hybrid formulation and the one developed by Antuono et al. [18] lead to a better assessment of load diagrams along the depth where a general compressive strain state occurs at the walls due to the preservation of the hydrostatic pressure, as expected. A substantial difference between SPH and experiments in evaluating can be observed near the holes of front walls where the results obtained by the three diffusive models highlight lower dynamic pressures than the experimental ones. Near these geometrical singularities, the reference experimental measurements are missing because of constructional constraints induced by the difficulty to place the pressure sensors near the edges of a structure. As solely carried out in the design purposes, a linear variation is imposed between the measurements deduced from a discrete battery of pressure sensors [15] and this approach was historically adopted in laboratory investigations addressed to study the interaction between waves and perforated breakwaters (see, e.g., [14]). The experimental load diagrams adopted for the evaluation of wave forces are therefore approximated with respect to the spatial discretization but, at the same time, they provide a safe dimensioning from the engineering point of view. The discrepancies between SPH and experimental values allow furnishing similar values of the resulting wave forces acting on the large inertia caisson and different results when the total force is applied to the slender structural components of the slotted wall. This result suggests the need to assess load diagrams on distinct elements of the perforated breakwaters for subsequent structural analyses.
In the reference experiments by Chen et al. [28], the pressure gauges were only placed along the solid contours of the vertical parts of the considered perforated breakwaters and, consequently, the experimental values of pressures were adopted to evaluate horizontal wave forces. As a result, the validation of the proposed hybrid SPH model refers to the available experimental values due to the lack of pressure measurements at the top cover of the caissons. Moreover, we observe that, in the present numerical simulations, the adopted wave trains interact with a small part of the top cover and for a short time interval just for the perforated breakwater number 1.
An application example of the adopted hybrid diffusive formulation is given in Figure 7, where the dynamic pressure distribution at the perforated breakwater number 1 is illustrated when the maximum wave crest impacts on the front wall. Strong variations in the values of are evident across the holes of the structures due to the loss of wave energy. In this case, negative dynamic pressure can occur during the passage of the water flux. The pressure drops are related to the strong horizontal velocity gradients occurring at the geometrical singularities, as expressed by the canonical Bernoulli’s theorem (see, e.g., [14]).
In order to assess the capabilities of the involved diffusive SPH forms in simulating the experimental values of wave pressures on the three faces of the caissons [28], Figure 8 shows the root mean square errors, , as a function of the magnitude of the diffusive parameter . In all cases the value of is less than 30 N/m2 and, particularly, for the second perforated breakwater, the use of the hybrid formulation gives the lowest errors.
(a) Perforated breakwater number 1
(b) Perforated breakwater number 2
4. Conclusions
The work presented an SPH analysis of dynamic pressures induced by the interaction between regular wave trains and vertical/perforated breakwaters. The research was addressed to overcome the limitations in the breakwater design, in particular, when different structural components are involved. The present phenomenon is characterized by the occurrence of different flow dynamics and the pressure field was evaluated through a hybrid formulation based on the coupling of the diffusive terms developed by Molteni and Colagrossi [17] (close to the breakwater walls) and Antuono et al. [18] (away from them) through a transition zone.
Numerical simulations of the wave pressures acting on a simple vertical breakwater and two types of perforated breakwaters using the adopted diffusive SPH schemes were compared with laboratory experiments. SPH analysis in time, frequency, and space domain allowed an assessment of the capabilities of the involved diffusive models to smooth out the high-frequency noise and furnish a suitable evaluation of dynamic loads. Application of Molteni and Colagrossi model provided the smoothest effect on the spurious frequencies in the underwater part of the breakwater but it failed on the estimation of the maximum wave loads close to the SWL and the horizontal structural parts. The present hybrid formulation gave less errors and a higher smoothing effect than the Antuono et al. model along the depth. However, these two models furnish a globally better evaluation of the load diagrams and, in particular, between the wave run-up and run-down where the highest wave pressures occur. With regard to the perforated structures, SPH simulations described the presence of pressure drops occurring across their holes. Further investigations will be addressed to analyze the performances of the proposed hybrid formulation considering breakwaters subjected to breaking wave loads.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work has been partially supported by European Commission, European Social Fund 2007/2013, and Calabria Region. The authors would like to thank V. S. O. Kirca (Istanbul Technical University) for providing some experimental data and A. Colagrossi and S. Marrone (CNR-INSEAN) for the fruitful discussions. Finally, the authors thank the anonymous reviewers for the useful comments.