A modified bond model for describing isotropic linear elastic material behaviour with the particle method
Rahav Gowtham Venkateswaran, Ursula Kowalsky, Dieter Dinkler
AA modified bond model for describing isotropic linear elasticmaterial behaviour with the particle method
R. G. Venkateswaran ∗ , U. Kowalsky, and D. DinklerInstitute of Structural AnalysisTechnische Universit¨at BraunschweigBraunschweig, Germany Abstract
Particle based methods such as the Discrete Element Method and the LatticeSpring Method may be used for describing the behaviour of isotropic linear elasticmaterials. However, the common bond models employed to describe the interactionbetween particles restrict the range of Poisson’s ratio that can be represented. Inthis paper, to overcome the restriction, a modified bond model that includes thecoupling of shear strain energy of neighbouring bonds is proposed. The coupling isdescribed by a multi-bond term that enables the model to distinguish between sheardeformations and rigid-body rotations. The positive definiteness of the strain energyfunction of the modified bond model is verified. To validate the model, uniaxial ten-sion, pure shear, pure bending and cantilever bending tests are performed. Compar-ison of the particle displacements with continuum mechanics solution demonstratesthe ability of the model to describe the behaviour of isotropic linear elastic materialfor values of Poisson’s ratio in the range 0 ≤ ν < . Keywords:
DEM, strain energy, multi-bond, shear, rigid-body rotation, Poisson’sratio
Regarding the atomic scale a material is made up of atoms which are known to interactwith each other through attractive and repulsive forces. Extrapolating this model tothe macroscopic scale is the fundamental idea of modelling continua with interactingdiscrete particles. Referring to the historical development of the theory of elasticity[14, 4] one realises that the mechanicians of the nineteenth century have made substantialcontribution to the development of this alternate discrete model of elasticity. Apart fromelasticity, it allows one to describe the failure of material as the absence of interactionbetween two previously interacting material points. Two of the most widely used particle ∗ Correspondence: [email protected] a r X i v : . [ c s . C E ] N ov ethods are the Lattice Spring Model (LSM) and the discrete element method (DEM)initially developed for modelling the motion of granular media. The LSM influencedby crystal elasticity discretises the domain with material points in an ordered manner,whereby each point interacts with each other by means of an elastic spring. The DEMinitially developed by Cundall [6] has been extended to model solids such as rocks bythe bonded-particle method (BPM) where two particles remain bonded together as longas a critical bond deformation is reached [17]. Although these methods are different intheir implementation they are similar in the fact that they use bond models to describethe interaction between two particles. Several bond models exist in the literature anda brief summary of the most commonly used models is provided here. For a detailedreview on the various bond models refer to [16].The fundamental bond model is a Hookean spring which penalises change from thereference length between two interacting particles. For a solid discretised with an infinitenumber of material points connected by such springs, the Poisson’s ratio is restrictedto 1 / / Cauchy symmetry are fulfilledand lead to the restriction on Poisson’s ratio [14, 4]. In order to overcome the restric-tion, Born [2] introduced a shear spring in addition to the longitudinal spring. The shearsprings can be understood as a penalisation of the change in orientation from the ref-erence configuration. This model has two stiffness parameters available for calibration.However, for an infinite number of such bonds, the Poisson’s ratio is restricted to therange 0 ≤ ν ≤ and 0 ≤ ν ≤ for the case of plane stress and plane strain respectively.Apart from these restrictions of Poisson’s ratio, the model cannot distinguish betweenrigid-body rotation and shear deformation [11, 3, 12]. Instead of using shear springs, theKirkwood-Keating spring model [13, 12] uses an angular spring to penalise the angularmotion. However, this model is known to be nonlinear due to the angular terms and itoffers no substantial advantage in comparison to the Born model [3]. The Lattice Beammodel (LBM) uses a Hookean spring to allow longitudinal forces and a beam to allowshear force and bending moment and is widely used in DEM. However, in DEM the purebending modes are neglected [18, 10]. The beam model has the similar restriction on thePoisson’s ratio as the
Born model but is capable of distinguishing between rigid-bodyrotation and shear deformation due to inclusion of the rotational degree of freedom.In order to overcome the limited range of Poisson’s ratio that can be represented bythese models, Zhao [21] evaluated the shear strains in the
Born model using a particlestrain tensor instead of particle displacements. This tensor is an approximate measureof the continuum strain tensor and is obtained by the least square method using theinformation of displacements of the particle under consideration and its neighbouringparticles. The calculation of the particle strain tensor is computationally expensive andfor a practical simulation with random particle arrangements there exists ambiguityregarding its existence. In such scenario only Hookean springs are used. Celigueta [5]proposed a nonlocal contact law for the DEM, where apart from the overlap of two2nteracting particles in contact also the influence of forces in the surrounding of thiscontact is included. By means of this term a better description of continuum elasticitywas obtained in comparison to the classical bond model of DEM. However, in orderto calculate this term information regarding the contact area between two particles incontact and their surroundings are required. Additionally, the nonlocal stress tensor ofeach particle and its neighbours are required which results in a computationally expensivemethod.In this work, a modification of the
Born model is proposed by introducing a multi-bond strain energy term that is capable of distinguishing between shear deformationand rigid-body rotation for the case of two-dimensional elasticity. Two new multi-bonds called
L-bond and
X-bond employing the proposed coupling of shear strain energyare introduced. Based on positive definiteness of the strain energy function, we alsointerpret that the restriction of the
Born model up to a certain Poisson’s ratio resultsfrom incapability to distinguish between shear deformation and rigid-body rotation.With the modified model it is shown that the strain energy function remains positivedefinite for values of Poisson’s ratio in the range 0 ≤ ν < .
5. To validate the modifiedmodel, the results of a thin plate under uniaxial tension, pure shear, pure bending andcantilever bending are compared with the respective plane stress continuum mechanicssolution.
Born model
An arbitrary rectangular domain discretised with rigid circular particles of radius r isshown in Figure 1(a). Each particle is bonded with its immediate neighbours (separatedby 2 r ) and also with its second neighbours (separated by 2 √ r ). Similar to the works [10,15], the fundamental building block called as a unit-cell is extracted from the discretiseddomain and is shown in Figure 1(b). This square unit-cell is the combination of two unit-cells: one which contains only the bonds with the immediate neighbours and anotherthat contains only the bonds between second neighbours. They are called first-neighbourand second-neighbour unit-cells as shown in Figure 1(c) and 1(d) respectively. Althoughother configurations such as a triangle or a hexagon can be used as a unit-cell, thesquare configuration is used here due to its simplicity. The modification proposed inthis paper is not restricted to a square configuration and it can be used also regardingother regular configurations. The Born model which is applied here employs springs innormal direction and tangential direction with stiffness parameters k n and k s as shownin Figure 2(a). A generic bond oriented at angle θ to the global coordinate system alongwith its local and global displacement components is shown in Figure 2(b).Regarding Figure 2(a), the strain energy stored in the bond Π b is given in terms ofquantities with respect to the bond coordinate system asΠ b = 12 k n ( u nB − u nA ) + 12 k s ( v sB − v sA ) (1)= l (cid:0) k n ( (cid:15) nn ) + k s ( (cid:15) sn ) (cid:1) , (2)3 a) A BCD (b)
A BCD (c)
A BCD (d)
Figure 1: (a) Domain discretised with the chosen unit-cell (b) Square unit-cell obtainedfrom the combination of (c) First neighbour unit-cell and (d) Second neighbour unit-cell
A Bk n k s (a) A Bx, uy, v n, u n s, v s l θ (b) Figure 2: (a) Springs used in a bond (b) Coordinate system of a generic bond4here (cid:15) nn = u nB − u nA l and (cid:15) sn = v sB − v sA l are the local longitudinal and shear strain respec-tively and l is the length of the undeformed bond.After substituting the local strain quantities in terms of the components of the globalstrain tensor by the transformation relations provided by Griffiths [10] the strain energyin a bond is obtained asΠ b = l (cid:18) k n ( c (cid:15) xx + sc ( (cid:15) xy + (cid:15) yx ) + s (cid:15) yy ) + k s ( c (cid:15) yx + sc ( (cid:15) yy − (cid:15) xx ) − s (cid:15) xy ) (cid:19) (3)where, c = cos θ and s = sin θ . By substituting the orientation θ of the individualconstituent bonds in Equation (3), the strain energy in the first and second neighbourunit-cells are obtained asΠ uc = Π normaluc + Π shearuc = k bn l (cid:20) ( (cid:15) xx ) + ( (cid:15) yy ) + ( (cid:15) xx ) + ( (cid:15) yy ) (cid:21) + k bs l (cid:20) ( (cid:15) yx ) + ( − (cid:15) xy ) + ( (cid:15) yx ) + ( − (cid:15) xy ) (cid:21) , (4)Π uc = Π normaluc + Π shearuc = k bn l (cid:20)(cid:18) (cid:15) xx (cid:15) xy (cid:15) yx (cid:15) yy (cid:19) + (cid:18) (cid:15) xx − (cid:15) xy − (cid:15) yx (cid:15) yy (cid:19) (cid:21) + k bs ( √ l ) (cid:20)(cid:18) (cid:15) yx − (cid:15) xy − (cid:15) xx (cid:15) yy (cid:19) + (cid:18) (cid:15) xx − (cid:15) xy (cid:15) yx − (cid:15) yy (cid:19) (cid:21) , (5)where k bn , k bs are the normal and shear stiffness of first neighbour bonds and similarly k bn , k bs are the stiffness parameters of second neighbour bonds. All bonds in the sameshear plane were assigned the same stiffness k bs . Every bond in the first neighbour unit-cell contributes only 1 / uc = Π uc + Π uc e uc = Π uc l t , (6)where t is the thickness of the domain.As described before, the Born model cannot distinguish between rigid-body rotationand shear [12, 11, 21, 3]. In Figure 3(a) an exploded view of each constituent bond of theunit-cell with first neighbour bonds is shown along with its local and global shear strains.The bonds are assembled back together in Figure 3(b) and one observes that although thefinal geometry describes pure rotation, the strain energy is not zero. Thus strain energy5
Bns(cid:15)
ABsn = (cid:15) yx BCns (cid:15) B C s n = − (cid:15) x y CD n s(cid:15)
CDsn = (cid:15) yx DA n s (cid:15) D A s n = − (cid:15) x y (a) A BCD (b)
A Cns BD ns (c)
A BCD (d)
Figure 3: (a) Split up of first neighbour unit-cell bonds (b) Assembly of first neighbourbonds (c) Split up of second neighbour unit-cell bonds (d) Assembly of second neighbourbondsis stored in the unit-cell due to shear of individual bonds although all angles subtendedbetween neighbour bonds are 90 ◦ . This can also be observed in Equation (4) where thecontribution of shear strain energy Π shearuc remains non-zero for any combination of (cid:15) xy and (cid:15) yx except for the trivial case where (cid:15) xy = (cid:15) yx = 0. This interpretation holds alsofor the unit-cell with the second neighbour bonds as shown in Figure 3(c), (d). To overcome the limitations explained in Section 2, a modified model is proposed wherethe shear strains of two neighbour bonds (multi-bond) are coupled. With this coupling,strain energy is stored only in the case of shear and not in the case of rotation. Hereagain, the unit-cell made up of only the first neighbour bonds is looked at first and theunit-cell made up of the second neighbour bonds afterwards.6 .1 L-bond
For a generic
L-bond
ABC made up of two first neighbour bonds AB and BC as shownin Figure 4(a), the strain energy due to normal strain remains unchanged as before.However the shear strains are combined in such a way that strain energy is stored onlyfor shear deformation and not for rigid-body rotation. The strain energy stored in ageneric
L-bond with length l AB = l BC = l is given byΠ ABC = 12 k mn (cid:20)(cid:0) (cid:15) ABnn l (cid:1) + (cid:0) (cid:15) BCnn l (cid:1) (cid:21) + 12 k ms (cid:20) − (cid:15) ABsn l + (cid:15) BCsn l (cid:21) , (7)where k mn and k ms are the normal stiffness and shear stiffness parameters employingthe modified model respectively. The first neighbour unit-cell is now made up of fourL-bonds as shown in Figure 4(b). AB Cnsn sxy (a)
DA BD CB A BCCDA (b)
Figure 4: (a) Generic
L-bond (b) First neighbour unit-cell made up with
L-bonds
As explained in Section 2, a first neighbour bond has a contribution factor of 1 / AB is obtained from DAB and
ABC . Therefore, in an L-bond the normalstiffness k mn and the shear stiffness k ms have a contribution factor of 1 /
4. With the trans-formation relations provided by Griffiths [10], the local strain components are writtenin terms of the global strain components and the strain energy of a generic
L-bond now7ields Π
ABC = k mn l (cid:20)(cid:18) (cid:15) xx c AB + (cid:0) (cid:15) xy + (cid:15) yx (cid:1) c AB s AB + (cid:15) yy s AB (cid:19) + (cid:18) (cid:15) xx c BC + (cid:0) (cid:15) xy + (cid:15) yx (cid:1) c BC s BC + (cid:15) yy s BC (cid:19) (cid:21) + k ms l (cid:20) − (cid:18) (cid:15) yx c AB + (cid:0) (cid:15) yy − (cid:15) xx (cid:1) c AB s AB − (cid:15) xy s AB (cid:19) + (cid:18) (cid:15) yx c BC + (cid:0) (cid:15) yy − (cid:15) xx (cid:1) c BC s BC − (cid:15) xy s BC (cid:19)(cid:21) , (8)where, c AB = cos θ AB , s AB = sin θ AB , c BC = cos θ BC and s BC = sin θ BC . Employingthe modified bond model, the strain energy stored in the first neighbour unit-cell isobtained from the sum of strain energy stored in individual L-bonds Π moduc = Π ABC + Π
BCD + Π
CDA + Π
DAB . (9)The strain energy stored in the first neighbour unit-cell employing the modified bondmodel can be evaluated with respect to the global coordinate system toΠ moduc = k mn l (cid:0) (cid:15) xx + (cid:15) yy (cid:1) + k ms l (cid:0) (cid:15) xy + (cid:15) yx (cid:1) = k mn l (cid:0) (cid:15) xx + (cid:15) yy (cid:1) + k ms l (cid:0) (cid:15) xy + (cid:15) yx (cid:1) + k ms l (cid:15) xy (cid:15) yx = Π normaluc + Π shearuc + k ms l (cid:15) xy (cid:15) yx . (10)Comparing the formulation to that of the first neighbour unit-cell employing the Bornmodel , see Equation (4), an extra term is added due to the coupling of the shear strainenergies of the neighbours. With this term the modified model is able to distinguishbetween shear (change in the angle subtended) and rigid-body rotation (no change inthe angle subtended), which can be interpreted as a modification to the
Born model . X-bond
Similar to
L-bonds , the unit-cell made up of the second neighbour bonds is to be modifiedsuch that strain energy is stored only in the case of shear deformations and not for purerotation. This is achieved with the
X-bonds . Again the shear strains of two neighbouringbonds are coupled and the normal strains of the bonds remain unchanged. The modifiedstrain energy of the unit-cell made up of one
X-bond isΠ moduc = Π AC + Π BD = k mn ( √ l ) (cid:20)(cid:0) (cid:15) ACnn (cid:1) + (cid:0) (cid:15) BCnn (cid:1) (cid:21) + k ms ( √ l ) (cid:20) (cid:15) ACsn − (cid:15) BDsn (cid:21) , (11)8here, k mn and k ms are the normal and shear stiffness parameters employing the modifiedmodel respectively. The shear stiffness parameters in the same shear plane are assumedto have the same stiffness parameter k ms . The strain energy stored in the second neigh-bour unit-cell with respect to the global coordinate system isΠ moduc = k mn ( √ l ) (cid:20)(cid:18) (cid:15) xx (cid:15) xy (cid:15) yx (cid:15) yy (cid:19) + (cid:18) (cid:15) xx − (cid:15) xy − (cid:15) yx (cid:15) yy (cid:19) (cid:21) + k ms ( √ l ) (cid:20)(cid:18) (cid:15) yx − (cid:15) xy − (cid:15) xx (cid:15) yy (cid:19) − (cid:18) (cid:15) xx − (cid:15) xy (cid:15) yx − (cid:15) yy (cid:19)(cid:21) . (12)Upon expansion of Equation (12) and comparison with the strain energy of the secondneighbour unit-cell, see Equation (5), one observes that similar to the first neighbourunit-cell with L-bonds , the coupling of shear strain energies results in an extra multi-bondterm to distinguish between shear and rotation.
The strain energy stored in a square unit-cell is obtained by summing up the strainenergy of an unit-cell made with
L-bonds , Equation (10), and of an unit-cell with
X-bonds , Equation (12) Π moduc = Π moduc + Π moduc . (13)The strain energy density e moduc is obtained similar to that of the unit-cell with the Bornmodel given in Equation (6).
Born model
In this section, the stiffness parameters of a square unit-cell employing the modifiedmodel and the
Born model are calibrated to the macroscopic elastic material parameters.The calibrated parameters are verified if they satisfy the condition of isotropy. The strainenergy function of the unit-cell employing these bond models is verified for its positivedefiniteness.
In order to model the behaviour of an isotropic elastic material, the stiffness parame-ters of the unit-cell employing the modified model must be calibrated with respect tothe macroscopic elastic material parameters. In the literature there exist broadly twoapproaches:1. A numerical calibration where the stiffness parameters are iteratively calibratedby solving an inverse problem to match the slope of the stress-strain curve of amaterial in the linear region under uniaxial loading. For more details refer [20, 8].9. An analytical calibration based on the equivalence of strain energy density of theunit-cell e uc with that of an equivalent elastic continuum e co = σ : (cid:15) . In partic-ular, the components of the elastic tensor C uc are derived from the strain energydensity of the unit-cell and are compared to the continuum description.In this work, the analytical calibration approach is used due to its independence ofthe fineness of discretisation. The components of the tensor of elasticity of the unit-cell are obtained by differentiating the strain energy density twice with respect to thecorresponding strain components. The individual components are summarised in thetensor of elasticity C moduc = ˆ C ˆ C C ˆ C
00 0 ˆ C , (14)where, ˆ C = k mn + 2 k ms + k mn , ˆ C = k mn − k ms and ˆ C = k mn + k ms are defined for themodified model and ˆ C = k bn + k bs + k bn , ˆ C = k bn − k bs and ˆ C = k bn + k bs forthe Born model . The elasticity tensor of the unit-cell has three components similar tothat of a planar continuum. By comparing the three components ˆ C , ˆ C and ˆ C withthose of the planar continuum elasticity tensor, the stiffness parameters are calibrated tothe macroscopic Young’s modulus E and the Poisson’s ratio ν . The calibrated stiffnessparameters for the case of plane stress and plane strain are summarised in Table 1.Table 1: Calibrated stiffness parameters of square unit-cellBond model Type Stiffness parameters k bn k bs k bn Born
Plane stress
E t (1+3 ν )3 (1+ ν ) (1 − ν ) E t (1 − ν )3 (1+ ν ) (1 − ν ) E t ν ) (1 − ν ) Plane strain
E t (1+2 ν )3 (1+ ν ) (1 − ν ) E t (1 − ν )3 (1+ ν ) (1 − ν ) E t (1 − ν )3 (1+ ν ) (1 − ν ) k mn k ms k mn Modified Plane stress
E t (1+3 ν )3 (1+ ν ) (1 − ν ) E t (1 − ν )6 (1+ ν ) (1 − ν ) E t ν ) (1 − ν ) Plane strain
E t (1+2 ν )3 (1+ ν ) (1 − ν ) E t (1 − ν )6 (1+ ν ) (1 − ν ) E t (1 − ν )3 (1+ ν ) (1 − ν ) As expected for a
Born model there exists a restriction on the Poisson’s ratio at ν = 1 / ν = 1 / Born model and the modified model10 . . . . . . ν [-] S t i ff n e ss p a r a m e t e r s [ - ] k bn k bn k bs k mn k mn k ms (a) Born and modified model . . . . . . . ν [-] ˆ C i j [ - ] ˆ C CM ˆ C CM ˆ C CM ˆ C uc ˆ C uc ˆ C uc (b) Born and modified model
Figure 5: Variation of the stiffness parameters (a) and two-dimensional elasticity tensorcomponents (b) as a function of ν for plane stressare identical since only the shear strain energy components were coupled in the modifiedbond model. The stiffness parameters and the elasticity constants of the unit-cell areplotted as a function of the Poisson’s ratio ν for both bond models in Figure 5 for thecase of plane stress. The stiffness parameters are normalised with respect to the Young’smodulus E and the thickness of the specimen t to obtain a qualitative behaviour.From the Figure 5(b) we observe that the components of the elasticity tensor of theunit-cell agree exactly with the equivalent continuum components in the range 0 ≤ ν ≤ .
5. In order to ensure the equality of the elasticity constants, the stiffness parametersof the unit-cell must take on the values as shown in Figure 5(a). Therefore, we interpretthe negative value of shear stiffnesses ( k bs , k ms ) for ν > / Apart from calibrating the stiffness parameters to the macroscopic isotropic elastic ma-terial parameters it must also be checked whether and in which case they satisfy thecondition for macroscopic isotropy. A material is considered to exhibit isotropy if theanisotropy factor Λ = 1. The anisotropic factor for the chosen unit-cell is calculated11ith the components of the elasticity tensor for both bond models as
Born model : Λ = 2 ˆ C ˆ C − ˆ C = 2 k bn + k bs k bn + 2 k bs (cid:40) = 1 , k bn = k bn + k bs (cid:54) = 1 , otherwise (15)Modified model: Λ = 2 ˆ C ˆ C − ˆ C = 2 k mn + 2 k ms k mn + 4 k ms (cid:40) = 1 , k mn = k mn + 2 k ms (cid:54) = 1 , otherwise . (16)This shows that for isotropic elasticity (Λ = 1) the unit-cell has only two independentstiffness parameters and a condition that ensures isotropy for all values of E and ν . Bysubstituting the stiffness parameters from Table 1 in the condition for isotropy it can beshown that this condition holds for both plane stress and plane strain. The square unit-cell shown in Figure 1(b) is made up of four particles with two de-grees of freedom per particle. The vector of displacements for one unit-cell is u =[ u A , v A , u B , v B , u C , v C , u D , v D ] T . With the displacement vector u and the stiffness ma-trix K moduc given in Appendix B, the strain energy of the unit-cell in Equation (13) canalternatively be written as Π moduc = 12 u T K moduc u . (17)Equation (17) shows that the strain energy of the unit-cell is a quadratic functionof displacements. In order to be thermodynamically stable, it must be positive definite.The necessary and sufficient condition for the positive definiteness of Equation (17) isthat all eigenvalues of the eigenvalue problem K moduc u = λ u are real and positive [9]. Therelating stiffness matrices are given in Appendix A and B respectively. The resultingeigenvalues and corresponding eigenforms are summarised in Table 2. After substitutingthe stiffness parameters of the unit cell from Table 1, the non-zero eigenvalues are plottedas a function of Poisson’s ratio ν for both plane stress and plane strain in Figure 6. Thestiffness parameters are normalised with respect to Young’s modulus E and the thickness t to obtain a qualitative behaviour. 12able 2: Eigenvalues and eigenforms of unit-cell with the Born model and modifiedmodel λ i Eigenform
Born model Modified bond model λ λ λ k bs λ k bn + k bs k mn + k ms λ k bn + k bs k mn + k ms λ k bn + k bs k mn + 2 k ms λ k bn + 2 k bs k mn + 4 k ms λ k bn + 2 k bn k mn + 2 k mn From Figure 6 the following observations are made for both plane stress and planestrain:1. Regarding the
Born model , the eigenvalue λ of the rigid-body rotation eigenformis non-zero everywhere expect at ν = 1 / ν = 1 / Born model deforms under rigid-bodyrotation and therefore it is proven that this model cannot distinguish between rigid-body rotation and shear deformations. Also, the eigenvalue becomes negative for ν > / ν > / Born model is not positive definiteafter this lower limit of the Poisson’s ratio. Therefore we expect the response tobe unstable for values of Poisson’s ratio above this limit. However, regarding themodified bond model, the eigenvalue of the rigid-body rotation eigenform remainszero for all values of ν and all eigenvalues remain positive. Hence, the strainenergy derived from the modified bond model is positive definite and therefore isstable for all values of ν although the shear stiffness k ms of the unit-cell takes ona negative value after a critical value of ν as shown in Figure 5(b). This impliesthat a negative stiffness does not necessarily lead to unstable results, rather it is13ue to negative eigenvalues. Similarly, in the work of Esin [7] it was concludedthat materials and structures with negative stiffness elements can exist when thenegative element energy is compensated by the energy of the rest of the system oran encompassing system that provides stabilisation.2. The eigenvalues ( λ , λ ) related to bending eigenforms result as expected for anisotropic material.3. Similarly, the eigenvalues of shear ( λ , λ ) eigenforms result as expected.4. As expected, the volumetric form λ has the highest eigenvalue. . . . . . − ν [-] λ i [ - ] λ λ λ λ λ λ (a) Born model - plane stress . . . . . − ν [-] λ i [ - ] λ λ λ λ λ λ (b) Modified model - plane stress . . . . . − ν [-] λ i [ - ] λ λ λ λ λ λ (c) Born model - plane strain . . . . . − ν [-] λ i [ - ] λ λ λ λ λ λ (d) Modified model - plane strain Figure 6: Variation of normalised eigenvalues of unconstrained unit-cell stiffness matrixas a function of ν emark: At ν = 1 / ν = 1 / Born model is rotationally invariant and hence we have λ = 0 at thisparticular value of ν . In order to validate the capability of the modified model and compare it to the
Bornmodel , four different examples are chosen for which closed-form continuum solutionsexist. With the help of these examples and based on the positive definiteness of the strainenergy function, the stability of the bond model will also be investigated. The discreteelement method results are obtained by solving the linear system of equations
K u = f for the unknown particle displacements u . The system stiffness matrix K is obtainedby assembling the individual unit-cell stiffness matrices derived in the Appendix A andB. External tractions are reformulated as forces and are added to the correspondingposition in the global force vector f . A thin square plate of length l = 0 . t = 0 .
01 [m] as shown in Figure7 is modelled with the chosen square unit-cell under plane stress conditions with E =2 × [N / m ]. The plate is constrained in y − direction along the bottom and in the x − direction along the left. A stress σ xx is applied on the right edge to model a uniaxialstate of deformation. σ xx lxy Figure 7: Thin plate under uniaxial tensionThe analytical solution for the displacement fields is u = σ xx E x and v = − ν σ xx E y . Theresulting displacements along u ( x = l, y ) and v ( x, y = l ) for three different values ofPoisson’s ratio ν are summarised in the Table 3 for both the bond models.We observe that both unit-cells employing the Born model and the modified bondmodel produce the exact results for all values of the Poisson’s ratio independent of the15able 3: Comparison of uniaxial tension test results (Displacements in [ × − m]) ν Analytical solution
Born model
Modified model u v u v u v . . . . . − .
03 0 . − .
03 0 . − . .
49 0 . − .
049 0 . − .
049 0 . − . ν > /
3, they still produce the expected result. In Section 4.3, based onthe eigenvalues of the stiffness matrix K uc of an unconstrained unit-cell employing the Born model it was observed that the strain energy was not positive definite (due tothe rigid-body rotation eigenform) for ν > / ν .This is shown in Figure 8, where the normalised eigenvalues of the constrained unit-cellemploying the Born model and the modified model are plotted as a function of ν . . . . . . . . ν [-] λ i [ - ] λ λ λ λ (a) Born model . . . . . . . ν [-] λ i [ - ] λ λ λ λ (b) Modified model Figure 8: Variation of eigenvalues of constrained unit-cell stiffness matrix as a functionof ν for uniaxial test As a second test, a thin square plate with edge length l = 0 . t =0 .
01 [m] as shown in Figure 9 is modelled under plane stress condition with E = 2 × [N / m ]16 xy σ xy σ xy σ yx Figure 9: Thin plate under pure shearapplying pure shear stress σ xy = σ yx = τ = 1 × [N / m ].The plate is constrained in both directions along the bottom edge to prevent rigid-body motion similar to that used by Ockelmann [15]. The analytical solution for thedisplacement fields is obtained as u = τ G y and v = 0, where G = E/ ν ) is theshear modulus. Here again, the test is conducted for three different values of ν and alsofor four different discretisations similar to the uniaxial test. The u and v displacementsof the domain discretised with 8 × Born model under pureshear are summarised in Figure 10(a), (b) and (c). In Figure 10(d), the results of 8 × ν = 0 .
49 are given. InFigure 10 the bonds connecting the particles are not visualised for clarity.We observe that the results for the unit-cell employing the modified model agreesexactly with the analytical solution, because similar to the case of uniaxial tension, theanalytical displacement field is a linear function of the position of particles. With themodified bond model, the exact solution was also observed for all the discretisations.However, for the unit-cell employing the
Born model the response of the plate is stifferin comparison with the analytical solution in the range 0 < ν ≤ / ν > /
3. In order to understand the reason for this behaviour, the plate is discretisedwith one unit-cell employing the
Born model and the modified model respectively. Thenormalised eigenvalues of the constrained stiffness matrix are plotted as function of ν asshown in Figure 11. For the stiffness matrix of an unit-cell employing the Born model ,we observe from Figure 11 that the eigenvalue λ b decreases with increasing values of ν and also becomes negative for ν > .
4. The eigenform corresponding to this eigenvalue isplotted as given as a quiver plot in Figure 12. However, as the plate is further discretised,the eigenvalues of the unit-cells that are not located at the bottom edge of the plate takeon the form as previously shown in Figure 6(a) for the case of an unconstrained unit-cellstiffness matrix employing the Born model. This unconstrained stiffness matrix alsoincludes the contribution of rigid-body rotation and is unstable for values of ν > / [m]1 e − xy v [m]4 . e − − . e − (a) Born model : ν = 0 u [m]1 . e − xy v [m]3 . e − − . e − (b) Born model : ν = 0 . u [m]3 e − − . e − xy v [m]1 . e − − . e − (c) Born model : ν = 0 . u [m]2 . e − xy v [m]4 e − − e − (d) Modified model: ν = 0 . Figure 10: Displacements u (left) and v (right) employing the Born model and themodified model 18he case of plane stress. Therefore, the upper limit on the Poisson’s ratio for the plateconstrained as shown in Figure 9 is ν = 1 / Born model . However, forthe unit-cell employing the modified model, due to the coupling of shear strain energyof neighbouring bonds, the eigenvalues remains positive for all values of ν as shown inFigure 11 for the constrained stiffness-matrix and in Figure 6(b) for the unconstrainedstiffness matrix respectively. . . . . . . . ν [-] λ i [ - ] λ b λ b λ b λ b λ m λ m λ m λ m Figure 11: Variation of eigenvalues of con-strained unit-cell stiffness matrix as a func-tion of ν for pure shear test Figure 12: Eigenform corresponding to λ Next, a thin rectangular plate of length l = a = 0 . h = 2 b = 0 .
125 [m] andthickness t = 0 .
01 [m] under pure bending is considered. The plate is constrained at twopoints in the y − direction and at one location in the middle along x − direction as shown inFigure 13. The applied bending moment M = 2604 .
17 [Nm] at the sides is modelled as alinearly varying load σ with opposite magnitude at both corners. Following Timoshenko[19], the analytical displacement fields are u = MyEI ( − x + l x ) and v = M EI ( νy + x − xl x ),where I = th / aM Mxy bb Figure 13: Thin plate under pure bending19 a) Born model (b) Modified model
Figure 14: Deformed configuration for ν = 0 . ×
2, 16 ×
4, 32 × ×
16. Also, the effect of ν on the accuracy has been checked by performing the testfor three different values of ν . The deformed configuration employing the Born model and the modified model for ν = 0 .
49 is given in Figure 14.The displacement u at the left edge of the plate and the displacement v of the axisof the plate employing the Born model and the modified model are compared with theanalytical solution for different discretisations and for different Poisson’s ratio and aresummarised in Figure 15. From the results we observe that regarding the
Born model for ν = 0 and ν = 0 .
3, the displacements u and v are much smaller in comparisonwith that of the analytical solution. In the case of pure bending, elements which arelocated faraway from the boundary undergo rigid body rotations as well. Since the Bornmodel cannot distinguish between rigid-body rotations and shear deformations, strainenergy is also stored for the rigid body rotation. As the rigid-body rotation eigenvaluetakes on a maximum value at ν = 0 and approaches zero at ν = 1 / Born model , the response of the plate is much stiffer for ν = 0than for ν = 1 /
3. And for ν > / Born model , that the displacements do notconverge. This is due to the non-zero eigenvalue of the rigid-body rotation eigenform.For the unit-cell employing the modified bond model, the coupling of the shear strainenergy of the neighbouring bonds allows the model to distinguish between rigid-bodyrotations and shear deformations. Therefore, satisfactory results are obtained employingthis modified bond model and also convergence to the analytical solution upon refinementis observed. 20 nalytical solution Born - 8 × × × ×
16 Modified - 8 × × × × − − . . · − − − − · − u [m] y [ m ] (a) ν = 0 . . . . . − − · − x [m] v [ m ] (b) ν = 0 − − . . · − − − − · − u [m] y [ m ] (c) ν = 0 . . . . . . − − · − x [m] v [ m ] (d) ν = 0 . − − . . · − − − − · − u [m] y [ m ] (e) ν = 0 . . . . . . − − · − x [m] v [ m ] (f) ν = 0 . Figure 15: Displacements u at the left edge of the plate (left) and v of the axis of theplate (right) employing the Born model and the modified model21 .4 Cantilever bending
As a final test, a thin rectangular plate under cantilever bending is considered. Theright edge of the plate is completely constrained and a constantly distributed load F isapplied on the left edge as shown in Figure 16. The plate is modelled with the chosenunit-cell under plane stress conditions. F axy bb a = l = 0 . b = h/ . t = 0 .
01 [m] E = 2 × [N / m ] F = 1 . × [N / m]Figure 16: Thin plate under cantilever bendingFrom [1], the analytical solution for the displacement fields are obtained as u = 3 F x y Eb + 3 F (1 + ν ) y Eb − F (2 + ν ) y Eb − F a y Eb − F a y Eb (cid:18) ν ) b a (cid:19) (18) v = − F νxy Eb − F x Eb − F a Eb (cid:18) ν ) b a (cid:19) + 3 F a x Eb (cid:18) ν ) b a (cid:19) . (19)For the derivation of the analytical solution the following weak boundary conditions wereapplied along the right edge ( x = a ): b (cid:90) − b u dy = 0 , b (cid:90) − b v dx = 0 and b (cid:90) − b y v dy = 0 . (20)Similar to the pure bending test, four different discretisations are used. The effectof ν on the results is studied as well. The deformed configuration employing the Bornmodel and the modified model are given in Figure 17.The distributions of the displacement u at the left edge of the plate and the displace-ment v of the axis of the plate regarding the Born model and the modified model aresummarised in Figure 18. From the results for ν < / Born model is stiffer in comparisonwith that of the analytical solution. This is because, for the case of cantilever bendingthere exists elements faraway from the boundary that undergo rigid-body rotation aswell. Since the
Born model cannot distinguish between rigid-body rotations and shear22 a) Born model (b) Modified model
Figure 17: Deformed configuration for ν = 0 . ν > / Bornmodel produces unstable results. The reason for this behaviour can be understood bylooking at the variation of normalised eigenvalues of an unconstrained unit-cell stiffnessmatrix given in Figure 6(a). In the range 0 < ν < /
3, the eigenvalue λ of the rigidbody rotation eigenform is non-zero and therefore a part of the work performed by theexternal force is stored as strain energy due to this rigid-body motion. Only the remain-ing part of the work performed is available for the bending and shear eigenforms whichare also included in the cantilever bending solution. For ν > / ν approaches the value 1 / λ reduces to zero and thisimplies that the strain energy stored due to rigid-body rotation decreases. Thereforethe response of the plate employing Born model becomes relatively softer for ν = 0 . ν = 0. Due to the presence of a non-zero eigenvaluefor the rigid-body rotation eigenform, the solution also does not converge towards theanalytical solution. Since the modified model is able to distinguish between rigid-bodyrotations and shear deformations, no strain energy is stored in the case of rigid-body ro-tation. Therefore the unit-cell employing the modified bond model produces satisfactoryresults. Also the solution converges towards the analytical solution upon refinement.23 nalytical solution Born - 8 × × × ×
16 Modified - 8 × × × × − − − · − − − − · − u [m] y [ m ] (a) ν = 0 . . . . . − . − − . · − x [m] v [ m ] (b) ν = 0 − − − · − − − − · − u [m] y [ m ] (c) ν = 0 . . . . . . − . − − . · − x [m] v [ m ] (d) ν = 0 . − − − · − − − − · − u [m] y [ m ] (e) ν = 0 . . . . . . − . − − . · − x [m] v [ m ] (f) ν = 0 . Figure 18: Displacements u at the left edge of the plate (left) and v of the axis of theplate (right) employing the Born model and the modified model24
Conclusion
Modelling of continuum isotropic elasticity with particle methods such as the discreteelement method or the lattice spring method is traditionally limited to a certain valueof Poisson’s ratio due to the bond model used. The modified bond model presented inthis paper for planar continuum overcomes this limitation by introducing a multi-bondterm that couples the shear strain energy of neighbour bonds. This term enables themodified bond model to distinguish between rigid-body rotation and shear deformation.Two bonds, namely the
L-bond and
X-bond that employ this proposed coupling wereintroduced. The positive definiteness of the strain energy function of the unit-cell em-ploying the modified bond model was ensured for values of Poisson’s ratio in the range0 ≤ ν < .
5. The results obtained under uniaxial, pure bending, cantilever bending andpure shear loadings were validated with the continuum mechanics solution. Moreover,the agreement of the numerical results with the continuum mechanics solution demon-strate the ability of the modified bond model to describe the behaviour of an isotropicelastic material. The concept of shear strain energy coupling of neighbour bonds providesan alternate method to those existing in the literature which are based on computationof a particle stress or strain tensor. The generalisation of this modified bond model forrandom discretisation and for three-dimensions are under progress. The implementationof this modified bond model within the framework of the DEM is also to be carried out.
A Stiffness matrix of the unit-cell with
Born model
In the work of Griffiths [10], the stiffness matrix of a unit-cell was obtained by assembling the stiffness matrix ofindividual constituent bonds in a procedure that is similar to that of the Finite Element Method (FEM). Here analternative approach of deriving the stiffness matrix from the strain energy of the unit-cell is taken. The strainenergy stored in a generic bond in terms of the local displacement components was given in Equation (1). Thelocal displacements of particles in a generic bond oriented at an angle θ to the global coordinate system can bewritten in terms of the global displacement with the transformation matrix Q as u lo = Q u gl , (cid:20) u nA v sA (cid:21) = (cid:20) cos θ sin θ − sin θ cos θ (cid:21) (cid:20) u A v A (cid:21) and (cid:20) u nB v sB (cid:21) = (cid:20) cos θ sin θ − sin θ cos θ (cid:21) (cid:20) u B v B (cid:21) . (21)With this transformation, the strain energy in a bond in terms of the global particle displacements is given byΠ b = 12 k n (cos θ ( u B − u A ) + sin θ ( v B − v A )) + 12 k s (cos θ ( v B − v A ) + sin θ ( u A − u B )) . (22)The strain energy of the unit-cell in terms of displacements are obtained by summing up the strain energy ofindividual constituent bonds given in equation (22) by taking in to account their orientation asΠ uc = Π b ( θ = 0 ◦ ) + Π b ( θ = 90 ◦ ) + Π b ( θ = 180 ◦ ) + Π b ( θ = 270 ◦ )+ Π b ( θ = 45 ◦ ) + Π b ( θ = 135 ◦ ) . (23)While summing up, the stiffness factor of individual bonds have to be considered as well. The stiffness factor forbonds with first neighbours is 1 /
2, since they are shared by two unit-cells (periodicity). However, the stiffnessfactor for bonds with second neighbours is 1 as they belong exclusively to each unit-cell. By using the appropriatestiffness factors and by substituting the orientation of individual bonds, the stiffness matrix of the unit-cell isobtained by differentiating equation (23) twice with respect to the appropriate global displacement as uc = ˆ K ˆ K ˆ K ˆ K K K − ˆ K ˆ K ˆ K − ˆ K ˆ K K − ˆ K ˆ K K ˆ K ˆ K ˆ K K ˆ K ˆ K K K ˆ K ˆ K K − ˆ K ˆ K ˆ K = 12 k bn + 12 k bn + k bs ˆ K = 12 k bn − k bs ˆ K = − k bn ˆ K = − k bn − k bs ˆ K = − k bs . (24) B Stiffness matrix of the unit-cell with the modified bondmodel
For a generic
L-bond as shown in Figure 4(a), the strain energy stored in terms of the local displacements of theparticles can be written asΠ
ABC = 12 k mn (cid:20) ( u nB − u nA ) + ( u nC − u nB ) (cid:21) + 12 k ms (cid:20) − ( v sB − v sA ) + ( v sC − v sB ) (cid:21) . (25)With the transformation matrix Q and Equation (22), the strain energy is now given in terms of the globalparticle displacements asΠ ABC = 12 k mn (cid:20)(cid:18) ( c AB u B + s AB v B ) − ( c AB u A + s AB v A ) (cid:19) + (cid:18) ( c BC u C + s BC v C ) − ( c BC u B + s BC v B ) (cid:19) (cid:21) + 12 k ms (cid:20) − s BC ( u C − u B ) + c BC ( v C − v B ) + s AB ( u B − u A ) − c AB ( v B − v A ) (cid:21) (26)where, c AB = cos θ AB , s AB = sin θ AB , c BC = cos θ BC and s BC = sin θ BC . The strain energy stored in theunit-cell made up of L-bonds Π moduc is obtained by summing up the strain energy of individual constituent L-bond as shown in Figure 4(b) and in Equation (10). Similarly the strain energy stored in the unit-cell made up of
X-bond is given in terms of the local displacements asΠ moduc = 12 k mn (cid:20) ( u nC − u nA ) + ( u nD − u nB ) (cid:21) + 12 k ms (cid:20) ( v sC − v sA ) − ( v sD − v sB ) (cid:21) . (27)Similar to L-bond , the local displacements are expressed in terms of the global displacements with the transfor-mation matrix Q . After this, the strain energy stored in the unit-cell with the modified bond model is givenby Π moduc = Π moduc + Π moduc . (28)By differentiating the strain energy twice with respect to the appropriate global displacements, the stiffness matrixis given by K moduc = ˆ K ˆ K ˆ K ˆ K − ˆ K ˆ K ˆ K − ˆ K ˆ K ˆ K − ˆ K − ˆ K ˆ K − ˆ K ˆ K ˆ K ˆ K ˆ K ˆ K K ˆ K ˆ K ˆ K − ˆ K ˆ K − ˆ K ˆ K ˆ K ˆ K ˆ K − ˆ K ˆ K ˆ K = 12 k mn + 12 k mn + k ms ˆ K = 12 k mn − k ms ˆ K = − k mn − k ms ˆ K = − k ms ˆ K = − k mn − k ms . (29) eferences [1] Barber, J.: Elasticity. Solid Mechanics and Its Applications. Springer Netherlands (2009)[2] Born, M., Huang, K.: Dynamical theory of crystal lattices. Oxford classic texts in the physical sciences.Clarendon Press, Oxford (1954)[3] Buxton, G.A., Care, C.M., Cleaver, D.J.: A lattice spring model of heterogeneous materials with plasticity.Modelling and Simulation in Materials Science and Engineering (6), 485–497 (2001)[4] Capecchi, D., Ruta, G., Trovalusci, P.: From classical to voigt’s molecular models in elasticity. Archive forHistory of Exact Sciences (5), 525–559 (2010)[5] Celigueta, M.A., Latorre, S., Arrufat, F., O˜nate, E.: Accurate modelling of the elastic behavior of a contin-uum with the discrete element method. Computational Mechanics (6), 997–1010 (2017)[6] Cundall, P.A., Strack, O.D.L.: A discrete numerical model for granular assemblies. G´eotechnique (1),47–65 (1979)[7] Esin, M., Pasternak, E., Dyskin, A.V.: Stability of 2d discrete mass-spring systems with negative stiffnesssprings. physica status solidi (b) (7), 1395–1409 (2016)[8] Flack, C.: Mehrfeldmodellierung von beton mit diskreten element methoden. Ph.D. thesis, TU Braunschweig(2020)[9] Fung, Y.: Foundations of solid mechanics. Prentice-Hall, Englewood Cliffs N.J. (1965)[10] Griffiths, D.V., Mustoe, G.G.W.: Modelling of elastic continua using a grillage of structural elements basedon discrete element concepts. International Journal for Numerical Methods in Engineering (7), 1759–1775(2001)[11] Hassold, G.N., Srolovitz, D.J.: Brittle fracture in materials with random defects. Phys. Rev. B , 9273–9281(1989)[12] Keating, P.N.: Effect of invariance requirements on the elastic strain energy of crystals with application tothe diamond structure. Phys. Rev. , 637–645 (1966)[13] Kirkwood, J.G.: The skeletal modes of vibration of long chain molecules. The Journal of Chemical Physics (7), 506–509 (1939)[14] Love, A.: A Treatise on the Mathematical Theory of Elasticity. Cambridge University Press (2013)[15] Ockelmann, F., Dinkler, D.: A discrete element model for the investigation of the geometrically nonlinearbehaviour of solids. Computational Particle Mechanics (3), 335–344 (2017)[16] Ostoja-Starzewski, M.: Lattice models in micromechanics . Applied Mechanics Reviews (1), 35–60 (2002)[17] Potyondy, D., Cundall, P.: A bonded-particle model for rock. International Journal of Rock Mechanics andMining Sciences (8), 1329 – 1364 (2004)[18] ˇSmilauer, V., Chareyre, B.: Dem formulation. The Yade Project, in yade documentation 2nd ed. edn. (2015)[19] Timoshenko, S., Goodier, J.: Theory of Elasticity. Engineering mechanics series. McGraw-Hill (1969)[20] Zhao, C., Hobbs, B.E., Ord, A., Hornby, P., Peng, S., Liu, L.: Particle simulation of spontaneous crackgeneration problems in large-scale quasi-static systems. International Journal for Numerical Methods inEngineering (11), 2302–2329 (2007)[21] Zhao, G.F., Fang, J., Zhao, J.: A 3d distinct lattice spring model for elasticity and dynamic failure. Inter-national Journal for Numerical and Analytical Methods in Geomechanics (8), 859–885 (2011)(8), 859–885 (2011)