A generalised multi-scale Peridynamics–DEM framework and its application to rigid–soft particle mixtures

The discrete element method (DEM) is the most dominant method for the numerical prediction of dynamic behaviour at grain or particle scale. Nevertheless, due to its discontinuous nature, the DEM is inherently unable to describe microscopic features of individual bodies which can be considered as continuous bodies. To incorporate microscopic features, efficient numerical coupling of the DEM with a continuous method is generally necessary. Thus, a generalised multi-scale PD–DEM framework is developed in this work. In the developed framework, meshfree discretised Peridynamics (PD) is used to describe intra-particle forces within bodies to capture microscopic features. The inter-particle forces of rigid bodies are defined by the DEM whereas a hybrid approach is applied at the PD–DEM interface. In addition, a staggered multi-scale time integration scheme is formulated to allow for an efficient numerical treatment of both methods. Validation examples are presented and the applicability of the developed framework to capture the characteristics mixtures with rigid and deformable bodies is shown.


Introduction
In the last decades, the computational mechanics community has gradually grown and various numerical methods have been developed. Thanks to the significant increase in computational resources and efficiency it is now possible to tackle highly complex problems. This involves the coupling of different numerical methods to overcome the drawbacks of individual methods, allowing to capture more complex phenomena. In this work a three-dimensional framework for the efficient numerical treatment of coupled continuous and discontinuous material behaviour is developed. The motivation behind this is the ability to describe the discontinuous dynamics of a particle system, whilst capturing the microscopic features of individual particles by a continuous description. A possible application is the simulation of grain collisions, in which microscopic phenomena of individual grains such as deformability or fragmentation are considered. Real-life examples include the compaction of powders [1], the compaction of powder mixtures [2], the identification of fracture origins in ceramics [3] and the characterisation of particle mixtures with hard and soft grains [4].
In the following, suitable numerical methods have to be picked to develop the desired coupling framework. The most widespread and well-established method applied in computational mechanics is the Finite Element Method (FEM) (e.g. [5][6][7]), which is based on the consideration of continuous media. Within the method, the media is discretised by finite elements and associated nodal points. The degrees of freedom, e.g. displacements or temperature, are only defined and solved for these discrete points within the continua. Up to now, the FEM is the method of choice for most standard solid mechanics applications. However, there are various problems for which the FEM is not applicable. One crucial shortcoming of the classical FEM is its inability to capture fracture and crack propagation. To con-sider this phenomenon the so-called eXtended Finite Element Method (XFEM) (cf. [8,9]), can be applied. It is based on the enrichment of ansatzfunctions with discontinuous functions within the FEM. Besides the additional numerical cost, the main disadvantage is that discontinuities are not naturally considered within the method. To overcome this problem Peridynamics (PD) was introduced as an alternative total Lagrangian formulation to classical continuum mechanics [10]. In contrast to classical methods, the underlying equations are integro-differential equations without spatial derivatives. Consequently, discontinuities in space are naturally considered within the non-local framework, even though the underlying integral equations describe continuous media. Moreover, the integro-differential equations allow for a direct meshfree discretisation using nodal integration. Up to now, a major shortcoming in PD is the treatment of contact. Initially, short-range force were applied to account for the contact between peridynamic bodies [11]. In [12] an overview of various possibilities to tread contact within PD is given. Additionally, a complex peridynamic specific contact model is introduced, conserving angular momentum during collision.
When it comes to the description of the dynamics of soil and granular materials the Discrete Element Method (DEM), introduced in [13], is commonly applied. Within the DEM, the underlying dynamics are captured on grain level on the basis of contact forces, formulated with respect to micromechanical parameters. Generally, the solid or granular material is represented by rigid particles with an associated particle size distribution. It is also possible to capture microscopic effects on grain level using the DEM. Examples are the modelling of progressive failure in fractured rock masses [14], the approximation of grain deformations under the consideration of compressible effects using the implicitly formulated deformable DEM [15] and the bonded sphere approach in which the microscopic behaviour of a rubber grain is described by a deformable agglomerate of rigid particles being able to move relative to each other [4]. Note that to describe the microscopic behaviour, these approaches are still based on discontinuous approaches relying on micro-mechanical parameters. Generally, it is numerically advantageous to describe grains as continua to capture their microscopic behaviour. The reason behind this is that the associated material models of continuum methods are formulated with respect to measurable material parameters, having a real physical meaning. In contrast to macroscopic material parameters, the micro-mechanical parameters represent fitting parameters which are not size independent and can depend on particle size distribution and particle shapes. Consequently, a problem-dependent complex numerical calibration on the basis of experimental measurements, e.g. triaxial compression tests, is always necessary (see e.g. [16]). In contrast, classical macroscopic material parameters, e.g. elastic, viscous and plastic properties, can be directly determined by standard experimental tests without the necessity of complex calibrations.
The most promising combination to capture the discontinuous dynamics of grains, whilst taking into account microscopic effects on the grain level, is the combination of DEM and PD. Firstly, the DEM is perfectly capable to capture the dynamics on the grain level and various contact models have been developed in the last decades. Secondly, PD allows to describe individual grains as continua and to capture their microscopic behaviour by applying peridynamic material models. Besides elasticity models, this includes fracture models as shown in [17][18][19], amongst others. An application of peridynamic fracture models for grain crushing is further considered in [20,21]. Thus, a PD-DEM coupling framework is superior to FEM-DEM coupling frameworks (see e.g. [22,23]) and to the 'meshfree numerical tool' developed for the simulation of mixtures of hard and soft grains (see e.g. [24,25]). The reason behind is that fracture is naturally included in the PD formulation and has been intensively studied (see e.g. [18,26]). Consequently, the discontinuouscontinuous approach is applicable for a wider range of applications. A different continuum approach is presented in [27], considering flexible DEM particles of arbitrary polyhedral shapes on the basis of the Virtual Element Method.
Recently a PD-DEM framework for the prediction of fracture of colliding grains in two-dimensional space has been implemented [28]. In the framework, the intra-particle forces within the arbitrarily shaped grains are computed on the basis of a peridynamic formulation. The inter-particle interactions, i.e. the forces between two grains coming into contact, are computed by DEM-like contact laws between particles of distinct discretised bodies. Moreover, in [29] the coupling of PD and DEM is considered from a computational point of view with respect to the software library ParticLS, in which meshfree methods and the DEM are considered. A shortcoming of the existing framework in [28], which is tackled in this work, is that all grains are treated as peridynamic bodies. However, it is not necessary to consider all grains on the microscopic scale for various applications. An example are deformablerigid mixtures of soft and hard grains, where all rigid grains can be treated as DEM grains, i.e. as a single discrete particle. Incorporating rigid DEM grains as well as microscopic phenomena in PD grains within the same framework constitutes a multi-scale approach. Consequently, it is desirable to have the possibility of a multi-scale time integration scheme. In the following, the term grain is replaced by body to allow the description of arbitrary objects in a generalised framework. Moreover, in discretised form bodies consists of particles.
This contribution targets an efficient numerical coupling of discontinuous and continuous material behaviour by developing a generalised multi-scale PD-DEM coupling framework. The fundamentals of the PD-DEM framework, i.e. conservation principles, PD and DEM, are recapitulated in Sect. 2. The proposed generalised multi-scale PD-DEM coupling framework is then introduced in Sect. 3. It is focused on the generalised force coupling as well as on the staggered solution scheme for the multi-scale time integration. In Sect. 4.1, the PD-DEM coupling is verified using a 3D Hertzian contact problem and the multi-scale time integration is verified in Sect. 4.2. Furthermore, numerical examples of coupled discontinuous-continuous material behaviour are presented in Sect. 5 in terms of simultaneous consideration of deformable and rigid bodies. Finally, the main findings are summarised and possible extensions of the current work are discussed in Sect. 6.

Conservation principles in mechanics
In numerical frameworks, it is crucial to obey fundamental laws of physics with respect to associated physical values to allow reliable predictions. In the case of dynamic systems, the fundamental laws are the conservation of linear and angular momentum.
To conserve linear momentum it is necessary to fulfil Newton's second law as well as Newton's third law f A = −f B . Newton's second law states that the resulting force f is directly proportional to the change of linear momentum p, i.e for a system with constant mass m, proportional to the change of velocity v and thus proportional to the second time derivative of the positional vector u which is the acceleration a. Newton's third law states that the force between two points A and B in contact is equal in magnitude with opposed direction. When fulfilling Newton's third law, the accumulated point-wise intra-particle forces within a generalised body B with volume v and boundary ∂B with area a vanish. Consequently, only traction forces t with associated normal direction n on ∂B have to be considered for the fulfillment of Eq. (1). Thus, including body forces m i b i , where b is the specific body force, it yields for a discontinuous particle system. For continuous bodies with infinitesimal particle volumes it yields (3) and locally where ρ is the density in the current configuration and divσ the divergence of the Cauchy stress tensor.
To conserve angular momentum L, it is necessary to fulfil where I is the moment of inertia andω the rate of angular velocity. Moreover, M act is the sum of acting moments, either directly applied as external moments or generated by interaction forces (e.g. contact forces or rolling resistance). Thus, Eq. (5) states that acting moments lead to a change of angular momentum.
In the case of the classical continuum description (Eq. 4), angular momentum is automatically conserved when the stress tensor is symmetric, i.e. σ = σ T . In contrast, for discontinuous particle systems it is required to consider the angular velocities in the degrees of freedom and to solve Eq. (5) with respect to the acting moments.

Peridynamics
In the theory of Peridynamics (PD), firstly introduced in [10], non-local particle interactions over a specific radius within the family H are considered. Considering a master particle I , its family represents the domain of influence and contains all neighbouring particles J whose distances are less equal than the horizon size δ, as depicted in Fig. 1.
In contrast to the updated Lagrangian type smoothedparticle hydrodynamics (SPH) method (cf. [30]), the original framework of PD is of total Lagrangian type. As a consequence, the neighbourhoods of particles do not change during the computation. Thus, scalars corresponding to material properties in the governing equations are defined with respect to the reference configuration and vectorial quantities are defined with respect to initial positions X and time t. Fig. 1 Family H of particle I with interacting neighbouring particles J Compared to classical continuum mechanics, the resulting equations in PD are integro-differential equations without spatial derivatives. In the context of the preservation of linear momentum, this leads to the replacement of the classical divergence of the stress tensor divσ , cf. Eq. (4), by its peridynamic counterpart L. Thus, the operator L represents the resulting force density from a peridynamic material model and the associated equation of motion is defined by Consequently, the peridynamic force is defined with respect to the volume V of the initial configuration by Similar to classical continuum mechanics, the conservation of angular momentum is directly handled via the material model. Thus, a fundamental requirement for peridynamic material models is that no angular momentum is generated due to deformation. To define peridynamic material models, the deformations of bonds ξ = X − X, defined as the vectors between the initial position of master particle X and all initial positions X of particles within its family H X , are considered. Generally, it is distinguished between bondbased and state-based PD, whereby the difference lays in the bond-force computation, cf. [31]. In bond-based models, the bond-force between two particles depends on the deformation of the associated bond only, whereas the bond-force between two particles depends on the collective deformation of bonds within the family for state-based peridynamic models. In the following, the force density computation for bondbased as well as state-based models is explained and specific elastic material models are introduced. Bond-based models are derived from a central potential and the force densities are computed with respect to a pairwise force function t, which generally depends on the relative displacement (displacement field w) and the bond vector itself. Thus, it yields Note, that bond-based models are automatically restricted to a Poisson's ratio of 0.25 within three-dimensional approaches due to the derivation from a central potential. Thus, there are generally not applicable for the modelling of nearly incompressible materials like rubber. The force vector between two particles is always parallel to its deformed bond vector η, i.e. the bond force densities are represented by where t is the scalar bond-force. Consequently, angular momentum is automatically preserved. In the further course, the micro-elastic brittle material model is utilised and the scalar bond-force is defined by Within Eq. (10) s = ξ +η − ξ ξ is the bond stretch, c = 18K πδ 4 a spring constant and φ a damage function. The spring constant can be defined with respect to the classical compression modulus K and is obtained by postulating the same stress power for the linear elastic continuum mechanical material model and the presented peridynamic model. Note that the correspondence between micro and macro material parameters is formulated with the assumption of compact support. Thus, the assumption is violated for particles whose initial distances to the surface are smaller than δ. In the scope of this work, no surface correction is applied. However, a study about different surface correction procedures is conducted in [32]. Moreover, in this work bond breakage is not considered and φ = 1 is used. Similar to [33,34], Eq. (6) is discretised in a meshfree fashion. Neglecting the body force term, it yields on particle level Thus, the peridynamic forces exerted on particle I due to its interactions within its family are defined by Note that the result of the meshfree discretised PD is a particle system. However, the system is describing a continuum body in a non-local manner.
As the name indicates, state-based PD are formulated with respect to states. From a mathematical point of view, states are generalised second-order tensors representing the mapping of bond vectors ξ ∈ H X to either a scalar or a vector. In the following, states are indicated by an underscore. The equation of motion in state-based PD is formulated with respect to force state densities T. This leads to the general state-based expression of for the non-local operator L. Applying nodal integration and multiplying by volume leads to the force definition on particle level Within the discretised version the force state density T I contains the collective information about forces exerted on/from particle I within its family. Consequently, T I ξ I J maps the bond vector ξ I J into a force vector per unit volume square. In the formulation of associated material laws the deformation state, which maps the bond vector into the deformed bond vector, is defined by Y[X, t] X − X = η.
In the following, Linear Peridynamic Solid (LPS) [31] will be used. In order to define the force state density, the scalar reference state the scalar deformation state as well as the deformed direction state are introduced. In this way, the bond extensions are defined with respect to the scalar extension state by Moreover, the extension state is split into an isotropic e i and a deviatoric part e d . For their definitions it is necessary to define the scalar weighted volume and the peridynamic dilatation with respect to the the weighting influence function w, which is further assumed to be 1. Note that the peridynamic dot product between two states is a scalar and is defined by Using Eqs. (18) and (20), the deviatoric scalar extension state is defined by Finally, the force state densities are defined with respect to the deformed direction state (Eq. 17) and the scalar force state t(Y) by Note that angular momentum is automatically conserved using the LPS formulation. The scalar force state in Eq. (23) consists of a co-isotropic and co-deviatoric part, i.e. t = t i + t d , whereby the parts are defined with respect to the peridynamic free energy function by The material parameter α is the micro shear modulus which can be related to the macroscopic shear modulus G by α = 15G m . Note that similar to the bond-based approach, the missing compact support leads to inaccurate α's near the surface. An extension of the described model to elasto-plasticity with von Mises plasticity is shown in [35].
In this work, fractions of the so-called PD correspondence formulation, which is based on a non-local stress tensor, is used. Based on the discretised shape tensor the non-local deformation gradient is defined by Based on this non-local deformation gradient, a classical material model can be utilised for the stress computation. The only requirement is the application of the same macro-mechanical material parameters. Thus, the material parameters of the bond-based as well as state-based elastic model needs to be transformed into the Lame constants λ and G. Applying the compressible Neo-Hookean free energy function with the compressible part of Ciarlet [36] W the second Piola-Kirchhoff stress tensor is computed by [37] where C = F T F is the right Cauchy-Green tensor, I C its first invariant and J = det F the Jacobean. As a result, the mechanically more meaningful non-local Cauchy-stresses are obtained by Moreover, in engineering practise the Green-Lagrange strain tensor E = 1 2 (C − 1) is commonly used for the evaluation of finite strains.

Discrete element method
The second method considered within the developed coupling scheme is the Discrete Element Method (DEM). In contrast to the continuum based PD, in the DEM discrete bodies are considered and thus, no discretisation is necessary.
In the following, spherical bodies are considered, whilst generally arbitrary shaped discrete elements can be used. Their equation of motion is defined with respect to Eq. (1) by where mb are body forces and f DEM inter-particle DEM forces. The DEM forces are classically composed of normal and shear forces of adjacent discrete elements being in contact via whilst the direct neighbours of bodies are progressively updated. The used rheological model for the DEM-force computation between two spherical bodies is depicted in Fig. 3. It consists of a divider, a spring in normal direction with normal stiffness k n , a spring in tangential direction with shear stiffness k s and a slider. In contrast to other numerical methods, the interpenetration of bodies is not prohibited within the DEM and the resulting overlap represents the basis for the force computation. The resulting forces of two overlapping spherical bodies A and B is visualised in Fig. 2. The normal forces are defined with respect to the overlap in normal direction d n and normal stiffness k n . The incremental shear forces between two colliding bodies A and B are defined by whereu s is the relative shear velocity and k s the shear stiffness. Note that the normal and shear stiffness can both be a function of the overlap in normal direction d n . The total shear force is limited by a Coulomb-like slip model, i.e.
where μ is the friction coefficient. The divider in Fig. 3 indicates the absence of forces when two adjunct particles do not overlap.
The no slip non-linear Hertz-Mindlin (HM) contact model, cf. [38], is used in this work where contact stiffnesses are defined as: with respect to the equivalent sphere radiusR = (1 The associated shear moduli G A c and G B c are computed by Besides contacts between two spherical bodies, contacts between spheres and kinematically constraint rigid walls are considered. The associated contact forces are defined analogously to the sphere-sphere interaction. Since walls do not possess a radius, their radii are assumed to be equivalent to the radii of their spherical contacts for the computation of the equivalent sphere radiusR. As mentioned in Sect. 2.1, the conservation of angular momentum is not automatically fulfilled in the DEM and it is required to consider the angular velocities as degrees of freedom. In this work no external moments are considered and thus, the acting moments M act result solely from interaction forces, cf. Fig. 2. As a result, it yields with respect to Eq. (5) where R c represents the radial vector from spheres origin to the contact point. Consequently, besides the equation of  motion, the first-order differential equation has to be solved for all DEM bodies as well.

Generalised force coupling
A generalised PD-DEM force coupling scheme is introduced based on the framework discussed in Sect. 2. All possible contact variations on the basis of three PD and three DEM bodies are depicted in Fig. 4. The associated bodies are colour-coded in green (DEM) and grey (PD), whilst the centers of discrete elements are depicted by squares and the discretised peridynamic particles by crosses. Moreover, hybrid particles acting as peridynamic particle as well as DEM body are highlighted in red. In the following, DEM bodies are also denoted by discrete particles. The inter-particle forces between two discrete elements (DEM-DEM) are obtained by the DEM (Eq. 32). The macroscopic intra-particle forces within a peridynamic body f PD are defined by Eqs. (12) and (14) for a bond-based and statebased peridynamic model, respectively.
The key of the developed formulation is the introduction of hybrid particles at the interface of PD and DEM. Indeed, all surface particles of the peridynamic bodies are also treated as discrete elements. This allows the computation of inter-particle forces between discrete elements and the surface of peridynamic bodies (PD-DEM), inter-particle forces between two peridynamic bodies (PD-PD), as well as self-contact of peridynamic bodies by DEM contact forces. An example where self-contact of a peridynamic body can take place is the compression of a highly deformable hollow sphere under compression when the top and bottom boundary of the hollow sphere come into contact. The resulting forces are further denoted by f Coupling and are computed by Eq. (32). A problem in the definition of hybrid particles is that the underlying peridynamic particles do not inherently have a radius. They represent integration points with associated volumes as integration weights accounting for the fraction of continuous volume they are representing. Therefore, the quasi radii R of the hybrid particles required for contact detection and force computation are defined on the basis of the associated volumes V of the integration points via Hence, the required DEM part of hybrid particles (i.e. contact forces) can be calculated. Note that the DEM contact model for f DEM and f Coupling can be different.
Since hybrid particles only represent integration points with their volume as weighting within the peridynamic framework, they do not possess the rotational degree of freedom of DEM bodies. Thus, the quasi radii R of hybrid particles are only used for contact detection and contact force computation, whereby the resulting contact forces are applied on the associated peridynamic integration points. Moreover, self-contact in peridynamic bodies is not considered for par-ticles within the same family. In the proposed formulation the drawback of the lacking kronecker-delta property resulting from the meshfree peridynamic discretisation diminished. Due to the hybrid approach at the surface of peridynamic bodies, boundaries are indirectly applied by the resulting coupling forces.
Summarising, the forces of the particle system are defined as the superposition of intra-particle peridynamic, interparticle DEM, inter-particle coupling and body forces by

Multi-scale time integration
The state-of-the-art for PD as well as DEM is to perform an inherently conditionally stable explicit time integration for the equation of motion. A detailed investigation of different stability criteria for classical peridynamic models is done in [39]. In this work, the CFL-criterion (cf. [40]) where h is the characteristic length and c the wave speed, defined by is applied. As discussed in [41], it is not obvious what the characteristic length within the peridynamic framework is. However, it is said that the time step estimation by the CFLcriterion is conservative when the particle spacing is assumed to be the characteristic length. The critical time step t DEM crit of the DEM is defined on the basis of the equivalent per-particle stiffnessk which is composed of the associated stiffnesses of all contacts (see e.g. [42]). Thus, each particle has an individual critical time step and the global critical time step is generally defined as the minimum of local critical time steps, i.e.
The reader is referred to [42] for the computation of the equivalent per-particle stiffnessk. An extension to include rotational degrees of freedom is presented in [43]. = v n− 1 2 + a n t u n+1 = u n + v n+ 1 2 t + a n t 2 .
Moreover, the rotational accelerations of pure DEM bodies, cf. Eq. (38), have to be integrated to fulfil the conservation of angular momentum. Thus, the same time integration as for the translational kinematics, cf. Eq. (44), is applied and it yields Based on the simultaneous use of PD and DEM with associated critical time steps defined in Eq. (41) and Eq. (43), two distinct critical time steps have to be considered within simulations. The straightforward approach is to define the critical time step as the minimum of both, i.e. t crit = min( t PD crit , t DEM crit ), compute all forces and update the kinematics monolithically using Eq. (44). However, the use of a single time step is computationally highly inefficient. The main issue is that the length scales of PD and DEM may vary in magnitudes due to the targeted multiscale approach. This may result in a significantly smaller critical time step t PD crit . This is why a staggered integration scheme is developed. Based on Eqs. (41) and (43), the maximal integral multiple m between t PD crit and t DEM crit or vise versa is determined. Even though t PD crit < t DEM crit is expected, both possibilities are considered to formulate the generalised multi-scale time integration scheme.
The key idea of the implemented staggered solution scheme is the successive update of kinematics resulting from peridynamic and DEM forces with respect to their associated critical time steps. The simplified update scheme for equivalent critical time steps is illustrated in Eq. (46): Based on the updated kinematics of the last time step u n , the inter-particle forces f DEM n+1 and f for discrete and hybrid particles. Afterwards, the peridynamic forces f PD n+1 are computed and the associated kinematics are updated. Thus, the hybrid particles are updated a second time, but now with respect to the intra-particle forces. Note, that it is essential to consider body forces of hybrid particles only once during the staggered time integration. Thus, they are always taken into account at the same time as the coupling forces. Consequently, body forces of peridynamic surface particles being in contact, i.e. they are hybrid particles, are always considered during the DEM time integration step. In order to simplify the treatment of body forces within the staggered time integration, the body forces of peridynamic particles are also considered in the first step.
For the staggered time integration of hybrid particles and t DEM = t PD this yields to v Coupling Note that the mean velocitiesV in the velocity updates are always defined with respect to their values of the last incremental update. Thus, for varying peridynamic and DEM time step sizes the symplecticity of the leapfrog integration is violated for the kinematic update of hybrid particles. The implementation of the generalised multi-scale time integration algorithm is explained in detail in Sect. 3.3.

Implementation
The open-source DEM framework Yade [44] is chosen as a basis for the implementation. This comes with a high variety of already implemented DEM models to choose from for DEM and hybrid force computation. In this paper, only the simple model introduced in Sect. 2.3 will be applied. An additional advantage of using Yade is that an optimised collider for contact detection as well as time integration schemes are already available. Model generation including definitions of DEM and peridynamic bodies with related discretisation are further handled via the available Python interface.
The peridynamic framework and the peridynamic material models described in Sect. 2.2 are implemented based on the pre-existing software architecture. Moreover, a coupling engine is implemented to account for the PD-DEM coupling formulation and the applied staggered time integration scheme outlined in Sect. 3.2. In the following, the extended simulation loop within Yade, cf. Algorithm 1, is explained first before focusing on the implementation of the multiscale time integration. Yade is run via a Python-interface and its classical simulation loop consists of the resetting of forces and the action of fundamental DEM engines. Note, that the fundamental DEM engines have to be defined in a specific order: approximate collision detection, exact collision detection with overlap computation, definition of physical properties of new interactions, DEM force computation via constitutive laws and time integration. The highlighted step Coupling Engine is the required extension for the generalised multi-scale PD-DEM coupling time integration. Since Yade is a DEM based framework, the global time step is set to t DEM crit and the iteration number n is incremented with respect to this time step. As described in Sect. 3.2, the maximum integral multiple m between both critical time steps is computed in the beginning and is further used in the simulation. Before the coupling engine is called, the classical DEM steps of neighbouring search and DEM force computation are performed. The kinematics of the associated particles are then updated, leading to new positions u DEM,Coupling n+1 for discrete and hybrid particles.
The pseudo code of the coupling engine is shown in Algorithm 2. It is distinguished between the two possible cases t PD crit >= t DEM crit and t PD crit < t DEM crit . In case of t PD crit >= t DEM crit it is checked if the current time t n is a m-mulitple of t PD crit . If this is the case, the peridynamic forces f PD n+1 are computed and the kinematics of corresponding particles are updated. The associated time integration is performed with respect to t PD crit = m t DEM crit . Otherwise, no time integration of peridynamic particles is performed.
The second and most probable case is t PD crit < t DEM crit . Similar to the first case, the time integration of peridynamic particles takes place in the coupling engine. However, an incremental time integration with t inc = t PD crit = t DEM m is necessary to account for the smaller peridynamic time step. Thus, the peridynamic forces and corresponding kinematic quantities are m times updated incrementally. After the last incremental update, t n+1 is reached.

PD-DEM coupling
In order to verify the implemented PD-DEM coupling scheme, the solution for the contact of a rigid sphere with an elastic half space obtained by Hertzian contact theory (cf. [45]) is considered. Dependent on the penetration depth d, whereĒ is the effective Young's modulus and R the radius of the sphere. Within the numerical PD-DEM approach, the rigid sphere is represented by a discrete DEM particle of radius R DEM = 5 cm while the elastic half-space is approximated by a peridynamic body of dimensions 45 cm × 45 cm × 22.5 cm, as illustrated in Fig. 5. The peridynamic body is discretised by 30 × 30 × 15 = 13500 particles based on a regular particle spacing of x = 1.5 cm. Using Eq. (39), the representative radius of the hybrid particles is R PD = 0.93 cm. The kinematics of surface/hybrid particles are constrained in associated normal directions whilst all internal peridynamic particles are unconstrained for the purpose of representing the elastic half space. An exception are the hybrid particles in the plane of the initial contact surface which are unconstrained as well.
To obtain the numerical force response, the described Hertz-Mindlin contact model with normal force definition in Eq. (33) is used for the DEM-like coupling forces between two bodies. Both the bond-based (Eq. 12) and state-based (Eq. 14) peridynamic models are used. The applied material parameters are listed in Table 1, whilst friction is not consid-ered in this example, i.e. μ = 0. In contrast to the analytical approach, the resulting coupling forces are DEM-like forces which generally depend on the micro-mechanical contact parameters between DEM bodies. When applying the Hertz-Mindlin contact model, the microscopic DEM parameters are represented by their macroscopic counterpart and the analytical Hertzian contact force is obtained for the contact of two single DEM bodies with associated radii R DEM A and R DEM B . Thus, the numerical contact forces are computed locally between spherical particles, whereas the analytical solution describes the global behaviour between a rigid sphere and an elastic half space. For this reason, a perfect match between predicted and analytical results is not expected. Nevertheless, the overall trend of the analytical result should be captured. It should be mentioned that a perfect match can be obtained by calibration, e.g. [16,46], however, this is not the aim of the current study.
In the simulations a prescribed velocity of v z = −1 cm s is applied for the discrete particle in the z-direction using a fixed time step of t = 10 −4 s . Thus, the penetration depth is indirectly used as control parameter and the acting normal forces are measured. The measurement is done incrementally by setting the velocity to zero at certain times, iterating until quasi-static conditions are reached and measuring the resulting normal force. This procedure is performed with an incremental penetration depth of d = 0.05 cm until d = 1 cm is reached. The comparison of analytical (green), PD bond-based (blue) and PD state-based (orange) normal force response with respect to the penetration depth is depicted in Fig. 6. Overall, the trend of the analytical Hertzian solution is captured well by the PD-DEM coupling framework applying the bond-based as well as state-based material model. In comparison to the analytical solution, the numerical contact force is underestimated for both the bond-based and the state-based PD material model. Moreover, the state-based solution is approaching the analytical solution with increasing penetration depth. Based on the fact that no calibration is performed for the micro-mechanical contact parameters, Hertzian contact is properly represented within the developed framework. Since the focus and goal of this contribution is the efficient numerical coupling of discontinuous and continuous material behaviour using the example of grain mixtures, it is relinquished to perform a calibration for the micro-mechanical contact parameters. Thus, a successful verification of the PD-DEM coupling framework has been performed and it is meaningful to use the framework for more complex applications.

Multi-scale time integration
In a next step the verification of the multi-scale time integration is performed. The verification example consists of   (Fig. 7). A velocity of v z = −1 cm s is continuously applied to all discrete particles until a compression of 0.2 cm is reached. Moreover, the displacements in the x-y plane are constraint for all peridynamic particles at the bottom. The material parameters listed in Table 1 are applied and the state-based PD formulation is used. Measured are the resulting normal forces on the wall with respect to the norm of applied displacements, u, of discrete particles in the z-direction. The normal forces are measurements under fully dynamic conditions (no damping) in order to investigate possible differences between monolithic and staggered time integration for the developed coupling approach.
In the following, all three different possibilities arising from the developed and implemented multi-scale time integration (Sects. 3.2 and 3.3) are considered. The first one is the classical monolithic solution of DEM and PD. The second possibility is the application of the staggered integration The resulting normal forces are plotted against the norm of enforced displacement of discrete particles in the z-direction in Fig. 8. Overall, the normal forces are oscillating due to the fully dynamic approach, whilst the peak values increase with respect to increasing displacements. On the one hand, there are no differences in normal forces between the staggered scheme with identical time steps compared to the staggered multi-scale scheme. Thus, even though the symplecticity of the integrator is violated, no negative effects are observable in the multi-scale time integration. On the other hand, there are differences between the results for the monolithic and the staggered scheme. The period length of the oscillations is slightly smaller for the staggered scheme and the magnitudes in normal forces are also lower. These effects result solely from the postulated successive integration of motion for hybrid particles and are expected. By applying the DEM like coupling forces on hybrid particles and updating their kinematics before the peridynamic forces are computed, a higher 'constraint' is applied on the peridynamic body. The reason behind is that the updated positions of hybrid particles are used to compute the peridynamic forces. Thus, the peridynamic force response is slightly higher than for the monolithic integration, leading to a reduced wave speed as well as to a reduced accumulated force within a total time step. In other words, the dynamic impact of discrete particles upon peridynamic bodies is slightly reduced in the staggered time integration scheme. Finally, the computational efficiency of the implemented multi-scale time integration is evaluated. Therefore, normalised computation times obtained on a single Intel Core i7-7560U (2.40 GHz) processor are compared. As a reference solution the monolithic approach with t = 10 −4 s is considered. The normalised computation times for the multiscale simulations with t PD = 10 −4 s and t DEM = m t PD with m = 1, 2, 4, 6 are evaluated and listed in Table 2. Overall, the normalised computation time decreases with increasing m. In case of equal time steps the staggered computation time is 27 % higher than the monolithic computation time. Equal normalised computation times are obtained for m = 2. A reduction of 16 % and 21 % is observed for the staggered multi-scale integration for m = 4 and m = 6, respectively. This proofs the computational advantage of the multi-scale time integration over the monolithic one. In the considered example the number of discrete particles is with 81 significantly smaller than the number of peridynamic particles (2000). On the one hand, it is obvious that the computational advantage increases with an increasing number of discrete DEM particles. On the other hand, when only a few DEM discrete particles are used and the applicable factor m between critical DEM and PD time step is not sufficiently high, the computational overhead might be higher than the gain obtained by the multi-scale integration.

General set-up
In the following, the potential of the developed multi-scale PD-DEM coupling formulation is presented on the basis of deformable-rigid mixture applications. A real-life example for this are rubber-sand mixtures. The rigid bodies are assumed to be perfect spheres and are represented by single discrete particles. In contrast, the deformable bodies are assumed to be of arbitrary shape. In the following, the simulations are performed on the cm scale and not on the length scale of real grain mixtures to save computational costs. This is sufficient for the present approach since it is the goal to reveal the potential of the multi-scale coupling framework and not to reproduce experimental based observations. For all numerical examples the inherent body forces including gravity are neglected and loads are induced by kinematic boundary conditions. Moreover, the state-based PD formulation (Eq. 14) is applied to account for incompressible behaviour of deformable bodies which cannot be covered by the bond-based approach with fixed Poisson's ratio of ν = 0.25. In the following, the deformable bodies are modelled as weakly compressible with ν = 0.48 to avoid additional kinematic constraints in the material model. The DEM and PD parameters used for the numerical examples are listed in Table 3.

Compression of deformable cube via generalised force coupling
In the first example, a three body interaction bounded by two walls, as depicted in Fig. 9, is considered. The centre body is a cubic rubber body of length 30 cm, discretised by 8000 peridynamic particles of volumes 3.375 cm 3 . The other two bodies are rigid and represented by discrete particles of radii R DEM = 15 cm. In this example friction is not considered, i.e. μ = 0. Both walls have a prescribed velocity of v wall = ±1 cm s in the z-direction to induce an indirect compression of the peridynamic cube. With a fixed time step of t = 10 −4 s the simulation is run for t = 12 s.
In this approach all forces of the developed generalised force coupling scheme, cf. Eq. (40), are coexisting. Interparticle DEM forces are acting between the walls and the adjacent discrete particles whereas intra-particle forces are acting between the peridynamic particles of the cube. Finally, PD-DEM coupling forces exist between the discrete particles and the surface of the cube (i.e. hybrid particles).
The final results at t = 12 s are depicted in Fig. 10. Note that particles on the edges are not visualised due to their missing compact support, cf. Sect.   subfigures, the elastic peridynamic cube is deformed due to the indirectly induced penetration of discrete particles. The associated strains in loading direction E zz are depicted in Fig. 10a. The extrema E zz = 0.134 are located in the centres of the free surfaces in the x-y plane, thus, at the centres of contact between the discrete particles and the peridynamic cube. In consequence of its penetration from top and bottom together with a nearly incompressible material, the elastic cube is squeezed perpendicular to the applied loading direction. This is reflected by the resulting strains E x x (Fig. 10b). Based on the introduced non-local deformation gradient in Eq. (27) and the associated Cauchy stresses Eq. (30), the von Mises stresses are computed and plotted in Fig. 10c and Fig. 10d. As expected, the highest von Mises stresses of σ V M = 1.292 N /cm 2 are obtained at the centre of the contact between the discrete particles and the peridynamic cube since it is the point of highest deformation. The capability of the developed framework to overcome the shortcoming of purely DEM based frameworks to predict stresses in soft particles is proven.

As shown in all
Overall, the induced deformations of the peridynamic cube are well captured using the developed generalised force coupling scheme. Meaningful elastic deformations and stresses of the considered deformable block, related to its weakly compressible material, are observed based on the penetration. Thus, an application to real life problems where the penetration of rubber objects by solid objects plays the superordinate role, is feasible.

Morphology changes in deformable-rigid mixtures
In the second example, the capability of the developed framework to account for major morphology changes in deformable-rigid mixtures is investigated on the basis of a simplified problem. Considered is a rectangular deformable body of dimensions 60 × 15 × 7.5 cm 3 sitting on two static discrete particles of radius R DEM = 5 cm and subjected to a induced loading by a discrete particle with same radius on top with a prescribed velocity of v = 1 cm s in the z-direction, as depicted in Fig. 11. The deformable body is equidistantly discretised by 2000 peridynamic particles with corresponding volumes of V PD = 3.375 cm 3 . In contrast to the previous example, friction is considered with a contact friction coefficient of μ = 0.546. The aim of this example is to investigate if it is possible to capture the squeezing of a lengthy deformable body through a network of rigid bodies thanks to their deformability induced via the elastic PD formulation. Performed is a completely dynamic simulation using t = 10 −4 s for t = 22.15 s. The resulting displacements in the x and z-direction of the deformable block are shown at three different time instances in Fig. 12. Over time, i.e. from t = 7.5 s (Fig. 12a and   Fig. 12b), t = 17.5 s (Fig. 12c and Fig. 12d) to t = 22.15 s (Fig. 12e and Fig. 12f) the initially straight block is gradually bend until it slips between the two static rigid discrete particles. Consequently, the block is compressed in the x-direction over time with maximum absolute values of 4.6 cm, 17.6 cm and 24.6 cm for the associated displacements. In contrast, the ends of the block are deflected in the z-direction, taking the maximum value of 7.6 cm at t = 17.5 s (Fig. 12d) before the bend block slips through the static kinematic boundary particles and the displacements in the z-direction continuously decrease.
Based on the captured dynamics, the capability of the developed coupling framework to account for major morphology changes in deformable-rigid mixtures is shown.

Compression of sphere pack
The last and most complex example is the constraint compression of a sphere pack consisting of a mix of 8 deformable and 8 rigid bodies, as illustrated in Fig. 13. All bodies have a radius of R = 15 cm, whereby each deformable body is uniformly discretised by 4166 peridynamic particles with V PD = 3.375 cm 3 . Note that by applying the uniform discretisation, the bodies of volume V = 14137 cm 3 are each represented by an accumulated PD volume of V dis = 4166 · 3.375 cm 3 = 14060 cm 3 . Consequently, a volumetric discretisation error of 0.55 % is induced.
In total, the problem consists of 8 discrete particles, 33328 peridynamic particles and 6 walls surrounding the deformable-rigid mixture. Similar to the compression of a single cube in Sect. 5.2, the loading is induced by the application of kinematic boundaries in terms of v wall = ±10 cm s for the walls in the x-z plane. The surrounding walls are modelled as static to induce constraints perpendicular to the compression direction. Similar to the previous example in Sect. 5.3, friction is considered with a contact friction coefficient of μ = 0.546. A simulation time of t = 3 s with t = 10 −4 s is considered. For the evaluation, the displacements of peridynamic particles in the z-direction as well as their norm of displacements in the x-y plane are considered at different times. The results at t = 1, t = 2 and t = 3 s are illustrated in Fig. 14. At the beginning, at t = 1 s (Fig. 14a, b) the peridynamic particles start to be compressed and a slight misalignment in the x-y plane is observable. However, the spherical bodies are still stacked on top of each other and the maximum norm of displacement in the x-y plane is 6.07 cm. After t = 2 s (Fig. 14c, d), a morphology change is observable and the bodies are not aligned on top of each other any more. On the one hand, PD bodies are compressed further and start clinging to the discrete element at the bottom. On the other hand, a discrete element is now located in the centre of the compressed packing. The morphology change is reflected in   Fig. 14c, d, the sphere pack appears to be a single clump of material, whereby the maximum norm of displacement in the x-y plane has slightly decreased to 33.79 cm. The reason for the slightly decreased maximum norm of displacement is the ongoing morphology change of individual bodies in the packing, i.e. peridynamic bodies are not only further compressed over time, but are still moving. This can lead to increasing or decreasing norms of displacement in the x-y plane of the associated peridynamic particles.
In order to quantify the results further, the mean of the norm of displacements in the x-y plane from all peridynamic particles are computed. The value gradually increases from 2.35 at t = 1 and 8.32 at t = 2 to 11.21 cm at t = 3 s. This increase of values reflects the continuous squeezing of deformable peridynamic bodies in the x-y plane. The final result of a more or less single clump of material where the deformable bodies cling to the discrete rigid particles is exactly as expected for deformable-rigid mixtures.
The results obtained for the sphere pack mixture at t = 3 s are now compared against sphere packs where all considered bodies are either rigid (DEM) or deformable (PD). The associated particle configurations including displacements in z-direction and the norm of displacements in the x-y plane are depicted in Fig. 15. In the compression of purely DEM bodies, none of the spheres is deflected in the x-y plane as shown in Fig. 15a, b. In contrast, the behaviour with PD bodies only, cf. Fig. 15c, d, is similar to the behaviour of the mixture (Fig. 14). Overall, morphology changes in the sphere pack are observable and the bodies start to form a single clump of material. The maximum norm of displacements in the x-y plane is with 32.21 cm slightly smaller than the one for the mixture (33.79 cm). Similar, the mean of the norm of displacements in the x-y plane is with 8.64 cm smaller in comparison to 11.21 cm. The reason why the mean deformations in the x-y plane for the purely deformable sphere pack are inferior to the mean deformations for the sphere mixture is related to the impact of the discrete particles in the mixture. Since the DEM bodies are rigid the entire deformation in the sphere pack mixture is carried by the deformable PD bodies and thus, by a total of eight bodies. In contrast, in the purely deformable sphere pack the deformation is distributed over all 16 bodies. Thus, the average deformation in a single deformable bodies is smaller.

Conclusion
In this contribution an efficient numerical coupling of discontinuous and continuous material behaviour based on an interface contact scheme is developed. In a first step, a generalised force coupling scheme is introduced, consisting of the superposition of inter-particle DEM, intra-particle peridynamic and inter-particle coupling forces. The key part in the formulation are surface particles of peridynamic bodies acting as hybrid particles. Based on the formulation, a multi-scale time integration scheme is proposed leading to the generalised multi-scale PD-DEM coupling framework.
It is shown that the trend of the analytical 3D Hertzian contact theory is captured reasonably well within the developed PD-DEM coupling framework, without calibrating the micro-mechanical DEM contact parameters. Thus, a proper force coupling of the discontinuous-continuous numerical approach is verified. The multi-scale time integration is verified on the basis of a fully dynamic example. It is shown that the staggered integration leads to a slightly smaller dynamic impact than the monolithic scheme due to subsequent kinematic updates of hybrid particles within a global time step. However, the general behaviour is similar. Comparing the computational efficiency it is possible to reduce the computation time up to 21 % when using the staggered scheme for the simple example considered. It is expected that the computational efficiency significantly increases for more complex simulations with higher number of peridynamic particles. The reason behind is the general non-linear scaling of the computation time with respect to the number of particles considered. Most of the computation time is spend for collision detection, which is considerably decreased when applying the multi-scale time integration in case of t DEM < t PD .
The coupling framework is then applied to the simulation of systems with deformable and rigid bodies. Firstly, the successful simulation of the penetration of rigid discrete elements into an elastic body is shown. Secondly, the ability of the framework to represent the major morphology change of squeezing a lengthy deformable body through rigid spheres is presented. Since it is essential to capture this phenomena in the simulation of deformable-rigid grain  mixtures under compressible loads, it verifies the applicability of the developed framework for this kind of problems. In the last numerical example the actual compression of a deformable-rigid sphere mixture is considered. Using this example, the capability of the framework to capture the effect of deformable bodies being squeezed and clinging to rigid bodies is proven. As a consequence, a successful application of the developed numerical continuous-discontinuous coupling framework for the simulation of grain mixtures is presented.
In future it is desirable to apply the approach for reallife grain mixtures, e.g. for rubber-sand mixtures. In order to predict their dynamic behaviour appropriately, only two steps are necessary. Firstly, it is required the calibrate the macroscopic material parameters of the rubber by laboratory experiments. In a subsequent step a classical calibration of required DEM contact parameters is necessary, whilst already applying the rubber parameters. Besides specific grain mixtures, the framework could also be applied to cap-ture the penetration of tires by stones on gravel road with stones being possibly trapped in grooves of the tire thread.
Overall, there are numerous possibilities to use the developed approach or to make it applicable to other problems with simple extensions or changes. One of these approaches is the simulation of grain crushing of specific stones in sand grain mixtures. Whilst the sand can still be modelled by DEM bodies it would only be necessary to model the crushable stones by PD bodies and to apply a material law including fracture.
The implementation of the coupling scheme into Yade offers new extensive additional possibilities. Generally, DEM bodies do not have to be spherical particles and more complex DEM shapes can be applied within the framework. Moreover, the applications are not limited to the presented Hertz-Mindlin contact model since various contact models have already been implemented. Generally, it is also possible to apply different contact models within a mixture.
Further, the PD-DEM force coupling scheme is implemented in a generalised object-orientated manner with respect to the applied multi-scale time integration scheme. Thus, the performed contact coupling in Yade is not limited to PD-DEM couplings and can be straight forward extended to other DEM coupling schemes.