Sahan Kuruneru

Ewen Maréchal

Michaël Deligant

Emilie Sauret

Sofiane Khelladi

Florent Ravelet

Chandra Suvash

Yuantong Saha

Gu

A Comparative Study of Mixed Resolved-Unresolved CFD-DEM and Unresolved CFD-DEM Methods for the Solution of Particle-Laden Liquid Flows

published or not. The documents may come

Various computational techniques have been developed to investigate the efficacy of heat exchangers. Ramagadia et al. [START_REF] Abhishek | Fully developed flow and heat transfer characteristics in a wavy passage: Effect of amplitude of waviness and reynolds number[END_REF] used finite volume method (FVM) together with a momentum interpolation method to investigate heat transfer characteristics of wavy channel (i.e. sinusoidal or arc-shaped walls) heat exchangers. Wavy channels are steadily gaining attention thanks to their manufacturing simplicity and potentially high energy savings and less power consumption. It was found that wavy-type heat exchangers exhibit higher heat transfer rates compared with straight channels due to the unsteady vortex shedding in the former device. Wang et al. [START_REF] Wang | Real-time fouling characteristics of a typical heat exchanger used in the waste heat recovery systems[END_REF] used a discrete phase model (DPM) coupled with a RNG k-ε turbulence model to investigate the real-time fouling characteristics of a H-type finned tube for waste heat recovery applications. It was found that fouling mainly occurs in the flow stagnation region specific regions such as in front of the tube and fins. The asymptotic fouling resistance decreases with increasing superficial inlet velocity, and fouling without removal increases linearly with time. However, the DPM model did not include full resolved particle-particle and particle-fluid interactions which is must be enabled in dense particle-fluid flows; moreover, the morphology variation caused by the foulants on the finned tubes was not considered. Bayomy et al. [START_REF] Bayomy | Electronic cooling using water flow in aluminum metal foam heat sink: Experimental and numerical approach[END_REF] studied the optimum design of an aluminium foam computer heat sink by taking into account the highest heat transfer and lowest pumping power. De Bellis & Catalano [START_REF] Bellis | Cfd optimization of an immersed particle heat exchanger[END_REF] used Reynolds averaged Navier Stokes equations (RANS) together with SIMPLEX and non-dominated sorting genetic algorithm II to maximize the heat transfer of an immersed particle based heat exchanger. Pierre et al. [START_REF] Pierre | Numerical computation of 3d heat transfer in complex parallel heat exchangers using generalized graetz modes[END_REF] developed an optimal weak-variational formulation in the form of a spectral method (i.e. generalized Graetz problem) for the numerical analysis of the temperature fields and effectiveness of parallel convective heat exchangers. In the case of a two inlet/outlet semi-infinite counter-current tubes, it was found that the heat exchanger effectivess saturates with the exchanger length and Péclet number. The final effectiveness is controlled by the thermal conditions though the dependence on the imposed hydrodynamics is diminutive. Gu et al. [START_REF] Gu | Airside heat transfer and pressure loss characteristics of bare and finned tube heat exchangers used for aero engine cooling considering variable air properties[END_REF] stated that the exclusion of the effects of the variability of air properties on the thermal-hydraulic characteristics of heat exchangers for aero-engine cooling can overestimate the heat transfer and pressure drop.

The major limitation of the cited publications is that the evaluation of heat exchanger performance is based on several oversimplified assumptions which make a systematic comprehension of complex multiphase transport and particulate fouling phenomena impossible. Firstly, the above cited publications that assess a heat exchangers efficiacy with the use of the Eulerian-Lagrangian or even with the use of CFD-DPM are based on assumptions such as neglecting fully resolved particleparticle interactions (i.e. zero particle volume) and neglecting the influence of the dispersed phase on the fluid continuum, and vice versa; which is incorrect considering the fact that heat exchangers in the chemical, oil & gas, and energy generation industries consist of dense multiphase (i.e. solid-liquid or solid-gas) flows [START_REF] Stw Kuruneru | Numerical investigation of the temporal evolution of particulate fouling in metal foams for air-cooled heat exchangers[END_REF]. Secondly, studies assume the system comprises single-phase flow which is not the norm in a myriad of engineering applications such as air-cooled heat exchangers [START_REF] Stw Kuruneru | Numerical investigation of the temporal evolution of particulate fouling in metal foams for air-cooled heat exchangers[END_REF][START_REF] Zf Tian | Numerical modelling and validation of gas-particle flow in an in-line tube bank[END_REF], or even fuel systems in turbojets [START_REF] Marchal | A hybrid finite volumediscrete elements for two-phase flows: application to snow showers in jet-engine fuel systems[END_REF]. Therefore, the development of robust numerical models is of paramount significance to decipher the mechanisms that govern multiphase transport and fouling in various engineering systems.

Numerical methods for particle-laden fluid flows

Studies have delved into the physics of multiphase transport using various advanced numerical techniques, each one having its own strengths and weaknesses. The standard Lagrangian based DPM (Discrete Phase Model) is similar to the Discrete Element Method (DEM) but the former neglects inter-particle collisions (i.e. zero particle volume); moreover, the DPM method neglects gas displacement by the particles. As such, the DPM is suited for dilute particulate suspensions [START_REF] Kloss | Comprehensive demdpm-cfd simulations-model synthesis, experimental validation and scalability[END_REF] where a larger time-step could be used to reduce computational effort [START_REF] Mezhericher | Modeling of particle pneumatic conveying using dem and dpm methods[END_REF]. The Multiphase Particlein-Cell (MPPIC) method is similar to the DEM methodology but particle movement and interactions are viewed statistically whilst excluding particle-particle and particle-wall interactions. The MPPIC method could also be deployed in dense solutions as it is impractical to use DEM method to simulate tens of thousands or millions of discrete particles. The Eulerian-Eulerian (i.e. two fluid model) method is not suitable to model dense (non-dilute) particle-fluid flows as the method treats both phases as an interpenetrating continua; moreover, the constitutive relations for solid particles and interphase interactions are generally not available [START_REF] Hp Zhu | Discrete particle simulation of particulate systems: theoretical developments[END_REF]. Monte-Carlo methods do not permit one to accurately resolve all particle-fluid interactions particle micromechanics; moreover, the visualization of particle accumulation is not achievable. Population balance approaches, which is originally based on the Smoluchowski equation, is found to provide accurate predictions of particle aggregate size distribution; the drawback of this method is that predicting or obtaining the actual micromechanics and microstructure of the aggregation phenomena and particle aggregate interactions with fluid cannot be achieved [START_REF] Marshall | Discrete-element modeling of particulate aerosol flows[END_REF]. The Lattice Boltzmann-Discrete Element Method (LB-DEM) could also be used to examine particle-fluid flows. However, the development of such models is still at its infancy. Moreover, it is significantly more computationally demanding than the widely used FVM-DEM method or the Two Fluid Model (TFM) method. It is also not well suited for specific engineering applications such as process modeling and control but it is suitable for fundamental research on particle physics. Moreover, numerical difficulties are inherent in LB-DEM simulations with strong particle-particle interactions [START_REF] Hp Zhu | Discrete particle simulation of particulate systems: theoretical developments[END_REF][START_REF] Yu | Powder processing models and simulations[END_REF]. The literature is devoid of material regarding the development of penalty methods. The development of these methods could then be used to compare against the results of the other multiphase numerical methods such as the CFD-DEM method based on [START_REF] Stw Kuruneru | A coupled finite volume and discrete element method to examine particulate foulant transport in metal foam heat exchangers[END_REF]. In particular, there are no studies on penalization on solid-liquid flows; moreover, investigation of solidliquid flows using a CFD-DEM is very rare and most of these CFD-DEM investigations on solid-liquid flows are based on fluidized beds. What is more, the literature is devoid of CFD-DEM investigation of solid-liquid flows in complex geometries. Clearly, the development of numerical algorithms to fully resolve multiphase transport in various engineering systems is extremely complex and challenging. [START_REF] Tryggvason | Virtual motion of real particles[END_REF].

Objective

Fouling phenomenon is indeed a complex multifaceted problem and it is a subject of debate and ongoing research. This provides the prime motivation to develop robust and advanced numerical methods to accurately unravel the mechanisms governing multiphase transport (i.e. solid-gas or solid-liquid) and particulate fouling. The development and successful implementation of advanced numerical methods permits engineers, for instance, to better optimize heat exchanger systems for the purposes of fouling alleviation and control. The objective of this investigation is to develop and compare two numerical method to assess the mechanisms that govern two-phase solid-liquid flows and particulate clogging on a filter. The two numerical methods is validated extensively against various benchmarks namely, experimental results.

Numerical methods

We investigate solid-liquid flows and solid particle accumulation on a filter based on two numerical methods. Section 2.1. covers the equations which is used in both methods. Section 2.2 covers the mathematical moedl based on penalization whereas the second method is based of the coupled CFD-DEM method developed in Open Field Operation and Manipulation (OpenFOAM), an opensource C++ CFD program. The results based on these two methods are compared against the experimental results.

Constitutive equations

The transport of incompressible and isothermal fluid is governed by the Navier-Stokes equations, and is given as

∂α d ∂t + ∇ • (α d v) = 0 (1)
∂(ρ f α d v) ∂t + ∇ • (ρ f α d vv) = -∇(α d p) + ∇ • (α d τ ) + (ρ f α d g) + F pf (2)
where α d is the fluid volume fraction, fluid velocity v, fluid density ρ f , fluid pressure p, gravitational acceleration g, fluid viscous stress tensor τ . The gas volume fraction is given as

α d = 1 - kc i=1 V i V c (3)
The interphase momentum transfer between the particles and fluid is denoted by F pf [START_REF] Stw Kuruneru | A coupled finite volume and discrete element method to examine particulate foulant transport in metal foam heat exchangers[END_REF] and is given as

F pf = 1 ∆V c kc i=1 f pf,i (4)
The presence of non-reactive dense particulate-fluid flows signifies the necessity to accurately resolve interactions between the individual particles and the hydrodynamic interactions between the particles and carrier fluid (two-way coupling) and particle and the walls of the domain (four-way coupling). The cohesive contacts (i.e. cohesion energy density) between the particles and walls have not been taken into account as the particle diameter is significantly greater than 1 mm. In this study, the density and diameter of a particle is set to 2400 kg/m 3 and 1-2 mm respectively.

Mesh cell size and DEM particle size

Although it is desirable to having a fine mesh in order to resolve details of the fluid flow field, numerical stability becomes concern if a DEM particle size approaches a FVM cell of identical size (or slightly smaller). Smoothing or approximation models could be deployed to circumvent this issue. The reader is referred to Goniva et al [START_REF] Kloss | Models, algorithms and validation for opensource dem and cfd-dem[END_REF] regarding other approximation methods (such as big particle void fraction method, divided void fraction method, etc) for use in the event a DEM particle is either very similar to the cell size or slightly exceeds the cell size. In this study an approximation method based on Wahyudi et al [START_REF] Hadi Wahyudi | 3d particle-scale modeling of gas-solids flow and heat transfer in fluidized beds with an immersed tube[END_REF] is developed in OpenFOAM. This approximation method is based on the designation of the maximum solid phase fraction per computational cell if the FVM cell is completely smeared with the DEM particle (around the vicinity of the filter); also, and the mass and momentum sources is distributed to neighbouring cells as a means to conserve mass and energy [START_REF] Hadi Wahyudi | 3d particle-scale modeling of gas-solids flow and heat transfer in fluidized beds with an immersed tube[END_REF]. The mesh cell size S c to particle size d p ratio is approximately 5:1 but this ratio is significantly lower at the filter region. Each cell is coloured according to its solid volume fraction χ s . The difference between the particle size and the resolution of the mesh has an impact on the choice of the coupling method. A) Most particles are smaller than the mesh cell size, therefore an unresolved approach is suitable. B) The particles and the mesh cells have similar size. Neither approach is appropriate but the hybrid method is applicable. C) The flow around the largest particle can be resolved accurately.

It is noteworthy that there is no concrete consensus among the research community regarding the preferred ratio of the mesh cell size to the particle diameter; values vary in literature. For instance, Li et al [START_REF] Li | Modeling of spout-fluidized beds and investigation of drag closures using openfoam[END_REF] suggested a ratio of 1.670 or above whereas both Geng & Che [START_REF] Geng | An extended dem-cfd model for char combustion in a bubbling fluidized bed combustor of inert sand[END_REF] and Wang et al [START_REF] Wang | Simulation of flow behavior of liquid and particles in a liquid-solid fluidized bed[END_REF] successfully used 1.00 and 1.33 respectively, whereas Feng & Yu [START_REF] Yq Feng | Effect of bed thickness on the segregation behavior of particle mixtures in a gas fluidized bed[END_REF] used 1.625. Kubicki & Lo [START_REF] Kubicki | Slurry transport in a pipeline-comparison of cfd and dem models[END_REF] numerically examined slurry transport with an S c /D p ratio of approximately one and it was found the numerical results are in good agreement with the experimental data. Kuruneru et al. [START_REF] Stw Kuruneru | Numerical investigation of the temporal evolution of particulate fouling in metal foams for air-cooled heat exchangers[END_REF] showed no difference in the particulate foulant distribution profiles in an idealized metal foam heat exchanger pertaining to a mesh cell size/particle diameter ratio of 1:1 or higher. What is more, Li et al [START_REF] Li | Modeling of spout-fluidized beds and investigation of drag closures using openfoam[END_REF] claimed that a ratio of 1.67 or higher is suitable (at least for their study) as the results corresponding to about 50 data points are all identical to the experimental results; interestingly though, about 68 % and 74 % of the data points corresponding to ratios of 1.02 and 1.12 respectively, closely match the experimental data points. All 48 data points for ratios of 1.67 or higher identically match the experimental data Wahyudi et al [START_REF] Hadi Wahyudi | 3d particle-scale modeling of gas-solids flow and heat transfer in fluidized beds with an immersed tube[END_REF] used very fine grids near the wall of a fluidized bed in order to thermally resolve the boundary layer. It is in this region that the particle diameter is smaller than the mesh cell size (half the particle diameter). However, in the study by Li et al [START_REF] Li | Modeling of spout-fluidized beds and investigation of drag closures using openfoam[END_REF], the minimum gas phase fraction was not specified; moreover, the 1.67 ratio is based on the Gidaspow drag law (ErgunWenYu) [START_REF] Gidaspow | Multiphase flow and fluidization: continuum and kinetic theory descriptions[END_REF], as such it remains to be seen whether the same ratio stands for other drag closures (i.e. Di Felice [START_REF] Di | The voidage function for fluid-particle interaction systems[END_REF], PlessisMasliyah [START_REF] Prieur | Flow through isotropic granular porous media[END_REF], Koch & Hill [START_REF] Donald | Inertial effects in suspension and porous-media flows[END_REF]) or a standard drag model. What is more, the Gidaspow model is not universally used in all dense particle-fluid systems. To be precise, the Gidaspow is ideally suited for packed beds whereas the Di Felice drag closure is derived for particle sedimentation. Thus, the 1:67 ratio ideally should not be used in all dilute or dense granular-fluid systems. The interested reader is referred to [START_REF] Stw Kuruneru | Coupled cfd-dem simulation of oscillatory particle-laden fluid flow through a porous metal foam heat exchanger: Mitigation of particualte fouling[END_REF] regarding grid size to DEM particle diameter ratio and validation of numerical model pertaining to particel bounce in a rectangular model compared with experimental data. In fact, there is no concrete consensus among the research arena regarding the modeling of solid-air drag closures [START_REF] Ebrahimi | Numerical and experimental study of horizontal pneumatic transportation of spherical and low-aspect-ratio cylindrical particles[END_REF][START_REF] Sturm | Coupled dem-cfd simulation of pneumatically conveyed granular media[END_REF][START_REF] Mjv Goldschmidt | Hydrodynamic modelling of dense gasfluidised beds: comparison and validation of 3d discrete particle and continuum models[END_REF].

Mixed resolved-unresolved FVM-DEM method

The methods presented in this section are implemented in a research C++ based CFD program. The present method has been developed to provide a tool capable of simulating these transient, complex flows in actual fuel system geometries [START_REF]Etude du colmatage des systmes carburant de turboracteurs par des suspensions denses de particules de glace[END_REF] Jet engines have complex fuel systems, involving several hydraulic components sensitive to clogging e.g. filters, valves or heat exchanger inlet screens. When a so-called "snow shower" occurs, ice particles settle in seconds to form a porous layer which is likely to occasion fuel flow restrictions, as shown in Figure 2. Little data were then available, which is why the model was designed to adapt to various time and space scales of this phenomenon. The fluid phase is modelled by equations ( 1) and ( 2) with α d = 1 in the entire domain. The particularity of this method comes from the asymmetry in the treatment of the interphase momentum transfer : The action of the flow on the particles F f p is modelled as in typical unresolved approach, while the sink term in the momentum equation F pf is calculated based on the local medium porosity, in the manner of resolved methods. For a complete description of the method, see [START_REF] Marchal | A hybrid finite volumediscrete elements for two-phase flows: application to snow showers in jet-engine fuel systems[END_REF] and [START_REF]Etude du colmatage des systmes carburant de turboracteurs par des suspensions denses de particules de glace[END_REF].

Flow through porous medium

The porous medium consists of a stack of ice crystals of various sizes. The feedback of the particles on the flow must be considered at two scales :

• Microscopic: ice particles are intrinsically porous. Experiments have shown that the fraction of fluid within the ice can reach up to 50

• Macroscopic: Even if the particles were impermeable, the stack of particles has interstices in which the fluid can flow.

The pressure drop related to the intrinsic porosity is modelled by Darcy's law. For a given finite volume cell of the mesh, the pressure gradient related to viscous loss is expressed:

∇p = - µ K (v f -v s ) ( 5
)
Where µ is the dynamic viscosity of the fluid, v f the superficial velocity of the fluid, v s is the (average) solid velocity within the cell, and and K is the intrinsic permeability. To account for the macroscopic effects, the equation ( 5) is weighted by a permeability function X , calculated for each mesh cell :

∇p = -X (χ s ) µ K (v f -v s ) (6)
Where χ s is the solid volume fraction, defined as the ratio of the total volume of particles entirely or partially located in a computational cell, to the volume of the FVM cell :

χ s = 1 -α d = kc i=1 V i ∆V c . (7)
The permeability function used in the present work is a power law of exponent η = 2, which acts as a tuning parameter for the macroscopic porosity. The final expression for the momentum sink term is :

F pf = -χ η s µ K (v f -v s ) (8)

Hydrodynamic forces

The forces considered in this work include the effects of pressure field, buoyancy and drag. For the sake of simplicity, the other hydrodynamic forces and the effects related to the particles rotation are not taken into account (see table 1). Since the particles accumulate, their interactions with each other and with the walls must also be considered. The contact force between a particle i and an other solid k in the simulation is denoted F c ik . The constitutive laws that govern the transport of the individual discrete particles are based on Newtons second law given as:

m i dv i dt = F g + F d + F p + k =i F c ik (9)
In this study, the particles are assumed smooth, rigid, and isothermal.

Particle interactions

There are two main approaches for modelling contact particle interactions. The first relies on the conservation of the momentum of binary and instantaneous shocks, and generally implements algorithms based on the management of collision events. It is well suited for inelastic collision of hard

Buoyancy F g (m p -ρ f V p ) g Drag F d -1 2 ρ f S p C D v r v r Pressure F p -V p (∇p) hydro
Table 1: Hydrodynamic forces considered in the present method.

spheres. In event-driven algorithms, the time step is determined by the smallest duration between two contacts, which tends to zero when the solid fraction increases. In the second approach, two particles are considered in contact when they interpenetrates slightly. Normal and tangential forces are then evaluated by a spring-dashpot model.

Industrial applications of clogging require simulations over periods of several seconds and involves locally high concentrations of particles. The very small time-step of conventional models can become problematic. A contact algorithm for inelastic collisions allowing for larger time steps was thus developed. The present method consider that the collisions between particles are perfectly inelastic. This choice is based on two arguments : 1) The coefficient of restitution for wet particles and particles in fluid is much lower than dry particles [START_REF] Liu | Cfd-dem simulation of liquid-solid fluidized bed with dynamic restitution coefficient[END_REF] 2) The experimental setup is designed to study static stacks of particles. The principle of the model is the following : let consider a pair of particles (i, r i , m i ) and (j, r j , m j ) where r and m refers to radius and mass. In absence of any interaction force, interpenetration between two particles at the end of a time step ∆t may occur (figure ). The principle of the contact handling algorithm is to compute the contact force F c ij required to prevent this behaviour for each pair of particles and to apply it before the position of the particles are updated.

The derivation of the method starts with the discrete equation of motion for a pair of particles :

           x n+1 i -x n i ∆t = v n i + ∆t m i F h i + k =i F c ik x n+1 j -x n j ∆t = v n j + ∆t m j F h j + k =j F c jk ( 10
)
Where F h i refers to the sum of all but contact forces. Here a 1st-order forward euler formulation is used but the method can be derived with more refined time schemes. The distance that would cover a particle if there was no contact divided by the time step dt is called predicted velocity and denoted v. The expression reads:

vn i = v n i + ∆t m i F h i (11)
Thereafter, let us denote any relative quantity φ ij = φ j -φ i . The two vector equation of motion are subtracted from one another and projected on the normal direction n ij to obtain a single scalar equation :

x n+1 ij • n ij = x n ij • n ij + ∆tv ij • n ij + k =i ∆t 2 m i F c ik + k =j ∆t 2 m j F c jk • n ij (12)
x ij = (x j -x i ) • n ij is the relative distance between two particle and vij = (v j -vi ) • n ij the relative predicted velocity projected along n ij . Decomposing contact forces by their amplitude f ik and direction n ik and introducing previous notations, equation [START_REF] Geng | An extended dem-cfd model for char combustion in a bubbling fluidized bed combustor of inert sand[END_REF] becomes:

x n+1 ij = x n ij + vij ∆t + k =i ∆t 2 m i f ik n ik + k =j ∆t 2 m j f jk n jk • n ij (13)
Avoiding the interpenetration at the end of the next time step is equivalent to verifying for each pair of particle (i, j) the inequality:

x n+1 ij ≥ r i + r j (14)
Combining eq [START_REF] Gidaspow | Multiphase flow and fluidization: continuum and kinetic theory descriptions[END_REF] and eq [START_REF] Mjv Goldschmidt | Hydrodynamic modelling of dense gasfluidised beds: comparison and validation of 3d discrete particle and continuum models[END_REF] finally gives:

x n ij + vij ∆t + k =i f ik ∆t 2 m i n ik • n ij + k =j f jk ∆t 2 m j n jk • n ij ≥ r i + r j (15)
The distance between the surface of two particles is usually denoted δ ij . For better readability, let us denote:

Λ j ik = ∆t 2 m i n ik • n ij (16)
Equation ( 15) holds for each pair of particles. A system of inequalities g(f ij ) is obtained:

g(f ij ) = δ ij -vij ∆t - k =i Λ j ik f k ik + k =j Λ i jk f k jk ≤ 0 (17)
Adding contact forces will impact the total kinetic energy of the particle system. Therefore contact forces must be computed so that the set of constraint is satisfied while minimizing the change of kinetic energy [? ]. This is achieved through an iterative procedure, fully described in [START_REF]Etude du colmatage des systmes carburant de turboracteurs par des suspensions denses de particules de glace[END_REF]. Starting from zero, contact forces are gradually increased in proportion of the value of the constraint functions. After convergence, the speed of a particle a time t + ∆t is computed with the resulting interactions forces:

v n+1 i = vn i - ∆t m i k =i f ik n ik (18)
In practice, the convergence of the method requires an under-relaxation of parameter ω. The iterative procedure is as follow :

1. Initialisation of forces and constraints:

f 0 ij = 0 and g 0 ij = δ ij -vij ∆t 2. Update forces : f k+1 ij = max 0; f k ij -ω.g k ij 3. Update constraints : g k+1 ij = δ ij -vij ∆t - k =i Λ j ik f k+1 ik + k =j Λ i jk f k+1 jk 4. Convergence test 2.2.

Numerical schemes and resolution algorithm

The complexity of solving the Navier-Stokes equations lies in the absence of an independent pressure equation, whose gradient plays a dominant role in the three momentum equations. Moreover, the continuity equation no longer behaves like a transport equation for the mass, but plays the role of constraint on the velocity field. Thereby the main difficulty in solving the incompressible Navier-Stokes equations is the calculation of the pressure field. The resolution algorithm adopted in this work is based on the SIMPLE algorithm, first developed by Caretto et al. in 1972, and adapted to unstructured meshes by Rhie and Chow in 1983. The reader is referred to Ferziger and Peric [START_REF] Joel | Computational methods for fluid dynamics[END_REF] regarding the standalone SIMPLE or PISO algorithms.

Due to the non-linearity of the equations, this class of algorithms involve an iterative procedure to achieve coupling between pressure and velocity fields. The deferred correction is a technique whose principle is to calculate higher order terms explicitly, and treat them as a source term in the second member of the equation. A low-order approximation of these terms is treated implicitly in the first member, and subtracted from the second member. As the iterative procedure converge, low-order terms tend to zero. In this work, we used a moving least square (MLS) interpolation scheme developped by Cueto-Felgueroso et al. for compressible flow [START_REF] Cueto-Felgueroso | Finite volume solvers and moving least-squares approximations for the compressible navierstokes equations on unstructured grids[END_REF] and applied to uncompressible flow by Ramrez et al. [START_REF] Ramrez | A new higher-order finite volume method based on moving least squares for the resolution of the incompressible navier-stokes equations on unstructured grids[END_REF]. A deferred correction is used for the calculation of convective flux. The low implicit term is a Linear Upwind Scheme, and the explicit correction is a second order MLS interpolation. The diffusive flux is directly calculated by a second order MLS interpolation. The volume integrals for the temporal term and the source term are approximated by a midpoint rule. A second order backward Euler scheme is used for time derivatives.

The linear algebra is handled through the Petsc library [START_REF] Balay | PETSc Web page[END_REF]. The three linear systems resulting from the (uncoupled) momentum equations are solved by the BiConjugate Gradient Stabilized (BiCGSTAB) method with block-Jacobi preconditioning. The pressure-correction equation is solved by Hypre's algebraic multi-grid (AMG) method.

Computational domain & mesh

The computational domain is meshed with quad cells of 1mm length. The total number of cells is 19,200. A schematic of the mesh is shown in Figure 4.

Unresolved FVM-DEM method

The numerical results obtained using the method in Section 2.1 is compared against the numerical results based on a coupled FVM-DEM method developed on the OpenFOAM platform, an opensource C++ based CFD program. In this study, the particles are smaller than the FVM mesh cell size, therefore the coupling between the finite volume method (FVM) and discrete element methodology (DEM) is achieved with an unresolved method, developed on the OpenFOAM platform.

Development of a coupled FVM-DEM

In this study, the sink term F pf accounts for gravity force, drag force, lift, and pressure gradient force, virtual mass, and Basset history forces. Brownian force is neglected as the size of particles in both numerical methods is significantly greater than 1 µm. A soft-sphere DEM method, similar to the Cundall and Strack model [START_REF] Peter | A discrete numerical model for granular assemblies[END_REF], in the form of spring-slider-dashpot model is employed to accurately resolve the trajectories of particle motion by integrating the Newtonian equations of motion. The soft-sphere approach permits one to explicitly define the properties of both particle and wall (i.e. density, Youngs modulus, Poisson ratio), and also the coefficient of restitution, coefficient of friction, and cohesion energy density between particle-particle and particle-wall interactions.The following equations are used to model the particle-particle and particle-wall interactions [START_REF] Stw Kuruneru | A coupled finite volume and discrete element method to examine particulate foulant transport in metal foam heat exchangers[END_REF].

The equations governing the motion of solid particles is given as:

m i dV ip dt = m i g + k i j=1 F C,ij + f pf,i , (19)
I i dω i dt = k i j=1 T i , (20)
where V ip is the translational velocity of a particle, and the number of particles in contact with particle i is denoted as k i , I i is the moment of intertia, rotational velocity ω i , and torque T t .

The normal contact force is given by

F cn,ij = (-k n δ 1.5 nij -η n u ij • n ij )n ij , (21)
The tangential force is given by

F ct,ij = (-k t δ tij -η t u sij ), ( 22
)
where k n is the non-linear normal spring stiffness, and k t is the tangential spring stiffness between particle contacts i and j. The following expressions connote the spring stiffness:

k n = 4 3 √ R * E 2(1 -σ 2 j ) , (23)
and

k t = 8 R * δ nij G 2(2 -σ) , (24)
where

R * = r i r j r i + r j and G = E (2(1 + σ)) , ( 25
)
where δ is the deplacement, unit vector n ij from the centres of particles i and j. The slip velocity of the contact point is given as:

u sij = u ij -(u ij • n ij )n ij + (r i ω i + r j ω j ) × n ij , (26)
where u ij is the relative velocity vector between the contact of particles i and j.

For particle sliding to occur, the following relation must be satistifed:

|F ctij | > µ f |F cnij | (27)
The resultant tangential force is expressed as a function of the friction coefficient µ f :

|F ctij | = -µ f |F cnij |u sij /|u sij |, (28)
The damping coefficient is given as

η n = λ(m * k n ) 0.5 δ 0.25 n , (29)
where m * is the effective mass and expressed as

m * = m i m j m i + m j , (30)
where m i and m j is, respectively, particles i and j. It is noteworthy that the emperical constant λ is related to the coefficient of restitution. Moreover, the damping coefficient η n is assumed to be identical to η t .

In order to accurately capture particle contacts, the DEM time-step must be around 10 -100 times smaller than the FVM time-step [START_REF] Xu | Numerical and experimental studies on transport properties of powder ejector based on double venturi effect[END_REF]; as such, the FVM time-step is set at 1 x 10 -5 s whereas the particle (discrete phase) collision resolution time-step is set to 20. The simulation is run from 0.00 s to 5.00 s and the particle injection commences at 0.20 s to allow for fluid flow development prior to the injection of particles. The particle velocity is initialized at the same velocity as the fluid velocity. The Youngs Modulus is assigned a slightly lower value than the actual to reduce computational effort. Trial numerical results show negligible difference in the particle distribution patterns irrespective of the Youngs Modulus value. The same observation is reached by Tsuiji et al. [START_REF] Tsuji | Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe[END_REF][START_REF] Tsuji | Discrete particle simulation of two-dimensional fluidized bed[END_REF].

The coupling between the two phases is achieved as follows: first, at each time-step, the DEM solver will relay the dynamic information such as positions and velocities of individual particles, in order to evaluate the porosity and the particle-fluid interaction force in a computational cell. Afterwards, the CFD solver will use this data to evaluate the gas flow field which computes the fluid forces acting on each DEM particle. Then all of these resultant forces are imported into the DEM in order to generate motion information of individual particles for the next time-step. The fluid force acting on each discrete particle will act in response on the carrier fluid from the DEM particles, thereby complying with Newtons third law of motion [START_REF] Qian | Cfd-dem simulation of the filtration performance for fibrous media based on the mimic structure[END_REF].

Numerical solution and algorithm control

A generalized (GAMG) solver and a Gauss-Seidel smoother is deployed to obtain discretized pressure equations whilst a smooth solver with a smoother symmetric Gauss-Seidel (sGS) is used to discretize the velocity equations while performing a single sweep smoothing iteration prior to re-calculating the residual in order to improve computational efficiency. These solvers operate on a LDU matrix class where the smooth solver is for symmetric and asymmetric matrices, and the smoother converges the solution to the required tolerance (or relative tolerance). The GAMG solver generates a solution on a mesh with a small number of smalls with minimal computational effort; afterwards, the solver maps the solution onto a finer mesh which uses it as a starting solution in order to generate an accurate solution on the fine mesh. This is achieved by geometrically coarsening the grid (geometric multi-grid) or directly harnessing the algebraic multi-grid irrespective of geometry. The mesh is coarsened or refined in steps. The agglomeration of cells is executed by a face area pair agglomerator. Merge levels is set to control the pace at which coarsening or refinement of the grids is performed. Typically in most situations, OpenFOAM coarsens/refines the grid one level at a time by making one cell out of four (i.e. mergeLevel 1). This level or merging generally yields optimal convergence. However, for cases with a simple mesh, coarsening (or refining) of the grids can be safely achieved at a rapid pace by coarsening (or refining) two levels at a time (i.e. mergeLevel 2).

A merged SIMPLE-PISO (PIMPLE) algorithm

The modular implementation and versatility of OpenFOAM permits one to implement a hybrid SIMPLE-PISO (PIMPLE) algorithm. This algorithm is used to couple the pressure-velocity equations for transient solutions with very large Courant numbers (1-10) which in turn assists in stabilizing the numerical convergence while preserving numerical accuracy. The PIMPLE algorithm consists of various important parameters that could be used depending on the case study. These parameters include the number of non-orthogonal correctors, number of correctors (inner loops/pressure correction), number of outer correctors (outer loops/pressure-momentum correction), momentum predictor, consistent (PIMPLEC), residual controls, under-relaxation. PIMPLEC (PISOC-SIMPLEC) is beneficial for cases with a large maximum Courant number (Co). Although the PIMPLEC algorithm is applicable for transient solutions with very large Courant number, simulations with this algorithm take longer to reach convergence compared with the SIMPLEC algorithm. This is due to the fact that the PIMPLEC algorithm involves both predictor and corrector steps. Unlike the standalone PISO algorithm, which is generally applicable for cases where Co ≤ 1, the PIMPLE algorithm permits the use of a high time-step or an adaptive time-step (i.e. assign a maximum Courant number) to numerically stabilize transient solutions which is beneficial for cases involving complex fluid flow patterns in complex geometries with skewed non-orthogonal meshes. It is also beneficial for particleladen gas flows based on unstructured meshes where the DEM particle size is on par with the FVM mesh cell size. Executing the pimple algorithm loops over the PISO algorithm in one time-step which permits some under-relaxation between these loops. This activity permits the use of larger timesteps, which is not possible in the standalone PISO algorithm. The PIMPLE algorithm attempts to solve the momentum equation one or more time at each time-step depending on the number of outer correctors assigned to the PIMPLE loop. In theory, the PIMPLE algorithm is identical to the PISO algorithm if one outer corrector is assigned. This is in theory identical to multiple PISO loops per time-step analogous to a transient SIMPLE algorithm. In other words, assigning only one outer correctors solves the momentum equation once only at each time-step which is the norm in a PISO algorithm. In the event the pressure-momentum coupling is calculated only once (one outer corrector), the standalone PISO algorithm is enforced irrespective of the number of correctors [START_REF] Balay | PETSc Web page[END_REF][START_REF] Bayomy | Electronic cooling using water flow in aluminum metal foam heat sink: Experimental and numerical approach[END_REF][START_REF] Cueto-Felgueroso | Finite volume solvers and moving least-squares approximations for the compressible navierstokes equations on unstructured grids[END_REF].

The number of correctors signifies the number of times the pressure field is corrected. For tetrahedral non-orthogonal mesh, a correction term is essential for the treatment of non-orthogonality. In other words, the number of non-orthogonal correctors corrects the solution Laplacian term of the pressure equation (surface normal gradient schemes). The value is ranges from 0 to 2; it is generally set to 0 for steady-state simulations and pure hexahedral mesh or 1 for transient and/or low-quality highly skewed meshes (i.e. max non-orthogonality angle is approximately 70 or higher). In this study, 2 correctors and 2 non-orthogonal corrector is assigned. The pressure is re-calculated based on the updated fluxes obtained from the outer loop correction. For example, in this study, Open-FOAM computes 50 SIMPLE outer loops and within one outer loop, the pressure is corrected twice. To ensure the robustness and stability of the PIMPLE algorithm, under-relaxation factors for outer iterations is enforced.

An under-relaxation factor is assigned to under relax the system of discretized equations An under-relaxation factor is assigned to under relax the system of discretized equations (eq. pp.115 [START_REF] Jasak | Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows[END_REF]); moreover, an under-relaxation factor is also used to relax the new pressure equations. In other words, solving the pressure equation again is executed in order to yield a better approximation of the correct pressure fields (eq. 3.145 pp.149 [START_REF] Jasak | Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows[END_REF]). The under-relaxation factor for the final equations and fields of velocity and pressure is set to 1 in order to comply with the conservation of mass. It is noteworthy that the PISO algorithm does not under-relax the fields and equations and the momentum corrector step is executed more than once. Likewise, the pressure-momentum correction (outer loops) is essentially a SIMPLE loop which requires under-relaxation to stabilize the solution. The number of outer correctors (outer loop correction) is set as the number of times the fluxes, pressure, and momentum are re-calculated within one time-step. However, slightly lower final residuals could be enforced only if a very high number of outer correctors (≥ 100) is assigned to find the correct solution within one time-step. Secondly, a residual sub-control is enforced for the PIMPLE algorithm in order to reduce the computational time whilst maintaining numerical stability. This allows OpenFOAM to exit the PIMPLE outer corrector loop once a solution fulfils the residual criteria during a time-step. A residual and tolerance criteria for both the PISO and PIMPLE loop/iteration is assigned, where OpenFOAM will escape the PISO loop/iteration when the final residuals within each PISO loop fall below the assigned final tolerance level. For instance, the geometric agglomerated algebraic multi-grid (GAMG) solver will iteratively solve the system of linear equations until the final residual for pressure falls below an allocated value (i.e. 1 x 10 -6 ) OpenFOAM will escape the PIMPLE loop and proceed to the next time-step if the initial residuals fall below the allocated tolerance value. The deployed PIMPLE algorithm is well suited for skewed complex geometries and meshes and multiphase transport (with a slightly high Courant number).

As the case study involves transient simulations, the relative tolerance is set to 0.0 to yield efficient PIMPLE simulations by forcing the solution to converge to the solver tolerance in each time step [START_REF]Relative tolerance openfoam userguide[END_REF]. The residuals for the pressure and velocity is set to 1 x 10 -6 and 1 x 10 -5 respectively whereas the residual control for the PIMPLE loop is assigned as 1 x 10 -3 for both pressure and velocity.

Numerical schemes

There is an array of numerical schemes available. For most practical engineering applications, the following schemes shown in Table 1 are used. Linear interpolation is widely used in a number of cases although a cubic interpolation could be deployed but it is rarely used except for very cases such as stress analysis. The corrected surface normal gradient schemes is generally used for most cases where the maximum mesh non-orthogonality does not exceed 70 • which is the case in our study.The uncorrected and orthogonal surface normal gradient schemes is normally deployed in the even the mesh exhibits very low non-orthogonality (i.e. ≤ 5 • ) The choice of the Laplacian scheme is also based on the maximum mesh non-orthogonality. A Gauss Linear corrected Laplacian scheme is deployed to obtain solution to the pressure poisson equation. Gauss Linear and Euler are deployed for the divergence and time derivatives respectively. For the numerical schemes for the particles, an Euler-implicit integration scheme is used. The coupling between the two phases permits the transfer of the corrected momentum from the discrete phase to the fluid continuum phase.

A recent study has found that the commercial software FLUENT-EDEM exhibits a miscalculation of the drag force (to a certain degree) which resulted in the overestimation of the mean particle velocity, which is attributable to the lack of an appropriate mesh interpolation scheme in their code [START_REF] Ebrahimi | Numerical and experimental study of horizontal pneumatic transportation of spherical and low-aspect-ratio cylindrical particles[END_REF].

In short, the velocity of gas is not interpolated to the particle location, and all of the solid particulates in the FVM grid encounter identical fluid velocity irrespective of the particle position within the fluid mesh cell [START_REF] Ebrahimi | Numerical and experimental study of horizontal pneumatic transportation of spherical and low-aspect-ratio cylindrical particles[END_REF]. OpenFOAM circumvents this issue by explicitly enforcing an interpolation scheme based on the carrier fluid bulk properties (i.e. fluid density, velocity, and dynamic viscosity). The reader is referred to Elghobashi (1994) [START_REF] Elghobashi | On predicting particle-laden turbulent flows[END_REF] and Xiao & Sun (2011) [START_REF] Xiao | Algorithms in a robust hybrid cfd-dem solver for particle-laden flows[END_REF] for additional details on these interpolation schemes. As the simulation is fully coupled, the momentum correction is transferred from the discrete solid phase to the fluid continuum phase. The discrete phase is coupled to the carrier phase (i.e. source terms are generated for the carrier phase via a semi-Implict scheme).

Computational domain & mesh

The FVM-DEM method is applied to the computational domain shown in Figure 1. The filter is composed of wires (300-500 µm diameter) and it forms square holes (500-800 µm). The geometry and mesh is designed in ANSYS. Afterwards, the mesh file is imported to OpenFOAM to run the simulations. Both numerical results (Section 2.1 & 2.2) are then compared against the experimental findings. The grid consists of unstructured tetrahedral mesh and the mesh is refined around the vicinity of the filter. A mesh sensitivity analysis is performed based on single-phase fluid flow. The inlet is set to 0.885 m/s. The results are shown in Table 2. The pressure drop is evaluated as the difference between the average pressure drop along the y-axis 40 mm from the inlet and the average pressure drop along the y-axis 40 mm from the rear face of the obstruction walls. The numerical simulations is executed in OpenFOAM.

Grid N umber of N odes P ressure Drop (P a) According to Table 2, less than 5 % difference is observed after every successive mesh refinement commencing Grid 3. As such, Grid 3 is used for all the simulations to ease computational burden. Open slightly vane air to allow flow in the test section 7. Wait until all the air exit the system 8. Adjust vane A to set the desired volume flow rate 9. Wait until the particles and the pressure drop stabilize 10. Record the pressure drop and take a picture of the particles

Comparison between Numerical and Experimental Findings

We compare the numerical and experimental particle distributions in three different cases as shown in Figure 2 & Figure 3. The inlet velocity is set to 0.885 m/s which corresponds to a volumetric flow rate of 4 m 3 /h. Four different solids concentration is investigated: 2 mL, 4 mL, 6 mL, 8 mL. The numerical pressure drop is evaluated and compared against the experimental values. Both the numerical and experimental pressure drop is calculated as difference between the average pressure drop along the y-axis 40 mm from the inlet and the average pressure drop along the y-axis 40 mm from the rear face of the obstruction walls at a particular time point. As time elapses, the solid particles aggregate and eventually clog the filter. Particle aggregate spread along the x-axis and y-axis is in good agreement with the experimental results. The pressure drop aligns well with the experimental observations based on a solid volume of 2 mL or 4 mL. However, there exists a large discrepancy between the experimental and numerical pressure drop at 6 mL and 8 mL. Additional simulations have shown that the solution greatly underestimates the pressure drop at a large solid volume (i.8 mL , 10 mL) irrespective of the inlet velocity. However, this is not the case for lower solid volume (i.e. 2 mL and 4 mL). The pressure drop miscalulation stems from the fact that total number of particles injected into the 2D configuration is not the same as a 3D configuration (experiments -Secion 3). Instead, a homogeneous repartition of the particle along depth of the experimental test channel is assumed and the equivalent surface of the particles in a 2D case is computed as (S s = V s /L) as specified in Section 2.1.

In 2D calculation, an artificial porosity is added to the particle in order to account the 3D effect. Although this assumption showed large discrepancy between the numerical and experimental pressure drop values in OpenFOAM, the penalization technique yield reasonably accurate numerical pressure drop results. However, it is noteworthy that although the penalization technique yields very accurate representation of the particle spread distribution profile, the penalization technique did show large discrepancy between the numerical and pressure drop distributions in several cases. This discrepancy stems from the fact that an ideal filter is used which actually blocks the particles but does not impede the incoming carrier fluid.

Secondly, the introduction of significantly fewer particles into the 2D system (i.e. 1911 particles for Case E1 compared with 96 particles for Case O1) results in a loosely packed particle aggregate bed formation which consists of more voids between the particle contacts (O1) than the 3D case (E1); moreover, the tightly packed particles for the 3D case (i.e. E1) remain motionless whereas the 2D case (i.e. O1) clearly shows that the particle undergo very faint unsteady sliding vibration-like movement as shown by the velocity contours (i.e O1) in Figure 2.

It is noteworthy that the additional validation cases (Sections 3.1.2 & 3.1.3) is based on solid-gas flows but the results in this section is solid-liquid flows. Although only a few numerical cases of solid-liquid flows align reasonably well with the experimental values, all the solid-gas flow numerical case studies are in excellent agreement with the three benchmark cases. As shown in Figure 7, at 6 mL and 8 mL solids volume, a slightly better agreement is observed between the experimental and the penalization method. However, the large discrepancy in the OpenFOAM results is linked to the 2D particle projection assumption in addition to the difference in velocities of the particles together with the compactness of the 2D and 3D cases. At 8 mL solids volume, both numerical methods greatly miscaluate the pressure drop values. Although both methods use the 2D particle injection assumption, it is interesting to note that at 6 mL and 8 mL OpenFOAM underestimates the pressure drop values, whereas the penalization method overestimates the pressure drop values, and the difference is more profound at 8 mL.

Conclusions and Perspectives

This study investigates the transient evolution of particle-laden liquid flow and particle accumulation on a filter by comparing the methodology and results of the FVM-DEM methodology developed in OpenFOAM and the the penalization technique. The results obtained by both numerical methods are then compared against experimental results.

Interestingly, many studies including this study here are based on engineering applications where granular media (i.e. particle-fluid flows) is either dense or dilute. However, certain applications may have regions of dilute or dense flows. As such, it may not be practical to use the DEM throughout the entire domain, rather, DPM could be used. In short, a FVM-DEM-DPM could be used to accurately and rapidly obtain the numerical solution (i.e. use DEM in areas of high particle (dense) concentrations and use DPM in areas of low particle (dilute) concentrations. This could be achieved by a straightforward algorithm which assigns which particle solver to use, DEM or DPM, depending on the porosity of the FVM mesh cell. The next phase of this project comprises the implementation of the energy equation to account for thermal transport between solid particles and fluid continuum in the two numerical approaches. The penalization technique is to be extended for 3D flows coupled with heat transfer and radiation effects. This topic will be investigated in the near future using OpenFOAM. This approximation method coupled with the PIMPLE algorithm will nullify numerical instability issues whilst preserving numerical accuracy which will be show in Section 4.3. Interestingly, Kuruneru et al [START_REF] Kuruneru | Cfd-dem of sandstone-laden air flow in a clear channel-preliminary investigation[END_REF] have shown that the use of both methods lead to a negligible change in particle deposition profiles and pressure drop characteristics in a 2D clear rectangular channel.

Figure 1 :

 1 Figure1: Each cell is coloured according to its solid volume fraction χ s . The difference between the particle size and the resolution of the mesh has an impact on the choice of the coupling method. A) Most particles are smaller than the mesh cell size, therefore an unresolved approach is suitable. B) The particles and the mesh cells have similar size. Neither approach is appropriate but the hybrid method is applicable. C) The flow around the largest particle can be resolved accurately.

Figure 2 :

 2 Figure 2: Clogging of a typical fuel filter [29].

Figure 3 :

 3 Figure 3: Particle position before and after a time step. In grey, the position if there is no interaction.

Figure 4 :

 4 Figure 4: Mesh

Figure 5 :

 5 Figure 5: Computational domain for OpenFOAM case

3 .

 3 Results & Discussions 3.1. Solid-Liquid Flow: Particle Accumulation and Clogging of a Filter 3.1.1. Experimental Method

Figure 6 :

 6 Figure 6: Diagram of experimental setup

Figure 7 :

 7 Figure 7: Test section dimensions

Figure 8 :

 8 Figure 8: Comparative Assessment between Penalization results (P1, P2), OpenFOAM results (O1, O2), and experimental results (E1, E2) based on three different solid concentration (P1/O1/E1: 4 mL, P2/O2/E2: 8 mL)

Figure 9 :

 9 Figure 9: Numerical (OpenFOAM) and experimental pressure drop values

Table 2 :

 2 3.95 & 3.96 OpenFOAM Numerical Schemes

	T erms

Table 3 :

 3 Mesh sensitivity analysis

Acknowledgements

The author acknowledges the computational resources provided by the Queensland University of Technologys High Performance Computing (HPC) facilities and Arts et Metiers Paris Tech. S.T.W. Kuruneru is grateful for the experimental results provided by Sofiane. S.T.W. Kuruneru acknowledges the financial support in the form of an Australian Postgraduate Award (APA) scholarship granted by the Australian government.

Conflict of Interest

None declared.