GPU-accelerated transient lattice Boltzmann simulation of bubble column reactors

Journal Pre-proofs GPU-Accelerated Transient Lattice Boltzmann Simulation of Bubble Column Reactors Shuli Shu, Jingchang Zhang, Ning Yang PII: DOI: Reference:

S0009-2509(19)30926-1 https://doi.org/10.1016/j.ces.2019.115436 CES 115436

To appear in:

Chemical Engineering Science

Received Date: Revised Date: Accepted Date:

22 April 2019 18 August 2019 11 December 2019

Please cite this article as: S. Shu, J. Zhang, N. Yang, GPU-Accelerated Transient Lattice Boltzmann Simulation of Bubble Column Reactors, Chemical Engineering Science (2019), doi: https://doi.org/10.1016/j.ces. 2019.115436

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2019 Published by Elsevier Ltd.

GPU-Accelerated Transient Lattice Boltzmann Simulation of Bubble Column Reactors

Shuli Shua,b, Jingchang Zhanga,b, Ning Yanga,b*

a

State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, P. O. Box 353, Beijing 100190, China b

School of Chemical Engineering, University of Chinese Academy of Sciences, Beijing 100049, China

Correspondence concerning this article should be addressed to Ning Yang at [email protected]

1

Abstract Design parameters for Bubble Column Reactors such as backmixing and life/dark cycles of microalgae cells are pertinent to transient flow dynamics and mesoscale coherent structures. Transient simulation is therefore of practical significance. Traditional CFD solvers are generally implemented on multi-core CPU workstations, which is time-consuming and difficult to meet the requirement of industry applications. We developed a GPU-accelerated transient LBM simulation of gas-liquid flow, achieving 5,000 times acceleration in comparison with the CFD simulation. A notable feature of this method, which allows a tunable and smaller dimensionless relaxation time, is the remarkable decrease of required mesh number in LBM simulation by several orders of magnitude. The solver is validated against experimental data, resolving two mechanisms of mesoscale structure, i.e., vortex rotation and radial migration, in the intensification of macromixing. The computation acceleration and resolution of mesoscale structures demonstrates the potential of this approach for fast simulation of industrial applications.

Keywords: Transient simulation; GPU acceleration; Multiphase flow; Lattice Boltzmann Method; Bubble Column Reactor

2

1.

Introduction The transient behavior of mesoscale coherent structures such as vortical cells is of great

significance and receiving more and more attention in rational design, scale-up and optimization of Bubble Column Reactors (BCRs). The transient coherent structures have direct impacts on backmixing of liquid phase, mass and heat transfer, reaction selectivity and yields (Buwa and Ranade, 2002; Delnoij et al., 1999; Shu et al., 2019). For example, the light/dark cycles of microalgae cells are governed by their exposure to light sources, which fundamentally depends on the transient mesoscale coherent flow structures (Kliphuis et al., 2010). Over several decades, Computational Fluid Dynamics (CFD) simulation has been extensively developed as a powerful tool to predict the mean flow profiles (Joshi, 2001; Laborde-Boutet et al., 2009; Tabib et al., 2008; Xiao et al., 2013; Yang et al., 2011). However, prediction of the dynamics of mesoscale coherent flow structures in large-scale BCRs is still challenging for traditional CFD simulation as both the computation of longer physical time (more than 1,000 seconds) and finer grid resolution (~6.5 mm in radial direction and ~30 mm in axial direction) are needed to fully resolve the dynamic behavior (Pfleger and Becker, 2001). Several weeks or months of computational time are still needed to complete a transient simulation of about several hundred seconds of physical time on relatively coarser grids (Masood et al., 2015; McClure et al., 2014; Silva et al., 2014), not to mention the simulation of longer physical time with finer grid. This lower computational efficiency results not only from the limited performance of traditional computer hardware, i.e., Central Processing Units (CPU), but the computational algorithm involving complex iterations in mathematical models and fluid solvers, and in particular, the incompatibility between computer hardware and computational algorithm (Li et al., 2013). The question is therefore how to coordinate 3

the above issues to effectively improve computation speed and meet the practical requirement of chemical and process industry. While Direct Numerical Simulation (DNS) (Deen and Kuipers, 2013; Dijkhuizen et al., 2010; Roghair et al., 2011; Shu and Yang, 2013) or Eulerian-Lagrangian model (Song et al., 2013; Witz et al., 2016; Xue et al., 2017a, b) are more accurate yet unaffordable, the phase-volume-averaged models such as mixture or two-fluid models are more feasible. While the mixture model is much simpler due to the relatively weak coupling among phase equations, it has been reported to give reasonable predictions for BCRs (Sanyal et al., 1999; Sokolichin et al., 1997). This simplicity in model structure makes it more suitable for computation acceleration with Graphic Processing Units (GPU). Consequently, the mixture model can serve as a suitable platform to accelerate the simulation of transient hydrodynamics of gas-liquid flows, in view of the reasonable simulation accuracy and the compatible model structure with GPU-accelerated parallel computation (Chen et al., 2005; Chen et al., 2004; Icardi et al., 2014; Swiderski et al., 2016). In recent years, GPU with stunning computational capability has become a mainstream hardware platform of high-performance computation (HPC). The theoretical peak performance of a single latest NVIDIA V100 GPU card has been faster than that of the world’s fastest CPU-based supercomputer system in 2001. The compatibility between the solver and hardware platform enables the acceleration of computation. Compared to traditional finite volume method (FVM), Lattice Boltzmann Method (LBM) is more suitable for GPU computation due to its natural parallelism. For example, more than 1500 folds of acceleration was achieved by a GPU-based LBM simulation of single-phase turbulent flow within a stirred tank, enabling the high-resolution simulation by using only a desktop computer without the need of much finer grid resolution and 4

supercomputers (Shu and Yang, 2018a). Witz et al. (2016) reported that the Eulerian-Lagrangain simulation of bubbly flow, when accelerated with GPU-based LBM simulations, can be used to simulate an industrial-scale stirred tank bioreactor. It is therefore very promising to accelerate the transient gas-liquid flows with the combination of LBM algorithm with GPUs. However, it is not trivial to directly apply the mixture model for the simulation of gas-liquid flows with LBM due to numerical instability. The explicit nature of LBM makes the computation difficult to converge for the gas-liquid flows of high density ratios, especially for the unsteady flows (Shu and Yang, 2013; Wang and Wang, 2005; Yuan and Schaefer, 2006). In this case, modification of the collision terms, a general method to enhance numerical stability in single-phase flows, may not be adequate for two-phase flow. Semi-implicit discretization of temporal terms (Sankaranarayanan et al., 2002; Sankaranarayanan and Sundaresan, 2008) or solving an additional pressure-Poisson equation (Inamuro et al., 2004) are two main alternatives to improve the convergence at the expense of computational efficiency. The former needs to solve additional nonlinear equations, while the latter only needs to solve linear equations and is therefore preferred. However, the pressure term recovered from the latter was divided by the dimensionless relaxation time τ. The parameter is related to the fluid kinematic viscosity v through v=(τ-0.5)Δx2/3Δt. This indicates that τ can only be set as unity to recover the correct pressure term in macroscale equations. In fact, the combination of grid size Δx and time step Δt in LBM should not only meet the requirement of numerical stability, i.e., the so-called incompressible limit: Δx/Δt should be larger than about 10 times of the local velocity u, (Δx/Δt>10u), but also need to give a correct fluid kinematic viscosity. For a given system of specific fluid kinematic viscosity, Δx should be less than 3v/(10u(τ-0.5)). Therefore, the required cells number N for the simulation of a given three5

dimensional system is related to the dimensionless relaxation time τ, i.e., N∝(τ-0.5)3, which implies that the grid nodes required by the fixed τ=1 is 1,000,000 times larger than that of τ=0.505. Moreover, the grid size Δx should be less than 0.006 mm and the corresponding time step Δt should be smaller than 0.006 ms for air-water systems, provided that τ is set as unity (Inamuro et al. (2004) and the local water velocity is at the level of 0.1 m/s or higher. In other words, about 4.6 billion grid nodes are needed for a fluid domain of 1 cm3. The LBM algorithm of Inamuro et al. (2004) is therefore only suitable for small-scale simulation. In contrast, the computational resource required could be theoretically decreased in several orders of magnitude, provided that a new LBM algorithm with tunable and smaller τ could be developed to use a smaller number of grid nodes and meanwhile recover the correct pressure term without any other constraints for the practical applications. The objective of this work is to develop a new GPU-accelerated LBM approach to accelerate the simulation of gas-liquid flows with mixture models. The macroscopic governing equations of the mixture model, especially the pressure term, can be correctly recovered from this new LBM algorithm which allows a tunable smaller τ. The algorithm was implemented on multiple GPUs, and the GPU-accelerated simulation was compared with CFD simulation. Then the simulation of a labscale BCR was conducted for validation. Finally, the dynamics and anisotropic fluctuations of low frequency coherent structures in gas-liquid flows in BCRs are investigated with the new model.

2.

Mathematical model

2.1 Governing equations In this work, the gas-liquid two-phase flow is simulated with a simplified mixture model (Manninen et al., 1996) in which only the conservation equations for gas-liquid mixture rather than 6

the individual phases are solved. The dispersion of gas phase is modeled by an additional transport equation for gas volume fraction, as listed in Table 1Tables Table 1. For simplicity, a constant slip velocity (ulg=0.2 m/s), which can only describe a

monodispersed bubble population, is assumed to describe the drag effect (parallel to the motion of bubbles). Lift force (perpendicular to the motion of bubbles) and virtual mass force (due to acceleration or deacceleration of bubbles) are omitted. As suggested by Sokolichin and Eigenberger (1999), the term  Dm illustrated in Table 1 is neglected in this work. In order to improve numerical stability, the gravity term in momentum equations of the mixture phase can be partly absorbed into the pressure term. Hence, the governing equations for the mixture phase can be reformulated as   um  0

(1)

um p*   umum     vm  vt  um T um   Fb   t m





(2)

where Fb=αgg(ρg-ρl)/ ρm and p*=p-ρlg·r. The standard mixture k-ε model with the standard wall function is used to model the effect of turbulence (Laborde-Boutet et al., 2009).

2.2 The algorithm based on Lattice Boltzmann Method

To eliminate the dimensionless relaxation time τ in the pressure term, we first segregate the pressure term from the equilibrium velocity distribution function, and introduce the term as an *

external force to predict the intermediate velocity of mixture phase um . The intermediate velocity is then corrected by other external forces of previous time step except for pressure. Then the pressure n1

*

term is updated by um . Finally, the velocity of mixture phase of the next time step um can be *

updated by correcting um with the updated pressure. In this case, the temporal scheme of pressure term is semi-implicit. 7

*

In this way, the intermediate velocity for mixture phase um is firstly calculated from the corresponding LBM equation Eq.(3) by treating the pressure term as an external force:



 



f  x  ci t, t  t   f  x, t   M 1S m  meq  I  12 S Fˆ

(3)

where f is the velocity function vector, ci is the discretized velocity vector, M is the transfer matrix, m is the velocity function vector in moment space, meq is the equilibrium moment, I is the unit matrix and F is the force term in the momentum space. The diagonal matrix S is related to the relaxation time, and S=diag(S1,S2,…,Sb) in which b denotes the number of discretized velocity vector. In the D3Q19 lattice model, S9=S11=S13-15=1/τ, S0=S3=S5=S7=0, S1=1.19, S2=S10=S12=1.4, S4=S6=S8=1.2, S16-18=1.98. τ is the non-dimensional relaxation time and related to the kinematic viscosity by v= (τ-0.5)Δx2/(3Δt) and v=(μm+μt)/ρm, and Δx is the width of grid. The main difference between the standard LBM equation and Eq. (3) in this work is the equilibrium moments meq. In this work, the moments relevant to the pressure term are modified, so that the reduced form of Eq.(2) neglecting the pressure term can be recovered from Eq. (3). Hence, meq can be rewritten as,

 0,19(u x2  u 2y  u z2 ), 475 / 63(u x2  u 2y  u z2 ), u x ,  2 u x , u y ,  2 u y , u z ,  2 u z ,  3 3 3  m eq    2u x2  u 2y  u z2 ,0, u 2y  u z2 ,0, u xu y , u y u z , u z u x ,0,0,0   

(4)

The transformation matrix M can be referred to Shu and Yang (2018a). The force terms in the equilibrium space Fˆ is written as,

2 2 2 Fˆ  (0, 0,38F  u, Fx ,  Fx , Fy ,  Fy , Fz ,  Fz , 6 Fx u x  2F  u, F  u  3Fx u x , 3 3 3 2 Fy u y  2 Fz u z , Fz u z  Fy u y , Fx u y  Fy u x , Fz u y  Fy u z , Fz u x  Fx u z , 0, 0, 0)

(5)

where F is  p n F   Fbn - n  m 

where n represents the last time step. 8

  t  

(6)

Following the implementation of collision and propagation process in Eq.(3), the intermediate *

velocity for mixture phase um can be calculated and corrected by , 1   um*    fi ci  Fb  / 0 2  

(7)

where  0 is unity. Then the pressure field can be obtained from 1 p *   um  t m

(8)

Finally, the mixture velocity is updated by umn 1  um* 

1 p t 2 m

(9)

The Chapman-Enskog expansion in the Supplementary Information shows that the reduced form of Eqs. (1)-(2) can be fully recovered from Eqs. (3)-(9). In other words, the pressure term in the macroscopic governing equations can be correctly recovered and is independent of τ. Therefore, the non-dimensional relaxation time τ in the proposed LBM algorithm is tunable and no longer constrained as a constant to recover a correct pressure term.

3. Flow systems and numerical setups

A centrally-aerated two-dimensional bubble column of Pfleger et al. (1999), as illustrated in Figure 1, was simulated to validate the new solver. The width (W), depth (D) and height (H) of the bubble column are 200 mm (W), 50 mm (D) and 1500 mm (H), respectively. The sparger (24mm ×12mm) located at the column bottom center introduces air into the column. The fluid properties of gas and liquid phases are listed in Table 2. The liquid level was set as 450 mm with a H/W ratio 2.25 in the model validation and computational speed tests. Then different liquid levels were simulated in the study of meso-scale coherent structure with the new model.

9

Velocity inlet boundary was applied for gas phase, and no-slip boundary for mixture phase Degassing boundary was imposed at outlet for gas phase and slip velocity boundary for mixture phase. The gas velocities and gas volume fraction at the inlet at different superficial gas velocities are listed in Table 2. No-slip boundary was set for wall, and the standard wall function for the nearwall treatment. 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 NonOscillatory (WENO) scheme (Jiang and Shu, 1996), and the fourth-order central difference method was used to discretize diffusion terms. The source terms in the k-ε model equations were linearized to avoid the unphysical solutions. In this work, GPU programming for both LBM model and FVM was carried out to accelerate computation. LBM and FVM are linked by exchanging the following variables: velocity (from LBM to FVM), gas volume fraction (from FVM to LBM), turbulent viscosity (from FVM to LBM), pressure corrections (from FVM to LBM). The coupling between LBM and FVM was implemented at each time-step. The flow domain was divided into several sub-domains, each of which corresponds to one GPU card, and one host CPU core corresponds to the GPU. All the parallel parts of the computation were implemented on GPUs, whereas CPU was used to launch the kernel functions on GPU, execute the serials parts of the code and complete the communications between each sub-domain. The computation was implemented on NVIDIA P100 GPU cards.

4. Computation speed and model validation

10

4.1 GPU acceleration To assess the computational speed, comparisons were made between the new solver implemented on GPU cards and the commercial CFD solver (Ansys Fluent) running on CPUs. Grid resolution and time step were kept the same in both the cases. In the default grid resolution, the total number of grid cells is 288,000 with a time step 0.001 s. Two NVIDIA P100 GPU cards were used for the LBM-mixture simulations, and two CPU cores were assigned for the CFD simulation. Figure 2 compares the required computational time for two NVIDIA P100 GPUs or two CPU cores to update one-second physical time of flow. 25,877.5 seconds were needed for the CFD simulation with two CPU cores. By contrast, it only took about 5 seconds to achieve one-second physical time of flow with the LBM-mixture simulation accelerated by two GPUs, suggesting a remarkable acceleration of more than 5,000 times. In terms of MLUPS (Million Lattice Updates Per Second), which is usually utilized as an index to assess the performance of LBM algorithm (Shu and Yang, 2018a), the performance of the LBM algorithm in this work is about 57.6 MLPUS. This also suggests that the ratio between the computational time and physical evolution time could reach 5:1, which enables the long-time transient calculation of gas-liquid flows with affordable timeframe, even for large-scale industrial reactors. This is of practical significance for the dynamic analysis of the meso-scale coherent structure of gas-liquid flows in BCRs.

4.2 Effects of grid resolution

Meandering bubble plume observed in a rectangle BCR at lower superficial gas velocity (Buwa and Ranade, 2002) was reported to cause the oscillation of liquid velocity. The radial liquid velocity at the BCR center (x=0.1m, y=0.025m and z=0.225m) was recorded for 600 seconds in the

11

simulation. The effects of grid resolution on the radial liquid velocity are presented in Figure 3. There was no oscillation in all the three cases at the start of calculation. The predicted radial liquid velocities started to oscillate around the 100th second. However, the oscillation of the coarse grid simulation (Δx=5 mm) was damped before the 300th second, whereas the oscillations of the other two cases were kept. The damped oscillation in the coarse grid simulation might arise from numerical dissipation. In this work, the lateral force, e.g. lift force, is neglected and the lateral motion of bubble plumes was captured in our simulation at relatively high-resolution mesh (2.5 mm). This may be relevant to the more accurate prediction of pressure gradient and liquid lateral flow with a relatively high-resolution mesh, so that the lateral motion of bubble plumes could be captured without the need of additional lift. Actually lift is not required when the mesh is fine enough to reach the scale of direct numerical simulation (DNS). The transient axial liquid velocities from 400-600 seconds were averaged to obtain the timeaveraged axial velocity. Effects of grid resolution on the prediction of time-averaged axial liquid velocities are presented in Figure 4. Grid independence is achieved for the grid finer than 2.5 mm. Coarse grid calculation (Δx=5 mm) even fails to give a symmetrical prediction of axial liquid velocity along the radial direction. This indicates that the grid refinement is significant to the simulation of gas-liquid two-phase flow. In the following simulations, the grid sizes were kept as 2.5 mm.

4.3 Model validation Figure 5 illustrates the evolution of gas hold-up and liquid velocity. Compared to the wellknown 'gulf-stream' or 'cooling-tower' pattern of long-time averaged flow field, two staggered rows of vortices can be found with the oscillating bubble plume, which is qualitatively consistent with 12

the experiments of Becker et al. (1999). The oscillation features of the bubble plume proved to be induced by these vortexes or meso-scale coherent structures moving downwards, as shown in Figure 5 (Delnoij et al., 1997). From the 300th to 307th second, the middle part of plume moves from right to left. Meanwhile, the rotation direction of the main vortices in the upper part reverses. From the 307th to 314th second, the flow tends to restore the initial flow pattern at 300th second. The transient simulation indicates that the radial mixing or the exchange induced by convection occurs not only at the both ends of the BCR, e.g., the inlet zone and gas degassing zone, but also between these two ends due to the meso-scale coherent vortices. In contrast to the transient observation, the timeaveraged 'gulf-stream' or 'cooling-tower' pattern only suggests the radial macromixing or the liquid exchange at both ends. Compared to the steady-state simulation, the transient simulation of longer physical time has the potential to serve as a better approach in revealing the mixing mechanisms or predicting the mixing time of liquid phase. Then the transient flow properties are quantitatively estimated by the spectral analysis on the temporal variation of lateral liquid velocity at the BCR center. The distinct peak at 0.075 Hz, as shown in Figure 6, corresponds to the oscillation period 13.3 seconds, which is comparable with the experimental measurement of oscillation period 13.4 seconds (Becker et al., 1999). This suggests that the evolution of these meso-scale coherent structures is of low frequency, and a longer physical time simulation might be indispensable to obtain the time-independent time-averaged results. The time-averaged results are then compared with the experimental data in literature. The simulation is averaged over 200th-600th seconds. Figure 7 (a)-(c) compares the simulated liquid velocity at three different levels with the experiment of Pfleger et al. (1999). In general, the predicted liquid velocities are in good agreement with the experimental data. However, there are some 13

deviations between the simulation and the experiments near the wall and lower parts of the bubble column. It is difficult to claim that the deviation comes from the numerical errors or model simplification as no error bars found in the literature and asymmetrical profile of the experimental liquid velocity. The asymmetrical properties of liquid velocity profile might suggest that the sampling time in the reported experiments might not be sufficient. Figure 7 (d) compares the simulated gas hold-up at z=0.37m with the experiment of Buwa and Ranade (2003). The simulation of gas holdup is essentially pertinent to the drag model. Since we only use a constant slip velocity and a simplified drag model, as shown in Table 1, we only concentrate on the effect of gas plume on liquid mixing.

5.

Dynamics of meso-scale coherent structure The trajectories of the local vortex cores are further resolved and magnified in Figure 8 to

analyze the dynamics of meso-scale coherent structure in the BCRs. The details on the identification of the vortex core can be referred to Holmén (2012). The vortex cores first appear in the upper corners and then move downwards to the bottom along the axial direction in the form of alternate clockwise or counter-clockwise liquid rotation, reflecting the effects of bubble plume on liquid mixing. Meanwhile, the three snapshots also indicate that these vortex cores first migrate to the center and then alternately turn back to the wall along the radial direction. Figure 9 illustrates the trajectory of vortex cores and shows a well-developed zone from H=1.75W(0.35m) to 0.5W(0.10m) . This indicates that only the zones within 0.5W from the two ends can affect the radial migration of the vortex cores, i.e., the inlet zone and degassing zone. Furthermore, this also suggests two mechanisms, namely, the vortex rotation and revolution (radial migration) in the intensification of radial macromixing at the two ends. By contrast, only the vortex 14

rotation contributes to the radial macromixing in the well-developed zone. It should also be noticed that the vortex cores never cross the symmetrical axis of the column. Then the characteristic parameters of the descending vortex cores generated at the left upper corner, such as size, intensity, and descending velocity, are depicted in Figure 10. The quantification approaches of vortex size and intensity were covered in Schielicke et al. (2016) and Holmén (2012). The vortex size first increases as the vortex departures from the upper corner, and then keeps constant and decreases as it approaches the bottom corner. The variation trend of intensity and descending velocities is analogous to that of the vortex size. Then the effects of column aspect ratio and superficial gas velocity (scaleup effect) on the mean descending velocity of the vortex cores are illustrated in Figure 11-12. With the increase of aspect ratio from 2.25 to 5.7, the descending velocities of the vortex core increase. However, the double increase of aspect ratio does not generate a corresponding increase of the descending velocity of vortex cores. The descending velocities approaches a constant when the aspect ratio is over 5.7, which might be one of the reasons why the increase in BCR height could lead to larger mixing time or poor back-mixing. The superficial gas velocity causes the increase in the descending velocity of vortexes, as shown in Figure 12, which can further generate smaller mixing time or better backmixing. The resolution of mesoscale vortex structure demonstrates the advantage of the new GPUaccelerated LBM algorithm which allows a coarse grid through a tunable and smaller relaxation time τ, for fast simulation and future application in industrial processes. Of course, further improvement or advanced model could be achieved via, e.g., better closure of interphase momentum exchange for higher gas holdup and consideration of polydispersed bubble size population. A 15

constant slip velocity is used and the term τDM is neglected in this work for simplicity. The basic assumption in Sokolichin and Eigenberger (1999) is the relatively lower volume fraction of gas phase and the higher density ratio between gas and liquid. In this paper, the model is simplified to meet the needs of GPU acceleration since less iteration and good linearity in model algorithm is favorable for GPU-accelerated computation. More model complexity could be involved to eliminate the assumptions, but first these complexities should be first evaluated and then new algorithms be devised to coordinate the model complexity, numerical stability and required computational resources for practical large-scale applications.

6.

Conclusions

A GPU-accelerated LBM simulation was proposed to boost the transient simulation of multiphase flow in BCR. The proposed algorithm allows a tunable and smaller τ and therefore enables the simulation of air-water systems with a coarse grid at the level of 2.5 mm. In contrast, the grid size in literature has to be at the level of 0.006 mm or even smaller for numerical stability. The new algorithm can therefore greatly reduce the required computational resources by several orders of magnitude (up to 7) , enabling the simulation of relatively large-scale gas-liquid reactors. With the integration of GPU-acceleration and the new algorithm, about 5,000-folds acceleration can be achieved by only two GPU cards, in contrast to the traditional CFD solver implemented on two CPU cores. In this case, it only needs about 5 seconds of computational time to update one second of physical time. The new solver resolves the low-frequency dynamics of mesoscale coherent structure, and reveals two mechanisms, e.g., vortex rotation and vortex revolution (radial migration), for the intensification of radial macromixing at both the inlet and degassing zones. By contrast, only 16

the vortex rotation contributes to radial macromixing in the rest zones. The approach demonstrates the potential of this new approach to analyze the transient dynamics and the effect of mesoscale coherent structure on macromixing in scaleup or optimization of practical BCRs.

Acknowledgement Financial support from National Key Research and Development Program of China (2017YFB0602500) and National Natural Science Foundation of China (Grant Nos. 91834303, 91634203, 21808222) are acknowledged.

Notation B

Model Parameter

b

Number of discrete velocities

c1

Model parameter

c2

Model parameter

cg

Mass fraction of gas phase

ci

Discretized velocity vector

CD

Drag force coefficient

Cμ

Model parameter

DMg

Turbulent diffusivity

db

Bubble size

f

the velocity function vector

Fb

Body force term

g

Gravity acceleration constant

I

The unit matrix 17

Gk

Generation term of turbulent kinetic energy

meq

the equilibrium moment

m

the velocity function vector in moment space

M

The transfer matrix

k

Turbulent kinetic energy

p

Pressure term

p*

Pressure term with static pressure

r

Location vector

S

a diagonal matrix related the relaxation time

Sct

Turbulent Schmidt number

S

Magnitude of strain rate

t

Time

u+

Non-dimensional velocity

u

Local velocity

ulg

Slip velocity between liquid and gas phase

v0

Kinematic viscosity

vt

Kinematic turbulent viscosity

x

The spatial location vector

y+

Non-dimensional distance

Greek letters 18

α

Volume fraction

Δt

Time step

σk0

Model parameter

σk

Model parameter

σε0

Model parameter

σε

Model parameter

ε

Turbulent dissipation rate

ρ

Density

μ

Viscosity

κ

Karman constant

τDm

Diffusion stress due to the phase slip

τm

Viscous stress of mixture

τmt

Turbulent stress of mixture

τw

Shear stress at wall

τ

The non-dimensional relaxation time

Abbreviations BCR

Bubble Column Reactor

CFD

Computational Fluid Dynamics

DNS

Direct Numerical Simulation

EMMS

Energy Minimization Multi-Scale

FVM

Finite Volume Method 19

GPU

Graphic Processing Unit

LBM

Lattice Boltzmann Method

MRT

Multiple Relaxation Time

SRT

Single Relaxation Time

Subscripts g

Gas phase

lg

Relative

P

Node P

l

Liquid phase

m

Mixture phase

t

Turbulence

References Becker, S., De Bie, H., Sweeney, J., 1999. Dynamic flow behaviour in bubble columns. Chem. Eng. Sci. 54, 4929-4935. Buwa, V.V., Ranade, V.V., 2002. Dynamics of gas–liquid flow in a rectangular bubble column: experiments and single/multi-group CFD simulations. Chem. Eng. Sci. 57, 4715-4736. Buwa, V.V., Ranade, V.V., 2003. Mixing in bubble column reactors: Role of unsteady flow structures. Canadian Journal of Chemical Engineering 81, 402-411. Chen, P., Duduković, M.P., Sanyal, J., 2005. Three-dimensional simulation of bubble column flows with bubble coalescence and breakup. AIChE J. 51, 696-712. Chen, P., Sanyal, J., Dudukovic, M.P., 2004. CFD modeling of bubble columns flows: implementation of population balance. Chem. Eng. Sci. 59, 5201-5207. Deen, N.G., Kuipers, J., 2013. Direct numerical simulation of wall-to liquid heat transfer in dispersed gas–liquid two-phase flow using a volume of fluid approach. Chem. Eng. Sci. 102, 268-282. Delnoij, E., Kuipers, J.A.M., van Swaaij, W.P.M., 1997. Dynamic simulation of gas-liquid two-phase flow: effect of column aspect ratio on the flow structure. Chem. Eng. Sci. 52, 3759-3772. Delnoij, E., Westerweel, J., Deen, N.G., Kuipers, J.A.M., van Swaaij, W.P.M., 1999. Ensemble correlation PIV applied to bubble plumes rising in a bubble column. Chem. Eng. Sci. 54, 5159-5171. 20

Dijkhuizen, W., Roghair, I., Van Sint Annaland, M., Kuipers, J.A.M., 2010. DNS of gas bubbles behaviour using an improved 3D front tracking model—Drag force on isolated bubbles and comparison with experiments. Chem. Eng. Sci. 65, 1415-1426. Holmén, V., 2012. Methods for vortex identification. Master’s Theses in Mathematical Sciences. Icardi, M., Ronco, G., Marchisio, D.L., Labois, M., 2014. Efficient simulation of gas–liquid pipe flows using a generalized population balance equation coupled with the algebraic slip model. Appl. Math. Model. 38, 4277-4290. Inamuro, T., Ogata, T., Tajima, S., Konishi, N., 2004. A lattice Boltzmann method for incompressible twophase flows with large density differences. J. Comput. Phys. 198, 628-644. Jiang, G., Shu, C., 1996. Efficient implementation of weighted ENO schemes. J. Comput. Phys. 126, 202228. Joshi, J.B., 2001. Computational flow modelling and design of bubble column reactors. Chem. Eng. Sci. 56, 5893-5933. Kliphuis, A.M., de Winter, L., Vejrazka, C., Martens, D.E., Janssen, M., Wijffels, R.H., 2010. Photosynthetic efficiency of Chlorella sorokiniana in a turbulently mixed short light-path photobioreactor. Biotechnol. Prog. 26, 687-696. Laborde-Boutet, C., Larachi, F., Dromard, N., Delsart, O., Schweich, D., 2009. CFD simulation of bubble column flows: Investigations on turbulence models in RANS approach. Chem. Eng. Sci. 64, 4399-4413. Li, J., Ge, W., Wang, W., Yang, N., Liu, X., Wang, L., He, X., Wang, X., Wang, J., Kwauk, M., 2013. From Multiscale Modeling to Meso-Science: A Chemical Engineering Perspective Principles. Modeling, Simulation, and Application: Springer. Manninen, M., Taivassalo, V., Kallio, S., 1996. On the mixture model for multiphase flow. Masood, R.M., Jovicic, V., Delgado, A., 2015. Numerical Simulation of Interfacial Closures for 3D Bubble Column Flows. Chem. Eng. Technol. 38, 777-786. McClure, D.D., Norris, H., Kavanagh, J.M., Fletcher, D.F., Barton, G.W., 2014. Validation of a Computationally Efficient Computational Fluid Dynamics (CFD) Model for Industrial Bubble Column Bioreactors. Ind. Eng. Chem. Res. 53, 14526-14543. Pfleger, D., Becker, S., 2001. Modelling and simulation of the dynamic flow behaviour in a bubble column. Chem. Eng. Sci. 56, 1737-1747. Pfleger, D., Gomes, S., Gilbert, N., Wagner, H.G., 1999. Hydrodynamic simulations of laboratory scale bubble columns fundamental studies of the Eulerian–Eulerian modelling approach. Chem. Eng. Sci. 54, 5091-5099. Roghair, I., Lau, Y.M., Deen, N.G., Slagter, H.M., Baltussen, M.W., Annaland, M.V., Kuipers, J.A.M., 2011. On the drag force of bubbles in bubble swarms at intermediate and high Reynolds numbers. Chem. Eng. Sci. 66, 3204-3211. Sankaranarayanan, K., Shan, X., Kevrekidis, I.G., Sundaresan, S., 2002. Analysis of drag and virtual mass forces in bubbly suspensions using an implicit formulation of the lattice Boltzmann method. J. Fluid Mech. 452, 61-96. Sankaranarayanan, K., Sundaresan, S., 2008. Lattice Boltzmann Simulation of Two-Fluid Model Equations. Ind. Eng. Chem. Res. 47, 9165-9173. Sanyal, J., Vásquez, S., Roy, S., Dudukovic, M.P., 1999. Numerical simulation of gas–liquid dynamics in cylindrical bubble column reactors. Chem. Eng. Sci. 54, 5071-5083. Schielicke, L., Névir, P., Ulbrich, U., 2016. Kinematic vorticity number–a tool for estimating vortex sizes and circulations. Tellus A: Dynamic Meteorology Oceanography 68, 29464. 21

Shu, S., Vidal, D., Bertrand, F., Chaouki, J., 2019. Multiscale Multiphase Phenomena in Bubble Column Reactors: A Review. Renew. Energ. 141, 613-631. Shu, S., Yang, N., 2013. Direct numerical simulation of bubble dynamics using phase-field model and lattice Boltzmann method. Ind. Eng. Chem. Res. 52, 11391-11403. Shu, S., Yang, N., 2018a. GPU-accelerated large eddy simulation of stirred tanks. Chem. Eng. Sci. 181, 132-145. Shu, S., Yang, N., 2018b. Numerical study and acceleration of LBM-RANS simulation of turbulent flow. Chin. J. Chem. Eng. 26, 31-42. Silva, M.K., Mochi, V.T., Mori, M., d’Ávila, M.A., 2014. Experimental and 3D computational fluid dynamics simulation of a cylindrical bubble column in the heterogeneous regime. Ind. Eng. Chem. Res. 53, 3353-3362. Sokolichin, A., Eigenberger, G., 1999. Applicability of the standard k–ε turbulence model to the dynamic simulation of bubble columns: Part I. Detailed numerical simulations. Chem. Eng. Sci. 54, 2273-2284. Sokolichin, A., Eigenberger, G., Lapin, A., Lübert, A., 1997. Dynamic numerical simulation of gas-liquid two-phase flows Euler/Euler versus Euler/Lagrange. Chem. Eng. Sci. 52, 611-626. Song, F.F., Wang, W., Li, J.H., 2013. A lattice Boltzmann method for particle-fluid two-phase flow. Chem. Eng. Sci. 102, 442-450. Swiderski, K., Narayanan, C., Lakehal, D., 2016. Application of N-phase algebraic slip model and direct quadrature method of moments to the simulation of air-water flow in vertical risers and bubble column reactor. Comput. Chem. Eng. 90, 151-160. Tabib, M.V., Roy, S.A., Joshi, J.B., 2008. CFD simulation of bubble column - An analysis of interphase forces and turbulence models. Chem. Eng. J. 139, 589-614. Wang, T., Wang, J., 2005. Two-fluid model based on the lattice Boltzmann equation. Phys. Rev. E 71, 045301. Witz, C., Treffer, D., Hardiman, T., Khinast, J., 2016. Local gas holdup simulation and validation of industrial-scale aerated bioreactors. Chem. Eng. Sci 152, 636-648. Xiao, Q., Yang, N., Li, J., 2013. Stability-constrained multi-fluid CFD models for gas-liquid flow in bubble columns. Chem. Eng. Sci. 100, 279-292. Xue, J., Chen, F., Yang, N., Ge, W., 2017a. Eulerian–Lagrangian simulation of bubble coalescence in bubbly flow using the spring-dashpot model. Chin. J. Chem. Eng. 25, 249-256. Xue, J., Chen, F., Yang, N., Ge, W., 2017b. A Study of the Soft-Sphere Model in Eulerian-Lagrangian Simulation of Gas-Liquid Flow. Int. J. Chem. React. Eng. 15. Yang, N., Wu, Z.Y., Chen, J.H., Wang, Y.H., Li, J.H., 2011. Multi-scale analysis of gas–liquid interaction and CFD simulation of gas–liquid flow in bubble columns. Chem. Eng. Sci. 66, 3212-3222. Yuan, P., Schaefer, L., 2006. Equations of state in a lattice Boltzmann model. Phys. Fluids 18, 042101.

22

Tables Table 1. Governing equations of the mixture model. Model formula For the mixture phase: Continuity equations:

Momentum

 m      m um   0 , ρm=αlρl+αgρg, um= (αlρlul+αgρgug)/ρm t

  mum 

equations:

t







t   mumum   p   m  m  Dm  mg ,  m  m 2 um ,



 mt   t  2 um ,  Dm  mcg 1  cg ulg ulg , cg= αgρg/ρm For the gas phase:

Continuity equations:

g t









u lg u lg 

4db 3C d

 g um  DMgg   g 1 cg ulg  , DMg=vt/Sct, vt= μt/ρm, Sct=1  

Momentum equations: ug  um 

(1g ) l

m

ulg , ulg=constant or

23



g

 l 

l

g

Table 2. Fluid properties and numerical steps. Fluid properties Property(unit)

Gas

Liquid

Density(kg/m3)

1.225

1000

Dynamic viscosity(kg/m.s)

1.8*10-5

1*10-3

Numerical steps Ug (cm/s)

Gas inlet velocity(cm/s)

Gas inlet volume fraction

0.14

9.96

0.50

0.16

11.38

0.50

0.3

15.24

0.70

0.59

29.97

0.70

24

Figures

Figure 1. Geometry of the centrally-aerated rectangle bubble column (Case 1) (unit: mm).

25

Computational time for one second physical time (s)

105 104 103 102 101 100

Figure 2.

2 NVIDIIA P100 GPU 2 Intel Xeon E5-2680-V2 CPU cores (Commercial software) card (LBM-mixture)

Comparison of GPU-accelerated LBM computation and CPU-based CFD simulation.

26

0.3

5 mm 2.5 mm 1.25 mm

Radial velocity (m/s)

0.2 0.1 0.0 -0.1 -0.2 -0.3

0

200

400

600

Time (s)

Figure 3. Effect of grid resolution (5 mm, 2.5 mm and 1.25 mm) on radial velocity at the BCE center (Ug=0.14 cm/s).

27

Axial liquid velocity (m/s)

0.2

0.1

5 mm 2.5 mm 1.25 mm

0.0

-0.1

-0.2 0.00

0.05

0.10

0.15

0.20

X(m)

Figure 4. Effect of grid size (Δx=1.25, 2.5, and 5 mm) on the time-averaged of axial liquid velocity predictions (superficial gas velocity 0.14 cm/s).

28

Figure 5. Evolution of gas hold-up and liquid flow field (y=D/2).

29

8

Power spectral density

fp=0.075 Hz 6

4

2

0 0.00

0.02

0.04

0.06

0.08

0.10

Frequency (Hz)

Figure 6. Spectral analysis of lateral liquid velocity component (x=W/2, y=D/2, z=H/2).

30

0.2

Experiment (Pfleger et al., 1999) z=0.37m Simulation

Axial liquid velocity (m/s)

Axial liquid velocity (m/s)

0.2

0.1

0.0

-0.1

(a) -0.2 0.00

0.05

0.10

0.15

0.1

0.0

-0.1

(b) -0.2 0.00

0.20

Experiment (Pfleger et al., 1999) z=0.25m Simulation

0.05

X(m)

0.020

Experiment (Pfleger et al., 1999) z=0.13m Simulation

Time Averaged gas hold-up

Axial liquid velocity (m/s)

0.2

0.1

0.0

-0.1

(c) -0.2 0.00

0.05

0.10

0.10

0.15

0.20

X(m)

0.15

0.015

0.010

0.005

(d) 0.000 0.00

0.20

Experiment (Buwa et al., 2001) Simulation

0.05

0.10

0.15

X(m)

X(m)

Figure 7. Comparison of time-averaged simulation with experiments at different levels.

31

0.20

Figure 8. Evolution of the coherent vortex structures.

32

Figure 9. Trajectories of vortex cores (solid symbols/solid line for the vortex on the right; empty symbols/dashed line for vortex on the left).

33

Figure 10. Characteristic parameters of vortex cores at the right-hand side during vortex descending: (a). vortex size; (b). vortex intensity; (c) vortex descending velocity.

34

Figure 11. Effect of aspect ratio (H/W) on the mean vortex descending velocity.

35

Figure 12. Effect of superficial gas velocity on the mean vortex descending velocity.

36

Highlights



A GPU-accelerated LBM solver for transient gas-liquid flow simulation is developed.



The proposed LBM allows a free tuning of dimensionless relaxation time.



The memory requirement can be reduced by up to 7 orders of magnitude.



About 5,000-folds acceleration for a lab-scale bubble column is achieved.



Effect of mesoscale vortex structure in macromixing intensification is resolved.

37

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

38