Entropic lattice Boltzmann method for turbulent flow simulations: Boundary conditions

Entropic lattice Boltzmann method for turbulent flow simulations: Boundary conditions

Physica A 392 (2013) 1925–1930 Contents lists available at SciVerse ScienceDirect Physica A journal homepage: www.elsevier.com/locate/physa Entropi...

579KB Sizes 2 Downloads 116 Views

Physica A 392 (2013) 1925–1930

Contents lists available at SciVerse ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Entropic lattice Boltzmann method for turbulent flow simulations: Boundary conditions S.S. Chikatamarla, I.V. Karlin ∗ Department of Mechanical and Process Engineering, ETH Zurich, 8092 Zurich, Switzerland



Article history: Received 16 July 2012 Received in revised form 13 December 2012 Available online 17 January 2013 Keywords: Entropic lattice Boltzmann method (ELBM) Complex boundaries Turbulent flow simulation

abstract We introduce a new concept of boundary conditions for realization of the lattice Boltzmann simulations of turbulent flows. The key innovation is the use of a universal distribution function for particles, analogous to the Tamm–Mott-Smith solution for the shock wave in the classical Boltzmann kinetic equation. Turbulent channel flow simulations demonstrate that the new boundary enables accurate results even with severely under-resolved grids. Generalization to complex boundary is illustrated with an example of turbulent flow past a circular cylinder. © 2013 Elsevier B.V. All rights reserved.

1. Introduction The entropic lattice Boltzmann method (ELBM) was conceived a decade ago as a means for computational fluid dynamics [1]. Unlike its predecessor, the conventional lattice Boltzmann method (LBM) [2], ELBM makes use of second law in the form of the discrete-time H-theorem to compute flows at high Reynolds numbers. Recent studies of ELBM in periodic geometry revealed its robustness, accuracy and efficiency for turbulent flow simulations [3,4]. However, it remains open how to perform high Reynolds number simulations for wall-bounded flows. To that end, there exist a large number of proposals for the boundary conditions for the LBM in the literature, including cases of generic curved walls [2,5,6]. Many of the boundary conditions for LB simulations were studied for laminar flow problems. Moreover, for the case of a flat walls, the candidate that outperforms in terms of accuracy, stability and simplicity is the half-way bounce back scheme. Many boundary conditions recover the accuracy and stability of the bounce back scheme also for the curved walls, with the limiting case accuracy (for flat walls aligned with the grid) still being the bounce back scheme [7]. What is much less discussed is their realizability for turbulent flows where the situation becomes dramatically different (extending the grid requirements from laminar flow simulations requires prohibitively large grids). A telling example is provided by recent LBM simulations of turbulent channel flow for small Reynolds numbers [8,9]. These simulations revealed that (i) the only boundary condition which is able to produce results in a fully resolved simulation is the ancient bounce back scheme, and (ii) at higher Reynolds numbers, application of the bounce back fails as it triggers discontinuities (shocks) which propagate from the walls into the bulk. Thus, the major problem of useful boundary conditions for LBM remains open for turbulent flows. In this paper, we revisit the problem of realization of boundary conditions for LBM, in general, and with a particular focus on the ELBM for high Reynolds number flow. The key difference of our approach is not a particular way to derive suitable variables at the boundary nodes (target values of the fluid density and velocity) but rather the way these fields are introduced to the lattice Boltzmann system. Direct manipulation of the population, either the equilibrium or non-equilibrium part or both, to match the wall conditions results in unstable configuration for the (LBM and the) ELBM scheme. A physical handling

Corresponding author. Tel.: +41 044 632 66 28. E-mail address: [email protected] (I.V. Karlin).

0378-4371/$ – see front matter © 2013 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2012.12.034


S.S. Chikatamarla, I.V. Karlin / Physica A 392 (2013) 1925–1930

Fig. 1. Bounce back rule for flat wall: ‘‘Missing’’ populations of the incoming velocities −vi (dashed arrows) are completed by the populations of the anti-parallel velocities of particles leaving the fluid domain vi (bold arrows).

of the density and velocity at the wall is needed for a stable scheme. To this end, we derive a distribution which bears similarity to the well-known Tamm–Mott-Smith approximation for a shock wave in the Boltzmann equation theory [10]. This key innovation resolves the issue of boundary conditions at high Reynolds numbers for wide applications of ELBM for complex flows. 2. Entropic lattice Boltzmann equation Let us remind that in the ELBM scheme, populations associated with the discrete velocities vi evolve according to the following kinetic equation, fi (x + vi , t + 1) − fi (x, t ) = αβ(fi


− fi ).



In the above, fi is the local equilibrium, which minimizes the entropy function, H (f ) = i fi ln(fi /Wi ), where the weights Wi are lattice-specific constants. In Eq. (1), α is the maximal over-relaxation, which is operationally available as the positive root of the entropy condition:

H (f + α(f eq − f )) = H (f ).


¯ This  entropy estimate is the key, as it assures the discrete-time H-theorem: For β ∈ [0, 1], the total entropy H (t ) = ¯ ¯ x H (f (x, t )) is not increasing, H (t + 1) ≤ H (t ). Whenever the simulation is resolved (populations stay close to the local equilibrium), the maximal over-relaxation parameter α becomes fixed automatically to the value α = 2. Then the ELBM (1), (2) self-consistently becomes equivalent to the standard LBM equation and recovers Navier–Stokes equations with the   kinematic viscosity ν = cs2

1 2β

1 2

, where cs is speed of sound (a O(1) lattice-dependent constant). When the grid is

coarsened, local instabilities due to lack of resolution typically lead to a collapse of LBM. It is in this situation that ELBM proceeds with the adaptive relaxation. To this regard, ELBM is a natural extension of LBM into sub-grid simulations; with the distinctive trait that the stabilization mechanism is directly informed through the second principle (H-theorem). This is well reflected by the specific way ELBM mends instabilities; most of the time during the simulation the relaxation parameter remains constant everywhere, so that indeed ELBM collapses to LBM. It is only at the onset of a local instability, that ELBM deploys its built-in entropic stabilization capability. These stabilization events maybe rare in time and very localized in space (so are incipient instabilities) but they make the whole difference. This local self-adaptive stabilization (no fine tuning of parameters) is what makes ELBM distinct from other LB methods. These features were verified in periodic flow geometries, and the question remains as to how to extend ELBM to wall-bounded flows. 3. Boundary conditions For the sake of presentation, we shall consider the example of a flat wall and the well-known bounce back (BB) boundary condition. The situation at a boundary node (B) is explained in Fig. 1: In order to complete the propagation step, the ‘‘missing’’ ¯ are identified with the populations of antipopulations of the discrete velocities (we denote these velocities as sub-set D) ¯ fibb = fk , where vk = −vi . Note that the bounce back rule implies that at parallel velocities impinging the wall: for i ∈ D, the boundary node, the density and the velocity will be as

{ρbb , ρbb ubb } =

 ¯ i∈D

fibb {1, vi } +

fi {1, vi }.


¯ i̸∈D

We shall now reverse the perspective: rather than considering ρbb and ubb as a mere implication of the bounce back populations, we shall view ρbb and ubb as the target values of the density and velocity at the boundary node, ρtgt =

S.S. Chikatamarla, I.V. Karlin / Physica A 392 (2013) 1925–1930


ρbb , utgt = ubb , and we shall realize these target values through a different set of populations. In general, populations which eq neq complete the propagation step at the boundary node can be written fi∗ = fi (ρtgt , utgt ) + fi , where the non-equilibrium  neq  neq neq = 0, i fi vi = 0. The non-equilibrium parts of contributions fi satisfy the mass and momentum conservation: i fi the populations are proposed to be of the following form (we shall explain why shortly): fi




= fieq (ρtgt , utgt ) − fieq (ρ, u), eq

= fi − fi (ρ, u),

¯; i∈D


¯. i ̸∈ D


Here the instantaneous density and velocity, ρ and u, are distinct from the target values and are uniquely defined by the mass and momentum conservation,

{ρ, ρ u} =


fi (ρtgt , utgt ){1, vi } +

¯ i∈D

fi {1, vi }.

¯ i̸∈D eq

Finally, these pre-collision populations fi∗ are modified according to the standard ELBM collision rule, fi′ = fi∗ +αβ(fi − fi∗ ). Summarizing, with the proposed populations (4) and (5), the propagation and the collision at a boundary node returns the following populations: eq


fi′ = (1 − [αβ − 1])fi (ρtgt , utgt ) + (αβ − 1)fi (ρ, u), eq



fi = fi + αβ[fi (ρ, u) − fi ] + [fi (ρtgt , utgt ) − fi (ρ, u)], ′

¯, i∈D ¯. i ̸∈ D

(6) (7)

¯ Eq. (7) is a result of the usual propagation Let us now explain populations (6) and (7). Firstly, for the set of velocities i ̸∈ D, eq eq and collision plus a forcing in the form [fi (ρtgt , utgt )− fi (ρ, u)] (cf. Eq. (8) and Ref. [11]) which compensates for the loss (or excess) of velocity and density with respect to the target values. Secondly, the ‘‘missing’’ populations (6) are represented as a linear combination (with a coefficient dependent on the viscosity through αβ ) between the two equilibria, the one at the instantaneous density and velocity, and the other one at their target values. In the context of the classical kinetic theory, this corresponds to a bi-modal distribution (a linear combination of two Boltzmann–Maxwell distribution functions) and is the well-known Tamm–Mott-Smith approximation for the shock wave solution of the Boltzmann equation [12,10]. Appearance of a similar expression is not surprising in the present context as the main problem is a discontinuity propagation (shock) in both LBM and ELBM simulations whenever the standard realization of bounce back is used [8,9]. Note that, while in the classical case, the Tamm–Mott-Smith distribution is an approximation, in our case it is constructed to be the exact result of the propagation–collision dynamics at the boundary node. For a resolved simulation (α = 2), populations (6) and eq (7) have the following limits: At β → 1/2 (high viscosity, or equilibration limit), we have fi′ = fi (ρtgt , utgt ) for all the populations; in other words, at high viscosity, all the populations are equilibrated at the target values of the density and the velocity, as expected. In the opposite limit, as β → 1 (low viscosity, or over-relaxation), the ‘‘missing’’ populations are eq ¯ Note that in this case the also set at the equilibrium but for the instant values of density and velocity, fi′ = fi (ρ, u), i ∈ D. state of the ‘‘missing’’ populations is highly non-equilibrium with respect to the target equilibrium. In other words, instead of a standard realization of the target values through the bounce back populations, we proposed their realization with a populations based on the Tamm–Mott-Smith distribution in order to avoid the otherwise uncontrolled shock development. For a general case of the boundary node, the target velocity utgt is computed using interpolation between the velocity at the adjacent fluid nodes and the velocity at the wall, while ρtgt is computed based on the bounce back condition to ensure mass conservation. 4. Simulations 4.1. Turbulent channel flow Moving on to simulation results; turbulent flow in a rectangular channel is the classical benchmark for assessing the new boundary condition. The ELBGK model with the proposed boundary conditions was realized for the standard D3Q15 lattice. We note that the present realization of the ELBM used an optimized root solver for the implicit equation (2), based on the methods suggested in Ref. [13], which makes each update about 2–4 times slower compared to the standard LBGK. The force which is universally used to drive the flow was implemented in the present ELBGK scheme following the method of Kuperstokh (so-called exact difference method, [11]); that is, the basic kinetic equation is fi (x + vi , t + 1) − fi (x, t ) = αβ(fi


− fi ) + [f eq (ρ, u + ∆u) − f eq (ρ, u)],


where the added term on the right hand side takes into account the change of the velocity due to the presence of the acceleration g per time step ∆t = 1 : ∆u = g ∆t. The flow is characterized by Reynolds number Reτ = uτ h/ν , with uτ the wall friction velocity and h the channel halfheight. In a resolved simulation, friction velocity can be directly measured from its definition, u2τ = τw /ρ , where τw is the stress at the wall. In a typical sub-grid simulation, the resolution near the wall may become insufficient to evaluate the stress


S.S. Chikatamarla, I.V. Karlin / Physica A 392 (2013) 1925–1930

Fig. 2. (Color online) Mean stream wise velocity as a function of a distance from the wall at Reτ = 180. Line: ELBM with present boundary condition, ∆+ ≈ 2; symbol: ELBM with ∆+ ≈ 4; dashed line: filtered DNS [14].

accurately. In that case it is evaluated from the averaged flow profile at the channel center: u2τ = −(dp¯ /dx)(h/ρ), where p¯ is the time-averaged pressure. This allows us to set the target value of Reτ at the beginning of the simulation by choosing accordingly the acceleration g = Re2τ ν 2 /h3 in the force term of the ELBGK equation (8). This pre-set value of Reτ is verified at the end of the simulation by comparing it to the log law of the wall. Let us remind that the log law suggests that the average velocity at high Reynolds numbers varies with the distance from the wall in a self-similar fashion; if the distance is measured in wall units, y+ = yuτ /ν , then the reduced average velocity u+ = u/uτ may be written as, u+ = κ −1 ln y+ + C + , where κ ≈ 0.41 is the von Karman constant, and C + ≈ 5.5 (the log law is typically expected to hold for y+ ≥ 30, with logarithmic accuracy). The aforementioned verification assumes that the boundary layer extends up to the middle of the channel, so that u0 = uτ (κ −1 ln(huτ /ν) + C + ), with u0 the average velocity at the channel center, gives implicit equation for uτ , and hence another estimate of Reτ . Agreement within five percent between the three ways of estimating Reτ was considered as an indication of consistency of the simulation. The lowest Reynolds number was Reτ = 180 which is a well-documented value for the resolved direct numerical simulation (DNS) [14]. A convenient way to characterize the resolution is to measure the grid spacing ∆ in wall units, ∆+ = uτ ∆/ν . For the lattice Boltzmann method, grid spacing is uniform ∆ = 1. In the DNS which uses non-uniform grids, resolution ∆+ varies from ∆+ ≈ 0.1 at the wall to ∆+ ≈ 5 at the channel center, while in the previous LBM simulations ∆+ ≈ 1.5 was used [8,9]. Note that LBM with the standard bounce back has a stability bound of ∆+ ≈ 2.3 [8] while the ELBM becomes contaminated with propagating shocks at ∆+ ≈ 3 [9]. With the present boundary condition, we have first verified that on fine grids (∆+ ≈ 1.5–2) the results are identical to the standard bounce back method. Then a series of simulations on a much coarser grid with the resolution ∆+ ≈ 4 were performed. Comparison of mean velocity remains excellent in Fig. 2. The same coarse grid was used to simulate channel at the well-documented higher value, Reτ = 380 [14]. As the number of nodes remains same, the resolution drops drastically, ∆+ ≈ 8 for Reτ = 380. Comparison of DNS data with the present simulation is even striking for the mean velocity in Fig. 3: while even the first node is located at the end of the buffer sub-layer, in the log law domain the comparison remains excellent. This obviously demonstrates that the ELBM with the present boundary condition is the parameter-free robust method for sub-grid simulations. 4.2. Turbulent flow past a circular cylinder The proposed change of perspective on the boundary conditions, illustrated above with the target density and velocity from the bounce back prescription allows for a rather straightforward generalization to curved and complex boundaries. Indeed, all that is required is a different input for the target ρtgt and utgt which would reflect the location of the wall while the realization through the Tamm–Mott-Smith populations (6) and (7) remains the same. To this end, we used interpolation for the velocity and conservation of mass at the boundary nodes to evaluate ρtgt and utgt . Note that, while the use of interpolation is a common practice in the case of curved boundaries [2], previous attempts to use interpolation on the populations rather than on the utgt failed to produce useful boundary conditions for high Reynolds number simulations. For the illustration, we present a simulation of a flow past a circular cylinder at Reynolds number Re = 3300 based on the far-field inlet velocity U0 and the cylinder diameter D, and compare it with a recent near wake DNS simulation [15]. The aspect ratios

S.S. Chikatamarla, I.V. Karlin / Physica A 392 (2013) 1925–1930


Fig. 3. (Color online) Mean stream wise velocity as a function of a distance from the wall at Reτ = 380. Symbol: ELBM with present boundary condition, ∆+ ≈ 8. Line: filtered DNS [14].

Table 1 Comparison of the Strouhal number and the drag coefficient by ELBM with reference DNS and LES simulations [15,16] and with experimental data [17,18]. Strouhal number ELBM DNS [15] DNS [16]

Drag coefficient 0.20 0.21 0.21

ELBM [17], experiment [18], experiment

0.91 0.92 0.85

Fig. 4. (Color online) Vorticity iso-surfaces in the ELBM simulation of a flow past a circular cylinder at Re = 3300.

of the simulation domain were chosen to be 4D, 21D and 30D in the span-wise, vertical and the stream-wise directions, respectively. The cylinder was placed 5D downstream of the inlet. The choice of the domain aspect ratios follows the DNS simulation of Ref. [15]. In Fig. 4, a snapshot of the vorticity is shown. The comparison of the time-averaged velocity to the DNS data [15] is excellent in Fig. 5; in particular, the ‘‘U-shape’’ streamwise velocity profile is correctly reproduced in the near wake (x/D = 1) which is possible only for well-resolved simulations, according to Ref. [15]. This demonstrates the robustness of the suggested approach to boundary conditions, as the number of points representing the circumference of the cylinder section is C ≈ 100 in the ELBM simulation (D = 30 grid points), whereas in the DNS [15] it was C = 1206. In Table 1, we present a comparison of the computed Strouhal number with the DNS [15] and LES [16] data, as well as the comparison of the drag coefficient with experimental data [17,18]. Both the Strouhal number and the drag coefficient demonstrate the same good quality already revealed by the comparison of the velocity profiles. This indicates that the proposed boundary conditions, combined with the ELBM, are suitable for sub-grid high Reynolds number simulations.


S.S. Chikatamarla, I.V. Karlin / Physica A 392 (2013) 1925–1930

Fig. 5. (Color online) Mean stream-wise (left) and vertical (right) velocity profiles at the locations x/D = 1; 2; . . . ; 6 in the wake of the circular cylinder. Line: ELBM simulation with D = 30 grid points; Symbol: DNS data [15].

Acknowledgment This work was supported by the ERC Advanced Grant ELBM. References [1] I.V. Karlin, A. Ferrante, H.C. Ottinger, Perfect entropy functions of the lattice Boltzmann method, Europhys. Lett. 47 (1999) 182–188. [2] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, Oxford, 2001. [3] S.S. Chikatamarla, C.E. Frouzakis, I.V. Karlin, A.G. Tomboulides, K.B. Boulouchos, Lattice Boltzmann method for direct numerical simulation of turbulence, J. Fluid Mech. 656 (2010) 298–308. [4] I.V. Karlin, S. Succi, S.S. Chikatamarla, Comment on ‘‘numerics of the lattice Boltzmann method: effects of collision models on the lattice Boltzmann simulations’’, Phys. Rev. E 84 (2011) 068701. [5] Z. Guo, C. Zheng, B. Shi, An extrapolation method for boundary conditions in lattice Boltzmann method, Phys. Fluids 14 (2002) 2007. [6] Q.S. Zou, X.Y. He, On pressure and velocity boundary conditions for the lattice Boltzmann BGK model, Phys. Fluids 9 (1997) 1591. [7] J. Latt, B. Chopard, O. Malaspinas, M. Deville, A. Michler, Straight velocity boundaries in the lattice Boltzmann method, Phys. Rev. E. 77 (2008) 056703. [8] P. Lammers, K.N. Beronov, R. Volkert, G. Brenner, F. Durst, Lattice BGK direct numerical simulation of fully developed turbulence in incompressible plane channel flow, Comput. & Fluids 35 (2006) 1137–1153. [9] M. Spasov, D. Rempfer, P. Mokhasi, Simulation of a turbulent channel flow with an entropic lattice Boltzmann method, Internat. J. Numer. Methods Fluids 60 (2009) 1241–1258. [10] C. Cercignani, Theory and Applications of the Boltzmann Equation, Scottish Academic Press, 1975. [11] A.L. Kuperstokh, New method of incorporating a body force term into the lattice Boltzmann equation, in: Proc. 5th International EDH Workshop, pp. 241–246. [12] H.M. Mott-Smith, The solution of the Boltzmann equation for a shock wave, Phys. Rev. 82 (1951) 885–892. [13] S.S. Chikatamarla, S. Ansumali, I.V. Karlin, Entropic lattice Boltzmann models for hydrodynamic in three dimensions, Phys. Rev. Lett. 97 (2006) 010201. [14] R. Moser, J. Kim, N.N. Mansour, Direct numerical simulation of turbulent channel flow up to Reτ = 590, Phys. Fluids 11 (1999) 943. [15] G.J. Wissink, W. Rodi, Numerical study of the near wake of a circular cylinder, Int. J. Heat Fluid Flow 29 (2008) 1060–1070. [16] A. Kravchenko, P. Moin, Numerical studies of flow over a circular cylinder at Re 3900, Phys. Fluids 12 (2000) 403. [17] C. Wieselberger, Neuere festellungen uber die gestze des flussigkeits und luftwiderstands, Phys. Z. 22 (1921) 321–328. [18] J.D. Anderson, Fundamentals of Aerodynamics, second ed., McGraw-Hill, 1991, 228–236.