A Hybrid Local/Nonlocal Continuum Mechanics Modeling and Simulation of Fracture in Brittle Materials

Classical continuum mechanics which leads to a local continuum model, encounters challenges when the discontinuity appears, while peridynamics that falls into the category of nonlocal continuum mechanics suffers from a high computational cost. A hybrid model coupling classical continuum mechanics with peridynamics can avoid both disadvantages. This paper describes the hybrid model and its adaptive coupling approach which dynamically updates the coupling domains according to crack propagations for brittle materials. Then this hybrid local/nonlocal continuum model is applied to fracture simulation. Some numerical examples like a plate with a hole, Brazilian disk, notched plate and beam, are performed for verification and validation. In addition, a peridynamic software is introduced, which was recently developed for the simulation of the hybrid local/nonlocal continuum model.


Introduction
The accurate simulation of discontinuities-induced failure is still a challenge. In order to understand the mechanism of fracture, Griffith [Griffith (1921)] proposed classical fracture mechanics in the 1920s for brittle materials. It is widely used to predict the mechanical behavior of structures within the framework of local continuum mechanics. However, classical fracture mechanics is performed based on the fact that the location of a preexisting crack is known. Furthermore, it is difficult to predict the crack nucleation.
To overcome the mentioned shortcomings in classical fracture mechanics, Francfort et al. [Francfort and Marigo (1998)] developed a phase-field model in 1998. Bourdin et al. and Borden et al. [Bourdin, Francfort and Marigo (2000); Borden, Verhoosel, Scott et al. (2012)] applied the phase-field model to fracture simulation with various computational methods. Zhou et al. [Zhou, Zhuang, Zhu et al. (2018) ;; Zhou and Xia (2018) ;Zhou, Zhuang and Rabczuk (2019)] studied the crack propagation, branching and coalescence with the phase-field model. However, the phase-field method does not model the development of discontinuities, but the propagation of highly localized damage, as it is still fundamentally a continuum-field based technique. By extending the local continuum framework, the extended finite element method (XFEM) is able to explicitly introduce discontinuities [Belytschko and Black (1999)]. Some researchers applied XFEM to crack propagations. However, the limitation will appear when XFEM encounters complex phenomena such as crack branching. Cundall et al. [Cundall and Strack (1979)] proposed the discrete element method (DEM), which is a mesh-free method, to analyze the crack propagation, branching and coalescence. However, fracture, obtained by DEM, has the size effect and it closely related to the size of particles.
In 2000, Silling [Silling (2000)] proposed a novel model, named peridynamics, based on integral equilibrium equations. Since integral equations are mathematically compatible with discontinuities, peridynamics can handle the crack initiation and propagation without introducing a complicated failure criterion [Kilic, Agwai and Madenci (2009)], and Ha et al. [Ha and Bobaru (2010); Wang, Oterkus and Oterkus (2018); Shou, Zhou and Berto (2019); Ren, Zhuang and Rabczuk (2016); Rabczuk, Ren and Zhuang (2019)] have successfully applied it to fracture simulations. However, peridynamics suffers from a high computational cost because a point in peridynamics is interacting with numerous points in a finite neighborhood, which results in the expensive computation of the resultant of forces.
To overcome the disadvantage of the expensive computation, some improvements have been developed. Ren et al. [Ren, Zhuang and Rabczuk (2017)] proposed a dual-horizon peridynamic model in which the variable horizon is used instead of a constant horizon.
Compared with other coupling methods, the morphing approach, proposed by Lubineau et al., established a constraint equation of stiffness which naturally guaranteed the equivalence of material parameters between peridynamics and classical continuum mechanics. The morphing technique is also easy to implement from a technical point of view because it fully relies on spatial modifications of the material parameters, which is easy to control.
Based on the morphing approach, a hybrid local/nonlocal continuum model was developed to couple bond-based peridynamics with classical continuum mechanics [Lubineau, Azdoud, Han et al. (2012)]. The hybrid model makes use of a priori function to indicate the nonlocal subdomain for peridynamics. This nonlocal subdomain is in charge of crack initiations and propagations. The hybrid model is compatible with the traditional boundary conditions, for instance, surface tractions. This hybrid model is also used to couple ordinary state-based peridynamics with classical continuum mechanics ]. Additionally, Han et al. ] also coupled bond-based peridynamics with damage mechanics by this hybrid model.
In 2014, Azdoud et al. [Azdoud, Han and Lubineau (2014)] presented an approach to adaptively update the peridynamic subdomain according to the evolution of a crack with a three-dimensional (3D) benchmark example. Based on it, this paper further describes this adaptive coupling algorithm for the hybrid model in detail. Then a validation example under a non-homogeneous deformation is investigated for the effectiveness of the hybrid model. After the validation, we simulate several two-dimensional (2D) numerical examples to display its efficiency and performance for fracture simulations. In addition, a peridynamic software was introduced. The software was recently developed to apply the hybrid model with the adaptive coupling approach to fracture simulations.
The remainder of this paper is organized as follows: Section 2 reviews bond-based peridynamics and the hybrid local/nonlocal continuum model; Section 3 describes a process of applying the hybrid model to the adaptive simulation of crack propagations; Section 4 introduces a software which was developed to implement the hybrid model, and validates the hybrid model with this software in a bi-dimensional example under non-homogeneous deformation; and some benchmark numerical examples are displayed to predict crack initiations and propagations in Section 5.

Bond-based peridynamics
In 2000, Silling [Silling (2000)] proposed peridynamics which is one kind of nonlocal model. In peridynamics, the equilibrium equation at point x for a quasi-static problem is expressed as where H δ (x) is the neighborhood of point x and δ is the horizon. b is the external body force. f (p → x) is the bond force such that wheref is the partial interaction due to the action of point p over point x (with respect to point x over point p) [Lubineau, Azdoud, Han et al. (2012)]. The bond force is considered as the superposition of two interactions (see Fig. 1). For linear elasticity and small deformations, a possible constitutive model [Silling (2000); Silling and Lehoucq (2010) where the relative position vector ξ = p − x is called bond. η is the relative displacement vector and C[x](ξ) is the micromodulus tensor, which are defined as follows [Silling (2000)], where c[x](ξ) is the coefficient function, and u is the displacement.

Bond failure
In peridynamics, when the bond stretch s is larger than a critical value, s crit , the bond will break in an irreversible manner [Silling and Askari (2005)]. This failure law has been widely used for fracture simulations [Ha and Bobaru (2010); Azdoud, Han and Lubineau (2014); Hu, Ha and Bobaru (2012)].
The failure law is implemented by introducing a history-dependent scalar-valued function µ(t, ξ) into the bond force [Silling and Askari (2005)]. µ(t, ξ) is defined as follows, where t and t denote the computational steps. Note that the critical bond stretch, s crit , is considered as an intrinsic material parameter.

A hybrid local/nonlocal continuum model
Ωm Figure 2: The subdomains Ω 1 , Ω 2 and Ω m One disadvantage of peridynamics is the high computational cost. In order to overcome this issue, Lubineau et al. [Lubineau, Azdoud, Han et al. (2012)] developed a hybrid model to couple peridynamics with classical continuum mechanics. In the hybrid model, let the domain Ω be the overall domain. It can be divided into three subdomains: Ω 1 , Ω 2 and Ω m [Lubineau, Azdoud, Han et al. (2012)]. Ω = Ω 1 ∪ Ω 2 ∪ Ω m , Ω 1 ∩ Ω 2 = ∅, Ω 1 ∩ Ω m = ∅ and Ω 2 ∩ Ω m = ∅ (see Fig. 2). The idea of the hybrid model is that there is only classical continuum mechanics in the subdomain Ω 1 while there is only peridynamics in the subdomain Ω 2 . However, both classical continuum mechanics and peridynamics are defined in the morphing subdomain Ω m . The displacement u is applied on the boundary Γ u of ∂Ω, and the external surface traction T is imposed on the boundary Γ T of ∂Ω. n is the outward direction that is normal to Γ T .
In this paper, we focus on homogeneous, isotropic and linearly elastic materials under small deformation. The hybrid model based on the unified governing equations over Ω can be built as follows: Kinematic admissibility and compatibility: Static admissibility: Constitutive equations: where ε is the infinitesimal strain tensor, σ is the Cauchy stress tensor. E(x) is the stiffness tensor.
Using c 0 (ξ) as a reference, let us redefine C[x](|ξ|) by introducing a priori function, α is also called morphing function. Comparing Eq. (5) with Eq. (15), we can obtain that According to the constraint of energy density conservation, with assumptions of a small perturbation and a homogeneous deformation over the neighborhood, the stiffness tensor, E(x), at any point x can be expressed with the morphing function as follows [Lubineau, Azdoud, Han et al. (2012)], When one point x with its neighborhood is located in the peridynamic subdomain, one can know that E(x) = 0 and α(p) = 1 ∀p ∈ H δ (x). Thus the relation between E 0 and c 0 (ξ) can be obtained according to Eq. (16).
Remark. Now the accurate definitions of Ω 1 , Ω 2 and Ω m can be given by α as follows, 3 The adaptive algorithm for simulation of crack propagations In this section, we will describe the process of implementing the hybrid model for the adaptive simulation of crack propagations.

Initially introducing peridynamic subdomain
To study the crack initiations and propagations, it is necessary to introduce a peridynamic subdomain, i.e., Ω 2 , before implementing the adaptive hybrid model. The peridynamic subdomain, which is used for fracture, should cover the zone where cracks will firstly initiate. One feasible approach for assigning the initial peridynamic subdomain is that one can take the zone undergoing strong deformation gradients [Lubineau, Azdoud, Han et al. (2012)] during the elastic stage as the potential zone for crack initiations. For complicated structures, the zone under strong gradients can be obtained by calculating the deformation gradient field in advance.
According to Eq. (18), we introduce the peridynamic subdomain by assigning a zone with α = 1. This paper takes advantage of the circular zone with a radius r 1 , centered on x 0 which is designated as a flag point. Over the circular zone, we assign the morphing function, Meanwhile, adjacent to the circular zone, we assign an annular zone with an inner radius r 1 and an outer radius r 2 , centered on The morphing function of the remaining part of the structure is set as α x 0 (x) = 0, ∀x ∈ {x| |x − x 0 | ≥ r 2 }. Note that r 1 and r 2 are scalars, and r 2 > r 1 .
After assigning the distribution of the morphing function, according to Eq. (18), we can obtain the initial distribution of peridynamic subdomain, Ω 2 . Meanwhile, Ω 1 and Ω m are also obtained.

Discretization of the hybrid model
After introducing the peridynamic subdomain, we now focus on numerically carrying out the hybrid model by the finite element method (FEM). The principle of minimum potential energy is used to solve the problem, and the potential energy of the overall domain combining classical continuum mechanics and peridynamics can be written as follows [Azdoud, Han and Lubineau (2014); ], where After we discretized and minimized the potential energy, a linear algebraic equation set can be derived for finite element computations as follows, where {d} is a vector of all the nodal displacements, where [·] and {·} denote a matrix and a vector. n is the number of the total elements, h(x) is the number of the relative elements of point x, [N ] is the matrix of shape function, [H] denotes a matrix of differential operators and [R] is a diagonal matrix in which the diagonal entries may be 0 or 1, depending on the nodes of one element. Note that all these definitions are the same to that of FEM.
At the discretization level, the subdomains Ω 1 and Ω m are meshed by the finite element (FE) while Ω 2 is meshed by the discontinuous Galerkin finite element (DGFE) which was successfully applied to peridynamics by Azdoud et al. [Azdoud, Han and Lubineau (2014); Chen and Gunzburger (2011)]. In DGFE, each element does not share nodes with other elements and the DGFEs are connected by bonds that are defined by the relative position between quadrature points of elements. Consequently, the crack initiations and propagations are driven by bond breaking through the boundaries of elements.

Adaptive expansion of the peridynamic subdomain
When bonds begin to break, the peridynamic subdomain adaptively expands by introducing new peridynamic zones. The new peridynamic zones should cover these broken bonds and ensure that crack tips always stay far away from the boundary of the peridynamic subdomain. This paper introduces new peridynamic zones at the centroids of damaged elements in which at least a bond breaks, by assigning zones with α = 1. Note that the centroids of damaged elements are designated as flag points. Fig. 3 is a diagrammatic sketch that helps illustrate the process of the adaptive expansion.
In the beginning, an initial peridynamic subdomain was assigned at some point x 0 . When one bond in the peridynamic subdomain breaks, two new circular peridynamic zones are introduced at the centroids of damaged elements, x 1 and x 2 . The new peridynamic zone, centered on the flag point x 1 , is introduced by assigning α x 1 (x) = 1, ∀x ∈ {x| |x − x 1 | ≤ r i }, and the morphing function at the adjacent zone is set as 0 Meanwhile, we can define the other peridynamic zone which is centered on the flag point x 2 , by assigning α x 2 (x) = 1, ∀x ∈ {x| |x − x 2 | ≤ r i }, and the morphing function at the adjacent zone is defined as 0 With this adaptive approach, the peridynamic subdomain, Ω 2 , can be driven by broken bonds automatically without identifying crack tips. At the discretization level, the update of the peridynamic subdomain results in an update of mesh which leads to some elements transferring from FE to DGFE by inserting nodes for the new transferred elements, and all DGFEs are available for crack initiations and propagations (see Fig. 4). Note that the freedom of solution will increase after inserting new nodes for DGFEs. For one given increment, the structure is meshed by FE and DGFE to assign the peridynamic subdomain firstly. Next, the linear equations are solved, and the displacement results are used to calculate the stretch of bonds which is used for updating broken bonds. If a bond is decided to fail, the stiffness matrix that this broken bond contributed previously needs to update. Then the linear equations will be solved iteratively until there is no new bond failure.

Flowchart of the adaptive algorithm
After that, the adaptive algorithm will choose and store the new flag points that are the centroids of new damage elements in which a bond was broken at the given increment, and there was no bond breakage in previous increments. Then the same increment is carried out again.
If no new damage element appears, the algorithm will go to the next increment until the last loading increment is completed. Note that in this algorithm, the update of the peridynamic subdomain is performed between increments rather than iterations in each increment.   (16)) of the morphing subdomain is built based on an assumption of a homogeneous deformation over the neighborhood. Thus, it is interesting and necessary to validate the hybrid model for non-homogeneous deformations, even if the stiffness constraints in the morphing subdomain are used far from strong gradients [Lubineau, Azdoud, Han et al. (2012)]. In this section, a 2D numerical example for which we assume plane stress, under non-homogeneous deformations, is carried out by a platform that is based on the hybrid model, and ANSYS that is based on classical continuum mechanics, respectively.

A brief introduction of PeriFem for implementing the hybrid model
To make the hybrid model with an adaptive algorithm available for implementation, a software with a friendly user interface, named PeriFem, including solving, pre-and post-processing modules is developed. In the pre-processing module, we can make use of the software to generate a geometric model, mesh the model, input material parameters, and apply boundary conditions. In addition, the software can also take advantage of mesh data from other software, such as ANSYS and GMSH. In the solving module, the algorithm is developed by the C++ language, and it can implement 2D and 3D simulations. In the post-processing module, the software can display the synchronous calculation results for every increment. It is convenient for the users to use the mouse and keyboard to check the contour results. Also, it is very simple to display the contour results of the inner parts of a 3D model by a slicing technique. Additionally, the users can also output all the calculation data and store them in files.

Case description
A homogeneous sample is studied to validate the hybrid model. The geometric parameters and boundary conditions of the structure are shown in Fig. 6(a). The left boundary of the plate is fixed in x direction, and the middle point of the left boundary is also fixed in y direction. The displacement condition is applied to the right boundary in x direction with u x = 2 mm. Both the top and bottom of the plate are free. This numerical example implements bilinear quadrilateral elements in a finite element framework. The grid size ∆x of elements is 1.5 mm.
The stiffness parameters for classical continuum mechanics including Poisson's ratio and Young's modulus are ν = 1/3 and E = 30 GPa, respectively. The horizon, δ, of the neighborhood in peridynamics is 4.5 mm. The coefficient of micromodulus, c 0 (ξ), is assumed to be an exponential function, where τ 0 is a constant, which is calculated according to the Poisson's ratio and Young's modulus. l is a characteristic length which is determined to guarantee e −δ/l approaching 0.
Here, l is set to be 0.5 mm in this example.
To implement the hybrid model in the simulation, an initial seed of peridynamic subdomain is necessary. The initial peridynamic subdomain is seeded in the center of the hole with r 1 = 25 mm and r 2 = 45 mm see Fig. 6(b). In the orange zone, α = 1 while in the blue zone, α = 0. α decreases from 1 to 0 in the rainbow zone.
Note that the span of the rainbow zone influences the error (i.e., ghost forces for the energy-based coupling). Actually, the error is driven by δ/L c [Lubineau, Azdoud, Han et al. (2012)], where L c is the span of the rainbow zone (i.e., L c = r 2 − r 1 in this paper). If δ is given, the error reduces as L c becomes larger. . Note that the legend sort of (a) and (b) is different from that of (c) and (d) The displacement fields calculated by PeriFem and ANSYS are shown in Fig. 7. The displacement fields in x direction are in good agreement by comparing Fig. 7(a) and Fig. 7(c), while Fig. 7(b) and Fig. 7(d) display a high similarity of displacements in y direction.
To exactly illustrate the accuracy of the hybrid model, the displacements by PeriFem and ANSYS on three cutting lines are depicted and compared in Fig. 8 and Fig. 9 which also exhibit a good agreement. After calculating the relative error of displacements in the three cutting lines, we can find that the maximum relative error of the displacement in x direction is 1.43% for cutting line 1, and the maximum relative error of the displacement in y direction is 3.96% for cutting line 2. The maximum relative error of cutting line 3 in x and y directions are 0.5% and 1.5%, respectively. Note that the relative error is defined by | u p −u a u a |, where u p is the displacement calculated by the PeriFem while u a is the displacement calculated by ANSYS.  The results from two different approaches are in good agreement, and the comparison illustrates the accuracy of the hybrid model for non-homogeneous deformations.

Application of the hybrid model to fracture simulations
This section makes use of the hybrid model to predict crack initiations and propagations.
The following 2D numerical examples aim at illustrating the performance of the hybrid model on fracture simulations.
All the samples are assumed to be linear and elastic, and the 2D samples under plane stress are meshed by bilinear quadrilateral elements. Poisson's ratios for 2D analysis is ν = 1/3. The coefficient of micromodulus, c 0 (ξ), is defined in Eq. (24). The horizon, δ, of the neighborhood in peridynamics is chosen as 3∆x.

Brazilian disk
The Brazilian disk containing a single crack was experimentally tested by Haeri et al. [Haeri, Shahriar, Marji et al. (2014)] to study the breakage process of the disk. In this example, we apply the hybrid model to studying the properties of Brazilian disk under compression. The boundary conditions and dimensions of the sample are illustrated in Fig. 10(a). The grid size of elements is ∆x = 0.5 mm. The Young's modulus E is 3.1 GPa and l = 0.1 mm. In the simulation, two seeds of the initial peridynamic subdomain are centered around the crack tips with r 1 = 4 mm and r 2 = 8 mm see Fig. 10(b). The critical stretch, s crit , is set to be 0.055. This simulation is performed in 100 progressive increments (steps) so that the increment of displacement yields ∆u y = 0.012 mm in each step.
From simulation results (see Fig. 11), it is validated that the hybrid model can simulate the whole process of the elastic deformation, crack initiations and propagations, and the structure failure. Fig. 11 including six loading steps, exhibits the failure process and the contour of damage which is defined by Silling et al. [Silling and Askari (2005)]. There is no crack at Step 83 ( Fig. 11(a)) while crack suddenly appears at step 84 ( Fig. 11(b)). Then the crack grows slowly until Step 89 at which the crack propagates further. A disastrous fracture suddenly appears at Step 92. From the figures, we can find that the fractures initiate from the pre-crack tips and propagate toward the direction of loading points. The failure path of the disk by the hybrid model is very similar to the experimental result (Fig. 6(c)) in Haeri et al. [Haeri, Shahriar, Marji et al. (2014)].
Note that the grids shown in Figs. 11, 14, 15, 18 and 21 are DGFEs where peridynamics is performed. As the crack propagates, the DGFEs adaptively update. From the simulation, we can find that peridynamics is always restricted to a relatively small zone, compared to 1. Introduction the whole structure, which can highly reduce the computational cost. Fig. 12(a) displays the load-displacement curve of the resultant force on the upper loading boundary. Fig. 12(b) is an enlargement of the curve in the red box of Fig. 12(a). The load-displacement curve is approximately linear before Step 83. Zooming in on the peak region shows that the resultant of the forces encounters some drops. With the help of Fig. 11, we can pinpoint that the force experiences a small drop when the crack appears at Step 84, and another small drop is encountered at Step 89 when the crack propagates further. The drastic drop of the force at Step 92 keeps in step with the disastrous fracture. From Fig. 11(f), we find that there are still connections in the disk, which accounts for the residual bearing capacity after Step 92.

A notched plate with an off-center circular hole
We apply the hybrid model to simulating a notched plate with an off-center circular hole under unilateral tension. This test has been simulated by various investigators [Dipasquale, Zaccariotto and Galvanetto (2014); Tabiei and Wu (2003); Ni, Zhu, Zhao et al. (2018)] for predicting the crack initiations and propagations. The boundary conditions and dimensions of the sample are illustrated in Fig. 13(a). The grid size of elements is ∆x = 0.25 mm. The Young's modulus is 71.4 GPa and l = 0.05 mm. In the simulation, one seed of the initial peridynamic subdomain, centered around the notch tip with r 1 = 4 mm and r 2 = 8 mm, Step 89 Step 84 Step 83 Step 88 Step 92 Step 91 Force ( [Tabiei and Wu (2003)] and Ni et al. [Ni, Zhu, Zhao et al. (2018)]. By comparing Fig. 14 and Fig. 15, we can conclude that different distances between the notch and the hole can induce different paths of crack propagations. Additionally, the simulation results also display the adaptive expansion of the peridynamic subdomain. Then the force drops, which corresponds to the crack appearance (see Fig. 14(a)). After that, the force encounters a drastic drop, which corresponds to the crack splitting the structure (see Fig. 14(c)). For the specimen with h = 15 mm, the force experiences a drop, which corresponds to the crack appearance (see Fig. 15(a)). Then the force grows up. The force encounters another drop, which corresponds to the crack reaching the hole (see Fig. 15(c)). Since the structure is not destroyed in this simulation, it still has a relatively large residual bearing capacity.

A notched beam under four-point bending
This numerical example is a four-point bending test of a single-edge notched beam. The boundary conditions and dimensions of the sample are illustrated in Fig. 17(a). The grid size of elements is ∆x = 0.5 mm. The Young's modulus is 30 GPa and l = 0.1 mm. In the simulation, one seed of the initial peridynamic subdomain, centered around the notch tip with r 1 = 4 mm and r 2 = 6 mm, is introduced (see Fig. 17(b)). The critical stretch, s crit , is set to be 0.05. This simulation is performed in 100 progressive increments (steps) so that the increment of displacement yields ∆u y = 0.02 mm in each step. Step 81, and after several steps, the crack grows rapidly to about 21 mm. The crack under four-point bending is almost straight. The small curve in the fracturing path is induced by the unstructured mesh. In Fig. 18, the adaptive expansion of the peridynamic subdomain always covers the crack tip and is consistent with the crack propagation. Fig. 19 displays the curves of the forces on the loading points of the beam versus the displacement, which reflect the brittle fracture in the test. The curves approximate linear before a drastic drop. The drastic drop corresponds to the sudden appearance of the long crack (see Fig. 18(b)). The stairs of the load-displacement curves are induced by crack propagations. The two load-displacement curves are consistent until the displacement reaches 1.88 mm. The branching of the two curves appears due to the asymmetry of the structure which is induced by the crack propagation (see Fig. 18(d)).   (2018)]. Its boundary conditions and dimensions are illustrated in Fig. 20(a). The grid size of elements is ∆x = 1 mm. The Young's modulus is 30 GPa and l = 0.2 mm. In the simulation, two seeds of the initial peridynamic subdomain are centered around the notch tips with r 1 = 8 mm and r 2 = 16 mm (see Fig. 20(b)). The critical stretch, s crit , is set to be 0.1. This simulation is performed in 50 progressive increments (steps). In each step, ∆u y = 0.08 mm and ∆u x = 0.02 mm. Fig. 21 is the simulation results that display the evolution of the damage and fracture. The damage and fracture initiate from each notch tip and propagate towards the opposite notch.
In the beginning, the crack which is initiated from the left notch tends to the bottom of the plate while another crack tends to the upper. This kind of fracture path is induced by the asymmetric boundary conditions. Later the crack tips form a neck-zone where shear deformations deflected the evolution of the cracks. The proceeding evolution finally brings the upper crack tip into contact with the lower crack, which eventually leads to the failure of the plate. In addition, the adaptive expansion of the peridynamic subdomain is always restricted to a relatively small zone as the cracks propagate, which can make the simulation much more efficient. Step 20 lead to a small residual bearing capacity of the structure. We can conclude that the sudden appearance of the crack at Step 20 (see Fig. 21) is a deadly fracture for this structure.

Conclusion
In this paper, the hybrid local/nonlocal continuum model is introduced, and the adaptive coupling approach is described in detail. Firstly, for the purpose of validation, a plate with a circular hole is simulated by the hybrid model and the classical continuum mechanics. The simulation shows that the results from both models are in good agreement, and the relative errors of the cutting lines are less than 3.96%. This validation illustrates the accuracy of the hybrid model under non-homogeneous deformations. Secondly, this paper applies this hybrid model to fracture simulations. Various numerical examples such as a Brazilian disk, a notched plate with a hole, a notched beam and a double-edge-notched plate are simulated to predict complicated crack initiations and propagations. The simulation results exhibit that (i) the hybrid model can simulate the whole process of the elastic deformation, crack initiations and propagations, and the structure failure, (ii) the process of the adaptive expansion of the peridynamic subdomain which is always restricted to a relatively small zone displays the efficiency of the hybrid model. Thirdly, the peridynamic software, PeriFem, also displays its performance through the numerical example.