Numerical models for fluid-grains interactions : opportunities and limitations

In the framework of a multi-scale approach, we develop numerical models for suspension flows. At the micro scale level, we perform particle-resolved numerical simulations using a Distributed Lagrange Multiplier/Fictitious Domain approach. At the meso scale level, we use a two-way Euler/Lagrange approach with a Gaussian filtering kernel to model fluid-solid momentum transfer. At both the micro and meso scale levels, particles are individually tracked in a Lagrangian way and all inter-particle collisions are computed by a Discrete Element/Soft-sphere method. The previous numerical models have been extended to handle particles of arbitrary shape (non-spherical, angular and even non-convex) as well as to treat heat and mass transfer. All simulation tools are fully-MPI parallel with standard domain decomposition and run on supercomputers with a satisfactory scalability on up to a few thousands of cores. The main asset of multi scale analysis is the ability to extend our comprehension of the dynamics of suspension flows based on the knowledge acquired from the high-fidelity micro scale simulations and to use that knowledge to improve the meso scale model. We illustrate how we can benefit from this strategy for a fluidized bed, where we introduce a stochastic drag force model derived from micro-scale simulations to recover the proper level of particle fluctuations. Conversely, we discuss the limitations of such modelling tools such as their limited ability to capture lubrication forces and boundary layers in highly inertial flows. We suggest ways to overcome these limitations in order to enhance further the capabilities of the numerical models.


Introduction
The understanding of the dynamics of suspension flows is still incomplete, although major advances have been achieved over the past 50 years.Complementary to experiments and theory, numerical simulation has become an indispensable tool to get some more insight into the intricate momentum transfer between the continuous fluid and dispersed solid particles.Suspension flows exhibit a large range of scales: from the scale of an individual particle to the scale of the process.It is hence natural to adopt a multiscale description of these flows.A nowadays worthwhile research path to follow is to adopt a micro/meso/macro separation of scales with corresponding numerical models at these different scales combined to a bottom-up strategy.The bottom-up strategy, also called upwards cascade of knowledge, involves transferring the knowledge acquired at the smallest scale with high-fidelity simulations to the larger scales as a tool to improve numerical models at these large scales.
At the micro scale, Particle Resolved Simulation (PRS) has been extensively used to study fundamental aspects of suspension flow dynamics.PRS models contain very few assumptions and are hence expected to supply high-fidelity simulation results on momentum, heat and mass transfer in suspension flows.Body-fitted Figure 1: A discrete view on a grid of a micro/meso/macro multi-scale approach techniques as, e.g., Arbitrary Lagrangian-Eulerian, are nowadays not favored as the remeshing step is computationally too expensive in 3D.Among fixed mesh techniques, the most popular PRS methods include Immersed Boundary method (IB), lattice-Boltzmann method (LB) and Distributed Lagrange Multiplier/Fictitious Domain (DLM/FD) method (see [2] for an extended discussion).Initially developed to model momentum transfer, they can easily be extended to model heat and mass transfer too.
At the meso scale, mass and momentum conservation equations are averaged in control volumes encompassing a few particles and momentum transfer between the two phases is modeled using appropriate drag laws.Meso-scale methods provide a good basis for intermediate size simulations of suspension flows with a relatively wide range of Reynolds number and porosity.They represent an attractive alternative to PRS but the accurate treatment of fluid/particle interactions is the key to reliable numerical predictions.
Finally, macro-scale models correspond to different variants of the Euler/Euler model as, e.g., the Two Fluid Model (TFM).In this family of models, the solid phase is somehow assumed to behave as a fluid and both phases are solved in an average way based on closure laws for momentum transfer.Although quite advanced, the main drawback of this type of approach is its deficiency to model particle/particle interactions very well due to the use of (non necessarily extremely reliable) correlations for solid phase pressure and viscosity.
In this paper we introduced our numerical models at the micro scale and the meso scale and discuss how to transfer knowledge from micro to meso in the case of a bubbling fluidized bed as an illustration of what can be achieved with such a multi-scale approach.Finally, we discuss the current limitations of the different models and the next stages in the development of more advanced models for suspension flows.

Numerical models
In the following section, the mathematical aspects of both micro-and meso-scale models are described.We assume that the fluid is Newtonian and that the flow is incompressible.In both models the solid phase is treated in a direct fashion by solving Newton-Euler equations for each individual particle.All parameters in this section are dimensionless, while dimensional variables are distinguished by a "star" symbol.To non-dimensionalize the set of equations, the following scales are introduced: for distributed Lagrange multiplier (for the micro-scale model only) and ρ * f V * 2 c L * 2 c for contact force.

A DLM/FD formulation as a micro scale model
Our micro-scale approach is based on the DLM/FD formulation of [1] and is implemented in a Finite Volume fashion (see [2] for more details).The dimensionless and non-variational form of combined momentum continuity and momentum equations are as follows: where u, U, ω stand for the fluid velocity vector, the particle translational velocity vector and the particle angular velocity vector, respectively.The solid domain and whole fluid/particle domain are denoted by P and Ω, respectively.λ represents the distributed Lagrange multiplier vector, F c the contact force, V p = M * p /(ρ * s L * d c ) the particle volume, M * p the particle mass, d the dimension of the system, the particle inertia tensor, r the position vector and R the vector connecting the gravity center and the contact point.The following dimensionless numbers are also introduced:

A two-way Euler/Lagrange (TWEL) meso scale model
Based on the principle of locally averaged Navier-Stokes equations proposed by [3], fluid phase variables are solved on grid cells the volume of which is an order of magnitude larger than that of a particle.The dimensionless fluid conservation equations read: where D and ε denote the strain rate tensor and the fluid volume fraction, respectively.The pressure gradient term contains the hydrodynamic pressure only and F f p represents the fluid/particle interaction force.The most general expression of F f p is: where τ is the point-wise fluid stress tensor and n is the normal vector to the solid boundary.Since point-wise variables (pressure and viscous stress) on the solid boundary are not resolved in meso-scale models, the surface integration over ∂P cannot be carried out and a closure law is used instead to evaluate F f p .Among the different contributions to the total hydrodynamic force (drag, Basset, added mass, Magnus, Saffman), the drag force is generally the order 1 contribution.Many drag correlations are available for spheres in the literature.They all predict roughly the same values.For instance, the correlation of [4]  where C d is the drag coefficient for a single unhindered particle for Re p less than 1000 ( [5]) and the pre-factor ε −χ accounts for the presence of neighboring particles.

Lagrangian tracking of particles with collisions using a DEM/soft-sphere approach
The Newton-Euler equations for each individual particle are as follows: where denotes the buoyancy force, F c the inter-particle or particle-wall contact force and F f p the fluid/particle interaction force.Likewise, T c and T f p represent the contact-induced and the fluid-induced torque acting on the particle.The contact force and torque are computed by a spring-dashpot-frictional slider DEM/softsphere model (see [6]).For the micro-scale model, since the flow field is resolved at the particle level, no closure law is used for the fluid/particle interaction force F f p and torque T f p , and these quantities can be computed a posteriori using the Lagrange multipliers as: For the meso-scale model, however, the fluid/particle interaction force contains several terms: For high density ratios (ρ r > 8), we can safely neglect the added mass force (F am ) and the Basset force (F ba ).The lift force (F li f t ) acting on a particle is due to the rotation induced pressure distribution around it.This rotation can be the result of either the fluid local shear rate (i.e.Saffman force) or the particle own rotation (i.e.Magnus force).Unfortunately, correlations for these two forces available in the literature are limited to dilute regimes.Thus, to our knowledge, all meso-scale models in the literature (including ours) neglect the lift force.Finally, the fluid-induced torque T f p is also assumed to be zero in the meso-scale model, in the absence of both an appropriate correlation and reliable fundamental understanding of the importance of particle rotation in suspension flows.

Results
We consider the fluidized bed shown in Fig. 2 We simulate this flow configuration with the micro scale model and the meso scale model.Assuming that microscale simulation results are highly reliable, we directly compare the two sets of data produced by the two models and identify what is flawed in the meso-scale model.
Integral measures as average expansion of the bed (Fig. 2(b) and Fig. 2(c)) and pressure drop across the bed are well predicted by both models.Based on an extensive analysis of both Eulerian and Lagrangian particles statistics, we observe that particles motion is underpredicted by the meso-scale TWEL model in particular in the transverse direction.We identify two causes for this underprediction.First, our simple Bounding Cube (BC) projection kernel (Fig. 3(a)) is not accurate enough and using a Gaussian kernel (GK, Fig. 3(b)) enables us to recover part of the proper particles motion.The GK also decouples particle size and grid size, giving rise to the opportunity to use grid size smaller than particle diameter.Second, filtering PRS data reveals that for a pair of (Re, ε) the corresponding drag coefficient is not a deterministic (unique) value but a distribution.This led to the definition of a new stochastic drag law C d that models the sub-grid fluctuations.The total drag coefficient C t d is hence defined as C t d = Cd + C d where Cd is the deterministic value given by, e.g., [4].The parameters of the stochastic drag model C d are extracted from the filtered PRS data, such that no parameter is empirical.In Fig. 4, we plot the granular temperature in the vertical z and transverse x directions as a measure of particles fluctuations, parametrized by the level of sophistication of the TWEL model.It is very visible that we obtain a hierarchy of solutions: (i) BC, (ii) GK, and (iii) GK + stochastic drag, where the most sophisticated version of the model, i.e., GK + stochastic drag, is almost able to re-

Discussion and Perspectives
Advanced numerical models embedded in a multi-scale approach create many opportunities to gain novel physical insight into fluid-grains interactions that govern suspension flows.We have shown an example in which using high-fidelity PRS data is extremely helpful in identifying and partly curing some defects of meso-scale TWEL models.This analysis can be conducted for heat (and mass) transfer as our DLM/FD method (like IB or LB) and TWEL model extend easily to heat (and mass) transfer too.This offers the opportunity to examine a broder range of flow problems.Although very promising and efficient, this type of multi-scale approach also exhibits limitations.The first limitation is computational.In fact, PRS are extremely expensive.The PRS of the fluidized bed presented above, that contains 2000 spheres only, requires a mesh comprises at least 50 million of cells (i.e. at least 24 cells over the particle diameter) and runs over a long physical time to acquire converged statistics.This translates into weeks of computing on, e.g., 256 cores of a supercomputer.Highly scalable codes with fast algorithms are hence sheer game-shifters.Very soon, computations with a billion of cells will not be uncommon anymore to gain even more understanding into the dynamics of suspension flows.The main limitation of micro-scale models is the local resolution.In other words, the smallest length scale resolved by the model is controlled by the grid size.This creates intrinsic limitations in terms of what can be accurately capture as, e.g., boundary layers or lubrication force.For dense suspensions, most PRS are not able to fully captured lubrication forces and hence rely on analytical corrections, that exist for spheres only.In dilute and semidilute flows, the upper bound in terms of particle Reynolds number is generally around 500, as above 500 the boundary layer is too thin to be properly captured by a constant grid size mesh, unless at least 100 cells are used over the particle diameter.Similarly, the extension to non-spherical and presumably complex shape puts even more stress on the fineness of the mesh.One remedy to this issue of local mesh size is to implement these micro-scale models with an Adaptive Mesh Refinement (AMR) technique as illustrated in Fig. 5.This would offer the opportunity to refine the mesh close to the particle surface without increasing the total number of cells.One interesting problem that one could attack with such a tool is a dense suspension of nonspherical particles, for which there is no lubrication force correction available in the literature and hence lubrication force would be fully computed/captured by the solution of the coupled Navier-Stokes/Newton's equations.
At the meso scale level, we have shown how to improve TWEL models based on PRS data.Our impression is that the accuracy of TWEL models relies mainly on the projection operator (BC, GK, others ?) and on drag laws conceptually novel (with respect to the deterministic laws introduced over the past 50 years).This is crucial to recover the right level of particles motion.This is a rewarding research path to follow.

Figure 2 :
Figure 2: Flow configuration and flow field with the two models.
(a).It contains 2000 monodisperse spheres and has a dimensionless size of 10 × 10 × 35.We impose periodic boundary conditions on the lateral walls and a homogeneous and constant in time upwards inlet velocity at the bottom wall.Dimensionless parameters are set to ρ r = 10, u * in /u * m f = 3 (where u * m f is the minimum fluidization velocity) and Re in = 6.