Improvement of the Viscous Penalty Method for Particle-Resolved Simulations

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 research area that is commonly found in the fields of environment and industrial processes. Among them, we can cite fluidized beds and chemical engineering, material 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 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 [13] [14]. We choose to use ITPM instead of Darcy penalty method [15] because ITPM was demonstrated to be second order convergence in space [13] whereas Darcy penalty is only first order [16]. In addition, ITPM is a more general approach allowing dealing with moving particles, which is our objective in future works.
The paper is organized as follows. In Section 2, the main features of the viscous penalty method are presented and discussed. In particular, a new definition of the solid phase function located at the off-diagonal viscosity coefficients is proposed. The uniform Stokes flow past a cylinder is considered in Section 3. Various numerical parameters such as the numerical diameter of the particle, the penalty viscosity, the augmented Lagrangian parameter or the solid phase function evaluation are studied. At the end, the best set of parameters is proposed for an improved ITPM method, whose convergence order is almost 2. Section 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 simulation approach. It is compared to Ergun correlation [17] for various solid volume fractions. Conclusions and perspectives are finally drawn in Section 5.

Fictitious Domain Approach
The simulation of solid particles interacting with a carrier fluid is difficult to im- 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 [18]. The concept that separates the particle interfaces and the mesh used to solve the conservation equations is called fictitious domain approach [15] [19]. 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 [13], incompressible two-phase flows involving a carrier fluid and a solid particle phase can be modeled on a fixed mesh with fictitious domain approaches by considering the incompressible Navier-Stokes equations together 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 0.5 C = . As explained by Kataoka [20] for fluid/fluid two-phase flows and Vincent [13] 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 read where u is the velocity in all phases (fluid and solid), p the pressure, t the time, g the gravity vector, ρ 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 si F [19] [21] which is not considered in the present work as only fixed particles are dealt with. The source term m F is used to impose a flow rate to the fluid. In the present work, only fixed particles are considered, so Equation (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 ( ρ and µ ) depend on C. They will be discussed later on in the present work.

Penalty Methods
As previously explained, the one-fluid model and the fictitious domain approach formulated 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, which belongs 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). The first publication that reports on this approach was by Saulev [22].
For fixed particles, various improvements were suggested based on Darcy and Volume penalty methods [15], [16], [23]. Concerning moving particles, the viscous penalty method of the first order of convergence in space was initially proposed by Ritz and Caltagirone [24]. The method was then improved by [12] [13] [25] [26] 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 1 C = requires defining a specific rheological law for the rigid fluid part without imposing the velocity.
As reported by [13] the solid constraint is intrinsically maintained if the deformation tensor is nullified in the solid sub-domain s Ω : For the resolution of the momentum conservation Equation (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 these 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 . Following the work of Caltagirone and Vincent [12], 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 decomposition is then to act distinctly on each term in order to strongly impose the associated stress. If we assume that the Navier-Stokes equa- Decomposing  according to the partial derivative of the velocity in Cartesian coordinates for the sake of simplicity, we obtain [12] 0 0 0 0 This decomposition is written in a compact form as where Λ is the tearing tensor, Θ is the shearing tensor and Γ is the rotation tensor.
Consequently, the divergence of the viscous stress tensor for a Newtonian fluid appearing in the one-fluid model (2) reads The main interest of formulation (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, (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 [13], a staggered grid (see ation. 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 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 [13]. In a given pressure cell, 10 test points are seeded in each direction, as illustrated in Figure 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% [13].
In our second order convergence penalty approach, a phase function C µ located at the viscous mesh nodes (white nodes in Figure 1) is introduced. As in [13], it can be interpolated from C: 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, instead of interpolating it. The effect of this improvement is studied in this paper.

Local Properties of the Equivalent Fluid
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 pioneering work of Fortin and Glowinski [27], an augmented Lagrangian method is applied to the unsteady Navier-Stokes equations dedicated   (9) where r is an augmented Lagrangian penalty parameter used to impose the incompressibility constraint, m is an iterative convergence index and AL  a nu- 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 conditioning of linear system is degraded while 0 r → does not act on the fluid constraint and keeps the conditioning of the matrix unchanged. As recommended by [27], 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 [27]. From numerical experiments, optimal values are found to be of the order of i ρ and i µ in each phase (fluid or solid) to accurately solve the motion equations in the related zones [26] [28]. Algebraic improvements have also been proposed by Vincent [29] 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 conditioning 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.

Discretization Schemes and Solvers
All the schemes and solvers utilized in the present work are presented and discussed in detail in [13]. 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 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 ( ) 1 n t + ∆ , except the non-linear inertial term that is linearized with a second order Adams-Bashforth scheme as follows ( ) The obtained linear system can be solved in three-dimensions with a BiCGSTAB II iterative solver [30], preconditioned under a Modified and Incomplete LU approach [31] to speed-up the convergence of the solver. In this work, direct MUMPS solver [32], [33] is preferred as it provides computer error residuals. All the code is working on massively parallel computers by using MPI devices and exchanges [13].

Uniform Stokes Flow past a Cylinder
A validation of the presented method and a numerical study of some of its parameters are conducted considering the steady uniform Stokes flow past an isolated cylinder. The analytical solution is illustrated in Figure 3. According to [34] [35], a uniform Stokes flow ( where

Simulations Setup
The computational domain used to simulate a uniform Stokes flow past a cylinder is a square of a Length   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 simulation 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 numerical radius as follows: [ ]     being modified to a larger extent. A first conclusion here is that choosing harmonic or discontinuous averages is more desirable as R n 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 [13]. 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. To answer this question, the same study is conducted on R n and average viscous laws by considering the C µ directly calculated on the viscous points without using the pressure nodes. Figure 6 shows the velocity L1 relative error in the whole domain (14) for Stokes flow past a cylinder for the three viscous laws discussed above but with the color function C µ computed on the viscous mesh instead of interpolating it from the C function on the pressure nodes as in previous simulations. One can observe that the minimum error is reached for a numerical radius 2 2 whereas the new C µ computation seems to have no influence on the discontinuous law results. On the other hand, for the harmonic average law, not only the minimum error is divided by 10 but also this error is reached for the physical diameter of the cylinder 2 n d R = . Therefore, and for the rest of this work, the color function C µ will always be computed by projecting the particle shape on the viscous mesh, together with the use of the harmonic average law to compute the viscosity in the Eulerian mesh containing the interface. This important conclusion is new and has never been obtained in previous penalty simulations of particle flows [13] [17] [36] [37].
A new simulation is carried out with a new set of parameter and the conclusion of the numerical study above. They are given by:  Viscous law: harmonic average law instead of arithmetic law.   The Lagrangian parameter 5 10 r = remains the same.
With this new set of numerical parameters, Figure 7 shows the huge improvement brought by the new set of numerical parameters on relative error for the velocity and the pressure. Indeed the error decreases from 100% to less than 5% for the pressure in the fluid domain, except in the cells containing the interface as illustrated in Figure 4(c). If we refer to Figure 4(a) and Figure 4(b), the error went from 50% to 10% in the fluid region near the cylinder and from 10% to less than 2% in the rest of the fluid domain for both velocity components.

Effect of the Viscosity Ratio and the Augmented Lagrangian Parameter r
The   The last numerical parameter to be studied in this work is the Lagrangian parameter r. Indeed, the augmented Lagrangian method is a kind of penalty technique, and the incompressibility is imposed when r → +∞ . Therefore, knowing from which value of r the solution does no longer depend on it is an important matter to be carefully studied. Indeed, the larger r is, the worse is the solving of the linear system. As a consequence, r has to be large to impose incompressibility and at the same time the smallest possible to keep the conditioning of the linear system as small as possible too. Figure 9 shows the velocity L1 relative error in the whole domain (14) for Stokes flow past a cylinder for a Lagrangian parameter r between 10 3 and 10 9 .
One can observe that the solution is augmented Lagrangian parameter independent for 5 10 r ≥ . This is the value that will be used in the rest of this work. Note that in this work, the resolution of the linear system is ensured by a direct solver, which allows us to use a large value of r. On the other hand, the use of an iterative solver can be difficult in the case of    The Lagrangian parameter is kept as 5 10 r = . Figure 10 shows the relative error for the velocity and the pressure between the results of the penalty simulation with the best set of parameters and the analytical solution. It can be seen that error is now lower than 1% for either velocity or pressure in the fluid area far from the particle, and about 10% in the region containing the interface. This is mainly due to the one-fluid model for which the

Order of Convergence
Given the fact that we have been able to find a satisfactory set of parameters to obtain an accurate result on velocity and pressure, as illustrated in Figure 10, a study of convergence order of the viscous penalty method is conducted by simulating a series of uniform flow past a cylinder using the best set of parameters and by changing the Eulerian mesh resolution using different d x ∆ .     ) is given in Figure 13.

Uniform Flow past a Square Configuration of Cylinders
The aim of this section is to validate the superficial mean fluid velocity where u is solution of the Navier-Stokes equation using the viscous penalty method (ITPM) with the best set of parameters proposed in the previous section.    (15) This relation is a generalization of the Forchheimer equation [39].
Given the non-dimensional friction factor p f is defined as: (16) and the Reynolds number Re is given by: Re The validation consists in  simulating a uniform flow past a square configuration of cylinders, for a given pressure drop P L ∆ and a solid volume fraction d α .
 extracting from the velocity field the superficial mean fluid velocity f u .  computing the friction factor p f using (16).
 comparing p f to the Ergun correlation [40] given by (17).

Conclusions and Suggestions
The Viscous Penalty Method ITPM [13] has been used to simulate two-dimensional fixed particulate flows. The first goal of the work was to set up the best linder. The viscosity ratio between the fluid and the penalty viscosity s µ inside the particle, the augmented Lagrangian parameter r, the viscous law, the solid fraction estimate C µ at the viscous nodes and the numerical radius of the particle were investigated. For the first time, we have been able to demonstrate that if C µ is directly calculated by projecting the real shape of the particle on the viscous nodes, the numerical radius of the particle n R does not have to be adapted compared to its real physical value. Moreover, the best accuracy is obtained when a harmonic law on the viscosity is used to build the equivalent properties of the one-fluid model in cells cut by the fluid/particle interface.
Concerning the penalty viscosity, imposing 1000 times the fluid velocity is the best compromise between error level and solving efficiency. To finish with setting of ITPM parameters, Ongoing works are developed in several directions:  The ITPM is used to extract the drag and lift force coefficient for various arrangements of spherical particles [14].  The ITPM is extended to heat transfers in particulate flows. As for the force coefficient, the heat transfer coefficient is extracted for any particle inside various arrangements of spheres [41].  The viscous penalty method is utilized to simulate the force exerted by an incompressible flow on ellipsoidal particles as well as heat transfer coefficients [42].