Journal of the Korean Society for Marine Environment and Energy

[ Original Article ]

Journal of the Korean Society for Marine Environment & Energy - Vol. 22, No. 4, pp.215-225

ISSN: 2288-0089 (Print) 2288-081X (Online)

Print publication date 25 Nov 2019

Received 08 Jul 2019 Revised 19 Aug 2019 Accepted 30 Oct 2019

DOI: https://doi.org/10.7846/JKOSMEE.2019.22.4.215

Numerical Simulation on a Sand Erosion under Wave in 2D by Lattice Boltzmann Method

Nyo Me Thet Naing¹

; Rho-Taek Jung²^{, †}

1Graduate Student, School of Naval Architecture and Ocean Engineering, University of Ulsan, Ulsan 44610, Korea
2Associate Professor, School of Naval Architecture and Ocean Engineering, University of Ulsan, Ulsan 44610, Korea

2차원 해양파내 모래침식을 위한 격자볼츠만 시뮬레이션

뇨메텟나잉¹

; 정노택²^{, †}

1울산대학교 조선해양공학부 대학원생
2울산대학교 조선해양공학부 부교수

Correspondence to: ^†rtjung@ulsan.ac.kr

Abstract

The phenomenon of sand erosion at slopping seabed induced by 2D nonlinear wave is studied numerically by using Lattice Boltzmann Method. Free surface motion is tracked by employing VOF (Volume of Fluid) single-phase free surface model and changes in seabed erosion is controlled by shear stress. Sediment transport and erosion mechanisms are predicted by stochastic process. In order to verify the code developed in this paper, wave generation is validated by wave theory and sediment model is verified by scour experiment of 2D cylinder. Moreover, when the nonlinear wave approaches the beach slope with vertical seawall at the end, it triggered scour and the comparison of numerical result with experiment values gave good agreement. It was found that this simulation is mainly related to the angle of repose of the sediments.

초록

수치 조파기(numerical wave maker)로 생성된 2차원 비선형파가 경사가 있는 모래바닥면을 침식시키는 양상에 대해서 격자 볼츠만 법(Lattice Boltzamann Method)을 이용하여 시뮬레이션 하였다. 자유수면의 형태변화는 VOF(Volume of Fluid)방식을 이용하였으며, 유체입자의 운동으로 인한 바닥면의 침식변화는 유체의 전단응력을 이용하고 침전물 수송 및 침식 메커니즘은 확률론적 과정(CA: Cellular Automata)을 포함시켰다. 본 논문에서 개발된 코드의 검증을 위해서 파이론에 의한 조파검증, 2차원 실린더의 세굴(scour)실험에 의한 검증을 수행하였다. 계속해서 본 논문은 비선형파가 모래경사면으로 진행시에 경사면 끝단에 위치한 수직구조물에 세굴이 발생하여 실험치와 유사하였고, 본 2차원 시뮬레이션은 퇴적물의 안식각(angle of repose)에 밀접한 영향을 받고 있음을 알 수 있었다.

Keywords:

Lattice boltzmann method, Wave maker, Immersed boundary, Scour, Seawall erosion

키워드:

격자볼츠만법, 조파기, 이멀스 경계, 세굴, 해양구조물 침식

1. Introduction

To predict sediment transport in coastal environment is of particularly importance for analyzing coast erosion and deposition accurately and solving the corresponding coast protection engineering problems further more. Coastal waves and nearshore currents are crucial hydrodynamics for coastal sediment transport. Scour results from any of the following hydrodynamics conditions, such as localized orbital velocity increases due to reflected waves, focusing of wave energy by structures that induces breaking, structure alignments that redirect currents and accelerate flows, flow constrictions that accelerate flow, downward directed breaking waves that mobilize sediment, flow separation and creation of vortices, transitions from hard bottom to erodible bed, wave pressure differentials and groundwater flow producing quick condition. The dynamics of sediment transport are complex under wave and current forcing.

The Numerical Wave Tank (NWT) which solves the Reynolds Averaged Navier-Stokes equation have proved to generate the real flow problems. Various kinds of wave makers, such as piston-type, flap-type wave makers, cylindrical wave makers, plunger wave makers or snake wave makers, can be deployed in the NWT. Navier-Stokes equation based internal mass-source functions was developed by Lin et al. [1999]. Wang et al. [2017] employed the mass source function in Lattice Boltzmann Method to generate waves. Zhao et al. [2012] proposed a method based on the combination of Lattice Boltzmann Equation for flow and Kinematic boundary condition for waves generation. Mao [1986] studied experimentally the critical condition of the onset of local scour in steady flow as pipeline model is directly placed on sand bed and described the role of vortices to initiate the scour. Sumer and Fredsøe [1992] carried out scour depth investigations for a comprehensive range of combined waves and currents. The prediction of scour formation under pipelines was done by Dupuis and Chopard [2002] by LBM D2Q7 model. The unique feature of this model is that it does not require any kind of empirical sediment transport formulae and the computer code can be parallelized with ease. Flower [1992] conducted a series of experiments to investigate seawall erosion under breaking waves. Ahmad et al. [2017] did the numerical modelling, using open source CFD model REEF3D, for breaking waves over a slope and the resulting erosion in the coastline.

As a recently developed CFD method, the Lattice Boltzmann Method has some advantages over conventional methods such as its simplicity, its scalability of parallel computation and the ease of handling complex geometries. It has been successfully applied to the analysis of a variety of complex physical phenomena, such as natural convection, turbulent flow, advection & dispersion problems, multi component flows and free surface flow problems. LBM simulations on only wave or only sediment model can be seen in literature so that this study is an effort of combining these phenomenon in order to estimate the flow problem around costal region. In this paper, the wave surface is treated by the single-phase free surface model from Thürey et al. [2007]. Wave generation is done by mass source function, Wang et al. [2017] and is also done by imposing immersed boundary at wave paddle. Scour formation model is estimated by employing Dupuis and Chopard [2002] sediment model.

2. Lattice Boltzmann Method

2.1 Single-Relaxation-Time Lattice Boltzmann Method (SRT LBM)

The Single-Relaxation-Time Bhatnagar-Gross-Krook (BKG) model without a force is given as follows: Bhatnagar et al. [1954].

f i x + e i Δ t, t + Δ t - f i x, t = - 1 τ f i x, t - f i e q x, t

(1)

where $f i x, t a n d f i e q x, t$ are the particle distribution function and equilibrium particle distribution function along the i direction respectively. The Δt is the advancing time step. The τ is the collision relaxation time parameter, and e_i is a discrete velocity vector. For two-dimensional problem considered in this work, the two-dimensional nine-velocity (D2Q9) lattice model is employed. The e_i= (0,0), (0, ±1), (±1, 0), (±1, ±1) for i = 0,..,8. The equilibrium distribution functions can be defined as

f i e q x, t = ρ w i 1 + e i . u c s 2 + 12 e i . u 2 c s 4 - 12 u 2 c s 2

(2)

The speed of sound in this model is $c s = c / 3 w h e r e c = ∆ x / ∆ t$ is the lattice speed with Δt and Δx being the discrete time step and lattice spacing respectively. The lattice speed c is set to be 1. The u is the fluid velocity. The w_i are weighting factors for different directions and define as w₀ = 4/9, w_1~4 = 1/9, w_5~8 = 1/36. The macroscopic density ρ of the fluid can be defined as a summation of 9 microscopic particle distribution function, $ρ = ∑ f i ⋅$ And the macroscopic velocity u which can be computed simply as an average of the microscopic velocities e_i weighted by the directional densities f_i, can be defined as, $ρ u = ∑ e i f i ⋅$ This equation transforms the discrete microscopic velocities of the LBM into a continuum of macroscopic velocities which represent the fluid’s motion. The pressure is related to density by $p = ρ c s 2 = ρ / 3$ . The relaxation time is related to the kinematics viscosity by υ = (2τ-1)(Δx)²/6Δt.

2.2 Immersed boundary method

The combination of IBM and LBM, called immersed-boundary-lattice-Boltzmann method (IB-LBM) was first proposed by Feng and Michelids [2005], well suited for the moving boundary simulations. There are several IBM variants, for example explicit or direct-forcing for rigid boundaries and explicit IBM for deformable boundaries. In this work, the multi-direct forcing method (MDFM) proposed by Wang et al. [2008] is applied. In order to include external force g(x,t) (body force, especially for the moving boundary), the evolution equation is split into the following two steps;

f i * x + e i Δ t, t + Δ t - f i x, t = - 1 τ f i x, t - f i e q x, t

(3)

f i x, t + Δ t = f i * x, t + Δ t + 3 w i ρ * x, t + Δ t e i ⋅ g x, t + Δ t Δ x

(4)

Supposing that f_i (x,t) and u(x,t) are known, the temporal f_i^*(x,t +Δt) and u^*(x,t + Δt) can be calculated. Let X_k(t + Δt) and U_k (t + Δt) (k = 1,...,N) be the Lagrangian points of the moving boundary and the boundary velocity at the points, respectively. The boundary velocity in here is the harmonic oscillatory motion of wave paddle. Note that the moving boundary is represented by N points, and the boundary Lagrangian points X_k generally differ from the Eulerian grid points x see Fig. 1. Then, the temporal velocities u^*(x,t + Δt) at the boundary Lagrangian points X_k are interpolated by

u * X k, t + Δ t = ∑ x u * x, t + Δ t W x - X k Δ x 2

(5)

Fig. 1.

Illustration of Immersed boundary method, blue points are Eulerian (mesh) points and green points are Lagrangian (moving boundary) points.

where $∑ x$ describes the summation over all lattice points x, W is a weighting function proposed by Peskin. The weighting function W is given by, W(x, y) = (1/ Δx)w(x / Δx) ⋅ (1 / Δy)w(y / Δy). Where w(r) can be defined depends on |r|, that is, when |r| ≤ 1, $w r = 3 - 2 r + 1 + 4 r + 4 r 2 / 8$ , when 1 ≤ |r| ≤ 2, $w r = 5 - 2 r - - 7 + 12 r - 4 r 2 / 8$ and w(r) = 0 elsewhere. The body force g(x,t + Δt) is determined by the following iterative procedure, Ota et al. [2011].

Step 0. Compute the initial value of the body force at the Lagrangian points by

g 0 X k, t + Δ t = s h U k t + Δ t - u * X k, t + Δ t Δ t

(6)

where it is noted that Sh/Δt = 1/Δx.

Step 1. Compute the body force at the Eulerian grid points of the l-th iteration by

g l x, t + Δ t = ∑ k = 1 N g l X k, t + Δ t W x - X k Δ S

(7)

where ΔS = Δx × s (s is the length of the Lagrangian grid).

Step 2. Correct the velocity at the Eulerian grid point by

u l x, t + Δ t = u * x, t + Δ t + g l x, t + Δ t Δ t / S h

(8)

Step 3. Interpolate the velocity at the Lagrangian point with

u l x, t + Δ t = u * x, t + Δ t + g l x, t + Δ t Δ t / S h

(9)

Step 4. Correct the body force with

g l + 1 X k, t + Δ t = g l X k, t + Δ t + S h U k t + Δ t - u l X k, t + Δ t Δ t

(10)

and go to step 1. From previous researches, g_l=5(x,t + Δt) is enough to keep the no-slip boundary condition at the boundary.

2.3 Mass source function

To include mass source S and body force F, equation (1) becomes, Wang et al. [2017].

f i x + e i Δ t, t + Δ t - f i x, t = - 1 τ f i x, t - f i e q x, t + Δ t F i

(11)

F i = 1 - 1 2 τ w i e i - u c s 2 + e i ⋅ u c s 4 e i F + e i e i 2 c s 2 + 2 S

(12)

Since the total mass of waves generated must be equal to the mass added by source function, the relation between desired wave velocity and source term in region Ω is

∫ 0 t ∫ Ω s x, y, t d Ω d t = ∫ 0 t ∫ 0 h u x, y, t d h d t = ∫ 0 t C η x, y d t

(13)

where, s(x, y, t) = S/ρ. u(x, y, t) is the horizontal particle velocity and C is the phase velocity of the desired wave. The detailed derivation can be seen in Wang et al. [2017].

3. Sediment Transport model

The sediment model mechanism is designed by implementing the concept of multi-particle cellular automata (CA). There are four rules to treat this model, namely, erosion, transportation, deposition, toppling, respectively. If the shields parameter or shear stress in the upper node is higher than the critical (threshold) shields parameter or threshold shear stress, then erosion condition is satisfied and the particle will be picked up and moved further away with probability p_erosion each particle belongs. Transportation of particles mainly depend on the combination of fluid velocity and falling speed at each lattice cell. The dispatching of sand particles is done stochastically. Particles less than N_thres are allowed to accumulate directly above the solid cell. But when N_thres limit is reached, the site solidifies and new incoming particles pile up on the site directly above. This solidification process implies a dynamically changing boundary for the fluid. In this model, particles do not have infinite cohesion, it is needed to consider the toppling rule, when a lattice site contains an excess amount of N_thres deposited particles with respect to its left or right neighbors, toppling rules is applied. The excess amount is calculated according to the value of angle of repose of pile. In real phenomenon, angle of repose depends on morphology of material, additions of solvents and coefficient of friction of material and so on. In case of sand, the values range from 45º (wet) to 34º (dry) but with water filled, it is between 15º to 30º. In contrast to real phenomenon, angle of repose in this study is simply the slope of sand particles collection in each lattice cells without accounting any cohesion or inter particle reaction. It is just the controlling parameter for toppling rule. So that for this simulation, the angle of repose value is set between 20º and 30º. The stable configuration may not be reached in one iteration and so that toppling rule is employed several times every 20 time steps of whole rules. This is a local rule and excess of sediments are not directly deposited on neighbor piles. They are actually directed in the direction of the piles. Moreover, erosion is favored when this local toppling process occurs. The bed shear velocity near the seabed is determined by considering the logarithmic velocity profile. The bed shear stress can be defined as, $τ 0 = ρ u * 2$ . Where u_* is the bed shear velocity and defined as, u / u_* = ln(30.2z / k_s) / κ. Where u is the velocity, κ = 0.4 is the von Karman constant, z is the vertical height above the seabed to nearest cell center, k_s = 3d₅₀ is the Nikuradse’s equivalent sand roughness, d₅₀ is the median grain size. The dimensionless form of the bed shear stress and its relationship to the sediment is named the Shields parameter, defined by θ_cr =τ_cr / g (ρ_s − ρ ) d. Where θ_cr is threshold Shields parameter, τ_cr is threshold bed shear stress, g is acceleration due to gravity, ρ_s is grain density, ρ is water density and d is grain diameter. The (critical) threshold bed shear stress calculated using the Shields diagram approach might lead to underestimation of the sediment transport because it does not account for the effects of the sloping bed. Dey [2003] proposed the modified critical shear formulation on sloping beds. The effect of a sloping bed is accounted for by considering the longitudinal bed slop θ, the transverse bed slope α, the angle of repose φ and drag and lift forces. The equation for the bed shear stress reduction factor r, when there is no angle of inclination of flow with respect to the longitudinal axis of the channel, is defined as follows;

r = 1 1 - η tan ⁡ φ tan ⁡ φ - sin ⁡ θ + η tan 2 ⁡ φ cos 2 ⁡ θ - sin 2 ⁡ α + sin ⁡ θ + η tan 2 ⁡ φ cos 2 ⁡ θ - sin 2 ⁡ α 2 + 1 - η 2 tan 2 ⁡ φ cos 2 ⁡ θ tan 2 ⁡ φ - sin 2 ⁡ α tan 2 ⁡ φ - sin 2 ⁡ θ - sin 2 α 0.5

(14)

where η is the ratio of lift force to drag force. For transverse bed slopes, θ = 0 and η = 0, so that equation (14) becomes, $r = cos ⁡ α 1 - tan 2 ⁡ α / tan 2 ⁡ φ$ . On the other hand, for longitudinal bed slopes, using α = 0, equation (14) becomes, r = cosθ (1− tanθ / tanϕ). The modified critical bed shear stress, τ_scr is calculated by multiplying the Shields critical bed shear stresses τ_cr with the reduction factor as, τ_scr = r·τ_scr.

4. Numerical Simulation

4.1 Wave generation

Firstly, the generation and propagation of linear wave with small amplitude on a constant water depth is simulated. The wave height H = 0.001m, wave length of L = 0.121 m in the water depth of h = 0.02 m with the wave period T = 1 s is generated. For intermediate water depth, kh is kept around 1.04, which is the same parameter as in Lin et al. [1999]. The computational domain is 0.5 × 0.022 m² and the grid size is set as 2500 × 110. Schematic diagram is shown in Fig. 2. In Fig. 2(a), the source region is at the left side of computational domain from the bottom to the free surface. In Fig. 2(b) and 2(c), the wave paddle is located inside the domain with some distance from the left boundary. Non-slip boundary condition is applied at the bottom and left boundaries. Open boundary condition is implemented at the top and right side of domain.

Fig. 2.

Schematic diagram for wave generation with (a) source function (b) flap-type (c) piston-type.

From Fig. 3 to Fig. 5, the phenomena of water particles velocity around wave paddle and source region at different time steps are displayed. Particles are traveling in the direction of wave propagation at leading part of wave crest while at the trough part; they are traveling in opposite direction. As a result, particles between crest and trough region are pushed downward whereas particles between trough and crest are forced upward. The snapshots shown here are for the linear waves.

Fig. 3.

The water particles velocity near flap type wave paddle motion at T/4, T/2, 3T/4 and T.

Fig. 4.

The water particles velocity near piston type wave paddle motion at T/4, T/2, 3T/4 and T.

Fig. 5.

The water particles velocity near source region (left of domain) at T/4, T/2, 3T/4 and T.

To validate the plane wave makers, several H/S ratios are simulated for both wave makers and investigated whether they agree with the wave maker theory proposed by Dean and Dalrymple [1991]. The comparisons can be found in Fig. 6 and a favorable agreement is achieved. Among all these tests, only one result by flap type is shown here. Fig. 7 depicts the comparisons at 2T, 3T, 4T, 5T, and 6T at the constant water depth of H/h = 0.2 with kh = 1.4.

Fig. 6.

Comparison of computed results for various H/S ratio with wave maker theory.

Fig. 7.

Wave generation by flap-type.

For weakly nonlinear long waves, the cnoidal wave theory is based on the shallow water approximation. A cnoidal wave trains of H/h = 0.2 and H/h = 0.3 are simulated and compared with the analytical result in Fig. 8. In case of H/h = 0.2, numerical results agree with analytical solution very well. As H/h approach to 0.3, the discrepancy between the numerical results and analytical solutions becomes significant, Fig. 8(d). This is not only due to the numerical inaccuracy but also due to the source function, which is based on the weak nonlinearity assumption of Boussinesq equations, is becoming less accurate for larger values of H/h.

Fig. 8.

Cnoidal wave generation.

4.2 Scour test

The numerical model is applied to the computation of the flow past a single circular cylinder on the sandy erodible bed of 1D deep. The flow conditions, circle diameter and the characteristic grain size resemble those of Mao’s experiment [1986]. In their experiment, they did with both vibrating and fixed pipe laid on the sand bed. In this study, the pipe is fixed. The Reynolds number is around 10⁴, so that this is treated as turbulent flow. The uniform sediment size of d₅₀ = 0.36 mm with Shields parameter of 0.1 is used. The diameter of the pipe is considered as 20 lattice units in numerical test. The computational domain is set as 1000 × 90. The schematic diagram is illustrated in Fig. 9. The free stream flow velocity at inlet is considered as 0.10 in lattice unit. The outlet is treated as open boundary and when sand particles reach the outlet, they will leave the computational domain. Then, probability of erosion p_erosion is chosen as 0.06 and N_thres is set as according to the sediment size. Since this is the turbulent flow with high Reynolds number, the Smagorinsky subgrid model is used with a constant C_smago = 0.4.

Fig. 9.

Computational domain with white and grey areas for fluid and sediment respectively.

Fig. 10.

Time evolution of the scour depth.

The scour profiles at 20000, 42000 and 135000 lattice time steps are shown in Fig. 11. They are fairly comparable with Mao’s test [1986] of 30, 80 and 217 minutes, respectively. Even though this model can approximately predict the scour profile of physical experiment, there is still drawback regarding development of scour hole with respect to real world time. In the study of Dupuis & Chopard [2002], it took much longer to achieve a characteristic time. This is because of the parameter they selected for simulation, such as falling velocity (gravitational force), p_erosion and number of threshold particles N_thres which is too small. Fig. 11. depicts the change of maximum depth with respect to lattice unit time for the entire period of scour progression. The empirical time development obtained from the relationship proposed by Sumer and Fredsoe [1992] is as, S_t = S(1−exp(−t/T)). The quantity S denotes maximum equilibrium scour depth and S_t is the scour depth at time t, the quantity T represents the time scale during which a substantial scour develops. Fig. 10 shows that the numerical model is able to capture the correct scouring rates at different stages of scour development. The scour rate is very high at the early stage of the time development of scour but at the later stage, it becomes much slower. This rate depends on the values chosen for the parameters such as u_fall, p_erosion and angle of repose and so on.

Fig. 11.

Comparison of bed profiles at different time steps.

Fig. 12. depicts the vortex shedding in the lee of circle at different significant stages of scour formation. At early stage of simulation, on the lee side of the circle, a pair of small attached vortex with large vortex behind them can be observed. The vortex shedding exits at the lee of the circle from the very beginning of scour development. It is believed that the continuous forming of vortex shedding is accounted for smearing out the ledge of deposits further downstream and changing the downstream slop of the scour hole in gentle nature.

Fig. 12.

Vortex shedding generation at the lee of the circle during scour development.

4.3 Sand erosion at seawall

To simulate the sediment erosion at the seawall, the cnoidal wave of wave height, H = 0.24 m, still water level h = 1.16 m with period T = 2 s is simulated. The bed material used for the sloping seabed consists of non-cohesive sand with median particle diameter d₅₀ = 0.13 mm. The computational domain is 46.66 × 2.67 m² and the grid size is set as 1400 × 80. In the left boundary, mass source function is applied to generate wave and the right boundary is placed as wall so that the fluid particles and sand particles will bounce back from this wall boundary. The schematic diagram is illustrated in Fig. 13. The flat bottom starts from the beginning x = 0 to 35m followed by a sloping seabed of m = 1:15. Then, N_thres is set as according to the sediment size.

Fig. 13.

Schematic diagram for sand erosion at seawall.

The parameter, p_erosion and repose angle ϕ are the crucial values to determine the erosion pattern during numerical simulation. Here, p_erosion of 0.03 and repose angle of 20º and 30º are used to compared the resultant eroded bed forms after maximum equilibrium scour depth is reached. The values of repose angle are selected as 20º and 30º because this value can be between 15º and 30º for water filled sand. And the degree of sloped seabed is merely around 3.814, so that very steep angle will give undesirable morphology of eroded bed. Fig. 14 shows that higher repose value gives steeper equilibrium depth and sediments are transported further offshore than smaller one. The eroded sediment mass is deposited approximately 1 m away from the seawall which creates a breaker bar. One of the contributing parameters in transported bed profile (morphology) is p_erosion, the erosion probability of each particle possesses. It is varied from types of materials (condition). Different values of these parameters will give slightly different morphology and faster or slower characteristic time to reach equilibrium condition.

Fig. 14.

Seabed profile sequence with repose angle of 20º and 30º at 40s and 120s without free surface line for simplicity.

The experimental data conducted by Fowler [1992] is used to validate the numerical result. The comparison of simulated and measured erosion profile is given in Fig. 15. In the experiment, the maximum scour depth at seawall is 0.151 m from the initial position. The result of p_erosion = 0.03 with repose angle of 20º and 30º both yield approximately 0.153 m of maximum equilibrium scour depth at seawall at t = 120s. Fig. 16 depicts the simulated bed profiles at 40s, 80s and 120s. At earlier time, the erosion rate is slow and it is gradually increase due to the high wave impacts on the vertical seawall which in turn creates velocities impacting towards seawall toe and mobilizes the sand particles from the toe. The eroded sediments are then carried offshore by the reflected waves and it continues until the equilibrium state is reached.

Fig. 15.

Comparison of simulated and measured erosion profile.

Fig. 16.

Simulated seabed profiles at 40s, 80s and 120s.

5. Conclusion

A code for simulating send-erosion under wave was developed and validated. The computational domain is designed as relatively small and only coarse grid system is used for all cases for this initial development of the code. Despite of using coarse grid system with single relaxation scheme, the results are fairly comparable with analytical solutions and experimental results. It is noted that when the wave height to water depth ratio is smaller than 0.05, the wave surface form is distorted and multigrid (multiblock) system or multi relaxation time scheme is needed for better outcomes. The important parameters for the scour prediction are number of threshold sand particles, probability of erosion each particle have and repose angle of sediment. Also, threshold shear stress gives big contribution in rate of erosion process. Moreover, very fine grain size is assumed to be used in all scour formation cases. So that influences given by porosity and permeability were not accounted. Each and every sediment particle is treated by CA method with discrete representation. Depending on the chosen parameters, such as p_erosion or repose angle, the evolution of the bed form may be slightly differing from each other. Nevertheless, the simulated results show favorable outcomes. The implementation of multiblock, installation of wave absorbing mechanisms and extending the domain to three-dimension is the future work of this study.

References

Ahmad, N., Bihs, H., Chella, M. A. and Arntsen, A., 2017, Numerical Modelling of Arctic Coastal Erosion due to Breaking Waves Impact Using REEF3D, Int. Ocean. Polar. Eng. Conf.
Dupuis, A. and Chopard, B., 2002, “Lattice gas modeling of scour formation under submarine pipelines,” Journal of Computational Physics, vol. 178, pp. 161-174. [https://doi.org/10.1006/jcph.2002.7025]
Bhatnagar, P.L., Gross, E.P. and Krook, M., 1954, “A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems”, Physical review, 94(3), 511. [https://doi.org/10.1103/PhysRev.94.511]
Dean, R.G. and Dalrymple, R.A, 1991, “Water Wave mechanics, for Engineers and Scientists”. World Scientific , River Edge, NJ. [https://doi.org/10.1142/9789812385512]
Dey, S., 2003, “Threshold of sediment motion on combined transverse and longitudinal sloping beds”, Journal of Hydraulic Research, 41(4), 405-415. [https://doi.org/10.1080/00221680309499985]
Feng, Z., Michaelides, E., 2004, “The immersed boundary-lattice boltzmann method for solving fluidparticles interaction problems”, Journal of Computational Physics 195(2), 602-628. [https://doi.org/10.1016/j.jcp.2003.10.013]
Fowler, J.E., 1992, “Scour problems and method of prediction of maximum scour at vertical seawalls”, Report, Department of the Army, Waterways Experiment Station, Corps of Engineers, Vicsburg, Mississippi U.S., Technical Report CERC-92-16.
Lin, P. and Liu, P.L.F., 1999, Internal wave-maker for navier-stokes equations models, Port, Coastal and Ocean Engineering, 125(4), 207-217. [https://doi.org/10.1061/(ASCE)0733-950X(1999)125:4(207)]
Mao, Y., 1986, “The Interaction Between a Pipeline and an Erodible Bed”, Technical Univ. of Denmark, Lyngby, Denmark.
Ota, K., Suzuki, K. and Inamuro, T., 2011, Lift generation by a two-dimensional symmetric flapping wing: immersed boundary-lattice Boltzmann simulations Fluid Dyn. Res. 44 045504. [https://doi.org/10.1088/0169-5983/44/4/045504]
Sumer, B.M. and Fredsøe, J., 1992, “A review of wave/current induced scour around pipelines”, 23^rd International Coastal Engineering Conference, vol.3, pp 2839-2852, Venice, Italy.
Thürey, N., 2007, “Physically based animated of free surface flows with the lattice boltzmann method”, Ph.D Thesis, University of Erlangen-Nuremberg, Germany.
Wang, K., Yu, Y., Yang, L. and Hou, G., 2017, Lattice Boltzmann based internal wave-maker, Int. Ocean and Polar. Eng. Conf, San Fransico, CA, USA.
Wang, Z., Fan, J. and Luo, K., 2008, “Combined multi-direct forcing and immersed boundary for simulating flows with moving particels”, Int. J. Multiph. Flow 34, 283-302. [https://doi.org/10.1016/j.ijmultiphaseflow.2007.10.004]
Zhao, Z., Huang, P., Li, Y. and Li, J., 2012, A lattice Boltzmann method for viscous free surface waves in two dimensions, Int. J. Numer. Meth. Fluids. [https://doi.org/10.1002/fld.3660]