# Numerical Modeling of Wave-Current Flow around Cylinders Using an Enhanced Equilibrium Bhatnagar-Gross-Krook Scheme

Authors: Liming Xing; Haifei Liu; Yu Ding; Wei Huang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Flow around cylinders is a classic issue of fluid mechanics and it has great significance in engineering fields. In this study, a two-dimensional hydrodynamic lattice Boltzmann numerical model is proposed, coupling wave radiation stress, bed shear stress, and wind shear stress, which is able to simulate wave propagation of flow around cylinders. It is based on shallow water equations and a weight factor is applied for the force term. An enhanced equilibrium Bhatnagar-Gross-Krook (BGK) scheme is developed to treat the wave radiation stress term in collision step. This model is tested and verified by two cases: the first case is the flow around a single circular cylinder, where the flow is driven by current, wave, or both wave and current, respectively, and the second case is the solitary waves moving around cylinders. The results illustrate the correctness of this model, which could be used to analyze the detailed flow pattern around a cylinder.


---

## Body

## 1. Introduction

The phenomena of flow around cylinders, which represent blunt bodies, widely exist in aviation, mechanical, and environmental engineering. In recent years, an increasing number of problems about complex flow around cylinders have been raised with the development of coastal engineering projects. Therefore, this topic attracts much attention among researchers.Flow around cylinders is a classic and complicated problem. The cross section is contracted, the velocity increases, and the pressure decreases along the path when the flow encounters cylinders. The separation of the boundary layer is formed around cylinders due to the viscous force, which is called the flow around cylinders. Additionally, cylinders are non-streamline objects, which influence the characteristics of flow around cylinders by many factors, such as the Reynolds number, the surface roughness, the turbulence intensity, and the cylinder size. All these lead to the complexity of flow around cylinders. The wave is one of the most common movement forms in water, and it is worth studying wave motion in shipping, coastal, and ocean engineering. Therefore, the research of wave propagation around cylinders is complicated, but significant.With the development of the fluid mechanics theory and the continuous updating of computer equipment, computational fluid dynamics has been greatly developed and numerical simulation became an important tool in research. Saiki and Biringen [1] introduced a virtual boundary technique to simulate uniform flows around cylinders, and the oscillations caused by this method can be attenuated by high-order finite differences. Based on this, Lima E Silva et al. [2] proposed the physical virtual model in which this immersed boundary was represented with a finite number of Lagrangian points, distributed over the solid-fluid interface. Ofengeim and Drikakis [3] presented numerical research on the interaction of plane blast waves and a cylinder, revealing that the blast-wave duration significantly influenced the unsteady flow around the cylinder. Breuer [4] computed the turbulent flow around a cylinder (Re = 3900) via large eddy simulation. Meneghini et al. [5] used a fractional step method to simulate laminar flows between two cylinders. Hu et al. [6] built a fully nonlinear potential model based on a finite element method to investigate the wave motion around a moving cylinder, and it provided certain important features that were absent in the linear theory. Wu and Shu [7] proposed a local domain-free discretization method that is able to simulate flow around an oscillating cylinder easier due to its advantage of handling the boundary. Claus and Phillips [8] used spectral/hp element methods to study the flow around a confined cylinder. The nonconforming spectral element method and adaptive meshes method were tested by Hsu et al. [9], demonstrating its feasibility on curve surfaces of cylinder.The lattice Boltzmann method (LBM) is a promising numerical simulation method of recent decades. Compared to traditional methods, LBM has many advantages: the algorithm is simple; it can deal with complicated boundary conditions; and it is suitable for parallel processing. These superiorities lead to wide usage of LBM in many research fields. Ginzburg and D’Humieres [10] introduced a new kind of boundary conditions, improving the accuracy close to the quasianalytical reference solution. Jiménez-Hornero et al. [11] used LBM to simulate the turbulent flow structure in an open channel with the influence of vegetation. Liu et al. [12] established a two-dimensional multiblock lattice Boltzmann model for solute transport in shallow water flows. Based on the Chapman-Enskog process, Liu and Zhou [13] proposed a lattice Boltzmann model to simulate the wetting-drying front in shallow flows.At the same time, many scholars have investigated the flow around cylinders based on the LBM. However, most studies are related to the heat transfer around cylinders. Yan and Zu [14] presented a numerical strategy to handle curved and moving boundaries for simulating viscous fluid around a rotating isothermal cylinder with heat transfer. Rabienataj Darzi et al. [15] used the LBM to analyze mixed convection flow and heat transfer between two hot cylinders. However, up to now, there is no LBM model for wave-current flow around cylinders.In this study, considering wave-current interaction, a two-dimensional hydrodynamic numerical model is developed based on the LBM. The model couples three types of stresses, including wave radiation stress, wind shear stress, and bed shear stress. Meanwhile, an enhanced local equilibrium function is developed to treat the wave radiation stress. It is used to simulate the propagation of waves in the flow around cylinders, and then two classic examples are used for validation, which can provide characteristics of flow around cylinders.


## 2. Methodology

### 2.1. Governing Equations

The two-dimensional shallow water equations including the continuity equation and momentum equation can be written in a tensor form as(1)∂h∂t+∂huj∂xj=0,∂hui∂t+∂huiuj∂xj=-g∂∂xih22+ν∂2hui∂xj∂xj-gh¯∂Zb∂xj+∂Sij∂xj+Fi,where the subscripts i and j represent the space direction indices and the Einstein summation convention is used; xj represents the Cartesian coordinate, taking x, y, and z in turn; uj represents the velocity component which takes u and v corresponding to that in x and y and directions, respectively. h represents the water depth; t represents the time; ν represents the kinematic viscosity; Zb represents the bed height of the datum plane and Fi represents the force term and defined as(2)Fi=τwiρ-τbiρ,where τwi represents the wind shear stress and τbi represents the bed shear stress.Wave Radiation Stress
(
S
i
j
). Longuet-Higgins and Stewart [16] defined the difference between the time-average momentum value and the static water pressure on the water column per unit area, known as the wave radiation stress.In (3), the wave radiation stresses Sxx, Sxy, Syx, and Syy are determined via local wave parameters. The wave radiation stress along the direction of wave propagation is Sx=E(2Cg/C-1/2), and the lateral one is Sy=E(Cg/C-1/2), where E=1/8ρgHw2, C is wave velocity, Cg represents the group velocity, and Hw represents the wave height. The conversion is conducted in the Cartesian coordinate system [17]:(3)Sxx=Sxcos2⁡θ-Sysin2⁡θ,Syy=Sxsin2⁡θ-Sycos2⁡θ,Sxy=Syx=Sxsin⁡2θcos⁡θ-Sycos⁡θsin⁡θ,where θ represents the angle between the wave direction and the x-axis.Bed Shear Stress
(
τ
b
i
). Bed shear stress (τbi) is generated by the wave-current interaction in the i direction, calculated as follows [18]:(4)τbi=ρCbuiujuj+πρ8fwuwjuwjuwj+FBρπ2Cbfw1/2uwjuwjuwj,in which Cb represents the bed friction coefficient, which may be either constant or calculated from Cb=g/Cz2, where Cz represents the Chezy coefficient given based on the Manning coefficient nb,(5)Cz=h1/6nb;uwi represents the wave bottom frictional velocity; FB represents the wave-current influence factor, which is equal to 0.917 for the waves and currents are in the same direction, −0.1983 for perpendicular relation and 0.359 for other angles [19]; and fw represents the wave friction factor, which is from 0.006 to 0.001 in practice [20].Wind Shear Stress
(
τ
w
i
). Wind shear stress (τwi) is usually expressed as(6)τwi=ρaCwuwiuwjwwj,where ρa is the density of air; Cw is the resistance coefficient; and uwi is the component of the wind velocity in i direction.


### 2.2. Lattice Boltzmann Method

On account of the lattice Boltzmann method with a D2Q9 lattice, an enhanced equilibrium BGK Scheme is developed in this paper. The wave radiation stressSij is treated in local equilibrium function at collision step.The discrete evolution process in the LBM with the enhanced force term [12, 21] can be written as(7)fαX+eαΔt,t+Δt-fαX,t=-1τfα-fαeq-3Δtωαeαjgh¯e2∂Zb∂xj+ΔtFα,where the external force term can be written as(8)Fα=3ωα1e2eαiτwiρ-τbiρ,where ωα represents the weight factor: ωα=4/9 for α=0; ωα=1/9 for α=1,3,5,7; ωα=1/36 for α=2,4,6,8. fα represents the distribution function of particles; fαeq represents the local equilibrium distribution function; Δt represents the time step; τ represents the single relaxation time; and eα represents the velocity vector of a particle in the α link.For the D2Q9 lattice shown in Figure1, each particle moves one lattice at its direction. The velocity of each particle is defined by (9)eα=0,0α=0,ecos⁡α-1π4,sin⁡α-1π4α=1,3,5,7,2ecos⁡α-1π4,sin⁡α-1π4α=2,4,6,8,where e=Δx/Δt and Δx is the lattice size.Figure 1
D2Q9 lattice.An equilibrium distribution functionfαeq can be expressed as(10)fαeq=Aα+Bαeαiui+Cαeαieαjuiuj+Dαuiui.Therefore, the equilibrium distribution function can be written as (11)fαeq=A0+D0uiuiα=0,A¯+B¯eαiui+C¯eαieαjuiuj+D¯uiuiα=1,3,5,7,A~+B~eαiui+C~eαieαjuiuj+D~uiuiα=2,4,6,8,where there must be(12)A1=A3=A5=A7=A¯,A2=A4=A6=A8=A~due to symmetry.Moreover, the local equilibrium distribution function must satisfy the following three conditions:(13)∑αfαeqX,t=hX,t,∑αeαifαeqX,t=hX,tuiX,t,∑αeαieαifαeqX,t=12gh2X,tδij-Sij+hX,tuiX,tujX,t.Hence, the relations among A0, A¯, and A~ are(14)A0+4A¯+4A~=h,2e2A¯+4e2A~=12gh2-Sij,A¯=4A~.We can obtain(15)B¯=h3e2,C¯=h4e2,D¯=-h6e2,B~=h12e2,C~=h8e2,D~=-h24e2.Therefore, the enhanced equilibrium distribution function fαeq is (16)fαeq=h-5gh26e2+5Sij3e2-2h3e2uiui,α=0,gh26e2-Sij3e2+h3e2eαiui+h2e4eαieαjuiuj-h6e2uiui,α=1,3,5,7,gh224e2-Sij12e2+h12e2eαiui+h8e4eαieαjuiuj-h24e2uiui,α=2,4,6,8.


### 2.3. Recovery of Wave-Current Coupling Equations

The recover deductions are following the Chapman-Enskog procedure.Based on (7), assuming Δt is small, taking Taylor expansion in time and space around point (X,t) leads to(17)fαX+eαΔt,t+Δt=fαX,t+Δt∂∂t+eαj∂∂xjfαX,t+12Δt2∂∂t+eαj∂∂xj2fαX,t+oΔt2.From Chapman-Enskog expansion, we have(18)fα=fα0+Δtfα1+Δt2fα2+oΔt2.Substitution of (17) and (18) into (7), one can obtain(19)Δt∂∂t+eαj∂∂xjfα0+Δtfα1+Δt2fα2+12Δt2∂∂t+eαj∂∂xj2fα0+Δtfα1+Δt2fα2=-1τΔtfα1+Δt2fα2-3Δtωαeαjgh¯e2∂Zb∂xj+ΔtFα.To order Δt, it is(20)∂∂t+eαj∂∂xjfα0=-1τfα1-3ωαeαjgh¯e2∂Zb∂xj+Fα.To order Δt2, it is(21)∂∂t+eαj∂∂xjfα1+12∂∂t+eαj∂∂xj2fα0=-1τfα2.Substitution of (20) into (21), we have(22)1-12τ∂∂t+eαj∂∂xjfα1=-12∂∂t+eαj∂∂xj-3ωαeαjgh¯e2∂Zb∂xj+Fα-1τfα2.Taking ∑[(20)+Δt×(22)] about α provides(23)∂∂t∑αfα0+∂∂xj∑αeαjfα0=-ε112e2∂∂xj∑αeαjeαkFk.Taking ∑eαi[(20)+Δt×(22)] about α provides(24)∂∂t∑αeαjfα0+∂∂xj∑αeαjeαjfα0+Δt1-12τ∂∂xj∑αeαjeαjfα1=-gh¯∂Zb∂xj+Fiδij.According to the law of conservation of mass, we know (25)∑αfαX,t=∑αfαeqX,t.If the center-scheme for the force term is applied, evaluation of the other terms in the above equations using (13) and (25) simplifies (23) and (24) and obtains (26)∂h∂t+∂huj∂xj=0∂hui∂t+∂huihuj∂xj=-g∂∂xih22-∂∂xjΛij-gh¯∂Zb∂xj+∂Sij∂xj+Fi,with(27)Λij=Δt2τ2τ-1∑αeαieαjfα1≈-ν∂hui∂xi+∂huj∂xi.Substitution of (27) into (26) leads to the following equations which were referred to as wave-current coupling equations (1).


## 2.1. Governing Equations

The two-dimensional shallow water equations including the continuity equation and momentum equation can be written in a tensor form as(1)∂h∂t+∂huj∂xj=0,∂hui∂t+∂huiuj∂xj=-g∂∂xih22+ν∂2hui∂xj∂xj-gh¯∂Zb∂xj+∂Sij∂xj+Fi,where the subscripts i and j represent the space direction indices and the Einstein summation convention is used; xj represents the Cartesian coordinate, taking x, y, and z in turn; uj represents the velocity component which takes u and v corresponding to that in x and y and directions, respectively. h represents the water depth; t represents the time; ν represents the kinematic viscosity; Zb represents the bed height of the datum plane and Fi represents the force term and defined as(2)Fi=τwiρ-τbiρ,where τwi represents the wind shear stress and τbi represents the bed shear stress.Wave Radiation Stress
(
S
i
j
). Longuet-Higgins and Stewart [16] defined the difference between the time-average momentum value and the static water pressure on the water column per unit area, known as the wave radiation stress.In (3), the wave radiation stresses Sxx, Sxy, Syx, and Syy are determined via local wave parameters. The wave radiation stress along the direction of wave propagation is Sx=E(2Cg/C-1/2), and the lateral one is Sy=E(Cg/C-1/2), where E=1/8ρgHw2, C is wave velocity, Cg represents the group velocity, and Hw represents the wave height. The conversion is conducted in the Cartesian coordinate system [17]:(3)Sxx=Sxcos2⁡θ-Sysin2⁡θ,Syy=Sxsin2⁡θ-Sycos2⁡θ,Sxy=Syx=Sxsin⁡2θcos⁡θ-Sycos⁡θsin⁡θ,where θ represents the angle between the wave direction and the x-axis.Bed Shear Stress
(
τ
b
i
). Bed shear stress (τbi) is generated by the wave-current interaction in the i direction, calculated as follows [18]:(4)τbi=ρCbuiujuj+πρ8fwuwjuwjuwj+FBρπ2Cbfw1/2uwjuwjuwj,in which Cb represents the bed friction coefficient, which may be either constant or calculated from Cb=g/Cz2, where Cz represents the Chezy coefficient given based on the Manning coefficient nb,(5)Cz=h1/6nb;uwi represents the wave bottom frictional velocity; FB represents the wave-current influence factor, which is equal to 0.917 for the waves and currents are in the same direction, −0.1983 for perpendicular relation and 0.359 for other angles [19]; and fw represents the wave friction factor, which is from 0.006 to 0.001 in practice [20].Wind Shear Stress
(
τ
w
i
). Wind shear stress (τwi) is usually expressed as(6)τwi=ρaCwuwiuwjwwj,where ρa is the density of air; Cw is the resistance coefficient; and uwi is the component of the wind velocity in i direction.


## 2.2. Lattice Boltzmann Method

On account of the lattice Boltzmann method with a D2Q9 lattice, an enhanced equilibrium BGK Scheme is developed in this paper. The wave radiation stressSij is treated in local equilibrium function at collision step.The discrete evolution process in the LBM with the enhanced force term [12, 21] can be written as(7)fαX+eαΔt,t+Δt-fαX,t=-1τfα-fαeq-3Δtωαeαjgh¯e2∂Zb∂xj+ΔtFα,where the external force term can be written as(8)Fα=3ωα1e2eαiτwiρ-τbiρ,where ωα represents the weight factor: ωα=4/9 for α=0; ωα=1/9 for α=1,3,5,7; ωα=1/36 for α=2,4,6,8. fα represents the distribution function of particles; fαeq represents the local equilibrium distribution function; Δt represents the time step; τ represents the single relaxation time; and eα represents the velocity vector of a particle in the α link.For the D2Q9 lattice shown in Figure1, each particle moves one lattice at its direction. The velocity of each particle is defined by (9)eα=0,0α=0,ecos⁡α-1π4,sin⁡α-1π4α=1,3,5,7,2ecos⁡α-1π4,sin⁡α-1π4α=2,4,6,8,where e=Δx/Δt and Δx is the lattice size.Figure 1
D2Q9 lattice.An equilibrium distribution functionfαeq can be expressed as(10)fαeq=Aα+Bαeαiui+Cαeαieαjuiuj+Dαuiui.Therefore, the equilibrium distribution function can be written as (11)fαeq=A0+D0uiuiα=0,A¯+B¯eαiui+C¯eαieαjuiuj+D¯uiuiα=1,3,5,7,A~+B~eαiui+C~eαieαjuiuj+D~uiuiα=2,4,6,8,where there must be(12)A1=A3=A5=A7=A¯,A2=A4=A6=A8=A~due to symmetry.Moreover, the local equilibrium distribution function must satisfy the following three conditions:(13)∑αfαeqX,t=hX,t,∑αeαifαeqX,t=hX,tuiX,t,∑αeαieαifαeqX,t=12gh2X,tδij-Sij+hX,tuiX,tujX,t.Hence, the relations among A0, A¯, and A~ are(14)A0+4A¯+4A~=h,2e2A¯+4e2A~=12gh2-Sij,A¯=4A~.We can obtain(15)B¯=h3e2,C¯=h4e2,D¯=-h6e2,B~=h12e2,C~=h8e2,D~=-h24e2.Therefore, the enhanced equilibrium distribution function fαeq is (16)fαeq=h-5gh26e2+5Sij3e2-2h3e2uiui,α=0,gh26e2-Sij3e2+h3e2eαiui+h2e4eαieαjuiuj-h6e2uiui,α=1,3,5,7,gh224e2-Sij12e2+h12e2eαiui+h8e4eαieαjuiuj-h24e2uiui,α=2,4,6,8.


## 2.3. Recovery of Wave-Current Coupling Equations

The recover deductions are following the Chapman-Enskog procedure.Based on (7), assuming Δt is small, taking Taylor expansion in time and space around point (X,t) leads to(17)fαX+eαΔt,t+Δt=fαX,t+Δt∂∂t+eαj∂∂xjfαX,t+12Δt2∂∂t+eαj∂∂xj2fαX,t+oΔt2.From Chapman-Enskog expansion, we have(18)fα=fα0+Δtfα1+Δt2fα2+oΔt2.Substitution of (17) and (18) into (7), one can obtain(19)Δt∂∂t+eαj∂∂xjfα0+Δtfα1+Δt2fα2+12Δt2∂∂t+eαj∂∂xj2fα0+Δtfα1+Δt2fα2=-1τΔtfα1+Δt2fα2-3Δtωαeαjgh¯e2∂Zb∂xj+ΔtFα.To order Δt, it is(20)∂∂t+eαj∂∂xjfα0=-1τfα1-3ωαeαjgh¯e2∂Zb∂xj+Fα.To order Δt2, it is(21)∂∂t+eαj∂∂xjfα1+12∂∂t+eαj∂∂xj2fα0=-1τfα2.Substitution of (20) into (21), we have(22)1-12τ∂∂t+eαj∂∂xjfα1=-12∂∂t+eαj∂∂xj-3ωαeαjgh¯e2∂Zb∂xj+Fα-1τfα2.Taking ∑[(20)+Δt×(22)] about α provides(23)∂∂t∑αfα0+∂∂xj∑αeαjfα0=-ε112e2∂∂xj∑αeαjeαkFk.Taking ∑eαi[(20)+Δt×(22)] about α provides(24)∂∂t∑αeαjfα0+∂∂xj∑αeαjeαjfα0+Δt1-12τ∂∂xj∑αeαjeαjfα1=-gh¯∂Zb∂xj+Fiδij.According to the law of conservation of mass, we know (25)∑αfαX,t=∑αfαeqX,t.If the center-scheme for the force term is applied, evaluation of the other terms in the above equations using (13) and (25) simplifies (23) and (24) and obtains (26)∂h∂t+∂huj∂xj=0∂hui∂t+∂huihuj∂xj=-g∂∂xih22-∂∂xjΛij-gh¯∂Zb∂xj+∂Sij∂xj+Fi,with(27)Λij=Δt2τ2τ-1∑αeαieαjfα1≈-ν∂hui∂xi+∂huj∂xi.Substitution of (27) into (26) leads to the following equations which were referred to as wave-current coupling equations (1).


## 3. Numerical Tests

### 3.1. Wave-Current Flow around a Circular Cylinder

This model is built based on the verified LBM hydrodynamic model [22]. The layout diagram of the channel is shown in Figure 2. The length is 7 m, and the width is 2 m. The bottom is flat and a solid cylinder with a 0.12 m radius is located at 2 m, 1 m. The initial water depth is 1 m and the flows go from the left to the right. The computational domain is divided by 140×40 computational grids. The time step is 0.01 s.Figure 2
The layout sketch of the channel.This case includes three different tests, which are driven by currents, by waves, and by both wave and current, respectively. The flow variables and wave parameters of three types situations are shown in Table1 (u0 is initial horizontal velocity and v0 is initial vertical velocity).Table 1
The flow variables and wave parameters.
Test
u
0 (m/s)
v
0 (m/s)
Wave period (s)
Wave amplitude (m)
1
1
0
-
-
2
0
0
0.5
0.1
3
1
0
0.5
0.1Test 1 (driven by the current). It can be seen that the water depth and flow velocity obviously varied due to the presence of the middle column (see Figure 3).Figure 3
x, y direction of the velocity (t=3s).
(a)
The magnitude ofu
(b)
The magnitude ofvWhen the flow encounters the cylinder, it passes around and a weak area emerges just behind the cylinder, where the circulation and a drop of water surface can be found.Test 2 (driven by the wave). The initial water is still and a wave maker is set at the inlet, where the incident waves are parallel generated in the x-axis. The water depth is intuitively depicted in Figure 4, where one can find regular wave propagation although there is a deformation caused by the cylinder.Figure 4
Three-dimensional water depth diagram (t=3s).In terms of the longitudinal velocityu, it is not always positive, as the flow is only driven by the wave (see Figures 4 and 5). This phenomenon is further described in Figure 6.Figure 5
Velocity vector diagram (t=3s).Figure 6
The magnitude ofu and v (t=3s).
(a)
u
(b)
vTest 3 (driven by both wave and current). Under the interaction effects of waves and currents, the wave run-up is pushed higher than before (see Figures 7 and 8), and the deformation process is more apparent (see Figure 9).Figure 7
The magnitude of water surface (t=3s).Figure 8
Three-dimensional water depth diagram (t=3s).Figure 9
The magnitude ofu and v (t=3s).
(a)
u
(b)
vTo illustrate the effects of currents and waves, the comparisons of the velocityu and the water depth h are plotted in Figures 10 and 11, respectively. It can be found that the wave-current interaction is not a simple superposition of waves and currents, and furthermore, wave-current interaction effects are greater than summation of these two effects separately.Figure 10
Comparison ofu (y=20, t=3s).Figure 11
Comparison ofh (y=20, t=3s).


### 3.2. Solitary Waves around Cylinders

This case is a classic cylinder model that has been simulated by many researchers before [23, 24]. In this section, a solitary wave around a cylinder is simulated first. The whole water channel is 60 m long and 30 m wide, and there is a circular cylinder with R=1.5875m in the center of the channel. The initial solitary wave with amplitude of 0.4 m is incident from left. Lattice size is 0.4 m, and the time step is 0.01 s.Figure12 shows the plots of three-dimensional perspective view of water surface at t = 8.7 s and 16 s. The solitary wave climbs up and a sequence of significant disperse waves after initial wave encountering the cylinder can be observed. At t=16s, the solitary wave is about to propagate out of the area. At the same time, disperse waves are fully developed to cover almost all the channel behind the frontal wave. The results are consistent with previous research [23].Figure 12
The plots of three-dimensional perspective view of water surface.
(a)
t = 8.7 s
(b)
t = 16 sFurthermore, a solitary wave around four cylinders is simulated. The simulation is conducted in an area of constant water depth (1 m), being 60 m long and 40 m wide. The distance between the centers of any two adjacent cylinders is 7.17 m, and the radius of four cylinders is the same with 2 m (see Figure13). The whole domain is divided by 150×100 computational grids. The time step is 0.01 s.Figure 13
The layout of the channel.Figure14 is the three-dimensional water depth of the wave around four cylinders at different times. The climbing up of water on the first cylinder can be observed at t = 7 s. At t = 9 s, the solitary wave encounters two middle cylinders and then runs up the front sides. Furthermore, a circular back disperse wave begins to turn up and propagates along the channel. The height of the middle part of the solitary wave decreases significantly due to obstruction of the frontal cylinder. The solitary wave encounters the rear cylinder at t = 10.5 s. The results show that the back disperse waves induced by the frontal cylinder form a circular wave pattern propagating towards the left open boundary. At the same time, the circular disperse waves, emerging from the two middle cylinders, are also expanded. Due to the complicated interactions between waves and cylinders, the diffracted wave patterns become fully irregular in the domain at t = 16 s. The results of the proposed model agree well with the work conducted by Zhong and Wang [23].Figure 14
The plots of three-dimensional perspective view of water surface.
(a)
t = 7 s
(b)
t = 9 s
(c)
t = 10.5 s
(d)
t = 16 s


## 3.1. Wave-Current Flow around a Circular Cylinder

This model is built based on the verified LBM hydrodynamic model [22]. The layout diagram of the channel is shown in Figure 2. The length is 7 m, and the width is 2 m. The bottom is flat and a solid cylinder with a 0.12 m radius is located at 2 m, 1 m. The initial water depth is 1 m and the flows go from the left to the right. The computational domain is divided by 140×40 computational grids. The time step is 0.01 s.Figure 2
The layout sketch of the channel.This case includes three different tests, which are driven by currents, by waves, and by both wave and current, respectively. The flow variables and wave parameters of three types situations are shown in Table1 (u0 is initial horizontal velocity and v0 is initial vertical velocity).Table 1
The flow variables and wave parameters.
Test
u
0 (m/s)
v
0 (m/s)
Wave period (s)
Wave amplitude (m)
1
1
0
-
-
2
0
0
0.5
0.1
3
1
0
0.5
0.1Test 1 (driven by the current). It can be seen that the water depth and flow velocity obviously varied due to the presence of the middle column (see Figure 3).Figure 3
x, y direction of the velocity (t=3s).
(a)
The magnitude ofu
(b)
The magnitude ofvWhen the flow encounters the cylinder, it passes around and a weak area emerges just behind the cylinder, where the circulation and a drop of water surface can be found.Test 2 (driven by the wave). The initial water is still and a wave maker is set at the inlet, where the incident waves are parallel generated in the x-axis. The water depth is intuitively depicted in Figure 4, where one can find regular wave propagation although there is a deformation caused by the cylinder.Figure 4
Three-dimensional water depth diagram (t=3s).In terms of the longitudinal velocityu, it is not always positive, as the flow is only driven by the wave (see Figures 4 and 5). This phenomenon is further described in Figure 6.Figure 5
Velocity vector diagram (t=3s).Figure 6
The magnitude ofu and v (t=3s).
(a)
u
(b)
vTest 3 (driven by both wave and current). Under the interaction effects of waves and currents, the wave run-up is pushed higher than before (see Figures 7 and 8), and the deformation process is more apparent (see Figure 9).Figure 7
The magnitude of water surface (t=3s).Figure 8
Three-dimensional water depth diagram (t=3s).Figure 9
The magnitude ofu and v (t=3s).
(a)
u
(b)
vTo illustrate the effects of currents and waves, the comparisons of the velocityu and the water depth h are plotted in Figures 10 and 11, respectively. It can be found that the wave-current interaction is not a simple superposition of waves and currents, and furthermore, wave-current interaction effects are greater than summation of these two effects separately.Figure 10
Comparison ofu (y=20, t=3s).Figure 11
Comparison ofh (y=20, t=3s).


## 3.2. Solitary Waves around Cylinders

This case is a classic cylinder model that has been simulated by many researchers before [23, 24]. In this section, a solitary wave around a cylinder is simulated first. The whole water channel is 60 m long and 30 m wide, and there is a circular cylinder with R=1.5875m in the center of the channel. The initial solitary wave with amplitude of 0.4 m is incident from left. Lattice size is 0.4 m, and the time step is 0.01 s.Figure12 shows the plots of three-dimensional perspective view of water surface at t = 8.7 s and 16 s. The solitary wave climbs up and a sequence of significant disperse waves after initial wave encountering the cylinder can be observed. At t=16s, the solitary wave is about to propagate out of the area. At the same time, disperse waves are fully developed to cover almost all the channel behind the frontal wave. The results are consistent with previous research [23].Figure 12
The plots of three-dimensional perspective view of water surface.
(a)
t = 8.7 s
(b)
t = 16 sFurthermore, a solitary wave around four cylinders is simulated. The simulation is conducted in an area of constant water depth (1 m), being 60 m long and 40 m wide. The distance between the centers of any two adjacent cylinders is 7.17 m, and the radius of four cylinders is the same with 2 m (see Figure13). The whole domain is divided by 150×100 computational grids. The time step is 0.01 s.Figure 13
The layout of the channel.Figure14 is the three-dimensional water depth of the wave around four cylinders at different times. The climbing up of water on the first cylinder can be observed at t = 7 s. At t = 9 s, the solitary wave encounters two middle cylinders and then runs up the front sides. Furthermore, a circular back disperse wave begins to turn up and propagates along the channel. The height of the middle part of the solitary wave decreases significantly due to obstruction of the frontal cylinder. The solitary wave encounters the rear cylinder at t = 10.5 s. The results show that the back disperse waves induced by the frontal cylinder form a circular wave pattern propagating towards the left open boundary. At the same time, the circular disperse waves, emerging from the two middle cylinders, are also expanded. Due to the complicated interactions between waves and cylinders, the diffracted wave patterns become fully irregular in the domain at t = 16 s. The results of the proposed model agree well with the work conducted by Zhong and Wang [23].Figure 14
The plots of three-dimensional perspective view of water surface.
(a)
t = 7 s
(b)
t = 9 s
(c)
t = 10.5 s
(d)
t = 16 s


## 4. Conclusion

This paper proposes a two-dimensional hydrodynamic model to investigate the wave-current interaction around cylinders. The lattice Boltzmann method was used to discretize the mathematical model in numerical simulation. A BKG scheme with an enhanced equilibrium is used to treat the wave radiation stress. The numerical results of both cases are in good agreement with practicalities and previous studies, demonstrating that this new model is able to produce reliable results for studying cylinders problems.



---

*Source: 8516879-2018-02-15.xml*