Project topics and materials
Get Complete Project Material File(s) Now! »
Improvement of the Viscous Penalty Method for particle-resolved simulations
This chapter is the article [26] authored by M.-A. Chadil, S. Vincent and J.-L. Estivalezes
Abstract
A numerical study of the parameters controlling the viscous penalty method is investigated to better set up Particle-Resolved Direct Numerical Simulations (PR-DNS) of particulate flows. Based on this analysis, improvements of the methods are proposed in order to reach an almost second order convergence in space. The viscous penalty method is validated in Stokes regime by simulating a uniform flow past a fixed isolated cylinder. Moreover, it is also utilized in moderate Reynolds number regime for a uniform flow past a square configuration of cylinder and compared in terms of friction factor to the well-known Ergun correlation.
Introduction
The motion of rigid particles interacting with a carrier fluid is a very active re-search area that is commonly found in the fields of environment and industrial processes. Among them, we can cite fluidized beds and chemical engineering, ma-terial manufacturing and design, sand dynamics, beach erosion under wave impact or nano-particle impact on human health. The simulation of such real problems is based on the use of Eulerian-Eulerian or Eulerian-Lagrangian models that require knowledge of constitutive laws for drag, lift, torque, collisions or heat transfers for the fluid-particle interactions. One way of designing these laws or validating them is to use resolved-scale particle approaches, in which all scales associated with the fluid flow and the hydrodynamic forces on the particle are directly simulated, unlike in point-particle or Eulerian-Eulerian approaches where drag and lift correlations are required a priori to simulate the problem.
The numerical simulation of resolved-scale particle motion is a highly developed field of research mainly based on fixed structured grids, as unstructured meshes adapted to the particle motion are difficult to design in three dimensions and CPU time consuming [62]. Among the wide variety of fictitious domain approaches, i.e. particles are treated as immersed interfaces on a fixed mesh, we can cite the nu-merical methods based on Lattice Boltzmann models [12, 51, 53, 64, 65] and the approaches that uses the Navier-Stokes equations, such as the Immersed Bound-ary Method (IBM) of Uhlmann [105, 106], the PURe-IBM approach of Tenneti et al. [103], the Distributed Lagrangian Method (DLM) of Glowinski and co-workers [78, 79] and the Implicit Tensorial Penalty Method (ITPM) of Vincent et. al. [20, 109], also called viscous penalty method.
In the present work, we choose to investigate viscous penalty methods on fixed Cartesian grids for fixed particles. Compared to other fictitious domains techniques, the main interest of penalty methods is to rely on fully coupled velocity solving with incompressible and solid constraints satisfaction instantly, thanks to an augmented Lagrangian method for the fluid and viscous penalty for the solid phase. Our main goal is first to characterize the accuracy and convergence order of the ITPM method on reference particle motion test cases but also to improve the numerical method and the setting of numerical penalty parameters, what has never been done. The reference method from which we left is published in [24, 109]. We choose to use ITPM instead of Darcy penalty method [60] because ITPM was demonstrated to be second order convergence in space [109] whereas Darcy penalty is only first order [6]. In addition, ITPM is a more general approach allowing to deal with moving particles, that is our objective in future works.
The paper is organized as follows. In Section 2.2, the main features of the vis-cous penalty method are presented and discussed. In particular, a new definition of the solid phase function located at the off-diagonal viscosity coefficients is pro-posed. The uniform Stokes flow past a cylinder is considered in section 2.3. Various numerical parameters such as the numerical diameter of the particle, the penalty viscosity, the augmented Lagrangian parameter or the solid phase function evalua-tion are studied. At the end, the best set of parameters is proposed for an improved ITPM method, whose convergence order is almost 2. Section 2.4 is devoted to the uniform flow past a square configuration of cylinders. With the previous best set of parameters of ITPM, the friction factor is calculated with our particle-resolved sim-ulation approach. It is compared to Ergun correlation [42] for various solid volume fractions. Conclusions and perspectives are finally drawn in section 2.5.
Model and numerical methods
Fictitious domain approach
The simulation of solid particles interacting with a carrier fluid is difficult to imple-ment with unstructured meshes in particular with 3D geometries. The commonly developed alternative approach consists in simulating this kind of flow on a fixed mesh not adapted to the shape of the particle, i.e. by considering a solid phase fraction, and to locate the fluid-solid interface thanks to an auxiliary phase function such as the Volume Of Fluid or the Level Set [95]. The concept that separates the particle interfaces and the mesh used to solve the conservation equations is called fictitious domain approach [60, 87]. Indeed, from the motion equation point of view, the interface is not known, only the presence of the solid phase is taken into account into the motion conservation equations thanks to a volume auxiliary function and associated specific forcing terms.
One-fluid model
As previously presented in [109], incompressible two-phase flows involving a car-rier fluid and a solid particle phase can be modeled on a fixed mesh with fictitious domain approaches by considering the incompressible Navier-Stokes equations to-gether with a phase function C describing the particle phase shape. By definition, the phase function C equals to 1 in the solid phase and 0 in the fluid medium. The fluid-solid interface is located by the isosurface C = 0.5. As explained by Kataoka [59] for fluid/fluid two-phase flows and Vincent [109] for particle flows, the resulting one-fluid model takes implicitly into account the coupling between different phases separated by resolved interfaces, i.e. the particles are larger than the mesh cell size. The motion equations reads
fl 3ˆˆut + (u • Ò) u4
Ò • u = 0 (2.1)
= ≠Òp + flg + Ò • ˵(Òu + Òtu)È + Fsi + Fm (2.2)
ˆC + u • ÒC = 0 (2.3)
where u is the velocity in all phases (fluid and solid), p the pressure, t the time, g the gravity vector, fl and µ respectively the density and the dynamic viscosity of the equivalent fluid. The four-way coupling between particles and fluid motions is ensured in the momentum equations by the presence of a solid interaction force Fsi [17, 74] which is not considered in the present work as only fixed particles are dealt with. The source term Fm is used to impose a flow rate to the fluid. In the present work, only fixed particles are considered, so equation (2.3) will be discarded.
The one-fluid model is almost identical to the classical incompressible Navier-Stokes equations, except that the local properties of the equivalent fluid (fl and µ) depend on C. They will be discussed later on in the present work. In the present form, equations (2.1-2.3) do not account for incompressibility and solid constraints. Satisfying these mechanical properties requires developing specific numerical meth-ods called penalty approaches. They are detailed in the next section 2.3.
As previously explained, the one-fluid model and the fictitious domain approach for-mulated to deal with particle flows require to consider each different phase (fluid, solid) as a fluid medium with specific material properties (density and viscosity for an isothermal flow). The domain is covered by a set of representative elementary volumes, i.e. the mesh cells on a numerical point of view, that are belonging to different sub-domains located by the phase function C. A way to satisfy fluid and solid constraints is to define penalty terms in the momentum equation (2.2). The first publication that reports on this approach was by Saulev [96]. For fixed par-ticles, various improvements were suggested based on Darcy and Volume penalty methods [6, 60, 63]. Concerning moving particles, the viscous penalty method of first order of convergence in space was initially proposed by Ritz and Caltagirone [91]. The method was then improved by [20, 87, 109, 110] to become a second order in space penalty method called ITPM. This method is detailed in the rest of this section and will be used in the present work.
Ensuring the solid behavior in the solid zones where C = 1 requires defining a specific rheological law for the rigid fluid part without imposing the velocity. As reported by [109] the solid constraint is intrinsically maintained if the deformation tensor is nullified in the solid sub-domain Ωs: ’P œ Ωs, Òu + ÒT u = 0 (2.4)
For the resolution of the momentum conservation equation (2.2) in the Navier-Stokes equations, this condition is asymptotically verified when µ æ +Œ. In other words, viscous penalty method consists in imposing large values of viscosity in the particles compared to the fluid viscosity to implicitly impose the solid behavior and also the coupling between fluid and solid. For fixed particles, the velocity of the Eulerian cells near the centroid of the particle is assumed to be zero. A Darcy penalty method is utilized to satisfy this conditions. The viscous penalty method is used in the rest of the solid particles. Indeed, it propagates the zero velocity in the whole solid medium. The effect of the ratio between the particles and the fluid viscosities will be studied in this work.
A specific model is designed for handling the solid particle behavior in the one-fluid Navier-Stokes equations. It is based on a decomposition of the strain tensor ‘ = Òu + ÒT u. Following the work of Caltagirone and Vincent [20], the strain tensor can be reformulated so as to distinguish several natural contributions of the strain tensor dealing with tearing, shearing and rotation. The interest of this de-composition is then to act distinctly on each term in order to strongly impose the associated stress. If we assume that the Navier-Stokes equations for a Newtonian fluid contain all physical contributions traducing shearing or pure rotation effects, the splitting of the viscous stress tensor allows to impose separately these contri-butions by modifying the orders of magnitude of each term, through the related viscosity coefficients. These penalty terms act directly in the motion equations and so ensure the coupling between the fluid and the solid part of the simulation domain
The main interest of formulation (2.7) is to dissociate stresses operating in a viscous flow and then to make the implementation of a numerical penalty method easier. For instance, in a solid phase, if µ is chosen larger than the surrounding fluid viscosity, (2.7) imposes that the local solid flow admits no shearing, no tearing and a constant rotation according to the surrounding flow constraints. These flow constraints are implicitly transmitted to the particle sub-domain as they are solved with the fluid motion at the same time. In the same way, the modifications of the flow motion by the particle movement are directly accounted for (two-way coupling).
For obtaining a second order convergence in space [109], a staggered grid (see Figure2.1) is needed to implement this strain tensor decomposition where the tear-ing viscosity µt = 2µ is located at the pressure nodes whereas the pure shearing µsh = 2µ and rotation µr = µ viscosities lie on a specific grid, at the center of the mesh grid cells. Defining µ in the solid 2 to 3 orders of magnitude larger than the fluid velocity is equivalent to having µt, µsh and µr tending to large values and so acting as viscous penalty terms in the motion equation. In these grid cells, the local medium will be almost solid.
Phase function
The phase function C located at pressure nodes is automatically built by projecting particles onto the pressure mesh (black nodes in Figure 2.1). The color function is defined as the amount of solid in a pressure cell, i.e. the local solid fraction. Therefore, in the cells containing the interface, C is computed thanks to virtual test points [109]. In a given pressure cell, 10 test points are seeded in each direction, as illustrated in Figure 2.2. By counting the number of test points belonging to the particle and dividing this number by the total number of test points, the solid fraction C is naturally obtained. It has been previously demonstrated that using 10 points by directions provides an error on C lower than 1% [109].
In our second order convergence penalty approach, a phase function Cµ located at the viscous mesh nodes (white nodes in Figure 2.1) is introduced. As in [109], it can be interpolated from C:
1 ÿ
Cµ = (2.8)
CN
where N denotes the indices of the pressure nodes located at the vertices of the cell to which Cµ belongs.
Alternatively, a projection of the particle on the viscous mesh is proposed in this work to provide the phase function Cµ by using test points, as presented in figure 2.2, instead of interpolating it. The effect of this improvement is studied in this paper.
Local properties of the equivalent fluid
On a discrete point of view, the flow grid cells cut by the fluid-solid interface must be distinguished compared to those entirely included in the particles or in the fluid. Different methods can be designed to define the homogenized viscosity µ in these mixed cells. Three different numerical viscous laws have been investigated according to the fluid and solid viscosities (µf and µs respectively), C for the diagonal viscous stress tensor terms, Cµ for the offdiagonal viscous contributions and also a conditional indicator function IC satisfying IC<0.5 = 1 if C < 0.5 or ICØ0.5 = 1 if C Ø 0.5:
1. Discontinuous law:
µ = [µf IC<0.5 + µsICØ0.5]
2. Arithmetic law:
µ = [(1 ≠ C)µf + Cµs]
3. Harmonic law: D
C = µf µs
µ Cµf + (1 ≠ C)µs
In the previous laws, C can be replaced by Cµ if the viscosity is located at shearing µsh or pure rotations µr nodes. Concerning the density, an arithmetic average is used whatever its location on the discretization grid. The effect of the choice of the viscosity average law is studied in this work.
Augmented Lagrangian method
Following the pioneer work of Fortin and Glowinski [45], an augmented Lagrangian method is applied to the unsteady Navier-Stokes equations dedicated to particulate flows. It allows dealing with the coupling between the velocity and pressure and to satisfy the fluid and solid constraints at the same time by solving a saddle point problem. Starting with uú,0 = un and pú,0 = pn, the augmented Lagrangian solu-tion reads
while ||Ò• uú,m|| > ‘AL , solve
A (uú,0, pú,0) = (un, pn) Ò • u ) (2.9)
Δt + u • Òu B≠Ò(
fl uú,m ≠ uú,0 ú,m≠1 ú,m r ú,m = ≠Òpú,m≠1 + flg + Ò • [µ(Òuú,m + ÒT uú,m)] + Fsi pú,m = pú,m≠1 ≠ rÒ • uú,m
where r is an augmented Lagrangian penalty parameter used to impose the in-compressibility constraint, m is an iterative convergence index and ‘AL a numerical threshold controlling the constraint. The augmented Lagrangian method is a kind of penalty technique: if r æ +Œ, the incompressibility is imposed but the solving of the linear system is difficult with iterative solvers as the conditionning of linear system is degraded while r æ 0 does not act on the fluid constraint and keeps the conditioning of the matrix unchanged. As recommended by [45], a constant value of r is used, for example equal to the average between the minimum and maximum eigenvalues of the linear system for Stokes flows [45]. From numerical experiments, optimal values are found to be of the order of fli and µi in each phase (fluid or solid) to accurately solve the motion equations in the related zones [108, 110]. Algebraic improvements have also been proposed by Vincent [111] to automatically estimate the local values of r . In the present work, an automatic algebraic estimate of r will be used to optimize as much as possible the conditionning of the linear system while maintaining expected incompressible and solid constraints in the related zones. The effect of the Lagrangian parameter r is considered in the following section 2.3.
Discretization schemes and solvers
All the schemes and solvers utilized in the present work are presented and discussed in detail in [109]. The mass and momentum conservation equations, containing the viscous and augmented Lagrangian penalty terms, are discretized with implicit Finite volumes on structured staggered meshes (see figure 2.1). The time derivative is approximated with a second order Euler scheme while the inertial, viscous and augmented Lagrangian terms are discretized with second-order centered schemes. All fluxes are written at time (n + 1)Δt, except the non-linear inertial term that is linearized with a second order Adams-Bashforth scheme as follows u • Òu ¥ 12un ≠ un≠12 • Òun+1 (2.10)
The obtained linear system can be solved in three-dimensions with a BiCGSTAB II iterative solver [48], preconditionned under a Modified and Incomplete LU approach [39] to speed-up the convergence of the solver. In this work, direct MUMPS solver [4, 5] is preferred as it provides computer error residuals. All the code is working on massively parallel computers by using MPI devices and exchanges [109].
Uniform Stokes flow past a cylinder
A validation of the presented method and a numerical study of some of its param-eters are conducted considering the steady uniform Stokes flow past an isolated cylinder. The analytical solution is illustrated in Figure 2.3. According to [13, 19], a uniform Stokes flow (Re = 10≠3) past a cylinder of diameter d = 2m, with the undisturbed velocity being noted UŒ = 1m/s, is solution of the Brinkman equation ≠Òp + µΔui ≠ Kµ ui = 0. The reference solution is given in polar coordinate frame (r,◊ ), centered on the particle, by:
uú (rú, ◊) = Y ú 1≠ 1 2K1(⁄2) ú 1 2 2 K1(⁄rú)
] 1 1 + 2K1 (⁄) 1 + rú + 2 K1(⁄rú) cos ◊
r ⁄K0(⁄) r ⁄K0 (⁄)
1 1 + 1 1 + 2 ≠ K0(⁄rú) + 22 sin ◊
⁄K0(⁄) (rú)2 K0(⁄) ⁄rú
[ ≠ 3 ≠ (⁄) 1 rú ≠ rú 4 cos ◊ (2.11)
pú (rú, ◊) = Re⁄2 3 1+ ⁄K0 4 (2.12)
2 2K1 (⁄) 1
where uú = u , pú = p , rú = 2r , fl = 1kg.m≠3 is the fluid density, ⁄ = d2 U flU 2 d 4K is the dimensionless permeability of the porous medium in Brinkman sens, K is the permeability of the inside and outside the porous cylinder, K0 and K1 are the modified Bessel functions of rank 0 and 1. For K æ 0, the porous cylinder can be likened to an impermeable solid particle whereas outside the cylinder, K æ +Œ to obtain a fluid behavior.
Simulations setup
The computational domain used to simulate a uniform Stokes flow past a cylinder is a square of a Length L = 2d, and the spatial discretization, using a regular Carte-sian grid called Eulerian mesh, is represented by the number of grid cells across the diameter of the particle Δdx = 20. The velocity and pressure exact solutions ((2.11), (2.12) respectively) for a Stokes flow past a cylinder were taken as initial condition, as illustrated in Figure 2.3. They were also implemented at boundary conditions as a Dirichlet condition to be able to simulate such a flow in a numerically small domain not extending to infinity as Stokes flow would require.
A first simulation of a uniform flow past a cylinder is carried out using a reference set of parameters presented bellow:
• Viscous law: Arithmetic average law is chosen for this simulation.
Figure 2.3: Exact solution for a Stokes flow past a fixed cylinder: the first velocity component field u1 (a), the second velocity component field u2 (b) , and the pressure field p (c), are plotted at each point of the domain (fluid and solid).
• Numerical radius: it has been previously mentioned [109] that on a numerical point of view, the cylinder radius has to be tuned according to its physical radius. Indeed, interpolations are used in the cells cut by the fluid-solid interface for viscous discrete nodes, inducing numerical variations of the solid phase compared to the real one. In our simulations, the numerical radius Rn of the cylinder is given by: Rn = d2 + eΔ16x e is a correction coefficient on Rn. It is imposed to be e = 0, i.e. Rn = d2 , so that Rn is the physical radius of the cylinder for this simulation.
• Computation of Cµ: For the first simulation, it is interpolated from C, known in the pressure mesh, on the viscous mesh.
• The viscosity ratio µs between the viscosity imposed in the Eulerian cells µf inside the cylinder µs and the fluid viscosity µf is chosen such that µs = 500 µf for this simulation.
• The Lagrangian parameter is r = 105.
Figure 2.4 shows the relative error in each point of the domain for the velocity
Error = I |uSimu≠uAnalytic| u
if uAnalytic ”= 0 (2.13)
| u |uAnalytic|
Simu| if Analytic = 0
and the pressure between the simulation results and the analytical solution given by (2.11) and (2.12). This error is about 100% for the pressure in the fluid domain as illustrated in Figure 2.4(c) and more than 50% in the fluid region near the cylin-der and about 10% in the rest of the fluid domain for both velocity components as illustrated in Figure 2.4(a), Figure 2.4(b).
Facing this huge error for both pressure and velocity, we decided to conduct a numerical study on the effect of previously listed numerical parameters on the sim-ulation results. Our main goal is to set up the selection of parameters to minimize these errors. At the end of each study, the simulation results obtained with new parameters will be given in order to show the improvement made.
Sensitivity of simulations to viscous law, numerical radius Rn and phase function computation Cµ on the viscous mesh
The viscous law and the numerical radius are first investigated. To do so, several simulations are carried out with discontinuous, arithmetic and harmonic average laws (for both C and Cµ which is interpolated at this state) and for different nu-merical radius as follows: Rn = + e e œ J≠16, 16K
All other parameters remain unchanged: = 500 and r = 10 µf
Figure 2.5 shows the velocity L1 relative error in the whole domain = q | q |uAnalytic| (2.14) Error uSimu ≠ uAnalytic| for the Stokes flow past a cylinder for different viscous laws. It can be observed that the minimum error for arithmetic average law is reached for Rn = d2 ≠ Δx whereas it is reached for a numerical radius Rn = d2 + Δ8x for harmonic and discontinuous average laws. This minimum error is about 1% for both harmonic and discontinu-ous law whereas it is 0.5% larger for the arithmetic law with Rn being modified in a larger extent. A first conclusion here is that choosing harmonic or discontinuous averages is more desirable as Rn is closer to the physical cylinder radius and the obtained error is smaller Until now, the color function on the viscous mesh Cµ was interpolated from C computed on the pressure mesh [109]. One interesting issue is how the error implied by the different average laws will change if Cµ is computed directly on the viscous mesh by projecting the cylinder shape with the virtual point procedure presented before in figure 2.2. To answer this question, the same study is conducted on Rn and average viscous laws by considering the Cµ directly calculated on the viscous points without using the pressure nodes.
Table of contents :
1 Introduction
2 Improvement of the Viscous Penalty Method for particle-resolved mulations
2.1 Introduction
2.2 Model and numerical methods
2.3 Uniform Stokes flow past a cylinder
2.4 Uniform flow past a square configuration of cylinders
2.5 Conclusions and Suggestions
3 Accurate estimate of drag forces using particle-resolved direct numerical simulations
3.1 Introduction
3.2 Model and numerical methods
3.2.1 Fictitious domain approach
3.2.2 Penalty methods
3.2.3 Discretization schemes and solvers
3.3 Lagrangian extrapolation of forces for immersed boundary methods
3.3.1 Low order naive approach
3.3.2 New high order method based on Lagrange extrapolation
3.4 Validation on flows interacting with an isolated particle
3.4.1 Drag coefficient
3.4.2 Simulations setup
3.4.3 Study of numerical parameters for the Lagrange extrapolation
3.4.4 Result on the drag coefficient
3.4.5 Pressure coefficient
3.5 Forces in fixed arrangements of spheres
3.5.1 Monodispersed arrangements of spheres
3.5.2 Bidisperse arrangements of spheres
3.6 Conclusions and perspectives
3.7 Appendix 1: Taylor Interpolation
4 Accurate calculation of heat transfer coefficients for motions around particles with a finite-size particle approach
4.1 Introduction
4.2 Numerical Methodology
4.3 Convective heat transfer forced by a uniform flow around a stationary sphere
4.4 Face-Centered Cubic periodic arrangement of spheres
4.5 Finite size random arrangement of spheres in a channel
4.6 Conclusion
5 Drag, lift and Nusselt coefficients for ellipsoidal particles using particle-resolved direct numerical simulations
5.1 Introduction
5.2 Numerical Methodology
5.3 Drag, Lift and Nusselt for an isolated stationary ellipsoid past by a uniform flow
5.4 Conclusion
6 Novel method to compute drag force and heat transfer for motions around spheres
6.1 Introduction
6.2 Numerical Methodology
6.2.1 Conservation equations
6.2.2 Fictitious domain approach and viscous penalty method
6.2.3 Drag force and heat flux computation using Aslam extension
6.3 Isolated stationary sphere past by a uniform flow
6.3.1 Drag force computation
6.3.2 Heat transfer computation
6.4 Face-Centered Cubic arrangement of stationary sphere past by a uniform flow
6.4.1 Monodispersed Face-Centred Cubic periodic arrangement of spheres
6.5 Conclusions
7 Conclusions and perspectives
Bibliography
