Dual-horizon Peridynamics
DDual-horizon Peridynamics
Huilong Ren a,b , Xiaoying Zhuang a,b , Yongchang Cai a , Timon Rabczuk b,c a State Key Laboratory of Disaster Reduction in Civil Engineering, College of Civil Engineering,TongjiUniversity, Shanghai 200092, China b Institute of Structural Mechanics, Bauhaus-Universit¨at Weimar, Germany c School of Civil, Environmental and Architectural Engineering, Korea University, South Korea
Abstract
In this paper we develop a new Peridynamic approach that naturally includes varying horizonsizes and completely solves the ”ghost force” issue. Therefore, the concept of dual-horizonis introduced to consider the unbalanced interactions between the particles with differenthorizon sizes. The present formulation is proved to fulfill both the balances of linear momen-tum and angular momentum. Neither the ”partial stress tensor” nor the ”‘slice” techniqueare needed to ameliorate the ghost force issue in [1]. The consistency of reaction forces isnaturally fulfilled by a unified simple formulation. The method can be easily implementedto any existing peridynamics code with minimal changes. A simple adaptive refinementprocedure is proposed minimizing the computational cost. The method is applied here tothe three Peridynamic formulations, namely bond based, ordinary state based and non-ordinary state based Peridynamics. Both two- and three- dimensional examples includingthe Kalthof-Winkler experiment and plate with branching cracks are tested to demonstratethe capability of the method in solving wave propagation, fracture and adaptive analysis .
Keywords:
Peridynamics, horizon variable, dual-horizon, ghost force, spurious wavereflection, adaptive refinement ∗ Corresponding author: Tel: +44 191 334 2504
Email addresses: [email protected]; (Xiaoying Zhuang), [email protected]; (Timon Rabczuk)
Preprint submitted to Elsevier June 23, 2015 a r X i v : . [ phy s i c s . c o m p - ph ] J un . Introduction Peridynamics has recently attracted wide interests for researchers in computational solidmechanics since it provides the possibility to model dynamic fracture with ease. The crack ispart of the solution from PD simulation instead of part of the problem and no representationof the crack topology is needed. The original PD method was proposed by Silling [2] in 2000and has been exploited onwards for extensive applications of mechanical problems includingimpact loading, fragmentation [3, 4], composites delamination [5], beam and plate structures[6, 7].In PD, the classical balance equations are formulated in an integral form instead ofpartial differential form. The particles interact with each other when the distances betweenthe particles are within a threshold value, called horizon . The equations of motion at anytime t is expressed as ρ ¨u ( x , t ) = (cid:90) H x f ( u ( x (cid:48) , t ) − u ( x , t ) , x (cid:48) − x ) − f ( u ( x , t ) − u ( x (cid:48) , t ) , x − x (cid:48) ) d V x (cid:48) + b ( x , t ) , (1)where H x denotes the horizon (spherical domain) belonging to x , u is the displacementvector, b denotes the body force, ρ is mass density in the reference configuration, and f is a pairwise force function that computes the force vector (per unit volume squared) [8].Eq.(1) shows the key idea of PD to unify continuous and discontinuous media within a singleconsistent set of equations. The governing equation is written in an integral form insteadof an partial differential form. The crack surfaces are formed as the outcome of motionand constitutive models, and there is no entanglement of additional crack kinematics orgeometry treatment. The fracture behaviour including crack branching and coalescence ofmultiple cracks is captured through the breakage of the bonds between particles. Therefore,the need of smoothing of crack surfaces or branching criterion in the extended finite elementmethod (XFEM) [9], meshless methods [10] or other partition of unity methods (PUM) [11]is completely removed. The extension of PD from 2D to 3D problems is greatly facilitatedfor the computer implementation, which is not always the case in other methods [12, 13].2he original PD started with the bond-based formulation (BB-PD) where the bondsbehave like springs and independent of each other. BB-PD can be regarded as a specialcase of a more general theory, the state-based peridynamics (SB-PD) [2, 14] which can besuited for, theoretically, any type of constitutive model and large deformation analysis. Thekey difference between SB-PD and BB-PD is that in the former, the bond deformationdepends on collective deformation of other bonds, whereas the bonds in the latter deformsindependently. The SB-PD was later extended into two types, namely the the ordinary statebased peridynamics (OSB-PD) and the non-ordinary state based peridynamics (NOSB-PD).In all the above types of PD, horizons sizes are commonly required to be constant to avoidspurious wave reflections and ghost forces between particles. However, in many applications,the spatial distribution of the particles with changing horizon sizes is necessary, e.g. adaptiverefinement, multiscale modelling and multibody analysis. In other words, in order to achieveacceptable accuracy, the entire numerical model has to be discretised with respect to thehighest particle resolution locally required, and the smallest horizon size used accordingly.This is computationally expensive and undesirable. The restriction of horizons being positionindependent practically reduces the efficiency of PD.In this paper, we aim to remove the issue of varying horizons and ghost force by de-veloping a new PD formulation. The new approach is based on the concept of horizonand dual-horizon. Though peridynamics has been developed for different types of physicalfields, e.g. thermal field and fluid field, we confine the present work to solving solid me-chanics problems. The content of the paper is outlined as follows. § §
3, the horizon and dual-horizon are introduced.New motion equations with varying horizons are derived based on the dual-horizon concept.The balance of linear momentum and the balance of angular momentum of the present PDformulation are proved. The applications of the new formulation for BB-PD, OSB-PD andNOSB-PD are described in details in §
4. In §
5, three numerical examples are presented tovalidate the present method. 3 . Ghost force and spurious wave reflection in peridynamics
In the original PD theory, the forces exerted on a particle is the summation of all thepairwise forces from the particles falling inside the horizon of that particle. When the horizonsizes are set constant for all particles, the bond force is always pairwise. On the other hand,it should be noted that the size of the particles that represent the mass quantity can vary[2, 15]. Spurious wave reflections emerge when the horizon sizes vary. We illustrate themechanism of the phenomenon as follows. Consider a particle x (cid:48) falling inside the horizon ofparticle x , see Fig.1. Let f xx (cid:48) denote the force vector acting on particle x due to particle x (cid:48) ,where the first subscript x indicates x being the object of force and the second subscript x (cid:48) indicates f xx (cid:48) being the source of force from x (cid:48) . Take the Fig. 1 for an example, as x (cid:48) ∈ H x ,the force f x (cid:48) x (cid:54) = 0, x / ∈ H x (cid:48) due to unequal size of horizons, hence f xx (cid:48) = 0. When computingthe reactive forces for particle x , the bond forces F xx (cid:48) = f xx (cid:48) − f x (cid:48) x = − f x (cid:48) x , which is addedto F x . Likewise, when computing forces for particle x (cid:48) , particle x exerts no force on x (cid:48) as x isnot inside the horizon of x (cid:48) . Consequently, the bond force F xx (cid:48) only exists unilaterally, whichis known as the “ghost force” [1] resulting in an unbalanced internal force. The balance oflinear momentum and balance of angular momentum in this case are violated, and henceyields spurious wave reflections in PD simulations.Efforts have been made by researchers to explore the possibility of making PD suitablefor nonuniform spatial discretisation. Bobaru et. al studied the convergence and adaptiverefinement in 1D peridynamics [16] and multi-scale modeling in 2D BB-PD [17]. Dipasqualeet al. recently [18] introduced a trigger based on the damage state of the material for 2Drefinement of BB-PD. Refined PD was developed in [19, 16, 17, 18, 20], however it is restrictedto BB-PD formulation and to the authors knowledge, has only been used for one- and two-dimensional problems. In these works, the spurious wave reflection or ghost force problem isnot solved, and the refinement is performed by checking the spurious reflection is within anacceptable range compared to the magnitude of the whole wave. Recently, the partial stress was proposed in [1] to remove the ghost force for varying horizons. However, the methodimposes certain restrictions to the deformation of the body and requires the computation of4 x H x' x x' -f x'x f x'x f xx' -f xx' Figure 1: Force vector in peridynamics with varying horizons. Reactive force − f x (cid:48) x on x due to x exertedon x (cid:48) as x (cid:48) ∈ H x . partial stress which is complicated. It to certain extent impairs the simplicity of the originalPD, especially for BB-PD and OSB-PD. Besides, the method does not completely removethe ghost force but with a small residual which is believed to be acceptable. Though theslice method was devised in the same work, it is applicable only to piecewise constant sizes ofhorizons, and requires additional computation to enforce the consistency between particlesclose to the interface.
3. Governing equations based on horizon and dual-horizon
In this section, the concept of horizon and dual-horizon will be introduced to formulatethe balance equations for particles with varying horizons. It is applied to all peridynamicformulations.
We will begin by restating the original concept of horizon in peridynamics. In peridy-namics theory, the particles interact with each other within a finite distance. A particle isconsidered to have an influence over other particles within a small domain centering thatparticle. The radius of the domain is known as “horizon”, see Fig.2. Particle x (thick solidline) is included in the horizons of particles x , x , x and x (thin solid line), however is5ot included in the horizons of particles x and x (dashed line). The concept of horizoncan be compared to the concept of “nodal support” in meshless methods. x H X x x x x x x Figure 2: The schematic diagram for horizon and dual-horizon, all circles above are horizons. The greenpoints { x , x , x , x } ∈ H (cid:48) x ,whose horizons denote by thin solid line; the red points { x , x } / ∈ H (cid:48) x ,whosehorizons denote by dashed line Horizon
The horizon H x is the domain where any particle falling inside will receive the forces exertedby x . Hence, x will undertake all the reactive or passive forces from particles in H x . Hence,the horizon can be viewed as passive force horizon. The reactive force acting on x followsNewton’s third law. As shown in Fig.1, the reactive force − f x (cid:48) x exerted on x by otherparticles is in opposite direction of the forces applied by x to other particles. Dual-horizon
Dual-horizon is defined as a union of points whose horizons include x , denoted by H (cid:48) x = { x (cid:48) : x ∈ H x (cid:48) } . In the notation of dual-horizon H (cid:48) x , the superscript prime indicates “dual”,and the subscript x denotes the object particle. It can be understood as the set of thehorizons that belong to the particles who can “see” x in their horizons. As shown in Fig.2,the dual-horizon with respect to x is the union of x , x , x and x , whose horizons aredenoted by thin solid circles. Particles x and x are not included in the dual-horizon of x since their horizons do not include x . In this case, x becomes the object “observed” by6he other particles. If x is within the horizon of x (cid:48) , then x (cid:48) has an active or direct effect on x , corresponding to the passive effect in horizon defined previously. Particle x receives theactive forces from other particles in the dual-horizon of x , and in this sense it is consideredas “dual” corresponding to horizon . For any point x , the shape of H (cid:48) x is arbitrary, anddepends on the sizes and shapes of horizons as well as the locations of the particles. Notethat the horizon can take other shapes other than circles or spheres.The bond inside the dual-horizon ( H (cid:48) x ) is termed as “dual-bond”, and this is correspond-ing to the bond in the horizon ( H x ). It can be seen that the bond of one particle can becomethe dual-bond for another particle interacting with it. Note that for each particle, the bondand the dual-bond are independent from each other; the same applies to the horizon anddual-horizon. It means the bond and dual-bond can break independently in the fracturemodels. For discretisations with varying horizons, often one particle is within the horizonsof other particles but not vice versa. In this case, a single horizon is not sufficient to definethe interactions between particles. The concept of two horizons proposed here can solvethis issue with a simple and direct physical meaning. It naturally takes into account theinteractions between particles of varying horizon sizes. For models with constant horizons,horizon and dual-horizon will degenerate to the horizon in original peridynamics. Reaction force by horizon and dual-horizon
In our dual-horizon peridynamic formulation, the horizons are differentiated between how aparticle receives and exerts forces with other particles. Under this new concept, computingthe force f xx (cid:48) between a pair of particles, denoted as particles x and x (cid:48) , is determined bywhether x is within the horizon of x (cid:48) , and whether x (cid:48) is inside the dual-horizon of x , and viceversa for f x (cid:48) x . It can be easily seen that the force density vector has the following property,if x ∈ H x (cid:48) or x (cid:48) ∈ H (cid:48) x , f xx (cid:48) (cid:54) = 0 , else f xx (cid:48) = 0 . (2)For any bond between two particles x and x (cid:48) belonging to a domain denoted as Ω in Fig.3,the direct force f xx (cid:48) acting on x due to x (cid:48) can be computed by two approaches as follows.7 pproach 1 computes the force in terms of x , f xx (cid:48) = (cid:54) = 0 if x (cid:48) ∈ H (cid:48) x x (cid:48) / ∈ H (cid:48) x , and Approach 2 is formulated with respect to x (cid:48) , f xx (cid:48) = (cid:54) = 0 if x ∈ H x (cid:48) x / ∈ H x (cid:48) . � x �� x ' x - x ' H X' H X f x'x f xx' � Figure 3: Double summation of force
For any domain under consideration e.g. the shaded area in Fig.3, the computation ofthe direct forces that take place between particles (no reactive force considered yet) shallundertake all forces from any particle that belongs to Ω . There are two approaches toachieve that. The first approach is by summing all the forces the particles undertake in thedual-horizon, (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H (cid:48) x f xx (cid:48) (cid:52) V x (cid:48) (cid:52) V x . (3)And the second is to add up all the forces of each particle that it applies to other particles, (cid:88) x (cid:48) ∈ Ω (cid:88) x ∈ H x (cid:48) f xx (cid:48) (cid:52) V x (cid:52) V x (cid:48) . (4)8ince the total force for any Ω is independent of the approach chosen to compute it, Eq. (3)and (4) shall be equal (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H (cid:48) x f xx (cid:48) (cid:52) V x (cid:48) (cid:52) V x = (cid:88) x (cid:48) ∈ Ω (cid:88) x ∈ H x (cid:48) f xx (cid:48) (cid:52) V x (cid:52) V x (cid:48) . (5)Eq. (5) indicates that changing the summation order from x → x (cid:48) shall be done together withchanging from H x (cid:48) → H (cid:48) x so as to keep all the variables consistent. When the discretisationis sufficiently fine, the summation shall approximate to the integral and Eq. (5) becomes (cid:90) x ∈ Ω (cid:90) x (cid:48) ∈ H (cid:48) x f xx (cid:48) d V x (cid:48) d V x = (cid:90) x (cid:48) ∈ Ω (cid:90) x ∈ H x (cid:48) f xx (cid:48) d V x d V x (cid:48) . (6) In the following derivation x (cid:48) − x will be denoted as ξ , and u ( x (cid:48) , t ) − u ( x , t ) as η , hence ξ + η represents the current relative position vector between the particles. The internalforces that are exerted at each particle from the other particles should include two parts,namely the forces from the horizon and the forces from the dual-horizon. The other forcesapplied to a particle include the body force and the inertia force. Let (cid:52) V x denote the volumeassociates to x . The body force for particle x can be expressed as b ( x , t ) (cid:52) V x , where b ( x , t )is the body force density. The inertia is denoted by ρ ¨u ( x , t ) (cid:52) V x , where ρ is the densityassociated to x . At any time t , for x (cid:48) in dual-horizon of x , the force vector of ˜ f xx (cid:48) is definedas ˜ f xx (cid:48) := f xx (cid:48) ( η , ξ ) · (cid:52) V x · (cid:52) V x (cid:48) , (7)where f xx (cid:48) ( η , ξ ) is the force density in the traditional peridynamics with unit of force pervolume squared; ˜ f xx (cid:48) is the force vector acting on particle x due to the attraction or repulsionfrom x (cid:48) . The forces from the dual-horizon are active forces as they are applied to x . Thetotal force applied to x from its dual-horizon, H (cid:48) x , can be computed by, (cid:88) x (cid:48) ∈ H (cid:48) x ˜ f xx (cid:48) = (cid:88) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) · (cid:52) V x · (cid:52) V x (cid:48) . (8)For any x (cid:48) inside the horizon of x , the force ˜ f x (cid:48) x acting on x (cid:48) due to x is defined as˜ f x (cid:48) x := f x (cid:48) x ( − η , − ξ ) · (cid:52) V x (cid:48) · (cid:52) V x (9)9here f x (cid:48) x ( − η , − ξ ) is the force density per volume squared; ˜ f x (cid:48) x is the force vector actingon particle x (cid:48) due to the attraction or repulsion from x ( x (cid:48) inside horizon of x ). The totalforce x exerted on other particles from horizon H x is the summation of ˜ f x (cid:48) x (cid:88) x (cid:48) ∈ H x ˜ f x (cid:48) x = (cid:88) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) · (cid:52) V x · (cid:52) V x (cid:48) (10)Based on Newton’s third law, the total force x undertaken in H x takes opposite direction − (cid:88) x (cid:48) ∈ H x ˜ f x (cid:48) x = − (cid:88) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) · (cid:52) V x · (cid:52) V x (cid:48) (11)By summing over all forces on particle x , including inertia force, body force, active force inEq. (8) and passive force in Eq. (11), we obtain the equations of motion ρ ¨u ( x , t ) (cid:52) V x = (cid:88) x (cid:48) ∈ H (cid:48) x ˜ f xx (cid:48) + (cid:32) − (cid:88) x (cid:48) ∈ H x ˜ f x (cid:48) x (cid:33) + b ( x , t ) (cid:52) V x . (12)Substituting Eqs. (8) and (11) into Eq. (12) leads to ρ ¨u ( x , t ) (cid:52) V x = (cid:88) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) (cid:52) V x (cid:48) (cid:52) V x − (cid:88) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) (cid:52) V x (cid:48) (cid:52) V x + b ( x , t ) (cid:52) V x . (13)As the volume (cid:52) V x associated to particle x is independent of the summation, we can elim-inate (cid:52) V x in Eq. (13), yielding the governing equation based on x : ρ ¨u ( x , t ) = (cid:88) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) (cid:52) V x (cid:48) − (cid:88) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) (cid:52) V x (cid:48) + b ( x , t ) . (14)When the discretisation is sufficiently fine, the summation is approximating to the integra-tion of the force on the dual-horizon and horizon,lim (cid:52) V x (cid:48) → (cid:88) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) (cid:52) V x (cid:48) = (cid:90) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) d V x (cid:48) (15)and lim (cid:52) V x (cid:48) → (cid:88) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) (cid:52) V x (cid:48) = (cid:90) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) d V x (cid:48) . (16)10hus the integration form of the equation of motion in peridynamics with dual horizon isgiven as ρ ¨u ( x , t ) = (cid:90) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) d V x (cid:48) − (cid:90) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) d V x (cid:48) + b ( x , t ) . (17)Eq. (17) is similar to Eq. (1) of the original peridynamics theory. When the horizons areset constant, i.e. both horizon and the dual-horizon are equal, the integrations in Eq. (17)degenerate to the original peridynamics theory. It means the traditional peridynamics canbe viewed as a special case of the present dual horizon peridynamics.About the implementation of the present peridynamic formulation, for any particle x ,the force density f xx (cid:48) ( η , ξ ) in H (cid:48) x can be determined when calculating the force in H x (cid:48) forparticle x (cid:48) . Therefore, it is not necessary to know exactly the dual-horizon geometry andthe formulation can be implemented with minor modification of any peridynamic codes. The internal forces shall satisfy the balance of linear momentum for any bounded body Ω given by (cid:90) Ω ( ρ ¨u ( x , t ) − b ( x , t )) d V x = (cid:90) Ω (cid:90) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) d V x (cid:48) d V x − (cid:90) Ω (cid:90) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) d V x (cid:48) d V x = . (18)For simplicities, let f xx (cid:48) represent f xx (cid:48) ( η , ξ ), and f x (cid:48) x the f x (cid:48) x ( − η , − ξ ). Proof : For convenience, we start with the discrete form of the Eq. (14) as (cid:88) x ∈ Ω ( ρ ¨u ( x , t ) − b ( x , t )) (cid:52) V x = (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H (cid:48) x f xx (cid:48) (cid:52) V x (cid:48) (cid:52) V x − (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x f x (cid:48) x (cid:52) V x (cid:48) (cid:52) V x (19)The first term on RHS of Eq. (19) can be substituted with Eq. (5) by changing thesummation domain and therefore we will get, (cid:88) x ∈ Ω ( ρ ¨u ( x , t ) − b ( x , t )) (cid:52) V x = (cid:88) x (cid:48) ∈ Ω (cid:88) x ∈ H x (cid:48) f xx (cid:48) (cid:52) V x (cid:48) (cid:52) V x − (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x f x (cid:48) x (cid:52) V x (cid:48) (cid:52) V x , (20)11he above transformation is based on the idea that f xx (cid:48) acting on x is computed by usingdual-horizon ( H (cid:48) x ) of x . And f xx (cid:48) is calculated by using horizon ( H x (cid:48) ) of x (cid:48) . Note that inboth calculations, f xx (cid:48) remains unchanged and only the definition domain changes. Sincethe dummy variables can be relabeled for x ↔ x (cid:48) in the first term on RHS, thus the firstterm has the same expression as the second term, i.e. (cid:88) x ∈ Ω ( ρ ¨u ( x , t ) − b ( x , t )) (cid:52) V x = (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x f x (cid:48) x (cid:52) V x (cid:52) V x (cid:48) − (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x f x (cid:48) x (cid:52) V x (cid:48) (cid:52) V x = 0As the forces are differentiated here between active and passive force, i.e. an activeforce from x corresponds to a passive force of x (cid:48) , the forces is always pairwise but withopposite direction and of the same magnitude. Hence, it is natural that the balance oflinear momentum is satisfied. Let y be the deformation vector state field defined by y [ x , t ] (cid:104) ξ (cid:105) = y ( x (cid:48) , t ) − y ( x , t ) ∀ x ∈ Ω , ξ = x (cid:48) − x , x (cid:48) ∈ H x , t ≥ u ( x , t ) = y ( x , t ) − x , (22)where y ( x , t ) refers to the current configuration coordinate for x in material configuration, u ( x , t ) is the displacement for x , H x is the horizon. Thus y (cid:104) x (cid:48) − x (cid:105) is the image of the bond x (cid:48) − x under the deformation.To satisfy the balance of angular momentum for any bounded body Ω ,it is required that: (cid:90) Ω y × ( ρ ¨u ( x , t ) − b ( x , t )) d V x = (cid:90) Ω y × (cid:18)(cid:90) x (cid:48) ∈ H (cid:48) x f xx (cid:48) ( η , ξ ) d V x (cid:48) − (cid:90) x (cid:48) ∈ H x f x (cid:48) x ( − η , − ξ ) d V x (cid:48) (cid:19) d V x = (cid:90) Ω (cid:90) x (cid:48) ∈ H (cid:48) x y × f xx (cid:48) ( η , ξ ) d V x (cid:48) d V x − (cid:90) Ω (cid:90) x (cid:48) ∈ H x y × f x (cid:48) x ( − η , − ξ ) d V x (cid:48) d V x = (23)For simplicities, let f xx (cid:48) represent f xx (cid:48) ( η , ξ ), and f x (cid:48) x present f x (cid:48) x ( − η , − ξ ).12 roposition : In the dual-horizon peridynamics, suppose a constitutive model of the form f = ˆ f ( y , Λ) (24)where ˆ f : V → V is bounded and Riemann-integrable on H and V is the vector state; Λdenotes all variables other than the current deformation vector state that may depend onfor some particular material.If (cid:90) x (cid:48) ∈ H x y (cid:104) x (cid:48) − x (cid:105) × f x (cid:48) x d V x (cid:48) = 0 ∀ y ∈ V , (25)or in discrete form (cid:88) x (cid:48) ∈ H x y (cid:104) x (cid:48) − x (cid:105) × f x (cid:48) x (cid:52) V x (cid:48) = 0 ∀ y ∈ V , (26)where f x (cid:48) x = f x (cid:48) x ( u ( x (cid:48) , t ) − u ( x , t ) , x (cid:48) − x ) , the force vector density acting on x (cid:48) , and H x is thehorizon of x . Then, the balance of angular momentum, Eq. (23) holds for any deformationof Ω for any given constitutive model. Proof : As the summation over the domain Ω is equivalent to the integrand of the sameexpression over that domain, we use the discrete form as (cid:88) x ∈ Ω y × ( ρ ¨u ( x , t ) − b ( x , t )) (cid:52) V x = (cid:88) x ∈ Ω ( x + u ) × ( ρ ¨u ( x , t ) − b ( x , t )) (cid:52) V x = (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H (cid:48) x ( x + u ) × f xx (cid:48) (cid:52) V x (cid:48) (cid:52) V x − (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x ( x + u ) × f x (cid:48) x (cid:52) V x (cid:48) (cid:52) V x (27)The first term on the RHS of Eq. (27) can be rewritten according to Eq. (5) by changingthe summation domain and thus it becomes (cid:88) x ∈ Ω y × ( ρ ¨u ( x , t ) − b ( x , t )) (cid:52) V x = (cid:88) x (cid:48) ∈ Ω (cid:88) x ∈ H x (cid:48) ( x + u ) × f xx (cid:48) (cid:52) V x (cid:48) (cid:52) V x − (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x ( x + u ) × f x (cid:48) x (cid:52) V x (cid:48) (cid:52) V x . (28)13wapping the dummy variables for the first term on RHS between x ↔ x (cid:48) in the doublesummation, then the first term takes the same summation to the second term on RHS (cid:88) x ∈ Ω y × ( ρ ¨u ( x , t ) − b ( x , t )) (cid:52) V x = (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x ( x (cid:48) + u (cid:48) ) × f x (cid:48) x (cid:52) V x (cid:52) V x (cid:48) − (cid:88) x ∈ Ω (cid:88) x (cid:48) ∈ H x ( x + u ) × f x (cid:48) x (cid:52) V x (cid:48) (cid:52) V x = (cid:88) x ∈ Ω ( (cid:88) x (cid:48) ∈ H x (( x (cid:48) + u (cid:48) ) − ( x + u )) × f x (cid:48) x (cid:52) V x (cid:48) ) (cid:52) V x = (cid:88) x ∈ Ω ( (cid:88) x (cid:48) ∈ H x ( y (cid:48) − y ) × f x (cid:48) x (cid:52) V x (cid:48) ) (cid:52) V x = (cid:88) x ∈ Ω ( (cid:88) x (cid:48) ∈ H x y (cid:104) x (cid:48) − x (cid:105) × f x (cid:48) x (cid:52) V x (cid:48) ) (cid:52) V x = (29)or in an integration form for sufficiently fine discretisation as (cid:90) x ∈ Ω y × ( ρ ¨u ( x , t ) − b ( x , t )) d V x = . (30)Therefore, the angular momentum over the entire analysis domain is satisfied. (cid:3) x x f xx' f x'x ' x x �� f xx' f x'x ' x - x ' �� x - x ' BB-PD and OSB-PD NOSB-PD
Figure 4: force vector of BB-PD,OSB-PD and NOSB-PD
It can be seen that the conservation of angular somehow depends only on the horizon;the dual-horizon is not involved. The latter is only needed in the Eq. (14) and Eq. (17).The bond-based, ordinary based and non-ordinary based peridynamics all satisfy the angular14omentum in the original horizon concept (see proof in [15]). This conclusion can be alsoillustrated for the BB-PD and OS-PD since the internal forces f xx (cid:48) and f x (cid:48) x are parallel tothe bond vector in the current configuration, see Fig.4.
4. Dual-horizon peridynamics
The dual-horizon formulation will now be applied to all existing peridynamics, namelyBB-PD, OSB-PD and NOSB-PD. The calibration of constitutive parameters with respectto the continuum model and some issues concerning the implementations will be discussedfor 2D and 3D problems.
In the bond-based peridynamics theory [4, 8], the pair force function is expressed as f ( η , ξ ) = ∂w ( η , ξ ) ∂ η ∀ η , ξ (31)For a micro-elastic homogenous and isotropic material, f ( η , ξ ) can be specified as f ( η , ξ ) = cs · η + ξ (cid:107) η + ξ (cid:107) , (32)where c is the micro-modulus, and s is the bond stretch calculated by s = (cid:107) η + ξ (cid:107) − (cid:107) ξ (cid:107)(cid:107) ξ (cid:107) . (33)The energy per unit volume in the body at a given time t is given by [8] W = 12 (cid:90) H x w ( η , ξ ) d V ξ . (34)For the BB-PD theory, by enforcing the strain energy density being equal to the strainenergy density in the classical theory of elasticity [4, 8], and by setting the influence function ω ( (cid:107) ξ (cid:107) ) = 1, we obtain C ( δ ) = 3 Eπδ (1 − ν ) plane stress C ( δ ) = 3 Eπδ (1 + ν )(1 − ν ) plane strain C ( δ ) = 3 Eπδ (1 − ν ) 3D . (35)15ote that the value of C ( δ ) takes half of the micro-modulus c used in the horizon-constantBB-PD since the bond energy for varying horizon is determined by both horizon and dual-horizon. It can be easily seen that when the horizon takes the same value as the dual-horizon,the bond energy between two particles is reduced to that of the original BB-PD.Let w ( ξ ) = C ( δ ) s ( δ ) ξ/ s ( δ )is the critical bond stretch. By breaking half of all the bonds connected to a given particlealong the fracture surface and equalizing the breaking bonds energy with the energy releaserate G [18, 4], we can get the expression between the energy release rate G and the criticalbond stretch s ( δ ): s ( δ ) = (cid:114) πG Eδ plane stress s ( δ ) = (cid:114) πG Eδ plane strain s ( δ ) = (cid:114) G Eδ . (36)Both the micro-modulus C ( δ ) and the critical stretch s ( δ ) are derived from the local con-tinuum mechanics theory, and they depend on the horizon radii for variable horizons.In the implementation of the bond-based peridynamics, fracture is introduced by remov-ing particles from the neighbour list once the bond stretch exceeds the critical bond stretch s . In order to specify whether a bond is broken or not, a history-dependent scalar valuedfunction µ is introduced [8], µ ( t, ξ ) = s ( t (cid:48) , ξ ) < s for all 0 ≤ t (cid:48) ≤ t, . (37)The local damage at x is defined as φ ( x , t ) = 1 − (cid:82) H x µ ( x , t, ξ ) d V ξ (cid:82) H x d V ξ . (38)The damage formulation of Eq. (38) is also applicable to OSB-PD. For any particle withdual-horizon ( H (cid:48) x ) and horizon ( H x ), the active force f xx (cid:48) and the passive force f x (cid:48) x in Eq.1614) or (17) are computed by the following expressions, respectively f xx (cid:48) = C ( δ x (cid:48) ) · s xx (cid:48) · η + ξ (cid:107) η + ξ (cid:107) , ∀ x (cid:48) ∈ H (cid:48) x (39) f x (cid:48) x = C ( δ x ) · s xx (cid:48) · − ( η + ξ ) (cid:107) η + ξ (cid:107) , ∀ x (cid:48) ∈ H x , (40)where C ( δ x (cid:48) ) and C ( δ x ) are the micro-modulus based on δ x (cid:48) and δ x computed from Eq. (35)respectively, and s xx (cid:48) is the stretch between particles x and x (cid:48) . The concept of “state” for peridynamics was firstly introduced in [15]. A state of order m is a function A : H → T m , ξ (cid:55)→ A (cid:104) ξ (cid:105) . (41)where H is the horizon domain, T m denotes the set of all tensors of order m , (cid:104) (cid:3) (cid:105) indicatethe vector to which a state operates. In this paper, the scalar state and vector state are useif not otherwise specified.The 2D OSB-PD formulation of bond force can be derived by following the similarprocedure as that was described in 3D. Both 2D and 3D OSB-PD can be written in aunified expression taking into account the dimension number as t (cid:104) ξ (cid:105) = nKθm ω (cid:104) ξ (cid:105) · ξ + n ( n + 2) Gm ω (cid:104) ξ (cid:105) e d (cid:104) ξ (cid:105) , (42)where n ∈ { , } is the dimensional number, K is the bulk modulus, G the shear modulus, ξ = (cid:107) ξ (cid:107) , e d (cid:104) ξ (cid:105) = e (cid:104) ξ (cid:105) − θξn , m = (cid:82) H ω (cid:104) ξ (cid:105) ξ · ξ d V ξ , θ = nm (cid:90) H ω (cid:104) ξ (cid:105) ξ · e (cid:104) ξ (cid:105) d V ξ , e (cid:104) ξ (cid:105) = (cid:107) ξ + η (cid:107) − (cid:107) ξ (cid:107) and d V ξ are the area in 2D, and volume in 3D, respectively. For any particle x , f xx (cid:48) and f x (cid:48) x in Eqs. (14) and (17) can be calculated by f xx (cid:48) = t (cid:104) ξ (cid:105) · η + ξ (cid:107) η + ξ (cid:107) , ∀ x (cid:48) ∈ H (cid:48) x , (43)and f x (cid:48) x = t (cid:48) (cid:104)− ξ (cid:105) · − η − ξ (cid:107) η + ξ (cid:107) , ∀ x (cid:48) ∈ H x . (44)Here ξ = x (cid:48) − x , and x and x (cid:48) are the coordinate vectors for x and particle x (cid:48) in the materialconfiguration, respectively. 17 .3. Dual-horizon non-ordinary state based peridynamics NOSB-PD uses the deformation gradient from classic continuum mechanics. It offers thepossibility to consistently include existing constitutive models. With the use of the shapetensor proposed in the NOSB-PD, the fundamental tensors in continuum mechanics canbe easily introduced in peridynamics, including the deformation gradient tensor or velocitygradient [21, 22]. The shape tensor for particle x is defined as K x = (cid:90) x (cid:48) ∈ H x ω ( ξ ) · ξ ⊗ ξ d V x (cid:48) (45)where ξ = x (cid:48) − x is the bond in the reference configuration, H x is the horizon for particle x and ω ( ξ ) is the influence function. The deformation tensor for particle x is given by F x = ∂ y ∂ x = (cid:90) x (cid:48) ∈ H x ω ( ξ ) y (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) · K − x , (46)where y := y ( x , t ) are the spatial coordinates, x are the material coordinates, and ξ = x (cid:48) − x .The spatial velocity gradient L x for particle x is defined in the current configuration usingthe chain rule, L x := ∂ v ∂ y = ∂ v ∂ x · ∂ x ∂ y = ˙ F x F − x , (47)where v := v ( x , t ) = ∂ y ∂t ( x , t ) is the velocity vector and ˙ F x is the rate of deformationgradient. The non-local form of the velocity gradient can be written as L x : = (cid:90) x (cid:48) ∈ H x ω ( ξ ) v (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) · K − x · (cid:18)(cid:90) x (cid:48) ∈ H x ω ( ξ ) y (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) · K − x (cid:19) − = (cid:90) x (cid:48) ∈ H x ω ( ξ ) v (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) · K − x · K x (cid:18)(cid:90) x (cid:48) ∈ H x ω ( ξ ) y (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) (cid:19) − · = (cid:90) x (cid:48) ∈ H x ω ( ξ ) v (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) · (cid:18)(cid:90) x (cid:48) ∈ H x ω ( ξ ) y (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) (cid:19) − (48)where ξ = x (cid:48) − x , v (cid:104) ξ (cid:105) = v ( x (cid:48) , t ) − v ( x , t ) and y (cid:104) ξ (cid:105) = y ( x (cid:48) , t ) − y ( x , t ). It can been seenthat the velocity gradient does not require shape tensor, which means the shape tensor isnot necessary to define velocity gradient. The incremental deformation gradient can also beobtained without shape tensor K . 18et (cid:52) u denotes the displacement increment with respect to the previous time step. Theincremental spatial deformation gradient C x in continuum mechanics is defined as C x := ∂ ( (cid:52) u ) ∂ y = ∂ ( (cid:52) u ) ∂ x · ∂ x ∂ y = ∂ ( (cid:52) u ) ∂ x · F − x , (49)The non-local counterpart is given by C x := (cid:90) x (cid:48) ∈ H x ω ( ξ ) (cid:52) u (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) · ( (cid:90) x (cid:48) ∈ H x ω ( ξ ) y (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) ) − . Many material models in non-ordinary state based peridynamic accounting for the geomet-rical and material non-linear problems are solved on the incremental spatial deformationgradient. The non-ordinary bond force is computed based on the Piola-Kirchhoff as T x (cid:48) x = ω ( ξ ) P x · K − x · ξ (50) P x = J x σ x F − T x , J x = det F x (51)where ξ = x (cid:48) − x , σ x is Cauchy stress tensor.Let us define the deformation gradient F (cid:48) x F (cid:48) x := F x · K x = (cid:90) x (cid:48) ∈ H x ω ( ξ ) y (cid:104) ξ (cid:105) ⊗ ξ d V x (cid:48) · K − x · K x = (cid:90) x (cid:48) ∈ H x ω ( ξ ) y (cid:104) ξ (cid:105) ⊗ ξ d V x . (52)Therefore J x = det F x = det F (cid:48) x det K x . (53)By substituting Eqs. (53), (52) and (51) to (50), we obtain T x (cid:48) x = ω ( ξ ) J x σ x F − T x · K − x · ξ = ω ( ξ ) det F (cid:48) x det K x · σ x · ( F (cid:48) x · K − x ) − T · K − x · ξ = ω ( ξ ) det F (cid:48) x det K x · σ x · F (cid:48)− T x · K x · K − x · ξ = ω ( ξ ) det F (cid:48) x det K x · σ x · F (cid:48)− T x · ξ . (54)19he shape tensor is only needed for J x = det F x , while the computation of the bond forcedensity does not require the shape tensor. For any particle with horizon H x and dual-horizon H (cid:48) x , the active force f xx (cid:48) and the passive force f x (cid:48) x in motion Eq. (14) or (17) canbe calculated by f xx (cid:48) = T xx (cid:48) ∀ x (cid:48) ∈ H (cid:48) x , (55)and f x (cid:48) x = T x (cid:48) x ∀ x (cid:48) ∈ H x , (56)where ξ = x (cid:48) − x , and x and x (cid:48) are the coordinate vectors for particle x and particle x (cid:48) inthe material configuration, respectively.The summary of three kinds of dual-horizon peridynamics see Table.1.original Peridynamics dual-horizon PeridynamicsModel (cid:82) H x ( f xx (cid:48) − f x (cid:48) x ) d V x (cid:48) (cid:82) H (cid:48) x f xx (cid:48) d V x (cid:48) − (cid:82) H x f x (cid:48) x d V x (cid:48) BB-PD f xx (cid:48) − f x (cid:48) x = c · s xx (cid:48) · n , ∀ x (cid:48) ∈ H x where n = η + ξ (cid:107) η + ξ (cid:107) f xx (cid:48) = C ( δ x (cid:48) ) · s xx (cid:48) · n , ∀ x (cid:48) ∈ H (cid:48) x ; f x (cid:48) x = C ( δ x ) · s xx (cid:48) · ( − n ) , ∀ x (cid:48) ∈ H x where n = η + ξ (cid:107) η + ξ (cid:107) OSB-PD f xx (cid:48) = t (cid:104) ξ (cid:105) · n , f x (cid:48) x = t (cid:48) (cid:104)− ξ (cid:105) · ( − n ) , ∀ x (cid:48) ∈ H x where n = η + ξ (cid:107) η + ξ (cid:107) f xx (cid:48) = t (cid:104) ξ (cid:105) · n , ∀ x (cid:48) ∈ H (cid:48) x ; f x (cid:48) x = t (cid:48) (cid:104)− ξ (cid:105) · ( − n ) , ∀ x (cid:48) ∈ H x where n = η + ξ (cid:107) η + ξ (cid:107) NOSB-PD f xx (cid:48) = T xx (cid:48) , f x (cid:48) x = T x (cid:48) x , ∀ x (cid:48) ∈ H x f xx (cid:48) = T xx (cid:48) , ∀ x (cid:48) ∈ H (cid:48) x ; f x (cid:48) x = T x (cid:48) x , ∀ x (cid:48) ∈ H x Table 1: Comparison of the original Peridynamics and dual-horizon Peridynamics . Numerical Examples Consider a rectangular plate with dimensions of 0.1 × (see Fig.5). The Young’smodulus, density and the Poisson’s ratio used for the plate are E = 1, ρ = 1 and ν = 0,respectively.Note that this is 1-D model. The initial displacement applied to the plate isdescribed by the following equation u ( x, y ) = 0 . − ( x .
01 ) ] , v ( x, y ) = 0 , x ∈ [0 , . , y ∈ [ − . , . , (57)where u and v denote the displacement in the x and y directions respectively. The wavespeed is v = (cid:112) E/ρ = 1 m/s. At any time step, the L error in the displacement is given by (cid:107) err (cid:107) L = (cid:107) u h − u analytic (cid:107)(cid:107) u analytic (cid:107) , (58)with (cid:107) u (cid:107) = (cid:18)(cid:90) Ω u · u d Ω (cid:19) . Three models for solving this problem, namely model A, B and C, are devised,see Figs.6(a) xy . m . m DABC coarse zonedense zone
Figure 5: setup of the plate (cid:52) x =5e-4m, yielding a critical time increment (cid:52) t max = (cid:52) x/v =5e-4s. The physical computational time step is set 0.5 second with the time increment (cid:52) t max =5e-4s.As BB-PD can not simulate the model with Poisson ratio ν = 0, the OSB-PD model isadopted. Model A was solved with the OSB-PD with constant horizons. Model B is theoriginal OSB-PD with variable horizon (without additional treatment for ghost force), andmodel C is our dual-horizon formulation for OSB-PD with variable horizon. The parametersused for the three models are listed in Table.2.Model (cid:52) x = (cid:52) y · δ Particle numbersA 0.001 3 4000B 0.001,0.0005 3/1.5 8628C 0.001,0.0005 3/1.5 8628
Table 2: Values of the Peridynamic parameters for the mixed model
Fig.7,8,9 show the displacement in the x-direction along the x-coordinate of all particles.The red star points are particles with ∆ x = 0 . x =0 . L errors of the wave profile along the line- CD (see Fig. 5). The L error (see Table. 3) of the present peridynamic formulation, can achieve the almost the same22 a) Particles distribution of model A(b) Particles distribution of model B and CFigure 6: The discretizations of model A,B and C. x d i s p l a c e m en t ( − m ) x coordinate (m) ∆ x=0.001 m Figure 7: Displacement wave u x versus x coordinate of model A at step 650 x d i s p l a c e m en t ( − m ) x coordinate (m) ∆ x=0.001 m ∆ x=0.0005 m Figure 8: Displacement wave u x versus x coordinate of model B at step 650 x d i s p l a c e m en t ( − m ) x coordinate (m) ∆ x=0.001 m ∆ x=0.0005 m Figure 9: Displacement wave u x versus x coordinate of model C at step 650 Table 3: L error of displacement for the initial condition (Gauss distribution of displacement) before andafter the wave reflection Displacement and velocity of monitor points
Two monitor points (A,B) are selected to evaluate the spurious wave reflection of ModelsB and C as shown in Fig. 5. The displacement and velocity curves show the spurious wavereflection exists in model B while not exists in model C, see Fig.10.
To study the performance of the present formulation for fracture modeling, the Kalthof-Winkler experiment is exploited. The geometry and model used for the test is depicted inFig.11, the thickness of the specimen is set as 0.01m. In the experiment, the evolution ofcrack pattern was observed to be dependent on the impact loading velocity. For a plate madeof steel 18Ni1900 subjected to an impact loading at the speed of v =32 m/s, brittle fracturewas observed [23]. The crack propagates from the end of the initial crack at an angle around70 ◦ vs. the initial crack direction. In [23], the BB-PD with constant horizon was used totest this example. For convenience of comparison, the present dual-horizon formulation isalso applied to the bond-based peridynamics to test the example. The material parametersused are the same as in [23], i.e. the elastic modulus E = 190 GPa, ρ = 7800 kg/m , ν = 0 . G = 6 . e . The impact loading was imposed by applyingan initial velocity at v = 22 m/s to the first three layers of particles in the domain as shownin Fig.11;all other boundaries are free. The plate is discretized with two different particlesizes, namely (cid:52) x coarse = 1 . (cid:52) x dense = 0 . (cid:52) x coarse =7 . a) x displacement of model C (b) y displacement of model C(c) x velocity of model C (d) y velocity of model C(e) x displacement of model B (f) y displacement of model B(g) x velocity of model B (h) y velocity of model BFigure 10: The displacement and velocity curve of monitor points A and B for model B and C
28f particles in the coarse subdomain and eight layers in the fine subdomain are employed.The total number of particles is 57968. The crack propagation speed is computed by . m . m . m . m . m . m v . m . m . m dense zone Figure 11: Kalthof-Winkler’s experimental setup V l − . = (cid:107) x l − x l − (cid:107) t l − t l − , (59)where x l and x l − are the positions of the crack tip at the times t l and t l − respectively. Anexpression for the Rayleigh wave speed c R [24] is given as c R c s = 0 .
87 + 1 . ν ν , (60)where ν is the Poisson’s ration, c s = (cid:112) µ/ρ is the shear wave speed and µ is the shearmodulus. The crack starts to propagate at 26.3 µ s. The highest crack speed reached is1530 m/s, about 54.4% of the Rayleigh speed(2799.2 m/s). The average crack speed is1077 m/s for the fine subdomain and 1094 m/s for the coarse subdomain. During the crackpropagation, the crack paths in the fine and coarse subdomains are nearly symmetrical, asshown in Figs. 12, 13 and 14. It can be seen from Fig.15 that the crack speed is very closeto that predicted by the PD formulation using constant horizons ( (cid:52) x uniform = 1 . e -3 m).The crack propagation initiated at an angle of 65.7 ◦ in the fine subdomain with respect tothe original crack and 65.8 ◦ in fine subdomain. Therefore, it can be concluded that with the29resent formulation, the transition from fine to coarse of horizons has almost no influenceon the crack propagation speed. Figure 12: The crack pattern of Kalthof-Winkler plate by the present dual horizon PD at step 350Figure 13: The crack pattern of Kalthof-Winkler simulation by the present dual horizon PD at step 650
The present dual-horizon Peridynamic formulation provides the possibility for adaptive-refinement within a simple and unified framework. Two examples of adaptively refinedperidynamics will be tested in this section: the Kalthof Winkler experiment in section 5.2and plate with pre-crack subjected to traction. The threshold for the adaptive refinementis determined by both the damage-state criteria and energy state, as proposed in [18]. Theadaptive refinement procedure consists of two steps.(1) Search for the particles that exceed the threshhold values. Note that the refinement isnot only applied to the particles above the threshold but also to the neighbouring particles.In this way, it can be ensure that the crack tip always remains inside the refined zone.30 igure 14: The crack pattern of Kalthof-Winkler simulation by the present dual horizon PD at step 875
Crack speed (m/s)
T i m e ( (cid:3)(cid:2) (cid:1)(cid:4) (cid:5) ) D e n s e s i d e C o a r s e s i d e R a y l e i g h s p e e d U n i f o r m h o r i z o n P D
Figure 15: The crack speed of Kalthof-Winkler simulation by the present dual horizon PD parent particle child particles Figure 16: Particle splitting for Adaptive refined peridynamics
We now repeat the Kalthof Winkler experiment from section 5.2. However, the initialdiscretization is based on uniform horizon size. The initial particle size is (cid:52) x initial =1.5625e-3 m. After refinement, the particle size is reduced to (cid:52) x refined =7.8125e-4 m around thecrack. The total particles is increased from 32,768 to 54,860 at the end of the simulation. Auniform refinement would result 32,768 × ◦ with respect to the initial crack direction.32 igure 17: The crack pattern of adaptive refined Kalthof-Winkler experiment at step 420Figure 18: The crack pattern of adaptive refined Kalthof-Winkler experiment at step 600Figure 19: The crack pattern of adaptive refined Kalthof-Winkler experiment at step 925 Crack speed (m/s)
T i m e ( (cid:3)(cid:2) (cid:1)(cid:4) (cid:5) ) A d a p t i v e r e f i n e d P D R a y l e i g h s p e e d U n i f o r m h o r i z o n P D
Figure 20: The crack speed of adaptive refined Kalthof-Winkler experiment
In the last example, we will test the capability of the new formulation in modelling thecrack branching. Hence, a plate with pre-crack subject to traction is considered as shownin Fig.21. The traction remains constant during the entire simulation. The same examplehas been well studied in [25] with BB-PD and also in [26] using the XFEM.The material parameters of the plate (soda-lime glass) are E = 72 GPa, ρ = 2440 kg/m and energy release rate being G = 135 J/m [25]. Two models, namely Case 1 and Case 2,are set up for comparison. Case 1 is solved by using the traditional BB-PD with constanthorizons, while Case 2 is modeled by an adaptively refined BB-PD using the present dual-horizon formulation. The plate is assumed under plane stress condition and all simulationsare carried out in 2D. The particle sizes chosen are 5 × − m for Case 1 and 1 × − m for Case 2, respectively. Case 2 uses an adaptively refined model, and once a particle isrefined, the minimal particle size will become the same as in Case 1. The particle numberare 16,000 for Case 1 and 4000-6424 for Case 2, respectively. The Rayleigh speed for thematerial is 3102 m/s according to Eq. (60). The max crack speed is 1881.4 m/s (60.6% of c R ) for Case 1 and 2184.1 m/s (70.4% of c R ) for Case 2, respectively. One possible reasonfor the different max crack speeds is the deviation of the mapping method.For Case 1, the crack starts to propagate at 11.9 µ s, and the first crack branching point34 µ s at the speed of 1043.6 m/s and the second crack branch point B µ s at the speed of 1147.1 m/s. It is observed once branching initiates, the crackpropagation speed decreases. This can be explained that the energy release to create morefracture surfaces decelerates the propagation speed. Afterward, the crack speed increases.For Case 2, the initial crack started to propagate at 10.3 µ s, first branching point B µ s with 1247.6 m/s, and second branching point B µ s with 1330.6 m/s. Thecrack pattern is similar compared with Case 1. . . �� pre-notch Figure 21: Setup of the pre-cracked plate under traction loadFigure 22: The crack pattern of pre-cracked plate with uniform horizons
6. Conclusions
This paper contributes to the development of peridynamics with varying horizons. There-fore, the interactions between particles are based on two independent horizons, passive forcefrom the horizon and the active force from the dual-horizon. The spurious wave reflectionand ghost force in the conventional peridynamics is naturally eliminated. Based on the35 igure 23: The crack pattern of adaptive refined pre-cracked plate
B 21 1 4 7 B 11 0 4 4 B 21 3 3 1 B 11 2 4 7
Crack speed (m/s)
T i m e ( (cid:3)(cid:2) (cid:1)(cid:4) (cid:5) ) U n i f o r m h o r i z o n P D A d a p t i v e l y r e f i n e d P D
Figure 24: The crack speed of pre-cracked plate under traction
Acknowledgements
The authors acknowledge the supports from FP7 Marie Curie Actions ITN-INSIST andIIF-HYDROFRAC (623667), the National Basic Research Program of China (973 Program:2011CB013800) and NSFC (51474157), the Ministry of Science and Technology of China(Grant No.SLDRCE14-B-28, SLDRCE14-B-31). The authors would like to express thegratitude to Dr. Jafar Dashlejeh and Dr. Erkan Oterkus for the beneficial discussion aboutPeridynamic theory.
ReferencesReferences [1] S.A. Silling, D.J. Littlewood, and P. Seleson. Variable horizon in a peridynamic medium. Technicalreport, Sandia National Laboratories (SNL-NM), Albuquerque, NM (United States), 2014.[2] S.A. Silling. Reformulation of elasticity theory for discontinuities and long-range forces.
Journal of theMechanics and Physics of Solids , 48(1):175–209, 2000.[3] Xin Lai, Bo Ren, Houfu Fan, Shaofan Li, CT Wu, Richard A Regueiro, and Lisheng Liu. Peridynamicssimulations of geomaterial fragmentation by impulse loads.
International Journal for Numerical andAnalytical Methods in Geomechanics , 2015.[4] W. Gerstle, N. Sau, and S. Silling. Peridynamic modeling of concrete structures.
Nuclear engineeringand design , 237(12):1250–1258, 2007.[5] E. Oterkus, E. Madenci, O. Weckner, S. Silling, P. Bogert, and A. Tessler. Combined finite elementand peridynamic analyses for predicting failure in a stiffened composite curved panel with a centralslot.
Composite Structures , 94:839–850, 2012.
6] E.T. Moyer and M.J. Miraglia. Peridynamic solutions for Timoshenko beams.
Engineering , 6:ArticleID:46262, 2014.[7] J. OGrady and J. Foster. Peridynamic plates and flat shells: A non-ordinary, state-based model.
International Journal of Solids and Structures , 51:4572–4579, 2014.[8] S.A. Silling and E. Askari. A meshfree method based on the peridynamic model of solid mechanics.
Computers & structures , 83(17):1526–1535, 2005.[9] T. Belytschko, N. Mo¨es, S. Usui, and C. Parimi. Arbitrary discontinuities in finite elements.
Interna-tional Journal for Numerical Methods in Engineering , 50(4):993–1013, 2001.[10] T. Rabczuk and T. Belytschko. Cracking particles: A simplified meshfree method for arbitrary evolvingcracks.
International Journal for Numerical Methods in Engineering , 61(13):2316–2343, 2004.[11] I. Babuˇska and J.M. Melenk. The partition of the unity finite element method. Technical ReportBN-1185, Institute for Physical Science and Technology, University of Maryland, Maryland, 1995.[12] K.W. Cheng and T.P Fries. Higher-order xfem for curved strong and weak discontinuities.
InternationalJournal for Numerical Methods in Engineering , DOI: 10.1002/nme.2768.[13] A. Gravouil, N. Moes, and T. Belytschko. Non-planar 3D crack growth by the extended finite elementand level sets - part ii: Level set update.
International Journal for Numerical Methods in Engineering ,53:2569–2586, 2002.[14] S.A. Silling and R.B. Lehoucq. Convergence of peridynamics to classical elasticity theory.
Journal ofElasticity , 93(1):13–37, 2008.[15] Stewart A Silling, M Epton, O Weckner, J Xu, and E Askari. Peridynamic states and constitutivemodeling.
Journal of Elasticity , 88(2):151–184, 2007.[16] F. Bobaru, M. Yang, L.F. Alves, S.A. Silling, E. Askari, and J. Xu. Convergence, adaptive refinement,and scaling in 1d peridynamics.
International Journal for Numerical Methods in Engineering , 77(6):852–877, 2009.[17] F. Bobaru and Y.D. Ha. Adaptive refinement and multiscale modeling in 2d peridynamics.
Journalfor Multiscale Computational Engineering , 9(635-659), 2011.[18] D. Dipasquale, M. Zaccariotto, and U. Galvanetto. Crack propagation with adaptive grid refinementin 2d peridynamics.
International Journal of Fracture , 190(1-2):1–22, 2014.[19] E. Askari, F. Bobaru, R.B. Lehoucq, M.L. Parks, S.A. Silling, and O. Weckner. Peridynamics formultiscale materials modeling. In
Journal of Physics: Conference Series , volume 125, page 012078.IOP Publishing, 2008.[20] K. Yu, X.J. Xin, and K.B. Lease. A new adaptive integration method for the peridynamic theory.
Modelling and Simulation in Materials Science and Engineering , 19(4):045003, 2011.[21] S.A. Silling and R.B. Lehoucq. Peridynamic theory of solid mechanics.
Advances in Applied Mechanics ,44:73–168, 2010.[22] T.L. Warren, S.A. Silling, A. Askari, O. Weckner, M.A. Epton, and J. Xu. A non-ordinary state-basedperidynamic method to model solid material deformation and fracture.
International Journal of Solidsand Structures , 46:1186–1195, 2008.[23] S.A. Silling. Peridynamic modeling of the Kalthoff–Winkler experiment.
Submission for the , 2001.[24] K.F. Graff.
Wave motion in elastic solids . Clarendon Oxford, 1975.[25] Y.D. Ha and F. Bobaru. Studies of dynamic crack propagation and crack branching with peridynamics.
International Journal of Fracture , 162:229–244, 2010.[26] T. Belytschko, H. Chen, J. Xu, and G. Zi. Dynamic crack propagation based on loss of hyperbolicityand a new discontinuous enrichment.
International Journal for Numerical Methods in Engineering ,58(12):1873–1905, 2003.,58(12):1873–1905, 2003.