A modified bond model for describing isotropic linear elastic material behaviour with the particle method

Particle based methods such as the Discrete Element Method and the Lattice Spring Method may be used for describing the behaviour of isotropic linear elastic materials. However, the common bond models employed to describe the interaction between particles restrict the range of Poisson's ratio that can be represented. In this paper, to overcome the restriction, a modified bond model that includes the coupling of shear strain energy of neighbouring bonds is proposed. The coupling is described by a multi-bond term that enables the model to distinguish between shear deformations and rigid-body rotations. The positive definiteness of the strain energy function of the modified bond model is verified. To validate the model, uniaxial tension, pure shear, pure bending and cantilever bending tests are performed. Comparison of the particle displacements with continuum mechanics solution demonstrates the ability of the model to describe the behaviour of isotropic linear elastic material for values of Poisson's ratio in the range $0 \leq \nu<0.5$.


Introduction
Regarding the atomic scale a material is made up of atoms which are known to interact with each other through attractive and repulsive forces. Extrapolating this model to the macroscopic scale is the fundamental idea of modelling continua with interacting discrete particles. Referring to the historical development of the theory of elasticity [14,4] one realises that the mechanicians of the nineteenth century have made substantial contribution to the development of this alternate discrete model of elasticity. Apart from elasticity, it allows one to describe the failure of material as the absence of interaction between two previously interacting material points. Two of the most widely used particle methods are the Lattice Spring Model (LSM) and the discrete element method (DEM) initially developed for modelling the motion of granular media. The LSM influenced by 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 DEM initially developed by Cundall [6] has been extended to model solids such as rocks by the bonded-particle method (BPM) where two particles remain bonded together as long as a critical bond deformation is reached [17]. Although these methods are different in their implementation they are similar in the fact that they use bond models to describe the interaction between two particles. Several bond models exist in the literature and a brief summary of the most commonly used models is provided here. For a detailed review on the various bond models refer to [16]. The fundamental bond model is a Hookean spring which penalises change from the reference length between two interacting particles. For a solid discretised with an infinite number of material points connected by such springs, the Poisson's ratio is restricted to 1/3 in the case of plane stress and 1/4 for plane strain and in three-dimensions [10,16,21]. This model has one stiffness parameter to describe macroscopic isotropic elasticity which is on the other hand described by two independent material parameters. Therefore, additional crystal symmetry conditions called Cauchy symmetry are fulfilled and lead to the restriction on Poisson's ratio [14,4]. In order to overcome the restriction, Born [2] introduced a shear spring in addition to the longitudinal spring. The shear springs can be understood as a penalisation of the change in orientation from the reference 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 the range 0 ≤ ν ≤ 1 3 and 0 ≤ ν ≤ 1 4 for the case of plane stress and plane strain respectively. Apart from these restrictions of Poisson's ratio, the model cannot distinguish between rigid-body rotation and shear deformation [11,3,12]. Instead of using shear springs, the Kirkwood-Keating spring model [13,12] uses an angular spring to penalise the angular motion. However, this model is known to be nonlinear due to the angular terms and it offers no substantial advantage in comparison to the Born model [3]. The Lattice Beam model (LBM) uses a Hookean spring to allow longitudinal forces and a beam to allow shear force and bending moment and is widely used in DEM. However, in DEM the pure bending modes are neglected [18,10]. The beam model has the similar restriction on the Poisson's ratio as the Born model but is capable of distinguishing between rigid-body rotation 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 by these models, Zhao [21] evaluated the shear strains in the Born model using a particle strain tensor instead of particle displacements. This tensor is an approximate measure of the continuum strain tensor and is obtained by the least square method using the information of displacements of the particle under consideration and its neighbouring particles. The calculation of the particle strain tensor is computationally expensive and for a practical simulation with random particle arrangements there exists ambiguity regarding 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 two interacting particles in contact also the influence of forces in the surrounding of this contact is included. By means of this term a better description of continuum elasticity was obtained in comparison to the classical bond model of DEM. However, in order to calculate this term information regarding the contact area between two particles in contact and their surroundings are required. Additionally, the nonlocal stress tensor of each particle and its neighbours are required which results in a computationally expensive method.
In this work, a modification of the Born model is proposed by introducing a multibond strain energy term that is capable of distinguishing between shear deformation and rigid-body rotation for the case of two-dimensional elasticity. Two new multibonds called L-bond and X-bond employing the proposed coupling of shear strain energy are introduced. Based on positive definiteness of the strain energy function, we also interpret that the restriction of the Born model up to a certain Poisson's ratio results from incapability to distinguish between shear deformation and rigid-body rotation. With the modified model it is shown that the strain energy function remains positive definite for values of Poisson's ratio in the range 0 ≤ ν < 0.5. To validate the modified model, the results of a thin plate under uniaxial tension, pure shear, pure bending and cantilever bending are compared with the respective plane stress continuum mechanics solution.

Unit-cells with Born model
An arbitrary rectangular domain discretised with rigid circular particles of radius r is shown in Figure 1(a). Each particle is bonded with its immediate neighbours (separated by 2r) and also with its second neighbours (separated by 2 √ 2r). Similar to the works [10,15], the fundamental building block called as a unit-cell is extracted from the discretised domain and is shown in Figure 1(b). This square unit-cell is the combination of two unitcells: one which contains only the bonds with the immediate neighbours and another that contains only the bonds between second neighbours. They are called first-neighbour and second-neighbour unit-cells as shown in Figure 1(c) and 1(d) respectively. Although other configurations such as a triangle or a hexagon can be used as a unit-cell, the square configuration is used here due to its simplicity. The modification proposed in this paper is not restricted to a square configuration and it can be used also regarding other regular configurations. The Born model which is applied here employs springs in normal direction and tangential direction with stiffness parameters k n and k s as shown in Figure 2(a). A generic bond oriented at angle θ to the global coordinate system along with 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 of quantities with respect to the bond coordinate system as  After substituting the local strain quantities in terms of the components of the global strain tensor by the transformation relations provided by Griffiths [10] the strain energy in a bond is obtained as where, c = cos θ and s = sin θ. By substituting the orientation θ of the individual constituent bonds in Equation (3), the strain energy in the first and second neighbour unit-cells are obtained as where k b n 1 , k b s 1 are the normal and shear stiffness of first neighbour bonds and similarly k b n 2 , k b s 2 are the stiffness parameters of second neighbour bonds. All bonds in the same shear plane were assigned the same stiffness k b s 1 . Every bond in the first neighbour unitcell contributes only 1/2 to the strain energy of the unit-cell due to periodicity (every first neighbour bond is shared by two unit-cells). However the bonds in the second neighbour unit-cell belong completely to a single unit-cell. The strain energy density of the square unit-cell is given by where t is the thickness of the domain. As described before, the Born model cannot distinguish between rigid-body rotation and shear [12,11,21,3]. In Figure 3(a) an exploded view of each constituent bond of the unit-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 the final geometry describes pure rotation, the strain energy is not zero. Thus strain energy

Modified bond model
To overcome the limitations explained in Section 2, a modified model is proposed where the 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. Here again, the unit-cell made up of only the first neighbour bonds is looked at first and the unit-cell made up of the second neighbour bonds afterwards. 6

L-bond
For a generic L-bond ABC made up of two first neighbour bonds AB and BC as shown in 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 only for shear deformation and not for rigid-body rotation. The strain energy stored in a generic L-bond with length l AB = l BC = l is given by where k m n 1 and k m s 1 are the normal stiffness and shear stiffness parameters employing the modified model respectively. The first neighbour unit-cell is now made up of four L-bonds as shown in Figure 4 As explained in Section 2, a first neighbour bond has a contribution factor of 1/2. Every bond in the unit-cell is obtained by the combination of two L-bonds. For example, the Bond AB is obtained from DAB and ABC. Therefore, in an L-bond the normal stiffness k m n 1 and the shear stiffness k m s 1 have a contribution factor of 1/4. With the transformation relations provided by Griffiths [10], the local strain components are written in terms of the global strain components and the strain energy of a generic L-bond now where, c AB = cos θ AB , s AB = sin θ AB , c BC = cos θ BC and s BC = sin θ BC . Employing the modified bond model, the strain energy stored in the first neighbour unit-cell is obtained from the sum of strain energy stored in individual L-bonds The strain energy stored in the first neighbour unit-cell employing the modified bond model can be evaluated with respect to the global coordinate system to Comparing the formulation to that of the first neighbour unit-cell employing the Born model, see Equation (4), an extra term is added due to the coupling of the shear strain energies of the neighbours. With this term the modified model is able to distinguish between shear (change in the angle subtended) and rigid-body rotation (no change in the 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 modified such that strain energy is stored only in the case of shear deformations and not for pure rotation. This is achieved with the X-bonds. Again the shear strains of two neighbouring bonds are coupled and the normal strains of the bonds remain unchanged. The modified strain energy of the unit-cell made up of one X-bond is where, k m n 2 and k m s 2 are the normal and shear stiffness parameters employing the modified model respectively. The shear stiffness parameters in the same shear plane are assumed to have the same stiffness parameter k m s 1 . The strain energy stored in the second neighbour unit-cell with respect to the global coordinate system is Upon expansion of Equation (12) and comparison with the strain energy of the second neighbour unit-cell, see Equation (5), one observes that similar to the first neighbour unit-cell with L-bonds, the coupling of shear strain energies results in an extra multi-bond term to distinguish between shear and rotation.

Unit-cell with the modified model
The strain energy stored in a square unit-cell is obtained by summing up the strain energy of an unit-cell made with L-bonds, Equation (10), and of an unit-cell with Xbonds, Equation (12) Π mod uc = Π mod The strain energy density e mod uc is obtained similar to that of the unit-cell with the Born model given in Equation (6).

Comparison of the modified model and the Born model
In this section, the stiffness parameters of a square unit-cell employing the modified model 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 strain energy function of the unit-cell employing these bond models is verified for its positive definiteness.

Calibration of stiffness parameters
In order to model the behaviour of an isotropic elastic material, the stiffness parameters of the unit-cell employing the modified model must be calibrated with respect to the macroscopic elastic material parameters. In the literature there exist broadly two approaches: 1. A numerical calibration where the stiffness parameters are iteratively calibrated by solving an inverse problem to match the slope of the stress-strain curve of a material in the linear region under uniaxial loading. For more details refer [20,8].
2. An analytical calibration based on the equivalence of strain energy density of the unit-cell e uc with that of an equivalent elastic continuum e co = 1 2 σ : . In particular, the components of the elastic tensor C uc are derived from the strain energy density of the unit-cell and are compared to the continuum description.
In this work, the analytical calibration approach is used due to its independence of the fineness of discretisation. The components of the tensor of elasticity of the unitcell are obtained by differentiating the strain energy density twice with respect to the corresponding strain components. The individual components are summarised in the tensor of elasticity where, The elasticity tensor of the unit-cell has three components similar to that of a planar continuum. By comparing the three componentsĈ 1 ,Ĉ 2 andĈ 3 with those of the planar continuum elasticity tensor, the stiffness parameters are calibrated to the macroscopic Young's modulus E and the Poisson's ratio ν. The calibrated stiffness parameters for the case of plane stress and plane strain are summarised in Table 1.

Modified
Plane stress As expected for a Born model there exists a restriction on the Poisson's ratio at ν = 1/4 for plane strain and ν = 1/3 for plane stress. For higher values negative shear stiffness is obtained. This holds true for the modified model as well. The normal stiffness parameters of the unit-cell employing the Born model and the modified model (a) Born and modified model (b) Born and modified model  Figure 5 for the case of plane stress. The stiffness parameters are normalised with respect to the Young's modulus 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 the unit-cell agree exactly with the equivalent continuum components in the range 0 ≤ ν ≤ 0.5. In order to ensure the equality of the elasticity constants, the stiffness parameters of the unit-cell must take on the values as shown in Figure 5(a). Therefore, we interpret the negative value of shear stiffnesses (k b s 1 , k m s 1 ) for ν > 1/3 as a necessity to ensure this equality. Although a negative stiffness parameter of an individual spring seems unintuitive, in the case of a response of the unit-cell the elasticity tensor components are crucial.

Check for isotropy
Apart from calibrating the stiffness parameters to the macroscopic isotropic elastic material parameters it must also be checked whether and in which case they satisfy the condition for macroscopic isotropy. A material is considered to exhibit isotropy if the anisotropy factor Λ = 1. The anisotropic factor for the chosen unit-cell is calculated with the components of the elasticity tensor for both bond models as Modified model: This shows that for isotropic elasticity (Λ = 1) the unit-cell has only two independent stiffness parameters and a condition that ensures isotropy for all values of E and ν. By substituting the stiffness parameters from Table 1 in the condition for isotropy it can be shown that this condition holds for both plane stress and plane strain.

Positive definiteness of strain energy
The square unit-cell shown in Figure 1(b) is made up of four particles with two degrees of freedom per particle. The vector of displacements for one unit-cell is With the displacement vector u and the stiffness matrix K mod uc given in Appendix B, the strain energy of the unit-cell in Equation (13) can alternatively be written as Equation (17) shows that the strain energy of the unit-cell is a quadratic function of displacements. In order to be thermodynamically stable, it must be positive definite. The necessary and sufficient condition for the positive definiteness of Equation (17) is that all eigenvalues of the eigenvalue problem K mod uc u = λu are real and positive [9]. The relating stiffness matrices are given in Appendix A and B respectively. The resulting eigenvalues and corresponding eigenforms are summarised in Table 2. After substituting the stiffness parameters of the unit cell from Table 1, the non-zero eigenvalues are plotted as a function of Poisson's ratio ν for both plane stress and plane strain in Figure 6. The stiffness parameters are normalised with respect to Young's modulus E and the thickness t to obtain a qualitative behaviour.
From Figure 6 the following observations are made for both plane stress and plane strain: 1. Regarding the Born model, the eigenvalue λ 3 of the rigid-body rotation eigenform is non-zero everywhere expect at ν = 1/3 and ν = 1/4 for plane stress and plane strain respectively. This implies that the Born model deforms under rigid-body rotation and therefore it is proven that this model cannot distinguish between rigidbody rotation and shear deformations. Also, the eigenvalue becomes negative for ν > 1/3 and ν > 1/4 for plane stress and plane strain respectively. This implies that the strain energy function derived from the Born model is not positive definite after this lower limit of the Poisson's ratio. Therefore we expect the response to be unstable for values of Poisson's ratio above this limit. However, regarding the modified bond model, the eigenvalue of the rigid-body rotation eigenform remains zero for all values of ν and all eigenvalues remain positive. Hence, the strain energy derived from the modified bond model is positive definite and therefore is stable for all values of ν although the shear stiffness k m s 1 of the unit-cell takes on a negative value after a critical value of ν as shown in Figure 5(b). This implies that a negative stiffness does not necessarily lead to unstable results, rather it is 13 due to negative eigenvalues. Similarly, in the work of Esin [7] it was concluded that materials and structures with negative stiffness elements can exist when the negative element energy is compensated by the energy of the rest of the system or an encompassing system that provides stabilisation.
4. As expected, the volumetric form λ 8 has the highest eigenvalue.   The analytical solution for the displacement fields is u = σxx E x and v = − ν σxx E y. The resulting displacements along u(x = l, y) and v(x, y = l) for three different values of Poisson's ratio ν are summarised in the Table 3 for both the bond models.

Uniaxial
We observe that both unit-cells employing the Born model and the modified bond model produce the exact results for all values of the Poisson's ratio independent of the   The plate is constrained in both directions along the bottom edge to prevent rigidbody motion similar to that used by Ockelmann [15]. The analytical solution for the displacement fields is obtained as u = τ 0 G y and v = 0, where G = E/2(1 + ν) is the shear modulus. Here again, the test is conducted for three different values of ν and also for four different discretisations similar to the uniaxial test. The u and v displacements of the domain discretised with 8 × 8 unit-cells employing the Born model under pure shear are summarised in Figure 10(a), (b) and (c). In Figure 10(d), the results of 8 × 8 unit-cells employing the modified model under pure shear for ν = 0.49 are given. In Figure 10 the bonds connecting the particles are not visualised for clarity.
We observe that the results for the unit-cell employing the modified model agrees exactly with the analytical solution, because similar to the case of uniaxial tension, the analytical displacement field is a linear function of the position of particles. With the modified 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 stiffer in comparison with the analytical solution in the range 0 < ν ≤ 1/3 and unstable for ν > 1/3. In order to understand the reason for this behaviour, the plate is discretised with one unit-cell employing the Born model and the modified model respectively. The normalised eigenvalues of the constrained stiffness matrix are plotted as function of ν as shown in Figure 11. For the stiffness matrix of an unit-cell employing the Born model, we observe from Figure 11 that the eigenvalue λ b 5 decreases with increasing values of ν and also becomes negative for ν > 0.4. The eigenform corresponding to this eigenvalue is plotted 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 take on the form as previously shown in Figure 6 the case of plane stress. Therefore, the upper limit on the Poisson's ratio for the plate constrained as shown in Figure 9 is ν = 1/3 employing the Born model. However, for the unit-cell employing the modified model, due to the coupling of shear strain energy of neighbouring bonds, the eigenvalues remains positive for all values of ν as shown in Figure 11 for the constrained stiffness-matrix and in Figure 6   The test was performed with four different discretisations: 8 × 2, 16 × 4, 32 × 8 and 64 × 16. Also, the effect of ν on the accuracy has been checked by performing the test for 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 axis of the plate employing the Born model and the modified model are compared with the analytical solution for different discretisations and for different Poisson's ratio and are summarised 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 comparison with that of the analytical solution. In the case of pure bending, elements which are located faraway from the boundary undergo rigid body rotations as well. Since the Born model cannot distinguish between rigid-body rotations and shear deformations, strain energy is also stored for the rigid body rotation. As the rigid-body rotation eigenvalue takes on a maximum value at ν = 0 and approaches zero at ν = 1/3 for the chosen unit-cell employing the Born model, the response of the plate is much stiffer for ν = 0 than for ν = 1/3. And for ν > 1/3 the results are unstable since the eigenvalue takes a negative value and hence the strain energy function is no longer positive definite. We also observe regarding the unit-cell employing the Born model, that the displacements do not converge. 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 strain energy of the neighbouring bonds allows the model to distinguish between rigid-body rotations and shear deformations. Therefore, satisfactory results are obtained employing this modified bond model and also convergence to the analytical solution upon refinement is observed.

Cantilever bending
As a final test, a thin rectangular plate under cantilever bending is considered. The right edge of the plate is completely constrained and a constantly distributed load F is applied on the left edge as shown in Figure 16. The plate is modelled with the chosen unit-cell under plane stress conditions. From [1], the analytical solution for the displacement fields are obtained as For the derivation of the analytical solution the following weak boundary conditions were applied along the right edge (x = a): Similar to the pure bending test, four different discretisations are used. The effect of ν on the results is studied as well. The deformed configuration employing the Born model and the modified model are given in Figure 17.
The distributions of the displacement u at the left edge of the plate and the displacement v of the axis of the plate regarding the Born model and the modified model are summarised in Figure 18. From the results for ν < 1/3 we observe that the response of the plate discretised with the unit-cell employing the Born model is stiffer in comparison with that of the analytical solution. This is because, for the case of cantilever bending there exists elements faraway from the boundary that undergo rigid-body rotation as well. Since the Born model cannot distinguish between rigid-body rotations and shear  deformations, strain energy is stored also for rigid-body rotation. For ν > 1/3 the Born model produces unstable results. The reason for this behaviour can be understood by looking at the variation of normalised eigenvalues of an unconstrained unit-cell stiffness matrix given in Figure 6(a). In the range 0 < ν < 1/3, the eigenvalue λ 3 of the rigid body rotation eigenform is non-zero and therefore a part of the work performed by the external force is stored as strain energy due to this rigid-body motion. Only the remaining part of the work performed is available for the bending and shear eigenforms which are also included in the cantilever bending solution. For ν > 1/3 the eigenvalue of the rigid-body rotation eigenform becomes negative and hence the strain energy function is no longer positive definite. As ν approaches the value 1/3, λ 3 reduces to zero and this implies that the strain energy stored due to rigid-body rotation decreases. Therefore the response of the plate employing Born model becomes relatively softer for ν = 0.3 in comparison with the results for ν = 0. Due to the presence of a non-zero eigenvalue for the rigid-body rotation eigenform, the solution also does not converge towards the analytical solution. Since the modified model is able to distinguish between rigid-body rotations and shear deformations, no strain energy is stored in the case of rigid-body rotation. Therefore the unit-cell employing the modified bond model produces satisfactory results. Also the solution converges towards the analytical solution upon refinement.

Conclusion
Modelling of continuum isotropic elasticity with particle methods such as the discrete element method or the lattice spring method is traditionally limited to a certain value of Poisson's ratio due to the bond model used. The modified bond model presented in this paper for planar continuum overcomes this limitation by introducing a multi-bond term that couples the shear strain energy of neighbour bonds. This term enables the modified 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 were introduced. The positive definiteness of the strain energy function of the unit-cell employing the modified bond model was ensured for values of Poisson's ratio in the range 0 ≤ ν < 0.5. The results obtained under uniaxial, pure bending, cantilever bending and pure shear loadings were validated with the continuum mechanics solution. Moreover, the agreement of the numerical results with the continuum mechanics solution demonstrate the ability of the modified bond model to describe the behaviour of an isotropic elastic material. The concept of shear strain energy coupling of neighbour bonds provides an alternate method to those existing in the literature which are based on computation of a particle stress or strain tensor. The generalisation of this modified bond model for random discretisation and for three-dimensions are under progress. The implementation of 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 of individual constituent bonds in a procedure that is similar to that of the Finite Element Method (FEM). Here an alternative approach of deriving the stiffness matrix from the strain energy of the unit-cell is taken. The strain energy stored in a generic bond in terms of the local displacement components was given in Equation (1). The local displacements of particles in a generic bond oriented at an angle θ to the global coordinate system can be written in terms of the global displacement with the transformation matrix Q as u lo = Q u gl , With this transformation, the strain energy in a bond in terms of the global particle displacements is given by The strain energy of the unit-cell in terms of displacements are obtained by summing up the strain energy of individual constituent bonds given in equation (22) by taking in to account their orientation as While summing up, the stiffness factor of individual bonds have to be considered as well. The stiffness factor for bonds with first neighbours is 1/2, since they are shared by two unit-cells (periodicity). However, the stiffness factor for bonds with second neighbours is 1 as they belong exclusively to each unit-cell. By using the appropriate stiffness factors and by substituting the orientation of individual bonds, the stiffness matrix of the unit-cell is obtained by differentiating equation (23) twice with respect to the appropriate global displacement as B Stiffness matrix of the unit-cell with the modified bond model For a generic L-bond as shown in Figure 4(a), the strain energy stored in terms of the local displacements of the particles can be written as With the transformation matrix Q and Equation (22), the strain energy is now given in terms of the global particle displacements as where, c AB = cos θ AB , s AB = sin θ AB , c BC = cos θ BC and s BC = sin θ BC . The strain energy stored in the unit-cell made up of L-bonds Π mod uc 1 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 Π mod uc 2 = 1 2 k m n 2 (u n C − u n A ) 2 + (u n D − u n B ) 2 + Similar to L-bond, the local displacements are expressed in terms of the global displacements with the transformation matrix Q. After this, the strain energy stored in the unit-cell with the modified bond model is given by Π mod uc = Π mod uc 1 + Π mod uc 2 .
By differentiating the strain energy twice with respect to the appropriate global displacements, the stiffness matrix is given by