Dual-horizon peridynamics: A stable solution to varying horizons
ScienceDirect
Dual-horizon peridynamics: A stable solution to varying horizons
Huilong Ren c , Xiaoying Zhuang d,e , Timon Rabczuk a,b,c, ∗ a Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam b Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City, Viet Nam c Institute of Structural Mechanics, Bauhaus-Universit¨at Weimar, 99423 Weimar, Germany d State Key Laboratory of Disaster Reduction in Civil Engineering, College of Civil Engineering, Tongji University, Shanghai 200092, China e Institute of Continuum Mechanics, Leibniz Universit¨at Hannover, Hannover, Germany
Received 11 April 2016; received in revised form 18 November 2016; accepted 19 December 2016Available online 7 February 2017
Highlights • Stable crack paths with variable horizon and particle sizes. • Dual-horizon peridynamics for multiple materials/material interfaces. • Fulfillment of balance of momentum and balance of angular momentum.
Abstract
In this paper, we present a dual-horizon peridynamics formulation which allows for simulations with dual-horizon with minimalspurious wave reflection. We prove the general dual property for dual-horizon peridynamics, based on which the balance ofmomentum and angular momentum in PD are naturally satisfied. We also analyze the crack pattern of random point distributionand the multiple materials issue in peridynamics. For selected benchmark problems, we show that DH-PD is less sensitive to thespatial than the original PD formulation.c ⃝ Keywords:
Horizon variable; Random point distribution; Crack stability; Multiple materials; Material interfaces
1. Introduction
Peridynamics (PD) has recently attracted broad interests for researchers in computational solid mechanics sinceit shows great flexibility in modeling dynamic fractures. The crack happens naturally, and no attention is paid todescribe crack topology in the PD simulation. Hence, complicated fracture patterns i.e. crack branching and coales-cence of multiple cracks, can be captured with ease by breaking the bonds between material points. No techniquessuch as smoothing the normals of the crack surfaces in order to avoid erratic crack paths is necessary in PD, whilesuch techniques are commonly used in the extended finite element method (XFEM) [1], meshless methods [2] or ∗ Corresponding author at: Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City, Viet Nam.
E-mail addresses: [email protected] (X. Zhuang), [email protected] (T. Rabczuk).http://dx.doi.org/10.1016/j.cma.2016.12.0310045-7825/ c ⃝ a r X i v : . [ phy s i c s . c o m p - ph ] M a r . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 other partition of unity methods (PUM) [3]. Meanwhile, PD can be easily extended from 2D to 3D, which cannotalways be done in other methods [4–6]. In addition, an important feature in PD is that it can now be easily coupledwith the more classical local continuum mechanics [7,8]. The traditional PD method was proposed by Silling [9]and has been applied onwards for a wide range of problems including impact loading, fragmentation [10,11],composites delamination [12], beam and plate/shell structures [13–15], thermal diffusion [16] and flow in porousmedia [17,18].In PD, the equation of motion is formulated in an integral form rather than the partial differential form. Thematerial point interacts with other material points in its horizon, a domain associated to the point. The horizonis often represented by circular area (2D) or spherical region (3D) with a radius. There are mainly 3 types of PDformulations, i.e. bond based peridynamics (BB-PD), ordinary state based peridynamics (OSB-PD) and non-ordinarystate based peridynamics (NOSB-PD). PD formulation starts with discretizing the solid domain into material points,where any two material points within horizon’s radius form a bond. The first version of PD is called the bond basedformulation (BB-PD) [9] as the bonds act like independent springs. One limitation of BB-PD is that only the materialwith Poisson’s ratio of 1/3 in 2D or 1/4 in 3D can be modeled. In order to eliminate the Poisson’s ratio restriction,the state based peridynamics (SB-PD) [19] is proposed, where “state” means the bond deformation depends not onlyon the deformation of the bond itself, but also on collective deformation of other bonds. The SB-PD was extendedinto two types, namely the OSB-PD and the NOSB-PD. After that, BB-PD is regarded as a special case of OSB-PD.NOSB-PD is capable of simulating the material with advanced constitute models, e.g. viscoplasticity model [20], theDrucker–Prager plasticity model [10] and Johnson–Cook (JC) constitutive model [21].Aside from the pure PD formulation above mentioned, much work was devoted to the coupling of peridynamicswith the continuum methods [7,8,22–24]. Such strategies enable to model the fracture domain with PD while otherdomain modeled with continuum methods. Seleson et al. developed a force-based coupling scheme for peridynamicsand classical elasticity in one dimension [8] and higher dimensions [22]. This method satisfies the Newton’s third lawand Patch test, and no ghost force exists in the blending domain. Lubineau et al. developed a morphing strategy tocouple non-local to local continuum mechanics. The strategy is based on the domain transition where the strain energyis divided into nonlocal part and local continuum part. In this method, the PD is discretized with finite element andthe local integration matrix is assembled.One basic constraint of traditional PD formulations is that the radii of horizons are required to be constant,otherwise spurious wave reflections shall emerge and ghost forces between material points will deteriorate the results.However, in many applications, for sake of computational efficiency, it is necessary to vary the horizon sizes withrespect to the spatial distribution of the material points, e.g. for adaptive refinement, and multiscale modeling. Inthe implementation of traditional PD discretized with particles, in order to achieve acceptable accuracy, the horizonradius has to be determined with respect to the lowest material point resolution locally required. It should be notedthat horizon is guided by the characteristic range of long-range forces [9,7], not the discretization of the computationaldomain.The dual-horizon peridynamics (DH-PD) is firstly proposed to solve the issue of spurious wave reflections whenvariable horizons are adopted. The key idea lies in the definition of dual-horizon, which is defined as the dual termof the horizons centering at each material points. Within the framework of DH-PD, the traditional peridynamics withconstant horizon can be derived. The derivation of DH-PD is based on the Newton’s third law and other additionaltechniques such as variational principal or Taylor expansion [25] are not required. The dual-horizon PD not onlyconsiders the variable horizons, but also reduce the computational cost, allowing the domain discretized like finiteelement method(FEM) or finite volume method (FVM), where the domain of interest utilizes dense node distributionand the other domain coarse node distribution.In this paper, we review the DH-PD and propose a universal principle in dual-horizon, the dual property of dual-horizon. We also discuss the influence of random material points distribution on crack pattern and multiple materialsissue in peridynamics. The content of the paper is outlined as follows. Section 2 begins with the equation of motionof dual-horizon peridynamics and then discuss three ways to calculate the bond force in DH-PD. Section 3 proves thegeneral dual property of dual-horizon, based on which the balance of linear momentum and the balance of angularmomentum of the present PD formulation are re-proved. In Section 4, the issue of multiple materials in peridynamicsis discussed with 2 numerical examples. In Section 5, the simulation of Kalthoff–Winkler experiment is carried out totest the influence of random material points distribution on the crack pattern in peridynamics. H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
Fig. 1. The configuration for deformed body.Fig. 2. Force vectors of BB-PD,OSB-PD and NOSB-PD.
2. Dual-horizon peridynamics
Consider a solid domain in the initial and current configuration as shown in Fig. 1. x is the material coordinate ininitial configuration Ω ; and y := y ( x , t ) and y ′ := y ( x ′ , t ) are the spatial coordinates in the current configuration Ω t ; u := u ( x , t ) and u ′ := u ( x ′ , t ) are the displacement vectors for associated points; v := v ( x , t ) and v ′ := v ( x ′ , t ) arethe velocity for associated points; ξ := x ′ − x denotes the initial bond vector; the relative displacement vector of bond ξ is η := u ′ − u , and the current bond vector is y ⟨ ξ ⟩ := y ( x ′ , t ) − y ( x , t ) = ξ + η .Since the bond ξ starts from x to x ′ , an alias xx ′ is used to denote ξ . Obviously, x ′ x is in opposite direction of xx ′ .Let f xx ′ := f xx ′ ( η , ξ ) denote the force vector density acting on point x due to bond xx ′ ; based on Newton’s third law, x ′ undertakes a reaction force − f xx ′ .As the conventional peridynamics formulation with constant horizon can be viewed as a special case of the DH-PD [26], we here mainly focus on the DH-PD. There are general three types of peridynamics, as shown in Fig. 2 andDH-PD will be applied to these three types. In conventional peridynamics, the horizon is required to be constant for a homogeneous body, where beinghomogeneous means that the strain energy density for any uniform deformation must be invariant with respect tochanges in δ [27]. When the computational domain is discretized with different material point sizes, as shown inFig. 3, in order to avoid the spurious wave reflection, the number of neighbors in the right side is dramaticallyincreased. The fact of too many neighbors not only increases the computational cost, but also contributes little tothe accuracy improvement, as the value in the point is smoothed out by its neighbors.When different horizon radii are used, spurious wave reflections occur in conventional peridynamics. The reasonis that the traditional peridynamics formulation with single horizon is not easy to account for the interaction betweentwo points with different horizon radii, which has been discussed in detail in [26]. In conventional PD, the horizon is a domain centered at each point with a radius of δ . In DH-PD formulation, the horizon H x is the domain where any material point x ′ falling inside forms bond xx ′ , where the bond xx ′ will exert . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 i , j and k with constant horizon radius δ = r i . r i = r j = r k , neighbor number N i ≈ , N j ≈ N i ≈ , N k ≈ N i ≈ { x , x , x , x } ∈ H ′ x , whose horizonsare denoted by thin solid line; the red points { x , x } ̸∈ H ′ x , whose horizons are denoted by dashed line. (For interpretation of the references tocolor in this figure legend, the reader is referred to the web version of this article.) direct force f xx ′ on x and whereas x ′ will undertake reaction force − f xx ′ . In this sense, the horizon in DH-PD can beviewed as direct force horizon. The dual-horizon in DH-PD is defined as a union of points whose horizons include x ,denoted by H ′ x = { x ′ | x ∈ H x ′ } . (1)The bond from the dual-horizon is named as dual-bond. For any point x ′ ∈ H ′ x , the dual-bond x ′ x is actually the bond x ′ x from H x ′ . Dual-bond x ′ x exerts reaction force on x . In this sense, the dual-horizon can be viewed as reaction forcehorizon. One simple example to illustrate the dual-horizon is shown in Fig. 4, where { x , x } ̸∈ H ′ x because H x and H x do not contain x . The direct force and reaction force from horizon and dual-horizon are shown in Fig. 5.Bond and dual-bond are broken independently in the fracture models. There are reasons accounting for it. Forexample, in BB-PD, the equivalence of energy release rate of homogeneous body in fracture models indicates thatdifferent horizon sizes correspond to different critical stretches. This means bond xx ′ and dual-bond x ′ x may havedifferent critical stretches. If they are broken at the same bond stretch, the one with larger critical bond stretch isbroken in advance, which indicates a lower material strength for that point and thus violates the equivalence of energyrelease rate. The internal forces that are exerted at each material point include two parts, the direct forces from the horizon andthe reaction forces from the dual-horizon. The other forces applied to a material point include the body force and theinertia force. Let ∆ V x denote the volume associated with x . The body force for x can be expressed as b ( x , t ) ∆ V x ,where b ( x , t ) is the body force density. The inertia force is denoted by ρ ¨ u ( x , t ) ∆ V x , where ρ is the density. H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
Fig. 5. The force vectors in dual-horizon peridynamics. Bond xx ′ exerts direct force f xx ′ on x and reaction force − f xx ′ on x ′ . In this case, f xx ′ ̸= , f x ′ x = The direct force due to bond xx ′ in H x is f xx ′ ( η , ξ ) · ∆ V x · ∆ V x ′ , ∀ x ′ ∈ H x . (2)The reaction force due to dual-bond x ′ x in H ′ x is − f x ′ x ( − η , − ξ ) · ∆ V x ′ · ∆ V x , ∀ x ′ ∈ H ′ x . (3)By summing over all forces on point x , including inertia force, body force, direct force in Eq. (2) and reaction forcein Eq. (3), the equation of motion for DH-PD is obtained. ρ ¨ u ( x , t ) ∆ V x = H x f xx ′ ( η , ξ ) ∆ V x ′ ∆ V x + H ′ x ( − f x ′ x ( − η , − ξ ) ∆ V x ′ ∆ V x ) + b ( x , t ) ∆ V x . (4)As the volume ∆ V x associated to x is independent of the summation, we can eliminate ∆ V x in Eq. (4), yielding thegoverning equation based on x : ρ ¨ u ( x , t ) = H x f xx ′ ( η , ξ ) ∆ V x ′ − H ′ x f x ′ x ( − η , − ξ ) ∆ V x ′ + b ( x , t ). (5)When the discretization is sufficiently fine, the discrete form converges to the integral form:lim ∆ V x ′ → H x f xx ′ ( η , ξ ) ∆ V x ′ = H x f xx ′ ( η , ξ ) d V x ′ (6)and lim ∆ V x ′ → H ′ x f x ′ x ( − η , − ξ ) ∆ V x ′ = H ′ x f x ′ x ( − η , − ξ ) d V x ′ . (7)Thus, the integral form of the equation of motion in dual-horizon peridynamics (DH-PD) is given as ρ ¨ u ( x , t ) = H x f xx ′ ( η , ξ ) d V x ′ − H ′ x f x ′ x ( − η , − ξ ) d V x ′ + b ( x , t ). (8)Note that Eq. (8) is also valid for points near the boundary area.When the horizons are set constant, as x ′ ∈ H ′ x ⇔ x ′ ∈ H x , we have H ′ x = H x . Then Eq. (8) degenerates totraditional constant horizon peridynamics. ρ ¨ u ( x , t ) = H x f xx ′ ( η , ξ ) − f x ′ x ( − η , − ξ ) d V x ′ + b ( x , t ). (9) We aim at solving Eq. (5), the discrete form of DH-PD with the Velocity-Verlet algorithm [28]. After discretizationof the domain, each material point’s neighbors based on its horizon radius are stored explicitly, while the neighbors indual-horizon are inferred from other material points’ horizons. Every time the bond force from one point’s horizon is . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 calculated, the dual-bond force from the other point’s dual-horizon is determined according to Newton’s third law. Forexample, assuming the material points with horizons are shown in Fig. 4, the horizons of x is H x = { x , x , x , x } ,and H ′ x = { x , x , x , x } . Direct forces on x from its horizon H x are summed as follows: x ∈ H x , calculate f xx , add f xx ∆ V x ∆ V x to x , add − f xx ∆ V x ∆ V x to x ; x ∈ H x , calculate f xx , add f xx ∆ V x ∆ V x to x , add − f xx ∆ V x ∆ V x to x ; x ∈ H x , calculate f xx , add f xx ∆ V x ∆ V x to x , add − f xx ∆ V x ∆ V x to x ; x ∈ H x , calculate f xx , add f xx ∆ V x ∆ V x to x , add − f xx ∆ V x ∆ V x to x .We do not calculate f x x , f x x , f x x , f x x from dual-horizon H ′ x . The force from H ′ x are determined when calculatingforces in H x , H x , H x , H x . x ∈ H x , calculate f x x , add f x x ∆ V x ∆ V x to x , add − f x x ∆ V x ∆ V x to x ; x ∈ H x , calculate f x x , add f x x ∆ V x ∆ V x to x , add − f x x ∆ V x ∆ V x to x ; x ∈ H x , calculate f x x , add f x x ∆ V x ∆ V x to x , add − f x x ∆ V x ∆ V x to x ; x ∈ H x , calculate f x x , add f x x ∆ V x ∆ V x to x , add − f x x ∆ V x ∆ V x to x .Hence, in the process of calculating forces in horizons, the forces from dual-horizons are automatically done.Compared with the algorithm for OSB-PD given in [28], force expression: f xx ′ − f x ′ x for point x with bond xx ′ and force expression: f x ′ x − f xx ′ for point x ′ with bond x ′ x are calculated independently. Compared with the forcesummation in traditional PD, the bond force for every bond in DH-PD is calculated only once. f xx ′ There mainly exist three types of peridynamics formulations [9,19], namely BB-PD, OSB-PD and NOSB-PD.Their key difference is the way they calculate the bond force density f x ′ x . The applications of DH-PD to the threetypes of peridynamics are direct [26]. The calibration of constitutive parameters with respect to the continuum modeland some issues concerning the implementations will be discussed, for example, in BB-PD.In BB-PD, the energy per unit volume in the body at a given time t is given by [29] W = H x w( η , ξ ) d V ξ , (10)where w( η , ξ ) is the strain energy for bond ξ . For the BB-PD theory, by enforcing the strain energy density beingequal to the strain energy density in the classical theory of elasticity [29], we obtain the expression of the microelasticmodulus C (δ) as C (δ) = E π δ ( − ν) , (11)where δ is the horizon radius associated to that material point, and ν is the Poisson’s ratio. Note that the value of C (δ) takes half of the microelastic modulus c used in the constant-horizon(CH) BB-PD since the bond energy for varyinghorizon is determined by both horizon and dual-horizon.Let w ( ξ ) = C (δ) s (δ)ξ / s (δ) is the critical bond stretch.By breaking half of all the bonds connected to a given material point along the fracture surface and equalizing thebreaking bonds energy with the critical energy release rate G [11], we can get the expression between the criticalenergy release rate G and the critical bond stretch s (δ) in three dimensions s (δ) = G E δ . (12)Both the microelastic modulus C (δ) and the critical stretch s (δ) depend on the horizon radii for variable horizons.In the implementation of the bond based peridynamics, fracture is introduced by removing points from the neighborlist once the bond stretch exceeds the critical bond stretch s . In order to specify whether a bond is broken or not, ahistory-dependent scalar valued function µ is introduced [29], µ( t , ξ ) = s ( t ′ , ξ ) < s for all 0 ≤ t ′ ≤ t , . (13) H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
The local damage at x is defined as φ( x , t ) = − H x µ( x , t , ξ ) d V ξ H x d V ξ . (14)For any material point with dual-horizon ( H ′ x ) and horizon ( H x ), the direct force f xx ′ and the reaction force f x ′ x inEq. (5) or (8) are computed by the following expressions, respectively f xx ′ = C (δ x ′ ) · s xx ′ · η + ξ ∥ η + ξ ∥ , ∀ x ′ ∈ H x (15) f x ′ x = C (δ x ) · s xx ′ · − ( η + ξ ) ∥ η + ξ ∥ , ∀ x ′ ∈ H ′ x , (16)where C (δ x ′ ) and C (δ x ) are the microelastic modulus based on δ x ′ and δ x computed from Eq. (11) respectively, and s xx ′ is the stretch between points x and x ′ .For the other two types of peridynamics [19], comparing the equation of motion in traditional state-based PD withEq. (9), the expression for bond force f xx ′ and f xx ′ can be obtained with ease. For more details please refer to [26].
3. The dual property of dual-horizon
Let F ( x , x ′ ) be any expression depending on bond xx ′ . The dual property of dual-horizon is that the double inte-gration of the term F ( x , x ′ ) in dual-horizon can be converted to the double integration of the term F ( x ′ , x ) in horizon,as shown in Eqs. (17) and (18). The key idea lies in that the term F ( x , x ′ ) can be both interpreted in H x and H ′ x ′ . Ω H ′ x F ( x , x ′ ) ∆ V x ′ ∆ V x = Ω H x F ( x ′ , x ) ∆ V x ′ ∆ V x Discrete form (17) Ω H ′ x F ( x , x ′ ) d V x ′ d V x = Ω H x F ( x ′ , x ) d V x ′ d V x Continuum form . (18) Proof.
Let Ω be discretized with N voronoi tessellations (or other shape), as shown in Fig. 6. Each polygon is denotedwith an index i ∈ { , . . . , N } , x i is the coordinate for i ’s center of gravity, ∆ V i is the volume associated to polygon i , H i and H ′ i are i ’s horizon and dual-horizon, respectively. So Ω = ≤ i ≤ N ∆ V i . Let F ( i , j ) := F ( x i , x j ) be any expression depending on bond x i x j . Consider the double summation of F ( i , j ) on Ω . ∆ V i ∈ Ω ∆ V j ∈ H ′ i F ( i , j ) ∆ V j ∆ V i = ≤ i ≤ N j ∈ H ′ i F ( i , j ) ∆ V j ∆ V i = j ∈ H ′ F ( , j ) ∆ V j ∆ V + j ∈ H ′ F ( , j ) ∆ V j ∆ V + · · · + j ∈ H ′ N F ( N , j ) ∆ V j ∆ V N . (19)In the third step, j ∈ { , . . . , N } means j belongs to that point’s dual-horizon. Each term F ( i , j ) ∆ V j ∆ V i in H ′ i canbe interpreted as term F ( i , j ) ∆ V i ∆ V j in H j . Let us sum all terms in a way based on point j ’s horizon H j , where j ∈ { , . . . , N } j ∈ H ′ F ( , j ) ∆ V j ∆ V + j ∈ H ′ F ( , j ) ∆ V j ∆ V + · · · + j ∈ H ′ N F ( N , j ) ∆ V j ∆ V N = i ∈ H F ( i , ) ∆ V ∆ V i + i ∈ H F ( i , ) ∆ V ∆ V i + · · · + i ∈ H N F ( i , N ) ∆ V N ∆ V i . (20) . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 Ω . In the second step of Eq. (20), for example, i ∈ H F ( i , ) ∆ V ∆ V i means gathering all terms j = j ∈ H ′ F ( , j ) ∆ V j ∆ V + j ∈ H ′ F ( , j ) ∆ V j ∆ V + · · · + j ∈ H ′ N F ( N , j ) ∆ V j ∆ V N i ∈ H F ( i , ) ∆ V ∆ V i + i ∈ H F ( i , ) ∆ V ∆ V i + · · · + i ∈ H N F ( i , N ) ∆ V N ∆ V i = j ∈ H F ( j , ) ∆ V j ∆ V + j ∈ H F ( j , ) ∆ V j ∆ V + · · · + j ∈ H N F ( j , N ) ∆ V j ∆ V N = ≤ i ≤ N j ∈ H i F ( j , i ) ∆ V j ∆ V i . (21)In the second step of Eq. (21), i and j are swapped, yielding ≤ i ≤ N j ∈ H ′ i F ( i , j ) ∆ V j ∆ V i = ≤ i ≤ N j ∈ H i F ( j , i ) ∆ V j ∆ V i . When i and j are replaced with x and x ′ , respectively, we have Ω H ′ x F ( x , x ′ ) ∆ V x ′ ∆ V x = Ω H x F ( x ′ , x ) ∆ V x ′ ∆ V x . (22)When N → ∞ so that ∆ V i →
0, we havelim N →∞ ≤ i ≤ N j ∈ H ′ i F ( i , j ) ∆ V j ∆ V i = i ∈ Ω j ∈ H ′ i F ( i , j ) d V j d V i . Hence, the dual property of dual-horizon in the integral form is Ω H ′ x F ( x , x ′ ) d V x ′ d V x = Ω H x F ( x ′ , x ) d V x ′ d V x . (23)Eq. (23) means the double integration of the term in dual-horizon can be converted to the double integration of theterm with x and x ′ swapped in horizon.When F ( x , x ′ ) = f xx ′ , we have Ω H x f xx ′ d V x ′ d V x = Ω H ′ x f x ′ x d V x ′ d V x . (24) H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
When F ( x , x ′ ) = ( x + u ) × f xx ′ , we have Ω H x ( x + u ) × f xx ′ d V x ′ d V x = Ω H ′ x ( x ′ + u ′ ) × f x ′ x d V x ′ d V x . (25) Two basic physical principles include the balance of the linear momentum and the balance of the angularmomentum. We showed in [26] that DH-PD fulfills these basic principles. In this manuscript, we provide a shorterproof based on the dual property of the dual-horizon.
Balance of linear momentum
The internal forces shall satisfy the balance of linear momentum for any bounded body Ω given in form by Ω (ρ ¨ u ( x , t ) − b ( x , t )) d V x = Ω H ′ x f xx ′ ( η , ξ ) d V x ′ d V x − Ω H x f x ′ x ( − η , − ξ ) d V x ′ d V x = . (26) Proof.
Based on Eq. (24), Eq. (26) is apparently satisfied, so the linear momentum is conserved. (cid:3)
Balance of angular momentum
To satisfy the balance of angular momentum for any bounded body Ω , it is required that Ω y × (ρ ¨ u ( x , t ) − b ( x , t )) d V x = Ω y × H ′ x f xx ′ ( η , ξ ) d V x ′ − H x f x ′ x ( − η , − ξ ) d V x ′ d V x = Ω H ′ x y × f xx ′ ( η , ξ ) d V x ′ d V x − Ω H x y × f x ′ x ( − η , − ξ ) d V x ′ d V x = . (27)For simplicities, let f xx ′ represent f xx ′ ( η , ξ ) , and f x ′ x represent f x ′ x ( − η , − ξ ) . Proposition.
In the dual-horizon peridynamics, suppose a constitutive model of the form f = ˆ f ( y , Λ ) (28) where ˆ f : V → V is bounded and Riemann-integrable on Hilbert function space and V is the vector state; Λ denotesall variables other than the current deformation vector state.If H x y ⟨ x ′ − x ⟩ × f x ′ x d V x ′ = ∀ y ∈ V , (29) then the balance of angular momentum, Eq. (27) holds for any deformation of Ω for any given constitutive model. Proof.
Based on Eq. (25), Eq. (27) can be written as Ω y × (ρ ¨ u ( x , t ) − b ( x , t )) d V x = Ω y × H ′ x f xx ′ ( η , ξ ) d V x ′ − H x f x ′ x ( − η , − ξ ) d V x ′ d V x . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 = Ω H ′ x y × f xx ′ ( η , ξ ) d V x ′ d V x − Ω H x y × f x ′ x ( − η , − ξ ) d V x ′ d V x = Ω H x y ′ × f x ′ x ( − η , − ξ ) d V x ′ d V x − Ω H x y × f x ′ x ( − η , − ξ ) d V x ′ d V x = Ω H x ( y ′ − y ) × f x ′ x ( − η , − ξ ) d V x ′ d V x = . (30)The expression in step 4 is equal to zero in all three types of peridynamics (see proof in [19]), so the angular momentumis conserved. (cid:3) It can be seen that the conservation of angular momentum somehow depends only on the horizon; the dual-horizonis not involved. The latter is only needed in Eqs. (5) and (8). This conclusion can be also illustrated for the BB-PDand OS-PD since the internal forces f xx ′ and f x ′ x are parallel to the bond vector in the current configuration, see Fig. 2.
4. Wave propagation in 1D homogeneous bar
Consider a 1D bar with length of L = ∆ x = . × − m and the other parts discretized with ∆ x = ∆ x . The horizon size is selected as δ i = . ∆ x i , as shown inFig. 7. Both ends of the bar are free. The initial displacement field is given by u ( x , ) = .
002 exp − ( x − . L ) . (31)The bar is solved with DH-PD (Case I) and traditional PD formulation (Case II). The wave profiles at different timesare shown in Fig. 8 for DH-PD and Fig. 9 for traditional PD. It can be seen the horizon size interfaces have limitedinfluence on the wave profiles in Case I, while the spurious wave emerged for traditional PD formulation. Two points x = .
32 m and x = .
34 m are selected to analyze the wave profile, and the displacements at different time are shownin Fig. 10. It shows that the incident wave passed the interface and the magnitude of reflected wave is smaller than3.17% of the magnitude of incident wave.On the other hand, for the traditional PD without any treatment on variable horizons, when the incident waveapproached the interface, the deformation of the fine material points caused the ghost forces among the coarse materialpoints. The ghost forces give rise to the spurious wave in the domain with coarse material points, as shown in Fig. 9( T = . H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
Fig. 8. The wave profiles for Case I at different times.Fig. 9. The wave profiles for Case II at different times. . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
5. Multiple materials
The material interface problem in peridynamics was formally discussed in [30], where the authors analyzed theconvergence of peridynamics for heterogeneous media, and found that the operator of linear peridynamics divergesin the presence of material interfaces. Some approaches have been proposed to deal with this issue. In Ref. [19], theinterfaces and mixtures between multiple materials are accounted for by using separate influence functions for eachmaterials. For example, if x is near an interface between two materials denoted by (1) and (2), two separate influencefunctions ω ( ) and ω ( ) can be defined such that each vanishes at points x + ξ outside its associated material. Theforce state is given by T xx ′ = ω ( ) ⟨ ξ ⟩ P ( ) x K − ξ + ω ( ) ⟨ ξ ⟩ P ( ) x K − ξ (32)where P ( ) x and P ( ) x are the first Piola–Kirchhoff stress from the respective constitution models and respectivedeformation gradient tensors F ( ) and F ( ) , K is the shape tensor. In Ref. [30], a peridynamic interface model wasproposed. However, such issue can be handled in a simpler way in DH-PD. In the DH-PD formulation, any point x has a force state T x (cid:3) , which acts like a force potential field applied to its nearby neighbors in H x . We assume x imposesreaction bond forces to its neighbors no matter what the neighbor’s material type is. The bond force can be determinedby one side or both sides. For examples, in the BB-PD, microelastic modulus C (δ) and critical material stretch s (δ) are the primary parameters. For two points i and j with initial coefficient [ C (δ i ), s (δ i ) ] (abbreviated as [ C i , s i ] ) and [ C (δ j ), s (δ j ) ] , respectively, when taking into consideration of the combination of the initial coefficients, the alternateforms can be formulated as the arithmetic average or geometric average of microelastic moduli or critical stretches.For the OSB-PD and NOSB-PD in multiple materials, the force state t xx ′ and T xx ′ can be calculated similarly to thesingle homogeneous material and hence, it is not necessary to distinguish between homogeneous and heterogeneousmaterials. The length of the bar is 1 m, discretized with grid spacing ∆ x = . , .
005 m. The elastic modulus and density are E = , E = ρ =
1, respectively. Boundary conditions are applied at 3 material points at each end, i.e. fixedboundary condition at left end (blue points in Fig. 11) and displacement boundary conditions of ∆ s = − α = − t = . f i = δ i j ∈ H i , j > i C i s i j ∆ V j (33)where s i j is the bond stretch for bond i j , δ i is the radius of horizon H i . Due to the error of discretization, Eq. (33)cannot sum the sectional force from material point i accurately; when right half of H i contains just 3 neighbors andvolume correction is considered, there are 2.5 material points in Eq. (33), which corresponds to 5/6 of the theoreticalsectional force. For Case I, the force discontinuity in the interface of different point spacing is observed in Fig. 12.The force solution diverts from the analytical solution. The reason for the result is that the microelastic modulus isdetermined by one point. Consider a point i close to the interface of two materials with microelastic modulus C i , theequilibrium of i requires the total force from left side being equal to that from right side. As a result, the displacementfield in left half and right half horizon are anti-symmetrical with respect to i . Hence, there is no discontinuity in thestrain field at the interface. At last, the spurious equilibrium displacement and force field are obtained.However, for Case II modeled with DH-PD, the force solution and displacement solution agree well with theanalytical solutions, respectively, see Fig. 13. Both the material interface and point spacing interface have minimalinfluence on the final solution, though small variation is observed at the interfaces of material or grid spacing. Thevariation of the sectional force arises from the nonlocality in peridynamics. In PD, the force can transmit over finitedistance, the equilibrium of a point means the forces from horizon and dual-horizon are canceled. In this sense, theequilibrium is determined by the domains of horizon and dual-horizon, two nearby material points may have differentdomains, therefore, the jumping of internal sectional force is reasonable. This is essentially different from the ghostforce which is unbalanced. H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
Fig. 10. The wave for points x = .
32 m and x = .
34 m.Fig. 11. The setup of 1D bar.Fig. 12. The displacement and sectional force for Case 1.
In conclusion, we show for a simple 1D example that DH-PD has the potential to effectively deal with multiplematerials.
Consider a plate subjected to constant velocity traction at both the top and down edges, as shown in Fig. 14. Thematrix is filled with many randomly distributed inclusions; all inclusions have the same material parameters. Thematerial parameters for matrix include E =
72 GPa, ν = / G =
40 J/m . The parameters for the inclusionsare E =
144 GPa, ν = / G =
80 J/m . The plate is discretized into 40,000 material points with point spacingof ∆ x = . − δ = . ∆ x . The material’s critical bond stretch is calculated basedon Eq. (12), i.e. s = s ≈ −
3. Three cases are conducted to show the capability of DH-PD in dealing withcomposite materials. The velocity of traction on boundary is v bc = . / s for Case I, II and v bc = / s for Case III.Two different heterogeneous structures are tested; the heterogeneous structures for Case II and III coincide. . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 ji depends on i . When calculating i , all neighbors in H i are pretended to have the same parametersas i . Bond i j depends on j , when calculating j , all neighbors in H j are pretended to have the same parameters as j . Bond i j and ji are brokenseparately. We use a simple damage rule to take into account the crack propagation in heterogeneous materials. In theframework of DH-PD, the bond and dual-bond are broken separately. For example, as shown in Fig. 15, bond ji is interpreted as i exerting force on j , and the strength of bond ji is assumed to depend only on i parameters, e.g. i ’scritical stretch; the same applies to bond i j . Figs. 16(b) and 16(c) show when initial crack tip is in the inclusion(Fig. 16(c)), the crack will firstly propagate in the inclusions; when initial crack tip is in the matrix, the crack firstlypropagates in the matrix. We also note two phenomena when the load rate is increased:(1) The number of crack branches is increased. Such a phenomenon is well known in dynamic fracture.(2) The crack propagates through the inclusions. This tendency has been also observed experimentally, i.e. the numberof aggregates/particle cracks in a weaker matrix material is increased with increasing strain rate/load rate. H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 (a) Initial crack. (b) Case I : v bc = . / s.(c) Case II : v bc = . / s. (d) Case III : v bc = / s.Fig. 16. The crack patterns in multiple-material plate.
6. Simulation of the Kalthoff–Winkler experiment
The Kalthoff–Winkler experiment is a classical benchmark problem studied by [31–37]. The setup is depicted inFig. 19; the thickness of the specimen is 0.01 m. In the experiment, the evolution of crack pattern was observed tobe dependent on the impact loading velocity. For a plate made of steel 18Ni1900 subjected to an impact loadingat the speed of v =
32 m/s, brittle fracture was observed [38]. The crack propagates from the end of the initialcrack at an angle around 70 ◦ with respect to the initial crack direction. In [38], the BB-PD with constant horizon wasused to test this example. For convenience of comparison, the present dual-horizon formulation is also applied to thebond based peridynamics to test the example. The material parameters used are the same as in [38], i.e. the elasticmodulus E =
190 GPa, ρ = , ν = .
25 and the energy release rate G = . × J/m . The impactloading was imposed by applying an initial velocity at v =
22 m/s to the first three layers of material points in thedomain as shown in Fig. 19; all other boundaries are free. For brevity, the surface correction method [39] or fictitiousboundary [40] are not employed.In order to create the irregular material points distribution in the computational domain, we utilize the Abaqussoftware to mesh the geometrical model with different mesh seeds setting, and then convert the finite element meshinto material points. For example, in the case of triangle element or tetrahedron element, each element’s area(orvolume) is evenly allocated to its nodes, as shown in Figs. 17 and 18. Every node (acted as material point) is associated . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 with a unique volume or mass. The material point size ∆ x is the equivalent diameter, which is calculated based on itsvolume by assuming the material point is circle or sphere. Six cases as shown in Table 1, are tested. Case I, II are intwo dimensions with the same material points distribution where the ratio ∆ x max / ∆ x min = .
1, while Case III–VIin three dimensions with particle size ratio ∆ x max / ∆ x min = .
4. Case I and III are based on the traditional constant-horizon peridynamics (CH-PD), while Case II and IV are carried out with the DH-PD. Case V is conducted to testthe influence of spurious wave reflections when variable horizons is used in traditional PD. Case VI is used to testthe effect of volume correction given in Appendix A. Two different discretizations (for I, II and III–VI, respectively)with a bias towards a ‘wrong’ crack path. In Case II and IV–VI, the radius of each material point’s horizon is selectedas three times of the material point size, e.g. δ i = . ∆ x i , while in Case I and III δ = . ∆ x max . The bondmicroelastic modulus C (δ x ) is calculated based on Eq. (11).The angle of the crack propagation with respect to the original crack direction is shown in Table 2. For relativelylow speed impact (e.g. ≤
30 m/s), the crack angle in [41,42] is approximately 70 ◦ . Table 2 shows that result byDH-PD is better than that by CH-PD. The results in Fig. 20 and Table 2 show that the crack pattern both in 2D and3D cases by the traditional constant-horizon bond based peridynamics differs from the experiment, while the crackpattern obtained by the DH-PD agrees well with the experiment besides the introduction of a mesh bias towards an‘incorrect’ crack path. The comparison of Case IV and Case V shows the crack pattern is affected by the spuriouswave reflections. The comparison of Case IV and Case VI shows the volume correction has some positive effect onthe crack propagation, therefore, the volume correction is recommended. H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 (a) Case I. (b) Case II.(c) Case III. (d) Case IV.(e) Case V. (f) Case VI.Fig. 20. The crack patterns in the simulations of Kalthoff–Winkler experiment.Table 1The parameters in 4 simulations of Kalthoff–Winkler experiment.Case δ max /δ min Type Volume correctionI 1 CH-PD YesII 3.1 DH-PD YesIII 1 CH-PD YesIV 4.4 DH-PD YesV 4.4 CH-PD YesVI 4.4 DH-PD NoTable 2The angle of crack propagation with respect to the original crack.Case I Case II Case III Case IV Case V Case VILeft 56.1 ◦ ◦ ◦ ◦ ◦ ◦ Right 57.6 ◦ ◦ ◦ ◦ ◦ ◦ . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 The crack propagation speed is computed by V l − . = ∥ x l − x l − ∥ t l − t l − , (34)where x l and x l − are the positions of the crack tip at the times t l and t l − , respectively. An expression for the Rayleighwave speed c R [43] is given as c R c s = . + . ν + ν , (35)where ν is the Poisson’s ratio, c s = √ µ/ρ is the shear wave speed and µ is the shear modulus. The crack starts topropagate at 20 µ s. The crack speed in all four cases simulation is shown in Fig. 21. The maximum crack propagationspeeds are 1669 m/s, 1658 m/s, 1699 m/s, 1664 m/s for Case I, II, III, IV, respectively. All of them are within thelimit of 75% of Rayleigh speed. The crack angle of Kalthoff–Winkler experiment by other methods [37,34] are 65.1 ◦ and 63.8 ◦ as shown in Figs. 22(a) and Fig. 22(b), respectively. Therefore, it can be concluded that with the presentformulation, the random material points distribution has limited influence on the crack propagation when simulatedby DH-PD. H. Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20
7. Conclusions
This paper contributes to the development of dual-horizon peridynamics. The dual property of dual-horizon, auniversal principle for dual-horizon, was proved. Based on the dual property of dual-horizon, we re-proved that theDH-PD fulfills both the balance of momentum and balance of angular momentum. Two issues of crack pattern onrandom material points distribution and multiple materials are discussed.Several numerical examples were presented to demonstrate the capabilities of DH-PD in dealing with issue of ghostforces, multiple materials problems, and crack stability problems on random particle arrangement. The first numericalexample shows that spurious wave reflection happened in the traditional bond based peridynamics, while the DH-PDcan handle it with ease. The second numerical example shows that DH-PD can handle the multiple materials problemwith minimal modification on traditional PD. The last numerical example shows that the crack patterns simulatedby bond based DH-PD are stable and agrees well with the experiment result despite the irregular material pointsdistribution, while that by traditional bond based PD are affected by the irregular material points distribution.There are several advantages of DH-PD over constant horizon PD.(1) The number of neighbors for each material point in DH-PD is less than that in constant-horizon PD when usingan inhomogeneous discretization of the computational domain, therefore reduced the computational cost greatly. Forcertain small material points, it has very limited neighbors but may has many dual-neighbors as it is contained innearby point’ horizons; in the constant-horizon peridynamics, the small material point may have significantly largeneighbor list, as discussed in Section 2.1.(2) The DH-PD improves the computational efficiency in force summations since the force density in DH-PD is onlycalculated once while the force density in traditional PD is calculated twice, as shown in Section 2.4.(3) The DH-PD can deal with the issue of multiple materials with minimal modification of the traditionalperidynamics.(4) The DH-PD can improve the stability of crack propagation and reduce the width of the smeared fractures comparedwith the constant-horizon peridynamics when using a inhomogeneous discretization of the computational domain.
Acknowledgments
The authors acknowledge the supports from the National Basic Research Program of China (973Program: 2011CB013800) and NSFC (51474157), Shanghai Municipal Commission of Science and Technology(16QA1404000), the Ministry of Science and Technology of China (SLDRCE14-B-31) and the ERC-CoG(Computational Modeling and Design of Lithium-ion Batteries (COMBAT)).
Appendix
A.1. Volume correction
In order to account for the effect of material points partly falling inside the horizon, some volume correctiontechniques are utilized to improve the accuracy of force summation in horizon. The commonly used volume correctionmethod in constant-horizon peridynamics is [28] ν x p − x i = − ∥ x p − x i ∥ r i + δ r i + if δ − r i ≤ ∥ x p − x i ∥ ≤ δ, ∥ x p − x i ∥ ≤ δ − r i , r i = ∆ / ∆ is the spacing between material points, or the material point size. In the dual-horizonperidynamics, as any material point has its unique horizon, the volume ratio of big material point with small materialpoint can be big. For example, one case as shown in Fig. 23, the center of material point j is falling outside of i ’s horizon but with certain nontrivial part (green domain in Fig. 23) of j belonging to H i . In order to improve theaccuracy of integration in H i , the green domain in Fig. 23 should be added to calculate the reaction force from j . Asthe material point is generated from finite element, the shape of each material point is irregular, the material point sizecan be estimated by its volume. Let r i denote material point i ’s radius in the shape of circle or sphere. Let two circles . Ren et al. / Comput. Methods Appl. Mech. Engrg. 318 (2017) 0–20 (or spheres) of radii R and r be located along the x -axis centered at ( , , ) and ( d , , ) , respectively, as shown inFig. 24. The intersection area A (or volume V ) between two circles (or spheres) is A = r cos − d + r − R dr + R cos − d + R − r d R −
12 [ ( − d + r + R )( d + r − R )( d − r + R )( d + r + R ) ] / . (A.2) V = π d ( R + r − d ) ( d + dr − r + d R + r R − R ). (A.3) References [1] T. Belytschko, N. Mo¨es, S. Usui, C. Parimi, Arbitrary discontinuities in finite elements, Internat. J. Numer. Methods Engrg. 50 (4) (2001)993–1013.[2] T. Rabczuk, G. Zi, S. Bordas, On three-dimensional modelling of crack growth using partition of unity methods, Comput. Struct. 88 (2010)1391–1411.[3] I. Babuˇska, J.M. Melenk, The Partition of the Unity Finite Element Method, Technical Report BN-1185, Institute for Physical Science andTechnology, University of Maryland, Maryland, 1995.[4] K.W. Cheng, T.P. Fries, Higher-order XFEM for curved strong and weak discontinuities, Internat. J. Numer. Methods Engrg..http://dx.doi.org/10.1002/nme.2768.[5] A. Gravouil, N. Moes, T. Belytschko, Non-planar 3D crack growth by the extended finite element and level sets - Part II: Level set update,Internat. J. Numer. Methods Engrg. 53 (2002) 2569–2586.[6] Y. Jia, G. Zhang, X. Xu, T. Zhuang, Rabczuk, Reproducing kernel triangular bspline-based FEM for solving PDEs, Comput. Methods Appl.Mech. Engrg. 267 (2013) 342–358.[7] Gilles Lubineau, Yan Azdoud, Fei Han, Christian Rey, Abe Askari, A morphing strategy to couple non-local to local continuum mechanics,J. Mech. Phys. Solids 60 (6) (2012) 1088–1102.[8] Pablo Seleson, Samir Beneddine, Serge Prudhomme, A force-based coupling scheme for peridynamics and classical elasticity, Comput. Mater.Sci. 66 (2013) 34–49.[9] S. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, J. Mech. Phys. Solids 48 (1) (2000) 175–209.[10] Xin Lai, Bo Ren, Houfu Fan, Shaofan Li, CT Wu, Richard A. Regueiro, Lisheng Liu, Peridynamics simulations of geomaterial fragmentationby impulse loads, Int. J. Numer. Anal. Methods Geomech. (2015).[11] W. Gerstle, N. Sau, S. Silling, Peridynamic modeling of concrete structures, Nucl. Eng. Des. 237 (12) (2007) 1250–1258.0