A Boundary Element Formulation for Crack Analyses Incorporating a Cohesive Zone Model

The surfaces of cracks behind the (fictitious) crack tip are not completely separated in some materials, such as concrete, brittle polymers, fiber-reinforced composites, tough ceramics and some alloys. The tractions can be transferred across the crack line along a relatively long extension of the crack, which is commonly called the wake zone, the bridging zone, or the cohesive zone. The main assumption is the occurrence of material softening beyond the peak load in a narrow layer with a negligible volume behind the fictitious crack tip, in which the action can be replaced by cohesive forces (Figure 1). Two types of constitutive laws for the material in the cohesive zone have been used in the literature. This study employs a trac-tion-displacement relationship in the cohesive zone instead of using a material constitutive law that is defined in terms of stresses and strains accompanied by a layer thickening law. Barenblatt [1] introduced a cohesive model that employed a fictitious crack. Hillerborg et al. [2] proposed a function for a softening model related to an opening (mode I) crack (Figure 2), which allowed finite element analyses of the problem to be performed, such as those by Petersson [3], Carpinteri [4], and Rots [5].

The boundary element method (BEM) is superior to the finite element method (FEM) in crack propagation analysis because remeshing of the domain is not necessary when the crack grows, and new elements can be introduced without affecting the existing elements. The coincidence of two crack surfaces requires two different boundary integral equations (BIEs) for the solution. The dual boundary element method (DBEM) is a technique that is employed for crack propagation analysis. The position of the collocation point to solve the traction BIE and the strategy to treat improper integrals are the essential features of the formulation. The tangential differential operator (TDO) in conjunction with integration by parts is an interesting procedure for reducing the order of singularity in fundamental solution kernels of traction BIEs when Kelvin-type fundamental solutions are employed. Kupradze [7] first presented an application using the TDO, and Sladek [8] employed the TDO in a curved crack solution. Bonnet [9] presented regularized formulations for the BEM employing the TDO for gradients in potential problems and in stress BIEs for elasticity problems, including fracture mechanics formulations. Additional care is necessary in formulations with TDO when non-conformal interpolations are used on boundary elements, as was shown in [10] for traction BIEs in plane elasticity problems, in [11] for traction BIEs in three-dimensional elasticity and in [12] for stress equations in plate bending with the effect of shear deformation.

This study employed isoparametric linear boundary elements with all of the nodal parameters positioned at the ends of the elements. The collocation points were placed or shifted to the interior of boundary elements and were always located along the boundary line. Continuous or discontinuous boundary-element types were used based on the geometry or loading conditions.

The dual equations of the method are the displacement and the traction BIEs, which are written below for a collocation point that is located on a smooth boundary:

with,

$$\begin{gathered}C_{a k i m}=\mu\left(\frac{2 v}{1-2 v} \delta_{a k} \delta_{i m}+\delta_{a i} \delta_{k m}+\delta_{a m} \delta_{k i}\right) \\ \mathrm{D}_{b m}[f(x)]=n_b(x) f_{, m}(x)-n_m(x) f_{, b}(x)\end{gathered}$$

where $\mathrm{U}_{\mathrm{ij}}\left(\mathrm{x}^{\prime}, \mathrm{x}\right)$ and $\mathrm{T}_{\mathrm{ij}}\left(\mathrm{x}^{\prime}, \mathrm{x}\right)$ are the displacement and traction, respectively, in the direction j at boundary point x due to a singular load in the direction i at collocation point $\mathrm{x}^{\prime}$ based on the Kelvin solution for two-dimensional problems, $\mathrm{u}_{\mathrm{i}}(\mathrm{x})$ and $\mathrm{t}_{\mathrm{i}}(\mathrm{x})$ are the displacement and traction at boundary point x , respectively, $\mathrm{n}_{\mathrm{a}}\left(\mathrm{x}^{\prime}\right)$ is the direction cosine of the outward normal at collocation point $\mathrm{x}^{\prime}, \mathrm{C}_{\mathrm{akim}}$ is the Hooke tensor for an isotropic material, $\sigma_{\mathrm{ibj}}$ is the stress $\sigma_{\mathrm{bj}}$ at boundary point $x$ due to a singular load in the direction i at collocation point $x^{\prime}, D_{b m}()$ is the tangential operator, $\mu$ is the shear modulus, $v$ is Poisson's ratio and $\delta_{i j}$ is the Kronecker delta.

In Eq. (1), the left member has a $1 / \mathrm{r}$ singularity, and the right member has a logarithmic singularity when the field point approaches the collocation point. Both integrals in Eq. (2) have $1 / \mathrm{r}$ singularities because the first employs the tangential differential operator. To solve general mixed-mode crack problems with a single domain formulation, the displacement BIE (Eq. (1)) is applied to one of the crack surfaces, and the traction BIE (Eq. (2)) is applied to the other surface. Although the integration path is still the same for coincident points on crack surfaces, the respective BIEs are now distinct. The collocation points must be positioned to satisfy the continuity requirements for each BIE. The continuity of the displacement function at $\mathrm{x}^{\prime}$ is the necessary condition for the displacement BIE, and it is satisfied when the collocation point is placed at the ends of the boundary element or inside the element. The continuity of the displacement derivative at $\mathrm{x}^{\prime}$ is required for the traction BIE, and it is satisfied when the collocation point is placed inside the boundary element [13, 14].

The cohesive zone is an extension of the crack in which material softening beyond the peak load is located in a narrow layer behind the fictitious crack tip (Figure 1). The cohesive zone is modeled as a crack region containing springs that connect coincidental surfaces. The points that were originally coincidental on opposite sides of the crack line separate into distinct points and are connected by the cohesive zone material (springs). Continued strain increases the separation between these two points and eventually leads to cracking. Several types of functions have been presented in the literature to represent the traction-displacement relation, which should be suitable for quasi-static loading. The linear function that was adopted in this study to represent the softening in the cohesive zone was chosen to allow comparisons with the results obtained in [15].

The $p_y$ value for the normal traction in Figure 3 corresponds to the limit without opening displacement in the cohesive zone, and $w_f$ is the opening displacement limit where the traction disappears. The simplified maximum principal stress criterion is used to determine the crack increment. When the maximum principal stress at a fictitious crack tip reaches the critical value of $\mathrm{p}_{\mathrm{y}}$, the tip will extend under further loading using the constitutive law shown in Figure 3. The tip advances in the direction perpendicular to the direction of the maximum principal stress at that point, and the extension is such that the maximum principal stress at the new tip position is equal to the critical $p_y$ value during continued loading. It is important to note that structural instability may occur while the crack is advancing; i.e. extension of the crack may increase the stress at the crack tip rather than release it.

The linear function shown in Figure 3 is given by:

Figure 4 shows a crack region that contains a cohesive zone. The external boundary line of the problem is $\Gamma_{\mathrm{e}}$, and the tangent direction that is employed to compute the boundary integrals is shown. $\Gamma_1$ and $\Gamma_4$ are portions of the crack surface without cohesive forces (open crack), and $\Gamma_2$ and $\Gamma_3$ are portions that contain cohesive forces. The cohesive zone has tractions on the crack surfaces with opposite signs and the same magnitudes. The dual equations for the problem shown in Fig. 4 can be simplified when the directions and magnitudes of the cohesive tractions and the directions of integration on crack surfaces $\Gamma_2$ and $\Gamma_3$ are taken into account, which yield:

There was introduced $\mathrm{t}^{\mathrm{c}}$ in Eqs. (4) and (5), which is the traction in the direction j on one of the crack surfaces in the cohesive zone. The index $c$ was used because surface $\Gamma_2$ or $\Gamma_3$ with tractions $\mathrm{t}^2$ or $\mathrm{t}^3$, respectively, can be used. $\Gamma=\Gamma_1+\Gamma_2+\Gamma_3+\Gamma_4+\Gamma_{\mathrm{e}}$, and $\Gamma_2=\Gamma_3=\Gamma_{\mathrm{c}}$.

The tractions on one of the crack surfaces in the cohesive zone ($\mathrm{t}_{\mathrm{c}}$) are the unknowns in the system of equations that is obtained from Eqs. (4) and (5), and the cohesive law is the additional equation that is used. Thus, the tractions in the cohesive zone are directly computed at each loading step in the incremental process.

Eq. (6) summarizes the system of equations where the dual equations, according to Eqs. (4) and (5), are converted to submatrices $\mathrm{H}_{\mathrm{ij}}$ and $\mathrm{G}_{\mathrm{ij}}$, whereas the cohesive law is converted to submatrices $\mathrm{A}, \mathrm{B}$ and C . The indices $\mathrm{e}, \mathrm{o}$ and z are related to the boundary portions and correspond to the external boundary, open crack boundary and cohesive zone, respectively. Submatrices A and B relate the displacements and tractions in the normal and tangent directions, respectively, on the crack surfaces in the cohesive zone.

This softening criterion works with tractions and openings in the normal direction, which means that the tractions in the tangent direction are zero in the cohesive zone. The numerical algorithm is summarized below considering the boundary portions shown in Figure 4:

(a) The normal tractions $\left(t_c\right)$ are less than or equal to $p_y$; i.e. no opening occurs in the cohesive zone, and opposite points on the crack surfaces have the same displacement. Submatrices B and C are zero, whereas submatrix A contains direction cosines that relate the displacements in the directions $\mathrm{x}_{\mathrm{i}}$; i.e. the openings in the normal and tangent directions are zero:

(b) Subsequent loads after the normal traction has reached $\mathrm{p}_{\mathrm{y}}$; opening occurs in the cohesive zone according to the cohesive law (Eq. (3)). The submatrices B and C are modified to introduce the cohesive law in the line that is related to the normal direction, whereas the tractions in the tangent direction become zero (the corresponding line in submatrix A becomes zero), and Eqs. (8) and (9) are replaced by:

(c) The opening reaches $\mathrm{w}_{\mathrm{f}}$ in subsequent loads; the tractions $\left(\mathrm{t}_{\mathrm{c}}\right)$ must be zero at this point in the cohesive zone according to Eq. (3). The absence of tractions in the normal and tangent directions causes the tractions in the direction $\mathrm{x}_{\mathrm{i}}$ to be zero. The tractions $\left(\mathrm{t}^{\mathrm{c}}{ }_{\mathrm{i}}\right)$ at this point on the crack surfaces are eliminated from the system of equations:

The boundary element code was derived from that used in [10]. Linear mapping functions were used to represent the displacements and tractions in the boundary elements. The same mapping function was used for conformal and non-conformal interpolations with nodal parameters positioned at the ends of the elements. The collocation points were always positioned on the boundary line at the position ( $\xi$ ') in the range $(-1,1)$ : (i) $\xi^{\prime}=-0.67$ for continuous elements, and (ii) $\xi^{\prime}=-0.67$ and $\xi^{\prime}=+0.67$ for discontinuous elements. Analytical expressions were used to evaluate the singular integrals in the Cauchy principal value sense, and the Gauss-Legendre scheme was used for regular integrals.

A single edge crack in a rectangular specimen was studied under uniform far-field tensile loading to consider displacement-controlled elongation (Figure 5). The specimen was a rectangle with a length of 1 and a height of $h$, and both were equal to 1 in the simulations. All the lengths are normalized by l ; while this is not a natural length scale for fracture problems, it is convenient and easy to interpret. The initial length of the single edge crack ($\mathrm{a}_0$) had values in the range ($0.03,0.58$), Poisson's ratio was 0.3 , and plane strain conditions were adopted according to [15].

As was discussed for the cohesive model, the physical crack tip was positioned at the end of the initial length $\mathrm{a}_0$, where the opening could be $\mathrm{w}_{\mathrm{f}}$, and the fictitious crack tip was positioned at the end of the cohesive zone ZC, where the normal traction value could be $p_y$ (or lower than $p_y$). The parameters $\mathrm{p}_{\mathrm{y}}$ and $\mathrm{w}_{\mathrm{f}}$ were equal to 0.01 and 0.001, respectively. It should be noted that all the stresses and tractions were normalized by the shear modulus $(\mu)$. The unloaded specimen was initially in a stress-free state and was loaded incrementally perpendicular to the top and bottom boundaries under displacement-controlled elongation. The loading displacement increment was taken to be 1E-5 in tension, and a limit of 2000 load steps was used. The element length of the rectangular boundary was 0.05 . Double nodes were introduced at the corners and at the fictitious crack tip. The element length at the initial crack length $\left(a_0\right)$ was 0.05, and the element length in the cohesive zone (ZC) was 0.025. The analyses were carried out until the first element on the cohesive zone was cracked; i.e. the initial crack length was increased by 0.025, which is the length of one element. Figure 6 compares the results with those obtained in [15].

The fictitious crack tip first starts to advance during the loading process, and the physical crack tip also starts to move with continued loading. The length of the initial crack ($a_0$) introduced instability in the results presented in [15], and the stresses at the point ahead of the fictitious crack tip were greater than $p_y$ even when many cohesive elements were used over the length where the crack extension procedure was applied. According to the authors of [15], the critical load for the onset of instability decreases with increasing initial crack length. The crack extension in [15] became stable with initial crack lengths greater than 0.48.

The issues noted in [15] did not appear with the formulation that was used in this study. The loads due to elongation (Figure 6) were different only slightly when the initial crack length $\left(\mathrm{a}_0\right)$ was less than 0.38 , and the curves nearly coincide for larger values. The behavior of the normal traction in the cohesive zone with this formulation was similar for all values of the initial crack length, and they did not differ as was described in [15]. Figure 7 shows the cohesive forces for two initial crack lengths.

The results that were obtained with the proposed formulation were similar to those in [15] and did not cause the issues that were noted in that study. The behavior of the proposed formulation was also similar for other problems; they are described in [16]. The direct computation of tractions in the cohesive zone required fewer loading steps than were required when the iterative procedure was used to obtain tractions based on the initial value for the same cohesive law [17].