The impact of intrinsic compressibility effects – changes in fluid volume due to pressure variations – on high-speed wall-bounded turbulence has often been overlooked or incorrectly attributed to mean property variations. To quantify these intrinsic compressibility effects unambiguously, we perform direct numerical simulations of compressible turbulent channel flows with nearly uniform mean properties. Our simulations reveal that intrinsic compressibility effects yield a significant upward shift in the logarithmic mean velocity profile that can be attributed to the reduction in the turbulent shear stress. This reduction stems from the weakening of the near-wall quasi-streamwise vortices. In turn, we attribute this weakening to the spontaneous opposition of sweeps and ejections from the near-wall expansions and contractions of the fluid, and provide a theoretical explanation for this mechanism. Our results also demonstrate that intrinsic compressibility effects play a crucial role in the increase in inner-scaled streamwise turbulence intensity in compressible flows, as compared with incompressible flows, which was previously regarded to be an effect of mean property variations alone.

We introduce CaLES, a GPU-accelerated finite-difference solver designed for large-eddy simulations (LES) of incompressible wall-bounded flows in massively parallel environments. Built upon the existing direct numerical simulation (DNS) solver CaNS, CaLES relies on low-storage, third-order Runge-Kutta schemes for temporal discretization, with the option to treat viscous terms via an implicit Crank-Nicolson scheme in one or three directions. A fast direct solver, based on eigenfunction expansions, is used to solve the discretized Poisson/Helmholtz equations. For turbulence modeling, the classical Smagorinsky model with van Driest near-wall damping and the dynamic Smagorinsky model are implemented, along with a logarithmic law wall model. GPU acceleration is achieved through OpenACC directives, following CaNS-2.3.0. Performance assessments were conducted on the Leonardo cluster at CINECA, Italy. Each node is equipped with one Intel Xeon Platinum 8358 CPU (2.60 GHz, 32 cores) and four NVIDIA A100 GPUs (64 GB HBM2e), interconnected via NVLink 3.0 (200 GB/s). The inter-node communication bandwidth is 25 GB/s, supported by a DragonFly+ network architecture with NVIDIA Mellanox InfiniBand HDR. Results indicate that the computational speed on a single GPU is equivalent to approximately 15 CPU nodes, depending on the treatment of viscous terms and the subgrid-scale model, and that the solver efficiently scales across multiple GPUs. The predictive capability of CaLES has been tested using multiple flow cases, including decaying isotropic turbulence, turbulent channel flow, and turbulent duct flow. The high computational efficiency of the solver enables grid convergence studies on extremely fine grids, pinpointing non-monotonic grid convergence for wall-modeled LES.

We present a massively parallel GPU-accelerated solver for direct numerical simulations of transitional and turbulent flat-plate boundary layers and channel flows involving fluids in non-ideal thermodynamic states. While several high-fidelity solvers are currently available as open source, all of them are restricted to the ideal-gas region. In contrast, the CUBic Equation of state Navier-Stokes solver (CUBENS) can accurately model and simulate the non-ideal thermodynamics of single-phase compressible fluids in the vicinity of the vapor-liquid saturation line or the thermodynamic critical point. By employing high-order finite-difference schemes and convective terms in split, kinetic-energy-, and entropy-preserving form, the solver is numerically stable, and robust with minimal numerical dissipation, enabling it to capture the steep variations of non-ideal thermodynamic properties. For cost-effective high-fidelity simulations, in addition to MPI parallelization, CUBENS is GPU-accelerated using OpenACC directives for computation offloading, and asynchronous GPU-aware MPI for efficient GPU-GPU communication. Moreover, CUBENS is compatible with both NVIDIA and AMD GPU architectures, achieving significant performance results while ensuring energy-efficient simulations. For instance, using 64 NVIDIA A100 GPUs compared to 8192 CPUs at the same computational cost results in a speedup of approximately 130×. In multi-node and multi-GPU configurations ranging from 2 to 128 compute nodes (8 to 512 GPUs), a strong scaling efficiency of around 52% and a weak scaling efficiency of 0.88 with 1024³ points per GPU, corresponding to approximately 5 billion degrees of freedom, are achieved. The CUBENS solver is validated against selected cases from the literature, covering transitional to turbulent ideal and non-ideal flows up to the transonic regime. In particular, we demonstrate the solver's suitability and applicability for direct numerical simulations of transitional boundary layers with fluids at supercritical pressure and with buoyancy effects. The development of this high-fidelity solver offers the potential for future fundamental research in non-ideal compressible fluid dynamics. Program summary: Program Title: CUBic Equation of state Navier-Stokes (CUBENS) CPC Library link to program files: https://doi.org/10.17632/6jfy758gyv.1 Developer's repository link: https://github.com/pcboldini/CUBENS Licensing provisions: MIT Programming language: Fortran 90, OpenACC, MPI, Python, MATLAB Nature of problem: This code solves the three-dimensional Navier-Stokes equations for non-ideal gas flows in a Cartesian domain, applicable to boundary layers and channels. Solution method: This code uses high-order central finite-differences with split-convective form, preserving kinetic energy and entropy (KEEP) and pressure-equilibrium-preserving (PEP) property, for spatial discretization. The time advancement is performed with a third-order Total Variation Diminishing low-storage Runge-Kutta scheme. Flow non-ideality is accounted for by cubic equations of state and complex transport-properties models. Alongside MPI parallelization, the solver is GPU-accelerated using OpenACC for computation offloading and CPU-GPU data transfer, along with GPU-aware MPI for GPU-GPU communication.

In this work we propose a novel method to ensure important entropy inequalities are satisfied semi-discretely when constructing reduced order models (ROMs) on nonlinear reduced manifolds. We are in particular interested in ROMs of systems of nonlinear hyperbolic conservation laws. The so-called entropy stability property endows the semi-discrete ROMs with physically admissible behaviour. The method generalizes earlier results on entropy-stable ROMs constructed on linear spaces. The ROM works by evaluating the projected system on a well-chosen approximation of the state that ensures entropy stability. To ensure accuracy of the ROM after this approximation we locally enrich the tangent space of the reduced manifold with important quantities. Using numerical experiments on some well-known equations (the inviscid Burgers equation, shallow water equations and compressible Euler equations) we show the improved structure-preserving properties of our ROM compared to standard approaches and that our approximations have minimal impact on the accuracy of the ROM. We additionally generalize the recently proposed polynomial reduced manifolds to rational polynomial manifolds and show that this leads to an increase in accuracy for our experiments.

Suspensions of finite-size solid particles in a turbulent pipe flow are found in many industrial and technical flows. Due to the ample parameter space consisting of particle size, concentration, density and Reynolds number, a complete picture of the particle–fluid interaction is still lacking. Pressure drop predictions are often made using viscosity models only considering the bulk solid volume fraction. For the case of turbulent pipe flow laden with neutrally buoyant spherical particles, we investigate the pressure drop and overall drag (friction factor), fluid velocity and particle distribution in the pipe. We use a combination of experimental (MRV) and numerical (DNS) techniques and a continuum flow model. We find that the particle size and the bulk flow rate influence the mean fluid velocity, velocity fluctuations and the particle distribution in the pipe for low flow rates. However, the effects of the added solid particles diminish as the flow rate increases. We created a master curve for drag change compared to single-phase flow for the particle-laden cases. This curve can be used to achieve more accurate friction factor predictions than the traditional modified viscosity approach that does not account for particle size.

In particle-laden turbulent wall flows, lift forces can influence the near-wall turbulence. This has been observed recently in particle-resolved simulations, which, however, are too expensive to be used in upscaled models. Instead, point-particle simulations have been the method of choice to simulate the dynamics of these flows during the last decades. While this approach is simpler, cheaper and physically sound for small inertial particles in turbulence, some issues remain. In the present work, we address challenges associated with lift force modelling in turbulent wall flows and the impact of lift forces in the near-wall flow. We performed direct numerical simulations of small inertial point particles in turbulent channel flow for fixed Stokes number and mass loading while varying the particle size. Our results show that the particle dynamics in the buffer region, causing the apparent particle-to-fluid slip velocity to vanish, raises major challenges for modelling lift forces accurately. While our results confirm that lift forces have little influence on particle dynamics for sufficiently small particle sizes, for inner-scaled diameters of order one and beyond, lift forces become quite important near the wall. The different particle dynamics under lift forces results in the modulation of streamwise momentum transport in the near-wall region. We analyse this lift-induced turbulence modulation for different lift force models, and the results indicate that realistic models are critical for particle-modelled simulations to correctly predict turbulence modulation by particles in the near-wall region.

Direct numerical simulations of channel flow and temporal boundary layer at a Reynolds number Reτ=1500 are used to assess the scale-by-scale mechanisms of wall turbulence. From the peak of turbulence production embedded at the small scales of the near-wall region, spatially ascending reverse cascades are generated that move through self-similar eddies growing in size with the wall distance. These fluxes are followed by spatially ascending forward cascades through detached eddies thus reaching sufficiently small scales where eventually scale energy is dissipated. This phenomenology is shared by both boundary layer and channel flow and is recognized as a robust physical feature characterizing wall turbulence in general. Specific features related to the flow configuration are indeed identified in the outer region. In particular, the central region of channels is characterized by a generalized Richardson energy cascade where large scales are in equilibrium with small scales at different wall distances through a combined forward cascade and spatial flux. On the contrary, the interface region of boundary layers is characterized by an almost two-dimensional physics where spatially ascending reverse cascades sustain long and wide interface structures with a forward cascade that survives only in the wall-normal scales. The overall scenario consists in a variety of scale motions that while protruding from the turbulent core towards the external region, squeeze at the interface thus sustaining vertical shear in a thin layer. The observed multidimensional physics sheds light on the complex interactions between outer entrainment and near-wall self-sustaining mechanisms with possible repercussions for theories.

We propose a numerical method tailored to perform interface-resolved simulations of evaporating multicomponent two-phase flows. The novelty of the method lies in the use of Robin boundary conditions to couple the transport equations for the vaporized species in the gas phase and the transport equations of the same species in the liquid phase. The Robin boundary condition is implemented with the cost-effective procedure proposed by Chai et al. [1] and consists of two steps: (1) calculating the normal derivative of the mass fraction fields in cells adjacent to the interface through the reconstruction of a linear polynomial system, and (2) extrapolating the normal derivative and the ghost value in the normal direction using a linear partial differential equation. This methodology yields a second-order accurate solution for the Poisson equation with a Robin boundary condition and a first-order accurate solution for the Stefan problem. The overall methodology is implemented in an efficient two-fluid solver, which includes a Volume-of-Fluid (VoF) approach for the interface representation, a divergence-free extension of the liquid velocity field onto the entire domain to transport the VoF, and the temperature equation to include thermal effects. We demonstrate the convergence of the numerical method to the analytical solution for multicomponent isothermal evaporation and observe good overall computational performance for simulating non-isothermal evaporating two-fluid flows in two and three dimensions.

This study presents direct numerical simulation results of two-layer Rayleigh–Bénard convection, investigating the previously unexplored Rayleigh–Weber parameter space 106≤Ra≤108 and 102≤We≤103. Global properties, such as the Nusselt and Reynolds numbers, are compared against the extended Grossmann–Lohse theory for two fluid layers, confirming a weak Weber number dependence for all global quantities and considerably larger Reynolds numbers in the lighter fluid. Statistics of the flow reveal that the interface fluctuates more intensely for larger Weber and smaller Rayleigh numbers, something also reflected in the increased temperature root mean square values next to the interface. The dynamics of the deformed two-fluid interface is further investigated using spectral analysis. Temporal and spatial spectrum distributions reveal a capillary wave range at small Weber and large Rayleigh numbers, and a secondary energy peak at smaller Rayleigh numbers. Furthermore, the maxima of the space–time spectra lie in an intermediate dispersion regime, between the theoretical predictions for capillary and gravity-capillary waves, showing that the gravitational energy of the interfacial waves is strongly altered by temperature gradients.

We consider suspensions of finite-size neutrally buoyant rigid spherical particles in channel flow and investigate the relevance of different momentum transfer mechanisms and the relation between the local particle dynamics and the bulk flow properties in the highly inertial regime. Interface-resolved simulations are performed in the range of Reynolds numbers and solid volume fractions. The Lagrangian particle statistics show that pair interactions are highly inhomogeneous and dependent on the distance from the wall: in their vicinity, the underlying mean shear drives the pair interactions, while a high degree of isotropy, dictated by more frequent collisions, characterizes the core region. Analysis of the momentum balance reveals that while the particle-induced stresses govern the dynamics in dense conditions, and moderate Reynolds numbers, <![CDATA[$Re, the turbulent stresses take over at higher Reynolds numbers. This behaviour is associated with a reduced particle migration toward the channel core, which decreases the importance of the particle-induced stress and increases the turbulent activity. Our results indicate that Reynolds stresses and the associated velocity fluctuations, characteristics of near-wall turbulence, prevail at high inertia over the resistance to deformation presented by the particles for volume fractions lower than 30 %.

The geometrically complex mechanisms of energy transfer in the compound space of scales and positions of wall turbulent flows are investigated in a temporally evolving boundary layer. The phenomena consist of spatially ascending reverse and forward cascades from the small production scales of the buffer layer to the small dissipative scales distributed among the entire boundary layer height. The observed qualitative behaviour conforms with previous results in turbulent channel flow, thus suggesting that the observed phenomenology is a robust statistical feature of wall turbulence in general. An interesting feature is the behaviour of energy transfer at the turbulent/non-turbulent interface, where forward energy cascade is found to be almost absent. In particular, the turbulent core is found to sustain a variety of large-scale wall-parallel motions at the turbulent interface through weak but persistent reverse energy cascades. This behaviour conforms with previous results in free shear flows, thus suggesting that the observed phenomenology is a robust statistical feature of turbulent shear flows featuring turbulent/non-turbulent interfaces in general.

We present the Fluid Transport Accelerated Solver, FluTAS, a scalable GPU code for multiphase flows with thermal effects. The code solves the incompressible Navier-Stokes equation for two-fluid systems, with a direct FFT-based Poisson solver for the pressure equation. The interface between the two fluids is represented with the Volume of Fluid (VoF) method, which is mass conserving and well suited for complex flows thanks to its capacity of handling topological changes. The energy equation is explicitly solved and coupled with the momentum equation through the Boussinesq approximation. The code is conceived in a modular fashion so that different numerical methods can be used independently, the existing routines can be modified, and new ones can be included in a straightforward and sustainable manner. FluTAS is written in modern Fortran and parallelized using hybrid MPI/OpenMP in the CPU-only version and accelerated with OpenACC directives in the GPU implementation. We present different benchmarks to validate the code, and two large-scale simulations of fundamental interest in turbulent multiphase flows: isothermal emulsions in HIT and two-layer Rayleigh-Bénard convection. FluTAS is distributed through a MIT license and arises from a collaborative effort of several scientists, aiming to become a flexible tool to study complex multiphase flows. Program summary: Program Title: : Fluid Transport Accelerated Solver, FluTAS. CPC Library link to program files: https://doi.org/10.17632/tp6k8wky8m.1 Developer's repository link: https://github.com/Multiphysics-Flow-Solvers/FluTAS.git. Licensing provisions: MIT License. Programming language: Fortran 90, parallelized using MPI and slab/pencil decomposition, GPU accelerated using OpenACC directives. External libraries/routines: FFTW, cuFFT. Nature of problem: FluTAS is a GPU-accelerated numerical code tailored to perform interface resolved simulations of incompressible multiphase flows, optionally with heat transfer. The code combines a standard pressure correction algorithm with an algebraic volume of fluid method, MTHINC [1]. Solution method: the code employs a second-order-finite difference discretization and solves the two-fluid Navier-Stokes equation using a projection method. It can be run both on CPU-architectures and GPU-architectures.

This paper presents a performance analysis of pencil domain decomposition methodologies for three-dimensional Computational Fluid Dynamics (CFD) codes for turbulence simulations, on several large GPU-accelerated clusters. The performance was assessed for the numerical solution of the Navier-Stokes equations in two codes which require the calculation of Fast-Fourier Transforms (FFT): a tri-periodic pseudo-spectral solver for isotropic turbulence, and a finite-difference solver for canonical turbulent flows, where the FFTs are used in its Poisson solver. Both codes use a newly developed transpose library that automatically determines the optimal domain decomposition and communication backend on each system. We compared the performance across systems with very different node topologies and available network bandwidth, to show how these characteristics impact decomposition selection for best performance. Additionally, we assessed the performance of several communication libraries available on these systems, such as Open-MPI, IBM Spectrum MPI, Cray MPI, the NVIDIA Collective Communication Library (NCCL), and NVSHMEM. Our results show that the optimal combination of communication backend and domain decomposition is highly system-dependent, and that the adaptive decomposition library is key in ensuring efficient resource usage with minimal user effort.

We present a multi-block finite-difference solver for massively parallel Direct Numerical Simulations (DNS) of incompressible flows. The algorithm combines the versatility of a multi-block solver with the method of eigenfunctions expansions, to speedup the solution of the pressure Poisson equation. This is achieved by employing FFT-based transforms along one homogeneous direction, which effectively reduce the problem complexity at a low cost. These FFT-based expansions are implemented in a framework that unifies all valid combinations of boundary conditions for this type of method. Subsequently, a geometric multigrid solver is employed to solve the reduced Poisson equation in a multi-block geometry. Particular care was taken here, to guarantee the parallel performance of the multigrid solver when solving the reduced linear systems equations. We have validated the overall numerical algorithm and assessed its performance in a massively parallel setting. The results show that 2- to 8-fold reductions in computational cost may be easily achieved when exploiting FFT-accelerated for the solution of the Poisson equation. The solver, SNaC, has been made freely available and open-source under the terms of an MIT license.

This work presents the GPU acceleration of the open-source code CaNS for very fast massively-parallel simulations of canonical fluid flows. The distinct feature of the many-CPU Navier–Stokes solver in CaNS is its fast direct solver for the second-order finite-difference Poisson equation, based on the method of eigenfunction expansions. The solver implements all the boundary conditions valid for this type of problems in a unified framework. Here, we extend the solver for GPU-accelerated clusters using CUDA Fortran. The porting makes extensive use of CUF kernels and has been greatly simplified by the unified memory feature of CUDA Fortran, which handles the data migration between host (CPU) and device (GPU) without defining new arrays in the source code. The overall implementation has been validated against benchmark data for turbulent channel flow and its performance assessed on a NVIDIA DGX-2 system (16 T V100 32Gb, connected with NVLink via NVSwitch). The wall-clock time per time step of the GPU-accelerated implementation is impressively small when compared to its CPU implementation on state-of-the-art many-CPU clusters, as long as the domain partitioning is sufficiently small that the data resides mostly on the GPUs. The implementation has been made freely available and open source under the terms of an MIT license.

Polymer additives are commonly utilized to manipulate bubbly flows in various applications. Here we investigate the effects of clean and contaminated bubbles driven upward (upflow) in Newtonian and viscoelastic turbulent channel flows. Interface-resolved direct numerical simulations are performed to examine sole and combined effects of soluble surfactant and viscoelasticity using an efficient three-dimensional finite-difference-front-tracking method. The incompressible flow equations are solved fully coupled with the FENE-P viscoelastic model and the equations governing interfacial and bulk surfactant concentrations. The latter coupling is accomplished by a nonlinear equation of state that relates the surface tension to the surfactant concentration. For Newtonian turbulent bubbly flows, the effects of Triton X-100 and 1-pentanol surfactant are examined. It is observed that the sorption kinetics highly affect the dynamics of bubbly flow. A minute amount of Triton X-100 is found to be sufficient to prevent the formation of bubble clusters restoring the single-phase behavior while even two orders of magnitude more 1-pentanol surfactant is not adequate to prevent the formation of layers. For viscoelastic turbulent flows, it is found that the viscoelasticity promotes formation of the bubble wall-layers and thus the polymer drag reduction is completely lost for the surfactant-free bubbly flows, while the addition of small amount of surfactant (Triton X-100) in this system restores the polymer drag reduction resulting in 25% drag reduction for the Wi=4 case.

We use interface-resolved simulations to study near-wall turbulence modulation by small inertial particles, much denser than the fluid, in dilute/semi-dilute conditions. We considered three bulk solid mass fractions, and, with only the latter two showing turbulence modulation. The increase of the drag is strong at, but mild in the densest case. Two distinct regimes of turbulence modulation emerge: for smaller mass fractions, the turbulence statistics are weakly affected and the near-wall particle accumulation increases the drag so the flow appears as a single-phase flow at slightly higher Reynolds number. Conversely, at higher mass fractions, the particles modulate the turbulent dynamics over the entire flow, and the interphase coupling becomes more complex. In this case, fluid Reynolds stresses are attenuated, but the inertial particle dynamics near the wall increases the drag via correlated velocity fluctuations, leading to an overall drag increase. Hence, we conclude that, although particles at high mass fractions reduce the fluid turbulent drag, the solid phase inertial dynamics still increases the overall drag. However, inspection of the streamwise momentum budget in the two-way coupling limit of vanishing volume fraction, but finite mass fraction, indicates that this trend could reverse at even higher particle load.

We present a numerical method for interface-resolved simulations of evaporating two-fluid flows based on the volume-of-fluid (VoF) method. The method has been implemented in an efficient FFT-based two-fluid Navier-Stokes solver, using an algebraic VoF method for the interface representation, and extended with the transport equations of thermal energy and vaporized liquid mass for the single-component evaporating liquid in an inert gas. The conservation of vaporizing liquid and computation of the interfacial mass flux are performed with the aid of a reconstructed signed-distance field, which enables the use of well-established methods for phase change solvers based on level-set methods. The interface velocity is computed with a novel approach that ensures accurate mass conservation, by constructing a divergence-free extension of the liquid velocity field onto the entire domain. The resulting approach does not depend on the type of interface reconstruction (i.e. can be employed in both algebraic and geometrical VoF methods). We extensively verified and validated the overall method against several benchmark cases, and demonstrated its excellent mass conservation and good overall performance for simulating evaporating two-fluid flows in two and three dimensions.

The rupture of a liquid film due to co-flowing turbulent shear flows in the gas phase is studied using a volume-of-fluid method. To simulate this multiphase problem, we use a simplified numerical setup where the liquid film is 'sandwiched' between two fully developed boundary layers from a turbulent channel simulation. The film deforms and eventually ruptures within the shear zone created by the co-flows. This efficient setup allows systematic variation of physical parameters to gauge their role in the aerodynamically-driven deformation and rupture of a liquid film under fully developed sheared turbulence. The developing pressure field over the deforming film and related aerodynamic effects is studied, the importance of which was previously suggested by other authors, and in particular the role of the inviscid lift and drag forces. A cumulative lift force is introduced to capture the effect of the alternating pressure minima and maxima forming over the film which amplify and eventually rupture the film. A velocity scale derived from the lift-induced drag force reflects the state of the turbulent boundary layer over the film and collapses the temporal development of this cumulative lift force as well as the amplitude of film deformation with some success for the different film thicknesses and Reynolds numbers.

Equation (2.10) in Costa, Brandt & Picano (2020) for the lift force model used in the point-particle direct numerical simulations (DNS), and which is derived from the classical lift force of Saffman (1965), (Equation presented) does not correspond to the force model actually used in the point-particle DNS with lift force presented in the manuscript. Instead, the following equation was used: (Equation presented) which replaces the first occurrence of the term |U_s| on the right-hand-side of (1) with |ω|D. We recall that two cases were considered in the manuscript depending on the value of J in the lift force equation: J = 1 in the case denoted PP-Saffman; and J given by (Equation presented) with ϵ = √|ω|ν/|U_s|, in the case denoted PP-McLaughlin. Also, equation (2.13) of the manuscript - describing the perfectly elastic hard-sphere rebound - is incorrect; the term D/2 should be D: (Equation presented) Despite the lapse in the manuscript, equation (4) was implemented correctly (Costa et al. 2020). The results from the point-particle DNS with the model reported in (2.10) of Costa et al. (2020) ((1) above) differ from those reported in the manuscript, and are shown (Figure presented) in figure 1 (cf. figures 7 and 8 of Costa et al. (2020)). The statistics presented here have been collected in the fully developed state from 600 samples over a time interval of 250h/U_b, which ensured statistical convergence of the results. The results from the point-particle cases presented in the original manuscript are also reproduced here with this (higher) statistical sampling, and show very minor differences with respect to figures 7 and 8 of Costa et al. (2020). In light of these results, the conclusions drawn from the results in the last section of § 3 of the manuscript must be therefore reformulated: (i) The Saffman lift model does not correctly predict the near-wall statistics of the interface-resolved DNS very close to the wall, including the near-wall concentration peak. (ii) The equation proposed by Mei (1992) that fits the model of McLaughlin (1991) shows results similar to those reported in the original manuscript for this model. That is, it predicts well the near-wall concentration peak, and fails to predict the other observables near the wall. (iii) Equation (2) for F_l presented above, with J = 1, predicts very well all the observables in figure 1. We have therefore accidentally discovered that the expression (2) for F_l predicts the observed particle statistics very well. Still, the reason for the strikingly good agreement remains elusive to us. We hope that this result can be further exploited for the improvement lift force models for point-particle simulations of wall-bounded turbulent flows.