Abstract

Breaking wave-induced scour is recognized as one of the major causes of coastal erosion and offshore structure failure, which involves in the full 3D water-air-sand interaction, raising a great challenge for the numerical simulation. To better understand this process, a nonlinear 3D numerical model based on the open-source CFD platform OpenFOAM® was self-developed in this study. The Navier–Stokes equations were used to compute the two-phase incompressible flow, combining with the finite volume method (FVM) to discretize calculation domain, a modified VOF method to track the free surface, and a model to closure the turbulence. The nearshore sediment transport process is reproduced in view of shear stress, suspended load, and bed load, in which the terms of shear stress and suspended load were updated by introducing volume fraction. The seabed morphology is updated based on Exner equation and implemented by dynamic mesh technique. The mass conservative sand slide algorithm was employed to avoid the incredible vary of the bed mesh. Importantly, a two-way coupling method connecting the hydrodynamic module with the beach morphodynamic module is implemented at each computation step to ensure the fluid-sediment interaction. The capabilities of this model were calibrated by laboratory data from some published references, and the advantages/disadvantages, as well as proper recommendations, were introduced. Finally, nonbreaking- and breaking wave-induced scour around the monopile, as well as breaking wave-induced beach evolution, were reproduced and discussed. This study would be significantly helpful to understand and evaluate the nearshore sediment transport.

1. Introduction

As ocean wave gets closer to the coastal region, the free surface deforms and wave steepens with increasing wave height and may induce unstable and break waves. Wave breaking and the resulting current are the causes of many complex phenomena in the surf zone, such as the coastland drift, beach erosion, sediment transport, and so on. They will erode the soil-supporting pile foundations [1], which were constructed within aforementioned coastal high-hazard areas for practical reasons [2]. Therefore, assessing the hydrodynamic and sediment transport processes are important for inferring the strength of potential destruction during the extreme environment (e.g., Wu et al. [3, 4]).

Over the last few decades, some research studies were performed to improve the understanding of above processes. Field surveys were performed by analyzing sediment diameter distribution, instrument data, and field videos (e.g., Paris et al. [5]; Szczuciński et al. [6]; FEMA [2]; and USGS [7]), which provided some substantial cognitions, such as the majority of the deposits are from beach and not from deep ocean floor (e.g., Morton et al. [1] and Jaffe et al. [8]). However, the field research studies raise considerably limited information due to frequent lack of the wave condition and topography, especially facing the complex three-dimension high-speed water roller. Additionally, theoretical approaches for studying this process are still inadequate [9]. Hence, numerical simulations and physical experiments, providing effective controlling factors, become the preferred methods.

Traditionally, physical experiments are clearly helpful to assess this process. Many research studies have been performed in the aspect of cross-shore sediment transport on sandy beach (e.g., Kobayashi and Lawrence [10]; Young et al. [11]; Jiang et al. [12]; and Daghighi et al. [13]) and the scour around offshore structure foundation (e.g., Kato et al. [14, 15]; Tonkin et al. [16]; and Kuswandi et al. [17]), which improve our understanding of beach evolution and its response to structure by wave. With the development of computer technology, the numerical model provided us with another way to understand the mechanisms, which solves the limited capability used in measurement devices. In the past decades, modeling sandy beach evolution has been studied with depth-averaged equations, such as shallow water equations (e.g., Simpson and Castelltort [18] and Pritchard and Dickinson [19]) and Boussinesq-type equations (e.g., Shimozono et al. [20] and Xiao et al. [21]). However, this process is highly nonlinear, as well as local but strong turbulence near the free surface, and the seabed needs to be considered. Recently, Nakamura and Yim [22], Li et al. [23], and Jacobsen et al. ([2426]), as well as Liu et al. [27], introduced a volume of fluid-type model to study this process, which was based on the generalized Navier–Stokes equations with a turbulence closure solver computing incompressible viscous multiphase flows. Aforementioned research studies not only substantially improved our understanding of breaking wave-induced longshore sediment transport mechanisms but also gave us some inspiration to explore the more complex problems, such as sediment-wave-structure interaction in some coastal areas. Pan and Huang [28] are the pioneers to numerically investigate the interactions among wave, monopile, and sandy beach to access the effectiveness of a linked 2D hydrodynamic and sediment scour model in this process modeling. However, some limitations existed, such as the 2D model could only simulate vertically averaged currents instead of 3D currents in this study, in particular about the vertical velocity component [28]. Huang et al. [29] pointed out that the approximation of the 3D scour phenomenon around the pile structure to a 2D problem may be a contribution to errors between the model predictions and observations. To the authors’ knowledge, the research of 3D modeling in this question is still limited so far and yet to be further investigated.

To address the aforementioned issues, a 3D numerical model was developed based on the open-source OpenFOAM platform in this study, supporting two-phased incompressible flow. The wave dynamic is achieved by its solver interFoam combined with the active wave generation and absorption boundary conditions, developed by Higuera et al. [30]. Unfortunately, the official version of OpenFOAM lacks the function of sediment transport simulation, and the third-party sediment module could not be open accessed. Additionally, similar sediment solving modules did not necessarily apply to current research cases due to the complexities of sediment transport mechanisms under different sea environments. Therefore, this study aimed to overcome those gaps by self-extending the sediment functionality based on an effective usage of the better parts from published sediment modules. This developed model will be further used to simulate the breaking wave-induced scour around the monopile, as well as investigate its difference with the scenario of monopile on a flat seabed and the potential scour mechanisms. The generation, advection, and dissipation of turbulence in the fluid model are also important for sediment transport. For breaking wave-driven sediment motion, several turbulence models have been employed, namely, the RANS model (Li et al. [23] and Jacobsen [25, 26]) and the LES model (Christensen and Deigaard [31] and Christensen [32]). Generally, LES was under consideration as the turbulence model because of a considerable part of the mixing is resolved. However, when using the LES, the stochastic nature of the flow reflects nonlinearly onto the sediment transport and results in larger spatial gradients; hence, the allowable calculation stability or time step will be lowered considerably. Thus, the RANS turbulence model was considered in this study because it is more preferable to predict the scour profile than the LES turbulence model [33, 34].

The rest of this paper is organized as follows. The numerical methods including the hydrodynamic model and beach morphodynamic model, as well as corresponding boundary conditions and numerical schemes are introduced in Sections 24, respectively. The validations of the main modules are conducted in Section 5. The model applications and discussions are given in Section 6. The main conclusions of this study are shown in Section 7.

2. Hydrodynamic Model

3D RANS equations were applied to simulate the two-phase incompressible viscous fluids, and the governing equations are as follows:Continuity equation:Momentum equation:where is the Cartesian coordinate system, is the velocity vector, ρ is the fluid density, is the pressure in excess of hydrostatic, and g is the gravitational acceleration.

The last term of equation (2) is used to describe the surface tension, in which is the surface tension coefficient, is the surface curvature, and is a scale field used to track the fluid movement. is the Reynolds stress tensor given bywhere is the eddy viscosity, is the rate of strain tensor, k is the turbulent kinetic energy, and is the Kronecker delta function confirmed as 1 if i equal to j and as 0 if i unequal to j.

The second order standard turbulent model was employed in this study to closure the set of RANS equations:where is the dimensionless coefficient, l is the turbulence length scale, and is the turbulent dissipation rate.

The transport equations of k and are as follows:where is the turbulence kinematic viscosity and , , , and are the empirical coefficients confirmed as 1.44, 1.92, 1, and 1.3, respectively.

The free-surface motions were captured by the modified VOF method and defined aswhere the air and water are described by volume fraction (α) ranging from 0 (air) to 1 (water). Compared with the traditional VOF method, the extra compression term was introduced in the last term of equation (6) to limit the excessive numerical diffusion and the smearing of the interface, where ur is the relative velocity.

3. Beach Morphodynamic Model

In this study, the details on each component of the beach morphodynamic module are selected and tested according to a sequence of classic examples, and meanwhile, some proper parameters are recommended, which is a comprehensive usage of some previous research studies.

The bed shear stress, linking the hydrodynamic features and sediment transport, is estimated by the method proposed by Arzani et al. [35]:where is the wall traction confirmed as , is the unit normal vector that is perpendicular to the bed surface mesh, and is the stress tensor given bywhere the first term is the pressure-induced stress term, is the pressure, and the second term is the viscous stress term controlled by the bottom fluid motion.

The volume fraction was introduced into equation (7) to identify the bed shear stress of the beach surface.

The bed load transport was described by the formula proposed by Engelund and Fredsøe [36]:where qb is the bed load transport rate per unit width, is the median sediment diameter, R is the relative density of the sediment, θ is the Shields number, and θc is the critical Shields number considering the effect of seabed slope (Allen [37]) and defined as

Equation (11) indicates that θc is adjusted to a lower value for the sediment moving down the slope (i.e., negative slope) and a higher value for the sediment moving up the slope (i.e., positive angle), and is the sediment repose angle. is threshold shields number under flat bed, which could be calculated by [38]where is the dimensionless sediment size.

The bed load transport rates in different directions arewhere is the bed elevation and C is used to reveal the effect of bed slope on the sediment flux and is specified as 1.5 in this study.

The suspended load transport can be described according to the classical convection-diffusion equation as follows:where c is the concentration of the suspended sediment, is the sediment diffusivity equaling to the turbulence eddy viscosity, is the turbulent Schmidt number, which is related to eddy viscosity and sediment diffusivity and is taken to be 0.8, and is the sediment setting velocity, which will be affected by the suspension concentrationwhere ζ is the suspended sediment size related constant and is specified to be 5.0 in the present research according to the study of Liang et al. [33] and is the settling velocity in clear water given by [39]

As far as we know, the coastal zone is the interface area of water, air, and sediment, resulting in a frequent interaction. To keep the sediment to drop out immediately when accidentally left in the air, the volume fraction was introduced in equation (15) as follows:where the velocity is multiplied by α, and the inclusion of in equation (18) is used for numerical stability. This stability problem occurs when becomes 0; hence, the presence of ensures that an advection-diffusion problem is solved locally rather than an advection problem, where the latter is more difficult to solve numerically [39].

According to sediment continuity, the Exner equation was used to describe the seabed evolution:where is the seabed elevation, n is the beach porosity, and qb is the bed load transport rate, and its components are given by equation (14). D is deposition rate solved bywhere is the sediment concentration at bed defined as the value at the nearest cell center of beach in this study.

Additionally, the sediment concentration at a reference level is calculated by the following formula:where T is the dimensionless excess shear stress calculated byin which is the critical shear stress and is the effective shear stress coefficient. The details could be found in van Rijn [40].

E is the entrainment rate, which could be calculated by the following formula:

Finally, the model coupling between water and sediment is achieved by moving the computational mesh in such a way that the bottom mesh is conformal with the seabed, which is mainly based on the method proposed by Jasak and Tukovic [41]. Meanwhile, to ensure that the bed slope does not exceed the sediment repose angle, the mass conservative sand slide algorithm given by Khosronejad et al. [42] was employed:where and are the bed elevations at point and its ith neighbor (see Figure 1), is the horizontal distance between the two cell centers, and and are the corresponding corrections imposed to satisfy the sediment repose angle, which could be obtained by balancing the cell mass as follows:where and are the projections of the cell and its ith neighbor.

4. Boundary Conditions and Numerical Schemes

4.1. Boundary Conditions

Proper boundary conditions are the key factors to copy laboratory experiments. In this study, the active wave generation and the absorption boundary developed by Higuera et al. [30] were employed. Compared with the traditional method proposed by Jacobsen et al. [24], this method does not need the relaxation zone so that it could significantly reduce the computational domain.

The other boundary conditions are set as follows. The top boundary is described as free to the atmosphere, and the seabed and the pile surface are set as the no-slip condition. The symmetric condition is applied at the two-side faces of the numerical tank to reduce the computational cost, in view of the symmetric distribution in the streamwise direction.

4.2. Numerical Schemes

In this study, the finite volume method (FVM) and Euler scheme are used to address the space and time term, respectively. The PIMPLE algorithm is employed for the pressure-velocity solver, which is a mixture between PISO (pressure implicit with the splitting of operators) and the SIMPLE (semi-implicit method for pressure-linked equations) method. The MULES (multi-dimensional universal limiter for explicit solution) scheme is applied to maintain the boundedness of volume fraction. The Gauss linear corrected scheme and Gauss linear scheme are used to solve the Laplacian term and gradient term, respectively. The sediment module is addressed as follows. The bed load is solved firstly by shear stress that is derived from the hydrodynamic model. Meanwhile, the suspended load is transported as the fluid motion based on the convection-diffusion equation, and this equation is also solved by FVM. Subsequently, the two-dimensional Exner equation is applied to update bed elevation, along with the finite area method (FAM) addressing the problem of fluid information mapped from the three-dimensional space to the two-dimensional space (Tuković [43]). The FAM is a variant of the FVM, operating on two-dimensional curved surfaces in the three-dimensional space, which makes the seabed morphology suitable for arbitrary curved surfaces. In the layout, the FAM follows the structure of the Finite Volume discretization library, sharing the basic field algebra, linear system support, boundary condition settings, and discretization techniques. Subsequently, the sand slide procedure is implemented to modify the problem of bed slopes exceeding the sediment repose angle.

The self-developed model in this study is a two-way coupling scheme between water motion and sediment transport, and the computational procedure is implemented in Figure 2. Firstly, the hydrodynamic parameters, such as , are computed using the hydrodynamic module. Second, above parameters are used to drive suspended load and bed load transport, resulting in the seabed change based on the Exner equation and the sand slide equation at the same th time step, i.e., . Third, the hydrodynamic parameters are also affected by the resulting seabed change. Fourthly, the values of hydrodynamic at the next nth time step are calculated by the hydrodynamic module from the previous results based on the updated bed elevation. Fifthly, the sediment transport and seabed change at the nth time step are performed using not only the previous sediment environment at th time step but also the hydrodynamic results at the nth time step. Finally, this process is repeated until the calculated time reaches the preassigned time. During the computation, the time step was automatically adjusted to ensure the Courant number (, in which and were the maximum velocity and the minimum grid size, respectively, and is the time step) is always less than one.

5. Validations of the Main Modules

5.1. Wave Module

Reliable hydrodynamic is the foundation to ensure the accuracy of sediment results. Therefore, the stability of wave propagation over a flat tank was first tested. The expression of the solitary wave proposed by Lee et al. [44] was used to reproduce the free surface:where is the free surface elevation, is the wave height, is the water depth, , and is the wave celerity. Based on one typical solitary wave condition ( = 0.3), seven wave gauges with a constant interval 0.48 m were placed along the tank to obtain the wave surface elevation, and the results are shown in Figure 3. The good agreement between numerical results and theoretical results was found. In addition, the wave surface was almost no attenuation in the numerical domain, indicating the numerical model could better model the wave dispersion and nonlinearity and could be employed in the following research.

5.2. Sand Slide Module

In this section, a constant beach slope of 45° with sediment repose angle of 30° was used to validate the robustness of the sand slide module. The beach morphology before and after employing the sand slide model is shown in Figure 4, respectively. We could find that the beach slope is less than or equal to the sediment repose angle, which indicates the sand slide module can work well once the beach slopes exceed the sediment repose angle.

5.3. Suspended Sediment Module

The transport capacity of the suspended sediment was tested by the net entrainment experiment of van Rijn [40]. The sand bed was laid after a rigid bed in this experiment, in which x = 0 is the starting position of the sand bed, and the inlet of the model is free of sediment (see Figure 5). The flow is stable with a water depth of 0.25 m and a velocity of 0.67 m/s, and the representative sediment diameter is 0.2 mm with the corresponding settling velocity of 0.022 m/s. The roughness height ks is taken to be 0.01 m, suggested in accordance with the research of van Rijn [45] and Wu et al. [46].

The concentration profiles at x = 4H, 10H, 20H, and 40H were chosen to validate the simulated results. In this model, the entire computational domain was discretized using a structured mesh, and the gird sizes of 0.002 m were kept constant in the x-direction, y-direction, and z-direction. Figure 6 shows that the calculated results match well with the experimental data, and the sediment concentration increases with the decreasing height above the bed surface.

5.4. Bed Morphodynamic Module

The experiment of Mao [47] was employed to validate the capacity of sediment transport with the morphodynamic module. The submarine pipeline with diameter of D = 0.1 m was laid on the sand bed, water depth is h = 0.35 m, sediment diameter is d50 = 0.00036 m, sediment density is ρsed = 2650 kg/m3, flow velocity is u = 0.35 m/s, and critical shields parameter is  = 0.048. To avoid the divergence of numerical results at the junction of the seabed and the pipeline caused by the severe deformation of grids, a sinusoidal hole with 0.1 D was used under the pipeline inspired by the research studies, such as Liang et al. [33], Brørs [48], and David et al. [49]. The mesh resolutions around pipelines were increased to solve the problem of mesh distortion as the scour hole broadens and deepens, as shown in Figure 7.

The scour profiles at 10 min and 30 min were chosen to verify the seabed scour morphologies, as shown in Figure 8. Overall, the good agreement between the simulation and experiment was found, in which the sediment accretion in the former behind the pipe was slightly overpredicted than that in the latter. This is mainly attributed to the fact that the RANS model tends to smooth out the fluctuations produced by the vortex shedding and therefore underestimates the interaction between the vortices shed with the seabed, and the similar conclusion was raised by Liang et al. [33]. Nevertheless, the sediment module established in this study could well reproduce the morphological change of the seabed.

6. Model Applications and Discussions

6.1. Nonbreaking Wave-Induced Scour around Monopile
6.1.1. Model Setup

The water-sediment transport capacity of the model had been verified well by the above calculation settings. Furthermore, this model would be used to show its robustness of copying the scouring process around a pile. The numerical setup for the model validation is the same as the experiment by Sumer et al. [50]. The flume is 28.0 m long, 2.0 m wide, and 1.0 m high, and the pile diameter D = 0.10 m fixed in the center of the sand bed, as shown in Figure 9. The sediment diameter is d50 = 0.00018 m, sediment density is ρsed = 2700 kg/m3, critical shields parameter is  = 0.047, and water depth is h = 0.40 m. A periodic wave with wave period of T = 4.5 s and wave height of H = 0.12 m was generated to propagate over the flat sand bed. Figure 9(c) shows the model grids, and the grids around the pile were refined.

6.1.2. Verification and Discussion

Figure 10 shows the computed and theoretical wave surface elevation proposed by Svendsen et al. [51], as well as the computed and measured dimensionless maximum scour depth (S/D) change around the pile, in which S was presented without distinguishing between specific locations such as behind or beside the pile. The simulated wave crests and troughs are almost equal and in of the phase compared to the theoretical data. Meanwhile, it could be seen that the scour depth developed with the periodic wave action is well captured. Nevertheless, a small difference between the simulated and the experimental scour topography could be found. The measured scour depth fluctuated frequently and even exhibited an obvious disparity at two adjacent moments, but the variation of simulated results is rather regular, which should be attributed to the turbulence model based on a RANS method used in this study that tends to smooth out the fluctuations produced by the vortex shedding.

Figure 11 shows the scouring process, and it is seen that the seabed at the upstream side of the pile was eroded after the wave impacting the pile, which is related to the strengthened circling flow. Then, the suspended sediment was deposited downstream. As the wave continue impacting the pile, the scour depth at the upstream side continuously increased and the sedimentation at the downstream side also exhibited the same tendency. Finally, the simulated maximum scour depth took place at the upstream side of the pile, which follows the popular rule. Therefore, the current model is capable of modeling the interaction between the periodic wave and the pile.

6.2. Breaking Wave-Induced Sandy Beach Evolution
6.2.1. Model Setup

The water-sediment transport capacity of the model had been verified well by the above calculation settings. In this section, the numerical evolution of sandy beach was calibrated by the experimental data conducted by Kobayashi and Lawrence [10]. The experiment setup is shown in Figures 12(a) and 12(b), with the beach slope of 1 : 12, sediment diameter of 0.18 mm, porosity of 0.4, and water depth of 0.8 m. In this experiment, a solitary wave with wave height of 0.216 m was generated to propagate over sandy beach. Another solitary wave with the identical height was generated on the evolved beach after this wave subsided, and that cycle repeats. Finally, the experimental results of beach profile, wave surface, and velocity at the end of the fourth wave were selected to test the robustness of the model. Figure 12(c) shows the model grids for sandy beach, and the grid near the seabed was refined.

6.2.2. Verification and Discussion

Figure 13 shows the comparison of wave surface elevations from the numerical model and physical experiment. Good agreements between them could be found, not only the free-surface elevations but also the arrival times.

Figures 14(a) and 14(b) show the evolved beach profile and profile change. We could found that the numerical results in the offshore zone agree well with the measured data. However, the beach profile near the coastline could not be predicted accurately, which may be due to an overestimation of the turbulent dissipation because of the underestimation of the backwash velocity. Meanwhile, the time series of streamwise velocity (u) at the ADV location is shown in Figure 14(c). The numerical model satisfactorily reproduced the streamwise velocity, despite there are still some issues to be solved in underestimating the backwash velocity and the associated beach profile evolution.

After completing the model validation, the wave field, suspended load transport, and the associated beach profile are shown and discussed in Figure 15. Figures 15(a)15(c) show that the ocean wave got steepened and skewed with close to the coastal zone, and then the wave broke as a plunging breaker due to wave shoaling. The breaking waves inject nearshore immediately forming a complex wave field entraining air on the free surface, picking up amounts of sediment from sandy beach into suspension. Subsequently, the suspended sediment further runup, along with the breaking wave-induced current. The beach profile have not obviously changed in this stage, which should be attributed to the upwash direction opposite with the sediment gravity and bottom friction, which can also be found in the laboratory by Young et al. [11].

After upwash reaching the maximum height along the beach, it will start to reverse, as shown in Figures 15(d)15(f). A hydraulic jump is formed when the retreating water tongue collided with the relatively still massive body of water, but there is no obvious sediment resuspension. As upstream retreating water continues to inject into the hydraulic jump zone, the still water body turns to a complex water-air mixing zone that led to a turbulent roll, forming a recirculation region. Meanwhile, an intensive sediment transport in the form of sheet flow was found on the landside of the hydraulic jump point, resulting in net erosion of the onshore region. In contrast, a mass of sediments were deposited on the seaside of the hydraulic jump point due to the sudden deceleration of sediment-rich seaward flow, as well as the long particle residence time induced by the recirculation region. Consequently, the bed morphology of net erosion in the nearshore zone and net deposition in the offshore zone were formed, which was consistent with the measurements by Kobayashi and Lawrence [10] and Young et al. [11].

6.3. Breaking Wave-Induced Scour around the Monopile
6.3.1. Model Setup

Considering a more complex condition, the response of beach evolution to the structure would be discussed. The experiment of local scour around a pile on sandy beach, given by Tonkin et al. [16], was applied to validate the model robustness. In this experiment, the water depth was 2.45 m, wave height was 0.22 m, beach slope was 1 : 20, sediment diameter was 0.35 mm, and sediment density was 2643 kg/m3. A single pile with the diameter of 0.5 m was placed upright on the location of still water level as shown in Figures 16(a) and 16(b). The wave gauges were applied to measure water surface elevation at the locations of stations A, C, and D. An electromagnetic flowmeter, placing at station B, was employed to measure the streamwise velocity. The scour depths at the front, side, and back of the pile, i.e., stations C, D, and E, were measured using digital camera technology. More details could be found in Tonkin et al. [16]. The numerical model was setup in the same dimension with this experiment, and the mesh is shown in Figure 16(c) with the refined grids around the pile.

6.3.2. Verification and Discussion

Figure 17 shows the time series of free-surface elevations at stations A, C, and D and the time series of streamwise velocity at station B from numerical results and measured data. Overall, the predicted results were in good agreement with the measured data, with the exception of the velocity since t = 13 s after the initial wave impact. Numerical and measured differences of Figure 17(b) after t = 13 s should be attributed to during the most turbulent part of the drawdown, at which the water level is too low to give a valid velocity measurement; the sensor is not fully submerged at this point (Tonkin et al. [16]). Fortunately, the current model is better to make up this deficiency and shows the subsequent velocity process until about t = 20 s. Nevertheless, due to the further reduction of water level, the time series of velocity start to exhibit the intense numerical oscillation.

Figure 18 shows the time series of the local scour depth around a single pile. We could find that the seabed morphology at the front and side of the pile (i.e., stations C and E) was sharply eroded of approximately 10 cm when the moment of wave impacting the pile, which is related to the stronger velocity (see Figure 17(b)). Then, the scour depth decreased as the net deposition occurred during the drawdown. The seabed morphology behind the pile (i.e., station D) did not exhibit an obvious change at the moment of the wave impacting the pile, and then the scour depth continuously increased from t = 12 s to t = 18 s. During the latter stage of the backwash, the scour depth decreased to some extent due to the sediment settlement. After about t = 20 s, the seabed around the pile was relatively in equilibrium. Finally, we noted that the predicted maximum scour depth is slight lesser than the measured data, which may be caused by seepage flow (or groundwater flow) as described in the Tonkin et al. [16], but the current model was incapable of reproducing it.

To more visually present the topography evolution of sandy beach during the interaction between the wave and the monopile, the three-dimensional scour patterns are shown in Figure 19. It is obvious that the remarkable erosion took place at the front of the pile when the wave impacting the pile, which is different with the aforementioned scenario of the pile placed on a horizontal seabed. As upstream retreating water, the scour depth on the seaside of the pile increased quickly. The vorticity fields, which are related to the sediment transport (Kobayashi and Oda [52]), were chosen to reveal the potential causes of above phenomena, as shown in Figure 20. At the runup stages, the flow separated and the vortex motion formed and developed mainly located at the area less than 0.5D behind the pile. Following the wave drawdown, the wake region appeared before the pile. The generated vortex was transported further upstream by advection, and the wake region was stretched with a length scale more than 2D. This may be related to the backwash velocity larger than the uprush velocity, which is showed in Section 6.2. Overall, it is acceptable that backwash-induced scour around a pile is larger than the uprush.

7. Conclusions

A nonlinear three-dimensional coupled model, which was composed of the hydrodynamic module and the self-developed beach morphodynamic module, was developed based on the CFD tool OpenFOAM in this study. To better simulate the breaking wave-induced seabed scour around the monopile, an innovative combination of sediment components was proposed which is different with published papers, and its robustness was verified by a series of tests. The results showed that this self-developed model was able to satisfactorily reproduce the wave transmission, suspended load and bed load motion, and the seabed morphology evolution, as well as proved the sand slide module works well.

The 3D model developed here was first applied to explore two preliminary scenarios, one with the nonbreaking wave-induced scour around the monopile, and the other with breaking wave-induced sandy beach evolution. In both cases, the calculated results were basically in agreement with the measurements and meanwhile showed us that the obvious scour occurred in the upstream side of the pile under the nonbreaking wave, as well as the backrush-induced beach deformation is larger than the uprush. Subsequently, the more complicated case of wave-induced scour around the monopile on sandy beach was simulated, and its main features were well reproduced by this model. The remarkable erosion took place at the front of the pile when the wave impacting on the pile during the wave runup, which is different with the scenario of the monopile placed on a horizontal seabed. Following the wave drawdown, scour depth at the upstream side of the pile increased quickly, and this should be attributed to the high-speed backwash and strengthen the wake region. Finally, we remark that further research studies are needed to investigate the effect of the offshore structures on beach morphological evolution based on the current study, which would be significantly helpful to evaluate the nearshore sediment transport.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This study was supported financially by the National Natural Science Foundation of China (grant no. 51779280), the Open Foundation of Key Laboratory of Water-Sediment Sciences and Water Disaster Prevention of Hunan Province (grant nos. 2020SS02 and 2020SS06), and the Guangdong Science and Technology Plan Program (grant no. 2013B020200008).