Accepted Manuscript
Lattice Boltzmann methods in porous media simulations: from laminar to turbulent flow Ehsan Fattahi, Christian Waluga, Barbara Wohlmuth, Ulrich Rude, ¨ Michael Manhart, Rainer Helmig PII: DOI: Reference:
S00457930(16)303061 10.1016/j.compfluid.2016.10.007 CAF 3295
To appear in:
Computers and Fluids
Received date: Revised date: Accepted date:
30 May 2016 9 September 2016 7 October 2016
Please cite this article as: Ehsan Fattahi, Christian Waluga, Barbara Wohlmuth, Ulrich Rude, ¨ Michael Manhart, Rainer Helmig, Lattice Boltzmann methods in porous media simulations: from laminar to turbulent flow, Computers and Fluids (2016), doi: 10.1016/j.compfluid.2016.10.007
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights • Systematic evaluation of lattice Boltzmann models for porous media
CR IP T
simulation, from creeping to turbulent flow regimes.
• Evaluation of collision models and boundary schemes in terms of parallel effciency and numerical accuracy.
• Effect of the lattice model on permeability and drag force.
AC
CE
PT
ED
M
AN US
• Parallel simulation of pore scale resolved flow scenarios.
1
ACCEPTED MANUSCRIPT
Lattice Boltzmann methods in porous media simulations: from laminar to turbulent flow
CR IP T
Ehsan Fattahia,b , Christian Walugaa , Barbara Wohlmutha,c , Ulrich Rüdeb , Michael Manhartd , Rainer Helmige
a M2  Zentrum Mathematik, Technische Universität München, Garching Department of Computer Science 1, Universität ErlangenNürnberg, Erlangen c Department of Mathematics, University of Bergen, Bergen, Norway d Fachgebiet Hydromechanik, Technische Universität München, Munich e Institut für Wasser und Umweltsystemmodellierung, Universität Stuttgart, Stuttgart
AN US
b
Abstract
The Lattice Boltzmann method has become a popular tool for determining the correlations for drag force and permeability in porous media for a wide range
M
of Reynolds numbers and solid volume fractions. In order to achieve accurate and relevant results, it is important not only to implement very efficient
ED
code, but also to choose the most appropriate simulation setup. Moreover, it is essential to accurately evaluate the boundary conditions and collision
PT
models that are effective from the Stokes regime to the inertial and turbulent flow regimes. In this paper, we compare various noslip boundary schemes
CE
and collision operators to assess their efficiency and accuracy. Instead of assuming a constant volume force driving the flow, a periodic pressure drop
AC
boundary condition is employed to mimic the pressuredriven flow through the simple sphere pack in a periodic domain. We first consider the convergence rates of various boundary conditions with different collision operators in the Stokes flow regime. Additionally, we choose different boundary conditions that are representatives of first to thirdorder schemes at curved boundaries Preprint submitted to Elsevier
October 7, 2016
ACCEPTED MANUSCRIPT
in order to evaluate their convergence rates numerically for both inertial and turbulent flow. We find that the multireflection boundary condition yields second order for inertial flow while it converges with third order in
CR IP T
the Stokes regime. Taking into account the both computational cost and
accuracy requirements, we choose the central linear interpolation bounce
back scheme in combination with the tworelaxationtime collision model. This combination is characterized by providing viscosity independent results
AN US
and secondorder spatial convergence. This method is applied to perform simulations of touching spheres arranged in a simple cubic array. Full and reducedstencil lattice models, i.e., the D3 Q27 and D3 Q19 , respectively, are
compared and the drag force and friction factor results are presented for Reynolds numbers in the range of 0.001 to 2, 477. The drag forces computed
M
using these two different lattice models have a relative difference below 3% for the highest Reynolds number considered in this study. Using the evaluation
ED
results, we demonstrate the flexibility of the models and software in two large scale computations, first a flow through an unstructured packing of spherical
PT
particles, and second for the turbulent flow over a permeable bed. Keywords: Lattice Boltzmann method, pore scale simulation, turbulence, Darcy flow, noslip boundary condition, collision operator, periodic pressure
CE
boundary
AC
1. Introduction The simulation of fluid flow and transport processes through porous media
is an important field that contributes to a number of scientific and engineering applications, e.g., oil extraction, groundwater pollution, nuclear waste storage, 3
ACCEPTED MANUSCRIPT
and chemical reactors. The application of porescale simulations is challenging in most practical situations because the system under study is often larger by orders of magnitude than the characteristic size of the pores. Thus, for
CR IP T
practical purposes, many computational techniques are based on macroscopic models that average over many pores and consider average flow rates.
While the rigorous, analytic upscaling of the porescale problem at lower Reynolds numbers has received much attention in the literature [1, 2], similar
AN US
approaches for higher Reynolds numbers have not been demonstrated yet primarily because of the immense mathematical difficulties that arise in
flow models if a moderate Reynolds number cannot be assumed. However, with the availability of evermore powerful supercomputers, direct numerical simulation (DNS) has been established as a third possibility for the analysis
M
of homogenized models that can complement the classical experimental and rigorous mathematical averaging approaches.
ED
In the past two decades, the lattice Boltzmann method (LBM) has attracted the interest of researchers in CFDrelated fields [3, 4]. In contrast
PT
to traditional CFD approaches based on the conservation of macroscopic quantities such as mass, momentum, and energy, the LBM models a fluid
CE
using the kinetics of discrete particles that propagate (streaming step) and collide (relaxation step) on a discrete lattice mesh. Owing to this kinetic nature, microscopic interactions in fluid flow can be handled even in com
AC
plex geometries such as those in microfluidic devices or porous media [5–7]. Moreover, the inherently local dynamics used in LBM afford efficient implementation and parallelization of both of the fundamental algorithmic stages.
This allows to harness the computational power of currently available and 4
ACCEPTED MANUSCRIPT
emerging supercomputing architectures [8–11]. In this study, we use the waLBerla framework (widely applicable LatticeBoltzmann from Erlangen) [10], which is specifically designed to be used for
CR IP T
massively parallel fluid flow simulations; this enables us to compute problems with resolutions of more than one trillion (1012 ) cells and up to 1.93 trillion
cell updates per second using 1.8 million threads [12]. waLBerla has already been used to study the flow through moderately dense fluidparticle systems
AN US
[13] and to simulate largescale particulate flows [14]. It is based on four main software concepts: blocks, sweeps, communication, and boundary handling [15]. The modular structure of this framework allows for the implementation of new LBM schemes and models.
However, having immense computational power at hand is not enough to
M
solve relevant problems. For a threedimensional LBM simulation, stencils that differ with respect to the velocity directions can be used. Lattice models
ED
generally require an exact evaluation of their velocity moments up to secondorder to consistently recover NavierStokes dynamics in the continuum limit.
PT
However, they may behave differently at a discrete level, which in turn can lead to violation of some important physical requirements. Investigating
CE
an intermediate Reynolds number flow (50 < Re < 500) through a nozzle, White and Chong [16] demonstrated that the D3 Q19 lattice model, a three dimensional model with a 19 directional stencil can cause a lack of isotropy
AC
in the flow field. This deficiency, which was found to be independent of the collision model, the grid resolution, and the Mach number, can be removed by using the D3 Q27 lattice model which has a 27directional stencil. It was found that lattice models with a plane having less than six velocity vectors 5
ACCEPTED MANUSCRIPT
are not fully isotropic and can produce qualitatively different results. Similar observations have been reported by other researchers [4, 17, 18]. Recently, Silva and Semiao [19] theoretically analyzed the truncation errors of different
CR IP T
lattice models and showed that reduced lattice schemes, i.e., D3 Q15 and
D3 Q19 , produce spurious angular velocities. They showed that, for convectiondominated flows the rotational invariance will be violated if a reduced lattice scheme is used. Nash et al. [20], on the other hand, investigated the Dean
AN US
flow at Re ≈ 400 with strong secondary flow and demonstrated that none
of the respective lattice models produces large artifacts in the velocity field. Their results suggest that using a greater number of velocity sets for a lattice model, produces better results for interpolating boundary schemes but worse results when using a simple bounceback boundary scheme.
M
Moreover, the explicitness of classical LBM means that the spatial and temporal discretization characteristics are strongly coupled. Hence, special
ED
care must be taken when performing porescale simulations to properly incorporate the physics at the boundaries and inside the domain without
PT
overresolving the problem. Several methods have been proposed for the implementation of the LBM on nonuniform grids to improve the geometrical
CE
flexibility [21–23], as well as the interpolating boundary conditions that can be used by the classical LBM [24–26]. As was already pointed out in the evaluation of Pan et al. [27], this requires a suitable combination of collision
AC
and boundary operators. The evaluation of different boundaries is mostly done for Stokes flow regimes [28, 29]; however, combinations of boundary conditions and collision models must also be evaluated for high Reynolds number flow. 6
ACCEPTED MANUSCRIPT
A common way to simulate pressuredriven flow in the LBM is to replace the pressure gradient with an equivalent body force and apply streamwise periodic boundary conditions. However, previous studies [30–33] have shown
CR IP T
that using this approach does not lead to correct flow fields for flow through
complex geometries, such as e.g., porous media. In this study, we drive flow by imposing a pressure gradient while applying a periodic boundary condition in the streamwise direction in order to allow the flow to develop based on the
AN US
geometry. We employ a simple periodic pressure scheme for an incompressible
flow and investigate it in combination with different types of collision operators and a number of first, second, and thirdorder bounceback schemes. We simulate flow through simple sphere packs using different LBM approaches, which we will briefly outline in Section 2. Their accuracy and
M
convergence rates for flow in the Stokes regime, are investigated in Section 3. In addition, the computational costs of different boundary schemes are as
ED
sessed in order to choose the best combination for highly resolved simulations in high Reynolds numbers flow. After finding a suitable configuration, we
PT
then examine the spatial convergence of the boundary schemes in the laminar steady and fluctuating flow regime. Using the results of this evaluation, in
CE
Section 3.2 we simulate the flow through a simple sphere pack for two different solid volume fractions and investigate the effect of lattice velocity sets on the drag force and permeability in high Reynolds number flow by comparing
AC
the full and the reduced stencil lattice models. By sampling over the regime Rep ∈ [0.001, 2, 477] for a regular packing of touching spheres, we numerically investigate the friction factor based on the Reynolds number. To highlight
the ability of the proposed LBM to deal with complex geometries Section 4 7
ACCEPTED MANUSCRIPT
reports on simulations with a more realistic sphere packing geometry that is generated using a rigid body dynamics simulator. In the first example we study the flow through such a random sphere pack. Finally, the flow in a models the turbulent flow over a permeable wall. 2. LBM approaches for porous media flow
CR IP T
model porous region is combined with a freeflow so that the coupled regime
AN US
The lattice Boltzmann equation (LBE) is a simplification of the Boltzmann
kinetic equation where we assume that particle velocities are restricted to a discrete set of values ek , k = 0, 1, . . . , q, i.e., that the particles can only move along a finite number of directions, connecting the nodes of a regular lattice [34, 35]. In the following, we consider a three dimensional setting with
M
q = 18, 26, the socalled D3 Q19 and D3 Q27 model, respectively. Discretizing in time using a timestep size of ∆t = tn+1 − tn , the semidiscrete LBE then
ED
becomes
k = 0, . . . , q,
(1)
PT
∆t−1 (fk (x + ek ∆t, tn+1 ) − fk (x, tn )) = gk (x, tn ),
where fk (x, tn ) represents the probability of finding a particle at position
CE
x and time tn with velocity ek . The lefthand side of (1) corresponds to a discrete representation of the Boltzmann streaming operator, while the
AC
righthand side g := [gk ]qk=0 is responsible for controlling the relaxation to a local equilibrium. It is generally split into g := ∆t−1 Ω(x, tn ) + F(x, tn ),
where Ω(x, tn ) is a collision term function of f := [fi ]qi=0 , and F is a forcing term that drives the flow. Provided that the modeled fluid is close to its equilibrium state, Bhatnagar, Gross, and Krook (BGK) [36] proposed that 8
ACCEPTED MANUSCRIPT
the discrete local equilibria in the collision processes can be modeled using a local Maxwell distribution. The collision operator takes the general form (2)
CR IP T
Ω(x, tn ) = −R(f (x, tn ) − f eq (x, tn )),
where R is a relaxation operator and f eq (x, tn ) is an equilibrium distribution function of f (x, tn ) that, for incompressible flow, is given by [37]
2 1 −2 1 −4 , fkeq (x, tn ) = wk ρ + ρ0 c−2 s ek · u + 2 cs (ek · u) − 2 cs u · u
(3)
AN US
where wk is a set of weights normalized to unity, ρ = ρ0 + δρ ( where δρ is the
density fluctuation, and ρ0 is the mean density which we set to ρ0 = 1), and √ cs = ∆x/( 3∆t) is the lattice speed of sound ( where ∆x is the lattice cell width). The macroscopic values of density ρ and velocity u can be calculated i.e.,
Xq
k=0
fk ,
u = ρ−1 0
ED
ρ=
M
from f as zeroth and first order moments with respect to the particle velocity, Xq
k=0
ek fk .
(4)
In a lattice Boltzmann scheme, the computation is typically split into a
PT
collision and streaming step, that are given as
CE
f˜k (x, tn ) − fk (x, tn ) = ∆t gk (x, tn ), fk (x + ek ∆t, tn+1 ) = f˜k (x, tn ),
(collision) (streaming)
respectively, for k = 0, . . . , q. The execution order of these two steps is
AC
arbitrary and may vary from code to code for implementation reasons. For instance, in waLBerla, the order is first streaming and then collision. This has the benefits that the stream and collision steps can be fused and that macroscopic values do not need to be stored and then later retrieved from memory. 9
ACCEPTED MANUSCRIPT
When studying the numerical properties (such as stability, convergence, and accuracy) and the computational costefficiency of LBMbased approaches, two important components must be investigated. First, the treatment of the
CR IP T
collision term Ω(x, t) in the relaxation step must be considered. Second, as the nodes of the lattice in the LBM approach are distinguished only into fluid and solid nodes, the choice of the bounceback operator for the fluidwall interaction is important for obtaining numerical accuracy.
AN US
Unfortunately, relaxation and bounceback cannot be treated as two independent components as their interplay is essential for finding an optimal
tradeoff between computational effort and physical accuracy. For instance, in the context of porous media flow, use of the BGK collision operator in an approach with a singlerelaxationtime (SRT) ω, i.e., RSRT := ωI has
M
been reported to suffer from numerical instabilities and to exhibit unwanted boundary effects [27, 38]. More precisely, the effective poresize modeled by
ED
the numerical scheme becomes viscositydependent. On a macroscopic level this has the effect of making the calculated permeability in the Stokes flow
PT
viscositydependent, even though it physically should be independent of fluid properties [27]. The limitations of the simple and computationally appealing
CE
approach outlined above motivate us to investigate more sophisticated schemes. Our goal, which is similar to that of Pan et al. [27], is to investigate different LBM approaches with the application of porous media flow in mind. However,
AC
instead of lumping all boundary conditions into a bodyforce, we place our focus on boundary conditions that are suitable for high Reynolds number flow driven by pressuregradients. Before we present a detailed evaluation of suitable combinations for pore10
ACCEPTED MANUSCRIPT
scale simulations, we will briefly summarize the various common approaches. 2.1. Collision operators
CR IP T
Pan et al. [27] reported that, by using multirelaxationtime [39, 40] (MRT) schemes, it is possible to overcome the deficiencies of the SRT approach when
applied to porous media flow. In the MRT model, the relaxationtimes
for different kinetic modes can be treated separately, which allows for finetuning of the free relaxationtimes in order to improve numerical stability
AN US
and accuracy. This is done using the following ansatz: RMRT := M−1 S M,
(MRT)
where the components of the collision matrix in the MRT are developed to
M
reflect the underlying physics of collision as a relaxation process. The rows of the orthogonal matrix M are obtained using the GramSchmidt orthogo
ED
nalization procedure applied to polynomials of the Cartesian components of the particle velocities ek , and S := diag[s0 , s1 , ..., sq ] is the diagonal matrix holding the relaxation rates sk : cf., e.g., d’Humières [40] for details. Clearly,
PT
within this class of relaxation schemes, SRT approaches are contained as a special case. In the D3 Q19 lattice model, the parameters s0 , s3 , s5 and s7
CE
are the relaxation parameters corresponding to the collision invariant ρ, and j := ρu, respectively, which are conserved quantities during a collision. These
AC
are set to zero in the absence of a forcing term, i.e., if Fk = 0. Pan et al. [27] proposed to choose the remaining modes as follows: (i) setting the viscous stress vectors s9 , s11 , s13 , s14 , and s15 equal to sν , which is related to the kinematic viscosity ν = µ/ρ as sν = (3ν + 0.5)−1 , and (ii), setting the kinetic modes s1 , s2 , s4 , s6 , s8 , s10 , s12 , s16 , s17 , and s18 to sζ = 8(2 − sν )(8 − sν )−1 , 11
ACCEPTED MANUSCRIPT
which is also related to the kinematic viscosity. However, as their evaluation showed, this does not lead to viscosityindependent results for the macroscopic quantities, except for the multireflection boundary scheme. As Khirevich
CR IP T
et al. [29] recently stated, this problem can be solved if one modifies the
aforementioned model in a way such that the symmetric energy modes, i.e.,
s1 and s2 , maintain a constant ratio (1/sν − 0.5)/(1/si − 0.5) while sν varies. Based on a symmetry argument, Ginzburg [38] proposed a model based on
AN US
tworelaxationtimes (TRT), which can be seen as the minimal configuration
that provides just enough free relaxation parameters to avoid nonlinear dependencies of the truncation errors on the viscosity in the context of porous media simulations [28, 41, 42]. The scheme was derived from the MRT approach by splitting the probability density functions fk into symmetric and
M
antisymmetric components fk+ := 12 (fk + fk¯ ) and fk− := 12 (fk − fk¯ ), where k¯
is the diametrically opposite direction to k [28, 41, 42]. Performing a separate
ED
relaxation using the two corresponding relaxation rates ω + and ω − yields the operator
(TRT)
PT
RTRT := ω + R+ + ω − R− ,
where R+ and R− are the tensorial representations of the operators extracting
CE
the symmetric and antisymmetric components, respectively. The eigenvalue of the symmetric components is again related to the kinematic viscosity as
AC
ω + = (3ν + 0.5)−1 , and the second eigenvalue ω − ∈ (0, 2) is a free parameter. For steady nonlinear flow situations, it has been demonstrated that most of the macroscopic errors depend on the socalled “magic” parameter 1 1 1 1 Λ= − − ω+ 2 ω− 2 12
ACCEPTED MANUSCRIPT
which must be determined for the specific flow setup. The choice Λ = given as an optimal value for stability. Another choice, namely Λ =
3 , 16
1 4
is
yields
the exact location of bounceback walls in case of Poiseuille flow in a straight
CR IP T
channel [28, 29]. Previous studies [26, 29] show that the accuracy depends on the magic number. In our studies, we choose different values in the range Λ ∈ (0, 34 ]. The exact numbers will be given in the corresponding section. 2.2. Boundary conditions
AN US
In porescale simulations, two types of boundary conditions are typically encountered: a solidwall interaction of fluid particles that come into contact
with the porous matrix (corresponding to a noslip condition in continuummechanics terminology); and a periodic pressure forcing that is applied to
M
drive the flow via a pressure gradient in order to compute fluid or matrix properties such as the Forchheimer coefficient or the (apparent) permeability. this study.
ED
In the following, we will outline the different schemes under consideration for
PT
2.2.1. Solidwall interaction
In a simple bounceback (SBB) scheme, the wall location is represented
CE
through a zeroth order interpolation (staircase approximation), and the collision of particles with a wall is incorporated by mimicking the bounceback
AC
phenomenon of a particle reflecting its momentum upon collision with a wall, as is expected to happen halfway between a solid and fluid node. Hence, the unknown distribution function is calculated as fk¯ (xf1 , tn+1 ) = f˜k (xf1 , tn ),
13
(5)
ACCEPTED MANUSCRIPT
where, as we recall, k¯ is the diametrically opposite direction to k, and we take the values fek after collision but before streaming on the right hand
side. However, in complex geometries, the wall position is not always located
CR IP T
halfway on a regular lattice. Hence, especially at coarse resolutions, this
boundary treatment introduces severe geometric errors, which in turn lead to
boundary layer effects that can be particularly problematic in lowporosity media in which a large part of the domain is occupied by solid nodes. In
AN US
this case, implementing a discretization that adequately resolves the solid boundaries is often computationally prohibitive as the meshes would become too large; consequently, several advanced schemes have been proposed to more accurately represent non mesh aligned boundaries using spatial interpolations [24–26, 43–45].
M
In our numerical study, we compare five different interpolating boundary condition schemes to model curved boundaries: the linear interpolation
ED
bounceback (LIBB) [24] scheme; the quadratic extension (QIBB) [44] scheme; the interpolation/extrapolation bounceback (IEBB) [25] scheme; the multi
PT
reflection (MR) [28] scheme; and the central linear interpolation (CLI) [28] scheme. For the sake of completeness, we briefly recapitulate these approaches.
CE
Let xfi denote a lattice node located at a distance of at most i > 0 cells from the boundary and let q = xf1 − xw /xf1 − xb  define a normalized wall distance; cf. Fig. 1. Bounceback schemes based on higherorder interpolations
AC
require more than one fluid node between nearby solid surfaces. For instance, the LIBB condition proposed by Bouzidi et al. [24] is given by (1 − 2q)f˜k (xf2 , tn ) + 2q f˜k (xf1 , tn ), q < 1/2, fk¯ (xf1 , tn+1 ) = 1 − 1 f˜¯ (xf , tn ) + 1 f˜k (xf , tn ), q ≥ 1/2. k 1 1 2q 2q 14
(6)
ACCEPTED MANUSCRIPT
𝑥𝑓3
𝑥𝑓3
𝑒𝑘
𝑥𝑓2
𝑒𝑘
𝑥𝑓1 𝑥𝑓1
𝑞 < 0.5 𝑞 > 0.5 𝑞 = 0.5
𝑥𝑓3
𝑥𝑓2
𝑥𝑓1
𝑥𝑤
𝑥𝑤
𝑥𝑤 𝑥𝑏
CR IP T
𝑥𝑓2
AN US
Figure 1: Example of different distances of the wall to the first fluid node. Simple bounceback occurs when a lower or higher value of q = 1/2, q emerges the need of interpolated bounceback schemes. For simplicity, we only display a twodimensional illustration.
It is also possible to use quadratic interpolation to obtain the QIBB scheme [44], where the unknown distribution function is calculated as
M
q(1 + 2q)f˜k (xf1 , tn ) + (1 − 4q 2 )f˜k (xf2 , tn ) − q(1 − 2q)f˜k (xf3 , tn ), q < 1/2, fk¯ (xf1 , tn+1 ) = 1−2q ˜ 1 2q−1 f˜¯ (x , t ) + ˜ q ≥ 1/2. ¯ (xf2 , tn ), k f1 n q q(2q+1) fk (xf1 , tn ) + 2q+1 fk
(7)
ED
While extrapolation can possibly give higher order consistency, it is sensitive to numerical instabilities. As an improvement to the scheme of Filippova and
PT
Hänel [46], an interpolation/extrapolation bounceback (IEBB) model was
CE
proposed by Mei et al. [25] in which a velocity at the wall is computed via ux , q < 1/2, f ubf = 2 1 − 3 ux , q ≥ 1/2 f1 2q
AC
to determine an auxiliary distribution function n h fk∗ (xb , tn ) = wk ρ + ρ0 c32 ek · ubf + 2c94 (ek · uxf1 )2 −
3 u 2c2 xf1
that is then used to linearly interpolate the unknown as fk¯ (xf1 , tn+1 ) = (1 − X)f˜k (xf1 , tn ) + Xfk∗ (xb , tn ). 15
· ux f 1
io
(8)
(9)
ACCEPTED MANUSCRIPT
Above we let X = (2q − 1)/(1/ω − 2) for q < 12 , while for q ≥
1 2
we let
X = (2q − 1)/(1/ω + 1/2), where ω = ω + in case of the TRT scheme and ω = sν for the MRT scheme.
CR IP T
An alternative approach to obtaining highly accurate bounceback conditions is the multireflection (MR) scheme [26, 28], which reads as fk¯ (xf1 , tn+1 ) =
1−2q−2q 2 ˜ fk (xf2 , tn ) (1+q)2
1−2q−2q 2 ˜ fk¯ (xf1 , tn ) (1+q)2
−
q2 f˜¯ (xf2 , tn ) (1+q)2 k
+ f˜k (xf1 , tn ) + Fkpc ,
AN US
−
q2 f˜ (x , t ) (1+q)2 k f3 n
+
(10)
where Fkpc is the postcollision correction term. The post correction term is related to the thirdorder moments and for the MRT we use the one (2)
4Λk presented in [26] for the MR1 model. This reads as Fkpc = − 3ν(1+q) 2 fk
M
with Λk = (1/sν − 0.5)/(1/sk − 0.5), and fk2 is obtained from the thirdorder nonequilibrium and can be calculated as
ED
f (2) = M−1 S[t − teq ]
(11)
PT
where t are the thirdorder moments. The correction term for the TRT collision is employed as it is introduced
CE
in [28] as
Fkpc =
1 1 4 ω − ( − − )(f − − f eq− ). 2 (1 + q) ω 2
(12)
AC
It is worth to mention that the nonequilibrium term should be calculated before the collision. Note that this scheme must access five distribution values at three fluid
nodes to effect the update. A computationally cheaper variant is given by the central linear interpolation scheme (CLI), which only requires three values at 16
ACCEPTED MANUSCRIPT
two fluid nodes, i.e., 1−2q ˜ f (x , t ) 1+2q k f2 n
−
1−2q ˜ f¯ (xf1 , tn ) 1+2q k
+ f˜k (xf1 , tn ).
(13)
CR IP T
fk¯ (xf1 , tn+1 ) =
It should be noted that neither of the two latter schemes involves a distinction
of cases for different values of q, which allows for an efficient implementation. We note that the CLI scheme can also be modified with a post collision
correction fkpc in order to eliminate the secondorder error for steady flows.
AN US
This additional correction is not considered in the present study, but we refer to [27–29] for a further discussion and results. 2.2.2. Periodic pressure boundary condition
In many applications, fluid flow is driven by a pressure differential. For can be written as
M
an incompressible flow, the corresponding periodicity boundary conditions
ED
u(x + L i, t) = u(x, t),
p(x + L i, t) = p(x, t) + ∆p,
(14)
where L is the length of the domain in the periodic direction i, and ∆p/L is domain.
PT
the pressure gradient applied between the inlet and outlet boundaries of the
CE
In LBMbased approaches, applying this type of boundary condition is
not straightforward. Simply adjusting the pressures at the inlet and outlet
AC
boundaries produces nonphysical mass defects at the periodic boundary, as
was reported in [47, 48]. The most commonly employed approaches replace the pressure gradient by incorporating an equivalent body force; see e.g. [49–53]; however, these approaches suffer from the inability to accurately predict the pressure gradient locally for flow situations in general geometries [30–33]. For 17
ACCEPTED MANUSCRIPT
porous media applications, in which there are complex pore geometries it is therefore desirable to apply boundary conditions that are not lumped into a volume forcing and hence does not rely on rough predictions of the pressure
CR IP T
field. In this work, we employ a pressure boundary condition that can be applied for incompressible periodic flows. Here we specify the equilibrium distribution function fieq and the nonequilibrium distribution fineq = fi − fieq
separately, as we shall describe below. As the extension to multiple periodic
AN US
boundaries is straightforward, we consider in our description an essentially onedimensional setting in which flow propagates from the left (L) to the right (R) boundary.
Since the pressure is related to the density via the equation of state p = ρc2s , we consider a density difference instead of a pressure difference in
M
the following. The density at the left boundary is obtained by (15)
ED
ρL = ρR + ∆ρ.
The relaxation dynamics of the nonequilibrium distribution, in presence of
PT
the periodic boundaries can be approximated as neq neq fi,L = fi,R .
(16)
CE
By using the above formulations, the unknown distribution functions can be
AC
computed as
neq eq fi,L = fi,R + fi,R (ρL , uR ),
(17)
neq eq fi,R = fi,L + fi,L (ρR , uL ).
(18)
As the momentum in a periodic channel does not change, the implementation of this approach is simple. For instance, by considering (3), updates (17) and 18
ACCEPTED MANUSCRIPT
(18) can be performed as (19)
fi,R = fi,L − wk ∆ρ.
(20)
3. Results and discussion
CR IP T
fi,L = fi,R + wk ∆ρ,
For steadystate flow at a low Reynolds number, a generally accepted and
experimentally confirmed macroscopic relation between the pressure gradient
AN US
and flow rate is given by Darcy’s law:
−1 ∇P = −µKD U,
(21)
where µ is the fluid dynamic viscosity, KD is the permeability, and U and P
M
are the volume averaged velocity and pressure, respectively. This model, which was first proposed by Darcy [54] based on experimental investigation, is valid
ED
in the regime Rep 1, where Rep := ρdp U/µ is the Reynolds number based
on a characteristic particle diameter dp (in this study, the sphere diameter),
PT
ρ denotes the density of the fluid, and U := U · i is the scalar velocity in
streamwise direction i. In the following sections, we use the volume averaged
CE
velocity, i.e., the Darcy velocity, to calculate the Reynolds number of the flow. 3.1. Porescale simulation at Re 1
AC
To verify the implementation and assess the quality of different schemes
before tackling more complex flow problems, in the following we will first evaluate different combinations of the LBM strategies discussed in Section
2 in the Darcy regime. As the flow is linear and irrotational in the Darcy (Stokes) regime, in this section, we conduct our simulations using the D3 Q19 19
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 2: 3D view of simple sphere array with equal radii.
lattice model. We consider flow through a periodic array of spheres arranged on a regular lattice as depicted in Fig. 2. To quantify the errors, we compute the dimensionless drag coefficient
FD , 6πµU r
M
CD =
(22)
ED
where r is the sphere radius, FD is the drag force acting on the sphere and U is the volume averaged velocity (Darcy velocity). To reduce computational cost, we compute the wall distances required for the higherorder boundary
PT
conditions only once before the timestepping loop begins and reuse it in each boundary handling step. This is possible because we are assuming a
CE
nondeformable porous medium. Moreover, we exploit the periodicity in the domain and consider only one cubic representative elementary volume (REV)
AC
containing a single sphere. For details on boundary handling, we refer to the previous Section 2.2. In the following sections, we shall summarize our results related to the spa
tial discretization effect, as well as the undesired effect of viscosity dependence of the schemes. 20
ACCEPTED MANUSCRIPT
3.1.1. Convergence analysis In this test, we change the size of the sphere and increase the computational
CR IP T
domain accordingly, such that the relative solid volume fraction is fixed to χ = D/L = 0.6. Although Sangani and Acrivos [55] prepared a semianalytical
solution for a simple cubic sphere pack based on the solid volume fraction, the correlated limits are not sufficient to evaluate the convergence rate of
the boundary conditions. Therefore, first we find the limit value of each
AN US
boundary condition by leastsquares curve fitting. Subsequently, we calculate
the convergence rate of the corresponding boundary scheme. We compute the dimensionless drag force CD and consider the relative error CD /CD,∞ − 1. Our results are plotted in Figs. 3 and 4.
It can be confirmed that the higher order boundary conditions generally
M
converge with higher order than the SBB boundary scheme which is expected to be asymptotically of first order. As Fig. 4(a) indicates, all of the interpola
ED
tion boundary schemes have almost secondorder accuracy when combined with the SRT model. Moreover, for this example, we observe that the MR
PT
scheme converges with third order. In Figs. 4(b) and 4(c), we list the results for the TRT and MRT models,
CE
respectively. Here we observe a difference between the CLI and the MR schemes, i.e., the CLI scheme converges with secondorder while the MR scheme converges with thirdorder. The results also exhibit that the QIBB
AC
and IEBB schemes converge with secondorder. Here, it is worth mentioning that the simulation results using the TRT and MRT collision can be altered by changing the magic number. In this section, the magic number is set to
3/16 in which the interpolation boundary schemes result in better accuracy 21
ACCEPTED MANUSCRIPT
[29]. The results of the convergence rate are summarized in Table 1. 5
5
SBB 4.02 LIBB QIBB 4 IEBB MR 3.98 CLI
5 4.02
75 100
25
50
75 100
CD
CD
50
4
25
50 75100
Radius
3.98 25
4
3.5
4
4.5
3.98
50
75 100
CR IP T
4.5
25
4.02
4
CD
4.5
In
4
3.5
25
50 75100
Radius
(a) SRT
3.5
25
50 75100
Radius
(b) TRT
(c) MRT
AN US
Figure 3: Dependence of the dimensionless drag CD on the resolution of the sphere for simple sphere pack with relative volume fraction of χ = 0.6. Three collision models are presented: (a) SRT, (b) TRT, and (c) MRT. Each figure depicts 6 boundary conditions scheme results. All simulations were conducted with Λ = 3/16 with the D3 Q19 lattice model.
N
2
102 103 10
4
10
5
10
6
N 25
3
50 75100
PT
Radius
10 10
0
1
N1
100
N2
102 103 10
4
10
5
10
6
N 25
Radius
(b) TRT
3
10
1
N1
N2
102 103 104
N3
105 50 75100
10
6
25
Radius
50 75100
(c) MRT
CE
(a) SRT
101
CD/CD, 1
N1
101
M
10
1
SBB LIBB QIBB IEBB MR CLI
ED
CD/CD, 1
10
0
CD/CD, 1
101
Figure 4: Logarithmic relative error of the dimensionless drag CD /CD,∞ − 1 as a function
of the sphere radius (in lattice units). Three collision models are presented, (a) SRT, (b)
AC
TRT, and (c) MRT. Each figure depicts 6 boundary conditions scheme results. The solid lines are eye showing convergence rate. All simulations were conducted with Λ = 3/16
with the D3 Q19 lattice model.
these experiments we also observe preasymptotic discretization error modes where the errors for most models do not smoothly decrease with respect to the 22
ACCEPTED MANUSCRIPT
Table 1: Limit value of the dimensionless drag force and the convergence rate α calculated from the last three results of each boundary condition in Fig. 3.
TRT
QIBB
IEBB
MR
CLI
α
1.29
1.74
2.05
2.60
2.95
1.87
CD,∞
3.98340812
3.976854186
3.97337407
3.97329727
3.97282669
3.97363293
α
1.14
1.26
2.41
2.06
2.94
2.05
CD,∞
3.98180682
3.97455861
3.97369797
3.97372456
3.97373716
3.97383257
α
1.15
1.26
2.40
2.54
2.95
1.98
CD,∞
3.98167626
3.97460279
3.97319284
3.97402321
3.97378283
3.97362
AN US
MRT
LIBB
CR IP T
SRT
SBB
resolution (r < 15). Comparing for instance the LIBB model with the QIBB model reveals that the higher order scheme is less affected by this discrepancy. To investigate the source of these fluctuations, we conduct a second series of tests where the center of the sphere is shifted streamwise inside the cell. In
M
all cases, the radius of the sphere is kept fixed at 4.5 lattice units and the relative volume fraction was set to χ = 0.6. In Fig. 5, we show the results for
ED
a displacement of 0.0 − 0.5 times a cell width. Although the results obtained
by the SRTSBB are relatively good for this setup, we observe that the SBB
PT
scheme for all considered collision schemes is more sensitive with respect to the positioning of the center. As this behavior stems only from the geometry
CE
representation, the use of interpolation models significantly increases the reliability of the results. Since the displacement along the spanwise showed
AC
a similar behavior, we do not list these results. 3.1.2. Effect of the viscosity on the computed permeability As many researchers have previously reported, the choice of boundary
conditions and collision models may have a strong impact on the simulation accuracy [13, 27, 38, 56]. To have a complete understanding of the effects 23
ACCEPTED MANUSCRIPT
0
0.05
0.1
0.05
0.05
0
0.05
0.1
0
0.05
0.1
CR IP T
CD/CDref1
0.05
0.1
CD/CDref1
SBB LIBB QIBB IEBB MR CLI
CD/CDref1
0.1
0.1
0.15 0 0.1 0.2 0.3 0.4 0.5 Displacement of the sphere center
0.15 0 0.1 0.2 0.3 0.4 0.5 Displacement of the sphere center
0.15 0 0.1 0.2 0.3 0.4 0.5 Displacement of the sphere center
(a) SRT
(b) TRT
(c) MRT
Figure 5: Relative error of the drag force (CD /CD,ref − 1) as a function of the sphere center
AN US
displacement (in a lattice cell). Three collision models are presented, (a) SRT, (b) TRT, (c) MRT. Each figure depicts the results of 6 boundary conditions. All simulations were conducted by the D3 Q19 lattice model, Λ = 3/16, r = 4.5, with a solid volume fraction of χ = 0.6.
of the boundary conditions and collision schemes on the results, we conduct
M
a viscosity dependency study for a setting where the flow is driven by a periodicpressuredrop boundary condition. This is especially important when
ED
the simulation across a wide range of Reynolds numbers is of interest, since the viscosity is one of the key parameters to achieve different Reynolds numbers
PT
in the LBM simulation. To exclude preasymptotic discretization error modes from our observations and only concentrate on the physical inconsistencies
CE
introduced by a specific combination of collision and boundary treatment, we use a sufficient resolution of 16.5 (in lattice units) for the sphere radius.
AC
The permeability in the Darcy regime is only a geometrical characteristic; therefore, it is used to show the effect of viscosity. As proposed by Adler [57], the average pressure gradient in the flow direction can be computed as ∇P · i = −FD L−3 , where L is the length of the domain in flow direction. By combining this equation with Eqs. (22) and (21), we find the reference 24
ACCEPTED MANUSCRIPT
permeability Kref . Here, we use the drag force provided by Sangani and Acrivos [55] to calculate the reference permeability. The plots in Fig. 6 display the ratio of the computed permeability to the
CR IP T
reference permeability, K/Kref , for a simple sphere pack with a relative solid
volume fraction of χ = 0.6. We consider viscosity in the range [0.029, 0.45] (LB unit) and compare different collision models and boundary conditions. In
Fig. 6(a) we show that for all boundary schemes, the SRT collision results in
AN US
a permeability that is strongly dependent on the viscosity. While the IEBB
scheme is less sensitive than the other schemes, the errors we can expect from an SRT collision model are still unacceptable for our purposes. The results plotted in Figs. 6(b) and 6(c) depict the same set of experiments conducted with the TRT and MRT models, respectively. The result achieved
M
by the TRT collision operator shows that the SBB, CLI, and MR schemes all produce a viscosityindependent permeability. The LIBB scheme results
ED
in a clear viscositydependency of the permeability. Also QIBB and IEBB result in a slightly viscositydependent permeability though the effect is less
PT
pronounced than for LIBB. The results for the MRT model show the same behavior as the TRT, however, we observe a small deviation in the MR
CE
scheme, in which the relative difference is below 0.02%. The CLI and SBB also reproduce viscosity independent results. In case of the SBB we note that though they are viscosityindependent, the results severely underpredict
AC
the permeability due to the staircase representation of the geometry. We point out that our results are in accordance with those found with the TRT collision operator in Pan et al. [27]. We further note that in the comparison between the two linear boundary 25
ACCEPTED MANUSCRIPT
relaxation time 1 1.2 1.4 1.6 1.8
relaxation time 1 1.2 1.4 1.6 1.8
0.6 0.8
1.04
1.04
1.02
1.02
1.02
0.98
0.98
1
SBB LIBB QIBB IEBB MR CLI
0.98 0.96 0.1
0.2 0.3 Viscosity
1
1
0.96
0.4
0.96 0.1
(a) SRT
relaxation time 1 1.2 1.4 1.6 1.8
CR IP T
K/Kref
1.04
K/Kref
0.6 0.8
K/Kref
0.6 0.8
0.2 0.3 Viscosity
(b) TRT
0.4
0.1
0.2 0.3 Viscosity
0.4
(c) MRT
Figure 6: Viscosity dependence of the permeability for the solid volume fraction of χ = 0.6.
AN US
The results represent the normalized permeability as K/Kref for different boundary conditions in the viscosity range [0.029, 0.45] in the LB unit corresponding to the relaxation time of the range [0.58, 1.85], which covers both the over and underrelaxed modes.
schemes, namely, the CLI and LIBB, we find that the LIBB does not yield
M
viscosityindependent results with any of the collision models under consideration, while the CLI exhibits better properties in this respect and provides
ED
viscosity independent results with the TRT and MRT collision models. It is worth mentioning that the SRT collision scheme should not be used
PT
when the flow behavior is being studied across a wide range of Reynolds numbers, i.e., correlation of the drag force or the permeability, because the
CE
viscosity in the LBM is one of the key parameters that can be adjusted to simulate a certain Reynolds number of the flow. There are also some relaxation combinations proposed for the MRT model for high Reynolds
AC
number flow, which fixes the ghost relaxation rates to achieve stability. This also leads to viscosity dependent results. For further details, we refer to [29].
26
ACCEPTED MANUSCRIPT
3.1.3. Computational cost Our results obtained in Sections 3.1.1 and 3.1.2 indicate that the TRT
CR IP T
and MRT collision operators in combination with the MR boundary scheme results in the most accurate solutions. Given that the highly optimized split
kernel of the TRT implementation in waLBerla is about as fast as the SRT model, there is no justification to use the MRT scheme, as in the optimized version the computation is about a factor of two more expensive. Although
AN US
control of the higher order moments in high Reynolds numbers is only possible with the MRT collision scheme, the TRT collision model can also be used as it can control antisymmetric moments.
While the identification of a favorable collision operator is clearly apparent, the choice of the right boundary treatment is less obvious. It should be noted
M
that, especially for a DNS of turbulent flow, where one has to provide enough resolution to resolve the smallest dissipative scales, massive parallelism cannot
ED
be avoided. Unfortunately, the more accurate boundary handling methods must access LBM nodes from up to three cells away from the particle boundary
PT
cell, see Section 2.2. In a massively parallel setting, such physical boundary points may be near a logical processor boundary that has been created by the
CE
domain partitioning for a distributed memory machine [10]. In waLBerla this situation is handled by extra ghostlayer exchanges, i.e., by communicating a significantly extended set of distribution functions between neighboring
AC
processors along the subdomain boundaries. The data dependencies also cause additional synchronization overhead in the parallel execution. When saving both the pre and postcollision PDF, the nodes necessary for the interpolation along boundaries can be reduced by one. However, this still 27
ACCEPTED MANUSCRIPT
results in an additional communication for the extra saved field which again reduces the efficiency of the code. Even in a weak scaling scenario with only few spheres embedded in the flow and where large subdomains can
CR IP T
be handled on each processor, we have measured that the more advanced
boundary schemes may already cause a slowdown of between 10% and 50% as compared to the SBB method.
However, for our computational objectives, more complex geometric con
AN US
figurations must be considered. In these cases, the communication of data
quickly becomes the critical bottleneck. For our application, it is particularly important to save on communication, as for production runs many time steps are necessary to find the desired solution. Thus, to keep the computing times acceptable, only moderately large blocks of LBM nodes can be assigned to
M
each processor. In such a strongscaling regime, the ratio between computation and communication can easily become a bottleneck despite the highly
ED
optimized waLBerla program. This holds, even when just one layer of ghost nodes is exchanged in every time step and it is essential that further
PT
communication is avoided wherever possible. The efficiency of the LBM both in weak and strong scaling situations is analyzed in detail in [12].
CE
As a good compromise between the cost (including communication on parallel computers) and numerical accuracy, we choose here the TRTCLI scheme for large porescale problems. It leads to a reasonably good accu
AC
racy, has no viscositydependence, and needs less communication than the MR scheme. Hence, using the slightly less accurate but significantly better parallelizable method results in a considerable reduction in runtime. A more concise quantitative analysis of the performance tradeoffs will require the 28
ACCEPTED MANUSCRIPT
Table 2: Dimensionless drag force of steady laminar flow at Rep = 46 and χ = 0.9 for the SBB, CLI, and MR boundary schemes. All of the simulations are conducted with the TRT collision operator and the D3 Q27 lattice model.
CLI
MR
36
26.546
22.7539
22.3431
72
24.6073
22.6939
22.5754
144
23.6916
22.6731
22.6392
288
23.2805
22.6681
22.6535
1.17
1.79
1.86
AN US
α
CR IP T
SBB
D/∆x
development of sophisticated parallel performance models analogous to the techniques of Godenschwager et al. [12], Feichtinger et al. [58], but this is beyond the scope of the current paper. These articles study the program op
M
timization for the waLBerla software framework on various supercomputer architectures for simulations with simple boundary conditions.
ED
3.2. Porescale simulation at Re > 1 3.2.1. Convergence rate
PT
In this section, we investigate the consistency of the boundary condition schemes under steady laminar and weakly turbulent flows. To examine the
CE
convergence rate, we use the low Reynolds results from the previous sections to choose a first, second and thirdorder boundary scheme, namely the SBB,
AC
CLI, and MR, respectively. For all of the evaluation in this section, we use the TRT collision operator with the D3 Q27 velocity set and fix the magic number to 3/16. For the laminar steady flow, we simulate the fluid flow through the simple sphere pack at Rep = 46. In this regime, inertia of the fluid flow plays an 29
ACCEPTED MANUSCRIPT
Table 3:
Dimensionless drag force of weakly turbulent flow at Rep ≈ 315 and χ = 0.9
for the SBB, CLI, and the MR boundary schemes without correction term. All of the simulations are conducted with the TRT collision operator and the D3 Q27 lattice model.
SBB
CLI
MR
102.6
64.8742
64.8312
64.8332
144
61.6376
63.9441
63.8481
202.5
62.6216
63.3753
63.2392
288
64.0017
63.0316
63.7219

1.58

AN US
α
CR IP T
D/∆x
important role and the flow is nonlinear. The results of the dimensionless drag is presented in Table 2 for different boundary conditions. The SBB and CLI boundary schemes converge with first and secondorder, respectively,
M
as in the Stokes regime, while the convergence rate of the MR drops to secondorder. This is because of the secondorder accuracy of the LBM in
ED
bulk fluids in the nonlinear regime. The results also show that the SBB boundary scheme with the highest resolution conducted here differs from the
PT
CLI and MR schemes by 2.7%.
Here we investigate the capability of the LBM in the unsteady weakly
CE
turbulent flow at Rep ≈ 315. We conduct the simulation for a sufficiently long time such that the relative difference of the timeaveraged drag force is
below 0.1%. The results of the dimensionless drag force for various boundary
AC
schemes are shown in Table 3. Since the flow with this Reynolds number is fluctuating, we averaged the dimensionless drag force for 50 L/UD for each case. We use the diffusive scaling to keep the compressibility error in the same scale as the discretization error. Increasing the resolution beyond the 30
ACCEPTED MANUSCRIPT
necessary one for a stable simulation, the drag force computed in different boundary conditions does not change significantly. The MR boundary scheme in this test exhibits an instability that is associated with the correction
CR IP T
term. Since the MR model without correction term already yields a second order boundary condition and is stable, we will thus present the result of the MR model in this section without the correction term. The computed
convergence rate of the CLI based on the three point values is 1.58, while for
AN US
other boundary schemes an asymptotic convergence rate cannot be computed. The results suggest that for the simulation of turbulent flow in porous media, where high enough resolution is needed to capture the small vortices effects and reach a stable simulation, increasing the resolution does not significantly change the computed drag force. However, the boundary layer
M
effects at high Reynolds numbers should be investigated closely for different boundary schemes. In the above tests, the SBB results differ by at most 2%
ED
from the CLI and MR results. 3.2.2. Lattice model effect
PT
In this section, we evaluate the effect of the lattice models on turbulent flow and compare the results of the simulations conducted with two different
CE
lattice models of the D3 Q19 and D3 Q27 . The flow dynamics are obtained for touching spheres arranged in a simple cubic array. The magic number for the
AC
simulation is fixed at 3/16, and the TRT collision operator in combination
with the CLI boundary scheme is used. To simulate the touching sphere array (χ = 1.0), the flow field is initialized
with the result of the simulation using the SBB scheme and is continued with the CLI boundary scheme. This method of initialization helps to avoid 31
ACCEPTED MANUSCRIPT
the instability arising from the narrow gaps of the pores. For the boundary nodes for which the neighbor fluid cell does not exist, the SBB scheme is used instead of the CLI. By handling the boundary in this way, the convergence
CR IP T
rate will decrease but the approach is still stable, especially for unsteady flow. To calculate the time average Darcy velocity in the turbulent flow, we average the volume averaged velocity for a period of 100–150 L/UD flow through times.
AN US
In order to resolve the boundary layer, Tenneti et al. [59] point out
that the grid resolution in the DNS approaches should increase properly with increasing Reynolds number. The boundary layer thickness can be p approximated as δb ∼ D/ Rep . To check the gridindependency of the
turbulent flow computations, we conducted the simulation at Rep = 1045
M
with different resolutions. The results are presented in Table 4 which depict that for δb > 5.5, the timeaveraged dimensionless drag force (Cd ) is not
ED
changing significantly. For the following simulation, we use 252 lattice cells as sphere diameter.
PT
Figure 7 shows the instantaneous velocity field for the highest Reynolds number of 2477 that is considered in this study. The time evolution of the
CE
dimensionless drag force is also presented in this figure which shows the chaotic behavior of the flow. A more detailed analysis of the timeaveraged results gives that the same pressure drop will lead to different Darcy velocities
AC
for different lattice models. For D3 Q27 , the average velocity is 2.8% less than the D3 Q19 model and leads to the Reynolds number Rep = 2477, while for the D3 Q19 model Rep = 2547. On the other hand, the drag force calculated by these two models differs by only 0.2%. Taking into account the dimensionless 32
ACCEPTED MANUSCRIPT
Table 4: Dimensionless drag force of turbulent flow at Rep ≈ 1045 and χ = 1.0, and also
the approximate boundary layer thickness δb in lattice unit. The simulations are conducted
Rep
Cd
δb
150
1045
356.1
4.63
180
1043
355.2
5.57
252
1045
350.2
7.79
360
1048
348.9
11.12
AN US
D/∆x
CR IP T
with the TRT collision operator, the CLI boundary scheme and the D3 Q19 lattice model.
drag coefficient, the relative difference is 2.98%.
F/( AU2D)
ED
M
4.5
4
3.5 3
2.5
0.46
0.47 2
t /r
0.48
CE
PT
0.45
(a)
(b)
AC
Figure 7: Instantaneous velocity field contour of turbulent flow of touching spheres in simple periodic array arrangement at Rep = 2477 simulated with D3 Q27 (left), time series of the drag force simulated with different Rep (right). The simulations were conducted
with the TRTCLI and Λ = 3/16 with the same simulation parameters.
Figure 8 shows the relative differences between the D3 Q19 and D3 Q27 33
0.49
1
10
0
101
102
10
3
10
4
10
4
10
3
10
2
10
1
10
Rep
0
10
1
10
2
10
3
CR IP T
10
AN US
Relative difference %
ACCEPTED MANUSCRIPT
Figure 8: Relative difference of the dimensionless drag force CD,reduced /CD,f ull − 1% of
the two lattice models, D3 Q19 and D3 Q27 , reduced stencil and full stencil, respectively. The simulations were conducted with the same setup using the TRTCLI and Λ = 3/16.
models for different Reynolds numbers of the flow. Here we can see that for
M
Reynolds numbers below 100, the violation of the rotational invariance by the D3 Q19 model does not significantly affect the dimensionless drag force.
ED
In the regime beyond that Reynolds number the rotational invariance starts
PT
to exert an effect on the physical behavior of the flow. 3.3. Friction factor for the touching spheres
CE
At Reynolds numbers beyond the Stokes flow, or Darcy regime, Forchheimer [60] observed a nonlinear deviation from this relation, for which he
AC
proposed the addition of a quadratic term to the Darcy equation as follows −1/2
−1 ∇P = −µKD U − CF KD
ρUU,
(23)
where CF is the constant inertial factor, which mainly depends on the flow path and is usually found experimentally. A nondimensional form of Eq. 34
ACCEPTED MANUSCRIPT
(23) can be provided in the form of a friction factor, i.e., FK =
1 + F, ReK
(24)
CR IP T
√ √ where FK = (∆p/L) KD /ρU 2 , and ReK = ρU KD /µ.
Figure 9(a) presents the results of the friction factor for different Reynolds numbers of flow through touching spheres. As can be seen, for ReK < 1 the
results are inversely proportional to ReK which indicates the Darcy regime.
AN US
The deviation starts at ReK = 0.5 (Rep = 11) where the inertial effect of the
flow is starting to become dominant. The unsteady laminar flow starts when ReK > 6.9 (Rep > 139) and the transition regime from unsteady to chaotic flow can be identified in the range 18 < ReK < 28 (369 < Rep < 582). By further increasing the Reynolds number, we can see the turbulent flow regime
M
in which the friction factor converges to a nearly constant value of 0.16 for this porous geometry. The obtained friction factor is in line with the recent
ED
experimental results presented by Huang et al. [61] for Rep < 1000. In figure 9(a), the results of the simulation with the D3 Q19 lattice models
PT
are also presented. It is not easy to distinguish the difference between these two models, as the application of both models results in a similar trend for the
CE
friction factor versus Reynolds number. Figure 9(b) shows the friction factor results in the transition regime. It can be seen that when Rep > 139, the same pressure gradient in the boundary can result in different Darcy velocities
AC
for the two models, with the D3 Q19 , giving the highest Darcy velocity. The maximum relative difference in our simulation is below 3% which occurs with the highest Reynolds number.
35
102
102
1 0.8 0.6
D3Q27 D3Q19 1/ReK
FK
FK=0.16
0
102 102 100 ReK
102
M
104
Rep 102
D3Q19 1/ReK
FK=0.16
0.2
100
101 ReK
102
(b)
ED
(a)
103
D3Q27
0.4
FK
4
102 10
104
AN US
10
Rep 100
CR IP T
ACCEPTED MANUSCRIPT
Figure 9: Permeabilitybased friction factor versus permeabilitybased Reynolds numbers
PT
and particle diameterbased Reynolds numbers. a) whole range of Reynolds number considered in this study, b) a zoom view representing the transition regime. The simulations were performed for touching spheres in simple array arrangement. The plot shows the results
CE
of the D3 Q19 and D3 Q27 models when the simulations were conducted with TRTCLI
AC
method and Λ = 3/16.
36
ACCEPTED MANUSCRIPT
4. Turbulent flow over a porous bed By using the results of the evaluation which is presented in the previous
CR IP T
sections, we present a porescale simulation for semirealistic porous structure consisting spherical particles arranged in a random sphere packing. To
construct such a semireal porous structure, the inhouse physics engine
framework [62, 63] is used to simulate the rigid body interaction between spherical bodies. In this simulation, several spherical particles with the same
AN US
diameter are initialized such that they fall under their own weight and collect
at the bottom of the bed. After the particles come to rest, the geometry is imported as a permeable bed in the waLBerla framework for fluid flow
PT
ED
M
simulation.
(b)
CE
(a)
Figure 10: (a) Packing structure with monosized spherical particles, (b) Instantaneous
AC
velocity field contour at Rep = 65. The simulation was conducted with the TRTCLI and Λ = 3/16.
Figure 10 illustrates the packing structure containing 592 monosized
spherical particles in a domain with periodic stream and spanwise flow. At the top and bottom of the domain noslip conditions are applied so that part 37
ACCEPTED MANUSCRIPT
0.0045 0.004 0.0035
0.0025 0.002 0.0015 0.001 0.0005 20
40
60
80 100
AN US
Rep
CR IP T
K/D
2
0.003
Figure 11: Dimensionless permeability versus Reynolds number.
of the spheres are cut out of the domain. The diameter of the spheres is set to 40 lattice units, and we apply the TRT model of collision and the CLI
M
scheme for the wall boundary condition using the D3 Q19 model. The flow field is shown in the right picture with two 2D slices. The permeability is
ED
calculated at different Reynolds number in laminar flow that is presented in Fig. 11.
A free flow over a sphere bed, similar to [11] is also simulated to show
PT
a large scale simulation of a complex geometry. To resolve all scales, the porous part of the domain is refined using three different levels of refinement.
CE
Periodic boundary conditions have been used in the streamwise and spanwise direction while the top boundary is a free slip boundary. The simulation is
AC
executed on the LIMA supercomputer [64], with 64 nodes and 24 cores on each node. With this configuration it takes roughly 24 hours of computation for the whole simulation when using a spatial resolution of 9 × 107 cells. Fig.
12 shows the contours of the flow velocity of the simulation after 5 × 106 38
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Figure 12: Turbulent flow over a porous bed: From left to right, the block structure, a 2D velocity contour plot and a 2D cut through the fine lattice.
ED
timesteps. The skeleton grid resulting from the block structure as it is used in waLBerla is shown on the left hand side of the figure to indicate the
PT
static refinement of the meshes. The grids are also displayed on the right hand side of the contour to show the grid resolution. We point out that
CE
the results are stored after coarsening the grids throughout the domain by a factor of 8 to keep the size of the outputs reasonable for postprocessing and
AC
visualization.
5. Conclusions In this study, the lattice Boltzmann method is evaluated for the three
dimensional simulation of the flow through periodic simple sphere pack. A 39
ACCEPTED MANUSCRIPT
periodicpressure drop boundary condition has been adapted to drive the flow. This article studies different collision operators and various types of boundary conditions. Two different lattice models, namely, D3 Q19 and D3 Q27
CR IP T
have been examined for a wide range of Reynolds numbers of the flow.
The evaluation at low Reynolds numbers shows that the periodic pressure
boundary can accurately simulate porous media flows. The convergence study in the Stokes regime is investigated by computing the drag force with
AN US
various sizes of spheres while the solid volume fraction is fixed. The results show that the simple bounceback (SBB) boundary scheme converges with
first order, the linear interpolation bounceback (LIBB) converges between first and secondorder, the quadratic interpolation bounceback (QIBB) and the interpolationextrapolation bounceback (IEBB) converge with second
M
order and the multireflection method (MR) converges with thirdorder. The SRT collision operator produces a viscositydependent permeability for all
ED
boundary schemes, while the MRT and the TRT models can produce a viscosity independent permeability for the SBB, CLI, and MR boundary
PT
conditions. We also investigate the effect of the displacement of the center of the sphere on the drag force to show the source of scattering results at
CE
low resolution. The results show that, the use of the interpolating boundary offers smooth results, while the SBB is highly dependent on the number of solid cells representing the sphere.
AC
Since the TRT in its optimized version results in a similar runtime to the
SRT collision model, and MRT, which is computationally more expensive, does not significantly improve the accuracy, the TRT collision scheme has been used for the simulation of high Re number flows. 40
ACCEPTED MANUSCRIPT
Investigating the convergence of various boundary conditions with the TRT collision scheme in inertial flow show that the convergence rate of the MR boundary scheme decreases to second order accuracy while the CLI
CR IP T
and SBB maintain the same secondorder convergence rate as in the Stokes
regime. The degradation in the convergence rate of the MR results from the
secondorder inheritance of the LBM for the bulk flow in the inertial regime. Evaluating the boundary schemes in weakly turbulent flow at Rep ≈ 315
AN US
shows that the CLI converges with between the first and second order, while
calculating the convergence rate of the other boundary schemes is not possible because of the scattering results. The drag coefficient calculated with the SBB boundary condition differs by at most 2% compared to the results calculated with the CLI and MR boundary schemes.
M
Based on this simulation, the application of the CLI combined with the TRT can be considered as a good balance between accuracy and efficiency.
ED
Comparing the two lattice models with full and reduced stencils in turbulent flow at Rep = 2477, reveals the clear differences in the flow field caused by the
PT
lack of isotropy in the reduced scheme. However, porous media properties such as permeability and the drag force show less than 3% difference between these
CE
two lattice models. For flow at Rep < 100, the difference is negligible, and that beyond this Reynolds number, the difference increases with increasing Reynolds number.
AC
We also examined the ability of the model and the framework with two
large scale simulations, such as flow through random packing of spherical particles in the laminar regime, and the turbulent flow over a permeable bed.
The results demonstrate the capability of the model for wide range of flow 41
ACCEPTED MANUSCRIPT
regimes in complex flows.
CR IP T
Acknowledgements Financial support from the German Research Foundation (DFG, Project WO 671/111) and also the International Graduate School of Science and
Engineering (IGSSE) of the Technische Universität München for research training group 6.03 are gratefully acknowledged. Our special thanks goes to
AN US
Simon Bogner and Dominik Bartuschat for many fruitful discussions and the
waLBerla primary authors Florian Schornbaum, Christian Godenschwager and Martin Bauer for their essential help with the optimization of the code. UR and BW wish to thank the Institute of Mathematical Sciences at National
M
University of Singapore for their hospitality. References
ED
[1] S. Whitaker, Flow in porous media I: A theoretical derivation of Darcy’s law, Transport in porous media 1 (1986) 3–25.
PT
[2] S. Whitaker, The Forchheimer equation: a theoretical development,
CE
Transport in Porous media 25 (1996) 27–61. [3] E. Fattahi, M. Farhadi, K. Sedighi, H. Nemati, Lattice boltzmann
AC
simulation of natural convection heat transfer in nanofluids, International
Journal of Thermal Sciences 52 (2012) 137 – 144.
[4] S. Geller, S. Uphoff, M. Krafczyk, Turbulent jet computations based on MRT and cascaded lattice Boltzmann models, Computers and Mathematics with Applications 65 (2013) 1956 – 1966. 42
ACCEPTED MANUSCRIPT
[5] M. Singh, K. Mohanty, Permeability of spatially correlated porous media, Chemical Engineering Science 55 (2000) 5393 – 5403.
CR IP T
[6] J. Bernsdorf, G. Brenner, F. Durst, Numerical analysis of the pressure
drop in porous media flow with lattice Boltzmann (BGK) automata, Computer Physics Communications 129 (2000) 247 – 255.
[7] J. Kim, J. Lee, K.C. Lee, Nonlinear correction to Darcy’s law for a flow through periodic arrays of elliptic cylinders, Physica A: Statistical
AN US
Mechanics and its Applications 293 (2001) 13 – 20.
[8] A. Peters, S. Melchionna, E. Kaxiras, J. Lätt, J. Sircar, M. Bernaschi, M. Bison, S. Succi, Multiscale simulation of cardiovascular flows on the IBM Bluegene/P: Full heartcirculation system at redblood cell resolu
M
tion, in: Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing, Networking, Storage and Analysis,
ED
IEEE Computer Society, Washington DC, USA, 2010, pp. 1–10. [9] M. Schönherr, K. Kucher, M. Geier, M. Stiebler, S. Freudiger,
PT
M. Krafczyk, Multithread implementations of the lattice Boltzmann method on nonuniform grids for CPUs and GPUs, Computers and
CE
Mathematics with Applications 61 (2011) 3730–3743.
AC
[10] C. Feichtinger, S. Donath, H. Köstler, J. Götz, U. Rüde, WaLBerla: HPC software design for computational engineering simulations, Journal of Computational Science 2 (2011) 105 – 112.
[11] E. Fattahi, C. Waluga, B. Wohlmuth, U. Rüde, Large scale lattice boltzmann simulation for the coupling of free and porous media flow, 43
ACCEPTED MANUSCRIPT
in: High Performance Computing in Science and Engineering: Second International Conference, HPCSE 2015, Soláň, Czech Republic, May 2528, 2015, Revised Selected Papers, Springer International Publishing,
CR IP T
Cham, 2016, pp. 1–18.
[12] C. Godenschwager, F. Schornbaum, M. Bauer, H. Köstler, U. Rüde, A framework for hybrid parallel flow simulations with a trillion cells in
complex geometries, in: Proceedings of SC13: International Conference
AN US
for High Performance Computing, Networking, Storage and Analysis, SC ’13, ACM, New York, NY, USA, 2013, pp. 35:1–35:12.
[13] S. Bogner, S. Mohanty, U. Rüde, Drag correlation for dilute and moderately dense fluidparticle systems using the lattice Boltzmann method,
M
International Journal of Multiphase Flow 68 (2015) 71 – 79. [14] J. Götz, K. Iglberger, M. Stürmer, U. Rüde, Direct numerical simulation
ED
of particulate flows on 294912 processor cores, in: Proceedings of the 2010 ACM/IEEE International Conference for High Performance Computing,
PT
Networking, Storage and Analysis, SC ’10, IEEE Computer Society, Washington, DC, USA, 2010, pp. 1–11.
CE
[15] D. Bartuschat, D. Ritter, U. Rüde, Parallel multigrid for electrokinetic simulation in particlefluid flows, in: High Performance Computing and
AC
Simulation (HPCS), 2012 International Conference on, pp. 374 –380.
[16] A. T. White, C. K. Chong, Rotational invariance in the threedimensional lattice Boltzmann method is dependent on the choice of lattice, Journal of Computational Physics 230 (2011) 6367–6378. 44
ACCEPTED MANUSCRIPT
[17] G. Mayer, G. Házi, Direct numerical and large eddy simulation of longitudinal flow along triangular array of rods using the lattice Boltzmann method, Mathematics and Computers in Simulation (MATCOM) 72
CR IP T
(2006) 173–178.
[18] S. K. Kang, Y. A. Hassan, The effect of lattice models within the lattice
Boltzmann method in the simulation of wallbounded turbulent flows, Journal of Computational Physics 232 (2013) 100 – 117.
AN US
[19] G. Silva, V. Semiao, Truncation errors and the rotational invariance
of threedimensional lattice models in the lattice Boltzmann method, Journal of Computational Physics 269 (2014) 259 – 279. [20] R. W. Nash, H. B. Carver, M. O. Bernabeu, J. Hetherington, D. Groen,
M
T. Krüger, P. V. Coveney, Choice of boundary condition for latticeBoltzmann simulation of moderatereynoldsnumber flow in complex
ED
domains, Physical Review E 89 (2014) 023303. [21] T. Lee, C.L. Lin, An eulerian description of the streaming process in
PT
the lattice Boltzmann equation, Journal of Computational Physics 185
CE
(2003) 445 – 471.
[22] G. EitelAmor, M. Meinke, W. Schröder, A latticeBoltzmann method
AC
with hierarchically refined meshes, Computers & Fluids 75 (2013) 127 –
139.
[23] A. Fakhari, T. Lee, Numerics of the lattice Boltzmann method on nonuniform grids: Standard LBM and finitedifference LBM, Computers & Fluids 107 (2015) 205 – 213. 45
ACCEPTED MANUSCRIPT
[24] M. Bouzidi, M. Firdaouss, P. Lallemand, Momentum transfer of a Boltzmannlattice fluid with boundaries, Physics of Fluids (1994present)
CR IP T
13 (2001) 3452–3459. [25] R. Mei, W. Shyy, D. Yu, L.S. Luo, Lattice Boltzmann method for
3D flows with curved boundary, Journal of Computational Physics 161 (2000) 680 – 699.
[26] I. Ginzburg, D. d’Humières, Multireflection boundary conditions for
AN US
lattice Boltzmann models, Physical Review E 12 (2003) 666–614.
[27] C. Pan, L.S. Luo, C. T. Miller, An evaluation of lattice Boltzmann schemes for porous medium flow simulation, Computers & Fluids 35 (2006) 898 – 909.
M
[28] I. Ginzburg, F. Verhaeghe, D. d’Humières, TwoRelaxationTime lattice Boltzmann Scheme: About Parametrization, Velocity, Pressure and
ED
Mixed Boundary Conditions, Communications in Computational Physics 3 (2008) 427+.
PT
[29] S. Khirevich, I. Ginzburg, U. Tallarek, Coarseand finegrid numerical behavior of MRT/TRT lattice Boltzmann schemes in regular and random
CE
sphere packings, Journal of Computational Physics 281 (2015) 708–742.
AC
[30] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annual Review of Fluid Mechanics 30 (1998) 329–364.
[31] J. Zhang, D. Y. Kwok, Pressure boundary condition of the lattice Boltzmann method for fully developed periodic flows, Physical Review E 73 (2006) 047702. 46
ACCEPTED MANUSCRIPT
[32] S. H. Kim, H. Pitsch, A generalized periodic boundary condition for lattice Boltzmann method simulation of a pressure driven flow in a
CR IP T
periodic geometry, Physics of Fluids 19 (2007) –. [33] O. Gräser, A. Grimm, Adaptive generalized periodic boundary conditions for lattice Boltzmann simulations of pressuredriven flows through confined repetitive geometries, Physical Review E 82 (2010) 016702.
[34] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann equation:
AN US
theory and applications, Physics Reports 222 (1992) 145–197.
[35] S. Succi, The LatticeBoltzmann Equation, Oxford university press, Oxford, 2001.
M
[36] P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral onecomponent
ED
systems, Physical Review 94 (1954) 511–525. [37] X. He, L.S. Luo, Lattice Boltzmann model for the incompressible Navier–
PT
Stokes equation, Journal of Statistical Physics 88 (1997) 927–944. [38] I. Ginzburg, Consistent lattice Boltzmann schemes for the Brinkman
CE
model of porous flow and infinite ChapmanEnskog expansion, Physical
AC
Review E 77 (2008) 066704.
[39] I. Ginzburg, P. M. Adler, Boundary flow condition analysis for the threedimensional lattice Boltzmann model, J. Phys. II France 4 (1994)
191–214.
47
ACCEPTED MANUSCRIPT
[40] D. d’Humières, Multiple–relaxation–time lattice Boltzmann models in three dimensions, Philosophical Transactions of the Royal Society of (2002) 437–451.
CR IP T
London. Series A: Mathematical, Physical and Engineering Sciences 360
[41] I. Ginzburg, Lattice Boltzmann modeling with discontinuous collision components: Hydrodynamic and advectiondiffusion equations, Journal of Statistical Physics 126 (2007) 157–206.
AN US
[42] I. Ginzburg, F. Verhaeghe, D. d’Humières, Study of simple hydrodynamic solutions with the tworelaxationtimes latticeBoltzmann scheme, Communications in Computational Physics 3 (2008) 519+. [43] D. Yu, R. Mei, L.S. Luo, W. Shyy, Viscous flow computations with the
ED
39 (2003) 329 – 367.
M
method of lattice Boltzmann equation, Progress in Aerospace Sciences
[44] P. Lallemand, L.S. Luo, Lattice Boltzmann method for moving bound
PT
aries, Journal of Computational Physics 184 (2003) 406 – 421. [45] B. Chun, A. J. C. Ladd, Interpolated boundary condition for lattice
CE
Boltzmann simulations of flows in narrow gaps, Physical Review E 75
(2007) 066705.
AC
[46] O. Filippova, D. Hänel, Grid refinement for latticeBGK models, Journal of Computational Physics 147 (1998) 219–228.
[47] L. Talon, D. Bauer, N. Gland, S. Youssef, H. Auradou, I. Ginzburg, Assessment of the two relaxation time latticeBoltzmann scheme to 48
ACCEPTED MANUSCRIPT
simulate stokes flow in porous media, Water Resources Research 48 (2012) n/a–n/a. W04526.
CR IP T
[48] A. Dupuis, From a Lattice Boltzmann model to a parallel and reusable
implementation of a virtual river, Ph.D. thesis, University of Geneva, 2002.
[49] N. Martys, X. Shan, H. Chen, Evaluation of the external force term in the
AN US
discrete Boltzmann equation, Physical Review E 58 (1998) 6855–6857.
[50] J. M. Buick, C. A. Greated, Gravity in a lattice Boltzmann model, Physical Review E 61 (2000) 5307–5320.
[51] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the forcing term in
M
the lattice Boltzmann method, Physical Review E 65 (2002) 046308. [52] A. Mohamad, A. Kuzmin, A critical evaluation of force term in lattice
ED
Boltzmann method, natural convection problem, International Journal of Heat and Mass Transfer 53 (2010) 990 – 996.
PT
[53] H. Huang, M. Krafczyk, X. Lu, Forcing term in singlephase and shanchentype multiphase lattice Boltzmann models, Physical Review E 84
CE
(2011) 046710.
AC
[54] H. Darcy, Recherches expérimentales relatives au mouvement de l’eau dans les tuyaux, MalletBachelier, Paris (1857).
[55] A. Sangani, A. Acrivos, Slow flow through a periodic array of spheres, International Journal of Multiphase Flow 8 (1982) 343 – 360.
49
ACCEPTED MANUSCRIPT
[56] D. Bartuschat, U. Rüde, Parallel multiphysics simulations of charged particles in microfluidic flows, Journal of Computational Science 8 (2015)
CR IP T
1 – 19. [57] P. M. Adler, Porous media : geometry and transports, ButterworthHeinemann Limited, 1992.
[58] C. Feichtinger, J. Habich, H. Köstler, U. Rüde, T. Aoki, Performance modeling and analysis of heterogeneous lattice boltzmann simulations
AN US
on cpu–gpu clusters, Parallel Computing 46 (2015) 1–13.
[59] S. Tenneti, R. Garg, S. Subramaniam, Drag law for monodisperse gas–solid systems using particleresolved direct numerical simulation of flow past fixed assemblies of spheres, International Journal of Multiphase
M
Flow 37 (2011) 1072 – 1092.
ED
[60] P. Forchheimer, Wasserbewegung durch Boden, Z. Ver. Deutsch. Ing 45 (1901) 1782–1788.
PT
[61] K. Huang, J. Wan, C. Chen, L. He, W. Mei, M. Zhang, Experimental investigation on water flow in cubic arrays of spheres, Journal of
CE
Hydrology 492 (2013) 61 – 68.
[62] K. Iglberger, U. Rüde, Massively parallel rigid body dynamics simulations,
AC
Computer ScienceResearch and Development 23 (2009) 159–167.
[63] T. Preclik, U. Rüde, Ultrascale simulations of nonsmooth granular dynamics, Computational Particle Mechanics (2015) 1–24.
50
ACCEPTED MANUSCRIPT
[64] R. R. Erlangen, Lima cluster, https://www.rrze.fau.de/dienste/
AC
CE
PT
ED
M
AN US
CR IP T
arbeitenrechnen/hpc/systeme/limacluster.shtml, 2016.
51