domain: spi phys

text:
Keywords: Point-force approximation, regularization, dispersed phase, Lagrangian particle tracking

Gaussian, an analytic expression from the literature is used to estimate the particle's self induced velocity. When combining source term regularization and removal of the particle's self induced velocity, good results are obtained for the terminal settling speed in the Stokes regime. Results obtained for horizontally separated particle pairs settling under gravity in the Stokes regime show equally good agreement with theoretical results. Because analytic expressions for the particle's self-induced velocity are no longer available at finite particle Reynolds numbers, correlations recently proposed in the literature are used to obtain correct settling velocities beyond the Stokes regime.

Regularization of the Lagrangian point force approximation for deterministic discrete particle simulations

Jean-François Poustis a , Jean-Mathieu Senoner a , Davide Zuzio a , Philippe Villedieu a a ONERA / DMPE, Université de Toulouse, F-31055 Toulouse -France

Abstract

The current article presents a regularization procedure of the Lagrangian point-force approach commonly used to account for the perturbation of a fluid phase by a dispersed particle phase. The regularization procedure is based on a nonlinear diffusion equation to naturally ensure parallel efficiency when the regularization length scale extends over several grid cells. The diffusion coefficient thus becomes a function of the particle source term gradient and expressions allowing to approximately adjust the regularization length scale according to the local particle to mesh size ratio are proposed, so that mesh refinement or polydisperse sprays may be handled. Elementary numerical test cases confirm the convergence of the present procedure under mesh refinement and its ability to locally adapt the regularization length scale. Furthermore, the chosen regularization length scale allows to match the leading order term of the perturbation flow field set by the particle beyond approximately two particle diameters in the Stokes regime. When applying the presented source term regularization procedure, the terminal velocity of a particle settling under gravity in the Stokes regime becomes relatively insensitive to mesh refinement. However, errors with respect to the theoretical settling velocity remain substantial and removal of the particle's self induced velocity appears necessary to recover the undisturbed fluid velocity at the particle location and correctly evaluate the drag force. As the current regularization procedure yields source terms that are close to c 2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/

Introduction

Dispersed two-phase flows are encountered in numerous natural phenomena and industrial applications such as cloud formation, volcanic eruptions, airframe icing, spray cooling, combustion chambers, fluidized beds, etc... Significant research efforts have helped improve our understanding of these flows via theoretical, experimental and numerical investigations. Accordingly, the literature on the subject is extremely vast and the interested reader is referred to the reviews of [START_REF] Drew | Mathematical modeling of two-phase flow[END_REF], [START_REF] Balachandar | Turbulent dispersed multiphase flow[END_REF], [START_REF] Jenny | Modeling of turbulent dilute spray combustion[END_REF] and [START_REF] Maxey | Simulation methods for particulate flows and concentrated suspensions[END_REF], among many others.

Due to the exponential increase of available computational resources in the last decades, a variety of numerical approaches have emerged to simulate dispersed two-phase flows. A major classification criterion pertains to the evaluation of the hydrodynamic forces acting on the dispersed inclusions, indifferently named particles in the following. These hydrodynamic forces may be directly evaluated by more advanced numerical methods or only approximately evaluated.

Arbitrary Lagrangian Eulerian (ALE) approaches deform the mesh of the flow domain so as to conform to the instantaneous particle configuration [START_REF] Hu | Direct numerical simulations of fluid-solid systems using the arbitrary Lagrangian-Eulerian technique[END_REF]. Direct numerical simulation (DNS) approaches for solid particles are based on the extension of Cartesian single-fluid solvers by the means of immersed boundary and penalty methods [START_REF] Peskin | Flow patterns around heart valves: A numerical method[END_REF]; [START_REF] Ritz | A numerical continuous model for the hydrodynamics of fluid particle systems[END_REF]; [START_REF] Mittal | Immersed boundary methods[END_REF]; [START_REF] Tenneti | Particle-resolved direct numerical simulation for gas-solid flow model development[END_REF]; [START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF]. These methods are fundamental in that they allow to increase the understanding of dispersed two-phase flow dynam-ics for configurations where analytic approaches currently fail. These approaches may therefore greatly help improve our understanding of two-phase flows, but they are often not applicable to realistic configurations because of the involved time and/ or length scale disparities. Therefore, simplified approaches remain of great interest.

An important class of methods deals with dispersed two-phase flows in the viscous Stokes regime. For such flow, the perturbation generated by an arbitrary shaped particle may be approximated via a Taylor expansion of the Green tensor, yielding the multipolar expansion [START_REF] Happel | Low Reynolds number hydrodynamics: with special applications to particulate media[END_REF]; [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF]. Due to the linear superposition principle applicable in the Stokes regime, the perturbation flow resulting from a set of particles may then be obtained via a sum of multipolar expansions, one for each particle. The problem then consists in finding the expansion coefficients that match the boundary conditions on the particle surfaces while taking into account the influence of a particle's flow perturbation on its neighbors. This may be done via methods that imply a resolution of the flow field at the particle surface or within the volume occupied by the particles, such as the Boundary Integral Method (BIM) [START_REF] Youngren | Stokes flow past a particle of arbitrary shape: a numerical method of solution[END_REF] or the Force Coupling Method (FCM) [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF]; [START_REF] Lomholt | Force-coupling method for particulate twophase flow: Stokes flow[END_REF]. Although these methods are very useful, they do not easily account for additional physical phenomena such as heating, melting, evaporation/ condensation, fragmentation, etc... at the particle scale. Additionally, dispersed two-phase flow approaches that may at least approximately reproduce convection dominated flows at the particle scale without explicitly resolving the latter are required for many practical applications [START_REF] Hervo | Large-eddy simulation of kerosene spray ignition in a simplified aeronautic combustor[END_REF]. Deterministic Lagrangian point-particle methods are an attractive choice to handle such cases. They allow for a straightforward implementation of physical models since each particle is tracked individually so that dispersed phase quantities are available at the particle scale. Moreover, the computational cost associated with the numerical resolution of the Lagrangian equations is moderate, although dependent on the number of particles to be simulated. One important limitation of deterministic Lagrangian point-particle methods resides in their coupling with the carrier fluid, which is generally ensured via the point-force (PF) approximation, also called Particle in Cell (PIC) method [START_REF] Harlow | The particle-in-cell method for numerical solution of problems in fluid dynamics[END_REF] [START_REF] Crowe | The particle-source-in cell (psi-cell) model for gas-droplet flows[END_REF]. In this approach, the coupling term is determined by its projection on the numerical grid, making the procedure non convergent under mesh refinement.

Different regularization techniques have been proposed in the literature to handle this issue. [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] proposed the Force Coupled Model (FCM), where the Dirac delta function is regularized as a Gaussian. By appropriately setting the regularization length scale, particle velocities may then be determined by a volumetric integral in the fictitious fluid volume occupied by the particle. Thus, the particle advancement step immediately follows from the fluid flow resolution step. However, this relation only holds in the viscous Stokes regimes and it implies that the fictitious flow field within the particle is at least approximately resolved. Thus, FCM may therefore be considered as a simplified immersed boundary/ penalization method. [START_REF] Gualtieri | Exact regularized point particle method for multiphase flows in the two-way coupling regime[END_REF] also used a Gaussian regularization in their Exact Regularized Point Particle (ERPP) method and applied it in the context of the deterministic Lagrangian point particle approach. They related their length scale to a characteristic diffusion time scale of vorticity. The authors split the perturbation flow field in a regular and singular component and showed that the contribution of the singular component could be neglected provided the particle Reynolds number was sufficiently small. [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF] derived and solved a filtered set of equations to explicitly account for the volume occupied by the Lagrangian particles so as to handle fluid flows in presence of a dense particle phase. Therefore, the regularization of the Dirac delta function is implicit within Capecelatro and Desjardins's modeling framework. As will be discussed shortly, the main concern then resides in the correct evaluation of unperturbed quantities at the particle location. This aspect was later addressed by [START_REF] Ireland | Improving particle drag predictions in Euler-Lagrange simulations with two-way coupling[END_REF], who referred to their extended version of Capecelatro and Desjardins's approach as volume-filtered Euler-Lagrange (VFEL).

The motivation of the present work lies in an extension of the regularization procedure proposed by [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF], albeit in the context of deterministic Lagrangian point particle approaches where the global volume fraction occupied by the particles is assumed negligible [START_REF] Saffman | On the settling speed of free and fixed suspensions[END_REF]. In order to locally adjust the regularization length scale and thus naturally handle varying particle to mesh size ratios, either due to polydispersity or spatial mesh size variations, the regularization is performed via the numerical resolution of a nonlinear diffusion equation. Using a regularization procedure based on the resolution of a transport equation allows to naturally ensure parallel efficiency, contrary to explicit regularization procedures which require great care in the parallel context [START_REF] Zwick | A scalable Euler-Lagrange approach for multiphase flow simulation on spectral elements[END_REF]. Moreover, the method appears general in the sense that it is equally suited for structured or unstructured numerical solvers. Finally, it appears that such equations have been successfully applied in image processing [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] and that conditions ensuring well-posedness and convergence to a unique solution are available, both in the continuous [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]; [START_REF] Walter | Differential and integral inequalities[END_REF] and the discrete case [START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF]; [START_REF] Drblíková | Convergence analysis of finite volume scheme for nonlinear tensor anisotropic diffusion in image processing[END_REF].

Because the point force approximation also involves the hydrodynamic force vector acting upon the particle, a regularization of the Dirac delta term alone may prove insufficient to ensure accurate coupling of the dispersed particle phase with the fluid phase. This is because the evaluation of the force terms acting on a given particle are based on the velocity field unperturbed by the current particle, but accounting for the perturbations of all other particles, in virtue of the Maxey-Riley Gatignol (MRG) equation [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF]; [START_REF] Gatignol | The Faxén formulas for a rigid particle in an unsteady non-uniform Stokes-flow[END_REF]. Errors associated with the evaluation of the unperturbed velocity seen by the particle may accumulate over time via the coupling term and yield substantial error levels. Current corrections are either based on an estimate of the flow perturbation for steady or unsteady viscous Stokes flows [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF]; [START_REF] Gualtieri | Exact regularized point particle method for multiphase flows in the two-way coupling regime[END_REF]; [START_REF] Ireland | Improving particle drag predictions in Euler-Lagrange simulations with two-way coupling[END_REF] or a truncated Taylor expansion to locally recover the unperturbed fluid velocity [START_REF] Horwitz | Accurate calculation of Stokes drag for pointparticle tracking in two-way coupled flows[END_REF]. In the present work, the analytic solution derived by [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] appears to yield a reasonably accurate estimate of the self-induced perturbation generated by the particle in the steady Stokes regime. At finite particle Reynolds number, a set of corrections recently proposed by [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] is used to estimate the velocity perturbation set by the particle according to the regularization length scale and particle Reynolds number. These correlations are compatible with the proposed methodology in a straightforward manner, allowing to retrieve correct settling speeds at finite particle Reynolds numbers.

The article is organized as follows. The first part briefly summarizes some fundamental results related to the point-force approximation and the FCM approach of Maxey and coworkers [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF]; [START_REF] Lomholt | Force-coupling method for particulate twophase flow: Stokes flow[END_REF], which will be used for validation and comparison purposes. Then, the regularization procedure based on the resolution of an unsteady nonlinear diffusion is presented. The third section discusses methods to approximate the unperturbed velocity field to improve the evaluation of the hydrodynamic force term in the point-force expression. Then, numerical results validating the proposed regularization procedure are presented. In particular, it is shown that the regularization procedure is mesh convergent.

The perturbation flow field resulting from the regularization procedure is also examined. Following these validations, the settling of an isolated particle and a pair of horizontally separated particles in presence of two-way coupling in the Stokes regime may be examined. Finally, the settling problem is also simulated at finite Reynolds numbers using the corrections proposed by [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF]. The main findings of the present work and possible strategies to achieve further improvements are summarized in the conclusion.

Point-force approximation

The present section briefly summarizes fundamental theoretical results on the point-force approximation. Although available in the classical textbooks of [START_REF] Happel | Low Reynolds number hydrodynamics: with special applications to particulate media[END_REF]; [START_REF] Kim | Microhydrodynamics: principles and selected applications[END_REF], they are rewritten for self-sufficiency of the present paper, as they will be used for validation and comparison in the forthcoming sections. In the Stokes regimes, the flow response to the application of a punctual force vector F located at position ξ:

-µ∇

2 u + ∇p = Fδ(x -ξ) (1)
is given by: u

(x -ξ) = G(x -ξ)F (2)
where G is the Oseen tensor:

G(x -ξ) = 1 8πµr I + (x -ξ) ⊗ (x -ξ) r 2 (3)
with r = |x -ξ| and x = ξ. Tensors will be written both in bold and underlined, whereas vectors will be written in bold only. The point force approximation correctly represents the leading order term of the perturbation induced by a translating rigid sphere submitted to the hydrodynamic force F, which decays proportionally to the distance from the sphere. Further analysis shows that this results from Lorentz's reciprocal theorem in the Stokes regime, which builds the starting point for the multipole expansion. Within this framework, the analytic solution for the flow field around a translating sphere in the Stokes regime may be described by a monopole and degenerate quadrupole. Because the degenerate quadrupole term is less straightforward to impose in a Navier-Stokes solver and because it inversely scales with the cube of the distance from the particle, the coupling term is simply approximated as the leading order monopole. This approximation is named point-force (PF) or Particle in Cell (PIC) method and will be abbreviated PF-PIC in the following [START_REF] Harlow | The particle-in-cell method for numerical solution of problems in fluid dynamics[END_REF][START_REF] Crowe | The particle-source-in cell (psi-cell) model for gas-droplet flows[END_REF]. The projection P < . > of the Dirac delta term on the numerical grid leads to an implicit regularization:

P < δ(x -x p ) >= 1 V C,i , if x p ∈ C i 0, otherwise. ( 4
)
where V C,i is a measure for the volume of i-th computational cell containing the center of gravity of the particle x p . Therefore, the discrete source term will become increasingly singular with mesh refinement. The source term may be projected over a set of neighbouring cells via conservative polynomial interpolations. Although such interpolations may efficiently remove the initial singularity, they are not mesh independent since the source term interpolation is generally based on a fixed set of neighbouring cells so that the regularization length scale is still dependent on a local characteristic mesh size. However, it is important to note that the stream function at large distances should remain convergent regardless of the mesh dependent regularization effect [START_REF] Batchelor | An introduction to fluid dynamics[END_REF].

In most applications, drag is the predominant hydrodynamic force acting upon the particle, and since it is proportional to the particle diameter d p in the Stokes regime, the singularity of the full coupling term depends on a ratio involving the particle diameter d p and the characteristic mesh length scale ∆x. This ratio will be denoted with the symbol ∆ in the following:

∆ = d p /∆x (5)
In order to remove the particle source term singularity, [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] proposed the Force Coupling Model (FCM). In this approach, the Dirac delta function is explicitly regularized as a Gaussian centered on each particle:

δ(x -x p ) → δ σ (x -x p ) = (2πσ 2 ) -n/2 exp -(x -x p ) 2 2σ 2 (6)
with n the number of spatial dimensions and σ the standard deviation. [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] then derived the fundamental solution of eq. 1Pointforce approximationequation.2.1 with the regularized source term of eq. 6Pointforce approximationequation.2.6. This solution is rewritten below because it will be used for comparison purposes in the following sections. The regularized Oseen tensor reads:

G δσ (x -ξ; σ) = 1 8πµ I r + (x -ξ) ⊗ (x -ξ) r 3 erf r √ 2σ + I r 3 - 3(x -ξ) ⊗ (x -ξ) r 5 σ 2 erf r √ 2σ - σ 2 2µ I - 3(x -ξ) ⊗ (x -ξ) r 2 σ r 2 δ σ (x -ξ) (7)
and the resulting flow field is given by eq. 2Point-force approximationequation.2.2 using relation ?? instead of 3Point-force approximationequation.2.3. Using σ as a free parameter, the authors set the regularization length scale to match the perturbation velocity at the particle's center with the terminal velocity of a particle settling under gravity in the Stokes regime. However, Lagrangian particle tracking approaches determine the particle velocity by an approximate force balance so that directly imposing the correct settling velocity is not viable for them, as will be further discussed in section 4Evaluation of the unperturbed fluid velocitysection.4. Moreover, if the characteristic stencil of the Gaussian function extends over several grid cells, great care is required to ensure the parallel efficiency of explicit regularization procedures [START_REF] Zwick | A scalable Euler-Lagrange approach for multiphase flow simulation on spectral elements[END_REF]. Therefore, the regularization procedure proposed in this paper is rather performed implicitly via the resolution of an unsteady nonlinear diffusion equation.

Regularization procedure via an unsteady nonlinear diffusion equation

In this section, an implicit source-term regularization procedure based on the numerical resolution of a diffusion equation is described. Such procedure should be inherently conservative and adapted to both structured as well as unstructured solvers. When standard numerical approaches are used for the resolution of such partial differential equation, only the first neighborhood of a computational cell is involved, ensuring parallel efficiency. First, the nonlinear diffusion equation will be presented. Then, the strategy to locally control the regularization length scale via this equation and the expression for the diffusion coefficient will be discussed.

Nonlinear diffusion equation

It was shown in the previous section that the singularity resulting from the PF-PIC source-term increases with the length scale ratio ∆ (see eq. 5Pointforce approximationequation.2.5). Thus, in Euler-Lagrange point-particle simulations with varying values of ∆, source terms associated with the largest values of ∆ may need to be regularized by spreading them over a finite and mesh independent length scale. This may be achieved with a linear diffusion equation and a constant diffusion coefficient [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF]. However, all source terms would then be regularized over the same characteristic length scale, either leading to excessive smoothing for source terms associated with smaller values of ∆ or insufficient smoothing of the source terms exhibiting the largest values of ∆. While insufficient smoothing may fail to efficiently remove the original singularity due to the point-force approximation and trigger numerical instabilities, excessive smoothing may result in poor prediction of the fluid flow perturbation induced by the particles. Therefore, the diffusion coefficient should be a function of space if polydisperse sprays and/or mesh refinement are to be handled accurately. In the present work, the diffusion coefficient is a function of the local source term gradient, leading to a non-linear formulation. Denoting the source term field to be regularized by φ, the problem to solve writes:

∂φ ∂τ -∇. D |∇φ| 2 ∇φ = 0, x ∈ Ω, τ ∈ [0; τ max ] (8) φ(0, x) = φ 0 (x), x ∈ Ω (9) ∂φ ∂n = 0, x ∈ ∂Ω ( 10
)
where Ω is the computational domain and τ max the final pseudo time of the regularization procedure. The initial condition φ 0 (x) is the projection of the Dirac delta functions related to the particle source terms on the numerical grid and writes:

φ 0 (x) = - p F p P < δ(x -x p ) > (11)
where x p and F p denote the position of the p-th particle and the force vector applied by the surrounding fluid to the particle p. The numerical projection operator P <> leads to a mesh dependent regularization, see eq. 4Point-force approximationequation.2.4. While eqs. 7Nonlinear diffusion equationequation.3.7-9Nonlinear diffusion equationequation.3.7 describe the regularization of a scalar field, the particle source term given by eq. 10Nonlinear diffusion equationequation.3.10 is a vector. Although a direct vector resolution with the definition of a proper tensor norm would have been possible, a component-wise resolution was chosen for simplicity. The problem defined by eqs. 7Nonlinear diffusion equationequation.3.7-9Nonlinear diffusion equationequation.3.7 is also known as Perona-Malik problem [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF]; [START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF]) and has been widely used and analyzed in the field of image processing.

The time variable τ and the diffusion coefficient D appearing in eq. 7Nonlinear diffusion equationequation.3.7 have no direct physical significance. For this reason, the maximum of the diffusion coefficient D is arbitrarily normalized to unity. However, both variables control the local regularization length scale. In the linear case, the diffusion of a Dirac delta function yields a Gaussian. From this analytic result, the diffusion length scale is found proportional to the square-root of the product between diffusion time and diffusion coefficient:

L D,s ≈ 10 σ (12)
with σ the standard deviation given by σ = 2Dτ max (13)

Therefore, the time required to diffuse a a source term over the regularization length scale L D,s writes:

τ max = L 2 D,s 200D (14)
These relations are a priori not valid in the nonlinear case. However, the nonlinear diffusion coefficient is designed in such way that it remains approximately constant in the vicinity of a given source term to be regularized. This aspect will be further discussed in section 5Numerical resultssection.5. For this reason, the above relations provide a sufficiently accurate estimate of the regularization length scale, even in the nonlinear case.

In order to make the regularization mesh independent, the regularization length scale L D,s is a function of the particle diameter. The choice for L D,s is motivated by two constraints. First, the singularity of the fluid velocity in the vicinity of the particle should be efficiently removed to ensure numerical stability. Second, the regularization should not alter the perturbation field set by the particle. As will be illustrated in section 5Numerical resultssection.5, an increase of the regularization length scale L D,s reduces the maximum fluid velocity levels in the vicinity of the particle. However, it also increases the distance at which the correct perturbation flow field is recovered, as may be deduced from eq. ??. Therefore, despite there being no precise quantitative criterion to define the regularization length scale L D,s , its value should result from a compromise between these constraints. The value proposed in the present work is:

L D,s = 6 d p (15)
Note that this value is comparable in magnitude to values proposed in the literature [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF]; [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF].

Expression of the diffusion coefficient

The current procedure aims at smoothing stiff source terms so that the diffusion coefficient should always be positive. As the magnitude of the source terms increases with the particle diameter to mesh size ratio, source term gradients may be used to estimate the local magnitude of this ratio and adjust the diffusion coefficient accordingly. Clearly, the diffusion coefficient should become largest for high source term gradients while vanishing for sufficiently low gradient values. Well-posedness of the initial value problem described by eqs. 7Nonlinear diffusion equationequation.3.7, 9Nonlinear diffusion equationequation.3.7 and 10Nonlinear diffusion equationequation.3.10 is guaranteed provided the diffusion flux function is positive and monotonously increasing. This condition is satisfied if the diffusion coefficient itself is described by a monotonously increasing function, see AppendixAProperties of the Perona-Malik equationappendix.A.

The constraints to be satisfied may be summarized as follows:

D(s) ≥ 0, ∀s ≥ 0 (16) ∂D(s) ∂s ≥ 0, ∀s ≥ 0 (17) lim s→0 D(s) = 0 (18) lim s→∞ D(s) = 1 (19)
It is easily verified that the following expression for the diffusion coefficient satisfies all these constraints:

D(|∇φ| 2 ) = 2 π arctan β |∇φ| 2 |∇φ| 2 max α (20)
with (α, β) ∈ R * + two adjustable parameters. The reasons motivating the choice for eq. 16Expression of the diffusion coefficientequation.3.16 and the values suggested for the parameters α and β are detailed in the following. The aim of the current procedure is to significantly smooth the largest source term gradients. The initial condition consists of a sum of Dirac delta functions, so that the initial numerical gradient values scale as ∆x -(n+1) , with ∆x the characteristic mesh size and n the number of spatial dimensions. Therefore, significant reductions of the largest gradient values, typically over an order of magnitude or more may be observed. In order to ease the adjustment of the local diffusion coefficient magnitudes, the magnitude of the diffusion coefficient is set via a normalization of the source term gradient by its global maximal value |∇φ| max at each iteration.

D(s 2 ) s/s max α = 0.5 α = 1 α = 2 α = 3 α = 4 α = 5
The arctan function is chosen because it allows for a relatively sharp separation between low normalized gradient values for which little to no diffusion should be applied and large values where the diffusion coefficient should reach an asymptotic value close to unity. An additional adjustment of the parameters α and β is necessary to control the magnitude of the diffusion coefficient according to the normalized source term gradient, as illustrated on figs. 1Evolution of the diffusion coefficient given by eq. 16Expression of the diffusion coefficientequation.3.16 for several values of the parameter α (at constant β).figure.caption.1 and 2Evolution of the diffusion coefficient given by eq. 16Expression of the diffusion coefficientequation.3.16 for several values of the parameter β with ∆ = 1 (at constant α).figure.caption.2. Fig. 1Evolution of the diffusion coefficient given by eq. 16Expression of the diffusion coefficientequation.3.16 for several values of the parameter α (at constant β).figure.caption.1 shows how a variation of the exponent α at constant β determines the abscissa at which asymptotic values close to unity (for large normalized gradient values) and close to zero (for low normalized gradient values) are obtained for the diffusion coefficient. Hence, choosing a low value for α implies an early saturation of the diffusion coefficient. As a consequence, significant smoothing would be applied even for relatively small normalized source term gradients and results would not significantly differ from those obtained with a linear diffusion equation. In combination with the parameter β, α also influences the average slope between the asymptotic regimes. Precisely controlling this slope is fundamental since it sets the magnitude of the diffusion coefficient, i.e. the length scale over which a source term is regularized. Furthermore, for the relations ??-?? to be approximately valid in the nonlinear context, the diffusion coefficient applied to a given particle should not significantly vary when several iterations are required to numerically solve the nonlinear diffusion equation. Therefore, the diffusion coefficient should not be too sensitive to local source term gradient variations either. Finally, the present work indicates that α is linked to the number of spatial dimensions. While preliminary one-dimensional tests yielded satisfactory results with α = 2, it was found that α = 4 was necessary to achieve equally good results in three dimensions. As previously mentioned, the parameter β controls the average slope of the diffusion coefficient in conjunction with the parameter α. Moreover, it also controls the abscissa for which the diffusion coefficient reaches its asymptotic value for large normalized gradient value. For this reason, this parameter is related to a global particle to mesh size ratio:

β = k(n) max Ω ∆ ( 21
)
where k is the actual constant to be adjusted.

Setting the particle to mesh size ratio value to ∆ = 1, fig. 2Evolution of the diffusion coefficient given by eq. 16Expression of the diffusion coefficientequation.3.16 for several values of the parameter β with ∆ = 1 (at constant α).figure.caption.2 shows diffusion coefficient profiles for different values of β obtained from variations of k. In order to obtain a pronounced upper asymptote, k must be sufficiently large and k = 150 was used in the present work.

Note that a change in the value of the global particle to mesh size ratio does not invalidate the choices made for these coefficients, as will be demonstrated in section 5Numerical resultssection.5. However, while the results obtained with the above formula and parameter values are satisfactory, they are probably not unique. This is because the influence of the parameters α and β on the diffusion coefficient profile is not clearly separable in the proposed formulation, in particular with respect to the slope between the two asymptotic regimes.

The present regularization procedure suffers from a lack of precision once the regularized source term becomes underresolved on the numerical grid, i.e. when it spans over a length scale comparable to the distance between the neighbouring cell centers so that mesh projected and regularized source terms become indistinguishable. Because of the normalization step by the largest source term gradient, the current regularization procedure will necessarily report a small source term contribution to the first set of neighbouring cells. Despite this contribution being small for low values of ∆, it results in a source term that is spread over an effective length scale larger than the target value L D,s . The regularized source term is underresolved for L D,s ≈ 2∆x, i.e. ∆ 0.33 for the present value L D,s = 6 d p . However, it appears that the effective regularization length scale tends to be overestimated starting from ∆ 0.5. For this reason, an additional constraint is applied to suppress diffusion when the effective regularization length scale L D,s is lower than 3∆x. This constraint is analogous to the one proposed by [START_REF] Capecelatro | An Euler-Lagrange strategy for simulating particle-laden flows[END_REF]:

2Dτ = L 2 D,s , if L 2 D,s > (3∆x) 2 0,
otherwise.

In the present work, eq. 7Nonlinear diffusion equationequation.3.7 was implemented in a structured Navier-Stokes solver using a standard centered spatial discretization of the diffusion operator. Time advancement is based on an explicit scheme. The resulting discrete equation is stable under the standard CFL condition:

δt ≤ ∆x 2 2n ( 22
)
where δt is the time step used for time integration and n the number of spatial dimensions. The properties of the fully discretized scheme are referenced in [START_REF] Weickert | Why the Perona-Malik filter works[END_REF]. The test cases indicate that the explicit CFL condition may become too restrictive for ∆ > 2, as about ten explicit iteration steps are required for ∆ = 2. Therefore, more efficient numerical resolution methods appear necessary to make the present methodology affordable for realistic applications in presence of large particle to mesh size ratios. Since the focus of the present work lies in methodological aspects, such performance considerations are postponed to future work.

Evaluation of the unperturbed fluid velocity

An additional difficulty of the regularization procedure lies in the accurate evaluation of the drag force, which sets the amplitude of the source term to be regularized, see eq. 10Nonlinear diffusion equationequation.3.10.

In virtue of the Maxey-Riley-Gatignol (MRG) equation [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF]; [START_REF] Gatignol | The Faxén formulas for a rigid particle in an unsteady non-uniform Stokes-flow[END_REF], the evaluation of the drag force of a given particle requires knowledge of the relative velocity between that particle u p and the undisturbed fluid at the particle location u f @p . The latter is defined as the fluid velocity that would be seen in absence of the considered particle at the same location but accounting for the fluid flow perturbations induced by all other particles present in the flow. Clearly, the undisturbed fluid velocity u f @p is not directly available in numerical simulations and the explicit evaluation of each particle's flow perturbation, also referred to as self-induced velocity in the literature, would be prohibitively costly in realistic applications [START_REF] Boivin | Direct numerical simulation of turbulence modulation by particles in isotropic turbulence[END_REF].

For this reason, u f @p is often approximated by the perturbed fluid velocity at the particle location u f @p . As may be seen from eqs. 2Point-force approximationequation.2.2 and 3Point-force approximationequation.2.3, the leading error term of this approximation scales as O(d p /∆x = ∆) in the Stokes regime. This is consistent with previous findings: the fluid flow singularity in the vicinity of the particle is implicitly regularized by the projection of the particle source term on the numerical grid only as long as the ratio between particle diameter and local mesh size is small compared to unity, motivating the previously described source term regularization procedure. The latter effectively removes the flow field singularity at the particle position so that the particle's self-induced velocity may be directly evaluated in a closed form under simplifying assumptions. Assuming a Gaussian source term regularization of variance σ, [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] derived an analytic expression for the particle's self-induced velocity u * ,St f @p in the steady Stokes regime by evaluating eq. ?? at x = ξ:

u * ,St f @p = u f @p -u f @p = 2 π 1/2 d p 2σ (u p -u f@p ) (23)
In the present case, the regularization length scale is set to σ = 0.6 d p in virtue of eq. ??. Therefore, despite efficiently removing the initial singularity and in turn mesh dependency of the point-force approximation, it appears from eq. ?? that source term regularization must be combined with a correction on the particle's self-induced velocity to warrant accurate results. Note that a simpler derivation of eq. ?? was recently proposed by [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF]. The latter derivation also accounts for the fluid volume occupied by the particles, yielding a more general expression. However, it is important to note that volume fraction corrections are small compared to corrections on the particle's self-induced velocity [START_REF] Ireland | Improving particle drag predictions in Euler-Lagrange simulations with two-way coupling[END_REF].

As most practical applications involve moderate but finite particle Reynolds numbers [START_REF] Hervo | Large-eddy simulation of kerosene spray ignition in a simplified aeronautic combustor[END_REF]; [START_REF] Norde | Eulerian and Lagrangian ice-crystal trajectory simulations in a generic turbofan compressor[END_REF], the range of validity of eq. ?? should be extended beyond the Stokes regime. Assuming a Gaussian source term regularization, [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] derived a set of correlations from numerical simulations to estimate u * f @p for particle Reynolds numbers up to 200 as a function of the regularization length scale. The previous expression of the particle's self-induced velocity in the steady Stokes regime is then multiplied by a corrective function Ψ:

u * f @p = u * ,St f @p Ψ(Re p,σ , F ) ( 24
)
where Ψ is decomposed in a corrective function derived from Oseen's linearization Ψ os and an empirical function χ:

Ψ(Re p,σ , F ) = Ψ os (Re p,σ ).χ(Re p,σ , F ) (25)
The Gaussian Reynolds number Re p,σ explicitly depends on the particle Reynolds number Re p and on the regularization length scale through the Gaussian variance σ:

Re p,σ = Re p σ d p (26)
whereas F represents an expression of the drag force non-dimensionalized with the inertial scaling:

F = 3π Φ(Re p ) Re p d p σ 2 (27)
Φ(Re p ) denotes the corrective function for the particle drag at finite particle Reynolds numbers. In the present case, the expression proposed by [START_REF] Schiller | A drag coefficient correlation[END_REF] is used:

Φ(Re p ) = (1 + 0.15Re 0.687 p ) (28)
For more details on the expressions of the functions Ψ os and χ, the reader is referred to [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF].

For the present methodology, the evaluation of the equivalent Gaussian variance σ appearing in eqs. ??-?? depends on the effective regularization length scale. If the regularized source term is sufficiently resolved on the numerical grid, σ may be directly deduced from eqs. ?? and 11Nonlinear diffusion equationequation.3.11. Once the regularized source term becomes underresolved, it is assumed that the equivalent Gaussian spans over the first neighbouring cell centers, i.e. L D,s = 2 ∆x. In the present case, source term regularization is suppressed starting from ∆ = 0.5 so that the first method is used for values of ∆ > 0.5 and the second for ∆ ≤ 0.5.

A different strategy to locally remove the particle's self-induced velocity at finite Reynolds number could consist in filtering the perturbed velocity field. As demonstrated earlier, the present non-linear diffusion equation allows to locally set a desired diffusion/ filtering length scale and could possibly serve this purpose. However, an approach based on filtering may only be successful if the resulting regularization length scale is similar in magnitude to the particle diameter. This is because a significant scale separation between the smallest fluid flow structures and the flow disturbances generated by the particles cannot be invoked in the general case. In practice, it turns out that the regularization length scale required to accurately predict the terminal settling velocity of a particle submitted to gravity is only weakly dependent on the particle Reynolds number and of the order of ten particle diameters, which seems too large. In a similar spirit, a method based on the subtraction of perturbed velocity fields filtered at two different scales, typically three and four particle diameters, was also tested in the present work. While this method yields excellent results for the settling of a single particle under gravity regardless of the particle Reynolds number, it fails to predict the settling of a pair of horizontally separated particles. This is because the perturbation flow field set by the neighbouring particle becomes relatively slowly decaying as the particles become separated by more than a few particle diameters. Hence, subtraction of two perturbed velocity fields filtered at similar length scales leads to a significant removal of the perturbation flow field set by the neighbouring particle. In turn, significant errors on the final settling speed are obtained. Therefore, our present numerical experiments tend to indicate that procedures based on the filtering of the perturbed flow field are not suited for the removal of the particle's self-induced velocity. Thus, the analytic expression originally proposed by [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] is combined with the corrections of [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] at finite particle Reynolds number to remove the particle's self-induced velocity / retrieve the undisturbed fluid velocity in the present work.

Numerical results

In the present section, test cases validating the proposed methodology are presented. First, elementary test cases related to the validation of the nonlinear diffusion equation for source term regularization are discussed. Then, the coupling with the fluid phase via the regularized source terms is assessed. Numerical results are compared to analytic solutions in the Stokes regime. Because the adjustment of the parameters defining the diffusion coefficient is more tedious as the number of dimensions increases, all results presented in the current section are obtained from three-dimensional simulations.

The present code solves the unsteady incompressible Navier-Stokes equations on structured Cartesian grids with staggered variable storage. It is is based on a Chorin projection method. First, the momentum equation is solved without pressure term to obtain a first prediction of the velocity field. The pressure field is then computed using a Poisson equation based on the predicted velocity field. Finally, a divergence-free velocity field is obtained by applying the pressure field correction to the predicted velocity field. Time integration is performed via a second order Runge-Kutta scheme. More details about the resolution of the numerical methods used for the resolution of the fluid phase may be found in [START_REF] Zuzio | An improved multiscale Eulerian-Lagrangian method for simulation of atomization process[END_REF].

Regarding the Lagrangian particle solver, fluid properties are evaluated at the particle position via a simple trilinear interpolation scheme. The source term is fully reported to the velocity cell containing the particle's center of gravity in virtue of eq. 10Nonlinear diffusion equationequation.3.10. Finally, the time integration scheme for the particles is identical to that of the fluid solver.

Numerical resolution of the non-linear diffusion equation

The first numerical test case illustrates the source term regularization resulting from the numerical resolution of the non-linear diffusion equation with a single point-force term as initial condition, i.e. a Dirac delta function projected on the numerical grid, see eq. 4Point-force approximationequation.2.4. The computational domain consists of a cubic box of size L = 10 d p . The numerical resolution is performed for a broad range of particle to mesh size ratios, namely ∆ = (1; 2; 4; 8; 16). Note that the largest values are not realistic and that they are only used to demonstrate numerical convergence under mesh refinement. Choosing the regularization length scale of L D,s = 6 d p , the diffusion time τ max is set according to eq. ??. The regularized source terms are displayed on fig. 3Solution of the non-linear diffusion equation under mesh refinement for particle to mesh size ratios ∆ ranging between 0.5 and 16. All source terms are normalized by Φ(0), the value of the Gaussian at r = 0.figure.caption.3, along with the Gaussian curve spread over the same characteristic length scale, i.e. σ = 0.6 d p .

-

3 -2 -1 0 1 2 3 r/d p 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Φ ∆ (r)/Φ(0) Gaussian Distribution ∆ =1 ∆ =2 ∆ =4 ∆ =16
Figure 3: Solution of the non-linear diffusion equation under mesh refinement for particle to mesh size ratios ∆ ranging between 0.5 and 16. All source terms are normalized by Φ(0), the value of the Gaussian at r = 0.

It appears that the regularization length scale is properly imposed, regardless of the value of ∆. Furthermore, the curves obtained for the largest particle to mesh size ratios indicate that the present method yields regularized source term distributions that closely match the Gaussian reference curve. This result was unexpected since the Gaussian distribution is not a solution of the present non-linear diffusion equation. However, the following test case will show that the diffusion coefficient prescribed via eq. 16Expression of the diffusion coefficientequation.3.16 yields approximately constant diffusion values in the vicinity of a given source term. In order to ensure that the current procedure is mesh convergent, L 2 and L ∞ errors of the different solutions are computed with respect to the numerical reference solution obtained for the largest refinement ∆ = 16:

L 2,∆ = i |φ ∆ (x i ) -φ 16 (x i )| 2 ∆V i 1/2 (29) L ∞,∆ = max (|φ ∆ (x i ) -φ 16 (x i ))|) ( 30
)
with ∆V i the volume of the i-th computational cell.

The results are shown in fig. 4L 2 and L ∞ errors of the regularized source terms under mesh refinement with respect to the solution obtained on the finest mesh, i. e. ∆ = 16.figure.caption.4. It appears that the error converges with mesh refinement for both norms and that the present numerical scheme is second-order accurate. Although the present implementation was only performed in a Cartesian framework, it is expected that second-order accuracy may be straightforwardly achieved in an unstructured finite-volume solver [START_REF] Coudière | Convergence rate of a finite volume scheme for a two dimensional convection-diffusion problem[END_REF]; [START_REF] Drblíková | Convergence analysis of finite volume scheme for nonlinear tensor anisotropic diffusion in image processing[END_REF].

In order to illustrate that the present methodology may naturally handle particle to mesh size ratio variations (either due to local mesh refinement or to polydispersion), the former test case is now performed for three point-force source terms spread on the same mesh and corresponding to three distinct particle to mesh size ratios, namely ∆ = 1, 2, 4 . Except for the largest value ∆ = 4, this range is assumed representative of the largest particle to mesh size ratios that may be encountered for dispersed two-phase flows in realistic applications. The aim is to evaluate the regularization procedure for particle to mesh size ratios where it will be active, i.e. ∆ > 0.5 and in presence of significant variations of this parameter. The distance between each particle is large enough to avoid interactions between the different source terms. The computational domain consists of a cubic box of size L = 30 d p . The regularized source term distributions are then compared to the respective Gaussian reference curves.

The result is displayed on fig. 5Regularized source term field and Gaussian reference curves with the same regularization length scale for multiple source terms. Both curves are normalized by the maximum of the Gaussian reference curves.figure.caption.5 and shows that the present methodology allows to locally adapt the regularization length scale. The adjustment of the local regularization length scale via the variation of the diffusion coefficient is illustrated in fig. 6Diffusion coefficient field for the multiple source term test case at five different non-dimensional instants τ * = τ /τ max .figure.caption.6, which displays the diffusion coefficient profiles at different instants for the multiple source term regularization problem. First, it appears that the diffusion coefficient magnitudes do not vary significantly over time in the vicinity of each source term, except in the computational cell where the source term is applied and its direct neighbors. Here, the source term gradients necessarily decrease as the source term gradient reaches a minimum. Due to the chosen normalization procedure by the largest gradient at each iteration, the source term corresponding to the largest value of ∆ displays a diffusion coefficient close to unity over the entire diffusion time τ max . The diffusion coefficients are adjusted for the smaller source term gradients, reducing the regularization length scale with respect to the source term with the largest value of ∆. Therefore, although approximate, the proposed expression for the non-linear diffusion coefficient allows to adjust the local regularization length scale in combination with an appropriate diffusion time scale. As the resulting diffusion coefficients are contained within a limited range in the vicinity of a given source term and since locally they do not significantly vary over time either, the resulting regularized source terms display shapes that are close to Gaussian curves. Note that figs. 3Solution of the non-linear diffusion equation under mesh refinement for particle to mesh size ratios ∆ ranging between 0.5 and 16. All source terms are normalized by Φ(0), the value of the Gaussian at r = 0.figure.caption.3 and 5Regularized source term field and Gaussian reference curves with the same regularization length scale for multiple source terms. Both curves are normalized by the maximum of the Gaussian reference curves.figure.caption.5 also indicate that for the present regularization length scale L D,S = 6 d p , explicit regularizations only involve the first set of neighbouring cells up to approximately ∆ = 1, thus involving no communication overhead for parallel simulations up to this value. However, several rows of neighbouring cells may be involved for ∆ values of order unity when considering anisotropic structured grids and unstructured grids. Moreover, wider regularization stencils, i.e. L D,s > 6 d p , may also be necessary for stability reasons if the hydrodynamic forces acting on the particle are large. Thus, the present procedure may typically be recommended for two-phase flow simulations involving important local refinements such as reactive two-phase LES's of aeronautic combustion chambers [START_REF] Hervo | Large-eddy simulation of kerosene spray ignition in a simplified aeronautic combustor[END_REF] and dispersed two-phase flows in resolved turbulent boundary layers.

Coupling with the fluid phase

Now that the approximate regularization procedure for the point-force source term is validated, the flow field perturbations resulting from these regularized source terms will be examined. First, the flow field response to the presence of a single point force written as a Dirac delta function in a fluid initially at rest is examined. A source term corresponding to the Stokes drag exerted upon a particle of diameter d p and moving with velocity u p,s = u p,s y in a fluid at rest, i.e. F = 6πr p µ f u p,s , is imposed at the center of a cubic domain of size L = 32 d p . For these simulations, the velocity boundary conditions are set to match Stokes's far field solution for the flow around a sphere of the same diameter and whose center is located at the same position [START_REF] Chatelin | Méthodes numériques pour l'écoulement de stokes 3d: fluides à viscosité variable en géométrie complexe mobile; application aux fluides biologiques[END_REF] The Cartesian mesh is chosen such that the particle to mesh size ratio ∆ is equal to 2. Results are displayed in a coordinate system (x, y, z) whose origin x 0 coincides with the source term position and the fluid velocity vector is denoted u = (u, v, w).

The simulation is continued until a steady state is reached for the fluid phase. The analytic solution of this problem corresponds to the well-known Oseen tensor when the source term is expressed as the product of the force vector with the Dirac delta function, see eq. 3Point-force approximationequation.2.3. This solution correctly represents the leading order term for the flow field around a rigid sphere in viscous uniform flow. In addition, [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] It appears that the analytic solution of [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] and the perturbation resulting from the present regularization procedure for the same length scale L D,s = 6 d p are almost superimposed for r y /d p > 1 in the direction of the force and r x /d p > 2 in the direction perpendicular to the force. On the contrary, differences become more and more significant close to the origin. This result is consistent with the regularization curves obtained in section 5.1Numerical resolution of the non-linear diffusion equationsubsection.5.1, which deviate most significantly from the Gaussian reference curves in the vicinity of the origin because the diffusion coefficient is necessarily lower in the cell where the source term is imposed. Therefore, errors are to be expected on the estimate of the particle's self-induced velocity: the analytic value obtained from the evaluation of eq. ?? at the origin deviates from the value obtained through the implicit numerical regularization by approximately 8%. The associated error levels vary only slightly with resolution as long as the regularized source term is discretized on at least 5 grid cells. Once the regularized source term becomes underresolved on the numerical grid, i.e. indistinguishable from the mesh projected source term, the regularization procedure is suppressed. In that case, the perturbation velocity at the particle position is estimated by assuming that the source term corresponds to a Gaussian spread over L D,s = 2 ∆x. This procedure yields comparable error levels, i.e. 5% for ∆ = 0.5.

In addition, the magnitude of the perturbation velocity at the origin decreases for increasing values of L D,s as the singularity originating from the point-force term is more efficiently damped. On the contrary, the distance at which the perturbation flow fields matches Stokes's solution increases for larger values of L D,s , in particular for the perturbation field in the direction of the force, see fig. 7Fluid For this reason, an increase of the regularization length scale is expected to deteriorate the accuracy of the predicted fluid perturbation flow field while it should alleviate numerical stability issues by reducing the peak of the perturbation velocity at the particle position. Therefore, the best choice for the regularization length scale seems to involve a compromise between stability and accuracy. Furthermore, it depends on the range of values of ∆ that are expected to be encountered in a given application. Therefore, while general recommendations on the best value for L D,s seem difficult, L D,s = 6 d p appears as a reasonable choice to validate the present methodology. In particular, this length scale yields a perturbation velocity that lies below the relative velocity between particle and fluid for the steady point-force problem, in loose analogy with a maximum principle.

In order to illustrate the aforementioned stability issues, results for the steady point-force problem without regularization are provided in AppendixB-Necessity of the regularization procedureappendix.B. In addition, the settling under gravity of an isolated particle in presence of two-way coupling described in section 5.3.1Stokes regimesubsubsection.5.3.1 is performed without regularization but with a correction on the perturbation velocity only in Ap-pendixBNecessity of the regularization procedureappendix.B. For ∆ = 1, pronounced oscillations on the settling velocity become visible in fig. B.15Settling velocity histories for ∆ = 1. Comparison of one-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations as well as two-way coupled simulations with regularized source terms without ("2W + RegST") and with removal ("2W + RegST + AC") of the particle's self-induced velocity. The final simulation ("2W + AC") concerns a PF-PIC two-way coupled simulation without source term regularization but with the removal of the particle's self-induced velocity.figure.caption.19. Clearly the occurrence of these oscillations depends on numerical details of the source term projection scheme as well as details of the considered Navier-Stokes solver so that general recommendations seem difficult to make. Nevertheless, these results tend to indicate that for large particle to mesh size ratios, source term regularization may become mandatory to ensure numerical stability.

Particle settling under gravity 5.3.1. Stokes regime

The current section examines the elementary test case of a particle settling under gravity. Despite its simplicity, this test case clearly demonstrates the convergence issues inherent to the Lagrangian point-force approximation in presence of two way-coupling [START_REF] Gualtieri | Exact regularized point particle method for multiphase flows in the two-way coupling regime[END_REF]; [START_REF] Zuzio | An improved multiscale Eulerian-Lagrangian method for simulation of atomization process[END_REF]; [START_REF] Ireland | Improving particle drag predictions in Euler-Lagrange simulations with two-way coupling[END_REF]; [START_REF] Horwitz | Correction scheme for point-particle models applied to a nonlinear drag law in simulations of particle-fluid interaction[END_REF]. All forces except drag and gravity are neglected so that particle dynamics are governed by the following equations:

dx p dt = u p (31)
du p dt = u f @p -u p τ p + 1 - ρ f ρ p g (32)
with u f @p is the fluid velocity at the particle position unperturbed by the latter, but accounting for the perturbations of other particles. The settling velocity of the isolated particle is simply obtained as:

u p,s = τ p 1 - ρ f ρ p g (33)
Assuming the particle Reynolds number Re p to be much smaller than unity, the particle relaxation time is analytically known as:

τ p = ρ p d 2 p 18µ f (34)
so that eq. 20Stokes regimeequation.5.20 may be evaluated explicitly. The perturbation of the fluid phase induced by the particle is simulated using the regularized source term ('2W + RegST') and the standard PF-PIC approximation ('2W') for comparison purposes. In both cases, the unperturbed fluid velocity at the particle location u f @p is approximated by its perturbed counterpart u f @p . A simulation with one-way coupling ('1W') is used as reference. Three different particle to mesh size ratios ∆ = (0.5, 1, 2) are simulated. The simulation domain is a cube of size L = 64 d p . In the transverse directions, slip boundary conditions are enforced while a periodic boundary condition is imposed in the direction of the gravitational force.

It was verified that doubling the size of the box did not affect results. All relevant physical parameters are summarized in table 1Summary of numerical and physical parameters for the particle settling under gravity in the Stokes regime. The particle Reynolds number is based on the terminal settling velocity.table.caption.9. Following [START_REF] Horwitz | Correction scheme for point-particle models applied to a nonlinear drag law in simulations of particle-fluid interaction[END_REF], a Stokes number St ∆ based on the ratio τ p /τ visc , with τ visc = ∆x 2 /ν f is used. 14,55,222) Table 1: Summary of numerical and physical parameters for the particle settling under gravity in the Stokes regime. The particle Reynolds number is based on the terminal settling velocity.

∆ (0.5, 1, 2) Re p 0.04 ρ p /ρ f 1000 St ∆ (
The evolution of particle velocity over time is shown in fig. 9Comparison of settling velocity histories in the Stokes regime for one-way ("1W"), standard PF-PIC two-way ("2W") coupled simulations and two-way coupled simulations with source term regularization ("2W + RegST").figure.caption.10. Both variables are made non-dimensional by dividing them respectively by the particle relaxation time τ p and the final settling velocity u p,s . It appears that the particle's settling velocity is much less sensitive to mesh refinement when the source term regularization is applied, although variations of about 10% are visible for the different values of ∆. Since the regularization procedure itself was shown to be mesh convergent, the remaining spread observed for the settling velocities is probably due to resolution effects. The residual error on the settling velocity remains significant with respect to the reference solution, although significantly lower than the solution obtained with the PF-PIC approach for the cases ∆ = 1 and ∆ = 2. It appears that the error levels on the terminal settling velocity are comparable for the case ∆ = 0.5 without regularization and ∆ = 1 and ∆ = 2 with regularization, which confirms that the regularization procedure should have little effect for ∆ ≤ 0.5. 1.8 1.9

Figure 9: Comparison of settling velocity histories in the Stokes regime for one-way ("1W"), standard PF-PIC two-way ("2W") coupled simulations and two-way coupled simulations with source term regularization ("2W + RegST").

These results highlight the necessity to remove the particle's self-induced velocity for a correct evaluation of the drag force. Therefore, the settling problem is simulated with the analytic correction given by eq. ?? for ∆ = (0.5, 1, 2). Despite some discrepancies, the results displayed in fig. 10Settling velocity histories in the Stokes regime for varying ∆. Comparison of one-way ("1W") with two-way coupled simulations with source term regularization and with analytic removal ("2W + RegST + AC") of the particle's selfinduced velocity.figure.caption.11 indicate that the correction is effective as error levels are significantly reduced compared to the results obtained with source term regularization only. The relative errors for the terminal settling velocities with respect to the reference solution amount to +0.8% for ∆ = 0.5, +3% for ∆ = 1 and -2.5% for ∆ = 2. These orders of magnitude are close to those reported by [START_REF] Zuzio | An improved multiscale Eulerian-Lagrangian method for simulation of atomization process[END_REF]; [START_REF] Ireland | Improving particle drag predictions in Euler-Lagrange simulations with two-way coupling[END_REF]; [START_REF] Horwitz | Correction scheme for point-particle models applied to a nonlinear drag law in simulations of particle-fluid interaction[END_REF]. Comparison of one-way ("1W") with two-way coupled simulations with source term regularization and with analytic removal ("2W + RegST + AC") of the particle's self-induced velocity.

Finite particle Reynolds number

The present section reproduces the same test case as described in the previous section, but at finite particle Reynolds numbers. The particle Reynolds numbers based on the terminal settling velocity are respectively Re p = 10 and Re p = 50. These Reynolds numbers are obtained by decreasing the particle to fluid density ratio and adjusting the fluid viscosity, see table 2Summary of numerical and physical parameters for the particle settling under gravity at finite Reynolds number. The particle Reynolds numbers are based on their respective terminal settling velocity.table.caption.12. The chosen particle to mesh size ratios are identical to those used in the Stokes regime, i.e. ∆ = (0.5, 1, 2). ∆ 0.5 1 2 Re p (10, 50) ρ p /ρ f 85.33 µ f

• 10 -3 , 8.9 • 10 -4 ) St ∆ (0.69, 0.37) (2.70, 1.47) (11, 5.9)

Table 2: Summary of numerical and physical parameters for the particle settling under gravity at finite Reynolds number. The particle Reynolds numbers are based on their respective terminal settling velocity.

For finite Reynolds number, the particle relaxation time is evaluated with the correlation of [START_REF] Schiller | A drag coefficient correlation[END_REF]

τ p = ρ p d 2 p 18µ f (1 + 0.15Re 0.687 p ) (35)
While still involving a single unknown, eq. 20Stokes regimeequation.5.20 may no longer be solved explicitly when using eq. 22Finite particle Reynolds numberequation.5.22. As discussed in the previous section, the analytic correction given eq. ?? is no longer valid and the corrections for finite particle Reynolds numbers proposed by [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF], summarized by eqs. ??-??, must be used instead.

As may be seen from figs. 11Settling velocity histories at Re p = 10 for ∆ = (0.5, 1, 2) with source term regularization ("RegST") and removal of the particle's self-induced velocity based on the corrections of [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] ("BC"). One-way ("1W") and standard PF-PIC two-way 2019) ("BC"). One-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations and two-way coupled simulations with regularized source terms only ("2W + RegST") for ∆ = 1 are also shown for comparison purposes.

("2W") coupled simulations and two-way coupled simulations with regularized source terms only ("2W + RegST") for ∆ = 1 are also shown for comparison purposes.figure.caption.13 and 12Settling velocity histories at Re p = 50 for ∆ = (0.5, 1, 2) with source term regularization ("RegST") and removal of the particle's self-induced velocity based on the corrections of Balachandar et al. ( 2019) ("BC"). One-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations and two-way coupled simulations with regularized source terms only ("2W + RegST") for ∆ = 1 are also shown for comparison purposes.figure.caption.14, the corrections proposed by [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] seem fully compatible with the present methodology, despite a slight increase of the maximum error on the terminal settling velocity to 6%

for ∆ = 0.5 at Re p = 10. In particular, errors with respect to the theoretical solution amount to only 0.6 % for Re p = 10 and 1.11% for Re p = 50 for ∆ = 1. Moreover, [START_REF] Horwitz | Correction scheme for point-particle models applied to a nonlinear drag law in simulations of particle-fluid interaction[END_REF] also simulated the settling problem under gravity at Re p = 10 with nondimensional parameter values similar to those of table 2Summary of numerical and physical parameters for the particle settling under gravity at finite Reynolds number. The particle Reynolds numbers are based on their respective terminal settling velocity.table.caption.12 and their error levels are comparable to the present ones. Therefore, while [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] raised some concerns about the compatibility of their corrections obtained from numerical simulations with a highly accurate spectral element code to solvers with less accurate and different numerical algorithms as used in the present tests, these differences do not seem to significantly affect the accuracy on the terminal settling velocity. This effect could possibly be explained by the significant decrease of the relative error levels on terminal settling velocities for the numerical simulations with source term regularization only as the problem becomes more dominated by convective effects.

Particle Pair

The present section considers the settling under gravity of a pair of identical particles in the Stokes regime. The previous section showed that the present approach is able to predict the terminal velocity of an isolated particle settling under gravity in the presence of two-way coupling, for particle Reynolds numbers up to 50. When considering a pair of particles, each particle alters the motion of the other. Depending on the distance between particle centers l s , both particles may experience reduced or increased drag, resulting in an increase or decrease in terminal velocity with respect to the settling of a single particle [START_REF] Akiki | Pairwise interaction extended point-particle model for a random array of monodisperse spheres[END_REF]. The settling velocity of the particle clearly also depends on the orientation of their separation, i.e. vertical or horizontal. Theoretical results for the terminal settling velocities as a function of the separation distance are available in [START_REF] Batchelor | Sedimentation in a dilute dispersion of spheres[END_REF][START_REF] Batchelor | Brownian diffusion of particles with hydrodynamic interaction[END_REF]. They are made non-dimensional by division through the settling speed of a single particle of same size in the Stokes regime:

r 2p s = u 2p p,s u p,s (36)
Note that the settling velocities are assumed aligned with a component of the coordinate system so that the previous relation involves only scalar quanti- Figure 12: Settling velocity histories at Re p = 50 for ∆ = (0.5, 1, 2) with source term regularization ("RegST") and removal of the particle's self-induced velocity based on the corrections of [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] ("BC"). One-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations and two-way coupled simulations with regularized source terms only ("2W + RegST") for ∆ = 1 are also shown for comparison purposes.

ties. Only the configuration with horizontal particle separation is studied in the present case. The simulation setup is exactly the same as for the settling of an isolated particle, see section 5.3.1Stokes regimesubsubsection.5.3.1. The simulations are performed for ∆ = (0.5, 1, 2). It is reminded that given the present choice of the regularization length scale L D,s = 6 d p , the analytic correction of the particle's self-induced velocity is applied for all three values of ∆ while the source term regularization is only enabled for ∆ = 1 and ∆ = 2. When simulating the particle pair problem, each particle experiences a repulsive force that increases the separation length l s in time and in turn modifies the terminal velocity [START_REF] Esmaily | A correction scheme for two-way coupled point-particle simulations on anisotropic grids[END_REF]. This increase in l s is not consistent with the analytic solution, which assumes fixed positions in the transverse direction. Therefore, in the simulations presented here, the particle separation in the transverse direction is held fix so that the particles are falling in straight lines. Three different separation lengths l s /d p = 2, 3, 4 are simulated, d p being the particle diameter. Simulation results are compared with theoretical results of the literature in table 3Settling velocity ratio of a pair of horizontally separated identical particles with respect to the settling velocity of a single isolated particle of the same diameter as a function of several nondimensional separation distances l s /d p . Numerical simulations are compared to theoretical results from [START_REF] Batchelor | Sedimentation in a dilute dispersion of spheres[END_REF][START_REF] Batchelor | Brownian diffusion of particles with hydrodynamic interaction[END_REF] (second column). Numerical simulation results are reported for ∆ = (0.5, 1, 2). [START_REF] Batchelor | Sedimentation in a dilute dispersion of spheres[END_REF][START_REF] Batchelor | Brownian diffusion of particles with hydrodynamic interaction[END_REF] (second column). Numerical simulation results are reported for ∆ = (0.5, 1, 2).

For the chosen regularization length scale of L D,s = 6 d p , the numerical perturbation flow field induced by the particle in the direction perpendicular to the applied force approximately matches the analytic solution beyond 2 diameters from the particle center, see fig. 8Fluid [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF], corresponding to L D,s = 6 d p . figure.caption.8, which seems consistent the overestimation of the settling velocity in the numerical simulations for this separation length, in particular for ∆ = 1 and ∆ = 2. However, when considering the evolution of error levels with respect to the theoretical solution for varying particle to mesh size ratios ∆ at fixed separation length scales, no clear trend appears. Also, error levels do not decrease uniformly with increasing separation length scale for ∆ = 0.5. These effects may be partly explained by the varying grid resolution of the perturbation flow field set by the particle. It is important to note that the present results were obtained with Faxén's corrective drag term [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF]; [START_REF] Gatignol | The Faxén formulas for a rigid particle in an unsteady non-uniform Stokes-flow[END_REF] based on the perturbed velocity field. The latter was found to have a significant impact on error levels for the largest particle to mesh size ratio ∆ = 2: while the average error level amounted to about -4% without Faxén's correction (not shown), it decreased to +2.5% with it. On the other hand, Faxén's correction did not significantly affect average error levels for ∆ = 0.5 and ∆ = 1. Resolution effects may again be partially invoked to explain this trend as curvature effects of the velocity field at the particle scale may only be captured with a sufficiently refined grid. [START_REF] Esmaily | A correction scheme for two-way coupled point-particle simulations on anisotropic grids[END_REF] also studied the settling of a horizontally aligned particle pair. Two of their simulations, namely for spherical particles with separation lengths l s /d p = (2, 4) with ∆ = 1 may be directly compared to those reported in the present work. For these two simulations, the error levels obtained in the present work are very similar to those obtained by [START_REF] Esmaily | A correction scheme for two-way coupled point-particle simulations on anisotropic grids[END_REF].

Conclusion

In the current work, a regularization procedure of the Lagrangian pointforce approach based on the numerical resolution of an unsteady nonlinear diffusion equation was presented. The proposed expressions for the diffusion coefficient allow to set the regularization length scale according to local particle diameter to mesh size ratios. Therefore, the present procedure may handle practical applications involving mesh refinement and polydisperse sprays. Moreover, it naturally ensures parallel efficiency when the particle source term regularization length scale spans over more than a few grid cells, as expected for anisotropic grids, particle to mesh size ratios of order unity and above as well as large particle source terms for which larger regularization length scales may be necessary to ensure numerical stability.

The regularization length scale is locally set by enforcing an approximately constant diffusion level in the vicinity of a considered particle / source term. In turn, the characteristic diffusion time scale may be approximately estimated from linear theory as the square root between diffusion coefficient and diffusion time. The proposed expression for the diffusion coefficient respects all constraints ensuring well-posedness and convergence to a unique solution of the partial differential equation. Using a standard centered discretization of the diffusion operator, convergence to second order accuracy under mesh refinement was demonstrated on structured grids and should be equally straightforward to obtain on unstructured grids.

Once implicitly regularized via the nonlinear diffusion equation, the source terms displayed shapes that were close to Gaussian curves. Therefore, the perturbation flow field obtained for the regularized source terms was close to the analytic solution obtained for Gaussian source term regularization by [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] in the steady Stokes regime. Their analytic solution was then used to retrieve the unperturbed fluid velocity by removal of the particle's self-induced velocity. Numerical simulations of the steady point-force problem showed that error levels on the estimate of this quantity due to deviation of the regularized source term from a Gaussian shape amounted to about 8%. Moreover, the leading fluid flow perturbation term was recovered beyond approximately two particle diameters from the particle's center of gravity for the chosen regularization length scale L D,s = 6 d p .

The source term regularization procedure and the removal of the particle's self-induced velocity were then used to simulate the settling of a particle under gravity in the Stokes regimes with two-way coupling enabled and particle diameter to mesh size ratios up to two. Using both ingredients, results appeared in good agreement for the terminal settling velocity, with error levels of about 3%. Note that source term regularization is disabled once the regularized source term becomes underresolved on the numerical grid, i.e. when it spans over a region comparable to the first neighbourhood of a computational cell, so that only removal of the particle's self-induced velocity is applied. This is because of a lack of precision to impose the regularization length scale in the context of underresolution. The settling problem was then also simulated for a pair of horizontally separated particles for three different particle to mesh size ratios, i.e. 0.5, 1 and 2. Even for the shortest separation length of two diameters between the particles centers of gravity, the resulting settling velocities were found in very good agreement with theoretical results of the literature [START_REF] Batchelor | Sedimentation in a dilute dispersion of spheres[END_REF][START_REF] Batchelor | Brownian diffusion of particles with hydrodynamic interaction[END_REF], regardless of the particle to mesh size ratio. For the largest value of the particle to mesh size ratio considered in the present work, i.e. 2, inclusion of Faxén's drag correction (based on the perturbed velocity field) lead to a reduction of error levels for the numerical simulations. On the other hand, the influence appeared much less pronounced for the two lower ratios, i.e. 0.5 and 1.

Although no precise quantitative recommendations can be given regarding the best choice of the regularization length scale, its value seems to involve a compromise between two opposing constraints. On the one hand, larger regularization length scales will decrease the magnitude of the perturbation velocity at the particle position, favoring numerical stability. On the other hand, the distance from the particle center at which the correct perturbation flow field is recovered will increase, particularly in the direction perpendicular to the force exerted upon the particle, deteriorating the prediction of neighbouring particle interactions or increasing the amount of unresolved interactions between particles that need to be modeled in turn. In the present work, the regularization length scale was then chosen as the smallest value which yielded a perturbation velocity at the particle position that lied below the relative velocity between fluid and particle for the steady point force problem, in loose analogy with a maximum principle.

Since the analytic correction of [START_REF] Maxey | Localized force representations for particles sedimenting in stokes flow[END_REF] is no longer valid at finite particle Reynolds numbers, corrections proposed by [START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] were used to simulate the settling problem at finite particle Reynolds number. The aforementioned corrections appeared compatible with the proposed methodology as very good agreement was observed for settling problems simulated at particle Reynolds number 10 and 50 for three different particle to mesh size ratios.

While these results seem encouraging, many issues need to be adressed. In order to make the present methodology affordable for realistic applications in the presence of large particle to mesh size ratios, typically larger than two, a more efficient and scalable resolution technique of the nonlinear diffusion equation needs to be found. However, it is unclear whether this equation may be solved more efficiently with implicit iterative methods. If no efficient means to solve this equation is found, it could be reformulated by decomposing the regularization procedure into two distinctive substeps: the first substep would be used to set the field of spatially varying diffusion coefficients, while the second would perform the actual regularization based on this diffusion coefficient field. While the first substep could be handled with a more general, not necessarily conservative transport equation, the second substep could be solved with standard implicit iterative schemes. In this context, a comparison of implementation efforts, accuracy and computational cost of the methodology in its present formulation with respect to a scalable explicit regularization procedure as proposed by Zwick and Bal-achandar (2019) on a case of intermediate complexity would certainly be of interest. From a general perspective, it appears that the present method solves a Eulerian conservation equation to perform the regularization step and that particle positions are only used to provide the initial conditions for the regularization problem. Therefore, the cost of the present procedure should not be significantly affected by the number of particles in the simulation, but only the maximum particle to mesh size ratio, contrary to explicit scalable regularization methods which are however computationally cheap at the particle scale [START_REF] Zwick | A scalable Euler-Lagrange approach for multiphase flow simulation on spectral elements[END_REF]. Therefore, similarly to Eulerian and Lagrangian approaches for dispersed phase modeling, both methods display somewhat complementary advantages. In addition to performance considerations, the accuracy of the present method should also be assessed on test cases of intermediate complexity such as homogeneous isotropic turbulence laden with particles. Changes in both particle and turbulence dynamics were reported in the literature for this test case, for instance regarding the rate of energy dissipation [START_REF] Esmaily | A correction scheme for two-way coupled point-particle simulations on anisotropic grids[END_REF], particle settling velocities [START_REF] Ireland | Improving particle drag predictions in Euler-Lagrange simulations with two-way coupling[END_REF] and particle preferential concentration [START_REF] Horwitz | Correction scheme for point-particle models applied to a nonlinear drag law in simulations of particle-fluid interaction[END_REF] when including corrections to recover the undisturbed fluid velocity at the particle position. It would be interesting to confirm some of these observations using the present approach. The final aim is to implement the proposed methodology in an unstructured Navier-Stokes solver. In such context, the accuracy levels than can be achieved with the proposed methodology on anisotropic grids need to be considered. Furthermore, performing numerical simulations of realistic gas particle systems such as solid propulsion rockets and aeronautic combustion chambers will require extending the regularization procedure to account for mass and energy source terms [START_REF] Rangel | An evaluation of the point-source approximation in spray calculations[END_REF]; [START_REF] Liu | Self-induced temperature correction for inter-phase heat transfer in Euler-Lagrange point-particle simulation[END_REF]. Finally, the regularization procedure also seems applicable to stochastic point particle simulations. In the latter, the regularization length scale is no longer imposed by physical but statistical considerations, typically the innerparticle distance set by the initial sampling of the spray and the evolution of this quantity during particle transport [START_REF] Norde | Eulerian and Lagrangian ice-crystal trajectory simulations in a generic turbofan compressor[END_REF]. With some modifications, the proposed regularization procedure should also be applicable in the statistical context, where it could possibly alleviate convergence issues of point-particle approaches [START_REF] Garg | A numerically convergent Lagrangian-Eulerian simulation method for dispersed two-phase flows[END_REF] at a reduced computational cost. as two-way coupled simulations with regularized source terms without ("2W + RegST") and with removal ("2W + RegST + AC") of the particle's selfinduced velocity. The final simulation ("2W + AC") concerns a PF-PIC two-way coupled simulation without source term regularization but with the removal of the particle's self-induced velocity.figure.caption.19 with an additional test case. The latter consists in a standard PF-PIC two-way coupled simulation for ∆ = 1 with application of the velocity correction procedure. Again, the regularization length scale of the Gaussian projection is obtained assuming that the associated regularization length scale is L D,s = 2∆x. As may be seen from fig. B.15Settling velocity histories for ∆ = 1. Comparison of one-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations as well as two-way coupled simulations with regularized source terms without ("2W + RegST") and with removal ("2W + RegST + AC") of the Stokes' analytic solution for the equivalent drag force is also plotted for comparison.

particle's self-induced velocity. The final simulation ("2W + AC") concerns a PF-PIC two-way coupled simulation without source term regularization but with the removal of the particle's self-induced velocity.figure.caption.19, pronounced oscillations on the settling velocity are visible for this test case. Moreover, the average error level on the settling velocity seems significant. Adjustments on the velocity correction were tested to reduce the error level, but they appeared to magnify the oscillations (not shown). Therefore, the main issue seems to lie in the increased sensitivity of the particle to the applied fluid velocity correction as the latter exceeds the relative velocity between fluid and particle in magnitude. Clearly, the occurrence of these oscillations depends on numerical details of the considered Navier-Stokes solver and employed particle source term projection scheme. Therefore, while quantitative recommendations on the best choice of the regularization length scale

Figure 1 :

 1 Figure 1: Evolution of the diffusion coefficient given by eq. 16Expression of the diffusion coefficientequation.3.16 for several values of the parameter α (at constant β).

Figure 2 :

 2 Figure2: Evolution of the diffusion coefficient given by eq. 16Expression of the diffusion coefficientequation.3.16 for several values of the parameter β with ∆ = 1 (at constant α).

Figure 4 :

 4 Figure 4: L 2 and L ∞ errors of the regularized source terms under mesh refinement with respect to the solution obtained on the finest mesh, i. e. ∆ = 16.

Figure 5 :

 5 Figure 5: Regularized source term field and Gaussian reference curves with the same regularization length scale for multiple source terms. Both curves are normalized by the maximum of the Gaussian reference curves.

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Diffusion coefficient field for the multiple source term test case at five different non-dimensional instants τ * = τ /τ max .

  perturbation in the direction of the force for ∆ = 2. Perturbations resulting from different regularization length scales L D,s = 4 d p , 6 d p , 8 d p are compared to Stokes' analytic solution and the solution of Maxey and Patel (2001), corresponding to L D,s = 6 d p .figure.caption.7.

Figure 10 :

 10 Figure10: Settling velocity histories in the Stokes regime for varying ∆. Comparison of one-way ("1W") with two-way coupled simulations with source term regularization and with analytic removal ("2W + RegST + AC") of the particle's self-induced velocity.

Figure 11 :

 11 Figure11: Settling velocity histories at Re p = 10 for ∆ = (0.5, 1, 2) with source term regularization ("RegST") and removal of the particle's self-induced velocity based on the corrections of[START_REF] Balachandar | Self-induced velocity correction for improved drag estimation in Euler-Lagrange point-particle simulations[END_REF] ("BC"). One-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations and two-way coupled simulations with regularized source terms only ("2W + RegST") for ∆ = 1 are also shown for comparison purposes.

  perturbation in the direction perpendicular to the force for ∆ = 2. Perturbations resulting from different regularization length scales L D,s = 4 d p , 6 d p , 8 d p are compared to Stokes' analytic solution and the solution of Maxey and Patel (2001), corresponding to L D,s = 6 d p . figure.caption.8. At l s /d p = 2, the perturbation seems slightly overestimated in fig. 8Fluid perturbation in the direction perpendicular to the force for ∆ = 2. Perturbations resulting from different regularization length scales L D,s = 4 d p , 6 d p , 8 d p are compared to Stokes' analytic solution and the solution of

  Figure B.13: Fluid perturbation in the direction of the force for the steady point-force problem using the standard PF-PIC approximation with ∆ = (0.5, 1, 2). Stokes' analytic solution for the equivalent drag force is also plotted for comparison.

  Figure B.14: Fluid perturbation in the direction perpendicular to the force for the steady point-force problem using the standard PF-PIC approximation with ∆ = (0.5, 1, 2). Stokes' analytic solution for the equivalent drag force is also plotted for comparison.

  table.caption.15.

	l s /d p r 2p s , theo.	∆ = 0.5	∆ = 1	∆ = 2
	2	1.1950	1.1974 (+0.2%) 1.2228 (+2.3%) 1.2303 (+2.9%)
	3	1.1273	1.1526 (+2.2%) 1.1382 (+0.96%) 1.1540 (+2.4%)
	4	1.0947	1.0788 (-1.5%) 1.1019 (+0.65%) 1.1204 (+2.3%)

Table 3 :

 3 Settling velocity ratio of a pair of horizontally separated identical particles with respect to the settling velocity of a single isolated particle of the same diameter as a function of several nondimensional separation distances l s /d p . Numerical simulations are compared to theoretical results from

Acknowledgements

The authors would like to warmly thank Jesse Capecelatro for fruitful discussions on volume filtering approaches in the context of Euler-Lagrange simulations. They are also grateful to Ghislain Blanchard for his insightful comments on the present work.

AppendixA. Properties of the Perona-Malik equation

The nonlinear diffusion problem proposed by [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF] writes: Well-posedness of this initial value problem is ensured if the flux function is monotonously increasing in s [START_REF] Brezis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] [START_REF] Walter | Differential and integral inequalities[END_REF].

Considering s to be a scalar quantity, the derivative of the flux function writes:

Therefore, well-posedness is ensured if the diffusion coefficient and its first order derivative need to be strictly positive functions. Note that proofs of existence and uniqueness for eqs. A.1Properties of the Perona-Malik equationequation.A.1-A.3Properties of the Perona-Malik equationequation.A.1 were even established for nonmonotone but regularized flux functions [START_REF] Catté | Image selective smoothing and edge detection by nonlinear diffusion[END_REF]; [START_REF] Alvarez | Image selective smoothing and edge detection by nonlinear diffusion[END_REF].

AppendixB. Necessity of the regularization procedure

The present appendix provides a few more details regarding stability issues of the standard PF-PIC approximation.

This first test case simulates the flow field response to the application of a punctual force described in section 5.2Coupling with the fluid phasesubsection.5.2. All numerical parameters are identical except that no regularization is performed. Instead, the standard PF-PIC approximation, where the entire source term is reported to the cell containing the particle's center, is used. The simulation is performed for the particle to mesh size ratios ∆ = (0.5, 1, 2). The flow field perturbations in the directions of the force and perpendicular to the force are respectively plotted on figs. B.13Fluid perturbation in the direction of the force for the steady point-force problem using the standard PF-PIC approximation with ∆ = (0.5, 1, 2). Stokes' analytic solution for the equivalent drag force is also plotted for comparison.figure.caption.17 and B.14Fluid perturbation in the direction perpendicular to the force for the steady point-force problem using the standard PF-PIC approximation with ∆ = (0.5, 1, 2). Stokes' analytic solution for the equivalent drag force is also plotted for comparison.figure.caption.18. Stokes' analytic solution is also shown for comparison.

As expected, the magnitude of the fluid velocity perturbation at the particle center r = 0 significantly increases with ∆ and exceeds three times the equivalent relative particle to fluid velocity for ∆ = 2. On the other hand, the far field solution is not affected by the divergence occurring at the origin, which is consistent with theoretical findings indicating that the convergence of the stream function at large distances should not be affected by the singularity at the origin [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]. In particular, the fluid perturbations obtained with the standard PF-PIC approximation seem to yield accurate solutions beyond approximately r y > 2d p in the direction of the force and r x > d p in the direction perpendicular to the force. Since regularization was shown to increase the distance from the particle for which the correct fluid perturbation flow field is recovered, this confirms that only the minimal amount of regularization ensuring stability should be applied.

In order to asses the influence of the increasing velocity perturbation at the particle position, results obtained for the single particle settling problem under gravity in section 5.3.1Stokes regimesubsubsection.5.3.1 are reported in fig. B.15Settling velocity histories for ∆ = 1. Comparison of one-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations as well seem difficult to make, the latter may be mandatory for reasons of numerical stability. Comparison of one-way ("1W") and standard PF-PIC two-way ("2W") coupled simulations as well as two-way coupled simulations with regularized source terms without ("2W + RegST") and with removal ("2W + RegST + AC") of the particle's self-induced velocity. The final simulation ("2W + AC") concerns a PF-PIC two-way coupled simulation without source term regularization but with the removal of the particle's self-induced velocity.