E Abriak

Johannes Joubert
email: johannes.joubert@imt-lille-douai.fr

Daniel N Wilke

Nicolin Govender

Patrick Pizette

Ugur Tuzun

Nor-Edine Abriak

3D gradient corrected SPH for fully resolved particle-fluid interactions

Keywords: Smoothed Particle Hydrodynamics (SPH), Solid-Fluid Coupling, Weakly Compressible, Gradient Correction, Particle Drag, Boundary Condition

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Processing of granular material is a very common problem in many different industries such as construction, pharmaceuticals or even food production and storage. This makes understanding the behaviour of these materials very important, especially from a design and optimization point of view. The discrete element method (DEM), as introduced in 1978 by Cundall and Strack [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], is a numerical tool used to describe the behaviour of a granular media by resolving interactions between individual particles. With its simplified contact laws, DEM allows itself to be cast into a massively parallel form which, by exploiting developments in multiprocessing, enables DEM to be applied to industrial scale problems [START_REF] Govender | Development of a convex polyhedral discrete element simulation framework for NVIDIA Kepler based GPUs[END_REF]. Unfortunately, granular materials are rarely found in isolation in their operating environments, with other physics (usually at a continuum scale) not only influencing the granular material, but also being of interest to the designer. As an example, moisture in geotechnical soil and pharmaceutical powders as well as particle-liquid interactions in food processing, chemical catalysis and comminution processes plays a large role in the overall system performance. In a more generalized sense, a few key examples include fluid-solid interaction, thermofluid-solid interaction and solid particle fracture. Due to coupling usually taking place at the particle-continuum interface, it is clear that the boundary condition implementations would have a significant effect on the quality of results obtained. Due to the parallel nature of DEM, it would be desirable to couple the DEM solver to other parallelizable solvers such as the lattice Boltzmann method (LBM) [START_REF] Ding | Study on the multiphase fluid-solid interaction in granular materials based on an LBM-DEM coupled method[END_REF][START_REF] Tran | Numerical modelling of backward front propagation in piping erosion by DEM-LBM coupling[END_REF] or smoothed particle hydrodynamics (SPH) [START_REF] Tang | A Lagrangian-based SPH-DEM model for fluid-solid interaction with free surface flow in two dimensions[END_REF][START_REF] He | A GPU-based coupled SPH-DEM method for particle-fluid flow with free surfaces[END_REF][START_REF] Cleary | Prediction of coupled particle and fluid flows using DEM and SPH[END_REF].

In recent years, smoothed particle hydrodynamics (SPH) has shown itself to be a unique and powerful technique to simulate many dynamic fluid systems. Due to its particle representation, it has been found that SPH very naturally handles multi-fluid and multi-physics environments with various groups having had success simulating multi-phase flow [START_REF] Colagrossi | Numerical simulation of interfacial flows by smoothed particle hydrodynamics[END_REF][START_REF] Chen | An SPH model for multiphase flows with complex interfaces and large density differences[END_REF][START_REF] Natsui | An SPH Study of Molten Matte-Slag Dispersion[END_REF], multi-species interactions [START_REF] Cummins | An SPH Projection Method[END_REF][START_REF] Hu | An incompressible multi-phase SPH method[END_REF][START_REF] Szewc | Modeling of natural convection with Smoothed Particle Hydrodynamics: Non-Boussinesq formulation[END_REF] and physical phenomena such as temperature driven buoyancy flows [START_REF] Szewc | Modeling of natural convection with Smoothed Particle Hydrodynamics: Non-Boussinesq formulation[END_REF][START_REF] Leroy | Buoyancy modelling with incompressible SPH for laminar and turbulent flows[END_REF].

However, due to the meshless Lagrangian formulation of SPH, boundary conditions have posed several problems in the past. Libersky and Petschek [START_REF] Libersky | Smooth particle hydrodynamics with strength of materials[END_REF] introduced the idea of generating a ghost particle outside the computational domain with properties based on a real fluid particle close to solid boundaries. The ghost particle positions are obtained by mirroring the corresponding fluid particles across the boundary while its other properties are fixed in such a way as to handle simple boundary conditions such as full-slip or no-slip. The relative simplicity required to impose boundary conditions has led to this method becoming a common technique used to handle boundary conditions while also inspiring many other variants related to this idea. Boundary conditions for simple geometries are handled elegantly using this implementation, however for more complex boundaries, it becomes increasingly difficult to determine a ghost particle's properties in such a way that consistently leads to accurate boundary condition representations. In order to resolve these issues, several studies [START_REF] Marrone | -SPH model for simulating violent impact flows[END_REF][START_REF] Bouscasse | Nonlinear water wave interaction with floating bodies in SPH[END_REF][START_REF] Chen | A corrected solid boundary treatment method for Smoothed Particle Hydrodynamics[END_REF][START_REF] Chow | Incompressible SPH (ISPH) with fast Poisson solver on a GPU[END_REF] have made use of pregenerated ghost-interior particle pairs whereby the internal particle obtain information from the actual SPH particles being updated via an interpolation scheme. The ghost particles then take on a specific value to enforce boundary conditions just as in the traditional ghost particle method. Again, while this does have several advantages over the classic ghost particle approach, there are still certain geometries that lead to complications when applying this boundary implementation, however methods have been proposed to alleviate these issues by using similar interpolation techniques that are now applied to the generated ghost particles to fill any voids that may be present [START_REF] Marrone | -SPH model for simulating violent impact flows[END_REF].

In order to avoid using ghost particles altogether, the methods proposed by Randles and Libersky [START_REF] Randles | Smoothed particle hydrodynamics: Some recent improvements and applications[END_REF], Bonnet and Lok [START_REF] Bonet | Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations[END_REF] and Vila [START_REF] Vila | On particle weighted methods and smooth particle hydrodynamics[END_REF] uses linear truncation operations to reduce the order of errors present in the generalized SPH gradient from a 0 ℎ order scheme to a first order scheme allowing for the recovery of a linear function's derivatives exactly. Furthermore, Fatehi and Manzari [START_REF] Fatehi | Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives[END_REF] extended this idea to the evaluation of the Laplacian thereby allowing for a full first order consistent approximations of all differential operators present in the incompressible Navier-Stokes equations (NSE). Using truncation methods, Fatehi and Manzari [START_REF] Fatehi | A consistent and fast weakly compressible smoothed particle hydrodynamics with a new wall boundary condition[END_REF] were able to recover accurate solutions to a lid driven cavity problem as well as a locally expanding channel problem while Hashemi et al. [START_REF] Hashemi | Evaluation of a pressure splitting formulation for Weakly Compressible SPH: Fluid flow around periodic array of cylinders[END_REF] successfully simulated flow around periodically arranged cylinders. Gradient truncation methods have also been used in a study by Hashemi et al. [26] as a way to handle solid-fluid interactions whereby fluid forces on the solid bodies were resolved by inverting the fluid's momentum equation at solid surfaces as a means to resolve surface pressures. At least conceptually, these studies present the possibility to allow SPH to deal with extremely arbitrary system geometries and boundary conditions. In the context of fluid-solid coupling, it becomes possible to couple fluid-solid interactions for a generalized particle shape in a fully resolved sense. This has an advantage over drag correlation based approaches such as those seen in [START_REF] Cleary | Prediction of slurry transport in SAG mills using SPH fluid flow in a dynamic DEM based porous media[END_REF][START_REF] Robinson | Fluid-particle flow simulations using two-way-coupled mesoscale SPH-DEM and validation[END_REF][START_REF] Sun | Three-dimensional simulation of a solid-liquid flow by the DEM-SPH method[END_REF] since a drag correlation strategy requires simple particle shapes with know drag correlations. This usually leads to these approaches being applied to spherical particles only. When compared to the strategy proposed in [START_REF] Canelas | SPH-DCDEM model for arbitrary geometries in free surface solid-fluid flows[END_REF][START_REF] Canelas | Resolved Simulation of a Granular-Fluid Flow with a Coupled SPH-DCDEM Model[END_REF], whereby a solid particle is represented by a set of SPH particles with rigid body constraints between each other, clearly then, this approach does take particles shape into account, however, with these particles being treated identically to fluid SPH particles during the forcing calculations, this approach lacks explicit enforcement of boundary conditions and again leads to a type of unresolved coupling strategy.

Keeping this in mind, the main goal of this paper is to present an extension to the ideas proposed in [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF] in the context of a solid-fluid coupled environment with the focus being on implementing generalized boundary conditions in a 3D environment. Three validation cases and four numerical experiments are presented with the intent of investigating the applicability of the presented method as a strategy to couple 3D fluid simulations with general purpose DEM. In order to critically analyse the solid particle interactions, this paper mainly deals with individual solid particles. Validation is looked at in the context of fluid only systems, single particle solid-fluid systems, and multi-particle solid-fluid systems. As for the numerical experiments, due to the expected low fidelity of the SPH model that is to be coupled to the DEM solvers, the first experiment investigates the behaviour of SPH fluid resolution on the forcing conditions of the solid particle. The second experiment investigates the behaviour of a low fidelity coupled system with different particle shapes specifically dealing with spheres, convex polyhedra and non-convex polyhedra. The third experiment explores boundary interactions with a solid. Finally, a multi-particle system with five solid particles are simulated for both a set of spheres and non-convex particles. The solver presented is monolithically coupled to the multi-GPU DEM framework Blaze-DEM, and also makes use of the NVIDIA CUDA toolkit to run a parallel SPH solver on a general purpose GPU.

Following this, the contributions of this paper can be summarised as:

• Extending the solid-fluid coupling pressure boundary condition implementation of [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF] from a 2D to 3D environment,

• Utilizing gradient corrected surface derivatives to apply Neumann boundary conditions to velocity fields and • Investigating low fidelity fully resolved flow behaviour and drag characteristics of gradient corrected SPH around moving solid particles.

Numerical methods

In this paper an incompressible Newtonian fluid will be used as the model of interest. As such the incompressible Navier-Stokes equations (NSE) will be the partial differential equations (PDE) that will be required to be discretized by the SPH methodology:

= 0 (1) = - ∇ + ∇ 2 + (2)
Here, is the fluid density, is the flow velocity vector, is the pressure, is the kinematic viscosity of the fluid and is the acceleration due to external forces. It should be noted that unlike many stencil or element methods, the SPH methodology handles advection by integrating the particle position, and so requires the NSE to be cast in a Lagrangian frame. Clearly, true incompressibility is assumed in the above continuity and momentum equations, however, a common strategy to handle this constraint with SPH specifically is to weaken the incompressibility requirement in (1) leading to an approximate enforcement of incompressibility via a stiff equation of state. This requires a generalized continuity equation since is now required to evolve as well:

= -∇ ⋅ (3)
It should be stated that momentum equation (2) remains in its incompressible form implying that it is assumed that the fluid stress induced by a non-zero ∇ ⋅ is negligible.

Equations of state

In order to couple the continuity and momentum equations, a common strategy used in SPH revolves around the uses of an equation of state that allows for the coupling of the pressure and density fields of a fluid. Throughout this paper, the weakly compressible (WC) fluid model will be used to approximate an incompressible fluid and is given as:

= 0 2 0 -1 , ( 4
)
where 0 is a reference density, is the coefficient of expansion and is the speed of sound.

It should be mentioned that does not pertain to any physical fluid property, but can rather be seen as a numerical tuning factor used to balance the degree of compressibility and the speed at which information can travel through the fluid. A larger value stiffens the constraint on the incompressibility, but also increases the rate at which information propagates through the medium and thus limits the maximum allowable time step of the numerical integration scheme. In order to ensure accurate representation of an incompressible fluid, should be at least 10 times the maximum expected fluid velocity [START_REF] Monaghan | Simulating free surface flows with SPH[END_REF].

With regards to the coefficient of expansion, it can be seen that for the case of = 1.0, a linear dependency between pressure and density is enforced. When choosing larger, it can be seen that positive pressure respond strongly to density variations while negative pressures grow at a slower rate. This behaviour plays well with the SPH methodology as it allows particle very close to each other to be strongly repelled while nearly incompressible particles are allowed to fill the void in low density regions with less resistance. It has been shown by Monaghan [32] that = 7.0 produces an acceptable response for the case of approximating incompressible flow.

SPH discretization

The original SPH methodology introduced by Gingold and Monaghan [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF] approximate exact function values using a kernel function ( , ℎ) that approaches the Dirac-function in the limit as ℎ → 0 allowing nearby function value information to be interpolated to the desired location. In its discrete form where the general field ( ) is represented by a set of "particles" this can be written as:

⟨ ( )⟩ = ∑ ∈ ( ) ( -, ℎ) = ∑ ∈ , ( 5
)
where is the set of particle indices of particles inside the support radius of the ℎ particle, is the field value at and is the volume associated to the ℎ particle. While has physical meaning strongly related to the SPH particle mass and density in the context of standard SPH, it should be noted that by using the truncation methods discussed in section 2.3.1, this link is severed with simply becoming a scaling factor and so only requires ∝ -1 , however by not keeping consistent with and leads to a scaled identity matrix being recovered as the truncation matrix for fully supported particles. It should also be noted that in general ⟨ ( )⟩ ≠ since there is no explicit constraint fixing the field value evaluated at the particle position to the particle value itself. This paper uses a quintic kernel given as: where = ‖ ‖ 2 ∕ℎ is the scaled distance and = 3∕359 ℎ 3 is the kernel normalization parameter [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF].

( , ℎ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (3 -) 5 + 6(2 -) 5 + 15(1 -) 5 for 0 ≤ ≤ 1 (3 -) 5 + 6(2 -) 5 for 1 < ≤ 2 (3 -) 5 for 2 < ≤ 3 0 otherwise (6)
Schematically, this process can be described as seen in Figure 1, where a fluid body is approximated by many individual SPH particles. The ℎ particle is represented by the black particle, while the blue particles represent the collection of particles associated with and the orange circle indicates the support of the kernel function.

Since this approximation is itself a continuous field, it becomes possible to take derivatives of (5) directly, however, in order to improve the accuracy of the derivative approximation, applying the chain rule to ∇ ⟨ ⟩ where is a constant, the following form of the derivative can be found:

⟨∇ ⟩ = ∑ ∈ ( -)∇ (7)
Finally, in order to approximate the second order derivates of , consider the scheme proposed by Morris et al. [START_REF] Morris | Modeling low Reynolds number incompressible flows using SPH[END_REF] whereby a hybrid SPH-finite difference approach is used to estimate the function Laplacian. This approximation has been found to be more accurate and robust than taking second derivates of (5) which has been shown to produce inconsistent results simply based on the particle distribution [START_REF] Basa | Robustness and accuracy of SPH formulations for viscous flow[END_REF]. Modifying the original proposed Laplacian, the Hessian is given as:

⟨∇ ⊗ ∇ ⟩ = 2 ∑ ∈ , ≠ - 2 ⊗ ∇ , (8)
where = -, = ‖ ‖ 2 and ⊗ indicates the outer product.

SPH modifications

The standard SPH scheme as discussed in Section 2.2 has several issues related to truncated support at boundaries, noisy pressure fields and numerical instability. The proposed methodology introduces modification to standard SPH to alleviate some of these issues, specifically gradient correction, mass correction and particle shifting to address each of the above mentioned issues.

Figure 2: Contours of entries for -1 in 2D space.

Gradient correction

Following the analysis in [START_REF] Fatehi | Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives[END_REF] a Taylor series expansion of [START_REF] Cleary | Prediction of coupled particle and fluid flows using DEM and SPH[END_REF] shows that the proposed SPH gradient recovers:

⟨∇ ⟩ = - ∑ ∈ ⊗ ∇ ⋅ ∇ + … (9)
Clearly then, by defining a inverted truncation matrix as:

-1 = - ∑ ∈ ⊗ ∇ , ( 10
)
it becomes possible to recover first order accurate gradient via the linear mapping of the gradient ⟨∇ ⟩ = ⋅ ⟨∇ ⟩ with the superscript indicating truncation.

A visual description of this process in 2D is supplied in Figure 2 where the contour of the inverse truncation matrix entries over a single particle's support domain is shown. An immediate observation that can be made is how the scaled identity matrix would emerge from a fully supported particle via integration over the domain. Due to the diagonal entries simply being rotations of each other, it is clear that they would have identical values. Furthermore, the off-diagonal terms are zero due to the anti-symmetries present over the x-or y-axis leading to identical cancellation. It also shows how at flat boundaries where the support is not complete, the matrix diagonals would clearly be smaller in the direction of the surface normal (leading to larger diagonals in its inverse) and as such scale the appropriate derivatives, while corners would couple the directional derivatives based on the corner's angle as would conceptually be necessary to extrapolate information at a corner, and so account for missing information in a manner consistent with linear extrapolation.

The idea applied in ( 9) can be expanded to the Laplacian by following a Taylor series expansion of (8) leading to:

2 ∑ ∈ , ≠ - - ⋅ ⟨∇ ⟩ ⊗ ∇ = -∇ ⊗ ∇ ∶ ̂ + … (11)
where = ∕ and ̂ is the fourth order tensor defined as:

̂ = ∑ ∈ ⊗ ⊗ ⊗ ∇ + ∑ ∈ ⊗ ⊗ ∇ ⋅ ⋅ ∑ ∈ ⊗ ⊗ ∇
It can be seen that if the Hessian approximation in ( 11) is contracted with a new truncation matrix ̂ , a first order consistent Laplacian can be obtained given that ̂ ∶ ̂ = -.

Using the corrected SPH gradients and Laplacians, ( 3) and ( 2) can be discretized into a collection of ODEs. More specifically, the equations of motion for the ℎ particle becomes:

= - ∑ ∈ ( -) ⋅ ⋅ ∇ (12) = - 1 ∑ ∈ ( -) ⋅ ∇ + 2 ∑ ∈ , ≠ - - ⋅ ⟨∇ ⊗ ⟩ ̂ ∶ ⊗ ∇ + (13)
When considering the symmetries present in ̂ , it can seen that the linear system required to be solved can be reduced from 2 to ( 2 + )∕2 linear equations where is the spatial dimensions of the PDE. This implies that for each particle, a system of linear equations need to be solved for gradient correction as well as a system with ( 2 + )∕2 linear equation for Laplacian correction. The inverted truncation matrix as well as the symmetric components of the three sub-tensors of ̂ are constructed and saved in GPU memory. Due to the small size of these systems, a direct linear solver is used to solve these systems directly on the GPU as well.

Mass correction

While weakly compressible SPH (WCSPH) is known to have a tendency of producing noisy pressure fields, a modification presented in [START_REF] Fatehi | A consistent and fast weakly compressible smoothed particle hydrodynamics with a new wall boundary condition[END_REF] has been shown to smooth out the pressure field allowing for a more realistic pressure field approximation to be obtained. This method adds an additional term to [START_REF] Hu | An incompressible multi-phase SPH method[END_REF] with the modified equation given as:

= -⟨∇ ⋅ ⟩ -ΓΔ ⟨ ∇ ⋅ ∇ ⟩ - ⟨ ∇ ⋅ ⟨ ∇ ⟩⟩ (14)
The SPH approximation of divergence of the pressure gradient is given by:

⟨ ∇ ⋅ ∇ ⟩ = 2 ∑ ∈ , ≠ 2( -) ( + ) - ⋅ ⟨∇ ⟩ ̂ ∶ ⊗ ∇ (15)
while the SPH approximation of divergence of the SPH approximation of the pressure gradient is given by:

⟨ ∇ ⋅ ⟨ ∇ ⟩⟩ = ∑ ∈ ⟨ ∇ ⟩ - ⟨ ∇ ⟩ ⋅ ⋅ ∇ (16)
Looking at the correction term, it should be noticed this scheme evaluates the difference between two different measures of the pressure Laplacian and with it accordingly modifies the density field until this difference is insignificantly small. Furthermore, it has been shown that ( 14) rapidly converges to (12) [START_REF] Fatehi | A consistent and fast weakly compressible smoothed particle hydrodynamics with a new wall boundary condition[END_REF]. It should also be noted that Γ in ( 14) is a smoothing factor set to be larger than or equal to 1.

Particle shifting

To counteract the numerical instabilities caused by the fluid particles having a tendency to form clusters, a particle shifting algorithm introduced by Xu et al. [START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF] is used to redistribute particles ensuring approximately uniformly distributed particles. The particle positions are shifted according to the following scheme:

= ‖ ‖ Δ ∑ ∈ ̄ 2 2 , ( 17
)
where is a decrowding factor set between 0.001 and 0.1, ‖ ‖ is the maximum fluid velocity, Δ is the time step size, ̄ is the average distance to all particles in the ℎ particle's support radius and being the shift applied to the particle position.

The pressure and velocity fields are adjusted via a simple first order Taylor series expansion. Writing this in terms of a directional derivative with a step gives the form:

= ⟨∇ ⟩ ⋅ , ( 18
)
where is either the pressure or velocity fields. The density is modified via a first order Taylor series expansion of density in terms of pressure:

= = 0 -1 2 (19)

Integration methods

A similar two step predictor-corrector integration scheme outlined in [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF] has been used throughout this work. More specifically, ( 13) and ( 14) are used to determine a new set of intermediate value * from the current variables at the ℎ time step for density and velocity. The pressure values are also updated using (4) to reflect the respective densities. The same process is followed again to determine * * from * , and finally the new variable +1 is set to the average between * * and * . It should be emphasised that the time derivative * is obtained by re-evaluating ( 13) and ( 14) again with respect to the predictor particle values * . This leads to the following integration scheme:

* = + Δ (20) * * = + Δ * (21) +1 = 1 2 ( * * + * ) (22)
Where is either velocity or density and is the corresponding time derivative as presented in ( 13) and ( 14). The particle positions are updated only once every time step using the new fluid velocities:

+1 = + Δ +1 (23)
The time step is limited by either a CFL condition based on the maximum fluid speed or the rate of viscous diffusion given as:

Δ = min + ‖ ‖ , 2 , ( 24
)
where is the smallest particle spacing, ‖ ‖ is the maximum fluid velocity, is the kinematic viscosity and is a user specified stability factor between 0 and 1 that, for the case of this study, was kept between 0.1 and 0.2. The maximum expected fluid velocity is used to evaluate the CFL condition initially, after which the time step size remains fixed for the remainder of the simulation.

Boundary conditions

Since truncation is used to correct for the differential operators' sensitivity to particle distribution, Dirichlet boundary conditions can now be applied by simply using 1 layer of boundary particles with fixed properties set to the required value. This boundary implementation is used for all no-slip boundary conditions throughout the paper.

While the no-slip boundary condition requires a Dirichlet boundary condition for velocity, the surface pressure is enforced via a Neumann boundary condition. As shown in [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF], by contracting ( 13) with a surface normal ̂ at solid boundaries it becomes possible to manipulate this expression in such a way to determine surface pressures:

= ̂ ⋅ ∑ ∈ ⋅ ∇ -̂ ⋅ ⟨ ∇ ⋅ ⟩ -̂ ⋅ + ̂ ⋅ ̂ ⋅ ∑ ∈ ⋅ ∇ (25)
This in turn allows one to obtain surface densities without the need for evolving it via [START_REF] Leroy | Buoyancy modelling with incompressible SPH for laminar and turbulent flows[END_REF].

In this paper, this idea has also been extended to the velocity fields allowing for Neumann boundary conditions to be implemented for these variables. To do this, consider a field . Its directional derivative along the surface normal ̂ can be approximated as:

⟨ ⟩ = ̂ ⋅ ⟨∇ ⟩ = ∑ ∈ ( -) ̂ ⋅ ⋅ ∇ , ( 26
)
from which the surface value can be determined as:

= ̂ ⋅ ∑ ∈ ⋅ ∇ - | | | | ̂ ⋅ ∑ ∈ ⋅ ∇ , ( 27
)
leading to a generalised handling of Neumann (or mixed) boundary condition.

Solid coupling

As mentioned during the introduction, the SPH solver is built upon the Blaze-DEM framework with the goal being to couple continuum physics to DEM via SPH. Since the majority of this paper deals with single DEM particles, a detailed discussion of the DEM solver is omitted, however a thorough overview of the solver can be found in [START_REF] Govender | Development of a convex polyhedral discrete element simulation framework for NVIDIA Kepler based GPUs[END_REF][START_REF] Govender | Collision detection of convex polyhedra on the NVIDIA GPU architecture for the discrete element method[END_REF][START_REF] Govender | Blaze-DEMGPU: Modular high performance DEM framework for the GPU architecture[END_REF]. Rather, this section focusses on the coupling between the fluid and solid particles.

Three categories of shapes are considered in this paper which can be classified as follows: Analytical geometries, convex polyhedra and non-convex polyhedra. Analytical geometries are fully characterised by a surfaces satisfying some analytical equation of the form ( ) = 0. Convex polyhedra are described by a collection of vertices, edges and planar surfaces. Nonconvex polyhedra are described by a collection of convex polyhedra. While this representation allows the DEM solver to determine particle collision effectively, the SPH solver needs some discretized representation of these solid particles in order to interface with it. Since only a single layer of particles is required to implement the fluid boundary conditions, only the surfaces of these solid particles are meshed with SPH particles. For the analytical case this simply means placing the SPH particles via a parameterisation of the surface. For the case of a convex polyhedron, a bilinear mesh is used for any face that is rectangular, while other faces are decomposed into triangular subsection that are then meshed accordingly. Non-convex polyhedra follow this same logic, however only the outer surfaces are meshed in this case. Figure 3a to 3c shows the meshed surface for the case of an analytical, convex polyhedral and non-convex polyhedral particle representation, respectively.

All simulations presented in this paper consider the influence of gravity as this is the driving force for solid particle motion, however gravity is never applied as a body force to the SPH fluid particles. Rather, the hydrostatic forces are accounted for with a buoyancy argument whereby the body force on the solid particle is determined as the difference between the gravitational force on the solid mass and the displaced fluid mass where is the solid volume:

= ( - ) (28)
Several methods exists to account for dynamic fluid forces and moments ranging from the commonly applied drag force correlations [START_REF] Cleary | Prediction of slurry transport in SAG mills using SPH fluid flow in a dynamic DEM based porous media[END_REF][START_REF] Robinson | Fluid-particle flow simulations using two-way-coupled mesoscale SPH-DEM and validation[END_REF][START_REF] Sun | Three-dimensional simulation of a solid-liquid flow by the DEM-SPH method[END_REF] to realising solids as a collection of SPH particles with a rigid body constraint between them [START_REF] Canelas | SPH-DCDEM model for arbitrary geometries in free surface solid-fluid flows[END_REF][START_REF] Canelas | Resolved Simulation of a Granular-Fluid Flow with a Coupled SPH-DCDEM Model[END_REF]. While drag correlation provide an efficient coupling strategy for the case of spherical particles, there clearly exists an issue when it comes to generalizing this approach to other particle shapes. In the strategy where the solid is treated in an SPH context, this idea closely resembles that of other under-resolved strategies since the solid boundary conditions are never truly enforced.

Due to the strong focus on particle shape throughout the paper, a fully resolved coupling strategy is used. This leads to the dynamic fluid forces and moments being computed via integration of the pressure and surface shear fields.

In its continuum form, the macroscopic forces and moments are written as:

= ∫ Ω -̂ + (∇ ⊗ + (∇ ⊗ ) ) ⋅ ̂ , ( 29
) = ∫ Ω - × -̂ + (∇ ⊗ + (∇ ⊗ ) ) ⋅ ̂ , ( 30
)
where Ω is the solid particle domain, is the position of the solid's center of gravity, ̂ is the surface normal and = is the dynamic viscosity. In its discrete form, ( 29) and ( 30) are written as:

= ∑ ∈ -̂ + ⟨∇ ⊗ ⟩ + ⟨∇ ⊗ ⟩ ⋅ ̂ Δ , ( 31
) = ∑ ∈ - × -̂ + ⟨∇ ⊗ ⟩ + ⟨∇ ⊗ ⟩ ⋅ ̂ Δ , ( 32
)
where is the set of particle indices of SPH particles associated to the solid particle's boundary, ̂ is the surface normal obtained from the DEM particle at the position of the ℎ SPH particle, ⟨∇ ⊗ ⟩ is the SPH approximation of velocity gradients and Δ is the surface area of the DEM particle associated with the surface SPH particle. The surface pressure is obtained from [START_REF] Hashemi | Evaluation of a pressure splitting formulation for Weakly Compressible SPH: Fluid flow around periodic array of cylinders[END_REF], while the surface velocities are obtained from the rigid body kinematic equations, specifically:

= + × ( -), (33)
where and are the linear and angular velocity of the solid particle, respectively. For the case of multiple particle systems, an additional linear and angular forcing condition and needs to be considered due to particle contact with these forces being a function of some measure of overlap between particles. Again the reader is referred to [START_REF] Govender | Development of a convex polyhedral discrete element simulation framework for NVIDIA Kepler based GPUs[END_REF][START_REF] Govender | Collision detection of convex polyhedra on the NVIDIA GPU architecture for the discrete element method[END_REF][START_REF] Govender | Blaze-DEMGPU: Modular high performance DEM framework for the GPU architecture[END_REF] for a detailed discussion. This leads to the following equations of motion, handled by the DEM solver:

= + + , (34)
⋅ = + - × ⋅ , (35)
where is the mass of the solid and is the inertial matrix of the solid.

Finally, when considering the effects of lubrication on multi-particle systems, it has been shown that fluid forces cannot be resolved accurately when solid interfaces get closer than the fluid resolution [START_REF] Bian | A splitting integration scheme for the SPH simulation of concentrated particle suspensions[END_REF]. While this does require correction in cases where a particle's own momentum is dwarfed by the fluid pressure gradient between two solids, such as neutrally buoyant systems, many studies have shown that if the particle is sufficiently energetic, this correction is not needed to capture the relevant physics [START_REF] Canelas | SPH-DCDEM model for arbitrary geometries in free surface solid-fluid flows[END_REF][START_REF] Canelas | Resolved Simulation of a Granular-Fluid Flow with a Coupled SPH-DCDEM Model[END_REF][START_REF] Feng | The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[END_REF][START_REF] Zhang | Numerical investigation on the role of discrete element method in combined LBM-IBM-DEM modeling[END_REF][START_REF] Tang | A new drag correlation from fully resolved simulations of flow past monodisperse static arrays of spheres[END_REF]. The cases considered in this study falls into the latter category, and so avoids the need for lubrication forces leading to contact being handled solely by the DEM solver.

Validation

In this section, the SPH methodology discussed in Section 2 is applied to three systems with each of the corresponding result used to systematically validate the response of the solver. The first case investigated addresses the classic lid driven cavity problem. These results are compared to literature as a means to validate the implemented solver as as a standalone fluid solver.

Next, in order to investigate the solid-fluid coupling behaviour, a single particle system representative of a submerged infinitely long cylinder is released under gravity. The results are compared to a similar 2D investigation.

Finally, the behaviour of multiple submerged solid particles are investigated by simulating a two particle system with the goal of showing the so-called drafting-kissing-tumbling (DKT) phenomenon frequently explored in multi-particle fluid simulations.

Lid driven cavity

In order to validate both the SPH solver and the proposed Neumann boundary conditions, a 3D lid driven cavity systems with symmetry boundary conditions on both faces aligned with the XY-plane is simulated. Specifically, the = 100 case is considered in this study.

A schematic representation of the system of can be seen in Figure 4 with the system being fully described by:

• Geometry of = = 1.0m, = 0.1m,

• Fluid density of = 1000kg∕m 3 ,

• Kinematic viscosity of = 1.0 × 10 -2 m 2 ∕s,

• Speed of sound of = 10.0m∕s,

• Lid velocity of = 1.0m∕s 2 .

The SPH particle diameter 0 is defined as the distance between SPH particles in the initial configuration. For this study an SPH particle diameter of 0 = 1×10 -2 m was used leading to a 100k particle simulation. The results presented are all evaluated at the z-mid plane.

Contours of x-velocity, y-velocity and velocity magnitude can be seen in Figures 5a to 5c, respectively.

The mid-plane velocities are compared, in Figure 6, to the classical results of Ghia et al. [START_REF] Ghia | High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[END_REF] as well as the ISPH results of Xu et al. [START_REF] Xu | Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach[END_REF].

Comparing the results to those of Ghia et al., a maximum difference of 1.29% and 1.11% was found for the x-and y-velocity, respectively, while compared to the results of Xu et al., a maximum difference of 2.53% and 2.58% was found for the x-and y-velocity, respectively.

Similarly, with respect to the results of Ghia et al., an average difference of 0.48% and 0.36% was found for the x-and yvelocity, respectively. The average difference with respect to Xu et al. was found to be 0.79% and 1.56% for the x-and y-velocity, respectively.

Infinitely long cylinder falling in a viscous fluid

In order to validate the fluid-solid forcing methods described above, a single particle system is simulated and compared to the results of Hashemi et al. [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF]. Specifically, the case of interest in this section deals with an infinitely long cylinder falling through a viscous fluid under gravity.

The simulation is carried out in 3D space with symmetric boundary conditions applied in the longitudinal direction of the cylinder aligned with the z-axis. Figure 7a shows the schematic representation of the geometry. For this study, is set to be 20 × 10 -3 m and is set to be 50 × 10 -3 m, while is chosen based on particle size to ensure that 10 SPH fluid particle layers are present in the z-direction. The cylinder has a radius of 1.25 × 10 -3 m. The cylinder is started at rest from the center of the cavity in the x-direction at a height of ℎ = 40 × 10 -3 m, while the caps of the cylinder lies exactly on both walls in the z-direction. A time step of Δ = 1 × 10 -5 s was used. It should be noted that since the solid particle's forcing and integration scheme is handled by the Blaze-DEM solver, the "cylinder" is actually only an approximation generated from a polyhedron with 12 planar faces as shown in 7b, which are then meshed with SPH particles.

An SPH particle diameter of = 1.0 × 10 -4 m was used in this case. The support radius of the simulation is set to be 3 0 , varying from that used in [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF]. This is justified by the fact that here a 3D simulation is conducted as opposed to the 2D simulation conducted in [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF].

The physical parameters used in these simulations are as follows:

• Fluid density of = 1000kg∕ 3 ,

• Solid density of = 1250kg∕m 3 ,

• Kinematic viscosity of = 1.0 × 10 -5 m 2 ∕s,

• Speed of sound of = 1.0m∕s,

• Gravitational acceleration of = 9.81m∕s 2 .

(a) (b) (c) Figures 8a to 8c shows the -and -velocities and velocity magnitude of the mid plane section, respectively.

Figure 9a shows the comparison between the solid particles vertical height at similar SPH particle resolutions while Figure 9b shows a comparison between the solid particle velocity and the results of Hashemi et al., along with the results presented by Glowinski et al. [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for the simulation of flow around moving rigid bodies: Application to particulate flow[END_REF] used for comparison in [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF].

Considering the fact that the geometry of the solid particle is not an ideal cylinder, but rather a polyhedral approximation, it can be seen that even with the additional complexity of a full 3D simulation requiring symmetric boundary conditions without the use of ghost particles, the results, although showing slight variations to other previous studies, clearly capture both the solid and fluid behaviour with a 2.8% and 5.5% variation from the results of Hashemi et al., in position and velocity, respectively.

Drafting-kissing-tumbling phenomenon

With section 3.2 showing that fluid forces are approximated well in the case where momentum transfer is only considered between the fluid and the solid particle, this section expands upon these results by introducing a second solid particle allowing for an investigation into the solver's behaviour under both particle-particle and particle-fluid momentum transfer modes. The physical parameters of the system are chosen as:

• Fluid density of = 1000kg∕ 3 ,

• Solid density of = 1100kg∕m 3 ,

• Kinematic viscosity of = 5.0 × 10 -6 m 2 ∕s,

• Speed of sound of = 1.0m∕s,

• Gravitational acceleration of = 9.81m∕s 2 .

When comparing the results found in this section to those in literature, it is important to take note of the remarks in [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF], where it was stated that this type of system's behaviour during the kissing stage is highly sensitive to its contact model used to the point where different contact models can lead to very distinct behaviour. Both cylinders are again approximated by the 12 sided polyhedra seen in Figure 7b, hence some variation in particular after contact is expected when compared to other studies that had more idealized cylinder representations.

Figures 11a, 11b and11c shows the velocity magnitude of the fluid during intervals representative of the drafting, kissing and tumbling phase, respectively.

Figures 12a and12b show the x-and y-positions of the solid particles, respectively, while Figures 13a and13b show the particle velocities in the x-and y-directions, respectively. The results are Compared to the results of Singh et al. [START_REF] Singh | Distributed Lagrange multiplier method for particulate flows with collisions[END_REF] and Hashemi et al. [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF].

As expected, it can be noticed that there is a very good agreement between all three simulations during the drafting phase, however this method see a response emerging once the kissing stage is initiated that is unique from either of the two other observed behaviours. In this case, it can be seen that the kissing phase transitioned into the tumbling phase very rapidly. This can again be attributed to the fact that the cylinder is approximated by a polyhedron where contact is made directly between 2.5 × 10 -4 141K Medium:

1.75 × 10 -4 290K Fine:

1.0 × 10 -4 894K

Table 1: Particle sizes and particle count for Section 4.1.

edges of the leading faces. This, in turn, leads to a highly unstable contact point where slight perturbations can cause highly dynamic behaviour. Of course, these perturbations do occur during the slow displacement of the leading particle to the side and, instead of the smooth force expected from idealized cylinders, the particles receives a discontinuous force leading to the kissing phase being cut short. However, due to the very similar kinematic conditions at the start of the tumbling phases, this still led to very good correlations between the tumbling phase of the current results and those of [START_REF] Hashemi | A modified SPH method for simulating motion of rigid bodies in Newtonian fluid flows[END_REF] albeit slightly dampened.

Numerical experiments

Four sets of results are presented in this section. The first set of results investigate the effects of SPH fluid resolution on the drag force experienced by a solid particle. This is done by considering a fixed geometry that is approximated using varying sizes of SPH particles.

The second investigation looks at a comparison between the flow behaviour and drag characteristics of different particle shapes at low SPH fluid resolution. This is done in order to determine whether even at low fidelity it is possible to resolve flow around various particle shapes to a degree where it is possible to capture the differences between particle shapes.

The next set of results deals with effects that boundary distance to the solid particles have on the drag experienced by the solid. These results give good indication of the sensitivity of the method's drag force predictions on pressure and shear. It is well known that these properties have a strong dependency to boundary effects when the infinite fluid approximation is sufficiently violated.

The last set of results presented deals with the case of multiple particle-fluid interactions and looks at the effects that particle shape have on the general behaviour of a cluster of solid particles.

Finally, it should be noted that all results presented are obtained from full 3D simulations.

Infinitely long cylinder falling in a viscous fluid

Building upon the results of section 3.2, this study looks at the same case considered during validation. Again, Figure 7a shows a schematic representation of the system of interest with the same geometry described before. Specifically, = 20 × 10 -3 m and = 50 × 10 -3 m, while is chosen in such a way to ensure that 10 SPH fluid particle layers are present in the z-direction. With this study focussed on the effect of varying the SPH particle diameter, the solid and fluid properties as well as other parameters for the numerical method remains the same as those used in 3.2.

Table 1 provides a summary of the particle sizes and total number of SPH particles for the three cases.

Figure 14a and 14b shows the cylinder height and velocity as a function of time for the three cases, respectively. From these results, it can be seen that although similar trends are present between the different simulations, the coarse resolution simulation differs by 16.33% from the fine resolution simulation with regards to the solid particle's terminal velocity. It can also be noted that a higher SPH resolution results in a slower moving solid particle (or equivalently in a larger drag force). For this reason, a comparison between particle drag characteristics are also provided. The non-dimensional drag coefficient is defined as:

= 1 2 2 , ( 36
)
where is the total drag force, is the projected area of the solid in the flow direction and is the vertical velocity of the cylinder. The value is averaged over the period where the velocity is seen to remain constant. The Reynolds number based on diameter at the terminal velocity is also provided. The results are summarised in Table 2. As mentioned above, these results show that coarser SPH fluid resolutions leads to an underestimation of drag forces which in turn leads to faster fall speed as captured by the larger . This is possibly due to the coarser fluid resolution not being able to capture the sharp gradient due to boundary layer flow around solid particles and so underestimate the surface shear on the particle. This is explored in Figure 15, where the tangential velocity along the cylinder is plotted at regular intervals over the cylinders surface for the fine case. The dots in the middle indicate the fluid resolution for the fine case. Clearly then, in the area where large gradient are present, even the fine resolution case only resolves this layer with 5-6 particles over a large portion of the leading face. Comparing this to the coarse and medium case, it is clear that these cases would have 2-3 and 3-4 particles in this region, respectively, leading to a much poorer estimation of the boundary layer effects and surface shear.

It should be noted that even when considering the lower fidelity models, while the results may not match the numerical values, the results give strong insight into the behaviour and trends of the system and so may possibly allow one to extract general Coarse Medium Fine 17.03 15.71 14.64 10. [START_REF] Singh | Distributed Lagrange multiplier method for particulate flows with collisions[END_REF] 12.18 14.02

Table 2: Coefficient of drag for an infinitely long cylinder at various resolutions.

Figure 16: Schematic representation for the various particles falling in an enclosed cavity behaviours of the system as would be of interest when dealing with a much larger number of solid DEM particles in situations such as design or optimization.

Drag characteristics of convex and non-convex particle shapes

The results presented in this section investigates the trends among the drag characteristic for various particle shapes under otherwise similar conditions. This is done to determine the degree to which the described fluid-solid coupling can represent various particle geometries and capture the differences due to particle shape, and so investigate the possibility of its applications to large scale DEM simulations.

With the trends between particle shapes being the main focus of this investigation, it was decided that low fidelity models are to be used to evaluate the drag behaviours of the particle shapes. Similarly to the previous section, a particle is started at rest at the height of ℎ = 40 × 10 -3 m and is then allowed to fall under the effects of gravity through the fluid. For this investigation, the particles are positioned in the center of the cavity while all walls have a no-slip boundary condition applied to them. The schematic representation of the system can be seen in Figure 16 with both and set to 15 × 10 -3 m and set to be 50 × 10 -3 m.

Four different particle types are investigated, namely a sphere, a cylinder, a ring and a cross as depicted in Figure 17 with the exact geometrical properties tabulated in Table 3. As with the previous set of results, the cylinder is approximated by a convex polyhedron while the ring and cross are both represented by non-convex polyhedra. The cross is built using five identical cuboids and the ring is built using 12 identical trapezoidal prisms. The sphere is described analytically. Arguably, it would be beneficial to keep the solid particles at a similar mass as well as projected area, however since it was desired to keep the material properties constant over all simulations and due to the surface of the solid particle being meshed with SPH particle, it was chosen that the area and mass of the solids be only approximately equivalent, but rather ensure that the SPH Sphere:

Cylinder:

• particle spacing remain as uniform as possible over the solid particles surface, due to the strong relationship between density and particle concentration under the SPH methodology.

The sphere was found to have a projected area of = 4.909 × 10 -6 m 2 and mass of = 1.022 × 10 -5 kg. Table 4 shows the comparison between the solid particle's projected area and mass with respect to that of the sphere. The characteristic length used to determine relative to sphere's diameter is also provided here. The solid and fluid properties used in these simulations are kept the same as discussed in section 4.1, while a time step size of Δ = 2 × 10 -5 s is used. An initial SPH particle spacing of 0 = 2.5 × 10 -4 m was used for each simulation while the support radius for this study was fixed at 2.4 0 . This led to each simulation requiring 693K SPH fluid particles.

Figure 18a and 18b shows the height and velocity profiles for the different solid particles. The velocity magnitude at = 0.7s for the various shapes for both the side and top view can be seen in Figures 19 and20, respectively, and, as with the previous section, to compare the drag characteristics, ( 36) is used to obtain the average coefficient of drag for each particle during the time when the velocity remains constant. These results can be seen in Table 5. From the results presented, it becomes clear that even under a low fidelity model, the fully resolved SPH flow around the different particle shapes leads to drastically different behaviours of the solid and fluid particles. As expected, the sphere has a considerably lower drag when compared to any of the other particle shapes. When comparing the cross and cylinder, it can be seen that the front section of the cross has the effect of breaking the flow and somewhat channelling the fluid around the larger flat faces. This leads to a more aerodynamic behaviour when compared to the purely flat face of the cylinder. Finally, it can also be seen that due to the viscous nature of fluid, the shear generated on the internal surface of the ring has a considerable effect on the drag response leading to the largest drag force of the four particle types under consideration.

As a supplement to the investigation above, the solid particle responses are also compared to well known results found in literature for the case of the sphere. As a way to explore the space surrounding the chosen design point, two additional simulations are preformed with varying solid particle density. Specifically, a density of = 1100kg∕m 3 and = 1500kg∕m 3 are used. The terminal velocity is given as the solution to:

( -) = 1 2 2 (37)
A correlation for based on as presented in [START_REF] Morrison | Data Correlation for Drag Coefficient for Sphere[END_REF]

The vertical positions and velocities for the three cases are show in Figures 21a and21b, respectively, with the red lines showing the corresponding theoretical terminal velocities of the particles. A clear trend is observed here with a somewhat expected behaviour where slower moving particles show a much closer correlation to the theoretical results. Specifically, a difference to the theoretical terminal velocity of 0.11%, 9.21% and 13.52% was observed for the case of = 1100kg∕m 3 , = 1250kg∕m 3 , and = 1500kg∕m 3 , respectively. This can again be explained by the simplified surface dynamics of the slower moving particle being approximated well by the coarse SPH resolution, while the sharper and more dominant boundary layer effects of faster moving particles require a much finer resolution to accurately describe the surface forcing conditions. Of course, this correlation does assume that the fluid body is infinite and so one would expect boundary effects to influence the solution to some degree (especially for the faster moving particles), however this does serve at the very least as a good way to estimate realistic particle responses.

Boundary effects on drag characteristics of rings

Due to the no-slip boundary conditions applied to the cavity walls in the previous section, it is expected that the results obtained differ from that of an infinite fluid approximation. In this section, the dependency of drag on the system geometry is investigated by varying the wall distances to the solid particle. This will have the effect of changing the flow area available to the fluid and thus lead to a pressure build up. Of course, it is also expected that the boundary positions will have an effect on the viscous momentum diffusion. Specifically, this investigation will only consider the slowest moving particle, namely the non-convex ring.

Due to these expected fluid behaviour it can be seen that this investigation serves as a good means of exploring the sensitivity of the fully resolved low fidelity coupling to the internal fluid forcing condition, specifically the surface pressure and shears.

Two additional simulations are compared against the results found in Section 4.2. The same geometry as seen Figure 16 is used for this study with = scaled to 2∕3 and 1∕2 of the original geometry.

Solid particle heights and velocities can be seen in Figure 22a and 22b, respectively. Furthermore, based on outer diameter and the for the three cases at the constant velocity can be found in Table 6. The influence of the choked flow due to the shrinking domain can be clearly seen with the expected trend of an increasing drag with the shrinking of the domain presenting itself. Quantitatively, the terminal velocity was reduced by 14.9% for the case of the 2/3 domain, while the 1/2 domain showed a decrease of 35.3%. A massive change in drag characteristics can be notices with the 2/3 and 1/2 domain case showing an increase of 38.19% and 139.1% to when compared to the original domain, respectively.

Finally, the y-velocity of the fluid for the three cases can be seen in Figure 23. It can be seen that the shrinking domain has the effects of constraining the flow development around the solid particle, leading to a much stronger interaction between the fluid travelling with the solid and that being displaced. Also, it can be noticed that the fluid is travelling faster around the outside of the solid, as expected with the reducing flow area, which in turn has a strong effect on surface shear as well.

Multi-solid particle simulation

Previous results allow one to probe the physics of a fluid reacting to a solid particle travelling within it, however with the essentially single parameter models based on terminal velocity that can be derived for the solid particle behaviours, these flow cases do not explore complex coupling in multiple spatial dimensions. For this reason, this section focus on the qualitative behaviour of a cluster of particles falling through a fluid with the particles starting initially in a staggered configuration. This leads to complex non-symmetric flow patterns influencing particle trajectories in a much more dynamic sense.

It was chosen to compare cases using the sphere and ring of section 4.2 as these shapes where found to have the extremal drag characteristics when considering the set of particles investigated.

The schematics as seen in Figure 24 show the initial configuration of the solid particles. A box with dimension = 50×10 -3 m, = 10 × 10 -3 m and = 8 × 10 -3 m houses the particles and fluid. The solid particles are initialised at the mid plane in the z-direction while the in-plane positions are characterised by ℎ = 33 × 10 -3 m, = 3 × 10 -3 m and = 4 × 10 -3 m. As with section 4.2, a time step of Δ = 2 × 10 -5 s and an initial spacing of 0 = 2.5 × 10 -4 m is used. This leads to 241K fluid particles taking part in the simulation. Fluid and solid properties as discussed in section 4.1 are used for these simulations as well.

No-slip boundary conditions are applied to the top and bottom planes and the planes normal to the z-direction. A symmetry boundary condition is applied to the planes normal to the x-direction essentially representing the case of many scattered particle settling in a large fluid bath. Snapshots of the fluid velocity magnitudes for the spheres and rings can be seen in Figures 25 and26, respectively.

As expected from the drag results in section 4.2, the set of spheres collectively travel faster through the fluid when compared to the set of rings. Furthermore, it can be noticed that the rings now have a strong rotational response due to the non-symmetric flow tripping the particle from its conditionally stable state. This further reduces descent rates due to energy now being stored and dissipated in a rotational mode as well. The sphere also experiences small rotational effects, but due to its rotationally symmetric nature this effect has a much less drastic consequence on the particle behaviours.

Finally, it should be noted that contact was not forbidden over the course of the simulations. Although these simulations deal with many fewer particles than a typical DEM simulation, contact did occur between solid particle and was handled effectively by the DEM solver with the SPH solver behaving as expected as well. This leads to higher degree of confidence in this method as a tool for dealing with larger scale solid-fluid coupled environments.

Conclusion

An investigation into the applicability of gradient corrected SPH for the use of solid-fluid coupling in 3D was presented. It was found that although lower fidelity models tend to present a lower than expected coupling force, the observed trends still allow information about the physics present in the system to be extracted. Furthermore, it was shown that high resolution models yield accurate representations of fluid-solid interactions. It was also shown that even at low fidelity scales for both the case of single and multiple particles, the presented fluid-solid coupling was able to represent the various fluid and solid particle behaviours associated to specific solid particle geometries which, at least for the case of the non-convex particles, are sufficiently complex to pose challenges for other boundary implementations. Finally, as a results of the strong comparison between the 2D and 3D results for the lid driven cavity and infinitely long cylinder, at least as an initial proof of concept, it has been found that by utilizing gradient truncation, it is possible to extended beyond only the recovery of surface pressures to also apply a variety of other Neumann boundary conditions for velocity fields as well.

With this in mind, it is proposed that future work focuses on an investigations into the behaviour of larger scale DEM-SPH coupling making use of low fidelity fully resolved flow as not only an analysis tool, but also as a design and optimization tool. Furthermore due to the flexibility of the boundary condition method, investigations into other multi-physics domains such as heat or elasticity may also be explored.

Figure 1 :

 1 Figure 1: Visual representation of the SPH discretization.

Figure 3 :

 3 Figure 3: SPH representation of solid particles described by an (a) analytical shape, (b) convex polyhedron and (c) non-convex polyhedron.

Figure 4 :

 4 Figure 4: Schematic representation of the lid driven cavity problem.

Figure 5 :

 5 Figure 5: Contours of (a) x-velocity, (b) y-velocity and (c) velocity magnitude for the lid driven cavity evaluated at the z-mid plane.

Figure 6 :

 6 Figure 6: Mid plane velocities of the lid driven cavity.

Figure 7 :

 7 Figure 7: System description with (a) the schematic representation for infinitely long cylinder falling in an enclosed cavity and (b) the polyhedral representation of the cylinder.

Figure 8 :Figure 9 :

 89 Figure 8: Fine scale mid plane profile for (a) fluid velocity x, (b) fluid velocity y, (c) fluid velocity magnitude at time = 0.63s.

Figure 10 :

 10 Figure 10: Schematic representation for the DKT problem.

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: Velocity magnitude of the fluid at (a) = 0.27s (b) = 0.54s and (c) = 0.66s.

Figure 14 :

 14 Figure 14: Response of an infinitely long cylinder's (a) height and (b) velocity as a function of time to a varying SPH resolution.

Figure 15 :

 15 Figure 15: Tangential velocity of the fine case over the from face of the cylinder.

  Radius of = 1.25 × 10 -3 m • Radius of = 1.25 × 10 -3 m • Height of ℎ = 1.75 × 10 -3 m Ring: Cross: • Inner radius of = 0.725 × 10 -3 m • Outer radius of = 1.45 × 10 -3 m • Height of ℎ = 1.75 × 10 -3 m • Cuboid height of ℎ = 1.0 × 10 -3 m • Cuboid base of = 2.25 × 10 -3 m • Cuboid width of = 0.75 × 10 -3 m

Figure 17 :Figure 18 :

 1718 Figure 17: Schematic representation of the particle geometry for the (a) sphere, (b) cylinder, (c) ring and (d) cross.

Figure 19 :Figure 20 :Figure 21 :

 192021 Figure 19: Fluid velocity magnitude of the (a) sphere, (b) cylinder, (c) ring and (d) cross at time = 0.7s.

Figure 22 :

 22 Figure 22: Comparison of ring (a) height and (b) velocity for the case of the original, 2/3, and 1/2 domain.

Figure 23 :

 23 Figure 23: Fluid y-velocity of the (a) 1/2, (b) 2/3 and (c) full domain cases at time = 0.7s.

Figure 24 :

 24 Figure 24: Schematic representation for the various particles falling in an enclosed cavity

Figure 25 :

 25 Figure 25: Fluid velocity magnitude for the case with multiple spheres at (a) = 0.175s, (b) = 0.35s, (c) = 0.525s and (d) = 0.7s.

Figure 26 :

 26 Figure 26: Fluid velocity magnitude for the case with multiple rings at (a) = 0.175s, (b) = 0.35s, (c) = 0.525s and (d) = 0.7s.

Table 3 :

 3 Geometric characteristics of the various particle type.

		Sphere Cylinder Ring Cross
	∕	1.000	1.000	1.009 1.031
	∕	1.000	1.050	1.060 1.031
	∕	1.000	1.000	1.160 0.900

Table 4 :

 4 Comparison between particle projected area and mass relative to the sphere.

Table 5 :

 5 Coefficient of drag for various particle shapes.

Table 6 :

 6 Coefficient of drag for Ring at different boundary geometries.

	1/2	2/3	Full
	8.12 10.68 12.55
	53.97 31.19 22.57