Modelling fracture of aged graphite bricks under radiation and temperature

The graphite bricks of the UK carbon dioxide gas cooled nuclear reactors are subjected to neutron ir- radiation and radiolytic oxidation during operation which will affect thermal and mechanical material properties and may lead to structural failure. In this paper, an empirical equation is obtained and used to represent the reduction in the thermal conductivity as a result of temperature and neutron dose. A 2D ﬁnite element thermal analysis was carried out using Abaqus to obtain temperature distribution across the graphite brick. Although thermal conductivity could be reduced by up to 75% under certain conditions of dose and temperature, analysis has shown that it has no signiﬁcant effect on the temperature distribution. It was found that the temperature distribution within the graphite brick is non-radial, different from the steady state temperature distribution used in the previous studies [1,2] . To investigate the signiﬁcance of this non-radial temperature distribution on the failure of graphite bricks, a subsequent mechanical analysis was also carried out with the nodal temperature information obtained from the ther- mal analysis. To predict the formation of cracks within the brick and the subsequent propagation, a linear traction–separation cohesive model in conjunction with the extended ﬁnite element method (XFEM) is used. Compared to the analysis with steady state radial temperature distribution, the crack initiation time for the model with non-radial temperature distribution is delayed by almost one year in service, and the maximum crack length is also shorter by around 20%.


Introduction
Nowadays, in the UK, the gas-cooled graphite moderated type dominates over 90% of the nuclear reactor [3] . Graphite is also considered to be used within the design of the new reactor generation (such as the Helium gas cooled very high temperature reactor or VHTR) as structural support and as moderator because of its extreme purity, ability to withstand extremely high temperatures [4] . The core structure of the reactor is formed from joining graphite bricks, with keys and keyways. During operation, graphite is exposed to neutron radiation and temperature gradients that could cause irradiation and radiolytic oxidation. These occurrences will affect material properties such as weight loss, porosity changes and thermal conductivity which can potentially lead to material and structural failure [5] .
Future reactors such as VHTR are expected to operate at higher temperatures than the current reactors and its core outlet temperature (COT) will be about 10 0 0 °C [5][6][7][8] . In order to enhance thermal efficiency of future nuclear reactors, extensive research has been carried out, aiming to provide an improved design for future reactors [6][7][8][9][10] . As the reactor is exposed to neutron dose and temperature, many material properties are changed by irradiation and radiolytic oxidation if air or carbon dioxide is used as a coolant. Some of the properties of nuclear graphite which change under service conditions include elastic properties and thermal expansion coefficient, irradiation creep, irradiation induced dimensional changes. The prediction of the nuclear graphite lifetime and its structural integrity are significant issues for the safety and reliability of reactor operation that could be affected by changing properties.
In previous studies [2] , a constitutive model has been developed for the effects of mechanical properties changes on stress states within the graphite brick. Resultant crack development from the stress state within the brick was also studied in [1] . For previous studies, steady state temperature distribution was used to calcu- late thermal stress assuming that thermal conductivity of graphite brick remains constant across the brick during the operation. It was also assumed that the temperature distribution is linear across the graphite brick. However, studies [11,12] have shown that thermal conductivity of nuclear graphite changes with both temperature and fast neutron irradiation. Hence, a new thermal analysis is required to calculate the temperature gradient across the brick by considering the changes in thermal conductivity due to dose and temperature.
In this paper, an empirical equation for the evolution of thermal conductivity is derived from available experimental data by Marsden et al. [12] . Using this relationship with relevant thermal properties, a thermal analysis is carried out. Resultant nodal temperatures are then implemented into mechanical finite element (FE) model with a constitutive model developed in [2] and fracture mechanics based damage model developed in [1] . Stress analyses and crack development are then carried out. Comparisons are also made to models with the linear temperature distribution across the brick.

Empirical relationship of thermal conductivity of graphite with neutron dose and irradiated temperature
The effects of temperature gradients and irradiation can induce damage to the graphite bricks and this leads to changes in its thermal conductivity. The common approach is based on relating thermal conductivity (K) to summation of thermal resistances (1/K) as a result of scattering hurdles [13,14] where the term α( x ) is a coefficient, which is function of orientation terms (basal plane, x ), 1

K GB
is the thermal resistance due to grain boundary scattering, 1 K u is the thermal resistance due to Umklapp scattering or the so-called phonon-phonon scattering [14] , and 1

K RD
is the thermal resistance due to lattice defect scattering. K GB depends on the graphite perfection and is taken into account only for lower temperatures, as it is insignificant at intermediate and high temperatures. K u can be scaled as a quadratic function of temperature [14] and it controls thermal conductivity at higher temperatures. Finally, K RD depends on the neutron irradiation, as the latter plays a significant role in producing various types of defects, which are more effective in scattering phonons. Therefore, many studies [5,13,15] emphasise that the change in conductivity is only caused by the phonon scattering in lattice defects. This phonon scattering increases with temperature. Hence it is necessary to consider both irradiation dose and temperature to estimate the change in thermal conductivity. For significant nuclear irradiated graphite grades e.g. Gilsocarbon, Kelly [11] proposed a formula as shown in Eq (2) , to calculate the thermal conductivity. It includes the terms, f and S k , obtained empirically as a function of dose and temperature [12] . 1 K( 30 ) − 1 and S k is a structural term which allows decrease in thermal conductivity at higher dose. Marsden et al. [12] used Eq (2) on irradiated and unirradiated graphite samples. The thermal conductivity is shown to be directly proportional to the operation temperature (T) for irradiated graphite but the opposite is true for unirradiated graphite. Without using an atomistic physical approach for calculating of the variation of thermal conductivity for the graphite such as the one proposed in [16] , the thermal conductivity dependence on dose and temperature can be deducted from experimental data curves or include empirical terms obtained from the experiments using Material Test Reactor (MTR). In the current paper, thermal conductivity variation results of MTR from Marsden et al. [12] are used to generate the empirical equation. The variation in thermal resistivity with temperature and dose is illustrated in Fig. 1 . The data was obtained from the nuclear experiments with equivalent DIDO nickel dose (EDND). The neutron dose is expressed in terms of equivalent DIDO nickel dose (EDND) which is calculated so that the damage rate experienced by the graphite at a given location is equivalent to that of graphite test samples irradiated in the DIDO reactor [17] . 1 n / c m −2 EDND is equal to 1 . 313 × 10 −21 displacement per atom (dpa) [12] .
Using the data of thermal resistivity from Fig. 1 , an empirical equation shown in Eq (3) for the dependence of the resistivity on dose (D) and temperature (T) is fitted.
where K is the thermal conductivity of the nuclear graphite, D is dose × 10 20 n / c m −2 EDND, T is the irradiation temperature in Kelvin and K o is the initial thermal conductivity of unirradiated graphite at 30 °C (about 130 W/m/K), and a, b and c are constants equal to 0.004185, 0.9716 and 209.4, respectively. Please note that the data for 25 °C was excluded from the fitting process to reduce the error sum of squares (SSE) from 3.742 to 1.858, with coefficients of 95% confidence bounds. The resultant Eq (3) can be used to calculate transient change in thermal conductivity for the thermal FE model for which changes in dose and temperature are expected. Finally, to illustrate the accuracy of the calculated thermal resistivity using Eq (5) , the values are plotted against experimental results as shown in Fig. 2 at 100 °C, 300 °C and 600 °C.

Geometry and mesh of finite element model of graphite brick
In this paper, failure of the graphite brick of Gilsocarbon under irradiation and temperature is estimated. In previous studies [1,2] , an assumed profile of temperature distribution across the graphite brick was applied to the model as one of the field variables to model the reactor conditions. The temperature is assumed to vary in the radial direction of the brick while it is kept constant in the circumferential direction and along the depth direction of the brick. Here, a separate thermal FE model is used to assess the effect of conductivity change on the temperature distribution within the brick. The empirical equation (presented in Section 2 ) of the evolution of the thermal conductivity as a result of irradiation and temperature change is used for the thermal model. Resultant nodal temperature distribution is imported to the mechanical FE model for stress and failure analyses.
The FE model is developed here based on the dimensions of a particular brick model used by Tsang and Marsden [2] . In order to save computation time, only one-eighth of hypothetical-brick design was modelled, as shown in Fig. 3 . The fillet radius at the corner of the keyway is 6 mm. Appropriate symmetry and circumferential periodic boundary conditions were applied along B1 and B2 (see. Fig. 3 (a)) as described by Zou et al. [24] . The FE mesh of the brick section used for both thermal and mechanical analyses is shown in Fig. 3 (b).The mesh is refined at the stress concentration region (at the keyway fillet), but larger element sizes were used for regions far from that interest, as shown in Fig. 3 (b).

Thermal finite element model
The heat transfer analysis from the thermal FE model provides the temperature distribution nodal field variables across the brick which can be implemented into the mechanical FE model as the predefined nodal fields. Fig. 4 shows a schematic illustration of a thermal FE model of the graphite brick. A steady-state heat transfer was applied with inner and outer temperature boundary conditions on the brick surfaces, T 1 = 500 °C and T 2 = 400 °C, respectively. Besides, the convection heat transfer from the inner and outer surfaces can be simulated by using the coefficient of heat transfer of the inner (h 1 ) and outer (h 2 ) surfaces assuming that the coolant is CO 2 . The initial thermal conductivity of the brick is taken as 130 W/m/K [18] . The plane strain elements of four-node linear diffusive heat transfer (DC2D4) were used.
For the particular section of the brick shown in Fig. 4 , the governing equations of the heat transfer are: where Q is the heat flux W/m 2 , T is the temperature, R Total is the total thermal resistance, R Outer is the convection resistance between the outer surface of the graphite brick and the coolant,  R Inner is the convection resistance between the inner surface of the graphite brick and the coolant, R Brick is the conduction resistance within the brick, A is the area of conduction or convection (assumed identical and unity) and X is the radial distance between the inner and outer surfaces. R inner and R outer are not considered for the current model to avoid a complex forced cooling model of flowing CO 2 coolant. In contrary, nodes at the inner and outer surfaces of the brick are fixed at 500 °C and 400 °C, respectively, by assuming that the CO 2 coolant will keep this temperature difference at the steady state. Therefore the change in the temperature within the brick will be given solely by the change in R Brick due to the change in its conductivity, K.
The dose is implemented as one of the field variables of the thermal model. It is ramped upwards from zero at the beginning of the reactor life. At the end of operation (30 years) the dose varies quadratically across the brick from 214 × 10 20 n/cm 2 at the inner radius to 106 × 10 20 n/cm 2 at the exterior. With the known information of dose and temperature, updated thermal conductivity is calculated using Eq (3) and stored as the field variables which are implemented in the mechanical FE model.

Mechanical finite element model
With regard to the mechanical model, a user-defined material model for the constitutive stress/strain behaviour of the graphite from Tsang and Marsden [2] was used. In the material model of graphite brick, different physical effects such as changes of graphite dimension, porosity and microstructure from the effect of irradiation, radiolytic oxidation and fast neutron damage are taken into account. This model was then adapted to investigate the behaviour of the brick fracture using a traction-separation cohesive model in [1] . The FE simulation was set for 30 years of plant operation. In this paper, temperature field is calculated from the separate thermal FE model and its effects on crack initiation and propagation within the graphite brick was investigated.
An ABAQUS solver with a user material subroutine (UMAT) was used to calculate the resultant stresses in the mechanical model. Fig. 5 summarises the mechanical FE model based on four field operating variables. In each analysis, four field-variables profiles (the operating and irradiation temperatures, neutron fluence and weight loss) are implemented to resemble the condition of reactor operation. The operating temperature is the same as the irradiation temperature during steady state whereas the operation temperature is ramped down to the room temperature at cooling down step after operating. The nodal temperature from the thermal analysis is implemented into the mechanical model. Fluence (radiation dose) variation is mentioned in Section 2 . An empirical relationship from Eason et al. [18] was used for the weight loss and it can be expressed as a function of the radiation dose as shown in Eq. (9) . However, the weight loss resulting from the thermal oxidation was neglected as thermal oxidation is not significant for lower temperature CO 2 cooled UK graphite moderators. It could be significant for new higher temperatures moderators. The coolant used in the simulation was CO 2 coolant, as radiolytic oxidation is more important than other types, but modifications to the material data gained from Tsang and Marsden [20] are needed for future generation reactors.
At the beginning of each increment, ABAQUS calculates the total strain within the model. The total strain is decomposed into estimated strain components related to different physical processes. There are seven strain components decomposed from the total strain ( ɛ T ), as shown in Eq (10) . ε T = ε e + ε pc + ε sc + ε dc + ε th + ε ith + ε idc (10) where ɛ e , ɛ pc , ɛ sc , ɛ dc , ɛ th , ɛ ith , ɛ idc are the elastic, primary creep (recoverable), secondary creep (unrecoverable), dimensional change, thermal, interaction thermal and interaction dimensional change strains, respectively.
From Hooke's law of linear elasticity, the ɛ e is calculated. The UKAEA model of Kelly and Brocklehurst [19] was used to model the creep behaviour. The ɛ pc and ɛ sc are irradiation induced creep strains and were assumed to be temperature independent. The ɛ pc leads to fast deformation and then gradually decreases with increasing irradiation fluence, followed by the ɛ sc , which is a steady state creep strain. The ɛ dc was assumed to be dependent on the irradiation fluence, temperature and weight loss, whereas the ɛ ith and ɛ idc were assumed to be dependent on irradiation creep in graphite. These interaction strains were proposed by Kelly and Burchell [20] as a correction factor to accommodate the change made by the total creep strain on CTE of graphite and dimensional change strain. These different strains create a discrepancy of various stresses that in turn creates a complex stress field. The UMAT then calculates the total stresses, updates them and prepares the Jacobin matrix (( ∂ σ ij )/( ∂ ɛ ij )) at the end of the increment to support convergence of solution. The ( ∂ σ ij ) and ( ∂ ɛ ij ) are changed in terms of stress and strain, respectively.
ABAQUS FE code then uses the current stresses in the linear damage model to predict initiation of cracks and then their propagation. Based on the stress field created, crack initiation and propagation can be estimated using a fracture mechanics based damage model.

Fracture mechanics based damage model
Graphite is a quasi-brittle material and the initiation of crack is assumed to be driven by the stress states within the material. Different failure criteria can be used for crack initiation (maximum principal stress, maximum nominal stress etc.). Maximum principal stress criterion is used in the current work and for this criterion, a crack initiates once the stress within the system exceeds the tensile strength of the graphite. For virgin graphite, tensile strength of 22.2 MPa is taken from Shi et al. [21] . Due to the lack of tensile strength evolution data with irradiation dose and weight loss, the evolution is assumed to be similar to the trend of flexural strength given in [22] . With this assumption, the tensile strength evolution with the weight loss is tabulated in Table 1 . Shi et al. [21] reasoned that the initial increase in strength is the result of hardening caused by neutron radiation although the overall reduction in strength is observed as the weight loss becomes higher than 5%. Numerically, extended finite element method (XFEM) is applied here for crack initiation for the present work. It allows crack initiation within arbitrary elements once the failure criterion is met. The direction of crack path is assumed to be perpendicular to the direction of the maximum principal stress and thus the current model simulates Mode I crack growth.
Graphite does not fail instantly when the maximum principal stress exceeds tensile strength like pure brittle materials. It has a strain softening region because of micro-cracking and grain bridging near the crack tip. A cohesive traction-separation damage model is applied to simulate this type of behaviour. It is a phenomenological model by combing complex damage mechanisms around the crack tip into a single cohesive zone. The approach has been used by other researchers [23,24] for fracture analysis of nuclear graphite. The traction stress (T) at the cohesive zone can be related to the crack separation ( δ) as shown in Fig. 6 . If graphite is a pure brittle material, it will fail when T exceeds the maximum traction stress ( T ), which is the same as the tensile strength. However graphite has a strain softening region once T is reached which is illustrated by the line XZ in Fig. 6 . The gradient of the line OX (K) represents the cohesive stiffness of an undamaged material. Once the crack initiates the stiffness reduces to K(1-D). A damage parameter (D) is set so that D = 0 and D = 1 represent zero damage at point X and a fully cracked condition at point Z respectively. It can be proved that D at point Y ( Fig. 6 ), can be linked to maximum crack separation ( δ z ), separation at maximum traction stress ( δ), and at the unloading point Y ( δ Y ) as shown in Eq (11) . The critical crack length opening is described in terms of the fracture stress ( T ) and the fracture toughness ( K c ) by Eq (12) . The area underneath the triangle OXZ ( Fig. 6 ) represents the critical strain energy release rate. Once this energy is achieved, the cohesive zone is completely damaged and the crack tip propagates further. There is no experimental data for the dependence of critical strain energy on the weight loss due to irradiation dose. Hence it is assumed here that the critical strain energy is dependent only on irradiation and the data from Ouagne et al. [25] are applied. (11)

Temperature distribution through the graphite brick
The thermal analysis model shown in Fig. 4 was implemented to analyse the heat flux and temperature distributions within the brick. Fig. 7 shows the heat flux through the brick and heat flux concentration can be observed at the fillet corner. Generally, it can be seen that heat is transferred from the hotter inner core side of the brick to the outer surface of the brick. Nodal temperature field from the thermal model is implemented into the mechanical model. In contrary to the assumption made by Tsang and Marsden [2] , a nonlinear change of the temperature along the radial direction of the brick is found. This is expected since the surface around the fillet region is assumed to have the same temperature (400 °C) as the outer surface of the brick at the other regions for the current model. Hence it has lower temperature than Tsang and Marsden model for which temperature varies linearly across the brick. Comparison of temperature distributions between the current model and the previous model by Tsang and Marsden is shown in Fig. 8 . Nodes along the path B of the brick (see Fig. 9 ), have been chosen to demonstrate the temperature variations under the irradiation for the current model and Tsang and Marsden model.
Although the thermal conductivity could be reduced by as much as 75%, the updated thermal conductivity does not have a noticeable effect on the temperature distribution. In other words, temperature distribution does not vary with the time chosen for the current model (30 years of irradiation). This may be because of a high initial thermal conductivity of Gilsocarbon (130 W/m/K) and also due to the short radial distance (150 mm) of the graphite brick (shown as X in Eq (7) ) through which heat is transferred. Nevertheless, the proposed method for the evolution of thermal conductivity could be important in other aspects of heat transfer analyses. This may include a detailed transient heat transfer analysis for start-stop load cycles or the case where heat is conducted through a larger distance of the graphite brick (axial length of the brick).

Stress distribution in the graphite brick
Since the maximum principal stress criterion is used for crack initiation within the graphite brick, it is worthwhile to investigate the distribution of stresses across the brick before implementing the damage model. Cracks are expected to initiate at the region where peak maximum principal stresses occurs. The current work simulated the stresses through the brick under two conditions; the radial and non-radial temperature distributions as shown in Fig. 8 . From both temperature distributions, maximum principal stress concentration can be found at the fillet corner of the graphite brick at the end of 30 year of service life as shown in Fig. 10 . Hence, the crack initiation is expected at the fillet corner. For both temperature distribution cases there is no significant difference in contours for the maximum principal stress (see Fig. 10 ).
The evolutions of principal stresses for both temperature distribution cases at the fillet corner during operation are plotted in Fig. 11 . There is an increase in the principal stress at the start of operation for both cases until the stress stabilises after 5 year of operation. Sudden increase in stresses after the steady state can be found at around 20 year life. The sudden increase in maximum principal stress can be related to the increase in the hoop strain related to dimensional change of graphite at higher dose as mentioned by Wadsworth et al. [1] . Final increase in the principal stress during the cooling stage after operating for 30 years is due to the mismatch in coefficient of thermal expansion (CTE) throughout the brick. In general, higher maximum principal stress is expected for a linear temperature distribution cases indicating the crack initiation time for the case could be shorter than the model with nonlinear temperature distribution.

Damage analysis of the graphite brick under temperature and radiation loads
According to the stress analysis from the previous section, cracks are expected to form at the fillet corner due to the concentration of maximum principal stress. The crack initiation time for linear temperature distribution (22 years and 7 months) is shorter than the cases with nonlinear temperature distribution (23 years and 2 months) because of higher maximum principal stress. The finding shows the importance of the temperature distribution for the lifetime estimation of the graphite bricks. The typical geometry of the crack in the brick after 30 year operation for linear and nonlinear temperature profiles are shown in Fig. 12 (a) and (b), respectively. Firstly, the crack propagates towards the inner surface and then changes its direction when reaches a distance of 6.5 mm. The direction of further crack growth is perpendicular to the hoop direction indicating that the maximum principal stress is in the direction of hoop when the crack reaches 6.5 mm. The total crack length for both cases at the end of operation is around 65.4 mm. When crack growth rate is observed, it is shown that the growth rates for both temperature distribution cases are similar as shown in Fig. 13 . Rapid crack growth occurs within 1 year from crack initiation time. Because of this crack, stresses are relaxed and the crack length reaches the plateau stage after one year from the initiation time. Further increase in crack length is found after cooling down stage. During the operation, the brick is exposed to different temperature and radiation dose. Hence, difference in CTE across the brick is expected and this CTE mismatch causes stresses when the cooling is applied. This mismatch stress further increases the crack length by 1.5 mm. The maximum principal stress distribution at the end of cooling down for both temperature distribution cases are shown in Fig. 14 .

Conclusions and future work
The empirical formula of varying thermal conductivity of common irradiated graphite form, Gilsocarbon, due to the change in temperature and dose has been developed. This formula is developed by using the experimental thermal conductivity of irradiated graphite from the literature. The resultant formula can be implemented into either thermal FE model or corresponding mechanical FE model and life estimation models.
Although, the thermal conductivity of graphite could be reduced by 65-75% under irradiation and temperature, the thermal analysis showed it does not have a significant effect on the temperature distribution because of the high initial thermal conductivity and the relatively small thickness of the graphite brick investigated. However, it could play a significant role in the applications when heat is conducted through much larger distance, for example   through the thickness of larger graphite bricks. The thermal analysis also showed that the temperature distribution does not vary radially within the graphite brick due to the effect of the thermal concentration in the fillet corner. Based on the temperature distribution obtained from the thermal analysis, stress and damage analysis of the graphite brick under neutron dose and temperature were carried out. The dose increases the weight loss resulting in a decrease in the strength of the brick ahead the maximum principal stress and thus leading to   crack initiation within the brick. The initiation of cracks occurs at the fillet corner of the graphite brick and it grows towards the inner bore of the brick. The maximum principal stress is lower at the expected crack initiation region for the case where the tem-perature distribution is nonlinear compared to the case with assumed linear temperature distribution. Hence the crack initiation for the case is delayed by about 7 months. Regardless of crack initiation time, both cases have similar crack growth rate and the crack growth stops after one year from initiation time. Further crack growth is only driven from CTE mismatch stresses during the cooling down stage.
The study shows the importance of temperature distribution on the lifetime estimation of graphite bricks. Subtle difference in temperature distribution (as shown by two different cases in the current work) could cause the difference in the crack initiation time. However, there are some areas of the current work that needs to be improved. A transient heat analysis considering the flow of the coolant may be implemented in future instead of using a simple temperature boundary conditions for inner and outer surfaces of the graphite brick. In addition, the usage of a more realistic 3D geometry of the brick with cut-out holes for the circulation of the coolant may be considered in order to predict more accurate transient heat flow and the temperature profiles. Moreover, changes in porosity and phases in graphite as described by Kyaw et al. [26] could affect its mechanical properties and this model could be coupled to the current mechanical model in the future work.