A hybrid wavelet-based adaptive immersed boundary finite-difference lattice Boltzmann method for two-dimensional fluid–structure interaction

Journal of Computational Physics 333 (2017) 24–48

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A hybrid wavelet-based adaptive immersed boundary finite-difference lattice Boltzmann method for two-dimensional fluid–structure interaction Xiongwei Cui, Xiongliang Yao ∗ , Zhikai Wang, Minghao Liu College of Shipbuilding and Engineering, Harbin Engineering University, Harbin, 150001 China

a r t i c l e

i n f o

Article history: Received 22 June 2016 Received in revised form 26 October 2016 Accepted 13 December 2016 Available online 19 December 2016 Keywords: Lattice Boltzmann Finite difference method Immersed boundary method Second generation wavelets Adaptive wavelet collocation method

a b s t r a c t A second generation wavelet-based adaptive finite-difference Lattice Boltzmann method (FD-LBM) is developed in this paper. In this approach, the adaptive wavelet collocation method (AWCM) is firstly, to the best of our knowledge, incorporated into the FD-LBM. According to the grid refinement criterion based on the wavelet amplitudes of density distribution functions, an adaptive sparse grid is generated by the omission and addition of collocation points. On the sparse grid, the finite differences are used to approximate the derivatives. To eliminate the special treatments in using the FD-based derivative approximation near boundaries, the immersed boundary method (IBM) is also introduced into FD-LBM. By using the adaptive technique, the adaptive code requires much less grid points as compared to the uniform-mesh code. As a consequence, the computational efficiency can be improved. To justify the proposed method, a series of test cases, including fixed boundary cases and moving boundary cases, are invested. A good agreement between the present results and the data in previous literatures is obtained, which demonstrates the accuracy and effectiveness of the present AWCM-IB-LBM. © 2016 Elsevier Inc. All rights reserved.

1. Introduction As an alternative to the traditional Navier–Stokes (N-S) equation solver, the lattice Boltzmann method (LBM) has been used popularly in the simulation of fluid flow problem over the last two decades [1]. Compared to the traditional numerical methods, the LBM is very new. Its noticeable advantages, such as simplicity, parallelizability and explicit calculation lead to its great popularity. With years of continuous improvement, the LBM has been successfully used in many different areas, such as flow dynamics, noise control problems and fluid–structure interaction problems [2–5]. Based on the kinetic theory, the LBM use the density distribution functions to describe the fictitious particles. The dynamic of the fictitious particles is studied by the evolution of the density distribution functions. For the standard LBM, the whole calculation process can be made up by two sub-processes, i.e. the collision and the streaming. The macroscopic variables, such as density and velocity, can be calculated by the distribution functions. In the standard LBM, it is required that particles should move from one mesh point to its neighbor points within a time interval in the streaming sub-process. So the calculation mesh is limited to uniform mesh. Because of this mesh requirement, the application of the LBM is limited

*

Corresponding author. E-mail address: [email protected] (X. Yao).

http://dx.doi.org/10.1016/j.jcp.2016.12.019 0021-9991/© 2016 Elsevier Inc. All rights reserved.

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

25

greatly. To eliminate this shortcoming, great improvements have been made by many researchers to implement the LBM on the non-uniform mesh. The strategies can be roughly classified into two types: interpolation type and differential type [6]. For the first type, He and his co-workers [7,8] firstly introduced an interpolation step into the LBM and proposed the interpolation-supplemented lattice Boltzmann method (ISLBM). Based on the ISLBM, some other interpolating type LBM were proposed. By proposing a generalized form of the interpolation supplemented LBM (GILBM), Imamura et al. [9] make the ISLBM applicable to arbitrarily structured mesh. And by using the local time step, the GILBM has a very high efficiency. Based on the Taylor series expansion and least squares theory, Shu et al. [10,11] proposed the Taylor series expansion and least squares-based lattice Boltzmann method (TLLBM). This method has been successfully used to simulate flows with complex geometry. For the second type, by using Taylor series expansion, the lattice Boltzmann equations (LBEs) are transformed into a set of differential equations. These differential equations can be solved by the conventional numerical methods such as finite difference (FD), finite element (FE), and finite volume (FV). Perhaps, the finite difference LBM (FDLBM) was firstly proposed by Reider and Sterling [12], and then was examined by Cao et al. [13]. Benefiting from the work of Mei and Shyy [14], the FDLBM can be implemented on the curvilinear coordinates. By using a predictor-corrector method, Lee and Lin [15] proposed the finite-element LBM (FELBM). And then, some other finite element-based lattice Boltzmann methods were proposed [16,17]. By introducing the FVM to solve the LBMs, Succi et al. [18] firstly proposed the finite volume LBM (FVLBM). After that, based on the cell-centered data structure-based FVM, the FVLBM were implied on the structured mesh [19]. And based on the vortex-centered data structure-based FVM, the FVLBM has been successfully implied on the unstructured mesh [20,21]. There are many adaptive-mesh-refinement (AMR) techniques used in the FD, FE and FV. By introducing these techniques into the LBM, an adaptive lattice Boltzmann can be proposed. J. Wu and C. Shu [6] proposed a solution-adaptive LBM by introducing an efficient stencil adaptive technique, originally proposed by Ding and Shu [22], into the ISLBM. In the method, by coupling 5-point orthogonal symmetric stencil and 5-point diagonal symmetric stencil, a 9-point symmetric structure can be generated, which is similar to the D2Q9 lattice model. Guzik et al. [23] proposed a novel adaptive FVLBM. In the adaptive FVLBM, the key point is transferring information across interfaces between different grid resolutions. This key point was solved by a new method which was developed following established techniques for finite-volume representations. On the basis of a block-structured AMR, A. Fakhari and T. Lee [24] proposed an adaptive finite-difference lattice Boltzmann method. By using pointers, there is no need to search all blocks to find the neighboring blocks of a given block in this method. Another highlight is that the use of a special finite-difference lattice Boltzmann scheme which is following the idea proposed in [25] ensures that all of the blocks at different refinement levels can be advanced in time with the same time step. More recently, based on a “bubble” function which can eliminate the need of time interpolation, Guo et al. [26] proposed a hybrid adaptive LBM to simulate two-dimensional flow. It is obvious that incorporating the AMR techniques, used in the FD, FE and FV, into the LBM for more efficient and effective adaptive method attracts increasing attention. In this paper, the Adaptive Wavelet Collocation Method (AWCM) is firstly, to the best of our knowledge, incorporated into FDLBM to develop a new adaptive LBM. The AWCM is based on the second generation wavelet transformation. By being wavelet transformed, a function’s local regularity can be reflected by wavelet amplitudes. In the rapidly varying regions, wavelet amplitudes are relatively large. While in the regions where the solution is well-behaved, the wavelet amplitudes are small. So in the regions where wavelet amplitudes are small, there is no need to use such a fine grid as in the regions where wavelet amplitudes are large. By omitting points with small wavelet amplitudes, as well as the wavelet transformation corresponding to these points, a sparse grid can be generated. And this sparse grid is dynamically adapted to evolving features of the solution [27]. Besides, the partial differential equations (PDEs) are solved in the physical space in the AWCM which is not the same as in the other transform-based methods such as Fourier transform-based method and adaptive wavelet Galerkin method. As a result, the derivations can be conveniently and efficiently approximated by using finite-differences, which is proposed by O.V. Vasilyev and C. Bowaman [28]. The details of this procedure and the error bound of this finitedifferences based derivation approximation are given in Ref. [28]. Since it was proposed, the AWCM has been a very useful solver for PDEs such as parabolic [28], hyperbolic [29] and elliptic [30]. It has been also applied to various engineering areas by many researchers [27,29,31–33]. In this paper, the AWCM is incorporated into the FDLBM proposed by T. Lee and C.-L. Lin [25]. Based on the grid refinement criterion which is proposed in this paper and explained in the following text, the adaptive omission and addition of computing points are conducted. As a result, a sparse grid is generated. The computing procedures are conducted on this sparse grid. For the computing grid to be compressed, the computing resources can be saved and computing efficiency can be improved. As is well known, some special treatments are needed in using FD-based derivation approximation near the boundaries. Even, special treatments will be more complex if there are one or more boundaries in the inner field of the calculation domain. And the AWCM is not suitable for the complex body-fitted mesh. So if the inner boundaries can be eliminated, the use of FD-based derivation approximation and AWCM will be more convenient and the computing efficiency will be also improved. Fortunately, the immersed boundary method (IBM) proposed by Peskin [34] can eliminate the need of these special treatments. In the IBM, a fixed Cartesian mesh is used to represent the flow field, and a Lagrangian mesh is used to represent the boundary. The discrete delta function interpolation is used to relate the variables on the two meshes to each other. As the brightest spot of the IBM, the effect of boundary is reflected by the restoring force acting on the Eulerian mesh in the vicinity of boundary. The N-S equations with the restoring force are solved over the whole Eulerian domain. So there is no inner boundary in solving the N-S equations, which is very suitable for the implementation of the AWCM and the FD-based derivation approximation. The IBM can be used to solve not only the fixed boundary problems but

26

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

also the problems involving moving boundaries. The IBM was firstly incorporated in LBM by Feng and Michaelides [35,36]. The inner boundaries in LBM can be also eliminated by using the IBM as well. However, as pointed above, the boundary effect is reflected by restoring force in the traditional IBM. There is no enforcement of non-slip boundary conditions. As a result, some streamlines may penetrate the solid body boundary. To overcome this drawback, Wu and Shu [37] proposed a new version of immersed boundary-lattice Boltzmann method (IB-LBM) in which the velocity correction is considered as unknown and is computed implicitly to ensure that the non-slip boundary condition is enforced. The restoring force can be computed from the obtained velocity correction. And more accurate numerical results can be also obtained. As same as the traditional IBM, the IBM proposed by Wu and Shu can be also used to solve the fixed and moving boundaries problems. The criterion for controlling grid refinement plays an extremely key role in the AMR techniques. The construction of the criterion can be based on many variables in LBM, such as macroscopic variables including the density, the velocity, the vorticity, etc. In the previous adaptive LBM, most of the criterions are based on the macroscopic variables either the velocity or the vorticity. However, in fact, the calculation process of the LBM is the evolution of the density distribution functions. The macroscopic variables i.e. the density, the velocity and the vorticity can be obtained by based on the density distribution functions. So why we don’t construct a new grid refinement criterion based on the density distribution functions? In this paper, the first attempt, to our best knowledge, to construct a new grid refinement criterion based on density distribution functions is conducted. Besides, in the previous adaptive LBM the relative differences [6,24] or magnitudes [23] of vorticity have been used as criterion successfully and good results were obtained. The first- or second-order derivatives of the velocity are also a very good choice for the criterion [26]. The criterions mentioned above are belonging to either the error estimate-based or the heuristic criterion. As pointed in the Ref. [38], the effectiveness of this type criterion requires a posteriori error estimates for the actual error introduced by AMR is not well understood. However, the order of accuracy of the algebra and finite-difference derivative approximation of the transformed function on the sparse grid can be obtained analytically, which has been given in Ref. [39] in detail. As a result, it is extremely reasonable to construct the grid refinement criterion based on the wavelet amplitudes of distribution functions. Inspired by the series success of the AWCM in solving partial difference equations and the IB-LBM in simulating the fluid–solid problems, these two strategies are firstly, to the best of our knowledge, coupled to develop a new adaptive lattice Boltzmann method in this paper. The rest of this paper is arranged as follows. In Section 2, the immersed boundary finite-difference lattice Boltzmann method is described. In Section 3, the AWCM algorithm is elaborated and the coupling of algorithm of the AWCM and the finite-difference LBM is also given. A series of test cases of flows past a fixed circular cylinder, a rigid circular cylinder oscillating in a viscous fluid at rest and flows past a rigid oscillating cylinder are studied and numerical results are compared with available data in previous literatures in Section 4. Some concluding remarks and recommendations for the future work are presented in Section 5. 2. Immersed boundary finite-difference lattice Boltzmann method The discrete Boltzmann equation (DBE) with single-relaxation time (SRT) without a forcing term [40] can be written as: eq

∂ fα fa − fα + eα · ∇ f α = − , ∂t λ

(1) eq

in which f α is the particle distribution function; t is the time; eα is the particle velocity in the α th direction; f α is the equilibrium distribution function and λ is the relaxation parameter. For the nine-velocity lattice in two dimension (D2Q9), the discrete velocity vectors can be defined as:

eα =

⎧ ⎨

α=0 (0, 0), ( cos θ , sin θ ), θ = α − 1 ) π /2, α = 1 − 4 , ( α a a ⎩√ 2(cosθα , sinθa ), θa = (2α − 9) π /4, α = 5 − 8

(2)

and the equilibrium distribution function can be expressed as:



eq

fα = ρ wα 1 +

eα · u c 2s

√

+

(ea · u)2 − c 2s |u|2 2c 4s



,

(3)

in which c s = 1/ 3 is the sonic speed and the weight factors are w 0 = 4/9, w 1−4 = 1/9 and w 5−8 = 1/36. The macroscopic density, ρ , and velocity, u are defined as follows:

ρ=



fα

(4)

α

ρu =



eα f a .

(5)

α

The Eq. (1) can be split into two sub-steps [25]: collision eq

gα = f α −

fα − fα

τ + 0.5

(6)

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

27

Fig. 1. Schematic of the IB-LBM (Lagrangian points: circular solid points; Eulerian points: square solid points).

and streaming

∂ gα + eα · ∇ g α = 0 , (7) ∂t where τ is the single relaxation parameter and τ is related to the kinematic viscosity υ through the equation: τ =  υ / c 2s t . In the streaming sub-step Eq. (7), by using the central differences to discretize the first and second derivatives along characteristic lines, the Lax–Wendroff scheme with second-order accuracy in both time and space can be obtained [24,25]:

g α (x, t + t ) =

σ (1+σ ) 2





g α (x − xa , t ) + 1 − σ 2 g α (x, t )

− σ (12−σ ) g α (x + xα , t )

,

(8)

where σ = |eα | t /| xα | is the Courant–Friedrichs–Lewy (CFL) number; t is the time step and xα is the grid spacing in the direction of eα . It is can be seen in Eq. (8) that when the CFL number is equal to 1 (σ = 1) the streaming process becomes equivalent to the perfect shift g α (x, t + t ) = g α (x − xα , t ), as is the case in the conventional LBM. When the CFL number is less than 1 (σ < 1), the streaming sub-step is solved with second-order accuracy in both time and space. It is should also be pointed out that when the time step t is defined, the CFL number changes only due to the grid spacing. In the IBM approach, a set of Eulerian points are used to represent the flow field, and a set of Lagrangian points are used to represent the boundary of a curved object which is immersed in the fluid field. For the case that the LBM is chosen as the N-S equation solver, the Eulerian points are the fixed Cartesian mesh nodes used in the LBM in fact. The two sets of points in the IB-LBM are showed in Fig. 1. Following the idea proposed by Wu and Shu [37], we incorporate the implicit velocity correction-based immersed boundary into the above finite difference lattice Boltzmann method. The external force term is introduced into the collision sub-step [41], so the collision sub-step with external force term can be written as eq

fα − fα

gα = f α − Fα = 1 −

ρu =



+ F α · t ,

 eα − u eα · u · f, wα + · e α 2 4

τ + 0.5



1



2λ

eα f α +

α

cs

1 2

f · t .

cs

(9) (10) (11)

The force density f is distributed from the Lagrangian points and w α are the weight factors used in the equilibrium distribution function. From Eq. (11), the fluid velocity is made up by two parts. One is contributed from the density distribution function, which is represented by the intermediate velocity u∗ . And the other is contributed from the force density f, which is represented by the velocity correction δ u. The u∗ can be expressed as

ρ u∗ =



eα f α

(12)

α

and the δ u as

1

ρ δ u = ft. 2

(13)

So Eq. (11) can be written as

u = u∗ + δ u.

(14)

28

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

The velocity correction δ u at Eulerian points is distributed from the velocity correction at the immersed boundary points. Through the Dirac delta function interpolation, the velocity correction δ u can be obtained from the following equation:



δ u ( x, t ) =

δ u B (X B , t ) δ (x − X B (s, t )) ds.

(15)

Here X B (sl , t ) , l = 1, 2, 3, · · · , S represents the Lagrangian points. δ (x − X B (s, t )) is smoothly approximated by a continuous Kernel distribution D i j , which was proposed by Peskin [34]:

(1 + cos (π |r |/2)) |r | ≤ 2 , |r | > 2 0       D i j xi j − XlB = δ xi j − X lB δ y i j − Y lB .

δ (r ) =

1 4

(16) (17)



From Eq. (17), it can be seen that if an Eulerian point xi j with xi j − X lB > 2 or y i j − Y lB > 2, D i j xi j − XlB



will is

equal to zero. This means that we can just calculate the Eulerian points xi j with xi j − X lB ≤ 2 or y i j − Y Bl ≤ 2. These l Eulerian points are called neighboring Eulerian points of XlB , which can be defined as a set of Eulerian points, P nei :







l S 1 2 P nei = xi j xi j − X lB ≤ 2, y i j − Y lB ≤ 2 . The union of all these sets of boundary points P nei = P nei ∪ P nei ∪ · · · ∪ P nei is called the boundary-neighboring-Eulerian-points (BNEP). So Eq. (15) can be approximated by

    l  l  δ u xi j , t = δ u B X B , t D i j xi j − XlB sl (l = 1, 2, · · · , S ),

(18)

l

where sl is the arc length of the boundary element. So the fluid velocity at Eulerian points xi j can be corrected as:

   δ u xi j , t = u∗ xi j , t + δ u xi j , t .

(19)

To satisfy the non-slip boundary condition, the fluid velocities at the boundary points XlB mush be equal to the wall veloci ties UlB XlB , t . The fluid velocities at the boundary points XlB can be also obtained through interpolation using the smooth delta function. The mathematical expression of this equality relation of velocities is





UlB XlB , t =

 







u xi j , t D i j xi j − XlB x y .

(20)

i, j

Substituting Eq. (19) into Eq. (20), the following equation can be obtained





UlB XlB , t =

 i, j

+









u∗ xi j , t D i j xi j − XlB x y

   i, j

l

δ ulB



      X l , t D i j xi j − X l sl D i j xi j − X l x y . B B B

(21)

In Eq. (21), the intermediate velocity u∗ can be obtained from Eq. (12). The unknowns in Eq. (21) are the velocity corrections, δ ulB . Now the key mission is solving Eq. (21). In Ref. [37], Eq. (21) was further re-written as the following matrix form

AX = B,   where X= δ u1B , δ u2B , · · · , δ u BS ; B = { u1 , u2 , · · · , u S } T with

     

ul = UlB XlB , t − u∗ xi j , t D i j xi j − XlB x y (l = 1, 2, · · · , S ).

(22)

(23)

i, j

The elements of matrix A are only related to the Lagrangian boundary points and the BNEP. After obtaining the velocity corrections at boundary points by solving Eq. (22), the velocity corrections and the corrected velocities can be calculated by Eq. (18) and Eq. (19). The calculation of force on BNEP can be simply calculated, which is an advantage of this version IBM. From Eq. (13), the force density can be calculated by

f = 2ρ δ u/ t .

(24)

Similarly, the force on the boundary points can be from







F XlB = 2ρ δ ulB t .

(25)

The solution procedure of the present immersed boundary finite-difference lattice Boltzmann method can be summarized as follows:

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

29

(1) Arrange initial values, and compute the matrix A and get A−1 . (2) Use Eq. (9) to obtain the density distribution function after the collision sub-step of the t = tn step (initially setting F α = 0). (3) Use Eq. (7) to obtain the density distribution function after the streaming sub-step of the t = tn step. (4) Use Eq. (4) and Eq. (12) to compute the macroscopic variables. (5) Solve equation Eq. (22) to obtain the velocity corrections at all boundary points and use Eq. (18) to get velocity corrections at BNEP. (6) Correct the fluid velocity at BNEP using Eq. (14) and obtain the force density using Eq. (24) (7) Compute the equilibrium distribution function using Eq. (3). (8) Repeat step (2) to step (7) until convergence is reached. 3. Adaptive wavelet collocation method In this section, we briefly explain the second wavelet transform, which is the basis of the AWCM. More details about wavelet transform are given in Refs. [28,39,42]. Then the AWCM and the scheme of incorporating AWCM into the finite difference LBM are given. 3.1. Wavelet transform The multi-dimension wavelet transform formulation can be accomplished by tensor products of the one-dimension wavelet basis [42]. So the one-dimension interpolating wavelet transform is expressed firstly, and the multi-dimension wavelet transform is expressed then. The 1D second generation wavelets are constructed on an . The performance of the construction is on a set of interpo  j

xk ∈  , which forms a set of nested spaces:

lating points,



j



j +1

j

V j = xk ∈  : xk = x2k , k ∈ K j ,

(26)

j

j

j

where xk are the collocation points of the j level of resolution. It can be seen that the restriction xk = x2k guarantees the nestedness of the spaces, i.e. V j ⊂ V j +1 . And let W j be the complement of V j in V j +1 . Following the construction of the second generation wavelets presented in Ref. [42] and Ref. [43], the scaling functions j j φk denoting the bases of the spaces V j and the wavelets ψk denoting the bases of the spaces W j are constructed such that a function u (x) can be decomposed as

u (x) =



j

j

ck0 φk 0 (x) +

j

j

dl ψl (x),

(27)

j = j 0 l∈L j

k∈K j 0 j

+∞  

j

where ck0 and dk are the scaling and wavelet coefficients, respectively. K j and L j are the index sets associated with scaling functions and wavelets on level j, respectively. j 0 is the coarsest level. The scaling functions φ (x) and wavelets ψ (x) are translations and dilatations of special functions. The interpolation subdivision scheme is used in the process of the scaling j +1 function construction. In the scheme, the scaling function values at even collocations x2k are equivalent to the function j +1

j

values at the collocations xk , and the values at the odd collocations x2k+1 are computed by means of Lagrange polynomial interpolation of degree p − 1 [44]. And the wavelets are defined as dilatations and translations of the special scaling function, ψ (x) = φ (2x − 1), with modified boundary wavelets used near the edges of a finite domain [27]. So the wavelet transform of the function u ( x) is the process of finding the scaling and wavelet coefficients. There is a fast wavelet transform (FWT) to perform this process [45]. In the FWT, the Lagrange polynomial interpolation is used hierarchically. The wavelet coefficients are calculated as the error between the function values and its predicted values by the local polynomial interpolation of the information on the next coarser level as showed in Fig. 2, j +1

j

dk = u 2k+1 −

j

r =max χk



p

j

hr u r



j

  j

uk = u xk

(28)

,

j r =min χk

p

j

where hr are the interpolation coefficients of degree p − 1, and χk is a special index set which is related to the polynomial degree. The details regarding the FWT, the inverse fast wavelet transform (IFWT), the computation of interpolation j coefficients and the definition of χk are given in Ref. [45]. The 1-D wavelet construction can be extended to multi-dimension by mean of tensor products [46,47]. For the 2-D case, the wavelet functions consist of the three different product functions as follows:

30

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

j Fig. 2. The calculation of wavelet coefficient dk for p = 4. Note that the filled circular points are in V j , and the filled square point is in V j +1 .

⎧ j j ψ (x1 ) φk2 (x2 ) η = 1 ⎪ ⎪ ⎨ k1 j ,η ψk1 ,k2 (x) = φkj (x1 ) ψkj (x2 ) η = 2 1 2 ⎪ ⎪ ⎩ j j ψk1 (x1 ) ψk2 (x2 ) η = 3 j

j

(x = (x1 , x2 )),

j

(29)

j

j

j

j

and the scaling function φk ,k (x) = φk (x1 ) φk (x2 ), where ψk (x1 ), ψk (x2 ), φk (x1 ), φk (x2 ) correspond to 1-D wavelets 1 2 1 2 1 2 1 2 and scaling functions. The n-dimensional wavelet functions consist of 2n − 1 distinctive product function families. Accordingly, the spaces which the multi-dimensional tensor product second generation wavelets are constructed on are also formed by a product of 1-D nested spaces:





j

V j = xk ∈  : k ∈ K j , k = (k1 , k2 , · · · , kn ) ,

(30) j



j

j

where K j is some index set associated with scaling functions of level j and the collocations points xk = x1,k , x2,k , · · · , j



1

2

xn,k [48]. Because of the nestedness of 1D spaces, the n-dimensional spaces are also nested, i.e. V j +1 = V j ⊕ W j , where W j n is the complement of V j in V j +1 . The FWT in n-dimensions consists of n sequential applications of the 1-D FWT from x1 direction, while the IFWT consists of the n sequential applications of 1-D IFWT in reverse order starting from xn direction [47,48]. Similarly to 1-D case, a function u (x) of level J in n-dimensions can be decomposed as J 2 −1    n



u (x) =

j

j

c k0 φk0 (x) +

j= j0

k∈K j 0 j ,η

j ,η

j ,η

d l ψl

(x),

(31)

η=1 l∈Lη j η

j

in which ψk (x) and φk (x) are respectively n-dimensional wavelets and scaling functions, L j is the index set associated with wavelets of family η of level j and j 0 is the coarsest level. 3.2. Collocations adaptation According to Eq. (27), the wavelets coefficients reflect how large function values deviates from the values calculated by the local interpolation polynomial. So the amplitudes of wavelets coefficients provide a perfect tool to measure the local j j +1 smoothness. If the coefficient dk is small enough, the function value u 2k+1 can be calculated by the local interpolation polynomial with a small enough error. The function u can be compressed through discarding the coefficients which are smaller than the threshold value ε . The function u can be compressed as



u ε (x) =

J 2 −1   n

j

j

c k0 φk0 (x) +

k∈K j 0



j ,η

d   l  η  j ,η  j = j 0 η =1 l∈L j :dl ≥ε 

j ,η

ψl

(x).

The approximation error resulting from neglecting the wavelet coefficients is of the order of

u (x) − u ε (x) ∞ ≤ C 1 , ε

(32)

ε [49]: (33)

where C 1 is a constant. Meanwhile, the points corresponding to the wavelets whose coefficients are smaller than the threshold value will be also removed, and a sparse collocation grid V E of essential points can be defined as



j

j

   j



V E = xk0 , xk : j > j 0 ∧ dk  ≥ ε .

(34)

It must be pointed out that the wavelets coefficients of the points associated with wavelets whose coefficients are lower than the threshold ε in the present iteration may be greater than the threshold ε in the next iteration. So it is suggested by Liandrat and Tchamitchian [50] that some points which do not belong to V E should be add to V E . And these points are j j called neighboring points. The set of neighboring points N k associated with xk is defined as: j



j +1

N k = xr

 : r − 2k ∞ ≤ 2P ,

(35)

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

31

where an integer P is a parameter determining the size of region to be add [27,39]. And there is no benefit for using P > 1 [27–30,39,47]. The one-to-one correspondence between grid points and wavelets or the scaling function allows the grid points to be j j j adaptive to the solution by adding or removing wavelets. From Eq. (28) and Fig. 2, the function values at xk−1 , xk , xk+1 and j

j

xk+2 are requested to calculate dk . This means that performing the FWT and IFWT requires each point to have complete interpolation stencils, which is called a minimum index set [27]. Usually, the sparse grid doesn’t meet the minimum index set requirements. Some points called reconstruction points whose corresponding wavelet amplitudes are smaller than threshold value will be added to the sparse grid, and these wavelet amplitudes will be assigned to be equal to zero. The FWT and IFWT in n-dimensions consist of d applications of n sequential applications of the 1-D FWT and IFWT. So it is also need to add some non-essential points to the sparse grid to ensure each point in every application of 1-D FWT has complete interpolation stencils. In this paper, the AWCM is incorporated into the finite-difference LBM. So it needs to calculate spatial derivatives. Based on the interpolating properties of second generation wavelets, a procedure of the derivative approximation using finitedifferences was discussed in detail by Vasilyev and Bowman [28] for the 1-D case and by Vasilyev [47] for multi-dimensions. j +1 j If the wavelets amplitudes of points x(2k ±1,2k ±1) and all the neighboring points of a grid point xk are smaller than the 1

2

j

threshold value, the actual function can be well approximated by a wavelet interpolant based on c m corresponding to the j points, called ghost points, belonging to some neighborhood of xk [30]. And the approximation error is of the order of ε . So the spatial derivatives can be calculated by using finite differences based on the associated Lagrange interpolation polynomial. The order of accuracy of the finite-difference derivative approximation can be estimated analytically, which has been given by O.V. Vasilyev et al. [28,47]. And the numerical tests results agree well with the analytical estimation. The details of the accuracy of the finite-difference derivative approximation are also given by Vasilyev et al. 3.3. The coupling algorithm of the AWCM and the finite-difference LBM It can be seen from the calculating procedure of the finite-difference LBM, the finite-difference operations only exist in the streaming sub-step. The calculation of the collision sub-step at a grid point just needs the information at this point, not relating to other points. So the AWCM is mainly used in the streaming sub-step. Meanwhile the omitting and adding points occur. The calculation of macroscopic variables and collision sub-step are performed on the adaptive grid points generated in the streaming sub-step. The first step of calculating will be performed on an uniform grid of the finest level J max . In the streaming sub-step of the first step, the density functions, f α , will be compressed based on the wavelet transform according to the threshold value, and the sparse grid V E of essential points will be also generated. To ensure the accuracy of the IBM, the neighboring Eulerian points of all Lagrange boundary points are not compressed, and the resolution levels of the neighboring Eulerian points remain the maximum value. Then the neighboring points of all points belonging to V E are added. The points of the union of V E and the set of the neighboring points may become essential points in the next step, which means that the calculating grid of the next step can only consist of these points. However, it must be pointed that the spares grid which the FWT or IFWT will be performed on must satisfy the minimum index set. Since the FWT and IFWT are used in every step, the calculating grid of every step must satisfy the minimum index set. The reconstruction points have to be also added to the calculating of the next step. And the calculating grid is called V C here. Lastly, the ghost points needed for differentiations are added and some other reconstruction points have to be also added. Here, it must be pointed that calculation of the spatial derivatives is just needed at points of V C . The calculation of macroscopic variables are also performed only on the V C . In summary, the general steps of the algorithm are as follows: (1) Arrange initial values, and compute the matrix A and get A−1 . (2) Use Eq. (9) to obtain the density distribution function after the collision sub-step of the t = tn step on the points of V C (initially setting F α , and assigning V C = V J max ). (3) Perform FWT on density distribution functions f α , and omit grid points according to threshold value (the neighboring Eulerian points are not omitted). (4) Add neighboring points and reconstruction points to generate V C for next step. (5) Add ghost points and some other reconstruction points and perform IFWT to obtain the reconstructed function values. (6) Calculate the spatial derivatives on the points of V C . (7) Use Eq. (7) to obtain the density distribution function after the streaming sub-step of the t = tn step on the points of VC. (8) Use Eq. (12) and Eq. (13) to compute the macroscopic variables on the points of V C . (9) Solve equation Eq. (22) to obtain the velocity corrections at all boundary points and use Eq. (18) to get velocity corrections at BNEP. (10) Correct the fluid velocity at BNEP using Eq. (14) and compute the equilibrium distribution function using Eq. (3). (11) Repeat step 2 to step 10 until the convergence is reached.

32

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Fig. 3. Flow past a fixed cylinder: computational domain and boundary conditions.

4. Numerical results and discussions To verify the present adaptive finite difference LBM based on the AWCM, test cases involving incompressible viscous flows past single fixed circular cylinder are chosen as numerical experiments to validate the present method for fixed boundary flows. Then simulations of a rigid circular cylinder oscillating in a viscous fluid at rest and flows past a rigid oscillating cylinder are conducted to validate the present method for moving boundary flows. Some existing results will serve as references for the comparison. The computational efficiency is also given in this section. In the present numerical experiments, the Reynolds number (Re) is defined as

Re =

U∞ D

υ

(36)

,

where U ∞ is the free stream velocity and D is the diameter of the cylinder and The drag coefficient of cylinder is defined as

Cd =

2F D 2 D ρU∞

υ is the kinematic viscosity of the fluid. (37)

,

where F D is the drag force, which can be calculated by

FD = −









F x XlB s.

(38)

l



F x X lB is the x component of boundary force in Eq. (25). The lift coefficient is defined as

Cl =

2F L 2 D ρU∞

(39)

,

where F L is the lift force, which can be calculated by

FL = −











F y XlB s.

(40)

l

F y X lB is the x component of the boundary force in Eq. (25). For the unsteady flow, the Strouhal number is defined as

St =

fq D U∞

(41)

,

where f q is the vortex shedding frequency. For boundary oscillating in a fluid at rest, the Keulegan–Carpenter number K c is defined as

Kc =

V max T D

,

(42)

where V max is the amplitude of the boundary velocity and T is the period of the oscillation. 4.1. Flow past a fixed cylinder In all simulations of flow past a fixed cylinder, the density ρ of the fluid is set as 1.0. The free stream velocity U ∞ is 0.1. The diameter D of the cylinder is 40. The finest resolution level J max is set as 10, and the coarsest level J min is set as 7. The simulations occupy a square domain with the length L of 1024, i.e. the finest resolution level J max is set as 10. The number of resolution levels is limited to 4, i.e. J min = 7. The finest grid space x and y are set as 1. The center of the cylinder is located at ( L /2, L /2), as shown in Fig. 3.

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

33

Table 1 Drag coefficient, length of bubbles and separation angle for flows past a cylinder at Re = 20 and Re = 40. Case

References

Cd

Ls

θs

Re = 20

Niu et al. [51] Shukla et al. [52] J. Wu and C. Shu [37] Present

2.144 2.07 2.091 2.11

0.945 0.92 0.93 0.93

42.9◦ 43.2◦ – 42.8◦

Re = 40

Niu et al. [51] Shukla et al. [52] J. Wu and C. Shu [37] Present

1.589 1.55 1.539 1.56

2.26 2.34 2.23 2.20

53.86◦ 52.7◦ – 53.0◦

Fig. 4. Steady flow past a cylinder at Re = 20 and Re = 40. The resolution level of the unshown computational area is the coarest level J min .

34

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Table 2 Drag coefficient, length of bubbles and separation angle for flows past a cylinder at Re = 100 and Re = 200. (Also Re = 185.) Case

References

Cd

Cl

St

Re = 100

Benson et al. [53] Ding et al. [54] J. Wu and C. Shu [37] Present

1.46 ± 0.01 1.325 ± 0.008 1.334 ± 0.012 1.358 ± 0.013

±0.38 ±0.28 ±0.37 ±0.342

0.17 0.164 0.163 0.167

Re = 200

Benson et al. [53] Ding et al. [54] J. Wu and C. Shu [37] Present

1.45 ± 0.04 1.327 ± 0.045 1.43 ± 0.051 1.340 ± 0.049

±0.65 ±0.60 ±0.75 ±0.655

0.193 0.196 0.195 0.196

Re = 185

Present

1.325 ± 0.044

±0.633

0.192

Fig. 5. Evolution of drag and lift coefficients at Re = 100 and Re = 200.

For the LBM, the time step t is set as 1. At the start, the density functions values are equal to the equilibrium distribution functions values. Also, the equilibrium distribution functions are used to implement the boundary conditions at the far field boundaries. For the AWCM, the order of the wavelet basis p is 4 and the threshold ε = 10−4 . Both steady and unsteady state experiments are conducted. For the steady state experiments, the Re is set as 20 and 40. For the unsteady state experiments, the Re is taken as 100 and 200. For the case of Re = 20 and 40, the flows present a steady state. Also, a pair of stationary recirculating eddies will develop behind the cylinder. When the Reynolds number increases, the length of the recirculation region from the rearmost point of the cylinder to the end of the wake will also increase. The comparison of drag coefficient C d , length of the recirculation region L (scaled by D /2) and separation angle θs with previous studies are also shown in Table 1. The contour plots of velocities, streamlines and the final collocation points are shown in Fig. 4. From the Table 1, it can be seen that the present numerical results agree well with those in the previous studies. From Fig. 4, the region of the recirculation becomes bigger when the Reynolds number grows from 20 to 40. Meanwhile, the finer collocation points region also becomes bigger according to the flow field, which means that the AMR can distribute the collocation points in accordance with the flow field for the steady state cases. Also, the refinement criterion based on the wavelet amplitudes of distribution functions is of efficiency for the steady state cases. For the case Re = 100 and Re = 200, the flow becomes unsteady. There will be a Kármán vortex street behind the cylinder. Also, there will be a lift force on the cylinder. Table 2 shows the comparison of drag coefficient, lift coefficient and Strouhal number for flows past a cylinder at Re = 100 and Re = 200 with previous studies. And the time evolution of the drag and lift coefficients are shown in the Fig. 5. Besides, the drag coefficient, lift coefficient and Strouhal number for flows past a cylinder at Re = 185 are also showed in Table 2 which will be used in the simulations of flows past a rigid oscillating cylinder. From Table 2, it can be found that a good agreement is obtained between the present results and the data in previous literatures. Because the flow will be in an unsteady state, the distribution of the collocation points should be dynamic according to the flows. The adaptive collocation point distributions in one cycle for Re = 100 and Re = 200 are plotted in Fig. 6 and Fig. 8 respectively. The corresponding instantaneous velocity contours are shown in Fig. 7 and Fig. 9. In the plots of collocation point distributions, the process of point omission and addition is clearly shown. From the plots of point distributions and velocity contours, it can be clearly obtained that the AWCM can distribute the collocation points

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

35

Fig. 6. Adaptive collocation point distribution in one cycle for Re = 100 at t 0 , t 0 + 1/(4 f s ), t 0 + 2/(4 f s ) and t 0 + 3/(4 f s ). The resolution level of the unshown computational area is the coarest level J min .

efficiently according to the flow field information. Since it is difficult, even impossible, to pre-distribute the mesh points for the calculation of continuously developing flows, it is of great use to distribute the mesh points according to the flow field information. It is well known that the vortices value and the area of vortices behind the cylinder will increase with augmentation of Reynolds number, which means that it needs more grid refinement for a higher Reynolds number. By comparing Fig. 8 with Fig. 6, it can be found that the area of higher resolution level collocation points for Re = 200 is larger than that for Re = 100. It can be also found that it needs some mesh refinement on the right boundary for Re = 200. This is caused by that the area of vortices for Re = 200 is large enough which leads to that there are rapidly varying regions near the right boundary for the density functions. 4.2. Oscillating cylinder in a viscous flow at rest In simulations of an oscillating cylinder in a viscous flow at rest, the density of the fluid, finest resolution level, the coarsest level, the finest grid space x and y, the time step and the threshold value are same as in the simulations of

36

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Fig. 7. Velocity contours and Streamlines in one cycle for Re = 100 at t 0 , t 0 + 1/(4 f s ), t 0 + 2/(4 f s ) and t 0 + 3/(4 f s ).

flows past a fixed cylinder. The diameter D of the cylinder is 20. The simulations occupy a square domain. Initially, the center of the cylinder is located at ( L /2, L /2), as shown in Fig. 10. The cylinder oscillates only in the vertical direction. And the velocity of the rigid cylinder is described as:



V (t ) = V max cos

2π t T



,

(43)

where V max is the amplitude of the velocity and T is the oscillating period. In the present simulation, the amplitude of the velocity is 0.04, and oscillating period is 2500. The Reynolds number Re is chosen as 100. Hence the Keugelan–Carpenter number is 5. Besides, it can be found that the fluid is filled both inside and outside the boundary in the LBM. The force, acting on the boundary, comes from the fluid both inside and outside the boundary. As pointed by K. Suzuki and T. Inamuro [55], the flow and the pressure outside the boundary are not affected by the internal fluid. However, when calculating the force and the torque acting on the boundary, the influence of the internal fluid is not negligible. The effect of inertial force can be calculated by the time derivative of the linear momentum of the internal fluid [5,55], which can be defined as

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

37

Fig. 8. Adaptive collocation point distribution in one cycle for Re = 200 at t 0 , t 0 + 1/(4 f s ), t 0 + 2/(4 f s ) and t 0 + 3/(4 f s ). The resolution level of the unshown computational area is the coarsest level J min .

Finer (t ) = ρ f V B

∂ U B (t ) , ∂t

(44)

where ρ f is the density of the fluid, V B is the volume of the boundary and U B (t ) is the velocity of the rigid boundary. From Eq. (44), when the boundary is fixed, the effect of the inertial force can be negligible. The evolutions of the drag coefficient and the drag coefficient accounting for the inertial fluid are showed in Fig. 11. Also, the reference values given by A. De Rosis et al. [56] are also showed in the figure. It can be seen that the present numerical results agree well with the reference results. From the evolutions of the drag coefficients, the drag coefficient decreases when the inertial fluid is taken into consideration. However, the oscillation period remains unchanged, which means the oscillation period is independent. Besides, whether or not the inertial fluid is taken into consideration, the drag coefficient is C d = 2 at t / Z = z (z is an integer) and C d = −2 at t / T = z + 1/2 (z is an integer), as highlighted by the black and red dashed lines in Fig. 11. The adaptive collocation point distributions in one cycle are shown in Fig. 12. And the corresponding instantaneous velocity contours are shown in Fig. 13. In Fig. 12 and Fig. 13, t 0 is the time when a oscillating cycle begins, i.e. t 0 = 2zTπ (z is

38

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Fig. 9. Velocity contours and Streamlines in one cycle for Re = 200 at t 0 , t 0 + 1/(4 f s ), t 0 + 2/(4 f s ) and t 0 + 3/(4 f s ).

Fig. 10. Oscillating cylinder in a viscous flow at rest: computational domain and boundary conditions.

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

39

Fig. 11. Oscillating cylinder in a viscous flow at rest: evolutions of the drag coefficient C d and the drag coefficient accounting for the inertial fluid effect C d∗ . The blue circular and red square points are the reference values given in the Ref. [56]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 12. Velocity contours and Streamlines in one cycle for oscillating cylinder in a viscous flow at rest at t 0 , t 0 + 0.25T , t 0 + 0.50T and t 0 + 0.75T . The resolution level of the unshown computational area is the coarest level J min .

40

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Fig. 13. Velocity contours and Streamlines in one cycle for oscillating cylinder in a viscous flow at rest at t 0 , t 0 + 0.25T , t 0 + 0.50T and t 0 + 0.75T .

Fig. 14. Flows past a rigid oscillating cylinder: computational domain and boundary conditions.

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

41

Fig. 15. Comparison of Cd(Present ), Cdrmse and Clrmse for flows past a rigid oscillating cylinder at Re = 185.

Fig. 16. Drag and lift coefficients for flows past a oscillating rigid cylinder at Re = 185.

an integer). From the plots of point distributions and velocity contours, it can be found that the AWCM can distribute the collocation points efficiently in the simulations involving moving boundaries. 4.3. Flows past an oscillating rigid cylinder In simulations of flows past a oscillating rigid cylinder, except that the cylinder is oscillating and the Reynolds number is 185, the computational domain, boundary conditions, the time step and the values used in the AWCM are as same as in

42

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Fig. 17. Adaptive collocation point distribution in one cycle for f os / f o = 1.0 at t 0 , t 0 + 1/(4 f os ), t 0 + 2/(4 f os ) and t 0 + 3/(4 f os ). The resolution level of the unshown computational area is the coarsest level J min .

the simulations of flows past a fixed cylinder in the present article (see Fig. 14). The motion of the oscillating cylinder is defined as:

y (t ) = A sin (2π f os t ) ,

(45)

where A is the amplitude of the oscillating and f os is the oscillating frequency. In the present simulation, the amplitude of the oscillating is A = 0.2D. Five cases with different oscillating frequency of f os / f o = 0.8, 0.9, 1.0, 1.1 and 1.2 are conducted, where f o is the natural vortex shedding frequency for flows past a fixed cylinder at Re = 185 which has been given in Table 2. The time-averaged drag coefficient Cd, the root-mean-square-error values of drag coefficient Cdrmse and lift Clrmse are shown in Fig. 15. As same as in the simulation of oscillating cylinder in a viscous flow at rest, the effect of inertial force is taken into account. The results of Y. Wang et al. [5] and E. Guilmineau and P. Queutey [57] are also given for comparison. It can be seen that the present results agree well with the reference data.

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

43

Fig. 18. Velocity contours and Streamlines in one cycle for oscillating cylinder in a viscous flow at rest for f os / f o = 1.0 at t 0 , t 0 + 1/(4 f os ), t 0 + 2/(4 f os ) and t 0 + 3/(4 f os ).

The evolutions of the Cd, Cdrmse and Clrmse at f os / f o = 0.9, 1.0, 1.1, 1.2 are given in Fig. 16. It can be seen that the drag and lift coefficients exhibit harmonic oscillations and the magnitudes grow as f os / f o increases, when f os / f o ≤ 1.0. When f os / f o > 1.0, the lift and drag coefficients exhibit a modulation phenomenon. The adaptive collocation point distributions in one cycle for f os / f o = 1.0 are plotted in Fig. 17. The corresponding instantaneous velocity contours is shown in Fig. 18. In Fig. 17 and Fig. 18, t 0 is the time when an oscillating cycle begins, i.e. t 0 = 2πzf (z is an integer). Also, the collocation point distributions when the oscillating cylinder reaches its upper extreme os position are plotted in Fig. 19 with f os / f o = 0.8, 0.9, 1.1, 1.2. And the corresponding streamlines and velocity contours are shown in Fig. 20. From the plots of point distributions and the corresponding velocity contours, the AWCM can distribute the collocation points efficiently according to the flows field information, which indicates that the present algorithm can handle the problems involving moving boundaries. To demonstrate the advantage of the adaptive LBM presented in this paper, we compare the computing time and the number of mesh points in our AWCM-IBM-FDLBM code with that in the uniform-mesh code. All the simulations are per-

44

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Fig. 19. Adaptive collocation point distribution with the cylinder locating at its extreme upper position at f os / f o = 0.8, f os / f o = 0.9, f os / f o = 1.1 and f os / f o = 1.2. The resolution level of the unshown computational area is the coarest level J min .

formed by a single core on a PC with Intel(R) Core(TM) i7-4790K 4.00GHz CPU and 16.0 GB RAM. In Table 3 and Table 4, the computing time for per step and the number of collocation points of our AWCM-IBM-FDLBM code are compared with the uniform-mesh code for fixed boundaries cases and moving boundaries cases. Both in the fixed boundary cases and the moving boundary cases, it can be seen that our proposed method is able to reduce the computation time for each single step and that the required points in AWCM-IBM-FDLBM code is much less than the uniform-mesh code from Table 3 and Table 4. Note that for the unsteady cases of the fixed boundary cases and the moving boundary cases, the distribution of collocation points is dynamic, so the computing time of per step and the number of the collocation points are changing all the time. It is reasonable to study the range of the computing time and the point number during one cycle of vortex shedding. For the fixed boundary cases, the computing time of per step and the number of required points become larger as the Reynolds number increases. For the flows past an oscillating rigid cylinder cases, the computing time of per step and the number of required points become slightly larger as the oscillating frequency increases.

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

45

Fig. 20. Instantaneous streamlines with the cylinder locating at its extreme upper position at f os / f o = 0.8, f os / f o = 0.9, f os / f o = 1.1 and f os / f o = 1.2.

Table 3 Comparison of the computing time of per step. Case study

Uniform-mesh

Adaptive-mesh

Speedup

Re = 20 Re = 40 Re = 100 Re = 200 Oscillating cylinder in a viscous flow at rest Flows past an f os / f o = 0.8 oscillating f os / f o = 0.9 rigid f os / f o = 1.0 f os / f o = 1.1 cylinder f os / f o = 1.2

1.0 1.0 1.0 1.0 1.01

0.2184 0.2418 0.2652–0.2808 0.3432–0.3743 0.2495–0.4055

4.579 4.136 3.561–3.771 2.672–2.914 2.466–4.008

1.01 1.01 1.01 1.01 1.01

0.3432–0.3744 0.3432–0.3587 0.3432–0.3588 0.3584–0.3896 0.3590–0.3900

2.698–2.943 2.816–2.943 2.815–2.943 2.592–2.818 2.564–2.786

Flows past a fixed cylinder

46

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

Table 4 Comparison of the number of collocation points. Case study

Uniform-mesh

Adaptive-mesh

Speedup

Re = 20 Re = 40 Re = 100 Re = 200 Oscillating cylinder in a viscous flow at rest Flows past an f os / f o = 0.8 f os / f o = 0.9 oscillating f os / f o = 1.0 rigid f os / f o = 1.1 cylinder f os / f o = 1.2

1025 × 1025 1025 × 1025 1025 × 1025 1025 × 1025 1025 × 1025

21533 23224 26758–27189 35952–37004 20008–20072

2.05% 2.23% 2.55%–2.59% 3.42%–3.52% 1.90%–1.91%

1025 × 1025 1025 × 1025 1025 × 1025 1025 × 1025 1025 × 1025

35461–36503 35817–36555 35820–36959 36236–36981 36320–37096

3.38%–3.47% 3.41%–3.48% 3.41%–3.52% 3.45%–3.52% 3.46%–3.53%

Flows past a fixed cylinder

5. Conclusions A novel adaptive collocation points lattice Boltzmann method is proposed in this work. In the proposed method, the adaptive wavelet collocation method is firstly, to the best of our knowledge, incorporated into the finite-difference lattice Boltzmann and the immersed boundary is also incorporated to eliminate the special treatments in using the FD-based derivative approximation near boundaries. A new grid refinement criterion based on the wavelet amplitudes of distribution functions is proposed. To justify the present AWCM-IB-LBM approach, a series of test cases are investigated, including fixed boundary cases and moving boundary cases. For the fixed boundary cases, flows past a fixed circular cylinder at Re = 20, Re = 40 for steady problems and Re = 100, Re = 200 for unsteady problems are carried out. For the moving boundary cases, a rigid circular cylinder oscillating with Re = 100 and K c = 5 in a viscous fluid at rest as well as flows past a rigid oscillating cylinder with five different oscillating frequency of f os / f o = 0.8, 0.9, 1.0, 1.1 and 1.2 at Re = 185 are chosen as the experiments. Numerical results obtained in this work are discussed and compared with available data in previous literatures. A good agreement between the present results and the data in previous literatures is achieved. The required points of the AWCMIBM-FDLBM code are much less than those of the uniform mesh code, which means the computing resources are much saved. As a consequence, the computing time of per step is reduced, which means that the present adaptive approach can improve the computational efficiency. It is reasonable to conclude that the proposed method is an efficient approach for simulation of fluid–solid interaction problems involving fixed boundaries and moving boundaries. Although the current AWCM-IB-LBM is developed in 2D, it is straightforward to extend this method to 3D based on the 3D second generation wavelet transform. Just replace D2Q9 lattice with D3Q19 or D3Q27. Also the current method solves the single phase flow problems. The extension to multiphase flow problems is another promising direction. Besides, the LBM has been used to treat other problems, such as shallow water flows [58], heat conduction problems [59] and wave propagations [60]. It is believed that this method can be extended to those fields. In a word, this work is still in the process. Acknowledgements The authors would like to thank Prof. Wu for his greatly useful advice about IB. The authors would also like to acknowledge the support of the National Natural Science Foundation of China (51479041, 51309365, 51279038, U1430236). References [1] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1) (1998) 329–364, http://dx.doi.org/10.1146/annurev.fluid. 30.1.329. [2] D. Yu, R. Mei, L.-S. Luo, W. Shyy, Viscous flow computations with the method of lattice Boltzmann equation, Prog. Aerosp. Sci. 39 (5) (2003) 329–367, http://dx.doi.org/10.1016/S0376-0421(03)00003-4, http://www.sciencedirect.com/science/article/pii/S0376042103000034. [3] O. Malaspinas, N. Fiétier, M. Deville, Lattice Boltzmann method for the simulation of viscoelastic fluid flows, J. Non-Newton. Fluid Mech. 165 (23–24) (2010) 1637–1653, http://dx.doi.org/10.1016/j.jnnfm.2010.09.001, http://www.sciencedirect.com/science/article/pii/S0377025710002429. [4] E. Vergnault, P. Sagaut, An adjoint-based lattice Boltzmann method for noise control problems, J. Comput. Phys. 276 (2014) 39–61, http://dx.doi.org/ 10.1016/j.jcp.2014.07.027, http://www.sciencedirect.com/science/article/pii/S0021999114005208. [5] Y. Wang, C. Shu, C. Teo, J. Wu, An immersed boundary-lattice Boltzmann flux solver and its applications to fluid structure interaction problems, J. Fluids Struct. 54 (2015) 440–465, http://dx.doi.org/10.1016/j.jfluidstructs.2014.12.003, http://www.sciencedirect.com/science/article/pii/S0889974614002709. [6] J. Wu, C. Shu, A solution-adaptive lattice Boltzmann method for two-dimensional incompressible viscous flows, J. Comput. Phys. 230 (6) (2011) 2246–2269, http://dx.doi.org/10.1016/j.jcp.2010.12.013, http://www.sciencedirect.com/science/article/pii/S0021999110006777. [7] X. He, L.-S. Luo, M. Dembo, Some progress in lattice Boltzmann method. Part I. Nonuniform mesh grids, J. Comput. Phys. 129 (2) (1996) 357–363, http://dx.doi.org/10.1006/jcph.1996.0255, http://www.sciencedirect.com/science/article/pii/S0021999196902557. [8] X. He, G. Doolen, Lattice Boltzmann method on curvilinear coordinates system: flow around a circular cylinder, J. Comput. Phys. 134 (2) (1997) 306–315, http://dx.doi.org/10.1006/jcph.1997.5709, http://www.sciencedirect.com/science/article/pii/S0021999197957090. [9] T. Imamura, K. Suzuki, T. Nakamura, M. Yoshida, Acceleration of steady-state lattice Boltzmann simulations on non-uniform mesh using local time step method, J. Comput. Phys. 202 (2) (2005) 645–663, http://dx.doi.org/10.1016/j.jcp.2004.08.001, http://www.sciencedirect.com/science/article/pii/ S0021999104003006.

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

47

[10] C. Shu, X.D. Niu, Y.T. Chew, Taylor series expansion and least squares-based lattice Boltzmann method: three-dimensional formulation and its applications, Int. J. Mod. Phys. C 14 (07) (2003) 925–944, http://dx.doi.org/10.1142/S0129183103005133, http://www.worldscientific.com/doi/pdf/10.1142/ S0129183103005133, http://www.worldscientific.com/doi/abs/10.1142/S0129183103005133. [11] X. Niu, Y. Chew, C. Shu, Simulation of flows around an impulsively started circular cylinder by Taylor series expansion- and least squares-based lattice Boltzmann method, J. Comput. Phys. 188 (1) (2003) 176–193, http://dx.doi.org/10.1016/S0021-9991(03)00161-X, http://www.sciencedirect.com/science/ article/pii/S002199910300161X. [12] M.B. Reider, J.D. Sterling, Accuracy of discrete-velocity BGK models for the simulation of the incompressible Navier–Stokes equations, Comput. Fluids 24 (4) (1995) 459–467, http://dx.doi.org/10.1016/0045-7930(94)00037-Y, http://www.sciencedirect.com/science/article/pii/004579309400037Y. [13] N. Cao, S. Chen, S. Jin, D. Martínez, Physical symmetry and lattice symmetry in the lattice Boltzmann method, Phys. Rev. E 55 (1997) R21–R24, http:// dx.doi.org/10.1103/PhysRevE.55.R21, http://link.aps.org/doi/10.1103/PhysRevE.55.R21. [14] R. Mei, W. Shyy, On the finite difference-based lattice Boltzmann method in curvilinear coordinates, J. Comput. Phys. 143 (2) (1998) 426–448, http:// dx.doi.org/10.1006/jcph.1998.5984, http://www.sciencedirect.com/science/article/pii/S0021999198959848. [15] T. Lee, C.-L. Lin, A characteristic Galerkin method for discrete Boltzmann equation, J. Comput. Phys. 171 (1) (2001) 336–356, http://dx.doi.org/10.1006/ jcph.2001.6791, http://www.sciencedirect.com/science/article/pii/S0021999101967919. [16] Y. Li, E.J. LeBoeuf, P.K. Basu, Least-squares finite-element scheme for the lattice Boltzmann method on an unstructured mesh, Phys. Rev. E 72 (2005) 046711, http://dx.doi.org/10.1103/PhysRevE.72.046711, http://link.aps.org/doi/10.1103/PhysRevE.72.046711. [17] A. Duster, L. Demkowicz, E. Rank, High-order finite elements applied to the discrete Boltzmann equation, Int. J. Numer. Methods Eng. 67 (8) (2006) 1094–1121, http://dx.doi.org/10.1002/nme.1657. [18] S. Succi, G. Amati, R. Benzi, Challenges in lattice Boltzmann computing, J. Stat. Phys. 81 (1) (1995) 5–16, http://dx.doi.org/10.1007/BF02179964. [19] G. Amati, S. Succi, R. Benzi, Turbulent channel flow simulations using a coarse-grained extension of the lattice Boltzmann method, Fluid Dyn. Res. 19 (5) (1997) 289–302, http://dx.doi.org/10.1016/S0169-5983(96)00026-3, http://www.sciencedirect.com/science/article/pii/S0169598396000263. [20] H. Xi, G. Peng, S.H. Chou, Finite-volume lattice Boltzmann method, Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat. Interdiscip. Topics 59 (5 Pt B) (1999) 6202–6205. [21] S. Ubertini, G. Bella, S. Succi, Unstructured lattice Boltzmann equation with memory, in: Discrete Simulation of Fluid Dynamics in Complex Systems, Math. Comput. Simul. 72 (2006) 237–241, http://dx.doi.org/10.1016/j.matcom.2006.05.009, http://www.sciencedirect.com/science/article/pii/ S0378475406001406. [22] H. Ding, C. Shu, A stencil adaptive algorithm for finite difference solution of incompressible viscous flows, J. Comput. Phys. 214 (1) (2006) 397–420, http://dx.doi.org/10.1016/j.jcp.2005.09.021, http://www.sciencedirect.com/science/article/pii/S0021999105004511. [23] S.M. Guzik, T.H. Weisgraber, P. Colella, B.J. Alder, Interpolation methods and the accuracy of lattice-Boltzmann mesh refinement, J. Comput. Phys. 259 (2014) 461–487, http://dx.doi.org/10.1016/j.jcp.2013.11.037, http://www.sciencedirect.com/science/article/pii/S0021999113008012. [24] A. Fakhari, T. Lee, Finite-difference lattice Boltzmann method with a block-structured adaptive-mesh-refinement technique, Phys. Rev. E, Stat. Nonlin. Soft Matter Phys. 89 (3) (2014) 033310. [25] T. Lee, C.-L. Lin, An Eulerian description of the streaming process in the lattice Boltzmann equation, J. Comput. Phys. 185 (2) (2003) 445–471, http://dx.doi.org/10.1016/S0021-9991(02)00065-7, http://www.sciencedirect.com/science/article/pii/S0021999102000657. [26] X. Guo, J. Yao, C. Zhong, J. Cao, A hybrid adaptive-gridding immersed-boundary lattice Boltzmann method for viscous flow simulations, in: The Fourth European Seminar on Computing (ESCO 2014), Appl. Math. Comput. 267 (2015) 529–553, http://dx.doi.org/10.1016/j.amc.2015.01.082, http://www. sciencedirect.com/science/article/pii/S0096300315001149. [27] S. Paolucci, Z.J. Zikoski, D. Wirasaet, WAMR: an adaptive wavelet method for the simulation of compressible reacting flow. Part I. Accuracy and efficiency of algorithm, J. Comput. Phys. 272 (2014) 814–841, http://dx.doi.org/10.1016/j.jcp.2014.01.025, http://www.sciencedirect.com/science/article/pii/ S0021999114000527. [28] O.V. Vasilyev, C. Bowman, Second-generation wavelet collocation method for the solution of partial differential equations, J. Comput. Phys. 165 (2) (2000) 660–693. [29] J. Regele, O.V. Vasilyev, An adaptive wavelet-collocation method for shock computations, Int. J. Comput. Fluid Dyn. 23 (7) (2009) 503–518. [30] O.V. Vasilyev, N.K.-R. Kevlahan, An adaptive multilevel wavelet collocation method for elliptic problems, J. Comput. Phys. 206 (2) (2005) 412–431, http://dx.doi.org/10.1016/j.jcp.2004.12.013, http://www.sciencedirect.com/science/article/pii/S0021999104005200. [31] H. Yousefi, A. Noorzad, J. Farjoodi, Simulating 2D waves propagation in elastic solid media using wavelet based adaptive method, J. Sci. Comput. 42 (3) (2010) 404–425, http://dx.doi.org/10.1007/s10915-009-9328-7. [32] M. de la Llave Plata, R. Cant, A wavelet-based multiresolution approach to large-eddy simulation of turbulence, J. Comput. Phys. 229 (20) (2010) 7715–7738, http://dx.doi.org/10.1016/j.jcp.2010.06.027, http://www.sciencedirect.com/science/article/pii/S0021999110003426. [33] N.K.R. Kevlahan, O.V. Vasilyev, An adaptive wavelet collocation method for fluid–structure interaction at high Reynolds numbers, SIAM J. Sci. Comput. 26 (6) (2005) 1894–1915, http://dx.doi.org/10.1137/S1064827503428503. [34] S. Charles Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517. [35] Z.-G. Feng, E.E. Michaelides, The immersed boundary-lattice Boltzmann method for solving fluid particles interaction problems, J. Comput. Phys. 195 (2) (2004) 602–628, http://dx.doi.org/10.1016/j.jcp.2003.10.013, http://www.sciencedirect.com/science/article/pii/S0021999103005758. [36] Z.-G. Feng, E.E. Michaelides, Proteus: a direct forcing method in the simulations of particulate flows, J. Comput. Phys. 202 (1) (2005) 20–51, http:// dx.doi.org/10.1016/j.jcp.2004.06.020, http://www.sciencedirect.com/science/article/pii/S0021999104002669. [37] J. Wu, C. Shu, Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications, J. Comput. Phys. 228 (6) (2009) 1963–1979, http://dx.doi.org/10.1016/j.jcp.2008.11.019, http://www.sciencedirect.com/science/article/pii/S0021999108006116. [38] P. Fernandes, P. Girdinio, P. Molfino, G. Molinari, A comparison of adaptive strategies for mesh refinement based on ‘a posteriori’ local error estimation procedures, IEEE Trans. Magn. 26 (2) (1990) 795–798. [39] D. Wirasaet, Numerical Solutions of Multi-Dimensional Partial Differential Equations Using an Adaptive Wavelet Method, Ph.D. thesis, University of Notre Dame, United States, 2007. [40] X. He, L.-S. Luo, Lattice Boltzmann model for the incompressible Navier–Stokes equation, J. Stat. Phys. 88 (3) (1997) 927–944, http://dx.doi.org/10.1023/ B:JOSS.0000015179.12689.e4. [41] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev., E Stat. Nonlin. Soft Matter Phys. 65 (4 Pt 2B) (2002) 046308. [42] W. Sweldens, The lifting scheme: a custom-design construction of biorthogonal wavelets, Appl. Comput. Harmon. Anal. 3 (2) (1996) 186–200. [43] W. Sweldens, The lifting scheme: a construction of second generation wavelets, SIAM J. Math. Anal. 29 (2) (1998) 511–546. [44] G. Deslauriers, S. Dubuc, Symmetric iterative interpolation processes, Constr. Approx. 5 (1) (1989) 49–68, http://dx.doi.org/10.1007/BF01889598. [45] Y.A. Rastigejev, S. Paolucci, Wavelet-based adaptive multiresolution computation of viscous reactive flows, Int. J. Numer. Methods Fluids 52 (7) (2006) 749–784, http://dx.doi.org/10.1002/fld.1202. [46] B.I. Daubechies, Ten lectures on wavelets, in: CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 61, Soc. for Industrial & Applied Math., 2010. [47] O.V. Vasilyev, Solving multi-dimensional evolution problems with localized structures using second generation wavelets, Int. J. Comput. Fluid Dyn. 17 (2) (2003) 151–168, http://dx.doi.org/10.1080/1061856021000011152.

48

X. Cui et al. / Journal of Computational Physics 333 (2017) 24–48

[48] A. Nejadmalayeri, A. Vezolainen, E. Brown-Dymkoski, O.V. Vasilyev, Parallel adaptive wavelet collocation method for PDEs, J. Comput. Phys. 298 (2015) 237–253, http://dx.doi.org/10.1016/j.jcp.2015.05.028, http://www.sciencedirect.com/science/article/pii/S0021999115003629. [49] D.L. Donoho, Interpolating Wavelet Transform, Tech. Rep. 408, Department of Statistics, Stanford University, 1992. [50] J. Liandrat, P. Tchamitchian, Resolution of the 1D Regularized Burgers Equation Using a Spatial Wavelet Approximation, Tech. Rep. NASA-CR-187480, ICASE-90-83, NAS 1.26:187480, AD-A231566, NASA Langley Research Center, 1990. [51] X. Niu, C. Shu, Y. Chew, Y. Peng, A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows, Phys. Lett. A 354 (3) (2006) 173–182, http://dx.doi.org/10.1016/j.physleta.2006.01.060, http://www.sciencedirect.com/science/article/pii/ S0375960106001472. [52] R.K. Shukla, M. Tatineni, X. Zhong, Very high-order compact finite difference schemes on non-uniform grids for incompressible Navier Stokes equations, J. Comput. Phys. 224 (2) (2007) 1064–1094, http://dx.doi.org/10.1016/j.jcp.2006.11.007, http://www.sciencedirect.com/science/article/pii/ S0021999106005663. [53] M. Benson, P. Bellamyknights, J. Gerrard, I. Gladwell, A viscous splitting algorithm applied to low Reynolds number flows round a circular cylinder, J. Fluids Struct. 3 (1989) 439–479, http://dx.doi.org/10.1016/S0889-9746(89)80026-X. [54] H. Ding, C. Shu, K. Yeo, D. Xu, Simulation of incompressible viscous flows past a circular cylinder by hybrid FD scheme and meshless least square-based finite difference method, Comput. Methods Appl. Mech. Eng. 193 (9–11) (2004) 727–744, http://dx.doi.org/10.1016/j.cma.2003.11.002, http://www. sciencedirect.com/science/article/pii/S0045782503005838. [55] K. Suzuki, T. Inamuro, Effect of internal mass in the simulation of a moving body by the immersed boundary method, Comput. Fluids 49 (1) (2011) 173–187, http://dx.doi.org/10.1016/j.compfluid.2011.05.011, http://www.sciencedirect.com/science/article/pii/S0045793011001708. [56] A.D. Rosis, S. Ubertini, F. Ubertini, A partitioned approach for two-dimensional fluid structure interaction problems by a coupled lattice Boltzmannfinite element method with immersed boundary, J. Fluids Struct. 45 (2014) 202–215, http://dx.doi.org/10.1016/j.jfluidstructs.2013.12.009, http://www. sciencedirect.com/science/article/pii/S0889974613002910. [57] E. Guilmineau, P. Queutey, A numerical simulation of vortex shedding from AN oscillating circular cylinder, J. Fluids Struct. 16 (2002) 773–794, http://dx.doi.org/10.1006/jfls.2002.0449. [58] H. Liu, J.G. Zhou, R. Burrows, Lattice Boltzmann simulations of the transient shallow water flows, Adv. Water Resour. 33 (4) (2010) 387–396, http://dx.doi.org/10.1016/j.advwatres.2010.01.005, http://www.sciencedirect.com/science/article/pii/S0309170810000175. [59] H. Rihab, N. Moudhaffar, B.N. Sassi, P. Patrick, Enthalpic lattice Boltzmann formulation for unsteady heat conduction in heterogeneous media, Int. J. Heat Mass Transf. 100 (2016) 728–736, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.05.001, http://www.sciencedirect.com/science/article/pii/ S0017931016307487. [60] S. Xiao, A lattice Boltzmann method for shock wave propagation in solids, Commun. Numer. Methods Eng. 23 (1) (2007) 71–84.