Aims and Scopes

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

where fix,t and fieqx,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

The speed of sound in this model is cs=c/3  where 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, ρ=∑fi⋅ 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=∑eifi⋅ 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=ρcs2=ρ/3. The relaxation time is related to the kinematics viscosity by υ = (2τ-1)(Δx)²/6Δt.

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;

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

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, wr=3-2r+1+4r+4r2/8, when 1 ≤ |r| ≤ 2, wr=5-2r--7+12r-4r2/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

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

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

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

Step 3. Interpolate the velocity at the Lagrangian point with

Step 4. Correct the body force with

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.

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

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

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].

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.

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.

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.