## Abstract

Super singular element method (SSEM) is an accurate method devoted to stress singularity analysis accomplished by stress analysts who have little or no experience of fracture mechanics. The basic idea of SSEM is to transform nodal displacements inside a superelement to global variables, which are the coefficients of the Williams eigenfunction series (ES) that accurately represents crack tip singularity. The stress singularity is modelled by the superelement and the Williams' ES, hence refining the finite-element meshes near the singular point and creation of crack faces can be avoided. The stress intensity factor, *T*-stress, and the higher order coefficients can be obtained directly from the global variables without any post-processing. The implementation of the SSEM is rather easy; since only matrix multiplication and transformation are required. The SSEM has been incorporated into standard finite-element program ABAQUS. Several two-dimensional crack problems have been performed to compare the SSEM with published results.

## 1. Introduction

It is well known that stress singularities occur at crack tips, corners and at non-isotropic multi-material interfaces. The computation of the stress intensity factors (SIFs) at such singular points has become an important consideration in fracture mechanics. One of the earlier workers in this field is Williams (1957). He obtained the displacement and stress fields at the vicinity of a two-dimensional crack tip by the eigenfunction expansion method. The first term of the eigenfunction expansion represents the stress singularity and its coefficient is called SIF. This coefficient is dependant on the applied load, the shape of the body and the crack length. The SIF connects with the fracture theory through the postulate that a crack will grow when the SIF reaches a critical value and this critical value is an independent material property. The second term of the eigenfunction expansion represents a constant stress acting parallel to the crack plane and is called *T*-stress. The sign and magnitude of *T*-stress substantially change the size and shape of the plane strain contour (Larsson & Carlsson 1973; Rice 1974). Cotterell & Rice (1980) have evaluated the *T*-stress in front of a crack and proposed a criterion that determined the direction of its propagation. Cracks with *T*<0 are stable in mode I, while those with *T*>0 kink. Tvergaard & Hutchinson (1994) have shown that *T*-stress affects mode I crack growth resistance in a ductile solid. Zhao & Chen (1998) have shown that *T*-stress affects microcrack toughening. Betegón & Hancock (1991) showed that *T*-stress can be used as a second parameter in conjunction with *J*-integral to characterize fully the crack-tip field in contained yielding. Similar to SIF, the values of *T*-stress are dependent on the type of loading, the relative crack length and overall geometry.

It was found that the higher order terms also contribute to the size effect of quasi-brittle materials (Dyskin 1997). Hence it is important to obtain the leading coefficients of the Williams expansion for the cracked geometries under consideration. Several analytical and numerical methods have been developed to evaluate the SIF and *T*-stress for cracked configurations (Fett 1997; Fleming *et al*. 1997; Tan & Wang 2003). However, only few methods are available in the literature which can be used to evaluate the higher order terms in the expansion series. Those methods include the higher order hybrid crack element (HCE; Karihaloo & Xiao 2001*a*) and the scaled boundary finite-element method (SBFEM; Wolf 2003). The HCE is formulated from a simplified variational principle using truncated asymptotic crack tip displacement and stress expansions. In the early development of the HCE, the coefficients in the Williams expansion that represent rigid body translations and rotation were neglected in the formulation of the HCE because they did not contribute to stresses. Consequently, the displacement field was naturally incompatible with the surrounding finite element (FE; Xiao *et al*. 2004). Recently, Xiao & Karihaloo (2007*a*) have recovered these coefficients by an indirect method that involves the application of a least-squares procedure to numerically generated data.

The SBFEM (Wolf 2003) reduces the governing partial differential equations to a set of second-order ordinary linear differential equations for the displacement, which can be solved analytically. The reduction is accomplished by introducing shape functions in one coordinate direction, while working analytically in the other direction. The SBFEM introduces a radial coordinate system based on a scaling centre and a defined curve, usually the boundary of the problem domain, which is discretized with elements based on the standard FE interpolation functions. If the scaling centre is located at the crack tip, the SBFEM semi-analytical series converges to the Williams expansion. However, the authors have stated that the rather complicated mathematics of the SBFEM in comparison with the finite-element method has obviated the uptake of the method by other researchers. A new virtual work formulation and model interpretation of the method (Deeks & Wolf 2002) have increased the transparency of the method considerably.

The super singular element method (SSEM) is a semi-analytical method, which combines the efficiency and accuracy of analytical solutions and the agility of standard finite-element method. The SSEM is an improved version of the fractal-like finite-element method (Leung *et al*. 2000; Tsang *et al*. 2004, 2007). The SSEM is designed to be easy to use for engineering purposes such that crack analysis can be done by a normal stress analyst. The basic idea of SSEM is to transform the nodal displacements inside a superelement to global variables, which are the coefficients of the Williams eigenfunction series (ES). Unlike some of the special singular element methods (Tong *et al*. 1973; Fleming *et al*. 1997) that were formulated to use a fixed number of leading terms in the eigenfunction expansion series, the number of eigenfunctions used in SSEM analyses is a parameter and can be varied. Also most of the special singular elements have compatibility problems, where the displacements do not match at the interelement boundary (Xiao *et al*. 2004). However, the SSEM does not seem to suffer from any displacement compatibility problem. The expansion coefficients are the output of the solution and depend on the far-field conditions. Crack problems in different geometry and subjected to different loading conditions will give a different set of expansion coefficients.

The conventional FEs are used within the SSEM; no new special elements have to be derived or created. The implementation and coding of the SSEM are rather easy as only matrix multiplication and transformation are needed; no cumbrous mathematics is involved as in SBFEM (Chidgzey & Deeks 2005). After the transformation, the order of unknowns is reduced from a very large value to a small value, which depends on the size of the transformation matrix. The singular point is modelled by a superelement; hence refining the finite-element meshes near the singular point and creation of crack faces can be avoided. Also the crack tip location can be anywhere within the SSEM, hence different crack length and different crack orientation within a structure can be analysed by using a single superelement. The SIF, *T*-stress and higher order terms can be obtained directly from the expansion coefficients without any post-processing.

## 2. Super singular element method

Superelement has been previously used by many authors (Dashevskii & Rotter 1984; Islam 1994; Yosibash & Schiff 1997) for crack and notch singularity analyses. One of the main advantages of using a superelement is that mesh refinement around a crack tip could be avoided, which can be very time consuming. The formation of the SSEM is more advanced than the traditional superelement technique (Tkachev 2000). In the traditional superelement technique mesh generation is not included; hence different superelements are required for different problems. However, the SSEM has automatic mesh generation capability (self-similar meshing) and hence a superelement can be applied to a number of different crack problems. Furthermore in traditional superelements, all the internal nodal displacements are eliminated, but in the SSEM all the internal nodal displacements are represented by a small set of variables.

Furthermore, an analytical representation of the crack tip singularity is embedded in the SSEM via the use of the Williams eigenfunction expansion series. However, other approaches reported in the literature such as the infinitely small element method (ISEM) of Go & Lin (1990) and the infinite-element method (IEM) of Liu & Chiou (2005), which also employ a self-similar method to create a theoretical infinite number of layers of elements around the crack tip, do not include any theoretical or numerical representation of the crack tip singularity. These approaches seem to rely solely on the use of infinitely small elements to represent the singularity. But Tong & Pian (1973) have shown that it is not possible to accurately model the crack tip singularity merely by the use of infinitely small elements; it is necessary to accurately represent the singularity using appropriate expressions. Therefore, the SSEM is considerably different from ISEM and IEM in that SSEM includes an accurate analytical representation of the crack tip singularity, whereas ISEM and IEM do not. Also, the SSEM provides accurate values of the SIFs directly from the first coefficient and the *T*-stress from the second coefficient of the eigenfunctions expansion series without the need for any post-processing of results. Similarly, the SSEM also provides the higher order terms of the series directly. It should be noted that both HCE and SBFEM also include accurate representations of the crack tip singularity.

In the implementation of the SSEM, infinite layers of self-similar meshes are created automatically within a super singular element, with the singular point as the self-similar centre. The large numbers of nodal displacement unknowns are transformed element by element to a smaller set of unknowns. The self-similar meshing and transformation are two of the main features of the SSEM and are described in more details in §§2*a*,*b*.

### (a) Self-similar meshing

Consider the two-dimensional domain *Ω* inside a boundary *Γ*_{0}, which contains a crack with crack tip singularity at *O*, shown in figure 1. Self-similar meshes can always be generated in the domain *Ω* by taking *O* as a centre of similarity. Let 0<*ξ*<1 be the similarity ratio. A set of curves {*Γ*_{1}, *Γ*_{2}, …, *Γ*_{n}, …}, which have similar shape as *Γ*_{0}, with proportionality constant {*ξ*, *ξ*^{2}, …, *ξ*^{n}, …} is generated in the domain *Ω*. The region between any two consecutive curves *Γ*_{n−1} and *Γ*_{n} is denoted as *n*th layer. Each layer is divided into the same number of elements such that every layer has a similar pattern. The number of elements in each layer is dependent on the element type and number of nodes on the boundary *Γ*_{0}. For a fixed value of similarity ratio *ξ*, the mesh profile can be uniquely determined in the domain *Ω* using self-similar meshing. A value of *ξ* close to 1 produces more fine meshes.

Mesh preparation can be a very difficult task in crack problems (e.g. to model a crack tip singularity and to create crack faces). By using the self-similar mesh, an infinite number of layers, and hence infinite number of meshes, can be generated for the crack problems. Thus complex mesh preparation can be avoided.

It is well known that stiffness matrices of any two-dimensional isoparametric elements of similar shapes are independent of their dimensions. Consequently, only the first layer of elements is enough to define all element stiffness matrices within the superelement.

### (b) Transformation

Using self-similar meshing, a layer of elements is enough to determine all the element stiffness matrix inside *Ω* and hence the global stiffness matrix. However, we end up with a system of algebraic equations with an infinite number of unknowns. The number of unknowns can be reduced significantly by executing a transformation procedure.

Let (*r*, *θ*) be the polar coordinate with centre at *O* (figure 1). The displacement at any location inside *Γ*_{0} (except *r*=0) can be determined using the Williams expansion series, that is(2.1)where *f*_{m} are the eigenfunctions, *λ*_{m} are the eigenvalues and *a*_{m} are the unknown coefficients in the series. The Williams expansion is an analytical solution of crack problems which satisfies all the boundary conditions on the crack faces and has the correct order of singularity around the crack tip. In two-dimensional cases, the eigenfunctions *f*_{m} and the eigenvalues *λ*_{m} can be determined analytically as elementary functions. The series coefficients *a*_{m} depend on the loading and geometry of the crack problems.

If the first *M* terms in the Williams series are used, equation (2.1) can be rewritten in matrix form, i.e.(2.2)where ** P** is called the transformation matrix and

**is called the generalized variable, which contains unknown coefficients in the expansion series.**

*a*For each self-similar element in the first layer, we have(2.3)where *K*_{e} is the elemental stiffness matrix; *u*_{e} is nodal displacement; and *R*_{e} is the corresponding nodal forces. Applying the transformation using equation (2.2), equation (2.3) becomes(2.4)The transformation matrix only needs to be determined once for every element in the first layer. There is a relationship between the transformation matrix of the first layer and that of the *n*-th layer, that is(2.5)where *G*_{n} is a diagonal matrix that depends on the similarity ratio *ξ* and the value of *n*.

After the transformation, the infinite numbers of unknowns reduce to the generalized variable ** a** and the size of

**is dependent on the number of Williams expansion terms used in the transformation. For each set of self-similar elements towards the centre, the assembly procedure in FE analysis becomes an infinity sum,(2.6)Each component in the infinity sum is a geometric series and hence the infinity sum can be found analytically, i.e.(2.7)where**

*a***is a square matrix, which is a function of**

*H**ξ*only; is a diagonal matrix; and the operator ⊗ is the component product of the two matrices. The final equilibrium equation inside the superelement becomes(2.8)where(2.9)and

*J*is the total number of self-similar elements. Now we have transformed an infinite number of nodal displacement variables to a generalized variable

**. Also all the calculations are done using the first layer only.**

*a*## 3. Numerical examples

SSEM has been implemented into ABAQUS (2006) via user element subroutine. Four crack problems will be studied here including single-edge crack under uniform traction, three-point bending, mixed-mode loading and wedge splitting. A total of six cases are investigated under these four crack problems.

### (a) Single-edge crack under uniform traction

Consider a single-edge cracked plate under uniform traction (figure 2*a*). The width and height of the plate are *w* and *h* units, respectively. The crack length is *a* units. A single 121-noded superelement and a single 6-noded triangle element are used in this study (figure 2*b*). The triangle element is used mainly for the purpose of displacement boundary conditions. Conventional second-order elements are used to generate the SSEM.

The SIFs and *T*-stress obtained by SSEM for aspect ratio *h*/*w*=2*w*=2 are presented in table 1. All the results are computed using a single mesh configuration as shown in figure 2*b* and they are in good agreement with the results determined by Larsson & Carlsson (1973), Yang & Ravi-Chandar (1999) and Tada *et al*. (2000). The displacement profiles determined from the finite-element analysis and the ES at the superelement boundary for *a*/*w*=0.5 are shown in figure 3. Both displacement results are in good agreement.

The leading five coefficients for the single-edge cracked plate subjected to uniaxial tension have been determined by HCE (Karihaloo & Xiao 2001*a*) and SBFEM (Chidgzey & Deeks 2005). The geometry of the plane is taken to be *a*=1 and *w*=*h*/2=4. Figure 2*c* shows the mesh configuration used by SSEM. A 49-noded superelement is surrounded by 764 conventional FEs. Table 2 shows the results obtained by SSEM and is compared with the coefficients computed in the above references. Very close agreement has been obtained between SSEM and SBFEM.

### (b) Three-point bending

Consider a single-edge cracked plate under three-point bending (figure 4). The geometry of the plate is 4*w*=*h*=4. *P* is the applied load per unit thickness. Figure 5 shows the meshes used in the analysis for crack length *a*/*w*=0.1 to 0.6. Each mesh has 1585 conventional FEs. Each crack tip is enclosed by a 33-noded superelement. The crack tip asymptotic field has been obtained by HCE (Karihaloo & Xiao 2001*b*). Table 3 shows non-dimensional results for the leading five terms of the coefficients in the eigenfunction expansion series obtained from SSEM and compared with HCE. Relative differences are generally less than 5%, except when the magnitude of a term is very small and close to zero.

To examine the size effect of superelement, now the crack tip is enclosed by 20 layers of conventional FE. The sizes of elements decrease towards the crack tip. To apply SSEM, some layers have to be deleted and replaced by a superelement. Five different meshes have been used to model the plate with crack length *a*/*w*=0.5 and they are shown in figure 6. The numbers of layers deleted in figure 6*a–e* are 5, 10, 15, 19 and 20 (i.e. without any layers of elements), respectively. Each superelement has 161 nodes. The leading five terms of coefficients in the series are shown in table 4. The coefficients shown in table 4 suggest that the size of the superelement has little effect on the results. Also shown in table 4 is the computational time for each case. As expected, the computational time is reduced with smaller number of layers.

There is an advantage in using a bigger size of SSEM such that different crack lengths can be analysed by using one mesh profile (i.e. without remeshing). Table 5 shows numerical results for various crack lengths using the mesh shown in figure 6*e*. When table 5 is compared with table 3, it shows that very accurate solutions can be obtained, up to the coefficient of the fifth term, for *a*/*w*=0.2 to 0.6. However, for *a*/*w*=0.1 only the first term can be determined accurately. This is because the crack length is too small when compared with the superelement length. However, as shown previously when a local approach is used as illustrated in figure 5*a*, the SIF value computed is very accurate (table 3).

### (c) Mixed-mode loading

Consider a central slant crack of length 2*a* in a rectangular plate of height 2*h* and width 2*w* with angle of inclination *θ* (figure 7*a*). The mixed-mode loading is achieved by a uniform tensile stress applied at the ends of the plate. Two different meshes, regular mesh and irregular mesh, have been depicted in figure 7*b*,*c* for the study of a crack with length 2*a*=1.2 in a plate of aspect ratio *h*/*w*=2*w*=2 and the angle of inclination *θ*=45°. There are 1380 conventional FEs in the regular mesh and 1216 conventional FEs in the irregular mesh. The numbers of nodes in the superelements used in the regular and the irregular meshes are 117 and 49, respectively. A higher number of nodes were used for the superelement of the regular mesh which was designed to enable the study of cracks with slant angles 0°<*θ*<90°. On the other hand, a small number of nodes were used for the superelement of the irregular mesh which was designed for the study of a single case of slant angle *θ*=45°. Table 6 shows the values of the SIF and the *T*-stress obtained from the regular and the irregular meshes. The table shows that for *θ*=45°, good agreements have been obtained for both meshes. The discrepancies are less than 0.1% for SIF results and 1.1% for *T*-stress results.

Table 7 presents the mixed-mode SIF values and *T*-stress values for various crack lengths and angles of inclination. The regular mesh shown in figure 7*b* is used for all the calculations. The SIF results have been compared with published results from Kitagawa & Yuuki (1977), who used a conformal mapping (CMAP) approach. The SSEM results agree well with published results, except for *K*_{II} result for *a*/*w*=0.8 and *θ*=15°. However, by considering the rate of change, it seems that the CMAP result for this case is in error.

The leading five coefficients for an angle cracked in an infinite plate under tension have been determined by HCE (Xiao & Karihaloo 2007*b*). The infinite plate is modelled by 2*w*=2*h*=100 and half crack length *a*=0.2. The analytical solutions of the SIFs for an infinite plate can be found in Sih & Liebowitz (1968). A uniform mesh of 100×100 square elements is used for all crack angles. The mesh configurations for the SSEM are the same as in Xiao & Karihaloo (2007*b*) and are shown in figure 8. The node number in the superelements enclosing the crack tips is varied from 25 to 33. A node on the middle of the right edge is fixed in both the vertical and horizontal directions, whereas a node on the middle of the left edge is fixed in the vertical direction only. Different boundary conditions were used in Xiao & Karihaloo (2007*b*).

Tables 8 and 9 show the leading five coefficients obtained by SSEM, compared with numerical results from HCE (Xiao & Karihaloo 2007*b*) and analytical solution (Sih & Liebowitz 1968). The SSEM gives identical results for both left and right tips. The SSEM results are in very good agreement with analytical solutions. Also the SSEM results are in close agreement with HCE, except for *a*_{3} coefficients. This is obviously because the HCE method does not give consistent results for the left and right crack tips as seen in table 8. This lack of consistency of the HCE for the left and right crack tips results for some coefficients has also been stated in Xiao & Karihaloo (2007*b*).

Tables 10 and 11 show the coefficients for *θ*=30° and 60°, respectively, using different mesh configurations. Table 10 shows the results obtained from meshes as in figure 8*b*,*c*,*f*, whereas table 11 shows the results obtained from meshes as in figure 8*c*,*d*,*f*. The SSEM results are insensitive to the mesh and the SSEM is less sensitive to the mesh configuration than HCE on *a*_{3} coefficient. Table 12 shows the coefficients for 15°≤*θ*≤75° using single mesh configuration as in figure 8*c*. Accurate coefficients can be determined using just a single mesh configuration in SSEM for various *θ*.

### (d) Wedge-splitting specimens

The leading coefficients for a wedge-splitting specimen with distributed and concentrated applied load have been determined by SBFEM (Chidgzey & Deeks 2005) and HCE (Karihaloo *et al*. 2003). The geometry of the wedge-splitting specimen is shown in figure 9*a*. The applied uniformly distributed load in figure 9*a* and point load in figure 9*b* are both assumed to be unity. A 49-noded superelement is used in the SSEM and the mesh configuration is shown in figure 9*c*.

Table 13 shows the leading coefficients for wedge splitting under distribution load obtained by SSEM, SBFEM (Chidgzey & Deeks 2005) and HCE (Karihaloo *et al*. 2003). The results from SSEM and SBFEM are in very good agreement. The greatest discrepancy is in the fourth coefficient, which has a difference of just over 1.4%. The HCE results seem less accurate when compared with those of the SSEM and SBFEM. For the first three terms, HCE gives results that are different from those of SSEM and SBFEM by approximately 4%. Furthermore, the sign of the fifth term obtained by HCE is opposite to that obtained by SSEM and SBFEM. Consequently, the difference between the HCE result and those of SSEM and SBFEM is very large.

The SSEM model for wedge-splitting specimen has been applied with point load and the leading coefficients are given in table 14. Again the results from SSEM and SBFEM are in excellent agreement. The greatest discrepancy is in fifth term and the difference is approximately 5%. However, the values of the fifth term are extremely small. Again, the order of magnitude and the sign of the fifth term obtained by HCE are inaccurate.

## 4. Conclusions

Super singular element has been proved very efficient to calculate the SIF and *T*-stress for various crack problems. The self-similar mesh together with the Williams ES reduces the number of unknowns significantly. The coefficients in eigenfunctions can be determined accurately without any post-processing. The mesh preparation for the crack problem can be simplified significantly by enclosing a crack tip by a superelement. Numerical results show very good agreement between the present method and published results.

## Footnotes

- Received October 12, 2007.
- Accepted April 23, 2008.

- © 2008 The Royal Society