Mathematical Problems in Engineering

Volume 2018, Article ID 4072758, 12 pages

https://doi.org/10.1155/2018/4072758

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

College of Ship-Building Engineering, Harbin Engineering University, Harbin 150001, China

Correspondence should be addressed to Zhikai Wang; moc.liamg@aw.iakihz

Received 11 January 2018; Revised 8 March 2018; Accepted 15 March 2018; Published 29 April 2018

Academic Editor: Manuel Pastor

Copyright © 2018 Zhikai Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. 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–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–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–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–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 single-relaxation-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 () to be higher according to . 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, , 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 , and the number of compute nodes can be reduced by letting [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.

#### 2. 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 is the particle distribution function at the position , 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 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 where is the discretized direction. The completely discretized equation, that is, LBGK equation, can be written aswhere stands for the particle distribution associated with the 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, , and is the particle velocity, . With the content of conservation laws,and the condition that the Knudsen number is , 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 diagonal relaxation matrix:Correspondingly, the collision term on the right hand of (3) is transformed into the moment space [25]:where is the velocity moment of the distribution function and is the transformation matrix.The moments can be written aswhere 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 asIn the present study, the relaxation parameters in are valued as , , , and . , are related to and the eddy viscosity.

##### 2.2. 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 are the coordinate values in the Lagrangian coordinate system and the Eulerian coordinate system, respectively. 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 . and are the force densities on the structure nodes and the fluid nodes. and 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 are completely dependent on calculating . 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 is the body force. Correspondingly, the evolution equation containing the discretized body force can be modified asCombining (3) and (7), the above evolution equation can be written in the MRT scheme asAnd the velocity can be computed byThe velocity in (18) can be divided into two parts, , where is the intermediate velocity and is the velocity correction. In the conventional IB-LB schemes, is explicitly calculated, and then the velocity correction can be decided in advance. Here, based on the implicit velocity correction, is kept unknown at first and is determined by the velocity correction at the structure nodes satisfying the no-slip boundary condition. Assuming the unknown velocity correction at the structure point is in the Lagrangian coordinate system, the velocity correction at the fluid point in the Eulerian coordinate system can be calculated byWith the continuous kernel distribution [12], (19) can be rewritten aswhere is the arc length of the structure boundary element.

Therefore, the corrected velocity can be expressed asBecause the intermediate velocity is obtained in advance, the only unknown variable in (22) is , 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 is still the only unknown variable, but it can be calculated directly now associated with the specific boundary velocity . After 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 . is the force term whose formation in the moment space can be rewritten as withwhere and are the components of macro velocity and macro force, so (17) can be rewritten as

##### 2.3. 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 diagonal relaxation matrix; that is, . Under the consideration of turbulence effects, should be replaced by which contains the influence of eddy viscosity:where is the original molecule viscosity similar to that in laminar flows and is the eddy viscosity. In the Smagorinsky model, can be written as [28]where is the Smagorinsky constant and is the filter size. , where is the strain-rate tensor. By using the resolved nonequilibrium momentum tensor, can be written aswhere . Combining (28) and (29), can be solved by [28]where . Then can be written asTherefore, in the simulating process of turbulent flows, the eighth and ninth parameters in the relaxation parameter matrix will be .

#### 3. Results and Discussion

##### 3.1. 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 . 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, and .