Abstract
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.
1. 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 body-fitted 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–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–16] and adaptive mesh [17–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 finite-difference LBM is that, with the help of using a special finite-difference 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 angle of attack (AOA) and 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.
2. 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.
2.1. 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 is the velocity of the fluid, is the dynamic viscosity, is the density, is the pressure, and is the body force. The no-slip boundary conditions on the boundary of the obstacle domain in the fluids can be described aswhere is the boundary of the obstacles and is the velocity of the obstacles. For problems involving only fixed obstacles are considered in the present method, the velocity of the obstacles 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 aswhereis the mask function used to describe the obstacles’ geometry and 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 [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 single-relaxation time (SRT) without a forcing term [23] can be written asin which is the particle distribution function; is the time; is the particle velocity in the direction; 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 asand the equilibrium distribution function can be expressed asin which is the sonic speed and the weight factors are , , and . The macroscopic density, , and velocity, , 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.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 . And the second part is contributed by the force density represented by . The can be written asand the asSo (17) can be written as
Substituting (20) into (14), the force density can be expressed asThen substituting (19) into (21), the second part of fluid velocity can be obtained asand the force density can be expressed as
2.3. 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 aswhere is the Courant-Friedrichs-Lewy (CFL) number; is the time step and is the grid spacing in the direction of . It is can be seen in (24) that when the CFL number is equal to 1 () the streaming process becomes equivalent to the perfect shift , that is, (13). When the CFL number is less than 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 asSimilarly, all other in-line hanging nodes can be obtained. Then the out-line hanging node 4 can be calculated asAfter 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 ).(3)Use (25), (26), and (27) to obtain the density distribution function on the hanging nodes.(4)Use (24) to obtain the density distribution function on the nodes of computational grid.(5)Use (10) and (18) to obtain the macroscopic density and intermediate velocity.(6)Use (22) and (23) to obtain the velocity corrections and the force density.(7)Correct the fluid velocity using (20) and compute the equilibrium distribution function using (9).(8)Repeat step to step until the convergence is reached.
3. 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 AOA and 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.
3.1. Flow Past a Fixed Cylinder
In the experiments of this article, the density of the fluid is set as . The free stream velocity is . The diameter of the cylinder is . The numerical simulations occupy a rectangle domain with the height of and the length of . The center of the fixed cylinder is located at , as shown in Figure 3. The computational meshes for these experiments are shown in Figure 4. For the cases of and , 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 ), 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.
(a) Computational domain set up for local refinement
(b) Three-level refined meshes for flows past a fixed circular
(a) Streamlines and velocity contours for
(b) Streamlines and velocity contours for
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 to . These mean that the computational meshes set up for the cases of are good enough to get correct results.
For the cases of and , 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 drag coefficient, the lift coefficient, and Strouhal number are compared with those of previous studies in Table 2. From Table 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 and are given in Figure 6, and the streamlines and velocity contours are shown in Figure 7.
(a)
(b)
(a) Streamlines and velocity contours for
(b) Streamlines and velocity contours for
3.2. 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–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 is set as . For different experiments conducted, .
(a) Computational domain set up for local refinement
(b) Three-level refined meshes for flows past two tandem arranged cylinders
From Figure 10, it can be seen that when is set as , 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 , 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 of upstream and downstream cylinders oscillate at the same frequency, which means the Strouhal numbers are equal. The comparison of the present results with the data in previous literatures is represented in Table 3, from which a good agreement can be seen.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
3.3. Flow Past a NACA-0012 Airfoil
For the practical application of the LBM-VP with local mesh refinement, the flows past a NACA-0012 airfoil at with and at with are selected for numerical experiments. Just like in the previous numerical experiments, the fluid density is set as and the free stream velocity is . The chord of the airfoil is .
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 with zero angle of attack will reach a steady state. And the drag coefficient of the present experiment is , which agrees well with the results given by Lockard et al. [32] with and Wu and Shu [5] with . For the case of with , 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 , while the value given by Suzuki et al. [33] is and that given by Kurtulus [34] is . 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.
4. 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 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 with and at with 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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors would like to acknowledge the support of the Fundamental Research Funds for the Central Universities (HEUCFP201744) and the National Nature Science Foundation of China (51779056).