Next Article in Journal
The INUS Platform: A Modular Solution for Object Detection and Tracking from UAVs and Terrestrial Surveillance Assets
Previous Article in Journal
The Performance of a Gradient-Based Method to Estimate the Discretization Error in Computational Fluid Dynamics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Revisiting the Homogenized Lattice Boltzmann Method with Applications on Particulate Flows

1
Lattice Boltzmann Research Group, Karlsruhe Institute of Technology, Straße am Forum 8, 76131 Karlsruhe, Germany
2
Institute for Mechanical Process Engineering and Mechanics, Karlsruhe Institute of Technology, Straße am Forum 8, 76131 Karlsruhe, Germany
3
Institute for Applied and Numerical Mathematics, Karlsruhe Institute of Technology, Englerstraße 2, 76131 Karlsruhe, Germany
*
Author to whom correspondence should be addressed.
Submission received: 22 December 2020 / Revised: 18 January 2021 / Accepted: 20 January 2021 / Published: 27 January 2021
(This article belongs to the Section Computational Engineering)

Abstract

:
The simulation of surface resolved particles is a valuable tool to gain more insights in the behaviour of particulate flows in engineering processes. In this work the homogenized lattice Boltzmann method as one approach for such direct numerical simulations is revisited and validated for different scenarios. Those include a 3D case of a settling sphere for various Reynolds numbers. On the basis of this dynamic case, different algorithms for the calculation of the momentum exchange between fluid and particle are evaluated along with different forcing schemes. The result is an updated version of the method, which is in good agreement with the benchmark values based on simulations and experiments. The method is then applied for the investigation of the tubular pinch effect discovered by Segré and Silberberg and the simulation of hindered settling. For the latter, the computational domain is equipped with periodic boundaries for both fluid and particles. The results are compared to the model by Richardson and Zaki and are found to be in good agreement. As no explicit contact treatment is applied, this leads to the assumption of sufficient momentum transfer between particles via the surrounding fluid. The implementations are based on the open-source C++ lattice Boltzmann library OpenLB.

1. Introduction

The simulation of particulate flows finds application in many fields as it permits a detailed examination of an engineering process. It allows access to data on particle behaviour, for which experimental measuring is complex or costly. It was already applied for the simulation of solid separation processes by Viduka et al. [1], testing different pulsation profiles of a jigging device, or by Li et al. [2] regarding the separation of soybeans from mustard seeds. The work presented here is the validation and improvement of a method for direct numerical simulation, also depicting the particle’s shape. The industrial relevance of this topic has been shown, e.g., by Champion et al. [3] who showed that particle shape significantly influences the performance of drug carriers. A discussion on the industrial relevance of particle shape for products and processes is also given by Davies [4], who references, e.g., the relevance for the production of rubber, as also stated by Scotti et al. [5].
On the most coarse level, particles can be depicted as a continuum by considering distributions with Euler–Euler approaches utilising an advection–diffusion equation [6,7,8]. Approaches which treat each particle as separate entity yield more details. With the discrete element method (DEM), e.g., separation processes in a cyclone can be studied [9] considering the particles as solid spheres. Furthermore, adsorption processes in a static mixer [10] or packing characteristics of powder [11] are of interest. The DEM can be used together with various approaches and models for the simulation and representation of the fluid. The coupling of DEM particles with a model for a non-Newtonian turbulent fluid is for example relevant in the simulation of anaerobic digestion [12]. A review of applications for DEM simulations is given by Zhu et al. [13].
For more sophisticated studies, hydrodynamic forces on a particle and its shape can be resolved by direct numerical simulations (DNS) [14]. This enables, e.g., the calculation of drag correlations [15]. Focusing on the shape, the investigation of the interaction between particles and a rough surface is possible [16]. At this point, however, proper collision treatment becomes important. For the absence of a significant influence of a fluid, the framework Grains3D [17,18] has been presented to investigate arbitrarily shaped convex particles under the influence of gravity and contact forces applying the Gilbert–Johnson–Keerthi algorithm for contact detection and distance calculation. The framework has later been extended to concave shapes by representing them as glued convex shapes [19]. Considering particles submersed in fluid, a simple approach is the one of glued spheres [20] to approximate various particle shapes and get more information on the distribution of acting forces. For this approximation, however, the accuracy is limited as the required number of spheres increases drastically for better approximations of a shape. To circumvent this problem, the object can be described by several Lagrange points on the surface, at which the hydrodynamic forces are calculated. This points can be independent of the underlying grid, which represents the fluid. This allows for a high accuracy in depiction of the shape, but requires frequent interpolation between fluid and particle points. This common approach is called immersed boundary method (IBM) [21]. One of its advantages is that it can be coupled with various approaches to solve the fluid system like the finite element method, the finite volume method or the lattice Boltzmann method (LBM) [22,23].
The coupled IBM–LBM has been frequently applied for the investigation of fluid structure interaction problems [24]. While many studies of, e.g., vortex induced vibrations [22,25] or the flow around a torus by Wu and Shu [26] only consider simpler geometries with analytical representations, mappings for more complex polygons exist, as discussed by Owen et al. [27]. Nevertheless, more complex structures can be simulated as well as described by Beny and Latt [28], who simulated multiple propellers on a GPU system. In IBM–LBM, it is possible to also consider simple strings as shown by Tian et al [22] who simulated filaments in a flow. A comparative study of IBM approaches was performed by Kang and Hassan [29], who investigated the flow past a cylinder and decaying Taylor–Green vortices.
The homogenized lattice Boltzmann method (HLBM) used in this work, in turn, has very few restrictions regarding the shape, since a voxel representation of almost any object can be chosen for simulations, as shown in [30]. Very thin objects, however, remain problematic, as they may not be captured by the lattice resolution. The main difference between IBM and HLBM is the way the objects are represented. The former employs Lagrangian points, which allows thin structures but requires interpolations and the creation of a distribution of those points on the surface of the object. The latter directly uses points on the fluid grid. The correct depiction of a structure, in this case, mainly depends on the grid resolution.
The LBM proved to be easy to parallelize due to the fact that most calculations are carried out in a strictly local collision step. This has been investigated for the use on computing clusters with GPU systems [31] as well as mixed CPU–GPU systems [32]. For LBM, additional approaches for the simulation of arbitrarily shaped particles are available, which are specific to their method. The most common one is the partially saturated cells method (PSM), in its original version proposed by Noble and Torczynski [33]. It depicts the object on the fluid grid by assigning a linear combination of a no-slip condition and bulk flow to a grid cell depending on its distance to the actual physical particle boundary. This is possible as LBM operates on a mesoscopic level considering the fluid by distributions of fluid-particles in a phase space. There are several aspects to this approach if moving objects are considered, e.g., the refilling of cells that were previously covered by a particle or the chosen no-slip condition. This has been investigated by Peng et al. [34].
Another topic is the way the hydrodynamic force is calculated. The stress integration method, e.g., described in [35,36], is also applicable for other simulation approaches and not specific to LBM but proved to be inefficient [37,38]. An alternative is the momentum exchange algorithm (MEA) [39,40], which calculates differences in the momentum via the mentioned particle distributions in the phase space. For this method, various formulations and improvements have been proposed [41,42,43,44] and are subject to current investigations [34].
In this study three applications are presented. For the first case of a settling sphere, the results are validated against ten Cate et al. [45], who performed simulations of a single sphere for different Reynolds numbers and compared the data to particle imaging velocimetry experiments. In addition, the results are compared to various drag correlations discussed in Section 2.1.
The second case simulates the tubular pinch effect first described by Segré and Silberberg [46]. They found that a neutrally buoyant particle in a tube flow equilibrates at a position between the tube’s center and its wall and that the results for single particles can be extended to mixtures by linear combination [47]. Later, Tachibana et al. [48] found that the equilibrium position depends on the ratio of sphere diameter to tube diameter, which should be larger than 0.2 for the effect to be clearly visible. Karnis et al. [49] studied further the influence of the Reynolds number and the distance between two spheres on this effect.
The case of hindered settling is also relevant in many applications of process engineering as usually collectives of particles are processed. This is reflected in many studies simulating this case [50,51,52,53]. The results are usually compared to an empirical correlation found by Richardson and Zaki [54]. However, more correlations exist as discussed in Section 2.2. The results in this study are compared to different correlations to give a better overview.
In this work, the HLBM proposed by Krause et al. [30,55] is applied. Contrary to other approaches for moving objects, HLBM does not represent objects by a sharp no-slip boundary but rather with interior fluid utilising a forcing scheme and a transition region to mimic such a condition. This is a similarity to the IBM and leads to the need for further research as existing MEA approaches are developed for no-slip boundaries in the first place. Therefore, a novel momentum loss algorithm (MLA) is tested in this study. The semi-locality (depending on the chosen approach for the calculation of hydrodynamic forces) of the HLBM allows for an easy and efficient implementation on parallel systems without the need of costly interpolations. On the other hand, no refilling of cells is required as in PSM approaches. Besides the new MLA the HLBM is revisited and for the first time evaluated for different forcing schemes and methods for the calculation of hydrodynamic forces. To the knowledge of the authors such a comparative study for HLBM has not been performed before. This finally leads to a new improved scheme, which is validated for the simple case of a settling sphere but also tested for the application in hindered settling simulations and the tubular pinch effect. The results are compared to common experimentally determined correlations and reference simulations. This shows the range of applicability of the presented method, e.g., regarding the Reynolds number. Such tests are currently not found in literature for HLBM, however, this information is of interest when choosing a proper simulation approach depending on an application. Therefore, additionally different cases are examined for the first time with HLBM, which is important as, e.g., the tubular pinch effect can’t be reproduced by all approaches for DNS particle simulations as stated by Li et al. [35] and Peng et al. [34]. Additionally, considering the case of hindered settling the investigation of particle interactions, first mentioned by Krause et al [55], is continued, which was not found in other publications. To the knowledge of the authors, the investigations and findings for HLBM presented in this work have not been shown before.
The remainder of this paper is organized as follows: in Section 2 the physical background and correlations found in the literature are discussed. Then, in Section 3 the methods used for the simulations are discussed along with various forcing and MEA schemes to be tested. Finally, in Section 4 the results of the numerical investigation are presented and the conclusion is drawn in Section 5. All simulations presented are performed with the open-source lattice Boltzmann C++ framework OpenLB [56].

2. Modelling

This study focuses on two component systems of rigid particles submersed in a water-like fluid. The latter is described by the incompressible Navier–Stokes equations, given by
u f t + ( u f · ) u f ν Δ u f + 1 ρ f p = F f in Ω × I , · u f = 0 in Ω × I .
The equations are defined on a spatial domain Ω R d of dimension d 2 , 3 and a time interval I R . The fluid velocity is denoted by u f : Ω × I R d , while p : Ω × I R describes the pressure, ρ f R > 0 the fluid’s density and ν R > 0 its kinematic viscosity. F f : Ω × I R d represents the total of external forces acting on the fluid.
The dynamics of the second component, the rigid particles, are governed by the equations based on Newton’s second law of motion
m p u p t = F p and J p ω p t = T p .
With the particle mass m p R > 0 , the particle velocity u p : Ω × I R d and the force acting on the particle F p : Ω × I R d , the trajectory of motion can be calculated. The rotational behaviour is modelled analogously with the moment of inertia J p R d ^ , the angular velocity ω p : Ω × I R d ^ and the torque T p : Ω × I R d ^ . Here d ^ = 1 for d = 2 and d ^ = 3 for d = 3 .
In this study, only the gravitational and buoyancy forces are considered as external forces. They are given by F G = ( 0 , 0 , m p g ) and F B = ( 0 , 0 , m p ρ f ρ p g ) for d = 3 ( F G = ( 0 , m p g ) and F B = ( 0 , m p ρ f ρ p g ) for d = 2 ) for a gravitational acceleration of g. This means, forces regarding particle-particle and particle-wall contact are neglected in this study and momentum between particles is only transferred via the fluid. As most cases consider only single particles with no relevant wall contact this is reasonable. The only exception is the investigation of hindered settling in Section 4.3. Here the system is discretized with a decent resolution to allow for sufficient momentum transfer via the fluid. Along with the hydrodynamic force F H : Ω × I R d accounting for the momentum transfer between fluid and particles, the combined forces are given by
F f = F H , F p = F G F B F H and T p = r × F H ,
for a vector r R d denoting the distance to the center of mass of a particle.

2.1. Drag Force

While in this study a DNS approach is chosen, for methods considering the particles as single points the hydrodynamic force has to be modelled. This is done via a drag force, accounting for the resistance an object experiences when moving relative to a surrounding fluid. It is defined by
F D = 1 2 ρ f ( u p u f ) 2 C D A ,
with A R > 0 being the area of the projection of the object in flow direction. While the determination of the drag coefficient C D R for various shapes is subject to current research [57,58], for spheres many correlations depending on the Reynolds number Re have been given. The latter is defined as
Re = d p u p ν ,
for the cases considered in this work, with a particle’s diameter d p R > 0 . For small Reynolds numbers the drag coefficient is described by Stokes’ law [59]
C D = 24 Re , for Re < 1 .
While in the Newton regime between Re = 10 3 and Re = 10 5 the coefficient is defined by C D = 0.44 , many correlations exist for the transition region of Reynolds numbers between 1 and 10 3 . Notable here is among others the contribution by Abraham [60], who took the boundary layer into account for a theoretical derivation, leading to
C D = 0.5407 + 24 Re 2 , for 0 < Re < 5000 .
Schiller and Naumann [61] used an empirical approach to obtain the commonly used approximation of
C D = 24 Re 1 + 0.15 Re 0.687 , for Re < 800
for the drag coefficient. As many more correlations can be found, the authors refer to Dey et al. [62] for a comprehensive overview.

2.2. Hindered Settling

The above mentioned drag correlations consider the case of free settling, more precisely the settling of a single sphere in an infinitely extended medium. In practical applications, collectives of particles are more common, which leads to the case of hindered settling caused by interactions between particles. The empirical correlations developed in many past studies proved to be sensible approaches in modeling such behaviour.
The oldest recapitulated here is the one by Steinour [63], who investigated spheres consisting of tapioca and glass, for a solid volume fraction ϕ of up to 49.8 % . The experimental conditions correspond to Re = 0.0025 and Re = 0.0026 , leading to
u Steinour ( ϕ ) = u S ( 1 ϕ ) 2 e 4.19 ϕ .
The hindered settling velocity is here expressed relative to the terminal settling velocity according to Stokes’ u S , obtained by inserting Equation (6) in Equation (4).
Another correlation was found by Oliver [64]. In the experiments with Kallodoc particles, Reynolds numbers up to 0.39 were reached with solid volume fractions up to 35 % . The results are reflected in the formula for the hindered settling velocity given by
u O ( ϕ ) = u S 1 0.75 ϕ 1 3 ( 1 2.15 ϕ ) .
One of the most commonly used expressions was found in investigations performed by Richardson and Zaki, who approached the topic from both, an analytical [54] and an experimental [65] perspective. Their formulation was subject to many adaptions regarding accuracy, convenience and extension of applicability [66,67,68,69]. In the original form it is given by u RZ ( ϕ ) = u S ( 1 ϕ ) n , with
n = 4.65 for Re < 0.2 4.35 Re 0.03 for 0.2 < Re < 1 4.45 Re 0.1 for 1 < Re < 500 2.39 for 500 < Re .
In a more comprehensive study taking also the results of other authors into account along with their findings, Barnea and Mizrahi [70] concluded with the general correlation
u BM ( ϕ ) = u S 1 ϕ ( 1 + ϕ 1 3 ) e 5 ϕ 3 ( 1 ϕ ) ,
which is independent of the Reynolds number. While this is the finding for the creeping flow regime ( Re 1 ), Barnea and Mizrahi suggest the application of drag correlations beyond Stokes’ law. Therefore in contrast to other publications, for higher Re the hindered settling velocity is not based on u S anymore, but on a terminal settling velocity determined, e.g., by Equation (8).

3. Methods

For the simulations, the lattice Boltzmann method was applied, which operates on a uniform cubical grid Ω h approximating the computational domain. This structure is, together with the discrete phase space, spanned by the d spatial dimensions and a discrete set of q velocities c i R d , i = 0 , , q 1 , denoted as lattice. The latter is usually labeled as DdQq. Choosing commonly used sets, the simulations in this study were performed on a D2Q9 and a D3Q19 lattice. As usual, the computations were conducted in lattice units meaning all values were scaled such that the grid spacing and time step size became δ x L = 1 and δ t L = 1 , respectively. The superscript L indicates that a value is given in lattice units. For reasons of readability, all values in this section were assumed to be in lattice units if not stated otherwise omitting the superscript.
The propagation of particle distributions f i : Ω h × I h R , i = 0 , , q 1 , for a discretized time interval I h , on the given lattice is described by the lattice Boltzmann equation. Together with the Bhatnagar—Gross—Krook (BGK) collision operator [71] it reads
f i ( x + c i , t + 1 ) f i ( x , t ) = 1 τ f i ( x , t ) f i eq ( ρ , u ) , for x Ω h , t I h , i = 0 , , q 1 .
The lattice relaxation time τ is related to the lattice kinematic viscosity by
ν = c s 2 τ 1 2 .
Beside a streaming of particle distributions, Equation (13) describes the relaxation of a system towards an equilibrium state given by the discrete Maxwell–Boltzmann distribution
f i eq ( ρ , u ) = w i ρ 1 + c i · u c s 2 + ( c i · u ) 2 2 c s 4 + u 2 2 c s 2 , for i = 0 , , q 1 .
Here, w i are weights originating from a Gauss–Hermite quadrature rule and c s = 1 3 is the lattice speed of sound, valid for both lattice configurations used in this work. The density and velocity on a grid point can be calculated by moments of the discrete distributions
ρ ( x , t ) = i = 0 q 1 f i ( x , t ) and u ( x , t ) = 1 ρ i = 0 q 1 c i f i ( x , t ) , for x Ω h , t I h ,
respectively. Furthermore, the pressure p is related to the density ρ in LBM by p ( x , t ) = ρ ( x , t ) c s 2 . In the rest of this section, the arguments for density and velocity are omitted for better comprehensibility. For a more comprehensive overview on this simulation approach, the authors refer to [23].

3.1. Homogenized Lattice Boltzmann Method

The DNS simulations in this study were performed utilizing the HLBM, which was first published by Krause et al. [55] and was extended to 3D by Trunk et al. [30]. The extensions of the method to the standard LBM with BGK collision could be divided into three parts. Namely, the representation of an object on the lattice, the forcing scheme applied to the fluid to account for the particles’ presence and a method to calculate the exchanged momentum between an object and fluid for the calculation of forces acting on the object. From this point, the trajectory of a particle is calculated according to Equation (2) in combination with a velocity-Verlet algorithm [72,73].
Object Representation. The depiction of arbitrary shapes within the HLBM is described in detail by Trunk et al. [30]. Here, only a brief overview is given. From a voxel representation of the object, which did not necessarily have the same spatial resolution as the lattice, parameters like the moment of inertia and the mass were calculated. For a smooth transition between fluid and particle domain, a Gaussian filter was applied. If an analytical representation of the shape was given, a smooth transition could also be defined utilizing trigonometric functions as described by Krause et al. [55]. This led to a mapping function d B : Ω × I [ 0 , 1 ] , which took on the value of 0 in the fluid domain and 1 in the area covered by the object. In the transition region it yielded values in between. Along with the current position of the object x p Ω this defined a domain B ( t ) = { x Ω : d B ( x , t ) 0 } which was covered by the body. To account for the presence of particles, the velocity used in the equilibrium distribution Equation (15) was replaced by the convex combination
u ˜ ( x , t ) = u ( x , t ) + d B ( x , t ) u p ( x , t ) u ( x , t ) for x B ( t ) , t I h .
The velocity difference Δ u = u ˜ u could be used in forcing schemes as stated in the next Section.
Forcing Schemes. As the effect on the fluid is given via a velocity difference, in previous publications a forcing scheme according to Shan and Chen [74] has been chosen, which only modifies the velocity inserted in the Maxwell–Boltzmann distribution. In this study, several forcing schemes (proposed by Shan and Chen [74], Guo et al. [75] and Kupershtokh et al. [76]) were considered (see Section 3.2) and then tested in Section 4.1.
Momentum Exchange. Various methods to calculate the exchanged momentum are given in Section 3.3 and are applied in Section 4.1. They were evaluated regarding the accuracy of the velocity profile of a falling sphere. In previous publications [30,55] an adapted version of the momentum exchange according to Ladd [39,40] was used.

3.2. Forcing Schemes

To incorporate a force in LBM, the definition of the velocity handed to the Maxwell–Boltzmann distribution, from now on labeled u eq , was modified and/or a source term S i was added to the right hand side of Equation (13). In all cases, the fluid velocity given in Equation (16) was redefined in the presence of a force F to
u = 1 ρ i = 0 q 1 c i f i ( x , t ) + F δ t 2 ρ , for x Ω h , t I h .
For a more comprehensive discussion on forcing schemes the authors again refer to Krüger et al. [23], however, the three forcing schemes used in this study are briefly summarized here.
The approach of Shan and Chen [74] only modifies the equilibrium velocity to
u eq = 1 ρ i = 0 q 1 c i f i ( x , t ) + τ F ρ , for x Ω h , t I h .
This scheme originates from methods for multi-component and multi-phase flows, but is also applicable to single-component systems. It is frequently used due to its simplicity and performance advantage, however, a dependency of the results on τ is found, e.g., by Huang et al. [77] examining various forcing schemes for multi-phase LBM. Along with this scheme, a velocity shift method is proposed by Shan and Chen [74], which states F = ρ τ Δ u . Guo et al. [75] proposed a scheme modifying both the velocity and the source term to
u eq = 1 ρ i = 0 q 1 c i f i ( x , t ) + F δ t 2 ρ , for x Ω h , t I h ,
S i ( x , t ) = 1 δ t 2 τ w i c i u c s 2 + ( c i · u ) c i c s 4 · F , for x Ω h , t I h .
At last, the exact difference method (EDM) by Kupershtokh et al. [76] only introduces a source term S i ( x , t ) = f i eq ( ρ , u + Δ u ) f i eq ( ρ , u ) to the system. The velocity difference Δ u it is based upon is calculated by Δ u = F δ t / ρ , since δ t = δ x = 1 in lattice units as stated before. In this approach the velocity u remains just as defined in Equation (16).

3.3. Methods for Momentum Exchange

In this study, different schemes for the momentum exchange between fluid and particle are also discussed. The first is an adaptation of the MEA by Ladd [39,40] which has been applied in [30,55]. As shown in [37,38], it can be written as
F H ( t ) = x B ( t ) i = 1 q 1 c i f i ( x + c i , t ) + f i ¯ ( x , t ) , for t I h .
Here f i ¯ is defined as distribution function according to the velocity c i ¯ , which is opposite to c i , i.e., c i ¯ = c i . Various Improvements have been suggested [41,78], as this approach fails, e.g., in the simulation of the tubular pinch effect described by Segré and Silberberg [46,47], as stated by Peng et al. [34]. In their study, they compared different MEA approaches and tested them along with other aspects of moving boundary implementations. One of the most common approaches is the one by Wen et al. [42]
F H ( t ) = x B ( t ) i = 1 q 1 c i u p f i ( x + c i , t ) + c i + u p f i ¯ ( x , t ) , for t I h .
It incorporates the particle velocity for higher accuracy and in order to achieve Galilean invariance.
Lastly, a new approach is proposed by the authors denoted as momentum loss algorithm. As in the HLBM objects are depicted with interior fluid rather than with a sharp no-slip boundary, the investigation for alternative approaches to calculate exchanged momentum is of interest. Examined in this work is the method of computing the momentum directly from the introduced forcing for the approaches by Kupershtokh et al. [76] and Shan and Chen [74]. For the latter this means
F H ( t ) = 1 τ x B ( t ) i = 1 q 1 c i f i eq ( ρ , u ˜ ) f i eq ( ρ , u ) , for t I h ,
while for the EDM it is slightly different namely,
F H ( t ) = x B ( t ) i = 1 q 1 c i f i eq ( ρ , u ˜ ) f i eq ( ρ , u ) , for t I h .
The forcing schemes and methods for momentum exchange discussed in this section allow for a comparative study. The previously used HLBM approach is based on the forcing introduced by Shan and Chen [74] and the momentum exchange by Ladd [39,40]. According to literature [34], there are methods available which yield a higher accuracy. This will be investigated in Section 4.1 leading to an updated version of HLBM which is expected to yield higher accuracy.

4. Results and Discussion of Numerical Experiments

In this Section, various schemes for forcing and momentum exchange are tested along with the remaining HLBM implementation. In Section 4.1 first, the approaches are compared regarding the settling velocity of a sphere and the best performing combination is selected for further validation against existing drag correlations for spheres. Afterwards, the implementation is applied to the cases of a neutrally buoyant particle in a pipe in Section 4.2 and hindered settling in Section 4.3, respectively.

4.1. Settling Sphere

In this section, simulation results are first compared to experimental data by ten Cate et al. [45] for evaluation and validation of the chosen methods and subsequently against empirical correlations discussed in Section 2.1.

4.1.1. Simulation Setup—Comparison to Literature

The settling of a single sphere has been investigated by ten Cate et al. [45] both, experimentally and numerically. Therefore a sphere of diameter d p = 0.015   m with density ρ p = 1120   kg   m 3 was placed 0.12   m from the bottom in a container of height 0.16   m and a width and depth of 0.1   m , as shown in Figure 1. This setup was used to examine the settling for different Reynolds numbers by varying the fluid density ρ f and dynamic viscosity μ f . The same setup was chosen for the simulations with different discretization parameters for each case. The values were given along with the simulation results of ten Cate et al. [45] for a resolution of four cells per particle radius in Table 1. The Reynolds number was given in this case according to the values in the table by Re = d p u ρ f / μ f . For the boundaries, a no-slip halfway bounce-back condition [79] was chosen.
In this study, the given setup is investigated for different combinations of forcing and momentum exchange schemes described in Section 3.2 and Section 3.3. Additionally, a convergence study is performed.

4.1.2. Results and Discussion—Comparison to Literature

For reasons of readability, further abbreviations are introduced. Here, GUO and SCF refer to the respective forcing schemes by Guo et al. [75] and Shan and Chen [74] (see Section 3.2). Additionally, MEA-L and MEA-W refer to the momentum exchange algorithms by Ladd [39,40] and Wen et al. [42]. As different resolutions are studied a factor N is introduced to scale the effective grid spacing to δ x / N .
From the simulations, the root mean squared error (RMSE) relative to the maximum absolute velocity of the respective case and the deviation of the maximum settling velocity compared to the one calculated according to the drag correlation by Abraham [60] are given in Table 2 for N = 1 and in Table 3 for N = 8 . To further compare the whole curves, similarity measures were calculated [80], i.e., a partial curve mapping (PCM) [81] and the area between curves according to Jekel et al. [80]. The latter was less influenced by noisy data or outliers. As reference solution in this calculation, the experimental data of ten Cate et al. [45] were used and approximated by a polynomial fit of 18-th order. The RMSE of the fit for case 1 to 4 was given by 0.01447 , 0.01481 , 0.01838 , and 0.02093 . The values were rounded to six decimal places to demonstrate that some differences (e.g., between MEA-L and MEA-W) were on a negligible level for this case.
For better comprehensibility, the values in Table 2 and Table 3 are averaged over the four cases for each setup in Table 4. The mean error regarding the maximum settling velocity in the reference was given by 6.35 % . For the same grid spacing, the combinations EDM and MEA-W, EDM and MEA-L and GUO and MEA-W performed best (about 5.5 % ) with the first one yielding slightly better results. The same holds for the area between the curves, PCM and RMSE. While for smaller grid spacing the results were predominantly better, the same conclusions could be drawn as for N = 1 , with the combination of EDM and MEA-W performing better or equal to all other combinations. Only for PCM, the combination of EDM and MLA yielded the best results while it performed worst in all other cases. Looking only at the results for the cases with lower Reynolds numbers however, this setup yielded the highest accuracy in all values except in the prediction of the maximum settling velocity. This might also be due to the fact that the latter was not calculated from experimental values, but from a given correlation.
A possible reason for this behaviour is that forcing terms for points in the interior of the particle domain did not cancel out in contrast to MEA approaches. While for static objects no difference in the results was observed for the MLA, the error increased with the Reynolds number. This can also be seen in Figure 2. Since the objects did not have strict no-slip boundaries, simply omitting contributions to the force of points beyond the physical boundary was no solution. The transition region between the cells relevant and the ones not relevant for the hydrodynamic force within the particle domain seemed to be dependent on both the discretization parameters and the velocity. Further research in the future is required at this point. In empirical studies the influence of different parameters like grid spacing and the velocity or Reynolds number need to be quantified to deduce further steps in improving the proposed MLA.
Some configurations are plotted in Figure 2 for N = 8 , along with the experimental results by ten Cate et al. [45]. At this point, it can be concluded that all combinations tested in this study yielded reasonable results with the RMSE being the same order of magnitude as the RMSE for the utilized fit. While the application of EDM and MLA should be restricted to low Reynolds numbers only, where it performed best among the considered schemes. The combination of SCF and MEA-L, which was used in previous publications of HLBM [30,55] was among the worst, therefore the overall best performing configuration of EDM and MEA-W was chosen for the remaining calculations.
To further validate the implementation, a convergence study was performed. The experimental order of convergence (EOC) given by
E O C ( N , N ^ ) = log ( err ( N ) ) log ( err ( N ^ ) ) log ( N ) log ( N ^ ) ,
was calculated for two different grid spacings δ x / N and δ x / N ^ with N < N ^ . The function err : [ 1 , 2 , 4 , 8 ] R gives the difference of a chosen error or similarity measure for a given N to the value in the same setup but with a grid spacing according to N = 8 . This way, for each case and setup, two values, i.e., E O C ( 1 , 2 ) and E O C ( 2 , 4 ) were obtained. Averaging all values for the experimental order of convergence of a given combination of schemes for a measure leads to Table 5.
The mean EOC values were between 1 and 2 for all schemes, thereby superlinear convergence can be assumed. Similar to the results in previous tables, combinations of EDM or GUO with MEA-W or MEA-L performed best while combinations with SCF were falling behind. The MLA yielded the highest convergence rate for PCM and the maximum settling velocity while it also yielded the worst for the other measures. The results for Re = 11.6 with EDM and MEA-W are depicted for different grid spacings according to the parameter N in Figure 3.
The combination of EDM and MEA-W was chosen for all further computations, as it proved best in overall performance. The applicability of the MLA, however, seemed to be limited to some low Reynolds number cases.
Thereby the HLBM was updated, while in previous publications a combination of SCF and MEA-L was used for the calculations the new approach consisting of EDM and MEA-W yielded better results and a higher accuracy.

4.1.3. Simulation Setup—Comparison to Correlations

To check if the drag correlations presented in Section 2.1 can be reproduced by the HLBM, again the setup of a single settling sphere was chosen. It was placed with the center 0.021   m above the bottom of a domain of size 0.0055   m × 0.0055   m × 0.022   m . The top and bottom of the domain were equipped with no-slip boundaries again, but in contrast to the previous case the sides were chosen to be periodic. The densities of the fluid and the sphere were ρ f = 1000   kg   m 3 and ρ p = 2500   kg   m 3 , respectively. The dynamic viscosity was varied to change the Reynolds number, the sphere had a diameter of d p = 0.0005   m and the gravitational acceleration was given by g = 9.81   m   s 2 . In this subsection Re is defined via the given ρ f and d p , while u is defined by the used drag correlation and μ f as shown in Table 6.
The grid spacing was given by δ x = 3.125 × 10 5   m , this means the particle diameter was equivalent to 8 grid cells. As stated by ten Cate et al. [45], DNS methods experience a non-physical dependency of the drag force on the kinematic viscosity. This has also been investigated by Rohde et al. [82] for a boundary taking the actual, non-voxelized, surface. This effect is described to be less significant for Re 1 and a scaling procedure is proposed to estimate the hydrodynamic radius for the simulation. In this study, the radius and grid spacing were kept fixed, only the temporal discretization was adapted. LBM simulations required the maximum occurring velocity in lattice units to be much smaller than c s to ensure incompressible flow conditions and stability. Estimating the maximum velocity u by the drag correlation in Equation (8) and choosing the discretization such that this velocity was equivalent to 0.01 in lattice units finally led to
δ t = 0.01 u δ x .

4.1.4. Results and Discussion—Comparison to Correlations

In all simulations, the sphere reached a stable settling velocity, which was used to calculate the drag coefficient using Equation (4) with the fluid velocity assumed to be zero. The results are presented in Table 6 along with Reynolds numbers according to the terminal velocity of the drag correlations by Stokes, Abraham and Schiller and Naumann, discussed in Section 2.1. Additionally, the deviation in % from the analytical value of the respective correlation is given.
Here, the closest correlation was the one by Schiller and Naumann with an average error of 7.78 % calculated from the absolute values of the ones given in Table 6. The results yielded a slightly worse error of 10.75 % comparing to Abraham. Additionally, as expected for Re < 1 the results were closer to the prediction using Stokes drag.
On the other hand, deviations in this region might also be due to the stronger viscosity dependency for low Re discussed in Section 4.1.3. Another reason for deviations might be the lattice relaxation time entering the region of under-relaxation [23] for Re < 1 and approaching the value of 0.5 for the high Re considered here. For a stable simulation τ > 0.5 is required [23]. This, however, is related to the chosen discretization. For reasons of comparability between the data points δ x was kept constant, however, adapting it to achieve a reasonable τ can presumably improve the results. Increasing the grid resolution and also the gravitational acceleration, Reynolds numbers up to 1591 were reached in tests. At this point, a further increase was limited by available computational resources, as not only did δ x need to be lowered, but also the domain size had to be increased to allow the particle to reach its terminal velocity. Overall, the results are in good agreement with the discussed correlations, this can be seen in Figure 4, where the Reynolds number was defined via the maximum settling velocity in the simulation.
For further analysis, δ t was varied by choosing a different lattice velocity to be equivalent to the one calculated by Schiller and Naumann in Equation (27). These results can be used in further simulations for an a-priori estimate of the error regarding the temporal resolution, for cases that allow an estimation of the maximum occurring velocity.
The results in Table 7 show that u L should be at least smaller or equal to 0.04 as this gave a 7.5 % deviation from the finest tested temporal resolution. A value of u L = 0.01 would be most desirable since from this point on the deviation was smaller than most errors measured in this work. Then again, this meant the computational cost quadruple compared to the just found minimum requirement regarding the lattice velocity, since u L was directly proportional to δ t . Using u L = 0.0025 as a reference was justified as halving u L = 0.005 only led to a change below 1 % regarding the maximum settling velocity. Overall, almost perfect linear convergence was observed considering the average E O C = 0.97 regarding the maximum settling velocity. The results are also depicted in Figure 5, where a convergence towards a value in between the predictions by Abraham [60] and Schiller and Naumann [61] is shown.

4.1.5. Further Discussion—Onset of Unsteadiness

Spheres settling under the influence of gravity experience a range of regimes of motion, depending on the density ratio between the spheres and the fluid as well as the Reynolds or Galileo number. This has been shown and experimentally investigated by Horowitz and Williamson [83]. Later, Rahmani and Wachs [84] investigated the same behaviour via simulations. Following their definition of the Galileo number as
Ga = 1 ρ p / ρ f g d p 3 ν ,
the simulations presented in Section 4.1.4 cover a range of Ga = 2.14 up to Ga = 548.97 . For this discussion, the size of the domain has been extended to 0.0375 m in z-direction. In this context, it is notable that the unsteadiness observed in experiments can only be reproduced in simulations if the setup is not perfectly symmetrical, i.e., if the sphere does not perfectly align with the grid for a domain free of disturbance. Therefore, the sphere is moved 0.3 δ x from the center of the domain in x- and y-direction respectively. It can be observed that for Ga = 137.24 , a sphere follows a slightly oblique path, while all spheres with lower Ga settle vertically straight, as depicted in Figure 6. This fits the results by Jenny et al. [85], who found the onset of oblique motion at about Ga = 150 , and Rahmani and Wachs [84].
For higher Galileo numbers, the motion seems to be characterized by an oblique oscillation. However, to ensure this classification, a closer evaluation with a longer settling path is required, especially since Rahmani and Wachs [84] as well as Jenny et al. [85] predicted a chaotic motion for this parameters. The vortex structure shown in Figure 7 also seems to be rather chaotic. A thorough investigation of the different settling regimes including the vortex shedding, which was found as depicted in Figure 7, however, exceeds the scope of the present work at this point.

4.2. Tubular Pinch Effect

In this section, 2D simulations regarding the tubular pinch effect are presented. It was first discovered [46] and then further investigated [47] by Segré and Silberberg. The effect describes the motion of neutrally buoyant particles in a tube flow towards an equilibrium position between the tube’s center and its wall. As reference simulation, results by Li et al. [35] and Inamuro et al. [36] are considered.

4.2.1. Simulation Setup

For all considered cases, the base setup depicted in Figure 8 was the same, only the parameters change. The tube with diameter D and length L was equipped with pressure boundary conditions at front and back. A flow was induced by a pressure difference Δ p L , considering the relation between pressure and density this was equivalent to setting ρ L = 1 ± p L / 6 , for an initial density of ρ L = 1 across the domain. For the particle, a periodic boundary was chosen. As contact with the given boundary implementation distorted the particle movement, the particle was always kept one particle diameter d p away. While pressure boundaries were applied to the fluid domain, the particle experienced the periodic boundary before the actual end of the domain and got cut off one particle diameter before the pressure outlet. The part of the particle having left the domain via this boundary entered the domain also one particle diameter from the pressure inlet. This is done similarly by Li et al. [35] while Inamuro et al. [36] present a boundary condition capable of handling particles passing through it. To account for this difference, the domain was enlarged by 2 d p in flow direction so that the domain length stayed the same from perspective of the circle. The top and bottom of the domain had a no-slip condition for the fluid as well as the submersed object. The periodicity for the particles was required to reduce computational load since, as depicted in Figure 9, the domain would require between 15 and 60 times the size, depending on the case. It also allowed us to study the influence of the distance between spheres by varying the domain length as presented by Inamuro et al. [36]. The application of periodicity as described by Li et al. [35] was required as the improved periodic particle boundary presented in Inamuro et al. [36] was formulated for particles represented by a bounce-back boundary and could therefore not be applied straightforward for HLBM.
The calculations in the reference literature were performed in lattice units, therefore the discretization parameters were chosen accordingly. Additionally, all values in this subsection and the one regarding the results of this case are therefore non-dimensionalized by those discretization parameters, just like shown in the references. As no value was affiliated with a unit, the superscript L , indicating lattice units, is also dropped in this section for reasons of readability. The viscosity is defined by Equation (14) for a given τ . Since, as also stated by Inamuro et al. [36], the error was proportional to the maximum occuring lattice velocity squared, it was monitored for the presented cases.
The trajectory of the submersed circle is described by the coordinates of its center ( x c , y c ) . All simulations are run for 6 × 10 2 time steps to ensure convergence according the the results shown in Figure 9.

4.2.2. Results and Discussion

The results of the simulations are printed in Table 8 for different cases with the respective parameters. For each case, the simulations were run for different starting points y start / D { 0.2 , 0.25 , 0.35 , 0.4 , 0.45 } of the circle and the simulations showed a convergence towards a single position for each case as depicted in Figure 9. The resulting y c / D was averaged over the last 20 , 000 time steps of each simulation and over the results of all starting points.
Inamuro et al. [36] investigated the influence of various parameter like Re , d p / D and the distance between two circles, which is reflected in the area covered by particles. They found the same relations as Karnis et al. in their experiments [49], namely that the equilibrium position of the circle gets closer to the wall for higher Re , but closer to the center for more distance between two objects and also gets closer to the center for increasing d p / D . Inamuro et al. [36] defined the Reynolds number for this case as the ratio of time- and space-averaged velocity times the distance between the walls to the kinematic viscosity. The values were chosen to approximately achieve a maximum lattice velocity of u L = 0.04 , this was also found in the current simulations with u L 0.041 across all runs. With an average absolute error of 4.75 % across the cases the results showed good agreement to the ones obtained by Inamuro et al. [36]. The most probable source of the deviation was the difference in boundary conditions at the in- and outflow.
With an error of 4.1 % , the last case is also in good agreement with Li et al. [35], however, they also stated an error of 5.16 % comparing with case 2 of Inamuro et al. [36]. As the parameter did not fully match any of those reference cases, a comparison between the two reference publications is complicated, especially since for the last case in Table 8, a maximum lattice velocity of u L 0.096 was measured.
Additionally, the applicability of the MEA-L was tested for HLBM in this case. As already stated by Li et al. [35] and Peng et al. [34], it was found that the MEA according to Ladd [39,40] is not able to reproduce the tubular pinch effect. In all cases, the circle reached its final position in the middle of the tube.

4.3. Hindered Settling

In this section, the phenomenon of hindered settling is simulated with HLBM. Although there was no explicit collision model implemented, the particles affected each other by momentum transfer via the fluid. This, of course, was only sufficient if the spatial and temporal resolution was chosen fine enough that the transferred momentum prevented an overlapping of particles by several cells. This effect has been observed in previous publications regarding HLBM, e.g., by Krause et al. [55], but has up to now only been described qualitatively. Therefore, the study of hindered settling was used to investigate the occurring error regarding the velocity and solid volume fraction. This case aimed at testing the quality of results without collision model and if the effect could be reproduced at all. Therefore, the simulations were compared to various correlations presented in Section 2.2.

4.3.1. Simulation Setup

For this simulations, random distributions of spheres were created in the bottom half of a domain of size 0.00625 m × 0.00625 m × 0.025 m. The latter was equipped with no-slip boundary conditions at the bottom and top and was periodic for fluid and particles at the sides. The considered solid volume fractions in the bottom half are 5 % , 10 % , 15 % , 20 % and 25 % , which required 373, 746, 1119, 1492 and 1865 spheres, respectively. The starting setup for the solid volume fraction ϕ = 0.2 is depicted in Figure 10.
The spheres had a diameter of d p = 0.0005   m and a density of ρ p = 2500   k g   m 3 , while the physical parameters of the fluid were given by ρ f = 1000   k g   m 3 and ν = 10 6   m 2   s 1 . To study the effect regarding different Reynolds numbers the gravitational acceleration was varied. For an a-priori estimation of Re , the terminal settling velocity u of a sphere according to Schiller and Naumann was used (see Section 2.1). Investigated are the cases of Re = 0.53 ( g = 0.056   m   s 2 ), Re = 5.29 ( g = 0.748   m   s 2 ) and Re = 49.64 ( g = 15.22   m   s 2 ).
For spatial discretization δ x = 2.27 × 10 5   m was chosen, thereby the diameter of a sphere was equivalent to 11 grid cells. Furthermore
δ t = u u L δ x
was chosen to ensure a sufficient temporal resolution, with u L = 0.005 , for all cases.

4.3.2. Results and Discussion

The hindered settling velocity of the particle collective was calculated by tracking the upper front to the clear water zone. Therefore, the settling velocity in direction of gravitation was averaged over the highest 5 % of particles. The velocity field during the process is visualized in Figure 11.
Since the number of simulated particles was low compared to an experimental setup, especially for low solid volume fractions, the calculated velocity is prone to oscillations. Therefore, it was further averaged over time, beginning with the first point, for which the absolute of the averaged velocity decreases. The endpoint was defined by the the time the first particle of the considered 5 % reached the bottom 15 % of the domain. The results of the averaging over the particles for Re = 49.64 are depicted in Figure 12 along with the part used for temporal averaging, represented by the dashed line. In the following, this latter average will be denoted as simulated hindered settling velocity u sim .
While a comparison to Steinour [63] resulted in throughout large errors, the deviations to different other correlations based on a single particle settling velocity according to Stokes [59] and to Schiller and Naumann [61] (here denoted by u S N ) are given in Table 9. The best agreement was achieved with the results by Barnea and Mizrahi [70] in combination with the drag correlation by Schiller and Naumann with an average error of 8.07 % . This was expectable since they already stated to use drag correlations beyond Stokes with their formula for higher Reynolds numbers. For the case of Re = 0.53 the used drag correlation had negligible influence as it yielded 4.77 % deviation to the results with Stokes and 4.82 % to the results in combination with u S N . Except for the correlation of Oliver built upon Schiller and Naumann yielding an average error of 7.05 % , mainly the case of Re = 49.46 contributed to the error. This was to be expected since the discretization remained the same, but the spheres moved faster, increasing the possibility of large particle volumes overlapping and thereby making further momentum transfer via the fluid impossible. Additionally, high solid volume fractions seemed to have a negative effect on the results, this is shown in the last column of Table 9. Similar as for the high Reynolds number cases the error likely originated from the missing explicit contact model. For this comparably large number of spheres the occurrence of multiple contacts in a short time was possible, which pushed spheres to overlap each other. Despite these errors it was shown that hindered setting could be simulated with this approach however, with limitations to the solid volume fraction and velocity in combination with spatial and temporal discretization.
In this setup the only boundary conditions applied were periodic boundaries and a bounce-back no-slip condition, which are both mass conserving. Since the collision step executed a standard BGK collision with forcing, which was also mass conserving [23], this was also true for the whole scheme. Additionally, during the simulations no notable accumulation or dilution of mass within the areas covered by particles were observed.
Independent of the used drag formula, the correlation according to Richardson and Zaki [54] yielded the worst results among the more precisely studied ones. The average error was 47.27 % for the combination with u S , intended by the originators. However, it dropped to 17.78 % using u S N . The varying deviations to the drag correlations and obviously also between those correlations show the complexity of this application case.
Many reasons for errors and deviations can be found for this application case. Besides the possible influence of the starting distribution of the spheres, it can also be deduced from Figure 12 that a larger domain with more particles was required for more reliable results. Since the clear water front was not necessarily sharp and easy to track, a broader domain with more particles would smooth oscillations in this chaotic top region. Lowering the percentage of tracked particles e.g., to 2 % , improved the results in comparison to the correlations, however, it also increased the oscillations. Using the top 5 % was a compromise between reliable results and a precise tracking of the front by only considering the uppermost particles. A broader domain and clear water front would therefore be beneficial. The extend of the simulation, however, was limited for DNS approaches due to high computational costs.
Furthermore, the probable main point is that no explicit particle-particle collision model is applied and the momentum is solely exchanged via the fluid. This, of course, can be a main source of error, especially if the particles start to overlap by several cells. The latter is caused by a resolution chosen too coarse in comparison to the velocity of the spheres, so that the repelling effect of momentum exchange via the fluid is diminished. The good agreement, e.g., with the results by Barnea and Mizrahi and the obvious dependency of the simulated hindered settling velocity on ϕ however, show that it is possible to sufficiently depict the effect of hindered settling with the chosen setup. While capable of properly simulating the settling phase, the presented method is not suitable for the study of bed formation, as particles start to overlap once sedimented due to the absence of an explicit collision model. No influence of this nonphysical behaviour back on the settling process was observed.

4.4. Computational Efficiency

The performance of the algorithm has been recently investigated by Bretl et al. [86]. In their study, also based on OpenLB [56], they investigated the speedup S, which was defined as ratio of execution times using a single process T 1 to using n p processes. Comparing to the theoretical value for a code which is 99.5 % parallelizable given by
S = T 1 T 1 ( 0.005 + 0.995 / n p ) ,
they found the results to be in good agreement. These, however, are results for the algorithm in the case of a single settling particle. Considering the case of hindered settling presented in Section 4.3, the dependency of the performance on the number of simulated objects was studied. As measure, the amount of million lattice site updates per second and core (MLUPps) was used. The cases for Re = 49.64 were evaluated with the simulations run on a system equipped with Intel Xeon E5-2660 v3 processors.
Each simulation used a lattice with 83,717,424 cells and took 176,000 time-steps utilizing 80 cores. The results are depicted in Figure 13 together with a fit to a rational function f with
f ( ϕ ) = a 1 ϕ + a 2 a 3 ϕ + a 4 ,
and fitting parameters a 1 , a 2 , a 3 , a 4 R . The function yielded 0.45 MLUPps for ϕ = 1 and approached 0.35 MLUPps asymptotically as ϕ . Despite the questionable physicality regarding such a solid volume fraction for a set of spherical particles, those numbers were useful to estimate the lower bounds of single core performance for this method. The decrease in MLUPps was to be expected, since for each cell covered by a particle, additional computations regarding the hydrodynamic force were required (see e.g., Equation (23)). The influence of the number of particles on the performance, however, did not only depend on the solid volume fraction, but also on the distribution of the objects. Since for the parallelization MPI was used, each particle was treated by the process responsible for the domain it resided in. As in this simulation the particles accumulated at the bottom of the domain, a small number of processors had to take care of the calculation of hydrodynamic forces for all particles. A discussion on distribution strategies of computational load among the processors (load optimal strategy vs. communication optimal strategy) is given by Henn et al. [87].

5. Conclusions

The homogenized lattice Boltzmann method has been revisited and various forcing schemes and approaches to calculate exchanged momentum have been evaluated. Among the latter is a new proposed algorithm based on the momentum loss on a cell. However, it showed substandard accuracy except in the case of low Reynolds numbers. It is assumed that the hydrodynamic force is overestimated due to effects correlated to the interior fluid. The amount of force contributed by the latter is hardly distinguishable from the amount related to momentum exchanged with the bulk fluid. As most MEAs are constructed for a sharp solid boundary, further research at this point should be considered in the future. Since issues related to the interior fluid also exist for the immersed boundary method, it might be possible to adapt improvements.
Besides this a combination of Kupershtokh forcing and the MEA by Wen et al. [38] gave the overall best results. This finally leads to an updated version of HLBM, for which the results are in good agreement to literature (e.g., regarding the velocity profile of a settling sphere).
Furthermore, the reproducibility of drag correlations was tested yielding an error of 7.78 % compared to the one by Schiller and Naumann in the range from Re = 0.24 up to Re = 948.67 . Therefore, the applicability of HLBM for Reynolds numbers up to the Newton regime is shown.
In the next application case, the tubular pinch effect was investigated in 2 dimensions. Here, the error to reference simulations in the literature regarding the influence of the Reynolds number, the ratio of circle to tube diameter and the distance between the circles, never exceeded 5.83 % .
The authors furthermore showed that the cases of hindered settling can be investigated with HLBM for a sufficient resolution, even without an explicit collision model. With a deviation of 8.07 % to the correlation by Barnea and Mizrahi, the results are in excellent agreement, especially considering the error regarding a single settling particle. For more reliable computations regarding hindered settling, particularly for higher Reynolds numbers, and for bed formation processes an additional collision model is mandatory. The implementation of such is planned for future research, especially to evaluate the actual improvement in the quality of results as an effect of the addition of an explicit contact model.
Overall, the updated version of HLBM performed well across all application cases with good agreement to existing literature regarding experimental and theoretical results.

Author Contributions

Conceptualization, R.T.; methodology, R.T.; software, R.T., N.H., T.W., M.J.K.; validation, R.T.; formal analysis, R.T., T.W.; investigation, R.T.; resources, M.J.K., H.N.; data curation, R.T.; writing—original draft preparation, R.T.; writing—review and editing, R.T., N.H., G.T., M.J.K.; visualization, R.T.; supervision, M.J.K., H.N.; project administration, M.J.K., R.T.; funding acquisition, M.J.K., R.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Deutsche Forschungsgemeinschaft (DFG) [grant KR 4259/8-2] within the priority program SPP2045 MehrDimPart “Highly specific and multidimensional fractionation of fine particle systems with technical relevance”. The authors also acknowledge support by the Open Access Publishing Fund of the Karlsruhe Institute of Technology (KIT).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Acknowledgments

The authors gratefully acknowledge the Steinbuch Centre for Computing at KIT for providing access to their high performance computer ForHLR II, where most computations were carried out as part of the CPE project. The authors also acknowledge DUG Technology for providing high-performance computing resources and expertise to support this research.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
BGKBhatnagar–Gross–Krook
DNSdirect numerical simulation
EDMexact difference method
EOCexperimental order of convergence
GUOGuo forcing
HLBMhomogenised lattice Boltzmann method
LBMlattice Boltzmann method
MEAmomentum exchange algorithm
MEA-Lmomentum exchange algorithm according to Ladd [39,40]
MEA-Wmomentum exchange algorithm according to Wen et al. [42]
MLAmomentum loss algorithm
MLUPpsmillion lattice site updates per second and per processor
PCMpartial curve mapping
RMSEroot mean squared error
SCFShan–Chen forcing scheme
Roman
Aprojected area of an object in flow direction
Bdomain covered by a particle
C D drag coefficient
c i i-th discrete lattice velocity
c s lattice speed of sound
d, d ^ spatial dimension
d B mapping function of an object on the lattice
d p particle diameter
Dtube diameter
f i particle distribution function in the phase space according to the i-th lattice velocity
f i eq Maxwell–Boltzmann distribution function according to the i-th lattice velocity
F force
F f force acting on the fluid
F p force acting on the particle
F B Buoyancy force
F D drag force
F G gravitational force
F H hydrodynamic force
ggravitational acceleration
Itime interval
I h discrete time interval
J p moment of inertia
Ltube length
m p particle mass
N, N ^ resolution parameter
n p number of processes
ppressure
p L pressure in lattice units
qdimension of the velocity space of a lattice
r distance to the center of mass of an object
Re Reynolds number
Sspeedup
S i source term in the LBM collision step according to the i-th lattice velocity
ttime
T p torque
u velocity
u ˜ redefined velocity in the presence of a particle
u BM hindered settling velocity according to Barnea and Mizrahi
u eq velocity used in the Maxwell–Boltzmann distribution
u f fluid velocity
u O hindered settling velocity according to Oliver
u p particle velocity
u RZ hindered settling velocity according to Richardson and Zaki
u S terminal settling velocity according to Stokes
u S N terminal settling velocity according to Schiller and Naumann
u sim simulated hindered settling velocity
u Steinour hindered settling velocity according to Steinour
u maximum settling velocity
u L velocity in lattice units
w i weighting function according to the i-th lattice velocity
x coordinates of a lattice point
x c x-coordinate of particle center
y c y-coordinate of particle center
y start initial y-coordinate of particle center
Greek
δ t time step size in SI units
δ x grid spacing in SI units
δ t L time step size in lattice units
δ x L grid spacing in lattice units
Δ u difference between fluid velocity and the velocity in presence of an object
μ f dynamic viscosity
ν kinematic viscosity
ρ density in lattice units
ρ f fluid density
ρ p particle density
τ lattice relaxation time
ϕ solid volume fraction
ω angular velocity
Ω spatial domain
Ω h discrete approximation of the spatial domain Ω

References

  1. Viduka, S.M.; Feng, Y.Q.; Hapgood, K.; Schwarz, M. Discrete particle simulation of solid separation in a jigging device. Int. J. Miner. Process. 2013, 123, 108–119. [Google Scholar] [CrossRef]
  2. Li, J.; Webb, C.; Pandiella, S.; Campbell, G. A Numerical Simulation of Separation of Crop Seeds by Screening—Effect of Particle Bed Depth. Food Bioprod. Process. 2002, 80, 109–117. [Google Scholar] [CrossRef] [Green Version]
  3. Champion, J.A.; Katare, Y.K.; Mitragotri, S. Particle shape: A new design parameter for micro- and nanoscale drug delivery carriers. J. Control. Release 2007, 121, 3–9. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Davies, R. A simple feature-space representation of particle shape. Powder Technol. 1975, 12, 111–124. [Google Scholar] [CrossRef]
  5. Scotti, R.; Wahba, L.; Crippa, M.; D’Arienzo, M.; Donetti, R.; Santo, N.; Morazzoni, F. Rubber-silica nanocomposites obtained by in situ sol–gel method: Particle shape influence on the filler-filler and filler-rubber interactions. Soft Matter 2012, 8, 2131. [Google Scholar] [CrossRef]
  6. Trunk, R.; Henn, T.; Dörfler, W.; Nirschl, H.; Krause, M.J. Inertial dilute particulate fluid flow simulations with an Euler–Euler lattice Boltzmann method. J. Comput. Sci. 2016, 17, 438–445. [Google Scholar] [CrossRef]
  7. Höcker, S.B.; Trunk, R.; Dörfler, W.; Krause, M.J. Towards the simulations of inertial dense particulate flows with a volume-averaged lattice Boltzmann method. Comput. Fluids 2018, 166, 152–162. [Google Scholar] [CrossRef]
  8. Li, Z.; Yazdani, A.; Tartakovsky, A.; Karniadakis, G.E. Transport dissipative particle dynamics model for mesoscopic advection-diffusion-reaction problems. J. Chem. Phys. 2015, 143, 014101. [Google Scholar] [CrossRef] [Green Version]
  9. Chu, K.; Wang, B.; Yu, A.; Vince, A. CFD-DEM modelling of multiphase flow in dense medium cyclones. Powder Technol. 2009, 193, 235–247. [Google Scholar] [CrossRef]
  10. Maier, M.L.; Milles, S.; Schuhmann, S.; Guthausen, G.; Nirschl, H.; Krause, M.J. Fluid flow simulations verified by measurements to investigate adsorption processes in a static mixer. Comput. Math. Appl. 2018, 76, 2744–2757. [Google Scholar] [CrossRef]
  11. Parteli, E.J.R.; Pöschel, T. Particle-based simulation of powder application in additive manufacturing. Powder Technol. 2016, 288, 96–102. [Google Scholar] [CrossRef]
  12. Dapelo, D.; Trunk, R.; Krause, M.J.; Bridgeman, J. Towards Lattice-Boltzmann modelling of unconfined gas mixing in anaerobic digestion. Comput. Fluids 2019, 180, 11–21. [Google Scholar] [CrossRef]
  13. Zhu, H.; Zhou, Z.; Yang, R.; Yu, A. Discrete particle simulation of particulate systems: A review of major applications and findings. Chem. Eng. Sci. 2008, 63, 5728–5770. [Google Scholar] [CrossRef]
  14. Shardt, O.; Derksen, J. Direct simulations of dense suspensions of non-spherical particles. Int. J. Multiph. Flow 2012, 47, 25–36. [Google Scholar] [CrossRef]
  15. Kravets, B.; Rosemann, T.; Reinecke, S.; Kruggel-Emden, H. A new drag force and heat transfer correlation derived from direct numerical LBM-simulations of flown through particle packings. Powder Technol. 2019, 345, 438–456. [Google Scholar] [CrossRef]
  16. Dolanský, J.; Chára, Z.; Vlasák, P.; Kysela, B. Lattice Boltzmann method used to simulate particle motion in a conduit. J. Hydrol. Hydromechanics 2017, 65, 105–113. [Google Scholar] [CrossRef] [Green Version]
  17. Wachs, A.; Girolami, L.; Vinay, G.; Ferrer, G. Grains3D, a flexible DEM approach for particles of arbitrary convex shape—Part I: Numerical model and validations. Powder Technol. 2012, 224, 374–389. [Google Scholar] [CrossRef]
  18. Rakotonirina, A.D.; Wachs, A. Grains3D, a flexible DEM approach for particles of arbitrary convex shape—Part II: Parallel implementation and scalable performance. Powder Technol. 2018, 324, 18–35. [Google Scholar] [CrossRef]
  19. Rakotonirina, A.D.; Delenne, J.Y.; Radjai, F.; Wachs, A. Grains3D, a flexible DEM approach for particles of arbitrary convex shape—Part III: Extension to non-convex particles modelled as glued convex particles. Comput. Part. Mech. 2019, 6, 55–84. [Google Scholar] [CrossRef] [Green Version]
  20. Nolan, G.; Kavanagh, P. Random packing of nonspherical particles. Powder Technol. 1995, 84, 199–205. [Google Scholar] [CrossRef]
  21. Uhlmann, M. An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 2005, 209, 448–476. [Google Scholar] [CrossRef] [Green Version]
  22. Tian, F.B.; Luo, H.; Zhu, L.; Liao, J.C.; Lu, X.Y. An efficient immersed boundary-lattice Boltzmann method for the hydrodynamic interaction of elastic filaments. J. Comput. Phys. 2011, 230, 7266–7283. [Google Scholar] [CrossRef] [Green Version]
  23. Krüger, T.; Kusumaatmaja, H.; Kuzmin, A.; Shardt, O.; Silva, G.; Viggen, E.M. The Lattice Boltzmann Method; Graduate Texts in Physics; Springer International Publishing: Cham, Switzerland, 2017. [Google Scholar]
  24. Feng, Z.G.; Michaelides, E.E. The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems. J. Comput. Phys. 2004, 195, 602–628. [Google Scholar] [CrossRef]
  25. Wu, J.; Wu, J.; Zhan, J.; Zhao, N.; Wang, T. A Robust Immersed Boundary-Lattice Boltzmann Method for Simulation of Fluid-Structure Interaction Problems. Commun. Comput. Phys. 2016, 20, 156–178. [Google Scholar] [CrossRef] [Green Version]
  26. Wu, J.; Shu, C. An improved immersed boundary-lattice Boltzmann method for simulating three-dimensional incompressible flows. J. Comput. Phys. 2010, 229, 5022–5042. [Google Scholar] [CrossRef]
  27. Owen, D.R.J.; Leonardi, C.R.; Feng, Y.T. An efficient framework for fluid-structure interaction using the lattice Boltzmann method and immersed moving boundaries. Int. J. Numer. Methods Eng. 2011, 87, 66–95. [Google Scholar] [CrossRef]
  28. Beny, J.; Latt, J. Efficient LBM on GPUs for dense moving objects using immersed boundary condition. arXiv 2019, arXiv:1904.02108. [Google Scholar]
  29. Kang, S.K.; Hassan, Y.A. A comparative study of direct-forcing immersed boundary-lattice Boltzmann methods for stationary complex boundaries. Int. J. Numer. Methods Fluids 2011, 66, 1132–1158. [Google Scholar] [CrossRef]
  30. Trunk, R.; Marquardt, J.; Thäter, G.; Nirschl, H.; Krause, M.J. Towards the simulation of arbitrarily shaped 3D particles using a homogenised lattice Boltzmann method. Comput. Fluids 2018, 172, 621–631. [Google Scholar] [CrossRef]
  31. Kuznik, F.; Obrecht, C.; Rusaouen, G.; Roux, J.J. LBM based flow simulation using GPU computing processor. Comput. Math. Appl. 2010, 59, 2380–2392. [Google Scholar] [CrossRef] [Green Version]
  32. Feichtinger, C.; Habich, J.; Köstler, H.; Rüde, U.; Aoki, T. Performance modeling and analysis of heterogeneous lattice Boltzmann simulations on CPU–GPU clusters. Parallel Comput. 2015, 46, 1–13. [Google Scholar] [CrossRef]
  33. Noble, D.R.; Torczynski, J.R. A Lattice-Boltzmann Method for Partially Saturated Computational Cells. Int. J. Mod. Phys. C 1998, 9, 1189–1201. [Google Scholar] [CrossRef]
  34. Peng, C.; Teng, Y.; Hwang, B.; Guo, Z.; Wang, L.P. Implementation issues and benchmarking of lattice Boltzmann method for moving rigid particle simulations in a viscous flow. Comput. Math. Appl. 2016, 72, 349–374. [Google Scholar] [CrossRef] [Green Version]
  35. Li, H.; Lu, X.; Fang, H.; Qian, Y. Force evaluations in lattice Boltzmann simulations with moving boundaries in two dimensions. Phys. Rev. E 2004, 70. [Google Scholar] [CrossRef]
  36. Inamuro, T.; Maeba, K.; Ogino, F. Flow between parallel walls containing the lines of neutrally buoyant circular cylinders. Int. J. Multiph. Flow 2000, 26, 1981–2004. [Google Scholar] [CrossRef]
  37. Mei, R.; Yu, D.; Shyy, W.; Luo, L.S. Force evaluation in the lattice Boltzmann method involving curved geometry. Phys. Rev. E 2002, 65. [Google Scholar] [CrossRef] [Green Version]
  38. Wen, B.; Li, H.; Zhang, C.; Fang, H. Lattice-type-dependent momentum-exchange method for moving boundaries. Phys. Rev. E 2012, 85. [Google Scholar] [CrossRef]
  39. Ladd, A.J.C. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation. J. Fluid Mech. 1994, 271, 285–309. [Google Scholar] [CrossRef] [Green Version]
  40. Ladd, A.J.C. Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results. J. Fluid Mech. 1994, 271, 311–339. [Google Scholar] [CrossRef] [Green Version]
  41. Clausen, J.R.; Aidun, C.K. Galilean invariance in the lattice-Boltzmann method and its effect on the calculation of rheological properties in suspensions. Int. J. Multiph. Flow 2009, 35, 307–311. [Google Scholar] [CrossRef]
  42. Wen, B.; Zhang, C.; Tu, Y.; Wang, C.; Fang, H. Galilean invariant fluid-solid interfacial dynamics in lattice Boltzmann simulations. J. Comput. Phys. 2014, 266, 161–170. [Google Scholar] [CrossRef] [Green Version]
  43. Chen, Y.; Cai, Q.; Xia, Z.; Wang, M.; Chen, S. Momentum-exchange method in lattice Boltzmann simulations of particle-fluid interactions. Phys. Rev. E 2013, 88. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Lorenz, E.; Caiazzo, A.; Hoekstra, A.G. Corrected momentum exchange method for lattice Boltzmann simulations of suspension flow. Phys. Rev. E 2009, 79. [Google Scholar] [CrossRef] [PubMed]
  45. Ten Cate, A.; Nieuwstad, C.H.; Derksen, J.J.; Van den Akker, H.E.A. Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity. Phys. Fluids 2002, 14, 4012–4025. [Google Scholar] [CrossRef]
  46. Segré, G.; Silberberg, A. Radial Particle Displacements in Poiseuille Flow of Suspensions. Nature 1961, 189, 209–210. [Google Scholar] [CrossRef]
  47. Segré, G.; Silberberg, A. Behaviour of macroscopic rigid spheres in Poiseuille flow Part 2. Experimental results and interpretation. J. Fluid Mech. 1962, 14, 136–157. [Google Scholar] [CrossRef]
  48. Tachibana, M. On the behaviour of a sphere in the laminar tube flows. Rheol. Acta 1973, 12, 58–69. [Google Scholar] [CrossRef]
  49. Karnis, A.; Goldsmith, H.L.; Mason, S.G. The flow of suspensions through tubes: V. Inertial effects. Can. J. Chem. Eng. 1966, 44, 181–193. [Google Scholar] [CrossRef]
  50. Baldock, T.; Tomkins, M.; Nielsen, P.; Hughes, M. Settling velocity of sediments at high concentrations. Coast. Eng. 2004, 51, 91–100. [Google Scholar] [CrossRef]
  51. Derksen, J.J. Eulerian-Lagrangian simulations of settling and agitated dense solid-liquid suspensions-achieving grid convergence. AIChE J. 2018, 64, 1147–1158. [Google Scholar] [CrossRef]
  52. Deshpande, R.; Antonyuk, S.; Iliev, O. DEM-CFD study of the filter cake formation process due to non-spherical particles. Particuology 2020. [Google Scholar] [CrossRef]
  53. Zaidi, A.A.; Tsuji, T.; Tanaka, T. Hindered Settling Velocity & Structure Formation during Particle Settling by Direct Numerical Simulation. Procedia Eng. 2015, 102, 1656–1666. [Google Scholar] [CrossRef] [Green Version]
  54. Richardson, J.; Zaki, W. The sedimentation of a suspension of uniform spheres under conditions of viscous flow. Chem. Eng. Sci. 1954, 3, 65–73. [Google Scholar] [CrossRef]
  55. Krause, M.J.; Klemens, F.; Henn, T.; Trunk, R.; Nirschl, H. Particle flow simulations with homogenised lattice Boltzmann methods. Particuology 2017, 34, 1–13. [Google Scholar] [CrossRef]
  56. Krause, M.J.; Kummerländer, A.; Avis, S.J.; Kusumaatmaja, H.; Dapelo, D.; Klemens, F.; Gaedtke, M.; Hafen, N.; Mink, A.; Trunk, R.; et al. OpenLB—Open source lattice Boltzmann code. Comput. Math. Appl. 2020. [Google Scholar] [CrossRef]
  57. Hölzer, A.; Sommerfeld, M. New simple correlation formula for the drag coefficient of non-spherical particles. Powder Technol. 2008, 184, 361–365. [Google Scholar] [CrossRef]
  58. Bagheri, G.; Bonadonna, C. On the drag of freely falling non-spherical particles. Powder Technol. 2016, 301, 526–544. [Google Scholar] [CrossRef] [Green Version]
  59. Stokes, G.G. On the Effect of Internal Friction of Fluids on the Motion of Pendulums. Trans. Camb. Philos. Soc. 1851, 9, 8–106. [Google Scholar] [CrossRef] [Green Version]
  60. Abraham, F.F. Functional Dependence of Drag Coefficient of a Sphere on Reynolds Number. Phys. Fluids 1970, 13, 2194. [Google Scholar] [CrossRef]
  61. Schiller, L.; Naumann, A. Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung. Z. Des Vereines Dtsch. Ingenieure 1933, 77, 318–320. [Google Scholar]
  62. Dey, S.; Ali, S.Z.; Padhi, E. Terminal fall velocity: The legacy of Stokes from the perspective of fluvial hydraulics. Proc. R. Soc. Math. Phys. Eng. Sci. 2019, 475, 20190277. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  63. Steinour, H.H. Rate of sedimentation.Nonflocculated Suspensions of Uniform Spheres. Ind. Eng. Chem. 1944, 36, 618–624. [Google Scholar] [CrossRef]
  64. Oliver, D. The sedimentation of suspensions of closely-sized spherical particles. Chem. Eng. Sci. 1961, 15, 230–242. [Google Scholar] [CrossRef]
  65. Richardson, J.; Zaki, W. Sedimentation and fluidisation: Part I. Chem. Eng. Res. Des. 1954, 75, S82–S100. [Google Scholar] [CrossRef]
  66. Rowe, P.N. A convenient empirical equation for estimation of the Richardson-Zaki exponent. Chem. Eng. Sci. 1987, 42, 2795–2796. [Google Scholar] [CrossRef]
  67. Di Felice, R.; Gibilaro, L.G.; Foscolo, P.U. On the hindered settling velocity of spheres in the inertial flow regime. Chem. Eng. Sci. 1995, 50, 3005–3006. [Google Scholar] [CrossRef]
  68. Di Felice, R. The sedimentation velocity of dilute suspensions of nearly monosized spheres. Int. J. Multiph. Flow 1999, 25, 559–574. [Google Scholar] [CrossRef]
  69. Lu, Y.; Wei, L.; Wei, J. A numerical study of bed expansion in supercritical water fluidized bed with a non-spherical particle drag model. Chem. Eng. Res. Des. 2015, 104, 164–173. [Google Scholar] [CrossRef]
  70. Barnea, E.; Mizrahi, J. A generalized approach to the fluid dynamics of particulate systems. Chem. Eng. J. 1973, 5, 171–189. [Google Scholar] [CrossRef]
  71. Bhatnagar, P.L.; Gross, E.P.; Krook, M. A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems. Phys. Rev. 1954, 94, 511–525. [Google Scholar] [CrossRef]
  72. Verlet, L. Computer ’Experiments’ on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98–103. [Google Scholar] [CrossRef]
  73. Swope, W.C.; Andersen, H.C.; Berens, P.H.; Wilson, K.R. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J. Chem. Phys. 1982, 76, 637–649. [Google Scholar] [CrossRef]
  74. Shan, X.; Chen, H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 1993, 47, 1815–1819. [Google Scholar] [CrossRef] [Green Version]
  75. Guo, Z.; Zheng, C.; Shi, B. Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys. Rev. E 2002, 65. [Google Scholar] [CrossRef] [PubMed]
  76. Kupershtokh, A.; Medvedev, D.; Karpov, D. On equations of state in a lattice Boltzmann method. Comput. Math. Appl. 2009, 58, 965–974. [Google Scholar] [CrossRef] [Green Version]
  77. Huang, H.; Krafczyk, M.; Lu, X. Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models. Phys. Rev. E 2011, 84. [Google Scholar] [CrossRef] [Green Version]
  78. Caiazzo, A.; Junk, M. Boundary forces in lattice Boltzmann: Analysis of Momentum Exchange algorithm. Comput. Math. Appl. 2008, 55, 1415–1423. [Google Scholar] [CrossRef]
  79. Sukop, M.C.; Thorne, D.T. Lattice Boltzmann Modeling; Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  80. Jekel, C.F.; Venter, G.; Venter, M.P.; Stander, N.; Haftka, R.T. Similarity measures for identifying material parameters from hysteresis loops using inverse analysis. Int. J. Mater. Form. 2019, 12, 355–378. [Google Scholar] [CrossRef]
  81. Witowski, K.; Stander, N. Parameter Identification of Hysteretic Models Using Partial Curve Mapping. In Proceedings of the 12th AIAA Aviation Technology, Integration, and Operations (ATIO) Conference and 14th AIAA/ISSMO Multidisciplinary Analysis and Optimization Conference, Indianapolis, IN, USA, 17–19 September 2012. [Google Scholar]
  82. Rohde, M.; Derksen, J.J.; Van den Akker, H.E.A. Volumetric method for calculating the flow around moving objects in lattice-Boltzmann schemes. Phys. Rev. E 2002, 65. [Google Scholar] [CrossRef] [Green Version]
  83. Horowitz, M.; Williamson, C.H.K. The effect of Reynolds number on the dynamics and wakes of freely rising and falling spheres. J. Fluid Mech. 2010, 651, 251–294. [Google Scholar] [CrossRef]
  84. Rahmani, M.; Wachs, A. Free falling and rising of spherical and angular particles. Phys. Fluids 2014, 26, 083301. [Google Scholar] [CrossRef]
  85. Jenny, M.; Dušek, J.; Bouchet, G. Instabilities and transition of a sphere falling or ascending freely in a Newtonian fluid. J. Fluid Mech. 2004, 508, 201–239. [Google Scholar] [CrossRef]
  86. Bretl, C.; Trunk, R.; Nirschl, H.; Thäter, G.; Dorn, M.; Krause, M.J. Preliminary study of particle settling behaviour by shape parameters via lattice Boltzmann simulations. In High Performance Computing in Science and Engineering ’20; Nagel, W., Kröner, D., Resch, M., Eds.; Springer: New York, NY, USA, 2020. [Google Scholar]
  87. Henn, T.; Thäter, G.; Dörfler, W.; Nirschl, H.; Krause, M.J. Parallel dilute particulate flow simulations in the human nasal cavity. Comput. Fluids 2016, 124, 197–207. [Google Scholar] [CrossRef]
Figure 1. Setup of the simulation according to ten Cate at al. [45].
Figure 1. Setup of the simulation according to ten Cate at al. [45].
Computation 09 00011 g001
Figure 2. Settling velocity in the four defined cases for different combinations of forcing and momentum exchange schemes along with the experimental results by ten Cate et al. [45].
Figure 2. Settling velocity in the four defined cases for different combinations of forcing and momentum exchange schemes along with the experimental results by ten Cate et al. [45].
Computation 09 00011 g002
Figure 3. Extract of the results for the case of Re = 11.6 with exact difference method (EDM) and momentum exchange algorithm (MEA)-W for different grid spacings along with the experimental data by ten Cate et al. [45].
Figure 3. Extract of the results for the case of Re = 11.6 with exact difference method (EDM) and momentum exchange algorithm (MEA)-W for different grid spacings along with the experimental data by ten Cate et al. [45].
Computation 09 00011 g003
Figure 4. Drag coefficient of correlations discussed in Section 2.1 along with the simulation results plotted against the Reynolds number, computed using the maximum settling velocity measured in the simulation.
Figure 4. Drag coefficient of correlations discussed in Section 2.1 along with the simulation results plotted against the Reynolds number, computed using the maximum settling velocity measured in the simulation.
Computation 09 00011 g004
Figure 5. Maximum settling velocity for μ f = 0.005   k g   m 1   s 1 over the chosen lattice velocity u L (see Equation (27)).
Figure 5. Maximum settling velocity for μ f = 0.005   k g   m 1   s 1 over the chosen lattice velocity u L (see Equation (27)).
Computation 09 00011 g005
Figure 6. Path in the x-z-plane for spheres with different Ga .
Figure 6. Path in the x-z-plane for spheres with different Ga .
Computation 09 00011 g006
Figure 7. Contours of the stream-wise vorticity for a sphere with Ga = 548.97 at t = 0.19 s.
Figure 7. Contours of the stream-wise vorticity for a sphere with Ga = 548.97 at t = 0.19 s.
Computation 09 00011 g007
Figure 8. Computational domain for the investigation of the tubular pinch effect.
Figure 8. Computational domain for the investigation of the tubular pinch effect.
Computation 09 00011 g008
Figure 9. Results of various starting positions of the circle according to case 2 of Inamuro et al. [36].
Figure 9. Results of various starting positions of the circle according to case 2 of Inamuro et al. [36].
Computation 09 00011 g009
Figure 10. Random distribution of spheres in the computational domain for the hindered settling simulations with a solid volume fraction of 20 % .
Figure 10. Random distribution of spheres in the computational domain for the hindered settling simulations with a solid volume fraction of 20 % .
Computation 09 00011 g010
Figure 11. Velocity field around the sphere at t = 0.23   s for the case of Re = 5.29 with a solid volume fraction of ϕ = 0.2 .
Figure 11. Velocity field around the sphere at t = 0.23   s for the case of Re = 5.29 with a solid volume fraction of ϕ = 0.2 .
Computation 09 00011 g011
Figure 12. Averages over the settling velocity of the top 5 % of particles for various ϕ and Re = 49.64 . The part used for temporal averaging is depicted as dashed line.
Figure 12. Averages over the settling velocity of the top 5 % of particles for various ϕ and Re = 49.64 . The part used for temporal averaging is depicted as dashed line.
Computation 09 00011 g012
Figure 13. Performance data as million lattice updates per second and core, regarding the number of simulated particles.
Figure 13. Performance data as million lattice updates per second and core, regarding the number of simulated particles.
Computation 09 00011 g013
Table 1. Values of the setup of the four cases for a settling sphere investigated by ten Cate et al. [45] along with the deviation of their results to the maximum settling velocity u predicted with the drag correlation by Abraham [60].
Table 1. Values of the setup of the four cases for a settling sphere investigated by ten Cate et al. [45] along with the deviation of their results to the maximum settling velocity u predicted with the drag correlation by Abraham [60].
Case ρ f in kg m 3 μ f in kg m 1 s 1 Re δ x in m δ t in s u in m s 1 (Abraham [60])Error in %
1970 0.373 1.5 1.671 × 10 3 3.891 × 10 4 0.038 10.6
2965 0.212 4.1 1.645 × 10 3 2.41 × 10 4 0.06 5.0
3962 0.113 11.6 1.61 × 10 3 1.526 × 10 4 0.091 4.5
4960 0.058 31.9 1.559 × 10 3 1.01 × 10 4 0.128 5.3
Table 2. Different error measures [80] for the four cases described by ten Cate et al. [45] for various combinations of momentum exchange and forcing schemes. The error in terminal velocity is given relative to the analytical values according to Abraham [60]. The results in the table refer to N = 1 .
Table 2. Different error measures [80] for the four cases described by ten Cate et al. [45] for various combinations of momentum exchange and forcing schemes. The error in terminal velocity is given relative to the analytical values according to Abraham [60]. The results in the table refer to N = 1 .
CaseForcingMomentum ExchangeRMSEPCMArea between CurvesError in % Maximum Velocity
1EDMMEA-W 0.04191 3.63347 0.00386 4.69880
1EDMMLA 0.02722 2.17069 0.00292 4.91581
1EDMMEA-L 0.04187 3.61223 0.00386 4.71524
1GUOMEA-W 0.04190 3.63317 0.00386 4.69881
1SCFMEA-W 0.02616 3.86421 0.00323 6.73646
1SCFMEA-L 0.02618 3.87909 0.00323 6.73834
2EDMMEA-W 0.07019 1.96464 0.00703 6.08454
2EDMMLA 0.11885 4.41254 0.01115 6.71691
2EDMMEA-L 0.07028 1.96609 0.00704 6.08627
2GUOMEA-W 0.07019 1.96473 0.00703 6.08455
2SCFMEA-W 0.10945 3.48661 0.01082 7.83986
2SCFMEA-L 0.10952 3.49102 0.01082 7.85094
3EDMMEA-W 0.04287 1.50668 0.00409 6.23779
3EDMMLA 0.08819 2.58684 0.01000 8.01712
3EDMMEA-L 0.04289 1.50718 0.00409 6.23545
3GUOMEA-W 0.04287 1.50659 0.00409 6.23780
3SCFMEA-W 0.05442 2.00959 0.00585 7.27864
3SCFMEA-L 0.05443 2.00904 0.00585 7.28794
4EDMMEA-W 0.02578 0.50550 0.00271 4.88694
4EDMMLA 0.10395 2.71061 0.01302 9.27789
4EDMMEA-L 0.02580 0.50635 0.00271 4.89684
4GUOMEA-W 0.02578 0.50551 0.00271 4.88695
4SCFMEA-W 0.03347 0.86398 0.00401 5.84097
4SCFMEA-L 0.03349 0.86637 0.00402 5.84066
Table 3. Different error measures [80] for the four cases described by ten Cate et al. [45] for various combinations of momentum exchange and forcing schemes. The error in terminal velocity is given relative to the analytical values according to Abraham [60]. The results in the table refer to N = 8 .
Table 3. Different error measures [80] for the four cases described by ten Cate et al. [45] for various combinations of momentum exchange and forcing schemes. The error in terminal velocity is given relative to the analytical values according to Abraham [60]. The results in the table refer to N = 8 .
CaseForcingMomentum ExchangeRMSEPCMArea between CurvesError in % Maximum Velocity
1EDMMEA-W 0.03241 4.35586 0.00422 8.24923
1EDMMLA 0.03399 2.84034 0.00448 8.27399
1EDMMEA-L 0.03241 4.35589 0.00422 8.24930
1GUOMEA-W 0.03241 4.35586 0.00422 8.24923
1SCFMEA-W 0.03442 4.41524 0.00457 8.49026
1SCFMEA-L 0.03442 4.41527 0.00457 8.49039
2EDMMEA-W 0.03224 1.11747 0.00347 4.32345
2EDMMLA 0.05467 1.03172 0.00550 4.44263
2EDMMEA-L 0.03224 1.11747 0.00347 4.32358
2GUOMEA-W 0.03224 1.11748 0.00347 4.32345
2SCFMEA-W 0.03326 1.10431 0.00359 4.59350
2SCFMEA-L 0.03326 1.10431 0.00359 4.59360
3EDMMEA-W 0.02349 0.64709 0.00162 4.43872
3EDMMLA 0.06483 1.49952 0.00563 4.83740
3EDMMEA-L 0.02349 0.64703 0.00162 4.43877
3GUOMEA-W 0.02349 0.64709 0.00162 4.43872
3SCFMEA-W 0.02079 0.76920 0.00152 4.67580
3SCFMEA-L 0.02079 0.76898 0.00152 4.67593
4EDMMEA-W 0.08260 1.06138 0.00223 4.96350
4EDMMLA 0.08638 1.67762 0.00927 6.09874
4EDMMEA-L 0.08259 1.06147 0.00223 4.96360
4GUOMEA-W 0.08260 1.06139 0.00223 4.96350
4SCFMEA-W 0.08877 1.03403 0.00251 5.16140
4SCFMEA-L 0.08875 1.03418 0.00251 5.16147
Table 4. Mean values of the results presented in Table 2 and Table 3. The averages are taken over the four application cases.
Table 4. Mean values of the results presented in Table 2 and Table 3. The averages are taken over the four application cases.
NForcingMomentum ExchangeRMSE (Mean)PCM (Mean)Area between Curves (Mean)Error in % Maximum (Mean)
1EDMMEA-W 0.04519 1.90257 0.00442 5.47702
1EDMMLA 0.08455 2.97017 0.00927 7.23193
1EDMMEA-L 0.04521 1.89796 0.00443 5.48345
1GUOMEA-W 0.04519 1.90250 0.00442 5.47703
1SCFMEA-W 0.05588 2.55610 0.00598 6.92398
1SCFMEA-L 0.05590 2.56138 0.00598 6.92947
8EDMMEA-W 0.04268 1.79545 0.00288 5.49373
8EDMMLA 0.05997 1.76230 0.00622 5.91319
8EDMMEA-L 0.04268 1.79546 0.00288 5.49381
8GUOMEA-W 0.04268 1.79545 0.00288 5.49373
8SCFMEA-W 0.04431 1.83069 0.00305 5.73024
8SCFMEA-L 0.04430 1.83069 0.00305 5.73035
Table 5. Mean experimental order of convergence (EOC) for different error and similarity measures regarding N = 8 . Results are averaged over all values for a given combination of forcing and momentum exchange scheme as well as all cases.
Table 5. Mean experimental order of convergence (EOC) for different error and similarity measures regarding N = 8 . Results are averaged over all values for a given combination of forcing and momentum exchange scheme as well as all cases.
ForcingMomentum ExchangeEOC (RMSE)EOC (RMSE)EOC (Area between Curves)EOC (Area between Curves)
EDMMEA-W 1.79 1.81 2.06 1.39
EDMMLA 1.37 2.00 1.49 1.66
EDMMEA-L 1.79 1.81 2.05 1.36
GUOMEA-W 1.79 1.81 2.06 1.39
SCFMEA-W 1.57 1.03 1.72 1.55
SCFMEA-L 1.57 1.03 1.72 1.55
Table 6. Deviation of of the drag coefficients calculated by homogenized lattice Boltzmann method (HLBM) simulations to different drag correlations.
Table 6. Deviation of of the drag coefficients calculated by homogenized lattice Boltzmann method (HLBM) simulations to different drag correlations.
μ f in kg m 1 s 1 C D (Simulated) StokesAbrahamSchiller and Naumann
ReError in %ReError in %ReError in %
2.0 × 10 2 86.6942 0.26 7.72 0.23 24.93 0.24 17.33
1.5 × 10 2 56.7155 0.45 7.33 0.40 17.98 0.42 8.42
1.0 × 10 2 29.6648 1.02 26.31 0.84 14.14 0.90 2.68
5.0 × 10 3 10.5937 4.09 80.42 2.90 9.38 3.08 2.74
2.5 × 10 3 4.5797 16.35 212.00 9.18 1.61 9.57 6.94
1.25 × 10 3 2.2965 65.40 525.82 26.57 3.29 26.83 5.34
6.25 × 10 4 1.2539 261.60 1266.86 70.48 0.79 69.50 3.52
3.125 × 10 4 0.7339 1046.40 3100.03 173.68 11.84 170.80 14.74
1.5625 × 10 4 0.5213 4185.60 8992.70 404.05 15.27 406.40 14.28
7.8125 × 10 5 0.4543 16,742.4031,595.60 900.47 8.31 948.67 1.76
Table 7. Maximum settling velocity for μ f = 0.005   k g   m 1   s 1 for a chosen lattice velocity u L (see Equation (27)). Additionally, given is the deviation to the velocity obtained with u L = 0.0025 and the change of maximum velocity to the value obtained with the next smaller u L .
Table 7. Maximum settling velocity for μ f = 0.005   k g   m 1   s 1 for a chosen lattice velocity u L (see Equation (27)). Additionally, given is the deviation to the velocity obtained with u L = 0.0025 and the change of maximum velocity to the value obtained with the next smaller u L .
u L 0.0025 0.005 0.01 0.02 0.04 0.08 0.16
u in m   s 1 0.02975 0.03002 0.03042 0.03102 0.031981 0.03375 0.03768
Deviation in %0 0.89 2.24 4.25 7.50 13.45 26.64
Change in %- 0.89 1.33 1.97 3.11 5.54 11.6
Table 8. Parameters of the simulations investigating the tubular pinch effect along with the results for multiple cases according to Inamuro et al. [36] and the one described by Li et al. [35].
Table 8. Parameters of the simulations investigating the tubular pinch effect along with the results for multiple cases according to Inamuro et al. [36] and the one described by Li et al. [35].
Case τ Δ p DL d p / D y c / D Error in %
Inamuro et al. case 2 1.4 8.167 × 10 4 200200 0.25 0.2642 3.33 %
Inamuro et al. case 6 0.757 2.337 × 10 4 200200 0.25 0.2548 5.83 %
Inamuro et al. case 11 1.4 1.633 × 10 3 200400 0.25 0.2726 4.35 %
Inamuro et al. case 14 0.95 3.207 × 10 3 100400 0.5 0.3590 5.48 %
Li et al. 0.75 2.670 × 10 4 100400 0.25 0.2756 4.10 %
Table 9. Simulated hindered settling velocity regarding different Re and ϕ in comparison to different correlations.
Table 9. Simulated hindered settling velocity regarding different Re and ϕ in comparison to different correlations.
Re ϕ usim in m s 1 Error in % (uRZ, uS)Error in % (uO, uS)Error in % (uBM, uS)Error in % (uRZ, uS-N)Error in % (uO, uS-N)Error in % (uBM, uS-N)
0.53 0.05 0.00071 23.38 10.17 3.97 15.90 1.46 5.35
0.53 0.1 0.00055 24.69 16.87 7.46 17.28 8.81 1.51
0.53 0.15 0.00045 20.36 19.03 6.15 12.46 11.18 2.95
0.53 0.2 0.00037 15.63 21.35 5.47 7.19 13.73 3.70
0.53 0.25 0.00031 5.06 18.20 0.80 4.52 10.27 10.57
5.29 0.05 0.00691 46.58 34.75 30.24 20.82 4.00 2.63
5.29 0.1 0.00551 48.15 37.75 30.71 22.56 8.41 1.95
5.29 0.15 0.00444 48.68 40.61 31.16 22.72 12.61 1.29
5.29 0.2 0.00367 47.08 41.07 29.17 19.62 13.29 4.22
5.29 0.25 0.00295 46.26 41.77 28.24 17.62 14.32 5.58
49.46 0.05 0.06962 74.81 67.71 65.47 18.17 3.13 10.26
49.46 0.1 0.05656 76.34 68.63 65.08 21.76 0.19 11.53
49.46 0.15 0.04757 76.80 68.71 63.73 21.83 0.06 15.84
49.46 0.2 0.03889 77.69 69.34 63.14 23.30 2.07 17.71
49.46 0.25 0.03300 77.49 68.01 60.58 20.96 2.17 25.90
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Trunk, R.; Weckerle, T.; Hafen, N.; Thäter, G.; Nirschl, H.; Krause, M.J. Revisiting the Homogenized Lattice Boltzmann Method with Applications on Particulate Flows. Computation 2021, 9, 11. https://0-doi-org.brum.beds.ac.uk/10.3390/computation9020011

AMA Style

Trunk R, Weckerle T, Hafen N, Thäter G, Nirschl H, Krause MJ. Revisiting the Homogenized Lattice Boltzmann Method with Applications on Particulate Flows. Computation. 2021; 9(2):11. https://0-doi-org.brum.beds.ac.uk/10.3390/computation9020011

Chicago/Turabian Style

Trunk, Robin, Timo Weckerle, Nicolas Hafen, Gudrun Thäter, Hermann Nirschl, and Mathias J. Krause. 2021. "Revisiting the Homogenized Lattice Boltzmann Method with Applications on Particulate Flows" Computation 9, no. 2: 11. https://0-doi-org.brum.beds.ac.uk/10.3390/computation9020011

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop