A SFBEM–FEM coupling method for solving crack problems based on Erdogan fundamental solutions

The boundary element method (BEM) has proven to be an efficient approach for crack analysis in fracture mechanics, while its versatility in application to crack problems of complex structures with irregular boundaries deserves further attention. In this study, to improve the applicability to complex crack analysis, a cracked superelement is first established with BEM to model the near-crack region, and the crack problem is then solved within the frame of finite element method (FEM). The stiffness matrix of the cracked superelement is formulated using the spline fictitious boundary element method (SFBEM) based on the Erdogan fundamental solutions for an infinite plane with a single crack. The proposed superelement is further incorporated into a finite element mesh to simulate the behaviour of the crack zone, and the governing equation of the crack problem is finally established and solved using the typical procedure of FEM. After obtaining the nodal displacements of the superelement, the stress intensity factors (SIFs) of crack tips can be obtained by a backward analysis with SFBEM. The accuracy and efficiency of the proposed SFBEM–FEM coupling method are demonstrated by two numerical examples involving a rectangular plate with a central crack and a square plate with 100 horizontal cracks. The present approach is further applied to the analysis of the SIFs of multiple cracks exposed in a steel anchorage box for a hanger of a suspension bridge, which indicates the merging of the cracked superelements to the commercial FEM software is computationally efficient for the analysis of complex structures with multiple cracks.


Introduction
The analysis of stress intensity factors (SIFs) is of great importance to the assessment of damage tolerance and the prediction of crack propagation for engineering structures. As analytical methods are generally limited to solving simple crack problems [1], various numerical methods, e.g. the finite element method (FEM) [2], the extended finite element method (XFEM) [3] and the boundary element method (BEM) [4], have been developed for analysis of SIFs in engineering practice. Due to its high flexibility in modelling complex structures and imposing complex boundary conditions, the FEM is the most commonly used numerical method for solving crack problems of engineering structures. However, because of the crack-tip singularity, mesh refinement is required at crack tips so as to capture the singular behaviour of stresses, leading to extra burden of mesh generation and large computational cost of FEM for crack problems. The XFEM alleviates the shortcomings of the traditional FEM associated with mesh refinement of crack regions, in which local enrichment functions consisting of a discontinuous function and the near-tip asymptotic fields are incorporated into finite element formulation to describe the jump in displacement across the crack surface and the singularity at the crack tip, resulting in extra degrees-of-freedom (DOFs) in XFEM for crack problems [5,6].
As an alternative to FEM, the BEM has proven to be a powerful numerical method to solve crack problems, in which only boundary rather than domain discretization is required for analysis of SIFs. The existing BEMs applied to crack problems can be roughly classified into two categories depending on the fundamental solutions used. For the first category, the Kelvin fundamental solutions [7] for noncrack problems are adopted in the formulation of BEM, in which the stress boundary conditions along crack surfaces need to be satisfied and special boundary elements are required at crack tips to reflect the singular behaviour [8,9]. For the second category, the Erdogan fundamental solutions [10] for crack problems are employed in the BEM formulation [11][12][13], in which the stress boundary conditions on the crack surfaces are automatically satisfied and the singular behaviour at crack tips can be naturally reflected. Therefore, there is no need to set boundary elements on the crack surfaces, and furthermore, the SIFs of the crack tips can be calculated directly using the specific fundamental solutions of SIFs.
However, for crack problems of complex structures with complex boundary conditions, it is not convenient to model the structures and apply the boundary conditions with BEM, and for multi-crack problems, multi-domain coupling techniques are further required in BEM formulation when the single-crack Erdogan fundamental solutions are used. Therefore, in order to make full use of the unique advantage of BEM in dealing with the singular behaviour at crack tips and the superior flexibility of FEM in modelling complex structures and boundary conditions, several researchers proposed the BEM-FEM coupling methods to solve complex crack problems [14][15][16][17][18].
Actually, the coupling method based on boundary discretization methods and FEM stems from the pioneering work of Zienkiewicz et al. [19] and Brebbia and Georgiou [20], and now have appeared in the areas of flexoelectricity [21], vibroacoustic response [22][23][24], fluid-structure interaction [25] and soil-structure interaction [26], among others. Note that, for crack problems thus far, the BEM formulation for the BEM subdomain containing cracks has been established using the noncrack Kelvin fundamental solutions, and therefore, one also needs to consider the stress boundary conditions on crack surfaces and employ special crack-tip boundary elements to reflect the singular behaviour, which makes it inconvenient to calculate the SIFs at crack tips, even if the coupled BEM-FEM procedure has been adopted.
In this study, a coupled BEM-FEM procedure is developed for the analysis of complex structures in the presence of cracks. The near-crack region is modelled as a cracked superelement, and its stiffness matrix is formulated with the spline fictitious boundary element method (SFBEM) [27][28][29], in which, to automatically satisfy the stress boundary conditions on crack surfaces and naturally capture the crack-tip singularity, the Erdogan fundamental solutions corresponding to an infinite cracked plane are employed instead of the traditional noncrack Kelvin fundamental solutions. For the other regions without cracks, the conventional finite elements are used to model the structure. Then, the proposed stiffness formulation of the cracked superelement can be assembled to the global finite element formulation, from which the nodal displacements of the structure, including those of the superelement, can be obtained following the typical FEM procedure. After obtaining the nodal displacements of the superelement, the SIFs of the crack tips can be readily obtained by a backward analysis with SFBEM focussing on the domain of the superelement. As the SIF fundamental solutions are known in advance as a part of the Erdogan fundamental solutions, the crack-tip SIFs can be calculated directly based on the superposition of the SIF fundamental solutions in the above backward analysis. The accuracy and efficiency of the proposed SFBEM-FEM coupling method are demonstrated by numerical examples involving a rectangular plate with a central crack and a square plate with 100 horizontal cracks. The present approach is further applied to the analysis of a multi-crack steel anchorage box for a hanger of a suspension bridge, indicating the wide applicability of the SFBEM-FEM coupling method to complex crack problems.

Formulation of stiffness matrix for cracked superelement by SFBEM
As the Erdogan fundamental solutions [10] for crack problems are derived by the method of complex variable function and expressed in a more complicated manner, it is more advantageous to employ SFBEM, a nonsingular indirect BEM, for the analysis of the cracked domain, in which singular integrals of the complex fundamental solutions are completely avoided [11][12][13]. In this section, the Erdogan fundamental solutions are first presented, and the stiffness matrix of the cracked superelement is then formulated using SFBEM.

Erdogan fundamental solutions for plane crack problems
Consider an infinite cracked plane subjected to concentrated forces Q and P at an arbitrary source point z 0 = x 0 + iy 0 , as shown in Fig. 1. The complex variable stress functions are taken as follows: where In the above equations, z=x + iy denotes an arbitrary field point; a is the half crack length; and υ is the Poisson's ratio of the material.
The stresses σ 0 x , σ 0 y and τ 0 xy and the displacements u 0 x and u 0 y at any point z=x + iy are derived from the stress functions and can be expressed as 2μ(u 0 where the superscript 0 denotes the responses are caused by Q and P applied at z 0 = x 0 + iy 0 ; the barred quantities are the complex conjugates; the primed quantities represent the derivatives with respect to z; and μ is the shear modulus of the material. It can be seen from Eqs. (1)-(3) that the stresses are singular at the point of concentrated forces, that is, the stresses become infinite as z approaches z 0 . The stresses also become infinite as z approaches ±a, having the required square root singularity at the crack tips.
Accordingly, the SIFs for the crack tips can be further derived as where K 0 I and K 0 II are the SIFs of model-I and model-II, respectively. Let P = 0, Q = 1 or P = 1, Q = 0 in Eqs. (1)- (5). Then, the above solutions become the Erdogan fundamental solutions for plane crack problems, which can be denoted as G K II (z 0 )(k = 1, 2), with k = 1 when P = 0, Q = 1 and k = 2 when P = 1, Q = 0.

Stiffness matrix of cracked superelement
Consider an arbitrary cracked superelement e with a contour L and with a total number of N nodes, as shown in Fig. 2, in which the crack is on the x-axis and its length is taken as 2a and a for an inner crack and an edge crack, respectively. In order to establish the stiffness matrix of the cracked superelement, assume that the nodal displacements U x,n and U y,n (n = 1, 2, . . . , N ) are known, and the nodal forces F x,n and F y,n (n = 1, 2, . . . , N ) are required to be determined from the nodal displacements, which is actually a displacement boundary value problem and can be solved by SFBEM.
Embed the cracked superelement into an infinite plane, and apply unknown fictitious loads X (k) (k = 1, 2) along a fictitious boundary S outside e , whose shape is similar to that of the real boundary L, as also shown in Fig. 2. Then, under the action of the fictitious loads X (k) (k = 1, 2), the displacements and stresses at any point z=x + iy in the infinite plane can be expressed as where z S ∈ S; and G (k) are the Erdogan fundamental solutions of plane crack problems presented in Sect. 2.1.
Because of using the Erdogan fundamental solutions, the stresses shown in Eq. (7) automatically satisfy both the governing differential equations within e and the stress boundary conditions on the crack surface. Therefore, only the displacement boundary conditions along the contour L for the cracked superelement need to be considered.
Substituting Eq. (6) into the displacement boundary conditions along L, one has where z L ∈ L; and u x (z L ) and u y (z L ) are the known displacement functions along L, which can be expressed as where U x,n and U y,n (n = 1, 2, . . . , N ) are the known nodal displacements; and η n (n = 1, 2, . . . , N ) are the displacement interpolation functions. In this study, the boundary displacements are determined by linear interpolation from two adjacent nodal displacements. As the source point z S and the field point z L are located on the fictitious boundary S and the real boundary L, respectively, the integral equations shown in Eq. (8) are nonsingular fictitious boundary integral equations, which need to be solved on a numerical basis. For this purpose, the unknown fictitious loads X (k) (k = 1, 2) are expressed in terms of a set of B-spline functions as follows [11,27]: where N S is the number of fictitious boundary elements; x (k) i (i = −1, 0, . . . , N S + 1) are the unknown spline node parameters; and ϕ i (z S )(i = −1, 0, . . . , N S + 1) are B-spline functions of the third order.
Substituting Eq. (10) into Eq. (8), one obtains the residual functions along boundary L as follows: To eliminate boundary residues R u x (z L ) and R u y (z L ), divide boundary L into N L segments, and let the integrals of the residues along each segment equal zero, namely, where L j is the jth segment on boundary L. Substitution of Eqs. (9) and (11) into Eq. (12) yields where X is the unknown vector consisting of the spline node parameters x (k) i (i = −1, 0, . . . , N S + 1) of the fictitious loads along S; A U is the influence matrix of X depending on the nonsingular integrals of the Erdogan fundamental solution multiplied by the B-spline function; U e = [U x,1 U y,1 U x,2 U y,2 · · · U x,N U y,N ] T is the nodal displacement vector of the cracked superelement; and H is the matrix depending on the integrals of the displacement interpolation functions.
Generally, overdeterminate collocation is employed to achieve a better solution with more boundary segments while keeping the number of fictitious boundary elements at a lower level. Therefore, Eq. (13) needs to be solved on a least-squares basis as follows [27]: where Once the unknown spline node parameters x (k) i (i = −1, 0, . . . , N S +1) are obtained from Eq. (14), the responses of the cracked superelement can be readily calculated from the fictitious loads X (k) (k = 1, 2) shown in Eq. (10). In particular, to obtain the nodal forces F x,n and F y,n (n = 1, 2, . . . , N ) shown in Fig. 2, the surface loads f x and f y along L need to be determined from Eq. (7), which can be expressed as where n x and n y are the cosine functions of normal vector for boundary L. Then, the nodal forces F x,n and F y,n (n = 1, 2, . . . , N ) can be obtained from the surface loads f x and f y using the static equivalent principle, which can be expressed as where η n (n = 1, 2, . . . , N ) are the displacement interpolation functions adopted in Eq. (9). Substituting Eqs. (10) and (16) into Eq. (17) yields where F e = [F x,1 F y,1 F x,2 F y,2 · · · F x,N F y,N ] T is the nodal force vector of the cracked superelement; and A F is the influence matrix of X depending on the nonsingular integrals of the Erdogan fundamental solution multiplied by the B-spline function and the displacement interpolation function. Substituting Eqs. (14) and (15) into Eq. (18), one obtains a relationship between the nodal force vector and the nodal displacement vector for the cracked superelement shown in Fig. 2, which can be written as where K e is the stiffness matrix of the cracked superelement, which can be expressed as It is worth noting that, due to the possible numerical errors induced by SFBEM, the stiffness matrix obtained in Eq. (20) may be slightly asymmetric. The matrix can be made symmetric by minimizing the square of the errors in the nonsymmetric off-diagonal terms [20], which yields  To reflect the degree of asymmetry of matrix K e , define the asymmetry index as follows: where · 2 denotes the L 2 norm of matrix. In general, when the asymmetry index ε is smaller than 5%, the stiffness matrix K e of the cracked superelement will sufficiently approach to a symmetric matrix, indicating the accuracy of SFBEM in the formulation of the superelement.

Analysis of crack problems by SFBEM-FEM coupling method
Now that the stiffness matrix of the cracked superelement has been formulated with SFBEM, it can be further assembled to the global FEM formulation for analysis of the entire structure. After global finite element analysis, the SIFs of the cracked superelement can be finally obtained by a backward analysis using SFBEM.

FEM formulation with cracked superelement
Consider a plane domain with an inner crack or an edge crack, as shown in Fig. 3, in which the near-crack region a is modelled with the cracked superelement e proposed in Sect. 2.2, and the rest region b is modelled using traditional finite elements. As the stiffness matrix of the cracked superelement is established under the local coordinate system x−y, it should be transformed to the form under the global coordinate systemx−ȳ shown in Fig. 3, which can be expressed as whereK e is the stiffness matrix of the cracked superelement under the local coordinate system x−y, which has been obtained in Eqs. (20) and (21) by SFBEM;K e is the stiffness matrix of the cracked superelement under the global coordinate systemx−ȳ; and T is the coordinate transformation matrix, which can be written as where where N is the number of the nodes of the cracked superelement; and β is the angle between the x-axis and thex-axis, as shown in Fig. 3. Suppose that, under the global coordinate systemx−ȳ, the displacement vector of the nodes of the cracked superelement is U a , while the displacement vector of the rest nodes of the finite element mesh is U b . The external force vectors are assumed to be F a and F b corresponding to U a and U b , respectively. Then, following the typical procedure of finite element analysis, the global stiffness equation can be obtained as where and K is the global stiffness matrix, which is expressed as whereK e is the stiffness matrix of the cracked superelement shown in Eq. (23); and K aa , K ab , K ba and K bb are the block stiffness matrices contributed by the traditional finite elements. Note that, for a multi-crack problem, multiple cracked superelements need to be employed to model different near-crack regions, and the global stiffness equation can be formulated as above by merging different superelements to the traditional finite elements.

SIF analysis of cracked superelement
After solving for the global displacement vector U = [U T a U T b ] T from Eq. (26), a backward analysis with SFBEM is required for the SIF analysis of the cracked superelement. For this purpose, the nodal displacement vector U e of the cracked superelement under local coordinate system x−y is first obtained from U a by coordinate transformation as follows: where T is the coordinate transformation matrix shown in Eq. (24). Substituting Eq. (28) into Eq. (14) and considering Eq. (15), the vector of spline node parameters of the fictitious loads can then be obtained as On the other hand, under the action of the fictitious loads X (k) (k = 1, 2), the mode-I and mode-II SIF at the crack tips of the cracked superelement shown in Fig. 2 can be directly obtained as where G (k) K I and G (k) K II (k = 1, 2) are the Erdogan SIF fundamental solutions presented in Sect. 2.1.
Using Eq. (10), the discretization form of Eq. (30) can be derived as where K SIF = [K I K II ] T is the SIF vector of the cracked superelement; X is the vector of spline node parameters of the fictitious loads; and A K is the influence matrix of X depending on the integrals of the Erdogan SIF fundamental solution multiplied by the B-spline function.

Substitution of Eq. (29) into Eq. (31) yields
where Note that, as the fundamental solutions of SIFs are available in the formulation, the SIFs of the crack tips can be directly obtained via Eq. (32) from the nodal displacements

Numerical examples
To validate the accuracy and efficiency of the proposed SFBEM-FEM coupling method, the cracked superelements formulated with SFBEM are incorporated into the commercial FEM software Abaqus with the help of the user-defined elements [30]. Two numerical examples for SIF analysis of a rectangular plate with a central inclined crack and a square plate with 100 horizontal cracks are considered in this section. The finite element meshes for these two examples are generated with the meshing techniques provided by Abaqus.

A rectangular plate with a central inclined crack
Consider a rectangular plate with one side fixed and the opposite side subjected to a uniform tensile load σ = 1 MPa, as shown in Fig. 4. The modulus of elasticity and Poison's ratio of the plate are taken as E = 1000 MPa and υ = 0.2, respectively, and the width and height of the plate are assumed to be 2W = 200 mm and 2H = 250 mm, respectively. An inclined crack AB is located at the centre of the plate, whose length and inclined angle are 2a and β, respectively, as also shown in Fig. 4. The crack problem shown in Fig. 4 is solved with the proposed SFBEM-FEM coupling method, in which the near-crack region is modelled using a 100 mm × 100 mm superelement with the crack located at the centre, and the other region is discretized with traditional Q4 plate elements. In particular, when β = 45 • , the mesh consists of a superelement and 192 Q4 elements, which is shown in Fig. 5. For the cracked superelement, a total of 40 nodes have been adopted, and 40 fictitious boundary elements and 80 boundary segments are employed in calculating the stiffness matrix of the 40-node superelement with the distance between the real boundary and the fictitious boundary being taken as the length of the boundary segment, i.e. d = 5 mm. It has been observed that the asymmetry index ε defined in Eq. (22) is just around 4.62% using the above parameters, and thus the accuracy of the proposed method can be guaranteed regarding the formulation of the superelement. For comparison, the crack problem is also analysed using the traditional FEM with Abaqus, in which, to obtain the reference solution, the near-crack region is discretized into a sufficiently fine mesh leading to a total number of 1748 Q4 elements for the entire problem (a/W = 0.3), as shown in Fig. 6. Comparing the finite element meshes shown in Figs. 5 and 6, it can be seen that the number of DOFs involved in the coupling method is much smaller than that of the traditional FEM due to the use of the cracked superelement. Furthermore, for a given value of β, the coupling mesh shown in Fig. 5 remains the same for different crack lengths, while for the domain shown in Fig. 6, it needs to be remeshed when the crack length changes.  When β = 0 • and β = 45 • , the SIFs K I and K II of crack tips A and B corresponding to different crack lengths are presented in Tables 1 and 2, respectively, from which it can be seen that the results of the SFBEM-FEM coupling method agree well with those of the traditional FEM with relative deviations being less than 1.5%, indicating the high accuracy of the present method. To investigate the variation of K I and K II with the inclined angle β under different ratios of the crack length to plate width, a/W , the results are depicted in Figs. 7 and 8 for crack tips A and B, respectively. It can be observed from the above figures that K I decreases with the increase of β for both crack tips. While for K II , it will first increase with the increase of β until β = 45 • ,   Tables 3 and 4, from which it can be observed that the distance between the fictitious boundary and the real boundary can be changed in a reasonable range with little influence on the results, indicating the numerical stability of the present method. As a general rule, the distance between the fictitious boundary and the real boundary in SFBEM should be taken as the same order of magnitude of the length of the real boundary segments. This conclusion has also been drawn in the previous study reported in [11,27]. Detailed investigation on the optimal location of the fictitious boundary can be found in [31].

A square plate with 100 horizontal cracks
Consider a square plate containing 10 × 10 = 100 equally spaced horizontal cracks, as shown in Fig. 9. The length of each crack is taken to be 2a = 40 mm. The plate is fixed at the bottom side and is subjected to a uniform tensile load σ = 1 MPa at the top side. The modulus of elasticity and Poison's ratio of the plate are assumed to be E = 200 GPa and υ = 0.25, respectively, and the width and height of the plate are assumed to be 2W = 3000 mm and 2H = 3000 mm, respectively. The crack problem is solved with the proposed SFBEM-FEM coupling method, in which the plate domain is discretized with 10 × 10 = 100 cracked superelements, as shown in Fig. 9. The width and height of each superelement are 2w = 300 mm and 2h = 300 mm, respectively, and the crack is located at the centre of the superelement. To obtain the SIF results with sufficient accuracy, a total of 40 nodes are adopted for each superelement, and 40 fictitious boundary elements and 80 boundary segments are employed in calculating the stiffness matrix of the superelement with the distance between the real boundary and the fictitious boundary being taken as the length of the boundary segment, i.e. d = 15 mm. The asymmetry index ε defined in Eq. (22) has been found to be just 3.73% with the above parameters, indicating the accuracy of the proposed method regarding the formulation of the superelement. Note that, for this particular problem, as the 10 × 10 = 100 cracked superelements have covered the whole domain of the plate, no traditional finite elements are required in the present approach. For comparison, the crack problem is also analysed using the traditional FEM with Abaqus, in which a total number of 140539 Q4 elements need to be used for obtaining the reference solution of the problem, as shown in Fig. 10. Obviously, the number of DOFs involved in the coupling method with the mesh shown in Fig. 9 is greatly reduced compared with that of the traditional FEM using the mesh shown in Fig. 10, and it is much easier to prepare for the mesh shown in Fig. 9 than that shown in Fig. 10.
To investigate the variation of the SIFs of the cracks located at the same horizontal section of the plate, sections A-A and B-B are taken into consideration with the cracks being numbered from 1 to 10, as shown in Fig. 9. Consider the left tips of the above cracks, and the values of K I for different cracks at sections A-A and B-B are shown in Figs. 11 and 12, respectively, from which it can be seen that the results obtained with the SFBEM-FEM coupling method and the traditional FEM are in good agreement with the relative deviations being less than 0.5%, indicating the high accuracy of the proposed method. It can be also observed from the above figures that, at section A-A, the values of K I range from 7.7905 to 8.0734 MPa · mm 1/2 with those of the cracks on both sides of the plate being a bit smaller, while at section B-B, the values of K I range from 7.6157 to 8.9169 MPa · mm 1/2 with those of the cracks on both sides of the plate being significantly larger. Note that, for this particular example, the computing CPU time required by the coupling method is only 9.1% of that required by the FEM alone for solving the final system of equations.

Engineering application
Crack propagation of the steel anchorage box for a hanger of a suspension bridge often occurs under the action of vehicle loads. The SIFs of cracks are important parameters for fatigue analysis of the steel anchorage box, and thus it is necessary to conduct the SIF analysis of the structure. A steel anchorage box is shown in Fig. 13. The back plate of the structure is assumed to be fixed at the upper and lower edge with two side edges being free. The bearing plate of the structure is subjected to a uniform pressure with an amplitude of σ = 1.4057 MPa, which is transformed from the tension amplitude of a hanger rod under the action of vehicle loads. The modulus of elasticity and Poison's ratio of the steel anchorage box are assumed to be E = 210 GPa and υ = 0.3, respectively. Suppose that, after field inspection, three cracks are detected at the back plate, including an inner crack and two edge cracks, as also shown in Fig 13. The lengths of the inner crack and the two edge cracks are assumed to be 2a 1 = 100 mm and a 2 = 100 mm, respectively. In this section, the cracked superelements formulated with SFBEM are incorporated into the commercial FEM software Abaqus with the user-defined elements, and the finite element mesh is generated with the meshing techniques provided by Abaqus.
The crack problem shown in Fig. 13 is first solved with the traditional FEM using Abaqus. To obtain the reference solution, a total of 6404 S4 shell elements are adopted for the analysis, and the finite element mesh is shown in Fig. 14. To validate the applicability of the present approach to engineering practice, the problem is also analysed with the SFBEM-FEM coupling method. The near-crack regions are modelled using one 400 mm × 400 mm centre-cracked superelement and two 300 mm × 300 mm edge-cracked superelements. To keep the superior flexibility of FEM in modelling the complex structure, the other region is discretized with 2601 traditional S4 elements. The coupling mesh is shown in Fig. 15. For the centre-cracked superelement, a total of 40 nodes have been adopted, and 40 fictitious boundary elements and 80 boundary segments are employed in calculating the stiffness matrix with the distance between the real boundary and the fictitious boundary being taken as d = 20 mm. For the edge-cracked superelement, a total of 41 nodes have been employed, and 60 fictitious boundary elements and 120 boundary segments are adopted in the formulation of the stiffness matrix, in which the distance between the real boundary and the fictitious boundary is taken as d = 6 mm. It can be found that, using the above parameters, the values of the asymmetry index ε defined in Eq. (22) are only 3.88% and 4.57% for the centre-cracked and the edge-cracked superelement, respectively, indicating the accuracy of SFBEM in calculating the stiffness matrices of the superelements.
The results of K I and K II for the crack tips A, B, C and D obtained by the coupling method and the traditional FEM are presented in Table 5. It can be seen from Table 5 that the results of the coupling method agree well with those of the traditional FEM. The relative deviations of the coupling method are less than 1.3%, indicating the high accuracy of the present approach. Note that, for this example, the elapsed CPU time of the coupling method is only 37.5% of that required by the FEM alone for solving the final system of equations.

Conclusions
To make full use of the unique advantage of SFBEM in dealing with the singular behaviour at crack tips while keeping the wide applicability of FEM in modelling complex structures and boundary conditions, this paper presents a SFBEM-FEM coupling method for the SIF analysis of complex crack problems. The stiffness matrix of the cracked superelement is first formulated by SFBEM based on the Erdogan fundamental solutions, in which the stress boundary conditions on the crack surfaces are automatically satisfied and the singular behaviour at crack tips can be naturally captured. Then, the proposed cracked superelements are incorporated into a finite element mesh to simulate the behaviour of different crack zones, and the global finite element equations are established and solved following the typical FEM procedure. Finally, the SIFs of the crack tips can be directly obtained by a backward analysis with SFBEM by superposition of the analytical SIF fundamental solutions, which are included as a part of the Erdogan fundamental solutions. Two numerical examples as well as an engineering application are presented to show the efficacy of the present approach by merging of the cracked superelements to the commercial FEM software. However, it should be pointed out that the present study is only limited to two-dimensional crack analysis. The coupling method for three-dimensional crack problems will be further considered in future study.