Wall-modeled large-eddy simulation integrated with synthetic turbulence generator for mutiple-relaxation-time lattice Boltzmann method

The synthetic turbulence generator (STG) lies at the interface of the Reynolds averaged Navier – Stokes (RANS) simulation and large-eddy simulation (LES). This paper presents an STG for the multiple-relaxation-time lattice Boltzmann method (LBM) framework at high friction Reynolds numbers, with consideration of near-wall modeling. The Reichardt wall law, in combination with a force-based method, is used to model the near-wall field. The STG wall-modeled LES results are compared with turbulent channel flow simulations at Re s ¼ 1000 ; 2000 ; 5200 at different resolutions. The results demonstrate good agreement with direct numerical simulation, with the adaptation length of 6 – 8 boundary layer thickness. This method has a wide range of potentials for hybrid RANS/LES-LBM related applications at high friction Reynolds numbers. V


I. INTRODUCTION
In industrial applications related to high Reynolds number flows, wall-bounded turbulence plays a crucial role in designing aircraft, cars, wind farms, and so on. 1 Direct numerical simulations (DNS) can accurately quantify the physics by solving the Navier-Stokes equations using high-order numerical approximations and grid refinement techniques. 2,3However, DNS is computationally too expensive for realworld applications, making it impractical for use in the design cycle.5][6] Nevertheless, LES still requires fine grid resolution in the near-wall region, requiring grid points to be proportional to OðRe n Þ with n ¼ 13/7. 2,3,7Wall-resolved LES (WRLES) is still far from an engineering tool.
To improve computational efficiency, studies have attempted to model the near-wall region by solving Reynolds-averaged Navier-Stokes (RANS) equations and using LES in the far field or by modeling the near-wall region with relatively few grids by reconstructing nearwall velocities.In past decades, the lattice Boltzmann method (LBM) has gained popularity for simulating fluid dynamic problems at a variety of scales, from micro-nanoscales [8][9][10][11][12] to macroscopic scales [13][14][15][16][17] at low Mach numbers.9][20] LBM's parallelization-friendly nature makes it attractive due to the local update of discrete particle distribution functions.Hou et al. 13 introduced the effective turbulent viscosity to model the subgrid-scale (SGS) turbulence in the LBM framework.This approach combines the advantages of LES techniques with the computational efficiency of LBM.
To further improve computational efficiency, the hybrid RANS-LES approach is becoming increasingly popular as it balances accuracy and computational cost.The RANS method is used to conserve macroscopic quantities in computationally less-demanding regions, while LES provides detailed flow information in computationally-intensive areas.Generating high-quality turbulence at the RANS/LES interface is critical for achieving accurate results, as highlighted by Wu. 21In the conventional computational fluid dynamics (CFD), studies have conducted wide range approaches by precursory DNS/LES data, 22 velocity field "recycling," [23][24][25] synthetic eddy method (SEM), [26][27][28] or involving control techniques 29,30 to generate turbulence.Most of the existing studies suffer from either relatively long turbulent-developed adaptation length, high computational cost, or hard to generalize for complex geometries.Shur et al. 31 successfully developed a synthetic turbulent generator (STG) for the detached eddy simulation (DES). 32The STG method creates velocity fluctuations based on the Fourier coefficients, which are given by the energy spectrum, whereas, for the SEM method, coherent structures are simulated by superimposing artificial eddies at the inlet plane.The results show both fast and robust with an adaption length of 2-4 boundary-layer widths.
Despite being an active area in the conventional CFD, studies using a turbulent generator in the LBM framework are still rare.Koda and Lien 33 generated turbulence for the channel flow by placing a "recycled" channel flow before the inlet, which requires a presimulated periodic turbulent channel flow.Nakayama et al. 34 used the similar recycling approach and additionally added a heat flux source from the realistic observation to simulate atmospheric flow over a city.Asmuth et al. 35 applied pre-synthetic turbulence inflow for the wind farm simulation.The above-mentioned methods need additional preparation work for the inlet turbulence before applying it to the applications.Buffa et al. 36 reconstructed turbulence by using the SEM method in LBM (SEM-LBM); however, this method may suffer from a relatively long adaptation length. 37Xue et al. 38 integrated the synthetic turbulent generator (STG) in the LBM framework at Re s ¼ 180 with wall-resolved LES-LBM simulation with an adaptation length of 2-4 boundary-layer widths.However, the friction Reynolds number is relatively low and the Bhatnagar-Gross-Krook (BGK) collision operator limited its applications for high Reynolds number cases.To reduce computational time and maintain STG accuracy, modeling the nearwall flow field with a wall-model (WM) in the LES-based LBM (WMLBM) is a non-straightforward task.The first WMLBM was proposed by Malaspinas and Sagaut 39 where they successfully reconstructed the first-layer near-wall velocity with the Musker wall function or log-law (see Ref. 40).2][43][44] However, the force-based method is rarely mentioned or described in detail in the LBM framework. 45n the present work, a synthetic turbulent generator is developed, integrated with a multiple-relaxation-time (MRT) LBM collision operator to tackle high friction Reynolds numbers.The near-wall region is modeled via a force-based wall model, using a wall law by Reichardt (see Ref. 46).The performance of the STG method is examined at three different friction Reynolds numbers: Re s ¼ 1000; Re s ¼ 2000, and Re s ¼ 5200 at various resolutions, and compared with DNS data from Hoyas and Lee (see Refs. 47 and 48).The proposed framework aims to produce high-quality turbulence at the RANS/LES-LBM interface to cater to high Reynolds number flow applications. 1,49,50

II. METHODOLOGY A. The multiple-relaxation-time lattice Boltzmann method
In this work, we utilize a three-dimensional (3D) lattice model with 19 discretized directions known as the D3Q19 model.The lattice cell is located at position x and time t, with a discretized velocity set c i for i 2 0; 1; The weight for the discretized directions is defined as The evolution equation for the distribution functions, accounting for collision and forcing, can be expressed as where X is a collision kernel and Dt is the lattice Boltzmann time step, which is set to unity.In this work, MRT collision kernel is chosen due to its higher numerical stability compared to the BGK model at high Reynolds numbers 51 where M is the transformation matrix from the population space to the moment space obtained via the Gram-Schmidt approach [M matrix is described in Eq. (A1)].S is the diagonal matrix with relaxation frequencies at different moments, S ¼ diagfx 0 ; x 1 ; …; x QÀ1 g, MRT collision operator will be equivalent to BGK with x i set to the same value x.The frequency x i is the inverse of the relaxation time s i .Note that we set s k ¼ s 9 ¼ s 11 ¼ s 13 ¼ s 14 ¼ s 15 , which are related to the kinematic viscosity , which is with c s being the speed of the sound, and c 2 s is equal to 1/3 lattice Boltzmann unit (LBU).Other relaxation parameters can be found in Eqs.(A2)-(A6).Instead of colliding in the population space, the MRT collides in the moment space; thus, Eq. (3) can be rewritten as where m is the moment space component, which is defined as mðx; tÞ ¼ Mfðx; tÞ: The moment space equilibrium m eq ðx; tÞ can be defined as Fðx; tÞ in Eq. ( 6) is the vector of F i ðx; tÞ, which is the volume force acting on the fluid cell 52 where g is the volume acceleration.The macro-scale quantities for the density, momentum, and momentum flux tensors can be calculated from the distribution function, the discrete velocities, and the volume force Note that the momentum flux, Pðx; tÞ, can be presented by the sum of the equilibrium and non-equilibrium parts, Pðx; tÞ ¼ P eq ðx; tÞ þP neq ðx; tÞ.

B. Smagorinsky subgrid-scale modeling
In this part, the lattice-Boltzmann-based Smagorinsky SGS largeeddy simulation techniques are summarized.In the LBM framework, the effective viscosity eff 3,13,33 is modeled as the sum of the molecular viscosity 0 and the turbulent viscosity where j Sj is the filtered strain rate tensor where C smag is the Smagorinsky constant, D represents the filter size, and s i is the relaxation time for the moment-space collision.
, with P neq being the non-equilibrium part of the momentum flux tensor shown in Eq. ( 12).With the help of Eq. ( 5), the total relaxation time s eff is obtained as Finally, s eff i is replaced by the related MRT collision operator relaxation s i to enclose the lattice-Boltzmann-based LES system.

C. Near-wall modeling for LBM
There have been various approaches in the LBM framework to model the near-wall field with a wall model.43][44]53 Reference 39 also coupled a RANS solver at the near-wall region with LBM; however, it is proven to be more time consuming than a monolithic LBM method.In the present work, we focus on developing a wall model that is based on forces in the nearwall region.
In this paper, the Reichardt wall law 46 will be used instead of the one from Musker. 40Similar to the Musker law, the Reichardt wall law covers the buffer zone and gives the full prediction of the u þ as a function of y þ , where the other wall law needs different expressions at different zones.Detailed discussions regarding wall function prediction errors can be found at Haussmann et al. 41 The Reichardt wall law shown in Fig. 1 is defined as where u þ and y þ are the dimensionless unit defined as where u s is the shear velocity, u s ¼ ffiffiffiffiffiffiffiffiffiffi s w =q p ; hÁi denotes ensemble average over space or time, s w is the wall shear stress, and u is the streamwise velocity.Figure 2 illustrates the use of a slip wall boundary condition for the wall treatment.Notice that, the force-based model was also mentioned in Kuwata and Suga 45 as the "wall function with specular reflection conditions."To generalize the wall model, a base vector e x is first needed to project the near-wall velocity u w on the wall-parallel direction where u 2 is the velocity in the second cell from the wall and n is the wall-normal vector.Then, we project the velocity in the first cell near the wall, u w , to obtain the scalar streamwise velocity, ûw , which is defined as Next, the aim is to compute the friction velocity u s by solving for u s in Eq. ( 16) By using the Newton method, we update the friction velocity locally.Note that, instead of using a plane averaged friction velocity, 41 this work uses the local value to make the algorithm more generalized.The wall shear stress can be estimated by s w ðx; tÞ ¼ u 2 s ðx; tÞqðx; tÞ.Finally, the force near the wall is defined as where F is the shear force acting on the wall.A notable advantage of the force-based method is that it does not require any reconstruction of populations to find the target velocity and density near the wall, making it easier to implement compared to other wall modeling methods.

III. SYNTHETIC TURBULENCE GENERATOR FOR THE WALL-MODELED LATTICE BOLTZMANN METHOD A. Synthetic turbulence generator formulation in LBM framework
In this study, a synthetic turbulence generator (STG) is positioned at the inlet of a channel flow.It needs velocity field from a k À x RANS simulation. 54The total velocity u in ðx; tÞ at the inlet is given by u in ðx; tÞ ¼ u RANS ðxÞ þ u 0 ðx; tÞ; (22)   where u RANS is the velocity vector obtained from a RANS simulation; then, the interpolated velocity will be applied on the LBM grid in the case of grid resolution differences. 38The STG generates the velocity fluctuations u 0 ðx; tÞ at the cell x at time t u 0 ðx; tÞ ¼ a ab v 0 ðx; tÞ: The time-averaged velocity fluctuation is zero, i.e., hu 0 ðx; tÞi ¼ 0. The Cholesky decomposition of the Reynolds stress tensor reads where R ab ¼ hu 0 a u 0 b i is taken from the Reynolds stress tensor using EARSM 55 when post-processing the 1D RANS data.v 0 ðx; tÞ in Eq. ( 23) is imposed by N Fourier modes given by where q n is the amplitude of a modified von Karman spectrum, n is the mode number, k n is the amplitude of the mode direction vectors d n with r n Á d n ¼ 0; / n is the random mode phase that is uniformly distributed in the interval of ½0; 2pÞ.A detailed description is found in Refs.31 and 38.The distribution function at the inlet of the channel flow can be defined as the sum of equilibrium part and the nonequilibrium part where f inðeqÞ i ðx; tÞ is the equilibrium part of the inlet distribution function, which can be calculated by FIG. 2. Sketch on the first layer near the wall.
ALGORITHM 1. STG-WMLBM: Implementation of synthetic turbulence generator for the MRT wall-modeled LBM.
1. Obtain the RANS velocity at the inlet u RANS ðxÞ in Eq. ( 22).2a.Read saved value from the RANS simulation: R ab , k, x field, etc. 2b.Compute the Reynolds stresses using EARSM 3. Calculate a ab with the help of Eq. ( 24). for all t from 0 to t end do for all cells do if cell x is at the RANS/LBM inlet then 4. Calculate v 0 ðx; tÞ from Eq. ( 25). 5. Compute the fluctuating velocity u 0 ðx; tÞ following Eqs.( 23) and ( 24), respectively.6. Compute boundary density q in ðx; tÞ thanks to Eq. (30).7. Reconstruct the particle's probability distribution function by combining Eqs. ( 26)-( 28).8. Update the LES-LBM relaxation time for the MRT framework t eff in Eq. ( 15) and replace in kinematic-viscosity-related relaxation time in Eq. ( 5).end if if cell x is at the wall function cell then 9. Compute u 2 .10. Compute the wall-parallel base vector e x with the help of Eq. ( 18).11.Compute the scalar streamwise velocity ûw using Eq.(19).12. Solve implicit function to obtain u s with the help of Eqs. ( 20) and (16).13.Compute s w with the help of s w ðx; tÞ ¼ u 2 s ðx; tÞqðx; tÞ.14.Update force on the cell F w ðx; tÞ by using Eq. ( 21) end if 15. Apply stream and collide with the consideration of forces to update the f i ðx; tÞ at each cell Eq.(3) end for end for Following the regularized scheme, 56 the non-equilibrium part of the inlet distribution function in Eq. ( 26) is obtained via where Q i ¼ c i c i À c 2 s I with I being the identity matrix.P in neq is the non-equilibrium part of the moment flux tensor, which is defined as The unknown variables in the ith direction at the inlet can be calculated via the known direction following Q i ¼ Q invðiÞ ; f inðneqÞ i ðx; tÞ ¼ f inðneqÞ invðiÞ ðx; tÞ, where the notation "inv" denotes the opposite direction of the unknown variable.For the density of the inlet boundary, we follow the idea from 57 q in ðx; tÞ ¼ 1 1 þ ûin ðx; tÞ ð2q ?ðx; tÞ þ q k ðx; tÞÞ; (30)   where ûin is the cross product with the normal unit vector n at the boundary ûin ¼ u in LB Á nðjû in j < 0:3c s Þ and u in LB is the velocity of the lattice Boltzmann domain at the interface.q ?and q k are the density calculated by where n 0 is the normal vector pointing toward the inlet boundary, f ?

B. Implementation summary
In Algorithm 1, we summarize the implementation details of STG in the WMLBM framework.

IV. TURBULENT CHANNEL FLOW SIMULATIONS
This section presents the turbulent channel flow simulations with the STG as the inlet at Re s ¼ 1000; Re s ¼ 2000, and Re s ¼ 5200.The present work employs three different resolutions, LBU, in reference to the height of the channel H, that is, LBU ¼ H=20; H=40, and H=60.

A. Numerical setup
The turbulent channel flow simulations use STG at the inlet and a pressure-free in the streamwise direction at the outlet.Periodic boundary conditions are employed in the spanwise direction (z) and wall functions are used at the first cell-layer near the top and bottom (y) planes, which are equipped with the slip boundary condition.To minimize the reflection wave's impact on the flow field, a sponge zone is placed near the outlet. 38The numerical setup for the channel flow is depicted in Fig. 3, where the boundary layer thickness (d) is defined as half of the channel height.The extent of the simulation domain is 20d Â 2d Â 1:6d in the x, y, and z directions, respectively.The sponge layer thickness is set to 1d.The Smagorinsky constant is set to C smag ¼ 0:01.The simulations are conducted at different friction Reynolds numbers and resolutions, and run for a total of 10 domainthrough times (10T).It is worth mentioning that the time convergence study of the STG has been addressed in our previous work (see Ref.  38).The present study finds that the turbulent statistics is stabilized after two domain-through times (2T).Thus, all statistical analyses begin after 2T.

B. Results
This work presents validation of the STG framework for high friction Reynolds numbers, i.e., Re s ¼ 1000; Re s ¼ 2000, and Re s ¼ 5200.The initial investigation focuses on the resolution with H ¼ 2d ¼ 20 LBU.The y þ values at H ¼ 20 LBU for Re s ¼ 1000, 2000, and 5200 are approximately 50, 100, and 260, respectively.The results, triggered by the STG inlet, are compared with DNS data. 47,48igure 4 shows the mean velocity field of u þ as a function of y þ .The STG-WMLBM results compare with the DNS data at different friction Reynolds numbers.The results of the STG-WMLBM at all three friction Reynolds numbers show good agreement with DNS reference starting from x=d ¼ 0. While the initial results are promising, further investigations are needed to analyze the Reynolds stresses in order to fully evaluate the effectiveness of the method.Figure 5 depicts the representation of hu 0 u 0 i þ as a function of y=d for different values of the friction Reynolds number (Re s ), namely, Re s ¼ 1000; 2000; 5200.The results obtained from STG-WMLBM have been compared against the DNS data to validate the accuracy of the LBM approach.At the initial stages of the flow, the LBM results display substantial deviations from the DNS data between the locations x=d ¼ 0 and x=d ¼ 4, since the turbulence in the flow has not yet fully developed.However, as the flow progresses downstream, these discrepancies reduce and the LBM results converge toward the DNS reference.By x=d ¼ 8, the LBM results have stabilized and show good agreement with the DNS data, demonstrating the effectiveness of the LBM approach in predicting the turbulent flow characteristics.
5. hu 0 u 0 i þ as a function of y=d at Re s ¼ 1000; 2000; 5200 for H ¼ 20 (LBU).Then, this work also examines the convergence of the algorithm by investigation on the comparison of hu 0 v 0 i þ between STG-WMLBM and DNS. Figure 6 shows hu 0 v 0 i þ as a function of y=d.The STG-WMLBM results compare with the DNS data at different friction Reynolds numbers.For Re s ¼ 2000 and 5200, the STG-WMLBM results exhibit excellent agreement with the DNS reference data after x=d ¼ 6 À 8.Meanwhile, for Re s ¼ 1000, the results converge to the DNS data at approximately x=d ¼ 8-10.However, discrepancies with the DNS data are observed near the wall, which can be attributed to the poor resolution.

Physics of Fluids
To further analyze the downstream development of the turbulence, the present work calculates the normalized mean absolute error (nMAE) between the STG-WMLBM and DNS, which is defined as nMAE :¼ R m i¼1 jy i;LBM À y i;DNS j R m i¼1 y i;DNS where m is the total number of the data points from the results, and y i;LBM and y i;DNS are the y axis data for LBM results and DNS reference, respectively.Figure 7 illustrates nMAE of hu 0 u 0 i þ along the downstream direction from the inlet to the outlet.The trend observed    to y þ % 260; 130; 86:7, respectively.Although the first few near-wall cells are slightly off the reference, the STG results match well with DNS further away from the wall; see Fig. 9.In Fig. 9, it can be observed that as the resolution increases, the u þ values are overestimated when compared to the DNS data at high y þ values.This phenomenon is a well-known effect in LES modeling named as the log-layer mismatch.The discrepancy may be attributed to the SGS turbulence modeling and the choice of shear estimation from either the first or second near-wall layer.Interested readers can refer to Asmuth et al. 44 for a discussion on the choice of shear estimation.Moreover, the results also reveal that at higher resolutions, the first off-wall layer exhibits a greater discrepancy with DNS than at lower resolutions.One possible explanation for this is the use of secondlayer shear stress to model the first layer.However, using first-layer information to estimate shear stress may increase the log-layer mismatch. 44,58Further investigation is required to fully understand this effect.The analysis of hu 0 v 0 i þ as a function of y=d is carried out at all three resolutions.Figure 10 displays the results of the STG data with the DNS reference.The results converge to the DNS reference around x=d ¼ 6 to x=d ¼ 8.

V. CONCLUSIONS
This paper presents a synthetic turbulent generator (STG) model based on the LES-LBM framework for high friction Reynolds number simulations (Re s ¼ 1000; 2000; 5200).We use wall function based on Reichardt's law 46 in combination with the force-based method, which simplifies the implementation.The STG simulations are examined at different resolutions [H ¼ 20, 40, 60 (LBU)] and compared with the DNS data, showing immediate convergence to the mean velocity field from the inlet of the channel flow.Further analysis of the Reynolds stress indicates that convergence to the DNS data occurs around x=d ¼ 6 to x=d ¼ 8.The presented STG model is computationally efficient and quickly converges to the DNS data even at low resolutions, which is promising for high Reynolds-number applications.One interesting future perspective for the presented framework is how to couple the LES-LBM turbulent flow back to RANS.Another future work can be to compare different turbulent generation methods for LBM.
The relaxation parameters that are not determined by the viscosity are set to x

i
are the probability density functions that point toward the boundary and are parallel to the boundary.
ARTICLE pubs.aip.org/aip/pofPhys.Fluids 35, 065115 (2023); doi: 10.1063/5.015352635, 065115-6 in the figure clearly demonstrates the level of agreement between the LBM results and the DNS data at different x=d.The error in the LBM results is observed to be in good agreement with the DNS data around x=d ¼ 8.This indicates that the LBM approach has effectively captured the turbulence characteristics in the flow near this location.The error then stabilizes until x=d ¼ 15, confirming the robustness of the LBM approach in the intermediate region of the flow.However, near the outlet, the error is observed to increase due to the presence of the buffer layer.This is a common phenomenon in turbulence simulations, as the buffer layer is known to have a significant impact on the accuracy of the simulation results.Nevertheless, the figure highlights the overall reliability of the STG-WMLBM approach in predicting the turbulence characteristics in the flow.Next, we will present the friction velocity u s =u theory at the three different friction Reynolds numbers Re s ¼ 1000; 2000; 5200.