Fem analysis oF pressure vessel with an investigation oF crack growth on cylindrical surFace analiza mes zbiornika ciśnieniowego z badaniem rozwoju pęknięć na powłoce walcowej

During the exploitation, pressure vessels are subjected to different load types (static, dynamic, thermal, etc.), and their failures usually occur in areas of geometrical discontinuity. The most common geometrical discontinuities are nozzles positioned on the cylindrical shell of the pressure vessel. The previous researches in the field [6, 15, 16] focused on this problem: stress and strain fields on pressure vessels with one nozzle were analysed, using analytical and numerical calculations, as well as experimental methods. Nozzles are usually positioned on cylindrical shells at angle bigger or less than 90o relative to longitudinal axis of the vessel [14, 19]. Each nozzle has its constructive characteristics, as well as specific effect on cylindrical shell strength. The aim of the research was to analyse the influence of the dimensions and positions of the two nozzles on the strain distribution on the vessel cylindrical shell subjected to internal pressure, and then to use obtained distributions for initial crack position prediction and consequent damage analysis. In general, damage to pressure vessels can be thought of as occurring in two stages: crack initiation and crack propagation. However, in conducting damage analyses on pressure vessel, it is often conservative to ignore the crack initiation process, assume that the component already contains a pre-existing crack or defect, and use the analysis technique to estimate the in-service extension of the defect [3, 12]. Hypothesizing the existence of an initial crack or defect seems especially appropriate for the case of large welded structures since the welds may contain defects (lack of fusion, voids, inclusions, etc.). Martina BAlAc Aleksandar GrBovic Aleksandar Petrovic vladimir PoPovic


Introduction
During the exploitation, pressure vessels are subjected to different load types (static, dynamic, thermal, etc.), and their failures usually occur in areas of geometrical discontinuity.The most common geometrical discontinuities are nozzles positioned on the cylindrical shell of the pressure vessel.The previous researches in the field [6,15,16] focused on this problem: stress and strain fields on pressure vessels with one nozzle were analysed, using analytical and numerical calculations, as well as experimental methods.Nozzles are usually positioned on cylindrical shells at angle bigger or less than 90º relative to longitudinal axis of the vessel [14,19].Each nozzle has its constructive characteristics, as well as specific effect on cylindrical shell strength.
The aim of the research was to analyse the influence of the dimensions and positions of the two nozzles on the strain distribution on the vessel cylindrical shell subjected to internal pressure, and then to use obtained distributions for initial crack position prediction and consequent damage analysis.In general, damage to pressure vessels can be thought of as occurring in two stages: crack initiation and crack propagation.However, in conducting damage analyses on pressure vessel, it is often conservative to ignore the crack initiation process, assume that the component already contains a pre-existing crack or defect, and use the analysis technique to estimate the in-service extension of the defect [3,12].
Hypothesizing the existence of an initial crack or defect seems especially appropriate for the case of large welded structures since the welds may contain defects (lack of fusion, voids, inclusions, etc.).

Fem analysis oF pressure vessel with an investigation oF crack growth on cylindrical surFace analiza mes zbiornika ciśnieniowego z badaniem rozwoju pęknięć na powłoce walcowej
To ensure reliability of pressure vessels during service it is necessary to (1) know properties of materials used in their design and ( 2) evaluate vessels' behaviour under different working conditions with satisfying accuracy.Due to various technical and/ or technological requirements, nozzles are usually welded on vessel's shell producing geometrical discontinuities that reduce the safety factor.To evaluate their influence, vessels with two different nozzles were experimentally studied and critical areas for crack initiation have been identified by 3D Digital Image Correlation (DIC) method.After that, the numerical analysis of equivalent 3D finite element model was performed and obtained results were compared with experimental values.In the most critical area, next to the one of the nozzles, crack was initiated and then growth of the damage was simulated using extended finite element method (XFEM).In this paper evaluation of stress intensity factors (SIFs)

sciENcE aNd tEchNology
Hence, characterization of the crack extension property of pressure vessel in terms of fracture mechanics parameters can be quite useful.
In recent years, for evaluation of fracture mechanics parameters like stress intensity factors (SIFs) finite element method (FEM) has been used.Many papers deal with the problems related to crack initiation and damage growth on pressure vessels.But, in the most of these investigations simple geometry was used in FEM calculations; researchers usually carry out simulations on 2D models of specimens [17,8] and very rare on 3D models of real structure [13].Moreover, even on 3D models of pressure vessels crack growth in plane was simulated, and only one crack direction was considered.In practice, cracks can grow in different directions (depending on the loads and constraints) and very often they change plane of propagation.This is hard to numerically simulate and therefore good evaluation of residual life of damaged structure is not easy to achieve.The best estimates are made when SIFs are known for real geometry; SIFs values obtained in FEM simulations on 2D models of specimens must be modified and adjusted in order to be used for residual life evaluation of the real structure.
In this research, numerical method known as extended finite element method (XFEM) was used in order to predict the in-service propagation of crack on the cylindrical shell of the pressure vessel.XFEM has caught a lot of attention since its inception in 1999.It is an alternative method to finite element method (FEM), which allows for the introduction of some knowledge (called enrichment) of the solution into the approximation space, using so-called the partition of unity method (PUM).Discontinuities may be incorporated into the approximation of the unknown field such that the faces and edges of the mesh do not need to match the discontinuity geometry.This way, the method permits the crack propagation without the need to remesh the domain between each step of the simulation.
The cracks are represented with the help of two signed distance functions that are discretized on the same mesh as the displacement field with first-order shape functions.After each step of the propagation simulation, the SIFs are computed from the numerical solution at several points along the crack fronts.Interaction integrals are used to extract the mixed-mode SIFs with the help of auxiliary fields.After that Paris-Erdogan crack growth model, for example, can be used for evaluation of the number of cycles that will grow crack to the critical length in the case of dynamical load.
This paper presents the application of the XFEM to the simulation of the crack propagation in pressure vessel with two nozzles under inservice loading.Firstly, 3D Digital Image Correlation Method (DIC) was used for determination of maximal strain values on the cylindrical shell of pressure vessel under different working conditions.Secondly, numerical model of the vessel was developed and values of stress and strain obtained by FEM were compared to experimental values to check the reliability of numerical model.Finally, based on the obtained results, i.e. identified areas of the greatest strain values between welded nozzles, crack was initiated in one of them and propagated using XFEM.Then, for each propagation step stress intensity factors (SIFs) have been calculated.

Experimental method
For obtaining reliable picture of pressure vessel stress-strain state, several pressure vessels with two adjacent nozzles were fabricated.Experiments were performed on horizontal pressure vessels each with the following dimensions: D = 378.4mm,1.5 a e = , L= 770 mm.Pressure vessels were made of X5CrNi1810 material (EN 10088).Two nozzles, nozzle 1 DN50 1 60.3  d = (m 1 2.9 e = ) and nozzle 2 DN32 ), were placed in the middle of cylindrical shell, so that the longitudinal axis of a nozzle was perpendicular to vessel's longitudinal axis (Figure 1).Minimum distance between two adja-cent nozzle centres was calculated according to standard EN 13445-3.Vessels were tested on the pressure testing installation, subjected to internal pressure of water at the temperature of 20°C (Figure 2).Experimental investigation of strain distribution between two nozzles was carried out using 3D optical method [11,2] considering that this method gives very reliable results.3D system Aramis uses two digital cameras for full field strain measurement and enables precise determination of critical area [21, 9], i.e. the highest strain values on cylindrical shell.Parameters for basic Aramis system setup were: measuring volume 100 x 75, measuring distance 800 mm, camera angle 26°, calibration object CP20 90 x 72.Before starting the measurement, specimen surface had to be prepared and free of grease and oil.At the clean measuring surface, white paint was applied by spraying as the base colour.After the base-colour dried, stochastic pattern of black dots was sprayed.Before using the system, it was necessary to adjust the sensor unit, adjust the angle between the lens focus and aperture.To ensure dimensional consistency of the measuring system, calibration was performed with the help of calibration panel, and the entire system was calibrated.Inner pressure was applied gradually and imaging by means of 3D system was performed manually.Pres-sure increase was in the accordance with procedure defined in standard EN 13445-3.
Initial calculations had shown that the vessel's material would behave within its elasticity limits if applied pressure is less than 1MPa.However, during the experiments the following pressures had been used -0.5MPa, 1MPa, and 1.5MPa in order to investigate full field strain distributions between two closely welded nozzles (elastic plus plastic strain).In this paper results for 1.5 MPa internal pressure are presented and Figure 3 shows strain field obtained in this case (maximum value 0.195%).

Finite element analysis of pressure vessel (numerical model)
Since the experiments with pressure vessels showed what was expected -maximum strain value was always near the welded nozzle DN50 -a numerical model was created with the aim of: a) checking obtained strain values and b) determining the direction of the crack growth after initiation in the area with the maximum strain.Pressure vessel model was designed in CATIA v5 and subsequently exported to Abaqus.Finite element (FE) model of the one half of the vessel was created in Abaqus (Figure 4), and initial numerical simulations (static nonlinear finite element analysis) were performed in order to verify the model fineness [17].Besides, Abaqus included definition of material properties (steel, Young's modulus of elasticity of 210000 MPa, Poisson's ratio of 0.3, UTS = 540 MPa, Yield strength = 230 MPa), loads (uniform pressure of 1.5 MPa) and boundary conditions (type ZSYMM, with restricted one translation and two rotations, Figure 5).The FE model used in simulation had 283605 nodes and 243487 linear hexahedral elements of type C3D8R, and was obtained by means of iteration process which consisted of comparison between the nu-merical results obtained for current, denser mesh and values obtained in calculations with coarser mesh.In each iteration process, mesh step was being refined until difference in stress and strain values in two consecutive steps was less than 5 %.When results obtained by finite element analysis were closed enough to the experimental values, numerical model was accepted as a satisfactory.Figure 6 shows the maximum strain value (0.209%) obtained in FE simulation.Comparing the results obtained by experimental method and numerical simulation, difference of about 7% was calculated confirming that values obtained by FE model can be considered relevant.At the same time, Figure 6 shows that area in which the crack is most likely to occur is around nozzle DN50.This was also confirmed by other authors [20] who showed that cracks appeared mostly in the areas with greatest stress/strain concentration, as can be seen in Figure 7.All of this justified the efforts aimed at defining sufficiently good FE model for a reliable simulation of a crack propagation in pressure vessel with two nozzles.

Stress intensity factors evaluation
To predict life of the component or assembly subjected to time dependent crack growth mechanisms, postulates of fracture mechanics (FM) must be used.The cracking rate can be described using FM parameters such as the stress intensity factors (SIFs).When these factors are known the critical crack size for failure can be computed for given fracture toughness.The fatigue crack growth rate in metals can usually be described by the empirical Paris-Erdogan relationship   size by a safety factor.The critical crack size is computed from the applied stress and fracture toughness and evaluation of service life of the structure can then be obtained by calculating the time (number of cycles) required for crack to grow from its initial size to allowable size.This is all impossible without evaluating SIFs first.There are three types of loading that a crack can experience: Mode I loading, where the principal load is applied normal to the crack plane, Mode II loading corresponds to in-plane shear loading, and Mode III refers to out-of-plane shear.A cracked body can be loaded in any one of these modes, or a combination of two or three modes.Each mode can be described by corresponding SIF.The stress intensity factor is usually given a subscript to denote the mode of loading, i.e., K I , K II , or K III [1].
Stress intensity solutions for many configurations have been published and most of them were obtained from numerical models [5,4,18,10].A variety of numerical techniques have been applied to problems in solid mechanics, including finite difference method, boundary integral equation methods and finite element method.In recent years, the latter two have been applied exclusively, but now it seems that XFEM offers better and easier approach.XFEM is still not fully recognized and needs to prove its practical value to be generally acknowledged.SIFs obtained by using XFEM for a complex 3D geometry are still not regarded as reliable without experimental verification.Here XFEM is applied to study of the crack propagation in cylindrical shell of pressure vessel with nozzles (like that shown in Figure 7) with the purpose of demonstrating its power and contributing to more objective judgment about method usefulness.
Similarly to strain verification method explained in section 3 of this paper, results obtained by XFEM and presented here were verified using experimental data, as described in [23].In brief, The XFEM was first used to the 3-point bending specimen (Figure 8a) to compare numerical values with the experimental results.After successful verification, the XFEM was used to simulate crack growth in casing pipe (Figure 8b), made of API J55 steel by high-frequency welding.
Residual life obtained in simulation was close to that observed in experiment with pipe.
In pressure vessel study, initial crack -defined as a semi-circle surface of radius 2 mm (which is crack size visible by eye) -emanates from the strain concentration location at the fillet between bigger nozzle and the body of the vessel.Abaqus includes defined hexahedron finite element mesh and, as Figure 9 shows, denser mesh was generated in the areas in which crack is expected to propagate (one half of bigger nozzle and along the wall of the vessel); the idea was to increase the accuracy of calculated values of SIFs along the crack front.Abaqus defines initial crack as a separate entity with no element mesh and the first step in 3D analysis of crack propagation is crack "opening" (Figure 10) followed by calculation of pressure vessel stresses which are then used for determination of SIFs in the nodes of the crack front.
It is important to emphasize that there are considerable differences between 2D simulation of crack propagations still dominant in papers [17,8,22] and 3D XFEM simulation shown here; the most significant difference is this: in 2D simulations, values for stress intensity factors are calculated in one point only -at the tip of the crack propagating in plane, whereas, in 3D simulations, the values are calculated in several points/nodes along the crack front that propagates in space.This way, it is possible to determine the stress intensity factors for all three modes, while 2D analysis determine K I and K II only.
Stress intensity factors Modes I, II and III were calculated by using Morfeo/Crack add-in for Abaqus [7].This add-in uses Abaqus solutions to calculate stress intensity factors in nodes of the crack front and generates a file with results.Then, the equivalent stress intensity factor eq K which combines all three SIF modes, was calculated as well as the kink angle (crack propagation angle) which defines the direction in which crack will be propagated in a next step.
After crack opening, Morfeo/Crack for Abaqus offers two choices: forced crack propagation in a plane and free crack propagation.Both options were used and since difference in cracks' paths was neg- ligible results with forced crack propagation in a plane (maximum displacement 0.2 mm per step) will be presented here.Simulation took about 72 hours on Intel CORE i7 CPU, 32GB RAM computer and was stopped after 200 propagation steps when crack reached critical size.

Results and discussion
As Figures 11, 12 and 13 show, the crack propagated in almost vertical plane all the time, with two propagation fronts: one along the wall of the cylindrical part of the vessel and the other along bigger nozzle.It can be clearly seen that crack "separates" finite elements, which is one of the most comprehensive XFEM features.
Table 1 shows values calculated by Morfeo/Crack for Abaqus after each step of crack propagation: curvilinear coordinate of each point along the crack front, coordinates of the crack front points in global xyz system, stress intensity factors for Modes I, II and III, as well as the values of eq K .The number of output values for each propagation step might be a large and depends on the number of points on the crack front, which, again, results from the density of the finite element mesh in propagation areas; this is why values obtained during simulation had to be processed and shown afterward (Tables 2 and 3).
Based on the selected results shown in Table 2 (crack front on the cylindrical part of the vessel) and Table 3 (crack front on the nozzle),

Table 2. Processed SIF values for crack front 1 on the cylindrical part of the vessel
Value of equivalent SIF K eq (MPa mm 0,5 ) Value of SIF Mode I K I (MPa mm 0,5 ) Step number

Min value along the front
Mean value along the front Step number it can be concluded that, during 200 propagation steps, number of points at which SIFs were calculated on both fronts varied a littlefrom 12 to 16 (front 1) and from 22 to 30 (front 2) -which implies that the fronts were formed in the areas with approximately constant mesh density.The number of nodes is larger on front 2 because the wall of the nozzle is thicker than the wall of the vessel body.Tables 2 and 3 also show minimum, maximum and mean values of eq K along the fronts in selected steps, as well as minimum, maximum and mean values of I K that show that mode I loading was the most dominant.It should be noted that eq K values noticeably varied along the crack front during many propagation steps which was expectable considering complex geometry of pressure vessel and 3D crack.For example, in step 30 which is chosen as the most illustrative, maximum eq K value for front 1 was 2804.57MPamm 0,5 whereas minimum value was 2661.2MPamm 0,5 (Figure 14), while for front 2 eq max K was 2032.42MPamm 0,5 and eq min K was 1882.05MPamm 0,5 (Figure 15), showing the differences of more than 8% along the front.Therefore, it was decided that for each step mean eq K would be used as a representative value.Mean I K values were then calculated in an analogous way.

Conclusions
At present, the safety performance of pressure vessels is arousing increasing attention.Many researches have been focusing on the plastic deformation of pressure vessels, but the macroscopic damage processes -including crack initiation and propagation -on these engineering structures haven't been studied yet.This paper brings one of the first attempts of real structure SIFs values estimation in the case when crack simultaneously grows in two perpendicular directions along pressure vessel cylindrical shell and nozzle.Loads measured in experiment and equivalent boundary conditions were used and crack growth obtained in simulation was analogous to cracks' expansions previously observed in vessels' exploitation.This led to conclusion that adequate SIFs values were obtained in calculations since crack growth directions and crack growth rate strictly depend on them.Based on calculated SIFs values, good predictions of residual life of damaged pressure vessel can be obtained.
The behaviour of cracked pressure vessel is explored using the XFEM, and based on the results of simulations and experiments a few conclusions can be drawn:-The initial crack was generated in the area with greatest strain concentration that was obtained during the experiment and subsequently confirmed by numerical analysis; since, in practice, damages occur in the same area (Figure 7), it can be claimed with a great deal of certainty that this is a critical area for a crack to appear on pressure vessels during their usage.
-In the beginning, the crack was positioned in such a way that it reached both cylindrical shell of the pressure vessel and bigger nozzle; therefore, during propagation, it formed two fronts propagating simultaneously with mutual angle of 90 0 (which again corresponds to practical situations).Because of different wall thicknesses of the cylindrical shell and bigger nozzle, these two fronts cannot extend uniformly, which was also shown by 3D simulation performed with the help of XFEM: Figure 13 clearly shows that the crack part on bigger nozzle is considerably shorter than crack part on the cylindrical shell of the pressure vessel.SIF values given in Tables 2 and 3 show that SIFs are considerably bigger for crack front 1 than for crack front 2, which is main reason why the crack propagated more on the cylindrical shell.
-Figure 16 shows an interesting observation: mean values of eq K as almost linear increase during propagation of both fronts (coefficients of determination are high: 0.9959 and 0.995 respectively).This means it is possible to establish correlation between crack length, vessel wall thickness and stress intensity factor values and then evaluate crack propagation speed in a certain vessel area.An effort was made in that direction and result is shown in Figure 17 cycles crack length is about 20 mm and after next 250 cycles its length is doubled meaning that crack reaches critical size after which complete failure occurs soon.Obtained small number of cycles indicates that vessel of this wall thickness under given pressure must be crack free because damage on the cylindrical shell grows extremely fast.
However, next investigations of pressure vessel integrity must consider the influence of the type of material and wall thickness on the rate of crack growth.
-In the end, it must be kept in mind that the main advantage of XFEM lies in possibility of SIFs values evaluation on complex cracked geometry (like pressure vessel with nozzles) which then can be used for predictions of crack propagation paths and evaluations of residual life, but -at the same time -XFEM results are mesh sensitive and depend on the mesh density in the fracture process region.Mesh size must be determined carefully to ensure the computational efficiency and accuracy; therefore, experimental verification of FE model is necessary, at least in the phase when the object is still undamaged.

Fig. 5 .
Fig. 5. Boundary condition of the type Symmetry (U z = ROT x = ROT y = 0) was applied in selected (red) area

n
da dN C K = here / da dN s the crack growth per cycle, K s the stress intensity range during the fatigue cycle max min K K K = − ()and C and n are material constants.Damage tolerance allows subcritical cracks to remain in a component, but when they grow an allowable flaw size must be defined, usually by dividing the critical

Table 3 .
Processed SIF values for crack front 2 on the nozzle