A Critical Heat Generation for Safe Nuclear Fuels after a LOCA

Abstract

This study applies a thermo-elasto-plastic-creep finite element procedure to the analysis of an accidental behavior of nuclear fuel as well as normal behavior. The result will be used as basic data for the robust design of nuclear power plant and fuels. We extended the range of mechanical strain from small or medium to large adopting the Hencky logarithmic strain measure in addition to the Green-Lagrange strain and Almansi strain measures, for the possible large strain situation in accidental environments. We found that there is a critical heat generation after LOCA without ECCS (event category 5), under which the cladding of fuel sustains the internal pressure and temperature for the time being for the rescue of the power plant. With the heat generation above the critical value caused by malfunctioning of the control rods, the stiffness of cladding becomes zero due to the softening by high temperature. The weak position of cladding along the length continuously bulges radially to burst and to discharge radioactive substances. This kind of cases should be avoid by any means.

1. Introduction

This study applies a thermo-elasto-plastic-creep finite element procedure to the analysis of an accidental behavior of nuclear fuels as well as normal behavior. The result will be used as a basic data for the robust design of nuclear power plant and fuels. We extended the range of mechanical strain from small or medium to large adopting the Hencky logarithmic strain measure in addition to the Green-Lagrange strain and Almansi strain measures for the possible large strain situation in accidental environments.

The accidental situation is one of the most severe cases where LOCA is involved simultaneously. In most LOCA simulation, the power is generated 100% at the initial state and assumed to be shut down suddenly through the reactor trip at the beginning of the accident. However, the case assumes an accident unusual to the normal LOCA without ECCS (event category 5) [1–3], and the power is not completely but partly shut down. Therefore, we may call it a kind of RIA compared with the nominal full shutdown. The coolant behaves also most severely to exhibit an instant loss of pressure and level as one might expect an explosive failure of reactor shell.

2. Finite Element Analysis of Nuclear Fuels

Finite element analysis procedures of nuclear fuels have been developed by many researchers [5–7] and mostly applied to basic designs of them. The nominal operation of the nuclear power plant is known to exhibit the highest pellet temperature between 1100~1400°C and the clad temperature 400~600°C [8]. In this study, we introduce the thermo-elasto-plastic-creep analysis procedure based on Green-Lagrange strain measure first. For the analysis of the fuels under accidental situation like LOCA, we extended the procedure to large strain capability adopting logarithmic strain measure and additive decomposition of plastic strain. Through the adoption of additive decomposition method, we could provide a continuation between small-medium strain analysis and large strain analysis including creep phenomena.

2.1. Thermo-Elasto-Plastic-Creep Analysis

Thermo-elasto-plastic-creep analysis procedure has been developed by many researchers [4, 6, 9–12]. Most of them adopted Green-Lagrange or Almansi strain measures adequate for small-medium strain. It is briefly introduced here. For the thermo-elasto-plastic-creep analysis of nuclear fuel, it is necessary to obtain time-dependent temperatures solving the following well- known transient equation of heat transfer: 𝑘𝑇_,𝑥𝑥+𝑇_,𝑦𝑦+𝑄−𝑐𝜌𝑇̇=0,𝑞=𝑓_𝐵onΓ_𝑓,𝑞=ℎ𝑇−𝑇_∞onΓ_ℎ, (1) where 𝑇= temperature, 𝑘 = coefficient of the thermal conductivity, 𝜌 = density, 𝑐 = specific heat, 𝑞 = heat flow rate, 𝑓_𝐵 = heat flux, and ℎ = convective heat transfer coefficient.

The partial differential equation of heat transfer (1) is expressed in other forms of functional to be minimized as in (2) [13].

Consider 1=2𝐓_𝜕^𝑇[𝜅]𝐓_𝜕𝑇̇+𝑓−𝑄𝑇+𝑐𝜌𝑇𝑑𝑉_𝐵𝑇+ℎ𝑇_∞1𝑇−2ℎ𝑇²𝑑𝑆, (2) where [𝜅] = conductivity matrix and {𝐓_𝜕} = temperature gradient.

After we introduce finite elements with proper shape functions and minimize the functional, the following finite element equation is derived: ^𝑡+Δ𝑡𝐂^{(𝑖)𝑡+Δ𝑡}𝐓̇^(𝑖)+^𝑡+Δ𝑡𝐊^𝐾+^𝑡+Δ𝑡𝐊^𝐶Δ𝐓^(𝑖)=^𝑡+Δ𝑡𝐐+^𝑡+Δ𝑡𝐐^{𝐶(𝑖−1)}−^𝑡+Δ𝑡𝐐^{𝐾(𝑖−1)},_𝑡+Δ𝑡𝐓^(𝑖)=^𝑡+Δ𝑡𝐓^(𝑖−1)+Δ𝐓^(𝑖). (3)(4)

Superscript 𝑖 indicates the quantity of iterative state 𝑖 at time 𝑡+Δ𝑡; therefore, ^𝑡+Δ𝑡𝐂 is the heat capacity matrix at time 𝑡+Δ𝑡, 𝐓̇^𝑡+Δ𝑡 is the rate of temperature at time 𝑡+Δ𝑡, ^𝑡+Δ𝑡𝐊^𝐾 is the conductivity matrix at time 𝑡+Δ𝑡, and ^𝑡+Δ𝑡𝐊^𝐶 is the convection matrix at time 𝑡+Δ𝑡. Δ𝐓 is the incremental nodal point temperature, ^𝑡+Δ𝑡𝐐 is the nodal point heat input vector, ^𝑡+Δ𝑡𝐐^𝐶 is the nodal point heat contribution due to convection, and ^𝑡+Δ𝑡𝐐^𝐾 is the equivalent nodal point heat contribution. Accordingly, to solve the unknown nodal point temperature, ^𝑡+Δ𝑡𝐓, the equilibrium iteration will be performed repeatedly until the difference in the right hand terms of (3) is smaller than the convergence tolerance.

After obtaining the temperature distributions of fuel, the incremental equilibrium equation for the quasi-static motion of the fuel can be solved to obtain the incremental displacements: ^𝑡+Δ𝑡𝐊^(𝑖−1)Δ^𝑡+Δ𝑡𝐔=^𝑡+Δ𝑡𝐏−^𝑡+Δ𝑡𝐑^(𝑖−1), (5) where ^𝑡+Δ𝑡𝐊^(𝑖−1): stiffness matrix, Δ𝐔^(𝑖): incremental displacement, Δ^𝑡+Δ𝑡𝐔^(𝑖): ^𝑡+Δ𝑡𝐔^(𝑖−1)+Δ^𝑡+Δ𝑡𝐔^(𝑖), ^𝑡+Δ𝑡𝐏: external load, ^𝑡+Δ𝑡𝐑^(𝑖−1): equivalent load, and all quantities ^𝑡+Δ𝑡𝐗⁽⁰⁾=^𝑡𝐗.

The stiffness ^𝑡+Δ𝑡𝐊^(𝑖−1) has linear term and geometrically nonlinear term. The temperatures are included in the external load resulting in the thermal deformation, and the thermal stain derived from the thermal deformation is considered in the evaluation of mechanical strain to compute the mechanical stresses.

The deviatoric stress is computed from the deviatoric elastic strain disregarding plastic and creep strain as ^𝑡+Δ𝑡𝐒=^𝑡+Δ𝑡𝐸1+^𝑡+Δ𝑡𝜈^𝑡+Δ𝑡𝐞^′−^𝑡+Δ𝑡𝐞^𝑃−^𝑡+Δ𝑡𝐞^𝐶. (6) The mean stress is computed from the mechanical mean strain as in (5).

Consider ^𝑡+Δ𝑡𝜎_𝑚=^𝑡+Δ𝑡𝐸1−2^𝑡+Δ𝑡𝜈^𝑡+Δ𝑡𝑒_𝑚−^𝑡+Δ𝑡𝑒_th, (7) where ^𝑡+Δ𝑡𝐒 = deviatoric stress tensor, ^𝑡+Δ𝑡𝐞^′ = deviatoric strain tensor, ^𝑡+Δ𝑡𝐞^𝑝 = plastic strain tensor, ^𝑡+Δ𝑡𝐞^𝑐 = creep strain tensor, ^𝑡+Δ𝑡𝜎_𝑚 = mean stress, ^𝑡+Δ𝑡𝑒_𝑚 = mean strain, ^𝑡+Δ𝑡𝑒_th = thermal strain, ^𝑡+Δ𝑡𝐸 = Young’s modulus, and ^𝑡+Δ𝑡𝜈 = Poisson's ratio.

If we split the plastic and creep stains at 𝑡+𝑑𝑡 into the terms at time 𝑡 and incremental terms, (7) becomes ^𝑡+Δ𝑡𝐒=𝐸1+𝜈^𝑡+Δ𝑡𝐞^′′−Δ𝐞^𝑝−Δ𝐞^𝑐, (8) where ^𝑡+Δ𝑡𝐞^′′=^𝑡+Δ𝑡𝐞^′−^𝑡𝐞^𝑝−^𝑡𝐞^𝑐. (9)

Creep strain may be represented by various creep laws, such as power creep law or exponential creep law. To confirm plastic stress based on the given displacement, ^𝑡+Δ𝑡𝑢^(𝑖), it must satisfy one of the yield criteria. Adopting the von Mises yield criterion, 𝑓_𝑦=^𝑡+Δ𝑡𝜎²−^𝑡+Δ𝑡𝜎²_𝑦=0, (10) where ^𝑡+Δ𝑡𝜎=1.5^𝑡+Δ𝑡𝐒⋅^𝑡+Δ𝑡𝐒:vonMisesequivalentstress,_𝑡+Δ𝑡𝜎_𝑦=yieldstress. (11)

We have to solve the equation to obtain the deviatoric stress tensor. By applying the following flow rules of plastic and creep deformations, Δ𝐞^𝑝=Δ𝜆^𝑡+Δ𝑡𝐒,Δ𝐞^𝑐=Δ𝑡^𝜏𝛾^𝜏𝐒,𝜏=(1−𝛼)𝑡+𝛼(𝑡+Δ𝑡):relaxedtime𝛼=relaxationparameter,Δ𝜆=1.5Δ𝑒^𝑝_𝑡+Δ𝑡𝜎:incrementalparameterforequivalentplasticstrainandequivalentstressat𝑡+Δ𝑡,𝛾=1.5Δ𝑒^𝑐_𝜏𝜎:parameterforrelatingincrementalequivalentcreepstrainandequivalentcreepat𝑡+Δ𝑡. (12)

We obtain ^𝑡+Δ𝑡𝐒=_𝑡+Δ𝑡1𝑎_𝐸+𝛼Δ𝑡^𝜏𝛾+Δ𝜆^𝑡+Δ𝑡𝐞^′′−(1−𝛼)Δ𝑡^𝜏𝛾^𝑡𝐒, (13) where 𝑎_𝐸=(1+𝜈)𝐸. (14)

Taking product of both sides of the equation, we obtain the equation 𝑓_𝑔=𝑎^2𝑡+Δ𝑡𝜎²+𝑏^𝜏𝛾−𝑐^2𝜏𝛾²−𝑑²=0, (15) where 𝑎=^𝑡+Δ𝑡𝑎_𝐸+𝛼Δ𝑡^𝜏1𝛾+Δ𝜆=2𝐺^′,𝑏=3(1−𝛼)Δ𝑡^𝑡+Δ𝑡𝐞^′′⋅^𝑡𝐒,𝑐=(1−𝛼)Δ𝑡^𝑡𝑑𝜎,₂=1.5^𝑡+Δ𝑡𝐞^′′⋅^𝑡+Δ𝑡𝐞^′′. (16)

The function, left side of (15), is known as the effective stress function and has to satisfy the condition of zero. Because 𝜎 is function of plastic strain and Δ𝜆 is also a function of plastic strain and yield stress, the equation has to be solved simultaneously for 𝜎 and Δ𝜆 using a root finding algorithm. We tried both Newton method and bisection method. The latter method is very reliable without being affected by the derivative of function. For the assumed ^𝑡+Δ𝑡𝜎=^𝑡+Δ𝑡𝜎_𝑦, we compute Δ𝜆 first and check whether the equation is satisfied. If not, follow the usual root finding algorithm to find the ^𝑡+Δ𝑡𝜎=^𝑡+Δ𝑡𝜎_𝑦 and Δ𝜆 simultaneously that satisfy the chosen type of yield criterion.

Here, we propose a first stage constitutive modulus as 𝐶^′′=_𝑡+Δ𝑡1𝑎_𝐸+𝜆+Δ𝑡^𝜏𝛾=2𝐺^′′; (17) see Figure 2.

The tangent constitutive matrix 𝐂^EPC can be obtained using the perturbation method [4, 14]. Each column of the matrix 𝐂^EPC is the column of perturbed stress tensor divided by the magnitude of perturbed element of corresponding strain. However, at the first stage of interaction, one can use quasi-elastic modulus of Figure 1(b) [4], where 𝐶^′′ is more close to the tangent modulus than 𝐶^′. We could find fast convergence with 𝐶^′′ than with 𝐶^′ in several examples.

(a) Elastic

(b) Quasi-elastic [4]

Figure 1 

First stage constitutive modulus.

Figure 2 

Proposed first stage quasi-elastic modulus.

2.2. Large Strain (Hencky Logarithmic) Measure

For large strain analysis, we need some basic relations that are briefly introduced here [9, 10, 15–18]. Deformation gradient is represented by ^𝑡₀𝐗=𝜕^𝑡𝐱𝜕⁰𝐱=^𝑡𝐉^𝑇0𝐉^−1𝑇, (18) where 𝐉 is the transformation between the cantesian coordinate and natural coordinate systems. It can be decomposed as ^𝑡₀𝐗=^𝑡₀𝐑^𝑡₀𝐔 (19) into rotation tensor ^𝑡₀𝐑 and right stretch tensor ^𝑡₀𝐔. Right Cauchy-Green gradients are expressed using the deformation gradient ^𝑡₀𝐂=^𝑡₀𝐗₀^𝑇𝑡𝐗. (20)

Solving the following equivalent problem, ^𝑡₀𝐂^𝑡₀𝐏=^𝑡₀𝚲₀^2𝑡𝐏. (21)

We obtain a diagonal matrix^𝑡₀𝚲 of principal stretches and their direction matrix. The right stretch matrix is obtained as ^𝑡₀𝐔=^𝑡₀𝐏^𝑡₀𝚲^𝑡₀𝐏^𝑇. (22)

Therefore, the rotation matrix is given by ^𝑡₀𝐑=^𝑡₀𝐗^𝑡₀𝐔⁻¹. (23)

Then, the Henckey logarithmic strain tensor is given by ^𝑡₀𝐄_𝐻=12^𝑡0𝐏^𝑡₀𝐇^𝑡₀𝐏^𝑇, (24) where ^𝑡₀𝐇=ln𝚲². (25)

We can also represent the Green-Lagrange strain tensor as ^𝑡₀𝐄_𝐺=12^𝑡0𝐂−𝐈 (26) and the Almansi strain tensor as ^𝑡₀𝐄_𝐴=12𝐈−⁰_𝑡𝐂. (27) The Henckey logarithmic strain tensor can also be written as ^𝑡₀𝐄_𝐻=ln^𝑡₀𝐔=12ln^𝑡₀𝐂. (28) We computed the above equations and compared them (three kinds of marks) with theoretical (29) in Figure 3: ^𝑡₀𝑒₁₁=12^𝑡0𝜆²−1:Green-Lagrangestrain,_𝑡𝑒₁₁=^𝑡₀𝜆:Smallstrain,_𝑡0𝐸^(𝐻)₁₁=ln^𝑡₀𝜆:Logarithmicstrain,_𝑡𝑡𝑒^𝐴₁₁=121−^𝑡₀𝜆⁻²:Almansistrain. (29)

Figure 3 

Strain-stretch relations for various strain measures.

Using the above procedure of thermo-elasto-plastic-creep with logarithmic strain measure with additive decomposition of plastic strain, our additive decomposition method follows the procedures in [9, 10, 15] except that we do not decompose plastic deformation gradient multiplicatively as 𝐗^𝐸_𝐴=₀^𝑡+Δ𝑡𝐗^𝑡₀𝐗^𝑝−1 (30) at the trial elastic state first but considered it additively as in (6) and (8). Therefore, we do not loose continuation between the later large strain procedure and the former small strain procedure including creep deformation.

The case without creep deformation in the following example compares the two procedures of large strain problems, additive decomposition and multiplicative decomposition. The first one is uniaxial isochoric extension with the following condition [15]:  𝐸 = 200 MPa: Young’s modulus;  𝜈 = 0.3: Poisson’s ratio;  𝜎_𝑦 = 0.75 MPa: yield stress;  𝐻 = 2 MPa: hardening modulus.

The results are compared in Figure 4 showing a good agreement.

Figure 4 

Normal stress 𝜎𝑥 along stretch 𝜆.

The second example is shear deformation; the result is compared with the result of the literature (solid square) obtained by the multiplicative decomposition in Figure 5.

Figure 5 

Shear stress 𝜏12 along shear strain 𝛾.

We also compared the following example (Figure 6) under shear load in Figure 7, which tests not only the stress integration process but also the equilibrium iteration procedure too.

Figure 6 

FEM model for shear load analyses of multiplicative and additive decomposition methods.

Figure 7 

Shear stresses under shear loads for different decompositions of plastic strain.

3. Behavior of Nuclear Fuels after Loss of Cooling Accidents

In the usual case of LOCA, the control rods are inserted into the fuel rods to shut down the reactor. The power is reduced successfully. However, such as the Fukushima accident, the control rods could not be inserted to the rods because of the malfunction of inserting system located at the lower part of the reactor. Therefore, we assumed that the heat generation of the fuels is partly reduced and partly working. In the literature concerning LOCA scenario, we saw a 100% power generation adopted in the analysis.

Before we analyze the behavior of nuclear fuels after LOCA (loss of coolant accident), we analyzed normal operation first using the following condition. The material properties of fuel components are the same as in the literature [6]. Mesh generation of the fuel rod was conducted using the quadrilateral 8-node elements, and the boundary conditions were applied as in Figure 8.

Figure 8 

Solid nuclear fuel model and boundary conditions.

The boundary conditions include the application of heat generation to the pellet and the application of coolant pressure and convection resulting from the coolant to cladding. By applying gas pressure to the gap between the pellet and cladding, the overall environment when actual nuclear fuel rods are activated is shown in Table 1.

Table 1 

Applied condition.

Figure 9 shows the temperature changes at pellet center and outer surface of the cladding along time, and temperatures go into steady state before 100 sec.

Figure 9 

Temperatures changes at pellet center and cladding outer surface along the time.

Suppose that there is a LOCA (loss of coolant pressure and relative heat transfer coefficient changed from 7500 mJ/s⋅mm²⋅^∘K to 75 mJ/s⋅mm²⋅^∘K) at time 70 sec. Moreover, there is an additional situation of the malfunctioning of control rod.

We assumed that one of four conditions of heat generation (complete shutdown = 0%, 55%, and 60% and full heat generation = 100%) is maintained after the LOCA. Figure 10 shows temperature changes at the two positions. As one may expect, shutdown makes the rapid drops of temperature at both positions of fuel. For the intermediate heat generations, the temperatures of pellet drop first and then increase afterwards. The temperature of clad increases fast and slowly increases later on. Finally, sharp increase of temperature is detected at the clad for full heat generation. The asterisk marks denote that the analyses have been stopped due to the instability of clad stiffness caused by the softening by temperature increase, surpassing the hardening by strain. The load applied to the clad due to the rod internal pressure is supported by the strength of clad. However, if the temperature of the clad is too high, Young’s modulus and strength of the clad are reduced, and it cannot sustain the applied load. Once it cannot resist the load and starts to deform, the situation is out of control and the clad becomes thinner and thinner causing larger and larger stresses resulting in fuel failure. This phenomenon is a kind of instability like a buckling.

Figure 10 

Temperatures changes at pellet center and cladding outer surface after LOCA at 70 sec.

The variations of von Mises stress of clad are shown in Figure 11. The heat generations above 55% induced instability at the clad to cause continuous expansion to burst at the weakest point along the length.

Figure 11 

Von Mises stresses changes at pellet center and cladding outer surface after LOCA at 70 sec.

In the early stage, there is a gap between the pellet and clad. After a long time of operation, the coolant pressure compresses the clad to cause a creep phenomenon and the gap is gone. Then the contact stress between pellet and clad raises the compressive stress at the inner part of the pellet. If this stress is above the ultimate strength of the pellet material, pellet may be fractured to pieces and may give some harm to the clad.

Hoop stresses at the center of pellet and clad surface are shown in Figure 12. We may have temporal time interval for the rescue of the plant in accident for the condition of heat generation under 55% of full generation. Over 60%, a catastrophic situation occurs when the clad of zero stiffness expands continuously to bust and discharge radioactive substances, a situation we have to avoid by any means.

Figure 12 

Hoop stresses changes at pellet center and cladding outer surface after LOCA at 70 sec.

Figure 13 exhibits the tendencies of thermal stresses in the pellet. Except the very early stage, outer portion of the pellet experiences tensile hoop stress, and the inner portion of it experiences compressive stresses. Every case including full heat generation shows lowering of thermal stresses in the pellet, which is thought to originate from smaller temperature difference between the extreme positions and insulation effect at the outer surface of the clad due to the loss of coolant. What is ironic is that the LOCA results in the reduction of temperature gradient in the pellet and causes more even distribution of temperature within the pellet to make it safer.

Figure 13 

Compressive and tensile stresses changes at pellet center and outer surface after LOCA at 70 sec.

4. Conclusions

We extended the finite element procedure of code TEPC (thermo-elasto-plastic-creep) from small and medium strain to Hencky logarithmic strain measure for the large strain analysis to handle the case that may happen in accidental situation. The TLPC (thermo-large strain-plastic-creep) which is the extended version of the TEPC has been applied to the analysis of nuclear fuels under LOCA without ECCS (event category 5) as well as normal operation.

A LOCA is concerned with the system pressure decrease, the rod internal pressure, the hoop stress on the cladding, the temperature of the cladding, strain, and so forth. We added an accidental condition of malfunctioning of the control rod system resulting in incomplete shutdown of the fuel power. And we checked a critical remaining heat generation that causes clad failure. The results are summarized as follows.(i)The fuels under consideration have reached quasi-steady state condition after 70 sec. The maximum core temperature was 1631°K, and the outer cladding temperature was 612.6°K when there was no contact between the pellet and cladding.(ii)We adopted a LOCA without ECCS (event category, 5), and the external pressure of the fuel rod has zero gauge pressure. After the LOCA (lower heat transfer coefficient of 75 mJ/s⋅mm²⋅^∘K) at time 70 sec with the heat generation of the pellet under the critical value, there was no instability at the cladding for the time being, which means one may have temporal time interval for the rescue of the system by recovery of the cooling system or shutdown of the system letting the heat generation be zero.(iii)After the LOCA at time 70 sec, there was an unstable condition that the stiffness of the cladding becomes zero which is followed by a local uncontrollable expansion of clad to burst, with above critical heat generation of the pellet.

For the analysis of this kind of unstable situation, we need certain special analyzing methods. One method is a postbuckling method usually applied to static problems, and it is expected to be improper for this kind of transient problem. The other one is a dynamic analysis, which usually is not necessary for this kind of quasi-static problem; however, it is expected to be the most effective technique which can handle the problem of zero stiffness.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2012R1A1A2008903). The authors are partially supported by BK21 plus project of the National Research Foundation of Korea.