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.@en

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 |Us| 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 ϵ = √|ω|ν/|Us|, 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/Ub, 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 Fl presented above, with J = 1, predicts very well all the observables in figure 1. We have therefore accidentally discovered that the expression (2) for Fl 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.@en

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 %. @en

The gravity-driven motion of rigid particles in a viscous fluid is relevant in many natural and industrial processes, yet this has mainly been investigated for spherical particles. We therefore consider the sedimentation of non-spherical (spheroidal) isolated and particle pairs in a viscous fluid via numerical simulations using the Immersed Boundary Method. The simulations performed here show that the critical Galileo number for the onset of secondary motions decreases as the spheroid aspect ratio departs from 1. Above this critical threshold, oblate particles perform a zigzagging motion whereas prolate particles rotate around the vertical axis while having their broad side facing the falling direction. Instabilities of the vortices in the wake follow when farther increasing the Galileo number. We also study the drafting-kissing-tumbling associated with the settling of particle pairs. We find that the interaction time increases significantly for non-spherical particles and, more interestingly, spheroidal particles are attracted from larger lateral displacements. This has important implications for the estimation of collision kernels and can result in increasing clustering in suspensions of sedimenting spheroids.@en

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. @en

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.@en

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.@en

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.@en

We present modal and non-modal linear stability analyses of Poiseuille flow through a plane channel with a porous substrate modeled using the Volume Averaged Navier–Stokes (VANS) equations. Modal stability analysis shows the destabilization of the flow with increasing porosity of the layer. The instability mode originates from the homogeneous fluid region of the channel for all the values of porosity considered but the governing mechanism changes. Perturbation kinetic energy analysis reveals the importance of viscous dissipation at low porosity values while dissipation in the porous substrate becomes significant at higher porosity. Scaling analysis highlights the invariance of the critical wavenumber with changing porosity. On the other hand, the critical Reynolds number remains invariant at low porosity and scales as Rec∼(H∕δ)1.4 at high porosity where δ is the typical thickness of the vorticity layer at the fluid–porous interface. This reveals the existence of a Tollmien–Schlichting-like viscous instability mechanism at low porosity values, and Rayleigh analysis shows the presence of an inviscid instability mechanism at high porosity. For the whole range of porosities considered, the non-modal analysis shows that the optimal mechanism responsible for transient energy amplification is the lift-up effect, giving rise to streaky structure as in single-phase plane Poiseuille flow. The present results strongly suggest that the transition to turbulence follows the same path as that of classical Poiseuille flow at low porosity values, while it is dictated by the modal instability for high porosity values.@en

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.@en

We present a direct comparison between interface-resolved and one-way-coupled point-particle direct numerical simulations (DNS) of gravity-free turbulent channel flow laden with small inertial particles, with high particle-to-fluid density ratio and diameter of approximately three viscous units. The most dilute flow considered, solid volume fraction , shows the particle feedback on the flow to be negligible, whereas differences with respect to the unladen case, notably a drag increase of approximately 10 %, are found for a volume fraction. This is attributed to a dense layer of particles at the wall, caused by turbophoresis, flowing with large particle-to-fluid apparent slip velocity. The most dilute case is therefore taken as the benchmark for assessing the validity of a widely used point-particle model, where the particle dynamics results only from inertial and nonlinear drag forces. In the bulk of the channel, the first- and second-order moments of the particle velocity from the point-particle DNS agree well with those from the interface-resolved DNS. Close to the wall, however, most of the statistics show major qualitative differences. We show that this difference originates from the strong shear-induced lift force acting on the particles in the near-wall region. This mechanism is well captured by the lift force model due to Saffman (J. Fluid Mech., vol. 22 (2), 1965, pp. 385-400), while other widely used, more elaborate, approaches aiming at extending the lift model for a wider range of particle Reynolds numbers can actually underpredict the magnitude of the near-wall particle velocity fluctuations for the cases analysed here.@en

We use interface-resolved direct numerical simulations to study the dynamics of a single sediment particle in a turbulent open channel flow over a fixed porous bed. The relative strength of the gravitational acceleration, quantified by the Galileo number, is varied so as to reproduce the different modes of sediment transport - resuspension, saltation and rolling. The results show that the sediment dynamics at lower Galileo numbers (i.e.resuspension and saltation) are mainly governed by the mean flow. Here, the regime of motion can be predicted by the ratio between the gravity and the shear-induced boundary force. In these cases, the sediment particle rapidly takes off when exposed to the flow, and proceeds with an oscillatory motion. Increasing the Galileo number, the frequency of these oscillations increases and their amplitude decreases, until the transport mode switches from resuspension to saltation. In this case, the sediment travels by short successive collisions with the bed. Further increasing the Galileo number, the particle rolls without detaching from the bed. Differently from the previous modes, the motion is triggered by extreme turbulent events, and the particle response depends on the specific initial conditions, at fixed Reynolds number. The results reveal that close to the onset of sediment motion, only turbulent sweeps can effectively trigger the particle motion by increasing the stagnation pressure upstream. We show that for the parameters in this study, a criterion based on the streamwise flow-induced force can successfully predict the incipient movement.@en

We report experimental observations of turbulent flow with spherical particles in a square duct. Three particle sizes, namely 2H/dp=40, 16 and 9 (2H being the duct full height and dp being the particle diameter), are investigated. The particles are nearly neutrally buoyant with a density ratio of 1.0035 and 1.01 with respect to the suspending fluid. Refractive index matched-particle image velocimetry (RIM-PIV) is used for fluid velocity measurement even at the highest particle volume fraction (20%) and particle tracking velocimetry (PTV) for the particle velocity statistics for the flows seeded with particles of the two largest sizes, whereas only pressure measurements are reported for the smallest particles. Settling effects are seen at the lowest bulk Reynolds number Re2H ≈ 000, whereas, at the highest Re2H ≈, particles are in almost full suspension. The friction factor of the suspensions is found to be significantly larger than that of single-phase duct flow at the lower Re2H investigated; however, the difference decreases when increasing the flow rate and the total drag approaches the values of the single-phase flow at the higher Reynolds number considered, Re2H=27 000. The pressure drop is found to decrease with the particle diameter for volume fractions lower than φ=10% for nearly all Re2H investigated. However, at the highest volume fraction φ=10%, we report a peculiar non-monotonic behaviour: the pressure drop first decreases and then increases with increasing particle size. The decrease of the turbulent drag with particle size at the lowest volume fractions is related to an attenuation of the turbulence. The drag increase for the two largest particle sizes at φ=10%, however, occurs despite this large reduction of the turbulent stresses, and it is therefore due to significant particle-induced stresses. At the lowest Reynolds number, the particles reside mostly in the bottom half of the duct, where the mean velocity significantly decreases; the flow is similar to that in a moving porous bed near the bottom wall and to turbulent duct flow with low particle concentration near the top wall.@en

We use interface-resolved numerical simulations to study finite-size effects in turbulent channel flow of neutrally buoyant spheres. Two cases with particle sizes differing by a factor of two, at the same solid volume fraction of 20 % and bulk Reynolds number are considered. These are complemented with two reference single-phase flows: the unladen case, and the flow of a Newtonian fluid with the effective suspension viscosity of the same mixture in the laminar regime. As recently highlighted in Costa et al. (Phys. Rev. Lett., vol. 117, 2016, 134501), a particle–wall layer is responsible for deviations of the mesoscale-averaged statistics from what is observed in the continuum limit where the suspension is modelled as a Newtonian fluid with (higher) effective viscosity. Here we investigate in detail the fluid and particle dynamics inside this layer and in the bulk. In the particle–wall layer, the near-wall inhomogeneity has an influence on the suspension microstructure over a distance proportional to the particle size. In this layer, particles have a significant (apparent) slip velocity that is reflected in the distribution of wall shear stresses. This is characterized by extreme events (both much higher and much lower than the mean). Based on these observations we provide a scaling for the particle-to-fluid apparent slip velocity as a function of the flow parameters. We also extend the scaling laws in Costa et al. (Phys. Rev. Lett., vol. 117, 2016, 134501) to second-order Eulerian statistics in the homogeneous suspension region away from the wall. The results show that finite-size effects in the bulk of the channel become important for larger particles, while negligible for lower-order statistics and smaller particles. Finally, we study the particle dynamics along the wall-normal direction. Our results suggest that single-point dispersion is dominated by particle–turbulence (and not particle–particle) interactions, while differences in two-point dispersion and collisional dynamics are consistent with a picture of shear-driven interactions.@en

Dense suspensions are usually investigated in the laminar limit where inertial effects are insignificant. In this regime, the main effect of the suspended phase is to alter the rheological behavior of the flow which always displays higher effective viscosity with respect to the carrier fluid. When the flow rate is high enough, i.e. at high Reynolds number, the flow may become turbulent and the interaction between solid and liquid phase modifies the turbulent dynamics that we know in single-phase fluids. In the present work, we study turbulent channel flows laden with finite-size particles at high volume fraction (F = 0:2) by means of Direct Numerical Simulations. A direct-forcing Immersed Boundary Method has been adopted to couple liquid and solid phases. The two-phase simulations have been performed fixing the bulk Reynolds number at Reb = Ub 2h=n = 12000 (Ub bulk velocity, h channel half-width and n the fluid kinematic viscosity). The particle size is relatively large with respect to the viscous length, i.e. 10 and 20 times, but smaller than large scales. We will present a detailed comparison of the statistical behavior of the particle-laden flow and the corresponding single-phase flow. The presence of the solid phase strongly alters the wall turbulence dynamics and its effect cannot be accounted only considering the higher rheological effective viscosity.@en

The macroscopic behavior of dense suspensions of neutrally buoyant spheres in turbulent plane channel flow is examined. We show that particles larger than the smallest turbulence scales cause the suspension to deviate from the continuum limit in which its dynamics is well described by an effective suspension viscosity. This deviation is caused by the formation of a particle layer close to the wall with significant slip velocity. By assuming two distinct transport mechanisms in the near-wall layer and the turbulence in the bulk, we define an effective wall location such that the flow in the bulk can still be accurately described by an effective suspension viscosity. We thus propose scaling laws for the mean velocity profile of the suspension flow, together with a master equation able to predict the increase in drag as a function of the particle size and volume fraction.@en

Suspensions of solid particles are frequently found in applications and environmental flows. Several studies concern the rheological properties of suspensions in laminar flows, but much less is known of turbulent suspensions. The present work fills this gap providing DNS data on dense suspensions of neutrally-buoyant rigid sphere in a turbulent channel flow at the bulk Reynolds number of Re = U0h/ν = 2800. We show that considering volume fractions Φ ≤ 0.1 the turbulent flow is similar to the unladen case with higher turbulence intensities. On the contrary, the flow behavior strongly changes at Φ = 0.2where the turbulence appears to be attenuated.@en

We present interface-resolved numerical simulations of turbulent channel flow laden with non-spherical rigid and neutrally-buoyant particles. We first focus on the case of oblate particles of aspect ratio 1/3 at volume fractions up to 15% and show that the turbulent drag is decreasing when increasing the particle volume fraction although the effective viscosity of the suspension actually increases. We relate the observed drag reduction to turbulence attenuation and to particle migration away from the near-wall region. Particles tend to align parallel to the wall with rotation rates significantly lower than those reported for spheres. In the second part of the study, we examine the effect of the particle slenderness on the observed drag reduction and show that the drag increases for flatter particles.@en

Inertial regimes in a channel flow of suspension of finite-size neutrally buoyant particles are studied for a wide range of Reynolds numbers: 500 Re 5000, and particle volume fractions: 0 0:3. The flow is classified in three different regimes according to the phase-averaged stress budget across the channel [2]. The laminar viscous regime at low Re and - where the viscous stress is the dominating term in the budget, the turbulent regime at high Re and relatively low where the momentum is mainly transferred by the action of the Reynolds stress and the inertial shear-thickening regime where the particle stress contributes the most to the significant enhancement of the wall shear stress. Particle distribution and dispersion properties provide additional evidence for the existence of the three different regimes.@en

Dense suspensions are usually investigated in the laminar limit where inertial effects are insignificant. When the flow rate is high enough, i.e. at high Reynolds number, the flow may become turbulent and the interaction between solid and liquid phases modifies the turbulence we know in single-phase fluids. In the present work, we study turbulent channel flows laden with finite-size particles at high volume fraction by means of Direct Numerical Simulations. A direct-forcing Immersed Boundary Method has been adopted to couple liquid and solid phases. We will show that the turbulence is attenuated in dense cases, even though the overall drag is increased because of the particle contribution to the total stress.@en