A nodal SN transport method for three-dimensional triangular-z geometry

Nuclear Engineering and Design 237 (2007) 830–839

A nodal SN transport method for three-dimensional triangular-z geometry Haoliang Lu ∗ , Hongchun Wu Department of Nuclear Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China Received 30 June 2006; received in revised form 23 October 2006; accepted 27 October 2006

Abstract A new nodal SN transport method has been developed to perform accurate transport calculation in three-dimensional triangular-z geometry, where arbitrary triangles are transformed into regular triangles via a coordinate transformation. The transverse integration procedure is applied to treat the neutron transport equation in the regular triangle. The neutron angular distributions of intra-node fluxes are represented using the SN quadrature set, and the spatial distributions of neutron fluxes and sources are approximated by a quadratic polynomial. The nodal-equivalent finite difference algorithm for 3D triangular geometry is applied to establish a stable and efficient iterative scheme. The present method was tested on four 3D Takeda benchmark problems published by the nuclear data agency (NEACRP), in which the first three problems are in XYZ geometry and the last one is in hexagonal-z geometry. The results of the present method agree well with those of the reference Monte-Carlo calculation method, the difference in keff being less than 0.1%. This shows that multi-group reactor core/criticality problems can be accurately and effectively solved using the present method. © 2006 Elsevier B.V. All rights reserved.

1. Introduction The nodal SN methods are widely accepted as accurate and efficient tools for the analysis of light water reactors (Badruzzaman, 1985) in Cartesian geometry and of fast reactors in hexagonal-z geometry (Wagner, 1989; Ikeda and Takeda, 1994). The regions of these methods treated with are structured meshes, such as rectangles and hexagons. For the problems with rapid variable fluxes in certain regions, like control rod positions, the accuracy of the nodal SN method is unsatisfactory (Takeda and Ikeda, 1991a). One of the ways is to use a method in triangular-z geometry so that the meshes of these regions can be refined. The finite methods, such as a finite element spherical harmonics method (McGhee et al., 1997; Cao and Wu, 2004) and a discrete ordinates finite element method (Ressel and Starke, 2000), can solve the problems with triangular meshes. However, the computational cost of these methods is great. So there is a need for developing an accurate and efficient nodal SN method for the solutions of three-dimensional reactor problems in triangular-z geometry. This paper describes a three-dimensional nodal SN method for triangular-z geometry. Unlike rectangular or hexagonal meshes ∗

Corresponding author. Tel.: +86 29 8266 3285; fax: +86 29 8266 7802. E-mail address: [email protected] (H. Lu).

0029-5493/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2006.10.025

with the same shapes and sizes, those of triangular meshes are different from each other. For simplicity, the coordinate transformation of the triangle (Gallagher, 1975) is necessary and described in Section 2.1. The transverse integration procedure is applied to the nodal SN theory scheme for 3D triangular-z geometry and this part is described in Section 2.2. The angular distributions of fluxes are represented using the SN quadrature set, and the spatial distributions of the fluxes and sources are approximated by an orthogonal quadratic polynomial. The radial leakage of the radial coupling equation is approximated by a second-order polynomial. To achieve a stable convergence, we apply the nodal-equivalent finite difference (NEFD) algorithm (Badruzzaman, 1985; Ikeda and Takeda, 1994) for 3D triangularz geometry and this part is described in Section 2.3. In order to investigate the applicability of this new method, we calculated all four 3D multi-group eigenvalue benchmark problems proposed from Osaka University to Nuclear Energy Agency committee on Reactor Physics Benchmarks (NEACRP) in 1991 (Takeda and Ikeda, 1991a,b) and this part is described in Section 3. Conclusions are summarized in Section 4. 2. Theory The nodal SN method for 3D triangular-z geometry described is similar to the discrete nodal transport method for Cartesian

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

831

The neutron transport equation with isotropic scattering is written in standard notation as μm

∂ψm (x, y, z) ∂ψm (x, y, z) ξ m ∂ψm (x, y, z) + ηm + ∂x ∂y hz ∂z

+Σt ψm (x, y, z) = Q(x, y, z)

Fig. 1. Regular triangular mesh.

geometry (DNTM) (Badruzzaman, 1985) and the nodal SN method for hexagonal geometry (NSHEX) (Ikeda and Takeda, 1994). In all these methods, the transverse integration procedure is used for deriving several 1D transport equations in the radial plane and one 1D equation in the axial plane. These equations relate the surface averaged fluxes to nodal averaged moments. To achieve a convergent iterative scheme, the NEFD algorithm is applied. 2.1. Coordinate transformation For a computational region, we can use triangular meshes to approximate it. However, unlike the rectangular or hexagonal meshes, these triangular meshes are not the same in the shapes and sizes with each other and we cannot apply the transverse integration procedure simply. In the present method, the area coordinate, which is usually used in the finite element methods is acceptable to transform an arbitrary triangle into a regular triangle. Suppose P is a point inside an arbitrary triangle, and this triangle could be divided into three triangles by connecting P with three vertexes. The areas of these four triangles are Δ, Δ1 , Δ2 , Δ3 , respectively. Then the area coordinates of the point P can be defined below λ1 =

Δ1 , Δ

λ2 =

Δ2 , Δ

λ3 =

Δ3 Δ

Thus the relationship between the orient physical coordinate and the area coordinate is set up. Using this relationship twice, we can transform an arbitrary triangle in physical coordinate (x, y) into an isoceles triangle in area coordinate (λ1 , λ2 , λ3 ), and then transform this new triangle into a regular triangle in computational coordinate (x , y ) (see Fig. 1). The final expression can be written as √   3(−x2 + x3 )  x2 + x3 x1 + x2 + x3  + −x1 + x+ y, x= 3 2 2   y1 + y2 + y3 y2 + y3 y= + −y1 + x 3 2 √ 3(−y2 + y3 )  + (1) y 2 where x1 , x2 , x3 , y1 , y2 , y3 are the coordinates of three vertexes.

(2)

Here, hz is the nodal axial width and −1/2 ≤ z ≤ 1/2. The values μm , ηm and ξ m are direction cosines relative to x, y and z directions. Using Eq. (1), one can transform Eq. (2) into a neutron transport equation in a regular triangle (see Fig. 1) as μm x

m   ξ m ∂ψm (x , y , z) ∂ψm (x , y , z) m ∂ψ (x , y , z) + η + x ∂x ∂y hz ∂z

+Σt ψm (x , y , z) = Q(x , y , z)

(3)

where (−y2 + y3 )μm + (x2 − x3 )ηm , 2Δ (−x1 + 1/2x2 + 1/2x3 )ηm + (y1 − 1/2y2 − 1/2y3 )μm √ ηm , x = 3Δ 2 x + 2/3 1 − ≤ x ≤ , −ys (x) ≤ y ≤ ys (x), ys (x) = √ 3 3 3 μm x =

After the coordinate transformation, the nodal averaged flux and surface averaged fluxes in the computational coordinate, which will be defined in Section 2.2, may be proved equal to those in the physical coordinate through a simple integral transformation. 2.2. Derivation of nodal coupling equation Omitted the superscripts of the computational coordinate (x , the Eq. (3) will be

y ), μm x

∂ψm (x, y, z) ∂ψm (x, y, z) ξ m ∂ψm (x, y, z) + ηm + x ∂x ∂y hz ∂z

+Σt ψm (x, y, z) = Q(x, y, z)

(4)

Integrating Eq. (4) over y and z, one obtains the onedimensional equation in x direction ¯ xm (x)} d{ys (x)ψ ¯ xm (x) + Σt ys (x)ψ dx m ¯ x (x) − Lzm = ys (x)(Q x (x)) − Lrx (x)

μm x

(5)

¯ xm (x) and Q ¯ x (x) represent the 1D averaged angular where ψ flux and averaged neutron source within the node, respectively. Lrxm (x) is a radial leakage, which is given by m m m μm u ψu (x) + μv ψv (x) √ (6) 3 √ m √ m m m m where μm u = −(μx − 3ηx )/2, μv = −(μx + 3ηx )/2,  1/2 and ψvm (x) = ψum (x) = −1/2 ψm (x, y, z)|y=ys (x) dz

Lrxm (x) =

832

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

 1/2

m

Lzx (x) is an axial leakage in the radial plant and written as

spatial distribution of the nodal flux is represented by a binary quadratic polynomial in the computational coordinate

ξm ¯ m m ¯ z− = (ψ (x) − ψ (x)) (7) hz z+  m (x) = ys ψ m (x, y, ±1/2) dy/(2y ). ¯ z± where ψ s −ys Solving Eq. (5) for μm x > 0, one obtains  x  m  ¯ xm (x) = 1 ¯ x (x ) − Lzm ys (x)ψ [ys (x )(Q x (x )) − Lrx (x )] μm x −2/3

ψm (x, y) = a(x2 + y2 ) + bx + cy + d

−1/2 ψ

m (x, y, z)| y=−ys (x) dz.

m Lzx (x)



× e−(Σt /μx )(x−x ) dx m

(8)

The outgoing surface averaged flux at the right (x = 1/3) from Eq. (8) is √  1/3 m m ¯x = 3 ¯ x (x) − Lzm ψ [ys (x)(Q x (x)) − Lrx (x)] μm −2/3 x × e−(Σt /μx )(1/3−x) dx m

(9)

Because there is only one surface in x direction, the equation for μm x < 0 differs from Eq. (8) and is written below  1/3 1 m  m  ¯ ¯ x (x ) − Lzm ys (x)ψx (x) = m [ys (x )(Q x (x )) − Lrx (x )] |μx | x 

× e−(Σt /|μx |)(x −x) dx √ 3 −(Σt /|μmx |)(1/3−x) ¯ m + ψx e 3 m

Similar expressions for nodal averaged surface fluxes are obtained for the u and v directions. Eqs. (8)–(10) are the analytical solutions in the x direction in the transverse integrated nodal method. In the DNTM, the nodal averaged flux is expanded in terms of Legendre polynomials (Badruzzaman, 1985). In the present method, the nodal averaged flux and source are expanded in terms of a quadratic polynomial as ¯ xm (x) = ψ

2 

m ¯ xi ψ hi (x)

(11)

¯ xi hi (x) Q

(12)

i=0

¯ x (x) = Q

2 

The coefficients a, b, c, d are solved using the nodal averaged flux and three surface averaged fluxes. The fluxes of three vertexes in a regular triangle is obtained ¯ vm ) ψ ¯ um + ψ ¯m 5(ψ ¯ m, − x − 2ψ 3 3 ¯ vm + ψ ¯m ¯ xm ) ψ 5(ψ m ¯ m, ψpu = − u − 2ψ 3 3 ¯ um ) ψ ¯ xm + ψ ¯m 5(ψ m ¯m ψpv = − v − 2ψ 3 3

m ψpx =

(14)

Using the surface averaged flux and two fluxes of its adjacent vertexes, one can approximate the spatial distribution of this surface flux by a quadratic polynomial. Then the spatial distribution of the radial leakage is expressed by m m m 2 Lrxm (x) = Lrx0 + Lrx1 x + Lrx2 x

(15)

where 1 1 m m ¯ um − ψpx ¯ vm − ψpx ) + μv (4ψ ), μu (4ψ 3 3 m m ¯ um ) + 2μv (ψpu ¯ vm ), = 2μu (ψpv −ψ −ψ

m Lrx0 = m Lrx1

(10)

(13)

m m m m m ¯ um ) + 3μv (ψpx ¯ vm ) Lrx2 = 3μu (ψpx + ψpv − 2ψ + ψpu − 2ψ

The nodal coupling equations in the axial direction are similar to those of the nodal transport method in Cartesian geometry. Integrating the Eq. (4) over x and y, one obtains the onedimensional equation in the z direction ¯ zm (z) ξ m dψ ¯ zm (z) = Q ¯ z (z) − L ¯m + Σt ψ xy (z) hz dz

(16)

m¯m m¯m m¯m ¯m where L xy (z) = 2[μx ψx (z) + μu ψu (z) + μv ψv (z)]. m And solving Eq. (16) for ξ > 0, one obtains  z  −(Σt /ξ m )(z−z ) ¯ z (z ) − L ¯ zm (z) = hz ¯m [Q dz ψ xy (z )] e m ξ −1/2

+e−(Σt /ξ

m )(z+1/2)

m ¯ z− ψ

(17)

i=0

where hi (x)={1, x, x2 + 2x/15 − 1/18}. The polynomials hi (x) form an orthogonal set with respect to the weight function ys (x)  1/3 ys (x)hi (x)hj (x) dx = 0, for i = j −2/3

The spatial distribution of the axial leakage in Eq. (7) is expanded by a flat approximation and a reasonable result can be achieved. However, the radial leakage in a triangle mesh cannot use a flat approximation due to the relatively large leakage contribution. Here, we develop a new spatial model for radial leakage in order to describe the spatial distribution more accurately. The

The outgoing surface averaged flux at the upper (z = 1/2) from Eq. (17) is  hz 1/2 ¯ m −(Σt /ξ m )(1/2−z) ¯ ¯m ψz+ = m [Qz (z) − L dz xy (z)] e ξ −1/2 m m ¯ z− +e−(Σt /ξ ) ψ

(18)

The spatial distribution of the flux and source is approximated by the quadratic polynomial ¯ zm (z) = ψ

2  i=0

m ¯ zi fi (z) ψ

(19)

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

¯ z (z) = Q

2 

¯ zi fi (z) Q

(20)

i=0

where the polynomial function series fi (z)={1, z, z2 − 1/12} are orthogonal over −1/2 ≤ z ≤ 1/2. To obtain a simple and conver¯m gent scheme, the transverse-leakage L xy (z) is approximated by a flat approximation. Substituting Eqs. (11), (12), (15) into Eqs. (8) and (9) and Eqs. (19) and (20) into Eqs. (17) and (18) with μm x > 0 and ξ m > 0, one obtains m ¯ x0 ψ

=

m ¯ Fx00 (Qx0

m − Lzx ) +

2 

m ¯ Fx0j Qxj

j=1

m m ¯ m ¯ xi ψ = Fxi0 (Qx0 − Lzx ) +

2 

+

2 

+

arrangements are complicated such as the KNK-II experiment fast reactor. To achieve a stable convergence and more efficient scheme of the nodal SN method, we apply the NEFD algorithm in 3D triangular-z geometry. To arrive at this algorithm, we reformulate Eqs. (21) through (24) and their analogues in u and v directions. Consider the x ¯ x0 − Lzm direction only. Eliminating Q x from Eqs. (21a), (21b) and (22), one obtains  2  m m  Fxi0 Fxi0 m m m m ¯ ¯ xj ¯ ψxi = m ψ + Fxij − m Fx0j Q Fx00 Fx00 j=1

m Gm x0j Lrxj

j=0

+ (21a)

j=0

(21b)

+

m¯ Pxi Qxi

i=1 m m m ¯ ¯ z0 ψ = Fz00 (Qz0 − Lxy ) +

+

2 

m Rm xi Lrxi

(22)

i=0

2 

2  

Rm xi

i=0

j=0 2 

Gm xij −

m ¯ m ¯m Fz0j Qzj + Hz0 ψz−

(23a)

 + Hzim −

m m m ¯ ¯ zi ψ = Fzi0 (Qz0 − Lxy ) +

m ¯ m ¯ z− Fzij , for i = 1, 2 Qzj + Hzim ψ

j=1

m m ¯ m ¯ z+ ψ = Pz0 (Qz0 − Lxy ) +

2 

(25)

(23b)

m ¯ zi + T m ψ ¯ z− Pzim Q

(24)

i=1

where the coefficients F, G, H, P, R, T are relative to Σ t , μm , ηm , ξ m and the properties of the triangular-z mesh. Eqs. (21)–(24) form the basic equations of the nodal SN method. 2.3. NEFD algorithm for 3D triangular-z geometry In the original DNTM, the standard algorithm of the nodal SN transport method is used. The outgoing surface averaged fluxes are solved using an initial guess of the incoming surface averaged fluxes and the angular flux moments are calculated to update the source moments. The nodal averaged flux is given by averaging three 0th-order angular flux moments in x, y and z directions. To reduce the computer memory and execution time, Badruzzaman (1985) developed the NEFD (nodal-equivalent finite difference) algorithm in which 0th-order angular source moments are eliminated and the nodal averaged flux is calculated using the neutron balance equation in rectangular-z mesh. And Ikeda and Takeda (1994) find that convergence is difficult to achieve using the standard nodal scheme in calculations for reactors whose core

 m Px0 m m − m Gx0i Lrxi Fx00

(26)

¯ z0 − L ¯m Eliminating Q xy from Eqs. (23a), (23b) and (24), one obtains  2  m Fm m  Fzi0 m m m ¯ ¯ zi ¯ + F ψ ψ = zi0 − F zij z0j Qzj m m Fz00 Fz00

j=1 2 

 m Fxi0 m m G x0j Lrxj , for i = 1, 2 m Fx00

i=1

m Gm xij Lrxj , for i = 1, 2

m m ¯ ¯ xm = Px0 ψ (Qx0 − Lzx ) +

2  

 2  m m  Px0 m m ¯ ¯ xm = Px0 ψ ¯m + P ψ − F xi x0i Qxi m m Fx00 Fx00

m ¯ Fxij Qxj

j=1 2 

833

m ¯ z+ = ψ

j=1

 m Fzi0 m ¯m H z0 ψz− , for i = 1, 2 m Fz00

 2  m m  Pz0 Pz0 m m ¯ P + − F zi z0i Qzi m m Fz00 Fz00 i=1   m Pz0 m m ¯ z− ψ Hz0 + Tm − m Fz00

(27)

(28)

The nodal balance equation is given by integrating Eq. (4) over −2/3 ≤ x ≤ 1/3 m¯m m¯m ¯m 2μm x ψx + 2μu ψu + 2μv ψv +

¯ ¯m = Q +Σt ψ

ξm ¯ m m ¯ z− (ψ − ψ ) hz z+ (29)

Considering a discrete direction m, seen in Fig. 2, we could find that there are two incoming surfaces for triangle 1 while only one incoming surface for triangle 2 in the radial plant. We need to treat these two cases separately. If there are two incoming surfaces in x and u directions, the surface in v direction is a outgoing surface and its surface averaged flux can be expressed in terms of the node-averaged angular flux, higher order averaged source moments and two incoming surface averaged fluxes in radial plant through the insertion of Eq. (14) into Eq. (26). Addition of Eqs. (28) and (29) to this equation for the outgoing surface, yields the set of ternary linear equations, and the averaged angular fluxes at the two outgoing surfaces and the node-averaged angular flux could be obtained.

834

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

Fig. 2. Distribution of triangular meshes.

In case of one incoming surface, such as x direction, the surfaces in u and v directions are outgoing surfaces and those surface averaged fluxes can be expressed in terms of the nodeaveraged angular flux, higher order averaged source moments and an incoming surface averaged fluxes in x direction through the insertion of Eq. (14) into Eq. (26) respectively. Once again, the addition of Eqs. (28) and (29) to the equations for the outgoing surfaces yield the set of quaternary linear equations and the averaged angular fluxes at the three outgoing surfaces and the node-averaged angular flux could be obtained. Using Eqs. (25) and (27), all the high-order moments are calculated and used to form the high-order source moments in every directions. The derivation and calculation of higher order m flux moments for μm x < 0 are similar to those for μx > 0. 3. Numerical results The nodal transport method described above has been implemented in the computer program DNTR3D. The accuracy and computational efficiency of DNTR3D are demonstrated by numerical results presented for four 3D neutron transport benchmark problems compiled by Takeda and Ikeda (1991a,b). These benchmark problems were set up to be used for checking the validity of 3D neutron transport codes. The benchmark problems consist of four core models, a small LWR core, a small FBR core, an axially heterogeneous FBR core (which are all in XYZ geometry), and a small FBR core with hexagonal-z geometry. The required cross-sections for all problems are taken from Takeda and Ikeda (1991a,b). The reference value is the average of Monte-Carlo calculations contributed to the benchmark, which was computed by the variance weighted procedure (Takeda and Ikeda, 1991b). Additional reference value is given by the 3D TDOT code which is developed from the 2D DOT 4.2 code. 3.1. Model 1 (small LWR core) This core consists of a fuel region, a void region (in place of control rod) and a reflector region. The core configuration is illustrated in Fig. 3. Calculations have been carried out in two cases:

Fig. 3. Core configuration of model 1.

• Case 1: The control rod position is empty (void). • Case 2: The control rod is inserted. S8 solutions are given by the DNTR3D for each case and the reference mesh numbers are 10 axial meshes and 383 radial triangular meshes in the 1/8 core model. The TDOT code is calculated with S8 discrete ordinates and the mesh size is 1 cm × 1 cm × 1 cm. Table 1 lists keff values for two control rod patterns (rod-out and rod-in) and the control rod reactivity worth. DNTR3D agrees with the Monte-Carlo method with the differences being less than 0.05%, which is almost the same as the standard deviation about the reference. Good agreement is also obtained between TDOT and DNTR3D by the S8 calculations. Now we discuss the spatial mesh effect and angular quadrature effect in the present nodal SN method. Several results for keff were calculated by independently varying the number axial meshes, radial meshes and order of the angular quadrature. These results are presented in Tables 2–4 respectively. According to these results, the spatial mesh effect in the axial plant is relatively larger than that of the radial plant. The reason for this is that the flat approximation of transverse leakage is used in the axial direction. However, the calculation with 10 axial layers is

Table 1 Comparison of keff and control rod worth for model 1 Method

Case 1

Case 2

CR-worth

Monte-Carlo TDOT DNTR3D

0.9780 ± 0.0006a

0.9624 ± 0.0006 0.9627, 0.03 0.9628, 0.04

1.66E−02 ± 0.09E−02 1.58E−02, −4.8 1.56E−02, −6.0

a b

0.9776, −0.04b 0.9775, −0.05

Relative standard deviation. Relative differences from reference (% k/k).

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839 Table 2 Axial mesh effect in keff for model 1 3 (5)a Rod-out Rod-in a b

0.21b 0.22

3 (7) 0.08 0.09

Table 6 Comparison of execution times for case 2 of model 1 6 (10)

15 (25)

Ref. (0.9775) Ref. (0.9628)

−0.04 −0.05

Numbers of fuel regional axial mesh (total axial mesh). Relative differences from reference (% k/k).

a b

1809 92

Relative differences from reference (%).

485

Method

Case 1

Case 2

CR-worth

0.11b 0.05

0.02 −0.01

Ref. (0.9775) Ref. (0.9628)

0.00 0.00

Monte-Carlo TDOT DNTR3D

0.9732 ± 0.0002a 0.9732, 0.00b 0.9736, 0.04

0.9594 ± 0.0002 0.9588, −0.06 0.9597, 0.03

1.47E−02 ± 0.03E−02 1.54E−02, 4.76 1.49E−02, 1.36

a b

S4

S6

S8

S10

−0.014a

−0.005 −0.007

Ref. (0.9775) Ref. (0.9628)

0.015 0.012

−0.022

Monte-Carlo

TDOT

4.9125E−03(0.10)a

−0.50b

8.6921E−04(0.13)

−0.15

−0.33 −0.18

Reflector 1G 2G

5.9109E−04 (0.21) 8.7897E−04 (0.23)

0.12 0.51

−0.23 −0.13

Control rod 1G 2G

1.2247E−03 (0.48) 2.4604E−04 (0.72)

0.43 0.07

0.36 −0.01

Relative standard deviation (%). Relative differences from reference (%).

Relative standard deviation. Relative differences from reference (% k/k).

3.2. Model 2 (small FBR core)

Table 5 Region averaged group fluxes and differences for rod-in case of model 1

a

0.9630 0.9627 (0.031) 0.9624

383

accurate enough. The angular quadrature effect is small and the differences are less than 0.03%. Region averaged group fluxes for case 2 are compared in Table 5. The region averaged fluxes are close to the MonteCarlo results and the difference is within 1% over all regions and all energy groups. To check the efficiency of the present nodal SN method, we compare the execution time of DNTR3D with that of TEPFEM for case 2. Here, the TEPFEM code is based on the finite element spherical-harmonics method and also applied to unstructured meshes (Cao and Wu, 2004). The results for the two codes, using comparative meshing schemes, are summarized in Table 6. Even though the keff values of the two methods are close to the reference (the difference is less than 0.1%), the execution time of DNTR3D is only one-twentieth of the TEPFEM’s.

b

TEPFEM (P5 ) DNTR3D (S6 ) Reference

232

Relative differences from reference (%).

Core 1G 2G

Execution time (s) (0.062)a

81a

Table 4 Angular quadrature effect in keff for model 1

a

Keff

Table 7 Comparison of keff and control rod worth for model 2

Numbers of radial mesh. Relative differences from reference (% k/k).

Rod-out Rod-in

Method

a

Table 3 Radial mesh effect in keff for model 1

Rod-out Rod-in

835

DNTR3D

This core is a small model of a FBR. The core configuration is illustrated in Fig. 4. Calculations have been carried out in two cases: • Case 1: The control rod is withdrawn (the control rod position is filled with Na). • Case 2: The control rod is half-inserted. For case 1, S8 solutions are given by the DNTR3D with 14 axial meshes and 420 radial triangular meshes and by the TDOT with the mesh size of 1 cm × 1 cm × 1 cm in the 1/8 core model. For case 2, S4 solutions are given by the DNTR3D with 32 Table 8 Region averaged group fluxes and differences for rod-in case of model 2 Monte-Carlo

TDOT

DNTR3D

Core 1G 2G 3G 4G

4.3482E−05 (0.06)a 2.4171E−04(0.05) 1.6200E−04(0.06) 6.0438E−06 (0.21)

−0.03b 0.03 −0.01 −0.35

−0.10 0.05 0.06 −0.24

Axial blanket 1G 2G 3G 4G

5.2209E−06 (0.27) 4.6772E−05(0.14) 4.6190E−05 (0.16) 3.6287E−06 (0.45)

1.36 1.03 0.55 −0.06

0.45 −0.21 −0.31 −0.20

Control rod position 1G 2.5902E−05 (0.46) 2G 1.6779E−04 (0.25) 3G 1.2551E−04 (0.31) 4G 7.0648E−06 (1.79)

0.93 0.07 −0.50 −5.81

0.24 −0.39 −0.73 −5.12

1.38 1.51 1.14 3.11

0.71 0.47 0.31 1.03

Control rod 1G 2G 3G 4G a b

1.6556E−05 (0.54) 9.1050E−05(0.26) 5.1815E−05 (0.30) 1.1073E−06(1.00)

Relative standard deviation (%). Relative differences from reference (%).

836

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

Fig. 4. Core configuration of model 2.

axial meshes and 420 radial triangular meshes and by the TDOT with the mesh size of 5 cm × 5 cm × 5 cm in the 1/4 core model. Table 7 lists keff values for the two cases, as well as the associated control rod worth. The eigenvalues of the DNTR3D agree well with those of the Monte-Carlo method with the differences being less than 0.05%. Once again good agreement is obtained between the TDOT and DNTR3D.

Table 8 shows the region averaged group fluxes for the rodin case. The difference in the fourth group in the Control rod position is −5.12%. Except for the above mentioned difference, the differences in the region averaged group fluxes are less than 1%. In the control rod, the differences of DNTR3D are less than those of TDOT. This is due to the triangular mesh refinement.

Fig. 5. Core configuration of model 3.

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

837

Table 9 Comparison of keff and control rod worth for model 3 Method

Case 1

Case 2

Case 3

CR-worth

CRP-worth

Monte-Carlo TDOT DNTR3D

0.9709 ± 0.0002a

1.0005 ± 0.0002 1.0003, −0.02 1.0015, 0.10

1.0214 ± 0.0002 1.0209, −0.05 1.0221, 0.07

3.05E−02 ± 0.03E−02 3.07E−02, 0.67 3.05E−02, 0.00

2.03E−02 ± 0.00E−02 2.02E−02, −0.49 2.01E−02, 1.00

a b

0.9705, −0.04b 0.9718, 0.09

Relative standard deviation. Relative differences from reference (% k/k).

Table 10 Region averaged group fluxes and differences for rod-in case of model 3 Monte-Carlo

TDOT

DNTR3D

Core 1G 2G 3G 4G

2.0334E−5 (0.06)a 1.1501E−4 (0.04) 8.2784E−5 (0.05) 3.4093E−6 (0.19)

0.19b −0.01 −0.11 −0.57

0.03 0.03 −0.07 −0.35

Internal blanket 1G 2G 3G 4G

1.2562E−5 (0.35) 1.1237E−4 (0.20) 1.1931E−4 (0.22) 7.4318E−6 (0.65)

0.97 0.46 0.51 0.11

0.39 −0.61 −0.42 −0.47

Radial blanket 1G 2G 3G 4G

1.5888E−6 (0.36) 1.6736E−5 (0.22) 2.1521E−5 (0.22) 1.7973E−6 (0.56)

1.61 0.74 0.73 0.98

1.38 0.16 0.04 0.35

Axial blanket 1G 2G 3G 4G

2.9405E−6 (0.23) 3.0115E−5 (0.12) 3.5408E−5 (0.13) 3.7516E−6 (0.32)

0.33 −0.97 −2.39 −4.37

0.68 −0.39 −0.47 −0.21

a b

Table 12 Region averaged group fluxes and differences for case 3 of model 4 GMVP

NSHEX

DNTR3D

Axial blanket 1G 2G 3G 4G

4.07E−5 (0.30)a 5.17E−5 (0.30) 2.36E−5 (0.42) 9.74E−6 (0.56)

−0.05b 0.27 −0.43 −0.02

−0.62 0.15 −0.57 0.30

Test zone 1G 2G 3G 4G

1.66E−4 (0.15) 1.17E−4 (0.17) 2.18E−5 (0.30) 1.18E−6 (0.62)

0.06 0.31 −0.15 −3.67

−0.30 −0.07 −0.14 −0.83

Driver without moderator 1G 1.11E−4 (0.13) 2G 7.79E−5 (0.15) 3G 2.12E−5 (0.25) 4G 3.91E−6 (0.31)

−0.56 −0.26 −0.32 −0.27

−0.23 0.07 −0.02 0.15

0.24 0.37 0.24 0.54

0.08 0.19 0.16 −0.54

Control rod 1G 2G 3G 4G a

Relative standard deviation (%). Relative differences from reference (%).

b

3.3. Model 3 (axially heterogeneous FBR core) This core has an internal blanket region and the core configuration is shown in Fig. 5. Calculations have been carried out in three cases: • Case 1: The control rods are inserted. • Case 2: The control rods are withdrawn. • Case 3: The control rods are replaced with core or blanket cells. S4 solutions are given by the DNTR3D and the mesh is chosen as with 11 axial meshes and 922 radial triangular meshes

9.53E−5 (0.14) 6.92E−5 (0.12) 1.34E−5 (0.21) 9.55E−7 (0.36)

Relative standard deviation (%). Relative differences from reference (%).

for all cases. With the mesh size of 1 cm × 1 cm × 1 cm, S4 and S8 solutions are given by the TDOT for case 1 and the others, respectively. Table 9 lists keff values for the three cases and the associated control rod worth. The eigenvalues of the DNTR3D agree well with those of the Monte-Carlo method with the differences being less than 0.10% and the differences between the TDOT and DNTR3D are of similar magnitude. The differences in region averaged group fluxes for the rodin case are shown in Table 10. The difference in the first group in the radial blanket is 1.38%. Except for this difference, the differences in the region averaged group fluxes are less than 1% for the present method. A −4.37% difference in the fourth group in the axial blanket is observed for the TDOT. This difference is

Table 11 Comparison of keff and control rod worth for model 4 Method

Case 1

Case 2

Case 3

CR-worth

GMVP NSHEX DNTR3D

1.0955 ± 0.0005a

0.9839 ± 0.0004 0.9842, 0.03 0.9844, 0.05

0.8802 ± 0.0004 0.8804, 0.02 0.8805, 0.03

22.3E−02 22.4E−02, 4.48 22.3E−02, 0.00

a b

1.0966, 0.10b 1.0960, 0.05

Relative standard deviation. Relative differences from reference (% k/k).

838

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

Fig. 6. Core configuration of model 4.

reduced to −0.21% for the present method. The same reduction is observed in the third group in the axial blanket. 3.4. Model 4 (small FBR core with hexagonal-z geometry) This core with the vacuum boundary conditions has a hexagonal-z geometry as shown in Fig. 6. Calculations have been carried out in three cases: • Case 1: The control rods are withdrawn. • Case 2: The control rods are half-inserted. • Case 3: The control rods are fully inserted. S4 solutions are given by the DNTR3D. For case 1 and case 2, the benchmark is calculated as a 1/4 core with 32 axial meshes and 404 radial triangular meshes. For case 3, the benchmark is calculated in the 1/8 core with 19 axial meshes and 404 radial triangular meshes. The reference value is given by the Monte-Carlo calculation code GMVP (Sugino and Takeda, 1996). Additional reference value is given by the NSHEX code (Sugino and Takeda, 1996) which is a nodal SN transport method

for 3D hexagonal geometry. NSHEX is calculated with S4 discrete ordinates and 10 axial meshes for case 3. Table 11 lists keff values for three cases, and the control rod worth. The eigenvalues and control rod worth of the DNTR3D agree well with those of the Monte-Carlo method with the differences being less than 0.05%. The differences between the NSHEX and DNTR3D are also small. The differences in region averaged group fluxes for case 3 are shown in Table 12. The differences in all region averaged group fluxes are less than 1% for the present method. A −3.67% difference in the fourth group in the test zone is observed for the NSHEX. This difference is reduced to −0.83% for the present method. However, a good agreement is obtained between the NSHEX and DNTR3D. 4. Conclusion A nodal SN method in the triangular-z geometry is developed for 3D neutron transport calculation, where arbitrary triangles are transformed into regular triangles via a coordinate transformation. The transverse procedure is used and the integral flux

H. Lu, H. Wu / Nuclear Engineering and Design 237 (2007) 830–839

is expanded by an orthogonal quadratic polynomial. The spatial distribution of the radial leakage in the radial coupling equation is approximated by a second-order polynomial. The NEFD algorithm is applied to establish a stable and efficient iterative scheme. Applicability of the DNTR3D is investigated with four NEACRP international benchmarks. The differences in keff between the DNTR3D and the Monte-Carlo method are less than 0.1% for all benchmarks. Good agreement of region averaged fluxes with the Monte-Carlo method shows that DNTR3D can calculate flux distributions accurately. Particularly, the local refined meshes can be obtained by using the triangles in the rod region and this can improve the calculation accuracy largely. Additionally, numerical results demonstrate the higher efficiency of the present method compared with the finite element spherical-harmonics method (TEPFEM), which is also developed for the neutron transport calculation with unstructured meshes. Acknowledgements This work is supported by National Natural Science Foundation of China grants No. 10475064 and National Key Laboratory of Reactor System Design Technology grants No. SYX-01-0509.

839

References Badruzzaman, A., 1985. An efficient algorithm for nodal-transport solutions in multidimensional geometry. Nucl. Sci. Eng. 89 (3), 281–290. Cao, L.Z., Wu, H.C., 2004. Spherical harmonics method for neutron transport equation based on unstructured-meshes. Nucl. Sci. Tech. 15 (6), 335–339. Gallagher, R.H., 1975. Finite Element Analysis Fundamentals. Prentice-hall, Englewood Cliffs, New Jersey, pp. 229–235. Ikeda, H., Takeda, T., 1994. A new nodal SN transport method for threedimensional hexagonal geometry. J. Nucl. Sci. Technol. 31 (6), 497–509. McGhee, J.M., Roberts, R.M., Morel, J.E., 1997. The DANTE Boltzmann transport solve: an unstructured mesh, 3D, spherical harmonics algorithm compatible with parallel computer architectures. In: Proceedings of Joint International Conference on Mathematical Methods and Supercomputing in Nuclear Applications, 5–10 October 1997, Saratoga Springs, New York, America. Ressel, K.J., Starke, G., 2000. A boundary functional for the least-squares finiteelement solution of neutron transport problems. SIAM J. Numer. Anal. 37 (2), 556–586. Sugino, K., Takeda, T., 1996. An improvement of the transverse leakage treatment for the nodal SN transport calculation method in hexagonal-z geometry. J. Nucl. Sci. Technol. 33 (8), 620–627. Takeda, T., Ikeda, H., 1991a. 3D neutron transport benchmarks. J. Nucl. Sci. Technol. 28 (7), 656–669. Takeda, T., Ikeda, H. 3D Neutron Transport Benchmarks, Technical Report OECD/NEA Committee on Reactor Physics (NEACRP), OSAKA University, NEACRP-L-330, March 1991b. Wagner, M.R., 1989. Three-dimensional nodal diffusion and transport theory methods for hexagonal-z geometry. Nucl. Sci. Eng. 103 (4), 377–391.