A disturbance corrected point-particle approach for two-way coupled particle-laden flows on arbitrary shaped grids
AA disturbance corrected point-particle approach for two-way coupledparticle-laden flows on arbitrary shaped grids
Pedram Pakseresht a , Sourabh V. Apte a, ∗ a School of Mechanical, Industrial and Manufacturing EngineeringOregon State UniversityCorvallis, OR 97331, USA
Abstract
A general, two-way coupled, point-particle formulation that accounts for the disturbance createdby the dispersed particles in obtaining the undisturbed fluid flow field needed for accurate com-putation of the force closure models is presented. Specifically, equations for the disturbance fieldcreated by the presence of particles are first derived based on the inter-phase momentum couplingforce in a finite-volume formulation. Solution to the disturbance field is obtained using two ap-proaches: (i) direct computation of the disturbance velocity and pressure using the reaction forcedue to particles at computational control volumes, and (ii) a linearized, approximate computationof the disturbance velocity field, specifically applicable for low Reynolds number flows. In both ap-proaches, the computed disturbance field is used to obtain the undisturbed fluid velocity necessaryto model the aerodynamic forces on the particle. The two approaches are thoroughly evaluatedfor a single particle in an unbounded and wall-bounded flow on uniform, anisotropic, as well asunstructured grids to show accurate computation of the particle motion and inter-phase coupling.The approach is straightforward and can be applied to any numerical formulation for particle-ladenflows including Euler-Lagrange as well as Euler-Euler formulations.
Keywords:
Undisturbed flow, Euler-Lagrange, Point-particle approach, Arbitrary grids.
1. Introduction
Particle-laden flows, wherein small size solid particles, liquid droplets or gaseous bubbles aredispersed in a fluid flow, are widely encountered in many engineering, biological and environmentalapplications. A wide range of numerical approaches resolving different scales of fluid and particlemotion have been developed for accurate and predictive simulation of these flows (van der Hoefet al., 2008; Balachandar and Eaton, 2010). The point-particle (PP) approach (Maxey and Ri-ley, 1983), in which particles are assumed spherical, subgrid and modeled as point sources, hasreceived much attention in modeling particle-laden flows owing to its simplicity and affordabilityin simulating motion of large number of particles ( O (10 − )). ∗ Corresponding author. 308 Rogers Hall, Corvallis, OR 97331, USA. Tel: +1 541 737 7335, Fax: +1 541 7372600.
Email addresses: [email protected] (Pedram Pakseresht), [email protected] (SourabhV. Apte)
Preprint submitted to Journal of Computational Physics December 4, 2020 a r X i v : . [ phy s i c s . f l u - dyn ] D ec his approach was originally developed for dilute particulate flows with particles smaller thanthe fluid length scale (or the computational grid size) wherein the presence of particles does not sig-nificantly perturb the characteristics of the flow, i.e., one-way coupled regime (Elghobashi, 1991).In modeling this regime, the fluid phase equations are solved irrespective of the presence of theparticles, while closures for fluid forces acting on particles such as drag, lift, added mass, pressuregradient, and history effect are employed to obtain the particle trajectories. Despite the origi-nal development of the PP approach for one-way coupled regimes, it has been widely applied tomodeling two-way coupled regimes wherein the fluid phase is indeed perturbed by the presence ofparticles (Elghobashi, 1991). Such a scenario may happen when the particle loading (or concen-tration) locally or globally becomes large, either due to a few large size particles or dense regimeof small particles. The effect of particles on the carrier phase is then modeled by applying theparticle reaction force to the fluid phase equations through a momentum source term. In addi-tion, for dense loading, the volume occupied by the particles is also removed when solving forthe fluid phase equations, by applying volume filtered equations, or volumetric coupling, whichresults in additional source terms in the continuity equation due to spatio-temporal variations inparticle volume fraction. Unlike the standard approaches, this formulation couples two phasesthrough both momentum and continuity equations (Ferrante and Elghobashi, 2004; Apte et al.,2008; Capecelatro and Desjardins, 2013; Pakseresht and Apte, 2019).Since particles in PP approach are assumed to be smaller than the grid and the local flow overthe particles is considered uniform, hence force closures based on uniform flow over a sphere aretypically employed. However, to extend the applicability range of PP approach to particles of sizeon the order of or slightly larger than the computational grid, that receive spatially varying flowfield, Faxen corrections have been developed (Annamalai and Balachandar, 2017). Typical forceclosures that are derived for a single particle rely on the undisturbed fluid flow, which is not readilyavailable in the two-way coupled simulations. By definition, the undisturbed flow is the velocityand pressure fields that would exist at the location of a particle if that particle was not there in theflow. For multi-particle systems, the undisturbed flow field seen by a particle only corresponds tothe flow field in the absence of that specific particle, however, includes the disturbed flow createdby the neighboring particles. In two-way coupled simulations, the fluid phase is altered by the self-induced disturbance of each individual particle through inter-phase momentum and mass exchange,and using the available disturbed flow field for force closures results in erroneous predictions. Theerror introduced by using the disturbed field may remain small when the particle size is muchsmaller than the grid, owing to the negligible disturbance of small particles. However, when theparticle size is on the order of or bigger than the grid resolution, such as those encountered indirect and large-eddy simulations, the disturbance due to particles becomes noticeable, hence theerrors can become significantly large.Pan and Banerjee (1996) first showed that the PP approach produces wrong prediction forvelocity of a single particle settling toward a no-slip wall, and in order to improve the predictions,they introduced a velocity-disturbance model, wherein the analytical Stokes solution at the locationof the particle is superimposed on the background flow to reflect the effect of the particle. UnlikePP approach, their model eliminates any dependency of the particle force computations to theundisturbed fluid velocity and results in more accurate predictions. Although their model can beapplied to both unbounded and wall-bounded regimes due to the availability of the Stokes solutionfor both, it is limited to small particle Reynolds numbers ( Re p ) and at steady state condition forwhich the analytical solution is available. Gualtieri et al. (2015) regularized the PP approach for2he unbounded flows by deriving analytical equations to remove the self-induced velocity distur-bance created by the particles, that is also extended to wall-bounded regimes (Battista et al., 2019).Their approach requires considerable computational resources to resolve the stencil over which theparticle force is distributed using a Gaussian filter function. Horwitz and Mani (2016, 2018) origi-nated a method to obtain the undisturbed field based on the enhanced curvature in the disturbedvelocity field for particle Reynolds numbers of Re p < .
0. A C-field library data was built usingreverse engineering technique that needs to be added to the current PP packages for recoveringthe undisturbed velocity. Their model is limited to (i) isotropic rectilinear computational grids,(ii) flows with particles of maximum size of the grid, and (iii) unbounded regimes. Ireland andDesjardins (2017) derived an analytical expression for recovering the undisturbed velocity in theunbounded regimes based on the steady state Stokes solution that was derived as the solution ofa feedback force distributed to the background flow using a Gaussian filter. Their model accountsfor the displaced fluid mass by the particles and is limited to unbounded regimes with small Re p .Using analytical and empirical expressions, Balachandar et al. (2019) developed a model thatcorrects the PP approach for cold particle-laden flows with Re p < Re p , and (iv)arbitrary interpolation and distribution functions, it is limited to unbounded flows. Paksereshtet al. (2020) extended this idea to wall-bounded flows by using empirical expressions as well aswall-modified Stokeslet solution. Their approach is applicable to large size particles and extremeanisotropic grids, typically employed in wall-bounded turbulent flows. Test cases performed onvelocity of a particle in both parallel and wall-normal motions showed that when the correctionschemes, that are developed for unbounded regimes, are employed for correcting particle force inwall-bounded regimes, errors on the same order of magnitude of the uncorrected PP scheme can beobtained. Their model is capable of recovering the undisturbed fluid velocity at any arbitrary walldistance and asymptotically approaches the regular unbounded correction schemes for particlestraveling sufficiently away from the no-slip wall. The above correction schemes by Esmaily andHorwitz (2018); Pakseresht et al. (2020) remove the self-induced disturbance for each individualparticle when correcting the particle forces, keeping the effects of all the other particles in theneighborhood unchanged. Thus, their approaches implicitly account for the effect of neighborson the individual particle force closures; however, in its present form is limited to Re p <
10 andtri-linear interpolation on rectilinear grids.Given the importance of the undisturbed field and the restrictions of the existing models, inthis work, a general formulation for estimating the disturbance field created by particles is derived.Solution to the disturbance field is obtained using two approaches: (i) direct computation ofthe disturbance velocity and pressure using the particle reaction forces at computational control3olumes, and (ii) an approximate computation of the disturbance velocity field, based on lowparticle Reynolds number assumption. Direct solution of the disturbance field is easily feasibleusing the same framework of the Navier-Stokes solver and is applicable to wall-bounded flows,complex geometries and boundary conditions, anisotropic as well as arbitrary, unstructured grids,and a wide range of Re p . The approach is straightforward and can be applied to any numericalformulation for particle-laden flows including Euler-Lagrange as well as Euler-Euler formulations.This direct solution approach does add additional computational cost, but its versatility, simplicity,and accuracy make it an attractive alternative. To reduce the computational cost, and yet keep thesame benefits mentioned above, an approximate solution, specifically applicable for low Re p , is alsopresented. Predictions from these two approaches are compared with existing uncorrected modelsfor motion of a particle in unbounded and wall bounded regimes. Furthermore, the effectivenessof these approaches for a range of grid types (structured or arbitrary shaped unstructured), gridanisotropy, and particle Reynolds numbers is evaluated.The rest of the paper is arranged as follows. Section 2 explains the existing issue in the forcecomputations of the standard two-way coupled PP simulations. The mathematical formulations forthe direct as well as approximate methods are derived in this section, as well. Section 3 validatesboth methods on a series of numerical test cases for a particle’s motion in unbounded and wall-bounded regimes using various grids and flows parameters. Section 4 concludes the paper withfinal remarks and summary of the work.
2. Mathematical Formulation
In this section, the existing issue in the force computations of the standard two-way coupled PPapproaches is first explained. Next, a general framework for correcting this issue is presented fol-lowed by a reduced order approximate, but computationally efficient method. For sake of simplicity,pressure-based incompressible fluid flow solvers are used here, however, the proposed frameworkcan be easily extended to any general flow solution techniques including those for compressibleflows. Moreover, for the simulations presented in this work, an Euler-Lagrange framework is em-ployed with that in mind that the present methods can be applied to Euler-Euler formulations, aswell.
Consider a particle-laden fluid flow as shown in Figure 1. In a typical point-particle approachin an Euler-Lagrange framework, the particles are assumed subgrid, and their motion is modeledby Newton’s second law as, dx i,p dt = u i,p , (1) du i,p dt = (cid:18) − ρ f ρ p (cid:19) g i + 1 V p ρ p F ti,p , (2)where x i,p and u i,p are particle centroid location and velocity, respectively, ρ p is the particle density, V p is the volume, and F ti,p represents the total fluid forces acting on each individual particle in the i direction. In the point-particle approach, since particles are assumed subgrid, the forces can notbe computed directly and instead are modeled using the available closures for drag, added mass,history effect (Maxey and Riley, 1983), lift force (Saffman, 1965) and Magnus effect (Rubinow andKeller, 1961), among others. 4 igure 1: Schematic diagram of flow over particles in a general domain with arbitrary boundary conditions and grids:the globally, undisturbed flow (left) and two-way coupled flow field (right). In the two-way coupling approach, the effect of the particles on the carrier fluid is modeledthrough an equal and opposite reaction force of the particles and using an appropriate interpolationkernel, modifying the fluid flow in the vicinity of the particles, ∂u i ∂t + u j ∂u i ∂x j = − ρ f ∂p∂x i + ν ∂ u i ∂x j − ρ f N p (cid:88) p =1 f ti,p ξ ∆ ( x cv − x p ) , (3) ∂u j ∂x j = 0 , (4)where ν is the fluid kinematic viscosity, ρ f is the fluid density, and f ti,p is the particle force perunit volume in the i direction. ξ ∆ ( x cv − x p ) is a kernel function to project the particles forces,that lie within the bandwidth (∆), to the computational cell center at x cv . N p is the total numberof particles that are located within the bandwidth of the projection function. The choice of thisprojection function is dependent upon the flow under consideration and accuracy needed.Typical force closure models used in computing the particle motion are based on the fluid flowfield that is undisturbed by the presence of the particle. As an example, the steady state Stokesdrag force over a particle with diameter of d p moving with velocity of u i,p in a fluid with dynamicviscosity of µ is F Stokesi,drag = 3 πµd p (cid:0) u uni, @ p − u i,p (cid:1) , (5)wherein u uni, @ p is the undisturbed fluid velocity at the location of the particle that is not influencedby the presence of the particle under consideration—that is, without the particle self-induceddisturbance. However, this undisturbed flow field is not readily available in a two-way coupledsimulation as the self-induced disturbance in the fluid flow created by the reaction force of theparticle alters the flow velocity and pressure fields. It should be noted that, the undisturbed fluidflow needed in the closure models for the motion of a particle refers to the velocity and pressurefields in the absence of that particular particle, however, accounts for the disturbance effect createdby any of the neighboring particles.In the present work, a general formulation is developed to compute the undisturbed flow fieldthat removes the disturbance created by all particles and is denoted as u ugi and p ug , a globally5ndisturbed velocity and pressure fields. This approach then allows formulating equations for thedisturbance field created by all particles in an Eulerian frame. For a particle under consideration,the undisturbed flow field at the particle location can then be written as, u uni, @ p = u ugi, @ p + (cid:88) nbr δu nbri, @ p , (6)where δu nbri, @ p is the velocity perturbation created by a neighboring particle at the location of theparticle p and in its absence. This neighboring effect can be substantial for regimes where inter-particle distance is comparable to the particle size. The mean inter-particle distance varies as φ − / , where φ is the local particle volume fraction. For example, with φ = 0 .
01, the nearestneighbor distance is about 3.7 times the particle diameter (Akiki et al., 2017). For systems withdilute to moderate volume loadings of φ ≤ − , the hydrodynamic inter-particle interactions areinsignificant, and thus the effect of neighboring particles on the undisturbed flow field is negligible.In this regime, u ugi ∼ u uni and p ug ∼ p un . In the present work, emphasis is placed on recovering theglobal, undisturbed flow field ( u ugi , p ug ). The formulation is thus directly applicable to dilute load-ings, and can be applied to moderate-to-dense loadings by explicitly incorporating the neighboringparticle effects in the future (Moore et al., 2019; Seyed-Ahmadi and Wachs, 2020).To recover the global undisturbed flow field, a general framework is proposed wherein governingequations for the disturbance created by all particles in the flow are first formulated. In the absenceof any particles, the fluid flow equations can be written as, ∂u ugi ∂t + u ugj ∂u ugi ∂x j = − ρ f ∂p ug ∂x i + ν ∂ u ugi ∂x j , (7) ∂u ugj ∂x j = 0 . (8)The velocity ( u di ) and pressure ( p d ) disturbance fields created by all the particles can be obtained bysubtracting the disturbed two-way coupled flow field, expressed by Eqs. 3-4, from the undisturbedflow field, given by Eqs. 7-8, ∂u di ∂t + u dj ∂u ugi ∂x j + u j ∂u di ∂x j = − ρ f ∂p d ∂x i + ν ∂ u di ∂x j + 1 ρ f N p (cid:88) p =1 f ti,p ξ ∆ ( x cv − x p ) (9) ∂u dj ∂x j = 0 , (10)where, u ugi = u i + u di , (11) p ug = p + p d . (12)Here, u i and p represent the standard two-way coupled velocity and pressure fields, u di and p d arethe disturbance fields created by all particles. The above equation has the boundary conditionof u di =0 far away from the particles. If the particle is near a no-slip wall, the disturbance fieldalso experiences the same no-slip condition of u di =0, making it a general formulation for any flowconfiguration, computational approach, and grid type. Solution to the above equations can beobtained by using two different approaches as described below.6 .2. Direct Method In order to solve the equations 9 and 10, the second, nonlinear term on the left hand side ofthe momentum equation, u dj ∂u ugi /∂x j , requires additional closure between the undisturbed anddisturbance fields. However, this term can be safely neglected in comparison to the third term byhypothesizing the following, ∂u ugi ∂x j (cid:28) ∂u di ∂x j , (13)which is approximately valid as the velocity gradient caused by the particle force in the disturbancefield is conjectured to be greater than that of the undisturbed field. Although this assumption isverified by the small errors for the studied cases reported in section 3, further investigations forcases with inherently large velocity gradient in the undisturbed field is left for future investigations.Knowing the particle forces, the disturbance field can then be obtained by directly solving thefollowing equations, ∂u di ∂t + u j ∂u di ∂x j = − ρ f ∂p d ∂x i + ν ∂ u di ∂x j + 1 ρ f N p (cid:88) p =1 f ti,p ξ ∆ ( x cv − x p ) , (14) ∂u dj ∂x j = 0 . (15)Note that the nonlinear, advective term contains the disturbance velocity ( u di ) and the two-way,coupled velocity ( u j ). The latter is readily available in a two-way coupled simulation.The same computational algorithm employed for the solution of the main two-way coupled flowfield (Eqs. 3-4) can be utilized to compute the disturbance field (equations above), which involvessolution of the disturbance momentum equations, and projection of the divergence-free disturbancecondition using a solution of a Poisson equation for the disturbance pressure. Direct solution ofthe disturbance field is then used to recover the undisturbed flow field from Eqs 11-12 and toaccurately compute the fluid forces acting on the particles. Compared to the existing correctionschemes, the direct method benefits from many advantages as explained below: • Direct solution is easily feasible using the same framework employed for solving the Navier-Stokes equations of the two-way coupled field. • Wall-bounded flows, complex geometries, and arbitrary boundary conditions can be auto-matically accounted for in the solution of the disturbance field. • Unlike the majority of the existing models, the direct method is capable of handling a widerange of Re p . • The disturbance velocity and pressure fields (and thus their gradients) are available, thus thecommon force closures for drag, lift, history effect, as well as pressure gradient forces can beaccurately computed. • The formulation is free of any tuning or empirical expressions, typically used for specificgrid configurations, and can be applied to both structured and arbitrary shaped, hybridunstructured grids with any grid aspect ratio.7
The formulation is free of any dependency on the interpolation and projection functionsemployed in the two-way coupled simulations. • The disturbance field is computed regardless of size of particles, hence it is adaptable forflows with any arbitrary size particles, particularly those with particles larger than grid.The main drawback of this approach is the additional computational cost for solving the dis-turbance field that requires full solution of the momentum as well as continuity equations, whichmakes the computations as nearly twice as expensive. The additional cost can still be tolerable fordirect numerical or large eddy simulations, as the approach is much more affordable than particle-resolved methods. However, for a faster computation, an approximate method is introduced in thefollowing part, which is shown to be reasonably accurate as compared to the direct method whilebeing significantly more cost efficient.
In this part, an approximate solution of the disturbance field is proposed that is derived basedon low particle Reynolds number assumption. In the limit of creeping flow ( Re p < . ∂u di ∂t = − ρ f ∂p d ∂x i + ν ∂ u di ∂x j + 1 ρ f N p (cid:88) p =1 f ti,p ξ ∆ ( x cv − x p ) , (16) ∂u dj ∂x j = 0 . (17)In order to further simplify the equations above, we conjecture that the fluid response to theparticle force, is approximately analogous to the flow that would be generated by the particle inthe real physics of the problem. This is in fact the main assumption employed in the two-waycoupled point-particle approach wherein it is assumed that the particle force can approximatelyproduce the same flow as the particle would do in the reality. Motivated by this analogy and inthe limit of steady state and Re p < .
1, we recall the Stokes solution that is the flow created aroundthe actual particle. In the Stokes regime, the drag on the particle that experiences slip velocity of u relp , consists of two terms (i) pressure and (ii) viscous forces as, F Stokesdrag = πµd p u relp (cid:124) (cid:123)(cid:122) (cid:125) Pressure force + 2 πµd p u relp (cid:124) (cid:123)(cid:122) (cid:125) Viscous force . (18)For low particle Reynolds numbers, the expression for these two forces are identical, with viscousforce being twice greater than the pressure force. Motivated by this, one can model the contributionof the pressure drag through an effective viscosity and rewrite the Stokes drag force as, F Stokes i, drag = 2 πµ eff d p u relp ; µ eff = K ν µ ; K ν = 1 . . (19)Rewriting the Stokes drag in this form facilitates the approximation that the effect of thepressure gradient term in Eq.16 can be modeled through an equivalent viscous term with aneffective viscosity of K ν =1 . K ν =1 .
5. The approximateequation then becomes, ∂u di ∂t = K ν ν ∂ u di ∂x j + 1 ρ f N p (cid:88) p =1 f ti,p ξ ∆ ( x cv − x p ) . (20)It is worth mentioning that the correction factor K ν can be Reynolds number dependent. Withincrease in Re p , the contribution of the pressure drag to the net drag is bound to increase (White,2006), and the value of K ν can potentially be changed. For the present study, however, the valueis kept fixed and equal to 1.5, even for higher Re p , and further adjustment for larger Re p casesare left for future. Similar to the direct method, the equation above is solved in the same Eulerianframe that is used for solving the two-way coupled flow field equations. This captures the resultantdisturbance field that is caused by all particles in the flow field. For dilute loadings, whereinthe disturbance of each particle is isolated from that of the neighbouring particles, this modelperfectly captures the self-induced disturbance of each particle required in force closures. However,for dense loadings, the neighbouring effect will be removed and such an effect should be addedexplicitly using the recently developed models by Moore et al. (2019); Seyed-Ahmadi and Wachs(2020). Depending on the application of interest, both unbounded and wall-bounded regimes canbe handled by this method since the boundary conditions are directly enforced for solving theequation above. Finally, owing to the linear, unsteady diffusion like equation with a source term,its solution is considerably faster than the direct method. Note that Eq. 20 is general and directlyapplicable to any arbitrary grid. Concerning the applicability of this method for Re p > .
1, it isshown later (section 3) that despite the fact that this method is constructed upon the assumptionof small Re p , it can reduce the errors for a wider range of particle Reynolds numbers of Re p ≤ The procedure in the present disturbance-corrected point-particle (DCPP) framework is simi-lar to the standard uncorrected point-particle approach with an additional step for recovering theundisturbed flow field. For the computations of the present methods, two sets of parameters andequations, corresponding to the disturbance as well as two-way coupled disturbed flow fields, aresolved separately yet on similar computational domains and identical boundary conditions. Notethat any interpolation and projection functions that are used for computations of the two-waycoupled flow field should be used for the computation of the disturbance field, to ensure that thepredicted disturbance is consistent with the one that particles actually sample in the disturbedtwo-way coupled flow field. Knowing the computed disturbance field, u di , and particles velocity, u i,p , from the previous time step, the following procedure is employed.1. Solve Eqs. 3 and 4 for the two-way coupled field to update the fluid velocity, u i , and pressure, p ,due to presence of particles. Note that this is the standard step in the uncorrected PP approaches.2. Knowing the disturbance field available from previous time step, and the updated disturbed9ow field from step 1, recover the undisturbed fluid velocity, u ugi , and pressure, p ug , fields at thelocation of particles using Eqs. 11 and 12.3. Use the undisturbed field to compute the net fluid forces acting on each particle, F ti,p .4. Update the velocity and location of each individual particle using Eqs. 1 and 2.5. Knowing the particles reaction forces from step 3, compute the disturbance field by solvingeither the direct method (Eqs. 7 and 8) or the approximate method (Eq.20).Present work is based on an energy-conserving scheme for unstructured, arbitrarily shapedgrid elements based on fractional time-stepping on a colocated mesh (Mahesh et al., 2004). Thevelocity and pressure are stored at the centroids of the volumes. The cell-centered velocitiesare advanced in a predictor step, the predicted velocities are interpolated to the faces and thenprojected. Projection yields the pressure potential at the cell-centers, and its gradient is usedto correct the cell and face-normal velocities, using an area weighted least-squares minimizationtechnique (Mahesh et al., 2004). Details of the algorithm on arbitrary shaped unstructured gridsfor particle-laden flows are given in Shams et al. (2011) and a brief description is presented inAppendix A for completeness. The same algorithm is used for the disturbance field in the directand approximate methods.
3. Results
In this section, the performance of the direct as well as approximate methods on recoveringthe undisturbed flow field and correcting the PP approach is verified in various scenarios. A singleparticle settling under gravity in an unbounded regime is investigated first. Settling velocity ofthe particle moving parallel and normal to a no-slip wall is performed next. As the final testcase, the unsteady motion of a single particle in an oscillatory field is examined, as well. Forsimplicity, drag force as the only fluid force acting on the particle is employed while other fluidforces such as lift, added mass, pressure gradient, and history effect are neglected. For each set,various grid configurations including isotropic and anisotropic rectilinear grids as well as tetrahedralunstructured grid are used to assess the accuracy of the present models on arbitrary shaped grids.A range of 0 . ≤ Re p ≤
100 is performed to evaluate the models for a wide range of scenarios thatmay happen in particle-laden flows. The grid resolution of 128 was chosen for all cases (with closeproximity to this resolution for the unstructured grid) as it was found to be sufficient to producethe grid-independent results.Three non-dimensional parameters are used to setup the cases: (i) particle-to-gird size ratio, Λ,(ii) particle Stokes number, St , and (iii) particle Reynolds number, Re p . The first dimensionlessparameter, Λ, is defined as Λ = d p d c , (21)where d p and d c are the particle diameter and the characteristic length of the grid, respectively.For rectilinear grids, d c can become a vector with three components each of which correspondingto the size of the grid in that direction, a i , hence three components for Λ, as well. However, forunstructured grids, finding an equivalent grid size for each direction is not trivial, therefore, a uni-fied d c based on the diameter of a sphere that has the equivalent volume of the grid, d c = (cid:112) V cv /π ,is defined. Particle Stokes number is defined as, St = τ p τ f , (22)10here, τ p = ρ p d p µ , (23)and, τ f = min ( d c ) ν , (24)are the respective particle relaxation time and fluid time scale in the Stokesian regime. The particleReynolds number in this regime is also defined as, Re Stkp = | u Stks | d p ν , (25)where, u Stks = (cid:18) − ρ f ρ p (cid:19) τ p g , (26)is the particle settling velocity under gravity vector of g . It is imperative to mention that for thestudied cases with Re Stkp > . Re p , differs from Eq. 25. For each of those cases, the proper expression isprovided, separately.The fluid velocity at the particle’s location, required for the drag force computation, is inter-polated using a three-point delta function with compact support that uses the nearest neighborsof a control volumes (Roma et al., 1999). For control volumes with resolution of ∆= √V cv , theinterpolation stencil utilizes only three points in one dimension and thus is easiest to implement: ξ ∆ ( x cv − x p ) = (5 − | r | − (cid:112) − − | r | ) + 1) , . ≤ | r | ≤ . , r = | x cv − x p | / ∆ (1 + √− r + 1) , | r | ≤ . . otherwise (27)Given the force balance acting over the particle, it is advanced using a first order Euler ap-proximation to solve Eqs. 1 and 2. Concerning the two-way coupled simulations, the particlereaction force is exerted to the nearby fluid control volumes using identical function as expressedabove. Although a simple three-point delta function is used for Euler-Lagrange interpolation andprojection purposes, the present methods can be easily adopted for any other functions. For thecomputations, we correct the PP approach using both direct and approximate methods and com-pare their results to those of the uncorrected PP as well as the corresponding reference for eachpart. The reference is obtained based on the one-way coupled simulations wherein the fluid phaseremains undisturbed, and drag force and particle motion are accurately computed. In this part, settling velocity of a single particle in an unbounded quiescent fluid is performed.A particle that is initially at rest settles under a gravity vector and experiences drag force, only.The particle equation of motion then becomes, du i,p dt = (cid:18) − ρ f ρ p (cid:19) g i − fτ p ( u i,p − u ugi, @ p ) , (28)11here u ugi, @ p is the interpolated fluid velocity at the location of the particle that is erroneouslynonzero in the uncorrected two-way coupled simulations, owing to the disturbance created by theparticle in the nearby computational cells. The direct as well as approximate methods, however,predict and remove this velocity as in the real physics of the problem (and one-way coupled sim-ulations) this velocity is zero. In general, the factor of f can correspond to any adjustment tothe Stokes drag to account for different effects. In this part, it follows the Schiller-Naumman ad-justment factor (Clift et al., 2005), as expressed below, to account for the finite Reynolds numbereffect of the particle on the Stokes drag in unbounded regime, f = 1 + 0 . Re . p ; Re p = | u p | d p ν . (29)Such an adjustment results in an effective particle relaxation time, τ ep , as, τ ep = τ p f . (30)Following Horwitz and Mani (2016), gravity vector of g =[1 , (1 + √ / , exp(1)] / | g | is chosen sothat particle sweeps through different positions among its adjacent computational cells, ensuringthat the present models are capable of handling any arbitrary positioning of the particle. The timestep, ∆ t , for the computations is also chosen as,∆ t = min (cid:0) . τ f , . τ ep , . τ cfl (cid:1) , (31)with τ cfl being the time scale associated with the fluid advection in high particle Reynolds numbercases based on Courant-Friedrichs-Lewy (CFL) condition less than one for time-accurate solutions.Here, a maximum CFL of 0 . u p ( t ), is decomposed into twocomponents (i) parallel, u || p , and (ii) normal, u ⊥ p to the reference settling velocity (terminal velocity)of u s , and are obtained respectively as, u || p ( t ) = u s · u p ( t ) | u s | u s , (32)and u ⊥ p ( t ) = u p ( t ) − u || p ( t ) . (33)The errors in these two velocity components as well as the total particle velocity in a time averagemanner, denoted by (), are then calculated using the respective following metrics, e (cid:107) = u (cid:107) p ( t ) . u s | u s | −
1; (34) e ⊥ = | u ⊥ p ( t ) || u s | ; (35) e = | u p ( t ) − u s || u s | . (36)12able 1 lists these errors for a particle with Re p =0 . d c over the entire grid. It isobserved that the errors associated with the uncorrected PP approach depend on Λ with biggerparticles producing stronger disturbances in the background flow. As an example, particle incase U2 produces five times larger errors compared to that of case U1 that has a particle fivetimes smaller. Similar comparison is observed between cases U8 and U9. Concerning the effect ofparticle Stokes number, results based on two different St =3 . St =10 (e.g., case U1 and U3,respectively) show small dependency of the errors on this parameter, consistent with the precedingworks. When the standard PP approach is corrected using the present methods, however, significanterror reduction is observed with nearly zero errors for the direct method across the board. Althoughthe approximate model yields slightly larger errors compared to the direct model, such as those inU2 and U9, the overall errors are still an order of magnitude lower than the uncorrected scheme.The affordability of the approximate method makes it an attractive scheme for recovering theundisturbed flow field, given the fact that the difference in the error reduction between these twomethods is still insignificant. The errors for approximate method in this case are also comparableto those reported by Pakseresht et al. (2020). Figure 2 qualitatively illustrates the performance ofthese models in predicting the time-dependent particle velocity for different grid resolutions. Theparticle relaxation time and settling velocity expressed by Eqs. 23 and 26, respectively, are usedfor normalizing the results.The performance of each method in capturing the effect of the particle on the fluid phase isillustrated in Fig.3 that pertains to cases U2 and U9 from Tab.1 with Λ=5 .
0. Contours of fluidvelocity magnitude normalized by the particle settling velocity of each case are shown at the timeinstance of 2 τ p , from the initial release of the particle. The uncorrected scheme is compared againstthe corrected results using both direct as well as approximate methods. For the uncorrected scheme,the fluid phase experiences smaller velocity compared to the corrected results, owing to the smallerslip velocity that the particle samples in this scheme, resulting in a smaller drag force exertedto the flow. Concerning the predictions of the corrected schemes, both direct and approximatemethods show nearly identical results in capturing the particle’s effect on the background flow.The observations here imply stronger inter-phase coupling when the PP approach is corrected,which could potentially yield more accurate predictions for two-way coupled particle-laden flows.The predictive capability of the present models for a wide range of particle Reynolds number,typically encountered in various applications, is investigated next. Table 2 provides settling inthe range of 1 100 and using the same three grid configurations employed earlier. The firstobservation from the results here is that the errors in the uncorrected scheme decreases as particleReynolds number increases, in line with the preceding works (Balachandar et al., 2019), suggestingthat the need for correction schemes becomes less important for Re p > . 67% reduces to 11 . 16% when Re p increases from 1 to 100 on the isotropic rectilinear grid. Thisis justified due to the fact that higher Re p particles move faster, and the residence time in theirown disturbance field, created in the previous time step, becomes smaller than that of the slowerparticles, hence the lower disturbance. Moreover, Balachandar et al. (2019) showed that as Re p Figure 2: Temporal evolution of particle velocity settling in an unbounded, quiescent flow at Re p =0 . St =10, isshown. Plotted are four different cases from Tab.1: (a) isotropic rectilinear grid (case U1), (b) anisotropic rectilineargrid with Λ =[4 . , . , . 0] (case U5), (c) anisotropic rectilinear grid with Λ =[0 . , . , . 6] (case U6) and (d) tetrahedralunstructured grid with Λ=1 . τ p , and settling velocity expressed by Eqs. 23 and 26, respectively, are used for normalizing theresults. XZ Y XZ Y XZ Y XZY XZ Y XZ Y XZ Figure 3: Contours of fluid velocity magnitude normalized by the particle settling speed at time instance of 2 τ p :particle and grid resolution (first column), uncorrected (second column), direct method correction (third column),and approximate method correction (fourth column). Results are based on case U2 (top row) and U9 (bottom row)from Tab.1 with Λ = 5 . ell shape case St Λ uncorrected e (cid:107) e ⊥ e corrected usingdirect method e (cid:107) e ⊥ e corrected usingapproximate method e (cid:107) e ⊥ e U1 10.0 [1.0,1.0,1.0] 54.2 0.077 54.2 -0.056 0.0006 0.056 0.95 0.0689 0.95U2 10.0 [5.0,5.0,5.0] 276.5 0.21 276.5 -0.28 0.007 0.28 5.0 0.54 5.04U3 3.0 [1.0,1.0,1.0] 52.13 0.06 52.13 -0.008 0.0003 0.008 0.08 0.1 0.17U4 10.0 [5.0,0.5,0.5] 35.42 2.86 35.54 0.0003 0.0003 0.0004 0.76 2.64 2.75U5 10.0 [4.0,2.0,0.2] 30.26 4.56 30.61 0.0003 0.0005 0.0007 1.83 4.31 4.7U6 10.0 [0.3,6.0,0.6] 21.22 3.53 21.51 0.0004 0.0002 0.0004 -2.31 3.42 4.13U7 3.0 [0.3,6.0,0.6] 14.52 3.82 15.02 0.0006 0.0001 0.0006 -2.28 3.85 4.48U8 10.0 1.0 25.55 0.44 25.56 -0.06 0.06 0.11 -3.01 0.48 3.06U9 10.0 5.0 129.66 2.17 129.68 0.65 0.39 0.81 -17.75 2.77 18.03U10 3.0 1.0 25.39 0.60 25.40 0.10 0.11 0.16 -3.63 0.62 3.70 Table 1: Listed are the percentage error in particle settling velocity e || , drifting velocity e ⊥ , and total velocity e , ofa particle in an unbounded, quiescent flow. Results compare the uncorrected scheme to the direct and approximatecorrection methods for Re p =0 . 1. Computational grids with different shapes and particle-to-gird size ratios are shown. increases, the region of maximum disturbance travels farther downstream so that the disturbedfluid velocity sampled at the particle location will be smaller for larger Re p . Nevertheless, whenthe PP approach gets corrected by either methods, more accurate predictions are achieved for thestudied range of Re p . It should be noted that even though the correction may not be necessary forlarge Re p cases (e.g., see error of 1 . 35 for uncorrected case of UR6), depending upon the particle-to-grid size ratio and grid anisotropy, the errors in uncorrected settling velocity could still be onthe order of 10% (see case UR3). For such cases, both methods are still effective in reducing theerrors of the uncorrected scheme by an order of magnitude. Such a capability of the present modelsfor reducing errors even for large Re p cases makes them more general to be employed without anyrestriction for a specific range of application. Figure 4 illustrates the time dependent velocity ofa single particle settling at different Re p predicted by the present correction schemes compared tothe uncorrected PP approach and the reference. Settling velocity of u s =(1 − ρ f /ρ p ) τ p | g | /f andparticle time scale of τ p are used for normalizing the results. In this part, the capability of the present models for wall-bounded regimes is evaluated. Fortest cases here, an additional non-dimensional parameter that is the normalized wall distance fromthe bottom of the particle is defined as, δ p = x ,p d p − . , (37)wherein x ,p is the wall-normal distance from the center of the particle (wall is assumed in x – x plane). Concerning complex geometries, computing this distance to the nearest wall might not bestraightforward and would have to be investigated in the future.Settling velocity of a particle at various wall distances is carried out using the present methodsand on the three aforementioned computational grids. For the studied cases, a particle that is16 Figure 4: Temporal evolution of particle settling velocity in an unbounded, quiescent flow at (a): Re p =0 . 1, (b): Re p =1 . 0, (c): Re p =10 . Re p =100 . 0. Results of the approximate method (solid blue), direct method (dashedred), uncorrected scheme (dash-dotted green) are compared against the reference velocity (dotted black). Resultsare based on cases U1 and UR1-UR3 from Tab. 1 and 2, respectively. The settling velocity of u s =(1 − ρ f /ρ p ) τ p | g | /f and partilce relaxation time of τ p are used for normalizing the results. ell shape case Re p Λ uncorrected e (cid:107) e ⊥ e Corrected usingDirect method e (cid:107) e ⊥ e corrected usingapproximate method e (cid:107) e ⊥ e UR1 1.0 [1.0,1.0,1.0] 44.67 0.12 44.67 0.16 0.007 0.16 2.18 0.03 2.18UR2 10.0 [1.0,1.0,1.0] 19.9 0.58 19.9 2.08 0.03 2.08 4.26 0.66 4.31UR3 100.0 [1.0,1.0,1.0] 11.15 0.21 11.16 2.82 0.03 2.82 3.44 0.24 3.45UR4 1.0 [0.3,6.0,0.6] 20.32 3.79 20.67 -0.003 0.001 0.003 -1.77 3.92 4.31UR5 10.0 [0.3,6.0,0.6] 7.53 1.47 7.67 -0.02 0.01 0.03 -0.27 1.62 1.65UR6 100.0 [0.3,6.0,0.6] 1.34 0.22 1.35 -0.56 0.02 0.56 -0.64 0.22 0.68UR7 1.0 1.0 15.32 0.47 15.33 -0.07 0.03 0.1 -3.07 0.42 3.10UR8 10.0 1.0 7.17 0.15 7.17 0.48 0.03 0.48 0.52 0.17 0.56UR9 100.0 1.0 2.33 0.03 2.33 0.49 0.01 0.49 0.59 0.03 0.59 Table 2: Listed are the effect of Re p on the errors in settling, drifting and total velocity of a particle in an unboundeddomain and predicted by direct, approximate and uncorrected schemes compared to the reference. Results are basedon isotropic rectilinear grid, anisotropic rectilinear grid as well as a tetrahedral unstructured grid. St = 10 . initially stationary and located at a given δ p , released to reach its settling velocity under a gravityvector of g =[exp(1) , , (1 + √ / / | g | that guarantees the particle’s motion on a plane parallelto the wall. In reality, the particle experiences a lateral force (Vasseur and Cox, 1977; Takemuraand Magnaudet, 2003) that is neglected in this study to isolate the parallel motion. The particle’sequation of motion in the presence of wall still follows Eq. 28 with the adjustment factor of f thataccounts for the wall effects on the particle’s drag coefficient. Concerning this factor, the empiricalexpression derived by Zeng et al. (2009) is employed that covers a wide range of δ p and Re p as, f || ( δ p , Re p ) = f || ( δ p ) f || ( δ p , Re p ) , (38)where, f || ( δ p ) = (cid:20) . − . 071 + 4 δ p − 815 log (cid:18) δ p 135 + 256 δ p (cid:19)(cid:21) ; (39) f || ( δ p , Re p ) = (cid:20) . (cid:16) − exp (cid:16) − (cid:112) δ p (cid:17)(cid:17) Re ( . . 313 exp ( − √ δ p )) p (cid:21) . (40) f || ( δ p ) captures wall effects on the Stokes drag for zero Re p , that approaches unity when δ p →∞ .The second term, f || ( δ p , Re p ), however, handles the wall-modified finite Reynolds number effecton the Stokes drag coefficient that converts to the Schiller-Naumman adjustment factor (Eq. 29)when particle travels sufficiently away from the wall.Using the correction factor expressed above, the particle’s equation of motion (Eq.28) is solvedusing the corrected and uncorrected PP approaches. Following the metrics presented in the pre-ceding subsection, the errors in settling, drifting and total velocity of the particle are measured incomparison to those of the one-way coupled simulations that serves as the reference. Table 3 showsthese errors for settling velocity of a particle with Re p =0 . St =10 on (i) isotropic rectilinear18 ell shape case δ p Λ uncorrected e (cid:107) e ⊥ e corrected usingdirect method e (cid:107) e ⊥ e corrected usingapproximate method e (cid:107) e ⊥ e A1 0.05 [1.0,1.0,1.0] 71.32 0.03 71.32 -1.38 0.005 1.38 3.99 0.09 3.99A2 0.5 [1.0,1.0,1.0] 42.53 0.02 42.53 -0.21 0.0008 0.21 0.89 0.04 0.9A3 1.0 [1.0,1.0,1.0] 45.81 0.02 45.81 -0.05 0.0005 0.05 0.27 0.04 0.27A4 2.0 [1.0,1.0,1.0] 49.73 0.02 49.73 -0.13 0.0006 0.13 -0.02 0.04 0.05A5 ∞ [1.0,1.0,1.0] 54.44 0.04 54.44 0.03 0.0008 0.06 1.14 0.05 1.15B1 0.05 [0.3,6.0,0.6] 29.28 1.62 29.33 0.08 0.00 0.08 7.6 1.62 7.77B2 0.5 [0.3,6.0,0.6] 24.18 1.71 24.24 0.14 0.00 0.14 5.64 1.72 5.9B3 1.0 [0.3,6.0,0.6] 26.73 2.26 26.82 0.27 0.00 0.27 5.32 2.28 5.79B4 2.0 [0.3,6.0,0.6] 28.8 2.61 28.91 0.53 0.001 0.53 4.44 2.64 5.17B5 ∞ [0.3,6.0,0.6] 28.53 2.37 28.63 1.25 0.02 1.25 3.59 2.46 4.36C1 0.05 1.0 30.28 0.50 30.28 0.53 0.29 0.64 -1.47 0.59 1.68C2 0.5 1.0 19.13 0.34 19.14 0.18 0.14 0.25 -1.21 0.33 1.28C3 1.0 1.0 18.60 0.29 18.60 -0.06 0.08 0.13 -2.17 0.29 2.21C4 2.0 1.0 21.63 0.38 21.63 -0.32 0.06 0.33 -2.75 0.45 2.8C5 ∞ Table 3: Percentage errors in the settling velocity, e || , drifting velocity, e ⊥ , and total velocity, e , of a single particlein parallel motion to a no-slip wall on different grids for St =10 and Re p =0 . δ p . grid (set A), (ii) anisotropic rectilinear grid (set B), and (iii) tetrahedral unstructured grid (setC). Each grid has six cases corresponding to settling at different δ p , that covers a wide range ofdistances from the wall.The first observation from Tab. 3 is that the uncorrected scheme produces significantly largeerrors in predicting particle velocity at all wall distances, with slightly larger values near the wall,consistent with the observations of Pakseresht et al. (2020). It is imperative to mention that thereported errors here in this work are slightly smaller than those of Pakseresht et al. (2020), poten-tially due to different Euler-Lagrange interpolation schemes employed in this study as comparedto their tri-linear interpolation. When the direct and approximate correction methods used toobtain undisturbed fluid velocity, the errors reduce by one or two orders of magnitude compared tothe uncorrected scheme. Of specific interest is the results obtained from the approximate methodthat shows same order of accuracy as those reported by Pakseresht et al. (2020) for rectilineargrids. Different grid configurations used in the present work, including unstructured grids, showthe applicability of the present methods for more complex geometries that are encountered in thereal world applications.In order to test the present methods for higher Re p in wall-bounded regimes, settling velocityof a particle near a no-slip wall is computed for 1 . 05, for which the errors in the settling velocity were observed to be more remarkablein the preceding part. Consistent with the unbounded cases (see Tab. 2), as Re p increases theerror in the particle velocity decreases and the need for correcting the PP approach diminishes. Asan example, for the particle settling on the unstructured grid, the total error of 25 . 95 for Re p =1 . ell shape case Re p Λ uncorrected e (cid:107) e ⊥ e corrected usingdirect method e (cid:107) e ⊥ e corrected usingapproximate method e (cid:107) e ⊥ e WR1 1.0 [1.0,1.0,1.0] 65.21 0.11 65.21 -0.92 0.005 0.92 3.86 0.01 3.86WR2 10.0 [1.0,1.0,1.0] 31.92 0.86 31.93 4.29 0.09 4.29 8.06 0.91 8.11WR3 100.0 [1.0,1.0,1.0] 16.48 0.37 16.48 10.72 0.24 10.72 11.08 0.4 11.09WR4 1.0 [0.3,6.0,0.6] 28.21 1.58 28.25 0.04 0.001 0.04 7.01 1.60 7.19WR5 10.0 [0.3,6.0,0.6] 19.28 0.84 19.30 0.11 0.004 0.11 3.86 0.95 3.98WR6 100.0 [0.3,6.0,0.6] 9.98 0.1 9.98 0.03 0.007 0.03 1.23 0.13 1.23WR7 1.0 1.0 25.94 0.6 25.95 0.43 0.38 0.59 -2.68 0.7 2.8WR8 10.0 1.0 12.09 0.18 12.09 0.21 0.05 0.22 -0.4 0.22 0.51WR9 100.0 1.0 3.72 0.04 3.72 -0.1 0.003 0.1 -0.04 0.03 0.06 Table 4: Effect of Re p on the particle settling velocity at δ p =0 . 05 predicted by the present models in comparisonwith uncorrected scheme showing errors in settling, e || , drifting, e ⊥ and total, e , velocities. (case WR7) decreases to 3 . 72 when Re p increases to 100 (case WR9). Nevertheless, both directand approximate models are able to reduce the errors even for large Re p by an order of magnitude.Not shown here, similar results were obtained for settling at other wall distances with large Re p . This section tests the ability of the present methods for recovering the undisturbed fluid velocityfor particles in wall-normal motion. Freely falling particle toward a no-slip wall is studied. In thisconfiguration, the drag coefficient of the particle increases as it approaches the wall, owing to thewall lubrication effect. Accordingly, for the wall adjustment factor to the particle’s drag coefficientof this part, the asymptotic expressions derived by Brenner (1961); Cox and Brenner (1967) as, f ⊥ ( δ p ) = (cid:16) . δ p (cid:17) for δ p > . 38 (Brenner, 1961) δ p (cid:16) . δ p log (cid:16) δ p (cid:17) + 1 . δ p (cid:17) for δ p < . 38 (Cox and Brenner, 1967), (41)are employed that include two parts depending on the wall normal distance of the particle. Thisadjustment factor is used in Eq. 28 to compute the particle’s equation of motion.Following the work of Pakseresht et al. (2020), a particle that is initially stationary and locatedat an arbitrary δ p =7, falls under gravity vector of g =[0 , − , . 100 and two particleStokes number of St = 3 . St = 10 . 0. For the computations of this part, the particle Reynoldsnumber is defined based on the unbounded Stokes regime, Re Stkp , expressed by Eq. 25. It shouldbe noted that the drag expression provided by Eq. 41 is only valid for Re p < . 1, however, it isstill employed for larger Re p cases just for numerical demonstration without advocating its use for Re p > . 1. The error for each method is measured based on the total time that the particle requiresto reach δ p =0 . t ref , that is obtained based on the one-way20 Figure 5: Normalized wall distances (left) and velocity (right) of a particle in moving normal to a no-slip wall.Results pertain to case N2 from Tab. 5 comparing direct as well as approximate methods to the uncorrected andthe reference solutions. The settling velocity, u s , based on Eq. 26 and the time scale of d p /u s are employed fornormalizing the results. coupled simulation. The deviation of each scheme from the reference is quantified based on therelative error as, e n = t − t ref t ref . (42)Table 5 shows that for the studied cases, the uncorrected scheme yields negative errors revealingthe fact that the particle in this scheme experiences smaller drag force, accelerates faster andreaches the wall distance of interest much quicker than it would in reality. When the PP approachis corrected using the direct method, however, small errors of O (0 . 1) are achieved that showsthe successful predictions of this method for the studied range of flow parameters and the gridconfigurations. For the approximate method, the errors are reduced to smaller values as well,however, for the highly skewed anisotropic rectilinear grid, this method results in errors on thesame order of magnitude of the uncorrected scheme. This is attributed to (i) the response ofa fluid to a source in a control volume, may differ based on the shape of the control volumeowing to the numerical approach used in solving the governing equation, (ii) even with anisotropicgrids, the distribution of the particle reaction force is done to the nearest neighbors of the controlvolume which could be asymmetric with high aspect ratio grids, and (iii) for particle motionnormal to a no-slip wall, the pressure distribution on the particle surface is asymmetric, and thusa simple approach to model the pressure gradient in the approximate method potentially needsto be modified. Concerning the effect of Re p , it is observed that for the studied computationalgrids, as Re p increases the error in the uncorrected scheme decreases, consistent with the previousobservations.Figure 5 shows qualitatively the prediction of the different methods on particle velocity andtrajectory of case N2 from Tab. 5. The settling velocity, u s , based on Eq. 26 and the time scaleof d p /u s are employed for normalizing the results. As illustrated, the direct method captures thetrajectory and velocity of the particle quite well in addition to the promising prediction of theapproximate method compared to the uncorrected approach.21 ell shape case Re Stkp St Λ uncorrected e n corrected usingdirect method e n corrected usingapproximate method e n N1 0.1 3.0 [1.0,1.0,1.0] -28.12 0.40 9.28N2 0.1 10 [1.0,1.0,1.0] -22.71 0.12 6.52N3 1.0 10 [1.0,1.0,1.0] -6.56 -0.28 0.56N4 10.0 10 [1.0,1.0,1.0] -1.29 -0.07 0.10N5 100.0 10 [1.0,1.0,1.0] -0.2 0.00 0.04N6 0.1 3.0 [0.3,6.0,0.6] -8.66 -0.96 20.08N7 0.1 10 [0.3,6.0,0.6] -8.65 -0.96 20.05N8 1.0 10 [0.3,6.0,0.6] -5.92 -0.26 14.28N9 10.0 10 [0.3,6.0,0.6] -0.74 0.00 1.36N10 100.0 10 [0.3,6.0,0.6] -0.04 0.00 0.09N11 0.1 3.0 1.0 -13.99 0.53 9.82N12 0.1 10 1.0 -12.04 0.37 7.90N13 1.0 10 1.0 -3.09 0.07 1.33N14 10.0 10 1.0 -0.54 -0.02 0.15N15 100.0 10 1.0 -0.07 0.00 0.02 Table 5: Percentage errors in the time that a particle initially located at δ p =7 . δ p =0 . Re Stkp and St with and without correction schemes. .4. Particle in oscillatory field As a final test case, the models are validated for unsteady motion of a single particle in anoscillatory flow field. A sinusoidal function is prescribed as a body force acting on the particle inan arbitrary direction of i , to resemble its unsteady motion in an oscillatory field. The followingequation of motion is used for the studied cases of this part, du i,p dt = sin( ωt ) − fτ p ( u i,p − u ugi, @ p ) , (43)with f being calculated based on Eq. 29 and ω being the frequency of the oscillation. The amplitudeof oscillation is set to be unity for the sake of brevity. Similar to the previous sections, u ugi, @ p is theinterpolated fluid velocity at the particle’s location, that is incorrectly nonzero in the uncorrectedtwo-way coupled simulations. For the reference, and consistent with other subsections, one-waycoupled results are used wherein u ugi, @ p =0. An additional non-dimensional parameter, Strouhalnumber, is also defined as, Str = τ ep τ w , (44)which expresses the ratio of the particle time scale (Eq. 30) and the time period of the oscillation, τ w =1 /ω . Table 6 lists the studied cases with the grid configurations employed in the previoussubsections as well as various flow parameters. Defining a constant particle Reynolds number forthis part might not be trivial due to the variation in the particle’s velocity. Therefore, in order to setup the cases of this part, one can choose the maximum particle Reynolds number, Re maxp , definedbased on the maximum particle’s velocity that occurs at the first crest point of its velocity profile.For the studied cases here, we perform two Strouhal number of Str =0 . Re maxp =0 . 097 and 99 . 87. Since finding an analytical expression for Re maxp might not be straightforward and setting up cases depends on this parameter as well, weprovide the dimensional parameters that are needed for reproducing the reported cases here. Thetime step for the computations of this part, ∆ t osc , follows the expression below which requires anadditional constraint to the ∆ t provided by Eq. 31, to accurately resolve the oscillation time scale,as well. ∆ t osc = min (∆ t, . τ w ) (45)Figure 6 shows excellent predictions of the present methods in capturing the time-dependentvelocity of the particle in the unsteady field with Re maxp =0 . 097 and different Strouhal numbersand using various grid configurations. As illustrated, the uncorrected scheme overshoots the crestand troughs of the particle velocity with significant deviation for small Strouhal number cases(left column of the figure), consistent with the observations of Horwitz and Mani (2016). Theperformance of the present models for higher particle Reynolds numbers of Re maxp =99 . 87, is shownin Fig. 7, signifying the capability of the present correction methods even for unsteady motions,as well.The test cases used in this study, underscore the applicability of the present models for a widerange of applications. Depending upon the accuracy necessary for simulation of particle-laden flows,either method can be chosen with the caveat that the direct method is more computationallyexpensive, but yields most accurate results. Concerning more sophisticated scenarios such asparticles close to two different walls, curved walls, corners or rough walls, both methods are still23 20 40 60 80-1.5-1-0.500.511.5 (a) referenceuncorrecteddirect methodapproximate method Figure 6: Time-dependent velocity [ m/s ] of a particle in an oscillatory motion is shown for the direct and approximatemethods in comparison to the standard uncorrected scheme and the reference solution. Cases (a) to (f) pertain toO1, O2, O5, O6, O7 and O8, respectively, from Tab. 6, that are all based on Re maxp =0 . Str =0 . Str =1 . 0, respectively. ell shape case Re maxp Str Λ ω ( s − ) µ ( N.s/m ) ρ p ( kg/m ) d p ( m ) ∆ t osc ( s )O1 0.097 0.1 [1.0,1.0,1.0] 0.1028 9.9700 180.0 1.0 0.0030O2 0.097 1.0 [1.0,1.0,1.0] 0.8968 8.7000 180.0 1.0 0.0034O3 99.87 0.1 [1.0,1.0,1.0] 0.0067 0.1480 180.0 1.0 0.0400O4 99.82 1.0 [1.0,1.0,1.0] 0.0588 0.1292 180.0 1.0 0.0400O5 0.097 0.1 [0.3,6.0,0.6] 1.1257 0.2730 5.0 0.3 0.0002O6 0.097 1.0 [0.3,6.0,0.6] 9.8136 0.2380 5.0 0.3 0.0003O7 0.097 0.1 1.0 0.0844 17.99 180.0 1.482 0.0037O8 0.097 1.0 1.0 0.7363 15.69 180.0 1.482 0.0042 Table 6: Listed are the cases performed for the unsteady motion of a particle in an oscillatory field. Two Str andvarious grid configurations are performed for Re maxp =0 . Re maxp =99 . 87 are investigated, as well. ρ f =1 . kg/m ) and St =10 . 0, are shared among all cases. referenceuncorrecteddirect methodapproximate method Figure 7: Time-dependent particle velocity [ m/s ] with Re maxp =99 . 87 in the oscillatory field predicted by the presentmethods and uncorrected scheme in comparison to the reference. Results in (a) and (b) pertain to cases O3 and O4from Tab. 6 with Str =0 . Str =1 . 0, respectively. 4. Conclusions A general, disturbance-corrected, point-particle (DCPP) formulation for two-way coupled com-putations of particle-laden flows is developed that recovers the undisturbed flow field necessaryfor accurate computation of the fluid forces acting on the dispersed particles The formulation isapplicable to arbitrary shaped grids, both structured or unstructured, and in complex configura-tions involving no slip walls. To this end, governing equations for the disturbance created by theparticle forces on the fluid flow are derived using the two-way coupled equations. Since the two-way coupled formulation and the disturbance flow equations are derived for any general boundaryconditions, the developed approach is applicable to any complex flow with or without no-slip walls.The formulation can be implemented in any fluid flow solver and is not limited to specific types ofgrids.Two models are developed to compute the disturbance field: (i) a direct method, and (ii)an approximate method. In the direct method, the non-linear disturbance momentum equationstogether with the continuity equation are solved using the same numerical formulation as the fluidflow solver for the two-way coupled field. This direct method provides the disturbance velocity andpressure fields, is free of any empirical or calibrated expressions, and makes it attractive for a widerange of computations including complex geometries, arbitrary shaped unstructured grids, as wellas particle-laden flows with any arbitrary particle size and particle Reynolds number. However,the cost associated with this approach is nearly doubled, as the disturbance momentum togetherwith the continuity constraint require additional Poisson solution for the disturbance pressure.Nevertheless, the accuracy gained by such computation warrants its use, and the cost is stillconsiderably lower than particle-resolved, direct numerical solutions wherein the grid resolutionsused are much finer than the particle size.In order to alleviate the computational cost associated with the direct method, a reduced order,approximate method was introduced, wherein a simplified momentum equation for the disturbancefield is solved. This approximate model is based on low Reynolds number assumption and neglectsthe non-linear, advective terms. In addition, in the steady, Stokes flow limit, the Stokes solutionover a spherical particle motivates an approximation for the pressure gradient term. The pressureand viscous terms are modeled by a modified viscous term with an effective increased viscositydetermined to match pressure and viscous forces in the Stokes flow limit. Since, the pressure fieldis not directly computed, the continuity constraint is only indirectly imposed, and the expensivestep of solving a Poisson equation for the disturbance field is not needed. Although the approximatemodel was constructed based on the assumption of small Re p , where inertial effects are negligible,the test cases for high particle Reynolds numbers, up to Re p =100, remarkably show good predictivecapability of this approach. In addition, the accuracy of this approximate method can be furtherimproved by making the effective viscosity a function of the particle Reynolds number, however,for majority of the cases studied, this was not necessary. The approximate method is shown to be26s accurate as the direct method and has considerable reduced computational cost that is on thesame order of the existing correction schemes.Both models were tested for various scenarios using isotropic and anisotropic rectilinear grids,tetrahedral unstructured grid, different particle sizes, a wide range of particle Reynolds numbers(0 ≤ Re p ≤ Re p with the fact that the need for correction diminishes whenparticle settling Reynolds number increases (Balachandar et al., 2019; Pakseresht et al., 2020).Prediction of particle settling near a no-slip wall was also evaluated for a single particle inparallel and normal motion to the wall. It was shown that both the direct and approximatemethods were capable of recovering the undisturbed field and reduced the errors to small valuesfor particle settling parallel to a no-slip wall. For particle motion normal to a wall, the direct methodshowed excellent prediction in recovering the undisturbed field and produced correct trajectory andvelocity of the particle. The approximate method also produced small errors for nearly isotropicgrids; even for particles larger than the grid resolution. However, for highly skewed anisotropicgrids, its overprediction in the undisturbed field yielded errors on the same order of magnitude asthe uncorrected scheme. An interpolation stencil that scales with the particle size, may alleviatethis issue. In addition, for particle motion normal to a no-slip wall, the pressure distribution on theparticle surface is asymmetric, and thus a simple approach to model the pressure gradient used inthe approximate method potentially needs to be modified. Nevertheless, the approximate methodis capable of capturing motion of a particle near a wall accurately, especially for nearly isotropicand arbitrary shaped grids.Finally, to test the models for unsteady motion of particles, an inevitable criterion in complexparticle-turbulence interaction, particle in oscillatory motion was investigated varying the Strouhalnumber (0 . . . . 0. Excellent predictions were achieved usingboth methods revealing their predictive capability even in unsteady motion.The present DCPP approach can be easily implemented in Euler-Lagrange and Euler-Eulerpackages as it leverages the identical algorithm, boundary conditions, and type of the computationalgird that are employed for solving the standard two-way coupled PP approaches. In its currentform, the DCPP approach is directly applicable to systems with dilute volume loading, whereinthe particle-particle hydrodynamic interactions are negligible. The approach can also be appliedto dense regimes, by explicitly modeling the hydrodynamic interaction and neighboring particleeffects. 5. Acknowledgements Financial support was provided under the NASA Contract Number NNX16AB07A monitoredby program manager Dr. Jeff Moder, NASA Glenn Research Center as well as the NationalScience Foundation (NSF) under Grant Number 1851389. In addition, the authors acknowledgethe Texas Advanced Computing Center (TACC) at the University of Texas at Austin for providingHPC resources that have contributed to the results reported here. The authors acknowledge Mr.Shashank Karra for generating the unstructured grids for various configurations used in the presentwork. 27 ppendix A. Numerical formulation The numerical approach used is based on fractional time-stepping on colocated, arbitraryshaped, unstructured, grid elements for constant density, incompressible flows. A semi-implicitscheme is used for the momentum equation solution, however, the inter-phase momentum exchangeterms are treated explicitly. The collocated grid arrangement is used for its easy application tostructured as well as arbitrary unstructured grids. (a) (b) cv1cv2 face Figure 8: Stencil used for the colocated grid, fractional time-stepping based algorithm for the two-way coupled flowfield ( u i , p ) and the disturbance field ( u di , p d ). Figure 8 shows the schematic of variable storage for fluid and particle phases. All variables arestored at the control volume (cv) center with the exception of the face-normal velocity u N (and u d N ), located at the face centers. The face-normal velocity is used to enforce continuity equation.Subscript ‘p’ is used to denote the disperse phase. Using these variable locations, integrating thegoverning equations over the control volume and applying Gauss’ divergence theorem to convertvolume integrals to surface integrals wherever possible, the discrete governing equations are derived.Accordingly, the continuity equation is1 V cv (cid:88) faces of cv u n +1N A face = 0 , (46)where ∆ t is the flow solver time-step, V cv is the volume of the cv, A face is the area of the face ofa cv and u N is the face-normal velocity. For the present colocated grid finite volume scheme theface-normal velocity u N is obtained through a projection scheme rather than interpolation of thecontrol volume based velocity to the faces. The discrete momentum equation for the i componentof velocity can be written as g n +1 i, cv − g ni, cv ∆ t + 1 V cv (cid:88) faces of cv g n +1 / i, face u n +1 / A face = − ∂∂x i p n +1cv +1 V cv (cid:88) faces of cv ( τ ij ) n +1 / N j, face A face + f n +1 i, cv , (47)where g i = ρ f u i represents the momentum per unit volume in the i direction, ρ f is the constantfluid density, ( τ ij ) face is the viscous stress at the faces of control volume, and N j, face represents thecomponents of the outward face-normal. Similarly, the velocity field ( u i, face ), and the momentum28 i, face = ρ f u i, face at the faces are obtained using arithmetic averages of the corresponding fieldsat the two control volumes associated with the face. The values at time level t n +1 / are obtainedby simple time-averaging. The interface coupling force is represented by f i, cv . The pressure field p n +1cv is unknown and is obtained using the best available guess at the current iteration. This getsupdated during the solution of the pressure Poisson equation. The above discretization is implicitand thus the time-steps are not limited by viscous stability limits. The use of symmetric centereddifferences makes the algorithm second order on uniform Cartesian grids. The main steps of thesolver are described below. • Step 1: Set the flow velocity at t n +1 using a second-order Adams-Bashforth predictor.Advance the particle positions and velocities using the undisturbed fluid velocity obtainedfrom the solution of the two-way coupled fluid flow equations ( u i and p ) and the disturbancevelocity and pressure fields ( u di and p d ) from the previous time step. • Step 2: Advance the two-way coupled fluid momentum equations using the fractional stepalgorithm, with the interphase force, f i , treated explicitly. ρ f u ∗ i − u ni ∆ t + ρ f V cv (cid:88) faces of cv (cid:2) u ∗ i, face + u ni, face (cid:3) u n +1 / A face = − δpδx i n + µ V cv (cid:88) faces of cv (cid:18) ∂u ∗ i, face ∂x j + ∂u nj, face ∂x i (cid:19) A face + f n +1 i , (48)where N is the face-normal component, and A face is the face area, µ is the fluid viscosity,and ρ f the density. The pressure gradient at the CV centers in the above equation is atthe old time-level and is obtained as described below in Step 6. The reaction force f n +1 i isobtained through Eulerian-Lagrangian interpolation (equation 27). In the above step, theviscous terms are treated implicitly, the three equations for the velocity components at theCV centers are solved using iterative scheme such as Gauss-Seidel or algebraic multigridsolvers. • Step 3: Remove the old pressure gradient to obtain the velocity field, (cid:98) u i : (cid:98) u i − u ∗ i ∆ t = + 1 ρ f δpδx i n (49) • Step 4: Interpolate the velocity fields to the faces of the control volumes and consider thecorrector step: ρ f u n +1N − (cid:98) u N ∆ t = − δpδx N n +1 , (50)where (cid:98) u N = (cid:98) u i, face N i, face is the approximation for face-normal velocity and N i, face are thecomponents of the face-normal. The face-based velocity is simply obtained as the average ofthe two control volumes that share the common face, (cid:98) u i, face = ( (cid:98) u i, cv1 + (cid:98) u i, cv2 ) / δp n +1 δx N = p n +1nbr − p n +1cv | S cv → nbr | , (51)29here the subscripts ‘cv’ and ‘nbr’ stand for the the control volume for which the velocity fieldis being solved and the neighboring control volumes sharing a common face, respectively and | S cv → nbr | represents the magnitude of the position vector connecting the two control volumes. • Step 5: The pressure field and the pressure gradients at t n +1 are unknown in the above step.A pressure Poisson equation is derived by taking a discrete divergence of the above equationsand solving for the pressure field at each control volume: (cid:88) face of cv δp n +1 δx N = ρ f ∆ t (cid:88) faces of cv (cid:98) u i, face A face (52) • Step 6: Reconstruct the pressure gradient at the cv-centers. The face-normal pressuregradient δp/δx N and the gradient in pressure at the cv-centers are related by the area-weighted least-squares interpolation (Mahesh et al., 2004): (cid:15) cv = (cid:88) faces of cv (cid:0) P (cid:48) i, cv N i, face − P (cid:48) face (cid:1) A face , (53)where P (cid:48) i, cv = δpδx i and P (cid:48) face = δpδx N . • Step 7: Compute new face-based velocities, and update the cv-velocities: u n +1N = (cid:98) u N − ∆ tρ f δp n +1 δx N (54) u n +1 i, cv = (cid:98) u i, cv − ∆ tρ f δp n +1 δx i, cv (55) • Step 8: Repeat steps 2–7 for the disturbance field, u di and p d , solving the momentum andcontinuity equations expressed by Eqs. 14 and 15 for the direct method. For the approximatemethod, repeat step 2 solving the approximate disturbance equation provided by Eq. 16. • Step 9: Using the disturbance field, compute the undisturbed flow velocity (direct andapproximate method) and pressure (direct method) and proceed to Step 1. References Akiki, G., Jackson, T., Balachandar, S., 2017. Pairwise interaction extended point-particle model for a random arrayof monodisperse spheres. Journal of Fluid Mechanics 813, 882–928.Annamalai, S., Balachandar, S., 2017. Fax´en form of time-domain force on a sphere in unsteady spatially varyingviscous compressible flows. Journal of Fluid Mechanics 816, 381–411.Apte, S., Mahesh, K., Lundgren, T., 2008. Accounting for finite-size effects in simulations of disperse particle-ladenflows. International Journal of Multiphase Flow 34 (3), 260–271.Balachandar, S., Eaton, J. K., 2010. Turbulent dispersed multiphase flow. Annual review of fluid mechanics 42,111–133.Balachandar, S., Liu, K., Lakhote, M., 2019. Self-induced velocity correction for improved drag estimation in Euler–Lagrange point-particle simulations. Journal of Computational Physics 376, 160–185.Battista, F., Mollicone, J.-P., Gualtieri, P., Messina, R., Casciola, C. M., 2019. Exact regularised point particle(erpp) method for particle-laden wall-bounded flows in the two-way coupling regime. Journal of Fluid Mechanics878, 420–444. renner, H., 1961. The slow motion of a sphere through a viscous fluid towards a plane surface. Chemical engineeringscience 16 (3-4), 242–251.Capecelatro, J., Desjardins, O., 2013. An Euler–Lagrange strategy for simulating particle-laden flows. Journal ofComputational Physics 238, 1–31.Clift, R., Grace, J. R., Weber, M. E., 2005. Bubbles, drops, and particles. Courier Corporation.Cox, R. G., Brenner, H., 1967. The slow motion of a sphere through a viscous fluid towards a plane surface—ii smallgap widths, including inertial effects. Chemical Engineering Science 22 (12), 1753–1777.Elghobashi, S., 1991. Particle-laden turbulent flows: direct simulation and closure models. Applied Scientific Research48 (3-4), 301–314.Esmaily, M., Horwitz, J., 2018. A correction scheme for two-way coupled point-particle simulations on anisotropicgrids. Journal of Computational Physics 375, 960–982.Evrard, F., Denner, F., van Wachem, B., 2020. Euler–Lagrange modelling of dilute particle-laden flows with arbitraryparticle-size to mesh-spacing ratio. Journal of Computational Physics: X, 100078.Ferrante, A., Elghobashi, S., 2004. On the physical mechanisms of drag reduction in a spatially developing turbulentboundary layer laden with microbubbles. Journal of Fluid Mechanics 503, 345–355.Gualtieri, P., Picano, F., Sardina, G., Casciola, C. M., 2015. Exact regularized point particle method for multiphaseflows in the two-way coupling regime. Journal of Fluid Mechanics 773, 520–561.Horwitz, J., Mani, A., 2016. Accurate calculation of stokes drag for point–particle tracking in two-way coupled flows.Journal of Computational Physics 318, 85–109.Horwitz, J., Mani, A., 2018. Correction scheme for point-particle models applied to a nonlinear drag law in simulationsof particle-fluid interaction. International Journal of Multiphase Flow 101, 74–84.Ireland, P. J., Desjardins, O., 2017. Improving particle drag predictions in Euler–Lagrange simulations with two-waycoupling. Journal of Computational Physics 338, 405–430.Liu, K., Lakhote, M., Balachandar, S., 2019. Self-induced temperature correction for inter-phase heat transfer inEuler–Lagrange point-particle simulation. Journal of Computational Physics 396, 596–615.Mahesh, K., Constantinescu, G., Moin, P., 2004. A numerical method for large-eddy simulation in complex geometries.Journal of Computational Physics 197 (1), 215–240.Maxey, M. R., Riley, J. J., 1983. Equation of motion for a small rigid sphere in a nonuniform flow. The Physics ofFluids 26 (4), 883–889.Moore, W., Balachandar, S., Akiki, G., 2019. A hybrid point-particle force model that combines physical and data-driven approaches. Journal of Computational Physics 385, 187–208.Pakseresht, P., Apte, S. V., 2019. Volumetric displacement effects in Euler–Lagrange LES of particle-laden jet flows.International Journal of Multiphase Flow 113, 16–32.Pakseresht, P., Esmaily, M., Apte, S. V., 2020. A correction scheme for wall-bounded two-way coupled point-particlesimulations. Journal of Computational Physics 420, 109711.Pan, Y., Banerjee, S., 1996. Numerical simulation of particle interactions with wall turbulence. Physics of Fluids8 (10), 2733–2755.Roma, A. M., Peskin, C. S., Berger, M. J., 1999. An adaptive version of the immersed boundary method. Journal ofcomputational physics 153 (2), 509–534.Rubinow, S., Keller, J. B., 1961. The transverse force on a spinning sphere moving in a viscous fluid. Journal of FluidMechanics 11 (03), 447–459.Saffman, P., 1965. The lift on a small sphere in a slow shear flow. Journal of Fluid Mechanics 22 (02), 385–400.Seyed-Ahmadi, A., Wachs, A., 2020. Microstructure-informed probability-driven point-particle model for hydrody-namic forces and torques in particle-laden flows. Journal of Fluid Mechanics 900.Shams, E., Finn, J., Apte, S., 2011. A numerical scheme for Euler–Lagrange simulation of bubbly flows in complexsystems. International Journal for Numerical Methods in Fluids 67 (12), 1865–1898.Takemura, F., Magnaudet, J., 2003. The transverse force on clean and contaminated bubbles rising near a verticalwall at moderate reynolds number. Journal of Fluid Mechanics 495, 235–253.van der Hoef, M. A., van Sint Annaland, M., Deen, N., Kuipers, J., 2008. Numerical simulation of dense gas-solidfluidized beds: a multiscale modeling strategy. Annu. Rev. Fluid Mech. 40, 47–70.Vasseur, P., Cox, R., 1977. The lateral migration of spherical particles sedimenting in a stagnant bounded fluid.Journal of Fluid Mechanics 80 (3), 561–591.White, F. M., 2006. Viscous fluid flow. Vol. 3. McGraw-Hill New York.Zeng, L., Najjar, F., Balachandar, S., Fischer, P., 2009. Forces on a finite-sized particle located close to a wall in alinear shear flow. Physics of fluids 21 (3), 033302.renner, H., 1961. The slow motion of a sphere through a viscous fluid towards a plane surface. Chemical engineeringscience 16 (3-4), 242–251.Capecelatro, J., Desjardins, O., 2013. An Euler–Lagrange strategy for simulating particle-laden flows. Journal ofComputational Physics 238, 1–31.Clift, R., Grace, J. R., Weber, M. E., 2005. Bubbles, drops, and particles. Courier Corporation.Cox, R. G., Brenner, H., 1967. The slow motion of a sphere through a viscous fluid towards a plane surface—ii smallgap widths, including inertial effects. Chemical Engineering Science 22 (12), 1753–1777.Elghobashi, S., 1991. Particle-laden turbulent flows: direct simulation and closure models. Applied Scientific Research48 (3-4), 301–314.Esmaily, M., Horwitz, J., 2018. A correction scheme for two-way coupled point-particle simulations on anisotropicgrids. Journal of Computational Physics 375, 960–982.Evrard, F., Denner, F., van Wachem, B., 2020. Euler–Lagrange modelling of dilute particle-laden flows with arbitraryparticle-size to mesh-spacing ratio. Journal of Computational Physics: X, 100078.Ferrante, A., Elghobashi, S., 2004. On the physical mechanisms of drag reduction in a spatially developing turbulentboundary layer laden with microbubbles. Journal of Fluid Mechanics 503, 345–355.Gualtieri, P., Picano, F., Sardina, G., Casciola, C. M., 2015. Exact regularized point particle method for multiphaseflows in the two-way coupling regime. Journal of Fluid Mechanics 773, 520–561.Horwitz, J., Mani, A., 2016. Accurate calculation of stokes drag for point–particle tracking in two-way coupled flows.Journal of Computational Physics 318, 85–109.Horwitz, J., Mani, A., 2018. Correction scheme for point-particle models applied to a nonlinear drag law in simulationsof particle-fluid interaction. International Journal of Multiphase Flow 101, 74–84.Ireland, P. J., Desjardins, O., 2017. Improving particle drag predictions in Euler–Lagrange simulations with two-waycoupling. Journal of Computational Physics 338, 405–430.Liu, K., Lakhote, M., Balachandar, S., 2019. Self-induced temperature correction for inter-phase heat transfer inEuler–Lagrange point-particle simulation. Journal of Computational Physics 396, 596–615.Mahesh, K., Constantinescu, G., Moin, P., 2004. A numerical method for large-eddy simulation in complex geometries.Journal of Computational Physics 197 (1), 215–240.Maxey, M. R., Riley, J. J., 1983. Equation of motion for a small rigid sphere in a nonuniform flow. The Physics ofFluids 26 (4), 883–889.Moore, W., Balachandar, S., Akiki, G., 2019. A hybrid point-particle force model that combines physical and data-driven approaches. Journal of Computational Physics 385, 187–208.Pakseresht, P., Apte, S. V., 2019. Volumetric displacement effects in Euler–Lagrange LES of particle-laden jet flows.International Journal of Multiphase Flow 113, 16–32.Pakseresht, P., Esmaily, M., Apte, S. V., 2020. A correction scheme for wall-bounded two-way coupled point-particlesimulations. Journal of Computational Physics 420, 109711.Pan, Y., Banerjee, S., 1996. Numerical simulation of particle interactions with wall turbulence. Physics of Fluids8 (10), 2733–2755.Roma, A. M., Peskin, C. S., Berger, M. J., 1999. An adaptive version of the immersed boundary method. Journal ofcomputational physics 153 (2), 509–534.Rubinow, S., Keller, J. B., 1961. The transverse force on a spinning sphere moving in a viscous fluid. Journal of FluidMechanics 11 (03), 447–459.Saffman, P., 1965. The lift on a small sphere in a slow shear flow. Journal of Fluid Mechanics 22 (02), 385–400.Seyed-Ahmadi, A., Wachs, A., 2020. Microstructure-informed probability-driven point-particle model for hydrody-namic forces and torques in particle-laden flows. Journal of Fluid Mechanics 900.Shams, E., Finn, J., Apte, S., 2011. A numerical scheme for Euler–Lagrange simulation of bubbly flows in complexsystems. International Journal for Numerical Methods in Fluids 67 (12), 1865–1898.Takemura, F., Magnaudet, J., 2003. The transverse force on clean and contaminated bubbles rising near a verticalwall at moderate reynolds number. Journal of Fluid Mechanics 495, 235–253.van der Hoef, M. A., van Sint Annaland, M., Deen, N., Kuipers, J., 2008. Numerical simulation of dense gas-solidfluidized beds: a multiscale modeling strategy. Annu. Rev. Fluid Mech. 40, 47–70.Vasseur, P., Cox, R., 1977. The lateral migration of spherical particles sedimenting in a stagnant bounded fluid.Journal of Fluid Mechanics 80 (3), 561–591.White, F. M., 2006. Viscous fluid flow. Vol. 3. McGraw-Hill New York.Zeng, L., Najjar, F., Balachandar, S., Fischer, P., 2009. Forces on a finite-sized particle located close to a wall in alinear shear flow. Physics of fluids 21 (3), 033302.