Simulation of Fluid and Structure Interface with Immersed Boundary–Lattice Boltzmann Method Involving Turbulence Models

The multiple-relaxation-time (MRT) version of the immersed boundary–lattice Boltzmann (IB-LB) method is developed to simulate fluid-structure interfaces. The innovations include the implicit velocity correction to ensure no-slip boundary conditions and the incorporated Smagorinsky’s algebraic eddy viscosity for simulating turbulent flows. Both straight and curved interfaces are investigated. The streamlines penetration can be well prevented, which means the no-slip boundary condition can be guaranteed. Due to the existence of two coordinate systems: the Lagrangian coordinate system and the Eulerian coordinate system, the velocity and force properties on the structure can be easily calculated. Several benchmark simulation cases are carried out to verify the correctness of the model, including flow around circular cylinder at Re = 20, 150, and 3900 and flow around square cylinder at Re = 150 and 1000. The results agree well with previous studies, especially in the events of lower Reynolds numbers. Due to the three-dimensional turbulence vortex effects, the discrepancy increases are associated with higher Reynolds numbers. In addition, the effect of rotating velocity on the interaction process of the square cylinder in flows is researched, coupled with the capability of dealing with the moving boundaries.


Introduction
Flow around bluff bodies is a classical and fundamental research topic in many engineering fields.Besides the experiments and theoretical studies, the topic has been widely studied by using various numerical methods, like the finite difference method, the finite element method, and the boundary element method.Most of the simulations are focused on the effective and efficient treatments of fluid-structure interfaces [1][2][3][4].
As a newly developed and potential numerical method, the LB method has obtained great successes in the areas of multiphase flows, porous flows, and turbulence flows [5][6][7].Due to the background of kinetic theory, the LB method provides a concise and efficient microscopic way dealing with complex fluid-structure interfaces [2,3,8].Generally, the popular LB schemes for boundary treatments include the traditional bounce-back scheme, the half-way bounce-back scheme, and the extrapolation scheme.Some novel strategies are also proposed, combining the advantages of other numerical methods among which the immersed boundary (IB) method is quite glamorous [9][10][11][12].As a great variability combining the IB method and the LB method, the IB-LB method owns two coordinate systems: the Eulerian coordinate system and the Lagrangian coordinate system.The Eulerian coordinate system is applied to simulate the fluid dynamics by solving the Navier-Stoke equations, and the Lagrangian coordinate system is introduced to describe the structure movements.Due to the separation strategy, the fluid-structure interaction problem can be investigated flexibly.However, Wu et al. proposed that the no-slip boundary condition could not be well guaranteed in the conventional IB-LB models [12][13][14].Then they developed a novel IB-LB model to enforce the no-slip boundary condition based on the Bhatnagar-Gross-Krook (BGK) scheme and the body force term proposed by Guo et al. [15] which can consider the effect of external forces on the momentum and momentum flux as well as the discrete lattice effect.The key issue in the IB-LB model is the implicit velocity correction process.
However, the IB-LB model proposed by Wu and Shu is based on the BGK scheme which contains only one singlerelaxation-time parameter (), so the inherited disadvantage is that the choice of the relaxation parameter must be high enough to keep the numerical stability, leading to the simulated kinematic viscosity (V) to be higher according to V = ( − 0.5)/3.However, as an advanced version of the single-relaxation-time (SRT) scheme, the simulated kinematic viscosity in the MRT scheme can be much smaller than that in the SRT.Additionally, the works by Guo and Zheng [16] show the advantage of the MRT compared with SRT in eliminating the numerical boundary slip.Then a novel immersed boundary-lattice Boltzmann model based on the MRT was proposed by Lu et al. [17], in which the numerical boundary slip reduced dramatically.In many real engineering fields, the work environment owns the attributes of low kinematic viscosity and high Reynolds number, Re = /V, where  is the characteristic length and  is the characteristic velocity.High Reynolds numbers are commonly accompanied by the turbulent flow containing eddies of a large spectrum of sizes [18].Particularly, in order to get sufficient information in simulating turbulent flows, quite a lot of computational resources are needed.The number of computational nodes necessary to resolve three-dimensional turbulent flow scales as (Re 9/4 ) = (() 9/4 V −9/4 ), and the number of compute nodes can be reduced by letting V → 0 [18].Therefore, based on the essence of the IB-LB model proposed by Wu and Shu, in which the no-slip boundary condition can be enforced by an implicit velocity correction, the MRT-IB-LB scheme is developed to extend the IB-LB model to be suitable for the turbulent flows.What is more, the Smagorinsky's algebraic eddy viscosity approach [19] is incorporated with the MRT-IB-LB scheme to improve the capability of simulating turbulent flows.
The rest of the paper is organized as follows.The MRT scheme of the IB-LB method with implicit velocity corrections is presented in Section 2. Then the developed MRT-IB-LB method in conjunction with the Smagorinsky subgrid scale (SGS) model [7] is introduced.In Section 3, numerical experiments are completed to test the performance of the present model.In addition, a rotating square cylinder in the laminar flow is simulated to confirm the inherited ability of IB method in dealing with moving boundaries.Finally, some conclusions are drawn in Section 4.

Numerical Implementation
2.1.MRT Lattice Boltzmann Equation.Different from the traditional numerical methods in which the governing equations are discretized, the LB method comes from the thought that the flows consist of particles, and the flow behavior is a statistical result of the particles.The motion of the particles can be expressed as the Boltzmann equation at molecular scales [20,21]: where  = (r, u, ) is the particle distribution function at the position r, u is the particle velocity, and  is the time variable.Ω() is the collision term containing various and complex collision forms between particles.Due to the fact that the direct calculation of Ω() is quite difficult, Bhatnagar et al. [22] proposed an effective collision operator: which is also named as the BGK collision term. is the relaxation time parameter and  eq is the equilibrium distribution function.In the LB simulation, ( 2) is taken into (1), and then the modified equation ( 1) is discretized in the velocity space with a finite set of velocity vectors {e  } where  is the discretized direction.The completely discretized equation, that is, LBGK equation, can be written as where   stands for the particle distribution associated with the th discretized direction.In present study, the D2Q9 model is used, so  ranges from 0 to 8. The equilibrium distribution function can be written as [23,24] where   is the weighting factors,  0 = 4/9,  1∼4 = 1/9,  5∼8 = 1/36, and  is the particle velocity,  = /.With the content of conservation laws, and the condition that the Knudsen number is Kn ≪ 1, the Navier-Stokes equations can be derived from the LB equation by using the multiscale Chapman-Enskog expression.
The main difference between the MRT scheme and the SRT scheme is that the single-relaxation-time parameter in (3) is replaced by a nonnegative 9 × 9 diagonal relaxation matrix: Correspondingly, the collision term on the right hand of (3) is transformed into the moment space [25]: where m is the velocity moment of the distribution function and  is the transformation matrix.
The moments m can be written as m = (, , ,   ,   ,   ,   ,   ,   )  , where  is the density,   =   , and   =   . and  are related to the energy and their equilibria can be written as , is related to energy flux, and   and   are related to the diagonal and off-diagonal components of the stress tensor.Their equilibria are, respectively, defined as In the present study, the relaxation parameters in Ŝ are valued as  0 =  3 =  5 = 1,  1 = 1.9,  2 = 1.54, and  4 =  6 = 1.9. 7 ,  8 are related to 1/ and the eddy viscosity.

MRT-IB Model with Implicit Velocity Correction.
In the IB method, the fluid-structure interaction can be simplified as the process that the force acting on the structure is spread on the nearby fluid nodes through a Dirac delta function ().On the contrary, the velocity of the structure nodes is distributed based on the surrounding fluid nodes.The interactive processing can be formulated as [10,12,13,26] where  and r are the coordinate values in the Lagrangian coordinate system and the Eulerian coordinate system, respectively.R(, ) is the Eulerian coordinate value converted from the Lagrangian coordinate value of the structure nodes.
The curve Γ stands for the immersed boundary in the computational domain Ω.F(, ) and f(r, ) are the force densities on the structure nodes and the fluid nodes.u(R(, ), ) and u(r, ) are the velocities of the structure nodes and the fluid nodes.The basic idea of IB method is that the presence of structures acts approximately as a body force on fluid motions.From ( 13), it can be clearly seen that, for a specific Dirac delta function (), the results of f(r, ) are completely dependent on calculating F(, ).However in the conventional IB-LB models, the force density is explicitly calculated in advance, and the no-slip boundary condition cannot be completely enforced [12].Thus, the idea considering the velocity correction at all structure nodes is adopted here, associated with the discretized body force form proposed by Guo et al. [15], which can be written as S is the body force.Correspondingly, the evolution equation containing the discretized body force can be modified as Combining ( 3) and ( 7), the above evolution equation can be written in the MRT scheme as And the velocity can be computed by The velocity u in ( 18) can be divided into two parts, u = û+Δu, where û = ∑   e  is the intermediate velocity and Δu = SΔ/2 is the velocity correction.In the conventional IB-LB schemes, S is explicitly calculated, and then the velocity correction can be decided in advance.Here, based on the implicit velocity correction, S is kept unknown at first and is determined by the velocity correction at the structure nodes satisfying the noslip boundary condition.Assuming the unknown velocity correction at the structure point is Δu  in the Lagrangian coordinate system, the velocity correction Δu at the fluid point in the Eulerian coordinate system can be calculated by With the continuous kernel distribution   [12], (19) can be rewritten as where Δ is the arc length of the structure boundary element.
Therefore, the corrected velocity can be expressed as Because the intermediate velocity û(r  , ) is obtained in advance, the only unknown variable in (22) is Δu  (R(, ), ), which is related to the velocity of the structure boundary points under the limit of the no-slip boundary condition.In other words, the velocity at the structure points must be equal to the velocity at the points obtained by interpolation using the smooth delta function: According to (22) and ( 23), we can see that Δu  (R(, ), ) is still the only unknown variable, but it can be calculated directly now associated with the specific boundary velocity u  (R(, ), ).After Δu  (R(, ), ) is calculated, the velocity correction at the fluid point can be obtained.Then the corrected velocity and force density can be all obtained.
Particularly, the discretized body force form in ( 15) is prepared for the SRT scheme and it must be modified before it is used in the MRT-IB model.Following the basic idea of MRT, the body force term is transformed into the moment space [6,27]: where  =  −1 Ŝ. 0∼8 is the force term whose formation in the moment space can be rewritten as  =  with where  , and  , are the components of macro velocity and macro force, so (17) can be rewritten as

Incorporating Eddy Viscosity into MRT-IB Model.
In accord with the fact that the turbulent flows exist in many real engineering fields, the Smagorinsky's algebraic eddy viscosity approach is incorporated with the MRT-IB-LB model to extend its ability to deal with the turbulent flow containing eddies of a large spectrum of sizes.In the LB simulation, the viscosity is only related to the relaxation time parameter  which is located in the nonnegative 9 × 9 diagonal relaxation matrix; that is,  7 =  8 = 1/.Under the consideration of turbulence effects,  should be replaced by  total which contains the influence of eddy viscosity: where V is the original molecule viscosity similar to that in laminar flows and V eddy is the eddy viscosity.In the Smagorinsky model, V eddy can be written as [28] V eddy = (  Δ) 2      Θ      , where   is the Smagorinsky constant and Δ is the filter size.|Θ| = √2Θ  Θ  , where Θ  is the strain-rate tensor.By using the resolved nonequilibrium momentum tensor, Θ  can be written as where Combining ( 28) and ( 29), V eddy can be solved by [28] where Ξ = Π  Π  .Then  total can be written as Therefore, in the simulating process of turbulent flows, the eighth and ninth parameters in the relaxation parameter matrix will be  7 =  8 = 2/( + √  2 + (2 √ 2  Δ) 2 √ 2Ξ/).

Domain Boundary Conditions.
In the following study, flow around stationary circular cylinder, flow around stationary square cylinder, and flow around rotating square cylinder are performed based on the above-mentioned large eddy MRT-IB-LB model.The schematic views of flow around circular and square cylinders are shown in Figure 1.The cylinder is placed at the horizontal axis in the channel, and the distance between the cylinder center and the inlet is 10 ( is the characteristic length of the cylinder).The computational domain covers 40×10.For the inlet, the parabolic velocity boundary condition is applied, and the pressure boundary condition is used for the outlet.To obtain the second order accuracy at the top and bottom boundaries, the half-way bounce-back method for the no-slip wall boundary condition is introduced.Then the interest of the study is focused on the proposed method dealing with the fluid-structure interface of the cylinders.In the present study, Δ, Δ are both the lattice units; that is, Δ = 1 and Δ = 1.

Effects of Capillary and Reynolds Number.
As a benchmark case, the flow around stationary circular cylinder is firstly simulated.Keeping  = 20, the number of Lagrangian points arranged on the cylinder surface is, respectively, set as 20, 30, 40, and 60.The pressure distribution on the surface The pressure and drag coefficients are defined as where the pressure on the solid boundary is calculated by It can be clearly seen that when the number is about half of the circumference associated with the arc length,  = 2Δ; the correctness can be guaranteed for both the pressure distribution and the hydrodynamic force.Comparison of streamlines around a stationary circular cylinder with different IB-LB models is depicted in Figure 3. Figure 3(a) shows the results of the traditional IB-LB model [12,29].And it can be clearly seen that the streamlines penetrate the circular cylinder surface, which means the no-slip boundary condition is not satisfied and the mass exchange happens.On the contrary, the penetration is well prevented in present model, as seen in Figure 3(b).
In the previous simulations, the eddy viscosity is not included.To test whether the penetration and the no-slip boundary condition can be satisfied or not in simulating fluid flows at high Reynolds numbers, the cases at Re = 150 and 3900 are both simulated.The eddy viscosity is introduced by modifying the Smagorinsky constant.The Smagorinsky constants are valued as   = 0.001 at Re = 150 and   = 0.16 at Re = 3900.Of course, the flow over a circular cylinder at Re = 150 is still laminar; the Smagorinsky constant valued here is used to the numerical model test.The focus is on whether the penetration and the no-slip boundary condition can be satisfied or not in simulating fluid flows with values of different Smagorinsky constant.  is valued as 0.01 when Re = 150, and the caused eddy viscosity is quite small compared with that induced when   = 0.16.It has been found that too small value of the Smagorinsky constant may cause numerical instability for turbulent flow simulation.However, in the present model, numerical stability can be ensured when the Smagorinsky constant is small enough to get a big range of   that can guarantee the penetration and the no-slip boundary condition.What is more, at high Reynolds numbers, mesh refining can improve the numerical stability.Therefore the cylinder diameter is taken as  = 40 and the rest simulation conditions remain unchanged.
Figure 4 shows the effects of Reynolds numbers on the flow around circular cylinder.The properties including the velocity color map (Figures 4(a) and 4(b)), vorticity contours (Figures 4(c) and 4(d)), and streamlines (Figures 4(e) and 4(f)) are drafted in order to describe the difference of flow states in detail.From the velocity color maps in Figure 4, it can be found that the regular and periodical vortex shedding phenomenon appears when the Reynolds number is 150, which is in general agreement with the results in [30].Obviously, when the Reynolds number rises to 3900, the vortex shedding becomes irregular and deviated from the horizontal axis.The phenomenon is also reported by Breuer [31]: the vortices shedding from the cylinder move downstream along an axis inclined with reference to the symmetry line in simulating the two-dimensional case at Re = 3900.Irregularly the axis of the vortex street changes from positive to negative angles and the other way around [31].With the increasing Reynolds number, the length of the vortex behind the cylinder is shortened, and vortex shedding occurs closer to the rear of the cylinder, associated with the higher Strouhal number recorded in Table 1.It should be stressed here that the streamline in Figure 4 is a good representation that the penetration and no-slip boundary condition can be well satisfied in present simulations.
Table 1 lists the drag coefficients and Strouhal numbers at Re = 150, 3900.It shows that the present numerical results agree well with the results in literature, especially for the case at Re = 150.For the case at higher Reynolds number (Re = 3900), the deviation of drag coefficient is relatively obvious; however, the Strouhal number matches well with that in literature.Two main reasons can cause this deviation: the computational domain whose effects especially focus on the hydrodynamic forces and the three-dimensional turbulence vortex effects.
According to the data recorded in Table 1, it can be found that the periodic time equals 3724 for the case at Re = 150.A nondimensional quantity,   =  0 ⋅ /, is defined [35],      where  0 is the free-stream velocity,  is the cylinder radius, and  is the time step.Then the periodic time can be rewritten as  150 = 11.172 at Re = 150 and  3900 = 9.602 at Re = 3900.Taking one-half cycle, for example, the pressure distributions around the cylinder surface from  = 0 to  = 3/8 ( is the periodic time) are depicted in Figure 5.The pressure distribution owns the same variation tread with that in literature.On both sides of the symmetry line, the pressure variations associated with time increasing have a reverse rule at Re = 150.However, the variation law seems to be more complex for Re = 3900.Compared with the pressure at the front and rear of the cylinder, the pressure on both sides has stronger fluctuations.For the case at Re = 3900, due to the complex vortex shedding forms as illustrated in Figure 4, in addition to both sides, the pressure at the rear also has obvious fluctuations.

Flow around Square
Cylinder.The simulation of flow around square cylinder is also a benchmark case which is often applied to test the newly developed numerical models.With the characteristic length as  = 40 and the computational domain as 40 × 10, the simulations at Re = 150 and 1000 are completed, in which the Smagorinsky constants are valued as   = 0.001 and   = 0.16, respectively.Table 2 lists the drag coefficients and Strouhal numbers at Re = 150 and 1000 compared with that in literature.Similar  with the situation of circular cylinder, the Strouhal numbers obtained at Re = 1000 match well the expected results, but the higher the Reynolds number is, the larger the deviation existing in the calculated hydrodynamic forces is.In addition, the present results have good agreements with the referenced results when the Reynolds number equals 150.The detailed predictions of drag and lift coefficients, as well as Strouhal numbers calculated through the Fast Fourier Transform, are displayed in Figure 6.
Furthermore, in order to investigate the capability of the developed model dealing with moving fluid-structure interfaces, the flow around a rotating square cylinder is   2, we can see the vortex shedding period of the flow around stationary square cylinder at Re = 150.Therefore, five cases with the rotating velocities,  = 8/ 150 , 4/ 150 , 2/ 150 , / 150 , /(2 150 ), are carried out.Additionally, to decrease the influence of initial conditions, 5000 steps are consummated in advance before the cylinder starts to rotate in all simulations.The results show that in the case with  = 8/ 150 , the running progress is stopped soon, due to the occurrence of instability induced by the highest velocity ( ≈ 1.414) at the cylinder corner.Due to the inherent restriction of LBM, the time step and the mesh spacing is coupled.To complete this simulation, only refined mesh is not enough.The values of the variable quantity should firstly be transformed from lattice units to physical units.Then finish the transformation process from physical units to lattice units under the limits in lattice Boltzmann models with considering the numerical instability.
In the case with  = / 150 , contours of the vorticity and pressure distributions around the rotating square cylinder at various stages through one-half of the rotation cycle are illustrated in Figure 7. Four sensor locations of which the attach angles separately equal  = 45, 84.6, 124.1, 163.7 are chosen to provide a subset that can capture some of the more dominant fluid dynamics associated with rotation.Due to the rotation, the traditional Karman vortex street behind the cylinder cannot be seen, but the vortex shedding pattern becomes more irregular.Several small vortices formed nearby the cylinder shed off together and turn into a big vortex.Compared with the vortex shedding pattern of a stationary square cylinder with upward angle of inclination, it can be found that the vortex moves downstream along an axis inclined with reference to the symmetry line, and the inclined angle is almost a constant associated with the referred rotating velocity.The corresponding pressure distributions around the cylinder are shown in Figure 7(b).The results show that the low pressures are just located at the centers of the formed vortices.The low pressures firstly occur nearby the corners of the cylinder where the boundary layer separation happens.Under the effects of the cylinder rotation and the vortex interaction, the low pressure area migrates on the rear surface of the cylinder.A detailed distribution law of the pressure on the cylinder over time is recorded in Figure 8.
Figure 8 displays the pressure distribution on the cylinder surface at four sensor locations just as that identified in Figure 7.The subgraph is a detailed description of the cylinder locations, in which A, B, C, and D represent the four corners and (a), (b), (c), and (d) stand for the four sensor moments.Wherever the four corners are, a large pressure fluctuation will happen near their local positions.The firstlargest pressure fluctuation occurs at the corner C at the sensor moment (c), and the second-largest pressure fluctuation occurs at the corner B at the sensor moment (a).From the subgraph, we can see that the two positions are both located Figure 9 shows the time histories of drag and lift coefficients of the rotating square cylinder with the velocity ranging from  = 0 to  = 2/ 150 .As shown in Figure 9, there is clear oscillation of force coefficients.Due to the inherent restriction, simulation of flows over moving objects by using IBM generally suffers the nonphysical temporal oscillation of pressure fields, which is the origin of spurious force oscillations produced.To overcome this deficiency, some remedies have been proposed [39][40][41].The results show that the   fluctuation amplitudes of the drag and lift coefficients are both enlarged, associated with a nonzero rotating velocity.When  = /(2 150 ), the drag coefficient exhibits the relative larger mean value and fluctuation period.As the velocity increases, the fluctuation period of the drag coefficient is seriously shortened and two main fluctuation frequencies exist.For the lift coefficient, as the rotating velocity increases, the mean lift coefficient value increases.The ways to suppress the force oscillation in the present model will be a direction for further study.

Conclusion
With the introduced Smagorinsky's algebraic eddy viscosity, a novel MRT-IB-LB model is developed for the study of fluid-structure interaction, in which the no-slip boundary condition can be guaranteed by an implicit velocity correction.What is more, the effects of turbulence are considered by resolving the nonequilibrium momentum tensor.As benchmark cases, the flows around cylinders at modified Reynolds numbers are simulated.The conclusions of the results are the following: (ii) The strength of the turbulence effect can be directly controlled through a Smagorinsky constant (  ).When Re = 150, the Smagorinsky constant is valued as 0.001, and when Re = 3900 for the circular cylinder and Re = 1000 for the square cylinder, the Smagorinsky constant equals 0.16.
(iii) The flow patterns and hydrodynamic properties of stationary cylinders are illustrated which are in good agreements with the results in literature.In addition, the results of flow around the rotating square cylinder verify the capability of the MRT-IB-LB model dealing with the moving boundaries.
(iv) At high Reynolds numbers, due to the three-dimensional turbulence effects, obvious relative error between the calculated hydrodynamic forces and that in literature can be found.However, taking advantage of the particle kinetic property of the LB method, the two-dimensional model can be easily transformed into a three-dimensional model.

Figure 1 :
Figure 1: Schematic view of flow around circular and square cylinders.(a) Stationary circular cylinder; (b) rotating square cylinder.

Figure 2 :
Figure 2: Effects of the number of Lagrangian points on the force profiles of the circular cylinder at Re = 20.(a) Pressure distribution on the surface; (b) time-history of drag coefficient.

Figure 3 :
Figure 3: Comparison of streamlines around a stationary circular cylinder at Re = 20.(a) Traditional IB-LB model; (b) present model.

Figure 6 :
Figure 6: Time histories of drag and lift coefficients and power spectra of Strouhal number.

Figure 7 :
Figure 7: Instantaneous vorticity contours and pressure color maps in the wake of the rotating square cylinder.(a) Vortex contours; (b) pressure color map.

Figure 8 :
Figure 8: Pressure distribution on the rotating square cylinder at four sensor locations.The subgraph in (b) is a detailed description of the cylinder locations.

Figure 9 :
Figure 9: Time histories of drag and lift coefficients of the cylinder with different rotating velocities.
(i) The proposed MRT-IB-LB model can well prevent the streamlines penetrating the cylinder surface and keep the no-slip boundary condition for steady flows and unsteady flows at high Reynolds numbers.

Table 1 :
Drag coefficient and Strouhal number of the stationary circular cylinder at Re = 150, 3900.

Table 2 :
Drag coefficient and Strouhal number of stationary square cylinder at Re = 150, 1000.