Article
Numerical study and acceleration of LBM-RANS simulation of turbulent flow

https://doi.org/10.1016/j.cjche.2017.05.013Get rights and content

Abstract

The coupled models of LBM (Lattice Boltzmann Method) and RANS (Reynolds-Averaged Navier–Stokes) are more practical for the transient simulation of mixing processes at large spatial and temporal scales such as crude oil mixing in large-diameter storage tanks. To keep the efficiency of parallel computation of LBM, the RANS model should also be explicitly solved; whereas to keep the numerical stability the implicit method should be better for RANS model. This article explores the numerical stability of explicit methods in 2D cases on one hand, and on the other hand how to accelerate the computation of the coupled model of LBM and an implicitly solved RANS model in 3D cases. To ensure the numerical stability and meanwhile avoid the use of empirical artificial limitations on turbulent quantities in 2D cases, we investigated the impacts of collision models in LBM (LBGK, MRT) and the numerical schemes for convection terms (WENO, TVD) and production terms (FDM, NEQM) in an explicitly solved standard kε model. The combination of MRT and TVD or MRT and NEQM can be screened out for the 2D simulation of backward-facing step flow even at Re = 107. This scheme combination, however, may still not guarantee the numerical stability in 3D cases and hence much finer grids are required, which is not suitable for the simulation of industrial-scale processes. Then we proposed a new method to accelerate the coupled model of LBM with RANS (implicitly solved). When implemented on multiple GPUs, this new method can achieve 13.5-fold acceleration relative to the original coupled model and 40-fold acceleration compared to the traditional CFD simulation based on Finite Volume (FV) method accelerated by multiple CPUs. This study provides the basis for the transient flow simulation of larger spatial and temporal scales in industrial applications with LBM–RANS methods.

Introduction

Fluid mixing is ubiquitous in chemical and process industries, and Computational Fluid Dynamics (CFD) is playing important roles in evaluating the mixing for various processes [1], [2], [3]. Traditional Finite Volume (FV)-based CFD models prove to be feasible for the steady-state simulation of mixing processes or the transient simulation of physically fast mixing processes. In the latter case, the flow can reach equilibrium within a relatively short time and can be effectively simulated by traditional CFD with current computational resources. For instance, the mixing of two or more miscible fluids of similar densities or viscosities can be treated as a single-phase flow of multiple components, and the FV-based CFD simulation was always decoupled into two procedures: the flow field was first obtained through steady-state simulation, and then the scalar transport equations, e.g., the concentration of components, were separately solved in a time-dependent manner with the calculated flow field [4], [5], [6]. The computational resource needed is relatively small and some critical parameters such as mixing time can be predicted with certain accuracy [7], [8]. On the other hand, in case of two-phase flow such as the mixing of gas bubbles and liquids, the two-fluid models may have to be applied to differentiate the movement tendencies of the mixing fluids at the governing equations level [9], [10], let alone Direct Numerical Simulation (DNS) [11] or Eulerian-Lagrangian simulation [12] of multiphase flow. In this case, the problems of higher computational cost and difficulties in numerical stability may become pronounced since the equations are strongly coupled, and the modeling accuracy is also sensitive to the turbulence models or interphase interaction models.

However, for some mixing processes with highly transient characteristics or the so-called active mixing [13], [14] in which the distribution of miscible solvent has non-negligible influences on flow field, even solving the equations of time-dependent single-phase turbulent flow of multiple components is time-consuming, let alone using the two-fluid models. A typical example is the mixing of crude oil in large–diameter storage tanks in which it takes hours or days to complete the mixing of two or multiple crude oils of different densities [15], [16] and only the transient flow simulation is relevant. To effectively simulate the mixing of two crude oils of different densities in storage tanks with acceptable computational cost, the Boussinesq approximation is generally applied, keeping only the density difference in the gravitational term as a compromise between computational cost and accuracy [17], [18] . Even though, it is still unfeasible to simulate transient mixing processes in industrial scale by solving FV-based CFD. For example, it took about 30 CPU hours to update only one second of physical time when using FV-based CFD and Detached Eddy Simulation (DES) to solve the single-phase flow and scalar transport equations for the passive mixing process (the solvent distribution has less impacts on flow field) [13], [14] in a cylindrical stirred tank (diameter T: 0.15 m and height H = T, cells number: 500000, time step: 0.001 s) [19]. For those active mixing processes, it may take about 432000 CPU hours or more to simulate only one day of physical time for crude oil mixing in a storage tank (more than 104 m3), even using only 500000 cells in total. It is therefore not affordable or a challenging issue for traditional FV-based CFD methods to simulate transient mixing processes in such large-scale industrial applications.

The parallel computation of Lattice Boltzmann Method (LBM) simulation accelerated by Graphic Processing Units (GPU) might be a suitable alternative. LBM is distinguished by the explicitly time-marching and highly localized computation properties, and proves to be more suitable for parallel computation and transient flow simulation [20], [21], in particular when utilizing the modern high-performance GPU cards. The computational speed of GPU-accelerated LBM is reported to be able to achieve nearly two-orders of magnitude faster than the LBM simulation using contemporary mainstream CPU [22]. In addition, the flow in crude oil mixing process in industrial scale is always turbulent. Reynolds-Averaged Navier–Stokes (RANS) turbulence modeling is more practical and economical than DNS [17], [18] or Large Eddy Simulation (LES) [23], [24], and may also be the only feasible approach for large-scale industrial turbulent flow simulation. The standard kε model, the most popular RANS model, solves the transport equations of turbulent kinetic energy or turbulent dissipation rate [25]. Compared to the FV-based CFD–RANS simulation, the GPU-accelerated LBM–RANS simulation is of practical interest for industrial applications.

Nevertheless, there are still several issues that need to be settled before exploiting the advantage of LBM–RANS simulation. First, to maintain the parallelism and computational efficiency of LBM, the RANS turbulence model may have to be explicitly solved. The LBM simulation using the Lattice Bhatnagar–Gross–Krook (LBGK) collision model has been coupled with the explicitly solved RANS model and validated in some benchmark problems of steady flow such as the turbulent pipe flows or separated flows [26], [27] or the turbulent flow around the aircraft wing of NACA 4412 profiled with surface grid refinement [28]. However, the equations in RANS models, which consist of a set of non-linear source terms, are strongly coupled with each other and the momentum conservation equations. This nature of RANS models makes the explicit computation of RANS models difficult to converge [29]. As a result, some artificial limitations were usually imposed for turbulent quantities to enhance the numerical stability. For example, the turbulent quantities were only allowed to be updated when the newly calculated turbulent viscosity was less than an upper threshold (e.g., 0.25 in lattice unit) [26], or the turbulent dissipation rate or turbulent kinetic energy was larger than a lower threshold [28]. The artificial numerical limitation may influence the accuracy of transient simulation, especially for the flow of longer evolution time in crude oil mixing. It is necessary to find a suitable combination of numerical treatments for the coupled model of LBM and explicitly solved RANS model, ensuring the numerical stability without resorting to the empirical artificial limitation. In this study, we first investigated the impacts of a series of factors such as the collision models in LBM, the discretization schemes of convection terms and the computational methods of production terms in standard kε model on numerical stability in a 2D case. These factors have not yet received a systematic study in the coupling of LBM and explicitly solved standard kε model. Rather than testing above-mentioned each single component which has been already reported in literature, we aim to find a suitable combination of these numerical schemes to ensure the numerical stability without any artificial limitations. A classical benchmark problem with complex flow separation, i.e. the turbulent flow over a backward-facing step was evaluated in the 2D simulations, and an optimized scheme combination was screened out for the 2D simulation cases.

Secondly, the optimized schemes for the 2D cases, however, still cannot guarantee the numerical stability in the 3D cases in our tests. Further refinement of the grids might improve the stability but at expense of much more computational resources. Hence it is still not suitable for the simulation of industrial-scale processes. To keep the numerical stability, the implicit method for RANS should be better, which may however greatly slow down the computational speed and compromise the advantage of explicit computation and parallelism of LBM. The key question is how to harmonize the explicitly solved LBM and implicitly solved RANS to accelerate the computation. We then proposed a new accelerating method to synchronize the LBM and RANS computations, or in other words, making the RANS simulation compatible with the LBM computation. This method shows the potential to be applied for accelerating the computation of long time active mixing processes in an industrial-scale mixing tank.

Section snippets

The Lattice Boltzmann Method

The governing equation of the Multiple Relaxation Time (MRT)–LBM [30], [31] isfx+eiΔtt+Δt=fxtM1Smmeq+I12SFwhere f is the velocity distribution function vector, ei is the discrete velocity vector, M is the transformation matrix, and S is the diagonal collision matrix with adjustable parameters. m, the velocity distribution function vector in moment space, is computed by m = Mf. meq is the equilibrium moment, I is the identity matrix, and F is the force term in moment space [32]. The D2Q9 model

Evaluation of Different Numerical Treatments

In this section, we investigated the influence of several key factors on the numerical stability of the coupled model of LBM and explicitly solved standard k–ε model, e.g., the collision models in LBM, the discretization schemes of convection terms and the computational methods of production terms in the standard k–ε model. The impacts of the artificial limits were also investigated. For the collision models in LBM, we compared the simulation using the LBGK model with that of the MRT model. For

A New Acceleration Method for Implicit 3D Simulation

When the RANS turbulence model is implicitly solved, the LBM–RANS simulation is more stable and does not require the artificial numerical limitation. However, there are still several open questions. The implicit method in the RANS simulation always needs many more iterations to ensure the computational convergence and the computational speed is relatively much slower. The LBM simulation, on the other hand, is an explicitly time-marching approach with much faster computational speed. Moreover,

Conclusions

With natural parallelism and explicit time-marching properties, the LBM–RANS models may hitherto be more practical and economical for LBM methods to be applied in industrial-scale turbulent flows, despite the active researches on DNS or LES in academia. To keep the efficiency of parallel computation of LBM, the RANS model should also be explicitly solved, which greatly affects the numerical stability; whereas to keep the numerical stability, the implicit method should be better for RANS, which

Nomenclature

    aP

    Coefficient of node P in discretization

    anb

    Coefficient of neighborhood nodes in discretization

    B

    Model parameter (= 5.0)

    Cμ

    Model parameter in RANS model (= 0.09)

    c1

    Model parameter in RANS model (= 1.44)

    c2

    Model parameter in RANS model (= 1.92)

    D2Q9

    Lattice model in 2D

    D3Q19

    Lattice model in 3D

    Dϕ

    the effective viscosity

    ER

    Expansion rate (= 1.5)

    ei

    Discrete velocity vector, m·s 1

    F

    Force term in moment space

    f

    Velocity distribution function vector

    Gk

    Generation of turbulent kinetic energy

    g

    Gravity acceleration, m·s 2

    Hin

References (41)

Cited by (12)

  • Regime mapping of multiple breakup of droplets in shear flow by phase-field lattice Boltzmann simulation

    2021, Chemical Engineering Science
    Citation Excerpt :

    c = 0 corresponds to the Langmuir adsorption isotherm, and e = 0 suggests the same solubility of soluble surfactants in both the bulk phases. To improve numerical accuracy and stability, the Multiple-Relaxation-Time (MRT) method was used to replace Eqs. (11)–(13) and (14)–(17) (Premnath and Abraham, 2007; Shu and Yang, 2013, 2018b). Fig. 1 illustrates a schematic map of droplet in shear flow within two parallel moving plates.

  • Utilizing lattice Boltzmann method for heat transfer analysis in solar thermal systems: A review

    2021, Sustainable Energy Technologies and Assessments
    Citation Excerpt :

    Another challenge dealt with when utilizing LBMs is reconciling the explicit nature of LBM with the implicit treatment needed for turbulent flow regimes e.g., when resolving non-linear source terms in Reynolds-averaged Navier–Stokes (RANS) models.Said source terms, are strongly coupled with both each other and the momentum conservation equations and makes the explicit computation of RANS models difficult to converge; As a result, some artificial limitations are usually imposed for turbulent quantities to enhance the numerical stability. Such artificial numerical limitation may have an adverse influence on the accuracy of transient simulation [56]. Similarly, difficulties arise when high thermal fluxes and boiling heat transfer phenomena need to be modeled.

  • High-resolution simulation of oscillating bubble plumes in a square cross-sectioned bubble column with an unsteady k-ε model

    2021, Chemical Engineering Science
    Citation Excerpt :

    The grid size, dx, was set as 2.5 mm unless otherwise mentioned. The LBM model was used to solve the mixture model and the classical Finite Volume Method (FVM) was used to solve the continuity equation of the gas phase as well as the transport equations of the turbulent quantities, i.e., the turbulent kinetic energy and the turbulent dissipation rate, as described in Shu and Yang (2018b). The Weighted Essentially Non-Oscillatory (WENO) scheme (Jiang and Shu, 1996) was used to discretize the convection terms of these transport equations, and the diffusion terms were discretized with a fourth-order central difference method.

  • Numerical study on flow field and pollutant dispersion in an ideal street canyon within a real tree model at different wind velocities

    2021, Computers and Mathematics with Applications
    Citation Excerpt :

    Researchers have studied the convection and diffusion of pollutants in street canyons under the influence of trees by numerical simulation [11–13], experimental techniques [14,15] and field observations [16]. Numerical simulations are widely used in different fields [17–20], including the numerical simulation of wind fields and pollutant concentration fields in street canyons, due to their advantages over experimental technology in resources consumption and the more abundant acquisition of physical information [21–25]. To adopt an appropriate turbulence model, mesh and Schmidt number, the prediction of pollutant concentrations, velocity and turbulent kinetic energy in a street canyon by numerical simulation must be consistent with wind tunnel test results [26].

  • GPU-accelerated transient lattice Boltzmann simulation of bubble column reactors

    2020, Chemical Engineering Science
    Citation Excerpt :

    A uniform grid in each direction (Δx = 2.5 mm) is set as the default configuration. The transport equations of the continuity equation of gas, turbulent kinetic energy and turbulent dissipation rate were solved by the traditional Finite Volume Method (Shu and Yang, 2018b). The convection term of these transport equations were discretized with the Weighted Essentially Non-Oscillatory (WENO) scheme (Jiang and Shu, 1996), and the fourth-order central difference method was used to discretize diffusion terms.

  • GPU-accelerated large eddy simulation of stirred tanks

    2018, Chemical Engineering Science
    Citation Excerpt :

    Lattice Boltzmann Method (LBM) is a novel approach with high parallelism in nature and explicit time-marching properties. Over the past two decades, LBM has been used as a fluid flow solver for the direct numerical simulation of immiscible multiphase flow such as bubble dynamics (Shu and Yang, 2013), high Reynolds number (up to 10 million) turbulent flows (Shu and Yang, 2018) and the simulation of stirred tanks. There are several pioneering studies about the parallel computation of LBM simulation of stirred tanks based on many CPUs (Central Processing Unit) (Derksen and Van den Akker, 1999; Derksen, 2003, 2011, 2012; Gillissen and Van den Akker, 2012; Hartmann et al., 2004a, 2004b, 2006a, 2006b).

View all citing articles on Scopus

Supported by the National Key Research and Development Program of China (2017YFB0602500), National Natural Science Foundation of China (91634203 and 91434121), and Chinese Academy of Sciences (122111KYSB20150003).

View full text