Boundary element simulation of fatigue crack growth in multi-site damage

This paper presents an eﬃcient and automatic scheme for modelling the growth of multiple cracks through a two-dimensional domain under fatigue loading based on linear elastic fracture mechanics. The dual boundary element method is applied to perform an analysis of the cracked domain and the J -integral technique is used to compute stress intensity factors. Incremental crack propagation directions are evaluated using the maximum principal stress criterion and a combined predictor-corrector algorithm implemented for propagation angle and increment length. Criteria are presented to control the mesh used on the slower growing cracks in the domain, improving computational eﬃciency and accuracy by the use of virtual crack tips to avoid the need for severe mesh grading. Results are presented for several geometries with multi-site damage, and sensitivity to incremental crack length is investigated. The scheme demonstrates considerable advantages over the ﬁ-nite element method for this application due to simplicity of meshing, and over other boundary element formulations for modelling domains with large ranges of crack growth rates.


Introduction
Fatigue is a common mode of failure in materials subjected to cyclical loading.Here, cracks propagating through a structure reach a critical length, upon which a sudden and often catastrophic fracture failure occurs.The growth of cracks can be difficult to detect and monitor, necessitating the requirement for methods to simulate the behaviour of such cracks.
Modern structures and components can contain many thousands of subcritical cracks and therefore an important design consideration is the extent of crack growth that can be sustained by the structure safely within its designed life.The damage tolerance approach to the fatigue assessment of a structure requires engineers to monitor cracks and also be able to calculate the remaining life.In order to make accurate calculations, detailed crack growth calculations must be performed.This is often done with the use of numerical elasticity calculations (e.g. with the finite element method (FEM)) in combination with a crack growth law.
A further consideration is multi-site damage, where a component contains multiple cracks of different sizes that can propagate at different rates and the interaction of cracks becomes significant.This phenomenon was brought into focus in 1988 when Aloha Airlines flight 243 experienced an explosive decompression following a structural failure of the fuselage in flight.The National Transportation Safety Board [1] determined in its accident report that the cause of the damage was failure of a fuselage lap joint from multi-site fatigue cracking.The consideration of multi-site damage has since become a critical consideration in aircraft design and maintenance; however, numerical simulation of multi-site fatigue crack growth remains challenging.
The increasingly complex geometries and crack interactions that are present in modern engineering structures require the development and use of numerical methods to simulate the propagation of cracks and to compute the resultant effect on stress fields.Linear Elastic Fracture Mechanics (LEFM) has long been used in damage tolerance assessments for cracked bodies.Here it is assumed that the crack tip plastic zone is small in comparison with the crack length.A complication arising in LEFM, and a particularly important one in performing numerical simulations, is that a stress singularity is present at the crack tip, so that values of local stress components become of limited use in assessing the stability of the crack and its propagation properties.Instead, stress intensity factors, K I and K II , are used to provide a convenient measure of crack stability, and also describe components, σ ij , of the stress tensor in the vicinity of the crack tip.For a pure mode I crack, for example, the stress field around the crack tip is given, in the usual polar coordinate system (r, θ) centered on the crack tip, by where K I is the mode I stress intensity factor and f ij is a dimensionless func-tion of θ, the angular coordinate measured from the crack axis (see Fig. 1).As r → 0, a singularity of order r −1/2 develops and the stress σ ij → ∞.The stress intensity factor is thus a scaling factor which indicates the severity of the singular crack tip stress field.In numerical simulations using any method based on LEFM, cracks are modelled in incremental analyses whereby a crack is incremented by a precalculated length and an analysis of the stress field in the structure is conducted and new stress intensity factors found.The crack propagation angle for each increment can be determined as a function of the stress intensity factors.
The FEM has been at the forefront of numerical simulation for decades and has been successfully used in fracture mechanics applications considering crack propagation [2].However, when applied to iterative crack growth schemes, the FEM requires extremely fine meshing around crack tips to resolve the high stress gradients, and each increment in the analysis requires a re-meshing of the domain.These problems were addressed by Moës et al. [3] who developed the extended finite element method (XFEM), which adds local enrichment to model the crack independently of the finite element mesh.
The scaled boundary finite element method [4] is an alternative approach which has benefits arising from (i) a reduced dimensionality of modelling, and (ii) a semi-analytical approach that can capture the leading order terms in the expansion (1) and yield stress intensity factors directly.However, the approach can become cumbersome for all but the simplest geometries.
The boundary element method (BEM) has become popular for modelling of fractures, and particularly for crack growth.This removes the need for remeshing, as only the boundary (including crack surfaces) are meshed, allowing subsequent increments to be considered simply by adding incremental elements to the existing mesh.The method is also well known for providing accurate boundary solutions to discontinuous fields.Using the classical form of the BEM, however, it is not possible to achieve a solvable system of equations in a single region, since collocation at coincident nodes on opposing crack surfaces causes the number of independent equations in the system to become insufficient.The Dual Boundary Element Method (DBEM) (see Hong & Chen [5], Portela et al. [6], and Chen & Hong [7]) addresses this by defining separate equations on collocated surfaces of the crack, one a function of displacement and the other of traction, giving a non-singular system of equations which can be solved.The method has been applied to both single and multiple crack problems [6,8,9,10,11].
This paper presents an algorithm for the application of the DBEM to the incremental analysis of multi-site damage crack growth problems.In the evaluation of crack tip stress intensity factors, the J-integral is used, and for calculating crack paths the predictor-corrector technique proposed by Portela et al. [6] has been incorporated into a multiple crack algorithm.Fatigue analysis based on Paris law is adopted for simplicity, although other crack growth laws could easily be substituted.Further correction for increment length is implemented based on Salgado et al. [9], and new growth criteria developed to control the extensions of cracks where a range of growth rates are present.Examples for multi-site damage fatigue problems are presented including single crack examples for validation and a series of multi-site damage applications with relevance to modern engineering structures.
The paper is organised in the following fashion.In section 2, we present the DBEM and its application for fatigue crack growth.In section 3, we extend the discussion to multi-site damage, presenting a new algorithm which includes the possibility of virtual crack extensions, and the approach is formalised in section 4. In section 5 we present validating examples, and discuss the performance of the algorithm.We close with some concluding remarks in section 6.

The dual BEM for fatigue crack growth
We start by defining our two-dimensional domain Ω ⊂ R 2 , having Lipschitz boundary Γ := ∂Ω.The static equilibrium state of an elastic body described by the domain can be expressed in the following way, where σ is the Cauchy stress tensor and b represents the body forces per unit volume.The stress tensor is related to components of the strain tensor, ε, through the constitutive relations and the strains are displacement derivatives according to the definition In the above λ is the Lamé constant, given by µ is the shear modulus, given by δ is the Kronecker delta, E is Young's modulus, ν denotes Poisson's ratio and u i are the displacements.In addition, the strains have to satisfy the compatibility equations We consider a boundary value problem in which the above differential equations are solved subject to some Dirichlet and Neumann boundary conditions, respectively where t i is a traction component and the overbars denote prescribed values.The Dirichlet and Neumann boundaries, Γ u and Γ t , form the entire boundary, so Γ = Γ u ∪ Γ t .Betti's reciprocal work theorem is applied to form the Boundary Integral Equation (BIE) for displacements at a point x ′ ∈ Ω\Γ, in terms of tractions and displacements at all points on Γ. Taking for simplicity b i = 0, i = 1, 2, this yields where T ij and U ij represent the Kelvin fundamental solutions for tractions and displacements respectively.Using the usual radial coordinate r := |x − x ′ |, the fundamental solutions become singular as r → 0; a strong 1/r singularity develops in T ij and a weak ln(1/r) singularity in U ij .Taking the limit of (10) as x ′ → Γ, gives where − ∫ indicates an integral taken in the Cauchy principal value sense, and the free coefficient c ij (x ′ ) arises as a result of the limiting procedure; it takes the value δ ij /2 for a smooth boundary at x ′ .
Repeated application of (11), taking x ′ at all nodes in turn, allows us to form a linear system which, upon prescription of a sufficient number of boundary conditions, can be solved to yield the nodal displacements and tractions.However in a single region analysis of a domain with cracks there are problems with the coincidence of the two crack surfaces giving rise to a singular system as a result of applying (11) at coincident nodes on the upper and lower crack surfaces, thereby yielding duplicate equations.
Several variations of the BEM have been proposed to address this problem, e.g., use of a fundamental solution that does not require the crack surfaces to be included in the BEM mesh [12] and a domain decomposition method [13].These suffer from lack of generality and difficulty of automating incremental crack growth.These problems were most effectively addressed in the development of the Dual Boundary Element Method (DBEM) [5,11,7], in which ( 11) is used when collocating on one crack surface and a hyper-singular traction integral equation on the opposing surface.This latter equation is derived through differentiation of (11) and application of Hooke's law to give an expression for traction components, t i , as where S ijk and D ijk contain derivatives of T ij and U ij and exhibit a 1/r 2 hyper-singularity and strong 1/r singularity respectively, and = ∫ indicates the Hadamard principal value integral.For piecewise flat cracks, it is possible to evaluate analytically the singular integrals required when integrating over the element containing x ′ [11].The DBEM linear system is then formed by repeated application of both ( 11) and ( 12), at a sufficient number of collocation points x ′ , as follows: where Γ = Γ T ∪ Γ D and the boundaries Γ T , Γ D are defined as in which Γ ext contains all portions of the domain boundary that do not describe a crack surface, and Γ + j and Γ − j are the two opposing surfaces of the j th crack of a total of M cracks.
A further advantage of the DBEM is that in an incremental analysis, the new elements added for the crack extensions form new columns and rows in the system matrix.The majority of the original system matrix can therefore be reused and techniques such as LU decomposition or Gauss elimination can be used to reduce this part of the matrix to optimise the re-analysis after each crack growth increment.
For the work in this paper the J-integral method [14] was used to evaluate the stress intensity factors.In the absence of body forces this is given by where S is an arbitrary counterclockwise contour surrounding the crack tip, W is the strain energy density, and the direction i = 1 relates to the direction of crack propagation, x 1 , as shown in Fig. 1.A relationship exists between J and the modal stress intensity factors, given by where E ′ is the elasticity modulus equal to E in plane stress and to E/(1−ν 2 ) in plane strain.In a mixed-mode two-dimensional problem, it is necessary to decouple the combination of modes in the J-integral, in order to compute K I and K II separately.The integral, J, in ( 17) can be given as the sum of two integrals using the method of Aliabadi [15], which involves taking a distribution of integration points on S that is symmetric about x 2 = 0, e.g. the points P and P ′ , shown in Fig. 1.This allows the displacements and stresses to be separated into their symmetric and antisymmetric parts, and the J-integral then to be decoupled, expressed in terms of the modal components J I , J II , giving The modal components of the J-integral can then be related to the modal stress intensity factors, giving Within the present analysis, discontinuous boundary elements are used throughout, satisfying the Hölder continuity requirement embedded in the hypersingular integral term in (12), and (as suggested by numerical tests) the J-integral path is taken as the circular contour defined by the radius of the fourth node from the crack tip.Both of these features are evident in Fig. 1.
Several methods of determining crack propagation directions in mixedmode problems have been presented, commonly those based on the orientation of the maximum principal stress (MPS) and also on maximum strain energy density [16].The latter finds more use in three-dimensional applications [17]; in the present work the MPS criterion is applied as proposed by Erdogan and Sih [18].Near the crack tip the maximum principal stress is taken as the tangential stress σ θθ , the singular part of which is given by The angle θ = θ t at which the crack will propagate is that maximising σ θθ , giving an expression which can be solved for the crack propagation angle, θ t , where only the negative result of the ± expression is of practical use in this application.
The MPS criterion is a continuous criterion; the crack propagation angle found using (23) is evaluated locally at the crack tip.In a discretised numerical scheme such as that proposed in the present work, this criterion does not take into account the variation in the stress field as the crack propagates; the incremental extension direction will always be defined in the same direction regardless of the crack extension length used.As a result the crack path calculated using (23) cannot be guaranteed to be unique.
A form of correction is required to consider the stress conditions at the extended crack tip, and to re-evaluate the crack propagation angle until the true crack path is achieved.Portela et al. [6] proposed an iterative, predictorcorrector scheme that achieves this by introducing a correction angle β.This approach is repeated here for completeness.Considering a crack increment, γ, the propagation direction θ t(γ) is computed from (23).After extending the crack by da along this angle, the actual crack path has deviated, therefore the crack propagation angle for the next increment is θ t(γ+1) .Fig. 2 shows that a correction angle β can be expressed as a geometric relationship, β = θ t(γ+1) /2.In an incremental analysis the following procedure is implemented at each crack tip in the domain: 1.For the first iteration of the predictor corrector algorithm, the crack propagation angle, θ t , is evaluated using the MPS criterion (23), 2. The crack is extended by da to the new crack tip, p i , 3. A full DBEM analysis of the domain is conducted and new crack tip stress intensity factors calculated, 4. A new crack propagation angle, θ i t(γ+1) , is computed from (23) with the new SIFs found in step 3, 5.The correction angle, β i , is calculated as β i = θ i t(γ+1) /2, 6.The crack increment is corrected as θ i+1 t(γ) = θ i t(γ) + β i , the previous extension in step 2, is replaced with an extension in the direction of the corrected angle, 7. The algorithm is repeated from step 2 until Notice that β → 0 as the crack increment da → 0; therefore in the limit the crack path tends to the direction of the tangent of the continuous crack path.
Cracks subjected to fatigue loading, characterised by a load with varying amplitude with time, will grow under certain conditions.Simulating the growth of cracks under fatigue loading allows fatigue life and residual strength to be determined, which can be used in damage tolerance assessments.In a multiply cracked domain, fatigue analysis must be conducted as part of an iterative algorithm; for each increment, all cracks in the domain will be subjected to a given number of load cycles resulting in different increment lengths for each crack.
The growth rate of a crack under fatigue loading, da/dN , is usually given as a function of the range of effective stress intensity factors, ∆K eff .Many formulae for this relationship exist but Paris Law [19] is used here as it is the most simple formulation, given by where a is the crack length, N the number of load cycles, and C and m are material constants.Tanaka [20] showed experimentally that for mixed-mode problems ∆K eff may be given by where ∆K = K max − K min and K min and K max are the minimum and maximum values of the stress intensity factor within a load cycle.In order to consider problems with non-zero stress ratio, R, ∆K can be expressed as follows: where, R = K min /K max .

Algorithm for simulation of multi-site damage
In a domain with multiple cracks it can be assumed that there will generally be a range of stress intensity factors at the different crack tips in the domain, and as a result a range of growth rates and associated crack increment lengths.In the present work the following procedure was adopted to determine initial growth conditions for all cracks in the domain.The number of load cycles required to increment the fastest growing crack tip by a given length is calculated by rearranging (24) in terms of dN and integrating, where da is a reference increment length, taken to be two times the length of the crack tip element on the fastest growing crack (i.e. that with the maximum value of ∆K eff in the domain), and a i is the current crack length.This choice of da is motivated by the desire to start and end the J-integral contour at x 2 = 0 (Fig. 1) at the fourth node from the crack tip.The integration assumes that the stress intensity factors are constant along the crack as it extends from a i to a i + da.
The crack growth increments for all remaining crack tips can then be calculated by integrating the inverse of ( 27), where N i is the current number of load cycles that have been simulated in the analysis, and dN is the number of load cycles calculated in (27); again the integration is performed assuming constant stress intensity factors.The stress intensity factors at the extended crack tip are not known prior to the extension.Once the crack has been extended and a re-analysis performed to correct the propagation angle, the increment length can be corrected as well to take into account the variation in crack tip stress intensity factors.The predictor-corrector scheme of Salgado and Aliabadi [9] has been implemented in the present work.
The procedure is adapted in the present work so that corrections for propagation angle and increment length are conducted in the same algorithm.The following steps are computed alongside the corresponding points in the propagation angle correction algorithm: 1.For the first iteration the number of load cycles for the increment and initial increment length, da, are calculated for each crack tip using ( 27) and (28) respectively, 2. All cracks are extended by da to the new crack tip along the propagation angle θ t calculated in the predictor corrector algorithm using (23), 3. A full DBEM analysis of the domain is conducted and new crack tip stress intensity factors calculated, 4. The number of load cycles for the increment are re-calculated considering variation in stress intensity factors by applying the trapezium rule to (27), 5. Crack increment lengths are re-calculated considering variation in stress intensity factors by applying the trapezium rule to (28), 6.The algorithm is repeated from step 2 until |da i − da i+1 | < ε a , where ε a is a pre-defined tolerance.
The scheme presented is simple and the combined algorithm for correction of propagation angle and increment length is computationally efficient.

Mesh grading and virtual crack extensions
In a multiple crack analysis the size of the crack increment da may become problematic for the more slowly growing cracks, i.e. those for which ∆K eff is considerably smaller than max∆K eff , since the new elements used to mesh the crack increment will become small.This is undesirable, not only because it is computationally inefficient to simulate crack growth using very small increments, but more importantly it may prevent a suitably accurate DBEM solution being achieved.Most implementations of the boundary element method require adjacent elements to be of similar length in order to compute accurately the near-singular integrals in ( 11) and ( 12).This is especially pertinent when using discontinuous elements.It is therefore necessary to ensure that the ratio between adjacent element lengths can be controlled.
A set of criteria are therefore proposed to enforce minimum adjacent element lengths when extending cracks to maintain suitable ratios of adjacent element lengths and ensure accuracy of the DBEM solution: • A given ratio between the lengths of adjacent elements must not be exceeded; numerical tests suggest that a target ratio of 1.5 is suitable.
• A crack extension must exceed a minimum increment length determined by the original geometry of the crack and the desired mesh density.
If either of these criteria is not met, then the crack is not extended for that increment of the crack growth scheme.That is, the boundary mesh will not be updated to include the new crack tip; the updated position of this crack tip is only recorded as a virtual crack tip.On subsequent iterations the analysis is performed at the meshed crack tip, incurring some approximation, however new crack tip positions are calculated by incrementing by da from the virtual crack tip with propagation angle θ t also calculated using stress intensity factors from the meshed crack tip.The distance from this tip to the meshed crack tip is then checked against the criteria again.The process of virtual increments is repeated until both criteria are met, the crack is then extended to the virtual crack tip and the mesh updated (i.e. the full crack is meshed and the tip is no longer virtual).The crack increment is then corrected for propagation angle and increment length.
This procedure enables cracks with slower growth rates (da/dN ) to be efficiently and accurately analysed using the DBEM.This is a significant advantage in terms of computational efficiency.

Incremental Crack Growth Scheme
The full structure of the incremental crack growth algorithm presented in this paper can be formalised in 13 steps as follows: 1.A DBEM analysis of the cracked domain is conducted and stress intensity factors are evaluated at each crack tip using the J-integral technique, 2. Effective stress intensity factors are calculated for all crack tips in the domain using (25), 3. The crack tip with the maximum range of effective stress intensity factor is identified, 4. The number of load cycles, dN , to increment this crack tip by an increment equal to two times the previous element is calculated using (27), 5.The extension length da is then calculated for all other crack tips in the domain using the number of load cycles, dN , evaluated in step 4 and (28), 6.The crack propagation direction, θ t is calculated for each crack tip using the MPS criterion (23), 7. The position of new crack tips is calculated using da and θ t , incrementing from the previous crack tip (or virtual crack tip as appropriate), 8.For each crack tip the distance, a ext , is evaluated from the meshed crack tip to the new crack tip (or virtual crack tip as appropriate): (a) If a ext does not satisfy either of the crack extension extension criteria, the crack increment is virtual.No further action for this crack tip on this increment is required and the crack tip is recorded as virtual (b) If a ext satisfies both extension criteria, the increment is applied to the crack to be meshed 9.The boundary element mesh is updated for all new crack increments 10.A new DBEM analysis is performed, and stress intensity factors are calculated at crack tips using the J-integral technique, 11.New effective stress intensity factors are calculated at all crack tips using (25), 12.The crack propagation angle, θ t , and crack increment length, da, are corrected considering the new stress intensity factors at each crack tip using the predictor corrector algorithms, 13.The new propagation direction and length for each crack tip are checked for convergence.For any unconverged crack tips, the algorithm is repeated from step 9, until all crack tips in the domain are converged.
The scheme presented is computationally efficient, and is especially appropriate for domains with multi-site damage and large ranges of crack growth rates.The criteria proposed to control the meshing and extension of cracks allow the simulation to be optimised for a desired level of solution accuracy.

Multi-Site Damage Crack Extension Applications
Three applications of the incremental crack-extension scheme are discussed.The first considers a singly cracked plate with three holes which is based on an experiment conducted by Bittencourt et al. [2] and is compared to a numerical BEM result by Leonel & Venturini [8], the second considers the case of offset cracks approaching across a plate under a uniaxial fatigue load, and the final example considers the effect of multiple cracks emanating from perforations in a panel.In all examples, the given Paris law constants, C and m are suitable when crack growth rates are considered in mm/cycle and SIFs in GPa √ m.

Cracked panel with three holes
To validate the crack growth scheme presented in this paper, an example presented by Leonel & Venturini [8] is demonstrated; a rectangular plate with three holes and a single edge crack with geometry as shown in Fig. 3 was tested with the incremental crack growth scheme presented in this paper.A cyclic point load of 4.448 kN was applied to the centre of the top edge of the plate.Young's modulus of E = 2.068 GPa and a Poisson's ratio of ν = 0.3 were used with Paris' law constants of C = 7.69 × 10 −12 and m = 3.0.The crack was modelled using 16 discontinuous, quadratic elements, and a further 36 elements (of the same type) were added over the nine increments considered.Fig. 4 shows the result of the nine increment analysis with the discussed DBEM model compared to the numerical result presented by Leonel & Venturini [8].
The DBEM result closely matches the alternative numerical result although with slightly different deviation around the lower hole; this result appears closer to the experimental result presented by Bittencourt et al. [2], though we note the likelihood of different discretistations used compared to [8].The method used by Leonel & Venturini uses the displacement correlation technique to calculate stress intensity factors; more investigation would be required to determine which method gives the best result in difficult geometries such as around the lower hole in this example.However, the comparison demonstrates, sufficiently to progress to the other examples including multi-site damage, that the proposed algorithm gives a satisfactory simulation of a single crack path under a cyclic fatigue load.

Panel with multiple approaching edge cracks
Consider a rectangular plate with two cracks of length 10 mm on opposing sides of the plate with a vertical offset of 5 mm, subjected to cyclic tension of 100 MPa to the top edge of the plate and fixed at the bottom edge.The geometry is shown in Fig. 5(a); a stress ratio of R = 0, Young's modulus of E = 30 GPa, and a Poisson's ratio of ν = 0.3 were used with Paris law constants of C = 10 −12 and m = 3.0.The cracks were modelled using a total of 52 discontinuous, quadratic elements, and a further 44 elements (of the same type) were added over the eleven increments considered.Fig. 5(b) shows the crack paths after eleven increments of the DBEM crack growth algorithm.Notice that the problem is asymmetric, and this results in somewhat greater stress intensity factors for crack 2, thus the growth is faster for this crack.As the fractures develop the rate of growth of crack 2 increases over that of crack 1.As a result, the increment length required to model the incremental extension of the slower crack 1 decreases.This analysis demonstrates the effectiveness of the present scheme in this situation; three virtual increments were made on crack 1, allowing the full path of crack growth to be modelled without the use of small elements that would have introduced numerical difficulties.
The cracks initially propagate without significant deviation, but as the mode II stress intensity factor increases the cracks begin to deviate away from each other before deviating back towards the opposing crack as is characteristic of this type of geometry.This example also demonstrates the requirement for the consideration of crack coalescence; further incrementation would result in crack 2 coalescing with crack 1.This has not been considered in the present work.
The predictor-corrector algorithms are next validated for path independence from increment length.Two values of the reference increment length, da = 1.71 and 2.64 mm, were tested to compare the effect of different increment lengths on the predicted crack propagation.In Fig. 6 we plot the crack paths for the two values of da.For the same geometry the crack paths correlate well for the different increment lengths, though small differences present around the middle of the plate, as a result of the coarseness in the meshing with the larger increment length, demonstrate the discretisation errors that can accrue.
A comparison of the structural life prediction is shown in Fig. 7 for the different increment lengths da = 1.71 and 2.64 mm; the structural life is estimated at 18000 cycles.The comparison shows that there is little difference in the number of cycles subjected to the structure, especially before the final acceleration to failure, a stage during which it is unlikely that the Paris law would remain valid.
Fig. 8 shows the stress intensity factors plotted against crack length for modes I and II for each crack.K II remains small for the first 10 mm of crack growth, resulting in the relatively straight crack paths for the initial growth shown in Fig. 5(b).As the cracks approach each other there is a rapid increase in K I at both crack tips, however the rate of growth of crack 2 does become visibly larger than that for crack 1 from 25 mm total crack length.

Multiply perforated panel with multiple cracks
Consider a perforated panel fixed at the left hand edge and subjected to a cyclic load to the right hand edge.The panel has three circular perforations of diameter 10 mm along the centreline, with cracks of length 3 mm emanating from the top and bottom of each perforation; the geometry and definitions of the crack numbering are shown in Fig. 9. Young's modulus of E = 30 GPa, a stress ratio of R = 0, and a Poisson's ratio of ν = 0.3 were used, with Paris law constants of C = 10 −13 and m = 2.6.Cyclic loads were applied to the right edge of the panel of 9 MPa in the x direction, and of 1 MPa in the y direction.The cracks were modelled using a total of 48 discontinuous, quadratic elements, and a further 188 elements (of the same type) were added over the incremental crack growth analysis.Fig. 10 shows an enlarged view of the crack paths after nine increments of the DBEM algorithm presented in this paper.It is clear that the bending stresses induced by the applied vertical load have resulted in larger crack growth rates for the cracks towards the bottom of the plate, particularly cracks 1 and 2. Growth of cracks 5 and 6 remains small and this shows the effectiveness of the crack extension criteria and use of virtual crack tips; over the nine increments of the algorithm made in the analysis only one new extension was made on crack 6, the remaining increments being accommodated through virtual extensions.
Fig. 11 shows the variation in stress intensity factors with crack length for each iteration of the DBEM algorithm, including those that resulted in a virtual increment (K I only shown for clarity).Some care should be taken in interpreting the figure, since the cracks all reach any prescribed length at different times.However, plotting the results in this form is a useful way of describing both the relationship K I vs. a for each crack and also the relationship between the different crack lengths in each iteration.The figure shows that, after an initial phase in which crack 1 exhibits the highest K I , cracks 2 and 3 become the most significant after approximately 2 mm of crack growth as the situation becomes more severe here (with two cracks approaching each other) than for crack 1 (which does not have another crack growing towards it).Cracks 2 and 3 also experience rapid changes in K II for the last increment as the cracks approach each other; a similar behaviour can be seen in the second example (Fig. 8).
This example demonstrates the ability of the algorithm to calculate increment lengths in cases in which the fastest growing crack changes during the analysis.The figure also shows how the virtual crack extension results in larger crack increments being used to model slower growing cracks; cracks 4, 5 and 6 all show some virtual increments.The result of the proposed extension criteria is that small increments are eliminated with no adverse effect on the results obtained.

Conclusions
This paper presents a computational scheme based on the Dual Boundary Element Method to perform incremental crack growth analyses of cracked domains with multi-site damage.The J-integral technique was used to compute stress intensity factors from which crack propagation angles were determined using the maximum principal stress criterion.A Paris law crack growth model was used to compute crack increment lengths based on a the number of fatigue load cycles required to propagate the fastest growing crack.More sophisticated crack growth laws could easily be substituted.Predictor-corrector algorithms for propagation angle and increment length were implemented in a combined fashion to ensure accurate modelling of crack paths.Criteria have been proposed to control the incrementation of slower growing cracks allowing multiply fractured domains to be simulated whilst maintaining the intrinsic benefits of the BEM in terms of mesh simplicity, numerical accuracy and computational efficiency.An essential feature of the new approach is the use of virtual crack tips to overcome difficulties associated with meshing slowly growing cracks.
Three examples were tested, demonstrating the performance of the model in simulating interacting mixed-mode crack growth, including cases in which a range of crack growth rates were present.The sensitivity of the results to the increment length, which is based on the initial mesh, was investigated.The two incremental lengths tested produced well correlated crack paths, with some small differences due to discretisation error in the coarser model.These differences made little difference to the fatigue life prediction.
Not included in the present work is consideration of crack coalescence, though since techniques have been developed for multiply cracked models (e.g.[8]) the extension to merging cracks is quite possible.This will form a natural extension to the present work, which focusses on the accommodation of different crack growth rates within the same model.

Figure 1 :
Figure 1: Circular J-integral contour path and coordinate reference system relative to crack tip.

Figure 2 :
Figure2: Geometric relationship of crack path correction angle to crack propagation angles calculated using the MPS criterion[6].

Figure 5 :
Figure 5: Panel with offset approaching edge cracks.

Figure 6 :
Figure 6: Crack paths for offset cracks for different increment length.

Figure 7 :
Figure 7: Structural life prediction for cracked plate shown in fig.5(b).Comparion of different increment lengths on total load cycles.

Figure 8 :
Figure 8: Plot of mode I and mode II stress intensity factors against crack length for cracked plate shown in fig.5(b).

Figure 9 :
Figure 9: Multi-perforated panel geometry and initial crack positions (dimensions in mm).

Figure 11 :
Figure 11: Plot of mode I and mode II stress intensity factors against crack length for the perforated panel shown in fig.10.