Enriched degree of freedom locking in crack analysis using the extended finite element method

Abstract

In this article, the enriched degree of freedom locking that can occur in a crack analysis with the extended finite element method is described. The discontinuous displacement field formulated by the enriched degree of freedom in the extended finite element method does not activate due to the enriched degree of freedom locking. Using the phantom node method, the occurrence of locking when two adjacent elements are simultaneously cracked in a loading step was verified. Two adjacent cracks can be determined to have developed simultaneously when an analysis model reveals a relatively uniform stress distribution on two adjacent elements. Numerical examples of a simply tensioned bar and a reinforced concrete beam are presented to demonstrate the erroneous analysis result due to the enriched degree of freedom locking. As a simple method to circumvent the enriched degree of freedom locking, the tensile strength of the neighboring elements was slightly increased in the numerical examples, and the effectiveness of the method was demonstrated. The proposed method is simple and easy for practicing engineers, and it can be easily applied to the three-dimensional crack propagation analysis.

Introduction

The fracture or cracking of concrete is a significant nonlinear phenomenon that dominates the overall behavior of the structural members. Owing to the development of the finite element method, analysis methods for the cracking of concrete have been available since late 1960s. The discrete crack method, which was introduced by Ngo and Scordelis, ¹ is aimed at predicting the initiation and propagation of dominant cracks by separating the pre-assumed cracks. One clear difficulty with such an approach is that the locations and orientations of the cracks are not known in advance. This drawback could be overcome by remeshing and refinement;^2,3 however, these approaches are extremely complex and computationally expensive. The smeared crack method, which was proposed by Rashid, ⁴ is based on the assumption that the cracked concrete remains a continuum, and cracks are smeared out over the continuum. ⁵ However, such continuum models are mesh-dependent, given that a finer discretization leads to a decrease in the fracture energy. ⁶

Effective numerical methods to overcome the drawbacks of the early solutions to the crack problems were developed, such as mesh-free methods,^7–9 boundary element methods,^10,11 extended finite element methods (XFEMs),^12,13 and generalized finite element methods.^14,15 The mesh-free method can simplify the crack analysis by extending the crack surface, as it does not use a structured mesh, but a scattered set of nodal points in the domain of interest. The boundary element method was used to solve crack problems, as it only requires a boundary and crack surface. It is more efficient, particularly in solving mixed mode crack propagation problems, due to its efficiency in stress concentration modeling and cohesive crack growth. ¹⁶ The generalized finite element method incorporates the analytically known or numerically computed handbook functions to enhance the local and global accuracies of the computed solution. ¹⁷

The XFEM, which was proposed by Belytschko and Black, ¹⁸ uses a partition of unity method ¹⁹ and an enriched function to calculate the discontinuous displacement field, without complex numerical computations. In addition, level set methods are used to track the shape and propagation direction of the cracks, so that the shape of the cracks, independent of the size or shape of the element, can be visualized.^20,21 Furthermore, given that the shape of the crack is expressed by a line or a surface inside the element, a remesh process due to crack propagation is not required. Therefore, the computation time due to crack propagation can be greatly reduced; furthermore, the method is effective in expressing multiple cracks. The XFEM describes the discontinuity of the displacement field using Heaviside functions and branch functions. The most important feature is that these enriched functions are multiplied by the shape functions of finite elements to easily describe the discontinuity. Thereafter, the cohesive crack method^22–24 was developed to improve the description of crack tip characteristics without the calculation of complex branch functions and stress intensity factors. The phantom node method proposed by Song et al. ²⁵ and Rabczuk et al. ²⁶ contributed to simplifying the calculation of the discontinuity displacement field of a fully cracked element. This method simplifies the calculation by adding a phantom node or phantom degree of freedom (DOF) to implement the additional displacement.

The XFEM proposed for solving simple discontinuity problems was used for arbitrary three-dimensional cracks. ²⁷ A study of cracks in various materials such as composite materials, ²⁸ incompressible solids, ²⁹ and metallic materials using a ductile damage model ³⁰ was then carried out. In addition, the application range of cracked plates, including holes, inclusions, minor cracks, ³¹ and multiple cracks of reinforced concrete flexural members, ³² is gradually expanding. Therefore, XFEM was developed by several researchers and is currently used as a means of crack analysis in various fields. For example, Abaqus, which is a popular commercial finite element analysis program, has provided XFEM functionality since 2009. ³³

However, problems with XFEM were observed in cases in which excessive stress is concentrated in the localized region, due to the unusual shape of the model, or the analysis of interaction behavior using two or more materials. In these cases, cracks occur simultaneously in adjacent elements, and the discontinuous behavior due to the cracks is not visible. This problem was reported by Yoo and Kim ³⁴ and was termed enriched DOF locking. When the load is increased excessively in a loading step, or interactions between different materials such as concrete and reinforcements are applied to the model, adjacent elements can simultaneously satisfy the crack initiation criterion. Under these circumstances, the XFEM determines and displays the cracks in all the elements satisfying the failure criterion. However, although the cracks are determined and displayed by the XFEM, the cracks are not actually open, and the stresses of the cracked elements are not released. Figure 1(a) presents an example of the enriched DOF locking phenomenon in the flexural analysis of a reinforced concrete beam. Adjacent cracks developed at the center of beam and the cracks would not propagate upward. Although the beam with the enriched DOF locking shows cracks, it behaves like a linear elastic beam. Figure 1(b) presents the deformed shape of beam without the locking. The crack at the center initiated at the bottom of the beam and propagated upward.

Figure 1.

Example of the enriched DOF locking in the XFEM analysis: (a) reinforced concrete beam with the enriched DOF locking and (b) reinforced concrete beam without the enriched DOF locking.

In this study, the characteristics of enriched DOF locking, where cracks occur simultaneously in two adjacent elements, were investigated. The XFEM formulation is described in section “XFEM.” Based on this, in section “Enriched DOF locking,” the verification of the continuity of the displacement field of the cracked element during the occurrence of the enriched DOF locking using two-node truss elements and four-node quadrilateral elements is described. In section “Numerical examples,” numerical examples are presented to demonstrate the effect of the enriched DOF locking. Methods to circumvent the locking are also presented considering these examples. Finally, the conclusions are presented.

XFEM

The discontinuous displacement field is specified using the enrichment functions, as given in equation (1), in the original XFEM

$$\mathit{u}\left(\mathit{X}\right)=\sum _{i=1}^N{N}_i\left(\mathit{X}\right)\{{\mathit{u}}_i+{\mathit{q}}_{i}H\left(\mathit{X}\right)\}$$

(1)

where $N$ is the number of nodes per element; $\mathit{X}$ denotes material coordinates; $N_i(\mathit{X})\mathrm{s}$ are the shape functions used in the conventional finite element method; ${\mathit{u}}_{\mathit{i}}$ and ${\mathit{q}}_{\mathit{i}}$ denote the standard and additional (enriched) DOFs at the ith node; and $H(\mathit{X})$ is an enriched function (Heaviside function) that represents discontinuous displacement, as shown in equation (2)

$$H\left(\mathit{x}\right)=\{\begin{array}{c}1:x>0\hfill \\ \multicolumn{1}{c}{0:x\le 0}\end{array}$$

(2)

However, the approximation given in equation (1) does not yield $\mathit{u}({\mathit{x}}_{\mathit{i}})={\mathit{u}}_{\mathit{i}}$ , which makes the imposition of essential boundary conditions difficult. Moreover, the interpretation of the results is difficult, as $\mathit{u}({\mathit{x}}_{\mathit{i}})$ is required to be constructed correctly by evaluating all the terms in the approximation. It is therefore desirable to have enrichment terms that vanish at all the nodes, and this is achieved by shifting the approximation. ³⁵ The phantom node method, which was proposed by Song et al., ²⁵ is a shifted version of the original XFEM, to simplify the implementation of the XFEM. It is important to note that the shifted formulation is still able to accurately reproduce the enrichment functions. The phantom node method forms overlapping elements by superimposing phantom nodes on cracked elements. Given that the overlapping element and the overlapped element are independent of each other, each element can be integrated without having to divide the cracked element into subdomains. As a result, it is possible to simplify the calculation required for the generation and propagation of cracks, when compared with the original XFEM.

The XFEM with the phantom node method was used for the basic formulation, to demonstrate the enriched DOF locking. The discontinuous displacement field $\mathit{u}(\mathit{X})$ of an enriched element without the consideration of the singularity of the crack tip is expressed by equation (3)

$$\mathit{u}\left(\mathit{X}\right)=\sum _{i=1}^N{N}_i\left(\mathit{X}\right)\{{\mathit{u}}_i+{\mathit{q}}_i[H\left(f\left(\mathit{X}\right)\right)-H\left(f\left({\mathit{X}}_i\right)\right)]\}$$

(3)

where ${\mathit{X}}_{\mathit{i}}$ denotes the coordinates of the ith node, and $f(\mathit{X})$ represents the location of cracks, as shown in Figure 2.

Figure 2.

Representation of a discontinuity in a truss element for XFEM with the phantom node method.

The truss element with two nodes is first discussed to illustrate the process of applying the phantom node method to XFEM. Equation (4) is obtained by expanding equation (3) for the truss element. In this case, $H=H(f(\mathit{X}))$ , and $u_1$ and $u_2$ are the standard DOFs at nodes 1 and 2, as shown in Figure 2. The additional DOFs $q_1$ and $q_2$ are nodal unknowns, corresponding to phantom nodes 3 and 4, as shown in Figure 2

$$\mathit{u}=u_1{N}_1+u_2{N}_2+q_1{N}_{1}H+q_2{N}_2(H-1)$$

(4)

When a crack occurs in element $A$ , it is divided into two regions, and the phantom elements $A_1$ and $A_2$ are formed by the phantom nodes. The DOF of each phantom element is given by equation (5)

$$\mathrm{element}{A}_1:\{\begin{array}{c}{u}_1^1=u_1\hfill \\ \multicolumn{1}{c}{u_2^1=u_2-q_2}\end{array}\mathrm{element}{A}_2:\{\begin{array}{c}{u}_1^2=u_1+q_1\hfill \\ \multicolumn{1}{c}{u_2^2=u_2}\end{array}$$

(5)

As shown in Figure 2, phantom element $A_1$ consists of real node 1 and phantom node 4, and phantom element $A_2$ consists of real node 2 and phantom node 3. To convert the displacement equation of the original XFEM for the phantom elements, it is necessary to replace it with the nodal displacements given in equation (5). For this, equation (4) is rewritten as follows

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{u}=(u_1+q_1)N_{1}H+u_1{N}_1(1-H)}\\ \multicolumn{1}{c}{+(u_2-q_2)N_2(1-H)+u_2{N}_{2}H}\end{array}$$

(6)

where $N_1=N_{1}H+N_1(1-H),N_2=N_{2}H+N_2(1-H)$ were used. This is summarized according to the DOF of the phantom element of equation (5)

$$\begin{array}{cc}\multicolumn{1}{c}{\mathit{u}=}& u_1^1{N}_1(1-H)+u_2^1{N}_2(1-H)\hfill \\ \multicolumn{1}{c}{}& +u_1^2{N}_{1}H+u_2^2{N}_{2}H={\mathit{u}}^1+{\mathit{u}}^2\hfill \end{array}$$

(7)

The region within the element in the range of $f(\mathit{X})\le 0$ has a displacement field governed by $u_1^1\mathrm{and}{u}_2^1$ of the phantom element $A_1$ , given that $H=0$ . On the contrary, the region in the range of $f(\mathit{X})>0$ is dominated by the DOFs $u_1^2\mathrm{and}{u}_2^2$ of the phantom element $A_2$ . Thus, this element has a discontinuous displacement field at $f(\mathit{X})=0$ . Figure 2 reveals that there is a displacement field different from the displacement field $u_{\mathrm{std}}^A$ of a conventional finite element without cracks. When a crack occurs, the displacement field through the phantom element is formed independently in each range, depending on the position of the crack $(f(\mathit{X})=0)$ . Therefore, the phantom elements independent of each other can be integrated without being divided into subdomains.

This process can be applied to quadrilateral elements with four nodes. Using the same method as before, the following equation is developed for the quadrangle element, by substituting $H_i=H(f({\mathit{X}}_{\mathit{i}}))$ and $N_i=N_{i}H+N_i(1-H)$ into equation (3)

$$\mathit{u}=\sum _{i=1}^4[{\mathit{u}}_i{N}_i(1-H)+{\mathit{u}}_i{N}_{i}H+{\mathit{q}}_i(H-H_i)N_i]$$

(8)

When a crack occurs in an element, a phantom node is added as shown in Figure 3, and a phantom element is formed according to the crack. Thereafter, the DOF of each phantom element is as follows

$${\mathbf{u}}_{\mathrm{i}}^1=\{\begin{array}{cc}{\mathit{u}}_i\hfill & :f\left({\mathit{X}}_i\right)\le 0\hfill \\ \multicolumn{1}{c}{{\mathit{u}}_i-{\mathit{q}}_i}& :f\left({\mathit{X}}_{\mathit{i}}\right)>0\hfill \end{array}{\mathbf{u}}_{\mathrm{i}}^2=\{\begin{array}{cc}{\mathit{u}}_i+{\mathit{q}}_i\hfill & :f\left({\mathit{X}}_i\right)\le 0\hfill \\ \multicolumn{1}{c}{{\mathit{u}}_i}& :f\left({\mathit{X}}_i\right)>0\hfill \end{array}$$

(9)

Figure 3.

Representation of phantom nodes in a quadrilateral element by crack type: (a) completely cracked elements and (b) partially cracked elements.

Equation (7) was then further expanded by duplicating it with the multipliers $H_i^{+}=H(f({\mathit{X}}_i)>0)$ and $H_i^{-}=H(-f({\mathit{X}}_i)\le 0)$ . When $H_i^{+}=1$ , $H_i^{-}=0$ is always satisfied. Therefore, equation (8) is expanded by substituting ${\mathit{u}}_{\mathit{i}}{N}_i={\mathit{u}}_{\mathit{i}}{N}_i{H}_i^{+}+{\mathit{u}}_{\mathit{i}}{N}_i{H}_i^{-}$

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{u}=\sum _{i=1}^4[{\mathit{u}}_i{H}_i^{+}{N}_i(1-H)+{\mathit{u}}_i{H}_i^{-}{N}_i(1-H)+{\mathit{u}}_i{H}_i^{+}{N}_{i}H+{\mathit{u}}_i{H}_i^{-}{N}_{i}H+{\mathit{q}}_i{H}_i^{+}{N}_i(H-1)+{\mathit{q}}_i{H}_i^{-}{N}_{i}H]}\end{array}$$

(10)

The aforementioned equation is then rewritten as

$$\begin{array}{cc}\multicolumn{1}{c}{\mathit{u}}& =\sum _{i=1}^4[{\mathit{u}}_i{H}_i^{+}{N}_{i}H+({\mathit{u}}_i-{\mathit{q}}_i)H_i^{+}{N}_i(1-H)+({\mathit{u}}_i+{\mathit{q}}_i)H_i^{-}{N}_{i}H+{\mathit{u}}_i{H}_i^{-}{N}_i(1-H)]\hfill \\ \multicolumn{1}{c}{}& =\sum _{i\in S_1}{\mathit{u}}_i^1{N}_i(1-H)+\sum _{i\in S_2}{\mathit{u}}_i^2{N}_{i}H\hfill \end{array}$$

(11)

where $S_1$ and $S_2$ denote the index sets of superposed elements 1 and 2, respectively.

In equation (11), the two terms in the second equation represent the displacement fields of phantom elements $A_1$ and $A_2$ , respectively. Like the truss element with two nodes, an independent displacement field is constructed, which corresponds to the range divided by the crack $(f(\mathit{X})=0)$ . In addition, it is possible to integrate without subdividing it into subdomains, which greatly reduces computational requirements.

Enriched DOF locking

Two-node truss elements

Figure 4 presents the relationship between the elements and cracks when two adjacent elements are cracked. When two neighboring elements A and B are cracked, the locations of the two cracks are given by $f(\mathit{X})=0$ and $g(\mathit{X})=0$ , respectively; and according to equation (4), the displacement field of each element is given by equations (12) and (13), respectively

$${\mathit{u}}^A=u_1{N}_1+u_2{N}_2+q_1{N}_{1}H+q_2{N}_2(H-1)$$

(12)

$${\mathit{u}}^B=u_2{N}_1+u_3{N}_2+q{\prime }_2{N}_{1}H+q_3{N}_2(H-1)$$

(13)

Figure 4.

The arrangement of phantom nodes when the enriched DOF locking occurs in truss elements.

The enriched DOFs added to each real node are denoted by $q_i$ , and the enriched DOFs of phantom nodes 5 and $5\prime$ are $q_2$ and $q{\prime }_2$ . In Figure 4, the enriched DOFs $q_2$ and $q{\prime }_2$ are presented as separated, to emphasize that the phantom elements of elements A and B are cracked. According to the phantom node method, ²⁵ only one phantom node is added to each real node when a crack develops in the element connected to the node. Therefore, the phantom node added to real node 2 is one. Moreover, nodes 5 and 5′ are identical, and $q_2$ and $q{\prime }_2$ are equal. If two phantom nodes were added to the real node, then the real node would be a free node when two adjacent elements are cracked, as shown in Figure 4. Using equation (7), equations (12) and (13) can be obtained from equations (14) and (15), respectively

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{u}}^A=}& u^{A_1}+u^{A_2}=u_1^{A_1}{N}_1(1-H)+u_2^{A_1}{N}_2(1-H)\hfill \\ \multicolumn{1}{c}{}& +u_1^{A_2}{N}_{1}H+u_2^{A_2}{N}_{2}H\hfill \end{array}$$

(14)

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{u}}^B=}& u^{B_1}+u^{B_2}=u_1^{B_1}{N}_1(1-H)+u_2^{B_1}{N}_2(1-H)\hfill \\ \multicolumn{1}{c}{}& +u_1^{B_2}{N}_{1}H+u_2^{B_2}{N}_{2}H\hfill \end{array}$$

(15)

The displacement field of phantom elements $A_1$ and $B_2$ is now expressed using equation (5) as equations (16) and (17). Furthermore $H=0$ is used in $A_1$ , and $H=1$ is used in $B_2$ . It should be noted that both ${\mathit{u}}^{A_1}$ and ${\mathit{u}}^{B_2}$ contain $q_2(=q{\prime }_2)$

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{u}}^{A_1}}& =u_1^{A_1}{N}_1+u_2^{A_1}{N}_2=u_1{N}_1+(u_2-q_2)N_2\hfill \\ \multicolumn{1}{c}{}& =u_{\mathrm{std}}^A-q_2{N}_2:f\left(\mathit{X}\right)\le 0\hfill \end{array}$$

(16)

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{u}}^{B_2}}& =u_1^{B_2}{N}_1+u_2^{B_2}{N}_2=(u_2+{q\prime }_2)N_1+u_3{N}_2\hfill \\ \multicolumn{1}{c}{}& =u_{\mathrm{std}}^B+q{\prime }_2{N}_1:g\left(\mathit{X}\right)>0\hfill \end{array}$$

(17)

where $u_{\mathrm{std}}^A=u_1{N}_1+u_2{N}_2$ and $u_{\mathrm{std}}^B=u_2{N}_1+u_3{N}_2$ represent the continuous displacement field of the two-node conventional truss element.

Using $N_1=1-x/L$ and $N_2=x/L$ , the nodal displacements at nodes 5 and 5′ in elements $A_1$ and $B_2$ are given as follows. Note that the displacement field in $A_1$ and $B_2$ can be evaluated over the entire domain, although the integrated domains are defined as $f(\mathit{X})\le 0$ and $g(\mathit{X})>0$ , respectively. Song et al. ²⁵ used this in the computation of the force matrix for a cracked element, by modifying it with the area fraction

$${\mathit{u}}^{{\mathit{A}}_1}\left(L\right)=u_{\mathrm{std}}^A-q_2{N}_2=u_2-q_2$$

(18)

$${\mathit{u}}^{B_2}\left(0\right)=u_{\mathrm{std}}^B+q{\prime }_2{N}_1=u_2+q_2$$

(19)

As mentioned earlier, nodes 5 and 5′ are identical, and $q_2$ and $q{\prime }_2$ are equal. Therefore, the nodal displacements at nodes 5 and 5′, given by equations (18) and (19), are equal, and $q_2=0$ . Substituting this into equations (16) and (17), the additional displacements $-q_2{N}_2$ and $q{\prime }_2{N}_1$ are eliminated. As a result, the phantom elements $A_1$ and $B_2$ have the same displacement field as the conventional finite element, despite the enriched DOF $q_2$ . Moreover, it should be noted that the nodal displacement at node 5 is equal to the nodal displacement at node 2 as follows

$${\mathit{u}}^{{\mathit{A}}_1}\left(L\right)={\mathit{u}}^{{\mathit{A}}_2}\left(L\right)={\mathit{u}}^{B_1}\left(0\right)={\mathit{u}}^{B_2}\left(0\right)=u_2$$

(20)

It can be said that phantom node 5 does not exist and that no cracks are present if two cracks develop on adjacent elements.

Figure 5 illustrates various cracking scenarios on the one-dimensional element mesh. Figure 5(a) presents the mesh before cracking, and Figure 5(b) presents the separated elements $A_1$ and $A_2$ after a crack develops on element $A$ . The displacement field on element $A$ is discontinuous, as intended. Figure 5(c) depicts the enriched DOF locking between two cracks on elements $A$ and $B$ . Real node 2 and phantom node 2′ are locked, and the displacement fields on elements $A$ and $B$ are continuous. Figure 5(d) presents the enriched DOF locking that occurred on nodes 2 and 3. Although three cracks developed, the structure behaved as though it were uncracked. Figure 5(e) presents two unlocked cracks that developed on elements $A$ and $C$ . These two cracks separated elements $A$ and $C$ into $A_1$ and $A_2$ , and $C_1$ and $C_2$ , respectively. The enriched DOF locking did not occur because the cracks did not develop on adjacent elements.

Figure 5.

Various cracking scenarios: (a) before cracking, (b) an unlocked crack on element A, (c) two locked cracks on elements A and B, (d) three locked cracks on elements A, B, and C, and (e) two unlocked cracks on elements A and C.

Four-node quadrilateral elements

To represent and interpret the cracks of the actual material, four-node quadrilateral elements were used rather than two-node truss elements. As shown in Figure 6, it was assumed that two cracks developed simultaneously on the adjacent elements $A$ and $B$ , and the position of each crack was given by $f(\mathit{X})=0$ and $g(\mathit{X})=0$ , respectively. If two cracks are assumed to develop at the bottom of the elements and the crack tips just reach the top edges of the elements, as shown in Figure 6, the opening displacements at the crack tips are assumed to be zero. Therefore, phantom nodes were not added to nodes 4, 5, and 6. Only phantom nodes 7, 8, and 9 were added to nodes 1, 2, and 3, respectively, where discontinuous displacements occur due to the cracks. It could therefore be shown that the nodal displacement of nodes 8, which is the same node as 8′, is equal to the nodal displacement of node 2, using the phantom elements $A_1$ and $B_2$ .

Figure 6.

The arrangement of phantom nodes when the enriched DOF locking occurs in quadrilateral elements.

Equation (10) or (11) was also used to demonstrate how the enriched DOF locking occurs in the quadrilateral element. Similar to the truss element, the aforementioned equation is developed for phantom elements $A_1$ and $B_2$ , which share a common phantom node. The displacement at phantom node 8 in the phantom element $A_1$ $(N_2=1)$ in the range of $f(\mathit{X})\le 0$ is given as follows

$${\mathit{u}}^{A_1}=\sum _{i\in S_1}{\mathit{u}}_i^1{N}_i(1-H)=u_2-q_2:f\left(\mathit{X}\right)\le 0$$

(21)

Likewise, the displacement at phantom node 8′ in the phantom element $B_2$ $(N_1=1)$ in the range of $g(\mathit{X})>0$ is given as follows

$${\mathit{u}}^{B_2}=\sum _{i\in S_2}{\mathit{u}}_i^2{N}_{i}H=u_2+q_2:g\left(\mathit{X}\right)>0$$

(22)

As can be seen in equations (18) and (19), $q_2=0$ because the displacements at phantom nodes 8 and 8′ calculated in each phantom element are equal.

Numerical examples

Simply tensioned bar

It was verified that the enriched DOF locking develops when two cracks occur simultaneously on the adjacent elements. In this section, the effect of the enriched DOF locking is demonstrated using Abaqus, ³⁶ which is a commercial finite element analysis program implementing the XFEM feature. First, three 10 mm × 10 mm elements were arranged adjacent to each other, and a pinned boundary condition was applied at nodes 1 and 5, as shown in Figure 7(a). The material is assumed to be linear elastic with the Young’s modulus $E$  = 30 GPa and Poisson’s ratio $v=0.2$ . The damage modeling that is required to simulate the cracking of the XFEM analysis consists of two ingredients: a crack initiation criterion and a crack evolution law. The maximum principal stress criterion with tensile strength of 3 MPa and the linear traction-separation with fracture energy $G_F$  = 0.1 N/mm were used for the crack initiation criterion and crack evolution law. The four-node bilinear plane stress element with reduced integration CPS4R was used. The thickness of the elements was 10 mm. The external load was applied at nodes 4 and 8 using the displacement control by 0.1 mm. Using simple statics, the cracking was expected to develop at a displacement of $0.003\mathrm{mm}$ .

Figure 7.

Analysis of a simply tensioned bar with the enriched DOF locking: (a) the undeformed shape, and boundary and loading conditions, (b) the deformed shape where the enriched DOF locking occurred, and (c) the principal stress–displacement graphs.

The deformed shape of the case with the enriched DOF locking is presented in Figure 7(b). Although the deformed shape had full penetrated cracks in all elements, the analysis continued like a linear elastic analysis. Figure 7(c) presents the principal stress at the center of each element, according to the magnitude of the applied displacement. Even after reaching a tensile strength of 3 MPa, the stress of elements 1 and 2 increases. The reason for this is that the discontinuity of the displacement field is not correctly calculated due to the locking of the phantom node. The calculated stress of element 3 gradually decreases because it represents the stress of one of the overlapped elements with a free edge, as shown in Figure 5(d). The stresses shown in Figure 7(c) correspond to the stresses in elements A1, B1, and C1 shown in Figure 5(d).

To solve the enriched DOF locking problem, it is necessary to ensure that the stresses or strains of adjacent elements do not simultaneously satisfy the crack initiation criterion. To achieve this, it is proposed that the tensile strength of the neighboring elements should be slightly increased. The tensile strength of the dark-shaded elements shown in Figure 8(a) is increased from 3 to 3.1 MPa. The tensile strength of the light-shaded element is retained as 3 MPa. Considering the heterogeneity of concrete, the suggested method to evade the enriched DOF locking is acceptable. Figure 8 presents the result of the case without the locking, which was solved using the proposed method. Due to the slightly increased tensile strength of elements 1 and 3 to circumvent the enriched DOF locking, only one crack developed in element 2, as shown in Figure 8(b). The locking did not occur at the phantom node between the adjacent elements. Thus, the discontinuous displacement in element 2 was calculated correctly. Moreover, as the applied displacement increased, the principal stress in element 2 decreased and reached zero after reaching the tensile strength. The stresses in elements 1 and 3 were also released since the entire structure was completely separated after cracking as shown in Figure 8(c).

Figure 8.

The results of a simply tensioned bar where the enriched DOF locking is circumvented: (a) assigning slightly higher tensile strength to the neighboring elements, (b) the deformed shape without the enriched DOF locking, and (c) the principal stress–displacement graphs.

Reinforced concrete beam

Next, the effect of the enriched DOF locking on the actual reinforced concrete structure was investigated. The shape of this specimen is presented in Figure 9. The material properties of the concrete were as follows: Young’s modulus = 33.1 GPa, Poisson’s ratio = 0.18, tensile strength = 3.67 MPa, and fracture energy $(G_F)$  = 0.112 N/mm. The material properties of the reinforcing bars are assumed to be: Young’s modulus = 200 GPa, Poisson’s ratio = 0.3, yield strength $(f_y)$  = 484 MPa, and ultimate strength $(f_k)$  = 630 MPa at the plastic strain of 0.1. The concrete was assumed linear elastic, except for the tensile fracture, which needs to be specified for XFEM. Tensile principal stress and linear damage evolution with fracture energy were used for the crack initiation criterion and the crack evolution law, respectively. The mechanical behavior of the reinforcing bars was modeled with isotropic hardening plasticity. The concrete and the reinforcing bars were modeled using 700 quadrilateral plane stress elements and 70 truss elements, respectively. Embedded region, which is one of the constraint methods in Abaqus, was applied to the models to implement perfect bonding between the concrete and steel bar. The vertical displacements were applied at one-third and two-third points using the displacement control by 2.0 mm. The center one-third portion is prone to the enriched DOF locking because of the uniform bending moment in the region.

Figure 9.

Dimensions of reinforced concrete beam (mm).

When the displacement was applied to the reinforced concrete beam, several cracks simultaneously developed at the center one-third portion of the beam because the concrete elements on the bottom side in this region were subjected to uniform tensile stress. As shown in Figure 10(a), the locking occurred because the cracks developed simultaneously in the adjacent elements located in this region. Figure 10(d) shows that the tensile stress at the bottom center element linearly increases as the applied displacement increases when the enriched DOF locking occurs. To evade the locking, the proposed method of assigning slightly higher tensile strength to the neighboring elements was also applied for this case as in the simply tensioned bar. By assigning slightly different material properties such as tensile strength, it is possible to ensure that the stresses or strains of adjacent elements do not simultaneously satisfy the crack initiation criterion. The location for the element with higher tensile strength can be easily determined by observing the crack pattern of the beam with locking as shown in Figure 10(a).

Figure 10.

The results of the reinforced concrete beam: (a) crack pattern of the reinforced concrete (RC) beam with locking, (b) analysis mesh to evade the enriched DOF locking, (c) crack pattern of the RC beam without locking, and (d) principal stress at the bottom center element versus applied vertical displacement.

Figure 10(b) shows the analysis mesh in which the tensile strength of the dark shaded elements was modified from 3.67 to 3.8 MPa. Only two elements—in the bottom layer at the center one-third portion—were modified because the cracks will propagate without locking once the initial locking is evaded. A trial was also conducted by modifying one element, but the enriched DOF locking still occurred. As a result, cracks initiated at the unmodified bottom elements and propagated upward as shown in Figure 10(c). The tensile principal stress at the bottom center element of the case without the locking is also shown in Figure 10(d). When the locking occurred, the stress in the locked element linearly increased beyond the tensile strength. On the contrary, the model analyzed without the locking revealed that the stresses were released after the tensile strength was reached. However, the stress did not drop to zero because the cracked concrete was completely constrained to the reinforcement.

Conclusion

In this study, the enriched DOF locking that can occur in XFEM analysis was investigated. It was verified that the inactivated discontinuous displacement field was due to the enriched DOF locking, which occurred when two adjacent elements shared the enriched DOF. Although the phantom node method was used in the verification, the enriched DOF locking could occur in all versions of the XFEM, unless another set of enriched DOFs was provided. The enriched DOF locking occurs when two adjacent elements are cracked simultaneously according to the crack initiation criterion. The crack initiation criterion determines the given element cracked if the analysis result of the element, such as the principal stress, exceeds the user-defined value. The analysis results can exceed the criterion at a load increment when the analysis model has a relatively uniform stress distribution over a certain range, due to the geometric characteristics or the stress transfer through the interconnected elements such as the reinforcing bars. Numerical examples of a simply tensioned bar and a reinforced concrete beam were presented to demonstrate the two structures prone to enriched DOF locking. Moreover, a simple method of increasing the tensile strength of the neighboring elements to circumvent the enriched DOF locking was proposed. By modifying the tensile strength of the neighboring elements, the simultaneous cracking of two adjacent elements can be evaded. Applying these methods to prevent enriched DOF locking enabled the achievement of appropriate cracking analyses, which revealed the discontinuous displacement fields and the release of the stress in the cracked elements after cracking. The number of modified elements can be minimized because the cracks will propagate properly once the initial locking is evaded. Considering the heterogeneity and variability of concrete, the proposed method of modifying the tensile strength by less than 5% is acceptable. The proposed method is simple and easy for practicing engineers since it does not involve computer programming, such as user subroutines, or complicated combination of model parameters. This method can be easily applied to the three-dimensional crack propagation analysis.