# Numerical Modeling of Wave Reflection from Sloped Impermeable Seawalls Using the SPH Method: Case Study of Chabahar Port

Authors: Behrouz Aghaei; Afshin Mohseni Arasteh; Kamran Lari; Masoud Torabi Azad
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

In this research, a comprehensive study is performed to investigate the interaction of regular waves with the impermeable seawall of the Chabahar port. First, a MIKE 21 SW model is used to transform the deep-water wave data to the nearshore zone. Then, the interaction of waves with the seawall is simulated using a well-known numerical smoothed particle hydrodynamics model named DualSPHysics. After validating the numerical results with the experimental data, a parametric study is performed to evaluate the effects of the wave height, wave period, and the slope of the seawall on the water level fluctuations and the wave reflection coefficient. The results showed that increasing the wave height slightly decreases the reflection coefficient. Meanwhile, a direct relationship was found between the wave height and the water level fluctuations near the wall. Generally, increasing the wave period resulted in higher reflection coefficients and water level fluctuations. Both the reflection coefficient and the water level fluctuations are greatly dependent on the slope of the seawall. Steeper slopes resulted in higher reflection coefficients and lower water level fluctuations near the seawall.


---

## Body

## 1. Introduction

The primary goal of constructing a coastal defense structure is to protect berths, maneuvering areas, and port facilities against waves. In a seawall, the energy of the incident wave is dissipated through several mechanisms, including the wave breaking and porous flow in the mound structure, partially reflecting wave to the sea, as well as penetration and overtopping into the port [1]. As one of the most critical mechanisms of wave dissipation in a seawall, wave reflection might cause erosion along the face of the structure [2, 3]. Vertical impermeable seawalls present a high reflection coefficient of around 90–100%. Although the ability of a seawall to absorb the incident wave energy is represented by the reflection coefficient [4], too much reflection might result in the formation of standing waves, increased velocities for water particles, beach erosion, and eventually undermining the structure [5]. To overcome the problems mentioned above, sloped and rubble mound seawalls can be used.Generally, there are two main approaches to model wave reflections from coastal structures [6, 7]. In the first method, the structure is assumed to be impermeable. Then, an empirical roughness factor is used to adjust for permeability and roughness [8–12]. In the second method, permeable/rough structures can be directly applied in the experiment. Each method has some positive and negative points. In the first case, the simulation is based on a “perfect” situation not usually encountered in real life, but then, the results are adjusted through an empirical factor. In the latter method, a dynamic similarity is assumed. Then, the form of turbulence and energy dissipation will remain at the smaller laboratory scales than the full scale.To date, the behavior of coastal structures under wave actions has been the issue of different experimental [13–15], analytical [16–19], and numerical studies [20–24]. Numerical modeling can significantly reduce the time and cost needed for studying the wave-seawall interaction problems. The most common numerical approach for simulating the interaction of waves with coastal structures is the volume-of-fluid (VOF) method [25–29]. Since this method can only track the averaged quantity of the fraction function in each cell, it has some difficulties in resolving complex free surface features or those that are smaller than the mesh size. This problem can be solved using a variable mesh grid. However, some problems such as false breaking might arise when various full fluid cells are adjacent to empty cells and empty cells receive unreal fluid flow due to VOF function convection. The problem of convection function has been analyzed and improved by Hirt and Lin [30, 31].In recent years, particle methods have appeared as an alternative to mesh-based methods. Smoothed particle hydrodynamics (SPH) is one of the most famous particle methods that has been extensively applied in fluid mechanics and coastal engineering applications, such as landslides [32–34], two phases flow [35–39], wave interaction with floating bodies [40–46], and sloshing [47–51]. This method does not require predefined connectivity between points/nodes. Also, the SPH method has the advantage of simple tracking of free surface and moving boundaries [52]. The most crucial drawback of this method is its high computing times for large domains that is mainly due to the implementation of the neighbor list in this method [53, 54]. Although this method has some drawbacks, it has become an inherent part of the numerical arsenal of industrial research and development laboratories and academic industries. The main reasons are the growing needs of industry and research for appropriate tools for complex hydrodynamics, recent advancements in SPH theory that resolved some problems with this method, and the emergence of GPUs that eases the use of SPH codes to study 3D flows in real-life scales while keeping the computational times manageable [52].In this regard, there are some SPH studies on wave interaction with coastal defense structures. Considering a rubble mound breakwater, Ren et al. (2014) simulated the hydraulic stability of the armor blocks using a 2D DEM-SPH model [55]. The DualSPHysics model was used by Altomare et al. (2014) to simulate the interaction between regular waves and a rubble mound breakwater and to calculate the wave run-up [8]. Altomare et al. (2015) investigated the applicability of the DualSPHysics model to assess the forces exerted by waves on two coastal structures in Belgian ports [56]. Solitary wave interaction with an impermeable seawall was addressed by Varnousfaaderani and Ketabdari (2015) using a modified turbulence SPH method [57]. A three-dimensional WCSPH model was used by Wen et al. (2016) to simulate the interaction of waves with a vertical breakwater [58]. Yeganeh-Bakhtiary et al. (2017) used a two-dimensional WCSPH model to simulate the wave reflection from a caisson breakwater [59]. Xu and Lin (2017) used a two-step projection ISPH method for simulating the propagation of regular and solitary waves in a flat wave flume and calculating the transformation of waves over a submerged breakwater [60]. Fathi and Ketabdari (2018) simulated the behavior of semicircular breakwater under regular waves using SPHysics software [61]. Ketabdari et al. (2019) used the WCSPH method to simulate the wave reflection from a multipurposed fixed breakwater [62]. Liu et al. (2019) used the SPH method to simulate wave-winged box type breakwater interaction [63]. Subramaniam et al. (2019) compared the results obtained by DualSPHysics and OpenFOAM on wave run-up on a curved dike [64]. Luthfi et al. (2019) used DualSPHysics to simulate the properties of the tsunami wave on a coastal slope [65]. Dai et al. (2021) used the SPH method to simulate an actual debris flow in Tibet, China [34]. Xu et al. (2021) performed a 3D SPH simulation on dam-break interaction with different obstacles [66]. Dang et al. (2021) used the DualSPHysics model to simulate the wave forces and overtopping over different seawall geometries [67]. The latest developments of meshfree methods used in coastal and ocean engineering problems can be found in a review study performed by Gotoh and Khayyer (2018) [68].One of the grand challenges of the SPH method is representing the ability of the method to deal with real-life problems. Performing SPH simulations of wave interaction with actual structures is an important step in the commercialization of the method. Studying the wave interaction with a real coastal structure can be considered a real-life problem. In this research, a comprehensive study is performed on wave interaction with an impermeable seawall. First, the DHI MIKE software is used to obtain wave characteristics at the location of the structure. Then, the open-source SPH-based DualSPHysics [69] software is used to study the effects of different wave heights and periods and the slope of the sea wall of the Chabahar port on the free surface fluctuations near the wall as well as the reflection coefficient. It should be noted that numerical simulations are performed assuming that the seawall is impermeable. However, as previously stated, the obtained reflection coefficient data can be adjusted for different types of seawall using empirical roughness factors.


## 2. Design Parameters

Design parameters of this study are presented in this section.


### 2.1. Study Area

The seawall is located at the Shahid Beheshti port of Chabahar (25° 17′ N, 60° 37′ E), Iran’s only oceanic port and the closest point to the Indian Ocean. Chabahar is a focal point of development in the southeast of Iran by enhancing transit routes among Central Asian countries and the northern part of the Indian Ocean [70]. It has the highest draft (16 meters) compared with other ports of Iran. With an area of over 4000 square meters and the capacity for simultaneous reception of more than 600 passengers, Chabahar has one of the largest international passenger terminals in the region [71]. To develop the Chabahar port, the construction of a seawall was planned at the right-hand side of the port. The location of this seawall is shown in Figure 1.Figure 1
Aerial photo from Chabahar port and the location of the seawall.


### 2.2. Hydrodynamic Data at the Location of the Sea Wall

To determine the characteristics of deep-water waves in the study area, the wave characteristics obtained from a comprehensive numerical hindcast of wave conditions in the Gulf of Oman for twenty-three years (1985–2007) are used [72]. The hindcast showed that in this period, the maximum wave height in deep water was about 4.3 m with a southwest direction. Meanwhile, the direction of the dominant waves was the south.In addition to local winds, the area may also be affected by tropical cyclones that generate relatively large sea states. At the Chabahar port, tropical cyclone waves are the dominant design wave condition as the wave heights can exceed those generated from other sources like the monsoon and local westerly winds [73]. Therefore, statistics and information considering rare but severe events should also be considered to determine the design waves.By performing statistical analysis on the tropical storm historical data from various sources, reviewing previous studies on tropical storms, numerical simulation of waves caused by historical tropical storms, numerical simulation of waves generated by tropical storms by the Monte Carlo method, and finally, the statistical analysis of extreme value analysis of maximum wave height, Dibajnia and Allahyar (2007) provided the design wave in deep water for different locations along the Oman Sea [73]. For a return period of 100 years, these values are shown in Figure 2. For Chabahar port, the wave heights for different return periods are given in Table 1.Figure 2
Spatial variations in the maximum significant wave height with a return period of 100-year in the Oman Sea [73].Table 1
Variation in maximum significant wave height by the return period in Chabahar [73].
LongitudeLatitudeReturn period (years)10 (m)25 (m)50 (m)100 (m)60.525.24.25.87.28.8


### 2.3. Sea-Level Variations

Sea-level changes should be considered in the seawall design as well as the wave transformation from offshore areas to the nearshore zone. Generally, sea-level changes arise from various factors such as tides and storm surges. Based on the report of the Integrated Coastal Zone Management (ICZM) plan of Iran [74], different tidal levels at the Chabahar port are given in Table 2. In this table, MHHW, MLHW, MSL, MHLW, and MLLW are the mean higher high water, mean lower high water, mean sea level, mean higher low water, and mean lower low water, respectively.Table 2
Mean tidal levels in Chabahar [74].
MHHW (m)MLHW (m)MSL (m)MHLW (m)MLLW (m)2.531.931.611.280.69In the Chabahar area, maximum storm surges are caused by tropical cyclones. Based on numerical simulations of the Gonu storm that was performed as a part of the monitoring studies of Iranian coasts, the maximum storm surge was 40 cm in the western part of Chabahar Bay (Figure3) [75]. A similar conclusion was also drawn by Allahdadi et al. (2018) [76]. According to Figure 4, the maximum storm surge in the coastal area of Chabahar and the location of the seawall has reached 10 cm. In the present research, the amount of the storm surge was considered to be 35 cm, and other changes in sea level were considered to be 30 cm. Therefore, the total increase of the seawater level (excluding tidal changes) was set to 65 cm.Figure 3
Seawater level variations due to the Gonu storm [75].Figure 4
Seawater level variations in the Chabahar area due to the Gonu storm [75].


### 2.4. Wave Transformation from Deep Waters to the Nearshore Region

To obtain accurate reflection from the seawall, the incident wave properties at the toe of the seawall are required. As mentioned earlier, during the monitoring studies of Chabahar Bay, wave characteristics data were recorded in a deep-water location. In this study, the DHI MIKE 21 SW software is used to transform deep-water waves to the nearshore region.The DHI MIKE 21 SW model is a third-generation spectral model based on unstructured grids. This model can simulate the growth, decay, and transformation of wind waves in deep water and coastal areas [77]. This model has two types of formulations: directional decoupled parametric (DDP) and full spectral (FS). In this study, the latter one is used. In this model, the action density Nσ,θ is related to the energy density, Eσ,θ.(1)N=Eσ,where σ is the relative frequency and θ is the wave direction (normal to the wave crest of each spectral component). The evolution of the wave spectrum in the position x,y and time t is described by the following spectral action balance equation:(2)∂∂tN+∂∂xCxN+∂∂yCyN+∂∂σCσN+∂∂θCθN=Sσ,where the first term on the left-hand side of the equation is the local rate of change of action density in time, and the second and third terms are the propagation of the action in geographical space with propagation velocities Cx and Cy in x and y spaces, respectively. The fourth term represents shifting of the relative frequency due to variations in depths and currents (with propagation velocity Cσ in σ space). The fifth term represents depth-induced and current-induced refraction and propagation in directional space (with propagation velocity Cθ in θ space). The term S=Sσ,θ on the right-hand side of the action balance equation represents the superposition of source functions describing various physical phenomena:(3)S=Sin+Snl+Sds+Sbot+Ssurf,where Sin represents the generation of energy by wind, Snl is the wave energy transfer by nonlinear wave-wave interaction, Sds is the wave energy dissipation due to whitecapping, Sbot is the dissipation due to bottom friction, and Ssurf is the dissipated energy due to the depth-induced breaking [78].To transform waves from deep water to the toe of the structure, the bathymetry data, the geometry of the model boundaries, the water level data (fixed or variable), the flow condition data (fixed or variable), the wind data, the initial conditions (stagnation condition and condition of using experimental spectra), and the boundary conditions at model boundaries have been used as input data in the SW model. As shown in Figure5, to achieve the highest possible accuracy, an unstructured triangulated mesh was generated with varying sizes of elements with finer triangles on the nearshore area and coarser triangles on the offshore area. To increase the accuracy of the model, the size of grids was reduced in several steps. In this simulation, the dimensions of the mesh at the toe of the seawall were about 10–20 m, which is a suitable mesh size for simulating the transformation of waves to the toe of the structure. To introduce the depth in the mesh file, an attempt was made to use the various available depth measurement data. The hydrographic file used in front of the seawall is obtained from hydrographic operations performed in recent years. For offshore areas, the depth data of the ETOPO1 database [79] are used. Figure 5 also shows the water depth in nearshore and offshore regions.Figure 5
Computational domain, grids, and water depth adjacent and far from the structure.To perform the simulations in the full spectral (FS) mode, the spectral space was divided into 25 frequency spaces and 32 directional segments. To apply the boundary conditions at the open boundaries of the model, the height, period, and wave direction were considered 7.2 m, 13 s, and 160°, respectively. Also, to make the simulation conditions more compatible with the actual conditions, the wind effect was considered by applying a wind speed of 20 m/s, whose direction was the same as the incoming waves at the open boundaries. The whitecapping coefficient was set to 0.2. After setting up the numerical model, the execution time of the model was adjusted, so that waves reached a stable state. Since the slope was relatively steep, a constant breaking coefficient of 0.8 was used to obtain the most critical wave conditions at the toe of the structure.Figure6 shows the wave height contour at the location of the sea wall in the condition of maximum water level. As this figure shows, the wave height at a distance of 300 m from the seawall where the critical section of the wall is located is 4.5 m. Instead of using irregular waves that require large runtimes, in this study, the significant wave is used, which is a well-defined and standardized method to denote the characteristics of the random waves in a sea state. The significant wave height, the peak period, and the mean period of the wave are 4.5 m, 12.9 s, and 11.15 s, respectively. To perform SPH simulations, these values will be downscaled (at a scale of 1/30) using the Froude similitude that is commonly used for wave-seawall interaction problems [80].Figure 6
Significant wave height adjacent to the structure.


## 2.1. Study Area

The seawall is located at the Shahid Beheshti port of Chabahar (25° 17′ N, 60° 37′ E), Iran’s only oceanic port and the closest point to the Indian Ocean. Chabahar is a focal point of development in the southeast of Iran by enhancing transit routes among Central Asian countries and the northern part of the Indian Ocean [70]. It has the highest draft (16 meters) compared with other ports of Iran. With an area of over 4000 square meters and the capacity for simultaneous reception of more than 600 passengers, Chabahar has one of the largest international passenger terminals in the region [71]. To develop the Chabahar port, the construction of a seawall was planned at the right-hand side of the port. The location of this seawall is shown in Figure 1.Figure 1
Aerial photo from Chabahar port and the location of the seawall.


## 2.2. Hydrodynamic Data at the Location of the Sea Wall

To determine the characteristics of deep-water waves in the study area, the wave characteristics obtained from a comprehensive numerical hindcast of wave conditions in the Gulf of Oman for twenty-three years (1985–2007) are used [72]. The hindcast showed that in this period, the maximum wave height in deep water was about 4.3 m with a southwest direction. Meanwhile, the direction of the dominant waves was the south.In addition to local winds, the area may also be affected by tropical cyclones that generate relatively large sea states. At the Chabahar port, tropical cyclone waves are the dominant design wave condition as the wave heights can exceed those generated from other sources like the monsoon and local westerly winds [73]. Therefore, statistics and information considering rare but severe events should also be considered to determine the design waves.By performing statistical analysis on the tropical storm historical data from various sources, reviewing previous studies on tropical storms, numerical simulation of waves caused by historical tropical storms, numerical simulation of waves generated by tropical storms by the Monte Carlo method, and finally, the statistical analysis of extreme value analysis of maximum wave height, Dibajnia and Allahyar (2007) provided the design wave in deep water for different locations along the Oman Sea [73]. For a return period of 100 years, these values are shown in Figure 2. For Chabahar port, the wave heights for different return periods are given in Table 1.Figure 2
Spatial variations in the maximum significant wave height with a return period of 100-year in the Oman Sea [73].Table 1
Variation in maximum significant wave height by the return period in Chabahar [73].
LongitudeLatitudeReturn period (years)10 (m)25 (m)50 (m)100 (m)60.525.24.25.87.28.8


## 2.3. Sea-Level Variations

Sea-level changes should be considered in the seawall design as well as the wave transformation from offshore areas to the nearshore zone. Generally, sea-level changes arise from various factors such as tides and storm surges. Based on the report of the Integrated Coastal Zone Management (ICZM) plan of Iran [74], different tidal levels at the Chabahar port are given in Table 2. In this table, MHHW, MLHW, MSL, MHLW, and MLLW are the mean higher high water, mean lower high water, mean sea level, mean higher low water, and mean lower low water, respectively.Table 2
Mean tidal levels in Chabahar [74].
MHHW (m)MLHW (m)MSL (m)MHLW (m)MLLW (m)2.531.931.611.280.69In the Chabahar area, maximum storm surges are caused by tropical cyclones. Based on numerical simulations of the Gonu storm that was performed as a part of the monitoring studies of Iranian coasts, the maximum storm surge was 40 cm in the western part of Chabahar Bay (Figure3) [75]. A similar conclusion was also drawn by Allahdadi et al. (2018) [76]. According to Figure 4, the maximum storm surge in the coastal area of Chabahar and the location of the seawall has reached 10 cm. In the present research, the amount of the storm surge was considered to be 35 cm, and other changes in sea level were considered to be 30 cm. Therefore, the total increase of the seawater level (excluding tidal changes) was set to 65 cm.Figure 3
Seawater level variations due to the Gonu storm [75].Figure 4
Seawater level variations in the Chabahar area due to the Gonu storm [75].


## 2.4. Wave Transformation from Deep Waters to the Nearshore Region

To obtain accurate reflection from the seawall, the incident wave properties at the toe of the seawall are required. As mentioned earlier, during the monitoring studies of Chabahar Bay, wave characteristics data were recorded in a deep-water location. In this study, the DHI MIKE 21 SW software is used to transform deep-water waves to the nearshore region.The DHI MIKE 21 SW model is a third-generation spectral model based on unstructured grids. This model can simulate the growth, decay, and transformation of wind waves in deep water and coastal areas [77]. This model has two types of formulations: directional decoupled parametric (DDP) and full spectral (FS). In this study, the latter one is used. In this model, the action density Nσ,θ is related to the energy density, Eσ,θ.(1)N=Eσ,where σ is the relative frequency and θ is the wave direction (normal to the wave crest of each spectral component). The evolution of the wave spectrum in the position x,y and time t is described by the following spectral action balance equation:(2)∂∂tN+∂∂xCxN+∂∂yCyN+∂∂σCσN+∂∂θCθN=Sσ,where the first term on the left-hand side of the equation is the local rate of change of action density in time, and the second and third terms are the propagation of the action in geographical space with propagation velocities Cx and Cy in x and y spaces, respectively. The fourth term represents shifting of the relative frequency due to variations in depths and currents (with propagation velocity Cσ in σ space). The fifth term represents depth-induced and current-induced refraction and propagation in directional space (with propagation velocity Cθ in θ space). The term S=Sσ,θ on the right-hand side of the action balance equation represents the superposition of source functions describing various physical phenomena:(3)S=Sin+Snl+Sds+Sbot+Ssurf,where Sin represents the generation of energy by wind, Snl is the wave energy transfer by nonlinear wave-wave interaction, Sds is the wave energy dissipation due to whitecapping, Sbot is the dissipation due to bottom friction, and Ssurf is the dissipated energy due to the depth-induced breaking [78].To transform waves from deep water to the toe of the structure, the bathymetry data, the geometry of the model boundaries, the water level data (fixed or variable), the flow condition data (fixed or variable), the wind data, the initial conditions (stagnation condition and condition of using experimental spectra), and the boundary conditions at model boundaries have been used as input data in the SW model. As shown in Figure5, to achieve the highest possible accuracy, an unstructured triangulated mesh was generated with varying sizes of elements with finer triangles on the nearshore area and coarser triangles on the offshore area. To increase the accuracy of the model, the size of grids was reduced in several steps. In this simulation, the dimensions of the mesh at the toe of the seawall were about 10–20 m, which is a suitable mesh size for simulating the transformation of waves to the toe of the structure. To introduce the depth in the mesh file, an attempt was made to use the various available depth measurement data. The hydrographic file used in front of the seawall is obtained from hydrographic operations performed in recent years. For offshore areas, the depth data of the ETOPO1 database [79] are used. Figure 5 also shows the water depth in nearshore and offshore regions.Figure 5
Computational domain, grids, and water depth adjacent and far from the structure.To perform the simulations in the full spectral (FS) mode, the spectral space was divided into 25 frequency spaces and 32 directional segments. To apply the boundary conditions at the open boundaries of the model, the height, period, and wave direction were considered 7.2 m, 13 s, and 160°, respectively. Also, to make the simulation conditions more compatible with the actual conditions, the wind effect was considered by applying a wind speed of 20 m/s, whose direction was the same as the incoming waves at the open boundaries. The whitecapping coefficient was set to 0.2. After setting up the numerical model, the execution time of the model was adjusted, so that waves reached a stable state. Since the slope was relatively steep, a constant breaking coefficient of 0.8 was used to obtain the most critical wave conditions at the toe of the structure.Figure6 shows the wave height contour at the location of the sea wall in the condition of maximum water level. As this figure shows, the wave height at a distance of 300 m from the seawall where the critical section of the wall is located is 4.5 m. Instead of using irregular waves that require large runtimes, in this study, the significant wave is used, which is a well-defined and standardized method to denote the characteristics of the random waves in a sea state. The significant wave height, the peak period, and the mean period of the wave are 4.5 m, 12.9 s, and 11.15 s, respectively. To perform SPH simulations, these values will be downscaled (at a scale of 1/30) using the Froude similitude that is commonly used for wave-seawall interaction problems [80].Figure 6
Significant wave height adjacent to the structure.


## 3. The SPH Method

The SPH formulation comprises two approximation steps, including kernel approximation and particle approximation [81]. Kernel approximation of a function is achieved by integrating the multiplication of the function and an interpolating kernel (W). The following integral interpolant equation is used to approximate any desired function Ax in x′:(4)Ax=∫Ax′Wx−x′,hdx′.In the particle approximation step, a finite number of particles are used to represent the entire system, and the continuous integral representation is replaced by a discretized form of summation over all the particles inside the compact support. A smoothing length is used to define the compact support domain:(5)Axa≈∑bAxbWxa−xb,hΔvb,where a is an individual particle, and Δvb=mb/ρb is the volume of the neighboring particle (b). The mass and density are denoted by m and ρ, respectively. Therefore, the above equation can be rewritten as(6)Axa=∑bAxbmbρbWxa−xb,h.Two kernel functions are provided in the DualSPHysics program. These include the cubic spline and Wendland kernels. Previous investigations show that the Wendland kernel is more computationally convenient and results in better numerical convergence at large neighbor numbers than the cubic spline kernel [82]. Hence, in this research, a Wendland kernel is used as the interpolation function:(7)Wx,h=αD1−q242q+1,0≤q≤2,where αD=7/4πh2 in 2D. The nondimensional distance between any two particles is defined by q=x/h, where x is the distance between particles a and b, and h represents the smoothing length.The momentum conservation equation in a continuum is given as follows:(8)dvdt=−1ρ∇P+g+Γ,where Γ refers to the dissipative terms, and g is the gravitational acceleration. Monaghan (1992) proposed using the artificial viscosity scheme [83](9)dνadt=−∑bmbPb+Paρb.ρa+Πab∇aWab+g,where Pk and ρk are the pressure and density of the particle k. The viscosity term Πab is given by the following equation:(10)Πab=−αcab¯μabρab¯,vab⋅rab<0,0,vab⋅rab>0.where rab=ra−rb, vab=va−vb, and rk and vk are the particle position and velocity, respectively. μab=hvab⋅rab/rab2+η2, cab¯=0.5ca+cb is the mean speed of sound, η2=0.01h2, and α is a tuning coefficient that is used to introduce the proper dissipation. For wave propagation problems, it is recommended to use α=0.01[56].While the density of each particle changes throughout the WCSPH simulation, its associated mass remains constant. The continuity equation is used to calculate these density changes, which has the following form in SPH:(11)dρadt=∑bmbVab∇aWab.In the WCSPH method, the fluid is treated as weakly compressible, and an equation of state is used to estimate the fluid pressure. To maintain the size of the time step at each moment at a reasonable value, the adjustment of the compressibility is performed to artificially lower the speed of sound [84]. Performing this adjustment limits the speed of sound to be at least ten times faster than the maximum fluid velocity, maintains the density fluctuations less than 1%, and does not introduce any significant deviations from an incompressible approach. The equation of state is written in the following form:(12)P=bρρ0γ−1,where γ=7, b=c02ρ0/γ, ρ0=1000kg/m3 is the reference density, and c0=cρ0=∂P/∂ρ|ρ0 is the speed of sound at the reference density.


### 3.1. Time Stepping

A second-order symplectic integration scheme is used in this study. This numerical integration algorithm preserves the geometrical properties such as energy time-reversal symmetry in the equations of motion. This characteristic that improves the long-term behavior of the solution, arises from the time-reversibility of this scheme in the absence of friction or viscous losses. The momentumνa, density ρa, and position xa are first written in the following form:(13)dνadt=Fa,dρadt=Da,dxadt=νa.First, the values of acceleration and density are approximated at the middle of the time step:(14)xan+1/2=xan+Δt2νan,ρan+1/2=ρan+Δt2Dan.Then,dνan+1/2/dt is used to calculate the modified velocity and the position of the particles at the end of the time step:(15)νan+1=νan+1/2+Δt2Fan+1/2,xan+1=xan+1/2+Δt2νan+1.Finally, the updated values ofνan+1 and xan+1 are used to calculate the corrected density dρan+1/dt=Dan+1.


### 3.2. Variable Time Step

In the symplectic scheme, the time step depends on the CFL condition, the forcing terms, and the viscous diffusion term. A variable time stepΔt is calculated from the following equations:(16)Δtf=CFL.minΔtf,Δtcv,Δtf=minah/fa,Δtcv=minahcs+maxbhνabxab/xab2+η2,where Δtf is based on the force per mass fa, and Δtcv combines the viscous and courant time step controls.


### 3.3. Boundary Conditions

In this study, the dynamic boundary condition is used. In this boundary condition, the same as the fluid particles, the dynamic boundary particles are governed by the conservation laws of mass and momentum. However, the dynamic boundary particles are kept fixed in their location or move according to the external function. When the distance between the fluid particle and boundary particles becomes less than2h, the density and the pressure of the affected particles increases. This in turn exerts a repulsive force on the fluid particles to maintain them within the domain [85].


### 3.4. Second-Order Wave Generation of Regular Waves

The nonhomogeneous boundary condition at the free surface leads to interactions between each possible pair of first-order wave components. Hence, the free surface of the water has a second-order effect known as the bound wave component. Also, the generated wave has another second-order effect known as the free wave component that arises from the first-order wavemaker motion and the boundary conditions at the wavemaker. These two components have the same frequency and different propagation velocities. Therefore, the wave field becomes spatially inhomogeneous. Hence, the free wave component is undesirable. An additional second-order bound effect can be added to the flap motion to prevent the generation of free wave components. The abovementioned process is recognized as “second-order steering” of the wavemaker motion [86]. The piston displacement for a second-order Stokes wave is obtained using the following equation [87]:(17)et=S02sinωt+δ+H232d.3coshkdsinh3kd−2m1sin2ωt+2δ.where et is the piston position, S0=Hsinhkd.coshkd+kd/2sinh2kd is the piston stroke, ω is the angular frequency of the wave, d is the water depth, H is the wave height, m1=H/S0, and δ is the phase shift.


### 3.5. Calculation of the Reflection Coefficients

To calculate the reflection coefficient from the impermeable seawall, the method of Goda and Suzuki (1976) [88] is used. As a result of the wave reflection from the seawall, the wave system becomes the superposition of several waves propagating in different directions. Since the wave frequency is constant, the waves propagating in each direction form a single train of progressive waves as a whole. If the amplitude of the superposed incident waves be aI and that of reflected waves be aR, the free surface will have the following form:(18)ηI=aIcoskx−ωt+εI,ηR=aRcoskx+ωt+εR,where ηI and ηR are the surface elevations of the incident and reflected waves, k=2π/L is the wavenumber, L is the wavelength, ω is the angular frequency of waves, εI is the phase angle of the incident wave, and εR is the phase angle of the reflected waves. If the surface elevations are recorded at two adjacent stations of x1 and x2=x1+Δl, the observed profiles of composite waves will have the following general form:(19)η1=ηI+ηRx=x1=A1cosωt+B1sinωt,η2=ηI+ηRx=x2=A2cosωt+B2sinωt,where(20)A1=a1cosϕI+aRcosϕR,B1=aIsinϕI−aRsinϕR,A2=aIcoskΔl+ϕI+aRcoskΔl+ϕR,B2=aIsinkΔl+ϕ1−aRsinkΔl+ϕR,ϕI=kx1+εI,ϕR=kx1+εR.Solving (21) yields the following estimations:(21)aI=12sinkΔl×A2−A1coskΔl−B1sinkΔl2+B2+A1sinkΔl−B1coskΔl2,aR=12sinkΔl×A2−A1coskΔl+B1sinkΔl2+B2+A1sinkΔl−B1coskΔl2.Actual wave profiles usually contain some higher harmonics. The use of the Fourier analysis enables to estimate of the amplitudes ofA1, B1, A2, and B2 for the fundamental frequency as well as for higher harmonics. The amplitudes of the incident and reflected waves, aI and aR, are then estimated by equation (22). This is the procedure to be taken for regular wave tests. Finally, the reflection coefficient can be calculated by the following equation:(22)Cr=aRaI.


## 3.1. Time Stepping

A second-order symplectic integration scheme is used in this study. This numerical integration algorithm preserves the geometrical properties such as energy time-reversal symmetry in the equations of motion. This characteristic that improves the long-term behavior of the solution, arises from the time-reversibility of this scheme in the absence of friction or viscous losses. The momentumνa, density ρa, and position xa are first written in the following form:(13)dνadt=Fa,dρadt=Da,dxadt=νa.First, the values of acceleration and density are approximated at the middle of the time step:(14)xan+1/2=xan+Δt2νan,ρan+1/2=ρan+Δt2Dan.Then,dνan+1/2/dt is used to calculate the modified velocity and the position of the particles at the end of the time step:(15)νan+1=νan+1/2+Δt2Fan+1/2,xan+1=xan+1/2+Δt2νan+1.Finally, the updated values ofνan+1 and xan+1 are used to calculate the corrected density dρan+1/dt=Dan+1.


## 3.2. Variable Time Step

In the symplectic scheme, the time step depends on the CFL condition, the forcing terms, and the viscous diffusion term. A variable time stepΔt is calculated from the following equations:(16)Δtf=CFL.minΔtf,Δtcv,Δtf=minah/fa,Δtcv=minahcs+maxbhνabxab/xab2+η2,where Δtf is based on the force per mass fa, and Δtcv combines the viscous and courant time step controls.


## 3.3. Boundary Conditions

In this study, the dynamic boundary condition is used. In this boundary condition, the same as the fluid particles, the dynamic boundary particles are governed by the conservation laws of mass and momentum. However, the dynamic boundary particles are kept fixed in their location or move according to the external function. When the distance between the fluid particle and boundary particles becomes less than2h, the density and the pressure of the affected particles increases. This in turn exerts a repulsive force on the fluid particles to maintain them within the domain [85].


## 3.4. Second-Order Wave Generation of Regular Waves

The nonhomogeneous boundary condition at the free surface leads to interactions between each possible pair of first-order wave components. Hence, the free surface of the water has a second-order effect known as the bound wave component. Also, the generated wave has another second-order effect known as the free wave component that arises from the first-order wavemaker motion and the boundary conditions at the wavemaker. These two components have the same frequency and different propagation velocities. Therefore, the wave field becomes spatially inhomogeneous. Hence, the free wave component is undesirable. An additional second-order bound effect can be added to the flap motion to prevent the generation of free wave components. The abovementioned process is recognized as “second-order steering” of the wavemaker motion [86]. The piston displacement for a second-order Stokes wave is obtained using the following equation [87]:(17)et=S02sinωt+δ+H232d.3coshkdsinh3kd−2m1sin2ωt+2δ.where et is the piston position, S0=Hsinhkd.coshkd+kd/2sinh2kd is the piston stroke, ω is the angular frequency of the wave, d is the water depth, H is the wave height, m1=H/S0, and δ is the phase shift.


## 3.5. Calculation of the Reflection Coefficients

To calculate the reflection coefficient from the impermeable seawall, the method of Goda and Suzuki (1976) [88] is used. As a result of the wave reflection from the seawall, the wave system becomes the superposition of several waves propagating in different directions. Since the wave frequency is constant, the waves propagating in each direction form a single train of progressive waves as a whole. If the amplitude of the superposed incident waves be aI and that of reflected waves be aR, the free surface will have the following form:(18)ηI=aIcoskx−ωt+εI,ηR=aRcoskx+ωt+εR,where ηI and ηR are the surface elevations of the incident and reflected waves, k=2π/L is the wavenumber, L is the wavelength, ω is the angular frequency of waves, εI is the phase angle of the incident wave, and εR is the phase angle of the reflected waves. If the surface elevations are recorded at two adjacent stations of x1 and x2=x1+Δl, the observed profiles of composite waves will have the following general form:(19)η1=ηI+ηRx=x1=A1cosωt+B1sinωt,η2=ηI+ηRx=x2=A2cosωt+B2sinωt,where(20)A1=a1cosϕI+aRcosϕR,B1=aIsinϕI−aRsinϕR,A2=aIcoskΔl+ϕI+aRcoskΔl+ϕR,B2=aIsinkΔl+ϕ1−aRsinkΔl+ϕR,ϕI=kx1+εI,ϕR=kx1+εR.Solving (21) yields the following estimations:(21)aI=12sinkΔl×A2−A1coskΔl−B1sinkΔl2+B2+A1sinkΔl−B1coskΔl2,aR=12sinkΔl×A2−A1coskΔl+B1sinkΔl2+B2+A1sinkΔl−B1coskΔl2.Actual wave profiles usually contain some higher harmonics. The use of the Fourier analysis enables to estimate of the amplitudes ofA1, B1, A2, and B2 for the fundamental frequency as well as for higher harmonics. The amplitudes of the incident and reflected waves, aI and aR, are then estimated by equation (22). This is the procedure to be taken for regular wave tests. Finally, the reflection coefficient can be calculated by the following equation:(22)Cr=aRaI.


## 4. Validation

The model is verified by comparing the water level fluctuations obtained by the present DualSPHysics model with available experimental data in the presence of a submerged breakwater [89]. The physical model tests were performed at the wave flume of the coastal and offshore engineering laboratory, Zhoushan Campus, Zhejiang University, China. The considered structure was a trapezoidal submerged breakwater placed at x = 32.0 from the left wall of the flume. The dimensions of the structure are shown in Figure 7. To compare the results of the present study with the experimental data, a regular wave with a height of 0.05 m and a period of 1.47 s is generated in a water depth of 0.5 m, and the water level fluctuations were compared at two points of X = 28.0 m and X = 33.54 m from the left wall of the flume, respectively. As Figures 8 and 9 show, there is generally a good agreement between the results of the present study with the available experimental results. Therefore, the model is accurate enough to be used to simulate wave interaction with the impermeable seawall of the Chabahar port.Figure 7
Experimental setup of the considered submerged breakwater [89].Figure 8
Comparison between the water level fluctuations obtained by experimental and the present study at the location ofX = 28.0 m (St 1).Figure 9
Comparison between the water level fluctuations obtained by experimental and the present study at the location ofX = 33.54 m (St 3).


## 5. Results

In this section, the effect of the wave characteristics and the slope of the sea wall on the water level and the reflection coefficient are studied. The geometry of the sea wall is shown in Figure10.Figure 10
Geometry of the impermeable seawall in prototype scale.To simulate the interaction of impermeable seawall with waves, the prototype was downscaled to a 1/30 scale. As shown in Figure11, the modeled seawall is placed after a ramp. Based on the Froude similitude, the primary model’s water depth, wave height, and period were considered to be 0.662 m, 0.15 m, and 2.35 s, respectively. To construct the model of the structure, the particle spacing of 0.005 m is used. This particle spacing resulted in 313875 particles.Figure 11
Initial condition.Although it is beyond the scope of this study, the wave overtopping from the sea is schematically shown in Figure12. It should be noted that the calculation of the wave reflection coefficient was performed based on the method of Guda and Suzuki (1976) [88].Figure 12
Overtopping from the seawall.


### 5.1. Effect of the Wave Height

In this section, the effect of the wave height on the reflection coefficient of the modeled impermeable seawall of the Chabahar port is studied. The water depth and wave period were 0.662 m and 2.35 s, respectively. Figures13 and 14 show that increasing the wave height increases the maximum water elevation in both X1 and X2, which are two wave gauges located at a distance of 1 m and 0.9 m upstream of the seawall, respectively. The harmonic analysis results on the water level variations are shown in Figure 15, where an almost direct relationship is found between the wave height and the wat surface variations at both wave gauges. Figure 16 shows that generally, the wave height has a limited effect on the reflection coefficient, such that the variations of the reflection coefficient are just between 0.6 and 0.75. A similar conclusion was also made by Mutray et al. (2006) [13]. However, increasing the wave height slightly decreases the wave reflection coefficient.Figure 13
Effect of the wave height on the water level variations at X1.Figure 14
Effect of the wave height on the water level variations at X2.Figure 15
Effect of the wave height on the height of the water surface elevations at two gauges.Figure 16
Effect of the wave height on the reflection coefficient.


### 5.2. Effect of the Wave Period

The dependency of the water level fluctuations and the reflection coefficient on the wave period are studied in this section. In this section, the water depth and wave height were 0.662 m and 0.15 m, respectively. As Figures17–19 show, the water level fluctuations at both wave gauges do not have a linear relationship with the wave period. First, the fluctuations increase by increasing the wave period. Then, from T = 2.0 s to T = 2.8 s, the water level fluctuations have an inverse relationship with the wave period. After that, a direct relationship between the wave period and the reflection coefficient is observed. Figure 20 shows that generally, increasing the wave period intensifies the wave reflection coefficient.Figure 17
Effect of the wave period on the water level variations at X1.Figure 18
Effect of the wave period on the water level variations at X2.Figure 19
Effect of the wave period on the height of the water surface elevations at two gauges.Figure 20
Effect of the wave period on the reflection coefficient.


### 5.3. Effect of the Slope

The effect of seawall slope on the water level changes, and the reflection coefficient is shown in this section. The water depth, wave height, and period were 0.662 m, 0.15 m, and 2.35 s, respectively. The water level fluctuations in X1 and X2 locations are shown in Figures21–23. Figure 23 shows that increasing the slope of the seawall decreases the water level fluctuations. Figure 24 shows an almost direct relationship between the wave reflection coefficient and the slope of the seawall. This discussion shows that although increasing the slope of the structure results in lower materials and improves the project’s economy, the maximum slope should be selected carefully to avoid high reflection coefficients and too large water surface fluctuations.Figure 21
Effect of the slope of the seawall on the water level variations at X1.Figure 22
Effect of the slope of the seawall on the water level variations at X2.Figure 23
Effect of the slope of the seawall on the height of the water surface elevations at two gauges.Figure 24
Effect of the slope on the reflection coefficient.


## 5.1. Effect of the Wave Height

In this section, the effect of the wave height on the reflection coefficient of the modeled impermeable seawall of the Chabahar port is studied. The water depth and wave period were 0.662 m and 2.35 s, respectively. Figures13 and 14 show that increasing the wave height increases the maximum water elevation in both X1 and X2, which are two wave gauges located at a distance of 1 m and 0.9 m upstream of the seawall, respectively. The harmonic analysis results on the water level variations are shown in Figure 15, where an almost direct relationship is found between the wave height and the wat surface variations at both wave gauges. Figure 16 shows that generally, the wave height has a limited effect on the reflection coefficient, such that the variations of the reflection coefficient are just between 0.6 and 0.75. A similar conclusion was also made by Mutray et al. (2006) [13]. However, increasing the wave height slightly decreases the wave reflection coefficient.Figure 13
Effect of the wave height on the water level variations at X1.Figure 14
Effect of the wave height on the water level variations at X2.Figure 15
Effect of the wave height on the height of the water surface elevations at two gauges.Figure 16
Effect of the wave height on the reflection coefficient.


## 5.2. Effect of the Wave Period

The dependency of the water level fluctuations and the reflection coefficient on the wave period are studied in this section. In this section, the water depth and wave height were 0.662 m and 0.15 m, respectively. As Figures17–19 show, the water level fluctuations at both wave gauges do not have a linear relationship with the wave period. First, the fluctuations increase by increasing the wave period. Then, from T = 2.0 s to T = 2.8 s, the water level fluctuations have an inverse relationship with the wave period. After that, a direct relationship between the wave period and the reflection coefficient is observed. Figure 20 shows that generally, increasing the wave period intensifies the wave reflection coefficient.Figure 17
Effect of the wave period on the water level variations at X1.Figure 18
Effect of the wave period on the water level variations at X2.Figure 19
Effect of the wave period on the height of the water surface elevations at two gauges.Figure 20
Effect of the wave period on the reflection coefficient.


## 5.3. Effect of the Slope

The effect of seawall slope on the water level changes, and the reflection coefficient is shown in this section. The water depth, wave height, and period were 0.662 m, 0.15 m, and 2.35 s, respectively. The water level fluctuations in X1 and X2 locations are shown in Figures21–23. Figure 23 shows that increasing the slope of the seawall decreases the water level fluctuations. Figure 24 shows an almost direct relationship between the wave reflection coefficient and the slope of the seawall. This discussion shows that although increasing the slope of the structure results in lower materials and improves the project’s economy, the maximum slope should be selected carefully to avoid high reflection coefficients and too large water surface fluctuations.Figure 21
Effect of the slope of the seawall on the water level variations at X1.Figure 22
Effect of the slope of the seawall on the water level variations at X2.Figure 23
Effect of the slope of the seawall on the height of the water surface elevations at two gauges.Figure 24
Effect of the slope on the reflection coefficient.


## 6. Conclusion

This study aimed to evaluate the effects of the wave height and period and the slope of the Chabahar port designed impermeable seawall on the reflection coefficient and the free surface variations near the seawall. First, the deep-water wave data were collected. Then, a MIKE 21-SW model was used to transform this deep-water wave data to the toe of the seawall. Using these values, an impermeable seawall was designed. Then, an open-source WCSPH model named DualSPHysics was used to simulate the wave interaction with the newly designed seawall. The WCSPH model was first validated by comparing the results of the water level fluctuations obtained in this study with available experimental data. Then, the model was used to simulate the regular wave interaction with the designed seawall. In light of this numerical study, the following conclusions can be obtained:(i)
The WCSPH method was truly capable of simulating the wave interaction with coastal structures. It was shown that both the phase and the height of the waves at the two stations were in accordance with the experimental results.(ii)
The wave height has a limited effect on the wave reflection coefficient such that the range of the variations of the wave reflection coefficient was just between 0.6 and 0.75. It was observed that the wave reflection coefficient decreases by increasing the wave height. However, an almost direct relationship was found between the wave height and water level fluctuations near the structure.(iii)
Generally, the wave period has a direct relationship with the wave reflection coefficient. However, a linear relationship was not found between the wave period and water level fluctuations near the structure.(iv)
Increasing the slope of the seawall slightly decreases the water level fluctuations near the wall. However, it can intensify the wave reflection coefficient. To avoid too much disturbance as well as erosion near the seawall, steep slopes should be avoided.



---

*Source: 7382416-2021-12-02.xml*