Microsoft Word - numero_32_art_1 N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 1 Influence of crack offset distance on the interaction of multiple cracks on the same side in a rectangular plate Neeraj Bisht, P. C. Gope, Kuldeep Panwar College of Technology, Govind Ballabh Pant University of Agriculture and Technology, Pantnagar (263145), Uttarakhand, India. neerajbisht30@gmail.com ABSTRACT. In the present work finite element method has been employed to study the interaction of multiple cracks in a finite rectangular plate of unit thickness with cracks on the same side under uniaxial loading conditions. The variation of the stress intensity factor and stress distribution around the crack tip with crack offset distance has been studied. Due to the presence of a neighbouring crack, two types of interactions viz. intensification and shielding effect have been observed. The interaction between the cracks is seen to be dependent on the crack offset distance. It is seen that the presence of a neighbouring crack results in the appearance of mode II stress intensity factor which was otherwise absent for a single edge crack. It can be said that the proximity of cracks is non-desirable for structural integrity. The von-Mises stress for different crack orientations has been computed. Linear elastic analysis of state of stress around the crack tip has also been done. KEYWORDS. Finite element method; Crack interaction; Von-Mises stress. INTRODUCTION he fracture mechanics theory can be used to analyse structures and machine components with cracks and to obtain an efficient design. The basic principles of fracture mechanics developed from studies of [1-3] are based on the concepts of linear elasticity. The interaction between multiple cracks has a major influence on crack growth behaviours. This influence is particularly significant in stress corrosion cracking (SCC), welding, riveting etc. because of the relatively large number of cracks initiated due to environmental effects. Pseudo – traction –electric – displacement –magnetic –induction method has been proposed [4] to solve the multiple crack interaction problems in the magneto elastic material. Most of the real life situations have the problem of multiple cracks and so it becomes imperative to study this interaction for an array of cracks, keeping this in mind interaction between two parallel cracks has been studied and a detailed analysis has been done in this regard. Since today, there have been over 20 approaches to calculate stress intensity factors. Some of these are the integral transform method [5], the Westergaard method [6],the complex variable function method [7], the singular equation integral method [8], conformal mapping [9], the Laurent series expansion [10], boundary collocation method [11], Green’s function method [12], the continuous distribution dislocation method [13], the finite element method [14], the boundary element method [15], the body force method [16] and the displacement discontinuity method [17]. The solutions of many of the fracture mechanics problems have been compiled in data hand books for stress intensity factors [18] and [19]. The configuration of multiple cracks is so complicated that a solution may not be available from the handbooks and literatures. The above mentioned methods with analytical features, which are usually suitable for special cases or very simple crack configurations, are not sufficient to obtain reasonable results for general orientations due to the multiple T N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 2 restrictions. In these cases numerical approaches are usually employed. In the numerical approaches proposed so far finite element method provides a very simple, effective and accurate technique for evaluation of fracture parameters. The historical development of computational fracture mechanics is found in the works of Ingraffea et al. [20] and Sinclair [21]. Sinclair presented an extensive review of numerical prediction models to determine stress intensity factors. The advantages and disadvantages of using finite element in computational fracture mechanics have been well addressed by Ingraffea [22]. The aspect of mesh refinement and associated error in computing stress intensity factors using finite element method has also been studied by Miranda et al. [23]. It has been reported that excessive mesh refinement may significantly degrade the calculation accuracy in crack problems. They also pointed out that the ratio between the longest and shortest element edge lengths should be kept below 1600 to avoid calculation errors in SIF calculations. For meshes with length ratios higher than 1600, improved numerical methods to deal with ill conditioned matrices would be necessary to not compromise the calculation accuracy of the calculated SIF. Many works on mesh generation algorithms and new methods to improve the numerical computation of SIF values have been found in the works of Miranda et al. [24, 25]. Recent studies have also shown that the coefficients of higher order terms can also play an important role in the fracture process in notched or cracked structures. It has been observed that in addition to the singular term, the higher order terms, in particular, the first non-singular stress term ( known as the T stress) may also have significant effects on the near notch tip stress field. The T-stress is considered in some studies as an auxiliary parameter for increasing the accuracy of the results for KI. Kim et al. [26], for instance, showed that this non-singular term has noticeable effects on the size and shape of plastic zone near the notch tip. It has been demonstrated that the first non-singular term may have considerable contributions to the stress components around the notch tip and also on the fracture resistance of notched components under mode I loading [27-29]. FINITE ELEMENT MODELLING he numerical simulations were run by means of the finite element (FE) software ANSYS to determine the stress intensity factors of two edge cracks on the same side of the specimen. The specimen is schematized by a 2D model. The specimen thickness in FE analysis was kept 1.0 mm. The other dimensions are length L= 200 mm, width W= 80 mm and the model was studied in plane strain condition, the specimen geometry and detailed dimensions are shown in Fig 1(a). The simulations have been run on full model. The stresses are applied at the two extremes of the specimen in the direction perpendicular to the crack plane (Fig. 1(b)). Figure 1(a): Specimen Geometry. Figure 1(b): Specimen with boundary conditions. Isoparametric quadrilateral element (PLANE 82) having 8 nodes with singularity elements at and around the crack tip has been used throughout the analysis and is shown in Fig. 2. The singularity phenomenon at the vicinity of the crack tips has been addressed by applying triangular element option of the PLANE 82 element. The radius of second row of elements is taken as a/8, where a is the half crack length and the radius ratio (second row/first row) is adjusted T N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 3 automatically. The number of elements around the circumference is taken as 32. The FE modelling parameters are optimized on the basis of the error analysis presented in Fig. 3 for two edge cracks. Figure 2: Plane 82 element with 8-nodes. Figure 3(a): Variation of stress intensity factor with number of crack tip elements. Figure 3(b): Variation of stress intensity factor with radius of first row of elements, KI,FEM and KI, THE are finite element and theoretical based computations of mode I stress intensity factor The error analysis has been done by varying the radius of the first row of the crack tip element and number of elements in the first row and computing the stress intensity factor for various parameters. The density of the FE mesh is modified by varying the number of the elements of the first row as 16, 20, 24, 32 and 40, keeping the radius of the first row as a/8 where a has been taken as 10 mm. Also, the radius of the first row (a/n) around the crack tip is varied, taking n= 8, 10, 12, 16 and 20 with 32 number of elements. For comparison KI is calculated theoretically from the relation given for two collinear edge cracks [30]. It is found that the KI calculation errors stay below 0.8% for all meshing strategies. Fig. 3(a-b) show variation of the normalized KI (FE based computation to the actual (theoretical) Mode I stress intensity factor) for different mesh configurations. It can be observed from Fig. 3 that, on average, the calculation with the radius of first element as 1.25 mm (i.e. a / 8=1.25) and number of elements around the crack tip as 32 yields least error i.e. closest to unity amongst the other meshing strategies. It is mentioned by Miranda et al. [32] that calculation error and associated standard deviation tend to increase for higher mesh density. These increasing errors are a result of an ill conditioned numerical problem [32]. The increasing error due to the decrease in first crack tip radius may be due to the stress singularity that is present at the crack tip and the linear elastic fracture mechanics concept fails to predict the stress intensity factor. The other FE modelling parameters are shown in Fig. 4 and 5. In Fig. 5 node 1 is the crack tip node and the nodes 2, 3, 4 and 5 are used to represent the crack path for evaluating the fracture parameters. The nodes are essentially taken in this order when a full crack model is employed. For the numerical simulation, a uniform pressure intensity of 1.0 MPa (80 N) is applied to the upper and lower edges in the vertical direction (y axis). The material properties are Young’s Modulus E=70 GPa, Poisson’s ratio=0.33. The mode I and mode II stress intensity factors are computed from the following relations using the displacement extrapolation method. N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 4 Figure 4: Crack tip meshing and different parameters. Figure 5: Nodes used for defining crack path. The analysis uses a fit of the nodal displacements in the vicinity of the crack. The actual displacements at and near a crack for linear elastic materials are given by Paris [31] as:     3 3 2 1 cos cos   2 3 sin sin 4 2 2 2 4 2 2 2 I IIK Kr ru k k G G                       (1)    θ 3 32 1 sin sin   2 3 cos cos 4 2 2 2 4 2 2 2 I IIK Kr rv k k G G                      (2) 2  sin   2 2 IIIK rw G    (3) where: u, v, w = displacements in a local Cartesian coordinate system, shown in Fig.6. r, θ = coordinates in a local cylindrical coordinate system, shown in Fig.6. G = shear modulus KI, KII, KIII = stress intensity factors relating to deformation shapes, shown in Fig.6. ν = Poisson’s ratio N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 5 Figure 6: Representation of crack with coordinate system. Evaluating Eq. (1-3) at 0180   and dropping the higher order terms Eq. (1-3) yields:  II K r u 1 k 2G 2π   (4)   2 2 IK rv k G   (5) 2     2 IIIK rw G   (6) For full crack models Eq. (4-6) can be reorganized to   | |   2I vG K k r    (7)   | | 2 1 II uG K k r     (8)   | |   2III w K G r    (9) where Δv, Δu and Δw are the motions of one crack face with respect to the other. k= 3−4ν, for plane strain or axisymmetric; (3−ν)/ (1+ ν) for plane stress; where ν is Poisson's ratio. The final factor    v r  is evaluated based on the nodal displacements and locations. For practical purposes the value of    v r  is approximated by limiting the value of    v r  by simply evaluating the following expression for a small fixed value of r (small in relation to the size of the crack) as: | |   v A Br r   (10) At point I shown in Fig.7, v=0 Hence, in the limiting condition N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 6 0 | |   lim r v A r  (11) Figure 7: Representation of nodes on crack face used for crack analysis. The Eq. (7) becomes 2   2  I GA K k  (12) Similarly | |   u C Dr r   (13) At point I of Fig. 7 u=0 Hence, in the limiting condition 0 | | lim r u C r  (14) and Eq. (8) becomes,   2 1 II GC K k   (15) Similarly for evaluating KIII, using | |w E Fr r   (16) and at point I of Fig. 7 w=0 and Eq. 9 in the limiting condition becomes, 0 | | lim   r w E r  (17) Thus, Eq. 9 becomes   2IIIK GE (18) N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 7 Mode III stress intensity factor in the present investigation is not considered because the thickness of the plate in FE analysis is taken as unity. However, the three dimensional effect or plate thickness effect on the stress intensity factors can be seen in the work of Kotousov et al. [33, 34]. In the present FE analysis, out of balance convergence and degree of freedom increment convergence criteria have been applied. The tolerance level has been taken as 0.001. These criteria are well documented in the ANSYS manual [35]. RESULTS AND DISCUSSION Stress Intensity Factor he stress intensity factors computed were normalized by K0 given for single edge cracked plate under uniform tension [30].  0  aK a f W  (19) where σ is the applied stress, a is the crack length and  af W is the geometric correction factor given as:   2 3 4 1.12 0.23 10.56 21.74 30.42 a a a aaf W W W W W                            (20) In the present investigation, a = 10 mm, W = 80 mm. Thus geometric correction factor becomes  af w  1.13 (21) and K0 = 400.57 MPa mm Fig. 8 shows the effect of crack offset distance H on the normalized stress intensity factors. From Fig. 8 it can be seen that normalized mode I SIF (KI/K0) for H=0.5 mm reduces drastically to about 66% value as compared to single edge crack. As the crack offset distance increases beyond H=0.5 mm the normalised mode I stress intensity factor starts increasing, and its value becomes 95% of single edge crack when H becomes 20 mm. So, it can be said that there is a shielding effect due to proximity of cracks due to which mode I SIF decreases as compared to single edge crack, however this shielding effect ceases to exist when the cracks move farther away. Figure 8: Variation of normalized stress intensity factors with crack offset distance H. The variation of normalized mode II SIF KII/K0 shown in Fig.8 is just opposite to that of mode I. Fig. 8 shows that when normalised mode I SIF attains the minimum value, associated mode II attains the maximum value which is almost 19% of mode I SIF for single edge crack. This behaviour is seen at H= 0.5 mm. Thereafter mode II SIF decreases with increase T N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 8 of crack offset distance H. There is an amplification effect as far as mode II SIF is concerned which also becomes insignificant for remotely placed cracks. This indicates that when two edge cracks are extremely close to each other (H ≤ 0.5 mm), shielding effect occurs and the normal stress components σxx and σyy reduces causing the shear component of the stress to increase. As a result, mode I SIF decreases and mode II SIF increases. Now, as the crack offset distance H increases from H=0.5 mm to higher value, mode I SIF increases and mode II SIF decreases. It is found that for about H ≥ 20 mm, the crack interaction becomes negligible and both cracks behave as a single crack. The variation of normalised KII/KI with H is shown in Fig. 9. This figure reveals that when the cracks are very close to each other mode II SIF plays a significant role in crack growth and contributes about 27% of mode I SIF. Beyond H= 20 mm, the contribution of mode II SIF remains below 2% only. These results indicate that special care should be taken while predicting growth from stress intensity factor, particularly when two or more edge cracks are very close to each other. Figure 9: Variation of KII/KI with crack offset distance H. Von-Mises Stresses Distribution of von-Mises stresses around the crack tip for various values of H at uniform stress σ = 1 MPa are shown in Fig. 10-11. von-Mises stresses are computed taking midpoint or element stresses. Identical boundary conditions, crack tip element size and full mesh model have been used for all crack configurations. The distribution of equivalent von Mises stresses for two edge cracks configuration for H = 2 mm and H = 20 mm have been shown in Fig.10 and 11. From Fig. 10 and 11 it can be observed that with increasing crack offset distance H the equivalent von Mises stresses increases. This indicates that crack offset distance have greater affect on the state of stress ahead of the crack tip. Analysis of the State of Stress around the crack tip The two dimensional finite element model presented in this work has been used to investigate the state of stress around the crack tip in plane strain condition. The state of stress at a radial distance r from the crack tip is schematically shown in Fig. 12. The results obtained for different crack configurations are presented in Fig. 13 to 15. The stresses computed for double edge cracks are normalized with σy0 which is the FE solution for single edge crack corresponding to θ=00 and r=0.5 mm along the crack plane. The variation of normalised σxx, σyy and τxy (normalised by σyy for single edge crack= 731.57 MPa) for two edge cracks on the same side of a rectangular plate and separated by offset distance H= 0, 2, 10 and 16 mm are shown in Fig. 13 – 15. The variation of these stresses around the crack tip at a radius of 0.5 mm from the crack tip is shown for the assumed plane strain condition. Identical loading, boundary conditions and mesh arrangements were used for all crack offset conditions. The stresses are taken at the midpoint of each crack tip element. It has been observed that symmetrical N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 9 distribution of σyy occurred on the upper and lower side of the crack, whereas σxx and τxy show asymmetrical distribution around the crack tip. Figure 10: Von Mises stress distribution for two offset edge crack geometry for H=2 mm. Figure 11: Von Mises stress distribution for two offset edge crack geometry for H=20. Figure 12: Schematic representation of state of stress ahead of the crack tip, r and θ are the polar coordinates. The σxx component is about 0.1 to 0.75 times of σyy at the crack tip element. Fig. 13 and 14 also reveal that σxx and σyy are lower in magnitude compared to a single edge crack. The normal stress component σxx increases with increasing offset distance, but its magnitude remains less than a single edge crack. It is also seen that σxx is least when two cracks are very close to each other (H ≤ 2 mm). It is well established that mode I SIF mainly depends upon variation of σyy at the crack tip. The present analysis shows that σyy is higher for single edge crack as compared to multiple cracks with different offset distance H. Hence, in Fig. 9 KI for H = 0 mm is found to be higher as compared to other multiple crack configurations. It is also seen in Fig.13 that σyy is minimum for H = 2 mm. Thus mode I SIF shown in Fig. 8 is found to attain minimum value at H = 2 mm, and thereafter it increases but remains less than a single edge crack. The variation of τxy for different crack offset distance H = 0, 2, 10, 16 mm is shown in Fig.15. The figure shows that shearing stress is zero for single edge crack (H = 0 mm) and maximum for H = 2 mm. It is seen that as H increases, the shearing stress decreases. Thus, mode II SIF which is mainly due to shear component attains maximum value at H = 2 mm and thereafter it reduces as H increases. The results shown in fig. 8 also confirms that for single edge crack (H = 0 mm), mode II SIF is zero. This is also in conformity with the theoretical results for single edge crack. N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 10 Figure 13: Variation of σxx/σy with H for specimen geometry A1. Figure 14: Variation of σyy/σy with H for specimen geometry A1. Figure 15: Variation of τxy/σy with H for specimen geometry A1. CONCLUSIONS 1. There is a shielding effect on mode I SIF due to the presence of a neighbouring crack. 2. The shielding becomes more pronounced as the cracks move towards each other. N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 11 3. For mode II SIF there is an amplification effect. Mode II SIF which is otherwise absent for single cracks develops due to the presence of neighbouring cracks. 4. The amplification is more significant for closer cracks. 5. Both the amplification and shielding effects cease to exist as cracks move farther away from each other. REFERENCES [1] Inglis, C.E., Stresses in a plate due to the presence of cracks and sharp corners, Transactions-Institute of Naval Architect, 55 (1913) 219-230. [2] Griffith, A.A., The phenomenon of rupture and flow in solids, Trans. Royal Soci. London, 221 (1920) 163-198. [3] Westgaard, H. M., Bearing pressure and cracks, J. Appl. Maths Mech., 6 (1939) 49-53. [4] Wen-Ye, T., Gabbert, U., Multiple crack interaction problems in magnetoelectrostatic solids, Europ. J. Mech. A/Solids, 23 (2004) 599. [5] Sneddon, I. N., Lowengrub, M., Crack Problems in Classical Theory of Elasticity, 1st Edn., Wiley, New York, (1969). [6] Westergaard, H. M., Bearing pressure and cracks, ASME J. Appl. Mech., 6 (1939) 49-53. [7] Muskhelishvili, N. I., Some basic problems of mathematical theory of elasticity, 4th Edn. Holland, Noordhoff, (1977). [8] Erdogan, F., Stress Intensity Factors, ASME J. Appl. Mech., 50 (1983) 992-1002. [9] Bowie, O. L., Analysis of an infinite plate containing radial cracks originating at the boundary of an internal circular hole, J. Math. Phys., 25 (1956) 60-71. [10] Isida, M., Stress intensity factors for the tension of an eccentrically cracked strip, ASME J. Appl. Mech., 33 (1966) 674-675. [11] Williams, M. L., On the stress distribution at the base of a stationary crack, ASME J Appl. Mech., 24 (1957) 104-114. [12] Sih, G. C., Paris, P. C., Erdogan, F., Crack tip stress intensity factors for plane extension and plate bending problems, ASME J. Appl. Mech., 29 (1962) 306-314. [13] Eshelby, J. D., Franck, F. C., Nabarro, F. R. N., The equilibrium of linear arrays of dislocations, Phil. Mag., 42 (1951) 351-364. [14] Watwood, V. B., The finite element method for prediction of crack behaviour, Nucl. Engg. Des., 11 (1969) 323-332. [15] Cruse, T. A., Numerical solutions in three dimensional elastostatics, Int. J. Solids Struct., 5 (1969) 1259-1274. [16] Nisitani, H., Solutions of notch problems by the body force method in fracture mechanics, Edited by G. C. Sih, fifth ed. Noordhoff International, (1978) 1-68. [17] Crouch, S. L., Starfield, A. M., Boundary element method in solid mechanics, with application in rock mechanics and geological mechanics, 1st Edn., George Allen & Unwin, (1983). [18] Pook, L. P., Stress intensity factor expressions for regular crack arrays in pressurised cylinders, Fatig. Fract. Eng. Mater. Struct., 13 (1990) 135–143. [19] Rooke, D. P., Cartwright, D. J., Compendium of stress intensity factors Great Britain, Ministry of Defence, Procurement Executive, (1976). [20] Ingraffea, A.R., Wawrzynek, P.A., Comprehensive structural integrity, 1st Edn., R.d.B.a.H. Mang (Ed.), Elsevier Science Ltd., Oxford, (2003). [21] Sinclair, G., Stress singularities in classical elasticity—II: Asymptotic identification, App. Mech. Rev., 57 (2004) 251– 298. [22] Ingraffea, A.R., Encyclopaedia of computational mechanics, John Wiley and Sons, (2004). [23] de Oliveira Miranda, A. C., Meggiolaro, M. A., Martha, L. F., de Castro, J. T. P., Stress intensity factor predictions: Comparison and roundoff error, Comput. Mater. Sci., 53 (2012) 354–358. [24] Miranda, A.C.O., Meggiolaro, M.A., Castro, J.T.P., Martha, L.F., Bittencourt, T.N., Fatigue life and crack path predictions in generic 2D structural components, Eng. Fract. Mech., 70 (2003) 1259–1279. [25] de Oliveira Miranda, A.C., Meggiolaro, M. A., de Castro, J. T. P., Martha, L. F., Fatigue life prediction of complex 2D components under mixed-mode variable amplitude loading, Int. J. Fatig., 25 (2003) 1157–1167. [26] Kim, J. K., Cho, S. B., Effect of second non-singular term of mode I near the tip of a V notched crack, Fatig. Fract. Eng. Mater. Struct., 2 (2009) 346–356. [27] Ayatollahi, M. R., Nejati, M., Determination of NSIFs and coefficients of higher order terms for sharp notches using finite element method, Int. J. Mech. Sci., 53 (2011) 164-177. [28] Ayatollahi, M. R., Dehghany, M., On T-stresses near V-notches, Int. J. Fract., 165 (2010) 121-126. N. Bisht et alii, Frattura ed Integrità Strutturale, 32 (2015) 1-12; DOI: 10.3221/IGF-ESIS.32.01 12 [29] Ayatollahi, M. R., Nejati, M., Experimental evaluation of stress field around the sharp notches using photoelasticity, Mater. Des., 32 (2011) 561-569. [30] Prashant, K., Elements of Fracture Mechanics, First Ed., McGraw-Hill Book Company, New-Delhi, (2009). [31] Paris, P. C., Sih, G. C., In ASTM Spec. Tech. Publ., STP 381, ASTM Philadelphia, STP 381 (1965) 30-83. [32] Rice, J. R., Matrix Computations and Mathematical software, 1st Edn., McGraw-Hill Computer Science Series, McGraw-Hill, (1981). [33] Kotousov, A., Berto, F., Lazzarin, P., Pegorin, F. Three dimensional finite element mixed fracture mode under anti- plane loading of a crack, Theo. App. Fract. Mech., 62 (2012) 26–33. [34] Kotousov, A., Effect of plate thickness on stress state at sharp notches and the strength paradox of thick plates, Int. J. Solids Struct., 47 (2010) 1916–1923. [35] ANSYS Manual Release 8.0.