UDK 621.039.536.2 Original scientific article/Izvirni znanstveni članek ISSN 1580-2949 MTAEC9, 47(6)825(2013) NUMERICAL CALCULATIONS OF STRESS-INTENSITY FACTORS FOR A WWER-1000 REACTOR PRESSURE VESSEL UNDER A PRESSURIZED THERMAL SHOCK NUMERIČNI IZRAČUNI FAKTORJA INTENZITETE NAPETOSTI V REAKTORSKI TLAČNI POSODI WWER-1000 PRI TOPLOTNEM ŠOKU POD TLAKOM Sergey Vladimirovich Kobelsky G. S. Pisarenko Institute for Problems of Strength, National Academy of Science of Ukraine, Kiev, Ukraine ksv@ipp.kiev.ua Prejem rokopisa — received: 2013-07-24; sprejem za objavo - accepted for publication: 2013-10-03 Two methods of determining stress-intensity factors based on the calculations of J- and G-integrals are analyzed for a simulated emergency-cooling regime of a WWER-1000 reactor pressure vessel. A finite-element analysis of the convergence of the results determining the stress-intensity factors using different approaches and different finite-element meshes was done. Keywords: reactor pressure vessel (RPV), finite-element method (FEM), J- and G-integral, crack, stress-intensity factor (SIF) Analizirani sta bili dve metodi za določanje faktorja intenzitete napetosti na osnovi izračuna J- in G-integrala pri simuliranem režimu hitrega nujnega ohlajanja tlačne posode reaktorja WWER-1000. Narejena je bila analiza raztrosa rezultatov s končnimi elementi pri določanju faktorja intenzitete napetosti in različnih mrežah končnih elementov. Ključne besede: tlačna posoda reaktorja (RPV), metoda končnih elementov (FEM), J- in G-integral, razpoka, faktor intenzitete napetosti (SIF) 1 INTRODUCTION The safety control of nuclear power plants and, in particular, of reactor pressure vessels (RPV) as their most critical elements, is one of the important theoretical and practical problems. One of the promising ways of solving this problem is mathematical modeling of the kinetics of the thermomechanical state of a reactor pressure vessel, which requires solutions of non-stationary non-linear problems of thermomechanics. The computational justification for the strength of the RPV is carried out on the basis of an analysis determining fracture-mechanics parameter values, namely, stress-intensity factors (SIF). The complexity of the problem under consideration, a large amount of computations and the necessity of performing multiple-choice calculations make it essentially impossible to justify the convergence of the obtained results when solving practical problems with the finite-element method (FEM).1 In the majority of cases, for the structures of a complex shape, calculations are performed on finite-element meshes with a spacing of 0.5 mm to 1 mm in the vicinity of crack front points. At the same time, in the finite-element method, the convergence of numerical-calculation results is a necessary but not a sufficient condition for their reliability. The degree of reliability of the results can be enhanced by comparing the values of the analyzed parameters obtained with the use of various hypotheses or criteria.2 It should be noted that the conditions for attaining the numerical convergence cannot automatically be extended to another class of problems. Thus, the numerical-convergence conditions for the problem to be solved using the deformation-plasticity theory, cannot be fulfilled for the problem to be solved on the basis of the flow theory. The goal of the present paper is to analyze the convergence of the results determining the SIF calculated using various concepts on different finite-element meshes. 2 PROBLEM STATEMENT For typical emergency-cooling conditions, known as the pressurized thermal shock, two computational models of the WWER-1000 reactor pressure vessel (RPV), with a built-in crack and without it, are analyzed (Figure 1). A circumferential semi-elliptical crack with the depth of a = 15 mm and the ratio a/c = 1/3 of the half-axis3 that is located at the level of welded joint No. 4 is postulated. The procedure of embedding the crack in the finite-element model of a pressure-vessel fragment is used. 3 METHOD OF CALCULATION AND THE SOFTWARE The problem is solved in a three-dimensional formulation using a mixed finite-element method scheme Figure 1: Computational model of the WWER-1000 reactor pressure vessel Slika 1: Model za izračune tlačne posode reaktorja WWER-1000 (MFEM).4 The main advantage of the MFEM, as compared to the classical FEM approach in the form of the displacement method, is the possibility of ensuring the continuity of the approximation not only for displacements but also for stresses and strains as well as the possibility of accurately satisfying the static boundary conditions on the body surface. The mixed scheme is implemented in the SPACE&RELAX software package created at the G. S. Pisarenko Institute for Problems of Strength of the National Academy of Sciences of Ukraine. The procedure of solving the problem includes three stages. At the first stage, the temperature fields are determined at the specified instants of time. At the second one, the fields of displacements, stresses and strains for the reactor-pressure-vessel (RPV) model are calculated taking into account the thermomechanical loading history. At the third stage, a refined calculation of the fields of displacements, stresses and strains is performed for the RPV fragment with a built-in crack, and the values of the stress-intensity factors (SIF) are determined. When solving the problem of the elasto-plastic formulation, the loading process is divided into separate time steps, within each of which the elasto-plastic problem is solved considering the residual plastic strains determined at the previous step. To solve the elasto-plastic problem at the current stage of loading, a two-step iteration process is used. The problem of the plasticity theory is solved with the method of variable parameters of elasticity, whereas at the internal steps pertaining to elasticity, the problems are solved using the conjugate-gradient method with the zeroth-order initial approximation. 4 METHODS OF DETERMINING THE SIF VALUES To determine the SIF values, the /-integral and G-integral concepts are used. In the problems of both the linear and nonlinear theories of elasticity as well as the deformation-plasticity theory, the /-integral value characterizes the energy-release rate for a virtual infinitesimal crack extension.5,6 The applicability of the /-integral concept is limited by the following requirements: the material is nonlinearly elastic and homogeneous; under elasto-plastic deformation, the load increases in proportion to one parameter, that is, the deformation-plasticity theory is applicable. The use of the /-integral in the flow theory has no theoretical justification. In this paper, the method of equivalent volume integration (EVI) proposed by G. P. Nikishkov is used for the calculation of the /-integral value.7 The G-integral of the crack closure is defined as the work of adhesive forces acting in a body ahead of the crack front and preventing a separation of its surfaces when the crack propagates negligibly. For the instant of the crack growth onset when the value of the crack propagation can be considered as infinitesimal, the G-integral is calculated with formula:8 G =1 a ■ A (1) Here, a is the stress calculated normally to the crack plane at the finite-element node located ahead of its front; A is the opening displacement calculated as a difference in the displacements normal to the crack plane for two mesh nodes adjacent to the crack front and lying on its different surfaces. In contrast to the /-integral, the G-integral value remains meaningful when used in the flow theory for an arbitrary loading history, because it defines the specific work necessary for the crack growth onset in an elastic-plastic body. The G-integral value is likely to be meaningful only on condition that tensile (positive) stresses are acting at the crack tip. Thus, it can be noted that when solving elasto-plastic problems with the flow theory for a complex law of loading, the use of the G-integral concept is more justified than the /-integral concept. 5 FINITE-ELEMENT MESH Figure 2 illustrates a fragment of the finite-element mesh being constructed in the vicinity of the crack front. The mesh extending along the front consists of three parts - the core, the transition zone and the coarsening zone. The core size is selected to be small in comparison with the size of the transition zone (5-10 %). A mesh that is regular in the plane perpendicular to the crack front is constructed in the core. The size of the mesh spacing along the front is adjusted by assigning the Figure 2: Finite-element mesh constructed in the vicinity of the crack front Slika 2: Mreža končnih elementov, konstruirana v bližini cela razpoke refinement factor whose value is selected so that, in the vicinity of the points of interest (e.g., the deepest point of the front), the ratio of the mesh spacing values along the front and in the plane directed to it does not exceed 10 x 1 x 1. The fragment of the finite-element mesh constructed in the crack plane is shown in Figure 3. The mesh shown in Figure 3a allows obtaining reliable data on the fracture parameters for the points located along the crack front (point 1 is the deepest point of the front) but does not allow evaluating their values at point 2 on the boundary between the base metal and the cladding. The mesh shown in Figure 3b makes it possible to obtain reliable data on the values of the stress-in- Figure 3: Finite-element mesh constructed in the crack plane: a) with an unclosed contour, b) with a closed contour Slika 3: Mreža končnih elementov, postavljena v ravnini razpoke: a) z nezaprtim obrisom, b) z zaprtim obrisom tensity factors and the corresponding integrals at points 1 and 2. Unless otherwise specified, the results below were obtained using the unclosed crack model. 6 SOLUTION OF THE PROBLEM IN THE ELASTO-PLASTIC FORMULATION TAKING INTO ACCOUNT THE LOADING HISTORY When solving the problem of the elasto-plastic formulation, taking into account the loading history, we analyzed the sensitivity of the SIF values determined with the following factors: • the use of the initial RPV model with a built-in crack and without it; • the use of different criteria for stopping the iterative process; • the use of the finite-element meshes of different densities. The SIF values presented in all the plots were normalized with the SIF value determined using the G-inte-gral for the crack with the mesh parameters of 21 pm x 5.5 pm. Relative temperature T - Tw is set along the horizontal axe, where Tw is the critical-brittleness temperature of the metal under the initial conditions. 6.1 Analysis of the influence of the crack presence in the initial RPV model and the used criteria for stopping the iteration process In the method of variable elasticity parameters, the use of two criteria for stopping the iteration process is analyzed - based on the discrepancy (2) and the G-inte-gral value (3): (r, r) < £, X (rj, rt) j=1 < (Gi -1, G) * 5*£ g (2) (3) where r is the discrepancy vector, G-k and Gi are the G-integral values at the i-k-th and i-th step iterations with respect to the non-linearity, k > 1 is the parameter, £r = 1 x 10-16, £o = 1 x 10-3. Figure 4 shows the plots of the results for fracture parameters - CTOD and relative SIF values calculated using the G-integral for the problems that were solved on three meshes (1000 pm x 600 pm, 400 pm x 275 pm, 47.5 pm x 27.5 pm) under the following conditions: • the initial RPV model with a built-in crack, the stopping criterion is (2), • the initial RPV model without a built-in crack, the stopping criterion is (3). The analysis of the results allows one to conclude that it is unreasonable to build a crack into the initial full RPV model. The use of criterion (3) instead of criterion 80 J 60 E d) ¡40 0J Cl O o 20 1,2 s. - 3,4 ■ ✓ 5,6 / =5= (a) 50 100 150 200 Relative temperature, °C 250 300 1 0,8 I 0,6 LL ^0,4 LL CO 0,2 0 1,2 s7 /7 /3,4 / 5,6 (b) 50 100 150 200 Relative temperature, °C 250 300 Figure 4: Plots showing the results of calculating the fracture parameters: a) is for the crack-tip opening displacement, b) is for the relative SIF values: 1 - 1000 ¡m x 600 ¡m (1), 2 - 1000 ¡m x 600 ¡m (2), 3 -400 ¡m x 275 ¡m (1), 4 - 400 ¡m x 275 ¡m (2), 5 - 47.5 ¡im x 27.5 ¡m (1), 6 - 47.5 ¡m x 27.5 ¡m (2) Slika 4: Prikaz rezultatov izračunov parametrov preloma: a) je za premik konice razpoke, b) je za relativne vrednosti SIF: 1 - 1000 ¡m x 600 ¡m (1), 2 - 1000 ¡m x 600 ¡m (2), 3 - 400 ¡m x 275 ¡m (1), 4 -400 ¡m x 275 ¡m (2), 5 - 47,5 ¡m x 27,5 ¡m (1), 6 - 47,5 ¡m x 27,5 ¡m (2) (2) did not lead to a noticeable difference in the results of the problem solution (the maximum difference in the SIF values was found to be less than 0.1 %) but made it possible to reduce considerably the total time of the solution (if the time of the problem solution with the use of criterion (2) can be taken to be a unity, then, with the use of criterion (3), the problem time was 0.55 for the mesh of 1000 pm x 600 pm, 0.4 for the mesh of 400 pm x 275 pm, and 0.7 for the mesh of 47.5 Mm x 27.5 pm). The character of the SIF variation curves (Figure 4b) shows that the effect of the mesh-spacing value is minimum at the stages of active loading. However, with an increase in the spacing, the descending branch of the curve departs increasingly more to the right, showing a tendency to be vertical in the limit. Particular attention should be given to the stress behavior during the mesh refinement - under unloading, a compressive stress zone occurs in the vicinity of the crack tip. The instant of time, at which this zone begins to be fixed, and its size are essentially dependent on the spacing value in the vicinity of the crack front. Thus, for the mesh of 1000 pm x600 pm, the compressive stresses did not occur at all, whereas for the mesh of 400 pm x 275 pm they appeared at the crack tip only, at the time instant of 6400 s. Figure 5: Plots of the time dependence of the compressive-stress-zone size in the vicinity of the crack tip: □ 21 ¡m x 5.5 ¡m, ♦ 47.5 ¡m x 27.5 ¡m Slika 5: Prikaz odvisnosti velikosti tlačne cone na čelu konice razpoke: □ 21 ¡m x 5,5 ¡m, ♦ 47,5 ¡m x 27,5 ¡m For the mesh of 21 pm x 5.5 pm, the last SIF value was obtained for the time instant of 2000 s, whereas for the mesh of 1000 ¡m x 600 ¡m, the SIF values were obtained up to the time instant of 6 400 s. The plots of the time dependence of the compressive-stress-zone size in the vicinity of the crack tip for the meshes of 47.5 ¡m x 27.5 ¡m and 21 ¡m x 5.5 ¡m are given in Figure 5. 6.2 Use of meshes oof different densities The solution of the problem with the finite-element method can be considered to be reliable only upon reaching the convergence. It is common practice to perform an analysis of the convergence by comparing the solutions obtained on two or more meshes, whose spacings differ by a factor of 2 to 5 times. Figure 6 compares the relative SIF values calculated on the meshes of 1000 ¡m x 600 ¡m, 47.5 ¡m x 27.5 ¡m and 21 ¡m x 5.5 ¡m using the J- and G-integral concepts. At the stages of active loading, the greatest difference in the SIF values is ~9 % for the meshes of 1000 ¡m x 600 ¡m and 47.5 ¡m x 27.5 ¡m, whereas for the mesh of 21 ¡m x 5.5 ¡m this differ- 1 0,8 cS 0,6 E LL 0,4 LL (0 0,2 0 1/7L Hi •3 /r6 /r5 50 100 150 200 Relative temperature, °C 250 300 Figure 6: Plots of the relative SIF values calculated on different meshes using the J- and G-integral concepts: 1 - 1000 ¡m x 600 ¡m (J), 2 - 1000 ¡m x 600 ¡m (G), 3 - 47.5 ¡m x 27.5 ¡m (J), 4 - 47.5 ¡m x 27.5 ¡m (G), 5 - 21 ¡m x 5.5 ¡m (J), 6-21 ¡m x 5.5 ¡m (G) Slika 6: Prikaz relativnih vrednosti SIF, izračunanih z različnimi mrežami, s konceptom J- in G-integrala: 1 - 1000 ¡m x 600 ¡m (J), 2 -1000 ¡m x 600 ¡m (G), 3 - 47,5 ¡m x 27,5 ¡m (J), 4 - 47,5 ¡m x 27,5 ¡m (G), 5 - 21 ¡m x 5,5 ¡m (J), 6-21 ¡m x 5,5 ¡m (G) ence is %. It can be noted that the SIF values determined using the G-integral concept turn out to be higher than those determined using the /-integral concepts, and also that with a decrease in the mesh-spacing value, the difference between the SIF values calculated with various methods decreases, too. 6.3 Determination of the allowable critical-brittleness temperature oof the metal The data of the SIF calculation for a circumferential under-the-cladding crack obtained using the /- and G-integral concepts is applied when determining the critical-brittleness temperature of the metal. A fracture-toughness curve is used, whose equation is given in3,9. Figure 7a shows the plots of the temperature dependence of the relative SIF values Kj(T)/Kmax and the fracture toughness KIc(T)/Kmax for the three meshes (1000 pm x 600 pm, 47 pm x 27.5 pm and 21 pm x 5.5 pm) obtained using the /-integral concept. On the other hand, Figure 7b illustrates the plots obtained using the G-in-tegral concept. Table 1 summarizes the values of increments in the critical-brittleness temperature of metal ATka compared for the meshes of 1000 pm x 600 pm, 47 pm x 27.5 pm and 47 pm x 27.5 pm , 21 pm x 5.5 pm. On the basis of the analysis of the results it is possible to draw a conclusion about their convergence when successively denser meshes are used. The critical-brittleness temperature values of the metal obtained using the /-integral concept turned out to be less conservative than those obtained using the G-in-tegral concept. The critical-brittleness temperature value of the metal increases and becomes determined with lesser conservatism when performing a refinement of a finite-element mesh. Table 1: Values of increments in the critical-brittleness temperature of the metal Tabela 1: Vrednost korakov kritične temperature krhkosti kovine Mesh, pm x pm A Tka, /- ATka, G- 47 x 27.5 ^ 1000 x 600 14 °C 14 °C 21 x 5.5 ^ 47 x 27.5 5 °C 3 °C 6.4 Analysis of the results obtained using the closed-crack model Figure 8 shows the results of determining the SIF values at a deep point of the crack and at the midpoint of Figure 7: Determining the critical-brittleness temperature value of the metal: a) the /-integral concept; b) the G-integral concept: 1 - 1000 pm x 600 pm, 2 - 47.5 pm x 27.5 pm, 3-21 pm x 5.5 pm Slika 7: Določanje vrednosti kritične temperature krhkosti kovine: a) koncept /- integrala, b) koncept G- integrala: 1 - 1000 pm x 600 pm, 2 - 47,5 pm x 27,5 pm, 3-21 pm x 5,5 pm Figure 8: Plots of the relative SIF values calculated on different meshes using the /- and G-integral concepts: a) is for the deep point of the crack, b) is for the point at the boundary between the base metal and the cladding: 1 - 1100 pm x 360 pm (/), 2 - 1100 pm x 360 pm (G), 3 - 100 pm x 100 pm (/), 4 - 100 pm x 100 pm (G), 5 - 12 pm x 12 pm (/), 6 - 12 pm x 12 pm (G) Slika 8: Relativne vrednosti SIF, izračunane z različnimi mrežami z uporabo koncepta /- in G-integrala: a) je za točko v globini razpoke, b) je za točko na meji med osnovnim materialom in oblogo: 1 - 1100 pm x 360 pm (/), 2 - 1100 pm x 360 pm (G), 3 - 100 pm x 100 pm (/), 4 - 100 pm x 100 pm (G), 5 - 12 pm x 12 pm (J), 6 - 12 pm x 12 pm (G) the boundary between the base metal and the cladding with the use of the J- and G-integral concepts. For the deep point, the results obtained using the G-integral concept, as before, turned out to be more conservative than those obtained using the J-integral concept. At the same time, for the point at the boundary between the base metal and the cladding, the results obtained using the J-integral concept turned out to be more conservative than those obtained with the G-integral concept. For this point, the maximum difference in the SIF values determined using both concepts was approximately ~5 %. 7 CONCLUSIONS The SIF determination methods based on the use of the J- and G-integral concepts were analyzed for a typical regime of the WWER-1000 RPV emergency cooling. When solving the problem of an elastic and elasto-plastic formulation without taking into account the loading history, the difference in the SIF values is negligible, corresponding to the theoretical statements about the equality of J and G integrals. The results of the SIF calculation obtained using the initial full model of the RPV - with a built-in crack and without it - practically coincide, which testifies to the uselessness of the crack embedding into the initial full model of the RPV. However, according to the requirements of the IAEA (International Atomic Energy Agency), a crack is to be built in the RPV fragment under consideration. It was shown that when solving the problem of the elasto-plastic formulation on the basis of the flow theory under non-proportional loading, the SIF values obtained using J- and G-integral concepts are in agreement even at the stages of unloading, although the use of the J-inte-gral concept under these conditions contradicts the theory statements. The results of the SIF calculation using the G-integral concept are more conservative than those obtained using the J-integral concept. The character of the temperature dependence of the fracture parameters - the stress at the crack tip and SIF - at the stages of unloading was shown to be considerably dependent on the density of the employed finite-element mesh. With a decrease in the size of the finite element in the vicinity of the crack front, an increase in the size of the compressive-stress zone is observed, and the descending branch of the SIF curve tends to become practically vertical. The analysis of the behavior of the fracture-parameter curves points to the convergence of the results, which is a necessary condition for solving the problem with the finite-element method. However, the use of the temperature dependence of SIF constructed from the solution of the problem on a dense finite-element mesh for the purpose of determining the brittleness temperature of the metal results in obtaining, as far as we know, less conservative results than those obtained using the curve constructed from the problem solution on a coarse finite-element mesh. 8 REFERENCES 1Z. Chen, The finite element method - Its fundamentals and applications in engineering, World Scientific, New Jersey 2011 2D. Kovacevic, I. Budak, A. Antic, B. Kosec, Special finite elements; theoretical background and application, Technical Gazette, 18 (2011) 4, 649-655 3 Authorizing Document of the Operating Organization 0606-2005, Procedure of calculating for brittle fracture resistance of WWER nuclear reactor pressure vessels (RPV - BFR - 2004), St. Petersburg, Moscow, 2004 4A.Yu. Chirkov, Mixed finite element method scheme for solving boundary-value problems of the elasticity theory and small elasto-plastic strains, Publishing House of the Institute for Problems of Strength, Kiev 2003 5G. P. Cherepanov, Brittle Fracture Mechanics, Nauka, Moscow 1974 6H. G. deLorenzi, On the energy release rate and the J-integral for 3-D crack configurations, Int. Journ. of Fracture, 19 (1982), 182-193 7 S. N. Atluri (ed.), Computational Methods in Fracture Mechanics, Mir, Moscow 1990 8 A. Yu. Chirkov, Development and implementation of the mixed finite element method in problems of vibration stability and structural element strength, PhD Thesis, Kiev, 2008 (in Russian) 9 Standard Document 0.03.391-09, Procedure for assessing the strength and in-service life of WWER RPV, Kiev, 2009