Numerical simulation of a turbulent pipe flow: FluidX3D LBM validation

. The validation of the computational LBM code FliudX3D is presented on the example of turbulent flow in a pipe at two Reynolds numbers: 5300 and 37700, built on the bulk velocity, pipe diameter and kinematic viscosity. Due to the LBM approach, the code performance allows massive calculations to be performed in a short period of time with a good agreement with the literature data for the lower Reynolds number. However, the lack of the possibility to refine the computational grid leads to insufficient resolution of the turbulent boundary layer for the higher Reynolds number.


Introduction
Over the past three decades, the Lattice Boltzmann method (LBM) has gained significance in the area of computational fluid dynamics (CFD) [1].Known for its outstanding adaptability and versatility, this method stands as an effective tool for simulating a wide spectrum of processes.Its applications cover various areas, such as immiscible multicomponent liquid flows [2], phase transitions [3], heat exchange and thermal transfer phenomena [4], turbulent flows [5] and others.A set of lattice Boltzmann equations describes the evolution of particles population moving in different directions.These equations describe interaction of populations within each cell and their interaction with an external volume force.These interactions are typically localized within one computational cell.To account for interactions between neighboring cells in the LBM, a linear, non-local transfer operator is employed.This natural locality of interactions in LBM provides an opportunity for high-performance parallel computing.Therefore, LBM is highly competitive in terms of speed compared to traditional discretization-based methods, such as finite difference and finite volume methods, which are rooted in conventional CFD principles [6][7][8][9][10].At the moment numerous successful implementations of LBM exist, with one of the most widely used being the open-source software package Palabos [11].However, despite its popularity, Palabos cannot provide a universal solution for cases of multiphase and multicomponent interactions.Furthermore, this code lacks support for GPU acceleration, which results in longer computation times, thus limiting the possibility to speed up calculations.
Currently, extensive efforts are being made to harness the computational powers of GPUs for LBM-based simulations.A noteworthy example is the open-source code FluidX3D [12].This code contains a basic set of LBM algorithms required for simulating fluid flow.However, FluidX3D achieves exceptional memory efficiency by employing a user-defined data structure in GPU memory.Furthermore, it achieves rapid calculations and low memory requirements through the innovative "Esoteric-Pull" [13] algorithm, for population transfer between neighboring cells.
Nevertheless, one of the primary challenges in LBM simulations is the need to implement specialized boundary-processing techniques, especially in scenarios involving curved or moving boundaries [14][15][16][17].This need arises from the limitation of LBM spatial cells to have any form other than cubic.The FluidX3D code does not incorporate these boundary algorithms and instead approximates all solid surfaces using a voxel representation.Nevertheless, the code's capacity to perform calculations with tens of billions of spatial cells calls into question the need for additional edgesmoothing algorithms.
It has been previously demonstrated that FluidX3D produces reliable results in laminar flow scenarios.This has been evidenced by comparisons with results obtained using the spectral elements method.In this paper, we extend our efforts to validate the code's performance by simulating turbulent flow within a cylindrical pipe.The simulation results obtained were compared to those obtained using the spectral element code Nek5000, as referenced in [18].This comparative analysis demonstrated an agreement in the results, particularly for scenarios characterized by high grid resolution and low Reynolds numbers.

Governing equations
In the LBM the observation is carried out over the distribution function (DF) of particles in the phase space , which evolves according to the Boltzmann equation: where � � is the particle velocity, F is the external body force, � is the density, and ���� is the collision operator.The coordinate space is divided into equal square or cubic cells, and a set of vectors �� � � � directed to the neighboring cells is selected from the velocity space.In general, these sets are described by the formula DnQm, where n corresponds to the number of the space dimensions, and m corresponds to the number of vectors in the set �� � � �.In calculations the observation is carried out over the values of the DF for the given vectors � � �� �, ���, called populations.In general, without taking into account external forces, populations are transformed according to the discrete Boltzmann equation ( 2) and the fluid density and velocity are restored as follows: 2) determines the change in populations in a cell, and then their movement to neighboring cells with the corresponding velocities � � � to the point � � � � � � �� which populations reach in time � � ��.
Algorithmically, this equation is solved in two stages: 1. Transformation of populations under the action of the collision operator within one cell -the collision step: where � � * -temporary values of populations; 2. Transfer of temporary populations to the neighboring cells according to directions of their vectors � � � is the streaming step: The most common form of the collision operator is the Bhatnagar-Gross-Krook operator (BGK).It is determined by the following formula: where � � �� are equilibrium populations, and � is the relaxation frequency associated with the viscosity according to the formula To determine equilibrium populations, it is considered that, on microscopic scales, the equilibrium velocity distribution function of particles has the form of a Maxwellian distribution (up to a shift by the total flow velocity).Further, this function is expanded in terms of Hermite polynomials, limited to the first three terms of the expansion.The final formula for equilibrium populations takes the following form: where � � are weight coefficients determined by a set of velocities and � � are components of the flow velocity at a given point in space.
When external forces are taken into account, equation ( 2) is transformed to the form: where � is the source of force.It is related to force populations as follows: where � � are the force populations.There are several approaches to defining � � , the most popular of which is the Guo approach: where � � are the components of the density of external forces.The flow velocity in this case is restored taking into account the forces: As a result, the complete time step algorithm in LBM consists of the following elements: 1. Calculation of force density based on the conditions of the problem 2. Recovery of moments according to (3) and ( 13) 3. Calculation of equilibrium populations according to (7) 4. Calculation of force sources according to (11) 5. Performing collision step according to (10) 6. Streaming step execution (5)

Returning to the first element
The processing of boundary conditions in LBM in the simplest case is performed using a voxel approximation of the boundary surface.In this case, the cells of the computational domain are divided into three types: solid media cells, boundary cells, and fluid cells.
In fluid cells, the flow step is performed according to the algorithm described above.In boundary cells, various methods are used to determine the populations reflected from the walls.In the simplest algorithm called Bounce Back (which is used in FluidX3D), at the streaming step, provided that the neighboring cell in the direction � � � corresponds to a solid medium, populations are reflected from the wall.

LES model
It is known that direct numerical simulation in LBM has a very low stability.To provide more stable simulations, a collision matrix operator is used or turbulence models are added.In FluidX3D software the Smagorinsky-Lilly LES turbulence model is implemented [19].In this model the additional stress viscosity is being calculated from the nonequilibrium stress tensor: where � � is the velocity vector corresponding to � � .After that the second variance of this tensor is being calculated: � � � � �,� � � �,� and then the intensity of the local filtered strain rate tensor for Navier-Stokes equation is calculated as follows: where � � is kinematic viscosity, � � is the Smagorinsky constant and � is the lattice step, which is commonly set to 1 in the LBM.And finally, total viscosity is calculated by adding stress viscosity: The obtained total viscosity is connected with relaxation rate in collision operator (6) the following:

Results
The simulation of the flow was performed for two different Reynolds numbers 5300 and 37700, defined as where D is the pipe diameter, � ���� is the average axial velocity and � is the kinematic viscosity.The length of the pipe was set 5 times larger than the diameter as shown at Fig. 1.The average axial velocity component was maintained equal to � ���� = 0.05 in the LBM units by adjusting external force directed along the pipe axis.The boundary conditions were set periodic for the inlet and outlet.Overall number of timesteps for averaging velocity fields was set to satisfy the condition: � ����� � ���� � ���� ⁄ .For the largest simulation with the total 2��� � 1� � spatial cells and ���� � 1� � timesteps overall simulation time was 25h.Simulations were performed on three NVIDIA A100 80Gb GPUs.
The velocity control module was added on top of FluidX3D software utilizing OpenMP library.Also, the additional GPU kernel and data field in GPU memory were added to collect the average velocity fields.The postprocessing of .vtkfiles was done with the help of python vtk library.Firstly, the DxDx5D numpy array was averaged over the z-axis.Then, the obtained DxD array was interpolated on the polar grid with scipy.interpolate.griddatatool and averaged over angle.Normalization of velocity and distance from the wall was performed over the dynamic velocity � � according to (19.1) - (19.3).
For the finest resolution of 800�800�3999 cells the plot of � � versus �1 � �� � for the �� � ����� shows an excellent agreement with the reference results (see Fig. 2).But for the �� � ������ FluidX3D significantly differ from the reference results.Due to the high Reynold number a very small boundary layer is affecting turbulence formation.Therefore, the deviation is caused by the poor resolution of the boundary layer.Comparison between reference data [18] and FluidX3D simulation for the Reynolds numbers 5300 and 37700.The requirement for a uniform distribution of cells affects the ability to resolve a turbulent boundary layer as it is shown at Fig. 3, where the average velocity was performed for different numbers of cells in diameter.It can be noticed that with a higher spatial resolution, the obtained velocity profile is getting closer to the literature data

Conclusions
Validation of the computational LBM code FluidX3D was carried out using the example of turbulent pipe flow at two Reynolds numbers: 5300 and 37700.For a smaller one, a good correspondence between the results and the literature data was obtained.At the same time, it is shown that with an increase in the Reynolds number, a problem arises associated with insufficient spatial resolution in the region of the turbulent boundary layer.This is primarily due to the fact that in this implementation of the code there is no possibility of refining the computational grid.At the same time, the performance of the code still makes it promising for computational fluid dynamics after expanding its functionality.
Numerical simulations and data analysis were conducted within the Russian Science Foundation grant no.21-15-00091.The development of the numerical code was supported by the state contract with IT SB RAS.The computational resources were provided by the Cascade Supercomputer of the Institute of Thermophysics SB RAS.

Fig. 1 .
Fig. 1. 3D visualization of vorticity isosurfaces and instantaneous velocity and vorticity fields in a cross section.Velocity and vorticity fields are normalized by � ���� and D.

Fig. 3 .
Fig. 3. Convergence of normalized velocity with the increase of the number of computational cells in pipe diameter.