The Boundary Element Method for Ordinary State-Based Peridynamics

Xue Liang and Linjuan Wang

1State Key Laboratory for Turbulence and Complex Systems,Department of Mechanics and Engineering Science,College of Engineering,Peking University,Beijing,100871,China

2CAPT-HEDPS,and IFSA Collaborative Innovation Center of MoE,College of Engineering,Peking University,Beijing,100871,China

3School of Astronautics,Beihang University,Beijing,100191,China

ABSTRACT The peridynamics (PD),as a promising nonlocal continuum mechanics theory,shines in solving discontinuous problems.Up to now,various numerical methods,such as the peridynamic mesh-free particle method(PD-MPM),peridynamic finite element method (PD-FEM),and peridynamic boundary element method (PD-BEM),have been proposed.PD-BEM,in particular,outperforms other methods by eliminating spurious boundary softening,efficiently handling infinite problems,and ensuring high computational accuracy.However,the existing PD-BEM is constructed exclusively for bond-based peridynamics(BBPD)with fixed Poisson’s ratio,limiting its applicability to crack propagation problems and scenarios involving infinite or semi-infinite problems.In this paper,we address these limitations by introducing the boundary element method (BEM) for ordinary state-based peridynamics(OSPD-BEM).Additionally,we present a crack propagation model embedded within the framework of OSPD-BEM to simulate crack propagations.To validate the effectiveness of OSPD-BEM,we conduct four numerical examples:deformation under uniaxial loading,crack initiation in a double-notched specimen,wedge-splitting test,and threepoint bending test.The results demonstrate the accuracy and efficiency of OSPD-BEM,highlighting its capability to successfully eliminate spurious boundary softening phenomena under varying Poisson’s ratios.Moreover,OSPDBEM significantly reduces computational time and exhibits greater consistency with experimental results compared to PD-MPM.

KEYWORDS Ordinary state-based peridynamics;boundary element method;crack propagation;fracture toughness

1 Introduction

Peridynamics,as a nonlocal continuum mechanics theory,has garnered increasing attention owing to its notable advantage in addressing discontinuous problems [1,2].This advantage arises from its unique approach of replacing spatial differential operators with integral operators in the equilibrium equation.Over time,peridynamic theory has seen continuous development and improvement,leading to the introduction of various theoretical models,e.g.,the dual horizon peridynamics[3,4],nonlocal operator methods[5],element-based peridynamic models[6,7],viscoelastic models[8],and elastoplastic theories [9].At the same time,a series of numerical approaches have emerged to leverage the advantages of Peridynamics (PD) in addressing various issues,including discontinuous problems[10,11],microscale problems[12,13]and multiscale problems[14–16].The one gaining the most attention is the Peridynamic Mesh-free Particle Method (PD-MPM),where the peridynamic equilibrium equation is discretized directly in terms of timing and spacing [17].Kilic et al.[18]explored a collocation point method drawing support from the Gaussian integral formula to optimize the nonlocal numerical integrals in the PD-MPM.Chen et al.[19] proposed the peridynamic finite element method (PD-FEM),drawing support from the energy principle.Tian et al.[20] adopted a central difference scheme for handling the PD equation.Liang et al.[21]put forward the bond-based peridynamics(BBPD)boundary element method(BBPD-BEM).

Several coupled numerical methods have been proposed based on the fundamental numerical approaches mentioned above.For example,to enhance computational efficiency,PD-MPM is coupled with numerical methods of classical continuum mechanics[22–26].These coupled methods essentially integrate two types of continuum mechanical media: the PD medium and the classical continuum medium.However,utilizing these coupled methods for investigating the responses of a single-phase material may result in a mismatch between numerical models and the actual research project.To address this issue,it is more advisable to couple different peridynamic numerical methods,such as PD-MPM with PD-FEM or PD-BEM.Therefore,the advancement of fundamental methods,including PD-MPM,PD-FEM,and PD-BEM,holds crucial importance for the practical applications of peridynamics.It is worth clarifying the term“PD-BEM”to avoid misunderstandings.PD-BEM is a numerical method constructed based on peridynamic theory,differing from numerical methods that combine the classical local continuum theory’s Boundary Element Method(BEM)with the mesh-free particle method based on peridynamic theory[25–28].

Our research is dedicated to PD-BEM,a focal point that offers distinct advantages.In comparison to alternative numerical methods,BEM stands out for its efficiency enhancement achieved through dimension reduction [21].Specifically,BBPD-BEM exhibits computational speeds two orders of magnitude faster than PD-MPM in computational domains without destruction.Furthermore,BBPD-BEM effectively circumvents spurious boundary-softening phenomena,facilitating PD calculations in infinite domains.However,it is essential to acknowledge that BBPD-BEM is built upon BBPD,which features only one independent material parameter for isotropic peridynamic materials.Additionally,the absence of a crack propagation model in BBPD-BEM restricts its ability to address crack propagation problems.To overcome these limitations,our paper aims to introduce BEM for ordinary state-based peridynamics(OSPD),one of the two typologies of state-based PD[2].Within the numerical framework of OSPD-BEM,we propose the crack propagation model inspired by the cohesive crack model[29–31]and the PD bilinear model[32,33].This approach addresses the identified shortcomings in BBPD-BEM,enhancing its versatility and applicability.

The process of outlining the BEM for OSPD,referred to as OSPD-BEM,unfolds through the following steps.Firstly,a boundary integral equation (BIE) for OSPD is derived,drawing support from Green’s function[34,35]and the nonlocal operator theory[36,37]in Section 2.Following this,a crack propagation model for the OSPD-BEM is proposed in Section 3.In Section 4,the accuracy and efficiency of OSPD-BEM are then demonstrated through the presentation of four numerical examples.The conclusions are given in Section 5.For ease of reading,the symbols in the paper are listed in Table 1 before the text begins.

Table 1:Nomenclature

2 The Boundary Integral Equation

Firstly,we briefly introduce the linear elastic OSPD[38,39]with the volume-constrained boundary[40,41].The equilibrium equation for the linear elastic OSPD is as follows[37]:

in whichQ(x,y)is a tensor function between pointxand pointy;P(y)is a vector function of pointy.β(x,y)is given as follows:

Γ(x,y)in Eqs.(1)and(4)is related to the nonlocal weighting functionϖ(x,y)in Eq.(1):

in whichψ(x)is dependent on the material parameters ϕ(x)andχ(x).In the three dimensional problem,material parameters ϕ(x)andχ(x)in Eq.(1)are expressed as

whereμdenotes shear modulus;κdenotes bulk modulus;andis

Figure 1:The diagram for volume-constrained boundary(color online)

For the two-dimensional problem,material parameters ϕ(x)andχ(x)in Eq.(1)are expressed as

It is noteworthy that the literature[5,43]also give the similar nonlocal operators.The connection of both nonlocal operators is in Appendix A.

Stemming from simplifying form,we introduce the following notations:

For two instances with different body forces and boundary constraints,Eq.(1)can be rewritten as

Following the derivations of the reciprocal theorem of the BBPD in [21],we can obtain the reciprocal theorem for OSPD

whose form is consistent with that of bond-based PD except for the operators defined in Eqs.(13)and(14).The integrals on the left side and right one respectively denote internal integrals and volumeconstrained boundary integrals.One can find that the boundary integrals are also volume integrals,which deprives the advantage of the dimensionality reduction for the BEM.Thus we introduce an extra condition to the volume-constrained boundary,which converts the integral on the right side of Eq.(15) from the one in the volume-constrained boundary to the one on the local boundary.The equivalency relationship to transform the integral in volume-constrained boundaryΩτinto the one on local boundary∂Ωis

whereuis authentic displacement state,it meets Eq.(1).usis a possible displacement state,it satisfies displacement constraint in Eq.(1):

whereOdenotes unit spherical surface;kdenotes unit normal vector of the unit spherical surfaceO;dΩmis differential solid angle on the surfaceOin directionk.f(x+yk,v(x+yk),x-zk,v(x-zk))denotes response function between pointx+ykand pointx-zk[34],which is given as follows:

whereKc(x+yk,x-zk)is the micromodulus tensor for the linear elastic OSPD[34].We know thatM (v)(x)on left side of Eq.(16)is analogous to body force density inΩτ,andon right side of Eq.(16)is analogous to surface traction on∂Ω.As a result,the equivalency relationship in Eq.(16)is perceived as a virtual work principle that is connected with possible deformed stateus.Applyingus=v1,u=v2andus=v2,u=v1to Eq.(16)respectively,one can obtain the following two equations:

Substituting Eqs.(20)and(21)into Eq.(15),the reciprocal theorem is rewritten as

which builds up a connection between internal integrals and classical boundary integrals,and is a precondition of BIE that will be derived as follows.

Imitating opinion in traditional BEM[46].We take into accountProblem1 corresponding to the actual problem in the reciprocal theorem,andProblem2 denotes infinite domain Green’s functionuG[35],Then,based on the Eq.(1),equilibrium equations corresponding toProblem1 andProblem2 are respectively described as

whereEidenotes a coordinate axis vector.meets

Applying first equations for Eqs.(23)and(24)to(26),we obtain

Figure 2:The diagram for the boundary limit process[21](color online)

The result for the integral on left side of Eq.(28)is-due to property of Dirac function.The result for the first item on the right side of Eq.(28)equals to Cauchy principal value(CPV)of integral on∂Ω.The result for the second item on the right side of Eq.(28)involves analyzing singularity for integrand.The singularity for integrand exists in Green’s functionuG.According to the literature[35],the logarithmic and Dirac singularity are involved.The integrable property of logarithmic singularity and the convolution property of the Dirac function cause the result for the second item on the right side of Eq.(28)to vanish.Then,Eq.(28)is given as below:

Adopting an approach in the literature [47,48] to build up a relation between PD boundary constraints and local boundary constraints,we obtain the following relationship:

which is the BIE of OSPD for the static case.For the dynamic case,we give BIE corresponding to the Laplace domain in Appendix B.Once a displacement and force flux on the local boundary∂Ωare obtained,one can calculate the displacement of any material point within the domainΩthrough the BIE Eqs.(33)and(B.13).Therefore,we just need to calculate the displacement and force on the local boundary through the discretization method in[21].The details will be repeated.So far,one can use the OSPD-BEM to simulate static and dynamic problems without fracture.For fracture problems,it is essential to research the crack propagation model for any BEM.Therefore,we propose a crack propagation model that is suitable for the present numerical method.

3 Crack Propagation Model

The crack propagation model for the OSPD-BEM proposed in this paper is inspired by the cohesive crack model [29–31] and the PD bilinear model [32,33].We take a continuum medium in Fig.3 as an example,and the construction of the model is given as follows.

Figure 3:The diagram of the crack propagation model(color online)

In Fig.3,S is the preset crack propagation path inΩ.Ωis divided into two subdivisionsΩ1andΩ2with a path S.e1and S1compose the boundary ofΩ1,whilee2and S2compose the boundary ofΩ2.We can convert the process of solvingΩto the process of solvingΩ1andΩ2.Therefore,the BIE of the OSPD Eq.(33)is applied to the boundarye1+S1and the boundarye2+S2,respectively.If there is no crack,S1and S2are the same surface,so the continuity condition should be satisfied on S1and S2.As shown in Fig.4,for any pair of associated points P1and P2on S1and S2,the continuity condition is expressed as

Figure 4:The diagram for the crack propagation model in the elastic stage(color online)

For the crack propagation model,we consider that the fracture at any material point P at the preset crack propagation path is divided into three stages: the elastic stage,the adhesive stage and the defunct stage.The elastic stage is given in Fig.4,and the calculation can be completed with the continuity condition Eqs.(34)and(35).The defunct stage is given in Fig.5,and the crack occurs at this time.At the defunct stage,the material point P is divided into P1and P2,and the force flux vectorsandvanish.The calculation is also easily completed at this stage.Next,we discuss the complicated adhesive stage in Fig.6.The adhesive stage needs to answer two problems.One is about the end of the elastic stage and the start of the defunct stage,and the other is about the constitutive relationship between the PD force flux vectors,i.e.,,and the displacement vectors,i.e.,u1andu2.We start with the problem two.In Fig.6,the PD force flux vectorssatisfy the equilibrium relationship

Figure 5:The diagram for the fracture calculation model in the defunct stage(color online)

Figure 6:The diagram of the crack propagation model in the adhesive stage(color online)

We assume that the force flux vectorand the relative displacementsatisfy the linear relationship shown in Fig.7.A more complex constitutive relationship,e.g.,those in[49,50],can be considered in the future.Next,let’s answer the problem one.In Fig.7,the force flux vectoris referred to as the elastic limit flux,and corresponds to the linear elastic stretch limit[32].According to the article[32,48,51],is expressed as

Figure 7:The diagram for the constitutive relationship in the adhesive stage(color online)

whereμdenotes shear modulus;κdenotes bulk modulus;hris horizon;G0is the elastic limit energy density [32].Whenin Eq.(34) achieves,the elastic stage is ended.uMin Fig.7 is referred to as the adhesive limit displacement,and it is determined by a Law of energy conservation for a process of the new crack growing[52].

whereGICis the critical energy release rate[51].Whenin Eq.(37)achievesuM,the defunct stage is started.The solving procedure of our crack propagation model can be summarized as Fig.8.

Figure 8:The flow chart for the solving process of the crack propagation model

4 Numerical Examples

We will take four numerical examples to confirm the accuracy and efficiency of OSPD-BEM in this section.Firstly,a two-dimensional square plate under uniaxial loading is investigated to estimate the numerical Poisson’s ratio as those in the literature[53,54].It displays that our numerical method can be more accurate than the PD-MPM.Secondly,we discuss the crack initiation in a double-notched specimen as the literature[21],which reveals the accuracy and efficiency of the OSPD-BEM compared with the PD-MPM.Thirdly,the wedge-splitting test is simulated to research a fracture toughness for the concrete specimen[55,56],and its result is compared to experimental results,which displays again the the accuracy and efficiency of OSPD-BEM.Finally,a three-point bending experiment is executed.

4.1 Two-Dimensional Square Plate under Uniaxial Loading

We consider a two-dimensional square plate withstanding the uniaxial tensile loading in the directionx,as displayed in Fig.9.The calculation parameters can be obtained in Table 2.adenotes length of the sides;hdenotes thickness;Edenotes elastic modulus;pdenotes tensile loading;hrdenotes horizon;Δis the grid spacing.For a given Poisson’s ratioν,we can obtain the deformation of the plate,and then calculate the numerical Poisson’s ratioswhich should be consistent with the given Poisson’s ratios.However,due to an error in the numerical method,the numerical Poisson’s ratio is probably different from the given Poisson’s ratio.The error between the two can be used to test the accuracy of one numerical method.Therefore,in this numerical example,we investigate the deformation state of the plate with different Poisson’s ratios,i.e.,ν=0.1,0.2,0.3.Then the numerical Poisson’s ratiosis calculated as

Table 2:Calculation parameters for the stretching square plate

whereΔydis the average absolute value of the displacement in the directionyfor all discrete material points.Δxdis the average absolute value of the displacement in the directionxfor all discrete material points.The numerical Poisson’s ratiosobtained by OSPD-BEM and PD-MPM are listed in Table 3.The errors between the numerical Poisson’s ratios and the given ones are also given in Table 3.

Table 3:Calculation results for the numerical Poisson’s ratios

Table 3:Calculation results for the numerical Poisson’s ratios

From Table 3,one can find that the errors of the OSPD-BEM are much smaller than those of PD-MPM.Besides,the numerical Poisson’s ratios predicted by the PD-MPM are always lower than the given Poisson’s ratio,which can be due to a spurious boundary softening phenomena [57].This phenomenon is observed from the results of PD-MPM in Fig.10.We find that Poisson’s ratio is variable for the OSPD in Table 3,while it is the fixed value for the BBPD.The BBPD-BEM is one order of magnitude more precise than the PD-MPM,while the magnitude is two orders for the OSPDBEM.The reason is that the spurious boundary softening phenomena [57] is more serious for the OSPD compared with the BBPD.

Figure 10:Displacement predicted by the PD-MPM and the OSPD-BEM for ν=0.3(color online)

4.2 Two-Dimensional Crack Initiation in a Double-Notched Specimen

In this example,we consider a crack initiation problem in the square plate with a double-notched edge crack exerted uniaxial tensile loading displayed in Fig.11.The calculation parameters are exhibited in Table 4.Edenotes elastic modulus;Crdenotes initial length of the crack;Gcdenotes critical energy release rate.According to a definition of critical energy release rate in PD[17],the load corresponding to the crack propagatingΔlength is defined as the utmost loadFm.Δis the grid spacing[17].For different Poisson ratios,the results predicted by the OSPD-BEM are compared with those of the PD-MPM.

Figure 11:The diagram for a double-notched specimen(color online)

Table 4:Calculation parameters for the double-notched specimen

For Poisson’s ratioν=0.3 and horizon hr=0.0125 m,the displacements predicted by the OSPD-BEM are compared with those predicted by the PD-MPM in Figs.12.One can find that both displacements in the directionxandyare consistent with those predicted by the PD-MPM.But utmost loadFmand calculation timetof both methods are different.For different horizons hr,the utmost loadFmand the calculation timetpredicted by OSPD-BEM and PD-MPM are displayed in Figs.13 and 14.

Figure 12:(Continued)

Figure 13:The load Fm vs.the reciprocal of the horizon 1/hr(color online)

Figure 14:The calculation time t vs.the reciprocal of the horizon 1/hr,ν=0.3(color online)

In Fig.13,we find the utmost loadFmacquired through OSPD-BEM can be higher than the one predicted through PD-MPM,and their results are consistent with each other when horizon hrapproaches zero.This phenomenon is consistent with the results in BBPD[21].The utmost loadFmdecreases with an increment of Poisson’s ratio,which agrees well with results in this literature [54].The efficiency is demonstrated in Fig.14.The calculation time for PD-MPM sharply improves as a reduction of horizon hr,while the one for OSPD-BEM does not change much.Meanwhile,the asymptotic compatibility[58]has also been demonstrated for the OSPD-BEM.

4.3 Wedge-Splitting Test

The wedge-splitting test is a significant experiment to investigate the fracture toughness of the concrete specimen[55,56].Many numerical methods based on classical continuum mechanics[59,60]and PD theory[32,61]have been applied to simulate this experiment.In this section,the OSPD-BEM is used to simulate this experiment process,and our result will be compared with those of the ordinary state-based peridynamic mesh-free particle method (OSPD-MPM) predicted by Bie et al.[61].The schematics of the wedge-splitting test specimen is displayed in Fig.15.The relative displacementthat is perpendicular to a crack is applied on the pre-cracked edge of a square plate,and hinge supports are applied to the symmetric points on the edge opposite the crack.The geometry and materials parameters are shown in Table 5.Edenotes elastic modulus;νdenotes Poisson’s ratio;hrdenotes horizon;G0denotes elastic limit energy density [32];GICdenotes critical energy release rate [51];Crdenotes initial length of a crack.

Figure 15:The schematics for the wedge-splitting test specimen(color online)

Table 5:Calculation parameters for the wedge-splitting test

To describe the propagation process of a crack,the curve of the crack lengthvs.the crack mouth opening displacement (CMOD) is given in Fig.16.Some specific crack propagation diagrams are shown in Fig.17.The main purpose of the wedge-splitting test is to obtain the curve about the external loadvs.CMOD.This curve is related to the fracture toughness of a concrete specimen[55].Drawing upon the crack propagation model of the OSPD-BEM in Section 3,the curve of the loadvs.the CMOD is given in Fig.18.In Fig.18,we find the result of the OSPD-BEM is closer to the experimental result[56]than the one of PD-MPM in[61].The reasons originate from two aspects.First,the OSPD-BEM is more accurate as a semi-analytical calculation method.Second,the softening process observed in the test[60,62]is introduced into the crack propagation model.

Figure 16:The crack length vs.the CMOD(color online)

Figure 17:The crack propagation diagram(color online)

The crack initiation loadFm(i.e.,the peak value in Fig.18)and its corresponding CMODUmare significant parameters for the fracture toughness of the concrete specimen.We further compare the fracture toughness predicted by our OSPD-BEM and Bie et al.[61]with the experimental results in Table 6.The error is calculated based on the experimental result[56].In Table 6,the accuracies ofFmandUmare respectively improved by one time and seven times when we compare the OSPD-BEM and Bie et al.[61]with the experimental result[56].

Figure 18:The load vs.the CMOD(color online)

Table 6:Fracture toughness comparison

Considering the calculation efficiency,we find that it takes 44000 s to simulate the wedge-splitting test in Fig.15 with the OSPD-MPM [61],while the time consumption has been reduced to 21000 s with the coupling method[61]between the OSPD-MPM and the node-based smoothed finite element method(NS-FEM).Nevertheless,the OSPD-BEM only takes 11000 s to deal with the same problem in Fig.15.The calculation efficiency is respectively improved by three times and one time compared with the OSPD-MPM[61]and the coupling method[61].By the way,200 boundary elements are used in the OSPD-BEM for the simulation result in Fig.15.The computer equipment is displayed as (1)CPU:12th Gen Intel(R)Core(TM)i5-1240P 1.70 GHz;(2)RAM:16.0 GB;(3)OS:Windows 11 Home 22H2 64 bit.

4.4 Three Point Bending Test

The three-point bending experiment is the main experimental approach to measure the mechanical properties of materials[63,64].The beam with a central crack is executed by the three-point bending experimental to investigate fracture toughness in Fig.19.As displayed in Fig.19,represents displacement loading,while the hinge support and the fixed support are applied to both ends of the beam.

The geometry and materials parameters are shown in Table 7.In Table 7,Edenotes elastic modulus;νdenotes Poisson’s ratio;hrdenotes horizon;G0denotes elastic limit energy density [32];GICdenotes critical energy release rate[51];Crdenotes initial length of crack.The three-point bending experiment in Fig.19 is simulated with a crack propagation model of OSPD-BEM in Section 3.To describe the propagation process of the crack,the curve of the crack lengthvs.the CMOD is given in Fig.20.Some specific crack propagation diagrams are shown in Fig.21.

Figure 19:The schematics for the three-point bending test specimen(color online)

Figure 21:The crack propagation diagram(color online)

Figure 22:The load vs.the CMOD(color online)

The curved line of external loadvs.CMOD is given in Fig.22,and we compare OSPDBEM results with FEM results predicted by Carpinteri et al.[62] and the bond-based peridynamic mesh-free particle method (BBPD-MPM) results predicted by Zaccariotto et al.[32,33] in Fig.22.We find the result of the OSPD-BEM is closer to those of Carpinteri et al.[62] than those of Zaccariotto et al.[32,33].The reason stems from the effect of Poisson’s ratio.The OSPD-BEM and Carpinteri et al.[62] can adopt a real Poisson’s ratio,while Zaccariotto et al.[32,33] is trapped in a fixed Poisson’s ratio.As displayed in [65],the effect of Poisson’s ratio becomes obvious,when a crack approaches boundary.The crack initiation loadFm(i.e.,the peak value in Fig.22) and its corresponding CMODUmin Fig.22 corresponding to Carpinteri et al.[62],Zaccariotto et al.[32,33]and the OSPD-BEM are shown in Table 8.The prediction of the fracture toughness parameters is close for these three numerical approaches.

Table 8:Fracture toughness comparison

5 Conclusions

In our research,we propose OSPD-BEM as an enhancement over prior efforts [21].This method inherits the advantages of the BBPD-BEM[21].Firstly,it eliminates the boundary softening phenomenon [57] resulting from boundary discretization.Secondly,it facilitates problem-solving in infinite domains.Thirdly,it exhibits advantages in terms of accuracy and efficiency compared to numerical methods that discretize the entire domain.In addition,this method is based on the OSPD,overcoming theoretical limitations present in previous works[21],such as fixed Poisson’s ratio.Furthermore,we introduce a fracture calculation model tailored for OSPD-BEM to investigate basic crack propagation problems.In future research,a viable approach to tackling complex crack growth problems,such as multiple cracks and crack branching,might involve the utilization of a coupling method[66,67].This method would integrate PD-MPM and OSPD-BEM in a complementary manner,with PD-MPM applied in cracked regions and OSPD-BEM in uncracked regions.

Acknowledgement:We acknowledge the High-Performance Computing Platform of Peking University for providing computational resources.

Funding Statement:The work is supported by the National Key R&D Program of China(2020YFA07 10500).

Author Contributions:The authors confirm contribution to the paper as follows:study conception and design:Xue Liang,Linjuan Wang;data collection:Xue Liang;analysis and interpretation of results:Xue Liang,Linjuan Wang;draft manuscript preparation: Xue Liang,Linjuan Wang.All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials:The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.

Appendix A.The Comparison of Different Nonlocal Operators

The literature[36]and the literature[5,43]provide two different nonlocal operators,and both can effectively represent the PD equilibrium equation and boundary condition.The nonlocal gradient,nonlocal curl and nonlocal divergence in both literature is a dual relationship.The nonlocal operators[5,43]is the weighted adjoint operator of the ones[36],if the following conditions are satisfied:

where the symbol on the left side of Eq.(A.1)appears in Eqs.(6)and(7),and the symbol on the right side of Eq.(A.1)appears in Eq.(6)of the literature[43].

In other words,there are

wherePis a vector,C∗is defined as follows:

The following conditions are satisfied:

Then

where the symbols are seen in Eq.(13)of the literature[43].

Appendix B.The Boundary Integral Equation in the Laplace Domain

It is convenient to deal with the dynamic problem in the Laplace domain for the BEM.It has two prominent advantages,compared with time domain.The first is eliminating time accumulation error;the second is that we can easily implement parallel computation.Thus,we derive the BIE of the OSPD for dynamic problems in Laplace domain as follows.First of all,we give the dynamic equation for the linear elastic OSPD in the Laplace domain

where the operatorsGandMare mentioned in Eq.(11);sis a variable corresponding to the Laplace transformation;ρis the mass density;is the Laplace transformation of the displacementuin the time domain;are the Laplace transformation of the boundary condition;is the Laplace transformation of the body force densityFbin the time domain;t0is the initial time;u0andare the initial conditions in the time domain.

The weighted residual method is adopted to derive the BIE in Laplace domain instead of the reciprocal theorem,following the derivation in the classical theory.We rewrite Eq.(B.1)in the form of weighted residual as follows:

Substituting Eq.(B.6)into Eq.(B.2),we obtain

Substituting Eqs.(B.1)and(B.3)into Eq.(B.7)and dividing the integral domainΩτin the third integral on the right hand side of Eq.(B.7)into two parts,we get

Considering the extra condition Eq.(16),which converts the integral in the volume constrained boundary to the one on the classical boundary,we simplify Eq.(B.8)as follows:

we get as follows:

whereD(x,s)andP(x,s)are the local boundary condition in the Laplace domain.Substituting Eqs.(B.11)and(B.12)into Eq.(B.10)yields

Eq.(B.13)is the BIE of the OSPD in the Laplace domain.