A micropolar peridynamics model with non-unified horizon for damage of solids with different non-local effects
aa r X i v : . [ c s . C E ] S e p A micropolar peridynamics model with non-unified horizon fordamage of solids with different non-local effects
Yiming Zhang a, , Xueqing Yang a , Xiaoying Zhuang b,c a School of Civil and Transportation Engineering, Hebei University of Technology, Xiping Road 5340,300401 Tianjin, P.R.China b Department of Geotechnical Engineering, Tongji University, Siping Road 1239,200092 Shanghai, P.R.China c State Key Laboratory for Disaster Reduction in Civil Engineering, Tongji University, Siping Road 1239,200092 Shanghai, P.R.China
Abstract
Most peridynamics models adopt regular point distribution and unified horizon, limitingtheir flexibility and engineering applications. In this work, a micropolar peridynamicsapproach with non-unified horizon (NHPD) is proposed. This approach is implemented ina conventional finite element framework, using element-based discretization. By modifyingthe dual horizon approach into the pre-processing part, point dependent horizon and non-unified beam-like bonds are built. By implementing a domain correction strategy, theequivalence of strain energy density is assured. Then, a novel energy density-based failurecriterion is presented which directly bridges the critical stretch to the mechanical strength.The numerical results indicate the weak mesh dependency of NHPD and the effectivenessof the new failure criterion. Moreover, it is proven that damage of solid with differentnon-local effects can lead to similar results by only adjusting the mechanical strength.
Keywords:
Bond-based peridynamic, Non-unified horizon, Novel failure criterion,Non-local effects, Finite element framework
1. Introduction
Cracks greatly degrade the durability of structures. Predicting the initiations andpropagations of cracks can assist the researchers to design better structures and maintaintheir performances. While partial differential equations are commonly used for describingthe balance relations of continuum physical fields, Cracks, i.e. discontinuities will introducesingularities in these equations. Hence it is a challenge to account in cracks in conventionalcontinuum-based frameworks.Despite of the difficulties, in last decades many researchers presented numerical ap-proaches for taking into account discontinuities in a continuous-discontinuous framework,such as remeshing and interface elements [1–3], numerical manifold method and extendedfinite element methods [4–8], cracking elements method [9–11], phase field method [12–14], and particle based methods [15–17]. Most of these methods deal with discontinuities
Preprint submitted to Theoretical and Applied Fracture Mechanics September 18, 2020 y moving boundaries, localizing strain, or smearing damaging. They mostly treat thecontinuous and discontinuous domains in different manners.Differently, peridynamics (PD) is a non-local theory solving the continuous-discontinuousproblems in the same framework [18, 19]. It uses the integral balance equations and avoidsthe singularity of partial differential equations when discontinuities appear. In the threetypes of PD formulation: bond-based, ordinary state-based, and non-ordinary state-based,bond-based peridynamics (BBPD) is the first proposed and the most popular one [20, 21].After introducing extra parameters or local/global rotational freedom degrees, the fixedPoisson’s ratio problem was solved [22–25]. Moreover, the numerical procedure of BBPDshows similarity to the other classical lattice element models [26, 27], easing its implemen-tations. From numerical point of view, when considering PD as a special type of meshfreemethod, the bond can be consider as a medium for accelerating the integration process[28]. Despite of the great success of PD, most PD formulations use regular grids of mate-rial points and unified horizon, which greatly limits the application and flexibility of PDmethods.In this work, inspired by the dual-horizon peridynamics proposed in [29, 30]. Thebond-based micropolar peridynamics with shear deformability [31, 32] was modified intoa peridynamics approach with non-unified horizon (NHPD). The approach is built in theframework of conventional finite element method (FEM), where the bonds are treated asbeam elements. By using the standard pre-processing step of FEM, it is proven that dual-horizon processing is a part of the modeling (pre-processing) procedure. The irregularmeshes built by Gmsh [33] are used for providing the peridynamics model with point-dependent horizon, the sizes of which can vary in space greatly. Then, an iterative domaincorrection strategy modified from the surface correction strategy is proposed, insuring theequivalence of strain energy density. Finally, a novel energy density based failure criterion isproposed, correlating the critical stretch to the experimental obtained mechanical strength.Comparing to the former PD methods, the proposed approach shows great flexibility re-garding discretization. Some benchmark tests are considered, indicating its reliability androbustness. Last but not least, results in the numerical tests show inspiring correlationsamong the size of horizon, stiffness, and strength of the structures. With the proposedframework, the damage of solid with different non-local effects can give similar results.The remaining parts of this paper are organized as follows: in Section 2, the NHPDis presented in details including the pre-processing, beam-like elemental matrix, iterativedomain correction strategy, and the damage criterion. In Section 3, several benchmarktests are used for demonstrating its robustness and reliability. And the relationship amongsize of horizon, stiffness and strength is revealed Finally, Section 4 contains concludingremarks.
2. The micropolar peridynamics with non-unified horizon (NHPD)
The theories of peridynamics, dual-horizon peridynamics, and correlated non-local op-erator can be found in such as [34–36], which will not be provided in details in this work. In2his Section, the NHPD will be proposed in a way like conducting a numerical simulationstep by step, from modeling to calculation.
For pre-precessing, standard FEM discretization is used in this approach. The domainis discretized into elements, then the nodes are transformed into material points, assuringthe equivalence of volume. For example, in Figure 1, the domain is discretized by lineartriangular elements. Then the volume of every triangular is divided equally to the threenodes, providing the material points. In this step, other types of elements can also be used,such as Voronoi diagrams [37]. (a) (b)
Figure 1: Discretization from elements to material points: (a) planar elements, (b) material points (sizesof the circle indicates their volumes, which are scaled for illustration)
Then a point to point distance checking is run for obtaining the distance between thepresent point to the nearest point, denoted as d min . The example illustrated in Figure 1is used again. The shortest distances between points are shown in Figure 2(a) where thearrows pointed from the present point to its nearest point. For every material point, thesize of its horizon equals to λ d min and λ ≥ λ = 3, thenon-unified horizon of the example is illustrated in Figure 2(b). Herein, λ can be consideredas a non-local parameter in NHPD model. With the increasing of λ , the non-local effectof PD model is enhanced. The influences of λ is studied in the numerical studies.After determining the sizes of horizon of every points, the bonds will be built. When apoint locates in the horizon of another point, a bond connecting these two points will beintroduced. The size of horizon will be used in the elemental matrix of the bond. For theconventional BBPD models, the horizon of bonds equal to the unified horizon. However inthis approach the horizon of bonds are non-unified as well, which are determined locally.Assuming a bond connects two points: A and B with the sizes of horizon h A and h B respectively and the length of this bond is l AB . Then there are two conditions: i) bothpoints locate inside the horizon of the other, ii) one point locates inside the horizon ofthe other, while the other point locate outside, see Figure 3 for example. The size of thehorizon of this bond is H AB = (cid:26) ( h A + h B ) / h B > h A > l AB ,h A h A > l AB > h B . (1)3 a) (b) Figure 2: Determination of the horizon: (a) Shortest distances d min between points, (b) Non-unifiedhorizon with λ =3 H AB will be used for building the stiffness matrix of the bond A − B . When building thebonds, the forces are pairwise introduced. There will be no ghost forces in the domain.The proposed procedure is consistent with the dual-horizon peridynamics [29, 30]. In ourmodel, this procedure is done in the pre-processing step, which is simpler. Figure 3: Two cases regarding a bond connecting two points with different sizes of horizon
The micropolar peridynamics model with shear deformability proposed in [31, 32] isused in this approach. Still, the bond A − B connecting points A and B is considered. For2D condition, the coordinations of A and B are [ x A , y A ] T and [ x B , y B ] T . The rotationalmatrix is defined as 4 AB = a b − b a a b − b a , (2)where a = ( x B − x A ) / l AB and b = ( y B − y A ) / l AB . x A and y A are the x and y coordinatesof point A . R AB is the same as the rotational matrix used in the beam element.Three degrees of freedom: two displacements along two axis and one rotational dis-placement are considered on each point, denoted as u A = [ ux A , uy A , m A ] T and u B =[ ux B , uy B , m B ] T respectively. By mimicking the beam element, the normal, shearing, androtational deformations of the bond are denotes as [ s, γ, ϑ ] T . For bond A − B , sγϑ AB = B TAB R AB (cid:20) u A u B (cid:21) , where B TAB = 1 l AB − − l AB / l AB /
20 0 − l AB l AB . (3)Then, during elastic loading the normal, shearing, and momentum forces of the bondare denoted as [ f n , f t , m ϑ ] T . For bond A − B , f n f t m ϑ AB = Ω AB α AB D AB sγϑ AB , (4)where D AB is the spring equivalent matrix as D AB = k n,AB k t,AB k ϑ,AB . (5)In Eq. 4, k n,AB , k t,AB , and k ϑ,AB are the bond normal, shearing, and rotational spring5quivalent stiffness factors, determined by k n,AB = c AB k t,AB = 12 d AB / l AB k ϑ,AB = d AB / l AB where c AB = Eπ t H AB (1 − ν ) plane stress6 Eπ t H AB (1 − ν ) (1 + ν ) plane strain , and d AB = E (1 − ν )6 π t H AB (1 − ν ) plane stress E (1 − ν )6 π t H AB (1 − ν ) (1 + ν ) plane strain , (6)where E is the elastic modulus and ν is the Poisson’s ratio. t is the thickness. For simplicity, t = 1 m is considered in this work.In Eq.4, α AB is the length correction coefficient. It is introduced for accounting theinfluences of bonds with different lengths. The short bonds are considered to have greaterinfluences on the mechanical responses than the long bonds. α AB is determined by takingthe mean value of the normalized values of l AB regarding points A and B as α AB = 12 (cid:20) exp (cid:18) l Amin − l AB l Amax − l Amin (cid:19) + exp (cid:18) l Bmin − l AB l Bmax − l Bmin (cid:19)(cid:21) , (7)where l Amax and l Amin are the maximum and minimum lengths of bonds connecting topoint A , and l Bmax and l Bmin are the maximum and minimum lengths of bonds connectingto point B . Ω AB is the domain correction coefficient which will be discussed in the nextsection.Correspondingly, the potential energy for the bond A − B , denoted as E AB , is deter-mined by E AB = Z Z f n s l AB f t γ l AB m ϑ ϑ dV A dV B , (8)in which, l AB appears only in the first two terms. Then, the beam-like elemental stiffnessmatrix K AB is determined by K AB = Ω AB α AB V A V B R TAB B AB L AB D AB B TAB R AB , where L AB = l AB l AB . (9)6n Eq. 9, V A and V B are the volumes of the material points A and B respectively. K AB willbe assembled into the global stiffness matrix one after another, just like the conventionalFEM models.Correspondingly E AB can be approximately determined by E AB = 12 (cid:18) [ u A u B ] K AB (cid:20) u A u B (cid:21)(cid:19) . (10) The iterative domain correction strategy is inspired by the energy-based surface cor-rection strategy [34, 38], which is used for correcting the stiffness of the surface materialpoints whose non-local effects are different from those of the inner points. When non-unified horizon inevitably introduce non-homogenized material points and more complexpoint to point bonds, the equivalence of the strain energy density cannot be insured auto-matically. Hence, all bonds in the domain need to be corrected. Comparing to the originalenergy-based correction method, another main difference of the proposed strategy is thatthe correction strategy will be run iteratively.Firstly, assuming a domain experiences unified normal strain ε along a specified direc-tion, the strain energy density e in the domain can be determined by e = E ε − ν ) plane stress E (1 − ν ) ε ν ) (1 − ν ) plane strain , (11)which shall be the true value for each material point.Furthermore, the trail value of the strain energy density on material point A , denotedby ˜ e A , can be determined by˜ e A = 14 n X B (cid:26) [ u A u B ] K AB,trail V A (cid:20) u A u B (cid:21)(cid:27) where K AB,trail = Ω
AB,j α AB V B R TAB B AB L AB D AB B TAB R AB , (12)where 1 / n is the set ofall points connecting to point A by bonds. Ω AB,j is the value of Ω AB at correction iterationstep j with Ω AB, = 1. Here we would like to mention that the correction iteration step j does not relate to the Newton-Raphson iteration step. The iterative domain correctionstrategy will be conducted before the main calculation starts. Once Ω AB is obtained, thevalues will not change during the calculation.Then, with Eq. 12, firstly applying ε along the x direction by setting u x = ε x , thenapplying ε along the y direction by setting u y = ε y , correspondingly the strain energy7ensity of points A and B along x and y direction: ˜ e A,x , ˜ e A,y and ˜ e B,x , ˜ e B,y will be obtained.With these values, the domain correction factor of bond A − B , Ω AB is determined byΩ AB,j +1 = Ω AB,j q ( a / p x ) + ( b / p y ) with p x = 12 (cid:18) ε ˜ e A,x + ε ˜ e B,x (cid:19) , p y = 12 (cid:18) ε ˜ e A,y + ε ˜ e B,y (cid:19) , (13)where a and b are the same as denoted in Eq. 2. When P AB {| Ω AB,j − Ω AB,j − |} < − , thedomain correction strategy will be stopped and Ω AB = Ω AB,j . The peridynamics theory shows differences from the conventional continuum-basedmethod the investigations of which are still undergoing. There are many different damagecriteria on the market, see [39–45] for example.In this work, a novel energy density based criterion is proposed. Based on Eq. 11, thecritical strain energy density under uni-axial tension with tensile stress equals to the tensilestrength F t is e = F t E (1 − ν ) plane stress F t (1 − ν ) (1 − ν ) E (1 − ν ) plane strain . (14)Focusing on the point A with bond A − B connecting to another point B , the balancerelations of the forces of all bonds connecting to A are fulfilled, see Figure 4. Whenconsidering only the stretch s of the bond A − B , based on Eq. 12, the dedication of thebond A − B to the strain energy density at uni-axial loading condition along the x -axiscan be considered as half of the total strain energy density e A as e A AB α AB V A V B c AB l AB s V A + V B ) , (15)where the stain energy density at point B is also considered.Setting e A = e , the corresponding stretch s is considered to be the critical stretch s .When using the same procedure on point B , the same result will be obtained. Finally thecritical stretch s of the bond A − B is obtained as s = s e ( V A + V B )Ω AB α AB V A V B c AB l AB , (16)8 a) (b) Figure 4: Bonds connecting to point A : (a) forces balance at point A , (b) decomposed forces along the x and y axes where, the subscript AB is ignored for s for simplicity. Eq. 16 directly correlates the crit-ical stretch to the experimentally obtained mechanical strength, bringing great flexibilityfor engineering practices.After determining the critical stretch s , the isotropic damage model presented forsome other PD formulations can also be used, such as the bilinear softening model [46, 47]and the exponential softening model [48]. However, this is beyond the topic of this work.Hence, conventional prototype microelastic brittle (PMB) is used here that once s AB ≥ s ,the bond A − B is assumed to break completely and its damage degree d AB is set to 1otherwise d AB = 0. Hence, for a material point A with some damaged bonds, its pointdamage degree ϕ A is defined as ϕ A = 1 − n P B { (1 − d AB ) Ω AB α AB } n P B { Ω AB α AB } . (17) ϕ is determined in the end of every load step, which can be considered as a post-processingstep.During the numerical iteration, in one step, damaging too many bonds may resultin numerical instability and overestimation of the damage zone. The implicit iterationprocedure is adopted for enhancing the numerical stability [49]. For convenience, thefollowing global matrix and vectors are defined: K = [ ( K AB ) and U = [ (cid:18)(cid:20) u A u B (cid:21)(cid:19) , (18)where S ( · ) denotes the assemblage of the beam-like elemental matrix or vector to theglobal form. According to the Newton-Raphson (N-R) method, for the iteration step j atthe load step i , the element-related incremental relation is U i,j = U i − + ∆ U j − | {z } + ∆∆ U | {z } known unknown , (19)9n which ∆ ( · ) denotes an increment of the corresponding value at the preceding load step, i −
1, while ∆∆ ( · ) stands for an increment of the value at the last N-R iteration step, j − K j ∆∆ U = F i − K j ( U i − + ∆ U j − ) (20)where F i is the loading forces at load step i . The total elastic energy of the system E = P AB { E AB } is used for checking whether the equilibrium iteration by means of the N-Rmethod converges. Thus, if if (cid:12)(cid:12)(cid:12)(cid:12) E j − E j − E j (cid:12)(cid:12)(cid:12)(cid:12) < ǫ, (21)then the N-R iteration converged at step j , where ǫ is a prescribed small value with ǫ = 10 − in all numerical examples. When the equilibrium iteration converges, the breakage of bondswill be checked. The value φ = s − s of every bond is obtained. Then, the indexes ofunbroken bonds with φ > φ . With a prescribed number o , the first o bonds in this list will be broken. And theN-R iteration will be rerun. When this list becomes empty in one iteration step, the N-Riteration of this load step converges. The algorithm of the described procedure is illustratedin Figure 5. The computing efficiency will be enhanced with the increasing of o while thenumerical stability will be reduced. o ≤
10 is recommended. It can be found that thisprocedure is similar to that of the cracking elements method [11, 52], which cracks theelement one after another.
3. Numerical investigations
Plane stress condition is considered for all the numerical examples provided in thissection.
The model, material and meshes of the intact disk test are shown in Figure 6. Threemeshes are considered. The analytical peak load per unit thickness is F peak = ( π D F t ) / λ are considered.The force-displacement curves are shown in Figure 7. From the results it can be found: • When λ ≥
2, the stiffness of the structure are generally similar with different λ ; • When λ ≥
2, the stiffness of the structure are generally similar with different meshes; • The values of the peak load increase considerably with the increasing of λ ; • The values of the peak load are slightly different with different meshes; • When λ = 3, the values of the peak load approach the analytical value.10 oading step i Convergent?N-R iteration NoBond id listempty?Yes Damage thefirst o bonds inthe id listNo Next step i = i + 1 Yes Calculate ϕ = s - s , then order the bondid list based on ϕ Figure 5: Calculation procedure within one N-R iteration step
Mesh IIIMesh II
Mesh I
F, d D = 10 cm xy E = 15 GPa (cid:1) = 0.21 F t = 3.81 MPaunit thickness Figure 6: Intact disk test: model, material and meshes λ , seeFigure 8. Because when λ increases, the global stiffness matrix becomes denser. F ( k N ) d (mm) λ = 4Mesh IMesh IIMesh III0.00150.00300.00450.00598.47750.00 0 0.04 0.08 0.12 0.16 0.2 F ( k N ) d (mm) λ = 3Mesh IMesh IIMesh III 0.00100.00200.00300.00400.00500.00562.50 0 0.03 0.06 0.09 0.12 0.15 F ( k N ) d (mm) λ = 2Mesh IMesh IIMesh III0.0090.00180.00270.00360.00450.00 0 0.02 0.04 0.06 0.08 0.1 0.12 F ( k N ) d (mm) λ = 1.5Mesh IMesh IIMesh III Figure 7: Intact disk test: force-displacement curves considering λ = 1 . λ = 2, λ = 3, and λ = 4 Furthermore, we follow the finding that the peak load changes with λ . Considering theresults with Mesh I, we obtain the equivalent tensile strength F eqt from F peak , dependingon λ . The results is illustrated in Figure 9 where the fitting curve is F eqt F t = 3 λ − . (22)Finally, the damage degree and deformation plots are shown in Figures 10 to 12. Gen-erally similar and reasonable failure patterns are found. The damage initiate from the12 .0050.00100.00150.00200.00250.00300.00 2 2.5 3 3.5 4 C P U ti m e f o r it e r a ti on ( s ) λ ( (cid:0) )Mesh IMesh IIMesh III Figure 8: Relationship between computing time and λ E qu i v a l e n t t e n s il e s t r e ng t h F t e q ( M P a ) F t e q / F t ( − ) λ (−) Figure 9: Relationship between equivalent tensile strength F eqt and λ (results of Mesh I) Figure 10: Intact disk test: damage degree and deformation plots (scale: 1:20) of Mesh IFigure 11: Intact disk test: damage degree and deformation plots (scale: 1:20) of Mesh II
Brazilian disk tests with a single slot and multiple slots were experimentally investigatedin [53, 54]. The models are shown in Figure 13. Mesh I shown in Figure 6 is used. Theslots are not explicitly modeled but the bonds intersect with the slots are removed. Aboutthe material properties, same values of E , ν , and η as taken in the last example are used.On the other hand, we consider different values of λ . Hence the values of F t is adjustedbased on Eq. 22, as: i) λ = 2 . , F t = 4 .
69 MPa, ii) λ = 3 , F t = 3 .
81 MPa, and iii)14 igure 12: Intact disk test: damage degree and deformation plots (scale: 1:20) of Mesh III λ = 3 . , F t = 3 .
21 MPa. 598 .
47 kN is used for obtaining normalized peak loads for allcases.For disk tests with an inclined slot, the force-displacement curves and normalized peakloads are shown in Figure 14, indicating agreeable results comparing to the results providedby phase field method [55] and Cracking Elements Method [11]. The damage degree anddeformation plots considering α = 30 ◦ and α = 60 ◦ are shown in Figures 15 and 16. (a) (b) Figure 13: Disk tests with slots: model (a) disk tests with a single slot (different angles of inclination α ),(b) disk tests with multiple slots For disk tests with multiple slots, the force-displacement curves and normalized peakloads are shown in Figure 17, indicating weak dependency between the results and λ after15 (cid:23)(cid:24)(cid:25)(cid:26)(cid:27)(cid:28)(cid:29)(cid:30)(cid:31) !" ;< => no r m a li ze d p ea k l o a d ( − ) α (°)phase field methodCEM λ = 2.5 λ = 3 λ = 3.5 F ( k N ) d (mm) λ = 2.5 α = 15° α = 30° α = 45° α = 60° α = 75° 0.0080.00160.00240.00320.00400.00 0 0.02 0.04 0.06 0.08 0.1 F ( k N ) d (mm) λ = 3 α = 15° α = 30° α = 45° α = 60° α = 75°0.0080.00160.00240.00320.00400.00 0 0.02 0.04 0.06 0.08 0.1 F ( k N ) d (mm) λ = 3.5 α = 15° α = 30° α = 45° α = 60° α = 75° Figure 14: Disk tests with an inclined slot: force-displacement curves and normalized peak loads comparingto the results obtained with phase field method provided in [55] and CEM in [11]Figure 15: Disk tests with an inclined slot: damage degree and deformation plots (scale: 1:20) with α = 30 ◦ igure 16: Disk tests with an inclined slot: damage degree and deformation plots (scale: 1:20) with α = 60 ◦ using adjusted F t . The damage degree and deformation plots are shown in Figures 18to 20, comparing to the experimental results provided in [54]. Generally the patterns ofcracking are similar. Plate made of PMMA with an inclined slot is a benchmark test for PD provided in [34],which was experimentally investigated in [56]. In the experiments, the cracks propagateaxis-symmetrically. The model, material and mesh of the test are shown in Figure 21. Thebonds intersect with the slot are removed for implicitly modeling the slot. This strategyinevitably makes the crack tips a little coarse, see Figure 22. This example is used fortesting the influences of mesh on crack propagation. There is a refined region on theleft side of the model for checking whether the crack will be attracted by this region.Different values of λ are considered and the values of F t is adjusted based on Eq. 22, as:i) λ = 2 . , F t = 14 .
286 MPa, ii) λ = 3 , F t = 10 MPa, and iii) λ = 4 , F t = 7 .
273 MPa.The relationship between the peak loads and the inclined angle is illustrated in Fig-ure 23. Considering different values of λ , after adjusting F t , the obtained peak loads aregenerally similar. The damage degree plots are shown in Figure 24. In most cases, similarto the experiments, axis-symmetrical cracks (damaged regions) are obtained. The crack isnot attracted by the refined region on the left side. Some unexpected branches are foundin some cases, such as in the case with λ = 2, θ = 62 . ◦ , which we attribute mainly to thecoarse modeling of the crack tips.
4. Conclusions
In this work, we present a micropolar peridynamics model with non-unified horizon(NHPD). The main features are summarized as17 @ABCDEFGHIJKLMNOPQR S no r m a li ze d p ea k l o a d ( − ) TUVWXY Z[ \]^_‘ a −) bcd λ e fgh λ i j λ k lmnopqrstuvwxyz{|}~(cid:127)(cid:128)(cid:129)(cid:130)(cid:131)(cid:132)(cid:133)(cid:134)(cid:135)(cid:136)(cid:137)(cid:138)(cid:139)(cid:140)(cid:141)(cid:142)(cid:143)(cid:144)(cid:145)(cid:146)(cid:147)(cid:148)(cid:149) (cid:150) (cid:151)(cid:152)(cid:153)(cid:154) (cid:155)(cid:156)(cid:157)(cid:158) (cid:159)(cid:160)¡¢ F £⁄¥ƒ d §¤'“ λ « ‹›fi α fl (cid:176)–† α ‡ ·(cid:181)¶ α • ‚„”»…‰(cid:190)¿(cid:192)`´ˆ˜¯˘˙¨(cid:201)˚¸(cid:204)˝˛ˇ—(cid:209)(cid:210)(cid:211)(cid:212)(cid:213)(cid:214)(cid:215)(cid:216)(cid:217)(cid:218)(cid:219)(cid:220)(cid:221)(cid:222)(cid:223)(cid:224)Æ (cid:226) ª(cid:228)(cid:229)(cid:230) (cid:231)ŁØŒ º(cid:236)(cid:237)(cid:238) F (cid:239)(cid:240)æ(cid:242) d (cid:243)(cid:244)ı(cid:246) λ (cid:247) łøœ α ß (cid:252)(cid:253)(cid:254) α (cid:255) =(cid:0)(cid:1) α (cid:2) (cid:3)(cid:4)(cid:5) 0(cid:6)(cid:7)(cid:8)5(cid:9)(cid:10)(cid:11)(cid:12)1(cid:13)(cid:14)(cid:15)(cid:16)(cid:17)(cid:18)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)2(cid:24)(cid:25)(cid:26)(cid:27)(cid:28)(cid:29)(cid:30)(cid:31) !"3 F :;<> d ?@AB λ C D α E FGH α I JKL α M NOP
Figure 17: Disk tests with multiple slots: force-displacement curves and normalized peak loads comparingto the results obtained with phase field method provided with CEM in [52]
QRSTUVWXYZ[
Figure 18: Disk tests with 2 slots: damage degree and deformation plots (scale: 1:20) comparing to theexperimental results provided in [54] ]^_‘abcdef Figure 19: Disk tests with 3 slots: damage degree and deformation plots (scale: 1:20) comparing to theexperimental results provided in [54] ghijklmnopq
Figure 20: Disk tests with 4 slots: damage degree and deformation plots (scale: 1:20) comparing to theexperimental results provided in [54]Figure 21: Plate with an inclined slot: model, material and mesh r −2−1 s t − u −2 −1 v wyxz{| } ~(cid:127)(cid:128)(cid:129) λ (cid:130) (cid:131)(cid:132)(cid:133)(cid:134) θ = 0° − − − − − − (cid:135)(cid:136)(cid:137)(cid:138)(cid:139) (cid:140) (cid:141)(cid:142)(cid:143)(cid:144) λ (cid:145) (cid:146)(cid:147)(cid:148)(cid:149) θ = 62.5° − − − − − − (cid:150)(cid:151)(cid:152)(cid:153)(cid:154) (cid:155) (cid:156)(cid:157)(cid:158)(cid:159) λ (cid:160) ¡¢£⁄ θ ¥ ƒ§ − − ¤ − ' − − “ − « ‹›fifl(cid:176) – †‡·(cid:181) λ ¶ •‚„” θ » …‰(cid:190)¿(cid:192) − − ` − ´ − − ˆ − ˜ ¯˘˙¨(cid:201) ˚ ¸(cid:204)˝˛ λ ˇ —(cid:209)(cid:210)(cid:211) θ (cid:212) (cid:213)(cid:214) − − (cid:215) − (cid:216) − − (cid:217) − (cid:218) (cid:219)(cid:220)(cid:221)(cid:222)(cid:223) (cid:224) Æ(cid:226)ª(cid:228) λ (cid:229) (cid:230)(cid:231)ŁØ θ Œ º(cid:236)(cid:237)(cid:238)(cid:239)
Figure 22: Plate with an inclined slot: bonds around the crack tips for cases θ = 0 ◦ and θ = 62 . ◦ (cid:240)æ(cid:242)(cid:243)(cid:244)ı(cid:246)(cid:247)łøœß(cid:252)(cid:253)(cid:254)(cid:255)2(cid:0)(cid:1)(cid:2)3(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)(cid:9)(cid:10)(cid:11)4(cid:12)(cid:13)(cid:14)(cid:15)(cid:16)(cid:17)(cid:18)(cid:19)(cid:20)5(cid:21)(cid:22)(cid:23)(cid:24) 0 (cid:25)(cid:26) F ((cid:30)(cid:31) θ !" λ = &’) λ * +,- λ . /17 Classical PDExperiments
Figure 23: Plate with an inclined slot: the relationship between the peak loads and the inclined angle,comparing to the experimental results provided in [56] and numerical results by classical PD in [34] = 2.2λ = 3.0Experimentsλ Figure 24: Plate with an inclined slot: the damage degree plots comparing to the experimental resultsprovided in [56] In the pre-processing step, normal FEM discretization is used for providing materialpoints. The horizon varies with different points the size of which depends on theshortest distance between neighboring points. The ratio of size of horizon to theshortest distance equals to a prescribed value λ . When λ increases, the non-localeffects are enhanced; • An iterative domain correction strategy is proposed for assuring the equivalence ofstrain energy density. Then, based on the maximum strain energy density, a novelfailure criterion is proposed for the NHPD which bridges the critical stretch to themechanical strength F t ; • Considering numerical studies regarding different values of λ and different meshes,the results indicate NHPD shows generally weak mesh dependency. Moreover, it isfound that if λ ≥ λ has weak influences on the stiffness of the structure while λ hasgreat influence on the equivalent strength of the structure F eqt . A linear relationshipbetween F eqt / F t and λ is obtained and F eqt = F t when λ = 3; • Considering the linear relationship between F eqt / F t and λ then adjusting the inputed F t , similar results can be obtained regarding different values of λ .The NHPD shows another routine for developing peridynamics models and the relationshipbetween equivalent strength and λ indicates correlations between strength and local/non-local damages.
5. Acknowledgement
The authors gratefully acknowledge financial support by the National Natural ScienceFoundation of China (NSFC) (51809069) and by the Hebei Province Natural Science FundE2019202441 and the 2019 Foreign Experts Plan of Hebei Province.22 eferences [1] P. Areias, J. Reinoso, P. Camanho, and T. Rabczuk, “A constitutive-based element-by-element crack propagation algorithm with local mesh refinement,”
ComputationalMechanics , vol. 56, pp. 291–315, 2015.[2] P. Areias, T. Rabczuk, and D. Dias-da-Costa, “Element-wise fracture algorithm basedon rotation of edges,”
Engineering Fracture Mechanics , vol. 110, pp. 113–137, 2013.[3] E. C. Mejia Sanchez, L. F. Paullo Muoz, and D. Roehl, “Discrete fracture propaga-tion analysis using a robust combined continuation method,”
International Journal ofSolids and Structures , vol. 193-194, pp. 405 – 417, 2020.[4] H. Zheng and D. Xu, “New strategies for some issues of numerical manifold methodin simulation of crack propagation,”
International Journal for Numerical Methods inEngineering , vol. 97, pp. 986–1010, 2014.[5] Y. Yang, G. Sun, H. Zheng, and X. Fu, “A four-node quadrilateral element fitted tonumerical manifold method with continuous nodal stress for crack analysis,”
Comput-ers and Structures , vol. 177, pp. 69–82, 2016.[6] Z. Wu, H. Sun, and L. N. Y. Wong, “A cohesive element-based numerical manifoldmethod for hydraulic fracturing modelling with voronoi grains,”
Rock Mechanics andRock Engineering , vol. 52, pp. 2335–2359, 2019.[7] J.-H. Song, P. Areias, and T. Belytschko, “A method for dynamic crack and shear bandpropagation with phantom nodes,”
International Journal for Numerical Methods inEngineering , vol. 67, pp. 868–893, 2006.[8] J.-Y. Wu and F.-B. Li, “An improved stable XFEM (Is-XFEM) with a novel enrich-ment function for the computational modeling of cohesive cracks,”
Computer Methodsin Applied Mechanics and Engineering , vol. 295, pp. 77–107, 2015.[9] Y. Zhang, R. Lackner, M. Zeiml, and H. Mang, “Strong discontinuity embedded ap-proach with standard SOS formulation: Element formulation, energy-based crack-tracking strategy, and validations,”
Computer Methods in Applied Mechanics and En-gineering , vol. 287, pp. 335–366, 2015.[10] Y. Zhang and X. Zhuang, “Cracking elements: a self-propagating strong discontinu-ity embedded approach for quasi-brittle fracture,”
Finite Elements in Analysis andDesign , vol. 144, pp. 84–100, 2018.[11] Y. Zhang and H. A. Mang, “Global cracking elements: a novel tool for Galerkin-based approaches simulating quasi-brittle fracture,”
International Journal for Nu-merical Methods in Engineering , vol. 121, pp. 2462–2480, 2020.2312] C. Miehe, L.-M. Sch¨anzel, and H. Ulmer, “Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crackpropagation in thermo-elastic solids,”
Computer Methods in Applied Mechanics andEngineering , vol. 294, pp. 449–485, 2015.[13] J.-Y. Wu and V. P. Nguyen, “A length scale insensitive phase-field damage model forbrittle fracture,”
Journal of the Mechanics and Physics of Solids , vol. 119, pp. 20–42,2018.[14] J.-Y. Wu, “A unified phase-field theory for the mechanics of damage and quasi-brittlefailure,”
Journal of the Mechanics and Physics of Solids , vol. 103, pp. 72–99, 2017.[15] T. Rabczuk and T. Belytschko, “Cracking particles: a simplified meshfree method forarbitrary evolving cracks,”
International Journal for Numerical Methods in Engineer-ing , vol. 61, pp. 2316–2343, 2004.[16] T. Rabczuk and T. Belytschko, “A three-dimensional large deformation meshfreemethod for arbitrary evolving cracks,”
Computer Methods in Applied Mechanics andEngineering , vol. 196, pp. 2777–2799, 2007.[17] T. Rabczuk, G. Zi, S. Bordas, and H. Nguyen-Xuan, “A simple and robust three-dimensional cracking-particle method without enrichment,”
Computer Methods in Ap-plied Mechanics and Engineering , vol. 199, pp. 2437–2455, 2010.[18] S. Silling, “Reformulation of elasticity theory for discontinuities and long-range force,”
Journal of the Mechanics and Physics of Solids , vol. 48, pp. 175–209, 2000.[19] H. Fei and L. Gilles, “Coupling of nonlocal and local continuum models by the ar-lequinapproach,”
International Journal for Numerical Methods in Engineering , vol. 89,no. 6, pp. 671–685, 2011.[20] D. Han, Y. Zhang, Q. Wang, W. Lu, and B. Jia, “The review of the bond-basedperidynamics modeling,”
Journal of Micromechanics and Molecular Physics , vol. 04,no. 01, p. 1830001, 2019.[21] J. Trageser and P. Seleson, “Bond-based peridynamics: a tale of two poisson’s ratios,”
Journal of Peridynamics and Nonlocal Modeling , Apr 2020.[22] Q.-Z. Zhu and T. Ni, “Peridynamic formulations enriched with bond rotation effects,”
International Journal of Engineering Science , vol. 121, pp. 118–129, 2017.[23] W.-J. Li, Q.-Z. Zhu, and T. Ni, “A local strain-based implementation strategy forthe extended peridynamic model with bond rotation,”
Computer Methods in AppliedMechanics and Engineering , vol. 358, p. 112625, 2020.[24] X. Gu and Q. Zhang, “A modified conjugated bond-based peridynamic analysis forimpact failure of concrete gravity dam,”
Meccanica , vol. 55, pp. 547–566, Mar 2020.2425] H. Yu, X. Chen, and Y. Sun, “A generalized bond-based peridynamic model for quasi-brittle materials enriched with bond tensionrotationshear coupling effects,”
ComputerMethods in Applied Mechanics and Engineering , vol. 372, p. 113405, 2020.[26] G.-F. Zhao, “Developing a four-dimensional lattice spring model for mechanical re-sponses of solids,”
Computer Methods in Applied Mechanics and Engineering , vol. 315,pp. 881–895, 2017.[27] M. Nikoli´c, E. Karaveli´c, A. Ibrahimbegovic, and P. Miˇsˇcevi´c, “Lattice element modelsand their peculiarities,”
Archives of Computational Methods in Engineering , vol. 25,pp. 753–784, Jul 2018.[28] M. A. Bessa, J. T. Foster, T. Belytschko, and W. K. Liu, “A meshfree unification:reproducing kernel peridynamics,”
Computational Mechanics , vol. 53, pp. 1251–1264,Jun 2014.[29] H. Ren, X. Zhuang, Y. Cai, and T. Rabczuk, “Dual-horizon peridynamics,”
Interna-tional Journal for Numerical Methods in Engineering , vol. 108, pp. 1451–1476, 2016.[30] H. Ren, X. Zhuang, and T. Rabczuk, “Dual-horizon peridynamics: A stable solution tovarying horizons,”
Computer Methods in Applied Mechanics and Engineering , vol. 318,pp. 762–782, 2017.[31] V. Diana and S. Casolo, “A bond-based micropolar peridynamic model with sheardeformability: Elasticity, failure properties and initial yield domains,”
InternationalJournal of Solids and Structures , vol. 160, pp. 201 – 231, 2019.[32] V. Diana and S. Casolo, “A full orthotropic micropolar peridynamic formulation forlinearly elastic solids,”
International Journal of Mechanical Sciences , vol. 160, pp. 140– 155, 2019.[33] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-d finite element mesh generator withbuilt-in pre- and post-processing facilities,”
International Journal for Numerical Meth-ods in Engineering , vol. 79, no. 11, pp. 1309–1331, 2009.[34] E. Madenci and E. Oterkus,
Peridynamic theory and its applications . Springer, 2014.[35] H. Ren, X. Zhuang, and T. Rabczuk, “Nonlocal operator method with numericalintegration for gradient solid,”
Computers & Structures , vol. 233, p. 106235, 2020.[36] H. Ren, X. Zhuang, and T. Rabczuk, “A higher order nonlocal operator method forsolving partial differential equations,”
Computer Methods in Applied Mechanics andEngineering , vol. 367, p. 113132, 2020.[37] X. Gu, Q. Zhang, and X. Xia, “Voronoi-based peridynamics and cracking analysis withadaptive refinement,”
International Journal for Numerical Methods in Engineering ,vol. 112, no. 13, pp. 2087–2109, 2017. 2538] Q. Le and F. Bobaru, “Surface corrections for peridynamic models in elasticity andfracture,”
Computational Mechanics , vol. 61, pp. 499–518, 2018.[39] D. Yang, X. He, X. Liu, Y. Deng, and X. Huang, “A peridynamics-based cohesive zonemodel (pd-czm) for predicting cohesive crack propagation,”
International Journal ofMechanical Sciences , vol. 184, p. 105830, 2020.[40] M. Zaccariotto, F. Luongo, G. sarego, and U. Galvanetto, “Examples of applicationsof the peridynamic theory to the solution of static equilibrium problems,”
The Aero-nautical Journal , vol. 119, no. 1216, p. 677700, 2015.[41] J. T. Foster, S. A. Silling, and W. Chen, “An energy based failure criterion for use withperidynamic states,”
International Journal for Multiscale Computational Engineering ,vol. 9, no. 6, pp. 675–688, 2011.[42] D. Huang, G. Lu, and Y. Liu, “Nonlocal peridynamic modeling and simulation oncrack propagation in concrete structures,”
Mathematical Problems in Engineering ,vol. 2015, p. 858723, Feb 2015.[43] Y. Zhang and P. Qiao, “A new bond failure criterion for ordinary state-based peridy-namic mode II fracture analysis,”
International Journal of Fracture , vol. 215, pp. 105–128, Jan 2019.[44] D. Dipasquale, G. Sarego, M. Zaccariotto, and U. Galvanetto, “A discussion on fail-ure criteria for ordinary state-based peridynamics,”
Engineering Fracture Mechanics ,vol. 186, pp. 378 – 398, 2017.[45] T. Rabczuk and H. Ren, “A peridynamics formulation for quasi-static fracture andcontact in rock,”
Engineering Geology , vol. 225, pp. 42–48, 2017.[46] V. Diana, J. F. Labuz, and L. Biolzi, “Simulating fracture in rock using a micropolarperidynamic formulation,”
Engineering Fracture Mechanics , vol. 230, p. 106985, 2020.[47] C. Xu, Y. Yuan, Y. Zhang, and Y. Xue, “Peridynamic modeling of prefabricatedbeams post-cast with steelfiber reinforced high-strength concrete,”
Structural Con-crete , vol. n/a, no. n/a.[48] Y. Tong, W. Shen, J. Shao, and J. Chen, “A new bond model in peridynamics theoryfor progressive failure in cohesive brittle materials,”
Engineering Fracture Mechanics ,vol. 223, p. 106767, 2020.[49] Y. Bie, S. Li, X. Hu, and X. Cui, “An implicit dual-based approach to couple peri-dynamics with classical continuum mechanics,”
International Journal for NumericalMethods in Engineering , vol. 120, no. 12, pp. 1349–1379, 2019.2650] Y. Zhang and X. Zhuang, “A softening-healing law for self-healing quasi-brittle mate-rials: analyzing with strong discontinuity embedded approach,”
Engineering FractureMechanics , vol. 192, pp. 290–306, 2018.[51] Y. Zhang and X. Zhuang, “Cracking elements method for dynamic brittle fracture,”
Theoretical and Applied Fracture Mechanics , vol. 102, pp. 1–9, 2019.[52] L. Mu and Y. Zhang, “Cracking elements method with 6-node triangular element,”
Finite Elements in Analysis and Design , vol. 177, p. 103421, 2020.[53] H. Haeri, K. Shahriar, M. F. Marji, and P. Moarefvand, “Experimental and numericalstudy of crack propagation and coalescence in pre-cracked rock-like disks,”
Interna-tional Journal of Rock Mechanics and Mining Sciences , vol. 67, pp. 20 – 28, 2014.[54] H. Haeri, A. Khaloo, and M. F. Marji, “Experimental and numerical analysis ofBrazilian discs with multiple parallel cracks,”
Arabian Journal of Geosciences , vol. 8,pp. 5897–5908, 2015.[55] S.-W. Zhou and C.-C. Xia, “Propagation and coalescence of quasi-static cracks inBrazilian disks: an insight from a phase field model,”
Acta Geotechnica , vol. 14,pp. 1195–1214, Aug 2019.[56] M. Ayatollahi and M. Aliha, “Analysis of a new specimen for mixed mode fracturetests on brittle materials,”