A Nonlinear Crack Model for Concrete Structure Based on an Extended Scaled Boundary Finite Element Method

Fracture mechanics is one of the most important approaches to structural safety analysis. Modeling the fracture process zone (FPZ) is critical to understand the nonlinear cracking behavior of heterogeneous quasi-brittle materials such as concrete. In this work, a nonlinear extended scaled boundary finite element method (X-SBFEM) was developed incorporating the cohesive fracture behavior of concrete. This newly developed model consists of an iterative procedure to accurately model the traction distribution within the FPZ accounting for the cohesive interactions between crack surfaces. Numerical validations were conducted on both of the concrete beam and dam structures with various loading conditions. The results show that the proposed nonlinear X-SBFEM is capable of modeling the nonlinear fracture propagation process considering the effect of cohesive interactions, thereby yielding higher precisions than the linear X-SBFEM approach.


Introduction
With the development of numerical analysis technology, structural fracture mechanics is an important approach to structural safety evaluation. Fracture process zone (FPZ) is defined as the intermediate space between cracked and uncracked portions of concrete [1]. Different from real cracks, the FPZ can still transmit stress, and the stress, σ, that FPZ transmits decreases with increasing crack open displacement, w. When the crack open displacement reaches a certain critical value, w c , the surface force of the crack surface becomes zero, as illustrated in Figure 1. The FPZ consists of microcracks, which are minute individual cracks; this gives rises to the cohesive tractions ahead of the crack tip, which comes from the aggregate interlocking and surface friction. Therefore, a nonlinear fracture-mechanics-based method needs to be applied to account for the effect of cohesive tractions during the fracture propagations. Two main approaches are often taken to model the FPZ, which are: the smeared crack models and the discrete crack models. The smeared crack models proposed by Rashid [2] are based on the continuum approach, where the computational mesh of the FPZ remains constant while the fracture propagation is modeled by the growth of a number of parallel cracks smeared over the elements within the FPZ. Using this approach, Bhattacharjee et al. [3] and Calayir et al. [4] successfully carried out dynamic cracking analyses of the Koyna gravity dam under the influence of nonorthogonal cracks. Cai et al. [5] also follow a similar approach incorporating the linear or bilinear softening dispersion crack models to predict the crack response of concrete gravity dams. The damage-based fracture mechanics was also developed sharing the similar concept by Bazant [6]. These methods have shown tremendous success in modeling concrete fractures under complex loading and boundary conditions. However, the strong mesh dependency often causes complications in determining the characteristic length scale, fracture strength, and fracture toughness to accurately describe the fracturing propagation process. The smeared crack models mainly combine with damage mechanics for crack propagation studies.
On the other hand, the discrete crack model proposed by Dugdale [7] and Barenblatt [8] utilizes a predefined fracture path as a part of the computational domain boundary, which reduces the mesh dependency comparing to the smeared crack model. Hillerborg et al. [9] firstly implemented the cohesive zone model (CZM) to describe the FPZ based on the discrete crack approach. Following this approach, numerous scholars have improved and implemented the cohesive zone method to model the fracture propagation process (Skrikerud et al. [10], Ayari et al. [11], Xie et al. and Yang et al. [12][13][14]). Among these improvements, the scaled boundary element method (SBFEM) was developed to model the complex fracture growth in terms of fracture branching and coalescing under complex loading conditions. However, the SBFEM fails to capture the nonlinearity brought by the cohesive interactions within the FPZ. Therefore, the mesoscale and atomic-scale-inserted CZMs were developed to model the FPZ [15][16][17][18]. However, the predetermined fracture paths often cause computational complexities and inaccuracies.
Recently, the extended scaled boundary finite element method (X-SBFEM) based on the level set method was developed on the basis of both the SBFEM [19,20] and the extended finite elements (XFEM) [21,22]. Capitalizing on the advantages of both methods, X-SBFEM [23,24] can make full use of XFEM to describe the discontinuous displacement field and SBFEM to solve the problem of the stress singularity with higher precision. Simulating the crack body section using XFEM and the crack tip using SBFEM, the method finally establishes the total equilibrium equation of the crack body and solves the equation, thereby overcoming the disadvantages of XFEM, such as obtaining the analytical form of the displacement and the stress asymptotic fields of the crack tip in advance and constructing the complex enhanced functions, which can express nonsmooth behaviors near the crack tip. In some special circumstances, the enrichment functions are discontinuous or have nonpolynomial forms to specially address the issue when the stiffness matrix is constructed using numerical integration.
Currently, the application of X-SBFEM in fractures [23] based on linear elastic fracture mechanics (LEFM) mainly focuses on the linear fracture. Due to the complexity of the three-dimensional analysis Two main approaches are often taken to model the FPZ, which are: the smeared crack models and the discrete crack models. The smeared crack models proposed by Rashid [2] are based on the continuum approach, where the computational mesh of the FPZ remains constant while the fracture propagation is modeled by the growth of a number of parallel cracks smeared over the elements within the FPZ. Using this approach, Bhattacharjee et al. [3] and Calayir et al. [4] successfully carried out dynamic cracking analyses of the Koyna gravity dam under the influence of nonorthogonal cracks. Cai et al. [5] also follow a similar approach incorporating the linear or bilinear softening dispersion crack models to predict the crack response of concrete gravity dams. The damage-based fracture mechanics was also developed sharing the similar concept by Bazant [6]. These methods have shown tremendous success in modeling concrete fractures under complex loading and boundary conditions. However, the strong mesh dependency often causes complications in determining the characteristic length scale, fracture strength, and fracture toughness to accurately describe the fracturing propagation process. The smeared crack models mainly combine with damage mechanics for crack propagation studies.
On the other hand, the discrete crack model proposed by Dugdale [7] and Barenblatt [8] utilizes a predefined fracture path as a part of the computational domain boundary, which reduces the mesh dependency comparing to the smeared crack model. Hillerborg et al. [9] firstly implemented the cohesive zone model (CZM) to describe the FPZ based on the discrete crack approach. Following this approach, numerous scholars have improved and implemented the cohesive zone method to model the fracture propagation process (Skrikerud et al. [10], Ayari et al. [11], Xie et al. and Yang et al. [12][13][14]). Among these improvements, the scaled boundary element method (SBFEM) was developed to model the complex fracture growth in terms of fracture branching and coalescing under complex loading conditions. However, the SBFEM fails to capture the nonlinearity brought by the cohesive interactions within the FPZ. Therefore, the mesoscale and atomic-scale-inserted CZMs were developed to model the FPZ [15][16][17][18]. However, the predetermined fracture paths often cause computational complexities and inaccuracies.
Recently, the extended scaled boundary finite element method (X-SBFEM) based on the level set method was developed on the basis of both the SBFEM [19,20] and the extended finite elements (XFEM) [21,22]. Capitalizing on the advantages of both methods, X-SBFEM [23,24] can make full use of XFEM to describe the discontinuous displacement field and SBFEM to solve the problem of the stress singularity with higher precision. Simulating the crack body section using XFEM and the crack tip using SBFEM, the method finally establishes the total equilibrium equation of the crack body and solves the equation, thereby overcoming the disadvantages of XFEM, such as obtaining the analytical form of the displacement and the stress asymptotic fields of the crack tip in advance and constructing the complex enhanced functions, which can express nonsmooth behaviors near the crack tip. In some special circumstances, the enrichment functions are discontinuous or have nonpolynomial forms to specially address the issue when the stiffness matrix is constructed using numerical integration.
Currently, the application of X-SBFEM in fractures [23] based on linear elastic fracture mechanics (LEFM) mainly focuses on the linear fracture. Due to the complexity of the three-dimensional analysis model, the theoretical researches of fracture mechanics still focus on the state of the two-dimensional analysis model. Therefore, based on the X-SBFEM algorithm, this paper introduces a nonlinear crack model which adopts the linear superposition of iterative methods to incorporate the cohesive interactions within the FPZ. The proposed approach was then implemented to model mixed-mode fracture of concrete beam and gravity dam structures. The results show improvements comparing to other methods. Close agreements were found between the numerical and experimental results.
The contents of this paper are arranged as follows. In Section 2, we explain the principle of X-SBFEM. In Section 3, a nonlinear crack model with iterative method for cohesive interactions in the FPZ is introduced and a flowchart for solving the cohesion is given. In Section 4, four numerical simulations (a three-point bending beam, a four-point shear beam, an experimental concrete gravity dam with single crack expansion, and a static cohesive crack propagation simulation of Koyna Dam) are modeled to validate the nonlinear model. In Section 5, we conclude that this paper has developed the X-SBFEM with the nonlinear model to improve the modeling of crack propagation.

Extended Scaled Boundary Finite Element Method
The core content of X-SBFEM based on the level set method focuses on the simulation of the nonsmooth behavior near the crack tip using semianalytical SBFEM in the form of super-elements and the simulation of the crack body using XFEM. The key lies in the way the algorithm addresses the boundary conditions at the joint. Figure 2 below shows the topological relationship in the model domain including a crack based on X-SBFEM.
Appl. Sci. 2017, 7, x FOR PEER REVIEW 3 of 20 model, the theoretical researches of fracture mechanics still focus on the state of the two-dimensional analysis model. Therefore, based on the X-SBFEM algorithm, this paper introduces a nonlinear crack model which adopts the linear superposition of iterative methods to incorporate the cohesive interactions within the FPZ. The proposed approach was then implemented to model mixed-mode fracture of concrete beam and gravity dam structures. The results show improvements comparing to other methods. Close agreements were found between the numerical and experimental results. The contents of this paper are arranged as follows. In Section 2, we explain the principle of X-SBFEM. In Section 3, a nonlinear crack model with iterative method for cohesive interactions in the FPZ is introduced and a flowchart for solving the cohesion is given. In Section 4, four numerical simulations (a three-point bending beam, a four-point shear beam, an experimental concrete gravity dam with single crack expansion, and a static cohesive crack propagation simulation of Koyna Dam) are modeled to validate the nonlinear model. In Section 5, we conclude that this paper has developed the X-SBFEM with the nonlinear model to improve the modeling of crack propagation.

Extended Scaled Boundary Finite Element Method
The core content of X-SBFEM based on the level set method focuses on the simulation of the nonsmooth behavior near the crack tip using semianalytical SBFEM in the form of super-elements and the simulation of the crack body using XFEM. The key lies in the way the algorithm addresses the boundary conditions at the joint. Figure 2 below shows the topological relationship in the model domain including a crack based on X-SBFEM.

Extended Finite Element Method
Based on the partition of unity methods simulating the crack body by XFEM, the general formula of a displacement field [25,26] is where N fem and N e respectively represent a node in a general element and an enriched node in an internal split crack. Correspondingly, is a normal degree-of-freedom, is a generalized degreeof-freedom related to (as shown in the square node in Figure 2), and ϑ(x) is the Heaviside step function.

Extended Finite Element Method
Based on the partition of unity methods simulating the crack body by XFEM, the general formula of a displacement field [25,26] is where N fem and N e respectively represent a node in a general element and an enriched node in an internal split crack. Correspondingly, q I is a normal degree-of-freedom, a J is a generalized degree-of-freedom related to ϑ (as shown in the square node in Figure 2), and ϑ(x) is the Heaviside step function. The equilibrium equation is where K e aa and K e bb represent the stress matrix respectively related to a normal degree-of-freedom and a generalized degree-of-freedom. Moreover, K e ab and K e ba are the coupling matrices, and P a and P b are the equivalent nodal forces in accordance to general degrees-of-freedom and generalized degrees-of-freedom, respectively.
where Γd and Γt respectively represent the crack face and the force interface. Moreover, p v and t respectively represent the body force in computational domain and the surface force on the force interface. According to the advantage of the description of discontinuous displacement field description, XFEM is used to simulate the main body of the crack.

Scaled Boundary Finite Element Method
X-SBFEM is used to simulate the crack tip for the high efficiency and high precision of stress singular field simulations. As shown in Figure 3, there are side-face forces at the face of the super-element at the crack tip in the case of SBFEM. Without taking the body force into account, the displacement field and the stress field given by the SBFEM are [21] {u(ξ, where N(η) is the interpolation shape function for one-dimensional line elements and {ϕ i } and λ i are the displacement modes and the eigenvalues, respectively. Furthermore, [B 1 (η)] and [B 2 (η)] are determined by the geometric shape of the boundary of the element. The formulas calculated by the virtual work principle are where {P} is the equivalent boundary nodal force and E 0 , E 1 , and E 2 are the coefficient matrices of the SBFEM governing equations. Equation (11) is a second-order homogeneous ordinary differential equation and its solution is where is the weight of this mode. Substitute Equation (12) into Equation (10) and get a quadratic eigenvalue problem.
Substitute Equation (12) into Equation (10) and get the equivalent node force on the boundary corresponding to the displacement mode.
Combining Equations (13) and (14) and introducing auxiliary variables can transform quadratic eigenvalue problems into standard linear eigenvalue problems. where and [ ] are modal matrices of the displacement modal matrix and the force modal matrix. For the square root singular problem for homogeneous materials, the stress intensity factors can be defined as where L0 is the distance between the scaling center and the point of the crack surface (the segment OA in Figure 3). From Equation (8), it can be deduced that if ≥ 1, the distribution of the stress mode tends to 0 when ξ → 0 [27]. Figure 4 shows a typical method simulating the crack with the coupling of the XFEM and the SBFEM [28] domain. In the process of XFEM and X-SBFEM coupling, there is a problem of the balance of virtual and real freedom, continuous displacement, and force. To solve the problem, the SBFEM simulates the crack surface using the boundary of the element, while the XFEM introduces an additional degree-of-freedom based on the use of the step function to describe the discontinuous displacement field. To coordinate the displacement of the two types of elements along the boundary, a transition matrix T is used to match the nodal displacements in SBFEM with those in XFEM. A transformation matrix, which can be derived from the previous equations that translates the Equation (11) is a second-order homogeneous ordinary differential equation and its solution is

Force Balance and Displacement Coordination
where c i is the weight of this mode. Substitute Equation (12) into Equation (10) and get a quadratic eigenvalue problem. [ Substitute Equation (12) into Equation (10) and get the equivalent node force q i on the boundary corresponding to the displacement mode.
Combining Equations (13) and (14) and introducing auxiliary variables can transform quadratic eigenvalue problems into standard linear eigenvalue problems.
where [Φ] and [Q] are modal matrices of the displacement modal matrix and the force modal matrix. For the square root singular problem for homogeneous materials, the stress intensity factors can be defined as where L 0 is the distance between the scaling center and the point of the crack surface (the segment OA in Figure 3). From Equation (8), it can be deduced that if λ i ≥ 1, the distribution of the stress mode tends to 0 when ξ → 0 [27].  Figure 4 shows a typical method simulating the crack with the coupling of the XFEM and the SBFEM [28] domain. In the process of XFEM and X-SBFEM coupling, there is a problem of the balance of virtual and real freedom, continuous displacement, and force. To solve the problem, the SBFEM simulates the crack surface using the boundary of the element, while the XFEM introduces an additional degree-of-freedom based on the use of the step function to describe the discontinuous displacement field. To coordinate the displacement of the two types of elements along the boundary, a transition matrix T is used to match the nodal displacements in SBFEM with those in XFEM. A transformation matrix, which can be derived from the previous equations that translates the unknown SBFEM nodal displacements (u E , u F , u A , and u B ) to the unknown XFEM nodal displacements (q 2 , q 3 , a 2 , and a 3 ) can be described as

Force Balance and Displacement Coordination
where I is a unit matrix. In order to ensure the compatibility of the displacement and to integrate the element stiffness matrix, it is necessary to rearrange the column displacement vectors of the SBFEM and the stiffness matrix according to whether the nodes are on the common boundary, where the subscript a represents the nodes on the common boundary of SBFEM and XFEM and subscript b represents the nodal degree-of-freedom of the noncommon boundary. The transformation matrix T is only related to the interpolation shape function at the crack opening at the common boundary of the SBFEM and the XFEM.
where I is a unit matrix. In order to ensure the compatibility of the displacement and to integrate the element stiffness matrix, it is necessary to rearrange the column displacement vectors of the SBFEM and the stiffness matrix according to whether the nodes are on the common boundary, = where the subscript a represents the nodes on the common boundary of SBFEM and XFEM and subscript b represents the nodal degree-of-freedom of the noncommon boundary. The transformation matrix T is only related to the interpolation shape function at the crack opening at the common boundary of the SBFEM and the XFEM.

A Nonlinear Crack Model with Iterative Method for Cohesive Interactions in the FPZ
The relative displacement of the crack surface, including the opening displacement (COD) and the sliding displacement (CSD) of the crack surface, does not exceed the limits shown in Figure 5. When the COD and CSD of the crack surface are not beyond the limits shown in Figure 5, the total loads generated by the structure include the external loads and the cohesive forces at the virtual crack surfaces. However, if the load exceeds the limit, the cohesive force is 0. In the case of cohesive forces, considering the stress intensity factor I as an example, the stress intensity factor consists of two parts [29], where is the total stress intensity factor and and are the components related to the external and cohesive forces, respectively. All three stress intensity factors can be calculated by the standard SBFEM solution stress intensity factor formula. Thus, > 0 when the crack opens as a

A Nonlinear Crack Model with Iterative Method for Cohesive Interactions in the FPZ
The relative displacement of the crack surface, including the opening displacement (COD) and the sliding displacement (CSD) of the crack surface, does not exceed the limits shown in Figure 5. When the COD and CSD of the crack surface are not beyond the limits shown in Figure 5, the total loads generated by the structure include the external loads and the cohesive forces at the virtual crack surfaces. However, if the load exceeds the limit, the cohesive force is 0. In the case of cohesive forces, considering the stress intensity factor I as an example, the stress intensity factor consists of two parts [29], where K I is the total stress intensity factor and K P I and K C I are the components related to the external and cohesive forces, respectively. All three stress intensity factors can be calculated by the standard SBFEM solution stress intensity factor formula. Thus, K P I > 0 when the crack opens as a result of the external force of the model, while K C I < 0 when the crack tends to close owing to the cohesive force. Equivalently, K I = 0 when force balance is achieved as a result of the roles of the external and cohesive forces. Therefore, K I ≥ 0 can be used as the criterion for judging whether the crack will continue to propagate or not [15]. . The area between the curves in Figure 5c is twice that of the mode II fracture energy, [30]. The key concept of this method is based on the relative displacement of the crack surface and its application to the linear superposition of an iterative scheme to solve and estimate the cohesive tractions on the crack surface. Specifically, A. Assume that the structure is only affected by the external force so that the relative displacement ∆ of the super-element crack surface ∆ can be obtained based on the linear elastic assumptions of X-SBFEM, and that the corresponding cohesive traction can be obtained according to Figure 5. B. As shown in Figure 6, the external force and the cohesive force obtained in the previous step are applied to the structure, wherein the cohesive traction is applied in the form of a side-face force, formulated in accordance to the following equation: where Assume that the structure is only affected by the external force so that the relative displacement ∆u 1 of the super-element crack surface ∆u 1 can be obtained based on the linear elastic assumptions of X-SBFEM, and that the corresponding cohesive traction t 1 can be obtained according to Figure 5. B. As shown in Figure 6, the external force and the cohesive force obtained in the previous step are applied to the structure, wherein the cohesive traction t 1 is applied in the form of a side-face force, formulated in accordance to the following equation: where Appl. Sci. 2017, 7, x FOR PEER REVIEW 8 of 20 Assume that the load can be expressed by the power series, The corresponding displacement mode is Substituting Equation (24) into Equations (23) and (11) Thus, the complete displacement of the boundary nodes and the equivalent nodal force are respectively, The SBFEM nonhomogeneous control equation can be easily obtained in accordance to Assume that the load can be expressed by the power series, The corresponding displacement mode is Substituting Equation (24) into Equations (23) and (11) Thus, the complete displacement of the boundary nodes and the equivalent nodal force are respectively, [Φ] and [Q] are modal matrices of the displacement modal matrix and the force modal matrix solved by Equation (15), respectively. The following equation can be obtained using Equations (27) and (28): The equivalent nodal force of the super-element boundary node generated by the distributed load of the crack surface is The displacement field of the crack tip element is where {ψ i } is the stress mode solved based on SBFEM. The relative displacement ∆u i+1 is solved using Equation (32).
C. Repeat the steps until the relationship between t i and ∆u i+1 becomes consistent with the pattern of variation plotted in Figure 5.
A simplified flowchart of the solution of the cohesive tractions is shown in Figure 7.
where { } is the stress mode solved based on SBFEM. The relative displacement ∆ is solved using Equation (32).
C. Repeat the steps until the relationship between and ∆ becomes consistent with the pattern of variation plotted in Figure 5.
A simplified flowchart of the solution of the cohesive tractions is shown in Figure 7.

A Three-Point Bending Beam
This model was first studied by Petersson (1981) as an experimental study of the mode I fracture propagation problem [31]. The geometry, boundary conditions, and material parameters of the beam are shown in Figure 8. The tensile strength f t is 3.33 MPa and the mode I fracture energy G f I is 137 N/m. The present example predicts the crack propagation path based on the LEFM maximum circumferential tensile stress criterion. Analyzing the single linear softening curve ( Figure 5b) and based on the mode I fracture energy G f I , the limit value of the linear softening curve w c is 0.0823 mm. The results of three crack propagation steps a = 10 mm, 20 mm, and 30 mm, for a 20 × 200 grid density, are calculated and compared with the results based on the linear elasticity method [24].

A Three-Point Bending Beam
This model was first studied by Petersson (1981) as an experimental study of the mode I fracture propagation problem [31]. The geometry, boundary conditions, and material parameters of the beam are shown in Figure 8. The tensile strength is 3.33 MPa and the mode I fracture energy is 137 N/m. The present example predicts the crack propagation path based on the LEFM maximum circumferential tensile stress criterion. Analyzing the single linear softening curve ( Figure 5b) and based on the mode I fracture energy , the limit value of the linear softening curve is 0.0823 mm. The results of three crack propagation steps a = 10 mm, 20 mm, and 30 mm, for a 20 × 200 grid density, are calculated and compared with the results based on the linear elasticity method [24].   [12] and experimental results published by Petersson [27], as shown in Figure 9. We can see that the elicited results based on LEFM are very different from the experimental data. This is particularly evident for the peak load, which is much higher than the experimental peak. This is because the LEFM-based method is incapable of simulating the energy dissipation of the fracture zone. Based on X-SBFEM, this study has used different methodologies to solve the cohesive tractions based on the iterative method of linear superposition by simulating the energy dissipation of FPZ. Figure 9b shows the Load-LPD curves for three different crack propagation steps considering FPZ nonlinearities. It can be seen from Figure 9 that the calculated results are in good agreement with the experimental results of Petersson, which shows that the method used in this study can simulate the energy dissipation of FPZ. Moreover, it can be seen from Figure 9 that the results of the three crack propagation steps are in good agreement with the experimental curve, which shows that different steps have minor effects on the calculated results.   [12] and experimental results published by Petersson [27], as shown in Figure 9. We can see that the elicited results based on LEFM are very different from the experimental data. This is particularly evident for the peak load, which is much higher than the experimental peak. This is because the LEFM-based method is incapable of simulating the energy dissipation of the fracture zone. Based on X-SBFEM, this study has used different methodologies to solve the cohesive tractions based on the iterative method of linear superposition by simulating the energy dissipation of FPZ. Figure 9b shows the Load-LPD curves for three different crack propagation steps considering FPZ nonlinearities. It can be seen from Figure 9 that the calculated results are in good agreement with the experimental results of Petersson, which shows that the method used in this study can simulate the energy dissipation of FPZ. Moreover, it can be seen from Figure 9 that the results of the three crack propagation steps are in good agreement with the experimental curve, which shows that different steps have minor effects on the calculated results.

A Four-Point Shear Beam
Arrea and Ingraffea first tested and analyzed the four-point unilateral shear beam [28]. The geometry and boundary conditions of the beam are shown in Figure 10. Assuming that the structure is in the plane stress state, the Young's modulus E is 24.8 GPa, Poisson's ratio is 0.18, tensile strength is 3.0 MPa, mode I fracture energy is 100 N/m, and mode II fracture energy is 10 N/m. The crack path is also predicted by using the LEFM-based maximum circumferential stress criterion. Analyzing the linear softening curve of Figure 5b, the limit value of CODs ᴡ obtained by the mode I fracture energy is 0.067 mm and the limit value of CSDs obtained by the mode II fracture energy is 0.02 mm. The results of the three crack propagation steps a = 10 mm, 20 mm, and 30 mm, using a 20 × 200 grid density, are calculated and compared with the results based on the linear elasticity method [24]. Figures 11a and 12a show the relationship between the load calculated by LEFM, the crack mouth sliding displacement (Load-CMSD), and the relationship between the load and the loading point displacement (Load-LPD). It can be seen from the figure that the calculated results are close to the numerical solutions of Yang et al. [12] and the effect of different steps on the calculation results is not considerable, which proves the applicability of the X-SBFEM algorithm to complex crack propagation problems. However, the FPZ energy dissipation of the crack tip makes the LEFM-based

A Four-Point Shear Beam
Arrea and Ingraffea first tested and analyzed the four-point unilateral shear beam [28]. The geometry and boundary conditions of the beam are shown in Figure 10. Assuming that the structure is in the plane stress state, the Young's modulus E is 24.8 GPa, Poisson's ratio υ is 0.18, tensile strength f t is 3.0 MPa, mode I fracture energy G f I is 100 N/m, and mode II fracture energy G f I I is 10 N/m. The crack path is also predicted by using the LEFM-based maximum circumferential stress criterion. Analyzing the linear softening curve of Figure 5b, the limit value of CODs w c obtained by the mode I fracture energy G f I is 0.067 mm and the limit value of CSDs s c obtained by the mode II fracture energy G f I I is 0.02 mm. The results of the three crack propagation steps a = 10 mm, 20 mm, and 30 mm, using a 20 × 200 grid density, are calculated and compared with the results based on the linear elasticity method [24]. Figures 11a and 12a show the relationship between the load calculated by LEFM, the crack mouth sliding displacement (Load-CMSD), and the relationship between the load and the loading point displacement (Load-LPD). It can be seen from the figure that the calculated results are close to the numerical solutions of Yang et al. [12] and the effect of different steps on the calculation results is not considerable, which proves the applicability of the X-SBFEM algorithm to complex crack propagation problems. However, the FPZ energy dissipation of the crack tip makes the LEFM-based method slightly different compared to previous calculations and experimental data [15]. In this study, the iterative method used to simulate the cohesive tractions is used to simulate the energy dissipation of the FPZ. The linear softening curve in Figure 5b and the curve of Figure 5c are used to solve the cohesive tractions of the vertical crack surface and the parallel crack surface. As shown in Figure 11b, it can be seen that the Load-CMSD curve calculated herein is in good agreement with Yang's experimental data (NFM-based). Figure 12b shows that the method considered herein can describe the snap-back phenomenon of the Load-LPD curve. It can be concluded that the iterative, linear superposition method based on X-SBFEM can yield high-precision results without coupling the interface unit (CIEs) near the crack tip.
Appl. Sci. 2017, 7, x FOR PEER REVIEW 12 of 20 method slightly different compared to previous calculations and experimental data [15]. In this study, the iterative method used to simulate the cohesive tractions is used to simulate the energy dissipation of the FPZ. The linear softening curve in Figure 5b and the curve of Figure 5c are used to solve the cohesive tractions of the vertical crack surface and the parallel crack surface. As shown in Figure 11b, it can be seen that the Load-CMSD curve calculated herein is in good agreement with Yang's experimental data (NFM-based). Figure 12b shows that the method considered herein can describe the snap-back phenomenon of the Load-LPD curve. It can be concluded that the iterative, linear superposition method based on X-SBFEM can yield high-precision results without coupling the interface unit (CIEs) near the crack tip. method slightly different compared to previous calculations and experimental data [15]. In this study, the iterative method used to simulate the cohesive tractions is used to simulate the energy dissipation of the FPZ. The linear softening curve in Figure 5b and the curve of Figure 5c are used to solve the cohesive tractions of the vertical crack surface and the parallel crack surface. As shown in Figure 11b, it can be seen that the Load-CMSD curve calculated herein is in good agreement with Yang's experimental data (NFM-based). Figure 12b shows that the method considered herein can describe the snap-back phenomenon of the Load-LPD curve. It can be concluded that the iterative, linear superposition method based on X-SBFEM can yield high-precision results without coupling the interface unit (CIEs) near the crack tip.

An Experimental Concrete Gravity Dam with Single-Crack Expansion
Carpinteri [32], Barpi [33], and Shi [34] tested and analyzed the single-crack expansion gravity dam. A 1:40 scale concrete gravity dam model is used herein. The geometry, boundary conditions, and material parameters of the gravity dam are shown in Figure 13. Consisting of concrete, the gravity dam is analyzed based on a plane strain assumption and on the bilinear softening curve, where ω c = 0.256 mm, G f = 184 N/m, and f t = 3.6 MPa. Suppose that the indentation length is 1/10 W (0.15 m), where W is the width of the dam at the elevation of the indentation. Hydrostatic pressure exists in the upstream face of the dam. This hydrostatic pressure acts on the upstream face and can be equivalently replaced by four concentrated loads, as shown in Figure 13. The hydrostatic pressure gradually increases until the dam breaks.

An Experimental Concrete Gravity Dam with Single-Crack Expansion
Carpinteri [32], Barpi [33], and Shi [34] tested and analyzed the single-crack expansion gravity dam. A 1:40 scale concrete gravity dam model is used herein. The geometry, boundary conditions, and material parameters of the gravity dam are shown in Figure 13. Consisting of concrete, the gravity dam is analyzed based on a plane strain assumption and on the bilinear softening curve, where = 0.256 mm, Gf = 184 N/m, and ft = 3.6 MPa. Suppose that the indentation length is 1/10 W (0.15 m), where W is the width of the dam at the elevation of the indentation. Hydrostatic pressure exists in the upstream face of the dam. This hydrostatic pressure acts on the upstream face and can be equivalently replaced by four concentrated loads, as shown in Figure 13. The hydrostatic pressure gradually increases until the dam breaks. In Figure 14, the results of the Load-COD curve in this study are compared with reference solutions, the experimental data of Carpinteri et al. [32], and the numerical simulation data of Barpi et al. [33] and Shi et al. [34]. It can be seen that the elicited results before the peak, load peak, and experimental data obtained by this method are in good agreement with all the other numerically elicited data. It can also be seen that the post-peak curve of the experiment was significantly higher than the numerical results, which means that the crack opening displacement (COD) of the experiment was greater under the same loading conditions. This phenomenon may be owing to the unanticipated rigid rotation in the experiment of Carpinteri et al. [32] which results in premature failure of the prefabricated crack. Compared with other numerical results shown in Figure 14a, the In Figure 14, the results of the Load-COD curve in this study are compared with reference solutions, the experimental data of Carpinteri et al. [32], and the numerical simulation data of Barpi et al. [33] and Shi et al. [34]. It can be seen that the elicited results before the peak, load peak, and experimental data obtained by this method are in good agreement with all the other numerically elicited data. It can also be seen that the post-peak curve of the experiment was significantly higher than the numerical results, which means that the crack opening displacement (COD) of the experiment was greater under the same loading conditions. This phenomenon may be owing to the unanticipated rigid rotation in the experiment of Carpinteri et al. [32] which results in premature failure of the prefabricated crack. Compared with other numerical results shown in Figure 14a, the initial stiffness obtained by the X-SBFEM method in this study is the closest to that obtained from the results of Carpinteri et al. [32]. The results obtained by Barpi et al. [34] were smoothened before the peak, while the large stiffness elicited in the pre-peak response obtained by Shi et al. [35] resulted in a crack opening displacement (COD) that was smaller in value compared to the experimental results. The crack trajectories obtained by the above experiment and numerical method are shown in Figure 14b. The crack trajectory obtained by the X-SBFEM method in this study is better matched with the experimental and with all the other numerical results. Among them, the crack trajectory obtained from the experiment directly developed towards the dam site before the dam ruptured. At the same time, the trajectories obtained by the numerical simulations of Barpi et al. [33] and Shi et al. [34] horizontally penetrated the dam body, which is quite different compared to the experimental results. The results generated in this study obviously match more closely the experimental results. Figure 15 shows the variation of the crack opening displacement and the distribution of cohesive traction with the crack path in this example model, based on the X-SBFEM method. The crack trajectories obtained by the above experiment and numerical method are shown in Figure 14b. The crack trajectory obtained by the X-SBFEM method in this study is better matched with the experimental and with all the other numerical results. Among them, the crack trajectory obtained from the experiment directly developed towards the dam site before the dam ruptured. At the same time, the trajectories obtained by the numerical simulations of Barpi et al. [33] and Shi et al. [34] horizontally penetrated the dam body, which is quite different compared to the experimental results. The results generated in this study obviously match more closely the experimental results. Figure 15 shows

Static Cohesive Crack Propagation Simulation of Koyna Dam
After a severe earthquake in 1967, the neck of the Koyna Dam suffered serious damages. Gioia et al. [35] (1992) simulated the generated crack based on linear elastic fracture mechanics and Bhattacharjee et al. [36] applied the smeared model to analyze the crack propagation. The geometry and material parameters of the Koyna Dam are shown in Figure 16

Static Cohesive Crack Propagation Simulation of Koyna Dam
After a severe earthquake in 1967, the neck of the Koyna Dam suffered serious damages. Gioia et al. [35] (1992) simulated the generated crack based on linear elastic fracture mechanics and Bhattacharjee et al. [36] applied the smeared model to analyze the crack propagation. The geometry and material parameters of the Koyna Dam are shown in Figure 16

Static Cohesive Crack Propagation Simulation of Koyna Dam
After a severe earthquake in 1967, the neck of the Koyna Dam suffered serious damages. Gioia et al. [35] (1992) simulated the generated crack based on linear elastic fracture mechanics and Bhattacharjee et al. [36] applied the smeared model to analyze the crack propagation. The geometry and material parameters of the Koyna Dam are shown in Figure 16   During the numerical simulation, the crack initially expanded horizontally. When the overflow increased gradually, the crack trajectory gradually expanded downward owing to the increase of the compressive stress in the downstream area of the dam. Figure 17  During the numerical simulation, the crack initially expanded horizontally. When the overflow increased gradually, the crack trajectory gradually expanded downward owing to the increase of the compressive stress in the downstream area of the dam. Figure 17 is the schematic diagram of the ultimate crack propagation path and corresponding cohesive tractions when the overflow reached 10.35 m.  Figure 18 compares the crack paths corresponding to an overflow of 10.35 m using the present method with the results of Zhong [34], Bhattacharjee and Léger [36], and Gioia and Bazant [35]. In the results reported in these published studies, the overflow was approximately 10.2 m, 10 m, and 14 m, respectively. The crack path predicted by Bhattacharjee and Léger [36] was initially nearly horizontal. It then turned downwards when the crack became equal to half of the width of the dam neck. The crack path based on X-SBFEM in this study is consistent with those predicted by Zhong [34] and Gioia and Bazant [35]. These studies reported the formation of very short horizontal extensions that initially curved downwards and towards the dam's heel.   Figure 18 compares the crack paths corresponding to an overflow of 10.35 m using the present method with the results of Zhong [34], Bhattacharjee and Léger [36], and Gioia and Bazant [35]. In the results reported in these published studies, the overflow was approximately 10.2 m, 10 m, and 14 m, respectively. The crack path predicted by Bhattacharjee and Léger [36] was initially nearly horizontal. It then turned downwards when the crack became equal to half of the width of the dam neck. The crack path based on X-SBFEM in this study is consistent with those predicted by Zhong [34] and Gioia and Bazant [35]. These studies reported the formation of very short horizontal extensions that initially curved downwards and towards the dam's heel. During the numerical simulation, the crack initially expanded horizontally. When the overflow increased gradually, the crack trajectory gradually expanded downward owing to the increase of the compressive stress in the downstream area of the dam. Figure 17 is the schematic diagram of the ultimate crack propagation path and corresponding cohesive tractions when the overflow reached 10.35 m.   [34], Bhattacharjee and Léger [36], and Gioia and Bazant [35]. In the results reported in these published studies, the overflow was approximately 10.2 m, 10 m, and 14 m, respectively. The crack path predicted by Bhattacharjee and Léger [36] was initially nearly horizontal. It then turned downwards when the crack became equal to half of the width of the dam neck. The crack path based on X-SBFEM in this study is consistent with those predicted by Zhong [34] and Gioia and Bazant [35]. These studies reported the formation of very short horizontal extensions that initially curved downwards and towards the dam's heel.   Figure 19 shows the relationships between the overflow and the crest displacement obtained herein and four previously published numerical simulations [34][35][36]. It can be seen from Figure 19 that the results of Gioia et al. [35] are consistent with those of Zhong et al. [34], Bhattacharjee and Leger et al. [36], and Li et al. [23] and the method presented herein during the initial loading, before the obvious nonlinearity took place. Subsequently, the results obtained by the method presented in this study are comparatively different from those reported by Gioia et al. [34] and Li et al. [23] based on the linear elastic fracture mechanics method. Nevertheless, they match closely to the results of Zhong et al. [35] and Bhattacharjee and Leger et al. [36]. In Figure 19, the resistance of the dam increases as a function of the crack length and there is no postpeak region owing to the stabilization effect of the self-weight of the dam. The results from these examples also reveal the minor difference encountered between the numerical simulation of the scaled-down dam model and the actual dam.
Appl. Sci. 2017, 7, x FOR PEER REVIEW 18 of 20 Figure 19 shows the relationships between the overflow and the crest displacement obtained herein and four previously published numerical simulations [34][35][36]. It can be seen from Figure 19 that the results of Gioia et al. [35] are consistent with those of Zhong et al. [34], Bhattacharjee and Leger et al. [36], and Li et al. [23] and the method presented herein during the initial loading, before the obvious nonlinearity took place. Subsequently, the results obtained by the method presented in this study are comparatively different from those reported by Gioia et al. [34] and Li et al. [23] based on the linear elastic fracture mechanics method. Nevertheless, they match closely to the results of Zhong et al. [35] and Bhattacharjee and Leger et al. [36]. In Figure 19, the resistance of the dam increases as a function of the crack length and there is no postpeak region owing to the stabilization effect of the self-weight of the dam. The results from these examples also reveal the minor difference encountered between the numerical simulation of the scaled-down dam model and the actual dam.
Meanwhile, it can be noted that the extended scaled boundary finite (NFEM-based) method is numerically stable and robust in modeling the nonlinear postpeak response of the dam up to a state where severe fractures and significant deformation start to occur.

Conclusions
The following conclusions are drawn based on the presented results: (1) A nonlinear X-SBFEM model using the linear superposition of the iterative method was developed and validated to include the cohesive tractions and the fracture energy from FPZ.  Meanwhile, it can be noted that the extended scaled boundary finite (NFEM-based) method is numerically stable and robust in modeling the nonlinear postpeak response of the dam up to a state where severe fractures and significant deformation start to occur.

Conclusions
The following conclusions are drawn based on the presented results: (1) A nonlinear X-SBFEM model using the linear superposition of the iterative method was developed and validated to include the cohesive tractions and the fracture energy from FPZ.