A General Numerical Method to Model Anisotropy in Discretized Bond-Based Peridynamics
AA Novel Numerical Method for Modeling Anisotropy in Discretized Bond-BasedPeridynamics
Naveen Prakash
Abstract
This work proposes a novel, general and robust method of determining bond micromoduli for anisotropiclinear elastic bond-based peridynamics. The problem of finding a discrete distribution of bond micromodulithat reproduces an anisotropic peridynamic sti ff ness tensor is cast as a least-squares problem. The proposednumerical method is able to find a distribution of bond micromoduli that is able to exactly reproduce adesired anisotropic sti ff ness tensor provided conditions of Cauchy’s relations are met. Examples of all eightpossible elastic material symmetries, from triclinic to isotropic are given and discussed in depth. Parametricstudies are conducted to demonstrate that the numerical method is robust enough to handle a variety ofhorizon sizes, neighborhood shapes, influence functions and lattice rotation e ff ects. Finally, an exampleproblem is presented to demonstrate that the proposed method is physically sound and that the solutionagrees with the analytical solution from classical elasticity. The proposed method has great potential formodeling of deformation and fracture in anisotropic materials with bond-based peridynamics. Keywords:
Peridynamics, Anisotropy, Bond-based Models, Least-Squares
1. Introduction and Motivation
Peridynamics modeling has been applied to a wide range of problems since it was first introduced bySilling [1]. The advantages of peridynamics, e.g. the ability to use a meshfree type of discretization to solvethe equations of motion have allowed the modeling of fracture initiation and propagation with relative ease[2]. Peridynamics has also been extended to include many multiphysics phenomena - heat conduction [3, 4],electrical conduction [5, 6], fluid flow [7] and corrosion [8]. Although there has been a significant amountof work in ordinary[9–11] and non-ordinary state-based peridynamics[12–17], bond-based peridynamicsstill attracts a lot of attention in literature.One of the areas of interest has been to include anisotropic response within the framework of Peridy-namics. For example, Hu and co-workers used bond-based peridynamics to model fracture in uni-directionfiber-reinforced composites using peridynamic bonds of di ff erent micromoduli in di ff erent directions to in-troduce anisotropy [18, 19]. While Hu and co-workers performed a 2D analysis, Oterkus and Madenci [20]used the same approach for a 3D analysis of composite laminate with multiple plies. Hu et al. used a similarapproach to model delamination growth and fatigue in composite structures [21, 22]. Kilic and Madenci onthe other hand modeled the fiber and matrix phases of a composite structure and retained the inhomogeneousnature of the material [23]. In more recent work, Mehrmashhadi et al. have presented a stochastic approach Email addresses:
[email protected] (Naveen Prakash), [email protected] (Naveen Prakash) a r X i v : . [ c s . C E ] N ov o modeling fiber reinforced composites using bond-based peridynamics. Another recent paper by Dianaand Casolo used the micropolar version of bond-based peridynamics to model orthotropic materials [24].Ghajari et al. has presented more fundamental work in using spherical harmonics expansion to representthe micromodulus as a function of the bond orientation [25]. A similar approach was followed in Zhou etal. [26] as well as in Azdoud et al. [27] for orthotropic and transversely isotropic materials. Trageser et al.have presented a detailed analysis of anisotropic bond-based peridynamics models in [28] with a specificfocus on 2D plane stress and plane strain models. A similar approach has been taken outside of the peridy-namics method as well with granular micromechanics in [29, 30] and virtual internal bond method by Gaoand Klein [31].Most approaches to model anisotropy so far have used either direct numerical simulations using di ff erentmicromoduli for di ff erent phases to induce anisotropy in a global sense, or made the micromoduli a functionof bond orientation and used a functional form for the micromoduli in terms of the bond orientation to thematerial’s axis to express this relationship. Many of these approaches start with some a priori assumptionabout the distribution of micromoduli, e.g. all bond micromoduli are constant in the case of isotropy ormicromoduli in one direction are sti ff er than other directions in case of orthotropy.This work presents a novel and general approach to calibrating micromoduli in terms for discretized neighborhoods for anisotropic bond-based peridynamics - without any assumptions on the distribution ofbond micromoduli a priori. Using a discretized neighborhood or a lattice as a starting point for peridynamicanalysis has not been commonly adopted yet. Liu and Hong adopted a discretized bond-based approach in[32] and [33] where as Nikravesh and Gerstle presented a state-based peridynamic lattice model (SPLM)using a hexagonal lattice as a starting point[34]. Gerstle argues that a discretized peridynamics might bemore realistic and computationally e ffi cient [35]. Again, discretized lattice particle methods have existedoutside of peridynamics as well, using both local and non-local models e.g. see [36–43]The main premise of the proposed method is that a collection of bonds in the neighborhood of a peri-dynamic material particle with a specific distribution of bond sti ff nesess or micromoduli results in some ef-fective peridynamic sti ff ness tensor. For example, if all micromoduli are equal, it will result in an isotropicsymmetry (note that this is not the only distribution that may lead to isotropic symmetry). The questionthen is, can one evaluate a collection of micromoduli such that they result in some desired arbitrary sti ff nesstensor?Most published literature deals with specific material classes such as orthotropy, but this work is afirst in that it presents a general, robust method for modeling any type of anisotropy in 3D bond-basedperidynamic models. No ad-hoc treatment is required for any specific material symmetry. The errorsaccrued in the calibration process are unambiguous and it is shown that if conditions of Cauchy’s relationsare satisfied, the bond micromoduli can be calibrated such that the e ff ective peridynamic sti ff ness is exactlyequal to the desired sti ff ness tensor regardless of the class of symmetry. Another significant advantage isthat this method takes in account volume corrections as well [44] thus removing the need for additionalcomputations. Although this method has the potential to address surface corrections in the future as well,this is not considered in the present work.Restrictions for each of the eight material classes arising from Cauchy’s relations and examples aregiven in section 4.1. More practical issues such as the e ff ect of horizon size, influence function, shape ofthe neighborhood and rotation of the neighborhood are also discussed with examples. Section 4.4 finallypresents an example simulation of a uniaxial tensile test with bond micromoduli distribution obtained fromthe proposed method. It is demonstrated that the proposed method is indeed mathematically and physicallysound, and the results are shown to agree well with the analytical solution.2 . Peridynamics Preliminaries Peridynamics is a continuum theory in which it is assumed that every material particle x interacts withevery other particle x (cid:48) within some interaction distance called the horizon H x through a peridynamic bond .The finite volume surrounding x within which these interactions take place is called the neighborhood of theparticle. Consider two particles x and x (cid:48) in the reference configuration, such that their undeformed relativeposition is given by ξ = x (cid:48) − x . If x and x (cid:48) are displaced by u and u (cid:48) respectively, their relative displacementis given by η = u (cid:48) − u such that their relative position in the deformed configuration is given by ξ + η . Corning Restricted 𝜹𝒙𝒙′
Materialparticles
Horizon
Body O ො𝒆 ො𝒆 ො𝒆 𝛏 𝛏 + 𝛈 𝒖𝒖′ Current configuration
Reference configuration
Body 𝑯 𝒙 Figure 1: Kinematics of peridynamic material particles, material particle x and x (cid:48) correspond to an interacting pair of particleswithin the horizon region of x , having a volume dV x (cid:48) [45]. The peridynamic equations of motion for a particle x at time t are given by, ρ ¨ u ( x , t ) = (cid:90) H x f ( η , ξ , t ) dV x (cid:48) + b ( x , t ) . (1)where ρ and ¨ u denote the density and acceleration of the material particle x , f is the pairwise force function(units of force per unit volume squared) of the bond between x and x (cid:48) and b is the body force per unit volumeat x . The net internal force per unit volume arising from non-local pairwise interactions between particlesis obtained from the integral of f over the horizon H x . Since the formulation does not involve spatialderivatives of displacement, the same governing equations can be applied in the presence of discontinuities.According to Silling [1], a micro-elastic material is defined as one in which the internal force is thegradient of a scalar potential w ( η , ξ ) i.e., f ( η , ξ , t ) = ∂ w ( η , ξ ) ∂ η . (2)For an elastic material, the scalar potential refers to a micro-elastic strain energy, i.e. strain energydensity per unit volume stored in the material. The material is said to be linear micro-elastic[1] if w ischosen to be, w ( ξ, η ) = c ( ξ ) s ( η , ξ ) | ξ | , (3)3here c ( ξ ) is a constant called the micromodulus, | ξ | is the magnitude of the reference bond vector ξ and sis the stretch of the bond which is given by, s = | ξ + η | − | ξ || ξ | . (4)The micromodulus c ( ξ ) can depend on the reference bond vector ξ , however for homogeneous isotropicmaterials the micromodulus is generally assumed to be independent of the direction of the bond. For thecurrent purposes, this dependence is retained. If Eq. (3) is substituted in Eq. (2), the expression for theinternal force density in a peridynamic bond for a homogeneous linear isotropic material can be derivedsuch that the internal force f is given by, f (cid:0) η, ξ (cid:1) = c ( ξ ) s ξ + η | ξ + η | , (5)where the force in the bond is a function of the stretch of the bond s , the micromodulus c ( ξ ), and acts in thedirection of the deformed bond as indicated by the unit vector (cid:0) ξ + η (cid:1) / (cid:0) | ξ + η | (cid:1) . This is commonly knownas the linear microelastic model of peridynamics and the reader is referred to Silling[1, 46] for a detailedderivation of f and Eq. (5). In this general form, bond stretch is a nonlinear function of u (cid:48) − u , the relativedisplacement of the bond and hence the force - displacement relationship is nonlinear.This constitutive model can be linearized by assuming small deformations, i.e. | η | (cid:28)
1. Rewriting Eq.(5), f (cid:0) η, ξ (cid:1) = c ( ξ ) (cid:32) | ξ + η | − | ξ || ξ | (cid:33) ξ + η | ξ + η | (6) ⇒ f (cid:0) η, ξ (cid:1) = c ( ξ ) (cid:32) | ξ | − | ξ + η | (cid:33) ξ + η (7)For | η | (cid:28)
1, using Taylor series expansion about η = and ignoring higher order terms of η , f (cid:0) η, ξ (cid:1) = f ( , ξ ) + ∂ f ∂ η (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) η = η (8)where the first term goes to zero, meaning that the internal force in the bond is zero in the undeformedconfiguration. Substituting for f in the second term of Eq. (8) from Eq. (7), f (cid:0) η, ξ (cid:1) = (cid:34) c ( ξ ) ∂∂ η (cid:32) | ξ | − | ξ + η | (cid:33) ⊗ (cid:0) ξ + η (cid:1)(cid:35) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) η = η + c ( ξ ) (cid:34)(cid:32) | ξ | − | ξ + η | (cid:33) ∂ (cid:0) ξ + η (cid:1) ∂ η (cid:35) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) η = η (9)The second term in Eq. (9) when evaluated at η = is zero and Eq. (9) becomes, f (cid:0) η, ξ (cid:1) = (cid:34) c ( ξ ) ∂∂ η (cid:32) − | ξ + η | (cid:33) ⊗ (cid:0) ξ + η (cid:1)(cid:35) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) η = η (10)where the partial derivative of 1 / | ξ | with respect to η is zero. Evaluating Eq. (10) and substituting η = , f (cid:0) η, ξ (cid:1) = (cid:34) c ( ξ ) ξ ⊗ ξ | ξ | (cid:35) η (11)where the quantity in brackets is a second order tensor with components that are a function of the micro-modulus c ( ξ ) and ξ . Rewriting Eq. (11), f (cid:0) η, ξ (cid:1) = (cid:34) c ( ξ ) ξ | ξ | ⊗ ξ | ξ | (cid:35) η | ξ | (12)4nd using the tensor product rule ( a ⊗ b ) c = ( b . c ) a , Eq. (12) can be written as, f (cid:0) η, ξ (cid:1) = c ( ξ ) ξ | ξ | · η | ξ | ξ | ξ | (13)Recognizing that ξ / | ξ | is nothing but the unit vector in the undeformed configuration of the bond andthat η · ξ / | ξ | is the component of relative displacement η in the direction of the undeformed bond, Eq. (13)can be written as , f (cid:0) η, ξ (cid:1) = c ( ξ ) η n | ξ | ξ | ξ | , (14)where η n denotes the component of relative displacement in the direction of the undeformed bond. Com-paring Eq. (14) to Eq. (4), it is found that the term η n / | ξ | is the linearized stretch s and the direction of theforce is changed from the deformed bond direction to the undeformed bond direction. Substituting for theinternal force density from Eq. (14) in Eq. (1), the peridynamic linear momentum equation with a linearizedforce - displacement relationship becomes, ρ ¨ u ( x , t ) = (cid:90) H x c ( ξ ) η n | ξ | ξ | ξ | dV x (cid:48) + b ( x , t ) . (15)The corresponding strain energy density for the linearized version of the linear microelastic model isgiven by, w ( ξ, η ) = c ( ξ ) η n | ξ | (16)The strain energy density stored at a material point can be written as the integral over all bonds connectedto the material point, W PD ( x ) = (cid:90) H x c ( ξ ) η n | ξ | dV x (cid:48) . (17)where the 1 / c ( ξ ) = c . Then the micromodulus c inEq. (15) can then be evaluated by comparing the stored peridynamic strain energy density to that obtainedusing classical continuum mechanics principles. The expression for the micromodulus c can be derived in3D as [2], c D = K πδ . (18)where, K is the bulk modulus of the material. Similar expressions for the micromodulus can be derived for2D plane stress and plane strain conditions as well.The general linearized bond-based model, with the strain energy density as shown in Eq. (17) is takenas the starting point for the proposed numerical method going forward and is used to derive the e ff ectiveperidynamic sti ff ness tensor as shown in the next section.
3. Least-Squares Method to Calibrating Micromoduli
Assume that under small strains and displacements, the macroelastic strain energy density under ahomogeneous strain field for linear elastic material is given by, W PD ( x ) = (cid:90) H x ω ( | ξ | ) c ( ξ ) η n | ξ | dV x (cid:48) , (19)5hich is the same as Eq. (17) but with an additional influence function ω ( | ξ | ) that weights contributions ofindividual bonds to the strain energy stored based on their undeformed bond lengths. Note that this equationcan be derived similarly starting with the inclusion of the influence function in Eq. (3), however is not doneto keep the equations brief.Assume that the relative displacements of the bonds η n is cause by some imposed homogeneous strainfield ε . Taking the derivative of Eq. (19) with respect to the second order strain tensor, ∂ W PD ∂ε kl = (cid:90) H x ω ( | ξ | ) c ( ξ ) η n | ξ | ∂η n ∂ε kl dV x (cid:48) , (20)where the scalar relative displacement can be written in terms of the bond length and the strain field as, η n = ξ i ε i j ξ j | ξ | , (21)and the derivative can be written as, ∂η n ∂ε kl = ξ k ξ l | ξ | . (22)Substituting into Eq. (20) and taking the derivative again, ∂ W PD ∂ε i j ∂ε kl = (cid:90) H x ω ( | ξ | ) c ( ξ ) ξ i ξ j ξ k ξ l | ξ | dV x (cid:48) . (23)which is nothing but the fourth order sti ff ness tensor, C PDi jkl = (cid:90) H x ω ( | ξ | ) c ( ξ ) ξ i ξ j ξ k ξ l | ξ | dV x (cid:48) , (24)or in vector notation can be written as, C PD = (cid:90) H x ω ( | ξ | ) c ( ξ ) ξ ⊗ ξ ⊗ ξ ⊗ ξ | ξ | dV x (cid:48) , (25)where the integral is over the horizon of particle x which covers an infinite number of it’s neighbors x (cid:48) andthe range of indices i jkl are from 1 to 3. Observe from the form of Eq. (24) and Eq. (25) that C PD hasthe following symmetries commonly known as the minor and major symmetries respectively which are alsopresent in the classical form of the fourth order material sti ff ness tensor, C PDi jkl = C PDjikl = C PDi jlk , (26)and, C PDi jkl = C PDkli j . (27)In addition, note from Eq. (25) that the peridynamics sti ff ness tensor also possesses an additional symmetryof the inner indices j and k , also known as Cauchy’s relations or Cauchy’s symmetry[47–49] , C PDi jkl = C PDik jl (28) We refer to this as Cauchy’s relations henceforth ff ness tensor fully symmetric. In other words, any two indices in the peridynamic sti ff nesstensor may be interchanged having no e ff ect on the sti ff ness tensor. There are no other restrictions on theperidynamic sti ff ness tensor other than those given in Eq. (26), Eq. (27) and Eq. (28), i.e. any material notnecessarily just isotropic can be represented using Eq. (25) provided the necessary symmetry conditions aremet.The discrete version for a discretized neighborhood at some particle in the domain can be written byreducing the integral to a finite summation following the method of Silling and Askari [2] in index notationas, C PDi jkl = M (cid:88) N = ω N c N ξ Ni ξ Nj ξ Nk ξ Nl (cid:12)(cid:12)(cid:12) ξ N (cid:12)(cid:12)(cid:12) ∆ V N . (29)This expression is written for a particle x as a summation over it’s Nth neighbor where N varies from 1 toM where M is the total number of neighbors, and therefore the total number of bonds. For example, thesti ff ness term C PD can be written as, C PD = M (cid:88) N = ω N c N ξ N ξ N ξ N ξ N (cid:12)(cid:12)(cid:12) ξ N (cid:12)(cid:12)(cid:12) ∆ V N . (30)Note that quantities for each bond appearing in Eq. (29) - e.g. ω N , c N etc. are identified with a superscript N except for the volume which is written with an N in the subscript due to convention.Voigt notation can be used to simplify, where by taking advantage of minor and major symmetries, thefourth order sti ff ness tensor can be written as a 6 × C PD αβ = M (cid:88) N = ω N c N ζ N α ζ N β (cid:12)(cid:12)(cid:12) ξ N (cid:12)(cid:12)(cid:12) ∆ V N , (31)where α and β are the two Voigt indices that vary from 1 to 6 such that they relate to cartesian indices i jkl as: 1 →
11, 2 →
22, 3 →
33, 4 →
23, 5 →
31 and 6 →
12. For example C in Voigt notation corresponds to C in conventional notation, or C due to minor symmetry. Similarly, ζ = ξ ξ , ζ = ξ ξ , ζ = ξ ξ , ζ = ξ ξ , ζ = ξ ξ and ζ = ξ ξ .For each one of the 21 independent constants, (31) can be written as a dot product of two vectors,˜ C PD αβ = M (cid:88) N = ω N ζ N α ζ N β (cid:12)(cid:12)(cid:12) ξ N (cid:12)(cid:12)(cid:12) ∆ V N . c N , (32)where the length of the vectors are equal to the number of neighbors or the number of bond micromoduli tobe evaluated - M. For example, for the 11 component, (cid:104) ˜ C PD (cid:105) = (cid:20) ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V · · · ω M ζ M ζ M | ξ M | ∆ V M (cid:21) c c c ... c M . (33) The collection of bonds here is referred to as the ‘neighborhood’, some works also refer to this as the ‘family’ of the particle.
7f the 21 independent sti ff ness components ˜ C PD αβ are arranged as components of a vector [ ˜ C PD ˜ C PD ˜ C PD ... ˜ C PD ] T then, ˜ C PD ˜ C PD ˜ C PD ... ˜ C PD = ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V · · · ω M ζ M ζ M | ξ M | ∆ V M ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V · · · ω M ζ M ζ M | ξ M | ∆ V M ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V · · · ω M ζ M ζ M | ξ M | ∆ V M ... ... ... ... ... ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V ω ζ ζ | ξ | ∆ V · · · ω M ζ M ζ M | ξ M | ∆ V M c c c ... c M , (34)or in matrix notation can be written as, ˜ C PD = X c , (35)where tilde is used to denote that the sti ff ness components are written in vector form (containing terms inthe voigt sti ff ness matrix).If the desired elastic sti ff ness of the material or the ‘reference’ sti ff ness is similarly denoted by ˜ C re f ,then ideally, the aim would be to find a solution c such that ˜ C PD = ˜ C re f . In general, for all practicalpurposes this coe ffi cient matrix is rectangular and wide. Therefore we seek to find a solution c that has theminimum error in a least-squares sense, that is, c = c : c i ≥ , min (cid:13)(cid:13)(cid:13) ˜ C re f − X c (cid:13)(cid:13)(cid:13) . (36) However, the least-squares solution does not guarantee non-negativity of all components of the solution c .This requirement of non-negativity of the micromoduli stems from the fact that the stored strain energydensity of the bond given in Eq. (16) has to be non-negative, w ( ξ, η ) ≥ c ≥
0. Therefore, an additional constraint needs to be imposed - all components of the solution vector c have to be non-negative.One way to find a unique solution to the problem is find that solution to the least-squares problemwhich also gives the lowest possible values of micromoduli. In other words, to find the minimum normleast-squares solution where in addition to minimizing the residual given by (cid:13)(cid:13)(cid:13) ˜ C re f − X c (cid:13)(cid:13)(cid:13) , the norm of thesolution (cid:107) c (cid:107) itself is minimized. This can be computed from the Moore-Penrose pseudoinverse of X (ina software such as MATLAB ® , the functions pinv or lsqminnorm can be used to compute this solution).However, once again, this does not guarantee non-negativity of the solution c . Therefore to find the solu-tion with the least norm and least residual within the space of positive solutions, the following scheme isproposed.Let v be a vector of size M × X such that, X v = . (37)where, X is the coe ffi cient matrix in Eq. (35). Then the desired solution is the vector c + v such that thenorm (cid:107) c + v (cid:107) is the minimum among all possible non-negative solutions, c + v ≥ . In other words, thedesired solution is obtained by a two-step process:1. First, evaluate the minimum norm least-squares solution c .8. Second, this solution c is perturbed by an amount v such that c + v is non-negative and such that theresidual (cid:13)(cid:13)(cid:13) ˜ C re f − X( c + v ) (cid:13)(cid:13)(cid:13) remains unchanged since v lies in the null space of X .This scheme can be posed as a quadratic program subject to linear equality and inequality constraintsas follows. Finding min (cid:107) c + v (cid:107) is equivalent to finding the min ( c + v ) T ( c + v ) which can be expanded as, min ( c + v ) T ( c + v ) = min (cid:16) c T c + v T c + c T v + v T v (cid:17) , (38) = min v T v + c T v , (39) = min v T Iv + c T v , (40)where I is an M × M identity matrix and the term c T c is dropped since it is a known positive constant anddoes not a ff ect the minimum. Therefore, the quadratic program can be now written as, Minimize v T Iv + c T v ; (41) S ub ject to c + v ≥ , (42) X v = . (43)Note that in the first step, other routines such as lsqlin and lsqnonneg can also be used to obtain theoriginal solution c , which may return one of multiple solutions for c . For example, using lsqnonneg returnsa solution of non-negative values in the solution vector but not necessarily of minimum norm. However,this is of no consequence as the second step aims to find the minimum norm solution ( c + v ) which is castas the objective function to be minimized in the quadratic program. The minor advantage of using say, lsqminnorm in the first step is that if for a particular case, this routine finds the minimum norm solutionwith all non-negative micromoduli, then that is the desired solution and the need for the second step iseliminated.Although a formal proof of existence of solutions is not attempted here, as will be seen in the forthcom-ing results, solutions are found for a very wide variety of practical cases thus demonstrating the generalityand robustness of the proposed method.
4. Results and Discussion
Results are given here to demonstrate the feasibility of this approach in producing a distribution of bondmicromoduli for anisotropic bond-based models. Anisotropy is generally defined by a symmetry group for amaterial, a group of orthogonal transformations under which the fourth order elasticity tensor for a materialis invariant. Simply put, suppose that a material is defined by a fourth order elasticity tensor C i jkl , such thatif under an orthogonal transformation defined by a matrix Q i j , C (cid:48) pqrs = Q pi Q q j Q rk Q sl C i jkl , (44)if C (cid:48) pqrs = C i jkl , then Q i j is said to be a symmetry transformation of C i jkl . In addition, a material is said topossess a symmetry plane if there exists a reflection transformation that satisfies Eq. (44) such that Q i j canbe defined as [50], Q i j = R ( n ) = I i j − n i n j (45)9here n is the unit normal to the plane of reflection. It turns out that all eight material symmetries can beobtained by successive introduction of symmetry planes starting from the most anisotropic - triclinic to theleast - isotropic [51, 52]. The group S of all possible orthogonal transformations Q i j for a material underwhich the material’s elasticity tensor remains invariant is known as the symmetry group for that material.Note that all anisotropic materials admit the identity transformation I and the central inversion transfor-mation − I [50] where, I = , − I = − − − . (46)The peridynamic sti ff ness tensor Eq. (25) of course admits these transformations as well. Notably, theinversion transformation can be written from Eq. (25) and Eq. (44) as, C PD = (cid:90) H x ω ( | ξ | ) c ξ ⊗ ξ ⊗ ξ ⊗ ξ | ξ | dV x (cid:48) = (cid:90) H x ω ( | ξ | ) c − ξ ⊗ − ξ ⊗ − ξ ⊗ − ξ | ξ | dV x (cid:48) . (47)In other words, the contribution of individual bonds to the peridynamics sti ff ness tensor is invariant if thebond is reflected about it’s own axis and therefore the micromoduli c ( ξ ) = c ( − ξ ). This can be observed forall cases in the forthcoming results.All eight possible symmetries for linear elastic materials are tested with material properties taken from[28] (except for triclinic which is obtained from [53]). Not all of these materials satisfy Cauchy’s relationsexactly, however are still used as inputs as is. This is done deliberately so that firstly, real materials may beused as examples, secondly so that errors in the calibration may be estimated when using realistic materialproperties which may not satisfy exactly Cauchy’s relations and finally so that the reader may observe thatCauchy’s relations emerge naturally in the calibrated elastic sti ff ness tensor.It is important to note that while any material sti ff ness tensor can be used to derive bond micromod-uli distributions using the proposed method, it is most e ff ective when using materials that exactly satisfyCauchy’s relations. This is of course because the bond-based peridynamics model itself, due to the assump-tion of pairwise interactions, generates an elasticity tensor which inherently satisfies Cauchy’s relations.The computational least-squares method is implemented in the software, MATLAB ® . The discretizedperidynamic neighborhood or ‘lattice’ is assumed to be a 3D regular grid with each particle imagined tobe a cube of unit size, unit volume, and therefore has unit lattice spacing ( ∆ = { e , e , e } unless otherwise stated. A spherical neighborhood with a horizon of 6 times the lattice spacing (i.e δ = ∆ = ω ( | ξ | ) = δ | ξ | is used for all examples unless otherwisestated. All elastic constants have units of GPa. For the sake of completeness, the degree of anisotropyfor each case is also evaluated and presented in the form of the universal anisotropy index [53, 54] (seeAppendix A). It is well known that there are eight possible symmetries for linear elastic material - triclinic, monoclinic,orthotropic, trigonal, tetragonal, transversely isotropic, cubic and isotropic [28, 50, 51]. Presented below isan example of each case starting with triclinic to isotropic in relatively decreasing order of anisotropy. Coordinates are also referred to as x,y,z in the text for ease and brevity .1.1. Triclinic Symmetry Triclinic, the most general material class has no planes of symmetry, i.e. all 21 constants in the uppertriangle of the Voigt sti ff ness matrix can be non-zero and are independent. The symmetry group for triclinicmaterials can be written quite simply as, S = { I , − I } . (48)With Cauchy’s restrictions imposed, the number of independent constants for a triclinic material reducesfrom 21 to 15, and can be visualized as follows, C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C → C (cid:12) ♦ ⊗ C C (cid:12) C (cid:3) C (cid:9) C ♦ (cid:3) C C C (cid:63) ⊗ C C (cid:3) (cid:63) (cid:9) C (cid:9) C (cid:63) ♦ ⊗ C C (cid:63) (cid:9) ⊗ (cid:12) . (49)Entries containing the same symbol are equal to each other. The other entries are left as is to di ff erentiatethem from the ones a ff ected by Cauchy’s relations. Note that for a triclinic material, 9 elastic constantsremain una ff ected while 6 out of the other 12 are left truly independent after Cauchy’s relations are imposed,leading to a total of 15 independent elastic constants. An example of a triclinic material, KIO , is givenbelow with properties are taken from [53]. The Voigt elastic sti ff ness matrix is given by,˜ C re f =
43 11 13 1 − − − − − − . (50)It appears that some terms a priori satisfy Cauchy’s relations exactly, for example C = C and C = C where as others are quite close resulting in an anisotropy coe ffi cient A U of 0.1965. (a) (b) (c) Figure 2: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for triclinic symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . The bond micromoduli are solved for using the least-squares method described in section 3. For thepurposes of illustration, the particles are depicted as colored spheres in Figure 2 with the color depicting11he magnitude of bond micromoduli between a particle within the horizon in the lattice and the referenceparticle at the center of the lattice. Figure 2a shows the complete neighborhood, Figures 2b and 2c showspartial views of the neighborhood with normal to the cut plane at n = − e and n = e respectively so thatthe reader may observe the distribution of bonds inside the volume of the neighborhood better. True to thetriclinic material class, there appear to be no symmetries present in the peridynamic neighborhood from avisual inspection of Figure 2. Bonds oriented in the x and z direction appear to be of higher sti ff ness thanthose oriented in the y direction in general, which appears to be qualitatively correct from looking at Eq.(50).The calibrated peridynamic sti ff ness tensor which can be obtained by ˜ C PD = X ( c + v ) is presented inVoigt notation in Eq. (51).˜ C PD =
43 11 .
67 13 1 − .
67 35 12 .
67 3 − .
33 313 12 .
67 43 2 − .
331 3 2 12 .
67 0 . − . − − . − .
33 13 12 3 0 . − .
33 1 11 . . (51)The 9 entries which remain una ff ected by Cauchy’s relations, namely C , C , C , C , C , C , C , C , C are left unchanged during the least-squares calibration and the other 12 are modified (compare thestructure of the matrices given in Eq. (51) and Eq. (49)).Note that the reference sti ff ness matrix in Eq. (50) is used as input to the model as is, and not modifiedin any way to satisfy Cauchy’s relations a priori. Cauchy’s relations are found to emerge naturally aspart of the least-squares calibration procedure. The relative error in the solution, which is defined here as (cid:13)(cid:13)(cid:13) ˜ C re f − ˜ C PD (cid:13)(cid:13)(cid:13) / (cid:13)(cid:13)(cid:13) ˜ C re f (cid:13)(cid:13)(cid:13) is 3.1873% where (cid:107)·(cid:107) represents the L norm. Monoclinic symmetry has a single plane of symmetry, which in the present case is the x − y plane or z = S = { I , − I , R ( e ) } . (52)This gives rise to 13 independent constants namely, C , C , C , C , C , C , C , C , C , C , C , C , C . (53)Cauchy’s relations reduce the number of independent constants from 13 to 9, C , C , C , C = C , C = C , C = C , C , C , C = C , (54)and can be visualized as follows, C C C C C C C C C C C C C C
00 0 0 C C C C C C → C (cid:12) ♦ C (cid:12) C (cid:3) C ♦ (cid:3) C (cid:63) (cid:3) (cid:63)
00 0 0 (cid:63) ♦ C C (cid:63) (cid:12) . (55)12he example material chosen is CoTeO , with a universal anisotropy index of 28.3169 and Voigt elasticsti ff ness matrix given by, ˜ C re f =
135 19 54 0 0 4219 13 15 0 0 654 15 269 0 0 180 0 0 14 25 00 0 0 25 66 042 6 18 0 0 18 . (56) (a) (b) (c) Figure 3: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for monoclinic symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . From Figure 3, it is seen that the neighborhood contains only one plane of symmetry - the x − y planewith the distribution of micromoduli reflected exactly across this plane. There are several bonds withmicromoduli equal to zero, which points to the wide range of values in the material sti ff ness tensor, e.g. C = GPa vs. C = GPa . The wide range of axial sti ff ness terms in orthogonal directions resultsin a lumped distribution of micromoduli. Bonds with the highest micromoduli are aligned preferentiallywith the z -axis followed by the x -axis, not surprising as C and C have the highest values in the materialsti ff ness tensor. Micromoduli for many bonds not oriented perfectly in the x − z plane are zero and most ofthe sti ff ness in the y direction is carried by a few bonds along the y -axis. While this is not ideal in a realisticmodeling scenario, one workaround could be to set a lower bound on the micromoduli in the second stepof the calibration procedure, i.e. specify c + v ≥ l in Eq. (42) where l is a vector of the minimum valueof micromoduli desired. However, this vector l will have to chosen arbitrarily by trial and error since anydesired value might not produce a solution that satisfies all the constraints in the calibration procedure. Evenif a solution is found, it may have a higher relative error than the original solution found.˜ C PD =
135 18 .
33 62 0 0 4218 .
33 13 14 .
33 0 0 662 14 .
33 269 0 0 22 .
670 0 0 14 .
33 22 .
67 00 0 0 22 .
67 62 042 6 22 .
67 0 0 18 . . (57)The corresponding peridynamic sti ff ness tensor given by Eq. (57) is calibrated to within 4.988% relative13rror. Note that once again that the 5 constants C , C , C , C and C that are not a ff ected by Cauchy’srelations remain unchanged while the others (denoted by symbols in Eq. (55)) which are modified. Orthotropic symmetry can be generated by specifying two additional symmetry planes such that thesymmetry group is given by, S = { I , − I , R ( e ) , R ( e ) , R ( e ) } . (58)Orthotropic symmetry is commonly encountered when modeling fiber reinforced composite materials,perhaps the class of anisotropy that is most commonly found in peridynamics literature. Orthotropic sym-metry specifies that the material sti ff ness tensor can have 9 independent constants given by, C , C , C , C , C , C , C , C , C , (59)with all others being zero. Cauchy’s relations reduce the number of independent constants from 9 to 6 (c.f.Azdoud et al. [27]), C , C , C , C = C , C = C , C = C , (60)which can be visualized as follows, C C C C C C C C C C C
00 0 0 0 0 C → C (cid:12) ♦ (cid:12) C (cid:3) ♦ (cid:3) C (cid:3) ♦
00 0 0 0 0 (cid:12) . (61)The example Voigt elastic sti ff ness matrix, for Te W, is given by,˜ C re f =
143 1 37 0 0 01 3 3 0 0 037 3 102 0 0 00 0 0 2 0 00 0 0 0 46 00 0 0 0 0 1 . (62)The distribution of calibrated bond micromoduli is given in Figure 4 with the corresponding peridy-namic sti ff ness tensor in Eq. (63). Figure 4b quite clearly shows the orthotropic nature of the neighborhoodalthough it is limited to a plane near the center. The bonds preferentially aligned with the x –axis have thehighest micromoduli followed by the z -axis.Similar to the monoclinic material considered, the elastic constants span a wide range of values, e.g.from 3 in C to 143 in C giving a high universal anisotropy constant of 54.0623. This leads to manybonds having a zero sti ff ness similar to the monoclinic case, especially bonds that are oriented towards the y -axis. For practical problems the strategy of specifying a lower bound outlined in section 4.1.2 might beuseful however there is no guarantee that a solution would be found or that the error would be the lowestpossible for any desired lower bound. 14 a) (b) (c) Figure 4: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for orthotropic symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . The calibrated peridynamics sti ff ness tensor is given in Eq. (63) for which the relative error is foundto be 5.0958%. Note that C = C a priori satisfy Cauchy’s relations exactly which help reduce the errorin this case. As with previous cases, the axial sti ff ness terms remain unchanged where as the 6 other termsnamely the o ff -diagonal axial terms as well as the shear sti ff nesses are modified during the calibration.˜ C PD =
143 1 43 0 0 01 3 2 .
33 0 0 043 2 .
33 102 0 0 00 0 0 2 .
33 0 00 0 0 0 43 00 0 0 0 0 1 . (63) Trigonal symmetry may be generated by choosing three planes of symmetry as well, similar to or-thotropic symmetry, with the normals to the planes lying in the x − z plane such that the symmetry groupcan be written as, S = (cid:40) I , − I , R ( e ) , R √ e ± e (cid:41) . (64)The restrictions on the non-zero elastic constants are, C = C , C = C , C , C , C = C , C = C = − C , C = ( C − C ) / . (65)Note that along with other restrictions, Cauchy’s relations notably impose C = C which combinedwith the last of Eq. (65) imposes, ( C − C ) / = C = ⇒ C = C . (66)Interestingly Cauchy’s symmetry impose C = C which trigonal symmetry by itself possesses.Therefore, combining Eq. (65) and Eq. (66) the number of independent constants reduces from 6 to 4,and can be written as follows, C = C = C = C , C , C = C = C = C , C = C = − C , (67)15nd can be visualized as follows , C C C C C C C C C C C C C C C C C C → C (cid:12) ♦ ⊗ (cid:12) C (cid:3) ⊗ ♦ (cid:3) C ⊗ ⊗ (cid:3) ♦ ⊗ ⊗ (cid:12) . (68)Note that the ⊗ symbol in the C position denotes that the entry is related but not equal (sign is oppositethat of C and C ). The Voigt elastic sti ff ness matrix, with a universal anisotropy index of A U = . C) is given by, ˜ C re f =
464 159 141 −
45 0 0159 464 141 45 0 0141 141 493 0 0 0 −
45 45 0 125 0 00 0 0 0 125 −
450 0 0 0 −
45 153 . (69) (a) (b) (c) Figure 5: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for trigonal symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . The resulting distribution of bonds is given in Figure 5. Note that from visual inspection of Figure 5 theperidynamic neighborhood is not symmetric about the y − z plane or the x − z plane but is still symmetricabout the x − y plane. The calibrated peridynamic sti ff ness tensor is given in Eq. (70). As in previouscases, the axial diagonal terms are captured exactly in addition to the axial-shear coupling terms C , C and shear coupling term C . The calibrated values of C and C of 155 are very close to their referencevalues as well, with the largest di ff erence stemming from the axial-coupling terms C , C and the shearterms C , C . Although there are 5 unique elastic constants appearing in e ff ective peridynamic sti ff ness It is implied that the matrices given in Eq. (55) are symmetric, therefore C JI = C IJ and the elastic constants C IJ that arenot denoted by symbols satisfy the constraints of the particular material class chosen. For example, it is implied in Eq. (67) that C = C and they are not truly independent of each other. C = C are related to C = C as noted before ( C ≈ C ). The overall relative errorin the calibration is 2.632%.˜ C PD =
464 155 130 . −
45 0 0155 464 130 .
33 45 0 0130 .
33 130 .
33 493 0 0 0 −
45 45 0 130 .
33 0 00 0 0 0 130 . −
450 0 0 0 −
45 155 . (70) Tetragonal symmetry can be imagined to be orthotropic symmetry with two additional planes of reflec-tion, with a symmetry group given by, S = (cid:40) I , − I , R ( e ) , R ( e ) , R ( e ) , R (cid:32) √ e ± √ e (cid:33) (cid:41) . (71)Tetragonal symmetry specifies that the elastic sti ff ness tensor can have 6 independent constants given by, C = C , C , C = C , C , C , C = C , (72)with all others being zero. Cauchy’s relations reduce the number of independent constants from 6 to 4 givenby, C = C , C , C = C = C = C , C = C , (73)and can be visualized as follows, C C C C C C C C C C C
00 0 0 0 0 C → C (cid:12) (cid:3) (cid:12) C (cid:3) (cid:3) (cid:3) C (cid:3) (cid:3)
00 0 0 0 0 (cid:12) . (74)The example Voigt elastic sti ff ness matrix, for Si, is given by,˜ C re f =
212 70 58 0 0 070 212 58 0 0 058 58 179 0 0 00 0 0 58 0 00 0 0 0 58 00 0 0 0 0 85 . (75)Interestingly, this material partially satisfies Cauchy’s relations as seen by the equality of C , C , C , C which will help reduce the error in the least-squares calibration process. The universal anisotropy index forthis material is given by A U = . a) (b) (c) Figure 6: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for tetragonal symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . The solution is shown graphically in Figure 6 and the e ff ective peridynamics sti ff ness tensor is given by,˜ C PD =
212 80 58 0 0 080 212 58 0 0 058 58 179 0 0 00 0 0 58 0 00 0 0 0 58 00 0 0 0 0 80 . (76)The structure of the sti ff ness matrix is quite similar to that of orthotropic symmetry in Eq. (63) whichis not surprising as the only di ff erence is the equality of C , C , C , C . Therefore it is evident that withCauchy’s relations imposed, orthotropic symmetry has 6 independent constants whereas tetragonal has 5.The error in the calibration for this particular material is found to be 3 . Transverse isotropy is a special case of orthotropy wherein there exists a plane of isotropy resulting in 5independent constants. It can also be thought of as being similar to tetragonal symmetry with an additionalconstraint on the C which consequently reduces the number of independent constants from 6 to 5. Alongwith orthotropy, transverse isotropy has also been modeled in peridynamics, usually for composite laminatesin 2D.The symmetry group for a transversely isotropic material can be written as, S = (cid:40) I , − I , R ( e ) , R (cid:32) √ a + b ( a e + b e ) (cid:33) : a , b ∈ R (cid:41) . (77)For example, if the x − y plane is the plane of isotropy, then the 5 independent constants can be written as, C = C , C , C = C , C C = C , C = ( C − C ) / . (78)Similar to the trigonal symmetry, the last of Eq. (78) along with C = C in addition to other restrictionsimposed by Cauchy’s relations reduce the number of independent constants from 5 to 3 given by, C = C = C = C , C , C = C = C = C , (79)18nd can be visualized as follows, C C C C C C C C C C C
00 0 0 0 0 C → C (cid:12) (cid:3) (cid:12) C (cid:3) (cid:3) (cid:3) C (cid:3) (cid:3)
00 0 0 0 0 (cid:12) . (80)Note that Cauchy’s relations render the structure of the sti ff ness matrix in Eq. (80) exactly the same asthat for a tetragonal material in Eq. (74). The only di ff erence being that for the tetragonal case, C can bedefined independent of C . In either case, Cauchy’s relations enforce C = C .The reference Voigt elastic sti ff ness tensor for transverse isotropy (MoN for this example) is given by,˜ C re f =
499 177 235 0 0 0177 499 235 0 0 0235 235 714 0 0 00 0 0 241 0 00 0 0 0 241 00 0 0 0 0 161 . (81)The resulting bond micromoduli from the least-squares method is shown in Figure 7, (a) (b) (c) Figure 7: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for transversely isotropic symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . and the resulting e ff ective peridynamics sti ff ness tensor is found to be,˜ C PD =
499 166 .
33 239 0 0 0166 .
33 499 239 0 0 0239 235 714 0 0 00 0 0 239 0 00 0 0 0 239 00 0 0 0 0 166 . . (82)A simple visual examination of Figure 7 reveals that the bond micromoduli as a function of spatialcoordinates resembles a transversely isotropic material. Higher micromoduli in the bonds with directions19ut of the x − y plane impart the higher sti ff ness in the C sti ff ness term in Eq. (82). It is also interesting tonote the di ff erence in the distribution of bond micromoduli in plane with normal n = − e , which is clearly afunction of the bond directions, and n = e in which there is a very low variation in bond micromoduli as afunction of the bond angle with respect to the x axis, as shown in Figures 7b and 7c. It’s easy to verify thatthe calibrated sti ff ness matrix of course still satisfies Cauchy’s relations. The relative error in the calibrationis found to be 1.5335%. Cubic symmetry can be generated using 9 planes of reflection symmetry, the symmetry group for whichcan be written as, S = (cid:40) I , − I , R ( e ) , R ( e ) , R ( e ) , R (cid:32) √ e ± √ e (cid:33) , R (cid:32) √ e ± √ e (cid:33) , R (cid:32) √ e ± √ e (cid:33) (cid:41) . (83)In the case of cubic symmetry, the number of independent constants are 3 given by, C = C = C , C = C = C , C = C = C . (84)Cauchy’s relations reduce the number of independent constants from 3 to 2 given by, C = C = C , C = C = C = C = C = C . (85)and can be visualized as follows, C C C C C C C C C C C
00 0 0 0 0 C → C (cid:3) (cid:3) (cid:3) C (cid:3) (cid:3) (cid:3) C (cid:3) (cid:3)
00 0 0 0 0 (cid:3) . (86)An example of a material with Cubic symmetry is MgAl O (Spinel) for which the Voigt elastic sti ff nessmatrix is given by, ˜ C re f =
252 145 145 0 0 0145 252 145 0 0 0145 145 252 0 0 00 0 0 142 0 00 0 0 0 142 00 0 0 0 0 142 . (87)As seen in Eq. (87), the material very nearly satisfies Cauchy’s relations as the o ff -diagonal axial terms C , C , C are numerically very similar to the shear terms C , C , C . The anisotropy index A U for thismaterial is evaluated to be 1.2372.The resulting bond micromoduli from the least-squares method is shown in Figure 8,20 a) (b) (c) Figure 8: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for cubic symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . and the resulting e ff ective peridynamics sti ff ness tensor is found to be,˜ C PD =
252 143 . . . . . . . . . . (88)The relative error in the solution for this particular material is found to be 0.8028%. Finally, the least anisotropic symmetry - isotropy may be generated by choosing any orthogonal trans-formation or alternatively any reflection symmetry plane. In other words, the symmetry group for isotropymay be defined as an infinite number of orthogonal transformations. In general, for isotropic materials, thefollowing relationships in terms of the Voigt notation will hold, C = C = C , C = C = C , C = C = C = ( C − C ) / , (89)Isotropic symmetry reduces the number of independent elastic constants to just 2. In addition, Cauchy’srelations also impose the following, C = C = C = C = C = C , (90)thereby reducing the number of elastic constants to only 1 given by, C = C = C = C = C = C = C = C = C . (91)The structure of the Voigt sti ff ness matrix can be visualized as follows, C C C C C C C C C C C
00 0 0 0 0 C → C (cid:3) (cid:3) (cid:3) C (cid:3) (cid:3) (cid:3) C (cid:3) (cid:3)
00 0 0 0 0 (cid:3) . (92)21t is easy to see from Eq. (91) that these conditions reduce to the well known limitation on the Poisson’sratio ( ν = / C = C = ⇒ E (1 − ν )(1 + ν )(1 − ν ) = E ν (1 + ν )(1 − ν ) = ⇒ ν = . (93)The isotropic material that is chosen as an example is Pyroceram 9608 which happens to satisfy exactlyCauchy’s relations [28]. The anisotropy index A U is of course, evaluated to be zero as it should be forisotropic materials. The sti ff ness matrix in Voigt notation is given as,˜ C re f = . . . . . . . . . . . . . (94)Unlike evaluating a single value for the micromodulus using the analytical expression given in Eq. (18),the least-squares approach gives multiple values of bond micromoduli, even in the isotropic case. For thepurposes of comparison, the micromodulus evaluated using Eq. (18) is 0.2535 N / m (with an influencefunction of 1) whereas the discrete calibration is a distribution bounded by 0.1726 N / m and 0.2076 N / m (note that the bond force has units of force per unit volume squared). (a) (b) (c) Figure 9: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for isotropic symmetry with horizon δ = ω ( | ξ | ) = δ | ξ | . The bond micromoduli thus obtained are shown graphically in Figure 9. An interesting result is that thecollinear bonds have the same micromoduli as seen for example in the partial view in Figure 9b along the y and the z -axes. The bond micromoduli do depend on bond orientation with respect to material axes howeverultimately results in an e ff ective isotropic material. The e ff ective peridynamic sti ff ness tensor obtained aftercalibration using the least-squares method in Eq. (35) is found to be exactly equal to the reference sti ff nesstensor. The relative error in the solution is e ff ectively zero ( e = . × − ).22 a) (b) (c) Figure 10: (a) Complete neighborhood, and partial views of the neighborhood with normal to the cut plane at (b) n = − e and (c) n = e for isotropic symmetry with horizon δ = ω ( | ξ | ) =
1. An additional lower bound on the solution of0.23 N / m is applied. One may find a non-constant micromoduli distribution unusual and counter-intuitive, understandably,as the most commonly adopted approach is to use a constant micromodulus. However, it can be shown thatan approximately constant micromodulus distribution can be recovered using the proposed method as well.This can be done by using a well informed choice of bounds in the quadratic program defined in Eqs. (41),(42) and (43). The value of the micromodulus using an influence function, ω ( | ξ | ) = / m . Therefore, choosing a lower bound of 0.23 N / m (after some numerical experimentation) results ina solution that is given in Figure 10. It is quite evident that the micromoduli are nearly constant within thehorizon except for some bonds near the edge of the spherical horizon. The mean value of this distributionis found to be 0.2435 N / m which agrees quite well with the analytical value of 0.2535 N / m . Note that, aperfect match is not expected since (by a rough estimation) the volume of the full spherical neighborhoodis approximately 905 m , compared to a value of 925 m for the discrete neighborhood. Accounting for theslight di ff erence in volume, the analytical value can be scaled down, according to the di ff erence in volume,to arrive at a true analytical estimate of 0.2480 N / m which is in excellent agreement with the mean valueof 0.2435 N / m for the discrete distribution.It is important to note that the calibrated peridynamic sti ff ness tensor always satisfies Cauchy’s relationsregardless of the material symmetry. This implies that if the reference sti ff ness matrix that is input to theproposed least-squares method satisfies Cauchy’s relations a priori, then the relative error of calibration ise ff ectively zero (close to machine precision) as seen for the isotropic case. This can be done by modifyingthe material’s sti ff ness matrix appropriately using an analytical approach such as that given in [28]. The factthat the relative error is zero also indicates that volume correction is accounted for automatically, which isa common problem encountered when using an analytical micromodulus. This method allows one to avoidadditional computation of volume correction factors for bonds.23 ymmetry Constants(Traditional) Constants(Cauchy’s rel. imposed) Example A U Error (%)(Original) Error (%)(Cauchy’s rel. imposed)
Triclinic 21 15 KIO W 54.0623 5.0958 9.2958e-14Trigonal 6 4 Ta C 0.5875 2.632 8.1990e-14Tetragonal 6 4 Si 0.1231 3.8634 1.1168e-11Trans. Isotropic 5 3 MoN 0.2551 1.5335 6.9451e-11Cubic 3 2 MgAl O (Spinel) 1.2372 0.8028 1.5723e-13Isotropic 2 1 Pyroceram 9608 0 2.7225e-11 2.7225e-11 Table 1: Summary of results for all eight symmetries tested.
Table 1 and Figure 11 summarize results obtained from all eight symmetries. The maximum error inthe calibration procedure was found to occur for the orthotropic symmetry, where as the lowest error wasfound in the case of isotropy. The last column also shows the relative error when the reference materialsti ff ness matrix is modified such that Cauchy’s relations are satisfied. These modified matrices are takenfrom [28], except for the triclinic case where the model is run with the matrix given in Eq. (51) as the input.Not surprisingly, it is found that the errors in calibration are e ff ectively zero.Most importantly, no evidence of correlation between the universal anisotropy constant and the error isfound, or in other words, the calibration error was found to depend on how closely the material satisfiesCauchy’s relations but not on the material symmetry itself . Triclinic Monoclinic OrthotropicTetragonal Transversely isotropic Cubic IsotropicTrigonal
Figure 11: All eight symmetries using a horizon ratio of 6 and an inverse influence function.
Since the method described in this work is a numerical method, various parametric studies are conductedto demonstrate the generality and robustness of the method. While combinations of all parameters areendless, a few of practical significance are investigated, namely the choice of horizon ratio, horizon shape,influence function and rotation of the lattice or the neighborhood. Numerical experiments with lattice24otations also serve as verification that the proposed method presents solutions that are meaningful andsatisfy conditions of the material’s symmetry group. ff ect of Horizon Ratio Naturally, one of questions that is of practical importance is how the horizon ratio a ff ects the calibrationprocedure, i.e. does the number of bonds in the neighborhood a ff ect the quality of the solution? m = 7 T r i c li n i c M ono c li n i c O r t ho t r op i c T r i gon a l T e t r a gon a l C ub i c I s o t r op i c T r a n s v e r s e l y i s o t r op i c m = 3 m = 5 m = 7 m = 3 m = 5 Figure 12: All eight symmetries using example materials from section 4.1 with horizon ratio m of 3, 5 and 7. In general it appears, from numerical experimentation, that the calibration procedure is largely unaf-fected by the size of the neighborhood, provided the horizon is not impractically small. Figure 12 showscalibrated neighborhoods assuming horizon ratios of 3, 5 and 7 for all eight material symmetries using thesame example materials as in section 4.1. Naturally the magnitude of a typical bond micromodulus for aparticular material changes with horizon size, decreasing with an increase in horizon in general due to alarger number of bonds, but the overall pattern simply from visual inspection is the same. Moreover, itis observed that the error in the calibration is una ff ected for most reasonable and practical values of thehorizon ratio (including all the cases shown in Figure 12).For a horizon as small as 2 times the lattice spacing, a convergent solution is found for most cases,some with an increase in the relative error. For example, the error remains exactly the same when using ahorizon of 2 for isotropic, cubic, transversely isotropic, tetragonal and orthotropic symmetries but the errorincreases to 6.732% for trigonal and 3.9036% for triclinic symmetry where as a converged solution is notfound for the monoclinic case. Notice that the anisotropy index for the monoclinic material in these cases isnot as high as the orthotropic material, however a solution is still found for the orthotropic material. Sincethe monoclinic material has only a single plane of symmetry, it appears that there simply are not enoughbonds that can contribute to the desired anisotropic elasticity tensor.25or all practical purposes since generally a horizon of 3 or more is used, the least-squares method isfound to be independent of horizon size. ff ect of Neighborhood Shape While a circular neighborhood in 2D and spherical neighborhood in 3D are the most common becauseof the ease of finding analytical expressions for micromoduli in polar or spherical coordinates respectively,the choice of neighborhood can in fact be arbitrary. For example, Ahadi and Krockmal recently used anellipsoidal influence function with a elastic-plastic ordinary state based model [55]. In this parametric study,three di ff erent types of neighborhoods are considered namely cubic, cuboid and ellisoidal in addition to thespherical neighborhood considered in the results section. The cubic neighborhood is assumed to have equalsides of 2 δ , the cuboid neighborhood has sides of δ , δ , δ in the x , y , z directions respectively and theellipsoidal neighborhood has semi-axes of δ , δ , δ in the x , y , z directions respectively, where δ =
6. Thechoice of edge lengths and semi-axes lengths for cuboid and ellipsoidal shapes are arbitrary and the sameprocedure could be applied to any reasonably sized cuboids or ellipsoids. T r i c li n i c M ono c li n i c O r t ho t r op i c T r i gon a l I s o t r op i c T e t r a gon a l C ub i c T r a n s v e r s e l y i s o t r op i c Figure 13: All eight symmetries using example materials from section 4.1 with a cubic, cuboid and ellipsoid neighborhoods.
Firstly, for all cases tested, the distribution of micromoduli could be successfully found using the least-squares method with no change in the errors from those given in Table 1. There appear to be some interestingconsequences for the micromoduli based on the type of neighborhood chosen. In most cases, the distributionof micromoduli for the cubic neighborhood is similar to it’s spherical counterpart. However, in some caseslike tetragonal symmetry, the distribution is notably di ff erent. While the bonds with the highest micromoduliare found to lie in the in the x − y plane at angles π/ , π/ , π/ π/ , π/ , π/ π/ x − y plane in the cubic neighborhood areof length 6 √ δ , they are only of length 6 δ in the spherical case. Quite naturally, there are more bonds inthe x − y plane for the cubic case that can contribute to the sti ff ness matrix. Note also that the magnitudeof micromoduli for the cubic neighborhood are slightly smaller than the spherical neighborhood. The sameprinciple applies for the cuboid and ellipsoid neighborhoods as well - since the neighborhood contains morebonds in the z direction, the magnitude of micromoduli for those bonds is lower than those oriented inthe other two directions. Nevertheless, all four neighborhoods lead to same e ff ective peridynamic sti ff nessmatrix. Similar behavior can be seen clearly for the transversely isotropic material as well, but it is not soobvious for some others like the trigonal or cubic symmetries. In e ff ect, it is demonstrated that the proposedmethod is robust enough to handle a wide range of discretized neighborhoods. ff ect of Influence Function Another issue of practical significance is the choice of influence function ω ( | ξ | ). The influence functioncan be chosen to weight the contribution of bonds based on the bond lengths. As Seleson and Parks note,the choice of influence function can permit a rich spectrum of dynamic behavior in non-local models [56].The examples in section 4.1 used an inverse weight function of δ/ | ξ | . In addition to the inverse weightfunction, three other choices are considered - a constant weight ω ( | ξ | ) =
1, a hat function ω ( | ξ | ) = − | ξ | δ and a power law type function ω ( | ξ | ) = − (cid:16) | ξ | δ (cid:17) − , for each of the example materials. While the inverseinfluence function is singular at the center for zero bond lengths, decreasing to a minimum of unity at thehorizon, the hat and the power law functions start at unity at the center and decrease to zero at the horizon.While the choice of influence function combined with a choice of the neighborhood shape can yield endlesscombinations, for this parametric study a spherical horizon is chosen. The results for all cases are given inFigure 14, the neighborhoods are sectioned so as to be able to observe the bond distributions better.27 r i c li n i c M ono c li n i c O r t ho t r op i c I s o t r op i c T e t r a gon a l (| |) 1 = ξ | |(| |) 1 = − ξξ (| |) 1 = ξ | |(| |) 1 = − ξξ | |(| |) 1 − = + ξξ C ub i c T r i gon a l T r a n s v e r s e l y i s o t r op i c | |(| |) 1 − = + ξξ Figure 14: All eight symmetries using example materials from section 4.1 with influence functions ω ( | ξ | ) = ω ( | ξ | ) = − | ξ | δ and ω ( | ξ | ) = − (cid:16) | ξ | δ (cid:17) − In general, the bond micromoduli appear to be weighted quite intuitively in accordance with the natureof the influence function. This is more obvious with certain symmetries than others. For example in thecase of tetragonal symmetry, the highest micromoduli are found for bonds that are longer and close to thehorizon at angles of π/ , π/ , π/ π/ x − y plane for a constant influence function of 1. Witha hat function that decreases from 1 for | ξ | = | ξ | = δ , the bonds with the highest micromoduli shifttowards the center of the neighborhood, roughly for bond lengths of 2 √ ∆ . For the power law, this shiftsfurther inwards to √ ∆ . Similar behavior can be observed very clearly in the isotropic case as well wherethe bonds with the highest micromoduli shift from 6 ∆ to 2 ∆ and further to ∆ .In all, this parametric study demonstrates that the least-squares method of calibration if general enoughto successfully evaluate the distribution of micromoduli for a wide range of influence functions. ff ect of Lattice Rotation Lattice rotations can be of two types, arbitrary rotations of the lattice that may not fall within thesymmetry group of a material and specific rotations that do. In the first case, the peridynamic elasticitytensor evaluated in the direction of the rotated lattice, in general, may not be equal to that in the unrotatedlattice. However, both lattices should ideally result in the same elasticity tensor in any fixed coordinatesystem. In the second case where the rotation imposed on the lattice falls within the symmetry group ofthe material, the first condition should be met in addition to the fact that the elasticity tensors evaluated in The neighborhood is also interchangeably referred to as the ‘lattice’. 𝒆 𝒆 𝒆 𝝃 𝒆 𝒆 𝒆 𝑸 𝑸
Figure 15: Schematic of the original unrotated configuration of the neighborhood / lattice and the rotated configuration along withthe original unrotated coordinate basis { e , e , e } and the coordinate basis { e (cid:48) , e (cid:48) , e (cid:48) } aligned with the rotated lattice. Assume that a peridynamic neighborhood / lattice is calibrated to some elasticity tensor where both thelattice and the elasticity tensor are in an orthogonal coordinate system { e , e , e } and the lattice is alignedwith the coordinate system such that, C PDi jkl = (cid:90) H x ω ( | ξ | ) c ( ξ ) ξ i ξ j ξ k ξ l | ξ | dV x (cid:48) . (95)Let us assume that another lattice is calibrated to the same reference elasticity tensor but the lattice is in adi ff erent orientation not necessary aligned with { e , e , e } such that the peridynamic sti ff ness tensor is thenexpressed as, ˜ C PDi jkl = (cid:90) H x ω ( | ξ | ) ˜ c (cid:16) ˜ ξ (cid:17) ˜ ξ i ˜ ξ j ˜ ξ k ˜ ξ l | ξ | dV x (cid:48) . (96)The lattice in both Eq. (95) and Eq. (96) are assumed to be exactly the same, in terms of the horizon, numberof bonds etc., except for the orientation such that every bond satisfies the transformation ˜ ξ i = Q i j ξ j . Thenthe condition for the numerical results from the least-squares method to be invariant with lattice rotationswould be, C PDi jkl = ˜ C PDi jkl (97)Let us also assume that the proposed least-squares method is able to calibrate the bond micromoduli inboth cases accurately such that Eq. (97) is assumed to be satisfied. Then, C PDi jkl = (cid:90) H x ω ( | ξ | ) c ( ξ ) ξ i ξ j ξ k ξ l | ξ | dV x (cid:48) = ˜ C PDi jkl = (cid:90) H x ω ( | ξ | ) ˜ c (cid:16) ˜ ξ (cid:17) ˜ ξ i ˜ ξ j ˜ ξ k ˜ ξ l | ξ | dV x (cid:48) . (98)Note that C PDi jkl and ˜ C PDi jkl are defined in the same coordinate system { e , e , e } even though the distributionof micromoduli and bond vectors are di ff erent. The peridynamic elasticity tensor can also be defined in therotated coordinate system { e (cid:48) , e (cid:48) , e (cid:48) } aligned with the lattice in Eq. (96) such that e (cid:48) i = Q i j e j . Then theperidynamic sti ff ness tensor in the rotated coordinate system can be written as, C (cid:48) PDi jkl = (cid:90) H x ω ( | ξ | ) ˜ c (cid:16) ˜ ξ (cid:17) ξ (cid:48) i ξ (cid:48) j ξ (cid:48) k ξ (cid:48) l | ξ | dV x (cid:48) . (99) For the sake of simplicity, the notiation in the influence function is left unchanged as the influence function depends only onthe length of bond which is invariant with respect to the change of coordinate basis. c (cid:16) ˜ ξ (cid:17) of the rotated lattice remains the same, only the sti ff ness tensor isexpressed in the rotated coordinate system. The primed bond vectors ξ (cid:48) and ˜ ξ are essentially the samevectors in di ff erent coordinate systems such that the transformation ˜ ξ i = Q i j ξ (cid:48) j and the transformation˜ ξ i = Q i j ξ j can be used to write, C (cid:48) PDi jkl = (cid:90) H x ω ( | ξ | ) ˜ c (cid:16) ˜ ξ (cid:17) ξ i ξ j ξ k ξ l | ξ | dV x (cid:48) , (100)Equation Eq. (100) retains the micromoduli distribution ˜ c (cid:16) ˜ ξ (cid:17) from Eq. (99). Now, if the following relationfor the micromoduli distribution holds, ˜ c (cid:16) ˜ ξ (cid:17) = c ( ξ ) , (101)then Eq. (100) can be re-written as, C (cid:48) PDi jkl = C PDi jkl , (102)which implies that the transformation Q i j must lie in the symmetry group for the material represented by C PDi jkl . In other words, if the bond distribution of an unrotated and rotated lattice remains the same, then thattransformation must lie in the symmetry group of that material. Therefore, this simple test can also be usedto verify the numerical results. Although there are infinitely many combinations of symmetries and latticerotations that can be tested, we choose three cases to demonstrate. For all cases, a cubic lattice of semi edgelength of 3 is chosen resulting in a total of 342 bonds.The first symmetry chosen is triclinic, the most anisotropic with no planes of symmetry. The mostelementary of transformations are chosen - reflections about the e , e , e axes. N o r m a li z ed m i c r o m odu l u s c / c m a x Bond Number
Triclinic - reflection transformations - normalized micromoduli
Ref ( e )OriginalRef ( e ) Ref ( e ) Figure 16: Reflection transformations about the e , e , e axes (left) and the distribution of normalized micromoduli c / c max where c max is the maximum value of the micromodulus in the original configuration(right) for all three transformations along with theoriginal unrotated configuration for triclinic symmetry. Figure 16 shows the original configuration along with the three reflection transformations imposed. Inaddition, the figure also shows a plot of the normalized micromoduli versus the bond number for all four30ases (original unrotated configuration included). In all cases, the calibration error of 3.1873% obtained insection 4.1 is unchanged. This is not surprising as the transformations in this case keep the cubic latticeorthogonal to the coordinate system, therefore the distribution of the micromoduli with respect to spatialcoordinates is exactly the same in all cases as seen in Figure 16 (left). However, the micromoduli corre-sponding to specific bonds in all cases are di ff erent as seen in the plot in Figure 16 (right), meaning that thetransformations do not fall under the symmetry group of the material according to Eq. (101) and Eq. (102).Next, a series of rotations given by the rotation transformation R are imposed such that R = R x R y R z ,where R x , R y and R z are equi-angle rotations of θ about the e , e , e axes given by, R x = θ − sin θ θ cos θ , R y = cos θ − sin θ θ θ , R z = cos θ − sin θ θ cos θ
00 0 1 . (103) 𝜋6 𝜋3 𝜋22𝜋3 5𝜋6 𝜋 N o r m a li z ed m i c r o m odu l u s c / c m a x Bond Number
Triclinic - reflection transformations - normalized micromoduli 𝜋 𝜋 Figure 17: Rotation transformations of π , π , π , π , π , π about all three axes (left) and distribution of bond micromoduli versusbond number (right) for triclinic symmetry. For this example, six rotations of π , π , π , π , π , π are chosen. These rotations are chosen arbitrarilyand do not hold any special significance for this symmetry. Figure 17 shows the size cases of rotations(left) for which the overall error in the calibration is once again unchanged at 3.1873%. Instead of plottingthe normalized micromoduli for all cases, only those for 0 , π , π are shown for brevity and to demonstratethat for these cases, the distribution of micromoduli are equal. It is quite obvious that a rotation of π isequivalent to the original configuration. As it turns out, an equi-angle rotation of π results in R being equalto the identity transformation I which indeed is within the symmetry group of triclinic materials. For allother rotations, the distribution of micromoduli are di ff erent from the original configuration and from eachother. Therefore, for the triclinic symmetry out of all cases tested, both the conditions in Eq. (97) and Eq.(101) are satisfied only for the identity transformation, however the first condition is satisfied for all cases.31 % D i ff e r en c e a ft e r r o t a t i on Bond Number
Trigonal - rotation in the x-y plane - % Difference
120 240 𝜋3 2𝜋3𝜋32𝜋30 𝜋
Figure 18: Rotation transformations of π , π and π about all three axes (left) and % di ff erence in the bond micromoduli from theunrotated configuration versus bond number (right) for trigonal symmetry. The second symmetry chosen is the trigonal material which has a few symmetry planes other than theubiquitous identity and inversion transformations in it’s symmetry group. Figure 18 shows the unrotatedconfiguration along with three rotations of π , π and π about the e axis, i.e. in the x − y plane. In allcases, the e ff ective peridynamic sti ff ness tensor is identical to the one found in section 4.1 with the error ofcalibration constant at 2.632% thus satisfying the first condition Eq. (97). Unlike in the triclinic case wherethe equi-angle rotation of π/ π and π fall within the symmetrygroup of the trigonal material. Note that the lattice after rotations is not orthogonal to the coordinate system.The plot in Figure 18 (right) shows the % di ff erence in the micromoduli resulting from a rotation of π and π from the original configuration. It is seen that the di ff erence is not identically zero however is negligible.The maximum absolute % di ff erence is found to be < .
04% which explains the negligible impact on thecalibration error. Thus the second condition given in Eq. (101) can be said to satisfied as well. In otherrotation transformations tested which did not fall under the symmetry group of the material (not presentedhere for the sake of brevity), only the first criteria given in Eq. (97) was satisfied as expected.32 % D i ff e r en c e a ft e r r o t a t i on Bond Number
Isotropic - rotation in 3D - % Difference
60 120 240 300 𝜋3 2𝜋3𝜋6 𝜋 𝜋 Figure 19: Rotation transformations of π , π , π , π , π , π about all three axes (left) and % di ff erence in the bond micromoduli fromthe unrotated configuration versus bond number (right) for isotropic symmetry. Lastly, isotropic symmetry, the least anisotropic with infinitely many transformations in it’s symmetrygroup, is chosen to demonstrate the e ff ect of lattice rotations. Similar to the triclinic case, equi-anglerotations of of π , π , π , π , π , π are chosen as shown in Figure 19. As in the triclinic case, rotations of π and π result in a return to the original unrotated configuration and therefore quite obviously, the distributionof micromoduli is identical. Note however, that even with a visual inspection it is easy to see that themicromoduli for all cases of rotations are identical, no matter the rotation imposed. The plot shows the %di ff erence in the micromoduli from the original unrotated configuration which demonstrates that the changein micromoduli is basically zero (close to machine precision). In addition, since the material (Pyroceram9608 in this case) exactly satisfies Cauchy’s relations, the error in calibration is also found to be zero for allcases thus satisfying both conditions Eq. (97) and Eq. (101). While the previous subsections focus on the foundations and robustness of the proposed numericalscheme, it still remains to be used in an actual simulation. In this section, a uniaxial tensile test case issimulated to show that:1. A simulation can indeed be run by using calibrated micromoduli from the proposed scheme.2. The results match with the analytical solution from classical elasticity.3. Volume corrections are inherently taken care of using this model.4. There is a computational advantage of using a discretized model.33 𝑥𝑦𝑧 (a) (b) (c) Figure 20: (a) Schematic of the uniaxial test, (b) discretized model with 13824 particles in total and (c) a cut section of modelwith boundary particles shown in red and interior particles shown in blue.
The schematic of the proposed numerical test is given in Figure 20. A cube of size 24m × × ∆ . A relatively low horizon ratioof 2 is deliberately chosen (a value of ≥ ff ectsat the surface due to an incomplete horizon of material particles. This region of particles with incompletehorizons is shown in red in Figure 20 which forms a ‘boundary layer’ around an inner cube of dimensions20m × × average of the evaluatedmicromoduli. This simple scheme at the surface is adopted for the current purposes although more work inthis area in required (see e.g. [11]). The simplest possible weight function, ω ( | ξ | ) = . × − .A force on the top surface of the model (to the top layer of particles) equivalent to a traction of 1 MPa isapplied. The bottom surface, at z = u z =
0. For computational purposes, zerodisplacements are prescribed, u x = x = u y = y = ff ness matrix positive definite (a static solution scheme is adopted wherein a sti ff ness matrix forthe linearized problem is assembled and inverted to obtain the solution, see Appendix B).34 a) (b) Figure 21: A cut section of the model showing (a) axial sti ff ness terms C = C = C and (b) the o ff -diagonal / shear sti ff nessterms C = C = C = C = C = C using the proposed calibration scheme for micromoduli. First, the calibrated elastic constants of the model are illustrated in Figure 21. As expected, the calibratedelastic sti ff ness constants within the interior of the cube, away from the surfaces are more accurate, and theydegrade near the surfaces, especially near the edges and corners. However, the average values of the elasticconstants are still found to be quite good agreement with actual material constants. The average axialsti ff ness ( C = C = C ) is found to be 102.7010 GPa such that the relative error is 0.48%. The averageo ff -diagonal and shear sti ff ness terms ( C = C = C = C = C = C ) is found to be 33.7470GPa such that the error is 1.89%. In the interior, not surprisingly, excellent agreement is found betweenthe numerically calibrated and actual material constants. The calibrated value of axial sti ff ness is found tobe 103.2009 GPa and the o ff -diagonal and shear sti ff ness is found to be 34.4002 GPa, which have relativeerrors of less than 0.001%. (a) (b) (c) Figure 22: (a) x-displacement, (b) y-displacement and (c) displacement contours for the uniaxial tensile test using the proposedcalibration scheme for micromoduli.
The displacement fields obtained from the force-controlled uniaxial tensile test are shown in Figure22. The maximum z-displacement at the top surface is found to be 2 . × − m and the minimum xand y-displacements are found to be − . × − m . Therefore, the average strain ε zz calculate fromthe maximum displacement is evaluated to be 1 . × − . The strain ε zz evaluated only for the innerparticles away from the edge is also evaluated to be 1 . × − as shown in Figure. Strain for individualparticles is calculated by means of non-local shape tensors (see Appendix C) which are then averaged in35he interior. The absolute relative error in the strain evaluated by the peridynamic solution is 1.47% overthe whole domain including surfaces and only 0.68% in the interior. The error in the interior of the domainshould be ideally be even lower - however note that the boundary layer created due to surface e ff ects carriessome of the applied load thereby introducing errors in the interior strain field as well. Figure 23: Cut section of the model showing axial strain component ε zz calculated only for the interior particles away from thesurfaces. Note that the the calibrated values of the elastic constants and the resulting solution prove also thatvolume correction is inherently included in the proposed model. No additional correction schemes arerequired to ‘correct’ for the partial volumes of the material particles included within the non-local horizon.However, surface e ff ects which are inherent to peridynamics model still persist, an avenue of research thatis to be pursued in more detail in the context of the proposed models. Many surface correction schemes fortraditional bond-based peridynamics have been investigated, e.g. see Le and Bobaru [57]. (a) (b) (c) Figure 24: A cut section of the model showing (a) axial sti ff ness terms C = C = C and (b) the o ff -diagonal / shear sti ff nessterms C = C = C = C = C = C and (c) axial strain component ε zz calculated only for the interior particles away fromthe surfaces using the analytical value of micromoduli. Next, the exact same simulation is repeated with the analytical micromodulus from Eq. (18) keepingall other parameters the same. Figure 24 shows the axial and o ff -diagonal / shear sti ff nesses, and the ε zz strain field. The average values of the axial and o ff -diagonal / shear sti ff ness (in the interior away from theedges) are evaluated to be 106.4309 GPa and 30.3218 GPa, resulting in a relative absolute error of 3.13%36 rrors in the solution using analytical and numerical calibration C (int.) C (int.) C (incl. surf. C (incl. surf.) ε zz (int.) ε zz (incl. surf.) Numerical < < Analytical
Table 2: Relative errors in axial sti ff ness terms C = C = C , o ff -diagonal / shear sti ff ness terms C = C = C = C = C = C and axial strain component ε zz averaged over the entire domain including surfaces and only in the interior using the analyticalmicromodulus and the numerical micromoduli distribution obtained from the proposed scheme. and 11.86% respectively. When considering the average over the entire domain, these are evaluated tobe 97.8902 GPa and 27.2428 GPa resulting in relative absolute errors of 5.15% and 20.81%. These errorare significantly higher than those obtained using the proposed numerical method. The average strain ε zz calculated from the maximum displacement at the top surface is found to be 1 . × − leading to anabsolute relative error of 2.93%. The average strain evaluated in the interior is 1 . × − leading to anabsolute relative error of 5.57%. Summary of the errors evaluated are given in Table 2.It is important to note that in a realistic engineering simulation, a horizon ratio of 3 or more is generallychosen which will no doubt reduce the errors when using an analytical micromodulus. However, a low valueof 2 is chosen for the current purposes to illustrate the di ff erence when using these two approaches. Alsonote that there is significant error in the elastic constants even in the bulk, away from the surface when usingthe analytical micromodulus. This points to the errors arising from not accounting for particle volumes ofparticles - i.e. not using volume corrections. Again, these errors should decrease with the use of a higherhorizon ratio and volume correction. (a) (b) Figure 25: Sparsity pattern of the sti ff ness matrix assembled for the example problem with horizon ratios of (a) δ = ∆ and (b) δ = ∆ . There is computational advantage in using smaller horizons though, especially in 3D. For a sphericalneighborhood, the number of bonds for a horizon ratio of 2 is 32 and increases rather rapidly to 122 bonds fora horizon ratio of 3. Figure 25 the sparsity pattern of the sti ff ness matrix assembled for the example problemwith horizon ratios of (a) δ = ∆ and (b) δ = ∆ . The total number of non-zeros in the sti ff ness matrixincreases from approximately 1.73 million to 10.91 million, an increase of 6x. Therefore, the proposedapproach holds some promise in reducing computational e ff ort of bond-based peridynamics models as well.37 . Summary and Conclusions A novel, general and robust method to evaluate micromoduli for anisotropic bond-based models hasbeen presented. Most literature has focused on specific models for certain material symmetries but littlework has been done for incorporating anisotropy in bond-based peridynamics in a general way. The numer-ical method presented in this work can be used for any of the eight elastic material symmetry classes. Theproposed method is based on a least-squares system of equations where the bond micromoduli are treatedas the unknowns to be solved for, such that the error in the peridynamic sti ff ness tensor and the referencematerial sti ff ness tensor is minimized in a least-squares sense.The peridynamic sti ff ness tensor, due to the limitation of bond-based peridynamics, inherits some con-ditions on the sti ff ness tensor known as Cauchy’s ‘symmetry’ or Cauchy’s ‘relations’ along with the moretraditional major and minor symmetries. The implications of Cauchy’s relations for each material symmetryare presented and an example material is chosen for each symmetry to demonstrate the numerical method.It is shown that the resulting peridynamic sti ff ness tensor always satisfies Cauchy’s relations, without anya priori modification. Cauchy’s relations emerge naturally from the proposed numerical method of cali-bration. Therefore, if the reference material satisfies Cauchy’s relations a priori, the error in calibration ofthe bond micromoduli is e ff ectively zero. This also means that volume corrections are fully accounted for,given the fact that the numerical method is for a discrete set of bonds.Extensive parametric studies are conducted to demonstrate that the numerical method is general enoughto cater to a large variety of horizon sizes, influence functions, neighborhood shapes and lattice rotations. Itis found that if the lattice rotation falls under the symmetry group of the material, then not only is the peridy-namic sti ff ness tensor unchanged from the original configuration, it is also seen that the rotated peridynamicsti ff ness tensor is invariant as well, thus serving as model verification.Finally an example simulation is conducted, using a uniaxial tensile test setup to demonstrate the e ff ec-tiveness of this method. The e ff ective elastic sti ff ness constants away from surface is found to be within0.001% of the actual constants and the resulting strain is within 1% of the analytical solution. It is shownthat using an analytical micromodulus will result in a much higher error and has a computational disadvan-tage as well, keeping all other parameters unchanged.This work opens up many avenues for future development of both, fundamental aspects of anisotropicbond-based peridynamics as well as application to real-world problems.
6. Acknowledgements
This work has been funded in part by Corning Incorporated. The author thanks Jason T. Harris andRoss. J. Stewart for discussions, suggestions and support. The author would also like to acknowledge thecontinued support for peridynamics from Sam S. Zoubi.
References [1] Stewart A Silling. Reformulation of elasticity theory for discontinuities and long-range forces.
Journal of the Mechanics andPhysics of Solids , 48(1):175–209, 2000.[2] Stewart A Silling and Ebrahim Askari. A meshfree method based on the peridynamic model of solid mechanics.
Computers & Structures , 83(17):1526–1535, 2005.[3] Florin Bobaru and Monchai Duangpanya. The peridynamic formulation for transient heat conduction.
International Journalof Heat and Mass Transfer , 53(19):4047–4059, 2010.[4] Bahattin Kilic and Erdogan Madenci. Peridynamic theory for thermomechanical analysis.
Advanced Packaging, IEEETransactions on , 33(1):97–105, 2010.
5] Naveen Prakash and Gary D Seidel. Electromechanical peridynamics modeling of piezoresistive response of carbon nanotubenanocomposites.
Computational Materials Science , 113:154–170, 2016.[6] Walter Gerstle, Stewart Silling, David Read, Vinod Tewary, and Richard Lehoucq. Peridynamic simulation of electromigra-tion.
Comput Mater Continua , 8(2):75–92, 2008.[7] Amit Katiyar, John T Foster, Hisanao Ouchi, and Mukul M Sharma. A peridynamic formulation of pressure driven convectivefluid transport in porous media.
Journal of Computational Physics , 261:209–229, 2014.[8] Ziguang Chen and Florin Bobaru. Peridynamic modeling of pitting corrosion damage.
Journal of the Mechanics and Physicsof Solids , 78:352–381, 2015.[9] Erdogan Madenci and Selda Oterkus. Ordinary state-based peridynamics for thermoviscoelastic deformation.
EngineeringFracture Mechanics , 175:31–45, 2017.[10] QV Le, WK Chan, and Justin Schwartz. A two-dimensional ordinary, state-based peridynamic model for linearly elasticsolids.
International Journal for Numerical Methods in Engineering , 98(8):547–561, 2014.[11] John Mitchell, Stewart Silling, and David Littlewood. A position-aware linear solid constitutive model for peridynamics.
Journal of Mechanics of Materials and Structures , 10(5):539–557, 2015.[12] Michael Breitenfeld.
Quasi-static non-ordinary state-based peridynamics for the modeling of 3D fracture . PhD thesis,University of Illinois at Urbana-Champaign, 2014.[13] John Foster, Stewart A Silling, and Weinong Chen. An energy based failure criterion for use with peridynamic states.
International Journal for Multiscale Computational Engineering , 9(6), 2011.[14] J. T. Foster, S. A. Silling, and W. W. Chen. Viscoplasticity using peridynamics.
International Journal for Numerical Methodsin Engineering , 81(10):1242–1258, 2010.[15] Yunteng Wang, Xiaoping Zhou, and Xiao Xu. Numerical simulation of propagation and coalescence of flaws in rock mate-rials under compressive loads using the extended non-ordinary state-based peridynamics.
Engineering Fracture Mechanics ,163:248–273, 2016.[16] Stewart A Silling, M Epton, O Weckner, J Xu, and E Askari. Peridynamic states and constitutive modeling.
Journal ofElasticity , 88(2):151–184, 2007.[17] Thomas L Warren, Stewart A Silling, Abe Askari, Olaf Weckner, Michael A Epton, and Jifeng Xu. A non-ordinary state-based peridynamic method to model solid material deformation and fracture.
International Journal of Solids and Structures ,46(5):1186–1195, 2009.[18] Wenke Hu, Youn Doh Ha, and Florin Bobaru. Peridynamic model for dynamic fracture in unidirectional fiber-reinforcedcomposites.
Computer Methods in Applied Mechanics and Engineering , 217:247–261, 2012.[19] Wenke Hu, Youn Doh Ha, and Florin Bobaru. Modeling dynamic fracture and damage in a fiber-reinforced composite laminawith peridynamics.
International Journal for Multiscale Computational Engineering , 9(6), 2011.[20] Erkan Oterkus and Erdogan Madenci. Peridynamic analysis of fiber-reinforced composite materials.
Journal of Mechanicsof Materials and Structures , 7(1):45–84, 2012.[21] YL Hu and Erdogan Madenci. Bond-based peridynamic modeling of composite laminates with arbitrary fiber orientation andstacking sequence.
Composite structures , 153:139–175, 2016.[22] YL Hu, NV De Carvalho, and Erdogan Madenci. Peridynamic modeling of delamination growth in composite laminates.
Composite Structures , 132:610–620, 2015.[23] B Kilic, A Agwai, and Erdogan Madenci. Peridynamic theory for progressive damage prediction in center-cracked compositelaminates.
Composite Structures , 90(2):141–151, 2009.[24] Vito Diana and Siro Casolo. A full orthotropic micropolar peridynamic formulation for linearly elastic solids.
InternationalJournal of Mechanical Sciences , 160:140–155, 2019.[25] M Ghajari, L Iannucci, and P Curtis. A peridynamic material model for the analysis of dynamic crack propagation inorthotropic media.
Computer Methods in Applied Mechanics and Engineering , 276:431–452, 2014.[26] Wu Zhou, Dahsin Liu, and Ning Liu. Analyzing dynamic fracture process in fiber-reinforced composite materials with aperidynamic model.
Engineering Fracture Mechanics , 178:60–76, 2017.[27] Yan Azdoud, Fei Han, and Gilles Lubineau. A morphing framework to couple non-local and local anisotropic continua.
International Journal of Solids and Structures , 50(9):1332–1341, 2013.[28] Jeremy Trageser and Pablo Seleson. Anisotropic two-dimensional, plane strain, and plane stress models in classical linearelasticity and bond-based peridynamics. arXiv preprint arXiv:1905.12761 , 2019.[29] Ching S Chang and Anil Misra. Packing structure and mechanical properties of granulates.
Journal of Engineering Mechan-ics , 116(5):1077–1093, 1990.[30] Anil Misra and Payam Poorsolhjouy. Granular micromechanics model of anisotropic elasticity derived from gibbs potential.
Acta Mechanica , 227(5):1393–1413, 2016.[31] Huajian Gao and Patrick Klein. Numerical simulation of crack growth in an isotropic solid with randomized internal cohesivebonds.
Journal of the Mechanics and Physics of Solids , 46(2):187–218, 1998.
32] Wenyang Liu and Jung-Wuk Hong. Discretized peridynamics for brittle and ductile solids.
International Journal for Numer-ical Methods in Engineering , 89(8):1028–1046, 2012.[33] Wenyang Liu and Jung Wuk Hong. Discretized peridynamics for linear elastic solids.
Computational Mechanics , 50(5):579–590, 2012.[34] Siavash Nikravesh and Walter Gerstle. Improved state-based peridynamic lattice model including elasticity, plasticity anddamage.
Computer Modeling in Engineering & Sciences , 116(3):323–347, 2018.[35] Walter Herbert Gerstle.
Introduction to practical peridynamics: computational solid mechanics without stress and strain ,volume 1. World Scientific Publishing Company, 2015.[36] Gavin A Buxton, Christopher M Care, and Douglas J Cleaver. A lattice spring model of heterogeneous materials withplasticity.
Modelling and simulation in materials science and engineering , 9(6):485, 2001.[37] E Schlangen and JGM Van Mier. Simple lattice model for numerical simulation of fracture of concrete materials and struc-tures.
Materials and Structures , 25(9):534–542, 1992.[38] Erik Schlangen and Edward J Garboczi. Fracture simulations of concrete using lattice models: computational aspects.
Engineering fracture mechanics , 57(2-3):319–332, 1997.[39] Anand Jagota and Stephen J Bennison. Spring-network and finite-element models for elasticity and fracture. In
Non-linearityand breakdown in soft condensed matter , pages 186–201. Springer, 1994.[40] WA Curtin and H Scher. Brittle fracture in disordered materials: A spring network model.
Journal of Materials Research ,5(3):535–553, 1990.[41] Martin Ostoja-Starzewski. Lattice models in micromechanics.
Appl. Mech. Rev. , 55(1):35–60, 2002.[42] Hailong Chen, Yang Jiao, and Yongming Liu. A nonlocal lattice particle model for fracture simulation of anisotropic materi-als.
Composites Part B: Engineering , 90:141–151, 2016.[43] M Grah, K Alzebdeh, PY Sheng, MD Vaudin, KJ Bowman, and M Ostoja-Starzewski. Brittle intergranular failure in 2dmicrostructures: experiments and computer simulations.
Acta materialia , 44(10):4003–4018, 1996.[44] Naveen Prakash. Calibrating bond-based peridynamic parameters using a novel least squares approach.
Journal of Peridy-namics and Nonlocal Modeling , 1(1):45–55, 2019.[45] Naveen Prakash.
Coupled Electromechanical Peridynamics Modeling of Strain and Damage Sensing in Carbon NanotubeReinforced Polymer Nanocomposites . PhD thesis, Virginia Tech, 2017.[46] SA Silling and RB Lehoucq. Peridynamic theory of solid mechanics.
Advances in Applied mechanics , 44:73–168, 2010.[47] Augustus Edward Hough Love.
A treatise on the mathematical theory of elasticity . Cambridge university press, 2013.[48] John D Clayton.
Nonlinear mechanics of crystals , volume 177. Springer Science & Business Media, 2010.[49] Jeremy Trageser and Pablo Seleson. Bond-based peridynamics: A tale of two poisson’s ratios. 2019.[50] Thomas Chi-tsai Ting and Thomas Chi-tsai Ting.
Anisotropic elasticity: theory and applications . Number 45. OxfordUniversity Press on Demand, 1996.[51] Peter Chadwick, Maurizio Vianello, and Stephen C Cowin. A new proof that the number of linear elastic symmetries is eight.
Journal of the Mechanics and Physics of Solids , 49(11):2471–2492, 2001.[52] Stephen C Cowin and Morteza M Mehrabadi. On the identification of material symmetry for anisotropic elastic materials.
The Quarterly Journal of Mechanics and Applied Mathematics , 40(4):451–476, 1987.[53] Maarten de Jong, Wei Chen, Thomas Angsten, Anubhav Jain, Randy Notestine, Anthony Gamst, Marcel Sluiter, ChaitanyaKrishna Ande, Sybrand van der Zwaag, Jose J Plata, Cormac Toher, Stefano Curtarolo, Gerbrand Ceder, Kristin A. Persson,and Mark Asta. Charting the complete elastic properties of inorganic crystalline compounds.
Scientific Data , 2, 03 2015.[54] Shivakumar I Ranganathan and Martin Ostoja-Starzewski. Universal elastic anisotropy index.
Physical Review Letters ,101(5):055504, 2008.[55] Aylin Ahadi and Jakob Krochmal. Anisotropic peridynamic model—formulation and implementation.
AIMS MaterialsScience , 5(4):742, 2018.[56] Pablo Seleson and Michael Parks. On the role of the influence function in the peridynamic theory.
International Journal forMultiscale Computational Engineering , 9(6), 2011.[57] QV Le and F Bobaru. Surface corrections for peridynamic models in elasticity and fracture.
Computational Mechanics ,pages 1–20, 2017.[58] Naveen Prakash and Stewart Ross J. A multi-threaded method to assemble a sparse sti ff ness matrix for quasi-static solutionsof linearized bond-based peridynamics. Journal of Peridynamics and Nonlocal Modeling . ppendix A. Universal Anisotropy Index Assuming that C is the sti ff ness tensor written in Voigt notation, the universal anisotropy index can beevaluated as follows [53, 54], s = C − , (A.1)where s is the compliance matrix in Voigt notation. The Voigt and Reuss estimates for the bulk and shearmodulus are, 9 K V = ( C + C + C ) + C + C + C ) (A.2)1 / K R = ( s + s + s ) + s + s + s ) (A.3)15 G V = ( C + C + C ) − ( C + C + C ) + C + C + C ) (A.4)15 / G R = s + s + s ) − s + s + s ) + s + s + s ) (A.5)The average of the Voigt and Reuss estimates are,2 K VRH = ( K V + K R ) (A.6)2 G VRH = ( G V + G R ) (A.7)The universal anisotropy index is, A U = G V G R + K V K R − ≥ Appendix B. Static Solution Scheme
The linearized equations of bond-based peridynamics are solved by assembling a sti ff ness matrix andinverting the sti ff ness matrix to find the displacement field, i.e. a static solution scheme is employed,ignoring inertial e ff ects. Complete details on the assembly of a sparse sti ff ness matrix in 2D is given in [58]which has been extended to 3D in this work. A brief explanation of the assembly of the sti ff ness matrixresulting from the global system of equations is presented.In the current work, the equations are implemented in a relatively straightforward program written inMATLAB ® . A coordinate (COO) description of the contents of the sti ff ness matrix is used to formulate asparse representation of the sti ff ness matrix thereby reducing memory and computational e ff ort. Informationabout the neighbors and the bond sti ff nesses are used to directly assemble the global sti ff ness matrix.Each particle in 3D has three degrees of freedom, therefore the global sti ff ness matrix is of size 3N × ff ness matrix- one of for each degree of freedom. Since the neighbor list of every particle is known, the contributions ofbonds to the sti ff ness of every degree of freedom can be directly accounted for.Every bond connects a pair of particles. The internal force acting on a particle due to a bond is a vectorof 3 dimensions that is a function of 3 components of displacements for each particle (total of 6, 3 for eachconnected particle), therefore, the sti ff ness matrix of each bond is a 3 × ff ness terms. The first 9 sti ff ness terms (11, 12, 13, 21, 22, 23, 31, 32, 33) represents a coupling between41he forces acting on a particle due to it’s own displacement and the other 9 sti ff ness terms (14, 15, 16, 24,25, 26, 34, 35, 36) couple the forces to the displacement of the pairing particle. By extension, if a particlehas M neighbors, then there are a total of 9 +
9M sti ff ness terms.In many ways, this is akin to a 3D truss element in traditional finite element analysis. The local sti ff nessmatrix of such a truss element can be described by a 6 × However , it isimportant to note that in this particular approach, the global sti ff ness matrix is assembled one particle (3rows) at a time. Therefore, for a pair of particles, say IJ, the sti ff ness matrix terms related to I total 18 terms,not 36.Algorithm 1 details the method to assemble the sti ff ness matrix.42 lgorithm 1 Evaluate sti ff ness matrix n = (cid:46) Count of the non-zero terms for i ← do K11 ←
0, K22 ←
0, K33 ← (cid:46) Diagonal terms K12 ←
0, K13 ←
0, K14 ←
0, K15 ←
0, K16 ← (cid:46) O ff -diagonal terms K21 ←
0, K23 ←
0, K24 ←
0, K25 ←
0, K26 ← (cid:46) O ff -diagonal terms K31 ←
0, K32 ←
0, K34 ←
0, K35 ←
0, K36 ← (cid:46) O ff -diagonal terms for j ← Neighbors of i do K11 ← K11 + ( ω ( | ξ | ) c ξ x ξ x ∆ V j ) / | ξ | K22 ← K22 + ( ω ( | ξ | ) c ξ y ξ y ∆ V j ) / | ξ | K33 ← K33 + ( ω ( | ξ | ) c ξ z ξ z ∆ V j ) / | ξ | K12 ← K12 + ( ω ( | ξ | ) c ξ x ξ y ∆ V j ) / | ξ | K13 ← K13 + ( ω ( | ξ | ) c ξ x ξ z ∆ V j ) / | ξ | K23 ← K23 + ( ω ( | ξ | ) c ξ y ξ z ∆ V j ) / | ξ | K21 ← K21 + ( ω ( | ξ | ) c ξ y ξ x ∆ V j ) / | ξ | K31 ← K31 + ( ω ( | ξ | ) c ξ z ξ x ∆ V j ) / | ξ | K32 ← K32 + ( ω ( | ξ | ) c ξ z ξ y ∆ V j ) / | ξ | K14 ← -( ω ( | ξ | ) c ξ x ξ x ∆ V j ) / | ξ | K15 ← -( ω ( | ξ | ) c ξ x ξ y ∆ V j ) / | ξ | K16 ← -( ω ( | ξ | ) c ξ x ξ z ∆ V j ) / | ξ | K24 ← -( ω ( | ξ | ) c ξ y ξ x ∆ V j ) / | ξ | K25 ← -( ω ( | ξ | ) c ξ y ξ y ∆ V j ) / | ξ | K26 ← -( ω ( | ξ | ) c ξ y ξ z ∆ V j ) / | ξ | K34 ← -( ω ( | ξ | ) c ξ z ξ x ∆ V j ) / | ξ | K35 ← -( ω ( | ξ | ) c ξ z ξ y ∆ V j ) / | ξ | K36 ← -( ω ( | ξ | ) c ξ z ξ z ∆ V j ) / | ξ | row(n + ←
3i - 2, col(n + ←
3j - 2, val(n + ← K14 row(n + ←
3i - 2, col(n + ←
3j - 1, val(n + ← K15 row(n + ←
3i - 2, col(n + ←
3j - 0, val(n + ← K16 row(n + ←
3i - 1, col(n + ←
3j - 2, val(n + ← K24 row(n + ←
3i - 1, col(n + ←
3j - 1, val(n + ← K25 row(n + ←
3i - 1, col(n + ←
3j - 0, val(n + ← K26 row(n + ←
3i - 0, col(n + ←
3j - 2, val(n + ← K34 row(n + ←
3i - 0, col(n + ←
3j - 1, val(n + ← K35 row(n + ←
3i - 0, col(n + ←
3j - 0, val(n + ← K36 n ← n + end for row(n + ←
3i - 2, col(n + ←
3i - 2, val(n + ← K11 row(n + ←
3i - 1, col(n + ←
3i - 1, val(n + ← K22 row(n + ←
3i - 0, col(n + ←
3i - 0, val(n + ← K33 row(n + ←
3i - 2, col(n + ←
3i - 1, val(n + ← K12 row(n + ←
3i - 2, col(n + ←
3i - 0, val(n + ← K13 row(n + ←
3i - 1, col(n + ←
3i - 0, val(n + ← K23 row(n + ←
3i - 1, col(n + ←
3i - 2, val(n + ← K21 row(n + ←
3i - 0, col(n + ←
3i - 2, val(n + ← K31 row(n + ←
3i - 0, col(n + ←
3i - 1, val(n + ← K32 n ← n + end for c shown in the algorithm is not constant but is taken from micromodulidistribution that is obtained from the proposed numerical method. After the row , col and val vectors arepopulated, it is relatively straightforward to form a sparse sti ff ness matrix using the sparse function inMATLAB. Appendix C. Non-local Shape Tensors
Non-local shape tensors, commonly used in state-based peridynamics are used to calculate the strain ε [14, 16]. First, a non-local shape tensor is evaluated as defined by K [ x ] = (cid:90) H x ω ( | ξ | )( ξ ⊗ ξ ) dV x (cid:48) (C.1)where ξ is the undeformed bond vector and ω ( | ξ | ) is the influence function. Using the shape tensor, theapproximate non-local deformation gradient is evaluated by F [ x , t ] = (cid:34)(cid:90) H x ω ( | ξ | )( Y (cid:104) ξ (cid:105) ⊗ ξ ) dV x (cid:48) (cid:35) K − (C.2)where Y (cid:104) ξ (cid:105) is the deformed image of the bond and is equal to ξ + η . This non-local deformation gradi-ent is consistent with the classical local definition of the deformation gradient for homogeneous deforma-tions. Using basic kinematics, the Lagrangian strain can be obtained from the deformation gradient with E = (cid:16) F T F − I (cid:17) . Note that for small deformations, the Lagrangian strain tensor E will be equal to theinfinitesimal strain tensor ε . Appendix D. Isotropic Symmetry with Varying Micromoduli
While a non-constant micromodulus is understandably counter-intuitive to the idea of an ‘isotropic’bond-based peridynamic material, we defer to the well-established and broader notion of material symme-try of the elasticity tensor rather than the isotropy of the pairwise force function [1]. In other words, wedefer to idea that the e ff ective peridynamic sti ffff