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.

1. 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 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 [1] . Among the wide variety of fictitious domain approaches, i.e. particles are treated as immersed interfaces on a fixed mesh, we can cite the numerical methods based on Lattice Boltzmann models [2] [3] [4] [5] [6] and the approaches that uses the Navier-Stokes equations, such as the Immersed Boundary Method (IBM) of Uhlmann [7] [8] , the PURe-IBM approach of Tenneti et al. [9] , the Distributed Lagrangian Method (DLM) of Glowinski and co-workers [10] [11] and the Implicit Tensorial Penalty Method (ITPM) of Vincent et al. [12] [13] , 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 [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.

2. Model and Numerical Methods

2.1. Fictitious Domain Approach

The simulation of solid particles interacting with a carrier fluid is difficult to implement 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 [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.

2.2. 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 $C=0.5$ . 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

$\nabla \cdot u=0$ (1)

$\rho \left(\frac{\partial u}{\partial t}+\left(u\cdot \nabla \right)u\right)=-\nabla p+\rho g+\nabla \cdot \left[\mu \left(\nabla u+{\nabla }^{t}u\right)\right]+{F}_{si}+{F}_{m}$ (2)

$\frac{\partial C}{\partial t}+u\cdot \nabla C=0$ (3)

where $u$ is the velocity in all phases (fluid and solid), p the pressure, t the time, $g$ the gravity vector, $\rho$ and $\mu$ 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 ${F}_{si}$ [19] [21] which is not considered in the present work as only fixed particles are dealt with. The source term ${F}_{m}$ 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 ( $\rho$ and $\mu$ ) depend on C. They will be discussed later on in the present work. In the present form, Equations (1)-(3) do not account for incompressibility and solid constraints. Satisfying these mechanical properties requires developing specific numerical methods called penalty approaches. They are detailed in the next section.

2.3. 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 $C=1$ 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 ${\Omega }^{s}$ :

$\forall P\in {\Omega }^{s},\nabla u+{\nabla }^{\text{T}}u=0$ (4)

For the resolution of the momentum conservation Equation (2) in the Navier-Stokes equations, this condition is asymptotically verified when $\mu \to +\infty$ . 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 tensor $\stackrel{¯}{\stackrel{¯}{ϵ}}=\nabla u+{\nabla }^{\text{T}}u$ . 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 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 contributions 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 instantaneously.

Decomposing $\stackrel{¯}{\stackrel{¯}{ϵ}}$ according to the partial derivative of the velocity in Cartesian coordinates for the sake of simplicity, we obtain [12]

$\stackrel{¯}{\stackrel{¯}{ϵ}}=2\left[\begin{array}{ccc}\frac{\partial u}{\partial x}& 0& 0\\ 0& \frac{\partial v}{\partial y}& 0\\ 0& 0& \frac{\partial w}{\partial z}\end{array}\right]+2\left[\begin{array}{ccc}0& \frac{\partial u}{\partial y}& \frac{\partial u}{\partial z}\\ \frac{\partial v}{\partial x}& 0& \frac{\partial v}{\partial z}\\ \frac{\partial w}{\partial x}& \frac{\partial w}{\partial y}& 0\end{array}\right]-\left[\begin{array}{ccc}0& \frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}& \frac{\partial u}{\partial z}-\frac{\partial w}{\partial x}\\ \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}& 0& \frac{\partial v}{\partial z}-\frac{\partial w}{\partial y}\\ \frac{\partial w}{\partial x}-\frac{\partial u}{\partial z}& \frac{\partial w}{\partial y}-\frac{\partial v}{\partial z}& 0\end{array}\right]$ (5)

This decomposition is written in a compact form as

${\stackrel{¯}{\stackrel{¯}{ϵ}}}_{ij}=2{\Lambda }_{ij}+2{\Theta }_{ij}-{\Gamma }_{ij}$ (6)

where $\Lambda$ is the tearing tensor, $\Theta$ is the shearing tensor and $\Gamma$ is the rotation tensor.

Consequently, the divergence of the viscous stress tensor for a Newtonian fluid appearing in the one-fluid model (2) reads

$\nabla \cdot \left(\mu \left(\nabla u+{\nabla }^{t}u\right)\right)=\nabla \cdot \left[{\mu }_{t}\Lambda \left(u\right)\right]+\nabla \cdot \left[{\mu }_{sh}\Theta \left(u\right)\right]-\nabla \cdot \left[{\mu }_{r}\Gamma \left(u\right)\right]$ (7)

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 $\mu$ 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 Figure 1) is needed to implement this strain tensor decomposition where the tearing viscosity ${\mu }_{t}=2\mu$ is located at the pressure nodes whereas the pure shearing ${\mu }_{sh}=2\mu$ and rotation ${\mu }_{r}=\mu$ viscosities lie on a specific grid, at the center of the mesh grid cells. Defining $\mu$ in the solid 2 to 3 orders of magnitude larger than the fluid velocity is equivalent to having ${\mu }_{t}$ , ${\mu }_{sh}$ and ${\mu }_{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.

2.4. 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

Figure 1. Discrete interpretation of the split viscous penalty approach on staggered grids: (●) pressure points, arrows (►, ▲) for velocity components and (δ) for pure shearing and rotations viscosities. The black line represents the interface between a particle and the carrier fluid.

Figure 2. Virtual points (♦) on a pressure cell in a staggered grid.

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}_{\mu }$ located at the viscous mesh nodes (white nodes in Figure 1) is introduced. As in [13] , it can be interpolated from C:

${C}_{\mu }=\frac{1}{4}\underset{N}{\sum }{C}_{N}$ (8)

where N denotes the indices of the pressure nodes located at the vertices of the cell to which ${C}_{\mu }$ belongs.

Alternatively, a projection of the particle on the viscous mesh is proposed in this work to provide the phase function ${C}_{\mu }$ by using test points, as presented in Figure 2, instead of interpolating it. The effect of this improvement is studied in this paper.

2.5. 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 $\mu$ in these mixed cells. Three different numerical viscous laws have been investigated according to the fluid and solid viscosities ( ${\mu }_{f}$ and ${\mu }_{s}$ respectively), C for the diagonal viscous stress tensor terms, ${C}_{\mu }$ for the off-diagonal viscous contributions and also a conditional indicator function ${I}_{C}$ satisfying ${I}_{C<0.5}=1$ if $C<0.5$ or ${I}_{C\ge 0.5}=1$ if $C\ge 0.5$ :

1) Discontinuous law:

$\mu =\left[{\mu }_{f}{I}_{C<0.5}+{\mu }_{s}{I}_{C\ge 0.5}\right]$

2) Arithmetic law:

$\mu =\left[\left(1-C\right){\mu }_{f}+C{\mu }_{s}\right]$

3) Harmonic law:

$\mu =\left[\frac{{\mu }_{f}{\mu }_{s}}{C{\mu }_{f}+\left(1-C\right){\mu }_{s}}\right]$

In the previous laws, C can be replaced by ${C}_{\mu }$ if the viscosity is located at shearing ${\mu }_{sh}$ or pure rotations ${\mu }_{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.

2.6. 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 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}^{\ast ,0}={u}^{n}$ and ${p}^{\ast ,0}={p}^{n}$ , the augmented Lagrangian solution reads while $‖\nabla \cdot {u}^{\ast ,m}‖>{ϵ}_{AL}$ , solve

$\begin{array}{l}\left({u}^{\ast ,0},{p}^{\ast ,0}\right)=\left({u}^{n},{p}^{n}\right)\\ \rho \left(\frac{{u}^{\ast ,m}-{u}^{\ast ,0}}{\Delta t}+{u}^{\ast ,m-1}\cdot \nabla {u}^{\ast ,m}\right)-\nabla \left(r\nabla \cdot {u}^{\ast ,m}\right)\\ =-\nabla {p}^{\ast ,m-1}+\rho g+\nabla \cdot \left[\mu \left(\nabla {u}^{\ast ,m}+{\nabla }^{\text{T}}{u}^{\ast ,m}\right)\right]+{F}_{si}\\ {p}^{\ast ,m}={p}^{\ast ,m-1}-r\nabla \cdot {u}^{\ast ,m}\end{array}$ (9)

where r is an augmented Lagrangian penalty parameter used to impose the incompressibility 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\to +\infty$ , 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 $r\to 0$ 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 ${\rho }_{i}$ and ${\mu }_{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.

2.7. 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 $\left(n+1\right)\Delta t$ , except the non-linear inertial term that is linearized with a second order Adams-Bashforth scheme as follows

$u\cdot \nabla u\approx \left(2{u}^{n}-{u}^{n-1}\right)\cdot \nabla {u}^{n+1}$ (10)

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] .

3. 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 ( $Re={10}^{-3}$ ) past a cylinder of diameter $d=2m$ , with the undisturbed velocity being noted ${U}_{\infty }=1\text{\hspace{0.17em}}\text{m}/\text{s}$ , is solution of the Brinkman equation $-\nabla p+\mu \Delta {u}_{i}-\frac{\mu }{K}{u}_{i}=0$ . The reference solution is given in polar coordinate frame $\left(r,\theta \right)$ , centered on the particle, by:

(a)(b)(c)

Figure 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).

${u}^{*}\left({r}^{*},\theta \right)=\left\{\begin{array}{l}\frac{1}{{r}^{*}}\left(-\left(1+\frac{2{K}_{1}\left(\lambda \right)}{\lambda {K}_{0}\left(\lambda \right)}\right)\frac{1}{{r}^{*}}+{r}^{*}+\frac{2}{\lambda {K}_{0}\left(\lambda \right)}{K}_{1}\left(\lambda {r}^{*}\right)\right)\mathrm{cos}\theta \hfill \\ -\left(1+\left(1+\frac{2{K}_{1}\left(\lambda \right)}{\lambda {K}_{0}\left(\lambda \right)}\right)\frac{1}{{\left({r}^{*}\right)}^{2}}-\frac{2}{{K}_{0}\left(\lambda \right)}\left({K}_{0}\left(\lambda {r}^{*}\right)+\frac{{K}_{1}\left(\lambda {r}^{*}\right)}{\lambda {r}^{*}}\right)\right)\mathrm{sin}\theta \hfill \end{array}$ (11)

${p}^{*}\left({r}^{*},\theta \right)=\frac{2}{Re}{\lambda }^{2}\left(-\left(1+\frac{2{K}_{1}\left(\lambda \right)}{\lambda {K}_{0}\left(\lambda \right)}\right)\frac{1}{{r}^{*}}-{r}^{*}\right)\mathrm{cos}\theta$ (12)

where ${u}^{*}=\frac{u}{{U}_{\infty }}$ , ${p}^{*}=\frac{p}{\rho {U}_{\infty }^{2}}$ , ${r}^{*}=\frac{2r}{d}$ , $\rho =1\text{\hspace{0.17em}}\text{kg}\cdot {\text{m}}^{-3}$ is the fluid density, $\lambda =\frac{{d}^{2}}{4K}$ is the dimensionless permeability of the porous medium in Brinkman sens, K is the permeability of the inside and outside the porous cylinder, ${K}_{0}$ and ${K}_{1}$ are the modified Bessel functions of rank 0 and 1. For $K\to 0$ , the porous cylinder can be likened to an impermeable solid particle whereas outside the cylinder, $K\to +\infty$ to obtain a fluid behavior.

3.1. 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 Cartesian grid called Eulerian mesh, is represented by the number of gridcells across the diameter of the particle $\frac{d}{\Delta x}=20$ . The velocity and pressure exact solutions ((11), (12) respectively) for a Stokes flow past a cylinder were taken as initial condition, as illustrated in Figure 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 below:

• Viscous law: Arithmetic average law is chosen for this simulation.

• Numerical radius: it has been previously mentioned [13] 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 ${R}_{n}$ of the cylinder is given by:

${R}_{n}=\frac{d}{2}+e\frac{\Delta x}{16}$

e is a correction coefficient on ${R}_{n}$ . It is imposed to be $e=0$ , i.e. ${R}_{n}=\frac{d}{2}$ , so that ${R}_{n}$ is the physical radius of the cylinder for this simulation.

• Computation of ${C}_{\mu }$ : For the first simulation, it is interpolated from C, known in the pressure mesh, on the viscous mesh.

• The viscosity ratio $\frac{{\mu }_{s}}{{\mu }_{f}}$ between the viscosity imposed in the Eulerian cells inside the cylinder ${\mu }_{s}$ and the fluid viscosity ${\mu }_{f}$ is chosen such that $\frac{{\mu }_{s}}{{\mu }_{f}}=500$ for this simulation.

• The Lagrangian parameter is $r={10}^{5}$ .

Figure 4 shows the relative error in each point of the domain for the velocity

$\text{Error}=\left\{\begin{array}{l}\frac{|{u}_{\text{Simu}}-{u}_{\text{Analytic}}|}{|{u}_{\text{Analytic}}|}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}\text{ }{u}_{\text{Analytic}}\ne 0\hfill \\ |{u}_{\text{Simu}}|\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{ }\text{if}\text{\hspace{0.17em}}\text{ }{u}_{\text{Analytic}}=0\hfill \end{array}$ (13)

and the pressure between the simulation results and the analytical solution given by (11) and (12). This error is about 100% for the pressure in the fluid domain as illustrated in Figure 4(c) and more than 50% in the fluid region near the cylinder and about 10% in the rest of the fluid domain for both velocity components as illustrated in Figure 4(a) and Figure 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 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.

3.2. 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}_{\mu }$ which is interpolated at this state) and for different numerical radius as follows:

${R}_{n}=\frac{d}{2}+e\frac{\Delta x}{16},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}e\in \left[-16,16\right]$

All other parameters remain unchanged: $\frac{{\mu }_{s}}{{\mu }_{f}}=500$ and $r={10}^{5}$ .

Figure 5 shows the velocity L1 relative error in the whole domain

$\text{Error}=\frac{\sum |{u}_{\text{Simu}}-{u}_{\text{Analytic}}|}{\sum |{u}_{\text{Analytic}}|}$ (14)

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 ${R}_{n}=\frac{d}{2}-\Delta x$ whereas it is reached for a numerical radius ${R}_{n}=\frac{d}{2}+\frac{\Delta x}{8}$ for harmonic and discontinuous average laws. This minimum error is about 1% for both harmonic and discontinuous law whereas it is 0.5% larger for the arithmetic law with Rn

(a) (b) (c)

Figure 4. Relative error (%) in the whole domain (fluid and solid zones) of the first component of velocity u1 (a), the second component of velocity u2 (b) and the pressure p (c) for Stokes flow past cylinder.

Figure 5. Relative error (%) of the first component of velocity u1 (left) and the second component of velocity u2 (right) for Stokes flow past a cylinder simulated using different viscous laws: (-∙ -∙) discontinuous, (◊) harmonic and (▲) arithmetic average. Cμ is interpolated from C computed by projecting the cylinder on the pressure mesh.

being modified to 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}_{\mu }$ 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}_{\mu }$ 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 Rn and average viscous laws by considering the ${C}_{\mu }$ 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}_{\mu }$ 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 ${R}_{n}=\frac{d}{2}-\frac{\Delta x}{2}$ for arithmetic average law instead of ${R}_{n}=\frac{d}{2}-\Delta x$ when ${C}_{\mu }$ was interpolated, whereas the new ${C}_{\mu }$ 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 ${R}_{n}=\frac{d}{2}$ . Therefore, and for the rest of this work, the color function ${C}_{\mu }$ 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.

• Numerical radius: it is kept unchanged i.e. equal to the physical cylinder radius ${R}_{n}=\frac{d}{2}$ .

Figure 6. Relative error (%) on the first component of velocity u1 (left) and the second component of velocity u2 (right) for Stokes flow past a cylinder simulated using different viscous laws: (-∙ -∙) discontinuous, (◊) harmonic and (▲) arithmetic average. Cμ is computed by projecting directly the cylinder shape on the viscous mesh.

• The color function ${C}_{\mu }$ is computed on the viscous mesh instead of being interpolated from C.

• The viscosity ratio is the same $\frac{{\mu }_{s}}{{\mu }_{f}}=500$ .

• The Lagrangian parameter $r={10}^{5}$ 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.

3.3. Effect of the Viscosity Ratio and the Augmented Lagrangian Parameter r

The viscous penalty method consists in imposing large values of viscosity in the Eulerian cells belonging to the solid phase, compared to the fluid viscosity. This penalty method allows ensuring the solid behavior in the particles.

Therefore, the viscosity ratio $\frac{{\mu }_{s}}{{\mu }_{f}}$ is to be carefully considered to simulate gas-solid flows as best as possible with the viscous penalty method. For this motivation, numerous simulations of a uniform Stokes flow past an isolated fixed cylinder were carried out, with different values of $\frac{{\mu }_{s}}{{\mu }_{f}}$ , to study the viscosity ratio effect on the viscous penalty method accuracy. Figure 8 shows the velocity L1 relative error in the whole domain (14) for Stokes flow past a cylinder for a viscosity ratio between 100 and 1000. It can be observed that the error of the second component of velocity seems to be viscosity ratio independent from $\frac{{\mu }_{s}}{{\mu }_{f}}\ge 600$ and to stabilize for the first component of velocity when $\frac{{\mu }_{s}}{{\mu }_{f}}\ge 900$ . Therefore, $\frac{{\mu }_{s}}{{\mu }_{f}}=1000$ seems to be a reasonable choice in order to get a viscosity ratio independent solution. This

(a)(b)(c)

Figure 7. Relative error (%) of the first component of velocity u1 (a); of the second component of velocity u2 (b) and of the pressure p (c) for Stokes flow past cylinder in the whole domain.

viscosity ratio will be used in the rest of this work.

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\to +\infty$ . 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 103 and 109. One can observe that the solution is augmented Lagrangian parameter independent for $r\ge {10}^{5}$ . 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 $r={10}^{5}$ . This point is not addressed in the present work.

A new simulation is carried out with the set of all most efficient parameters, summarized below:

• Viscous law: harmonic average law.

• Numerical radius: physical cylinder radius ${R}_{n}=\frac{d}{2}$ .

Figure 8. Relative error (%) on the first component of velocity u1 (▲) and the second component of velocity u2 (Δ) for Stokes flow past a cylinder simulated with different viscosity ratio.

Figure 9. Relative error (%) on the first component of velocity u1 (▲) and the second component of velocity u2 (Δ) for Stokes flow past a cylinder simulated with different Lagrangian parameter values.

• The color function ${C}_{\mu }$ is computed on the viscous mesh.

• The viscosity ratio is $\frac{{\mu }_{s}}{{\mu }_{f}}=1000$ instead of $\frac{{\mu }_{s}}{{\mu }_{f}}=500$ .

• The Lagrangian parameter is kept as $r={10}^{5}$ .

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

(a)(b)(c)

Figure 10. Relative error (%) on the first component of velocity u1 (a); the second component of velocity u2 (b) and the pressure p (c) for Stokes flow past cylinder in the whole domain.

physical proprieties of the equivalent fluid in the mixed cells are neither fluid nor solid but an average of them, consequently the velocity and the pressure in these cells are less accurate.

3.4. 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 $\frac{d}{\Delta x}$ .

Figure 11 shows L1 relative error in the whole domain (14) for both component of velocity with respect to the Eulerian mesh resolution given by $\frac{d}{\Delta x}$ . The order of convergence computed from these error between the simulation results and (11), (12) is 1.67 based on a logarithmic data fit. It can be observed that some oscillations appear when refining the Eulerian mesh. A possible reason could be the effect of the particle interface position with respect to the Eulerian mesh. To assess this assumption, we have conducted different simulation by changing only the position of the cylinder inside the same Eulerian mesh: the cylinder center coordinates $\left({x}_{c},{y}_{c}\right)$ are:

${x}_{c}=i\frac{\Delta x}{10},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{where}\text{\hspace{0.17em}}\text{ }i\in \left[0,10\right]$

${y}_{c}=j\frac{\Delta y}{10},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{where}\text{\hspace{0.17em}}j\in \left[0,10\right]$

Figure 12 shows the effect of the position of the Lagrangian mesh with respect to the Eulerian mesh. It is observed that the way that the interface intersects the Eulerian mesh clearly affects the velocity results. In the convergence order study, the Eulerian mesh refinement changes the way the interface cuts the Eulerian cells, which explains the oscillations.

Figure 11. Relative error (%) on both velocity components for Stokes flow past a cylinder with respect to Eulerian mesh refinement.

Figure 12. Relative error (%) on the first component of velocity u1 (a) and the second component of velocity u2 (b) for Stokes flow past cylinder with respect to cylinder center position.

4. Uniform Flow past a Square Configuration of Cylinders

To validate the viscous penalty method outside the Stokes regime and with the new set of parameters prescribed in the previous section, an additional test is investigated: a uniform flow past a square configuration of cylinders. This configuration consists in putting a cylinder of a diameter d in a periodic square of length L. This configuration is equivalent to an infinite array of cylinders equidistant from each other in each direction. The Eulerian mesh refinement respects the condition given in [14] : $\Delta x=\frac{d}{5\sqrt{Re}}$ which ensure the boundary layer resolution if $Re>16$ and $\Delta x=\frac{d}{20}$ if $Re<16$ . The fluid is accelerated using a pressure drop ${F}_{m}=\frac{\Delta P}{L}$ as a source term in the momentum equation. The domain length L is fixed, given a solid volume fraction ${\alpha }_{d}$ , by:

$\frac{L}{d}=\frac{1}{2}\sqrt{\frac{\text{π}}{{\alpha }_{d}}}$

An illustration of a uniform flow past a square configuration of cylinders for different solid volume fraction ( ${\alpha }_{d}=0.2$ , ${\alpha }_{d}=0.4$ and ${\alpha }_{d}=0.6$ ) is given in Figure 13.

The aim of this section is to validate the superficial mean fluid velocity

$〈{u}_{f}〉=\left(1-{\alpha }_{d}\right)\frac{{\int }_{V}\left(1-C\right)u\text{d}V}{{\int }_{V}\left(1-C\right)\text{d}V}$

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. Numerous correlations have been proposed for predicting $\frac{\Delta P}{L}$ from the $〈{u}_{f}〉$ : Darcy [38] was the earlier pioneer in the subject by proposing in the Stokes limit the linear relation $\frac{\Delta P}{L}=\frac{\mu }{K}〈{u}_{f}〉$ . At higher Reynolds number, this relation is no longer linear due to inertial effects. Ergun [17] established a semi-empirical relation given by:

(a)(b)(c)

Figure 13. Streamlines and u1 component of velotcity field for a steady uniform flow along the x-axis of a square configuration of cylinders for different solid volume fractions: ${\alpha }_{d}=0.2$ (a); ${\alpha }_{d}=0.4$ (b); ${\alpha }_{d}=0.6$ (c).

$\frac{\Delta P}{L}=150\frac{{\alpha }_{d}^{2}}{{\left(1-{\alpha }_{d}\right)}^{3}}\frac{\mu 〈{u}_{f}〉}{{d}^{2}}+1.75\frac{{\alpha }_{d}}{{\left(1-{\alpha }_{d}\right)}^{3}}\frac{\rho {〈{u}_{f}〉}^{2}}{d}$ (15)

This relation is a generalization of the Forchheimer equation [39] .

Given the non-dimensional friction factor ${f}_{p}$ is defined as:

${f}_{p}=\frac{\Delta P}{L}\frac{d}{\rho {〈{u}_{f}〉}^{2}}$ (16)

and the Reynolds number Re is given by:

$Re=\frac{\rho 〈{u}_{f}〉d}{\mu {\alpha }_{d}}$

the Ergun Equation (15) can be written as:

${f}_{p}=\frac{{\alpha }_{d}}{{\left(1-{\alpha }_{d}\right)}^{3}}\left(\frac{150}{Re}+1.75\right)$ (17)

The validation consists in

• simulating a uniform flow past a square configuration of cylinders, for a given pressure drop $\frac{\Delta P}{L}$ and a solid volume fraction ${\alpha }_{d}$ .

• extracting from the velocity field the superficial mean fluid velocity $〈{u}_{f}〉$ .

• computing the friction factor ${f}_{p}$ using (16).

• comparing ${f}_{p}$ to the Ergun correlation [40] given by (17).

Figure 14 shows the good agreement of the friction factor deduced from the superficial mean fluid velocity $〈{u}_{f}〉$ using (16) and Ergun's correlation [40] . This validates the viscous penalty method at Higher Reynolds number with the best set of parameters found in the previous section: harmonic average for

viscous laws, ${R}_{n}=\frac{d}{2}$ , ${C}_{\mu }$ computed on the viscous mesh, $\frac{{\mu }_{s}}{{\mu }_{f}}=1000$ and $r={10}^{5}$ .

5. Conclusions and Suggestions

The Viscous Penalty Method ITPM [13] has been used to simulate two-dimen- sional fixed particulate flows. The first goal of the work was to set up the best

Figure 14. Friction factor for a uniform flow past a square configuration of cylinders.

numerical parameters in order to obtain lower errors as possible when the simulations are compared to the analytical solution for the Stokes flow around a cylinder. The viscosity ratio between the fluid and the penalty viscosity ${\mu }_{s}$ inside the particle, the augmented Lagrangian parameter r, the viscous law, the solid fraction estimate ${C}_{\mu }$ 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}_{\mu }$ is directly calculated by projecting the real shape of the particle on the viscous nodes, the numerical radius of the particle ${R}_{n}$ 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, $r={10}^{5}$ allows satisfying the incompressibility and solid constraints with lower errors as possible. Using larger values of ${\mu }_{s}$ and r does not improve the accuracy of ITPM, due to numerical errors coming from the rest of the numerical methods and solver efficiency. A convergence study was conducted with respect to mesh refinement. An order of 1.67 was obtained for all velocities inside the fluid.

A second problem was considered at larger particle Reynolds number: the uniform flow past a square arrangement of cylinders. With the best set of ITPM parameters, comparisons of simulations with reference correlations of Ergun allowed us to demonstrate that for various solid fractions ranging from 0.2 to 0.6, the simulations were in very good agreement with the expected values.

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] .

Acknowledgements

This work was granted access to the HPC resources of CINES under the allocation A0032b06115 made by GENCI (Grand Equipement National de Calcul Intensif) and to the resources of CALMIP supercomputing center under the allocation 2017-P1529.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Chadil, M. , Vincent, S. and Estivalèzes, J. (2019) Improvement of the Viscous Penalty Method for Particle-Resolved Simulations. Open Journal of Fluid Dynamics, 9, 168-192. doi: 10.4236/ojfd.2019.92012.