- Email: [email protected]
Assessment of three-dimensional crack growth in ductile layered material systems
Engineering Fracture Mechanics 88 (2012) 15–27
Contents lists available at SciVerse ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Assessment of three-dimensional crack growth in ductile layered material systems A. Burke-Veliz a,⇑, S. Syngellakis b, P.A.S. Reed b a b
CIMA-Automotive Engineering Research Center, Tecnologico de Monterrey Campus Toluca, Toluca, Estado de Mexico 50110, Mexico School of Engineering Sciences, University of Southampton, Southampton SO17 1BJ, United Kingdom
a r t i c l e
i n f o
Article history: Received 9 February 2011 Received in revised form 9 December 2011 Accepted 23 March 2012
Keywords: FE Fatigue crack growth Shielding Ductile
a b s t r a c t Fatigue crack growth in three-dimensions may follow complex patterns that depend on the local and global conditions of the material around the crack front. In layered material systems, crack growth may be enhanced or delayed over portions of the crack front as the latter approaches a dissimilar layer. Three-dimensional elasto-plastic finite element analyses were developed for the study of the effects of dissimilar layers ahead of the crack front on the crack driving force (crack tip opening displacement) and crack growth; the simulation of the latter was based on previous experimental data obtained from a fatigued specimen with identical tri-layer architecture. The efficiency of the model under elastic conditions was first assessed by comparison of its predictions with an exact analytical solution. Crack growth in the top layer was simulated up to a depth of around 95% of its thickness. The modelling was based on a re-meshing scheme and was applied to bi-layer and tri-layer architectures subjected to three-point bending. The crack front was positioned as close as possible to the dissimilar layer so that shielding and anti-shielding effects would be clearly observed. These three-dimensional simulations revealed that such effects are less intense than those predicted by previous 2D analyses that assumed through-width cracks. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction The search for more efficient and lighter designs has promoted the development of tougher materials and more complex structures, which, along with more demanding service conditions, have resulted in less predictable patterns of damage. Due to such conditions, stress and crack propagation analyses of modern engineering material systems require the use of efficient and reliable numerical models capturing the geometric and material complexity of components and, in the case of fatigue, damage tolerance approaches properly accounting for the development of plasticity. The combination of different material properties presents a wide range of possibilities in creating new products. For such multifunctional components, a compromise needs to be found between material strength and damage tolerance on the one hand and performance requirements such as improved wear resistance and environmental protection, on the other. The layered architecture discussed in this paper is representative of a typical plain bearing design used in the automotive industry. Variation of growth rate due to mechanical properties mismatch ahead of the crack path has been the subject of many studies. This mismatch may be due to the presence of secondary phase particles or, as in the present case, of layers. Pioneering studies by Suresh and co-workers [1–3] on shielding and anti-shielding of crack growth through plastically mismatched layers were based on the J integral as the crack driving force (CDF). The crack tip opening displacement (CTOD) has also been used to assess the growth potential of stationary [4,5] and fatigue cracks [6] in diverse materials under elasto-plastic ⇑ Corresponding author. Tel.: +52 722 279 3193; fax: +52 722 2799990x2107. E-mail address: [email protected] (A. Burke-Veliz). 0013-7944/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2012.03.007
16
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
Nomenclature CDF CTOD K tb tL1 ti tL2 w dm ebr ebc esr esc Damax CDFmax Dai (da/dN)i (da/dN)max C and n CTODi De P a CTODmin CTODmax J
crack driving force crack tip opening displacement stress intensity factor backing layer thickness lining layer of thickness of tri-layer architecture interlayer thickness lining layer of thickness of bi-layer architecture (ti + tL1) strip width distance from crack tip to nodes where CTOD are measured element size in the radial direction for blunted crack tips element size in the circumferential direction for blunted crack tips element size in the radial direction for sipiderweb crack tips element size in the circumferential direction for sipiderweb crack tips maximum crack extension step maximum Crack driving force value along the front crack front advance at position i per load cycle crack growth at position i per load cycle maximum crack growth at position i per load cycle fitting parameters for crack growth crack tip opening displacement at position i strain range load applied to three point bending test crack length minimum crack tip opening displacement along the front maximum crack tip opening displacement along the front J integral estimate
conditions. Previous work in the literature has shown that growth behaviour depends on local material conditions nearby the crack, especially when heterogeneous materials are analysed. The study of shielding and anti-shielding phenomena found in the literature is based on two-dimensional models of growth and little work has been found on modelling elliptical or quasi-elliptical crack growth, which would reveal the three-dimensional shielding effect on such cracks. Stress and displacement fields around elliptical cracks were initially investigated through analytical solutions for infinite elastic and isotropic solids [8–10]. Later work [7,8] provided expressions for the CDF (represented by K) for finite solids and specific loading configurations based on finite element (FE) analyses using a structured mesh resembling an extrusion of the original 2D crack tip mesh configuration [9]. A major issue arising in these studies is the constraint on the crack front shape to be perfectly elliptical. Crack surfaces were thus defined by only the ellipse major and minor axis and the CDF values at the ends of these axes dictated the crack surface extension. Lin and Smith [10], using elastic FE analysis, compared predictions of semi-elliptical and quasi semi-elliptical crack fronts under identical loading conditions and concluded that the crack aspect ratio could differ by up to 20%, especially under bending. This work also showed that loading conditions could affect locally the crack front leading to uneven growth and changes on the crack surface. One major challenge in the development of crack propagation analyses in 3D is the extension or evolution of the surface crack and the numerical model representing it. The extension of the surface crack depends on the growth conditions that the crack front experiences locally leading to refined meshes for a more accurate local estimation at higher computational costs. The crack front evolution can be forecast through re-meshing techniques [11,12], at a high computational cost, or through the implementation of X-FEM [13], which is still not commercially available for complex architectures and under development for a broader range of applications such as interacting discontinuities (crack and material interface). X-FEM and remeshing techniques do not depend on any assumptions regarding crack growth trends but rather estimate the crack front evolution according to the chosen form of crack growth law as described by Riddell et al. [9]. The final shape of the crack surface is gradually formed by individual extension steps leading to possibly complex crack fronts and surfaces. The objective of the work described in this paper is the study of the evolution of a three-dimensional crack front in a multi-layer material system, representative of automotive plain bearing architecture. Bending loading conditions are analysed to a greater extent in this work due to its similarity to the actual working conditions experienced by the bearing. Tensile loading has also been analysed, but to a lesser extent. As a preliminary step, the model efficiency and accuracy is assessed by comparing its predictions with the analytical solution for an embedded crack in an infinite elastic solid [14]. Then, crack driving force estimations in three dimensions are compared to those from previous two-dimensional analyses of the same layered architecture under the same loading and boundary conditions. Thus, differences in shielding and anti-shielding
17
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
predictions are identified and the validity of using a simpler two-dimensional model to analyse a complex three-dimensional problem is re-assessed. 2. Materials and specimen The layered flat strip modelled in this paper was previously tested and characterised in [15,16], shown in Fig. 1. It consists of three different materials: a backing steel layer of thickness tb = 1.80 mm, a lining aluminium alloy layer of thickness tL1 = 0.38 mm and an aluminium foil interlayer of thickness ti = 0.04 mm. In addition, a hypothetical bi-layer system was also analysed; this comprises the same backing layer as that of the tri-layer model and an aluminium alloy lining layer of thickness tL2 = tL1 + ti = 0.42 mm. The strip has a width w = 19.5 mm and length between supports L = 40 mm. The material characterisation of each layer material was carried out through monotonic tensile tests with specially prepared single-layer specimens [16]. The stiffness and strength of the layer materials is given in Table 1 relative to those of the lining alloy, which has a Young’s modulus of 70 GPa and a yield stress of 57 MPa. The mechanical properties of each layer material are entered into the finite element code as the true stress–strain curves of Fig. 2; thus, a realistic stiffness mismatch, at any level of plastic deformation, is accounted for in the modelling. This becomes increasingly important as the crack tip or front approaches an interface between dissimilar layers. The crack causes high stress concentration in the material ahead of the crack tip, which may reduce or increase the actual stiffness discrepancy between two adjacent layers as different levels of strain hardening develop. The developed FE models led to the assessment of crack growth initiated at the middle of the lining surface, as observed experimentally in the strip subjected to high bending stresses [16]. The simulation of crack growth in such a multilayer
Fig. 1. Three-point bending and test specimen configuration.
Table 1 Elasto-plastic mechanical properties normalised with respect to those of the lining. Layer
Material description
Yield stress ratio
Young modulus ratio
Lining Interlayer Backing
Aluminium alloy Pure aluminium Medium-carbon steel alloy
1 0.66 8.05
1 0.96 2.82
Fig. 2. True stress–true strain curves for the three-point bending model.
18
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
system was simplified in previous work by assuming a width-through crack, plane strain conditions and a blunted crack tip [15]. The three-dimensional crack growth in this specimen is here simulated assuming an initial semi-circular flaw with a radius of 5 lm, which is consistent with earlier experimental observations [16]. This specimen was also subjected to bending and tension at a loading ratio equal to 0.1 for the sake of comparison. The numerical analyses are based on the assumption of homogeneous materials for all the layers and a perfectly defined, flat interface without material diffusion between adjacent layers. It should also be noted that the element size close to the crack tip is comparable to that of typical microstructural features of the lining of the tested flat strip. The adopted minimum element size was however only aimed at simulating the non-linear deformation process and the positioning of crack tip at small distances from the layers’ interface where shielding and anti-shielding is most intense. 3. Methods In this work, the methodology described by Lin and Smith [17], Riddell et al. [9] and Carter et al. [11] is adopted, in general terms; this is based on re-meshing schemes following the extension of the crack and the use of DCTOD as a crack driving force (CDF) parameter accounting for significant elasto-plastic behaviour around the crack tip. The numerical analysis involves the generation of crack front elements that are suitable for the extraction of CDF values (standard crack mesh configurations to calculate J-integral, K or CTOD), estimating the crack front advance according to a crack growth law and, finally, re-meshing. These steps are repeated until a critical or pre-defined crack size is reached. 3.1. Crack -element generation For the development of FE models to predict crack growth, many concepts explained in great detail by Carter et al. [11] were adopted. Key issues identified by Carter et al. are the volume decomposition into smaller entities and modification of the volume containing the crack elements, so that the crack front and the boundaries of the volume containing it are sufficiently separated and thus a good quality mesh could be generated. A similar decomposition into volumes comprising three main element groups was carried out here according to the element functionality: elements with low-stress gradients, crack tip elements and the transition elements between the first two groups. This decomposition is shown in Fig. 3. The crack tip elements, comprising collapsed bricks at the crack front and patterned surrounding elements, are created in a direct manner by specifying nodal coordinates on a plane oriented perpendicularly to the crack front. Such planes contain blunted or spider web crack tip patterns to incorporate the desired degree of strain and stress singularity that facilitates the estimation of the crack driving force. The crack front mesh is then defined by a set of planes on which crack driving force estimates are obtained. In other words, the crack front elements are created so that they form an extrusion of the crack tip pattern used on the planes. The section with the desired crack pattern is created perpendicular to the crack front to improve the computation accuracy of CDF, usually based on 2D methods [9,11,17]. The elements with low-stress gradients and the transition elements can be generated through the available automatic meshing functionality so that they fill the volumes formed by the boundaries of the flat strip (component outer surfaces and interfaces) and the outer surface of the volume comprising the crack tip elements. A mixed-approach involving direct generation and automatic meshing provides the best results for the creation of models containing quasi semi-elliptical shapes or more complex crack surface patterns such as coalescing cracks. The creation of the crack tip elements relied on ANSYS [18] modelling tools that allow the use of local coordinate systems. Such local coordinate systems were positioned and oriented according to the CDF results at the previous step, which were
Fig. 3. Typical meshes of: (a) the lining surface and (b) a cross section containing the crack (with the von Mises stress contours MPa).
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
19
post-processed through MATLAB routines developed for that purpose. These MATLAB routines included modelling, solving, post-processing and crack extension steps and were used to control the crack growth simulation in a continuing loop until a desired crack length is obtained. 3.2. Crack driving force post-processing The post-processing operations for generating CDF values was carried out following the approach proposed by Lin and Smith. Blunted tip and spider web configurations have pairs of directly opposite nodes at the crack faces, thus allowing efficient calculation of the opening and sliding displacement between the faces. The CTOD estimation was based on the relative movement of these opposite nodes at a particular distance away from the crack tip. The selection of DCTOD over J within the elasto-plastic fracture mechanics (EPFMs) framework was due to the possibility of assessing crack fronts located closer to the material interface and calculation simplicity. J estimates in 2D problems at slightly greater distances from the tip to the material interface have required additional meshing features and complexity [19]. The eventually selected pair of nodes was positioned at a distance dm = 960 nm and dm = 900 nm away from the tip for the blunted and spider web models, respectively; this is schematically shown in Figs. 4 and 5. These locations were considered sufficiently close to crack tip to provide reliable estimates of the CTOD. Typical minimum element sizes for the blunted tip in the radial and circumferential direction were ebr = 300 nm and ebc = 26.1 nm, respectively; while for the spider-web elements, esr = 300 nm and esc = 82.1 nm. The element size in the present 3D analyses was slightly larger than that used in the earlier 2D analysis developed on the same material architecture [15]. 3.3. Crack front advance and remeshing The crack front advance was based on the growth around the crack front relative to a specified maximum extension step Damax, corresponding to the CDFmax, according to the relation [9]
da Dai ¼ dadN
i
dN max
Damax ¼
f ðCDFi Þ Damax f ðCDFmax Þ
ð1Þ
where Dai refers to the crack advance at location i, (da/dN)i to the expected growth rate at location i, and (da/dN)max to the maximum growth rate along the crack front. The use of small values of crack advance is recommended for problems with sudden crack front changes, deflections that cause out of plane crack growth and other factors that contribute to the formation of irregular crack surfaces. This condition also applies to changing or uneven loading conditions existing along
Fig. 4. CTOD measurement in blunted tip mesh configuration.
Fig. 5. CTOD measurement in spider-web mesh configuration.
20
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
Fig. 6. Relative crack extension along the front.
the material on which the crack front is embedded, such as bending where the material stress depends on its position to the neutral axis. The estimations of da/dN were based on a crack growth law in terms of the CTOD as the crack driving force as proposed by Tanaka et al. [20] for small cracks. This law was based on a study of ductile materials under cyclic loading and led to satisfactory results when compared to experimental data [20]. The form of the relationship proposed by Tanaka et al. is
da ¼ CðDCTODÞn dN
ð2Þ
where fitting parameters C and n account for the material behaviour. The n exponent values obtained for materials such as copper and steel range between 1.17 and 1.13. Crack growth laws of this form can also be found in the work of Tomkins [21] where the value n = 1 is adopted. The present analyses were only concerned with the proportional growth of the crack front according to the crack driving force estimations. Therefore, the selected value of C would not affect the prediction of the crack front growth process as shown when Eq. (1) is substituted in Eq. (2) resulting in
da Dai ¼ dadN
i
dN max
Damax ¼
DCTODni Damax DCTODnmax
ð3Þ
The value of n = 1.129 for the Al alloy lining was estimated through a combination of experimental data and numerical simulation [22] of crack growth in the lining layer of the tri-layered flat strip measured during interrupted three-point bending fatigue tests [16] under a loading ratio 0.1 using acetate strips to replicate the cracked surface. These tests were performed at a maximum strain range De = 0.0063 ± 0.0002, which corresponds to a maximum load P = 920 N. Re-meshing operations with this methodology only require the generation of a new crack front and transition elements. The new crack elements along the crack front were created according to the CDF estimates at the previous step. As the crack front length grows, the number of elements is automatically increased, as shown in Fig. 6, so that a fairly constant element size is maintained throughout the crack propagation analysis. 4. Comparison of crack tip modelling The generation of FE models for the study of fracture mechanics has used distinct crack tip mesh patterns that facilitate the computation of the crack driving force. Previous 2D work on multi-layered systems used refined blunted tips [15] for the study of shielding and anti-shielding effects under extensive elasto-plastic deformation. The development of 3D models with a similar element size was expected to result in considerable computational cost. Therefore, the creation of spider web crack tip meshes, adopting quarter-point and collapsed brick elements [23] for elastic and elasto-plastic analyses, respectively, was also considered. The suitability of blunted or initially sharp tips in spider web configurations was numerically tested as an initial validation step before the simulation of crack growth and the assessment of shielding and anti-shielding effects. This initial linear elastic analysis was carried out to determine the accuracy and numerical stability of the CTOD estimated from both models and its sensitivity to mesh refinement along the crack front; it was applied to a specimen containing an embedded crack, subjected to tension as shown schematically in Fig. 7. The dimensions of the analysed plate are the same as those of the flat strip specimen shown in Fig. 1 apart from its thickness which is taken equal to 2tL2 = 0.76 mm, for the purpose of assessing the effectiveness of a mesh representative of those subsequently used for analysing multi-layered systems of equal size. In this validation test, the surface crack is replaced by a penny crack and the free-of-traction conditions on the face of the lining are replaced by symmetry conditions. The embedded crack has a radius of 0.266tL2, is oriented perpendicular to the plane of
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
21
Fig. 7. Geometry of embedded crack under tension for numerical (left) and exact solution (right).
symmetry and subjected to a remote uniform tensile stress of 20 MPa. The material used for these analyses was assumed homogeneous and elastic, having the lining properties, that is, a Young’s modulus of 70 GPa and a Poisson’s ratio of 0.33. In order to assess the accuracy and numerical stability of the developed meshes a comparison with the exact solution for a flat ellipsoidal crack was obtained by Green and Sneddon [14]. According to this solution, the crack opening for a penny shaped flaw is given by
COD ¼
4r1 ð1 m2 Þ Ep
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dm ð2a dm Þ
ð4Þ
where a is the crack radius and dm the normal (radial) distance from the crack front. This solution and the problem’s axis symmetric nature indicate that for a given crack of radius a the crack opening at a distance dm from the tip would be constant irrespective of its position along the crack front. This refinement varied the number of elements from 18 to 54 providing an initial estimate of the computational cost expected from subsequent elasto-plastic analyses. Each elastic analysis using the ANSYS Sparse Direct Solver required from 1 to 3 h of computing time using four Intel Xeon processors with 8 GB of RAM memory analysing a model comprising between 60,000 and 120,000 nodes. The CTOD estimates along the crack front in both crack tip configurations showed an inherent variability generated by the numerical approximation of the model. Such variability is in contrast to a constant value yielded by the exact solution along the whole front. The standard deviation calculation of CTOD values along the crack front was developed for seven levels of refinement in both crack tip configurations and is shown in Table 2. The standard deviation obtained at the different levels of refinement showed that models using spider web crack tip configurations behaved more closely to the analytical solution and values up to 20 times smaller to those of blunted models were found. Mesh refinement also proved to reduce the standard deviation estimations, as expected. These studies also revealed that CTOD estimates obtained from spider web configurations are closer to the analytical solution, as expected. According to the mesh sensitivity analysis, the degree of accuracy achieved using the spider web crack tip configuration showed error estimates between 1% and 5%. The best compromise between accuracy en computational costs was achieved with 30 elements along the front and an element size of 5.85 lm with error estimates of 2% and only 1 h of simulation time.
Table 2 Standard deviation of normalised CTOD along the crack front. Number of elements along the front
Blunted configuration
Spider web configuration
54 46 38 30 26 22 18
2.3E03 3.6E03 4.0E03 4.8E03 5.1E03 6.6E03 7.2E03
4.1E04 5.2E04 8.2E04 2.0E03 1.4E03 3.2E03 6.5E03
22
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
5. Crack growth modelling The proportional growth along the crack front was modelled here according to Eq. (4). The crack growth study was carried out at tension and bending loading conditions with a maximum tensile load Pmax = 2684 N and a maximum bending load Pmax = 200 N and loading ratio of 0.1, which causes substantial yielding around the crack. The current loading and modelling conditions allow the analysis of crack fronts placed closely to the layers’ interface, where shielding and anti-shielding effects are most intense. The models developed to account for various crack sizes comprised between 60,000 and 1,20,000 elements with 1,40,000 to 2,40,000 nodes. This significant increase in the number of DOFs during the solution process arises from the need to keep the element size constant as the crack grows. The processing time for the analysis of each growth step required between 12 and 26 h and result files for both layered architectures that exceeded 50 GB of hard disc per loadcase. The individual computational cost for the two modelled architectures at each loadcase was comparable despite the greater propagation depths of 0.4 mm for the bilinear one in comparison to 0.36 mm for the tri-linear. Crack front proximity to the lining-interlayer interface in the tri-layer resulted in slow convergence rates due to the greater compliance of the interlayer. The crack front evolution in bi-layer and tri-layer architectures was almost identical for crack depths less than 0.27 mm. Small differences could be observed at greater depths, which were consistent with the CDF trends in respective two-dimensional analyses but yielded almost identical outcomes until the crack tip came sufficiently close to the layers interface. The difference between the predicted crack fronts in tri-layer and bi-layer became obvious at depths greater than 0.3 mm.
5.1. Crack driving force evolution and shielding The evolution of the estimated values of CTOD and their pattern along the crack front in both layered systems subjected to bending are shown in Fig. 8. The similarities of the two plots may be noted. The CTOD pattern had maximum values at the surface while the minimum was located at the deepest point of the crack front as observed by previous authors when bending stresses are applied [10,17,24]. The sampling points at which the CTOD was estimated along the front were positioned at equal distances; in Fig. 9 the front contour coordinate has been normalised with respect to the length of the front. The evolution of CTOD with crack depth at two key locations along the crack front, for both layered systems, is shown in Fig. 9. In both cases, the CTOD estimates at the surface showed very similar values and rising trends. At the deepest point, they were almost identical for crack depths less than, approximately, 0.31 mm; they were also rising but at a lower rate than at the surface. The CTOD values for the tri-layer and bi-layer models became significantly different as the crack fronts approached the respective interfaces. Shielding and anti-shielding effects were therefore observed at the deepest point of the crack front at depths greater than 0.31 mm (82% of the lining thickness). Crack growth beyond this point showed accelerating and decelerating trends according to the adjacent layer stiffness. However, it should be noted that these trends did not affect to a great extent the magnitude of the estimated CTOD. The disproportional growth along the front, shown through DCTOD in Fig. 9, can be better represented by the CDF ratio of points A and B, which respectively represent the points located nearest and furthest to the materials interface. Fig. 10 demonstrates more clearly the evolution of these two parameters for the bending load case. It is thus possible to identify three different growth phases: initial transition from semi-circular to quasi semi-elliptical (I), crack growth under combined
(a)
(b)
Fig. 8. Crack front evolution in: (a) tri-layer and (b) bi-layer architectures under P = 200 N.
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
23
Fig. 9. Crack driving force evolution in bi-layer and tri-layer architectures subjected to three point bending.
Fig. 10. Ratio of CDF obtained at the deepest point (CTODd) to that at the surface (CTODs) in bi-layer and tri-layer architectures subjected to three-point bending.
tension-bending conditions (II) and the start of shielding and anti-shielding mechanisms (III). The initial transition phase corresponds to the changes that the crack front must undergo to achieve a more ‘‘natural’’ shape consistent with the loading conditions. The rate of change depended on the difference between the current and ‘‘natural’’ shape. The second phase in the plot of Fig. 10 shows the typical results of crack growth in a monolithic specimen subjected to tension-bending with maximum CDF values at the free surface where the tensile bending stress also attains its maximum value. As the crack front grew progressively deeper into the strip, lower CDF ratios were estimated. This behaviour has been previously observed [10] in the crack growth in monolithic strips and it is associated with significant change of the crack aspect ratio as the crack front approaches the neutral axis of the strip. However, the multi-layer architecture produced a different stress distribution to that observed in a monolithic strip. The stress distribution through the lining thickness of an undamaged specimen is equivalent to a combined uniform tension-bending loading on a monolithic strip. Under a load P = 200 N, the obtained stress distribution was equivalent to a uniform tension of about 53.5 MPa and maximum bending stress of about 8.5 MPa, thus uniform tension was dominant. The ratio between uniform and bending tensile stresses is a key indicator of how stress and deformation are distributed along the crack front; however, it is important to consider that these conditions may vary as the crack extends and plasticity spreads throughout the material. The third phase in Fig. 10 shows clearly the appearance of shielding and anti-shielding. Shielding is associated with increasing rate at which CTODB/CTODA, drops in the bi-layer architecture. In contrast, due to anti-shielding, the rate at which CTODB/CTODA decreases remains almost constant. If the numerical simulation had allowed the crack to grow further, closer to the layers’ interface, the anti-shielding effect may have become more obvious. However, the small separation between the
24
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
front and the interface could not allow a significant shape change before the crack penetrates into the interlayer. Pursuing the crack growth within the softer interlayer, as in previous two-dimensional analyses [15], would have elucidated the local behaviour further. This however was beyond the available computational capacity required by the adopted model but offers prospects for future interesting work. The 3D modelling challenges of a front placed within two different materials and the computational challenges arising from even larger plastic deformations in the softer interlayer point to the complexity of such a simulation. The crack growth analyses applied to the bi-layer architecture produced sudden drops and rises in the evolution of CDF when the crack depth exceeded 0.34 mm (90% of the lining thickness), as shown in the enlarged view of the inset in Fig. 9a (Points 1, 2 and 3). The existence of these gradient discontinuities may be due to the extension steps used here not being sufficiently small to capture the smooth changes of CDF. It should be noted that a front shape is produced here as a consequence of combining the shielding effect with that of the overall bending deformation state governing the crack front extension. Despite the existence of this local non-smoothness, the overall shielding trend is clear. Therefore, model refinement to achieve further crack extension was not considered necessary. On the other hand, crack growth under tensile loading showed similar trends to bending in areas such as the existence of three evolution areas and overall shielding and anti-shielding behaviour, shown in Fig. 11. In this particular case, the CDF ratio between points A and B indicated more clearly how the existence of a more compliant material would increase the CDF ratio between A and B as the crack tip approaches to the interface while a stiffer material would cause an opposite effect. These changes become more evident as the distance between the material interface and the crack tip is reduced. The gradual shift towards greater or smaller CDF ratios demonstrates the need of developing consecutive elasto-plastic analysis for a reliable estimation of CDF.
1.25 1.2 1.15 1.1 1.05 1 0.95 0.9 0.85 0.8
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Fig. 11. Ratio of CDF obtained at the deepest point (CTODd) to that at the surface (CTODs) in bi-layer and tri-layer architectures subjected to pure tension.
-4
x 10 4
3
Bi-layer 3D deepest point Tri-layer 3D deepest point Bi-layer 2D Tri-layer 2D
2
1 0.05
0.1
0.15
0.2
0.25
Fig. 12. Comparison between CDF estimates obtained from 2D and 3D models.
0.3
0.35
0.4
0.45
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
25
5.2. Comparison to 2D analyses Shielding or amplification arising from crack growth through dissimilar materials has been studied in the past through two-dimensional analyses, Despite their simplified representation of the actual problem, two-dimensional models may provide a reliable initial assessment before the development of more rigorous but also more complicated three-dimensional modelling or a good approximation to the correct answer when a quick output is required. It is evident from Fig. 12 that the two-dimensional simplification of the three-point bending simulation reduces the stiffness of the specimen more dramatically by modelling a width-through crack instead of half-penny type crack. The CDF estimates from 2D and 3D analyses of the flat strip subjected to a load P = 200 N were clearly different, as shown in Fig. 12b. The higher CDF estimates from the 2D model are expected since a through crack along the strip width is assumed. The acceleration and deceleration effects on crack growth caused by shielding and anti-shielding in 3D models did not appear as strong as predicted by two-dimensional analyses. This difference can be attributed to the more constrained growth of the surface crack shape as it expands as a whole. Shielding and anti-shielding are localised phenomena that only affect small segments of the crack front. The small separation between the crack front and the layers’ interface would affect the CDF values locally causing disproportional growth that modifies the surface crack profile. At the same time, the neighbouring material discourages the disproportional growth of the most affected areas. In other words, shielding and anti-shielding effects are averaged or diluted into a greater portion of material based on the crack front growth in its entirety. 5.3. Plastic zone The plastic zone around the crack tip has been extensively studied in test specimens with width-through cracks. Tri-axiality effects have been of interest in these specimens, often accounted for by assuming plane strain or plane stress conditions according to the location of the crack tip in the specimen. The shape of the plastic zone around a half penny-type crack is more complicated, especially when the loading conditions differ from uniform remote tension normal to the crack surface within an infinite solid. Fig. 13 shows typical von Mises plastic strain contours around the crack front for the bi-layer architecture subjected to bending and how shielding occurs by a stiffer and stronger backing. This condition can be observed where the crack front is the closest to the interface through the development of a broader plastic zone in this region with higher levels of plasticity. It is expected that the aforementioned broadening may facilitate the crack initiation process in the backing layer by increasing the stress concentration area within the backing. In contrast, as the crack front approaches the interlayer, in the tri-layer architecture, high values of plastic strain can be observed in this intermediate layer. For clarity, only deformations above 5e4 are shown since plastic deformation covers around 20% of the top surface. The lower yield strength of the interlayer produced a protuberance in the plastic zone shape as shown in Fig. 14. This protuberance expands through the interlayer forming a secondary plastic zone that is expected to promote faster crack propagation.
Fig. 13. Plastic zone broadening at the deepest point of the crack due to backing shielding in bi-layer architecture (Von Mises plastic strain contour).
Fig. 14. Plastic zone shape and low strain area next to the crack in the tri-layer architecture (Von Mises plastic strain contour showing only strains greater than 5e4 for clarity).
26
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
6. Discussion and conclusions The adopted methodology put emphasis on the evolution of the crack front as it approaches a dissimilar layer to assess shielding and anti-shielding effects on half-penny cracks. Modelling crack growth in systems comprising dissimilar layers that exhibit both elastic and plastic mismatch made the analysis complex and time consuming. The shielding and anti-shielding effects were predicted to be less significant in three-dimensional crack growth in comparison to those arising from width-through (or through-thickness) cracks in two-dimensional specimens. The crack front segments affected most intensively by shielding and anti-shielding in 3D are expected to accelerate or decelerate, respectively, crack growth in the neighbouring material. At the same time, the neighbouring material discourages the disproportional growth of the parts along the crack front most affected by shielding or anti-shielding. Crack growth processes are defined by the local conditions of material and loading around the crack front. Local loading and material conditions may change as the front extends leading to shapes that cannot be forecasted a priori, in contrast to homogeneous materials subjected to typical loading cases [10,17,24] and their local loading profile on the front. The study of shielding and anti-shielding in this work displayed loading and material conditions that vary significantly as the crack extends towards a dissimilar layer. At the same time, the elasto-plastic nature of the multi-layered architecture displays a stiffness gradient from the crack front to the interface which will change according to the stress concentration conditions, loading, and stress strain curve of the involved materials. An appropriate estimation of crack driving forces requires updated data regarding CDF and the state of the material which leaded to computationally expensive analysis and sophisticated modelling routines. The effects of shielding and anti-shielding were studied in layered architectures subjected to bending and tension. The crack growth under bending loading revealed localised speed reductions as the front approached the neutral axis. Additionally, shielding and anti-shielding effects caused less significant modifications to the front profile. The study of such layered systems under tensile loading revealed more clearly accelerations or reductions of CDF as the front approached a dissimilar material. Occurrences of crack deflection and bifurcation observed experimentally still remain to be investigated numerically; such phenomena lead to interesting crack development, especially through the interlayer. The simulation of surface crack extension into the interlayer offers great modelling challenges: from the prediction of three-dimensional deflected and bifurcated forms to the high level of mesh refinement necessitated by the thinness of the interlayer. The development of such analyses with elastic materials would significantly simplify the task; it is however important, from a practical point of view, to identify deflections and bifurcations using elasto-plastic material models, as in previous two-dimensional studies [15].
Acknowledgments The financial support of the Mexican Council of Science and Technology and the School of Engineering Sciences is gratefully acknowledged. The training provided by and discussions with MAHLE Engine Systems was also invaluable in carrying out this research.
References [1] Sugimura Y et al. Fracture normal to a bimaterial interface: Effects of plasticity on crack-tip shielding and amplification. Acta Metall Mater 1995;43(3):1157–69. [2] Sugimura Y, Grondin L, Suresh S. Fatigue crack growth at arbitrary angles to bimaterial interfaces. Scripta Metall Mater 1995;33(12):2007–12. [3] Suresh S, Sugimura Y, Tschegg EK. Growth of a fatigue crack approaching a perpendicularly oriented bimaterial interface. Scripta Metall Mater 1992;27:1189–94. [4] Pippan R, Flechsig K, Riemelmoser FO. Fatigue crack propagation behavior in the vicinity of an interface between materials with different yield stresses. Mater Sci Engng, A 2000;283(1–2):225–33. [5] Wang B, Siegmund T. Simulation of fatigue crack growth at plastically mismatched bi-material interfaces. Int J Plast 2006;22(9):1586–609. [6] Groh S et al. Fatigue crack growth from a cracked elastic particle into a ductile matrix. Phil Mag 2008;88(30):3565–83. [7] Newman Jr JC, Raju IS. An empirical stress-intensity factor equation for the surface crack. Engng Fract Mech 1981;15(1-2):185–92. [8] Raju IS, Newman JC. Stress-intensity factors for a wide range of semi-elliptical surface cracks in finite-thickness plates. Engng Fract Mech 1979;11:817–29. [9] Riddell WT, Ingraffea AR, Wawrzynek PA. Experimental observations and numerical predictions of three-dimensional fatigue crack propagation. Engng Fract Mech 1997;58(4):293–310. [10] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates. Part II: crack shape change. Engng Fract Mech 1999;63(5):523–40. [11] Carter BJ, Wawrzynek PA, Ingraffea AR. Automated 3-D crack growth simulation. Int J Numer Method Engng 2000;47(1–3):229–53. [12] Schollmann M, Fulland M, Richard HA. Development of a new software for adaptive crack growth simulations in 3D structures. Engng Fract Mech 2003;70(2):249–68. [13] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engng Fract Mech 2008;75(5):943–60. [14] Green AE, Sneddon IN. The distribution of stress in the neighbourhood of a flat elliptical crack in an elastic solid. Math Proc Cambridge Philos Soc 1950;46(1):159–63. [15] Burke-Veliz A, Reed PAS, Syngellakis S. A numerical study of crack shielding and deflection under extensive plasticity. Engng Fract Mech 2009;76(9):1345–56. [16] Ali M, Syngellakis S, Reed PAS. A comparison of fatigue performance of HVOF spray coated and conventional roll bonded aluminium bearing alloys. Mater Sci Technol 2008;25(5):575–81.
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
27
[17] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates. Part I: the numerical technique. Engng Fract Mech 1999;63(5):503–22. [18] Ansys, Ansys reference manual; 2005. [19] Kolednik O et al. On the fracture behavior of inhomogeneous materials – a case study for elastically inhomogeneous bimaterials. Int J Solids Struct 2005;42(2):605–20. [20] Tanaka K, Hoshide T, Sakai N. Mechanics of fatigue crack propagation by crack tip plastic blunting. Engng Fract Mech 1984;19:805–25. [21] Tomkins B. Fatigue crack propagation: an analysis. Phil Mag 1968;18(155):1041–66. [22] Burke-Veliz, A., Finite element modelling of crack growth in multi-layered architectures. In: Computational engineering and design group. PhD Thesis. Southatmpton: University of Southampton; 2009. [23] Barsoum R. Triangular quarter point elements as elastic and perfectly plastic crack tip elements. Int J Numer Meth Engng 1977;11:85–98. [24] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates. Part III: stress intensity factor and fatigue crack growth life. Engng Fract Mech 1999;63(5):541–56.
Contents lists available at SciVerse ScienceDirect
Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech
Assessment of three-dimensional crack growth in ductile layered material systems A. Burke-Veliz a,⇑, S. Syngellakis b, P.A.S. Reed b a b
CIMA-Automotive Engineering Research Center, Tecnologico de Monterrey Campus Toluca, Toluca, Estado de Mexico 50110, Mexico School of Engineering Sciences, University of Southampton, Southampton SO17 1BJ, United Kingdom
a r t i c l e
i n f o
Article history: Received 9 February 2011 Received in revised form 9 December 2011 Accepted 23 March 2012
Keywords: FE Fatigue crack growth Shielding Ductile
a b s t r a c t Fatigue crack growth in three-dimensions may follow complex patterns that depend on the local and global conditions of the material around the crack front. In layered material systems, crack growth may be enhanced or delayed over portions of the crack front as the latter approaches a dissimilar layer. Three-dimensional elasto-plastic finite element analyses were developed for the study of the effects of dissimilar layers ahead of the crack front on the crack driving force (crack tip opening displacement) and crack growth; the simulation of the latter was based on previous experimental data obtained from a fatigued specimen with identical tri-layer architecture. The efficiency of the model under elastic conditions was first assessed by comparison of its predictions with an exact analytical solution. Crack growth in the top layer was simulated up to a depth of around 95% of its thickness. The modelling was based on a re-meshing scheme and was applied to bi-layer and tri-layer architectures subjected to three-point bending. The crack front was positioned as close as possible to the dissimilar layer so that shielding and anti-shielding effects would be clearly observed. These three-dimensional simulations revealed that such effects are less intense than those predicted by previous 2D analyses that assumed through-width cracks. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction The search for more efficient and lighter designs has promoted the development of tougher materials and more complex structures, which, along with more demanding service conditions, have resulted in less predictable patterns of damage. Due to such conditions, stress and crack propagation analyses of modern engineering material systems require the use of efficient and reliable numerical models capturing the geometric and material complexity of components and, in the case of fatigue, damage tolerance approaches properly accounting for the development of plasticity. The combination of different material properties presents a wide range of possibilities in creating new products. For such multifunctional components, a compromise needs to be found between material strength and damage tolerance on the one hand and performance requirements such as improved wear resistance and environmental protection, on the other. The layered architecture discussed in this paper is representative of a typical plain bearing design used in the automotive industry. Variation of growth rate due to mechanical properties mismatch ahead of the crack path has been the subject of many studies. This mismatch may be due to the presence of secondary phase particles or, as in the present case, of layers. Pioneering studies by Suresh and co-workers [1–3] on shielding and anti-shielding of crack growth through plastically mismatched layers were based on the J integral as the crack driving force (CDF). The crack tip opening displacement (CTOD) has also been used to assess the growth potential of stationary [4,5] and fatigue cracks [6] in diverse materials under elasto-plastic ⇑ Corresponding author. Tel.: +52 722 279 3193; fax: +52 722 2799990x2107. E-mail address: [email protected] (A. Burke-Veliz). 0013-7944/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2012.03.007
16
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
Nomenclature CDF CTOD K tb tL1 ti tL2 w dm ebr ebc esr esc Damax CDFmax Dai (da/dN)i (da/dN)max C and n CTODi De P a CTODmin CTODmax J
crack driving force crack tip opening displacement stress intensity factor backing layer thickness lining layer of thickness of tri-layer architecture interlayer thickness lining layer of thickness of bi-layer architecture (ti + tL1) strip width distance from crack tip to nodes where CTOD are measured element size in the radial direction for blunted crack tips element size in the circumferential direction for blunted crack tips element size in the radial direction for sipiderweb crack tips element size in the circumferential direction for sipiderweb crack tips maximum crack extension step maximum Crack driving force value along the front crack front advance at position i per load cycle crack growth at position i per load cycle maximum crack growth at position i per load cycle fitting parameters for crack growth crack tip opening displacement at position i strain range load applied to three point bending test crack length minimum crack tip opening displacement along the front maximum crack tip opening displacement along the front J integral estimate
conditions. Previous work in the literature has shown that growth behaviour depends on local material conditions nearby the crack, especially when heterogeneous materials are analysed. The study of shielding and anti-shielding phenomena found in the literature is based on two-dimensional models of growth and little work has been found on modelling elliptical or quasi-elliptical crack growth, which would reveal the three-dimensional shielding effect on such cracks. Stress and displacement fields around elliptical cracks were initially investigated through analytical solutions for infinite elastic and isotropic solids [8–10]. Later work [7,8] provided expressions for the CDF (represented by K) for finite solids and specific loading configurations based on finite element (FE) analyses using a structured mesh resembling an extrusion of the original 2D crack tip mesh configuration [9]. A major issue arising in these studies is the constraint on the crack front shape to be perfectly elliptical. Crack surfaces were thus defined by only the ellipse major and minor axis and the CDF values at the ends of these axes dictated the crack surface extension. Lin and Smith [10], using elastic FE analysis, compared predictions of semi-elliptical and quasi semi-elliptical crack fronts under identical loading conditions and concluded that the crack aspect ratio could differ by up to 20%, especially under bending. This work also showed that loading conditions could affect locally the crack front leading to uneven growth and changes on the crack surface. One major challenge in the development of crack propagation analyses in 3D is the extension or evolution of the surface crack and the numerical model representing it. The extension of the surface crack depends on the growth conditions that the crack front experiences locally leading to refined meshes for a more accurate local estimation at higher computational costs. The crack front evolution can be forecast through re-meshing techniques [11,12], at a high computational cost, or through the implementation of X-FEM [13], which is still not commercially available for complex architectures and under development for a broader range of applications such as interacting discontinuities (crack and material interface). X-FEM and remeshing techniques do not depend on any assumptions regarding crack growth trends but rather estimate the crack front evolution according to the chosen form of crack growth law as described by Riddell et al. [9]. The final shape of the crack surface is gradually formed by individual extension steps leading to possibly complex crack fronts and surfaces. The objective of the work described in this paper is the study of the evolution of a three-dimensional crack front in a multi-layer material system, representative of automotive plain bearing architecture. Bending loading conditions are analysed to a greater extent in this work due to its similarity to the actual working conditions experienced by the bearing. Tensile loading has also been analysed, but to a lesser extent. As a preliminary step, the model efficiency and accuracy is assessed by comparing its predictions with the analytical solution for an embedded crack in an infinite elastic solid [14]. Then, crack driving force estimations in three dimensions are compared to those from previous two-dimensional analyses of the same layered architecture under the same loading and boundary conditions. Thus, differences in shielding and anti-shielding
17
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
predictions are identified and the validity of using a simpler two-dimensional model to analyse a complex three-dimensional problem is re-assessed. 2. Materials and specimen The layered flat strip modelled in this paper was previously tested and characterised in [15,16], shown in Fig. 1. It consists of three different materials: a backing steel layer of thickness tb = 1.80 mm, a lining aluminium alloy layer of thickness tL1 = 0.38 mm and an aluminium foil interlayer of thickness ti = 0.04 mm. In addition, a hypothetical bi-layer system was also analysed; this comprises the same backing layer as that of the tri-layer model and an aluminium alloy lining layer of thickness tL2 = tL1 + ti = 0.42 mm. The strip has a width w = 19.5 mm and length between supports L = 40 mm. The material characterisation of each layer material was carried out through monotonic tensile tests with specially prepared single-layer specimens [16]. The stiffness and strength of the layer materials is given in Table 1 relative to those of the lining alloy, which has a Young’s modulus of 70 GPa and a yield stress of 57 MPa. The mechanical properties of each layer material are entered into the finite element code as the true stress–strain curves of Fig. 2; thus, a realistic stiffness mismatch, at any level of plastic deformation, is accounted for in the modelling. This becomes increasingly important as the crack tip or front approaches an interface between dissimilar layers. The crack causes high stress concentration in the material ahead of the crack tip, which may reduce or increase the actual stiffness discrepancy between two adjacent layers as different levels of strain hardening develop. The developed FE models led to the assessment of crack growth initiated at the middle of the lining surface, as observed experimentally in the strip subjected to high bending stresses [16]. The simulation of crack growth in such a multilayer
Fig. 1. Three-point bending and test specimen configuration.
Table 1 Elasto-plastic mechanical properties normalised with respect to those of the lining. Layer
Material description
Yield stress ratio
Young modulus ratio
Lining Interlayer Backing
Aluminium alloy Pure aluminium Medium-carbon steel alloy
1 0.66 8.05
1 0.96 2.82
Fig. 2. True stress–true strain curves for the three-point bending model.
18
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
system was simplified in previous work by assuming a width-through crack, plane strain conditions and a blunted crack tip [15]. The three-dimensional crack growth in this specimen is here simulated assuming an initial semi-circular flaw with a radius of 5 lm, which is consistent with earlier experimental observations [16]. This specimen was also subjected to bending and tension at a loading ratio equal to 0.1 for the sake of comparison. The numerical analyses are based on the assumption of homogeneous materials for all the layers and a perfectly defined, flat interface without material diffusion between adjacent layers. It should also be noted that the element size close to the crack tip is comparable to that of typical microstructural features of the lining of the tested flat strip. The adopted minimum element size was however only aimed at simulating the non-linear deformation process and the positioning of crack tip at small distances from the layers’ interface where shielding and anti-shielding is most intense. 3. Methods In this work, the methodology described by Lin and Smith [17], Riddell et al. [9] and Carter et al. [11] is adopted, in general terms; this is based on re-meshing schemes following the extension of the crack and the use of DCTOD as a crack driving force (CDF) parameter accounting for significant elasto-plastic behaviour around the crack tip. The numerical analysis involves the generation of crack front elements that are suitable for the extraction of CDF values (standard crack mesh configurations to calculate J-integral, K or CTOD), estimating the crack front advance according to a crack growth law and, finally, re-meshing. These steps are repeated until a critical or pre-defined crack size is reached. 3.1. Crack -element generation For the development of FE models to predict crack growth, many concepts explained in great detail by Carter et al. [11] were adopted. Key issues identified by Carter et al. are the volume decomposition into smaller entities and modification of the volume containing the crack elements, so that the crack front and the boundaries of the volume containing it are sufficiently separated and thus a good quality mesh could be generated. A similar decomposition into volumes comprising three main element groups was carried out here according to the element functionality: elements with low-stress gradients, crack tip elements and the transition elements between the first two groups. This decomposition is shown in Fig. 3. The crack tip elements, comprising collapsed bricks at the crack front and patterned surrounding elements, are created in a direct manner by specifying nodal coordinates on a plane oriented perpendicularly to the crack front. Such planes contain blunted or spider web crack tip patterns to incorporate the desired degree of strain and stress singularity that facilitates the estimation of the crack driving force. The crack front mesh is then defined by a set of planes on which crack driving force estimates are obtained. In other words, the crack front elements are created so that they form an extrusion of the crack tip pattern used on the planes. The section with the desired crack pattern is created perpendicular to the crack front to improve the computation accuracy of CDF, usually based on 2D methods [9,11,17]. The elements with low-stress gradients and the transition elements can be generated through the available automatic meshing functionality so that they fill the volumes formed by the boundaries of the flat strip (component outer surfaces and interfaces) and the outer surface of the volume comprising the crack tip elements. A mixed-approach involving direct generation and automatic meshing provides the best results for the creation of models containing quasi semi-elliptical shapes or more complex crack surface patterns such as coalescing cracks. The creation of the crack tip elements relied on ANSYS [18] modelling tools that allow the use of local coordinate systems. Such local coordinate systems were positioned and oriented according to the CDF results at the previous step, which were
Fig. 3. Typical meshes of: (a) the lining surface and (b) a cross section containing the crack (with the von Mises stress contours MPa).
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
19
post-processed through MATLAB routines developed for that purpose. These MATLAB routines included modelling, solving, post-processing and crack extension steps and were used to control the crack growth simulation in a continuing loop until a desired crack length is obtained. 3.2. Crack driving force post-processing The post-processing operations for generating CDF values was carried out following the approach proposed by Lin and Smith. Blunted tip and spider web configurations have pairs of directly opposite nodes at the crack faces, thus allowing efficient calculation of the opening and sliding displacement between the faces. The CTOD estimation was based on the relative movement of these opposite nodes at a particular distance away from the crack tip. The selection of DCTOD over J within the elasto-plastic fracture mechanics (EPFMs) framework was due to the possibility of assessing crack fronts located closer to the material interface and calculation simplicity. J estimates in 2D problems at slightly greater distances from the tip to the material interface have required additional meshing features and complexity [19]. The eventually selected pair of nodes was positioned at a distance dm = 960 nm and dm = 900 nm away from the tip for the blunted and spider web models, respectively; this is schematically shown in Figs. 4 and 5. These locations were considered sufficiently close to crack tip to provide reliable estimates of the CTOD. Typical minimum element sizes for the blunted tip in the radial and circumferential direction were ebr = 300 nm and ebc = 26.1 nm, respectively; while for the spider-web elements, esr = 300 nm and esc = 82.1 nm. The element size in the present 3D analyses was slightly larger than that used in the earlier 2D analysis developed on the same material architecture [15]. 3.3. Crack front advance and remeshing The crack front advance was based on the growth around the crack front relative to a specified maximum extension step Damax, corresponding to the CDFmax, according to the relation [9]
da Dai ¼ dadN
i
dN max
Damax ¼
f ðCDFi Þ Damax f ðCDFmax Þ
ð1Þ
where Dai refers to the crack advance at location i, (da/dN)i to the expected growth rate at location i, and (da/dN)max to the maximum growth rate along the crack front. The use of small values of crack advance is recommended for problems with sudden crack front changes, deflections that cause out of plane crack growth and other factors that contribute to the formation of irregular crack surfaces. This condition also applies to changing or uneven loading conditions existing along
Fig. 4. CTOD measurement in blunted tip mesh configuration.
Fig. 5. CTOD measurement in spider-web mesh configuration.
20
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
Fig. 6. Relative crack extension along the front.
the material on which the crack front is embedded, such as bending where the material stress depends on its position to the neutral axis. The estimations of da/dN were based on a crack growth law in terms of the CTOD as the crack driving force as proposed by Tanaka et al. [20] for small cracks. This law was based on a study of ductile materials under cyclic loading and led to satisfactory results when compared to experimental data [20]. The form of the relationship proposed by Tanaka et al. is
da ¼ CðDCTODÞn dN
ð2Þ
where fitting parameters C and n account for the material behaviour. The n exponent values obtained for materials such as copper and steel range between 1.17 and 1.13. Crack growth laws of this form can also be found in the work of Tomkins [21] where the value n = 1 is adopted. The present analyses were only concerned with the proportional growth of the crack front according to the crack driving force estimations. Therefore, the selected value of C would not affect the prediction of the crack front growth process as shown when Eq. (1) is substituted in Eq. (2) resulting in
da Dai ¼ dadN
i
dN max
Damax ¼
DCTODni Damax DCTODnmax
ð3Þ
The value of n = 1.129 for the Al alloy lining was estimated through a combination of experimental data and numerical simulation [22] of crack growth in the lining layer of the tri-layered flat strip measured during interrupted three-point bending fatigue tests [16] under a loading ratio 0.1 using acetate strips to replicate the cracked surface. These tests were performed at a maximum strain range De = 0.0063 ± 0.0002, which corresponds to a maximum load P = 920 N. Re-meshing operations with this methodology only require the generation of a new crack front and transition elements. The new crack elements along the crack front were created according to the CDF estimates at the previous step. As the crack front length grows, the number of elements is automatically increased, as shown in Fig. 6, so that a fairly constant element size is maintained throughout the crack propagation analysis. 4. Comparison of crack tip modelling The generation of FE models for the study of fracture mechanics has used distinct crack tip mesh patterns that facilitate the computation of the crack driving force. Previous 2D work on multi-layered systems used refined blunted tips [15] for the study of shielding and anti-shielding effects under extensive elasto-plastic deformation. The development of 3D models with a similar element size was expected to result in considerable computational cost. Therefore, the creation of spider web crack tip meshes, adopting quarter-point and collapsed brick elements [23] for elastic and elasto-plastic analyses, respectively, was also considered. The suitability of blunted or initially sharp tips in spider web configurations was numerically tested as an initial validation step before the simulation of crack growth and the assessment of shielding and anti-shielding effects. This initial linear elastic analysis was carried out to determine the accuracy and numerical stability of the CTOD estimated from both models and its sensitivity to mesh refinement along the crack front; it was applied to a specimen containing an embedded crack, subjected to tension as shown schematically in Fig. 7. The dimensions of the analysed plate are the same as those of the flat strip specimen shown in Fig. 1 apart from its thickness which is taken equal to 2tL2 = 0.76 mm, for the purpose of assessing the effectiveness of a mesh representative of those subsequently used for analysing multi-layered systems of equal size. In this validation test, the surface crack is replaced by a penny crack and the free-of-traction conditions on the face of the lining are replaced by symmetry conditions. The embedded crack has a radius of 0.266tL2, is oriented perpendicular to the plane of
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
21
Fig. 7. Geometry of embedded crack under tension for numerical (left) and exact solution (right).
symmetry and subjected to a remote uniform tensile stress of 20 MPa. The material used for these analyses was assumed homogeneous and elastic, having the lining properties, that is, a Young’s modulus of 70 GPa and a Poisson’s ratio of 0.33. In order to assess the accuracy and numerical stability of the developed meshes a comparison with the exact solution for a flat ellipsoidal crack was obtained by Green and Sneddon [14]. According to this solution, the crack opening for a penny shaped flaw is given by
COD ¼
4r1 ð1 m2 Þ Ep
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dm ð2a dm Þ
ð4Þ
where a is the crack radius and dm the normal (radial) distance from the crack front. This solution and the problem’s axis symmetric nature indicate that for a given crack of radius a the crack opening at a distance dm from the tip would be constant irrespective of its position along the crack front. This refinement varied the number of elements from 18 to 54 providing an initial estimate of the computational cost expected from subsequent elasto-plastic analyses. Each elastic analysis using the ANSYS Sparse Direct Solver required from 1 to 3 h of computing time using four Intel Xeon processors with 8 GB of RAM memory analysing a model comprising between 60,000 and 120,000 nodes. The CTOD estimates along the crack front in both crack tip configurations showed an inherent variability generated by the numerical approximation of the model. Such variability is in contrast to a constant value yielded by the exact solution along the whole front. The standard deviation calculation of CTOD values along the crack front was developed for seven levels of refinement in both crack tip configurations and is shown in Table 2. The standard deviation obtained at the different levels of refinement showed that models using spider web crack tip configurations behaved more closely to the analytical solution and values up to 20 times smaller to those of blunted models were found. Mesh refinement also proved to reduce the standard deviation estimations, as expected. These studies also revealed that CTOD estimates obtained from spider web configurations are closer to the analytical solution, as expected. According to the mesh sensitivity analysis, the degree of accuracy achieved using the spider web crack tip configuration showed error estimates between 1% and 5%. The best compromise between accuracy en computational costs was achieved with 30 elements along the front and an element size of 5.85 lm with error estimates of 2% and only 1 h of simulation time.
Table 2 Standard deviation of normalised CTOD along the crack front. Number of elements along the front
Blunted configuration
Spider web configuration
54 46 38 30 26 22 18
2.3E03 3.6E03 4.0E03 4.8E03 5.1E03 6.6E03 7.2E03
4.1E04 5.2E04 8.2E04 2.0E03 1.4E03 3.2E03 6.5E03
22
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
5. Crack growth modelling The proportional growth along the crack front was modelled here according to Eq. (4). The crack growth study was carried out at tension and bending loading conditions with a maximum tensile load Pmax = 2684 N and a maximum bending load Pmax = 200 N and loading ratio of 0.1, which causes substantial yielding around the crack. The current loading and modelling conditions allow the analysis of crack fronts placed closely to the layers’ interface, where shielding and anti-shielding effects are most intense. The models developed to account for various crack sizes comprised between 60,000 and 1,20,000 elements with 1,40,000 to 2,40,000 nodes. This significant increase in the number of DOFs during the solution process arises from the need to keep the element size constant as the crack grows. The processing time for the analysis of each growth step required between 12 and 26 h and result files for both layered architectures that exceeded 50 GB of hard disc per loadcase. The individual computational cost for the two modelled architectures at each loadcase was comparable despite the greater propagation depths of 0.4 mm for the bilinear one in comparison to 0.36 mm for the tri-linear. Crack front proximity to the lining-interlayer interface in the tri-layer resulted in slow convergence rates due to the greater compliance of the interlayer. The crack front evolution in bi-layer and tri-layer architectures was almost identical for crack depths less than 0.27 mm. Small differences could be observed at greater depths, which were consistent with the CDF trends in respective two-dimensional analyses but yielded almost identical outcomes until the crack tip came sufficiently close to the layers interface. The difference between the predicted crack fronts in tri-layer and bi-layer became obvious at depths greater than 0.3 mm.
5.1. Crack driving force evolution and shielding The evolution of the estimated values of CTOD and their pattern along the crack front in both layered systems subjected to bending are shown in Fig. 8. The similarities of the two plots may be noted. The CTOD pattern had maximum values at the surface while the minimum was located at the deepest point of the crack front as observed by previous authors when bending stresses are applied [10,17,24]. The sampling points at which the CTOD was estimated along the front were positioned at equal distances; in Fig. 9 the front contour coordinate has been normalised with respect to the length of the front. The evolution of CTOD with crack depth at two key locations along the crack front, for both layered systems, is shown in Fig. 9. In both cases, the CTOD estimates at the surface showed very similar values and rising trends. At the deepest point, they were almost identical for crack depths less than, approximately, 0.31 mm; they were also rising but at a lower rate than at the surface. The CTOD values for the tri-layer and bi-layer models became significantly different as the crack fronts approached the respective interfaces. Shielding and anti-shielding effects were therefore observed at the deepest point of the crack front at depths greater than 0.31 mm (82% of the lining thickness). Crack growth beyond this point showed accelerating and decelerating trends according to the adjacent layer stiffness. However, it should be noted that these trends did not affect to a great extent the magnitude of the estimated CTOD. The disproportional growth along the front, shown through DCTOD in Fig. 9, can be better represented by the CDF ratio of points A and B, which respectively represent the points located nearest and furthest to the materials interface. Fig. 10 demonstrates more clearly the evolution of these two parameters for the bending load case. It is thus possible to identify three different growth phases: initial transition from semi-circular to quasi semi-elliptical (I), crack growth under combined
(a)
(b)
Fig. 8. Crack front evolution in: (a) tri-layer and (b) bi-layer architectures under P = 200 N.
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
23
Fig. 9. Crack driving force evolution in bi-layer and tri-layer architectures subjected to three point bending.
Fig. 10. Ratio of CDF obtained at the deepest point (CTODd) to that at the surface (CTODs) in bi-layer and tri-layer architectures subjected to three-point bending.
tension-bending conditions (II) and the start of shielding and anti-shielding mechanisms (III). The initial transition phase corresponds to the changes that the crack front must undergo to achieve a more ‘‘natural’’ shape consistent with the loading conditions. The rate of change depended on the difference between the current and ‘‘natural’’ shape. The second phase in the plot of Fig. 10 shows the typical results of crack growth in a monolithic specimen subjected to tension-bending with maximum CDF values at the free surface where the tensile bending stress also attains its maximum value. As the crack front grew progressively deeper into the strip, lower CDF ratios were estimated. This behaviour has been previously observed [10] in the crack growth in monolithic strips and it is associated with significant change of the crack aspect ratio as the crack front approaches the neutral axis of the strip. However, the multi-layer architecture produced a different stress distribution to that observed in a monolithic strip. The stress distribution through the lining thickness of an undamaged specimen is equivalent to a combined uniform tension-bending loading on a monolithic strip. Under a load P = 200 N, the obtained stress distribution was equivalent to a uniform tension of about 53.5 MPa and maximum bending stress of about 8.5 MPa, thus uniform tension was dominant. The ratio between uniform and bending tensile stresses is a key indicator of how stress and deformation are distributed along the crack front; however, it is important to consider that these conditions may vary as the crack extends and plasticity spreads throughout the material. The third phase in Fig. 10 shows clearly the appearance of shielding and anti-shielding. Shielding is associated with increasing rate at which CTODB/CTODA, drops in the bi-layer architecture. In contrast, due to anti-shielding, the rate at which CTODB/CTODA decreases remains almost constant. If the numerical simulation had allowed the crack to grow further, closer to the layers’ interface, the anti-shielding effect may have become more obvious. However, the small separation between the
24
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
front and the interface could not allow a significant shape change before the crack penetrates into the interlayer. Pursuing the crack growth within the softer interlayer, as in previous two-dimensional analyses [15], would have elucidated the local behaviour further. This however was beyond the available computational capacity required by the adopted model but offers prospects for future interesting work. The 3D modelling challenges of a front placed within two different materials and the computational challenges arising from even larger plastic deformations in the softer interlayer point to the complexity of such a simulation. The crack growth analyses applied to the bi-layer architecture produced sudden drops and rises in the evolution of CDF when the crack depth exceeded 0.34 mm (90% of the lining thickness), as shown in the enlarged view of the inset in Fig. 9a (Points 1, 2 and 3). The existence of these gradient discontinuities may be due to the extension steps used here not being sufficiently small to capture the smooth changes of CDF. It should be noted that a front shape is produced here as a consequence of combining the shielding effect with that of the overall bending deformation state governing the crack front extension. Despite the existence of this local non-smoothness, the overall shielding trend is clear. Therefore, model refinement to achieve further crack extension was not considered necessary. On the other hand, crack growth under tensile loading showed similar trends to bending in areas such as the existence of three evolution areas and overall shielding and anti-shielding behaviour, shown in Fig. 11. In this particular case, the CDF ratio between points A and B indicated more clearly how the existence of a more compliant material would increase the CDF ratio between A and B as the crack tip approaches to the interface while a stiffer material would cause an opposite effect. These changes become more evident as the distance between the material interface and the crack tip is reduced. The gradual shift towards greater or smaller CDF ratios demonstrates the need of developing consecutive elasto-plastic analysis for a reliable estimation of CDF.
1.25 1.2 1.15 1.1 1.05 1 0.95 0.9 0.85 0.8
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
Fig. 11. Ratio of CDF obtained at the deepest point (CTODd) to that at the surface (CTODs) in bi-layer and tri-layer architectures subjected to pure tension.
-4
x 10 4
3
Bi-layer 3D deepest point Tri-layer 3D deepest point Bi-layer 2D Tri-layer 2D
2
1 0.05
0.1
0.15
0.2
0.25
Fig. 12. Comparison between CDF estimates obtained from 2D and 3D models.
0.3
0.35
0.4
0.45
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
25
5.2. Comparison to 2D analyses Shielding or amplification arising from crack growth through dissimilar materials has been studied in the past through two-dimensional analyses, Despite their simplified representation of the actual problem, two-dimensional models may provide a reliable initial assessment before the development of more rigorous but also more complicated three-dimensional modelling or a good approximation to the correct answer when a quick output is required. It is evident from Fig. 12 that the two-dimensional simplification of the three-point bending simulation reduces the stiffness of the specimen more dramatically by modelling a width-through crack instead of half-penny type crack. The CDF estimates from 2D and 3D analyses of the flat strip subjected to a load P = 200 N were clearly different, as shown in Fig. 12b. The higher CDF estimates from the 2D model are expected since a through crack along the strip width is assumed. The acceleration and deceleration effects on crack growth caused by shielding and anti-shielding in 3D models did not appear as strong as predicted by two-dimensional analyses. This difference can be attributed to the more constrained growth of the surface crack shape as it expands as a whole. Shielding and anti-shielding are localised phenomena that only affect small segments of the crack front. The small separation between the crack front and the layers’ interface would affect the CDF values locally causing disproportional growth that modifies the surface crack profile. At the same time, the neighbouring material discourages the disproportional growth of the most affected areas. In other words, shielding and anti-shielding effects are averaged or diluted into a greater portion of material based on the crack front growth in its entirety. 5.3. Plastic zone The plastic zone around the crack tip has been extensively studied in test specimens with width-through cracks. Tri-axiality effects have been of interest in these specimens, often accounted for by assuming plane strain or plane stress conditions according to the location of the crack tip in the specimen. The shape of the plastic zone around a half penny-type crack is more complicated, especially when the loading conditions differ from uniform remote tension normal to the crack surface within an infinite solid. Fig. 13 shows typical von Mises plastic strain contours around the crack front for the bi-layer architecture subjected to bending and how shielding occurs by a stiffer and stronger backing. This condition can be observed where the crack front is the closest to the interface through the development of a broader plastic zone in this region with higher levels of plasticity. It is expected that the aforementioned broadening may facilitate the crack initiation process in the backing layer by increasing the stress concentration area within the backing. In contrast, as the crack front approaches the interlayer, in the tri-layer architecture, high values of plastic strain can be observed in this intermediate layer. For clarity, only deformations above 5e4 are shown since plastic deformation covers around 20% of the top surface. The lower yield strength of the interlayer produced a protuberance in the plastic zone shape as shown in Fig. 14. This protuberance expands through the interlayer forming a secondary plastic zone that is expected to promote faster crack propagation.
Fig. 13. Plastic zone broadening at the deepest point of the crack due to backing shielding in bi-layer architecture (Von Mises plastic strain contour).
Fig. 14. Plastic zone shape and low strain area next to the crack in the tri-layer architecture (Von Mises plastic strain contour showing only strains greater than 5e4 for clarity).
26
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
6. Discussion and conclusions The adopted methodology put emphasis on the evolution of the crack front as it approaches a dissimilar layer to assess shielding and anti-shielding effects on half-penny cracks. Modelling crack growth in systems comprising dissimilar layers that exhibit both elastic and plastic mismatch made the analysis complex and time consuming. The shielding and anti-shielding effects were predicted to be less significant in three-dimensional crack growth in comparison to those arising from width-through (or through-thickness) cracks in two-dimensional specimens. The crack front segments affected most intensively by shielding and anti-shielding in 3D are expected to accelerate or decelerate, respectively, crack growth in the neighbouring material. At the same time, the neighbouring material discourages the disproportional growth of the parts along the crack front most affected by shielding or anti-shielding. Crack growth processes are defined by the local conditions of material and loading around the crack front. Local loading and material conditions may change as the front extends leading to shapes that cannot be forecasted a priori, in contrast to homogeneous materials subjected to typical loading cases [10,17,24] and their local loading profile on the front. The study of shielding and anti-shielding in this work displayed loading and material conditions that vary significantly as the crack extends towards a dissimilar layer. At the same time, the elasto-plastic nature of the multi-layered architecture displays a stiffness gradient from the crack front to the interface which will change according to the stress concentration conditions, loading, and stress strain curve of the involved materials. An appropriate estimation of crack driving forces requires updated data regarding CDF and the state of the material which leaded to computationally expensive analysis and sophisticated modelling routines. The effects of shielding and anti-shielding were studied in layered architectures subjected to bending and tension. The crack growth under bending loading revealed localised speed reductions as the front approached the neutral axis. Additionally, shielding and anti-shielding effects caused less significant modifications to the front profile. The study of such layered systems under tensile loading revealed more clearly accelerations or reductions of CDF as the front approached a dissimilar material. Occurrences of crack deflection and bifurcation observed experimentally still remain to be investigated numerically; such phenomena lead to interesting crack development, especially through the interlayer. The simulation of surface crack extension into the interlayer offers great modelling challenges: from the prediction of three-dimensional deflected and bifurcated forms to the high level of mesh refinement necessitated by the thinness of the interlayer. The development of such analyses with elastic materials would significantly simplify the task; it is however important, from a practical point of view, to identify deflections and bifurcations using elasto-plastic material models, as in previous two-dimensional studies [15].
Acknowledgments The financial support of the Mexican Council of Science and Technology and the School of Engineering Sciences is gratefully acknowledged. The training provided by and discussions with MAHLE Engine Systems was also invaluable in carrying out this research.
References [1] Sugimura Y et al. Fracture normal to a bimaterial interface: Effects of plasticity on crack-tip shielding and amplification. Acta Metall Mater 1995;43(3):1157–69. [2] Sugimura Y, Grondin L, Suresh S. Fatigue crack growth at arbitrary angles to bimaterial interfaces. Scripta Metall Mater 1995;33(12):2007–12. [3] Suresh S, Sugimura Y, Tschegg EK. Growth of a fatigue crack approaching a perpendicularly oriented bimaterial interface. Scripta Metall Mater 1992;27:1189–94. [4] Pippan R, Flechsig K, Riemelmoser FO. Fatigue crack propagation behavior in the vicinity of an interface between materials with different yield stresses. Mater Sci Engng, A 2000;283(1–2):225–33. [5] Wang B, Siegmund T. Simulation of fatigue crack growth at plastically mismatched bi-material interfaces. Int J Plast 2006;22(9):1586–609. [6] Groh S et al. Fatigue crack growth from a cracked elastic particle into a ductile matrix. Phil Mag 2008;88(30):3565–83. [7] Newman Jr JC, Raju IS. An empirical stress-intensity factor equation for the surface crack. Engng Fract Mech 1981;15(1-2):185–92. [8] Raju IS, Newman JC. Stress-intensity factors for a wide range of semi-elliptical surface cracks in finite-thickness plates. Engng Fract Mech 1979;11:817–29. [9] Riddell WT, Ingraffea AR, Wawrzynek PA. Experimental observations and numerical predictions of three-dimensional fatigue crack propagation. Engng Fract Mech 1997;58(4):293–310. [10] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates. Part II: crack shape change. Engng Fract Mech 1999;63(5):523–40. [11] Carter BJ, Wawrzynek PA, Ingraffea AR. Automated 3-D crack growth simulation. Int J Numer Method Engng 2000;47(1–3):229–53. [12] Schollmann M, Fulland M, Richard HA. Development of a new software for adaptive crack growth simulations in 3D structures. Engng Fract Mech 2003;70(2):249–68. [13] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engng Fract Mech 2008;75(5):943–60. [14] Green AE, Sneddon IN. The distribution of stress in the neighbourhood of a flat elliptical crack in an elastic solid. Math Proc Cambridge Philos Soc 1950;46(1):159–63. [15] Burke-Veliz A, Reed PAS, Syngellakis S. A numerical study of crack shielding and deflection under extensive plasticity. Engng Fract Mech 2009;76(9):1345–56. [16] Ali M, Syngellakis S, Reed PAS. A comparison of fatigue performance of HVOF spray coated and conventional roll bonded aluminium bearing alloys. Mater Sci Technol 2008;25(5):575–81.
A. Burke-Veliz et al. / Engineering Fracture Mechanics 88 (2012) 15–27
27
[17] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates. Part I: the numerical technique. Engng Fract Mech 1999;63(5):503–22. [18] Ansys, Ansys reference manual; 2005. [19] Kolednik O et al. On the fracture behavior of inhomogeneous materials – a case study for elastically inhomogeneous bimaterials. Int J Solids Struct 2005;42(2):605–20. [20] Tanaka K, Hoshide T, Sakai N. Mechanics of fatigue crack propagation by crack tip plastic blunting. Engng Fract Mech 1984;19:805–25. [21] Tomkins B. Fatigue crack propagation: an analysis. Phil Mag 1968;18(155):1041–66. [22] Burke-Veliz, A., Finite element modelling of crack growth in multi-layered architectures. In: Computational engineering and design group. PhD Thesis. Southatmpton: University of Southampton; 2009. [23] Barsoum R. Triangular quarter point elements as elastic and perfectly plastic crack tip elements. Int J Numer Meth Engng 1977;11:85–98. [24] Lin XB, Smith RA. Finite element modelling of fatigue crack growth of surface cracked plates. Part III: stress intensity factor and fatigue crack growth life. Engng Fract Mech 1999;63(5):541–56.