Evaluation Method for Cohesive Crack Propagation in Fragile Locations of RCC Dam Using XFEM

: Roller compacted concrete (RCC) dams own a large number of horizontal construction layers, which can easily lead to weak joints among layers and generate interlayer joints with different scales to reduce the dam bearing capacity. In this study, extended ﬁnite element method (XFEM) is used to simulate crack propagation, the ﬁnite element description is ﬁrst taken on the strong discontinuity. Subsequently, the displacement function of the crack-tip in the quadrilateral element and the geometric determination method of the crack-tip strengthening region are established. Afterwards, the discrete form of the governing equation is derived and the XFEM increment discretization method of the cohesive crack with the crack-tip reinforcement is proposed using the virtual node method to represent the discontinuity of the fracture element. These methods are validated through simulating mixed-mode cracking of one-sided notched asymmetric four-point bending beam. Eventually, the proposed methods are applied to RCC gravity dam to study the development rule and propagation path of the interlayer joints, so as to evaluate the effect of different lengths of the interlayer joints on the dam structural performance. The estimated critical values of dam deformation are helpful to prevent the dam failure during long term operation.


Introduction
Roller compacted concrete (RCC) dams are usually built in thin and horizontal lifts, which allow rapid construction and gain recognition for building new dams and rehabilitating existed dams. Roller compacted layer thickness is usually only a few tens of centimeters, resulting in a large number of horizontal construction layers. Over the past four decades, a number of design technologies and construction methods have been adopted to enhance the final product while maintaining the construction speed. However, the weak bonding surfaces among compacted layers still easily crack with different scales, due to material properties, concrete mixed ratio, construction technology, climate conditions and other uncertain factors. For example, 5 of the 10 horizontal cracks on the RCC dam surface of the Fenhe reservoir two were located at the interlaminar junction (including a transmissible horizontal crack with 170 m length appeared at the 874 m elevation of the shutdown surface in the winter of 1998) [1]. Cracks in horizontal construction joints near the dam crest of the Guanyingge occurred when concrete pouring stopped in the upper and lower reaches from late October to early April of the following year [2]. The interlayer joints could reduce dam bearing capacity and durability and their further expansion would endanger dam safety operation. It is critical to study the development law of the interlayer joints under environmental loads, so as to evaluate the effect of the interlayer joints on dam structural performance.
Fracture mechanics is primarily used to study the fracture mechanism of the weak surfaces among RCC layers, such as traditional finite element and adaptive mesh [3], joint the crack tip in the quadrilateral element and the geometric determination method of the crack-tip strengthening region are proposed. In Section 3, the variation of the continuous component and the discontinuous component of the displacement mode is derived to obtain the discrete form of the governing equation, the increment discretization simulation approach using the virtual node method is proposed as well. In Section 4, the proposed methods are firstly validated by composite cracking simulation on an asymmetric fourpoint beam, then they are applied to the RCC gravity dam to study the propagation path of typical interlaminar joints and evaluate their effect. In Section 5, concluding remarks complete the paper. In comparison with the aforementioned works concerning crack propagation in concrete dams, there are limited efforts have been done for the case in which cracks probably appear at the weakest surfaces of the RCC dam upstream and the dam heel. In this study, the nonlinear fracture behavior of typical weakest interlayers of RCC dam will be studied by the XFEM under excess water level conditions made in ABAQUS.

XFEM Description of the Strong Discontinuity
The XFEM allows the discontinuous surface to pass through the element [9,10], i.e., the mesh is independent on the discontinuity, hence it is necessary to describe the discontinuities geometrically. The primary method is the level set method. The level set function is taken to construct an extended shape function in the XFEM. Here, the discontinuity on both sides of the fracture is defined as the strong discontinuity and the displacement field is discontinuous. In the calculation, a certain form of the extended shape function can be used to describe the discontinuity. As shown in Figure 1, a crack in the computational domain terminates an element through the mesh. Penetration element nodes, crack-tip element nodes are recorded as L (represented by " ") and K (represented by " "), respectively. Different extended shape functions are used for the node sets L and K. (1) Crack penetration element: As the displacement field on both sides of the fracture surface jumps, the extended shape function ψ j (x) [10] is where H(x) is the Heaviside function, The level set function f (x) [10] is: where n + is the element outside normal vector of the fracture atx. For points not on the crack, f (x) is the shortest distance from x to the crack. If the position of x is consistent with the direction n + , take positive; conversely, take negative.
(2) Crack-tip element: For elements around the crack tip, the extended shape function ψ j (x) [10] is where Φ l (x) is used to simulate the discontinuous displacement field at the crack tip. It is the basis of the Westergaard field at the asymptotic crack tip in the plane composite fracture of linear elastic fracture mechanics.
In Equation (4), only the first function is discontinuous when the crack traverses and the others are used to improve the singularity of the solution near the crack tip, where (r, θ) represents a local polar coordinate. The crack-tip function cannot only make the crack terminate precisely in the element, but also obtains better accuracy in the relatively rough 2D finite element mesh [29,30]. When Equation (4) is used to construct the crack-tip shape function, it can represent the geometric discontinuity and the displacement discontinuity of the elements near the fracture. The distribution state of the crack-tip displacement field can also be accurately obtained.
Based on the above two shape functions, the node displacement vector can be divided into the conventional displacement and reinforced displacement, caused by the cracks. When the Heaviside function and the asymptotic crack-tip displacement function are used to strengthen the nodes of the crack penetration element and the crack-tip element, the approximate form of the displacement vector in the XFEM discontinuous displacement mode [10] can be rewritten as follows: (5) where N i (x) and N j (x) are the finite element shape functions. u i , a j and b lk are the node normal freedom degree and the strengthened freedom degree. I represents a set of all nodes. J represents a set of penetration element nodes. K represents a set of nodes of the crack-tip element.
In Equation (5), the first item on the right is the same as the conventional finite element method, while the second and the third are unique to the XFEM, called enhancements, to describe the existence of cracks. If there is no crack penetrating the element, the displacement field of the element will be the first term, or the first two terms if the element is penetrated, or the sum of the first and third terms when the element locates at the crack-tip.
and Φ l (x) can improve the original extended finite element model, since the calculated node displacement does not need to be modified by the displacement formula.
For the strengthening node of the crack-tip branch function, if only the node of the crack-tip element is selected, the size region of the element tends to 0 and the strengthening region of the crack-tip branch function also tends to 0. In order to improve the convergence rate, the element-based reinforcement area is replaced by the geometry-based reinforcement area. Due to the close relationship between the calculation accuracy and the crack-tip position, the crack tip is suggested to be the center of a circle with radius r to eliminate the influence of the crack-tip position on the accuracy, as shown in Figure 2. The nodes in the circle are strengthened by the crack-tip function.
where S is the element area where the crack tip is located, r enr = 3 [31]. The thick red line in Figure 2 is the crack. " " represents the strengthening node, increasing one freedom degree in the direction of each normal freedom degree. " " indicates the crack-tip strengthening node, increasing four freedom degrees in each direction of the normal freedom degree. "×" is the double strengthening node, increasing five freedom degrees in each direction of the normal freedom degree. Obviously, the nodes are strengthened in a large circle centered with the crack tip.

Governing Equation Derivation
The cohesive fracture zone is shown in Figure 3, whose calculation domain is Ω, including the crack zone Γ c , making the displacement discontinuity on both sides. According to the stress characteristics of Γ c , it can be divided into a fracture zone Γ d (the displacement on both sides is independent and does not transfer the stress) and a cohesive zone Γ coh (the displacement is discontinuous but the stress is transmitted). In Figure 3, Γ F and Γ u are the external force boundary and the displacement boundary. The stress transfer follows the constitutive law of the cohesive zone. The governing equations and the boundary conditions in the region Ω are: where ∇ is the gradient operator, σ is the stress vector, Ω is the domain of the body, n is the unit outward normal vector on Γ F , F is the external force, t + and t − are cohesive forces in Γ + coh and Γ − coh , n + and n − are outside normal vectors on Γ coh . Under the conditions of small strain and displacement, the kinematic equations are: where ε is the strain vector, u is the displacement vector and w is the crack opening displacement. The constitutive equation is: where t is the cohesive force on Γ coh and C is the 4th order Hooke tensor. In order to avoid the simulation singularity, the crack tip is set at the element boundary. Assuming that the crack propagation is carried out in the element, each propagation would penetrate the crack-tip element. There are only two kinds of elements in the model: one is an ordinary element; the other is penetrated by the crack. When the mesh is subdivided and the length-width ratio of the mesh is close to 1.0, the simplified method has little effect on the calculation results. Since the displacement field of the cohesive crack tip is no longer singular, Φ l (x) can be: The displacement of the cohesive cracks can be: where x is the coordinate vector. u, a and b are the element normal freedom degree vector, H function enhances the freedom vector and the crack-tip strengthened freedom degree vector, respectively. u T is the real displacement vector after synthesis. N is the function array. H and φ are the Heaviside function array and the crack-tip strengthening function array of the element nodes, respectively. The weak form description of the problem is given to facilitate the derivation and the displacement u T belongs to the kinematic allowable displacement space µ.
The virtual work equation under the action of no physical force is: where t is an external force acting on the boundary. According to the principle of the virtual work, the variation of the continuous component and the discontinuous component of the displacement mode can be obtained.
Let w(v) denote the crack opening in the global coordinate system, w x−y = (N + coh − N − coh )v, N + coh and N − coh are corresponding shape functions of the cracks Γ + coh and Γ − coh , for crack penetration elements, N + coh = N, N − coh = -N. The N + coh and N − coh of the crack-tip elements are expressed as Nφ, φ on both sides of the crack is equivalent, but reverse. In order to obtain the virtual work by the cohesion force, w x−y is converted to w l in the local coordinate through the transformation matrix M, that is w l = Mw x−y .
where s is the area of one-sided crack surface on the unit length during the crack propagation, t represents the crack surface cohesion and is a function of the crack opening. The incremental formt can be expressed as: where T is the material matrix and M is the transformation matrix that converts the global coordinate to the local. The incremental form of the stress and the strain is: where B u , B a and B b represents strain matrices corresponding to the normal freedom degree. Substituting Equations (17)- (19) into Equation (16), written as an incremental form, the discrete form of the governing equation is: where: where u and u represents the general part and the strengthening part of the nodes, f int u and f int u are the nodal internal forces, hence the right side of Equation (20) denotes the unequilibrium force of the nodes. Since the internal tensile stress of the crack is continuous, For the crack penetration elements: For the crack-tip elements: where B r i (r = u, a, b) represents corresponding strain matrix and i is the node number.
The integral load array is integrated by the element load arrays.
where (f i ext ) u is the load of the node normal freedom degrees, (f i ext ) a and (f i ext ) b are additional loads of the H function strengthening node and the crack-tip branch function strengthening node, respectively.
where F b , F s and F are the distributed body force, the distributed surface force and the concentrated force, respectively. The partial derivative of the shape function is: In order to obtain the partial derivative of the crack-tip strengthening function in the global coordinate (x, y), the partial derivative of the crack-tip local coordinate (x 1 , y 1 ) must be firstly obtained, then the coordinate transformation is carried out as shown in Figure 4. According to Equations (33) and (34), the partial derivative of the crack-tip strengthening function to (x 1 , y 1 ) can be obtained as follows.
Take coordinate conversion on Equation (35) to obtain

XFEM Increment Discretization Simulation
Since concrete has obvious nonlinear fracture characteristics, the crack propagation can be simulated by the virtual crack model [32]. The nonlinear constitutive model can be used to reflect the strain softening and the extended finite element can be constructed by the presupposition virtual node method. As shown in Figure 5, the virtual node method superposes the nodes on the initial real nodes to represent the discontinuity of the fracture element. In the cohesive crack model, the nonlinear closure force near the crack tip is equivalent to the virtual crack line, hence the nonlinear deformation near the crack tip can be transformed into the whole crack line. The constitutive model of the cohesive cracks describes the relationship between the nonlinear closure force in the virtual fracture zone and the relative displacement of the fracture surface. The model assumes that the macrocrack occurs when the cohesive force of the crack surface equals to 0 and the area enclosed by the distribution curve of the cohesive force in the virtual fracture zone and the coordinate axis is equal to the fracture toughness.
Due to the asymmetry of the practical structures and the bearing loads, the crack forms are diverse. Since the stress of the crack tip cannot be accurately known, the stress of the Gaussian point around the crack tip is obtained by weighted average. The weight function W is where r represents the distance between the known node and the crack tip and l is used to control the speed of W fading from the crack tip to the known node, usually 3 times the size of the typical element. The traction separation model is used to simulate the constitutive relationship, i.e., assuming that concrete is initially linear elastic, the initial damage and the evolution will occur after the crack is initiated. The extended element will not be damaged under the compression. According to concrete tensile tests, the nonlinear softening constitutive relationships are more consistent with the actual situations and can be used to simulate the traction separation response. When the maximum main stress at the crack tip equals the tensile strength, the introduced crack is always perpendicular to the direction of the maximum main stress. However, the crack propagation is not bound to the element boundary, it is assumed that each crack propagation needs to break through the complete element, i.e., the crack tip coincides with the element boundary to avoid the singularity simulation.
An incremental iteration process for the cohesive fracture expansion is shown in Figure 6. The key of the iterative algorithm is the expansion of new crack-tip element. The Newton-Raphson orthogonal residual method is used in the incremental load step and the energy criterion is regarded as the convergence criterion, where the strain energy of the initial elastic load step E ref stays as the convergence parameter. In order to make the crack expand smoothly, the stress iteration is added in each load increment step to ensure that the maximum main stress at the crack tip is equal to the tensile strength.
The approach to determine the length of the virtual crack zone is explained as follows. 1 The initial length of the crack (including the virtual crack length) is given and the cohesive force of the crack surface is not included before initial cracking. The displacement and the opening degree of each point on the crack surface are solved. 2 The iteration starts from the second step and the cohesive force at the corresponding position is calculated according to the opening degree of each point. Keep iterating to obtain the true opening degree on the crack surface until the energy convergence criterion is satisfied. 3 The exponential softening constitutive model is adopted for the relationship between the cohesive force and the crack opening, i.e., the cohesive force on the crack surface decreases exponentially with the increase of the crack opening. When the cohesive force is less than a certain value, the increase of the crack opening has little effect on it, thus a critical value of the normal cohesive force can be set. When the cohesive force is less than the value at a point, it is considered to belong to a macrocrack zone. Then the opening value of the crack tip will be the critical value of the cohesive force. Through the trial calculation, the critical value was set as 6.5 × 10 2 N and corresponding critical opening value was 5.25 × 10 −4 mm. 4 Coordinates of the node on the crack surface corresponding to the critical opening value can be calculated through the interpolation of the opening degree of known integral nodes. The distance between the node and the virtual crack tip will be the length of the virtual fracture zone.

Validation by the Composite Cracking Simulation
The proposed methods were validated by one-sided notched asymmetric four-point bending beam, experimental test completed in 1998 by Gálvez, etc. [33]. The numerical simulation of cracking under mixed loading of the I-II type was carried out in ABAQUS. The geometric form, loads and boundary conditions of the bending beam were shown in Figure 7. The thickness of the specimen was 50 mm, the compression strength f c was 57 MPa, the tension strength f t was 3.0 MPa, the failure energy G f was 69 N/m and the elastic modulus E was 38 GPa. The numerical model was built as shown in Figure 8. The concrete beam adopted C3D8R element and the mesh was refined in crack propagation, support and loading point regions. The initial crack surface was preset according to the crack location in the test and there was no need to mesh the crack. The initial crack module should be kept at a small distance from the beam element boundary when the model was assembled as the initial crack included in Figure 8. In order to obtain a more accurate crack propagation path, the local mesh refinement was carried out for the crack propagation region through selecting the greater stress position and the lengthwidth ratio of the mesh was close to 1.0 to reduce the dependence of results on the mesh size. The model was loaded according to the test, controlled by the crack mouth opening displacement (CMOD), i.e., 0.004 mm/min to 40% of the peak value and 0.08 mm/min to the end of loading and the conversion of load-CMOD and load-displacement curves. Uniform load 16 MPa (P = 8 kN) was applied to the model, then the load-crack mouth opening displacement (P-CMOD) curve and the crack development path were obtained using the linear and exponential strain-softening relations on the crack surface, as shown in Figure 9. The crack propagation process was shown in Figure 10, the four-point beam notch deformed along the horizontal and vertical directions, where the I-II failure pattern was more obvious.  The displacements of the beam notch were in agreement with test results, as shown in Figure 11. In Ref. [34], 100 pairs of nonlinear spring elements were adopted to simulate the crack along the default crack surface and multiplied the element stress by the element area to output the load according to the crack development path. In this study, the XFEM did not need to presuppose the crack path and would be more effective in practical structures.

Application on RCC Dam
The RCC gravity dam, as shown in Figure 12, located in the main stream of Tingjiang River, whose maximum dam height is 111 m, officially started construction in April 1998.

XFEM Model and Material Parameters
Due to various influence factors such as construction measures, equipment capacity, climate conditions, etc., there were defects in the plane of 40 m away from the foundation during construction, while the layer at 120 m elevation was the interface between two materials R 180 150 and R 180 100 and turned into the weakest surface among the most unfavorable layers [35,36]. The interlayer joints also easily cracked in the interfaces between the dam concrete and the foundation. Strictly speaking, each RCC layer (0.3 m thickness) should be used as a block element mesh to simulate the layered characteristics of the RCC dam, but it was bound to greatly increase the computational workload. Hence, this study only focused on the two weakest surfaces, as shown in Figure 13. The modeling range of the foundation is shown in Figure 13, where the influence of the impervious curtain and the drainage holes was not considered. Fixed constraints were imposed at the bottom boundary and normal hinge constraints were imposed on the front and back boundaries. The model was primarily composed of four-node isoparametric elements. The mesh size of the local refined area was 1 m × 1 m. Convergence tests were implemented to determine the number of elements and nodes. The whole model consisted of 6762 elements and 6944 nodes, including 2533 elements in the dam body and 4229 elements in the foundation.
Mechanical parameters of dam concrete and rock mass are listed in Table 1, according to the experimental results [35]. The fracture energy is given according to the test results [37], G f = 100 N/m.

Material Constitutive Relationships
The linear elastic traction-separation model was used in the dam crack area, i.e., the concrete was linear elastic in the compression and would maintain in the initial tensile stage. When the tensile stress exceeded the admissible strength, the initial damage occurred and the virtual crack model based on the XFEM would be used. The contact between crack surfaces was rigid and the separation after contact was allowed, as shown in Figure 14. The strain softening curve of the crack surface is shown in Figure 15.  The appearance of cracks means that the cohesive force of the element begins to decline, i.e., the initial damage criterion would be the initial criterion of the crack. When the fracture criterion f reaches 1.0 (within a certain error range), a new crack will be introduced and the existed crack length will be enlarged in the next incremental step. The error range of f is: where f tol means the tolerance, f tol = 0.05. If f > 1.0 + f tol , the increment step will be reduced to satisfy the initial criterion. When the maximum main stress criterion is adopted, the fracture criterion f can be: where σ 0 max represents the maximum admissible main stress. The symbol is the Macauley bracket to ensure that the pure compressive strain does not damage the material. When f satisfies Equation (39), the damage occurs.
The overall average damage degree of the crack through the element boundary is denoted by D and determined by the linear criterion. For the linear softening process: where δ T 0 e f f is the effective tractive force, G c indicates the fracture energy, δ 0 m is the effective deformation, δ max m is the maximum effective deformation during loading and δ f m is the effective deformation of a complete fracture. The normal stress and the tangential stress affected by the damage are where T n and T s are the normal stress and the tangential stress in the case of the linear elasticity.

Excess Water Level Method
The imposed loads contain concrete gravity, upstream water weight, upstream water pressure, uplift pressure and seepage pressure in the joints. In order to calculate the propagation of the interlayer cracks, the excess water level method K H (trapezoidal overload) was used, as shown in Figure 16, increasing the water level step by step until the crack expands unsteadily.
where H represents the loading water level and H 0 is the design water level.

Propagation Paths of the Interlayer Joints
The model calculation consists of four steps to ensure the accuracy and prevent the non-convergence. The initial iteration step length was 0.001, the minimum step size was 1 × 10 −20 , the maximum step length was 0.1, the maximum number of the iterations was 10,000 and the nonlinear control switch was turned on in the solver.
Step 1 was the initial in-situ stress equilibrium analysis to eliminate the deformation caused by the gravity action of the rock mass. The deformation was about 10 −30 order of the magnitude at the end of the step. The dam gravity was imposed in step 2 and the upstream water pressure was imposed in step 3. In step 4, the water pressure was imposed on the crack surface with the magnitude of the water pressure p = γh at one end (γ is the bulk density and h is the depth of the water) and the other end 0, assuming a triangular distribution of the water pressure in the initial section of the crack surfaces. The water pressure was transformed into the equivalent joint force of the crack penetration element by the Gaussian integral [38], while the change of the seepage pressure caused by the water flow in the short crack during the propagation was ignored.
In the output of the field variables, the PHILSM (Signed distance function to describe the cracksurface, normal level set ϕ), the PSILSM (Signed distance function to describe the initial crack front, tangential level set ψ) and the STATUSXFEM (Status of the XFEM element) variables were selected. The XFEM element state: the STATUSXFM value is 1.0 when the element is completely penetrated, while the value is 0.0 if there is no crack. When the element is partially broken through, the STATUSXFE value is (0.0, 1.0). Hence, the crack location can be traced in the postprocessing and whether the element is strengthened or not can be checked.
1. The weakest interlayer surface located at 120 m elevation: The initial depth of the interlayer joint was 4 m (10% of the plane width). The XFEM calculation results show that the initiation cracking load of the weak surface was 116 m water head, about 1.17 times of the dam height, K H = 1.17. Figure 17 shows the morphology of the interlayer crack. When the crack occurred, the fracture surface became a free surface and only the part connected to the dam body was constrained. 2. The weak interlayer surface located in the dam heel: The initial depth of the interlayer joint was 2 m (2.74% of the plane width). The initiation load was 95 m water head. The crack growth pattern of the dam heel is shown in Figure 18. Figure 19 shows the gradual process of the crack propagation. The overall crack propagation trend was similar, i.e., the crack approached the horizontal at the initial stage and inclined downward to the dam bottom with the increase of the water head.

Effects of the Interlayer Joint Propagation
In order to assess the effect of the interlayer cracking on the dam structural behavior, taking the weakest surface located at the dam heel (initial crack length 2 m) for illustration, during the crack propagation, the variation of the horizontal displacement and the vertical displacement of the dam crest, the shear displacements of the interlayer joint and the opening displacement of the dam heel interlayer crack are listed in Table 2. The horizontal displacement increased obviously when the initial length of the interlayer joints extended from 2 to 3 m. Subsequently, both the horizontal displacement and the vertical displacement of the dam crest increased gradually. The shear displacement and the opening displacement increased linearly with the step-by-step propagation of the interlayer crack, basically synchronous. The measured time series of the dam crest in the recent 7 years are shown in Figure 20, where positive values represent that the dam deformed downstream and rose, otherwise it deformed upstream and sinks. The measured deformation actually includes the influence of environmental load factors such as water pressure and temperature, complex geological conditions, the influence of reservoir water permeability, etc. Although they indicate annual periodic variation of the dam horizontal and vertical deformations during normal operation, the slightly increasing trends to the downstream and the uplift should include the influence of interlayer cracks. Additionally, they will become more and more significant with the increase of the dam age.  The relationships between the horizontal and the vertical displacement and the length of the interlayer joint were constructed through the polynomial regression, as shown in Figure 21. The deformation of the dam crest can be estimated from the length of the interlayer joint. When the length was 9 m, the interlayer crack extended to the curtain, at this time, the horizontal displacement was 8.3 mm and the vertical displacement was 1.29 mm. For RCC gravity dams, the further propagation of the cohesive cracks at the dam heel would damage the grouting curtain. Generally speaking, most of curtain centerlines were located at the distance from the upstream dam heel 0.06 H 0 -0.10 H 0 (H 0 is the maximum water depth, calculated from the bedrock at the normal water level), or about 0.1 times the width of the dam bottom. Therefore, the center line of the curtain was 9 m from the dam heel in the RCC dam. When the crack depth in the dam heel reached 9.73 m, it had extended to the centerline of the curtain. In this case, the penetrated water in the joint would increase the uplift pressure on the foundation surface to weaken the dam stability, thus the remedial measures must be taken to prevent the RCC dam from losing its stability during long term operation.

Conclusions
This study was devoted to evaluate the effect of cohesive crack propagation in the fragile locations of RCC dams with the XFEM. The following conclusions can be drawn.

•
The established displacement function of the crack tip in the quadrilateral element and the geometric determination method of the crack-tip strengthening region were demonstrated to be more suitable for solving the cohesive cracks. The generalized Heaviside function used in the XFEM, similar to the exact local mode, could improve the adaptability of the displacement approximation function. Moreover, the adopted virtual node method could overcome the deficiency of distinguishing the type of inclined crack element by the level set.

•
The computational results for the cohesive cracks through the XFEM increment discretization simulation were in agreement with the one-sided notched asymmetric four-point bending beam test results by Gálvez. In contrast to the modeling results by Cendón, the load-crack opening displacement curve and the crack development path in the loading stage were very similar. Hence, the proposed method simulating the cracking process without presupposing the crack path was proved to be feasible. The drawback is that a variable number of the freedom degrees per node are needed at the cohesive cracks in the presented method.

•
The proposed methods were applied to simulate the cohesive crack propagation of the two typical weakest interlayer surfaces in the RCC gravity dam. In order to consider the influence of the plastic zone near the crack tip, the crack end was limited to the element edge and different models were used to simulate the constitutive relations. The localization caused by the softening of the crack surface could be simulated well by an overall average damage model. The gradual propagation and the morphology of the interlayer cracks obtained under the excess water level conditions indicate the crack approached the horizontal at the initial stage and inclined downward with the increase of the water head, which can make the dam deform downstream and rise.
Due to the construction technology of RCC dams, the interlayer defects are objective and unavoidable. It is vital to predict their possible propagation paths through numerical simulation for controlling the development of harmful cracks in time. Based on the regression results presented, the estimated critical values of the dam crest deformation can provide scientific support for the remedial measures to prevent the dam instability and failure during long term operation. In addition to studying the influence of single crack, the future study will focus on the stability assessment of multiple random interlayer cracks for RCC dams.

Data Availability Statement:
The data presented in this study are available in [Evaluation method for cohesive crack propagation in fragile locations of RCC dam using XFEM].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: RCC Roller Compacted Concrete XFEM Extended Finite Element Method