A family of non-conforming crack front elements of quadrilateral and triangular types for 3D crack problems using the boundary element method

Abstract

This study focuses on establishing non- conforming crack front elements of quadrilateral and triangular types for 3D crack problems when the dual boundary element method is applied. The asymptotic behavior of the physical variables in the area near the crack front is fully considered in the construction of the shape function. In the developed quadrilateral and triangular crack front elements, the asymptotic term, which captures the asymptotic behavior of the physical variable, is multiplied directly by the conventional Lagrange shape function to form a new crack front shape function. Several benchmark numerical examples that consider penny-shaped cracks and straight-edge crack problems are presented to illustrate the validity and efficiency of the developed crack front elements.

Key words： Taylor expansion; crack front elements; stress intensity factors; dual boundary element method

Introduction

Experimental [¹–⁴] and numerical [⁵–⁹] are two of the most popular methods to solve engineering problems, such as solid mechanical, heat transfer, acoustic, and electro-magnetic problems. In structural mechanical problems with cracks or V-notch interfaces, physical quantities vary singularly in the areas near the crack front and interface [⁵–⁹]. Accurate approximation of such singular variations is important but difficult in numerical methods for these problems. Many additional methods have been developed to circumvent these difficulties [⁵–⁹]. The finite element method (FEM) [¹⁰–¹²] and the boundary element method (BEM) [¹³–¹⁵] are two of the most extensive methods for structural mechanical problems. Special elements with specially defined shape functions are typically utilized to approximate the variations of physical quantities in both methods and eventually achieve high accuracy and efficiency. Many specially defined crack elements have been developed in FEM and BEM implementations to compute the stress intensity factors (SIFs) along the crack front accurately and efficiently. The most efficient and valuable crack front elements are quarter-point triangular and quadrilateral elements, which can be integrated into standard FEM [¹⁶,¹⁷] and BEM [¹⁸–²¹] codes. In BEM application in crack problems, the dual boundary integral equation method, in which the traction boundary integral equation is involved, is the most widely used scheme [⁵,¹⁹–²⁴].

Non-conforming elements, inside which interpolation nodes are arranged, are frequently used due to the existence of hyper-singular integrals. The calculation of hyper-singular integrals depends largely on the transformation of traditional conforming quarter-point elements into nonconforming ones. However, quarter-point triangular and quadrilateral elements are usually regarded as conforming elements, which are employed in the multidomain BEM. The shape functions of non-conforming crack front elements are obtained by solving linear equations. The derivation process is tedious, and the special shape functions are not closely related with traditional shape functions [²⁰,²⁵,²⁶]. This study proposes a family of non-conforming crack front elements of quadrilateral and triangular types for 3D crack problems.

In the deduction of the presented crack front elements, the shape functions are constructed to capture the asymptotic behavior of displacement fields in parametric space. The function of the analytical distribution of displacement is expressed by Taylor series up to the second order to capture the asymptotic behavior. Variable transformation is performed to generate a square-root term that captures the variation of the physical quantities. The shape functions are constructed in consideration of the square-root term. Then, the conforming element is transformed into a nonconforming one by a multiplication of the transforming matrix. Several benchmark examples are presented to illustrate the accuracy of the new crack front elements. The advantages of the non-conforming crack front elements are as follows. The special shape functions are obtained using the traditional shape functions and displacement approximate terms in the intrinsic coordinate system. Unlike the traditional construction method for crack front elements [²⁰,²⁵,²⁶], the method proposed in this study can be easily extended to the construction of special shape functions for V-notch interfaces, whose displacements may possess a property of order r^{\alpha } (\mathrm{1}>\alpha >\mathrm{0}).

The rest of this paper is outlined as follows. In Section 2, crack front elements within a quadrilateral mesh are proposed. Then, crack front elements within a triangular mesh are obtained in the same manner in Section 3. The process of obtaining SIFs is described in Section 4, and several numerical examples are presented in Section 5.

Crack front elements within a quadrilateral mesh

We provide a detailed derivation of a three-node, non-conforming crack tip element. The method for the three-node, non-conforming crack tip element is then extended to the 3D boundary surface element.

Construction of a three-node, non-conforming crack tip element

As shown in Fig. 1, the coordinates of the distance vector from arbitrary point P(x, y) in the element to the crack tip can be expressed as

(1a)x-x_{\mathrm{0}}=N_{\mathrm{0}}{x}_{\mathrm{0}}+N_{\mathrm{1}}{x}_{\mathrm{1}}+N_{\mathrm{2}}{x}_{\mathrm{2}}-x_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\xi +\mathrm{1})\left[(\xi -\mathrm{2})x_{\mathrm{0}}+\mathrm{2}(\mathrm{1}-\xi )x_{\mathrm{2}}+\xi x_{\mathrm{1}}\right],

(1b)y-y_{\mathrm{0}}=N_{\mathrm{0}}{y}_{\mathrm{0}}+N_{\mathrm{1}}{y}_{\mathrm{1}}+N_{\mathrm{2}}{y}_{\mathrm{2}}-y_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\xi +\mathrm{1})\left[(\xi -\mathrm{2})y_{\mathrm{0}}+\mathrm{2}(\mathrm{1}-\xi )y_{\mathrm{2}}+\xi y_{\mathrm{1}}\right],

where \xi is the parametric coordinate of the point \left(x,y\right), \left(x,y\right) is the coordinates of any point in the element, \left(x_{\mathrm{0}},y_{\mathrm{0}}\right), \left(x_{\mathrm{1}},y_{\mathrm{1}}\right), and \left(x_{\mathrm{2}},y_{\mathrm{2}}\right) are the coordinates of the three nodes in Fig. 1 (Points 0, 1 and 2), respectively, N_{\mathrm{0}}\left(\xi \right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\xi (\xi -\mathrm{1}), N_{\mathrm{1}}\left(\xi \right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\xi (\mathrm{1}+\xi ) and N_{\mathrm{2}}\left(\xi \right)=(\mathrm{1}+\xi )(\mathrm{1}-\xi ) are shape functions that are defined on these nodes, and (\xi +\mathrm{1})/\mathrm{2} represents the distance property in the parametric space, which also reflects the distance in the Cartesian space. The source points in the parametric space are -\lambda, 0, and \lambda.

Fig.1 Three-node, non-conforming quadratic element

Full size|PPT slide

The displacement \mathit{u} and traction fields \mathit{t} near the crack tip behave as square root singularity or inverse square root singularity, which can be written as [²⁷]

(2)\mathit{u}={\mathit{a}}_{\mathrm{0}}+{\mathit{a}}_{\mathrm{1}}\sqrt{r}+{\mathit{a}}_{\mathrm{2}}r+\cdot \cdot \cdot ,

(3)\mathit{t}={\mathit{b}}_{\mathrm{0}}/\sqrt{r}+{\mathit{b}}_{\mathrm{1}}+{\mathit{b}}_{\mathrm{2}}r+\cdot \cdot \cdot .

where {\mathit{a}}_i and {\mathit{b}}_i (i=0, 1, …) represent the coefficient vectors, respectively.

For the three-node, non-conforming quadratic element, the shape functions that are defined based on these collocation nodes can be written as

N_{\mathrm{coll}}^{\mathrm{0}}\left(\xi \right)={\displaystyle \frac{\xi (\xi -\lambda )}{\mathrm{2}{\lambda }^{\mathrm{2}}}},\hspace{0.17em}{N}_{\mathrm{coll}}^{\mathrm{1}}\left(\xi \right)={\displaystyle \frac{(\xi +\lambda )\xi }{\mathrm{2}{\lambda }^{\mathrm{2}}}},

(4)N_{\mathrm{coll}}^{\mathrm{2}}\left(\xi \right)=-{\displaystyle \frac{(\xi +\lambda )(\xi -\lambda )}{{\lambda }^{\mathrm{2}}}}.

where \lambda is the offset parameter, and N_{\mathrm{coll}}^{\mathrm{0}}\left(\xi \right), N_{\mathrm{coll}}^{\mathrm{1}}\left(\xi \right), and N_{\mathrm{coll}}^{\mathrm{2}}\left(\xi \right) are the non-conforming shape functions of the three collocation nodes in Fig. 1 (Points 0, 1 and 2). Using a two-order Taylor expansion of u near the crack front, we can obtain the following:

(5)\mathit{u}={\mathit{u}}_{\mathrm{0}}^{\mathrm{0}}+{\mathit{u}}_{\mathrm{0}}^1(\xi +\mathrm{1})+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}!}}{\mathit{u}}_{\mathrm{0}}^2{(\xi +\mathrm{1})}^{\mathrm{2}}+\cdot \cdot \cdot ,

where

{\mathit{u}}_0^0={{\displaystyle \sum _{m=\mathrm{0}}^{\mathrm{2}}\frac{{\mathrm{d}}^{\left(\mathrm{0}\right)}{N}_{\mathrm{coll}}^m\left(\xi \right)}{\mathrm{d}{\xi }^{\left(\mathrm{0}\right)}}}}_{\xi =-\mathrm{1}}{\mathit{u}}_m

=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_0+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right){\mathit{u}}_1+\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_2,

{\mathit{u}}_0^1={{\displaystyle \sum _{m=\mathrm{0}}^{\mathrm{2}}\frac{\mathrm{d}{N}_{\mathrm{coll}}^m\left(\xi \right)}{\mathrm{d}\xi }}}_{\xi =-\mathrm{1}}{\mathit{u}}_m

=\left(-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right){\mathit{u}}_0+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_1-{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_2,

{\mathit{u}}_0^2={{\displaystyle \sum _{m=\mathrm{0}}^{\mathrm{2}}\frac{{\mathrm{d}}^{\mathrm{2}}{N}_{\mathrm{coll}}^m\left(\xi \right)}{\mathrm{d}{\xi }^{\mathrm{2}}}}}_{\xi =-\mathrm{1}}{\mathit{u}}_m={\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_0+{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_1+{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_2,

and {\mathit{u}}_0, {\mathit{u}}_1, and {\mathit{u}}_2 are the displacements at three collocation nodes. {\mathit{u}}_0^0, {\mathit{u}}_0^1, and {\mathit{u}}_0^2 represent the zero- first- and second-order derivatives at Node 0.

Compared with Eq. (2), Eqs. (4) and (5) show that the shape functions of the non-conforming quadratic element cannot reflect the square-root asymptotic behavior of the displacements and tractions near the crack front. To introduce the asymptotic behavior of the displacements, a substitution of (\xi +\mathrm{1})/\mathrm{2} by \sqrt{(\xi +\mathrm{1})/\mathrm{2}} is introduced in the Taylor expansion. Thus, the Taylor expansion in Eq. (5) can be transformed into

(6)\mathit{u}={\mathit{u}}_0^0+{\mathit{u}}_0^1\sqrt{\mathrm{2}(\xi +\mathrm{1})}+{\mathit{u}}_0^2(\xi +\mathrm{1})+\cdot \cdot \cdot .

In Eq. (6), {\mathit{u}}_0^0, {\mathit{u}}_0^1, and {\mathit{u}}_0^2 are unchanged. After substituting expressions for {\mathit{u}}_0^0, {\mathit{u}}_0^1, and {\mathit{u}}_0^2 into Eq. (6), the following equation can be obtained:

(7)\begin{array}{l}\mathit{u}=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_0+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right){\mathit{u}}_1+\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_2\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\sqrt{\mathrm{2}}\left(-\left({\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right){\mathit{u}}_0+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_1-{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_2\right)\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\sqrt{\xi +\mathrm{1}}+\left({\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_0+{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_1+{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_2\right)\left(\xi +\mathrm{1}\right)\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\left[\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}\right)-\sqrt{\mathrm{2}}\left({\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right)\sqrt{\xi +\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\left(\xi +\mathrm{1}\right)\right]{\mathit{u}}_0\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left[\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right)+\sqrt{\mathrm{2}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right)\sqrt{\xi +\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\left(\xi +\mathrm{1}\right)\right]{\mathit{u}}_1\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left[\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right)-{\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{2}}}{{\lambda }^{\mathrm{2}}}}\sqrt{\xi +\mathrm{1}}+{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}\left(\xi +\mathrm{1}\right)\right]{\mathit{u}}_2.\hfill \end{array}

New shape functions for collocation points can be defined as

(8)\begin{array}{l}\mathit{u}=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_0+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right){\mathit{u}}_1+\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_2\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\sqrt{\mathrm{2}}\left(-\left({\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right){\mathit{u}}_0+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right){\mathit{u}}_1-{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_2\right)\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \sqrt{\xi +\mathrm{1}}+\left({\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_0+{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_1+{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}{\mathit{u}}_2\right)\left(\xi +\mathrm{1}\right)\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\left[\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}\right)-\sqrt{\mathrm{2}}\left({\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right)\sqrt{\xi +\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\left(\xi +\mathrm{1}\right)\right]{\mathit{u}}_0\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left[\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\lambda }^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}\right)+\sqrt{\mathrm{2}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\lambda }}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right)\sqrt{\xi +\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\left(\xi +\mathrm{1}\right)\right]{\mathit{u}}_1\hfill \\ \hfill \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\left[\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{{\lambda }^{\mathrm{2}}}}\right)-{\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{2}}}{{\lambda }^{\mathrm{2}}}}\sqrt{\xi +\mathrm{1}}+{\displaystyle \frac{\mathrm{2}}{{\lambda }^{\mathrm{2}}}}\left(\xi +\mathrm{1}\right)\right]{\mathit{u}}_2.\hfill \end{array}

where N_{\mathrm{scoll}}^{\mathrm{0}}\left(\xi \right), N_{\mathrm{scoll}}^{\mathrm{1}}\left(\xi \right), and N_{\mathrm{scoll}}^{\mathrm{2}}\left(\xi \right) are the new shape functions at collocation points for crack tip elements.

A three-node non-conforming crack tip element can be successfully constructed with the shape functions of Eq. (8). In the next section, we extend this crack tip element to a nine-node quadrilateral crack front element.

Construction of a nine-node non-conforming crack front element

A nine-node non-conforming quadratic element is illustrated in Fig. 2. The relations between the three-node and nine-node non-conforming quadratic elements can be established as follows. In the nine-node non-conforming quadratic element, interpolation nodes do not coincide with geometric nodes, a condition that is similar to that in the three-node non-conforming quadratic element. The geometric shape functions in the nine-node non-conforming quadratic element are constructed through the products of shape functions along two different directions:

(9)\begin{array}{l}{N}_{\mathrm{0}}(\xi ,\eta )=N_{\mathrm{1}}\left(\xi \right)N_{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{1}}(\xi ,\eta )=N_{\mathrm{0}}\left(\xi \right)N_{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{2}}(\xi ,\eta )=N_{\mathrm{0}}\left(\xi \right)N_{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{3}}(\xi ,\eta )=N_{\mathrm{1}}\left(\xi \right)N_{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{4}}(\xi ,\eta )=N_{\mathrm{2}}\left(\xi \right)N_{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{5}}(\xi ,\eta )=N_{\mathrm{0}}\left(\xi \right)N_{\mathrm{2}}\left(\eta \right),\hfill \\ N_{\mathrm{6}}(\xi ,\eta )=N_{\mathrm{2}}\left(\xi \right)N_{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{7}}(\xi ,\eta )=N_{\mathrm{1}}\left(\xi \right)N_{\mathrm{2}}\left(\eta \right),\hfill \\ N_{\mathrm{8}}(\xi ,\eta )=N_{\mathrm{2}}\left(\xi \right)N_{\mathrm{2}}\left(\eta \right),\hfill \end{array}

where N_i(\xi ,\eta ) (i=\mathrm{0},\hspace{0.17em}\mathrm{1},\hspace{0.17em}\mathrm{...},\hspace{0.17em}\mathrm{8}) are the shape functions that are defined on the eight geometric nodes (from 0 to 8), N_{\mathrm{0}}\left(\xi \right), N_{\mathrm{1}}\left(\xi \right) and N_{\mathrm{2}}\left(\xi \right) have been defined before, and \xi and \eta are the parametric coordinates.

Fig.2 Nine-node non-conforming quadratic element

Full size|PPT slide

Similarly, the non-conforming shape functions of the nine-node non-conforming quadratic element can be obtained through a similar method, which can be expressed as follows:

(10)\begin{array}{l}{N}_{\mathrm{coll}}^{\mathrm{0}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{1}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{1}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{0}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{2}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{0}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{3}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{1}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{4}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{2}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{5}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{0}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{2}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{6}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{2}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{7}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{1}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{2}}\left(\eta \right),\hfill \\ N_{\mathrm{coll}}^{\mathrm{8}}(\xi ,\eta )=N_{\mathrm{coll}}^{\mathrm{2}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{2}}\left(\eta \right),\hfill \end{array}

where N_{\mathrm{coll}}^i(\xi ,\eta ) (i=\mathrm{0},\hspace{0.17em}\mathrm{1},\hspace{0.17em}\mathrm{...},\hspace{0.17em}\mathrm{8}) are the non-conforming shape functions that are defined on the eight collocation nodes (from 0 to 8).

In Fig. 2, given that \xi =-\mathrm{1} is the crack front, the special shape functions can be expressed as

\begin{array}{l}{N}_{\mathrm{Nscoll}}^{\mathrm{0}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{1}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{1}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{0}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{2}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{0}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{3}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{1}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{0}}\left(\eta \right),\hfill \end{array}

(11)\begin{array}{l}{N}_{\mathrm{Nscoll}}^{\mathrm{4}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{2}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{1}}\left(\eta \right),\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{5}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{0}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{2}}\left(\eta \right),\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{6}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{2}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{0}}\left(\eta \right),\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{7}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{1}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{2}}\left(\eta \right),\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{8}}(\xi ,\eta )=N_{\mathrm{scoll}}^{\mathrm{2}}\left(\xi \right)N_{\mathrm{coll}}^{\mathrm{2}}\left(\eta \right),\hfill \end{array}

where N_{\mathrm{Nscoll}}^i(\xi ,\eta ) (i=\mathrm{0},\hspace{0.17em}\mathrm{1},\hspace{0.17em}\mathrm{...},\hspace{0.17em}\mathrm{8}) are the new shape functions on the eight collocation nodes for crack front elements (from 0 to 8).

The displacement distribution near the crack front can be accurately approximated using Eq. (11). The results are demonstrated in the following numerical examples.

Crack front elements within a triangular mesh

In this section, we construct a triangular crack front element based on a six-node conforming quadratic element in the first step. The six-node conforming crack front element is then converted into a six-node non-conforming crack front element in the second step. The six-node non-conforming quadratic element is shown in Fig. 3, where k is the position parameter. From the geometric meaning of the area coordinate, (\mathrm{1}-\xi -\eta ) represents the distance in the parametric space, which also reflects the distance in the Cartesian space. The geometric shape functions of the element are

(12)\begin{array}{l}{N}_{\mathrm{0}}=\xi (\mathrm{2}\xi -\mathrm{1}),\hspace{0.17em}\hspace{0.17em}{N}_{\mathrm{1}}=\eta (\mathrm{2}\eta -\mathrm{1}),\hfill \\ N_{\mathrm{2}}=(\mathrm{1}-\xi -\eta )\left[\mathrm{2}(\mathrm{1}-\xi -\eta )-\mathrm{1}\right],\hspace{0.17em}\hspace{0.17em}{N}_{\mathrm{3}}=\mathrm{4}\xi \eta ,\hfill \\ N_{\mathrm{4}}=\mathrm{4}\eta (\mathrm{1}-\xi -\eta ),\hspace{0.17em}\hspace{0.17em}{N}_{\mathrm{5}}=\mathrm{4}\xi (\mathrm{1}-\xi -\eta ).\hfill \end{array}

Fig.3 Six-node non-conforming quadratic element

Full size|PPT slide

ƒTo perform Taylor expansion of the displacement fields near the neighborhood of the crack front (\xi +\eta =\mathrm{1}), the following coordinate transformation is applied:

(13){\xi }_{\mathrm{1}}=\mathrm{1}-\xi -\eta ,\hspace{0.17em}\hspace{0.17em}{\eta }_{\mathrm{1}}=\xi -\eta .

Then, Taylor expansion of the displacement fields in the neighborhood of \xi +\eta =\mathrm{1} is performed similar to that in the neighborhood of {\xi }_{\mathrm{1}}=\mathrm{0}. We can derive the following expression:

(14)\mathit{u}={\mathit{u}}_0^0+{\mathit{u}}_0^1({\xi }_{\mathrm{1}}-\mathrm{0})+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}!}}{\mathit{u}}_0^2{({\xi }_{\mathrm{1}}-\mathrm{0})}^{\mathrm{2}}+\cdot \cdot \cdot .

where {\mathit{u}}_{\mathrm{0}}^{\mathrm{0}}={{\displaystyle \sum _{m=\mathrm{0}}^{\mathrm{5}}\frac{{\partial }^{\left(\mathrm{0}\right)}{N}_m}{\partial {\xi }_{\mathrm{1}}^{\left(\mathrm{0}\right)}}}}_{{\xi }_{\mathrm{1}}=\mathrm{0}}{\mathit{u}}_m, {\mathit{u}}_0^1={{\displaystyle \sum _{m=\mathrm{0}}^{\mathrm{5}}\frac{\partial N_m}{\partial {\xi }_{\mathrm{1}}}}}_{{\xi }_{\mathrm{1}}=\mathrm{0}}{\mathit{u}}_m, and {\mathit{u}}_0^2={{\displaystyle \sum _{m=\mathrm{0}}^{\mathrm{5}}\frac{{\partial }^{\mathrm{2}}{N}_m}{\partial {\xi }_{\mathrm{1}}^{\mathrm{2}}}}}_{{\xi }_{\mathrm{1}}=\mathrm{0}}{\mathit{u}}_m. {\mathit{u}}_m\left(m=\mathrm{0},\hspace{0.17em}\mathrm{1},\hspace{0.17em}\mathrm{...},\hspace{0.17em}\mathrm{5}\right) denotes the displacements at the six interpolation nodes in the element. From Eq. (13), we obtain

(15)\xi ={\displaystyle \frac{\mathrm{1}-{\xi }_{\mathrm{1}}+{\eta }_{\mathrm{1}}}{\mathrm{2}}},\hspace{0.17em}\eta ={\displaystyle \frac{\mathrm{1}-{\xi }_{\mathrm{1}}-{\eta }_{\mathrm{1}}}{\mathrm{2}}}.

Considering Eq. (15), the following relations are obtained:

(16a){\displaystyle \frac{\partial N_m\left(\xi \right)}{\partial {\xi }_{\mathrm{1}}}}={\displaystyle \frac{\partial N_m\left(\xi \right)}{\partial \xi }}{\displaystyle \frac{\partial \xi }{\partial {\xi }_{\mathrm{1}}}}+{\displaystyle \frac{\partial N_m\left(\xi \right)}{\partial \eta }}{\displaystyle \frac{\partial \eta }{\partial {\xi }_{\mathrm{1}}}},

(16b){\displaystyle \frac{{\partial }^{\mathrm{2}}N}{\partial {\xi }_{\mathrm{1}}^{\mathrm{2}}}}={\displaystyle \frac{{\partial }^{\mathrm{2}}N}{\partial {\xi }^{\mathrm{2}}}}{\displaystyle \frac{\partial \xi }{\partial {\xi }_{\mathrm{1}}}}{\displaystyle \frac{\partial \xi }{\partial {\xi }_{\mathrm{1}}}}+\mathrm{2}{\displaystyle \frac{{\partial }^{\mathrm{2}}N}{\partial \xi \partial \eta }}{\displaystyle \frac{\partial \xi }{\partial {\xi }_{\mathrm{1}}}}{\displaystyle \frac{\partial \eta }{\partial {\xi }_{\mathrm{1}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}N}{\partial {\eta }^{\mathrm{2}}}}{\displaystyle \frac{\partial \eta }{\partial {\xi }_{\mathrm{1}}}}{\displaystyle \frac{\partial \eta }{\partial {\xi }_{\mathrm{1}}}}.

To introduce the asymptotic behavior of the displacements, we correct the Taylor expansion expression by replacing {\xi }_{\mathrm{1}} with \sqrt{{\xi }_{\mathrm{1}}}. The Taylor expansion in Eq. (14) can be transformed into

(17)\mathit{u}={\mathit{u}}_0^0+{\mathit{u}}_0^1\sqrt{{\xi }_{\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}!}}{\mathit{u}}_0^2{\xi }_{\mathrm{1}}+\cdot \cdot \cdot ,

where

\{\begin{array}{l}{\mathit{u}}_0^0=\xi (\mathrm{2}\xi -\mathrm{1}){\mathit{u}}_{\mathrm{0}}+(\mathrm{1}-\xi )(\mathrm{1}-\mathrm{2}\xi ){\mathit{u}}_{\mathrm{1}}+\mathrm{4}\xi (\mathrm{1}-\xi ){\mathit{u}}_{\mathrm{3}},\hfill \\ {\mathit{u}}_0^1=-\mathrm{0.5}[(\mathrm{4}\xi -\mathrm{1}){\mathit{u}}_{\mathrm{0}}+(\mathrm{3}-\mathrm{4}\xi ){\mathit{u}}_{\mathrm{1}}+\mathrm{2}{\mathit{u}}_{\mathrm{2}}+\mathrm{4}{\mathit{u}}_{\mathrm{3}}\hfill \\ \hspace{0.17em}-\mathrm{8}(\mathrm{1}-\xi ){\mathit{u}}_{\mathrm{4}}-\mathrm{8}\xi {\mathit{u}}_{\mathrm{5}}],\hfill \\ {\mathit{u}}_0^2={\mathit{u}}_{\mathrm{0}}+{\mathit{u}}_{\mathrm{1}}+\mathrm{4}{\mathit{u}}_{\mathrm{2}}+\mathrm{2}{\mathit{u}}_{\mathrm{3}}-\mathrm{4}{\mathit{u}}_{\mathrm{4}}-\mathrm{4}{\mathit{u}}_{\mathrm{5}}.\hfill \end{array}

From Eq. (17), we obtain the shape functions for the conforming crack front element with the crack front at edge \xi +\eta =\mathrm{1} as follows:

(18)\{\begin{array}{l}{N}_{\mathrm{scoll}}^{\mathrm{0}}(\xi ,\eta )=\xi (\mathrm{2}\xi -\mathrm{1})-\mathrm{0.5}(\mathrm{4}\xi -\mathrm{1})\sqrt{\mathrm{1}-\xi -\eta }\hfill \\ \hspace{0.17em}+\mathrm{0.5}(\mathrm{1}-\xi -\eta ),\hfill \\ N_{\mathrm{scoll}}^{\mathrm{1}}(\xi ,\eta )=(\mathrm{1}-\xi )(\mathrm{1}-\mathrm{2}\xi )-\mathrm{0.5}(\mathrm{3}-\mathrm{4}\xi )\sqrt{\mathrm{1}-\xi -\eta }\hfill \\ \hspace{0.17em}+\mathrm{0.5}(\mathrm{1}-\xi -\eta ),\hfill \\ N_{\mathrm{scoll}}^{\mathrm{2}}(\xi ,\eta )=-\sqrt{\mathrm{1}-\xi -\eta }+\mathrm{2}(\mathrm{1}-\xi -\eta ),\hfill \\ N_{\mathrm{scoll}}^{\mathrm{3}}(\xi ,\eta )=\mathrm{4}\xi (\mathrm{1}-\xi )-\mathrm{2}\sqrt{\mathrm{1}-\xi -\eta }+\mathrm{2}(\mathrm{1}-\xi -\eta ),\hfill \\ N_{\mathrm{scoll}}^{\mathrm{4}}(\xi ,\eta )=\mathrm{4}(\mathrm{1}-\xi )\sqrt{\mathrm{1}-\xi -\eta }-\mathrm{2}(\mathrm{1}-\xi -\eta ),\hfill \\ N_{\mathrm{scoll}}^{\mathrm{5}}(\xi ,\eta )=\mathrm{4}\xi \sqrt{\mathrm{1}-\xi -\eta }-\mathrm{2}(\mathrm{1}-\xi -\eta ).\hfill \end{array}

In the application of the dual BEM for crack problems, due to the existence of hyper-singular integrals, the conforming crack front element should be transformed into a non-conforming one to retain Hölder’s continuity condition through the following linear transformation:

(19){\left[\begin{array}{l}{N}_{\mathrm{Nscoll}}^{\mathrm{0}}\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{1}}\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{2}}\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{3}}\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{4}}\hfill \\ N_{\mathrm{Nscoll}}^{\mathrm{5}}\hfill \end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{l}{N}_{\mathrm{scoll}}^{\mathrm{0}}\hfill \\ N_{\mathrm{scoll}}^{\mathrm{1}}\hfill \\ N_{\mathrm{scoll}}^{\mathrm{2}}\hfill \\ N_{\mathrm{scoll}}^{\mathrm{3}}\hfill \\ N_{\mathrm{scoll}}^{\mathrm{4}}\hfill \\ N_{\mathrm{scoll}}^{\mathrm{5}}\hfill \end{array}\right]}^{\mathrm{T}}{\left[\begin{array}{llllll}{N}_{\mathrm{scoll}}^{\mathrm{0}}\left(P_{\mathrm{0}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{1}}\left(P_{\mathrm{0}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{2}}\left(P_{\mathrm{0}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{3}}\left(P_{\mathrm{0}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{4}}\left(P_{\mathrm{0}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{5}}\left(P_{\mathrm{0}}\right)\hfill \\ N_{\mathrm{scoll}}^{\mathrm{0}}\left(P_{\mathrm{1}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{1}}\left(P_{\mathrm{1}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{2}}\left(P_{\mathrm{1}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{3}}\left(P_{\mathrm{1}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{4}}\left(P_{\mathrm{1}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{5}}\left(P_{\mathrm{1}}\right)\hfill \\ N_{\mathrm{scoll}}^{\mathrm{0}}\left(P_{\mathrm{2}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{1}}\left(P_{\mathrm{2}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{2}}\left(P_{\mathrm{2}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{3}}\left(P_{\mathrm{2}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{4}}\left(P_{\mathrm{2}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{5}}\left(P_{\mathrm{2}}\right)\hfill \\ N_{\mathrm{scoll}}^{\mathrm{0}}\left(P_{\mathrm{3}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{1}}\left(P_{\mathrm{3}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{2}}\left(P_{\mathrm{3}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{3}}\left(P_{\mathrm{3}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{4}}\left(P_{\mathrm{3}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{5}}\left(P_{\mathrm{3}}\right)\hfill \\ N_{\mathrm{scoll}}^{\mathrm{0}}\left(P_{\mathrm{4}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{1}}\left(P_{\mathrm{4}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{2}}\left(P_{\mathrm{4}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{3}}\left(P_{\mathrm{4}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{4}}\left(P_{\mathrm{4}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{5}}\left(P_{\mathrm{4}}\right)\hfill \\ N_{\mathrm{scoll}}^{\mathrm{0}}\left(P_{\mathrm{5}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{1}}\left(P_{\mathrm{5}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{2}}\left(P_{\mathrm{5}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{3}}\left(P_{\mathrm{5}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{4}}\left(P_{\mathrm{5}}\right)\hfill & N_{\mathrm{scoll}}^{\mathrm{5}}\left(P_{\mathrm{5}}\right)\hfill \end{array}\right]}^{-\mathrm{1}},

where N_{\mathrm{Nscoll}}^m\left(m=\mathrm{0},\hspace{0.17em}\mathrm{1},\hspace{0.17em}\mathrm{...},\hspace{0.17em}\mathrm{5}\right) denotes the shape functions of the non-conforming crack front element, P_m\left(m=\mathrm{0},\hspace{0.17em}\mathrm{1},\hspace{0.17em}\mathrm{...},\hspace{0.17em}\mathrm{5}\right) denotes the parameter coordinate of the interpolation nodes in (\xi ,\eta )coordinate systems, which can be obtained by the position parameter. In this work, Eqs. (11) and (19) are applied to approximate the distribution of displacements near the crack front.

The crack front element may lie on another edge. In this case, we can modify the order of the nodes or construct another type of shape function in the same manner.

Calculation of SIFs

Quadrilateral and triangular non-conforming crack front elements can be conveniently integrated into the dual boundary element method (DBEM). To compute the SIFs along the crack front, we construct a local coordinate system for Point O on the crack front, as shown in Fig. 4. The classic one-point formula is used to compute the SIFs, which is similar to that in Eq. (20) [⁵]:

(20)\{\begin{array}{l}{K}_{\mathrm{I}}={\displaystyle \frac{E}{\mathrm{4}(\mathrm{1}-{\nu }^{\mathrm{2}})}}\sqrt{{\displaystyle \frac{\mathrm{\pi }}{\mathrm{2}r\cos \phi }}}\left(u_b^{P^{+}}-u_b^{P-}\right),\hfill \\ K_{\mathrm{II}}={\displaystyle \frac{E}{\mathrm{4}(\mathrm{1}-{\nu }^{\mathrm{2}})}}\sqrt{{\displaystyle \frac{\mathrm{\pi }}{\mathrm{2}r\cos \phi }}}\left(u_n^{P^{+}}-u_n^{P-}\right),\hfill \\ K_{\mathrm{III}}={\displaystyle \frac{E}{\mathrm{4}(\mathrm{1}+\nu )}}\sqrt{{\displaystyle \frac{\mathrm{\pi }}{\mathrm{2}r\cos \phi }}}\left(u_t^{P^{+}}-u_t^{P-}\right),\hfill \end{array}

where E is Young’s modulus, \nu is Poisson’s ratio, and \phi is the angle between P^{+}O (or P^{-}O) and \mathit{n}. P is the point on the crack surface, which may be P₁, P₂, or P₃. By using \phi, the distance r^{P_{\mathrm{2}}}=P_{\mathrm{2}}O, r^{P_{\mathrm{3}}}=P_{\mathrm{3}}O, K^{p_{\mathrm{2}}} and K^{p_{\mathrm{3}}} in Eq. (21) can be obtained. Then, a two-point formula considering the whole element can be derived after performing linear extrapolation between \left(r^{P_{\mathrm{2}}},K^{p_{\mathrm{2}}}\right) and \left(r^{P_{\mathrm{3}}},K^{p_{\mathrm{3}}}\right) at location r=\mathrm{0}. We can compute the SIFs at the crack front by using Eq. (21):

(21)K^O={\displaystyle \frac{r^{p_{\mathrm{3}}}{K}^{p_{\mathrm{2}}}-r^{p_{\mathrm{2}}}{K}^{p_{\mathrm{3}}}}{r^{p_{\mathrm{3}}}-r^{p_{\mathrm{2}}}}}.

Fig.4 Local coordinate system at Point O

Full size|PPT slide

Numerical examples

Several numerical examples are presented in this section to verify the accuracy and efficiency of the newly developed crack front elements. In the first example, the two constructed non-conforming crack front elements are tested in the DBEM analysis of a penny-shaped crack under uniform or linear polynomial tractions. The crack open displacements, \mathrm{\Delta }u, are obtained for a comparison with existing results. Efficiency is also verified. In the second example, a penny-shaped crack under uniform inclined traction is analyzed with DBEM. Three types of SIFs, which contain K_I, K_II, and K_III, are computed, and their definitions can be found in Ref. [²⁷]. The numerical results obtained by the two presented elements are compared with the exact results. In the final example, straight edge crack problems are compared with existing finite element analysis and other results.

Penny-shaped crack under uniform or linear polynomial tractions

This example concerns a penny-shaped crack embedded in an infinite space, as illustrated in Fig. 5. The boundary of the penny-shaped crack is described as

(22)x=a\cos \omega ,\hspace{0.17em}\hspace{0.17em}y=a\sin \omega ,\hspace{0.17em}\hspace{0.17em}z=\mathrm{0},

where \omega varies from 0 to \mathrm{2}\mathrm{\pi }, and a is the radius. The crack is opened by imposing traction {\sigma }_{zz}(x,y) along the z-direction symmetrically to the upper and lower surfaces of the penny-shaped crack. In this example, a=1.0, Young’s modulus E=1.0, and Poisson’s ratio \nu =\mathrm{0.25} are evaluated. Three traction boundary conditions, namely, {\sigma }_{zz}(x,y)={\sigma }_{\mathrm{0}}, {\sigma }_{zz}(x,y)=-{\sigma }_{\mathrm{0}}x, and {\sigma }_{zz}(x,y)=-{\sigma }_{\mathrm{0}}y with constant {\sigma }_{\mathrm{0}}, are imposed on the penny-shaped crack. The reference solutions to the considered problems are \mathrm{2}{\sigma }_{\mathrm{0}}/\sqrt{\mathrm{\pi }}, \mathrm{4}{\sigma }_{\mathrm{0}}x/\left(\mathrm{3}a\sqrt{\mathrm{\pi }}\right), or \mathrm{4}{\sigma }_{\mathrm{0}}y/\left(\mathrm{3}a\sqrt{\mathrm{\pi }}\right), respectively [²⁸]. In this example, 96 nine-node non-conforming elements are used inside the crack surface and 32 nine-node non-conforming crack front elements are employed near the crack front, as shown in Fig. 5(a).

Fig.5 Mesh of the crack surface. (a) Non-conforming quadrilateral element; (b) non-conforming triangular element

Full size|PPT slide

Figure 6 shows that the K_I calculated by non-conforming crack front elements of the quadrilateral type agrees well with the reference solutions. The largest errors are approximately 0.2%, 1.7%, and 2%, respectively.

Fig.6 Results of the normalized SIF K_{\mathrm{I}} calculated by non-conforming crack front elements of the quadrilateral type

Full size|PPT slide

Then, we perform an analysis using the triangular elements. As shown in Fig. 5(b), on the crack surface, 162 six-node non-conforming elements, including 26 six-node non-conforming crack front elements, are involved in the analysis.

Figure 7 indicates that the K_I computed by employing non-conforming crack front elements of the triangular type agrees well with the reference solutions. The largest errors are about 0.9%, 2%, and 2%, respectively.

Fig.7 Results of the normalized SIF K_{\mathrm{I}} calculated by non-conforming crack front elements of the triangular type

Full size|PPT slide

Figure 8 indicates that the \mathrm{\Delta }u calculated by non-conforming crack front elements of the quadrilateral and triangular type agrees well with the reference solutions [²⁹]. The largest errors are about 1.5% in the two cases.

Fig.8 Results of \mathrm{\Delta }u over the crack surface obtained using our method

Full size|PPT slide

Figure 9 shows that even when position parameter k varies from 0.1 to 0.3, the SIFs still retain considerable accuracy with the analytical solution. Notably, k=\mathrm{1}-\lambda for crack front elements of the quadrilateral type. The largest errors are about 1.2% in the two cases.

Fig.9 Influence of position parameter k

Full size|PPT slide

When the two types of crack front elements are used, high accuracy can be achieved within 1% by using fewer elements than Liu’s work [³⁰]. In Liu’s work, 4704 elements are employed and the error is about 1.66%, while less than 200 elements are employed in our work, and the errors are within 1% stably.

Penny-shaped crack under uniform inclined traction

This example concerns the same crack as that in the first example. The load in this example, however, is different from that in the first numerical example. In this example, the crack is under a uniform inclined traction, as illustrated in Fig. 10. The parameters illustrated in Fig. 10 are evaluated by a=1.0, \sigma =\mathrm{1.0}, E=1.0, \mathit{\gamma }=\mathrm{\pi }/\mathrm{4}, and \nu=0.25. Here, \mathit{\gamma } is the angle between the uniform inclined traction and coordinate plane xoy. The same meshes as those in the first example are applied. The exact SIFs for this problem can be found in Tada et al.’s work [³¹]. The numerical and exact results are illustrated and compared in Figs. 11 and 12. The results obtained by our method agree well with the exact solutions. The largest errors are about 0.3%, 0.7%, 0.7%, 1%, 2%, and 2%, respectively.

Fig.10 Penny-shaped crack under uniform inclined traction

Full size|PPT slide

Fig.11 Results of three normalized SIFs K_{\mathrm{I}}, K_{\mathrm{II}}, and K_{\mathrm{III}} calculated by non-conforming crack front elements of the quadrilateral type

Full size|PPT slide

Fig.12 Results of three normalized SIFs K_{\mathrm{I}}, K_{\mathrm{II}}, and K_{\mathrm{III}} calculated by non-conforming crack front elements of the triangular type

Full size|PPT slide

Straight-edge crack compared with existing finite element analysis

In this example, a single edge crack, as illustrated in Fig. 13, is considered. The geometry parameters are c/b=1, t/c=3,and h/b=6. Uniform tension is imposed on two ends. The results are compared with those of the plane strain solution obtained by Mi and Aliabadi [⁵], Murakami and Keer [³²], and Raju and Newman [³³]. A total of 96 and 88 elements are used for the crack faces. The crack front elements proposed in this study are employed. In the work of Mi and Aliabadi [⁵], BEM was employed, whereas in the work of Raju and Newman [³³], FEM was employed. The plane strain solution is 2.827, which can be found in the work of Murakami and Keer [³²]. In our work, when z=0, the result is 2.2834 when crack front elements of the quadrilateral type are used; the result is 2.828 when crack front elements of the triangular type are employed. In Raju and Newman’s work [³³], the SIF result is 2.776. Figure 14 shows that the results in the two cases agree well with those of Mi and Aliabadi [⁵] when z/ω changes from 0 to 0.4. When z/ω approaches 0.5, the results in the two cases agree well with those of Raju and Newman [³³]. This example demonstrates the credibility of our method.

Fig.13 Geometry model of the singular edge crack and mesh

Full size|PPT slide

Fig.14 Normalized stress intensity factors along the crack front for the single edge crack

Full size|PPT slide

Conclusions

This study developed a family of non-conforming crack front elements of quadrilateral and triangular types in DBEM for 3D crack problems. In the deduction of the special shape functions, the asymptotic behavior of the distribution of the displacements and tractions near the crack front is considered. In the neighborhood of the crack front, the traditional shape function is expanded into Taylor expansion. The expanded shape function is then modified by considering the asymptotic behavior of the physical quantities. For the non-conforming quadrilateral crack front elements, the relations between the 2D and 3D Lagrange-type elements are determined. For the non-conforming triangular crack front elements, the transformation matrix between the non-conforming and conforming elements is introduced in consideration of the positions of the collocation nodes. Then, numerical examples using these non-conforming crack front elements are provided to verify the efficiency of our method. The numerical results are in good agreement with existing results. The proposed method may serve as a reference for addressing V-notch problems.

Acknowledgements

This work was supported in part by the National Natural Science Foundation of China (Grant Nos. 11602229 and 11602082), Hunan Provincial Natural Science Foundation of China (Grant No. 2017JJ3061), and Key Scientific and Technological Project of Henan Province (Grant No. 192102210227).

References

1 Zhang D H, Li Y Q, Xie G Z, Digital image correlation method for measuring deformations of vinyl chloride-coated metal multilayer sheets. Modern Physics Letters B, 2019, 33(5): 1950050

2 Zhang D H, Li Y Q, Liu J X, A novel 3D optical method for measuring and evaluating springback in sheet metal forming process. Measurement, 2016, 92: 303–317

3 Zhang D H, Xie G Z, Li Y Q, Strain and mechanical properties of the VCM multilayer sheet and their composites using digital speckle correlation method. Applied Optics, 2015, 54(25): 7534–7541

4 Zhang D H, Bai D P, Liu J B, Formability behaviors of 2A12 thin-wall part based on DYNAFORM and stamping experiment. Composites Part B: Engineering, 2013, 55: 591–598

5 Mi Y, Aliabadi M H. Dual boundary element method for three-dimensional fracture mechanics analysis. Engineering Analysis with Boundary Elements, 1992, 10(2): 161–171

6 Sladek V, Sladek J, Tanaka M. Nonsingular BEM formulations for thin-walled structures and elastostatic crack problems. Acta Mechanica, 1993, 99(1–4): 173–190

7 Zhang Y M, Gu Y, Chen J T. Boundary element analysis of the thermal behaviour in thin-coated cutting tools. Engineering Analysis with Boundary Elements, 2010, 34(9): 775–784

8 Yao Z H, Wang H T. Some benchmark problems and basic ideas on the accuracy of boundary element analysis. Engineering Analysis with Boundary Elements, 2013, 37(12): 1674–1692

9 Zhang J M, Lin W C, Dong Y, A double-layer interpolation method for implementation of BEM analysis of problems in potential theory. Applied Mathematical Modelling, 2017, 51: 250–269

10 Feng S Z, Han X, Wang G. An efficient on-line algorithm for the optimal design of multi-material structures under thermal loads. International Journal of Thermal Sciences, 2018, 132: 567–577

11 Feng S Z, Cheng Y H. An element decomposition method for heat transfer analysis. International Journal of Heat and Mass Transfer, 2018, 123: 437–444

12 Feng S Z, Bordas S P A, Han X, A gradient weighted extended finite element method (GW-XFEM) for fracture mechanics. Acta Mechanica, 2019 (in press)

13 Cheng C, Niu Z, Recho N. Analysis of the stress singularity for a bi-material V-notch by the boundary element method. Applied Mathematical Modelling, 2013, 37(22): 9398–9408

14 Cheng C, Ge S, Yao S, et al. Thermal stress singularity analysis for V-notches by natural boundary element method. Applied Mathematical Modelling, 2016, 40(19–20): 8552–8563

15 Zhu B J, Qin T Y. Application of hypersingular integral equation method to three-dimensional crack in electromagnetothermoelastic multiphase composites. International Journal of Solids and Structures, 2007, 44(18–19): 5994–6012

16 Henshell R D, Shaw K G. Crack tip finite elements are unnecessary. International Journal for Numerical Methods in Engineering, 1975, 9(3): 495–507

17 Lv J H, Jiao Y Y, Wriggers P, Efficient integration of crack singularities in the extended finite element method: Duffy-distance transformation and conformal preconditioning strategy. Computer Methods in Applied Mechanics and Engineering, 2018, 340: 559–576

18 Sáez A, Gallego R, Dominguez J. Hypersingular quarter-point boundary elements for crack problems. International Journal for Numerical Methods in Engineering, 1995, 38(10): 1681–1701

19 Hong H K, Chen J T. Derivations of integral equations of elasticity. Journal of Engineering Mechanics, 1988, 114(6): 1028–1044

20 Mi Y, Aliabadi M H. Discontinuous crack-tip elements: Application to 3D boundary element method. International Journal of Fracture, 1994, 67(3): R67–R71

21 Pan E, Yuan F G. Boundary element analysis of three-dimensional cracks in anisotropic solids. International Journal for Numerical Methods in Engineering, 2000, 48(2): 211–237

22 Xie G Z, Zhang J, Huang C, A direct traction boundary integral equation method for three-dimension crack problems in infinite and finite domains. Computational Mechanics, 2014, 53(4): 575–586

23 Li J, Feng W Z, Gao X W. A dual boundary integral equation method based on direct evaluation of higher order singular integral for crack problems. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(2): 387–398 (in Chinese)

24 Xie G Z, Zhang D, Meng F, Calculation of stress intensity factor along the 3D crack front by dual BIE with new crack front elements. Acta Mechanica, 2017, 228(9): 3135–3153

25 Ariza M P, Saez A, Dominguez J. A singular element for three-dimensional fracture mechanics analysis. Engineering Analysis with Boundary Elements, 1997, 20(4): 275–285

26 Ariza M P, Dominguez J. Boundary element formulation for 3D transversely isotropic cracked bodies. International Journal for Numerical Methods in Engineering, 2004, 60(4): 719–753

27 Aliabadi M H, Rooke D P. Numerical Fracture Mechanics. Boston: Kluwer Academic Publishers, 1991

28 Shah R C, Kobayashi A S. Stress intensity factor for an elliptical crack under arbitrary normal loading. Engineering Fracture Mechanics, 1971, 3(1): 71–96

29 Chen L S, Kuang J H. A displacement extrapolation method for determining the stress intensity factors along flaw border. International Journal of Fracture, 1992, 57(4): R51–R58

30 Liu Y. On the displacement discontinuity method and the boundary element method for solving 3-D crack problems. Engineering Fracture Mechanics, 2016, 164: 35–45

31 Tada H, Paris P C, Irwin G R. The Stress Analysis of Cracks Handbook. Vol. 130. New York: ASME, 2000

32 Murakami Y, Keer L M. Stress intensity factors handbook, Vol. 3. Journal of Applied Mechanics, 1993, 60(4): 1063

33 Raju I S, Newman J C Jr. Three Dimensional Finite-Element Analysis of Finite-Thickness Fracture Specimens. NASA Technical Note, NASA TN D-8414. 1977