Modelling of Non-Spherical Particles in Dilute Non-Colloidal Suspensions Using SPH

Suspensions and their applications can be found in many engineering, environmental or medical fields. Considering the special field of dilute suspensions, possible applications are cement paste or procedural processes in the production of medication or food. While the homogenized behavior of these applications is well understood, contributions in the field of pore-scale fully resolved numerical simulations with non-spherical particles are rare. Using Smoothed Particle Hydrodynamics as a simulation framework we therefore present a model for Direct Numerical Simulations of single-phase fluid containing non-spherically formed solid aggregates. Notable and discussed model specifications are the surface-coupled fluid-solid interaction forces as well as the contact forces between solid aggregates. Moreover we simulate and analyze the behavior of dilute non-colloidal suspensions of non-spherical solid particles in Newtonian fluids. The focus of this contribution is the numerical model for suspensions and its implementation in SPH. Therefore shown numerical examples present application examples for a first numerical analysis of influence factors in suspension flow. Results show that direct numerical simulations reproduce known phenomena like shear induced migration very well. Moreover the present investigation exemplifies the influence of concentration and form of particles on the flow processes in greater detail.


Introduction
Suspensions and their applications can be found in many fields of mechanical, civil and environmental engineering. Coarse-graining and therefore simplyfing the multi-phase response of the mixture, suspensions are often constitutively described on a continuum-scale as non-Newtonian fluids with either shear-thinning or shear-thicking behavior [1,2]. Such continuum models present a good approximation when it comes to general investigations of fluid flow of suspensions on length scales significantly larger than the characteristic particle diameter. However there is no information about the microscopical influence of solid-solid contacts as well as hydro-mechanical fluid-solid interactions which could lead to well-known phenomena like de-mixing [3], particle formation like hydroclusters [4,5] and shear-wall migration [6,7]. To overcome limitations of classical macroscopic continuum formulations, we present a multi-scale modeling approach for fully resolved suspensions. The model takes into account a Newtonian carrier fluid coupled with a repulsive force for the solid-solid contact (calculated from a potential). Contrary to standard continuum approaches for suspensions, both the fluid and the solid phase are fully resolved which leads to the possibility of investigating the physical properties of the flow and the solid contacts individually.
Various numerical methods to investigate suspension flow in Computational Fluid Dyamics (CFD) exist. Often either a fully or coupled Discrete Element Method (DEM) is used. Thus contributions exhibit Lattice-Bolzmann-DEM approaches [8], SPH-DEM approaches [9,10], Finite-Element-DEM approaches [11,12] or Particle-Finite-Element-Method approaches (PFEM) [13] to name only a few. All of these methods have in common, that they are either Eulerian grid-based methods (FEM, LB) which are limited when large deformations are invloved or free surface flow is dominant or else are not able to discretize the surrounding fluid without a coupling algorithm (DEM). To overcome these limitations we present a quasi-incompressible Smoothed Particle Hydrodynamics (SPH) approach. SPH as a meshless Lagrangian method presents remarkable advantages when it comes to fully resolved solid-fluid interfaces and interactions [5,14].
Based on that, SPH is the ideal choice to model suspensions since the carrier fluid as well as the solid particles can be discretized with SPH particles (see also Fig. 2). The fully resolved model allows to build rigid solid aggregates of arbitrary forms using solid SPH-particles and to construct the fluid phase with fluid SPH-particles. Therefore each material point P(x, t) is either occupied by the fluid or solid phase. This renders usage of surface coupled hydromechanical forces easy and allows to simulate aggregates with complex non-sperical forms. The presented SPH model [15] is implemented into the highly scalable numerical tool HOOMD-blue [16,17].
There are various contributions considering the rheological behavior of suspensions [1,18,19]. Considering non-colloidal suspensions, Tanner [20] as well as Hofmann [21] give a detailed introduction and comparison between theory and experiments. As described there, colloidal suspensions contain small solid particles (d < 1µm) and therefore can be described using Brownian forces. Different to that, non-colloidal suspensions are defined by larger particle diameters (d ≥ 1µm) and a low aspect ratio between particle volume and particle surface compared to colloids. Thus a larger Péclet number which represents the ratio of mechanical to Brownian forces is characteristic for non-colloidal suspensions which results in the assumption that fluid mechanical forces dominate the flow and Brownian forces can be neglected. Normally, non-Brownian suspensions at small Péclet number show a Newtonian behavior with a viscosity independent from their shear rate [20]. However, Bian et al. [4] recently found a shear rate dependency including the formation of hydroclusters when considering these suspensions under confinement. Similar to Tanner [20,22], most publications consider high concentrated, i.e. dense suspensions where other effects dominate the effect and impact of contact between single solid particles. Moreover the named publications only consider the ideal case of spherical solid particles. Therefore this work focuses on the influence and flow behavior of dilute suspensions with only a few particles with arbitrary and specially non-spherical forms.

Solid particle motion in single-phase Newtonian fluid
The rheological behavior of suspensions strongly depends on the concentration C s , i.e. the amount of solid particles in the fluid. It can be computed as number concentration C s , using the number of entities (e.g. particles) of a constituent n s divided by the total number of SPH particles N or as volume concentration, using the volume of a constituent dV s divided by the total volume of the mixture dV In general, a suspension can be discretized as a Newtonian or non-Newtonian fluid with embedded solid particles which influence the flow behavior [1]. Considering only dilute suspensions (C s = 1-5 %) with a Newtonian solvent, the mechanical behavior can be described using the mass and momentum conservation with a constitutive equation for the Cauchy stress tensor. After discussing underlying conservation equations, the numerical model is introduced.

Balance equations for single-phase fluid flow
Suspensions consist of a solvent (liquid phase) and a solute (solid particles). In contrast to a true solution, the solute in a suspension will not be dissolved and can float freely. The solvent can be discretized as a Newtonian or non-Newtonian fluid. In literature, the two phases of a suspensions with Newtonian solvents are often described by non-Newtonian fluid formulations [23]. These formulations present homogenized models for fluid and solid phases. Since in this case contact forces can not be considered at all, we use fully resolved direct numerical simulations based on a coupled formulation of the momentum conservation including an extra part for solid-solid contact instead.
Underlying conservation equations are the balance of mass and the balance of momentum, i.e. the Navier-Stokes equations for quasi-incompressible single-phase bulk fluid flow in local forṁ where ρ f is the fluid density, v f is the flow velocity and ρ f g are the body forces. T denotes the Cauchy stress tensor which is split into an extra or non-equilibrium part and into a pressure or equilibrium part T = T eq + T neq . While the equilibrium part is given by the non-equilibrium part of the stress tensor is referred to as viscous extra stress assuming a divergence-free velocity field v f . This leads to the updated form of the balance of linear momentum for a homogeneous quasi-incompressible fluid with constant dynamic viscosity µ f Since the fluid phase is considered to be barotropic, the fluid pressure p can be computed as a function of the fluid density ρ f using an equation of state in form of Tait's equation where ρ f 0 is the initial fluid density, c the (numerical) speed of sound and p 0 is a constant background pressure to avoid negative pressures in the system. Moreover γ = 7 is a common choice for quasi-incompressible flow [24,25].

Introduction to suspensions
As introduced, additional to the solvent part, the solute part of a solution has to be described. In general the underlying conservation equations (Eqs. (3) and (4)) are also valid for the fluid-solid interaction since we consider a surface-coupled approach where no mass exchange between interfaces and bulk phases exist (Fig. 2). Solid-solid interaction within the solute part of the suspension is described by a repulsive force approach [26] This commonly used law computes a force acting alongside the vector x ab that connects the center of masses of solid bodies P a and P b . F 0 denotes the magnitude of the repulsive force (dependent on the considered boundary value problem), τ −1 sets the active range of the force and s is the distance between the surfaces of solid bodies P a and P b . The repulsive force F rep ab is independent of the surrounding fluid and considered by an additional term in the balance equation of linear momentum (Eq. (7)) which therefore can be expressed by means of a whole-domain formulation with α = {s, f} for the solid and fluid phase respectively.

Numerical modeling of suspensions using SPH
For the purpose of studying the effects of heterogeneous particles in a single-phase fluid, we employ pore-scale resolved Direct Numerical Simulations (DNS). Our method of choice is Smoothed Particle Hydrodynamics (SPH) using an implementation that has previously been shown to accurately reproduce effective hydraulic properties of porous media [27][28][29]. In particular, our implementation incorporates the SPH scheme proposed in [14] together with the boundary conditions proposed in [30] and presents an extension of the highly optimized Molecular Dynamics software HOOMD-blue [16,17]. The implementation targets both CPU and GPU computation and was verified for several representative test cases performed on the supercomputers Hazel Hen (HLRS, Germany) and BinAC (Tübingen, Germany) [15]. The use of SPH is motivated by the fact that we consider the solvent consisting of arbitrary formed particles instead of spherical solid particles. The meshless nature of SPH renders the discretization of heterogeneous solid particles comparatively trivial. Given the Lagrangian nature of SPH, implying that nonlinear convective terms are not required to be modelled, SPH is rather stable at finite Reynolds numbers. Since a detailed study of the numerical scheme is beyond the scope of this contribution, the reader is referred to [31,32] for technical details concerning details of the implementation. The numerical implementation and therefore discretization of introduced equations results in four cases of SPH particle interaction, which will be dealt with further below: A. Fluid particle with fluid particle, B. Fluid particle with not moving solid particle, C. fluid particle with (moving) solid particle and D. Moving solid particle with moving solid particle.

SPH for single-phase fluid (case A and B)
In SPH the discretization of governing partial differential equations gives rise to a set of interacting collocation points (referred to as SPH particles). At each unstructured integration point, every field function f (x) can be represented by local interpolation using convolution with the Dirac delta function The Dirac delta function δ(x − x ) is replaced by a kernel function W with compact-supported and smoothing length h (Fig. 1) leading to an approximation of the field function f (x) Since SPH is a particle method, the kernel representation is converted into a spatial discretization and therefore every field function transforms into a particle property f i = f (x i ). At point P i the propertiy f i is computed by summation over function values of all neighboring particles P j with particle volume V j Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 15 March 2020 doi:10.20944/preprints202003.0240.v1 Figure 1. Representation of the kernel function W with compact support κh and dependent on particle distance r ij .
Moreover continuous differential operators can be discretized in the same way and hence transform into viscous and pressure short-range interaction forces F V ij and F P ij resprectively. Following this, the local force balance resulting from the discretization of the local balance of momentum (Navier-Stokes equations (7)) can be expressed as The shown discretization applies to discrete integration points x i , referred to as SPH particles P i , of mass m i and subject to the advection particle velocity v i . All forces acting on SPH particle P i can be considered as a combined force F SPH i . This contains the volumetric force on a particle F B i = m i g as well as the pairwise dissipative viscous interaction forces F V ij and the pairwise conservative pressure interaction forces F P ij both acting between particle P i and its nearest neighbor particles P j . They are discretized by and In Eq. (15) and (16) the particle number density n i , which satisfies n i = 1/V i , is used. The effective viscosity is represented by the harmonic mean of two neighboring particles P i and P j , while the pressure interaction forces, which approximate the pressure gradient, are calculated by an antisymmetric gradient stencil to satisfy the conservation criterion F P ij = − F P ij . The equations are evaluated for different cases in which the neighbor particle P j either is a fluid particle (x j ∈ Ω f , case A) or ghost particle (x j ∈ Ω G , case B) and therefore part of a periodic boundary or part of a (solid) Dirichlet boundary. For the latter case (x j ∈ Ω G ), a replicated viscosity µ * j , a fictitious particle velocity u * j (both Eq. (15)) and a fictitious particle pressure p * j (Eq. (16)) is used to ensure no-slip and no-penetration boundary conditions. Further details are again out of the scope of the presented publication and can be found in [30,33,34].

SPH for suspension flow (case C)
After introducing the discretization for single-phase fluid flow, the solid-fluid interactions (case C) will be derived for the usage of SPH. Every solid particle P s can be discretized by a rigid collection of solid SPH particles with a no-slip and no-penetration boundary condition at all fluid-solid interfaces Γ fs (Fig. 2). Doing so, earlier introduced particle forces (F B i ) and short-range interaction forces (F V i , F P i ) can be calculated for all solid SPH particles with fluid neighbors. Since solid bodies are rigid, containing solid SPH particles P i can not move independently. Therefore one can examine a total force F s (based on the balance of momentum in Eq. (14)) and a total torque M s acting on every solid particle P s arose by the surrounding fluid particles. and Viscous and pressure interaction forces, F V i and F P i respectively, are summed up for the amount of containing SPH particles n s inside a solid body and r i = x i − x M is the distance of each solid SPH particle to the bodys center of mass (see Fig. 3). Using the moment of inertia J M of a solid body the translational and angular acceleration (v a ,ω a ) and velocity (ṽ a ,ω a ) of a solid particle P s a can be calculated by and finally the total velocity v i of each SPH particle P i of the solid body can be updated to v i =ṽ a +ω a × r i .
The implementation generally provides no-slip and no-penetration boundary conditions at all fluid-solid interfaces Γ fs . Thus fluid velocity components (v f = 0) as well as the relative velocity between solid and fluid phase states to be zero (v f = v s ) on Γ fs .

Solid-solid interactions in SPH (case D)
While calculating suspensions also solid bodies P s a and P s b will interact with each other. To regard this, an additional term was included into the momentum balance equation (10). The repulsive force is discretized as presented in Eq. (9) taking into account x ab = (x b − x a )/x ab . Based on that, the discrete balance of momentum (Eq. (14)) is updated to Moreover the computation of the total force F s is extended by the addition term of the repulsive force while the calculation of the total torque remains the same as in Eq. (18).

Numerical analysis of non-colloidal suspensions
The presented SPH approach is implemented as SPH extention in the highly optimized Molecular Dynamics software HOOMD-blue [16,17]. It is used to perform 3-dimensional fully resolved Direct Numerical Simulations of a boundary value problem (BVP) of suspension flow in a channel as shown in Fig. 4. The scope of this contribution is the presentation of a numerical model for suspensions and its implemementation in SPH. Therefore here shown BVPs present application examples for a first numerical analysis of influence factors (as concentration and form) in suspension flow.
The simulation domain consists of the fluid domain which is periodic in e x -and e y -direction and is limited by solid particle layers (walls) in direction of e z . The channel has a height h = 0.02m which is chosen to be the characteristic reference length L ref . We select a forth order Wendland kernel representation with compact support κh = 3.4 dx [35] and a resolution of 60 particles over the channel height. Thus this results in an initial particle distance dx = L ref /60. For our investigations, wall particles are accelerated with a constant velocity v x = 0.005 m/s such that occurring Reynolds numbers stay constant as Re = 100. Initial viscosity of the fluid phase is chosen to be µ f = 0.001 Pa/ s and initial density to be ρ f 0 = 1000 kg/m 3 . Since we are interested in investigating the flow behavior dependent of the suspended particles, we perform simulations with various numbers of solid aggregates (two or four solid bodies) as well as with various aggregate forms (spherical and ellipsoidal). As presented in Tab  particle) and therefore the concentration varies between 1.2 % (two spherical aggregates) and 1.8 % (four ellipsoidal aggregates) but always stays in a dilute regime.

Discussion of relevant model parameters
The simulation results are analyzed in terms of resulting (steady state) velocities and motion of the solid bodies. The intention is to use the fully resolved SPH simulations to reproduce and investigate phenomena like de-mixing and migration processes. Nevertheless we do not aim to perform a full analysis of the flow regimes. Fig. 5 shows the velocity vector as well as structure of containing solid bodies for the steady state. Velocity parts in e y -and e z -direction remain to be negligible since they are sufficient small compared to the velocity in direction of e x which results to be linear over the height h.  analytical solution represents the flow of a fully Newtonian fluid, the best approximation is given in the simulation without solid bodies (red line). Simulations containing solid bodies exhibit a plateau of zero velocity in the centerline region of the simulation domain (z = ±0.003m) due to the there accumulated solid aggregates. This region remains to be wider for a higher number of spherical aggregates compared to ellipsoidal aggregates. Simulation results show that the numerical model is capable of reproducing the effect of shear wall migration, where solid aggregates move to regions with lower shear rates [6,7].  and with four solid aggregates (spherical vs. ellipsoidal) where the start position remains the same and only the form changes. Simulations are performed until suspensions are almost in a stationary regime (approx. 20 s), while simulations without solid particles reach a steady state much faster. We observe an oscillation around the centerline region (z = 0) in the motion of all solids until they arrange themselves in this line. Due to their form, ellipsoidal bodies (dashed line) perform a non-linear rotational motion and therefore take more time to arrange in the centerline than spherical bodies (straight line). Moreover the higher concentration ensures a faster arrangement since solid bodies are slowed down in in their motion due to particle collisions between the aggregates. Considering Fig. 7 compared to the motion of surface points in Fig. 8 one can perceive that rotation and motion in direction of e x proceeds parallel. It can be observed that the rotation stabilizes faster than oscillation around centerline (z=0). Since spherical solids are rotationally symmetric, their surface gives less windage than surfaces of ellipsoids and therefore rotation stops relatively fast. Nevertheless the results show that even with only a small number of solid bodies effects like shear-induced particle migration (as observed in experiments [36]) become visible.

Conclusion
The current contribution presented a mesh-less approach for Direct Numerical Simulations of the motion of solid particles (grains) in a Newtonian carrier fluid. Underlying conservation equations for linear momentum of single phase flow were derived and a contact force was included into the model. Fluid-solid motion interaction is captured via interfacial forces between solid and fluid SPH particles. The model was implemented in a quasi-incompressible SPH scheme and an example of a dilute suspension under shear impact was investigated. The resulting velocity profiles of the dominating velocity v x were compared to numerical results of single-phase fluid flow (without solid aggregates) of Newtonian fluid as well as the analytical solution for this case. All in all the numerical results show that suspension flow differs to single-phase fluid flow mainly in the centerline region (z = 0) where the flow velocity is zero. Due to hydro-mechanical forces the aggregated particles are accelerated and migrate to this region where they stay. Dependent on their form this process happens faster or slower. Spherical particles therefore move much faster to the center region than elliptical particles do. Due to the non-linear rotation non-spherical particles also enforce a non-zero velocity in the centerline region as it can be observed in Fig. 6. This reproduces well known effects where solid aggregates move to regions with lower shear rates, reported by [6,7] as well as general shear-induced particle migration, as reported by [36]. However, the number of solid particles is to small to investigate phenomena like de-mixing as observed in [3]. Therefore future investigations will include simulations with more solid particles to analyze these effects more closely.
To come closer to real-life applications, like for example concrete pumping processes, ongoing investigations are BVP where the fluid flow is induced by a volumetrical force corresponding to a pressure difference. Interesting details are the migration of particles dependent on their form and size as well as the evolution of flow profiles dependent on the concentration. Ongoing research moreover includes the implementation of a lubrication correction to stabilize the model by correcting hydro-dynamic forces acting between two solid particles as proposed in [4,5].