A new procedure for solving steady-state and transient-state nonlinear radial conduction problems of nuclear fuel rods

Annals of Nuclear Energy 110 (2017) 492–500

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

A new procedure for solving steady-state and transient-state nonlinear radial conduction problems of nuclear fuel rods Jiannan Tang a, Mei Huang a,b,⇑, Yuanyuan Zhao a, Xiaoping Ouyang b, Jijuan Huang a a b

North China Electric Power University, Beijing 102206, China Northwest Institute of Nuclear Technology, Xi’an 71002, China

a r t i c l e

i n f o

Article history: Received 13 December 2016 Received in revised form 25 May 2017 Accepted 26 May 2017 Available online 19 July 2017 Keywords: Half-boundary method Nuclear fuel rod Nonlinearity Steady-state Transient-state

a b s t r a c t Half-boundary method (HBM) is extended to solve one-dimensional radial conduction problems in nuclear fuel rods. The accuracy and efficiency of HBM are first tested by solving three one-dimensional radial conduction problems in cylindrical coordinates. They are steady-state and transient-state problems with linear and nonlinear material properties. HBM shows a comparable accuracy to the analytical solutions, and a higher accuracy than finite volume method (FVM) especially with fewer elements in the discretized model. Besides, the computation time ratio of HBM to FVM decreases as the element numbers increase, and it is only 0.16% when 1000 elements are involved. In order to further validate the practical applicability of HBM, nonlinear heat conduction problems of nuclear fuel rods (NFRs) are solved by HBM. The numerical results of NFRs are carried out both in steady-state and transient-state cases. HBM has a matched result with ANSYS. In addition, the size effect of pellet, helium gap and zircaloy cladding on the highest temperature in NFRs is also studied using HBM. The pellet dimension affects the highest temperature mostly. Ó 2017 Published by Elsevier Ltd.

1. Introduction The temperature distribution analysis of the nuclear fuel rod (NFR) is one of the most important parts of nuclear safeties. With the sustained requirements of accurate and efficient numerical solutions to the heat transfer problems, more and more attention has been attracted to the innovation of the new computational methods. Currently, finite difference method (FDM) (Wang et al., 2009; Zhang et al., 2014), finite volume method (FVM) (Krishnaprakas and Narayana, 2001; Lyra et al., 2005; Todd et al., 2001; Chen et al., 2008; Versteeg and Malalasekera, 1995) and finite element method (FEM) (Duda, 2015; Loh et al., 2007) are the major methods to solve heat transfer problems. But all these methods have drawbacks in balancing calculation accuracy and efficiency. Besides, most of the practical problems are complex nonlinear problems which require the solving of nonlinear heat transfer equations where nonlinearity arises with the dependence on the physical properties associated with an unknown temperature. Particularly, nonlinear material properties influence the temperature distribution in the fuel rod a lot due to the giant temperature ⇑ Corresponding author at: North China Electric Power University, Beijing 102206, China. E-mail address: [email protected] (M. Huang). http://dx.doi.org/10.1016/j.anucene.2017.05.061 0306-4549/Ó 2017 Published by Elsevier Ltd.

gradient. To solve the nonlinear problems, the common method is known as Picard iteration. Nonlinearity is treated iteratively by repeatedly updating physical properties for the newest temperature solutions until achieving the convergence in both temperature and physical properties. Some other methods were also utilized to solve nonlinearity. Tomatis (2013) applied the Kirchhoff transformation method which integrated thermal conductivity to avoid the iterative scheme, and the produced temperature distributions were comparable to the analytical solutions. In our current work, half-boundary method (HBM) coupled with iterative scheme is adopted to solve the radial conduction problems of NFRs with nonlinear material properties. HBM is a method aiming at reducing the order of the partial differential equations. It handles the higher-order derivatives in the partial differential equations as independent variables. For example, for a one-dimensional problem, only the second-order matrices are involved in HBM, and the matrix order is free from the influence of the boundary condition types and the discretized models. The building of the complex matrices and the arduous matrix inversion required in FVM (Krishnaprakas and Narayana, 2001; Lyra et al., 2005; Todd et al., 2001; Chen et al., 2008) are hence successfully avoided. Also, the temperature and heat flux fields can be solved simultaneously. In addition, HBM does not have accuracy limit on few element numbers, which is a drawback in FEM (Zhang et al.,

493

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500

Nomenclature g T, X1 V, X2 cp

q r Dr t

Dt k

heat source term temperature (K) heat flux, J/(m2 s) heat capacity, J/(kgK) density, kg/m3 space coordinate, m distance between adjacent nodes, m time, s

Subscripts i grid node j iteration step Superscripts m time moment

2014; Duda, 2015; Loh et al., 2007) and FVM. Different from FDM (Wang et al., 2009), HBM ensures the local and global energy conservation due to an integral conservation over each element, similar as FVM. Therefore, HBM provides much higher calculation accuracy and efficiency as compared with FVM, FEM and FDM. In addition, due to the feature of HBM, HBM can solve the problems which have the two boundaries (temperature and heat flux) on only one side of the geometric models because the calculation marches from one side to the other. Moreover, no uniform grid nor certain element length is required. HBM can calculate the variable fields at any point of interest in the domain without solving the results for all nodes, which makes it more flexible. In fact, HBM was first proposed to conduct free vibration analyses (Li, 2005; Sakeyama and Huang, 2000; Huang et al., 2015). HBM was then extended to solve heat transfer problems in Cartesian coordinates (Huang et al., 2017). In this paper, HBM is first applied to solve heat conduction problems of NFRs in cylindrical coordinates. Three simple one-dimensional radial conduction problems were first studied to verify the accuracy and efficiency of HBM in cylindrical coordinates. The results were compared with those obtained by FVM and analytical solutions. The NFR model was then built, and HBM was utilized to solve the radial temperature distributions in NFRs in steady and transient states with linear and nonlinear material properties. All numerical calculations were carried out using Python 2.7 codes with double precision in which the Digits environment variable was assigned to be 64, and all the calculations were performed on a PC with a 3.60 GHz CPU and 8.0 GB of RAM. ANSYS was utilized to test the accuracy of HBM in the NFR model. Compared with FVM and ANSYS (FEM), HBM shows superiority in both efficiency and accuracy. This study lays a foundation for the further application of HBM in solving intricate models coupled with complex geometries, huge mesh numbers and multi-physics interactions.


 1 @ @T þ gðrÞ ¼ 0 rkðT; rÞ r @r @r

(


 @T 1 @ @T ¼ rkðT; rÞ þ gðr; tÞ @t r @r @r

ð1Þ

where qðT; rÞ, cp ðT; rÞ and kðT; rÞ are the parameters of material properties which depend on the temperature and location, while the source term gðr; tÞ depends on the location and time.

@ ðrVÞ @r

þ rgðrÞ ¼ 0

ð3Þ

V ¼ kðT; rÞ @T @r

To solve the steady-state heat conduction problems of NFRs, the second-type boundary condition (heat flux equals to zero) was adopted on the pellet central axis of NFR, and the third-type boundary condition (heat convection) was used between the outer-surface of the cladding and the coolant. The one-dimensional discretized domain is illustrated in Fig. 1, where Dr i indicates element i. In consideration of energy conservation, Eq. (3) is integrated from the arbitrary node i over Dri to the adjacent node i + 1, and thus Eq. (3) becomes:

8R < r i þDr i ri : R r i þDr i ri

R r þDr þ rii i rgdr ¼ 0 R r þDr V dr ¼ rii i @T dr kðT;rÞ @r

@ ðrVÞdr @r

ð4Þ

The derivatives are integrated as

R ri þDri ri

@V @r

R ri þDri ri

@T @r

dr ¼ T ðiþ1Þ  T i and

dr ¼ V ðiþ1Þ  V i , where T i and V i are the temperature and

heat flux at node i, respectively. Eq. (4) can be hence written as:

8 < r ðiþ1Þ V ðiþ1Þ  r i V i þ 12 ðr i þ r ðiþ1Þ Þg i Dr i ¼ 0   Vi : 12 Vkðiþ1Þ þ Dr i ¼ T ðiþ1Þ  T i ðTÞ k ðTÞ i i

ð5Þ

where ki ðTÞ is the thermal conductivity between node i and node i + 1, and is a function of temperature. By rearranging the variables at node i + 1 at next step to the left side of Eq. (5), Eq. (5) becomes:

(

In this work, the one-dimensional HBM formulae are deduced in cylindrical coordinates. The heat conduction equation in cylindrical coordinates can be written as:

ð2Þ

Reduce the order of Eq. (2) by setting V ¼ kðT; rÞ @T (the variable @r V is heat flux) and then Eq. (2) is changed to:

2. Half-boundary method

qðT; rÞcp ðT; rÞ

time step, s Coefficient of thermal conductivity, W/(mK)

rðiþ1Þ V ðiþ1Þ ¼ r i V i  12 ðr i þ r ðiþ1Þ Þg i Dr i T ðiþ1Þ þ 12

V ðiþ1Þ ki ðTÞ

Dr i ¼  12

Vi Dr i ki ðTÞ

ð6Þ

 Ti

Further simplify Eq. (6) into a matrix form:

"

0 rðiþ1Þ ri 1  12 kDðTÞ i

#

T ðiþ1Þ V ðiþ1Þ

"

 ¼

0

#

ri

ri 1  12 kDðTÞ i

Ti Vi

"

 þ

 12 ðr i þ r ðiþ1Þ Þg i Dr i

#

0 ð7Þ

Then Eq. (7) can be transformed to the expression below:

X ðiþ1Þp ¼

2 X apki  X ki þ bip

p ¼ 1; 2

ð8Þ

k¼1

2.1. Steady-state ¼ 0, For steady-state heat conduction problems, qðT; rÞcp ðT; rÞ @T @t the governing Eq. (1) becomes:

r i Dr i r ðiþ1Þ ki ðTÞ

where X p1 ¼ T i , X i2 ¼ V i , a11i ¼ 1, a12i ¼ 12 ri , b1i ¼ a22i ¼ rðiþ1Þ

Dr 2 ðrðiþ1Þ þr i Þ  4ri ðiþ1Þ gi , ki ðTÞ

b2i ¼

r þri  ðiþ1Þ 2rðiþ1Þ

þ 12

Dr i , ki ðTÞ

a21i ¼ 0,

g i Dr i . Eq. (8) shows the

relationship between node i and node i + 1.

494

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500

Fig. 1. The schematic of one-dimensional discretized domain.

Based on Eq. (8), we can also obtain the relationship between the arbitrary node i and the adjacent node i-1 as well as the relationship between node i-1 and node i-2. By iteration, we can get the relationship between node i and node i-2. By repeating this process, the relationship between the boundary node 1 and the arbitrary node i can be obtained as below:

X ip ¼

2 X cpki  X 1k þ dpi

ðp ¼ 1; 2

1 6 iÞ

ð9Þ

k¼1

 i P h r n D r n r1 rn , c21i ¼ 0, c22i ¼ rr1i , where c11i ¼ 1, c12i ¼ in¼1 12 rnþ1 þ 12 kDn ðTÞ knþ1 rn n h io P Pi1 Dr 2m ðrðmþ1Þ þr m Þ rðmþ1Þ þrm r n Dr n Dr n r 1 d1i ¼ i1 , m¼1  4r ðmþ1Þ km ðTÞ g m  4rðmþ1Þ g m Dr m n¼m ðr ðnþ1Þ kn ðTÞ þ kn ðTÞÞ rn Pi1 rnþ1 þrn d2i ¼ n¼1  rnþ1 g n Drn .

For transient-state condition, the time term is included. Similar to the steady-state formulae derivation, the order of Eq. (1) is also , and then Eq. (1) is simplified as: decreased by setting V ¼ kðT; rÞ @T @r

qðT; rÞcp ðT; rÞ @T ¼ 1r @r@ ðrVÞ þ gðr; tÞ @t V ¼ kðT; rÞ @T @r

ð10Þ

8 R r þDr R tm þDtm R r þDr R tm þDtm @V i i > qðTÞcp ðTÞri @T dtdr ¼ rii i tm r @x dtdr > tm @t > ri < R ri þDri R tm þDtm þ ri rgðtÞdtdr tm > > > : R ri þDri @T dr ¼ R ri þDri V dr @x

4

qi ðTÞci ðTÞ

2

ri

t ðmþ1Þ

tm

Z

t ðmþ1Þ

tm

" þ

qi ðTÞci ðTÞ 2

54 3" 5

ðmþ1Þ

T ðiþ1Þ

ðmþ1Þ

V ðiþ1Þ

ðmþ1Þ

Ti

3 5

# ð14Þ

ðmþ1Þ

Vi

m m 1 m Dri ðr iþ1Þ T m ðiþ1Þ þ r i T i Þ þ 2 r i þ r iþ1Þ Þg i Dr i Dt

#

In this paper, the relationship between node i and node i + 1 at the next time step t + D t is expressed as below:

X ðiþ1Þp ¼

2 X m ðmþ1Þ apki  X ik þ bpi

ðp ¼ 1; 2

i P 1Þ

k¼1

where

Dr 2i ðr ðiþ1Þ þ r i Þqi ðTÞci ðTÞ ðr ðiþ1Þ þ ri Þqi ðTÞci ðTÞDr 2i  ki ðTÞDtm rðiþ1Þ þ

a12i ¼

8Dtm rðiþ1Þ ki ðTÞ ; ðr ðiþ1Þ þ ri Þqi ðTÞci ðTÞDr 2  8ki ðTÞDt m r ðiþ1Þ

4Dri Dtm r ðiþ1Þ ðr ðiþ1Þ þ ri Þqi ðTÞci ðTÞDr 2i  ki ðTÞDtm rðiþ1Þ 

4Dtm rðiþ1Þ Dr i ; ðr ðiþ1Þ þ ri Þqi ðTÞci ðTÞDr 2i  8ki ðTÞDt m r ðiþ1Þ

a21i ¼ 

2ki ðTÞðrðiþ1Þ þ ri Þqi ðTÞci ðTÞDr i ðr ðiþ1Þ þ r i Þqi ðTÞci ðTÞDr 2i  8ki ðTÞDt m r ðiþ1Þ



2kðiþ1Þ qi ðTÞci ðTÞDr i ðr ðiþ1Þ þ ri Þ ; m ðr ðiþ1Þ þ ri Þqi ðTÞci ðTÞDr m i  ki ðTÞDt r ðiþ1Þ

a22i ¼

ð13Þ

The different values of h account for the different weights at tðmþ1Þ and tm moments. When h ¼ 0, the average variables only depend on current moment t m which is called the pure explicit

32

0

8ki Dtm ri ðr ðiþ1Þ þ ri Þqi ðTÞci ðTÞDr 2i  8ki Dtm r iþ1 

Tdt ¼ ½hT ðmþ1Þ þ ð1  hÞT m Dt m

Dxi kðiþ1Þ ðTÞ

1 Dr i 2 ki ðTÞ

ð11Þ

ð12Þ

 12

1

k

Vdt ¼ ½hV ðmþ1Þ þ ð1  hÞV m Dt m

1

r i Dtm

ues as: T ¼ hT ðmþ1Þ þ ð1  hÞT m , V ¼ hV ðmþ1Þ þ ð1  hÞV m where h 2 ½0; 1. Therefore, the integral of T and V over time can be written as:

Z

riþ1Þ Dt m

i ðTÞ r i Dr i  qi ðTÞc 2

a11i ¼

To solve the transient-state heat conduction problems of NFRs, the second and third-type boundaries mentioned in the steadystate part are involved. In addition, the initial condition when time = 0 s is obtained based on the steady-state calculations. Similar to the steady-state formulae derivation, Eq. (10) is integrated from location ri to r(i+1) and from time tm to t(m+1). In the proR tðmþ1Þ @T dt ¼ T mþ1  T m . Besides, V and T are cess of integrating, tm @t functions of time. A weighted interpolation of variables at two adjacent time moments is adopted to approximate the average val-

r iþ1Þ Dri

2

ðmþ1Þ

By integrating Eq. (10) over the element Dri and the time interval Dtm which is the time step from t m to t ðmþ1Þ , the following expression is obtained:

ri

2

¼4

2.2. Transient-state

(

Euler method. As for h ¼ 1, the average variables equal to the values at next time moment tðmþ1Þ which is called the pure implicit Euler method. For h ¼ 1=2; it is the Crank-Nicolson method (Crank and Nicolson, 1947). The explicit Euler method has strict time step restriction (He and Li, 2010), while the implicit Euler method has no time step restriction. In this paper, implicit Euler method (h ¼ 1) is used in HBM. Substituting Eqs. (12) and (13) into Eq. (11), and then discretizing Eq. (11), and further simplifying it into a matrix form, similar to the progress for the steady-state problems, we obtain the following matrix form:

m

b1i ¼

qi ðTÞci ðTÞDrmi ðrðiþ1Þ þ ri Þ ; ðr ðiþ1Þ þ ri Þqi ðTÞci ðTÞDr 2i  ki Dt m r ðiþ1Þ Dr

i ðq ðTÞci ðTÞDr i T m ðiþ1Þ qi ðTÞci ðTÞDr2i  ki Dtm rðiþ1Þ i m m þ qi ðTÞci ðTÞDr i T m i þ g i Dt Dr i Þ;

ð15Þ

495

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500 m

b2i ¼

2ki ðq ðTÞci ðTÞDr i T m ðiþ1Þ qi ðTÞci ðTÞDr2i  8ki Dtm rðiþ1Þ i þq

m i ðTÞci ðTÞDr i T i

þ

m gm i Dt Dr i Þ:

Similar to the steady-state process, based on Eq. (15), the relationship between the arbitrary node i and the boundary node 1 at next time step t ðmþ1Þ can be obtained.

ðmþ1Þ

X ðiþ1Þp ¼

2 X ðmþ1Þ cm þ qm pki  X 1k pi

p ¼ 1; 2

ð16Þ

p ¼ 1; 2 ; d ¼ 1; 2

ð17Þ

k¼1

where

cm pdi ¼

2 X m am pki ckdði1Þ k¼1

qm pi ¼

2 2 X X m m am bipk  g m pki qkði1Þ þ ik k¼1

p ¼ 1; 2

ð18Þ

k¼1

3. Verification results and discussion In this section, the accuracy and efficiency of HBM were verified before we actually solved the temperature distribution in NFRs. Three typical one-dimensional radial heat conduction problems were hence first solved. The numerical results of FVM and the analytical solutions (exact solution) were also obtained to make a comparison with those of HBM. The computation time of HBM and FVM were also recorded. The geometric model of the three one-dimensional heat conduction problems is presented in Fig. 2. For problem 1, the cylinder is hollow in the center, and is constant in the coefficient of thermal conductivity (k = 0.5 W/(mK)). The temperature of the inner face is 200 K, while that of the outside surface is 100 K. The radii of inner and outer surfaces are rin = 0.05 m and rout = 0.25 m, respectively. 3.1. Steady-state linearity

Fig. 2. The geometric sketch for one-dimensional radial heat conduction problems.

Fig. 3 presents the numerical results obtained by HBM and FVM as compared with the analytical solution for different element numbers. The numerical result obtained by HBM matches the

Fig. 3. A comparison between the numerical results of HBM and FVM, as well as the analytical solution for different element numbers of problem 1.

496

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500

Table 1 Numerical results of HBM and FVM for different element numbers. Element number

Location

Analytical solution

FVM

HBM

Result

Error

Result

Error

4

0.02 0.06 0.10

156.93 131.74 113.86

151.61 124.73 105.91

3.39% 5.32% 6.98%

155.45 130.69 113.37

0.95% 0.79% 0.44%

12

0.02 0.06 0.10

156.93 131.74 113.86

155.35 129.41 111.03

1.01% 1.77% 2.49%

156.74 131.61 113.80

0.12% 0.10% 0.05%

100

0.02 0.06 0.10

156.93 131.74 113.86

156.76 131.47 113.52

0.11% 0.21% 0.30%

156.93 131.74 113.86

0.00% 0.00% 0.00%

1000

0.02 0.06 0.10

156.93 131.74 113.86

156.90 131.68 113.80

0.02% 0.04% 0.06%

156.93 131.74 113.86

0.00% 0.00% 0.00%

Table 2 Computation time ratios of HBM to FVM for element numbers changing from 4 to 10,000. Element number

4 20 100 400 1000 5000 10,000

Matrix Order

Time ratio (HBM/FVM)

FVM

HBM

4 20 100 400 1000 5000 10,000

2 2 2 2 2 2 2

95.37% 74.81% 39.51% 11.82% 4.16% 0.46% 0.16%

analytical solution perfectly even when few elements are involved in the calculation. Moreover, HBM shows a much higher accuracy than FVM although the accuracy of FVM increases with the increased element numbers. The numerical results and computation errors of HBM and FVM for 4, 12, 100 and 1000 elements are listed in Table 1. The errors produced by HBM are much smaller than those by FVM, particularly when the model has few elements. When 100 elements are involved in the calculation, the error of HBM has become zero, while for FVM the error still exists even though the element number has increases to 1000. In fact, the error of FVM cannot be eliminated due to the boundary approximation from the adjacent nodes. Different from FVM, the heat flux value at the boundary of HBM is accurate, and thus HBM has good accuracy even few elements are included. The computation time ratio of HBM to FVM for different element numbers of example 1 is presented in Table 2. As the element number is small, the computation time of HBM and FVM is similar. However, with the increase of elements HBM obviously requires less time in carrying out the calculation than FVM. When the element number is above 5000, the computation time of HBM is less than 0.5% of that of FVM. With the increase of elements, the numerical calculations of HBM and FVM both require matrix inversion. As can been seen from Table 2, the matrix order of FVM keeps the same as the element number, and thus the CPU cost needed for building the complex matrix and for solving the inverse matrix increases correspondingly. On the contrary, the calculation time of HBM is much less since its matrix order is constant 2. This indicates that the computation time of FVM increases proportionally to the element numbers while that of HBM does not increase much. Therefore, in addition to the high accuracy, HBM has high efficiency especially for complex problems involving a great number of grids.

Fig. 4. The calculation flow chart of HBM with nonlinear material properties.

3.2. Steady-state nonlinearity In practice, the thermal conductivity coefficients of most materials are varied, especially when the temperature gradient is big. The inconsistent material properties influence the temperature distribution a lot. Hence, the nonlinearity of material properties is considered in this section. The same geometric model illustrated in Fig. 2. is used. But the temperature of the inner surface is changed to 1000 K, and that of the outside surface is still 100 K. In this problem, the thermal conductivity coefficient is assumed to be a linear function of temperature:

kðTÞ ¼ 0:5 þ 0:001T

ð19Þ

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500

Fig. 5. Numerical results obtained by HBM and FVM with linear and nonlinear material properties.

Table 3 Numerical results of HBM and FVM for nonlinear material properties. Material properties

Nonlinear

Linear

Location

FVM result

HBM Result

Relative error

0.06 0.10 0.14 0.18 0.22

926.84 698.33 520.23 363.57 214.22

926.84 698.32 520.24 363.56 214.23

0.00% 0.14% 0.19% 0.27% 0.00%

0.06 0.10 0.14 0.18 0.22

898.04 612.38 424.23 283.70 171.48

898.04 612.37 424.22 283.70 171.48

0.00% 0.16% 0.24% 0.00% 0.00%

Fig. 6. Nonlinear transient-state temperature distributions obtained by HBM with 80 elements and by FVM with 2000 elements.

In addition, the temperature distribution along the radial direction changes as the thermal conductivity coefficient is changed from linearity to nonlinearity. The maximum difference is approximately 96.4 K, and appears at the location of r = 0.135 m. Such big difference influences the accuracy of computational results a lot, and hence cannot be ignored in the practical engineering problems. 3.3. Transient-state nonlinearity with heat source

To solve the nonlinearity, the physical properties at a trial temperature (Tj) are first evaluated, and then these physical property values are used to calculate a new temperature Tj+1, and the physical properties at T1 are updated correspondingly. The relative error of temperature is calculated using the equation below:

Error ¼jT jþ1  T j j=T j

497

ð20Þ

In this paper, the precision limit is set to be 1.0e-8. If the error is smaller than the precision limit, the iteration is converged, otherwise the physical properties are updated at each iteration with the latest solution until the convergence is achieved. The flow chart for solving the one-dimensional temperature distributions with nonlinear material properties using HBM is illustrated in Fig. 4. Since the analytical solution to the nonlinear problem is difficult to obtain, the numerical results of FVM are treated as a reference to examine the accuracy of HBM. The accuracy of FVM increases with the element numbers, so 2000 elements are built in the numerical model to ensure the reasonable accuracy of FVM results, while the results of HBM are obtained with 20 elements. The calculated results of HBM and FVM with linear (average) and nonlinear thermal conductivity coefficients are presented in Fig. 5. The data of partial results are listed in Table 3. The results obtained by HBM with 20 elements match those by FVM with 2000 elements perfectly for both linear and nonlinear conditions. Table 3 shows that the relative errors of HBM results to FVM results are less than 0.27%. Considering the very few elements involved in the HBM calculation, the high accuracy of HBM can be expected.

In reality, the temperature may change with time. Therefore, in this section, HBM is utilized to solve transient-state problem with nonlinear material properties. In order to couple the practical conditions and increase the complexity of the problem, a constant heat source is considered to trigger the transient state. The 2e7 W/m3 heat source is located between r = 0.05 m and r = 0.15 m. The previous geometric model in Fig. 2 and the nonlinear thermal conductivity coefficient shown in Eq. (19) as well as cp = 460 J/(kgK) and q = 7860 kg/m3 are applied. The initial temperature distribution along the radius from the inner surface to the outer surface is constant 100 K. Assuming that at time t > 0 s, the heat source begins to work. The numerical results of FVM with 2000 elements and of HBM with 80 elements at different times are ploted in Fig. 6. The time step is 0.05 s. Fig. 6 shows the radial temperature distributions obtained by HBM and FVM with nonlinear material properties at every 60 s. The HBM results still match FVM results very well. Although only 80 elements are involved in HBM while 2000 elements are used in FVM, the maximum relative error of HBM results to those of FVM is below 0.2%, which confirms the high accuracy of HBM in solving the transient-state problem with nonlinear material properties. The highest temperature increases steadily for every 60 s in the area where the heat source is located, but then decreases sharply in the region without heat source. 4. HBM on nuclear fuel rod model The NFRs, which serve as the energy source of nuclear power plants, are one of the most important parts. Heat conduction in NFRs dissipates heat and prevents the pellet-cladding thermalmechanical interaction (Jiang et al., 2011; Imani et al., 2015) and core meltdown (An et al., 2014; U.S. NRC, 1979) safety issues caused by the high temperature. Therefore, it is significant to obtain the temperature distribution in NFRs. Actually, many researchers have developed some new methods to calculate the

498

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500

The heat transfer in NFRs is simplified as one-dimensional radial conduction problems. The geometric model for the NFR is illustrated in Fig. 7. The simplified model consists of a cylindrical UO2 fuel pellet and a thin layer of cladding, between which is the helium gap. The radius of the pellet is 4.3 mm, and the thicknesses of the helium gap and the cladding are 0.03 mm and 0.5 mm, respectively. The temperature of the cooling water is kept constant as Tf = 573 K. In the normal working condition, the volumetric heat generation rate is set as a constant, 3.064e8 W/m3. The central axis of the pellet endures the highest temperature, about 1000 K. The temperature changes along the NFR radius. Given the big temperature gradient experienced by the NFRs, the change in material properties caused by the temperature change plays a significant role. Therefore, in this part the effect of material nonlinearity on heat conduction in NFRs is studied. The property values of different materials in NFRs are obtained from the published data, and are summarized in Table 4. Fig. 7. The geometric sketch for the nuclear fuel rod model.

4.1. Steady-state

temperature distribution (Eskandari et al., 2012; Yu et al., 2012; Geelhood et al., 2014; Roostaii et al., 2016; Regis et al., 2000). In this section, HBM is used to solve the heat conduction problems of NFRs located in a pressurized water reactor. The transientstate heat conduction problems with linear and nonlinear material properties are carried out. The results can help optimize the design of NFRs and ensure the safe operation of nuclear power plants.

In this section, the radial temperature distributions in NFRs with linear and nonlinear material properties in the normal working condition are simulated using HBM. The numerical results obtained by ANSYS are also provided to verify the accuracy of HBM. Fig. 8 shows the cross-section view of the mesh model and the radial temperature distribution in NFRs with nonlinear material properties in the normal working condition obtained by ANSYS.

Table 4 A summary of material properties corresponding to the nuclear fuel rod. Material

Property

UO2

q

k

Dependence on temperature T (K) 2

10; 960  ða þ bT þ cT þ dT Þ   106 T > 923 K a ¼ 0:997; b ¼ 9:082 c ¼ 2:705  1010 ; d ¼ 4:391  1013   105 T > 923K a ¼ 0:997; b ¼ 1:179 c ¼ 2:429  109 ; d ¼ 1:219  1012     1 4:715109 exp  16;361 4 þ T T2 0:0375þ2:16510

He

264:256 þ 0:047T

q

0:0818  8:275  105 ðT  600Þ 4

T  6:79  10

8 2

Cp

0:0468 þ 3:81  10 5190.0

q

6490.0

k

7:51 þ 2:09  102 T  1:45  105 T 2 þ 7:67  109 T 3 255:66 þ 0:1024T ð273K < T < 1100KÞ

Cp

Units

Ref.

kg/m3

Ramirez et al. (2006)

W/(mK)

Lucuta et al. (1996)

J/(kgK)

Newman (2009)

kg/m3

Ramirez et al. (2006)

T

Cp k

Zircaloy

3 3

T

W/(mK) J/(kgK) kg/m3 W/(mK)

Hagrman (1995)

J/(kgK)

IAEA (2006)

Fig. 8. The mesh model of ANSYS and the temperature distribution in nuclear fuel rods obtained by ANSYS with nonlinear material properties in the normal working condition.

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500

Fig. 9. The radial temperature distributions in nuclear fuel rods simulated by ANSYS and HBM with linear and nonlinear material properties in the normal working condition.

Fig. 10. The highest temperatures in nuclear fuel rods associated with the dimension changes of the pellet, the helium layer and the zirconium cladding.

The mesh sensitivity is verified, and the number of elements is large enough to get an accurate result. Fig. 9 presents the radial temperature distributions in NFRs obtained by HBM and ANSYS with linear and nonlinear material properties in the normal working condition. The results of HBM coincide with those of ANSYS, indicating the high accuracy of HBM. On the other hand, the temperature difference associated with material properties is significant. In the center of NFR, the nonlinear material property leads to a 70 K higher temperature than linear material property. Hence, the nonlinear material properties really cannot be ignored in this problem because the material properties change a lot due to the sharp temperature gradient. As can be seen in Fig. 9, the central part of NFR is hottest. The temperature decreases along the radius to the outer surface of the cladding. Due to the small coefficient of thermal conductivity of the helium gap, the temperature decreases significantly in this part. In order to avoid melting of NFRs, the dimension of each part is changed. The influence of the geometric dimension change on the highest temperature in NFRs is studied. Each time, only one of the three radii of pellet, helium gap and cladding is increased from 50% to 150% of the original dimension. The highest temperature is calculated for every 10% increment. The corresponding

499

Fig. 11. The time-dependent temperature distributions at different locations along the radius of the nuclear fuel rod with 130% and 70% heat powers of the steadystate value at the first and last 60 s. P is the pellet, and I indicates the interfaces between the pellet and the helium gap, and between the helium gap and the cladding.

Fig. 12. The radial time-dependent temperature distributions in the nuclear fuel rod at different moments with 130% and 70% heat powers of the steady-state value at the first and last 60 s.

highest temperatures in NFRs associated with the dimension changes are recorded in Fig. 10. As Fig. 10 shows, the highest temperature in the NFR increases with the increase of pellet radius and the helium gap and cladding thicknesses. The temperature increase resulting from the pellet elongation is most conspicuous, followed by the helium gap and the cladding. The thick helium gap and cladding impede the heat flux from being transferred to the coolant quickly, while the elongated pellet generates more heat energy.

4.2. Transient-state The nuclear power plant is stable in most of time, but the unstable situation should be considered in case of safety issues. In this section, the HBM is applied to solve the one-dimensional transient-state radial conduction of NFRs with nonlinear material properties. The model is the same one as depicted in Fig. 4. The steady-state temperature obtained in the previous section (Fig. 9) is used as the initial temperature. The heat power of the pellet is assumed to suddenly increase by 30% (g = 3.9832e8 W/m3) and is kept for 60 s.

500

J. Tang et al. / Annals of Nuclear Energy 110 (2017) 492–500

After that, the heat power decrease to 70% (g = 2.1448e8 W/m3) of the original power value. The temperature changes at different locations along the radius of NFRs are displayed in Fig. 9. The temperatures in NFRs respond quickly to the change of heat power. When the power is suddenly increased by 30%, the temperature first increases. After 18 s, the temperature becomes stable. The highest temperature rises to 1119.40 K, and maintains for about 42 s. At 60th second, the temperature decreases quickly as the power reduces to 70% of the original power. After 15.7 s, the temperature becomes stable again. The highest temperature changes to 835.22 K. As can be seen in Fig. 11, the temperature changes more drastically as the distance is closer to the NFR center. Fig. 12 presents the radial temperature distributions in the NFR at different moments. In the first 60 s, the heat power is maintained as 130% of the steady-state value while in the last 60 s, the heat power turns to 70%. At each moment, the temperature decreases smoothly in the pellet and the cladding, while a big temperature gradient is produced in the helium gap. This is mainly caused by the small value of qcp of the helium although the low k of the helium leads to a poor heat conduction. In addition, in the first 60 s when the heat power is 130% of the normal state power, the overall temperature increases with a declined growth rate as the time lasts. However, the temperature growth rate declines. In the last 60 s, the decrease in the temperature presents a similar trend.

5. Conclusion In this paper, a new half boundary method (HBM) was proposed to solve the heat conduction problems of NFRs. The equations were presented with detailed derivation. Before the simulation of heat conduction of NFRs, three typical one-dimensional radial conduction problems were solved in cylindrical coordinates to test the accuracy and efficiency of HBM. The results of HBM were compared with those obtained by FVM and the analytical solution. HBM shows a comparable result to the analytical solution, and a much higher accuracy than FVM even with small element numbers. The computation time ratio of HBM to FVM decreases with the increase of element numbers, and it is only 0.16% when 1000 elements are involved in the calculation. In the simulation of the temperature conduction in NFRs, HBM was coupled with iterative scheme to solve one-dimensional radial conduction problems in NFRs with linear and nonlinear material properties. The calculation results were carried out in both steady-state and transient-state cases. The results show the significant influence of nonlinearity on the temperature distribution. The ANSYS comparison results also confirm the high accuracy of HBM. In addition, the size effect of pellet, helium gap and cladding on the highest temperature in NFRs was also studied. The pellet dimension affects the highest temperature mostly, followed by the helium gap and the cladding. NFRs can fast respond to the power change.

Acknowledgments The authors would like to appreciate the financial support from the State Key Laboratory of Intense Pulsed Radiation Simulation and Effect (No. SKLIPR1402X) of China.

References An, C., Moreira, F.C., Su, J., 2014. Thermal analysis of the melting process in a nuclear fuel rod. Appl. Therm. Eng. 68, 133–143. Chen, Z., He, C., Wu, B., 2008. High order finite volume methods for singular perturbation problems. Sci. China, Ser. A Math. 51, 1391–1400. Crank, J., Nicolson, P., 1947. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Math. Proc. Cambridge Philos. Soc. 43 (1), 50–67. Duda, P., 2015. Numerical and experimental verification of two methods for solving an inverse heat conduction problem. Int. J. Heat Mass Transfer 84, 1101–1112. Eskandari, Bavandi, Mihandoost, 2012. Studies on nuclear fuel rod thermal performance. Energy Proc. 14, 142–147. Geelhood, K.J., Luscher, W.G., Cuta, J.M., 2014. FRAPTRAN-1.5: A Computer Code for the Transient Analysis of Oxide Fuel Rods. U.S. Nuclear Regulatory Commission. Hagrman, D.T., 1995. SCDAP/RELAP5/MOD 3.1 code manual: MATPRO, A library of materials properties for Light-Water-Reactor accident analysis. Office of Scientific & Technical Information Technical Reports, Volume 4. He, Y., Li, J., 2010. A penalty finite element method based on the Euler implicit/explicit scheme for the time-dependent Navier-Stokes equations. J. Comput. Appl. Math. 235, 708–725. Huang, M., Sakiyama, T., Matsuda, H., Morita, C., 2015. Free vibration analysis of stepped rectangular plates resting on non-homogeneous elastic foundations. Eng. Anal. Boundary 50, 180–187. Huang, M., Tang, J., Zhao, Y., Ouyang, X., 2017. A new efficient and accurate procedure for solving heat conduction problems. Int. J. Heat Mass Transfer 111, 508–519. IAEA, 2006. Thermophysical Properties Database of Materials for Light Water Reactors and Heavy Water Reactors. International Atomic Energy Agency. Imani, M., Aghaie, M., Zolfaghari, A., Minuchehr, A., 2015. Numerical study of fuelclad mechanical interaction during long-term burnup of wwer1000. Ann. Nucl. Energy 80, 267–278. Jiang, Y., Cui, Y., Huo, Y., Ding, S., 2011. Three-dimensional FE analysis of the thermal–mechanical behaviors in the nuclear fuel rods. Ann. Nucl. Energy 38, 2581–2593. Krishnaprakas, C.K., Narayana, K.B., 2001. Use of GCG methods for the efficient solution of matrix problems arising from the FVM formulation of radiative transfer. Numer. Heat Transfer 40, 515–533. Li, R., 2005. Fundamental of Finite Volume Method. National Defense Industry Press. Loh, J.S., Azid, I.A., Seetharamu, K.N., Quadir, G.A., 2007. Fast transient thermal analysis of Fourier and non-Fourier heat conduction. Int. J. Heat Mass Transfer 50, 4400–4408. Lucuta, P.G., Matzke, H., Hastings, I.J., 1996. A pragmatic approach to modelling thermal conductivity of irradiated UO2 fuel: review and recommendations. J. Nucl. Mater. 232, 166–180. Lyra, P.R.M., Lima, R.d.C.F.d., Carvalho, D.K.E.d., Silva, G.M.L.L.d., 2005. An axisymmetric finite volume formulation for the solution of heat conduction problems using unstructured meshes. J. Braz. Soc. Mech. Sci. Eng. XXVII, 407– 414. Newman, C., 2009. Three dimensional coupled simulation of thermomechanics, heat, and oxygen diffusion in UO2 nuclear fuel rods. J. Nucl. Mater. 392, 6–15. Ramirez, J.C., Stan, M., Cristea, P., 2006. Simulations of heat and oxygen diffusion in UO2 nuclear fuel rods. J. Nucl. Mater. 359, 174–184. Regis, C.R., Cotta, R.M., Su, J., 2000. Improved lumped analysis of transient heat conduction in a nuclear fuel rod. Int. Commun. Heat Mass Transfer 27, 357–366. Roostaii, B., Kazeminejad, H., Khakshournia, S., 2016. Influence of porosity formation on irradiated UO2 fuel thermal conductivity at high burnup. J. Nucl. Mater. 479, 374–381. Sakeyama, T., Huang, M., 2000. Free-vibration analysis of right triangular plates with variable thickness. J. Sound Vib. 234, 841–858. Todd, R.S., He, J., Webley, P.A., Beh, C., Wilson, S., Lloyd, M.A., 2001. Fast finitevolume method for PSA/VSA cycle simulations experimental validation. Ind. Eng. Chem. Res. 40, 3217–3224. Tomatis, D., 2013. Heat conduction in nuclear fuel by the Kirchhoff transformation. Ann. Nucl. Energy 5, 100–105. U.S. NRC, 1979. Backgrounder on the Three Mile Island Accident. United State Nuclear Regulatory Commission, 1–7. Versteeg, H., Malalasekera, W., 1995. An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Pearson Education Limited. Wang, C.-C., Chen, C.o.-K., Chiang, C.-Y., 2009. Analysis on error bounds of heat transfer problems using the fast residual correction method. Numer. Heat Transfer 56, 455–477. Yu, H., Tian, W., Yang, Z., Su, G.H., Qiu, S., 2012. Development of Fuel ROd Behavior Analysis code (FROBA) and its application to AP1000. Ann. Nucl. Energy 50, 8– 17. Zhang, Y., Xu, Y., Liu, S., 2014. A fast numerical method for the analysis of the heat transfer performance of graded metallic honeycomb materials. Int. J. Heat Mass Transfer 79, 507–517.