A Coupled Lattice Boltzmann-Volume Penalization for Flows Past Fixed Solid Obstacles with Local Mesh Refinement

A coupled Lattice Boltzmann-Volume Penalization (LBM-VP) with local mesh refinement is presented to simulate flows past obstacles in this article. Based on the finite-difference LBM, the local mesh refinement is incorporated into the LBM to improve computing efficiency. The volume penalization method is introduced into the LBM by an external forcing term. In the LBM-VP method, the processes of interpolating velocities on the boundaries points and distributing the force density to the Eulerian points near the boundaries are unnecessary. Performing the LBM-VP on a certain point, only the variables of this point are needed, which means the whole procedure can be conducted parallelly. As a consequence, the whole computing efficiency can be improved. To verify the presented method, flows past a single circular cylinder, a pair of cylinders in tandem arrangement, and a NACA-0012 are investigated. A good agreement between the present results and the data in the previous literatures is achieved, which demonstrates the accuracy and effectiveness of the present method to solve the flows past obstacle problems.


Introduction
Flows past obstacles are related to many applications, especially in marine and offshore engineering, aerospace engineering, current or wind turbines, and so on.Numerical studies and simulations of flows past obstacles are of great use in these areas.As an alternative to the traditional Navier-Stokes (N-S) equation solver, the Lattice Boltzmann method (LBM) has been used widely in the simulations of flows past obstacles [1].The reason why LBM is used so popularly is that it has many noticeable advantages like simplicity in coding, parallel computation, and explicit calculation.Just like in the traditional N-S equation solvers, two main methods are used in LBM when dealing with complex boundary: body-fitted grid methods and immersed boundary methods (IBM).
For the body-fitted grid methods, generating a bodyfitted grid is an extremely expensive process when complex boundaries are involved.Still, it is also difficult to create a high quality body-fitted grid with simple boundaries.Besides, in the process of creating body-fitted grids, the structured and unstructured grids are frequently used.However on the structured and unstructured grids, the order of accuracy is lower than that on the Cartesian grids [2].Compared with the body-fitted grid methods, the immersed boundary methods, which are proposed by Peskin [3], can be implemented easily.In IBM, there is no need to create body-fitted grids.The flow field is represented on a fixed Cartesian mesh on which the N-S equations are solved.And the boundary is represented by a Lagrangian grid.The variables on these two grids are related by a discrete delta function interpolation.The brightest spot of IBM is that it uses a restoring force to reflect the boundary effect on the flow.So the N-S equations can be solved without considering the boundaries, which means that the N-S equations can be solved on the whole Cartesian mesh.Feng and Michaelides [4] incorporated the IBM into the LBM firstly to solve fluid particles interaction problems.As mentioned above, the boundary effect on the flows is reflected by restoring force.However, the nonslip boundary conditions are not always guaranteed.As a result, in the simulations of flows past boundaries some streamlines may penetrate the solid body boundary.Wu and Shu [5] proposed a new version of IBM, in which the velocity correction near the boundary is considered as unknown and is computed implicitly to guarantee the nonslip boundary conditions.In the IBM, a Vandermonde matrix should be recomputed at every time step, which means more computing time is needed [2].
Recently, Benamour et al. [6] incorporated another type of IBM, Volume Penalization (VP), which is firstly proposed by Arquis and Caltagirone [7], into LBM.In the VP, the solid boundary is considered as a porous medium with very small permeability.By using a mask function, the solid boundaries are modeled on the Euler grids.So the Lagrangian grids of the VP are just part of the Euler grids which are marked by a mask function.This method has been successfully used by several researchers to simulate and study flows past obstacles [8][9][10].By the forcing term proposed by Guo et al. [11], the force induced by boundaries in the VP method is incorporated into LBM.Compared with the LBM-IBM, there is no need to use a delta function to interpolate the velocity at boundaries in the LBM-VP.Also, the process of distributing the force density into the Eulerian points near the boundaries is unnecessary.Performing the VP procedure on a certain Lagrangian point just needs the variables of this certain Euler point, which means that the VP can be performed parallelly.Considering the parallelizability of the LBM, the whole LBM-VP calculating procedure can be conducted parallelly.
For the standard LBM, the whole calculation process consists of two subprocesses: the collision and the streaming.In the streaming subprocess, particles move from one mesh point to its neighbor points within a time step.So the calculation mesh is limited to uniform mesh, which limits the application of the LBM greatly.To eliminate this disadvantage, many researchers have proposed several improvements to implement the LBM on the nonuniform mesh [12][13][14][15][16] and adaptive mesh [17][18][19].In this article, the finite-difference LBM, proposed by Lee and Lin [20], is adopted to perform the LBM-VP on the nonuniform grids aiming to improve the computing efficiency.One of the highlights of this finitedifference LBM is that, with the help of using a special finitedifference Lattice Boltzmann scheme, all of the blocks at different refinement levels can be advanced in time with the same time step.So it is very easy to conduct the LBM on the nonuniform mesh, especially when the VP or other IBM are incorporated into LBM.
In this article, to validate the proposed LBM-VP with local mesh refinement and its computing efficiency, flows past a circular cylinder, flows past a pair of circular cylinders in tandem arrangement, and flows past the NACA-0012 airfoil with 10 ∘ angle of attack (AOA) and 15 ∘ AOA are studied.And the results are compared with available data in previous literatures.The computing time of the LBM-VP performed on nonuniform mesh is also compared with that of the LBM-VP performed on uniform mesh to show the computing efficiency.The rest of this paper is arranged as follows.In Section 2, the volume penalization method, the Lattice Boltzmann method, and the local mesh refinement are described.Also the whole computing procedure is given in this section.The numerical experiments and the comparison of results are given in Section 3. In Section 4, some concluding remarks and recommendations for the future work are presented.

Mathematical and Numerical Formulation
In this section, the fluid-solid interaction (FSI) technique based on the volume penalization method is briefly introduced firstly.Then the incorporation of volume penalization method into the Lattice Boltzmann method and the LMR technique used in the Lattice Boltzmann are introduced in detail.Finally, a whole computing procedure is given.

Brief Review of the Volume Penalization Method and LBM.
Let us consider the fluid-structure interaction (FSI) between incompressible viscous fluid and rigid obstacles in the fluid.The motion of the fluid is governed by the incompressible Navier-Stokes equations: where u is the velocity of the fluid,  is the dynamic viscosity,  is the density,  is the pressure, and F  is the body force.The no-slip boundary conditions on the boundary of the obstacle domain Ω  in the fluids can be described as where Ω  is the boundary of the obstacles and U  is the velocity of the obstacles.For problems involving only fixed obstacles are considered in the present method, the velocity of the obstacles U  is equal to zero.The computational domain is shown in Figure 1.Ω  is the fluid domain.The union of these two domains Ω = Ω  ∪ Ω  is the entire domain.
The Dirichlet problem (1)-( 3) can be solved by the volume penalization method [8,21].In the volume penalization method, the solid obstacles are modeled as porous media.By adding a penalization term on the velocity, the momentum equation ( 1) is modified as where is the mask function used to describe the obstacles' geometry and  ≪ 1 is the penalization parameter.It can be seen that there is no Dirichlet boundary conditions in (4).The solution of the penalized N-S equation ( 4) tends towards the exact solution of N-S equation imposing no-slip boundary conditions with  → 0 [8,21,22].The hydrodynamic forces acting on the fixed obstacle can be obtained through The N-S equation can be recovered by the Lattice Boltzmann equation through the multiscale Chapman-Enskog analysis.The discrete Boltzmann equation (DBE) with singlerelaxation time (SRT) without a forcing term [23] can be written as in which   is the particle distribution function;  is the time; e  is the particle velocity in the th direction;  eq  is the equilibrium distribution function, and  is the relaxation parameter.For the nine-velocity lattice in two dimensions (D2Q9), the discrete velocity vectors can be defined as and the equilibrium distribution function can be expressed as in which   = 1/ √ 3 is the sonic speed and the weight factors are  0 = 4/9,  1-4 = 1/9, and  5-8 = 1/36.The macroscopic density, , and velocity, u, are defined as follows: Equation ( 7) can be split into two substeps [20]: collision: and streaming: where  is the single-relaxation parameter and  is related to the kinematic viscosity  through the equation  = /( 2  Δ).

The Incorporation of VP into LBM.
The added penalization term on the velocity in the modified N-S equation can been regarded as an external force density, acting on the fluid phase.Similarly, the DBE should be also modified by adding an external forcing term.In our article, the form of DBE with an external force term proposed by Guo et al. [11] is adopted.And the external force term is introduced into the collision substep.The collision substep is modified as From ( 17), the fluid velocity is made up by two parts [5].The first part is contributed by the density distribution function represented by the intermediate velocity u * .And the second part is contributed by the force density F VP represented by u.The u * can be written as and the u as So ( 17) can be written as Substituting ( 20) into ( 14), the force density can be expressed as Then substituting ( 19) into ( 21), the second part of fluid velocity can be obtained as and the force density can be expressed as

The Local Mesh Refinement Technique in LBM.
As proposed in [20], the streaming substep ( 13) is solved by using a Lax-Wendroff scheme with second-order accuracy in both time and space.The streaming substep can be rewritten as where  = |e  |Δ/|Δx  | is the Courant-Friedrichs-Lewy (CFL) number; Δ is the time step and Δx  is the grid spacing in the direction of e  .It is can be seen in ( 24) that when the CFL number is equal to 1 ( = 1) the streaming process becomes equivalent to the perfect shift   (x,  + Δ) =   (x − Δx  , ), that is, (13).When the CFL number is less than 1 ( < 1), the streaming substep is solved with second-order accuracy in both time and space.It should also be pointed out that when the time step Δ is defined, the CFL number changes only due to the grid spacing.Without loss of generality, a local fine-grid domain which is surrounded by a coarse-grid is considered, as shown in Figure 2.For the collision substep's localization feature, there is no need to transfer information between two different refinement levels grids.However, in the streaming substep, the transfer of information is required.To conduct this information transfer, some hanging nodes are needed.These hanging nodes can be classified into three types: in-line hanging node (ILHN), out-line hanging node (OLHN), and corner hanging node (CHN).The in-line hanging nodes will be calculated firstly.By using an interpolation scheme, the density distribution function values of in-line hanging node 3 can be calculated as Similarly, all other in-line hanging nodes can be obtained.
Then the out-line hanging node 4 can be calculated as After all other out-line hanging nodes are calculated, similarly, the out-line hanging node 2 is gotten.The corner hanging node 2 can be calculated: When density distribution function values at all hanging nodes are obtained, the streaming substep can be conducted.In summary, the general steps of the algorithm are as follows.
(1) Design the computational grid, and arrange initial values on the computational grid.
(2) Use ( 15) to obtain the density distribution function after the collision substep of the  =   step on the points of different refinement levels (initially setting   = 0).
(4) Use (24) to obtain the density distribution function on the nodes of computational grid.
(6) Use ( 22) and ( 23) to obtain the velocity corrections and the force density.

Numerical Results and Discussions
To verify the present Lattice Boltzmann-Volume Penalization with local mesh refinement, test cases involving incompressible viscous flows past single fixed circular cylinder are chosen as numerical experiments.Then flows past a pair of circular cylinders in tandem arrangement and flows past NACA-0012 airfoil with 10 ∘ AOA and 15 ∘ AOA are conducted to validate the present method for complex boundaries.Some existing results are chosen as references for comparison.The computational efficiency is given as well.
In the following numerical experiments, the Reynolds number (Re) is defined as where  ∞ is the free stream velocity and  is the diameter of the cylinder and  is the kinematic viscosity of the fluid.The drag coefficient of cylinder is defined as where   is the drag force.The lift coefficient is defined as where   is the lift force.For the unsteady flow, the Strouhal number is defined as where   is the vortex shedding frequency.

Flow Past a Fixed Cylinder.
In the experiments of this article, the density  of the fluid is set as 1.0.The free stream velocity  ∞ is 0.1.The diameter  of the cylinder is 40.
The numerical simulations occupy a rectangle domain with the height  of 1024 and the length  of 2048.The center of the fixed cylinder is located at (/2, /2), as shown in Figure 3.The computational meshes for these experiments are shown in Figure 4.For the cases of Re = 20 and 40, the flows will reach a steady state.Behind the fixed cylinder, a pair of stationary recirculating eddies will develop.
With the Reynolds number increasing, the end of the wake will also move farther away from the rearmost point of the cylinder.The comparison of drag coefficient   , length of the recirculation region  (scaled by /2), and separation angle   with the results of previous literatures is shown in Table 1.
The contour plots of velocities and streamlines are shown in Figure 5.
From Table 1, the present numerical results agree well with those in the previous studies.It can be seen in Figure 5 that the region of the recirculation becomes bigger when the Reynolds number grows from 20 to 40.These mean that the computational meshes set up for the cases of Re = 20, 40 are good enough to get correct results.
For the cases of Re = 100 and 200, the flows will reach a unsteady state.A Kármán vortex street will develop behind the cylinder.Also, a lift force will act on the cylinder.With the augmentation of Reynolds number, the vortices value and the area of vortices behind the cylinder will increase.The   2, it can be found that a good agreement is gotten between the presents results and those in previous studies.The evolution of drag and lift coefficients for the cases of Re = 100 and 200 are given in Figure 6, and the streamlines and velocity contours are shown in Figure 7.

Flow Past Two Tandem Arranged Cylinders.
To verify the capability of the present method to simulate complex flows, flows past a pair of circular cylinders tandem arranged, which have been studied by many researchers [28][29][30][31], are chosen as experiments.Compared with the flows past a single circular cylinder, the flows past a pair of cylinders tandem arranged are more complicated.For this problem, the distance between these two cylinders represented by /, as shown in Figure 8, plays an extremely important role in the development of vertex structures of the flows, as well as the evolution of drag and lift coefficients.In Figure 8, the computational domain and boundary conditions are presented, and the computational meshes are shown in Figure 9.For the presented experiments, the Reynolds number Re is set as 200.For different / experiments conducted, / = 1.5, 2.0, 3.0, 4.0.
From Figure 10, it can be seen that when / is set as 1.5, 2.0, 3.0, vortex shedding develops behind the downstream cylinder.But between the upstream cylinder      and the downstream cylinder, there is no vortex shedding.For these three experiments, the drag coefficients of the upstream cylinders are positive, as shown in Figure 11, but the downstream cylinders' are negative.The drag coefficient amplitudes of both upstream and downstream cylinders reduced further with the augment of the distance between two cylinders.When / turns to 4.0, great changes happen.The vortex shedding develops not only behind the downstream cylinder but also between the upstream cylinder and the downstream cylinder.The drag coefficient becomes positive.For all the four experiments, the drag coefficients    3, from which a good agreement can be seen.The streamlines and the velocity contours for flows past the NACA-0012 airfoil are shown in Figures 12 and 13.The flow behind the NACA-0012 airfoil at Re = 500 with zero angle of attack will reach a steady state.And the drag coefficient of the present experiment is 0.179, which agrees well with the results given by Lockard et al. [32] with 0.1762 and Wu and Shu [5] with 0.1759.For the case of Re = 1000 with AOA = 15 ∘ , the flow becomes unsteady, and a lift force will act on the airfoil.In the present experiment, the Strouhal number for the drag coefficient is 0.72, while the value given by Suzuki et al. [33] is 0.73 and that given by Kurtulus [34] is 0.712.Obviously, a good agreement is gotten between our result and those in the previous literatures.From the two experiments, it can be concluded that the present LBM-VP with local mesh refinement can solve the practical problems with complex boundaries.
To demonstrate the advantage of the LBM-VP with local mesh refinement, we compare the computing time and the number of the used mesh points of the LBM-VP on the local refined mesh with those on the uniform mesh code.All the experiments are conducted by eight threads on a PC with Intel(R) Core(TM) i7-4790K 4.00 GHz CPU and 16.0 GB RAM.In Table 4, the computing time per step and the number of the used points of the LBM-VP on local refined mesh are compared with the uniform mesh code for all the experiments in this article.For all the experiments, it can be seen that the local mesh refinement technique used in the LBM-VP can reduce the computing time per step.

Conclusions
In this article, the coupled Lattice Boltzmann-Volume penalization method (LBM-VP) with local mesh refinement is presented to simulate flows past obstacles.Compared to the direct-forcing IBM and the velocity IBM, there is no need to use a delta function to interpolate the velocity at boundaries and to distribute the force density into the Eulerian points near the boundaries.On a certain Lagrangian point, only the variables of this certain Euler point are needed to perform the VP procedure, which means that the LBM-VP can be conducted parallelly.By incorporating the local mesh refinement technique based on the finite-difference LBM, the LBM-VP are performed on the local refined mesh.As a result, the number of used mesh points is reduced, thus improving the computing efficiency.To justify the present LBM-VP with local mesh refinement, flows past a single circular cylinder with Re = 20, 40, 100, 200 are carried out firstly.Then slightly more complex experiments, flows past a pair of circular cylinders in tandem arrangement, are studied.For the practical application of the present method, the flows past a NACA-0012 airfoil at Re = 500 with AOA = 0 ∘ and at Re = 1000 with AOA = 15 ∘ are chosen as experiments.
By comparing the numerical experiment results with the data in previous literatures, our simulations show good agreement with available data.And the local mesh refinement technique reduces the number of used points of the LBM-VP.As a result, the computing time of per step is reduced, meaning that the computing efficiency improved.Although the current method is developed in 2D, it is easy to extend it to 3D by replacing D2Q9 lattice with D3Q15 or D3Q19.

Figure 1 :
Figure 1: The computational domain of obstacle and fluid.

Figure 2 :
Figure 2: Grid structure with difference refinement levels.

Figure 3 :
Figure 3: Computational domain and boundary conditions for flows past a fixed cylinder.

( a )Figure 4 :
Figure 4: Computational meshes for flows past a cylinder.

Figure 5 :
Figure 5: Steady flow past a cylinder at Re = 20 and Re = 40.

Figure 6 :
Figure 6: Evolution of drag and lift coefficients at Re = 100 and Re = 200.

( a )
Computational domain set up for local refinement (b) Three-level refined meshes for flows past two tandem arranged cylinders

Figure 9 :
Figure 9: Computational meshes for flows past a pair of circular cylinders in tandem arrangement.

Figure 10 :
Figure 10: Streamlines and velocity contours for flows past two cylinders in tandem arrangement.

Table 1 :
Drag coefficient, length of bubbles, and separation angle for flows past a cylinder at Re = 20 and Re = 40.
∘drag coefficient, the lift coefficient, and Strouhal number are compared with those of previous studies in Table2.From Table

Table 2 :
Drag coefficient, length of bubbles, and separation angle for flows past a cylinder at Re = 100 and Re = 200.

Table 3 :
Drag coefficient and Strouhal number for flows past a pair of circular cylinders in tandem arrangement at Re = 200.

Table 4 :
Comparisons of the computing time per step and number of the used points.