One-group steady state diffusion in a heterogeneous slab

Available online at www.sciencedirect.com

Mathematics and Computers in Simulation 80 (2010) 2142–2158

One-group steady state diffusion in a heterogeneous slab B. Ganapol Department of Aerospace and Mechanical Engineering, University of Arizona, Tucson, AZ, United States Available online 28 April 2010 Dedicated to Professor Ernest Mund, University of Brussels.

Abstract We derive a customized solution to the diffusion equation in a heterogeneous slab with an arbitrary contiguous source. While aspects of solution have previously been reported in the literature, its full explicit form as presented here by solving a 3-term recurrence, has not. We conclude with a series of exercises demonstrating the utility of the solution. © 2010 IMACS. Published by Elsevier B.V. All rights reserved. Keywords: Neutron diffusion; Criticality; Analytical solution; Recurrence

Ernest H. Mund is the quintessential professor. His knowledge of mathematical reactor physics and, more importantly, his ability to give that knowledge away is unsurpassed in academic circles. Moreover, as his colleague over these many years, I have greatly benefited from his inquiring mind and uncanny know-how. Ernest—May the years ahead of you and Monique be continuous in inspiration, discovery and happiness.

1. Introduction At PHYSOR’02, Ernest and I presented a new method of solution to the one-group neutron diffusion equation in highly heterogeneous 1D slab media for both steady state and transient applications [1]. Since that presentation, the solution has evolved to become arguably the most efficient solution to the 1D diffusion equation in heterogeneous media to date. Here, based on that work, I will investigate the method further with our focus confined to the steady state algorithm. We cast the solution in terms of a three-point boundary value recurrence relation, which we solve analytically in lieu of a tridiagonal solver. The algorithm makes use of the homogeneous and particular solutions to this difference equation to arrive at the solution for critical and subcritical heterogeneous systems (see also [3, Chapter 10] for the PN equations). We investigate the consequences of the numerical behavior of the algorithm primarily in anticipation of application to the time-dependent case [2]. We conclude however, that a tridiagonal solver provides the most robust algorithm.

E-mail address: [email protected]. 0378-4754/$36.00 © 2010 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.matcom.2010.04.008

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

2143

Fig. 1. 1D heterogeneous plane geometry.

2. Solution to the 1D-one-group diffusion equation 2.1. Diffusion theory in a heterogeneous medium In the following, we consider the one-group characterization of a heterogeneous reactor configuration shown in Fig. 1 where the appropriate steady state diffusion equation in slab j is Dj

d 2 φj (x) + (νf j − a j ) φj (x) = −Qj (x), dx2

(1)

and a spatially dependent source is included in each region. The boundary conditions at the free surface xa (a = 0 or n) can either be vacuum as characterized by the extrapolated endpoint φa (xa ) = 0

(2)

or zero current  dφa (x)  = 0. dx x=xa

(3)

The region interfaces are numbered from the left, beginning at zero. Continuity of the flux and current at the internal interfaces require  (φj−1 (x) − φj (x))x=x = 0, 

j−1

 dφj−1 (x) dφj (x)  −Dj−1 = 0. + Dj  dx dx x=xj−1

(4)

A more convenient form for (1) is d 2 φj (x) Qj (x) + Bj2 φj (x) = − , dx2 Dj

(5)

with Bj2 ≡

νf j − a j . Dj

(6)

Leading reactor physics texts give the solution to this, the most fundamental problem of reactor physics, but rarely in the generality of the heterogeneous system that we consider here. Most treatments are for homogeneous or reflected slabs, with full heterogeneity, relegated to later chapters and written in a matrix form that, as discussed in the following section, is considered to be the final result. The novelty here is that we treat full heterogeneity with the analytical simplicity of the single homogeneous slab.

2144

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

2.2. The standard solution The most straightforward solution to (5) is in terms of the homogeneous and particular solutions of the ODE in each region j φj (x) = aj eiBj x + bj e−iBj x + φp j (x).

(7)

A convenient alternative (as will be discussed) is to substitute hyperbolic or trigonometric functions for the exponentials depending upon the sign of Bj2 . When Eqs. (2)–(4) are satisfied, there result 2n-equations for the 2n-unknowns aj , bj that satisfy: ⎡ ⎤ a1 ⎢ ⎥ ⎢ b1 ⎥ ⎢ ⎥ ⎢ a2 ⎥ ⎢ ⎥ ⎢ ⎥ b2 ⎥ A·⎢ ⎢ ⎥=h ⎢ .. ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ an ⎦ bn where the banded A-matrix has elements e±iBj xj ,

e±iBj xj−1

and the h vector is in terms of the particular solution and its derivatives. The coupling coefficients aj and bj therefore come from matrix inversion of a 2n by 2n matrix. We obtain the largest eigenvalue of the system, keff , by replacing all fission cross sections by f j /keff and solving the homogeneous system giving the critical condition det A(keff ) = 0. For a single moderated fuel region and a two-region fuel/reflector system, the conditions for criticality are cos(B1 1 ) = 0 and D 2 B2 cot(B1 1 ) = 0 D 1 B1  respectively, where Bj = Bj2 , and x0 is the reactor centerline. For n large – of the order of 700 – this solution strategy can be rather costly for fixed source problems and is prohibitively so for time-dependent scenarios using numerical Laplace transform inversion. For this reason, we now propose a variant of the above using a more convenient customized solution methodology. tanh(B2 2 ) −

2.3. Fully customized solution The analysis begins with the development of the solution to (5) as − φj (x) = φj h+ j (x) + φj−1 hj (x) + φp j (x),

where the interfacial fluxes are φj ≡ φj (xj ),

φj−1 ≡ φj (xj−1 ),

(8)

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

2145

and h+ j (x)

≡

χj+ (x) χj

h− j (x)

,

≡

χj− (x) χj

(9)

.

Considering the two branches of the sine function depending upon the sign of Bj2 , we have ⎧ ⎨ sinh(Bj (x − xj−1 )), Bj2 < 0 + χj (x) ≡ ⎩ sin(B (x − x )), Bj2 > 0 j j−1 ⎧ ⎨ sinh(Bj (xj − x)), Bj2 < 0 − χj (x) ≡ ⎩ sin(B (x − x)), Bj2 > 0 j j

(10)

and χj = χ(βj j ),

(11)

where χ is either sin or sinh as appropriate. The first two terms of (8) are the homogenous solutions and the third term is the particular solution. The proposed solution will automatically satisfy flux continuity at slab interfaces by requiring the particular solution to satisfy φp j (xj ) ≡ 0, since

 h+ j (x)

≡

φp j (xj−1 ) ≡ 0,

x = xj , x = xj−1

1, 0,

(12) 

h− j (x)

≡

0, x = xj . 1, x = xj−1

The solution approach employed is the method of customized solutions and takes advantage of the fact that the homogeneous and particular solutions can take on a variety of forms and still be legitimate. Here, we choose a form that leads to numerical and analytical simplicity. When we apply the method of variation of parameters, the following particular solution results:    xj  x χj + −   −    +  φp j (x) = hj (x) dx Qj (x )hj (x ) + hj (x) dx Qj (x )hj (x ) . (13) D j Bj xj−1 x The lower limits of integration come from Eqs. (12). The advantage of the customized form of the solution is apparent when we satisfy the interfacial conditions of (4) to give the following three-point recurrence relation for the coupling coefficients between adjacent regions: φj + bj φj−1 + γj φj−2 = fj ,

2 ≤ j ≤ n,

(14)

with fj ≡ −[gp−j + γj gp+j−1 ], where gp+j ≡ and

 γj ≡

χj Dj Bj



xj

xj−1

Dj−1 Bj−1 χj−1

(15)

 dx Qj (x )h+ j (x ),



 χj , Dj Bj

gp−j ≡

χj D j Bj



xj

xj−1

bj ≡ −(Yj + γj Yj−1 ),

 dx Qj (x )h− j (x ),

(16)

(17)

2146

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

Yj ≡

⎧ ⎨ cosh(Bj j ),

Bj2 < 0

⎩ cos(B  ), j j

Bj2 > 0

(18)

For vacuum boundary conditions, we augment (14) with φ0 = φn = 0 making a recurrent boundary value problem. For symmetry at x0 , we have φ0 =

gp+1 φ1 + , Y1 Y1

which when introduced into (14) gives a modified equation for j = 2 φ2 + b2 φ1 = f2 , with the following replacements: b2 → b2 +

γ2 , Y1

f2 → f2 − γ2

gp−1 Y1

,

and now φ0 ≡ 0 becomes the natural condition. Similarly, for a zero current condition at xn modified n th equation is bn−1 φn−1 + γn−1 φn−2 = fn , with the modifications bn → bn +

1 , Yn

fn → fn −

gp+1 Yn

.

Again, φn ≡ 0 is now the natural boundary condition. Thus, for both vacuum and zero current conditions at the free surfaces, the recurrence of (14) is closed with φ0 ≡ 0 and φn ≡ 0. 2.4. Solution of the recurrence relation The recurrence given by (14) is one of the boundary value type, i.e., φ0 and φn are specified rather than starting values of φ0 and φ1 . We find the solution to this finite difference equation most conveniently through a tridiagonal solver, which indeed is the most robust and efficient numerical strategy. However, it is informative to seek a little known rather unconventional solution to the recurrence that leads to an explicit expression. Let φj be specified by φ j ≡ ρ j φ1 + g j φ 0 +

j 

αj,l fl .

(19)

l=2

Noting that φ0 ≡ 0, this equation simplifies to φj ≡ ρj φ1 +

j 

αj,l fl .

(20)

l=2

The summation is vacant when j = 0, 1. Therefore, for consistency, from (19), we must have g0 ≡ 1, g1 ≡ 0

and ρ0 ≡ 0, ρ1 ≡ 1.

(21)

From substitution of (19) into (14), we find that gj + bj gj−1 + γj gj−2 = 0,

(22)

ρj + bj ρj−1 + γj ρj−2 = 0,

(23)

αj,l + bj αj−1,l + γj αj−2,l = 0,

(24)

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

2147

and αj,j fj + αj,j−1 fj−1 + bj αj−1,j−1 fj−1 = fj ,

(25)

to maintain equality. Therefore, to satisfy (24), we have αj,j = 1,

αj,j−1 = −bj ,

(26)

which serve as the starting values for the recurrence of (25). Three initial value recurrence relations now represent the original boundary value recurrence relation. The expression given by (19) consists of two homogeneous solutions to the recurrence relation (first two terms) and a particular solution. Noting that if αj,l = νl gj + εl ρj ,

(27)

the recurrence relation is satisfied and from (25), we can show that ρl+1 + bl+1 ρl γl+1 ρl−1 = − , gl ρl+1 − gl+1 ρl gl ρl+1 − gl+1 ρl

νl =

gl+1 + bl+1 gl γl+1 gl−1 εl = − = . gl ρl+1 − gl+1 ρl gl ρl+1 − gl+1 ρl

(28)

This is nothing more than variation of parameters applied to a finite difference equation. The denominator in these expressions has a convenient explicit form according to the following analysis. When we multiply Eqs. (22) and (23) by ρl and gl , respectively, and subtract, there results the simple two-point recurrence relation Tl − γl+1 Tl−1 = 0, where Tl ≡ gl ρl+1 − ρl gl+1 . Upon solution of this simple recurrence, there results Tl =

l+1 

γj =

j=2

τ1 , τl+1

where τl ≡

D l Bl , χl

and (27) becomes τl αj,l = [ρj gl−1 − gj ρl−1 ]. τ1

(29)

Finally, for the zero boundary condition at xn , we find from (20) evaluated at j = n 1 αn,l fl , ρn n

φ1 = −

(30)

l=2

and (20) becomes φj ≡

j  l=2

with εj l ≡

τl τ1

ρj  αn,l fl , ρn n

αj,l fl +

l=2

  ρj [ρj gl−1 − gj ρl−1 ] (j − l) − [ρn gl−1 − gn ρl−1 ] . ρn

(31)

2148

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

Now the explicit expression for the interfacial flux is  xl n  1 dx [εj l χl− (x ) − γl εj l+1 χl+ (x )] Ql (x ). φj = − Dl Bl xl−1

(32)

l=1

This elegant solution has a shortcoming however. The recurrence can become unstable and lead to incorrect results as we discuss below. In addition, we can obtain explicit expressions for the solution to the homogeneous recurrence relations. If   kj yj ≡ , (33) kj+1 where k is ρ or g respectively, then the recurrences become yj + Aj yj+1 = 0 with

⎡b

j+2

⎢ A ≡ ⎣ γj+2

1 ⎤ γj+2 ⎥ ⎦

−1 and

  1 y0 ≡ 0

(34)

(35)

0   0 y0 ≡ . 1

or

(36)

We easily find the solution to the recurrence of (34) in terms of the fundamental matrix solution ⎤−1 ⎡ l−1  kl = (−1)l [1 0]⎣ Aj ⎦ y0 .

(37)

j=0

Thus, the inner products give ⎤−1   ⎡ l−1 1  gl = (−1)l [1 0]⎣ Aj ⎦ 0 j=0 (38)

⎤−1   . 0 ρl = (−1)l [1 0]⎣ Aj ⎦ 1 j=0 ⎡

l−1 

2.5. Summary of the analytical solution in 1D heterogeneous plane geometry For completeness, we summarize all the elements of the analytical solution in the following: (i) Customized solution χj φj (x) = Bj Dj

 x

xj

 +  dx h− j (x )hj (x)Qj (x ) +



x

xj−1

  −  dx h+ j (x )hj (x)Qj (x )

⎧ − −  ⎨ [εj l h+ j (x) + εj−1 l hj (x)]hl (x )

 xl n  1 dx − ⎩ +[ε Bl Dl xl−1 l=1

⎫ ⎬

+ − +  ⎭ j l+1 hj (x) + εj−1 l+1 hj (x)]hl (x )

(39) Ql (x )

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

2149

(ii) Homogeneous solutions

h+ j (x) ≡

χj+ (x) χj

 χj+ (x)

≡ 

χj− (x) ≡

χj− (x)

h− j (x) ≡

,

χj

sinh(Bj (x − xj−1 )),

Bj2 < 0

sin(Bj (x − xj−1 )),

Bj2 > 0

sinh(Bj (xj − x)),

Bj2 < 0

sin(Bj (xj − x)),

Bj2 > 0

(iii) Coupling coefficients

εj l

τl ≡ τ1

τl ≡



 ρj [ρj gl−1 − gj ρl−1 ] (j − l) − [ρn gl−1 − gn ρl−1 ] , ρn

D l Bl . χl

(iv) Solutions to homogeneous recurrences ⎤−1   1 gl = (−1)l [1 0]⎣ Aj ⎦ 0 j=0 ⎤−1   ⎡ l−1  0 ρl = (−1)l [1 0]⎣ Aj ⎦ 1 j=0 ⎡

l−1 

Fig. 2. Configuration and nuclear properties for the MYRRHA ADS concept used in the numerical comparison of tridiagonal and recurrence solvers.

2150

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

⎡b

j+2

A ≡ ⎣ γj+2 −1

1 ⎤ γj+2 ⎦ . 0

While a fully explicit analytical solution is elegant, it only has numerical relevance if evaluated by a precision enhanced computer algebra package. It is of little limited value otherwise. 3. Determination of keff We obtain the multiplication factor keff most easily by requiring the reactor to support a flux in the absence of any sources. This is realized if φ1 is required to be non-zero when fl is zero in (20). Thus, for criticality ρn (keff ) = 0.

(40)

This is the identical relation to (14) with no sources and φ1 normalized to unity. We therefore conduct a search for the largest value of keff producing a zero of ρn , where we have divided each fission cross section by keff . The specification of ρn can most easily be accomplished via (23). 4. Numerical implementation We have implemented both the tridiagonal and recurrence methods of solution for the special case of uniform regional sources only. The most convenient form of the solution for this case is − φj (x) = (φj − φp j ) h+ j (x) + (φj−1 − φp j ) hj (x) + φp j ,

(41)

where the particular solution is now φp j =

Qj . Dj βj2

(42)

In the following sections, we will test the analytical solution and its implementation through a series of exercises. 4.1. Internal numerical implementation checks Several internal checks of the numerical implementation are possible. We apply the one-group model to the MYRRHA accelerator driven system for waste disposal and energy production shown in Fig. 2 for this purpose. Table 1 specifies the dimensions and nuclear properties of each region. We use the tridiagonal solver exclusively in this section. The first check is for preservation of horizontal symmetry. Fig. 3 shows the fluxes for reversed MYRRHA half-cores. For this case, we reversed the zero current and zero flux boundary conditions as well as the regions. The results are mirror flux images as expected. The second check is for full and half-core symmetries as shown in Fig. 4, where we have superimposed the half core flux with zero current at the centerline on the full core configuration. Over the region of coincidence, there is no difference. We therefore conclude that the boundary conditions are implemented correctly. Table 1 Nuclear properties of the MYRRHA ADS concept used in the numerical comparison of tridiagonal and recurrence solvers. Region

j (cm) D (cm) νΣf (cm−1 ) Σa (cm−1 ) Q (n/cm3 )

1

2

3

4

6.405 3.1714 0 0.012441 0.7119 × 102

10.52 2.2438 0.049371 0.041823 0

10.974 2.2455 0.049701 0.041632 0

77.079 2.6501 0 0.010655 0

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

Fig. 3. Mirror symmetry of the MYRRHA 1-GP ADS.

Fig. 4. Comparison of full and half-core symmetries.

Fig. 5. Artificial subdivision of homogeneous regions.

2151

2152

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

Fig. 6. Introduction of central void.

Fig. 7. Comparison of tridiagonal and recurrence solvers for MYRRHA ADS concept.

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

Fig. 8. Degradation of recurrence solver with many regions.

Fig. 9. Control plate between fuel bundles.

Fig. 10. Control rod between two 8-plate assemblies (39 regions).

2153

2154

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

The penultimate check is for proper programming of the internal boundary conditions. We artificially divide each homogeneous region into 10 sub-regions, which should leave the flux distribution unaffected. Fig. 5 verifies our expectation. The final check is for the introduction of a near void at the core center. We simulated a void by a large diffusion coefficient (1000 cm) and a small absorption cross section (10−14 cm−1 ). As shown in Fig. 6, the flux in the void remains constant at the centerline value after which the flux behaves as it should in the original core. All internal checks are satisfied indicating that the implementation is ready for trial.

4.2. Comparison of solution strategies We compared the solutions given by the tridiagonal solver and recurrence relation for MYRRHA core. Rather than using the recurrence relation as specified by (14) directly or the analytical result of (39), we found a more numerically stable solution by the following analysis to reduce the recurrent boundary value problem to a recurrent initial value one. Starting with (14), we seek the solution in the form φj = uj + ρj φ1 .

(43)

Fig. 11. Failure of recurrence relation.

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

2155

Substitution into (14) leads to uj + bj uj−1 + γj uj−2 = fj

(44)

for 2 ≤ j ≤ n with u0 = u1 = 0 and (22) for φj . From the vacuum boundary condition introduced into (43) un φ1 = − ρn

(45)

and therefore φj = uj −

un ρj . ρn

(46)

This formulation is more numerically stable than either (14) or (20). Fig. 7 shows a direct comparison of the tridiagonal and recurrence solvers indicating nearly perfect agreement. As one can anticipate however, the recurrence relation will eventually degrade with the number of iterations as a result of round off error. This is indeed shown to be the case in Fig. 8 where each of the homogeneous regions was artificially partitioned into 100 regions for a critical reactor. We clearly observe the degradation of the recurrence with the error increasing from 10−16 to 10−12 . While for this case the degradation is insignificant, it does indicate that for more sensitive configurations, the recurrence will eventually fail as will soon be demonstrated. In conclusion, for a modest number of regions of the order of less than 400, either solver seems appropriate. Thus, for more regions, the tridiagonal solver is the best choice.

Fig. 12. Comparison of plate and homogeneous geometry.

2156

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

4.3. Selected demonstrations The first demonstration is for a control rod (plate) situated between two fuel bundles each containing 8-fuel plates as shown in Fig. 9. Water fills the channels between the elements as well as the control rod channel. Each channel also

Fig. 13. Control rod insertion.

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

2157

contains a uniform source representing the source of neutrons slowing to thermal energies. The rod contains about a dollar’s worth of reactivity and the “rod out keff ” of the system is 0.11515. Fig. 10 shows the resulting fluxes for four control rod positions and 2 fuel rod bundles. The control rod effectively suppresses the flux peak as required. Looking closely, one observes the individual detail of flux peaks within the bundles. The total number of regions for this case was 39. Fig. 11 shows a close-up of one and two bundles without control rods as calculated by both solvers. Now we see

Fig. 14. Increase in the fission cross section.

2158

B. Ganapol / Mathematics and Computers in Simulation 80 (2010) 2142–2158

the fuel and water channels more clearly. In addition, the failure of the recurrence relation is particularly spectacular. This is the recurrence of (14) which is less stable than (23) or (44). The same fate awaits these equations but at about three times the number of regions. This demonstration ends with a comparison of four and eight bundles (see Fig. 12) indicating the effect of the fuel homogenization in the channel. Apparently, homogenizing the fuel, underpredicts absorption. The final demonstration is to show that the method remains faithful for the extreme case of 20 assemblies and a total of 723 regions. For this case, we place a control rod between each assembly, which is completely inserted or out in the first two cases and alternatively inserted or out in the last case. From Fig. 13, the simulation seems to be quite representative even for the large number of regions considered. Finally, Fig. 14 shows the effect of increasing the fission cross section, which generates a higher keff . The flux in the fuel increases with fission cross section with the last plate showing a nearly critical flux. Again, the simulation seems entirely reasonable. 5. A final thought The one-group, diffusion equation has been solved analytically in a 1D heterogeneous medium. The solution can find use in Monte Carlo, diffusion acceleration of transport computations, testing of simple numerical algorithms and most importantly in the classroom. We have finally eliminated the drudgery of past solutions for heterogeneous systems. Extension of this work to other 1D geometries, multigroup and time-dependence will be subjects of future presentations with my good friend and colleague Ernest Mund. References [1] B. Ganapol, E.H. Mund, P. Ravetto, M.M. Rostagno, Multigroup diffusion kinetics benchmark of an ADS system in slab geometry, in: Proceedings of PHYSOR’02 7B-03, 2002. [2] B. Ganapol, One-group time dependent diffusion with delayed neutrons in a heterogeneous medium, Transactions of the American Nuclear Society 87 (2002) 162–165. [3] B. Davison, Neutron Transport Theory, Oxford University Press, London, 1957.