Analysis of surface crack in cylinder by finite element alternating method

Masayuki Kamaya1 Institute of Nuclear Safety System,Inc.(INSS), 64Sata,Mihama-cho,Fukui919-1205,Japane-mail:kamaya@inss.co.jpToshihisa NishiokaKobe University, 5-1-1Fukaeminamimachi,Higashinada-ku,Kobe685-0022,Japan e-mail:nishioka@maritime.kobe-u.ac.jp Analysis of Surface Crack in Cylinder by Finite Element Alternating MethodThefinite element alternating method(FEAM),in conjunction with thefinite element analysis(FEA)and the analytical solution for an elliptical crack in an infinite solid subject to arbitrary crack-face traction,is used for evaluating the stress intensity factor (SIF)of surface cracks in a cylinder.The major advantage of this method is that the SIF can be calculated by using the FEA results for an uncracked body.A newly developed system allows the FEAM to be performed by a simple method,which consists of the conventional FEA for an uncracked body and a subroutine for the FEAM alternating procedure.It is shown that the system can derive the precise SIF of circumferential, longitudinal,and inclined surface cracks in a cylinder.The crack growth predictions are performed for an inclined crack and projected longitudinal and circumferential crack in a cylinder.The results suggests that the crack characterizing procedure prescribed in Sec. XI may cause an unconservative evaluation in the crack growth prediction,and that the FEAM is valid for complex problems,to which the SIF evaluation by the FEA cannot be adopted easily.͓DOI:10.1115/1.1903003͔1IntroductionThefinite element alternating method͑FEAM͒is an alternating technique in conjunction with thefinite element analysis͑FEA͒and the analytical solution,which was derived by Nishioka et al.͓1,2͔and cited as the VNA solution,for an elliptical crack in aninfinite solid subject to arbitrary crack-face traction.The FEAM can be used for evaluating the stress intensity factors͑SIFs͒of elliptical embedded and semielliptical surface cracks in afinite body,and has been utilized for various problems,such as surface cracks in pressure vessels͓3͔,multiple interacting crack problems ͓4–7͔,etc.͓8,9͔.The major advantage of this method is that the SIF can be calculated by using the FEA results for an uncracked body.Therefore,the SIF of cracks in complex structural compo-nents or complex stress state can be evaluated relatively easily as far as the conventional FEA result of the uncracked body is avail-able.However,much effort is required in order to perform this method due to the complexity of the VNA solution and the alter-nating procedure.In the present study,an analysis system using the FEAM for evaluating the SIFs of surface cracks is developed.In order to make the most of the advantage of the FEAM,the conventional FEA for the uncracked body with initial loading condition is sepa-rated from the FEAM alternating routine.This separation hides the complex alternating procedure in the subroutine which has a simple interface.The subroutine requires the geometric data of cracks,such as size,location,and inclination,in addition to the FEA results of the original problem,in which the crack is not included.The commercial code,ABAQUS version6.3͓10͔,is used for FEA.The analysis system is applied to circumferential,longi-tudinal,and inclined surface cracks in a cylinder.The crack growth analyses are also performed in order to investigate the crack characterizing procedure for an inclined crack.2Finite Element Alternating Method2.1Basic Procedure.The FEAM is based on the principle of superposition,the basic procedure of which is shown in Fig.1.At first,the FEA is solved for the uncracked body under the initial loading,and the stress around the crack is calculated by the VNA solutions based on the crack-face traction obtained by the FEA.In order to satisfy the free surface conditions,the second FEA is performed under the surface load that is obtained by reversing the stresses at the surface derived by the VNA solutions.The FEA result is quoted as the traction for the VNA solution again.The alternation between the FEA and stress calculation by the VNA solutions are repeated until the traction becomes negligible com-pared to the initial load.2.2Solution for an Elliptical Crack in an Infinite Solid Subject to Arbitrary Crack-Face Traction(the VNA Solution). In this section,only the Mode I problem is considered.The solu-tion for Modes II and III is given in Ref.͓1͔.Let x3be the coor-dinate normal to the elliptical crack face,and let the tractions along the crack surface be expressed in the form␴33͑0͒=͚i=01͚j=01͚m=0M͚n=0m A3,m−n,n͑i,j͒x12m−2n+i x22n+j,͑1͒where A’s are undetermined coefficients and the parameters i and j specify the symmetries of the load with respect to the axes of the ellipse,x1and x2,wherein the commas are introduced for conve-nience only.M is an arbitrary integer which is related to the order of the polynomial.The solution corresponding to the arbitrarily distributed crack-face tractions expressed by Eq.͑1͒can be assumed in terms of the harmonic potential functions:f3=͚i=01͚j=01͚k=0M͚l=0kC3,k−l,l͑i,j͒F2k−2l+i,2l+j,͑2͒where C’s are undetermined coefficients and F’s are basic poten-tials wherein again the commas are used for convenience only. The components of stress␴ij as well as displacement can be ex-pressed by the potential function f3͓1͔.Thus,stress is given in a matrix form͕␴͖=͓P͔͕C͖,6ϫ16ϫN Nϫ1͑3͒where͓P͔is a function of the coordinates͑x1,x2,x3͒and N is the1Corresponding author.Contributed by the Pressure Vessels and Piping Division for publication in the J OURNAL OF P RESSURE V ESSEL T ECHNOLOGY.Manuscript received October28,2003;final manuscript received January31,2005.Review conducted by:Sharif Rahman.total number of coefficients A and C .Satisfying the boundary condition on the crack surface,the re-lation between the parameters A and C in Eqs.͑1͒and ͑2͒can be summarized in a matrix form͕A ͖=͓B ͔͕C ͖.N ϫN ϫN N ϫ1͑4͒The detailed expression of the components of ͓B ͔is given in Ref.͓1͔.For a complete polynominal loading expressed by Eq.͑1͒,the maximum degree of polynomial M c and the number of coeffi-cients N can be expressed,respectively,by M c =2M +1and N =͑M +1͒͑2M +3͒.For any load given in Eq.͑1͒,the coefficients C can be deter-mined by solving Eq.͑4͒.Once the coefficients C are determined,the Mode I stress intensity factors corresponding to the given load are evaluated from the following equation:K I =8␮ͩ␲a 1a 2ͪ1/2Q 1/4a 1a 2H ,͑5a ͒whereQ =a 12sin 2␪+a 22cos 2␪,͑5b ͒H =͚i =01͚j =01͚k =0M͚l =0k͑−2͒2k +i +j͑2k +i +j+1͒!ͩcos ␪a 1ͪ2k −2l +iͩsin ␪a 2ͪ2l +jC 3,k −l ,l ͑i ,j ͒.͑5c ͒␪denotes the elliptical angle as shown in Fig.2and ␮is the shearmodulus.The elliptic angle concerning point P is measured from the x 1axis.The elliptical crack size is defined by major semiaxis a 1and minor semiaxis a 2.2.3Procedure of Analysis.The analysis system is intended to simplify performing the FEAM.For this purpose,the complex alternating routine is placed in a subroutine with a simple inter-face.The subroutine is carried out using the data about crack geometry and the FEA result of the original problem.The flow chart of the analysis procedure is shown in Fig.3.Step ͑0͒Perform the FEA for the uncracked body under the given external loads.The information about the model,such as model geometry,location of nodes,and integration points as well as the stresses is utilized in subsequent FEAM.In ABAQUS ,these data are stored in the “result file”and “restart file.”Step ͑1͒The geometric data of the crack,such as size,position and inclination,is stored in this step.Step ͑2͒In order to evaluate the crack-face traction expressed by Eq.͑1͒,the interpolation points ͑IPs ͒,at which the stress is calculated,are determined.In the present study,the IPs are se-lected as shown in Fig.4.Step ͑3͒For efficient calculation,some basic parameters are calculated before the iteration process.In order to obtain the stress at the IPs,the interpolations are carried out using the stress at adjacent integration points of the finite element model.In this step,the related integration points used for the interpolations are determined for each IP.The surface information,such as element number at surfaces and distance from crack center,which is re-lated to the matrix ͓P ͔of Eq.͑3͒,is also calculated in this step.Step ͑4͒The initial stress at crack-face is evaluated at the IPs by interpolation using the volume coordinate.Step ͑5͒To satisfy the stress boundary condition at the crack-face,reverse the stress evaluated in Step ͑4͒.Then evaluate the coefficients A in Eq.͑1͒for applied stress by using the least squares method.The maximum degree of polynomial M c =5is used in the present study.Step ͑6͒The coefficients C of Eq.͑2͒can be derived by solving Eq.͑4͒with coefficients A evaluated in the previous step.Step ͑7͒Calculate the SIF for the current iteration by substitut-ing coefficients C in Eq.͑5͒.Step ͑8͒Calculate the residual stresses on external surfaces of the body due to the loads applied in Step ͑5͒.To satisfy the stress boundary condition at the free surface,reverse the residual stresses and calculated equivalent nodal forces.The stresses at the external surfaces can be calculated using Eq.͑3͒,in which the matrix ͓P ͔was determined previously in Step ͑3͒.Although the matrix ͓P ͔has stress singularity at the crack front,the magnitude of the matrix decays rapidly with the distance from the crack front.Thus,at the surfaces that are far from the crack,the residual stresses can be neglected in order to reduce the calculation time.Step ͑9͒Conduct the FEA for the uncracked body under exter-nal loads consisting of the nodal forces calculated in Step ͑8͒.Step ͑10͒Again,calculate the stress at IPs according to the same procedure as Step ͑4͒.Step ͑11͒Repeat all steps in the iteration process,from Step ͑5͒to Step ͑11͒,until the stresses acting on the crack face,which is calculated in Step ͑10͒,become negligible.In the present study,the threshold for the end of the iteration was set to 1%of the mean initial stress obtained in Step ͑4͒.Step ͑12͒To obtain the final SIF,add the SIF calculated inStepFig.1Finite element alternating method for finite cracked body under remoteloadingFig.2Elliptic angle ␪for point P on the edge of an elliptical crack͑7͒of all iteration processes.The VNA solution is valid for entire elliptical crack.Therefore, in order to calculate the SIF of a semielliptical crack,it is neces-sary to assume the stress at the portion existing outside the body in Step͑5͒.In Ref.͓1͔,numerical experimentation was conducted for arriving at an optimum stress distribution on a crack surface extended into thefictitious region outside of the body.From the conclusion of the numerical experimentation,in the present study,thefictitious stress for the region x2Ͻ0is assumed to vary in the x1direction and remains constant in the x2direction as shown in Fig.5.3Analyses of Surface Crack in Cylinder3.1Circumferential and Longitudinal Surface Cracks.The FEAM is applied to circumferential and longitudinalsurface Fig.3Flow chart of thefinite element alternating methodcracks in a finite cylinder,respectively.The cylinder is subjected to an axial stress P z and internal pressure P i .The model geometry and crack shape are shown in Fig.6.The thickness,t ,and half-length of the cylinder,L ,are t /R i =0.1and L /t =15,where R i is inside radius of the cylinder.The finite element mesh for the uncracked cylinder is shown in Fig.7.Due to the symmetries of the problem,only one quarter of the cylinder was modeled by finite elements.Fine elements were adopted around the crack portion for accurate analysis,although the mesh is not included in the crack.The finite element mesh shown in Fig.7consists of 18,0008-noded solid elements.The relationship between the accuracy of the analysis and the density of the mesh is discussed elsewhere ͓11͔.The crack size and loading conditions used for each problem are summarized as follows:Circumferential crack:b /a =0.5,b /t =0.4P z =1.0,P i =0.0Longitudinal crack:b /a =0.4,b /t =0.5P z =0.0,P i =1.0,where Poisson’s ratio of the material is taken to be 0.3and this value is used for all calculations of the present study.According to the procedure shown in Fig.3,the FEA for the uncracked body shown in Fig.7is performed under the initial loading conditions of P z and P i for each crack.Then the FEAM alternating subroutine is carried out.The data required in the sub-routine,such as the crack shape,location and the FEA results,are automatically taken into the subroutine through the data readpro-Fig.6Geometry of crackedcylinderFig.7Finite element mesh for an uncrackedcylinderFig.4Location of interpolation points …total 709points are employed…Fig.5Residual stress distribution over the entire crack surfacegram.The obtained SIF is shown in Fig.8together with ones obtained by conventional FEA given in Refs.͓12,13͔.The SIFs are normalized by the following equation:F T =K IP z ͱb ͑circumferential crack ͒,͑6a ͒F L =K I P i ͱ␲b ͑longitudinal crack ͒,͑6b ͒where K I denotes the SIF of mode I.The present results agree well with those obtained in the references.Based on the reference re-sults,the calculation error of the present study is about 1.2%and 3.7%in the circumferential and longitudinal cracks,respectively,although it is mentioned that there is about a 4%error in the reference data of the longitudinal crack.Figure 9shows the variation of residual stress on the crack surface with each cycle of iteration in the alternating technique ͓see Step ͑10͒in the flow chart ͔.It is seen that the residual stress decreases rapidly and monotonically with the number of itera-tions.After the third iteration,the residual stress becomes 1.8%and 5.7%of the initial stress in the circumferential and longitudi-nal cracks,respectively.3.2Inclined Surface Crack.The SIF of surface cracks,which rotated the x 2-axis by angle ␤as shown in Fig.10,is evaluated by the FEAM.The SIF of arbitrarily inclined surface cracks can be obtained by just changing the input data at Step ͑2͒,while the same FEA result for the uncracked body is used.Calcu-lations have been performed for all the possible combinations of inclined angle of ␤=0°,15°,30°,and 45°and crack aspect ratio b /a =0.25,0.5,and 0.8,where the depth of crack is taken to be b /t =0.4.The cracks exist at the center of the longitudinal direc-tion of the cylinder that is subjected to an axial stress P z .The radius and length of the cylinder are t /R i =0.1and L /t =15,which is the same geometry as the previous model.The finite element mesh is shown in Fig.11.In order to analyze the inclined cracks,the mesh is generated for whole of the model and consists of 58,0008-nodeed solid elements.The crack front is under a mixed mode deformation represented by K I ,K II ,and K III ,the SIFs for three fundamental deformation modes.They are normalized by the following equations,respec-tively,F 1=K IP z ͱb ,F 2=K IIP z ͱb ,F 3=K IIIP z ͱb .͑7͒The obtained values of F 1,F 2,and F 3are shown in Fig.12for ␤=0°,15°,30°,and 45°,respectively.F 1decreases as theincli-Fig.8Normalized stress intensity factor for a surface crack in a cylinder:…a …circumferential crack;and …b …longitudinalcrackFig.9Change in the residual stress on the cracksurfaceFig.10Definition of inclination angle ␤nation angle␤increases,although F2and F3show the oppositetendency.However,F1is shown to dominate F2and F3,especiallyfor small inclination.4Crack Growth Prediction of Inclined Crack in Cyl-inderGrowth prediction for a surface crack inside a cylinder is per-formed by using the FEAM.According to the ASME boiler andpressure vessel code Sec.XI͑Sec.XI͓͒14͔,an inclined crack,which exists in a cylinder and that does not lie in either the lon-gitudinal or circumferential plane,is projected onto these planesas shown in Fig.13,and the growth prediction is performed foreach projected crack separately.In order to examine the validity ofthis rule,the crack growth predictions are performed for the pro-jected longitudinal crack,circumferential crack and original in-clined crack,and the results are compared.The crack exists at the center of the cylinder and the orientationis inclined by␤=45°to the longitudinal direction.The geometriesof the initial crack and cylinder are summarized as follows:a o=5mm,a T=a L=a o cos͑45°͒b o=b T=b L=2.5mmR i=100mm,t=10mmFig.11Finite element mesh for an uncracked cylinderFig.12Normalized stress intensity factor for a inclined surface crack in a cylinder…b/t=0.4,t/Ri=0.1,L/t=15…:…a…inclination angle␤=0°;…b…Inclination angle␤=15°;…c…inclination angle␤=30°;and…d…inclination angle␤=45°.L =150mm.The cylinder is subjected to the axial stress of 200MPa and in-ternal pressure of 20MPa.The internal pressure of 20MPa yields about 210MPa hoop stress at the internal surface in the current problem.The finite element mesh,which is shown in Fig.11,is used for the analysis.The procedure of the crack growth prediction is summarized in Fig.14.The SIF is calculated by the FEAM based on the cracksize at each step.The crack extension length in the depth direction ⌬b and surface direction ⌬a are assumed to be calculated by:⌬a =DK ͑␪=0°͒m⌬t ⌬b =DK ͑␪=90°͒m ⌬t ,͑8͒where D and m are material constants,and D =5.24ϫ10−14and m =3.24,which are values of the stress corrosion cracking are adopted ͓15͔.K ͑␪=90°͒and K ͑␪=0°͒are the SIF of mode I obtained by the FEAM where ␪denotes the crack tip position defined in Fig.2.⌬t is the interval of the step,which is controlled so that ⌬b and ⌬a do not exceed 5%of the crack size.The extension lengths ⌬b and ⌬a are given in m ,the SIF in MPa ͱm and ⌬t in seconds.Then,the FEAM is again performed based on the extended crack size.The cracks are assumed to maintain the semielliptical shape.Figure 15shows the relationship between the time and crack size.The crack growth rate of the inclined crack is faster than those of the projected circumferential and longitudinal cracks.The depth of cracks b T and b L becomes 4.7and 5.6mm at 837h,respectively,while b o becomes 8.0mm.This implies that the pro-jection procedure prescribed in Sec.XI may cause an unconser-vative evaluation in the crack growth prediction.By using the FEAM,we can evaluate the SIF of the inclined crack without difficulty and can reasonably predict the growth.5ConclusionIn the present study,a FEAM analysis system was developed.This system allows the FEAM to be performed by a simple method,which consists of the conventional FEA for theuncrackedFig.14Procedure of the crack growthpredictionFig.13Inclined crack and projected circumferential and axial crack incylinderFig.15Relationship between time and crack size obtained by the crack growth prediction using the FEAM:…a …surface direc-tion;and …b …depth directionbody and the subroutine for the FEAM alternating procedure.The system was successfully used to obtain accurate SIFs of circum-ferential,longitudinal,and inclined surface cracks in a cylinder. Then,the crack growth predictions were performed for an inclined crack and projected longitudinal and circumferential cracks in a cylinder.The results showed that the crack characterizing proce-dure prescribed in Sec.XI may cause an unconservative evalua-tion in the crack growth prediction.References͓1͔Nishioka,T.,and Atluri,S.N.,1983,“An Analytical Solution for Embedded Elliptical Cracks,and Finite Element Alternating Method for Elliptical Surface Cracks,Subjected to Arbitrary Loadings,”Eng.Fract.Mech.,17,pp.247–268.͓2͔Vijayakumar,K.,and Atluri,S.N.,1981,“An Embedded Elliptical Crack,in an Infinite Solid,Subject to Arbitrary Crack-Face Tractions,”ASME J.Appl.Mech.,48,pp.88–96.͓3͔Nishioka,T.,and Atluri,S.N.,1982,“Analysis of Surface Flaw in Pressure Vessels by a New3-Dimensional Alternating Method,”ASME J.Pressure Vessel Technol.,104,pp.299–307.͓4͔O’Donoghue,P.E.,Nishioka,T.,and Atluri,S.N.,1984,“Multiple Surface Cracks in Pressure Vessels,”Eng.Fract.Mech.,20,pp.545–560.͓5͔Nishioka,T.,Tokunaga,T.,and Akashi,T.,1994,“Alternating Method for Interaction Analysis of a Group of Micro-Elliptical Cracks,”J.Soc.Mater.Sci.Jpn.,43,pp.1271–1277.͓6͔Nishioka,T.,Akashi,T.,and Tokunaga,T.,1994,“On the General Solution for Mixed-Mode Elliptical Cracks and Their Applications,”Tran.JSME,60,pp.364–371.͓7͔Nishioka,T.,and Kato,T.,1999,“An Alternating Method Based on the VNA Solution for Analysis of Damaged Solid Containing Arbitrarily Distributed Elliptical Microcracks,”Int.J.Fract.,97,pp.137–170.͓8͔Raju,I.S.,Atluri,S.N.,and Newman,Jr.J.C.,1989,“Stress-Intensity Factors for Small Surface and Corner Cracks in Plates,”ASTM STP1020,Philadel-phia,USA,pp.297–316.͓9͔Krishnamurthy,T.,and Raju,I.S.,1990,“A Finite-Element Alternating Method for Two-Dimensional Mixed-Mode Crack Configurations,”Eng.Fract.Mech.,36,pp.297–311.͓10͔ABAQUS Inc.,2002,“ABAQUS/Standard User’s Manual Ver.6.3,”ABAQUS Inc., USA.͓11͔Kamaya,M.,and Nishioka,T.,2004,“Evaluation of Stress Intensity Factors by Finite Element Alternating Method,”PVP-481,pp.113–120.͓12͔Bergman,M.,1995,“Stress Intensity Factors for Circumferential Surface Cracks in Pipes,”Fatigue Fract.Eng.Mater.Struct.,18,pp.1155–1172.͓13͔Raju,I.S.,and Newman,Jr.,J.C.,1982,“Stress Intensity Factors for Internal and External Surface Cracks in Cylindrical Vessesls,”ASME J.Pressure Ves-sel Technol.,104,pp.293–298.͓14͔ASME,2001,“ASME Boiler and Pressure Vessel Code Section XI,”New York,USA,Fig.FC-2400-3.͓15͔Kamaya,M,and Totsuka,N.,2002,“Influence of Interaction between Multiple Cracks on Stress Corrosion Crack Propagation,”Corros.Sci.,44, pp.2333–2352.。