Microsoft Word - numero_34_art_25 T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 237 Focussed on Crack Paths Modelling 3D crack propagation in ageing graphite bricks of Advanced Gas-cooled Reactor power plant Thi-Tuyet-Giang Vo University of Manchester – The Modelling and Simulation Centre, UK. Philippe Martinuzzi, Van-Xuan Tran EDF Energy R&D UK Centre – The Modelling and Simulation Centre, George Begg Building, University of Manchester, M139PL Manchester, UK. Philippe.Martinuzzi@edfenergy.com Neil McLachlan, Alan Steer EDF Energy Generation ABSTRACT. In this paper, crack propagation in Advanced Gas-cooled Reactor (AGR) graphite bricks with ageing properties is studied using the eXtended Finite Element Method (X-FEM). A parametric study for crack propagation, including the influence of different initial crack shapes and propagation criteria, is conducted. The results obtained in the benchmark study show that the crack paths from X-FEM are similar to the experimental ones. The accuracy of the strain energy release rate computation in a heterogeneous material is also evaluated using a finite difference approach. Planar and non-planar 3D crack growth simulations are presented to demonstrate the robustness and the versatility of the method utilized. Finally, this work contributes to the better understanding of crack propagation behaviour in AGR graphite bricks and so contributes to the extension of the AGR plants’ lifetimes in the UK by reducing uncertainties. KEYWORDS. 3D crack propagation; Advanced Gas-cooled Reactor; Graphite brick; eXtended Finite Element Method; Code_Aster; Strain energy release rate. INTRODUCTION GR cores are composed of thousands of graphite bricks. After years of service, the heterogeneous irradiation and temperature fields gradually modify both the microstructure and internal stresses of the bricks. This modification may lead to crack initiation and propagation. One goal of the Plant Life EXtension (PLEX) program, developed by EDF Energy, aims at understanding the crack propagation in graphite bricks for the coming years. Part of the research program conducted at EDF Energy R&D UK Centre relies on studying the capabilities to model graphite brick crack propagation using the X-FEM in Code_Aster, an open source finite element software package developed by EDF R&D since 1989. X-FEM’s major advantage is that it represents the crack without explicitly meshing it. This method of crack propagation allows many steps of crack propagation to be performed using a single mesh and enables a more accurate calculation of the stress intensity factors and the strain energy release rate around the crack front. Recent developments in A T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 238 Code_Aster enable modelling of the behaviour of ageing graphite (spatial heterogeneity of the Young’s modulus and internal stresses). In this paper, crack propagation in AGR graphite bricks with ageing properties is studied using the X-FEM. The accuracy of the strain energy release rate computation in a heterogeneous material is evaluated using a finite difference approach. A parametric study for crack propagation, including the influence of different initial crack shapes and propagation criteria, is then conducted. The robustness and the consistency of the propagation are also presented. In Section 2, the implementation of X-FEM model is given. A description of the methodology for crack propagation in ageing graphite bricks is presented in Section 3. The benchmark study, the numerical simulation results for planar and non-planar crack problem on the graphite bricks appear in Section 4 with the discussions. And finally, in Section 5, some conclusions and prospects are provided. CRACK PROPAGATION USING X-FEM Basics of Level Set Method and X-FEM he Level Set Method (LSM) is the numerical scheme developed by Osher and Sethian [1] to track the evolution of interfaces and shapes. In the LSM the interface is represented as the zero level set of a function of one higher dimension. To describe a crack, two level set functions are required [2]. The normal level set (lsn) represents the signed distance to the crack surface and the tangent level set (lst) represents the distance to the crack front. The X-FEM allows modelling the crack growth without having to modify the original mesh. In cracking, the discontinuity of displacement due to the crack is introduced by a generalized Heaviside function and the addition of the asymptotic fields near the crack tip in order to improve the accuracy in elastic fracture mechanics. The formula used for the displacement approximation uh(x) at point x is given in Eq. (1). The first term is the standard one (continuous). The second term brings the discontinuity of displacement for the elements crossed by the crack (Heaviside enrichment) and the third term brings the singularities for the elements around the crack tip (asymptotic enrichment). Fig. 1a shows a representation of a crack with additional degrees of freedom for enriched elements. Note that the crack passes through the elements and that the crack is not explicitly meshed [4]. 4 ( ) ( ) ( ) 1 ( ) ( ) ( ) ( ( )) ( ) ( ( ), ( )) n n n h i i j j k k i N x j N x K k N x L u x a x b x H lsn x c x F lsn x lst x                 (1) where ai are the traditional displacement degrees of freedom associated with node i, Φ i is the linear shape function of node i, bj are the enriched degrees of freedom associated with node j, kc  are the enriched degrees of freedom associated with node k. The function H(x) is a Heaviside function, which is discontinuous across the crack surface. Consider an element that is crossed by the crack, H(lsn(x))=1 for the nodes located on the side where lsn(x)>0, whereas H(lsn(x))=-1 for those located on the other side where lsn(x)<0. Only one layer of elements is enriched. Note that a geometrical (fixed area) enrichment can also be done. The expressions of the function F  , α = 1.4, are the following [3]:  sin , cos , sin sin , cos sin2 2 2 2F          (2) The local polar co-ordinates (  ,  ) can be expressed in terms of the level sets at the crack tip as (see Fig. 1b): 2 2 1; tan lst lsn lst lsn            (3) The propagation of crack thus consists in updating the level sets and determining the new crack front. The mesh is refined in the area near the crack front (the distance is evaluated by the values of the level sets) to improve the accuracy of the stress intensity factors and the strain energy release rate computation. The flowchart for a Code_Aster model with X- FEM crack propagation is presented in Fig. 1c. T T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 239 (a) (b) (c) Figure 1: A representation of a crack with additional degrees of freedom for enriched elements (circle symbols represent the elements enriched by the Heaviside function and square symbols represent the elements enriched by the singular functions) (a), global Cartesian coordinate system (r,z), local Cartesian coordinate system (n,t), local polar coordinate system (  , ) and representation of the level sets (b) and flowchart for Code_Aster model with X-FEM crack propagation (c). [4] Criteria of propagation – The maximum hoop stress criterion Erdogan and Sih [5] proposed the maximum hoop stress criterion. The maximum hoop stress criterion assumes that the crack will propagate where the hoop stress  is maximised. Near the crack front, a good approximation of the hoop stress field is given by Eq. (4). 2 1 1 1 3 cos cos sin 2 2 22 I IIK K r                    (4) The angle of propagation is deduced by differentiating the Eq. (4) and is given in Eq. (5).   2 1 12 tan 8 4 I I II II II K K sign K K K                  (5) where KI, KII are the values of the stress intensity factors. For non-planar crack propagation, at each point along the crack front, a plane perpendicular to the crack front at this point will be considered. The directions of the non-planar crack paths are also computed using Eqs. (4) and (5). CRACK PROPAGATION IN AGEING GRAPHITE BRICKS Algorithms of crack propagation in ageing graphite bricks using X-FEM he analysis is done in two steps (see Fig. 2). The first one is the ageing analysis. In order to perform this calculation, the complex ABAQUS User MATerial (UMAT) routine, developed by Amec Foster Wheeler since the 1990s, was used (Code_Aster being able to directly read UMAT routines). From this calculation, it is possible to extract the different material properties at any instant tcrack which also vary in space. For this study, the assumption is made that the crack will propagate instantly, therefore the properties of graphite are not evolving during crack propagation and the behaviour can be considered elastic. The important values from the UMAT calculation are the Young’s modulus and T T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 240 the stress state at instant tcrack that will be considered as the initial stress state in the crack propagation analysis. These values are extracted and passed into the second step. This step is based on an algorithm where the LSM is coupled with the X-FEM to introduce and represent the crack and perform the crack propagation (see Fig. 2). Step 1 Step 2 Irradiation and oxidation of graphite Crack propagation (XFEM) Introduction of the crack (XFEM) Extraction of: Material properties: E(x,y,z,tcrack) Stress state: σ_ij(x,y,z,tcrack) Non linear behaviour (UMAT) Crack propagates when G=Gc Automatic mesh refinement Figure 2: Description of the analysis done for crack propagation in ageing AGR graphite bricks. The computation of the energy release rate G is done by using the G-theta method [7]. When G is smaller than Gc (the critical strain energy release rate), there is no crack propagation. When G reaches Gc, there is propagation. Heterogeneous Young’s modulus in space Regarding the heterogeneity of the Young’s Modulus, developments have been made in Code_Aster in order to enable the calculation despite some theoretical approximations: indeed, the spatial variation of Young’s Modulus in the component is taken into account but there is a missing term, in gradient of E, which is not calculated. In fact, the generic definition of G which is derived from the evolution of the potential energy (Ep) with the increase of the crack surface (A) with propagation is described in Eq. (6). pE G A     (6) Therefore, in a 2D model, by simply calculating the potential energy for a crack of length l and then for a crack of length l+dl (where dl is small compared to l), it is possible to approximate G with Eq. (7). This value will be called Gdiff. ( ) ( )p p diff E l dl E l G dl     (7) The previous results of Martinuzzi and al. [6] had shown that the missing term in the calculation of G appeared to be negligible in the considered case of AGR bricks. RESULTS AND DISCUSSIONS Benchmark of crack propagation tools on virgin graphite bricks ig. 3 represents the geometry of graphite bricks used in AGR. Depending on the reactors, the dimensions and the geometry may slightly vary. A benchmark was conducted in order to test the ability of four potentially interesting methods – X-FEM with Griffith adapted criteria, Cohesive Zone Models (CZM) [8], Damage modelling (these three first methods are developed in Code_Aster) and Configurational Forces (developed at the University of Glasgow) – to model crack propagation. Fig. 4 shows the results of the crack paths obtained from these four methods on the first and simplest case: a 100 mm slice of F T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 241 brick is slotted and externally loaded in order to replicate the internal stress state of graphite bricks after years of service and exposition to oxidation and irradiation. The crack path is two-dimensional and planar. Figure 3: Geometry of the graphite bricks. Section (a), lateral view (b), and representation of the methane holes (c). Figure 4: Crack propagation methods tested for AGR graphite bricks: X-FEM, CZM, Damage modelling and Configurational Forces. From the results obtained in the benchmark study, the most promising one developed in Code_Aster was the use of the X-FEM. In this case, the crack paths were similar to the experimental ones. Incremental crack growth was based on energetic criteria relevant to quasi brittle materials: the direction of propagation along the crack front was based on the maximum hoop stress criterion, and the length of propagation for each step of calculation depends on the value of the strain energy release rate along the crack front. The crack paths are presented in Fig. 5 for different loading scenarios using X-FEM, and compared with their experimental counterparts. Planar crack propagation in ageing AGR graphite brick The behaviour of these ageing graphite bricks has been studied for many years. Fast neutron irradiation causes graphite to shrink initially before expanding again to exceed its original volume. The attenuation in dose between the bore and the outside causes differential rates of graphite shrinkage that generate internal stresses that vary through life. Fast neutron irradiation also causes the elastic modulus of the material to increase significantly early on, while oxidation causes   a)  b) c)  T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 242 progressively increasing weight loss that reduces the elastic modulus. These changes occur more rapidly close to the fuel channel walls, so there is a significant variation in elastic modulus across the thickness of the brick. Many numerical analyses were performed using UMAT routines in ABAQUS. These UMAT routines use material properties that depend on three different field variables (temperature, irradiation dose, and weight loss due to oxidation). Developments in Code_Aster have been made to enable stress analysis to be run by directly reading the UMAT routines. Figure 5: Crack paths obtained for the different test cases of the benchmark using the X-FEM (left) and experimental method (right). Simulation of internal loading on slotted sliced bricks (a), simulation of external loading on sliced bricks (b), external loading on full- size bricks (c), and the Brokenshire test [9] (d). Figure 6: Evolution of the strain energy release rate (G) with crack propagation obtained by Code_Aster and ABAQUS. The values of G have been calculated in three different ways in Code_Aster, two local approaches (G0 and Gav) and a global approach (Gdiff). ABAQUS calculation is based on a global approach T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 243 A parametric analysis has also been conducted using X-FEM in Code_Aster on a full-size brick with initial stresses representing post stress-reversal state of an AGR graphite brick when there are tensile stresses at the keyway corners. The influences of initial crack shape, propagation criteria, and the choice of the fracture mechanics post-processing tools are studied. A first comparison with an ABAQUS model gives good agreement between both codes. Particularly, the evolution of the strain energy release rate with crack propagation is similar (see Fig. 6). Additionally, a good agreement between the local approach and the global approach is seen in Fig. 7. 0 2 4 6 8 10 12 14 16 0 0.01 0.02 0.03 0.04 0.05 0.06 E m e rg y  r e le a se  ( J) Crack area (m2) Code_Aster ABAQUS Figure 7: Evolution of the strain energy release rate with crack propagation obtained by Code_Aster and ABAQUS. Fig. 8 shows the full-size brick and the evolution of a crack path under planar propagation. Simulations are performed for 94 steps. The study stopped when the crack propagated fully along the brick height. All crack propagations steps were automatically performed in Code_Aster. It took 7 hours to perform the entire study. Step 2 Step 60 Step 94 Figure 8: Photograph of a fuel brick (a) and evolution of a planar crack path using X-FEM on a full-size graphite brick with stresses coming from UMAT analysis (b). Non- planar crack propagation in virgin graphite brick The non-planar crack propagation model is also conducted by using X-FEM. An initial elliptic crack, representing a large, postulated, initial defect, is introduced with initial stresses representing the post stress-reversal state of an AGR graphite a) b) T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 244 brick. The influences of mesh size, propagation criteria, and the choice of the fracture mechanics post-processing tools are studied. Crack propagation simulations are conducted, and the different stages in the crack growth paths are depicted in Fig. 9. The angle of propagation is defined by Eq. (5). Simulations are performed over 16 steps, until the crack reaches the interstitial keyway and reaches the outer boundary of graphite brick. The simulations reveal the non-planar character of propagation and the results clearly demonstrate the capabilities of X-FEM. Ongoing work will be to improve the robustness of 3D non-planar crack propagation in ageing AGR graphite bricks. Step 1 Step 4 Step 8 Step 12 Step16 Figure 9: Evolution of a non-planar crack path using X-FEM on a full-size graphite brick with stresses coming from UMAT analysis. CONCLUSION n this paper, numerical applications of X-FEM are presented to demonstrate the crack propagation modelling method in Code_Aster and to show the capabilities of solving the challenging problem in computational fracture mechanics in ageing AGR graphite bricks. The results obtained in the benchmark study shows that the crack paths from X-FEM were similar to the experimental ones. The accuracy of the strain energy release rate computation in a heterogeneous material is evaluated using a finite difference approach. The evolution of the planar crack path on an ageing AGR graphite brick and the non-planar crack path on virgin graphite brick reveal the robustness and the consistency of the propagation. Further work to improve the robustness of 3D non-planar crack propagation in AGR graphite bricks is ongoing. This work contributes to the better understanding of crack propagation behaviour in service, and so contributes to the extension of the AGR plants’ lifetimes in the UK by reducing uncertainties. ACKNOWLEDGEMENTS his work was done as part of a TSB collaborative project between EDF Energy R&D UK Centre, EDF Energy Generation, EDF R&D in France and the University of Manchester. The authors gratefully acknowledge EDF Energy Generation for their valuable input data and support. The views expressed in this paper are those of the I T T.-T.-G. Vo et alii, Frattura ed Integrità Strutturale, 34 (2015) 237-245; DOI: 10.3221/IGF-ESIS.34.25 245 authors, and not necessarily those of EDF Energy Generation. Technical support for this work from EDF R&D is greatly appreciated. REFERENCES [1] Osher, S., Sethian, J.A., Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics, 79 (1988) 12–49. [2] Stolarska, M., Chopp, D.L., Moës, N., Belytschko, T., Modelling crack growth by level sets in the extended finite element method, Int. J. Numer. Methods Eng., 51 (2001) 943–960. [3] Geniaut, S., Galenne, E., A simple method for crack growth in mixed mode with X-FEM, International Journal of Solids and Structures, 49 (2012) 2094–2106. [4] Tran, V-X, Geniaut, S., Development and industrial applications of X-FEM axisymmetric model for fracture mechanics, Engineering Fracture Mechanics, 82 (2012) 135–157. [5] Erdogan, F., Sih, G.C., On the crack extension in plates under plane loading and transverse shear, ASME Journal of Basic Engineering, 85 (1963) 519–527. [6] Martinuzzi, P., Pellet, L., Geoffroy, D., Modelling crack propagation in irradiated graphite, The 4th EDF Energy Nuclear Graphite Symposium: Engineering Challenges Associated with the Life of Graphite Reactor Cores, (2014). [7] Code_Aster, Analysis of structures and thermomechanics for surveys and research. [8] Park, K., Paulino, G. H., Cohesive zone models: a critical review of traction-separation relationships across fracture surfaces, Applied Mechanics Reviews, (2011) 64(6) 060802. [9] Brokenshire, D.R., A Study of Torsion Fracture Tests, Ph.D. Thesis, Cardiff University, (1996).