Coupling of peridynamics with finite elements for brittle crack propagation problems

Journal Pre-proofs Coupling of peridynamics with finite elements for brittle crack propagation problems Dong Yang, Xiaoqiao He, Shenghui Yi, Yajie Deng, Xuefeng Liu PII: DOI: Reference:

S0167-8442(19)30535-X https://doi.org/10.1016/j.tafmec.2020.102505 TAFMEC 102505

To appear in:

Theoretical and Applied Fracture Mechanics

Received Date: Revised Date: Accepted Date:

23 September 2019 19 December 2019 20 January 2020

Please cite this article as: D. Yang, X. He, S. Yi, Y. Deng, X. Liu, Coupling of peridynamics with finite elements for brittle crack propagation problems, Theoretical and Applied Fracture Mechanics (2020), doi: https://doi.org/ 10.1016/j.tafmec.2020.102505

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2020 Published by Elsevier Ltd.

Coupling of peridynamics with finite elements for brittle crack propagation problems Dong Yang1, Xiaoqiao He1,2*, Shenghui Yi2, Yajie Deng3, Xuefeng Liu1 1Department

of Architecture and Civil Engineering, City University of Hong Kong, Tae Chee Avenue, Kowloon, Hong Kong, China. 2Centre for Advanced Structural Materials, City University of Hong Kong Shenzhen Research Institute, 8 Yuexing 1st Road, Shenzhen Hi-Tech Industrial Park, Nanshan District, Shenzhen, China. 3School of Aerospace Engineering and Applied Mechanics, Tongji University, 1239 Siping Road, Shanghai, China.

Abstract In this work, a new coupling model of finite element method (FEM) and ordinary state-based peridynamics (OSBPD) is proposed to accurately predict brittle fracture under quasi-static and dynamic loading. PD grids are directly coupled with FEM meshes without an overlapped zone. The transferred forces between PD grids and FEM meshes are reciprocal and balanced inherently, which avoids the complexity of the transitional zone-based coupling models and the non-equilibrium phenomenon of the structural force reported in the sharing nodes-based coupling model. More importantly, the proposed coupling model is free of significant influence of ghost forces and significant wave reflection at the coupling interface. Also, the surface effect of PD can be approximately removed. The accuracy of the proposed coupling model is verified through comparing with the FEM results and pure OSBPD results.

Keywords: Coupling model; Peridynamics; Brittle fracture; Crack branching; Dynamic crack propagation

1

1．Introduction Accurate modeling of crack propagation problems is still an active challenge in computational mechanics, especially for complex crack problems, such as crack branch and three-dimensional curve cracks. This challenge originates from the limitation of classical continuum mechanics (CCM), in which governing equations are partial differential equations (PDEs) and hence spatial derivatives cannot be obtained when cracks are involved. In order to overcome the limitation of CCM, Silling [1] introduced a nonlocal form of CCM called ‘peridynamics (PD)’ where the governing equation is in an integro-differential form. Cracks are spontaneously predicted in PD through terminating bond interactions without crack tracking approaches [2, 3] in extended finite element method (XFEM) and remeshing approaches in cohesive zone model (CZM). This implicit nature enables PD to suitably predict complex crack problems [4, 5]. It is noted that the implicit nature of PD is different from the one of the smeared crack approaches such as gradient damage model (GDM) [6] and phase filed model (PFM) [7], in which cracks are the result of damage accumulation. The implicit nature of smeared crack approaches makes it hard to tackle crack problems when the detailed crack characteristics are needed [8]. On the contrary, the detailed crack characteristics can be explicitly obtained in PD, which makes it possible to handle the interactions between crack surfaces, such as cohesive tractions [9, 10] and hydraulic pressure [11, 12]. So far, there are three main PD models such as bond-based PD (BBPD) [1, 13], ordinary state-based PD (OSBPD) [14, 15] and non-ordinary state-based PD (NOSBPD) [16]. Poisson’s ratio is fixed at 1/3 for plane stress case and 1/4 for plane strain and three-dimensional cases in BBPD. State-based PD (SBPD) removes the restriction of Poisson’s ratio, and hence it can describe more general materials. PD has been attracting an increasing interest among researchers [17] due to its powerful ability to simulate various complex crack problems. However, due to the fact that PD is a nonlocal method, this nonlocal nature brings it with some inherent drawbacks, such as the surface effect, low computational efficiency and inflexible boundary conditions. An efficient way to overcome these drawbacks is to couple PD with FEM. The external boundaries and linear elastic parts in the whole computational domain are discretized into FEM zones, while the parts of damage or crack being occurred are discretized into PD zones. This coupling strategy can utilize the respective advantages of both two methods and at the same time avoid the respective drawbacks. To date, considerable efforts have been made to couple PD with FEM. Macek and Silling 2

[18] made the first attempt to implement the PD model into FEM framework using truss elements. The coupling approach proposed by Kilic and Madenci [19] introduced an overlapped region in which both PD and FEM equations are used simultaneously. Agwai et al. [20] proposed a sub-modeling coupling approach which involves FEM for global analysis and PD for sub modeling to predict failure. Liu and Hong [21] introduced interface elements between FEM and PD regions to transfer coupling forces. Also, a morphing approach was proposed by Lubineau et al. [22] to couple local and nonlocal theories, which was further developed to couple state-based peridynamics with FEM [23]. Seleson et al. [24, 25] introduced a blending-based coupling approach to couple PD with classical continuum mechanics. Compared with the transitional zone-based approaches mentioned above, a sharing node-based coupling approach without an overlapped zone was proposed by Galvanetto and coworkers [26, 27]. The FEM elements in the sharing zone which cover PD particles at their corners exert forces only to FEM nodes while bonds connecting PD particles and FEM nodes only apply forces to PD particles. Also, without an overlapped zone, Bie et al. [28] proposed a coupling approach to couple state-based peridynamics with node-based smoothed FEM. Other recent coupling works can be found in Refs. [2932]. It is worth noting that in the transitional zone-based coupling approaches, complex techniques such as mapping, morphing functions and blending functions are needed to determine the transferred force values. For the sharing node-based approach, it is relatively simple, but the non-equilibrium phenomenon of the structural force needs to be carefully reduced in order to achieve higher accuracy and better convergence. How to solve the PD-FEM coupling system is another important issue in order to obtain the better computational efficiency. There are two main ways to solve the PD-FEM coupling system, such as implicit methods [26, 30, 33] and explicit methods [19, 21, 28, 34] respectively. Generally, computation cost at every time step is much less by using explicit methods than by using implicit methods [35], due to the fact that there are no large-matrix operations in explicit methods. For dynamic crack propagation problems of short duration under blast or high velocity impact, explicit methods are more suitable and hence commonly adopted in PD [13, 36, 37]. However, small time step is required in explicit methods to insure numerical stability while implicit methods are normally unconditionally stable. The restriction of small time step limits the effectiveness of explicit methods when it comes to structural problems of moderate duration, especially under static or quasi-static conditions. Fortunately, the adaptive dynamic relaxation (ADR) method [38] was extended in PD to effectively simulate static 3

or quasi-static crack propagation problems [19, 21, 28, 39]. On the other hand, implicit methods are also utilized [26, 30, 33] to obtain static solution of crack propagation problems in PD. But the comparison of numerical efficiency between explicit and implicit methods in PD is still less. A comparison of numerical efficiency between explicit and implicit methods was conducted in the investigation of crystal plasticity [40], which states that the explicit method is faster than the implicit method in the plastic stage, but slower in the elastic stage. It is noted that crack propagation problems are not considered in the comparison. Overall, the comparison of numerical efficiency between explicit methods and implicit methods needs to be further explored, especially when cracks and complex material models involved. In this work, a new simple coupling model of FEM and OSBPD is proposed to effectively investigate dynamic and quasi-static crack propagation problems. Similar to the sharing node-based coupling approach, there is no overlapped zone in the proposed coupling model. The transferred force values between PD grids and FEM meshes can be easily determined through directly computing the interactions between virtual PD material particles and true PD material particles, which avoids complex techniques such as mapping, morphing and blending functions in the transitional zone-based approaches. On the other hand, due to the fact that the determined transferred forces are inherently reciprocal and balanced, the proposed coupling model inherently satisfies the Newton’s third law, and is free of non-equilibrium phenomenon of the structural force reported in the sharing-node based coupling approach. And it is worth noting that the grid size of FEM elements near PD grids is the same as the PD grid size to ensure that the nodes in these FEM elements can be simultaneously regarded as virtual PD material particles. The introduction of these virtual PD material particles can also remove the surface effect of PD zone. Moreover, the influence of ghost forces is limited to some extent and comparable accuracy can be achieved in the proposed coupling model. Explicit difference methods and the adaptive dynamic relaxation method are adopted in the proposed coupling model for dynamic and quasi-static analysis, respectively. First, the accuracy of the proposed coupling model is verified compared with FEM results. And then, some benchmarks for quasi-static and dynamic brittle fracture are further conducted. It is expected that the proposed coupling model can provide an accurate tool to predict the fracture process of brittle materials under static or dynamic loading. The rest of the paper is organized as follows. The theory of OSBPD is reviewed in Section 2, together with the failure model for brittle fracture. Section 3 presents the proposed coupling model of 4

OSBPD and FEM and its numerical implementations. In Section 4, the validity of the proposed coupling model is examined by analyzing three linear elastic numerical examples. More benchmarks are simulated in Section 5, including one quasi-static brittle fracture case and two dynamic brittle fracture cases. Some key comments are concluded in Section 6. All the numerical examples are written using Fortran90, and an open source IDE platform called “Code::Blocks”, together with GNU Fortran compiler, is adopted to compile and run the Fortran90 codes. This computing environment is similar to Visual Studio with Intel Fortran compiler. For all the Fortran codes in this paper (including FEM model, PD model and PD-FEM coupling model), the same computational procedure or program structure is adopted, which is the same as the codes provided in Ref. [41]. It is worth noting that all the codes are serial without applying any parallel computing approaches. ANSYS is used to generate the meshes. The features of the computer are: Intel(R) Core(TM) i7-4790 CPU @ 3.60GHz, 16GB RAM and Windows 7 OS.

2. The theory of OSBPD In this work, the OSBPD model is adopted to be coupled with FEM due to the generality of OSBPD without the limitation of Poisson’s ratio. This section consists of two main parts, i.e. linear elastic constitutive model of OSBPD and the corresponding failure model for describing brittle fracture process. 2.1 Linear elastic constitutive model of OSBPD As shown in Fig. 1, the position vector of a material particle in the reference and deformed configurations is denoted as particles x and

x

and y , respectively. During this deformation process, material

x experience displacements, u and u  , respectively. According to the

OSBPD theory, the displacement field of material particle other material particles

is influenced by its interactions with

x within its horizon H x with an internal length or radius  , as shown

in Fig. 1. Material particles particle

x

x within the horizon H x are called the family members of material

x . Similarly, material particles x

interacts with their family members within

5

H x .

Fig. 1 Diagram of the OSBPD model

As proposed by Silling [15], the motion equation of material particle

x

at time t can be

written as,  u  x , t  

 T  x, t 

x   x  T  x , t  x  x   dV x   b  x , t  ,

(1)

Hx

where  is the density, u is the acceleration of material particle volume of material particle

is the differential

x . b is the body force density vector. x  x is the bond vector in the

reference configuration. T  x, t  x  x  T  x, t  x  x particle

x , and dVx

is the bond force vector that material

x exerts on material particle x . T  x,t  x  x is the force density vector with respect

to x  x in the force state T x,t  and T  x,t  x  x

is the force density vector with respect to

x  x in the force state T  x ,t  . State is a mathematical object defined in state-based peridynamics, operating like a mapping. For instance, the force state T x,t  is basically an array that stores the force density vectors of PD bonds pertaining to material particle

x , and the angle bracket ,

, is a

notation used to extract the force density vector with respect to a certain PD bond. For linear elastic isotropic materials, different constitutive models [14, 41-43] have been proposed in OSBPD to determine the force density vector. The constitutive model proposed in [41] is adopted in this study. The advantage of this OSBPD model is that it can be reduced to the BBPD model. The force density vector T  x ,t  x   x

for three-dimensional case can be expressed as, 6

4v  1  15  9   4 T  x, t  x  x =  3  4 1  2v x   x  where 

is the shear modulus,

deformed configuration. 

v

 y  y  x  x   x  x

y  y ,   y  y

(2)

is the Poisson’s ratio and y   y is the bond vector in the

is the dilation coefficient, which is defined as =

9 4

3



Hx

y  y  x  x dV x  . x  x

(3)

Similarly, the force density vector for two-dimensional case can also be obtained as

 2 3v  1  12    3 T  x, t  x  x =  2 1  v x   x   h

 y  y  x  x   x  x

y  y   y  y

for plane stress,

(4)

for plane strain,

(5)

and

 y  y  x  x  

 2 4v  1  12   T  x, t  x  x =   3 2 1  2v x   x   h

x  x

y  y   y  y

where h is the thickness, and the dilation coefficient is defined as =

2 h

2

y  y  x  x dV x  . x  x



Hx

(6)

It’s noted that the above parameters are derived based on an assumption that the material particle is infinitesimal which means that the complete horizon is a disk for 2D case and a sphere for 3D case. But in fact, this assumption cannot be achieved after the whole domain is discretized. This leads to the fact that the strain energy density at a material particle is less than the one calculated from CCM, which is called ‘discretization error’. Another error is called ‘surface effect’ [44], which is due to the fact that the horizon of the material particle near a boundary is not complete. Surface effect is approximately removed in the proposed coupling model in Section 3. In order to reduce the ‘discretization error’, the effective and direct strain energy density correction method [41] is adopted. 2.2 Failure models Apart from the constitutive model of OSBPD, failure models are essential to accurately describe the failure or fracture behavior of brittle materials. To describe brittle fracture, many different failure models [39, 41, 45] have been proposed in PD. In this work, the critical bond stretch criterion proposed in [41] is adopted due to its simplicity. This criterion can be stated that once bond stretch exceeds the critical bond stretch, the interaction of this bond will be terminated. After this bond is removed, the 7

deformation of this bond is not considered when computing the dilation coefficient. Bond stretch

s

is defined as, s

And the critical bond stretch



(7)

sc is determined by computing the critical energy release rate Gc as,

     sc       

in which,

y  y  x  x . x  x

Gc  5   3   3          4 3      4

Gc 16 6     2   2     9  

3D case

,

(8)

2D case

is the bulk modulus.

Using the above critical bond stretch criterion, the fracture process of brittle materials can be well predicted in OSBPD.

3. Coupling model of OSBPD and FEM In this section, a new simple coupling model of OSBPD and FEM is proposed. First, the detailed coupling strategy is explained. And then, explicit solution methods are described to solve the governing equations of the coupling model for the dynamic and quasi-static analysis. 3.1 Coupling strategy As shown in Fig. 2, a coupling model of a 2D plate is established to describe the proposed coupling strategy in detail. The whole body is discretized into FEM meshes with rectangular elements and uniform PD grids, and PD grids are directly coupled with FEM meshes. The mesh size of FEM elements near PD grids is the same as the PD grid size x , and the nodes in these FEM elements are simultaneously regarded as virtual PD material particles. The zone containing these virtual PD material particles is called the coupling zone. The zone containing the remaining FEM meshes is defined as the FEM zone. It is noted that the element type and mesh size in the FEM zone is not restricted, such as triangular elements or rectangular elements. Assuming the horizon size 𝛿 of the PD particles in the PD zone is equal to 3x , the size of the coupling zone is 2 . This is due to the fact that the maximum influence range of a OSBPD bond is 2 , as the bond between material particles 8

xa and xb shown

in Fig. 2(b). In order to distinguish different nodes or particles, the PD particles in the PD zone are highlighted in red, FEM nodes in the FEM zone are highlighted in blue, and FEM nodes in the coupling zone are highlighted in green. It is noted that green nodes in the coupling zone are in fact FEM nodes, and at the same time these green nodes are regarded as virtual PD particles. The introduction of these virtual PD particles has two roles. The first one is to approximately remove the surface effect of PD zone. The surface effect is an inherent deficiency of PD, which is due to the fact that the horizon of the PD particles close to a boundary is not complete. Virtual particle method is an effective way to ensure that all PD material particles in PD domain have a fully integral horizon [44], and hence the surface effect can be removed. This approach has also been used in the sharing node-based PD-FEM coupling technique. It is worth noting that in the proposed coupling model, when the PD zone is fully wrapped in the FEM zone, the surface effect can be completely removed, because each PD particle in the PD zone has a complete horizon due to the existence of the virtual PD particles. However, when the PD zone is not fully wrapped in the FEM zone, as shown in Fig. 2, the surface effect will be partly removed.

(a) Complete model (b) Enlarged zone Fig. 2 Coupling model of a 2D plate: (a) complete model; (b) enlarged zone

Another role is to easily determine the transferred force values between FEM meshes and PD grids. In the proposed coupling model, the transferred force values between FEM meshes and PD grids can be determined directly through computing the interactions between true PD particles (red) with virtual PD particles (green). For simplicity, let fij denote the bond force vector that material particle

xj exerts on material particle xi . V represents the volume of any material particle. According to 9

the constitutive model of OSBPD in Section 2, the bond force that virtual PD particle true PD particle

xa exerts on

xb can be calculated as fba = T xb , t  xa  xb  T xa , t  xb  xa  V . And the

bond force that true PD particle

xb exerts on virtual PD particle

xa is equal to

fab   fba = T xa , t  xb  xa  T xb , t  xa  xb  V . It is noted that the interaction exerted on a PD particle is in the form of force density (N/m3) while the interaction exerted on a FEM node is in the form of force (N). Therefore, the bond force exerts on PD particle

f ba is directly the interaction that FEM node xa

xb , while the interaction that PD particle xb exerts on FEM node xa can be

obtained through multiplying the bond force

f ab and the volume of virtual particle xa , that is,

V  fab . Actually, the core theoretical basis of the proposed coupling model is the Newton’s third law, in which the action and reaction forces are reciprocal. To ensure that the determined coupling forces between FEM meshes and PD grids are reciprocal and balanced inherently, FEM nodes are directly regarded as virtual PD particles in the proposed coupling model. According to this theoretical basis, the coupling zone established in the proposed model has two roles. For PD particles (red), the coupling zone is utilized like a ‘fictitious material layer’ with depth 2 to apply external displacement (or force) constraints on the PD zone. On the other hand, For FEM nodes (blue and green), the coupling zone is like a fictitious boundary to which external forces are applied by the PD zone. 3.2 Numerical implementation Since there are no large-matrix operations in explicit methods, computation cost every time step is much less for explicit methods than for implicit methods. And the explicit methods are also easily implemented to solve the PD-FEM coupling models. Therefore, explicit methods are adopted in this study to solve the proposed coupling model for dynamic problems. For dynamic problems of short duration, many different explicit methods, such as the central-difference time integration method [13], the explicit Velocity-Verlet scheme [5, 28], and explicit forward and backward difference techniques [41], have been utilized in PD. In this study, explicit forward and backward difference techniques [41] are utilized in this work to solve the proposed coupling model for dynamic problems. let  be any bounded body, and let  be divided into two subdomains  FEM and  PD . If 10

the lumped mass matrix is adopted, the governing equations for FEM modeling can be decoupled. At any FEM node x in  FEM , the discretized form of motion equation at the nth time step can be expressed as, ( x , t n )  mx u



x H x  x  PD

V + F ext (t n )  F int (t n ) f xxFEM-PD 

 x  FEM  ,

(9)

where mx is the lumped mass, which can be obtained through operating the lumped mass matrixes of FEM elements pertaining to the FEM node.

f xxFEM-PD = T  x , t n  x   x  T  x , t n  x  x   V , 

which is the bond force on FEM node x (virtual PD particle) exerted by PD particles

x in PD ,

like fab in Fig. 2. f xxFEM-PD exists only when x   H x  x    PD is not empty. F ext (t n ) is the  external force exerted on FEM node x , and F int (t n ) is the internal nodal force resulting from the deformation of the FEM elements pertaining to FEM node x . The detailed calculation of these variables is described in Ref. [19]. For any PD particle x in  PD , the discretized form of motion equation can also be expressed as,



 u( x, t n )V 

x H x  x  FEM

f xxPD-FEM V + L( x , t n )V +b( x , t n )V 



= T x , t n where  is the density, f xxPD-FEM  on PD particle x exerted by FEM nodes



 x  PD 

,

(10)

x   x  T  x , t n  x  x   V , which is the bond force

x (virtual PD particles) in  FEM , like fba in Fig. 2.

f xxPD-FEM exists only when x   H x  x    FEM is not empty. L( x , t n ) is the accumulated bond  force

on

L  x, t n  =

PD

particle

 T x, t  n

xH x

x

exerted

by

other

PD

particles

within

its

horizon

as

x  x  T  x, t n  x  x  V .

After calculating the acceleration in Eqs. (9) and (10), the velocity and the displacement of the PD particle or the FEM node at the next time step can be obtained by employing explicit forward and backward difference techniques in two steps. Taking PD particle x as an example, the velocity and the displacement at the next time step can be calculated as, 11

( x, t n )  u ( x, t n ) u ( x, t n1 )  t  u

u( x, t n1 )  t  u ( x, t n1 )  u( x, t n )

,

(11)

in which, t is the stable time step which is described in detail in literature [13]. The above explicit forward and backward difference techniques in Eq. (11) can also be used to solve the dynamic solution of the FEM node. For static or quasi-static problems, the adaptive dynamic relaxation method (ADR) is commonly utilized in PD studies [34]. The ADR method is adopted in this work to obtain the quasi-static solution of the coupling model. In the ADR method, an adaptive damping coefficient is introduced into the governing equations to obtain the quasi-static steady solution. Given that the virtual mass matrix is diagonal in the ADR method, the governing equations for quasi-static problems can also be decoupled. The governing equations of the FEM node and the PD particle at the nth time step for quasi-static problems can be re-expressed as

 f u( x, t n )  c nf  f u ( x, t n )  F fn   p u( x, t n )+c np  p u ( x, t n )=Fpn =



x H x  x  PD



x H x  x  FEM

f xxFEM-PD V + F ext (t n )  F int (t n ) 

+ L( x , t n )+b( x , t n ) f xxPD-FEM 

 x  FEM  (12)

 x  PD 

in which,  p is the virtual density of the PD particle and  f is the virtual lumped mass of the FEM node.  p can be approximately set to be  3c  3x  for each PD particle, where c is the micro-



elastic modulus in BBPD. This micro-elastic modulus c can be calculated as c  9E h 3



plane stress case and c  48 E 5h 3





for 2D

for 2D plane strain case.  f can be determined based on

Greschgorin’ s theorem [38] through operating the element stiffness matrixes of FEM elements pertaining to the FEM node. c np and c nf are adaptive damping coefficients for PD particles and FEM nodes. The determination of these adaptive damping coefficients is explained in detail in Refs. [34, 41]. Taking PD particle x as an example, the velocity and the displacement at the (n+1/2)th time step can be calculated using the central difference technique as, u ( x , t n +1/2 )   2  c np t  u ( x , t n -1/2 )  2tFpn  p  u( x , t n +1 )  u( x , t n )  t  u ( x , t n +1/2 ) 12

 2  c t  . n p

(13)

The time step t in the ADR method is set to be 1s. In order to start the explicit central difference algorithm, u( x , t 0 )  0 and u ( x , t 0 )  0 are assumed in [38], which yields,

u ( x , t1/2 )  t  Fp0

 2  . p

(14)

The above central difference technique can also be utilized to solve quasi-static solution of the FEM node. The flowchart of the computational procedure of the proposed coupling model for quasi-static problems is shown in Fig. 3. Besides the above explicit solution methods for the proposed coupling model, the proof of the balance of linear momentum and angular momentum has also been derived in Appendix A.

13

Start Input computational parameters： Information of geometry model, neighbour information, maximum time step nmax , maximum load step kmax , material properties

Set k =0 Loop for load step : k  k 1 Set n =0

Loop for iteration step : n  n 1

Calculate coupling forces and internal nodal forces  for x   FEM 

Yes Remove the broken bond

s  sc

No

Calculate coupling forces and bond forces  for x   PD 

No Explicit time integration (ADR): Calculate new velocities

n  nmax

and displacements  for x   

Yes

Stop Satisfy the convergence criterion

No

Yes No

k  kmax

Yes

End Fig. 3 Flowchart of the computational procedure of the proposed coupling model

In all the following cases studied in Sections 4 and 5, the horizon size  in the PD zone and the 14

pure OSBPD is set to be 3x , and the effect of horizon size on the numerical results is discussed in Appendix B. The coupling zone is discretized using uniform rectangular FEM elements. The mesh size of the coupling zone is the same as the grid size of the PD zone.

4. Verification of the proposed coupling model “Ghost forces” are a well-known concept reported in some local-nonlocal coupling models [23, 29, 46], which result from fictitious boundary conditions when addressing the local/non-local interface in the coupling zone. “Ghost forces” may break the Newton’s third law, leading to the overall structural force non-equilibrium phenomenon. They can also lead to significant wave reflection at the coupling interface and deteriorate the numerical results [19, 29]. Before applying the coupling models to simulate different crack propagation problems, the accuracy of the PD-FEM coupling models needs to be verified and evaluated. In order to further check the accuracy of the proposed coupling model, three different numerical cases are performed, including a pre-notched beam under three-point bending (TPB), an elastic wave propagation case, and a 2D plate under uniform traction and pure shear conditions. 4.1 Quasi-static analysis of a three-point bending beam Some PD-FEM coupling models may suffer from the overall structural force non-equilibrium phenomenon, which can be easily revealed as the residual force in a static equilibrium structure, especially under an unsymmetrical loading condition. In order to verify the proposed coupling model, quasi-static analysis of a three-point bending beam is performed. The geometry and boundary condition are illustrated in Fig. 4(a). The elastic modulus E is 40GPa, and Poisson’s ratios v is 0.2. The thickness of the beam is 50mm, and the length of the center crack is 20mm. The 100mm off-center load is set to be 3000N. The coupling model is established in Fig. 4(b). The outermost part of the beam is discretized as FEM zone (in blue) using rectangular elements with an approximate grid size of 10mm. The middle part of the beam is PD zone (in red) with the uniform grid size x of 1mm. The whole coupling model consists of 3572 PD particles (in red), 2478 FEM nodes (in green and blue), 1834 rectangular FEM elements and 876 triangular FEM elements in total. The coupling model is assumed to be under plane stress condition and solved using the ADR method. The convergence criterion for this case is defined as

f n 1  f n  1.0 103 N , where f n1 and f n are the resultant force of the 15

reaction forces at the supports at (n  1)th and nth time iteration steps, respectively.

(a)

(b)

(c) Fig. 4 Description of the three-point bending beam: (a) geometry and boundary conditions (units: mm); (b) the proposed coupling model; (c) the FEM model

4.5 Reaction force at left support Reaction force at right support Resultant force

4.0 3.5

External force

Force (kN)

3.0 2.5 2.0 1.5 1.0 0.5 0.0 -0.5

0

100

200

300

400

500

Iteration steps (x100)

Fig. 5 Reaction forces predicted by the coupling model for the three-point bending beam

16

The predicted reaction forces by the coupling model for this case are shown in Fig. 5. The predicted reaction forces at left and right supports converge to 750.002N and 2249.995N, respectively. The resultant force of reaction forces is 2999.997N, which coincides with the applied external load (3000N). The residual force is not found in the established coupling model, which further indicates that the proposed coupling model is free of the overall structural force non-equilibrium phenomenon and satisfies the Newton’s third law. In addition, the FEM model is also established for comparison, as shown in Fig. 4(c). The FEM model is established through replacing the PD zone of the coupling model with uniform rectangular FEM elements. The stress fields around the crack tip predicted by the coupling model and FEM model are compared in Fig. 6 and Fig. 7. Overall, the profiles of the stress distributions predicted by the coupling model agree with the ones predicted by the FEM model. In Fig. 7, the stress  xx distribution in front of the crack tip predicted by the coupling model is very consistent with the FEM result, except the points near the crack tip. This is because that the stress singularity at the crack tip is relatively reduced in the PD analysis based on the nonlocal nature. Besides, the crack mouth opening

(a1)

(b1)

17

(a2)

(b2)

Fig. 6 Predicted stress fields around the crack tip: (a1)  xx and (a2)  yy by the coupling model; (b1)  xx and (b2)  yy by the FEM model 14 Coupling FEM

12 10

xx (MPa)

8 6 4 2 0 -2 -4

0

20

40

60

80

y (mm)

Fig. 7 Stress  xx distributions along the x-axis of a distance length in front of the crack tip

displacement (CMOD) is also obtained in both the two models. The obtained CMOD is equal to 5.5638e-6m and 5.7413e-6m, respectively for the coupling model and FEM model. The relative error is less than 3.09%, which may be caused by the surface effect. The surface effect cannot be completely removed in this case because no virtual PD particles are around crack surfaces. Overall, the difference between the coupling model and FEM model is acceptable. 4.2 Elastic wave propagation in a rectangular plate

18

(a)

(b) (c) Fig. 8 Description of the rectangular plate: (a) geometry and boundary conditions (units: m); (b) the proposed coupling model; (c) the FEM model

The second case is to check whether the proposed coupling model is free of significant wave reflection phenomenon through analyzing elastic wave propagation in a rectangular plate. The geometry and boundary conditions are shown in Fig. 8(a). The elastic modulus, density and Poisson’s ratio are 1Pa, 1kg/m3 and 0, respectively. Actually, this problem is a one-dimensional problem solved in 2D [47]. The initial displacement conditions can be expressed as,

u  x, y   0.0002  e

 100 x   , u x, y  0, x  0, 0.1 , y  0.02, 0.02 .      

The wave speed can be determined as

2

(15)

E   1m/s. The OSBPD-FEM coupling model for this

case is established as shown in Fig. 8(b). The middle part of the rectangular plate is discretized into the PD zone with red PD particles. In the PD zone, the uniform grid size x is 0.0005m. The outermost part is discretized with triangular FEM elements with blue FEM nodes. The whole coupling model consists of 4141 PD particles (in red), 4689 FEM nodes (in green and blue), 1540 rectangular FEM elements and 5870 triangular FEM elements in total. In addition, the FEM model is established for comparison, as shown in Fig. 8(c). The FEM model is established through replacing the PD zone of the coupling model with uniform rectangular FEM elements. Both two models are assumed to be under plane stress condition. The total solution time is set to be 0.2s, and the stable time step t is set to be 1.0  10 4 s. Both the coupling model and FEM model are solved by the explicit time 19

integration techniques mentioned in Section 3.

0.12

0.12 FEM Coupling

0.10

0.10 0.08

0.06

ux (mm)

ux (mm)

0.08

FEM Coupling

Incident wave

0.04 0.02

0.06

Reflected wave

0.04 0.02

0.00

0.00 Difference

Difference

-0.02 0.00

0.02

0.04

0.06

0.08

-0.02 0.00

0.10

x (m)

0.02

0.04

0.06

0.08

0.10

x (m)

(a) 300th time step (b) 1300th time step Fig. 9 Predicted displacements in x direction by two models: (a) at 300th time step; (b) at 1300th time step

(a) by the FEM model

(b) by the coupling model Fig. 10 Predicted displacement fields in x direction by the two models at 300th time step

The predicted displacements in x direction at the 300th and 1300th time step of both the coupling model and FEM model are scattered in Fig. 9. As observed in Fig. 6, the predicted wave profiles at both two time steps by the coupling model are in an excellent agreement with the ones by the FEM model. Compared with the similar results reported in Refs. [48] and [28], the results predicted by the proposed coupling model are smoother. The differences between the displacements in x direction predicted by the coupling model and the FEM model are also plotted in Fig. 9. The difference is defined 20

as u xc  u xf , where uxc and uxf represent the displacement in x direction predicted by the coupling model and the FEM model respectively. The difference fluctuates at the left coupling interface and the right coupling interface, but fortunately the difference relatively remains small. For 300th time step, the predicted displacement fields in x direction of both the two models are also compared in Fig. 10, which also shows a good agreement. This case indicates the proposed coupling model is free of significant wave reflection phenomenon at the interface between the PD zone and FEM zone. 4.3 Quasi-static analysis of a 2D plate

(a)

(b)

(c)

(d)

Fig. 11 Description of the 2D plate: geometry (units: mm) and pure traction (a) and pure shear (b); (c) the proposed coupling model; (d) the FEM model

To further evaluate the accuracy of the proposed coupling model, a 2D plate under pure traction and pure shear conditions is simulated as the third case. The geometry and boundary conditions are shown in Fig. 11(a) and (b). In Fig. 11(a), the plate is under pure traction, and the boundary conditions are u x  0.01mm at the left edge and u x  0.01mm at the right edge, respectively. In Fig. 11(b), the plate is under pure shear, and the boundary conditions are set as follows, 21

u x  0.02  y .  u y  0.02  x

(a)

(16)

(b)

(c) (d) Fig. 12 Predicted displacement fields in y direction: by the FEM model (a) and by the coupling model (b) under axial tension loading; vertical displacements by the FEM model (c) and by the coupling model (d) under simple shear loading

The elastic modulus E is 192GPa, and Poisson’s ratios v is 0.2. The proposed coupling model for this case is established, as shown in Fig. 11(c). The innermost part of the coupling model is discretized as the PD zone with red PD particles. The edge length of the PD zone is 15mm. In the PD zone, the uniform grid size x is 0.5mm. The outermost part is discretized using triangular FEM elements with blue FEM nodes for simplicity. The whole coupling model consists of 961 PD particles (in red), 9240 FEM nodes (in green and blue), 8976 rectangular FEM elements. In addition, a FEM model is also established in order to provide the approximate accurate solution for comparison, as 22

shown in Fig. 11(d). The FEM model is established through replacing the PD zone of the coupling model with uniform rectangular FEM elements. Both the coupling model and the FEM model are solved using the ADR method. The convergence criterion for this case is defined as th uAn 1  uAn  1.0 108 m , where uAn 1 and uAn are the displacement of point A at (n  1)th and n

time iteration steps, respectively.

0.010

0.002

FEM Coupling

0.001

uy (mm)

ux (mm)

0.005 0.000

0.000

-0.005

-0.001

-0.010

-0.002

-0.03

-0.02

-0.01

0.00

0.01

0.02

-0.03

0.03

-0.01

0.00

y (m)

(a)

(b)

0.01

0.02

0.03

0.01

0.02

0.03

0.6 FEM Coupling

0.4

FEM Coupling

0.4 0.2

uy (mm)

0.2

ux (mm)

-0.02

x (m)

0.6

0.0

0.0

-0.2

-0.2

-0.4

-0.4

-0.6 -0.03

FEM Coupling

-0.02

-0.01

0.00

0.01

0.02

-0.6 -0.03

0.03

-0.02

-0.01

0.00

y (m)

x (m)

(c)

(d)

Fig. 13 Displacement variations along centerlines: (a) ux (x, y = 0) and (b) uy (y, x = 0) under axial tension loading; (c) ux (y, x = 0) and (d) uy (x, y = 0) under simple shear loading

The displacement fields in y direction under pure traction and pure shear predicted by both two models are shown in Fig. 12. The displacement fields predicted by the proposed coupling model are in a perfect agreement with the counterparts predicted by the FEM model. In addition, the predicted displacement variation along centerlines under pure traction and pure shear by both the coupling model and the FEM model is compared in detail in Fig. 13, which shows a very good agreement. To 23

quantitatively evaluate the accuracy of the proposed coupling model, the relative error on the displacement along the horizontal centerline is shown in Fig 14. The relative error ex is defined as

u

xc

 u xf  u xf , in which uxc and uxf are the displacements in x direction predicted by the coupling

model and the FEM model, respectively. And the relative error ey is defined as

u

yc

 u yf  u yf , in

which uyc and uyf are the displacements in y direction predicted by the coupling model and the FEM model, respectively. It can be observed from Fig. 14 that the perturbations occur at the coupling boundaries. The absolute value of the maximum error under pure traction and pure shear is less than 2%, which is slightly less than the result predicted by the linear morphing function, but greater than the result predicted by the cubic morphing function [23]. Overall, these errors are relatively small and

2.5

1.5

2.0

1.0

1.5

0.5

Relative error ey (%)

Relative error ex (%)

acceptable.

1.0 0.5 0.0 -0.5 -1.0 -1.5 -0.03

0.0 -0.5 -1.0 -1.5 -2.0

-0.02

-0.01

0.00

0.01

0.02

-2.5 -0.03

0.03

x (m)

-0.02

-0.01

0.00

0.01

0.02

0.03

x (m)

(a)

(b)

Fig. 14 Relative errors on the displacement along centerlines: (a) ex (x, y = 0) under pure traction and (b) ey (x, y = 0) under pure shear

Besides the predicted displacements, the horizontal reaction forces under pure traction are also obtained through both two models, and the comparison is shown in Fig. 15. The reaction force predicted by the FEM model is 1920.025N, while the reaction force predicted by the coupling model converges to 1919.89N. The relative error is about 0.007%, which is negligible.

24

9

Reaction force (kN)

8

FEM Coupling

7 6 5 4 3 2 1 0

0

100

200

300

400

500

Iteration steps

Fig. 15 Predicted horizontal reaction forces along the left edge under axial tension loading

On the whole, based on the above three numerical cases, even though spurious coupling effects or errors still occur in the coupling interface, they are relatively small and negligible. It can therefore be concluded that the proposed coupling model is free of the overall structural force non-equilibrium phenomenon, significant wave reflection at the interface and ghost forces to some extent. Moreover, the coupling model can be utilized to accurately conduct structural analysis.

5. Numerical examples In this section, the proposed OSBPD-FEM coupling model is applied to study several typical crack propagation examples, including a quasi-static brittle crack propagation case and two dynamic crack propagation cases. 5.1 Splitting test of a 2D plate

(a) (b) Fig. 16 Description of splitting test: (a) the geometry and conditions (in mm); (b) the proposed coupling model

As shown in Fig. 16, we consider a splitting test on a pre-cracked specimen. The elastic modulus, 25

Poisson’s ratio and the fracture energy are 210GPa, 0.3 and 2700N/m, respectively. The loading procedure of this splitting test is controlled by crack mouth opening displacement (CMOD). In the coupling model, the middle part of the specimen is discretized into the PD zone with PD particles. In the PD zone, the uniform grid size x is 0.01mm. The outermost part is discretized by triangular FEM elements with FEM nodes. The discretization for the coupling model is shown in Fig. 16(b). The coupling model consists of 2020 PD particles (in red), 2668 FEM nodes (in green and blue), 1000 rectangular FEM elements and 3024 triangular FEM elements in total. In addition, a pure OSBPD model is established for comparison. The pure OSBPD model consists of 10000 PD particles with a uniform grid size of 0.01mm. Both the coupling models and the pure OSBPD model are assumed to be under plane stress conditions. The initial displacement is 0.015mm, and then an increment of 0.0002mm is applied step by step. 150 load steps are considered in total. Both the coupling model and the pure PD model are solved using the ADR method. The convergence criterion for this case is defined as

f n 1  f n  1.0 102 N , where

f n1 and

f n are the reaction force at the loading point at

(n  1)th and nth time iteration steps, respectively.

50

OSBPD Coupling

Load (N)

40 30 20 10 0 0.00

0.01

0.02

0.03

0.04

0.05

CMOD (mm)

Fig. 17 Predicted load-CMOD curves

The predicted load-CMOD curves by the coupling model and the pure OSBPD model are compared in Fig. 17, which shows a good agreement even though the post-peak part of the predicted load-CMOD curves is not smooth. A smaller displacement increment may be adopted in order to obtain a smoother load-CMOD curve. The damage maps at the 60th and 150th loading steps obtained by the coupling model and the pure OSBPD model are compared in Fig. 18, which also show an excellent 26

agreement.

(a) 60th step

(b) 150th step

(c) 60th step

(d) 150th step

Fig. 18 Damage maps at two steps: (a) and (b) by the coupling model; (c) and (d) by the pure OSBPD model

5.2 Traction of a pre-cracked plate

(a)

27

(b) Fig. 19 Description of the pre-crack plate: (a) the geometry and conditions (in mm); (b) the coupling model

A crack branching in a pre-cracked plate with an initial crack is studied in the second example. The geometry and boundary conditions are shown in Fig. 19(a). The upper and lower edges are subject to a uniform traction of  =12MPa . This case has been simulated using BBPD [5, 49] and using XFEM [50]. The plate is made of soda-lime glass, and the elastic modulus, Poisson’s ratio, density, and the fracture energy are 72GPa, 0.22, 2440kg/m3 and 135N/m, respectively. To simulate this case, both the OSBPD-FEM coupling model and the pure OSBPD model are established. The coupling model is shown in Fig. 19(b). The domain where crack will branch is discretized into the PD zone with red PD particles. In the PD zone, the uniform grid size x is 0.5mm. The rest part is non-uniformly meshed with triangular FEM elements. The whole coupling model consists of 6890 PD particles (in red), 5042 FEM nodes (in green and blue), 3120 rectangular FEM elements and 3070 triangular FEM elements. The pure OSBPD model consists of 16100 PD particles with a uniform grid size of 0.5mm, and the horizon size  is also set to be 3x . In the two models, the time step is set to be 0.02 s , and the total time is 42 s with 2100 steps.

(a) 1200th step

28

(b) 1200th step

(c) 2100th step

(d) 2100th step Fig. 20 Damage maps at two time steps: (a) and (c) by the coupling model; (b) and (d) by the OSBPD model

Crack propagation process can be observed through damage maps. Fig. 20 shows the damage maps obtained by the two models at the 1200th and 2100th steps. The crack paths predicted by the two models are compared in Fig. 21. It can be seen that the crack paths predicted by the coupling model are in an excellent agreement with the ones predicted by the pure OSBPD model. Apart from predicting crack paths, PD can also approximately estimate the crack propagation speed. As proposed in [49], the crack tip in PD can be identified at each time step by the farthest particle with a damage level above 0.35. The crack speed can be approximately calculated as n 1 n V n 1/2  xtip  xtip

t

n 1

 tn  ,

(17)

n 1 n where xtip and x tip are the positions of the crack tip at the time steps t n 1 and t n , respectively. In

this case, the crack speed is calculated by considering a time interval of 2 s , or per 100 time steps. The calculated crack speeds obtained by the two models are compared in Fig. 22. The predicted trend 29

of the crack speed is in a good agreement with the ones reported in Refs. [28] and [49]. It is noted that the trend of crack speed before branching experimentally observed [51] is successfully captured in both the two models here. The maximum crack speed is reached before crack branching at t  22μs , which is close to the maximum crack speed experimentally-measured in Ref. [52].

0.020 0.015

OSBPD Coupling

0.010

y (m)

0.005 0.000

-0.005 -0.010 -0.015 -0.020

0.05

0.06

0.07

0.08

0.09

0.10

x (m)

Crack propagation speed (m/s)

Fig. 21 Numerical crack paths of the pre-cracked plate predicted by the two models

2000 1600 1200 800 Coupling OSBPD

400 0 0

5

10

15

20

25

30

35

40

45

Time (s)

Fig. 22 Comparison of the crack propagation speed calculated by the two models (1580m/s is the maximum speed measured experimentally in [52])

5.3 Kalthoff-Winkler experiment

30

Fig. 23 Description of the Kalthoff-Winkler experiment (units: mm)

An important study of dynamic crack propagation is the Kalthoff-Winkler experiment [53], which has been simulated by different numerical methods [28, 50]. The geometry and loading setup are described in Fig. 23. It was observed in the experiment that at lower impact velocities, crack will propagate forward along an angle of about 70°. In this section, a low impact velocity v0 = 22m/s is considered. The plate is made of soda-lime glass, and the elastic modulus, Poisson’s ratio, density, and the fracture energy are 190GPa, 0.25, 7800kg/m3 and 69000N/m, respectively. Because the loading and geometry of the experiment are symmetric, as shown in Fig. 24(a), only the lower part of the plate is used to establish numerical models. For this case, both the OSBPD-FEM coupling model and the pure OSBPD model are established. The coupling model is shown in Fig. 24(b). The domain where crack will propagate is discretized into the PD zone with red PD particles. In the PD zone, the uniform grid size x is 0.5mm. The rest part is non-uniformly meshed by triangular FEM elements and rectangular FEM elements. The whole coupling model consists of 15831 PD particles (in red), 13175 FEM nodes (in green and blue), 7434 rectangular FEM elements and 10490 triangular FEM elements. The pure OSBPD model consists of 40301 PD particles with a uniform grid size of 0.5mm, and the horizon size  is also set to be 3x . In the two models, the time step is set to be 0.02 s , and the total time is 88 s with 4400 steps.

31

(a) (b) Fig. 24 Description of the Kalthoff-Winkler experiment: (a) the geometry and boundary conditions of the plate (units: mm); (b) the coupling model

Fig. 25 shows the damage maps obtained by the two models at the 2400th and 4400th steps. The crack paths predicted by different numerical methods are compared with the experimental one in Fig. 26. It can be seen that the crack paths predicted by the coupling model are in an excellent agreement with the ones predicted by the pure OSBPD model. The predicted crack propagation angle is slightly less than the experimental one, but better than the one predicted by XFEM [50]. In addition, according to Eq. (16), the crack speeds obtained by the two models are compared in Fig. 27. In this case, the crack speed is calculated by considering a time interval of 2 s , or per 100 time steps. It was stated in [54] that the crack speed cannot be higher than the Rayleigh speed. The Rayleigh speed for this case is equal to 2871.8m/s. The maximum predicted crack speeds are 1684m/s and 1587m/s predicted by the coupling model and the pure OSBPD model, respectively, which are less than the limit of 75% of the Rayleigh speed.

32

(a) 2400th step

(b) 2400th step

(c) 4400th step

(d) 4400th step

Fig. 25 Damage maps at two steps: (b) and (d) by the coupling model; (a) and (c) by the pure OSBPD model

0.10 Coupling OSBPD XFEM Experiment

y (m)

0.08 0.06 0.04 0.02 0.00 0.00

0.02

0.04

0.06

0.08

0.10

x (m)

Fig. 26 Experimental and numerical crack paths of the Kalthoff-Winkler experiment

33

Crack propagation speed (m/s)

3000 Rayleigh Coupling OSBPD

2500 2000 1500 1000 500 0 0

10

20

30

40

50

60

70

80

90

Time (s)

Fig. 27 Comparison of the crack propagation speed calculated by the two models (2871.8m/s is the Rayleigh speed)

6. Conclusion In this work, a new simple coupling model of FEM and OSBPD is proposed to accurately predict dynamic and quasi-static brittle fracture. The accuracy of the proposed coupling model has been verified by several linear elastic cases and brittle fracture cases under dynamic and quasi-static loading. Some key remarks can be made as follows: (1) In the proposed coupling model, PD grids are directly coupled with FEM meshes without a overlapped zone. The transferred force values can be easily determined through computing the interactions between virtual PD particles in the coupling zone and true PD particles in the PD zone. Hence, the transferred forces between PD grids and FEM meshes are reciprocal and balanced inherently. The proposed coupling model avoids the complexity of the transitional zone-based coupling models and also the non-equilibrium phenomenon of the structural force reported in the sharing nodes-based coupling model. (2) The proposed coupling model approximately removes the surface effect of PD zone, and at the same time, it is free of significant influence of ghost forces and significant wave reflection at the interface between PD grids and FEM meshes. (3) The coupling model takes advantage of the generality of the OSBPD, and hence it can be utilized to describe more general materials without the limitation of Poisson’ ratio.

Acknowledgement The financial supports from the Science and Technology Innovation Commission of Shenzhen 34

Municipality (Ref: JCYJ20160229165310679) and the Research Grants Council of the Hong Kong Special Administrative Region, China (Ref: CityU 9680086) are gratefully acknowledged.

Appendix A A.1 Balance of linear momentum To satisfy the balance of linear momentum for the bounded body  described above, it is required that



x FEM



( x , t )-F ext (t )  +  mx u

 

x FEM x H x  x  PD

V + f xxFEM-PD 

   u( x, t )V -b( x, t )V 

x PD

 

x PD x H x  x  FEM

V + f xxPD-FEM 



L( x , t )V 

x PD



F int (t ) .

(A.1)

x FEM

=0

Proof Due to the symmetrical characteristic of stiffness matrix of FEM elements,



F int (t ) =0 is

x FEM



automatically satisfied.

L( x , t )V =

 

x PD x H x

x PD

T  x , t  x   x  T  x , t  x  x   V V =0

is

also automatically satisfied for ordinary state-based peridynamics, as derived in [41]. Hence, only

 

V + f xxFEM-PD 

 

V + f xxFEM-PD 

x FEM x H x  x  PD

x FEM x H x  x  PD

 

=

x FEM x H x  x  PD

=

 

V =0 is needed to be proofed. f xxPD-FEM 

 

V f xxPD-FEM 

x PD x H x  x  FEM

x PD x H x  x  FEM

T  x , t  x   x  T  x , t  x  x   V V +

  T  x, t 

x   x V V -

x FEM x H x  x  PD

+

  T  x, t 

  T  x, t 

 

x PD x H x  x  FEM

T  x , t  x   x  T  x , t  x  x   V V

x  x  V V

x FEM x H x  x  PD

x   x V V -

x PD x H x  x  FEM

  T  x, t 

x  x  V V

x PD x H x  x  FEM

(A.2) In Eq. (A.2), if the parameters x and

x in the second and fourth terms are exchanged, the second

and fourth terms become 35

  T  x, t 

x  x  V V =

  T  x, t 

x  x  V V =

x FEM x H x  x  PD

  T  x, t 

x   x V V ,

(A.3)

  T  x, t 

x   x V V .

(A.4)

x PD x H x  x  FEM

x PD x H x  x  FEM

x FEM x H x  x  PD

 

Combining Eqs. (A.2), (A.3) and (A.4),

x FEM x H x  x  PD

V + f xxFEM-PD 

 

x PD x H x  x  FEM

V =0 can hold. f xxPD-FEM 

And hence, the balance of linear momentum is satisfied for the proposed FEM-PD coupling system. A.2 Balance of angular momentum To satisfy the balance of angular momentum for the bounded body  described above, it is required that



x FEM

( x , t )-F ext (t )   y +  mx u

 



x FEM x H x  x  PD

 yV + f xxFEM-PD 

   u( x, t )V -b( x, t )V   y

x PD

 

x PD x H x  x  FEM

 yV + f xxPD-FEM 



L( x , t )  yV 

x PD



F int (t )  y (A.5)

x FEM

=0

Proof



L( x , t )  yV =

 

x PD x H x

x PD

T  x , t  x   x  T  x , t  x  x    yV V =0 is also satisfied

for ordinary state-based peridynamics, as derived in [41]. Also, due to the fact that FEM is based on the classical continuum mechanics, the symmetry of stress tensor naturally ensure the balance of angular momentum. Hence, the summation of the angular momentum generated by all the internal nodal forces must vanish, that is



F int (t )  y =0 .

x FEM

 

Hence, only

x FEM x H x  x  PD

 

x FEM x H x  x  PD

=

f xxFEM-PD  yV + 

f xxFEM-PD  yV + 

 

x PD x H x  x  FEM

 

x PD x H x  x  FEM

f xxPD-FEM  yV 

  T  x, t 

x   x  yV V -

  T  x, t 

x   x  yV V -

x FEM x H x  x  PD

+

f xxPD-FEM  yV =0 is needed to be proofed. 

  T  x, t 

x  x   yV V

  T  x, t 

x  x   yV V

(A.6)

x FEM x H x  x  PD

x PD x H x  x  FEM

In Eq. (A.6), if the parameters x and

x PD x H x  x  FEM

x in the second and fourth terms are exchanged, the second 36

and fourth terms become

  T  x, t 

x  x   yV V =

  T  x, t 

x  x   yV V =

x FEM x H x  x  PD

  T  x, t 

x   x  yV V ,

(A.7)

  T  x, t 

x   x  yV V .

(A.8)

x PD x H x  x  FEM

x PD x H x  x  FEM

x FEM x H x  x  PD

And then, Eq. (A.6) can be re-expressed as

 

x FEM x H x  x  PD

=

f xxFEM-PD  yV + 

  T  x, t 

 

x PD x H x  x  FEM

f xxPD-FEM  yV 

x   x  ( y  y)V V +

x FEM x H x  x  PD

  T  x, t 

x   x  ( y  y)V V

(A.9)

x PD x H x  x  FEM

For ordinary state-based peridynamics, Eq. (A.9) is automatically satisfied, and hence the balance of angular momentum can hold.

Appendix B The horizon size of PD is generally set as   mx , and the size of the coupling zone in the proposed coupling strategy is equal to 2 . The effect of change in horizon size on the numerical results needs to be considered. The 2D plate under pure tension in Section 4.3 is used to study the effect of change in horizon size. The geometry and boundary conditions are shown in Fig. 11(a). In this study, the value of m is chosen as 3, 4 and 5, respectively. When the value of m is fixed, we consider the effect of horizon size

 by changing grid size x . At the same time, the results predicted by the FEM model are adopted as the accurate solution to compare with the ones predicted by the coupling model. Fig. B.1 shows the comparisons of displacement variations along vertical centerline predicted by the FEM model and the coupling models. It can be seen from Fig. B.1 that when the value of m is fixed, with the decrease of the horizon size, the spurious coupling effect at the coupling interface is reduced gradually, and the displacement distribution predicted by the coupling model also gradually converges to the accurate FEM solution. This study indicates that an appropriate horizon size is needed to suppress the coupling effect at the coupling interface. Under the consideration of both the accuracy and computational cost, “m=3” is recommended in the proposed coupling strategy.

37

FEM x=0.5mm,  =1.5mm, m=3 x=1.0mm,  =3.0mm, m=3 x=1.5mm,  =4.5mm, m=3

0.002

0.001

uy (mm)

uy (mm)

0.001 0.000

0.000

-0.001

-0.001

-0.002

-0.002

-0.03

-0.02

-0.01

0.00

0.01

0.02

FEM x=0.5mm,  =2.0mm, m=4 x=1.0mm,  =4.0mm, m=4 x=1.5mm,  =6.0mm, m=4

0.002

0.03

-0.03

-0.02

-0.01

y (m)

0.01

0.02

0.03

y (m)

(a)

(b) FEM x=0.25mm,  =1.25mm, m=5 x=0.5mm,  =2.5mm, m=5 x=1.0mm,  =5.0mm, m=5

0.002 0.001

uy (mm)

0.00

0.000

-0.001 -0.002 -0.03

-0.02

-0.01

0.00

0.01

0.02

0.03

y (m)

(c) Fig. B.1 Displacement variations along centerline uy (y, x = 0) under axial tension loading: (a) m=3; (b) m=4 and (c) m=5

References [1] Silling SA. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids 2000;48(1):175-209. [2] Ferté G, Massin P, Moës N. 3D crack propagation with cohesive elements in the extended finite element method. Computer Methods in Applied Mechanics and Engineering 2016;300:347-374. [3] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics 2002;69(7):813-833. [4] Fan HF, Bergel GL, Li SF. A hybrid peridynamics–SPH simulation of soil fragmentation by blast loads of buried explosive. International Journal of Impact Engineering 2016;87:14-27. [5] Ha YD, Bobaru F. Studies of dynamic crack propagation and crack branching with peridynamics. International Journal of Fracture 2010;162(1-2):229-244. [6] Kuhl E, Ramm E, De Borst R. An anisotropic gradient damage model for quasi-brittle materials. Computer Methods in Applied Mechanics and Engineering 2000;183(1-2):87-103. [7] Nguyen VP, Wu JY. Modeling dynamic fracture of solids with a phase-field regularized cohesive zone model. Computer Methods in Applied Mechanics and Engineering 2018;340:1000-1022. [8] Cuvilliez S, Feyel F, Lorentz E, Michel-Ponnelle S. A finite element approach coupling a continuous gradient damage 38

model and a cohesive zone model within the framework of quasi-brittle failure. Computer Methods in Applied Mechanics and Engineering 2012;237:244-259. [9] Yang D, He XQ, Yi SH, Liu XF. An improved ordinary state-based peridynamic model for cohesive crack growth in quasi-brittle materials. International Journal of Mechanical Sciences 2019;153:402-415. [10] Yang D, Dong W, Liu XF, Yi SH, He XQ. Investigation on mode-I crack propagation in concrete using bond-based peridynamics with a new damage model. Engineering Fracture Mechanics 2018;199:567-581. [11] Ouchi H, Katiyar A, York J, Foster JT, Sharma MM. A fully coupled porous flow and geomechanics model for fluid driven cracks: a peridynamics approach. Computational Mechanics 2015;55(3):561-576. [12] Nadimi S, Miscovic I, McLennan J. A 3D peridynamic simulation of hydraulic fracture process in a heterogeneous medium. Journal of Petroleum Science and Engineering 2016;145:444-452. [13] Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Computers & Structures 2005;83(17-18):1526-1535. [14] Le QV, Chan WK, Schwartz J. A two‐dimensional ordinary, state‐based peridynamic model for linearly elastic solids. International Journal for Numerical Methods in Engineering 2014;98(8):547-561. [15] Silling SA, Epton M, Weckner O, Xu J, Askari E. Peridynamic states and constitutive modeling. Journal of Elasticity 2007;88(2):151-184. [16] Wu CT, Ren B. A stabilized non-ordinary state-based peridynamics for the nonlocal ductile material failure analysis in metal machining process. Computer Methods in Applied Mechanics and Engineering 2015;291:197-215. [17] Javili A, Morasata R, Oterkus E, Oterkus S. Peridynamics review. Mathematics and Mechanics of Solids 2018:1081286518803411. [18] Macek RW, Silling SA. Peridynamics via finite element analysis. Finite Elements in Analysis and Design 2007;43(15):1169-1178. [19] Kilic B, Madenci E. Coupling of peridynamic theory and the finite element method. Journal of Mechanics of Materials and Structures 2010;5(5):707-733. [20] Agwai A, Guven I, Madenci E, Damage prediction for electronic package drop test using finite element method and peridynamic theory, in: 2009 59th Electronic Components and Technology Conference, IEEE, 2009, pp. 565-569. [21] Liu WY, Hong J-W. A coupling approach of discretized peridynamics with finite element method. Computer Methods in Applied Mechanics and Engineering 2012;245:163-175. [22] Lubineau G, Azdoud Y, Han F, Rey C, Askari A. A morphing strategy to couple non-local to local continuum mechanics. Journal of the Mechanics and Physics of Solids 2012;60(6):1088-1102. [23] Han F, Lubineau G, Azdoud Y, Askari A. A morphing approach to couple state-based peridynamics with classical continuum mechanics. Computer methods in applied mechanics and engineering 2016;301:336-358. [24] Seleson P, Beneddine S, Prudhomme S. A force-based coupling scheme for peridynamics and classical elasticity. Computational Materials Science 2013;66:34-49. [25] Seleson P, Ha YD, Beneddine S. Concurrent coupling of bond-based peridynamics and the Navier equation of classical elasticity by blending. International Journal for Multiscale Computational Engineering 2015;13(2). [26] Galvanetto U, Mudric T, Shojaei A, Zaccariotto M. An effective way to couple FEM meshes and Peridynamics grids for the solution of static equilibrium problems. Mechanics Research Communications 2016;76:41-47. [27] Zaccariotto M, Mudric T, Tomasi D, Shojaei A, Galvanetto U. Coupling of FEM meshes with Peridynamic grids. Computer Methods in Applied Mechanics and Engineering 2018;330:471-497. [28] Bie YH, Cui XY, Li ZC. A coupling approach of state-based peridynamics with node-based smoothed finite element method. Computer Methods in Applied Mechanics and Engineering 2018;331:675-700. [29] Giannakeas IN, Papathanasiou TK, Bahai H. Wave reflection and cut‐off frequencies in coupled FE‐Peridynamic grids. International Journal for Numerical Methods in Engineering 2019:1-27. 39

[30] Ni T, Zaccariotto M, Zhu QZ, Galvanetto U. Static solution of crack propagation problems in Peridynamics. Computer Methods in Applied Mechanics and Engineering 2019;346:126-151. [31] Sun W, Fish J. Superposition-based coupling of peridynamics and finite element method. Computational Mechanics 2019:1-18. [32] Wang XN, Kulkarni SS, Tabarraei A. Concurrent coupling of peridynamics and classical elasticity for elastodynamic problems. Computer Methods in Applied Mechanics and Engineering 2019;344:251-275. [33] Huang XH, Bie ZW, Wang LF, Jin YL, Liu XF, Su GS, He XQ. Finite element method of bond-based peridynamics and its ABAQUS implementation. Engineering Fracture Mechanics 2019;206:408-426. [34] Kilic B, Madenci E. An adaptive dynamic relaxation method for quasi-static simulations using the peridynamic theory. Theoretical and Applied Fracture Mechanics 2010;53(3):194-204. [35] Dokainish M, Subbaraj K. A survey of direct time-integration methods in computational structural dynamics—I. Explicit methods. Computers & Structures 1989;32(6):1371-1386. [36] Fan H, Bergel GL, Li S. A hybrid peridynamics–SPH simulation of soil fragmentation by blast loads of buried explosive. International Journal of Impact Engineering 2016;87:14-27. [37] Diyaroglu C, Oterkus E, Madenci E, Rabczuk T, Siddiq A. Peridynamic modeling of composite laminates under explosive loading. Composite Structures 2016;144:14-23. [38] Underwood P. Dynamic Relaxation, in:

Computational Methods for Transient Analysis. North-Holland; 1986.

[39] Zhang H, Qiao PZ. A state-based peridynamic model for quantitative fracture analysis. International Journal of Fracture 2018;211(1-2):217-235. [40] Luo J, Ramazani A, Sundararaghavan V. Simulation of micro-scale shear bands using peridynamics with an adaptive dynamic relaxation method. International Journal of Solids and Structures 2018;130:36-48. [41] Madenci E, Oterkus E. Peridynamic Theory, in:

Peridynamic Theory and Its Applications. Springer; 2014.

[42] Madenci E. Peridynamic integrals for strain invariants of homogeneous deformation. ZAMM‐Journal of Applied Mathematics and Mechanics 2017;97(10):1236-1251. [43] Silling SA. Linearized theory of peridynamic states. Journal of Elasticity 2010;99(1):85-111. [44] Le QV, Bobaru F. Surface corrections for peridynamic models in elasticity and fracture. Computational Mechanics 2018:1-20. [45] Foster JT, Silling SA, Chen WN. An energy based failure criterion for use with peridynamic states. International Journal for Multiscale Computational Engineering 2011;9(6). [46] Badia S, Parks M, Bochev P, Gunzburger M, Lehoucq R. On atomistic-to-continuum coupling by blending. Multiscale Modeling & Simulation 2008;7(1):381-406. [47] Ren HL, Zhuang XY, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering 2017;318:762-782. [48] Ren HL, Zhuang XY, Cai YC, Rabczuk T. Dual‐horizon peridynamics. International Journal for Numerical Methods in Engineering 2016;108(12):1451-1476. [49] Dipasquale D, Zaccariotto M, Galvanetto U. Crack propagation with adaptive grid refinement in 2D peridynamics. International Journal of Fracture 2014;190(1-2):1-22. [50] Belytschko T, Chen H, Xu JX, Zi G. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. International Journal for Numerical Methods in Engineering 2003;58(12):1873-1905. [51] Field JE. Brittle fracture: its study and application. Contemporary Physics 1971;12(1):1-31. [52] Bowden FP, Brunton JH, Field JE, Heyes AD. Controlled fracture of brittle solids and interruption of electrical current. Nature 1967;216(5110):38. [53] Kalthoff JF, Winkler S. Failure mode transition at high rates of shear loading. DGM Informationsgesellschaft mbH, Impact Loading and Dynamic Behavior of Materials 1988;1:185-195. 40

[54] Ravi-Chandar K. Dynamic fracture of nominally brittle materials. International Journal of Fracture 1998;90(1-2):83102.

41

Highlights 

A simple coupling model of FEM and OSBPD is proposed



The proposed coupling model is free of significant influence of ghost forces, significant wave reflection and non-equilibrium phenomenon of the overall structural force



The surface effect of OSBPD can be approximately removed in the proposed coupling model

42

Declaration of interests The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

43

Credit Author Statement:

Dong Yang: Conceptualization， Methodology, Formal analysis, Data curation, Computation, Writing – Original Draft. Xiaoqiao He: Formal analysis, Conceptualization, Project Administration, Supervision, Funding Acquisition. Shenghui Yi: Formal analysis, Resources, Validation, Writing- Reviewing and Editing, Funding Acquisition. Yajie Deng: Formal analysis, Software, Validation, Writing- Reviewing and Editing. Xuefeng Liu: Formal analysis, Software, Validation, Writing- Reviewing and Editing.

44