Numerical Analysis of Vertical Breakwater Stability under Extreme Waves

Abstract

The purpose of this study is to perform a numerical simulation of caisson breakwater stability concerning the effect of wave overtopping under extreme waves. A numerical model, which solves two-dimensional Reynolds-averaged Navier–Stokes equations with the k−ε turbulence closure and uses the volume of fluid method for surface capturing, is validated with the laboratory observations. The numerical model is shown to accurately predict the measured free-surface profiles and the wave pressures around a caisson breakwater. Considering the dynamic loading on caisson breakwaters during overtopping waves, not only landward force and lift force but also the seaward force are calculated. Model results suggest that the forces induced by the wave overtopping on the back side of vertical breakwater and the phase lag of surface elevations have to be considered for calculating the breakwater stability. The numerical results also show that the failure of sliding is more dangerous than the failure of overturning in the vertical breakwater. Under extreme waves with more than 100 year return period, the caisson breakwater is sliding unstable, whereas it is safe in overturning stability. The influence of wave overtopping on the stability analysis is dominated by the force on the rear side of the caisson and the phase difference on the two ends of caisson. For the case of extreme conditions, if the impulse force happens at the moment of the minimum of load in the rear side, the safety factor might decrease significantly and the failure of sliding might cause breakwater damage. This paper demonstrates the potential stability failure of coastal structures under extreme sea states and provides adapted formulations of safety factors in dynamic form to involve the influence of overtopping waves.

1. Introduction

It is important to provide guidelines of structure stability for designing coastal protection structures or harbor breakwaters. Many previous studies on wave-structure interaction have discussed the seaward and lifting forces including those induced by the effects of wave overtopping and breaking from the physical or numerical models. Among these studies, the most important historical failures of the vertical structure have been documented [1]. Large-scale hydraulic experiments have been performed over the past several decades, leading to different empirical formulations that allow for the calculation of loads induced by wave breaking on vertical breakwaters [2,3,4,5]. Some semi-empirical equations for impulsive loads on vertical breakwaters are based on experimental data and, in some cases, with prototype-measured data [6]. Recently, the most commonly semi-empirical formulas for calculating pressures on a vertical structure wave have been developed [7,8,9]. When the wave obliquely enters the elongated structure, the maximum wave force decreases [10,11]. Lately, research results have reported that an increase in wave pressure caused by the diffracted wave can be one of the reasons for destruction of the conventional caisson breakwater [12].

In the theoretical analysis, there are some idealized analytical models which have been presented to study the forces exerted on the vertical structures [13,14]. However, these analytical approaches are complicated and inconvenient due to the complex geometries. In order to overcome the limitations, the numerical models have been developed since they are more flexible and efficient. The nonlinear shallow water or Boussinesq equations are developed to simulate the wave–structure interaction in coastal processes, especially for porous structures [15,16,17,18,19,20,21,22]. Besides, a boundary element method was used to solve the unsteady Forchheimer equations and described the wave transmission and reflection by a multilayer breakwater with arbitrary shapes [23]. Recently, the development of more efficient Navier–Stokes solvers have allowed a more sophisticated modelling of the wave-structure interaction problem [24,25,26,27]. For instance, the Reynolds-averaged Navier–Stokes (RANS) approach employing the volume of fluid (VOF) method for the modelling of complex (turbulent) flows has become popular in the past two decades [28,29,30]. A robust model was presented to investigate the functionality of rubble mound breakwaters with special attention focused on wave overtopping processes [31]. Then a numerical analysis considering wave loads corresponding to a low-mound and a rubble-mound breakwater with both regular and irregular incident wave conditions was carried out [32]. The results showed that this RANS equation model had a high potential to become a complementary tool to analyze the hydraulic response of caisson structures.

However, the existing research described above only considered the forces on the seaward and bottom of the caisson breakwater. The harbor-side loads induced by wave overtopping on a caisson breakwater have been investigated [33] and concluded that the backward loads should be considered as failure mode when designing the caisson breakwater. Very little information has been reported on wave forces acting on the backward side of caisson breakwater, which is an important factor in calculating the sliding and overturning stabilities for the engineering design. In this study, we shall extend the numerical model based on the Reynolds-averaged Navier–Stokes (RANS) equations to simulate the extreme wave overtopping induced loads on the front, bottom and rear sides of caisson breakwater. The numerical results are validated with the experimental data of surface elevations and wave pressures acting on the caisson breakwater [34]. The numerical model is applied to actual caisson breakwaters under various extreme wave conditions and the results are compared to the empirical formulas to discuss the stability of sliding and overturning [35]. Finally, some conclusions are made.

2. The Numerical Model

In this paper we used a numerical model (COBRAS model) which is a depth and time resolving 2DV numerical model solving the Reynolds-averaged Navier–Stokes (RANS) equations and the $k-\epsilon$ turbulence closure model to check experimental data. The RANS equations for the ensemble-averaged velocity, $<u_i>$, and the ensemble-averaged pressure, <p>, are well-known and can be expressed as

$$\frac{\partial 〈u_i〉}{\partial x_i}=0$$

(1)

$$\frac{\partial 〈u_i〉}{\partial t}+〈u_j〉\frac{\partial 〈u_i〉}{\partial x_j}=-\frac{1}{\rho }\frac{\partial 〈p〉}{\partial x_i}+g_i+\frac{1}{\rho }\frac{\partial {\tau }_{ij}}{\partial x_j}-\frac{\partial 〈{u^{\prime }}_i^{}{u^{\prime }}_j^{}〉}{\partial x_j}$$

(2)

in which i, j = 1, 2 for two-dimensional flow and ${\tau }_{ij}$ is the viscous stress. $\rho$, $p$ and t denote the water density, pressure and time. The Reynolds stress $-\rho 〈{u^{\prime }}_i^{}{u^{\prime }}_j^{}〉$ is calculated with a nonlinear eddy viscosity relationship [36]. This nonlinear closure relationship for Reynolds stress mimics the more complete Reynolds stress closure without solving the transport equations of Reynolds stress. A demonstratation showed that standard linear eddy viscosity closure tends to over predict the diffusion of turbulence under breaking waves [29], and the nonlinear relationship suggested above predicts better turbulence statistics when compared with laboratory data of breaking wave over sloping beaches [37]. The governing equations are solved by a finite difference scheme. A two-step projection method is adopted to calculate the RANS equation. The evolution of free surface is calculated by the VOF method. The functionality of rubble mound breakwaters with special attention focused on wave overtopping processes uses the COBRAS model, and the numerical result comparison with experimental data has good agreement [38]. The mean overtopping discharge at the root of the South Breakwater of Póvoa de Varzim Harbour (Portugal) modelling by COBRAS model has good agreement with physical model tests [38].

A small-scale (1:36) physical model was carried out in a wave flume (100 m × 1.5 m × 2 m) to study the wave forces on a composite vertical breakwater [34], which was composed of an impermeable vertical caisson and permeable rubble foundation as shown in Figure 1. The specifications of the breakwater were as follows: the water depth was 0.526 m, the freeboard height was 0.153 m, the mound height was 0.138 m and the caisson weight was 780.66 kg. There were four wave gauges marked g1–g4 to measure the surface elevations and nine pressure transducers to measure wave pressures where five pressure transducers were along the front face of the caisson (U1–U5) and four pressure transducers were set on the bottom of the caisson (V1–V4). For the permeable rubble foundation, the median diameter was 0.7 cm and the porosity was 0.49. The incident wave height and period was experimental measured at wave gauge g1. Wave gauges (g2–g4) were used to calculate the reflection of the breakwater. A numerical flume was implemented with the same characteristics in the experimental setup. The computational grid system was discretized with non-uniform meshes in the x and z directions ($\Delta x$ = 0.1 m, $\Delta z$ = 0.05 m). Thus, the total number of cell meshes was 300 × 34. Additionally, the total simulation time was 64 s and the Courant number was 0.3 for all cases. The corresponding time step was automatically adjusted during calculations to satisfy stability constrains by both advection and diffusion processes, in which the maximum time step was adequately chosen as ${10}^{-2}$ s compared with experimental measurements. The aforementioned computational conditions were used throughout this study.

The comparison of wave surface elevations for the four wave gauges between the numerical results and experimental data for wave overtopping is shown in Figure 2. The agreement between numerical results and experimental data was good except the wave gauge g2, where the vicinity of a node of the partial standing wave effected the interaction of the incident wave and reflection wave. Figure 3 and Figure 4 show the comparison of wave dynamic pressure time series between the numerical simulations and experimental data for the regular overtopping wave condition (H = 0.25 m and T = 2 s). The pressure transducers U1~U2 were above the still water label. In Figure 3, the pressure transducers U3~U4 are exposed to air under the wave trough in that wave condition. The comparison includes five pressure transducers along the vertical front face of the caisson (U1–U5) and four pressure transducers underneath the caisson (V1–V4). The solid line represents the numerical simulations while the dashed line represents the laboratory measurements. From Figure 2, Figure 3 and Figure 4, good qualitative comparisons between numerical simulation and experimental data are investigated. It indicates that the numerical model could simulate the waves propagating through the porous medium. In general, the COBRAS numerical model has a good prediction for the dynamic wave pressure time series at the vertical face of the caisson than those beneath the caisson. To quantify the comparisons, a validation method proposed by Wilmott was used as Equation (3) [39]

$$\phi =1-\frac{{\displaystyle \sum {\left|{\mathrm{X}}_{num}-{\mathrm{X}}_{\exp }\right|}^2}}{{\displaystyle \sum \left({\left|{\mathrm{X}}_{num}-\overline{{\mathrm{X}}_{num}}\right|}^2+{\left|{\mathrm{X}}_{\exp }-\overline{{\mathrm{X}}_{\exp }}\right|}^2\right)}}$$

(3)

In Equation (3), X_num is the numerical result, X_exp is the experimental data and the “$\overline{}$” is an averaging operator. The index of validation (φ), ranging from 0 to 1, represents the degree of agreement between numerical predictions and experimental data, and the unit is dimensionless. The validated results are shown as

Table 1. The index of model validation.

Case	φ	Case	φ	Case	φ
g1	0.960	U1	0.723	V1	0.714
g2	0.638	U2	0.777	V2	0.741
g3	0.873	U3	0.772	V3	0.739
g4	0.920	U4	0.888	V4	0.708
--	--	U5	0.895	--	--

. The value of validated index, ranging from 0.71 to 0.96, is acceptable. However, it has to be noted that a smaller index is found for the wave dynamic pressure beneath the caisson breakwater. In general, the average value of all validated indexes is 0.8, which indicates that numerical predictions from the COBRAS model and experimental data are consistent. Consequently, the comparisons results provide confidence in the subsequent application of the model to the prototype scale, and which is described in the next section. In this study, we focus on the load on the rear sides of caisson breakwater induced by wave overtopping, so the overtopping discharge won’t be concerned.

3. Stability Analysis of Prototype Vertical Breakwater

3.1. Prototype Wave Load Analysis

The validated numerical model was used to evaluate the total wave forces on the caisson breakwater in Taiwan at a prototype scale. There are two parts of the caisson breakwater: the porous rubble foundation and vertical caisson with parapet [40]. The vertical caisson is installed on a rubble foundation, which has a thin filter layer between caisson and foundation and is placed on a foreshore slope of 1 on 100. The design wave condition for the breakwater of Taichung harbor is 50 year return period wave condition. The specifications of the breakwater is as follows: the water depth is 22 m, the freeboard height is 11 m and the mound height is 6 m. Figure 5 shows a vertical breakwater cross-section at Taichung harbor in Taiwan. Tests were simulated using regular waves for various values of wave conditions as shown in

Table 2. The wave conditions (including return period of typhoon wave and extreme condition).

Case	H (m)	T (s)	Return Period (Year)
1	7.77	11.30	50
2	8.46	11.70	100
3	9.24	12.20	200
4	9.45	12.30	250
5	12.87	14.13	50 *
6	13.70	14.63	100 *
7	14.47	15.25	200 *
8	14.63	15.38	250 *

mean water depth = 28.86 m at x = 0 m. * represents the climate change effected extreme wave condition [41].

. In order to understand the force contribution clearly, the regular wave was applied. Besides the condition of different return period typhoon waves, extreme waves due to climate change were also used [41]. According to [41], the wave height will increase up to 65% and the corresponding wave period might increase to 25% in 2039 around Taiwan. The numerical domain is 1125 m long and 80 m height. The mesh resolution at the structure is 0.2 m in the x-direction and 0.5 m in the y-direction. The rubble mounted layers beneath the caisson are defined using the porous flow parameter in

Table 3. The porous flow parameter for rubble mounted layers.

Layer	Porosity	The Median Diameter (m)
1	0.49	0.118
2	0.53	0.387
3	impermeable materials

.

In the past, the existing researches only considered the forces on the seaward and bottom of the caisson breakwater. Very little information has been reported on overtopping wave forces acting on the backward side of caisson breakwater, which is an important factor in calculating the sliding and overturning stabilities for the engineering design. In order to bring more insight on the wave forces induced by the overtopping waves, the numerical model is extended for simulating the wave pressure acting on the front, bottom and back sides of the caisson breakwater. A free body diagram shown in Figure 6 was used to analyze the overtopping wave induced forces and moments acting on the caisson breakwater. In Figure 6, W denotes the weight of caisson, f₁ denotes the horizontal force acting on the front face of the caisson, f₂ denotes the horizontal force acting on the rear face of the caisson and f₃ denotes the uplift force acting on the bottom of the caisson. Buoyancy force is included in still water. The fixed reference point of the moment for the landward overturning is denoted by the point s₁, while that for seaward overturning is denoted by the point s₂. The moment of f₁ is m₁, and the moments of f₂ is m₂. The moment of f₃ at fixed point s₁ is denoted by m₃, and the moment of f₃ at fixed point s₂ is denoted by m₄. The moment of W at fixed point s₁ is m_w1, and the moment of W at fixed point s₂ is m_w2.

The relative forces and moments acting on a vertical structure can be determined by integrating the numerically calculated pressure over a portion of vertical structure. Figure 7 presents the time histories of horizontal forces on the seaward-side (upper panel), harbor-side (middle panel) and uplift force for the studied case 1. In case 1, wave overtopping happens, indicating that all cases in Table 2 are necessary to investigate the influence of wave overtopping (IWO) on the structure stability. The force calculated included the behavior of flow transmitting through the rubble mound. In case 1, the overtopping wave happened after t = 75 s, so the transmitting flow disturbed a fluctuation of f₂ before overtopping. The calculated moments induced by the three forces are shown in Figure 7. From Figure 7, the oscillation of the rear force f₂ and the induced moment m₂ occurs after wave overtopping. Moreover, the peak values occasionally appear at some instantaneous time which might be caused by the jet of the overtopping wave. In addition, the differences between the water level on the front face and that on the rear side of caisson breakwater will influence the pressure distribution underneath the caisson. Therefore, the corresponding upward force f₃ and the induced moments (m₃ and m₄) will vary with the surface elevation difference for wave overtopping. The negative maximum of f₃ is larger than the positive maximum of f₃, which is the same as the experimental results in [34].

3.2. Stability Analysis

It is important to safely design a vertical breakwater with respect to the stability criteria of sliding and overturning. The movements of a caisson breakwater are influenced by the total active forces exerted on the seaward, harbor side and underneath the vertical structure. The minimum safety factors for sliding and overturning are suggested by [35] for engineering design and these values must not be less than 1.2 in practice. These safety factors for sliding and overturning are reviewed as the following:

Using an empirical formula to calculate the safety factors for sliding and overturning, the force acting on the rear side of the caisson was usually taken to be static for simplicity’s sake. Few literatures studied the horizontal force acting on the rear side along the caisson for overturning wave. Therefore, the total force determination requires a good understanding of the overtopping wave behavior in the front and rear sides of the caisson as well as beneath the caisson. Here, the sliding safety factor and overturning safety factor without the force on the backward side can be presented as Equations (4) and (5).

$$f_s=\mu \frac{W-f_3}{f_p}$$

(4)

$$f_0=\frac{m_W-m_3}{m_p}$$

(5)

where $\mu$ is the coefficient of friction between the caisson and the foundation, W denotes the weight of caisson, f_p denotes the maximum horizontal force along the front side of caisson, m_p denotes the moment of f_p at the point s₁, and m_w denotes the moment of W at the point s₁.

Due to the rear side force being considered, the safety factor against seaward sliding and overturning should be accounted. Consequently, the corresponding safety factor against landward sliding and seaward are presented as Equations (6)–(9), where superscript “+” and “−” represent landward and seaward, respectively. The subscript “s” and “o” mean the sliding and overturning.

$$f_s^{+}=\mu \frac{W-f_3}{f_1-f_2}$$

(6)

$$f_s^{-}=\mu \frac{W-f_3}{f_2-f_1}$$

(7)

$$f_o^{+}=\frac{m_{W1}-m_3}{m_1-m_2}$$

(8)

$$f_o^{-}=\frac{m_{W1}-m_3}{m_2-m_1}$$

(9)

For example, the time histories of the safety factors against sliding including sliding safety factor (f_s), landward sliding safety factor (f_s⁺) and seaward sliding safety factor (f_s⁻) in case 1 are shown as Figure 8, and that against overturning including overturning safety factor $f_0$, landward overturning safety factor $f_o^{+}$ and seaward overturning safety factor $f_o^{-}$ in case 1 are shown as Figure 9. The water wave in the rear side caused the sliding safety factor and overturning safety factor to oscillate, which means Equations (6)–(9) provide a time series variation value as shown in Figure 8 and Figure 9. Overall, the caisson breakwater in Taichung harbor has sufficient stability for sliding and overturning on the IWO in the 50 year return period wave condition.

4. Results and Discussions

To investigate the structure stability on the IWO, the sliding and overturning safety factors considering the water level fluctuates on the rear side of caisson were calculated under 50 to 250 year return period typhoon wave conditions as well as extreme wave conditions. The represented safety factors for sliding on the IWO determined by the minimum of landward sliding safety factor and seaward sliding safety factor was denoted by f_s^±. In the same way, the overtopping safety factors determined by the landward overturning safety factor and seaward overturning safety factor, was denoted by f_o^±. The comparisons of the safety factors against sliding are shown as Figure 10. Generally, the safety factors following severe wave conditions get lower, and cases 6–8 are less than 1.2, indicating the risk of sliding failure might occur in extreme wave conditions. Although f_s^± is closed to f_s, f_s^± varies due to the IWO.

The safety factors against overturning are shown as Figure 11. All of the cases are greater than 1.2. Therefore, sliding failure is more dangerous than the overturning failure in Taichung harbor breakwater. Comparing between Figure 10 and Figure 11, the tendency of safety factors between sliding and overturning are similar. This study shows the impact of overtopping waves in respect of the phase difference between the front and rear side. Not only wave condition but also the configuration of caisson such as its width and crown height might also affect the safety factor of breakwater. The impact on the fluctuation of the safety factor in case 2 is smaller than in case 1.

To focus on the IWO on the stability analysis of vertical breakwaters, the load analysis of case 4, an example of f_o^± < f_o, is shown as Figure 12. When time is 189.5 s, the minimum safety factor against sliding on the IWO f_s^± = 1.32 is greater than the safety factor against sliding in static f_s = 1.43. At that time, the force f₁ in Figure 12b is 4395 kN and the force f₂ in Figure 12c in dynamic is 1840 kN, while that in static is 2084 kN, and the force f₃ in Figure 12d is 5614 kN. Because the force on the rear side of caisson considers the IWO is less than the force in still water, that caused the corresponding safety factor against sliding to be lower. Therefore, the impulse force happens at the moment of the minimum of load in the rear side, the safety factor drops significantly and the failure of sliding might cause breakwater damage. Furthermore, overtopping safety factors in case 4, an example of f_o^± < f_o, are shown at Figure 13. At 189.6 s, the minimum safety factor against overturning on the IWO f_o^± = 1.921 is less than the safety factor against overturning in static f_o = 2.048. At the time point, the moment m₁ in Figure 13b is 42,230 kNm, the moment m₂ in Figure 13c in dynamic is 12,035 kNm; while that in static is 14,317 kNm and the moment m₃ in Figure 13d is 69,279 kNm. Since the moment m₂ on the IWO is less than that in static, the corresponding safety factor against overturning is lower.

In summary, the risk of sliding is greater than that of overturning in Taichung harbor breakwater. Under the extreme wave more than 100 year return period, the breakwater is unstable due to sliding, but safe from overturning. The IWO on the stability analysis is dominated by the force on the rear side of the caisson and the phase difference on the two ends of caisson.

5. Conclusions

To evaluate the stability of caisson breakwater concerning the effect of wave overtopping, a 2D VOF-type RANS model using a k−ε turbulence closure (COBRAS) has been demonstrated to be suitable for simulating the complex hydrodynamics induced by overtopping waves. Experimental wave profiles and dynamic wave pressures on the seaward, backward and the bottom of the caisson breakwater can be suitably simulated by the present model. The horizontal forces on the front and rear sides and the uplift force as well as the moments on the caisson could also be analyzed. According to the numerical simulations, the risk of sliding is greater than that of overturning in the case of caisson breakwater in Taichung harbor. Under the extreme wave more than 100 year return period, the caisson breakwater is unstable due to sliding, but that is safe for overturning. Therefore, it is important to calculate the force on the rear side of the caisson and the phase difference of surface elevation on the two sides of the caisson as we analyze the stability analysis for wave overtopping for engineering design.