Viscous Flow Fields Induced by the Run-Up of Periodic Waves on Vertical and Sloping Seawalls

Abstract

Unsteady two-dimensional Reynolds Averaged Navier-Stokes (RANS) equations coupled with the $k-\epsilon$ turbulence model was developed to simulate the viscous flow field near seawalls. The complex free-surface configuration was captured using the particle level set method (PLSM). An innovative solid-fluid coupling method was employed to mimic the solid-fluid interaction on fixed Cartesian grids. The proposed numerical model was applied to investigate the viscous flow fields induced by periodic waves propagating on three different seawalls, including vertical and steep seawalls, which are commonly designed in Taiwan with slopes of 1:2 and 1:5. The results reveal that the wave run-up height is closely related to Iribarren number and breaking wave. Partial standing wave or asymmetric recirculating cells in front of a steep seawall lose symmetry with increasing slope. The wave pressure on the seawall face increases as the Ursell number or seawall slope increases. The undertow caused by breaking waves has a significant effect on the bottom shear stress in front of the seawall.

1. Introduction

Seawall structures built in coastal regions are widely used to protect areas of human habitation, conservation, and leisure activities from the action of waves, tides, and currents. Wave propagation on a sloping seawall (seabed) often leads to wave run-up, run-down, breaking, and overtopping. The reflected waves and incident waves together may form standing waves in front of the seawall. Severe scouring, which is related to the (partial) standing wave and breaking waves, may cause damage to the seawall and reduce its preventive ability. Therefore, wave run-up, wave pressure on seawalls, and toe scouring are important issues concerning the stability of seawalls.

Wave run-up and run-down on seawalls have been commonly studied in the past, being affected by many parameters, such as the incident wave condition, seawall slope, and roughness. Early studies focused on the wave run-up of non-breaking solitary waves due to the simpler wave and flow fields. Ref. [1] conducted experiments to investigate the wave run-up of solitary waves on impermeable slopes. Empirical formulas for maximum wave run-up heights as a function of beach slope and incident wave height were presented on the basis of the laboratory tests. Ref. [2] experimentally studied the evolution of a breaking solitary wave on a 1:60 sloping bottom. They proposed a simple formula for predicting the maximum run-up height of a breaking solitary wave on a uniform beach with a wide beach slope range (1:15–1:60). Wave run-up and run-down are accompanied by complex changes in the level of free liquid, viscous fluid flow, bed shear effects, and breaking wave bubbles. If relying solely on theoretical solutions, it is not easy to accurately describe all physical phenomena.

Laboratory experiments with high-precision equipment were also adopted to investigate the characteristics of wave and flow fields induced by waves propagating over different slopes. Ref. [3] performed experiments for wave shoaling on a 1:20 slope. A Laser Doppler Anemometer (LDA) was used to investigate the structure of the turbulent flow induced by breaking waves. Observation of large oblique vortex structures with descending trajectories were seen to play a dominant role in the mass transportation. Refs. [4,5,6] used the LDA to study the characteristics of flow fields under different breaker types. The undertow was observed in the surf zone under the spilling and plunging breakers. Kinetic turbulent transportation in the surf zone was also studied in the experiments. The experimental results showed that fluctuating kinetic energy is transported shoreward by plunging breakers and seaward by spilling breakers. Ref. [7] employed Laser Doppler Velocimetry (LDV) to measure the flow fields in front of a caisson breakwater with wave-dissipating blocks. This study demonstrated that turbulent is generated not only by the wave-breaking process but also by the flow through the porous armor layer on the caisson breakwater. Using Particle Image Velocimetry (PIV), Refs. [8,9] experimentally study on the turbulence and wave energy dissipation of spilling breakers in a surf zone. The measurements showed that a notable number of intermittent turbulent eddies penetrated into the bottom boundary layer. Ref. [10] conducted both laboratory experiments and numerical simulation to investigate the wave run-up and run-down processes in detail. Laboratory experiments were performed to measure the flow velocity fields by using high-speed PIV (HSPIV). The results were elaborated for providing insights into swash flow dynamics, generated by a non-breaking solitary wave on a steep slope.

Some researchers studied run-up by solving the non-linear shallow-water equations. Ref. [11] determined the maximum run-up height of solitary waves by solving linear shallow-water equations. The comparison of the analytical solution and the experimental data [12] showed fairly good agreement, and so then the regression law for solitary wave run-up on a slope was proposed. Some researchers studied run-up by solving the non-linear shallow-water equations. Ref. [13] solved the Boussinesq equations by incorporating the eddy viscosity and bottom friction to determine the maximum run-up height. Comparisons of the solution with experimental data show good agreement.

Field measurements of wave run-up on seawalls have also been conducted. Ref. [14] investigated wave run-up on a seawall slope in the field and compared the test data with those obtained using an empirical formula. Ref. [15] gave an overview of 20 years (1994–2013) of field measurements at the Petten site in the Netherlands and cataloged the main research findings obtained from these measurements. Ref. [16] conducted physical model tests to investigate wave run-up on berm coastal structures, and proposed a formula for predicting wave run-up. The new formula was verified by comparison with the field measurements obtained by [17]. A wave run-up monitoring system installed on seawalls to measure the wave run-up heights were conducted by [18]. The empirical formulas recommended in the [19,20] were adopted to estimate the run-up height. Influential parameters including Iribarren number based on an equivalent slope, the roughness of the slope, berm, shallow foreshore and oblique wave attack are taken into account to make the empirical formula more in line with the actual application environment. The good agreement between the forecasted and measured wave run-up heights under storm conditions indicated that the proposed approach could be used for forecasting of wave run-up on real seawalls.

In recent years, much effort has been made to improve numerical models in order to study wave-structure interactions. Refs. [21,22] solved the Reynolds Averaged Navier-Stokes (RANS) equations to simulate the complex free surface of breaking waves in the surf zone. Ref. [23] applied this numerical model and conducted laboratory experiments to study the characteristics of flow field and solitary wave run-up. Comparisons of numerical results and experimental data showed fairly good agreement in terms of free-surface elevation. Ref. [24] investigated the effects of bed shear stress in an idealized swash flow by solving large-eddy simulation (LES). Ref. [25] focused on wave breaking and wave circulation processes, as well as turbulence exchange in the surf zone.

In this study, we focus on periodic wave attacks on the most common slope-type seawall in Taiwan, and explore the distribution of pressure and undertow on the seawall surface, as well as the bottom shear stress caused by the flow field. A two-dimensional numerical model was developed to investigate the run-up and run-down of periodic waves on a seawall, involving a vertical wall and seawalls with slopes of 1:2 and 1:5. The unsteady two-dimensional Reynolds Averaged Navier-Stokes (RANS) equations and $k-\epsilon$ turbulence model were solved to simulate the viscous flow field near the seawalls. The complex free-surface configuration was captured using the particle level set method. The solid-fluid coupling method developed by [26] was employed to simulate the solid-fluid interaction on fixed Cartesian grids. The free-surface evolutions, the wave pressure on the seawalls, the bottom shear stress prior to the seawalls and the time-averaged flow field are discussed. Under the simulation conditions of this study, it is possible to clarify the causes of the reverse flow of the sloping bottom and its influence on the bottom shear stress, and to explore the standing wave characteristics in front of the seawall.

2. Governing Equations and Boundary Conditions for the Numerical Model

For an incompressible, viscous fluid, the continuity equation and the unsteady Reynolds Averaged Navier-Stokes (RANS) equations can be written as

$$\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0,$$

(1)

$$\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}=\frac{-1}{\rho }\frac{\partial p}{\partial x}+\nu \left(\frac{{\partial }^{2}U}{\partial x^2}+\frac{{\partial }^{2}U}{\partial y^2}\right)+\frac{\partial }{\partial x}\left(2{\nu }_t\frac{\partial U}{\partial x}-\frac{2}{3}k\right)+\frac{\partial }{\partial y}\left[{\nu }_t\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)\right],$$

(2)

$$\frac{\partial V}{\partial t}+U\frac{\partial V}{\partial x}+V\frac{\partial V}{\partial y}=\frac{-1}{\rho }\frac{\partial p}{\partial y}+\nu \left(\frac{{\partial }^{2}V}{\partial x^2}+\frac{{\partial }^{2}V}{\partial y^2}\right)+\frac{\partial }{\partial x}\left[{\nu }_t\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)\right]+\frac{\partial }{\partial y}\left(2{\nu }_t\frac{\partial V}{\partial y}-\frac{2}{3}k\right),$$

(3)

where $U$ and $V$ are the time-averaged horizontal and vertical velocity components in the Cartesian coordinates ($x$,$y$); t is the time; $p$ is the hydrodynamic pressure; $\rho$ is the fluid density; $\nu$ is the kinematic viscosity; ${\nu }_t$ is the eddy viscosity; and $k$ is the turbulent kinetic energy.

In the present study, to take the wall damping effect into account, low-Reynolds-number $k-\epsilon$ models were adopted, involving empirical constants and additional terms expressed as follows:

$$\begin{array}{c}\frac{\partial k}{\partial t}+U\frac{\partial k}{\partial x}+V\frac{\partial k}{\partial y}=\frac{\partial }{\partial x}\left(\frac{{\nu }_t}{{\sigma }_k}\frac{\partial k}{\partial x}+\nu \frac{\partial k}{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{{\nu }_t}{{\sigma }_k}\frac{\partial k}{\partial y}+\nu \frac{\partial k}{\partial y}\right)+[2{\nu }_t{\left(\frac{\partial U}{\partial x}\right)}^2+\\ 2{\nu }_t{\left(\frac{\partial V}{\partial y}\right)}^2+{\nu }_t\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)]-\epsilon ,\end{array}$$

(4)

$$\begin{array}{c}\frac{\partial \epsilon }{\partial t}+U\frac{\partial \epsilon }{\partial x}+V\frac{\partial \epsilon }{\partial y}=\frac{\partial }{\partial x}\left(\frac{{\nu }_t}{{\sigma }_{\epsilon }}\frac{\partial \epsilon }{\partial x}+\nu \frac{\partial \epsilon }{\partial x}\right)+\frac{\partial }{\partial y}\left(\frac{{\nu }_t}{{\sigma }_{\epsilon }}\frac{\partial \epsilon }{\partial y}+\nu \frac{\partial \epsilon }{\partial y}\right)+C_{\epsilon 1}{f}_1\frac{\epsilon }{k}[2{\nu }_t{\left(\frac{\partial U}{\partial x}\right)}^2+\\ 2{\nu }_t{\left(\frac{\partial V}{\partial y}\right)}^2+{\nu }_t\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)]-C_{\epsilon 2}{f}_2\frac{{\epsilon }^2}{k}+E,\end{array}$$

(5)

where $C_{\epsilon 1}$, $C_{\epsilon 2}$, ${\sigma }_k$, and ${\sigma }_{\epsilon }$ are the dimensionless empirical constants. $f_1$ and $f_2$ are the empirical constants. Ref. [27] recommended the following empirical constants in a full turbulent flow field for attached boundary-layer calculations:

$$C_{\epsilon 1}=1.14; C_{\epsilon 2}=1.92; {\sigma }_k=1.0; {\sigma }_{\epsilon }=1.3,$$

(6)

$$f_1=1.0; f_2=1-0.3\times exp\left(-R_T^2\right); E=2\nu {\nu }_t{\left(\frac{{\partial }^{2}U}{\partial y^2}\right)}^2; R_T=\frac{k^2}{\nu \epsilon }$$

(7)

Ref. [28] proposed the following terms in a low-Reynolds-number $k-\epsilon$ model to modify the general turbulent transport equation:

$${\nu }_t=C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon }$$

(8)

where

$$C_{\mu }=0.09; f_{\mu }=exp\left[\frac{-3.4}{{\left(1+R_T/50\right)}^2}\right]$$

(9)

The present model describes the turbulent transport with completed terms, including turbulence production, advection, diffusion, and dissipation, which are meaningful for the evolution of turbulent kinetic energy.

For numerical computation, it is easier to solve the governing equations in their dimensionless form. In this work, the still water depth $h_o$ was used to non-dimensionalize $x$ and $y$, and the characteristic velocity $u_o$, defined as $c{H}_i/h_o$, is used to non-dimensionalize u and v, where $c$ is the phase velocity, and $H_i$ is the incident wave height. Consequently, $h_o/u_o$, $\rho u_o^2$, $u_o^2$, $u_o^3/h_o$, and $u_o{h}_o$ are used to non-dimensionalize $t$, $p$, $k$, $\epsilon$, and ${\nu }_t$, respectively. The Reynolds number is then defined as $Re=u_o{h}_o/\nu$.

To obtain solutions for the RANS equations and the turbulent transport equations, it is necessary to provide appropriate boundary conditions at all boundaries of the solution domain, as well as the initial conditions at $t=0$ in the entire domain. The boundary conditions that must be satisfied are: (1) the kinetic and the dynamic free-surface boundary conditions at the interface $\mathsf{\Gamma }$ between the air and water phase, (2) the upstream boundary on the wavemaker plate, and (3) the fixed boundary condition including the boundaries of a rigid, impermeable breakwater and wave tank, where the velocities, turbulent kinetic energy, and dissipation rate are set to zero. The initial conditions of the velocities, hydrodynamic pressure, and the surface displacement are also set to zero.

3. Numerical Method

In the proposed model, the governing equations were discretized using the finite analytical method [29], and the staggered grid system, developed by [30], was adopted to prevent odd-even decoupling between pressure and velocity. The level set method was utilized to determine the unknown free surface or to identify the location of the liquid-solid interface, and the projection method was adopted to solve the coupling of velocity and pressure.

The present numerical model was developed on Cartesian coordinates. Since the sloping seawall surfaces do not always fit neatly into such a coordinate system, it is necessary to discuss which method can be applied to irregularly shaped domains. The iterative calculation near a solid boundary is not the same as it is in a fluid region. The proposed study applies an innovative method for simulating the interaction of the flow with stationary structures of irregular shape based on Cartesian grids. Two boundary velocities were imposed, referred to as the solid and mass-conserving boundary velocities, to satisfy the no-slip boundary condition and mass conservation in the ghost cells around the immersed solid boundary. For more details of used techniques for the solid-liquid interface, please refer to [26].

The evolution of the level-set function was solved using the fourth-order TVD Runge-Kutta method [31] and the fifth-order WENO scheme [32].

The level set method is a numerical technique originally proposed by [31] to elucidate the evolution of interfaces and shapes. The advantage of the level set method is that one can perform numerical computations involving curves and surfaces with an Eulerian approach (with a fixed Cartesian grid). Additionally, the level set method makes it easy to follow shapes that alter topology, e.g., when a shape splits in two, develops holes, or the reverse of these operations. Generally, the function $\varphi$ is defined as the signed normal distance from the interface $\mathsf{\Gamma }$ and satisfies $\left|\nabla \varphi \right|=1$. In this work, the interface $\mathsf{\Gamma }$ between the air and water is the zero level set of a smoothed distance function $\varphi \left(\underset{\_}{x},\mathrm{t}\right)$, where $\underset{\_}{x}$ is specified by $\left(x, y\right)$. $\varphi >0$ and $0<\varphi$ denote the air region and water region, respectively. However, numerical diffusion occurs after a finite amount of computational time, i.e., the level-set function $\varphi$ becomes irregular and ceases to be a distance function. Thus, the level-set function $\varphi$ must be re-initialized at each time step to ensure that the level-set function $\varphi$ remains a smooth distance function. This re-initialization can be performed by iterating the following partial differential equation to a steady state and then replacing $\varphi \left(\underset{\_}{x},\mathrm{t}\right)$ with $d\left(\underset{\_}{x},\mathrm{t}\right)$ [33]. Hence, the level-set function remains a distance function for which $\left|\nabla \varphi \right|$ converges to unity without changing its zero-level set. In the numerical implementation, however, mass error may arise during the re-distancing procedure. To overcome the drawbacks of the level set method, Ref. [34] proposed the hybrid particle level set method to improve the mass conservation of the traditional level set method and reduce numerical diffusion. The underlying concept is to use Lagrangian marker particles to adjust the level-set function in the vicinity of the free surface.

4. Verification of the Numerical Model

This section examines the accuracy of the numerical scheme as well as the suitability of the model equations through a series of numerical experiments. The involvement of spatial variations in water level on the sloping bed and the distribution of dynamic pressure on the front of seawall models are demonstrated with experimental data.

4.1. Free-Surface Elevation

Ref. [4] performed experiments to investigate variations in the free surface, including spilling breakers as a periodic wave propagated toward a slope seawall. The experiments were conducted in a 40 m-long, 0.6 m-wide and 1.0 m-high wave tank, and a 1:35 mildly sloping was located downstream. The experimental data were used to test the validation of the present model for a wave-seawall interaction. The Cartesian coordinate system was used, and the origin ($x$,$y$) = (0 m, 0 m) was set up on the mild slope, where the water depth was 0.38 m. The positive $x$ was in the horizontal direction of the wave motion, and the positive $y$ was defined in the vertical upward direction, as measured from the still water level. Figure 1 is a diagram of the experimental setup and its coordinate system. The wave gauges were deployed to measure the wave height on the sloping seawall. The still water depth was $h_o$ = 0.45 m, the incident wave period was $T$ = 2 s, and the wave height was $H_i$ = 0.125 m. The numerical coordinate system was the same as that of the experimental setup, and cnoidal wave theory [35] was used as the incident wave condition at the upstream boundary.

The mesh sizes in the computational domain are arranged, and the time step is set to 0.002 s. The finer grids are arranged in the area close to sloping bed, where the waves will be relatively complex due to wave shallowing, reflections, and even breaking waves. However, high-resolution mesh arrangement means more computing resources are required. Thus, it is necessary to try different mesh resolutions to achieve the expected simulation results. The final mesh size used in this study is as follows, and the simulation results are in good agreement with experimental results. There were 30 uniform grids from $y/h_o=0$ to $y/h_o=1.0-1.5H_i/{\mathrm{h}}_o$. In order to accurately describe the complex wave evolution, the vertical uniform grid spacing was 1/10 times the wave height from $y/h_o\ge 1.0-1.5H/h_o$. The horizontal uniform grid spacing was 1/200 and 1/320 times the wavelength in the numerical wave tank, where the finer grid spacing was from $x/L=-0.2$ to the end of the numerical tank.

Figure 2 presents the spatial variations in the maximum water level (${\mathsf{\zeta }}_{\max } - \overline{\zeta }$), mean water level ($\overline{\zeta }$), and minimum water level (${\zeta }_{\min } - \overline{\zeta }$) in the wave tank. The present numerical results were compared with laboratory data [4] and numerical results [36,37]. As shown in the figure, the present model yields more accurate crest elevations before entering the breaking wave region, but results in slight overprediction in the inner surf zone compared with [36,37]. A comparison between the numerical results and the experimental data revealed that the proposed numerical model can accurately simulate the evolution of a free surface during wave propagation on a sloping bottom.

It should be mentioned that the results of [4] are used to demonstrate the accuracy of present numerical model for simulating the interaction between waves and seawalls with gentle slopes. The wave breaking patterns can be different on very mild and steep slopes what we discuss in this study.

4.2. Wave Pressure

In coastal engineering, the determination of the wave force on an ocean structure provides important information for the design of the structures. The proposed numerical model could be applied to calculate the wave pressure on a seawall. The accuracy of the model was verified by comparing the numerical results with experimental data [38]. Figure 3 provides a diagram of the experimental setup and its coordinate system. The experiment was carried out in a 72.5 m-long and 2 m-wide wave flume. Two different seawall models were arranged at the furthest end from the wavemaker over a rigid sloping bed, where the height of the flaring shaped seawall (FSS) and vertical seawall (VW) was $B$ = 0.81 m. The initial still water depth was $h_o$ = 1.0 m, and the water depth at the toe of the seawall was $d$ = 0.54 m. The incident wave was described using cnoidal wave theory [35]; the wave height was $H_i$ = 0.05 m, and the wave period was $T$ = 3.0 s. The definition of the shape of the FSS is shown as:

$$\begin{gathered} \frac{x}{B}=\frac{-exp\left[cot\left(\frac{\pi }{3}\right)\tau \right]}{B\times 81\times {10}^6}\times cos\tau -0.12, \\ \frac{y}{B}=\frac{-exp\left[cot\left(\frac{\pi }{3}\right)\tau \right]}{B\times 105\times {10}^6}\times sin\tau +0.135, \end{gathered}$$

(10)

where $28.8\le \tau \le 31.2$.

The mesh sizes in the computational domain are arranged. There were 20 uniform grids from $y=0 \mathrm{m}$ to $y=0.5 \mathrm{m}$. In order to accurately describe the complex wave evolution, the vertical uniform grid spacing was 0.01 m from $y\ge 0.5 \mathrm{m}$. The horizontal uniform grid spacing was 0.01 m and 0.005 m, where the finer grid spacing was from the toe of the seawall to the end of the numerical tank. The time step is 0.003 s.

Figure 4 presents a comparison between the experimental data and the numerical results for the maximum wave pressure distribution on the seawalls [38]. In Figure 4, $P_w$ is the maximum wave dynamic pressure, $H$ is the wave height, and $z=y-h_o$ is relative depth of the wall. The circles and triangles correspond to the experimental data for the FSS and VW, respectively, and the solid blue line and dashed red line correspond to the numerical results for the FSS and VW, respectively. The comparison shows that the numerical results are in good agreement with the experimental data obtained by [38]. The numerical results show that the present model can be used to properly simulate the interaction between viscous flow and an irregular structure.

5. Numerical Results and Discussion

As waves hit a seawall, the evolution of the wave field in front of the seawall is affected by the water depth, wave conditions, and the shape of the seawalls. The wave fields become more complex due to wave reflection, diffraction, shoaling, and breaking. In most of the previous studies on this topic, solitary waves were commonly used to study the wave-structure interaction, because the physical mechanism was relatively easy to grasp. However, it is to be expected that flow fields with periodic waves near a seawall are more in line with actual physical problems than with solitary waves. In this section, the proposed numerical model was applied to investigate the propagation of a periodic wave over a flat bottom and the impact on seawalls. The physical characteristics, including free-surface evolution, wave pressure on seawalls, bottom shear stress prior to seawalls, and time-averaged flow fields, are discussed.

Two different periodic waves and three types of seawalls were specified, as listed in

Table 1. Numerical conditions of the incident periodic wave and the numerical results of the wave run-up.

Case	Slope (S)	${\mathit{H}}_{\mathit{i}}$ (m)	${\mathit{h}}_{\mathit{o}}$ (m)	$\mathit{T}$ (s)	$\mathit{L}$ (m)	${\mathit{H}}_{\mathit{i}}/\mathit{L}$	$\mathit{U}\mathit{r}$	$\mathit{R}/{\mathit{H}}_{\mathit{i}}$	${\mathit{\xi }}_{\mathit{o}}$
1	Vertical	0.05	0.35	1.266	2.00	0.025	4.7	0.88
2	1:2							2.02	3.16
3 *	1:5							1.36	1.26
4	Vertical			1.854	3.20	0.016	11.9	0.91
5	1:2							2.24	4.00
6 *	1:5							1.96	1.60

* breaking wave case.

. $T$ is the wave period, $H_i$ is the wave height, $h_o$ is the still water depth, $L$ is the wavelength, and $S$ is the slope of the seawall. When wave conditions are determined, the wave characteristic parameters are also determined at the same time, including $H_i/L$ (wave steepness), $R/H_i$ (the ratio of maximum run-up height to wave height), $Ur$ (Ursell number), and ${\xi }_o$ (Iribarren number). The Ursell number and the Iribarren number are defined as follows, respectively:

$$Ur=\frac{H_i{L}^2}{h_o^3}$$

(11)

$${\xi }_o=\frac{tan\alpha }{\sqrt{H_{i}L}}$$

(12)

The Ursell number represents the nonlinearity of the incident waves. Higher Ursell number can be regarded as the incident wave has higher harmonic components. After having verified the accuracy of the numerical scheme, the effects of the Ursell number and the Iribarren number are systematically discussed in terms of the run-up height, instantaneous and time-averaged flow fields, wave pressure on seawalls and bottom shear stress.

The Iribarren number is used to represent the effects of seawall slope and wave steepness. It means that if the seawall slope is larger or wave steepness smaller, the Iribarren number increase. The effects of the Iribarren number are considered to discuss in terms of the maximum wave run-up height for each case.

In this study, vertical wall and slopes of 1:2 and 1:5 for common seawalls along the coast of Taiwan were specially selected. In Case 3 and Case 6, breaking waves can be observed during wave run-up and run-down processes.

Figure 5 schematically illustrates a seawall located 15 wave lengths away from the wavemaker of a two-dimensional numerical wave tank. The second-order Stokes wave theory was used as the incident wave condition at the upstream boundary. The Cartesian coordinate system was used with the origin ($x$,$y$) = (0 m, 0 m) at the toe of the seawall. In the figure, the positive $x$ is in the horizontal direction of wave motion, and the positive $y$ is defined in the vertical upward direction, as measured from the bottom.

Appropriate simulation resolution can be obtained with different grid spacings while reducing computation time. As shown in Figure 6, non-uniform and staggered grids were adopted in the computation domain. To precisely capture the variations in the boundary layer flow, mesh refinement was performed near the bottom for a 16-times boundary-layer thickness ($\mathsf{\delta }$). There were 30 grids from $y/h_o=0$ to $y/h_o=3\delta /h_o$, 15 grids from $y/h_o=3\delta /h_o$ to $y/h_o=6\delta /h_o$, and 10 grids from $y/h_o=6\delta /h_o$ to $y/h_o=16\delta /h_o$. The boundary-layer thickness, $\delta$, was defined as:

$$\delta =\sqrt{2\nu /\omega }$$

(13)

where $\nu$ represents kinematic viscosity, and $\omega$ represents angular frequency. There were 30 uniform grids from $y/h_o=16\delta /h_o$ to $y/h_o<0.79$. In order to accurately describe the complex wave evolution, the vertical uniform grid spacing was set to 1/10 times the wave height from $0.79\le y/h_o$.The horizontal uniform grid spacing was set to 1/200 and 1/320 times the wave length in the numerical wave tank, where the finer grid spacing was from $x/L=-0.2$ to the end of the numerical tank. The mesh arrangement in the computational domain is based on [39].

To verify the accuracy of the proposed numerical scheme, the incident periodic wave in this numerical wave flume was compared with the theoretical solutions. Figure 7 shows a comparison between the periodic wave forms generated from the numerical results and the theoretical solutions. The black circles (○) and solid red line (-) correspond to the numerical results and theoretical solutions, respectively. A comparison between the numerical results and the theoretical solutions revealed that the proposed numerical model can accurately simulate the incident periodic wave. In addition, the analysis time $t^{*}$ = 0 identified the instant at which the periodic wave is stable.

5.1. Wave Run-Up on the Seawall

As a wave leaves a flat-bottom area and travels along a sloping area, the wave amplitude subsequently increases due to the decreasing water depth. The ratio of the local wave height to the wavelength, $H/L$, increases as the wave propagates onshore. As a result, the wave profile gradually loses its symmetry, and the wave becomes a breaking wave. When waves propagate towards and interact with the sloping seawall, the waves run-up on the face of the slope. As the waves reach the maximum run-up, the waves run-down the face of the slope due to the fact that the gravitational force is greater than the inertial force. Next, the waves run-up on the face of the slope again when the next rushing wave arrives. In this go-around and begin-again process, the reflected waves and incident waves together form a standing wave, and the waves’ energy simultaneously dissipates during the occurrence of the breaking waves.

In this section, the proposed numerical model was applied to investigate the wave run-up on seawalls induced by different types of periodic waves ($Ur$ = 4.7 and $Ur$ = 11.9) propagating on three different types of seawall, involving a vertical seawall and steep seawall with slopes of 1:2 and 1:5. Figure 8 shows the snapshot of free-surface elevation as the waves reached the maximum run-up in each numerical conditions. Table 1 shows the relative maximum wave run-up height ($R/H_i$). As shown in Figure 8, on the 1:2 sloping seawall the maximum wave run-up height was higher than on the 1:5 sloping seawall for the same incident wave conditions (Case 2 and Case 3; Case 5 and Case 6). The higher wave run-up on the 1:2 sloping seawall is likely related to the breaking wave phenomena in Case 3 and Case 6 (1:5 sloping seawall). The wave energy dissipates as the breaking wave occurs, so the wave run-up height elevation is lower. Under the same incident wave conditions (Case 1~Case 3; Case 4~Case 6), the maximum wave run-up height elevation on the vertical wall was the lowest. Since most of the wave energy was reflected back offshore, the maximum wave run-up height elevation on the vertical wall was the lowest in these cases.

Figure 8 and Table 1 also show that the relative wave run-up height ($R/H_i$) decreases as wave steepness ($H_i/L$) increases, by comparing the different incident wave condition propagations on the same type of seawall (Case 5 > Case 2 and Case 6 > Case 3). The comparison between Case 2 and Case 5 shows that the wave run-up height elevation increased with decreases in wave steepness. The same results were obtained in the comparison between Case 3 and Case 6. This indicates that waves easily deform and break with increases in wave steepness. As a result, the wave run-up height elevation on the seawall is lower due to the energy dissipation as the breaking wave occurs.

In addition, from the results in Table 1, we can also observe that the relative wave run-up height ($R/H_i$) increases as the Iribarren number (${\xi }_0$) increases (Case 3 < Case 6 < Case 2 < Case 5). The results are consistent with the empirical formula of EurOtop manual as follows:

$${\mathrm{R}}_{u2\%}/H_{m0}=1.65·{\gamma }_b·{\gamma }_f·{\gamma }_{\beta }·{\xi }_{m-1,0}$$

(14)

where ${\mathrm{R}}_{u2\%}$ is wave run-up height, which is exceeded by 2% of the number of incident waves; $H_{m0}$ is spectral wave height; where ${\xi }_{m-1,0}$ denotes the Iribarren number; ${\gamma }_b$, ${\gamma }_f$, ${\gamma }_{\beta }$ denote the influence factor for berm, roughness of the slope and oblique wave attack, respectively. It should be mentioned that the empirical formula is used to evaluate the wave run-up height caused by irregular waves.

5.2. Instantaneous Flow Fields

In coastal engineering, the most relevant hydraulic processes to be considered in the wave-structure interaction encompass the free-surface evolutions and the characteristics of flow fields, the wave pressure on seawalls, and the bottom shear stress prior to seawalls. In this section, the present numerical model was applied to simulate and investigate the viscous flow fields induced by periodic wave propagation on three different seawalls. In the following, the simulation results for the Case 2 and Case 3 are presented, and breaking waves can be observed in Case 3 during wave run-up and run-down processes. Both the wave run-up with the effects of breaking and non-breaking wave on the flow fields were investigated. In Figure 9 and Figure 10a the free-surface and flow field velocity, (b) the boundary-layer flow, (c) the bottom shear stress prior to the seawalls, and (d) the wave pressure on the seawalls is presented. Figure 9 and Figure 10 show the time history of instantaneous flow fields induced by periodic wave ($Ur$ = 4.7) propagation on the steep seawalls with slopes of 1:2 and 1:5, and Figure 9 and Figure 10I–IV show the corresponding flow fields at four-time instances.

When the periodic waves propagate on a sloping seawall, part of the waves run-up on the face of the seawall, and others reflect back offshore. As a result, the partially reflected wave and the incident wave together form a partial standing wave. The vertical velocity reaches the maximum at the anti-node ($x/L\approx 0.05,-0.45$) of the partial standing wave, while the horizontal velocity reaches the maximum at the node ($x/L\approx 0.21,-0.18$) of the partial standing wave, as shown in Figure 9I. When the anti-node of the partial standing wave closest to the seawall reaches its maximum elevation, the kinetic wave energy totally converts to potential energy. As a result, the flow velocity in front of the seawall is close to zero, and the maximum negative wave pressure is simultaneously caused by this partial standing wave, as shown in Figure 9II. On the contrary, when the anti-node of the partial standing wave closest to the seawall returns to its minimum elevation, the wave once again converts kinetic energy to potential energy. At the same time, the wave reaches the maximum wave run-up height elevation, and the maximum positive wave pressure is simultaneously caused by this partial standing wave, as shown in Figure 9IV. In Figure 9II,IV, we can observe that the reverse flow phenomenon occurred due to the differences in the free-surface elevations.

The wave shoaling process becomes longer as the incident wave propagates on a mildly sloping seawall, and the wave height and steepness increase due to wave shoaling. As a result, the wave profile gradually loses its symmetry, and the wave becomes a breaking wave. The reflected wave energy decreases due to the dissipation that occurs in the breaking wave. Therefore, the standing wave in Figure 10 cannot be observed. The evolution of the flow field on the slope varies rapidly due to the effect of fluid viscosity and non-linear effects. As a result, the wave profile gradually loses its symmetry, and the wave becomes a breaking wave, as shown in Figure 10I. In Figure 10II,III, it can be observed that the wave amplitude increases until the wave collapses and becomes a plunging breaking wave. As this moment, the wave releases a large amount of wave energy and the creates significant pressure on the seawall. Afterwards, the kinetic wave energy dissipates and causes part of the wave to run-up the face of the sloping seawall. Figure 10IV shows the maximum wave run-up height elevation.

5.3. Time-Averaged Flow Fields

In order to study the characteristics of flow fields’ long-term effects on periodic waves near seawalls, this section discusses the numerical results of the flow fields under steady state. The time-averaged flow field, the maximum wave pressure on seawalls, and the horizontal velocity profile in the boundary layer are discussed. The results for the non-breaking wave (Case 2) and the breaking wave (Case 3) are presented, respectively. As shown in Figure 11 and Figure 12, (a) the time-averaged flow field, (b) the averaged horizontal velocity profile, (c) the averaged horizontal velocity profile in the boundary layer, and (d) the maximum wave pressure on the seawall is presented.

Figure 11 shows the time-averaged simulation results for the periodic wave ($Ur$ = 4.7) propagating on the 1:2 sloping seawall. The partial standing wave motion and the steady recirculating cell can be observed in front of the steep seawall. In this case, the steady recirculating cell becomes asymmetric due to a free-surface limitation and the extrusion of the seawall, which leads to the steady recirculating cell in front of the steep seawall being less clear than that in front of the vertical wall. The key mechanism of scouring is the steady streaming that occurs in the vertical plane caused by partial standing waves in front of a seawall. A moderate understanding of the characteristics of the steady recirculating cell formed by partial standing waves would be helpful with regard to preventing severe scouring in front of seawalls, which could affect the safety of seawalls and structures.

Figure 12 shows the time-averaged simulation results for the periodic wave ($Ur$= 4.7) propagation on the 1:5 sloping seawall. As the wave propagates on a mildly sloping seawall, the wave steepness rapidly increases due to the wave shoaling process. This process results in wave profile asymmetry, which then leads to wave breaking. The breaking wave phenomenon on the mildly sloping seawall can be observed in Figure 12. Waves exert enormous force on the seawall because the wave largely releases its energy during breaking. A comparison of Figure 12d with Figure 11d shows that the maximum positive wave pressure on the 1:5 sloping seawall is five times larger than that on the 1:2 sloping seawall. The time-averaged flow field becomes more complex as the breaking wave occurs. There are no standing wave phenomena or steady recirculating cells, as shown in Figure 12. Once the wave form has been destroyed, the remaining water moves up the shore as swash and returns under the force of gravity as ‘undertow’ (indicated in Figure 12a by the red arrow). Undertow is a subsurface flow of water returning seaward from the shore due to wave action. Undertow is considered one of the dominant mechanisms in the erosion of beaches. It is important for coastal protection and seawall safety to understand the characteristics of undertow in surf zones.

In this section, the present numerical model was applied to investigate the characteristics of time-averaged flow fields near the different types of seawall. We can observe that the same periodic incident wave propagation on different sloping seawalls (seabed) causes separately specified physical phenomena. These phenomena may cause damage to the seawall or the erosion of beaches. Therefore, a deeper understanding of the characteristics of flow fields near seawalls and structures is useful for coastal protection.

5.4. Wave Pressure on Seawalls

In this section is presented an attempt to study the characteristics of wave pressure on seawalls, using a vertical seawall and steep seawall with slopes of 1:2 and 1:5, respectively. The present numerical model was applied to simulate two periodic incident waves ($Ur$ = 4.7 and $Ur$ = 11.9) propagating on each type of seawall. Figure 13a–c shows a comparison of the maximum wave pressure distributions with different wave conditions on a vertical wall, a 1:2 sloping seawall, and a 1:5 sloping seawall, respectively.

As shown in Figure 13a, the wave pressure distributions on the vertical wall concur with the theory of [40]. By comparing Figure 13b with c, one can observe slightly different characteristics of the wave pressure distributions between the vertical wall and the 1:2 sloping seawall. The positive and negative wave pressure on the 1:2 sloping seawall during the wave run-up processes can be simultaneously observed in Figure 9IV. It is speculated that this phenomenon is related to the evolution of the free-surface elevation on the sloping seawall. Positive wave pressure occurs on the sloping seawall when the free-surface elevation is higher than the still water level; on the contrary, negative wave pressure occurs on the sloping seawall when the free-surface elevation is lower than the still water level. Figure 12c shows that the incident waves create a huge amount of wave pressure on the 1:5 sloping seawall due to the occurrence of wave breaking. Figure 13 shows that the incident wave with $Ur$ = 11.9 (red mark) acts on the three seawalls, resulting in a wave pressure greater than that of the incident wave with $Ur$ = 4.7 (blue mark). In particular, on the 1:5 sloping seawall the maximum positive wave pressure was 2–6 times larger than that on the 1:2 sloping seawall, as shown in Figure 13c. Therefore, it is important to carefully strengthen seawalls to prevent long waves with breaking waves from acting on mildly sloping seawalls (seabeds).

5.5. Bottom Shear Stress

This section discusses variations in the bottom shear stress prior to seawall involvement, using a vertical seawall and steep seawall with slopes of 1:2 and 1:5, respectively. To show each calculation point horizontally with respect to the toe of the seawall, a Cartesian coordinate system is used with the origin $x/L$ = 0 at the toe of the seawall. The positive $x$ is in the horizontal direction of the wave motion. The dimensionless bottom shear stress ($\mathsf{\tau }$) is defined as the ratio of the dimensional bottom shear stress to ${\mathsf{\tau }}^{*}$ [41]:

$${\mathsf{\tau }}^{*}=\rho \sqrt{\frac{\nu }{\omega }}\frac{g\left(H_i/2\right)k}{coshk{h}_o}$$

(15)

where $\nu$ is the kinematic viscosity; $g$ is referred to as the acceleration of gravity; $k$ and $\omega$ are the wave number and circular frequency and denoted as $2\pi /L$ and $2\pi /T$, respectively.

Figure 14a,b show a comparison of the bottom shear stress distributions prior to different seawall types with wave conditions of $Ur$ = 4.7 and $Ur$ = 11.9, respectively. From Figure 14, it can be seen that the absolute value of the maximum positive and negative bottom shear stress at the toe of seawall is vertical wall (VW) < S = 1:2 < S = 1:5. The horizontal velocity in the positive x-direction decreases due to the (partial) reflected wave near the seawall. The offsetting effect is more obvious when the seawall slope is steeper, while the maximum bottom shear stress at the toe of the seawall is smaller. As the periodic waves propagate on the 1:5 sloping seawall, the undertow is a subsurface flow of water returning seaward from the shore due to occurrence of wave breaking. Therefore, the maximum negative bottom shear stress occurred at the toe of the 1:5 sloping seawall. The bottom shear stress was zero at the toe of the vertical wall due to the offsetting effect of the incident wave and reflected wave. Observing the entire bottom shear stress distribution, the bottom shear stress prior to the 1:2 sloping seawall is the smallest. It is speculated that there is undertow generation as the breaking wave occurs on the 1:5 sloping seawall and the total reflected wave occurs in front of the vertical wall.

6. Conclusions

The unsteady two-dimensional Reynolds Averaged Navier-Stokes (RANS) equations coupled with the $k-\epsilon$ turbulence model were numerically solved to investigate the viscous flow field induced by periodic waves near seawalls. Comparing the numerical results with those of previous studies, the validity of this model for simulating wave-structure interaction is demonstrated. The proposed numerical model was applied to investigate the viscous flow fields induced by periodic waves propagating on three different types of seawalls commonly designed in Taiwan, including a vertical seawall and steep seawall with slopes of 1:2 and 1:5. Based on the numerical results, we draw the following conclusions:

• Comparison between the numerical results and the experimental data revealed that the proposed numerical model can accurately simulate the evolution of the free surface and the maximum wave pressure on the structure during wave propagation on a sloping bottom.
• By comparing the different incident wave condition propagations on the same type of seawall, it can be seen that the wave run-up height elevation increases with decreases in wave steepness. As a result, the wave run-up height elevation on the seawall is lower due to the energy dissipation as the breaking wave occurs. The wave run-up height elevation on the vertical wall was the lowest. Since most of the wave energy was reflected back offshore.
• Numerical results reveal that the relative wave run-up height ($R/H_i$) increases as the Iribarren number increases and vice versa. The results are consistent with the empirical formula of EurOtop manual.
• For periodic waves propagating on a steep-slope seawall, it is found that the reflected wave and the incident wave can together form a partial standing wave or asymmetric recirculating cells. As the incident wave propagates on a mildly sloping seawall, the wave profile gradually loses its symmetry, and the wave becomes a breaking wave. The reflected wave energy decreases due to the dissipation that occurs in the breaking wave. Therefore, the standing wave phenomena or the steady recirculating cells cannot be observed.
• Positive wave pressure occurs on the sloping seawall when the free-surface elevation is higher than the still water level and vice versa. The incident wave with $Ur$ = 11.9 acts on the three seawalls, resulting in a wave pressure greater than that of the incident wave with $Ur$ = 4.7. The maximum positive wave pressure on the 1:5 sloping seawall is 2–6 times larger than that on the 1:2 sloping seawall.
• As the flow velocity on the slope seawall becomes smaller, the onshore flow reverses near the bottom and in the boundary layer. This reverse flow is primarily caused by an adverse pressure gradient imposed on the boundary.
• The absolute value of the maximum positive and negative bottom shear stress at the toe of the seawall is the vertical wall (VW) < S = 1:2 < S = 1:5. The maximum negative bottom shear stress occurs at the toe of the 1:5 sloping seawall due to undertow as wave breaking occurs.

Understanding of wave run-up process is import for protecting coastal structures. In present study, we show that the undertow plays a key role in wave run-up, especially when wave breaking occurrs. The boundary-layer flow on the sloping bed is closely related to the undertow and requires further study.