Comp. Part. Mech. DOI 10.1007/s40571-015-0097-9

Development of a coupled discrete element (DEM)–smoothed particle hydrodynamics (SPH) simulation method for polyhedral particles Benjamin Nassauer1 · Thomas Liedke1 · Meinhard Kuna1

Received: 14 April 2015 / Revised: 11 December 2015 / Accepted: 17 December 2015 © OWZ 2016

Abstract In the present paper, the direct coupling of a discrete element method (DEM) with polyhedral particles and smoothed particle hydrodynamics (SPH) is presented. The two simulation techniques are fully coupled in both ways through interaction forces between the solid DEM particles and the fluid SPH particles. Thus this simulation method provides the possibility to simulate the individual movement of polyhedral, sharp-edged particles as well as the flow field around these particles in fluid-saturated granular matter which occurs in many technical processes e.g. wire sawing, grinding or lapping. The coupled method is exemplified and validated by the simulation of a particle in a shear flow, which shows good agreement with analytical solutions. Keywords DEM · SPH · Fluid particle coupling · Liquid-solid flow · Polyhedral particles

1 Introduction In many technical processes, like lapping, grinding or wire sawing, we deal with mixtures of fluids and solid particles. While the solid particles are used as abrasives to machine the surface of a given workpiece, the fluid is needed to transport these particles into the kerf and to remove debris and heat from the process zone.

B

Thomas Liedke [email protected] Benjamin Nassauer [email protected] Meinhard Kuna [email protected]

1

Institut für Mechanik und Fluiddynamik, TU Bergakademie Freiberg, Lampadiusstr. 4, 09596 Freiberg, Germany

The actual removal process is the sum of uncountable single indentation procedures of abrasive particles and depend on the shape, the hardness and size distribution of the particles as well as the properties of the carrier fluid. In order to gain a deeper understanding of the process and to predict macroscopic machining parameters like the removal rate or the surface roughness, a numerical method is needed which is capable of simulating sharp-edged particles as well as the carrier fluid. For the first part, a discrete element method [6,7] with polyhedral particles was developed and presented in [24], which models the influence of the fluid in terms of drag forces due to a static shear flow without modelling the fluid itself. A more realistic approach needs to incorporate the interaction of the fluid with the particles i.e. by means of a coupled fluid– particle simulation. In general, there are two methods for that [35]. In the partitioned approach, the equations for the fluid flow and the equations describing the displacement of the particles are solved separately with two different solvers. The interaction is done by interpolating the fluid forces due to pressure or shear of the particles. The solution is obtained iteratively using the result of one simulation as boundary condition for the other one until the solution converges. The second approach is called monolithic approach or direct coupling. Here, as the name already suggests, the equations of the fluid flow and the motion of the solid particles are solved simultaneously by constructing a set of equations for the total system, which can be solved at once. The partitioned method provides a certain modularity because several kinds of already existing software can be used. But for the monolithic approach specialized software is needed which makes the development more complex but leads usually to faster code because no iterative procedure is needed.

123

Comp. Part. Mech.

Several approaches for the partitioned coupling of DEM with fluids can be found in the literature. E.g. Grof et al. [10] described a coupling of DEM with finite volume method and Galindo-Torres [8] presented a coupling of DEM with sphero-polyhedral particles and the lattice boltzmann method. The direct coupling of DEM with a fluid can be achieved using the dissipative particle dynamics e.g. Bierwisch et al. [2] or the smoothed particle hydrodynamics (SPH) method. Examples of these can be found in [12,17,18,27,28]. Cleary and Morrison [5] described a DEM–SPH coupling as well, but used a coupling via drag forces instead of a direct numerical simulation. The main advantage of combining DEM and SPH is that both methods are particle based and the underlying time integration scheme is quite similar. In both methods, the particles are driven by forces from neighbouring particles and their equations of motion are integrated numerically. The only difference is that the inter-particle forces in SPH are derived from the Navier–Stokes equations and those forces in DEM are due to particle contact process [12]. Usually the particles used in DEM are spherical and thus the methods mentioned above deal with the coupling of SPH to spherical DEM particles. For polyhedral DEM particles, the coupling with SPH particles is more complex and is related to the problem of defining wall boundaries to a SPH simulation. For this, two basic concepts exist. The first one is to fill the walls with SPH dummy particles to ensure the support of the kernel in the wall region. The other concept considers the non-vanishing surface integral when smoothing the flow quantities close to the boundary or to introduce artificial repulsion forces to prevent particles from entering the wall volume [1]. Regarding the first concept, [14] introduced the usage of ‘ghost’ particles to mirror real fluid particles at the wall surface. This method can be used to model symmetry and periodic boundary conditions but needs simple-shaped wall geometries. The ghost particles have to be recreated with every timestep. The second method has the advantage that only the wall surface needs to be described by a layer of SPH particles [20] and the disadvantage of the missing kernel supporting the area inside the wall. In [34], different concepts are compared and it was found that the usage of dummy particles for the kernel support combined with a pressure response from the boundary surface as proposed by [1] currently matches experimental results best. For the coupling of SPH particles to polyhedral DEM particles, we also use dummy particles (virtual SPH particles) located on the polyhedron surface combined with a pressure-based repulsive force. But our proposed method uses independent virtual particles for each SPH particle near the wall and the missing kernel support volume is compensated by using a one-dimensional kernel function for the virtual particles.

123

The described method is used in a software that has been previously developed by the authors [23,24] which itself is based on the open source software LIGGGHTS [13].

2 Representation of polyhedral particles For the discrete element method, usually spherical particles are used. If non-spherical particles are needed, they are composed of smaller spheres which are glued together [2,15,33]. With this method the precise reproduction of a sharp-edged particle is not possible and a realistic approximation needs a great amount of sub-particles, which makes this procedure inefficient. Another way is the use of polyhedral particles which are usually defined by a hierarchical set of connected vertices, edges and facets. If the shape of the particles is restricted to convex polyhedra, which fits very well with most granular material, it is possible to use a set of half spaces to define the particles [26]. This representation reduces the memory consumption and provides the possibility for new contact detection algorithm as presented in [3] or [23] and [24]. In the following, a short summary of the half-space representation of convex polyhedra, which is then further used for the SPH coupling, is given. For a more detailed explanation please see [23,24]. Each half space can be described by the inequality (x − aα ) · aα ≤ 0,

(1)

where x is the vector to an arbitrary point and aα is the normal vector from the center of mass to plane α. With α ranging from one to the number of boundary faces, the complete set of half spaces is defined. A 2D example for of this is given in Fig. 1, which shows a triangle defined by three half spaces. More complex polyhedra can be defined by a greater set of plane vectors. Randomly shaped polyhedrons, which are used, for example, to simulate abrasive grains, can be created by defining a set of randomly oriented plane vectors with varying lengths. An example for such a particle, constructed with eighteen planes is given in Fig. 2. In theory, these particles can be arbitrarily complex, but is useful to reduce the number of planes (facets) to a low value in order to reduce the computational costs. In our case, a maximum of twenty-one planes was needed to create realistic random particles. Based on the geometric representation, which is equivalent to the usual vertex-face structure, the physical properties such as volume, mass, inertia tensor and so on can be calculated using standard procedures. The contact force calculation between colliding particles is based on the intersection volume of the overlapping geometries of both particles. To calculate the overlapping volume, a

Comp. Part. Mech.

aα=1

aα=2 aα=3

x2 x1

c

Fig. 1 Representation of particle geometry by a set of half spaces, 2D-Example. Each half space is indicated by a shaded hatching. The combination of three half spaces lead to the red marked triangularshaped polygon. c is the vector from origin to the centre of the particle. α1 , α2 and α3 are vectors defining the planes

Fig. 2 Example of an arbitrarily shaped convex polyhedron. The arrows indicate the plane normal direction

Fig. 3 Determination of overlap, 2D-simplification. The overlapping region is a subset of half spaces and can be calculated using the same techniques as for a single polyhedron. n¯ f is the direction of force, derived from the surface normals ni of single half spaces and applied at the centre of mass Co of the overlap volume Vo . d indicates the indentation depth of the overlapping region. Ai is the surface Area of face i

The direction n¯ f of the normal force is defined as the normalized vector obtained by integrating the surface normal niA over the part A of the surface of the overlapping region that belongs to one particle (see Fig. 3). For polyhedra, this means, n f is obtained as the sum of the surface normals ni over all parts Ai , weighted by the area: 

kind of a boolean intersection is performed. This again results in a polyhedra which consists of a subset of half-spaces from both colliding particles (see Fig. 3). Volume and centroid of the intersection can therefore be calculated using the same technique as used for a usual polyhedral particle. A generalized contact law as proposed in [23] is used to calculate the contact forces. It is based on the Hertz contact law and depends on the volume of the overlapping region and on the indentation depth and includes a viscous damping term which is capable to account for material damping. The contact force in normal direction is given by    Fn = E ∗ k Vo d 1 + cd˙ .

(2)

With 4 k = √ ≈ 0.752 3 π

(3)

and E being the Young’s modulus, Vo the overlap volume, d the indentation depth, d˙ the magnitude of the relative velocity between both particles in direction of the normal force and c the constant for the damping. The point of application of the force is the centre of mass of the overlap polyhedra.

ni ds 1  = n¯ f =  A Ai ni .  ni ds  i Ai A

(4)

i

The friction force is modelled by a smoothed Coulomb friction model and is correlated with the relative velocity between the particles. No information about the contact history or other forces on the particle is needed. To avoid numerical instabilities, a smooth algebraic function is used for the transition between static and kinetic friction. For more information on this see [23]. To reduce the computational effort, a spherical dummy particle is defined for each polyhedral DEM particle. The radius of this particle equals the distance from centre to the most distant vertex of the according polyhedron. The physical properties like mass and inertia tensor of the polyhedron are mapped to the spherical dummy. This has the great advantage that almost all simulations steps ranging from (preliminary) contact detection to the integration of translational motion can be performed by standard DEM algorithms. Only the calculation of the rotational motion differs due to the the nonspherical contribution of the inertia tensor. Therefore, a more advanced integrator is used which is capable of handling rigid body dynamics. For this, we used the module “fix-rigid”, which is part of the LIGGGHTS DEM software package [13]

123

Comp. Part. Mech.

and is based on a rotational Verlet integration scheme (see e.g. [30]). The previous explained contact force calculation comes into use only after a contact between two dummy particles is detected. In this case, the underlying polyhedral particles are checked for collision in a second contact search step, and repulsive forces are calculated only if the contact is verified. For more details see [23]. The equivalent forces due to a ‘contact’ with SPH particles will be described in the following.

3 Smoothed particle hydrodynamics Smoothed particle hydrodynamics (SPH) is a meshless simulation method developed in 1977 by Lucy [16] and Gingold and Monaghan [9] for the simulation of astrophysical problems. Since then it has been advanced to various other fields like e.g. fluid mechanics. In the following, a short introduction to SPH is given in order to lay out the fundamentals before the coupling to our polyhedral discrete elements is described. For a deeper understanding of SPH see e.g. [9,16–18] or others. The basic principle of the method is that a physical property Φ(r) at a point r of a continuum can be approximated ˜ by a smoothed quantity Φ(r) which is build by integral interpolation of the surrounding volume Ω [17], see Fig. 4, ˜ Φ(r) ≈

Ω

Φ(r )W (r − r , h)dr .

In a further approximation, the integral interpolant in equation (5) is replaced by a summation interpolant over b discrete points in space. Φ(rb )W (ra − rb , h).

(6)

Each of these points is associated with a finite volume which is expressed by the fraction of mass m b and density ρb . The value Φb is the quantity of Φ at position rb . At this stage of approximation it becomes clear that these points can be interpreted as particles with certain physical properties. For the choice of the kernel there are several possibilities like a Gaussian kernel [9,18], a cubic spline kernel [4,17, 19,31] or a quintic spline kernel [21,27,36]. We use a cubic spline kernel of the form [17]: ⎧ 3 2 3 3 ⎪ ⎨1 − 2 q + 4 q , if 0 ≤ q ≤ 1 W (q) = σ 41 (2 − q)3 , if 1 < q ≤ 2 . ⎪ ⎩ 0, else

(7)

With the dimensionless distance q=

r |ra − rb | = h h

(8)

and the normalizing constant σ which depends on the dimension of the integration:

σ =

Ω

ρb

b

(5)

In this equation, r is any point in Ω and W (r − r , h) is the kernel function. Its width parameter h is often referred to as ‘smoothing length’ and its value determines the influence radius of the kernel function.

 mb

˜ a) ≈ Φ(r

⎧2 ⎪ ⎨ 3h ,

if 1D if 2D

10 2, ⎪ ⎩ 7π1 h , π h3

(9)

if 3D

Because mass and density are fixed for particle b the gra˜ a ) can be expressed dient of the summation interpolant Φ(r in terms of the gradient of the kernel function itself: ˜ a) = ∇ Φ(r

 mb b

=

r−r

ρb

 mb b

ρb

Φ(rb )∇W (ra − rb , h) Φb ∇Wab

(10)

The gradient of the kernel, found by chain rule, is

x2

r r x1

V, m, ρ

Fig. 4 Schematic view of SPH particles representing a continuous volume Ω

123

1 ra − rb ∂ W (q) h |ra − rb | ∂q ⎧ ⎪−3q + 49 q 2 , if 0 ≤ q ≤ 1 ra − rb σ ⎨ 3 = − 4 (2 − q)2 , if 1 < q ≤ 2 |ra − rb | h ⎪ ⎩ 0, else.

∇Wab =

(11)

The flow of a fluid can be described by the Navier–Stokes equation. For an incompressible fluid the Lagrangian form of the Navier–Stokes equation is

Comp. Part. Mech.

∇p η ∂v =− + ∇ 2 v + g. ∂t ρ ρ

(12)

With velocity vector v, time t, pressure p, density ρ, dynamic viscosity η and gravity vector g, the terms on the right hand side of the equation are associated with pressure, friction and gravity. All of them can be expressed in SPH formalism as follows [11]. Using chain rule, the pressure term in (12) can be expressed as ∇p =∇ ρ

  p p + 2 ∇ρ. ρ ρ

(13)

Using (10) this is transformed into      m b pb pa  m b ∇p ∇Wab + 2 ρb ∇Wab . (14) = ρ ρb ρb ρa ρb b

b

From this equation, the symmetric expression according to Monaghan [18] follows: ∇p  = mb ρ b



pb pa + 2 2 ρa ρb

 ∇Wab .

(15)

Mass and density of the particles are defined at the beginning of the simulation. While the mass is constant throughout the simulation and thus conservation of mass is satisfied automatically, the density can change. Following [18], the change of density is calculated by  dρa = m b (va − vb ) · ∇Wab . dt

(16)

From this follows that if two particles approach each other the density of the particles rises, while if they move away from each other the density of the particles diminishes. Using the Tait equation of state, the pressure pa can be calculated from the density by pa = B

ρa ρ0

χ

 −1 ,

 m b (ηa + ηb )(ra − rb ) · ∇Wab η 2 ∇ v= (va − vb ) ρ ρa ρb (ra − rb ) · (ra − rb ) b  m b (ηa + ηb )(va − vb ) ∂ Wab . = ρa ρb h|ra − rb | ∂q b

(18) Here, ηb is the dynamic viscosity of particle b. The viscosity model consists of a combination of a finite difference approximation of a first derivative and the SPH form of a first derivative. By this model, second derivatives of the kernel, which could lead to instabilities, are avoided [21]. Inserting Eqs. (15) and (18) in Eq. (12) gives the full Navier–Stokes equation in SPH formalism    pa pb dv  + 2 ∇W ab −m b = dt ρa2 ρb b  m b (ηa + ηb )(va − vb ) ∂ Wab + g. + ρa ρb h|ra − rb | ∂q

(19)

For the time integration of equation (19), we use the Velocity Verlet method.

4 Interaction between solid particles and fluid 4.1 General considerations

b



For the viscous term in the Navier–Stokes equation, Morris et al. [21] presented the expression

(17)

where ρa is the density of the particle and ρ0 is the reference density. The model parameter B is a kind of a compression module and χ the adiabatic exponent. Because the Reynolds number of the considered fluid flow is low, the value χ = 1 is used to obtain stable simulations [21]. To get a large stable timestep the model parameter B, which is associated with the fluid speed of sound and the density by B = c2f ρ, can be set to small values. According to [18], it should be ensured that the Mach number, which is the fraction of fluid velocity v f to speed of sound c f , is still below 0.1.

The numerical procedure of a DEM simulation and a SPH simulation is almost the same. Therefore, DEM particles and SPH particles can be handled in the same simulation framework. The interaction between two DEM particles or two SPH particles is the same as in separate simulations and requires no special treatment. But additionally, the interaction between a DEM particle and a SPH particle has to be defined. Several methods for this solid–fluid coupling can be found in the literature. For example, the mirroring of SPH particles at the solid surface. Hereby, the velocity of the created imaginary particles has to be chosen in a way that the no slip boundary condition is satisfied at the surface (see e.g. [14,32]). In a similar approach, SPH particles are placed on fixed positions inside the solid (see e.g. [21,27,36]). Another approach is to place boundary particles on fixed positions of the solid surface (see e.g. [18,22,28]). A drawback of these methods is that additional particles are needed. In simulations with large solid–fluid boundary interfaces, this leads to significant increase of numerical effort. In the following, an approach for the solid–fluid interaction is presented which avoids additional permanent SPH particles.

123

Comp. Part. Mech.

4.2 Determination of the virtual contact point on the polyhedron surface In order to apply forces from a SPH particle to a polyhedral particle, a contact point has to be defined. Therefore, the polyhedron point which is closest to the SPH particle is identified and, in case it lies within the sphere of influence of the SPH particle, this point is defined as the virtual contact point. Hence, forces and moments acting on the DEM particle due to a SPH particle are applied to this point. To reduce the computational effort, for both, DEM and SPH particles, a surrounding sphere is defined, which enables the software to deal with them like with ordinary DEM particles (see Sect. 2). For SPH particles, usually the influence radius, which is in our case two times the smoothing length, is used. For the polyhedral particles, we use the most distance vertex to define the surrounding sphere. The great advantage of this is that all conventional contact algorithms can be applied before the specialized contact algorithm for SPH particles and polyhedral particles comes into use (see Fig. 5). Only after a contact between two surrounding spheres is detected, additional steps are required. The exact position of the contact point may be located on a face, edge or vertex of the polyhedron. In the following, the algorithm used to determine the contact point pc as well as the corresponding vector rc will be described in more detail. Figures 6 and 7 illustrate the procedure. Both are 2D simplifications of three dimensional objects. Depending on the

Fig. 6 Part of a polyhedron with the intersection of two half spaces (2D simplification), exemplifying different distances for the contact algorithm between SPH particle and polyhedron. The grey shaded area indicates the space occupied by the half spaces j defined by the plane vectors a j . pci are the vectors from centre c to the contact point of SPH particle i and ri j are vectors pointing from SPH particle i to half space j. rv is a vector pointing to the vertex defined by the intersection of the half spaces. (Color figure online)

Fig. 7 Schematic view to the determination of the contact point between polyhedron and SPH particle. Because of the 2D simplification the figure shows both, the SPH to-face and the SPH to-edge contact scenario

context, the 2D vertices correspond to edges in 3D and 2D edges correspond to planes.

Fig. 5 Schematic view of different relative positions between SPH particles and polyhedral DEM particle (2D simplification). For each Particle a surrounding sphere is used for the contact search. Only if these spheres overlap, further considerations are necessary. In the example, particle 1 is disregarded in the first step, while particle 2 and 3 are checked for contact with the polyhedron. Particle 3 is considered to be in contact with the polyhedron

123

– In the first step, the contact to a polyhedron face is checked. For this, a loop over all particle half spaces is performed and the distance of each half space to the SPH particle is calculated by the following procedure. – Given are: The connection vector rd between the centre of the polyhedron and the centre of the SPH particle, the defining vector a for the half space and the vector p1 to one vertex of the polyhedron which is part of the half space

Comp. Part. Mech.

– An auxiliary vector p∗1 directed from SPH particle to vertex P1 is defined by p∗1 = p1 − rd

(20)

– This vector p∗1 is projected on a, resulting in vector rc , which is directed from SPH particle perpendicular to the boundary of the half space rc =

p∗1 · a a |a|2

p12 = p2 − p1

(22)

If the value is positive the SPH particle is inside the half space, which obviously means that this half space is not a valid contact surface and the next half space in the loop is considered. Such a invalid case is exemplified in Fig. 6 by the pair particle p S P H 1 and half space 2. – The range of influence of the SPH particles equals the doubled smoothing length 2h. So, if the half space is valid and the distance rc is greater than the doubled smoothing length 2h, the SPH particle is out of range and no contact is established. In this case, the contact search can be stopped completely because this also means that no valid contact point exists at all. if |rc | > 2h: no face contact.

(25)

– The projection of p∗1 on p12 leads to vector b, which is directed from the possible contact point to vertex P1 . b=

p∗1 · p12 p12 |p12 |2

(26)

– The point is a valid contact point if it is between the vertices P1 and P2 , i.e. the scale factor of b in equation (26) has to be between −1 and 0: −1 ≤

p∗1 · p12 ≤ 0. |p12 |2

(27)

If not, the point is not a valid contact point and the search is continued with the next edge. – If the contact point is valid, a vector rc is calculated: rc = p∗1 − b.

(28)

– As for the face contact, the distance between SPH particle and contact point is checked by

(23) |rc | > 2h,

– If the SPH particle is closer, then 2h the vector pc to a preliminary contact point can be calculated by pc = rd + rc

– Given are: the connection vector rd between centre of polyhedron and the SPH particle, the vertices p1 and p2 of the considered edge. – The auxiliary vector p∗1 is calculated according to Eq. (20). – The vector between vertex p2 and p1 reads

(21)

– From the direction of vector rc we can deduce if the SPH particle is inside or outside of the considered half space. If the SPH particle is outside, the scaling factor in equation (21) has to be negative: p∗1 · a < 0. |a|2

is performed, and for each face, another loop over all edges is performed. For each edge, the following steps have to be taken:

(24)

– As the contact point is on the surface of the polyhedron it must be part of all defining half spaces. That is why the inequality (1) is checked for all remaining half spaces. – If the point is part of all half spaces, the point is considered to be a valid face contact point and the contact search is finished. Otherwise, the contact search continues with the main loop and checks for a face contact with the next half space. – If no face contact was found, the edges of the polyhedron have to be checked (see particle p S P H 2 in Fig. 6). Again a loop over all half spaces, respectively, polyhedron faces,

(29)

and if the distance is greater than the doubled smoothing length, the point is discarded. – If the point is a valid contact point, the according vector is calculated and the search is finished: pc = rd + rc

(30)

– If neither a face nor an edge contact was found, at last a check for a vertex contact is performed. Therefore a loop over all vertices is performed and the following steps are executed: – Given are: The connection vector rd between centre of polyhedron and the SPH particle and the vertex p1 which is also the possible contact point: pc = p1 .

(31)

123

Comp. Part. Mech.

– Vector rc is calculated by rc = p1 − rd .

(32)

– The vertex loop continues with the next vertex if |rc | > 2h

(33)

– else a valid contact point was found and the search is finished. 4.3 Force calculation between SPH particle and polyhedron In normal direction an elastic force is used, which hinders the SPH particles from entering the polyhedron interior (Fig. 8). To relate the normal stiffness of this elastic force to the stiffness of the fluid the following considerations are made: The SPH particle is simplified as a cube with edge length l and mass m. The density ρ0 of the SPH particle in the initial state is m ρ0 = 3 . l

(34)

If the particle approaches the surface and the deformation is only perpendicular to the surface, then the density ρ of the particle is m ρ= , 2|rc |l 2

(35)

where rc is the vector from the centre of the SPH particle onto the surface of the DEM particle. Using the equations (17), (34) and (35), the pressure p can be calculated as  p=B

l 2|rc |

χ

 −1 .

(36)

Fig. 8 Contact force calculation in normal direction between rigid DEM particle and SPH particle

123

Fig. 9 Contact between SPH and polyhedral particle (2D simplification). For SPH particles with no polyhedral contacts (SPH 1), the usual scheme from Sect. 3 is used. For particles which are in contact with a polyhedron (SPH 2 and SPH 3), a virtual particle is created (indicated blue and green) and used for force calculation. The dotted lines indicate which neighbour particles are used for the summation. Note that SPH 2 and SPH 3 cannot ‘see’ the virtual particles from each other. (Color figure online)

This pressure acts on an area of size l 2 . Thus the force F on the solid is ⎧ χ   rc ⎨l 2 B l − 1 2|r | |rc | , if 2|rc | < l c F= (37) ⎩0, else The equilibrium of forces between SPH particle and polyhedron is ensured by using Newton’s third law i.e. applying the opposite force to the SPH particle. For the purpose of great generality, the described procedure is used for all contacts although this is not fully correct for the edge and vertex contacts. Besides the normal force a shear force occurs on a solid– fluid interface. In order to calculate the shear force, a virtual dummy particle is created by projecting the SPH particle perpendicular onto the surface of the DEM particle. Due to the no slip condition, the virtual particle gets the local velocity of the DEM particle at this point. For mass and density, the values from the considered SPH particle are used. Then the shear force is obtained by means of the viscous fluid Eqs. (18) and (19), and applied to the pair of virtual particle and SPH particle. The virtual particles are created and removed ‘on the fly’ during the contact force calculation for each single SPH particle which is in contact with a polyhedron. The virtual particle is not stored and is not part of the neighbour search algorithm. No internal particle lists have to be altered for it. This means that for one SPH particle only the virtual particle created by itself is part of the summation interpolant of this particle. All virtual particles from other contact pairs are not part of it (see Fig. 9). To correct this, we adjust the normalizing constant [see eq. (9)] of the virtual particle, as proposed in [25], to

Comp. Part. Mech.

2 , 3hl 2

10

(38)

8

which is derived by the integration of the one-dimensional kernel in combination with the projected area of the SPH particle.

force [µN]

σV P =

s = 10µm, l = 1µm s = 5µm, l = 1µm s = 2µm, l = 1µm s = 1µm, l = 1µm s = 10µm, l = 2µm theoretical

6 4 2

5 Validation of the simulation method

0

5.1 Pure Couette flow

1

2 3 time [µs]

4

5

Fig. 10 Couette flow between two plates: force on fixed plate versus time 5 4 velocity [m/s]

The presented method was implemented in the DEM code LIGGGHTS [13]. In the following, it is validated by three simple examples. The first one is a Couette flow, where the relative velocity of a moving plate made of polyhedrons leads to a shear flow of SPH particles. The second one shows an icosahedron inside a Couette flow. For both the material, parameters of Polyethylene glycol PEG300 with a density of ρ = 1.13 g/cm3 and a dynamic viscosity of η = 90 mPas are used. The model parameters are h = 1.5l, χ = 1 and B = 0.01 GPa. All solid parts, like e.g. walls and moving plates, are composed of polyhedral particles which use the interaction forces explained in Sect. 4.3.

0

s = 10µm, l = 1µm s = 5µm, l = 1µm s = 2µm, l = 1µm s = 1µm, l = 1µm s = 10µm, l = 2µm theoretical

3 2 1 0

0

2

4 6 height [µm]

8

10

Fig. 11 Couette flow between two plates: velocity profile after 5 µs

In the first example, a pure Couette shear flow between a fixed and a moving plate is simulated. In this case, a linear velocity profile is expected and the force F on the plate is v F = ηA . s

(39)

A, v and s are the surface area of the plate, the velocity of the moving plate and the distance between the plates, respectively. Simulations are performed with SPH particles of size l = 1 µm and a gap s of 1, 2, 5, and 10 µm as well as a particle size l = 2 µm and a gap s = 10 µm. In all simulations, a ratio v/s = 0.5 µs−1 and a surface area of A = 200 µm2 are used. The plates are composed of cubical-shaped polyhedral particles. While these assembled walls are moved with a constant speed, the SPH particles are accelerated due to DEM–SPH interaction forces until they have adopted the velocity of the wall particles and the typical shear profile is formed. Figure 10 shows the force on the fixed plate. The simulations show good agreement with the analytical value of the force F = 9.0 µN. The range of the results is between 5 and 10 % above the analytical value, depending on the number of the SPH particle layers used to fill the gap. The simulation with s = 1 µm shows the largest difference between simulation and theoretical value because in this simulation only one layer of SPH particles is placed between the plates. Thus, the volume integral cannot be solved properly. From Fig. 11, it can be seen that the corresponding velocity pro-

file between the plates agrees well with the theoretical linear velocity profile. 5.2 Particle in a Couette shear flow In the second example, a solid particle is placed in the centre of a shear flow, see Fig. 12. The distance between the plates is s = 20 µm and the velocity of the moving plate is vw = 10 m/s, which leads to a shear rate γ˙ = vw /s = 0.5 µ s −1 . The DEM particle approximates the shape of a sphere and is a icosahedron with an outer radius of R = 2 µm. The size of the SPH particles is l = 1 µm. This value is based on the rough approximation that the spatial resolution of the SPH particles should be twice the size of the DEM particles [29]. Thus, we choose a resolution, which should be sufficient for the equivalent spherical particle. To take the finer substructure of the icosahedron into account, this value has to be set much smaller, which includes higher computational costs. As in the previous example, all walls are assembled by cubical-shaped polyhedral DEM particles. Figure 13 shows the total force on the particle for the case that the particle is fixed in space. The torque acting from the fluid flow on the particle is shown in Fig. 14. For comparison, the analytical solutions for force (40) and torque (41) on a spherical particle in a shear flow according to Stokes law are used.

123

Comp. Part. Mech.

vw

6

velocity [m/s]

5

s/2

2

0

0

1

2 3 time [µs]

4

5

Fig. 15 Velocity of a free particle over time and theoretical local fluid velocity in a Couette shear flow

Fig. 12 Fixed DEM particle in an SPH particle shear flow

0,8

20 DEM-SPH outer sphere vol. eq. sphere inner sphere

15 10 5 0

1

2 3 time [µs]

4

Fig. 14 Torque on a fixed particle over time. For comparison the analytical solutions according to Stokes law are shown for different spheres (see previous force plot)

F = 6π η Rvr el

(40)

M = 8π η R (γ˙ − ω).

(41)

With vr el being the relative velocity between particle and fluid, γ˙ is the shear rate of the fluid flow and ω the angular velocity of the particle. Because the polyhedral particle is only a rough approximation to a sphere, three types of ideal

123

0,7 0,6 0,5 DEM-SPH theoretical

0,4 0,3 0,2 0,1 0,0

5

Fig. 13 Total force on a fixed particle over time together with analytical solutions for different spheres. ‘inner’ denotes a sphere which fits inside the polyhedron and ‘outer’ a sphere which encloses the polyhedron. ‘vol. eq. sphere’ is a sphere with the equivalent volume of the polyhedron

3

angular velocity [µs−1 ]

25

force [µN]

DEM-SPH theoretical

3

1

s/2

0

4

0

1

2 3 time [µs]

4

5

Fig. 16 Angular velocity of free particle over time and theoretical fluid shear rate in simple shear flow

spheres are used for comparison: First, a sphere with a diameter that suits exactly inside the polyhedron (inner sphere); second a sphere which includes the polyhedron (outer sphere) and third a sphere with equivalent volume of the polyhedron. While the steady state drag force in Fig. 13 matches the analytical solution of a volume equivalent sphere within the range of 5 %, the torque (see Fig. 14) is about 30 % below the analytical value. We assume this is due to perturbations caused by the periodic boundary conditions. The reason for the scatter observed in Figs. 13 and 14 is presumably a result of a too coarse spatial SPH resolution or is caused by the generalized pressure calculation of the edge and vertex contacts.

5.3 Free particle in a Couette flow If the particle is free to move, the velocity approaches the theoretical local fluid velocity of a simple shear flow, which is in our case v = 5.0 m/s (see Fig. 15), and the angular velocity approaches the theoretical shear rate of the shear flow, which is γ˙ = 0.5 µs−1 (see Fig. 16).

Comp. Part. Mech.

6 Conclusion The coupling of a discrete element method (DEM) with polyhedral particles and smoothed particle hydrodynamics (SPH) was presented and exemplified. Because both methods are particle based and the underlying integration scheme is quite similar, the two simulation techniques are directly coupled in both ways through interaction forces between the solid DEM particles and the fluid SPH particles. The presented method was implemented into the opensource DEM software LIGGGHTS that has been previously extended by polyhedral particles [24]. The coupled method provides the possibility to simulate the behaviour of particle systems with many discrete particles inside a fluid flow. The method is validated by means of three examples for which exact solutions are known. The results show that the movement of individual DEM particles inside a flow and the flow field around individual particles is modelled with sufficient accuracy. This allows for the simulation of mixtures of fluid and sharp-edged solid particle systems. Such applications are important in many technical processes e.g. wire sawing, grinding or lapping. Acknowledgments This work was funded by the German Research Foundation in the project KU 929/19-2.

References 1. Adami S, Hu XY, Adams NA (2012) A generalized wall boundary condition for smoothed particle hydrodynamics. J Comput Phys 231(21):7057–7075 2. Bierwisch C, Kübler R, Kleer G, Moseler M (2011) Modelling of contact regimes in wire sawing with dissipative particle dynamics. Philos Trans R Soc A 369(1945):2422–2430 3. Boon CW, Houlsby GT, Utili S (2012) A new algorithm for contact detection between convex polygonal and polyhedral particles in the discrete element method. Comput Geotech 44:73–82 4. Cleary PW, Monaghan JJ (1999) Conduction modelling using smoothed particle hydrodynamics. J Comput Phys 148(1):227–264 5. Cleary PW, Morrison RD (2012) Prediction of 3d slurry flow within the grinding chamber and discharge from a pilot scale sag mill. Miner Eng 39:184–195 6. Cundall PA (1988) Formulation of a three-dimensional distinct element model–part i. a scheme to detect and represent contacts in a system composed of many polyhedral blocks. Int J Rock Mech Min Sci Geomech Abstr 25(3):107–116 7. Cundall PA, Strack OD (1979) Discrete numerical model for granular assemblies. Int J Rock Mech Min Sci Geomech Abstr 16(4):47–65 8. Galindo-Torres SA (2013) A coupled discrete element lattice boltzmann method for the simulation of fluid-solid interaction with particles of general shapes. Comput Methods Appl Mech Eng 265:107–119 9. Gingold RA, Monaghan JJ (1977) Smoothed particle hydrodynamics—theory and application to non-spherical stars. Mon Not R Astron Soc 181:375–389 10. Grof Z, Cook J, Lawrence CJ, Štˇepánek F (2009) The interaction between small clusters of cohesive particles and laminar flow: coupled dem/cfd approach. J Petrol Sci Eng 66(1–2):24–32

11. Höfler C (2013) Entwicklung eines Smoothed Particle Hydrodynamics (SPH) Codes zur numerischen Vorhersage des Primärzerfalls an Brennstoffeinspritzdüsen. PhD thesis, Karlsruher Institut für Technologie 12. Huang YJ, Nydal OJ (2012) Coupling of discrete-element method and smoothed particle hydrodynamics for liquid-solid flows. Theor Appl Mech Lett 2(1):012002 13. Kloss C, Goniva C, Hager A (2012) Models, algorithms and validation for opensource dem and cfd-dem. Prog Comput Fluid Dyn Int J 12(2/3):140–152 14. Libersky LD, Petschek AG, Carney TC, Hipp JR, Allahdadi FA (1993) High strain lagrangian hydrodynamics: a three-dimensional SPH code for dynamic material response. J Comput Phys 109(1):67–75 15. Lu G, Third JR, Müller CR (2015) Discrete element models for non-spherical particle systems: from theoretical developments to applications. Chem Eng Sci 127:425–465 16. Lucy LB (1977) A numerical approach to the testing of the fission hypothesis. Astron J 82(12):1013–1024 17. Monaghan JJ (1992) Smoothed particle hydrodynamics. Ann Rev Astron Astrophys 30(1):543–574 18. Monaghan JJ (1994) Simulating free surface flows with SPH. J Comput Phys 110(2):399–406 19. Monaghan JJ (2000) Sph without a tensile instability. J Comput Phys 159(2):290–311 20. Monaghan JJ, Kajtar JB (2009) SPH particle boundary forces for arbitrary boundaries. Comput Phys Commun 180(10):1811–1820 21. Morris JP, Fox PJ, Zhu Y (1997) Modeling low reynolds number incompressible flows using SPH. J Comput Phys 136(1):214–226 22. Müller M, Schirm S, Teschner M, Heidelberger B, Gross M (2004) Interaction of fluids with deformable solids. Comput Anim Virtual Worlds 15(3–4):159–171 23. Nassauer B, Kuna M (2013) Contact forces of polyhedral particles in discrete element method. Granul Matter 15(3):349–355 24. Nassauer B, Liedke T, Kuna M (2013) Polyhedral particles for the discrete element method: geometry representation, contact detection and particle generation. Granul Matter 15(1):85–93 25. Nassauer B (2015) Entwicklung einer 3D Diskrete Elemente Methode mit polyederförmigen Partikeln zur Simulation der fluidgekoppelten Prozesse beim Drahtsägen. PhD thesis, TU Bergakademie Freiberg, Freiberg 26. De Pellegrin DV, Stachowiak GW (2005) Simulation of threedimensional abrasive particles. Wear 258(1–4):208–216 Second international conference on erosive and abrasive wear 27. Potapov AV, Hunt ML, Campbell CS (2001) Liquid-solid flows using smoothed particle hydrodynamics and the discrete element method. Powder Technol 116(2–3):204–213 28. Qiu L-C (2013) Numerical modeling of liquid-particle flows by combining SPH and DEM. Ind Eng Chem Res 52(33):11313– 11318 29. Robinson M, Ramaioli M, Luding S (2014) Fluid-particle flow simulations using two-way-coupled mesoscale SPH-DEM and validation. Int J Multiph Flow 59:121–134 30. Rozmanov D, Kusalik PG (2010) Robust rotational-velocity-verlet integration methods. Phys Rev 81(5):056706 31. Shao S, Lo EYM (2003) Incompressible sph method for simulating newtonian and non-newtonian flows with a free surface. Adv Water Resour 26(7):787–800 32. Takeda H, Miyama SM, Sekiya M (1994) Numerical simulation of viscous flow by smoothed particle hydrodynamics. Prog Theoret Phys 92(5):939–960 33. Tao H, Jin B, Zhong W, Wang X, Ren B, Zhang Y, Xiao R (2010) Discrete element method modeling of non-spherical granular flow in rectangular hopper. Chem Eng Process 49(2):151–158 34. Valizadeh Alireza, Monaghan Joseph J (2015) A study of solid wall models for weakly compressible SPH. J Comput Phys 300:5–19

123

Comp. Part. Mech. 35. van Zuijlen AH, Bijl H (2005) Implicit and explicit higher order time integration schemes for structural dynamics and fluidstructure interaction computations. Comput Struct 83(2–3):93–105 Advances in Analysis of Fluid Structure InteractionAdvances in Analysis of Fluid Structure Interaction

123

36. Zhu YI, Fox PJ, Morris JP (1999) A pore-scale numerical model for flow through porous media. Int J Numer Anal Meth Geomech 23(9):881–904