An alternative solution of the neutron diffusion equation in cylindrical symmetry

Annals of Nuclear Energy 38 (2011) 1140–1143

Contents lists available at ScienceDirect

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

An alternative solution of the neutron diffusion equation in cylindrical symmetry Saed Dababneh a,b,⇑, Kafa Khasawneh a, Zaid Odibat a a b

Prince Abdullah Bin Ghazi Faculty of Science and IT, Al-Balqa Applied University, Salt 19117, Jordan Jordan Nuclear Regulatory Commission, P.O. Box 2587, Amman 11941, Jordan

a r t i c l e

i n f o

Article history: Received 15 August 2010 Received in revised form 27 November 2010 Accepted 9 December 2010 Available online 3 January 2011 Keywords: Neutron diffusion Homotopy perturbation method Flux calculation Critical system

a b s t r a c t Alternative analytical solutions of the neutron diffusion equation for both infinite and finite cylinders of fissile material are formulated using the homotopy perturbation method. Zero flux boundary conditions are investigated on boundary as well as on extrapolated boundary. Numerical results are provided for one-speed fast neutrons in 235U. The results reveal that the homotopy perturbation method provides an accurate alternative to the Bessel function based solutions for these geometries. Ó 2010 Elsevier Ltd. All rights reserved.

In a recent work (Khasawneh et al., 2009), the nonlinear analytical homotopy perturbation method (HPM) has been introduced as a means to solve nuclear reactor physics problems; namely the one-group neutron diffusion equation in hemispherical symmetry. In tandem with that, this work provides a description for the application of HPM to solve time-independent neutron diffusion equations for cylindrical geometries at different boundary conditions. These boundary conditions are the zero flux at boundary ZF and the zero flux at extrapolated boundary EBC. The numerical results obtained using HPM are compared to the outcome of canonical Bessel-function based solutions, both computed in this work. The general theory including the fundamental HPM mathematical derivations, along with its application to infinite and finite cylinders, are given in Section 2. Numerical results are subsequently provided for critical radius calculations, as well as for flux distribution in the critical infinite and finite cylinders (Section 3).

which is easily solved. This method does not require a small parameter in an equation, nor it incorporates unjustified assumptions. It has a significant advantage in that it provides an analytical approximate solution to a wide range of nonlinear problems in applied sciences, see for example Ganji et al. (2007), Chowdhury and Hashim (2008), Marinca and Heriçanu (2008), Cveticanin (2009) and references therein. The HPM yields the solution in terms of a rapid convergent series with easily computable components. For linear equations the method gives exact solution, and for nonlinear equations it provides approximate solution with fairly good accuracy (He, 2006a). Analysis of the method and different illustrative examples (He, 1999, 2000, 2003, 2006b; Odibat, 2007) present HPM as a promissing means to solve general differential equations. The method was first introduced and used to solve nuclear reactor physics problems in a recent work by the authors (Khasawneh et al., 2009). The HPM is applied to neutron diffusion in a bare infinite and finite cylinders, in Sections 2.2 and 2.3, respectively.

2. Theory

2.2. Bare infinite cylinder

2.1. The homotopy perturbation method (HPM)

In this section we consider an infinite cylindrical reactor of radius a. Since the cylinder is infinite, the flux in this reactor is a function of r only, and the time-independent diffusion equation can be written as

1. Introduction

The HPM, which is a coupling of a homotopy technique and a perturbation technique, deforms continuously to a simple problem

r2 /ðrÞ þ B2 /ðrÞ ¼ 0; ⇑ Corresponding author at: Prince Abdullah Bin Ghazi Faculty of Science and IT, Al-Balqa Applied University, Salt 19117, Jordan. E-mail addresses: [email protected], [email protected] (S. Dababneh), [email protected] (K. Khasawneh). 0306-4549/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.anucene.2010.12.011

ð1Þ

where B2, the buckling of the reactor, is given by

B2 ¼

m Rf  Ra D

:

ð2Þ

S. Dababneh et al. / Annals of Nuclear Energy 38 (2011) 1140–1143

1141

Substituting the Laplacian in cylindrical coordinates and using x = Br, Eq. (1) reduces to

x2 /00 ðxÞ þ x/0 ðxÞ þ x2 /ðxÞ ¼ 0:

ð3Þ

In order to solve this equation using HPM, we construct the following homotopy

Hðu; pÞ ¼ x2 u00 ðxÞ þ xu0 ðxÞ þ px2 uðxÞ ¼ 0;

p 2 ½0; 1:

ð4Þ

The variation of p from zero to unity corresponds to the variation of Hðu; pÞ form u0(x); which is an initial approximation obtained by solving the simple equation x2u00 (x) + xu0 (x) = 0, to /(x); which yields the solution of Eq. (3). The basic assumption of HPM is that the solution of Eq. (4) can be expressed as a power series in p,

uðxÞ ¼ u0 ðxÞ þ pu1 ðxÞ þ p2 u2 ðxÞ þ    :

Substituting Eq. (5) into the homotopy given in Eq. (4) and solving the resulting equations, then considering the terms with identical powers of p, we obtain a set of equations of the form

8 0 p : x2 u000 ðxÞ þ xu00 ðxÞ ¼ 0; > > > > > 1 2 00 0 2 > > > p : x u1 ðxÞ þ xu1 ðxÞ ¼ x u0 ðxÞ; > < 2 2 00 p : x u2 ðxÞ þ xu02 ðxÞ ¼ x2 u1 ðxÞ; > > > .. > > > . > > > : k 2 00 p : x uk ðxÞ þ xu0k ðxÞ ¼ x2 uk1 ðxÞ;

/1 ð0Þ : 0; /2 ð0Þ : 0;

ð6Þ

  2 1 d d/ðr; zÞ d /ðr; zÞ r þ þ B2 /ðr; zÞ ¼ 0: 2 r dr dr dz /ðr; zÞ ¼ RðrÞZðzÞ;

/k ð0Þ : 0:

ð15Þ

ð7Þ

8 u1 ðxÞ ¼  A4 x2 ; > > > > A 4 > > < u2 ðxÞ ¼ ð16Þð4Þ x ;

ð8Þ

.. .

ð16Þ

we can rewrite the diffusion equation (Eq. (15)) as

where A is a constant. Similarly, we obtain

  2 1 d dR 1d Z r þ B2 ¼  : rR dr dr Z dz2

ð17Þ

The left hand side of Eq. (17), which is a function of r only, appears to depend on the function at the right hand side, which is a function of z alone. We resolve this by setting each side equal to the same constant, thus arriving at the eigenvalue equations 2

k

uk ðxÞ ¼ Að1Þ x2k : 4k k!k!

The final solution of Eq. (4) can be expressed in the form of a power series in p

A 4

uðxÞ ¼ A  p x2 þ p2

A 1 x4  p3 x6 þ    : 4  16 4  16  36

ð9Þ

Then the solution of Eq. (3) is

  x2 x4 x6 /ðxÞ ¼ lim A 1  p þ p2  p3 þ  ; p!1 4 4  16 4  16  36   x2 x4 x6  þ  ; ¼A 1 þ 4 4  16 4  16  36

ð10Þ ð11Þ

1 X Að1Þk k¼0

4k k!k!

ð12Þ

ðBrÞ2k :

ð13Þ

or 1 X Að1Þk k¼0

4k k!k!

ð18Þ ð19Þ

The solution of the axial part given in Eq. (18) is

Z m ðzÞ ¼ cosðmzÞ:

ð20Þ

For this finite cylinder, and in order to solve the radial part we can rewrite Eq. (19) as

  1 d dR r þ k2m R ¼ 0; r dr dr

ð21Þ

k2m ¼ B2  m2 :

2k

x ;

1d Z ¼ m2 ; Z dz2   1 d dR þ B2 ¼ m2 : r rR dr dr

where

and in a closed form as

/ðBrÞ ¼

ð14Þ

2

Then applying separation of variables

u0 ðxÞ ¼ A;

/ðxÞ ¼

r2 /ðr; zÞ þ B2 /ðr; zÞ ¼ 0;

the buckling B is again given by Eq. (2). Writing the diffusion equation (Eq. (14)) using the Laplacian in cylindrical coordinates, we get

/0 ð0Þ : finite;

The solution of the first component of Eq. (6) is

> > > > > > :

Fig. 1. Circular cylinder of radius a and height H.

ð5Þ

Clearly, this flux should satisfy the boundary conditions. 2.3. Bare finite cylinder It is worth now investigating the critical radius and flux distribution in a finite cylinder of height H and radius a (Fig 1). The symmetry in this reactor implies that the flux is a function of both r and z. Then the time-independent diffusion equation is written as

ð22Þ

Since there is no h and no z dependence in Eq. (21), then the first term in the equation is simply the Laplacian (with only r dependence), making Eq. (21) similar to Eq. (1), with R replacing / and km replacing B. Then considering xm = kmr, and in analogy with Eq. (3), we get

x2m /00 ðxm Þ þ xm /0 ðxm Þ þ x2m /ðxm Þ ¼ 0:

ð23Þ

Hence the solution of Eq. (23) can be obtained using the HPM technique. The same procedure of Subsection 2.2 yielding Eq. (13) has been followed and the result is obviously (compare Eq. (13))

Rm ðkm rÞ ¼

1 X Am ð1Þk k¼0

4k k!k!

ðkm rÞ2k :

ð24Þ

The final solution of the diffusion equation (Eq. (14)) is therefore

1142

/ðr; zÞ ¼

S. Dababneh et al. / Annals of Nuclear Energy 38 (2011) 1140–1143 1 X

am Rm ðkm rÞZ m ðzÞ;

ð25Þ

m¼

m¼0

where Rm(kmr) and Zm(z) are given by Eqs. (24) and (20), respectively. Again, this flux should satisfy the appropriate boundary conditions. 3. Numerical results In order to provide numerical results, we follow the example proposed originally by Cassell and Williams (2004) and later adopted by us for the solution of the diffusion equation in hemispherical symmetry (Khasawneh et al., 2009); namely 1 MeV neutrons diffusing in pure 235U. The following data will correspondingly be used:

m ¼ 2:42; rf ¼ 1:336b; rs ¼ 5:959b; rc ¼ 0:153b; Nc ¼ 0:0478  1024 atoms cm3 :

np ; Hc

n ¼ 1; 2; 3; . . . :

ð30Þ

In order to account for the first zero of the solution of the axial part, i.e. n = 1, we may write

ZðzÞ ¼ cos



p

Hc

 z :

ð31Þ

On the other hand, applying the ZF boundary condition on the radial part

Rm ðkm ac Þ ¼ 0;

ð32Þ

the first zero calculated using Eq. (24) is km ac = 2.40483. It is worthwhile here to notice that though this dimensionless value is similar to the one obtained for the infinite cylinder, but the critical radius itself depends on the critical height through the values of m and consequently km. The corresponding buckling B2, according to Eq. (22), is



2  2 2:40483 p þ : ac Hc

The diffusion coefficient D = 1/(3Ncrtotal) is then equal to 0.9363 cm.

B2 ¼

3.1. Critical radius calculations

This is the same expression for the buckling obtained using standard Bessel functions (Duderstadt and Hamilton, 1976), establishing HPM as an efficient alternative for calculating the critical dimensions of fissile systems. Table 2 lists calculated critical radii of a finite cylinder for a set of chosen critical height values. It is worthwhile to mention that, with accuracy to five significant figures, the smallest critical dimensions would be Hc = 10.529 cm and a corresponding ac,ZF of 665.253 cm. On the other hand, for relatively large critical height

In order to calculate the critical radius for an infinite cylinder, a computer code based on the HPM derived Eq. (13) has been written to calculate, by iteration, the critical value of Ba which was found to be 2.40483. If the ZF boundary condition is applied, the critical radius ac,ZF is 8.05886 cm, while applying the EBC yields a critical radius ac,EBC of 6.18626 cm. These are the same critical radii calculated using Bessel functions (Duderstadt and Hamilton, 1976), as illustrated in Table 1. This is the first indication of how perfect the HPM can reproduce numerical results when compared to other canonical mathematical techniques. It is worth mentioning here that the analytical form of the HPM solution (Eq. (13)) is identical to Bessel polynomials, in this particular case. On the other hand, and for a finite cylinder, both critical radius and height need to be calculated. In order to do that, boundary conditions should be applied on both eigenfunction equations (Eqs. (20) and (24)). According to ZF boundary condition, for the flux to vanish on the top and bottom surfaces of the cylinder, i.e.

  Hc ¼ 0; / r;  2

ð26Þ

we need to apply boundary conditions on the axial part

  Hc Zm  ¼ 0; 2   Hc ¼ 0; cos m 2

ð27Þ ð28Þ

ð33Þ

Table 2 Chosen heights of a finite cylinder and the corresponding critical radii. Hc (cm)

ac,ZF (cm)

ac,EBC (cm)

10.529 10.530 11.000 12.000 13.000 14.000 15.000 20.000 30.000 40.000 50.000 100.00 200.00 400.00 600.00 800.00 1000.0

665.253 439.235 27.8171 16.7953 13.7383 12.2267 11.3144 9.47880 8.60656 8.35371 8.24399 8.10420 8.07035 8.06195 8.06040 8.05986 8.05961

663.381 437.362 25.9445 14.9227 11.8657 10.3541 9.44178 7.60620 6.73396 6.48111 6.37139 6.23160 6.19775 6.18935 6.18780 6.18726 6.18701

which is valid only if

m

Hc p ¼n ; 2 2

n ¼ 1; 2; 3 . . . :

ð29Þ Table 3 Normalized flux across an infinite cylinder of fissile material.

Then the eigenvalues are

Table 1 Critical radius ac (in cm) for infinite cylinder using different boundary conditions. Identical results are obtained using HPM and Bessel functions. BC

ac (cm)

ZF EBC

8.05886 6.18626

r/ac,ZF

Normalized flux

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

1.000000 0.985594 0.942999 0.874050 0.781711 0.669929 0.543446 0.407578 0.267958 0.130272 0.000000

S. Dababneh et al. / Annals of Nuclear Energy 38 (2011) 1140–1143

1143

where C is a normalization constant. A computer code based on Eq. (34) has been written. Fig. 3 illustrates the general behaviour of the calculated flux for a finite cylinder of a given set of critical dimensions. As expected from the symmetry, the maximum value of the normalized flux is at the center of the cylinder (r = 0, z = 0) and starts to decrease when moving towards any surface. Evidently, the flux vanishes at r = ac,ZF and z ¼ H2c . 4. Conclusion The homotopy perturbation method, applied to the neutron diffusion problem in cylindrical symmetry, can exactly reproduce the results obtained using Bessel functions. This conclusion supports the earlier findings of Khasawneh et al. (2009) regarding the application of HPM to solve nuclear reactor physics problems in different geometries. In the particular case of cylindrical symmetry, the analytical form of the solution is identical to that obtained using canonical approaches.

Fig. 2. Flux distribution of an infinite cylinder.

Acknowledgment The financial support provided by King Abdullah II Fund for Development KAFD is gratefully acknowledged by K. Kh. References

Fig. 3. Flux distribution across a core of finite cylinder.

values, the corresponding critical radius approaches 8.05886 cm, which is the value obtained for an infinite system. 3.2. Flux distribution In order to calculate the flux distribution for an infinite cylinder, we wrote a computer code based on Eq. (13), where the constant A has been normalized to unity and ZF boundary condition was applied. Table 3 lists the normalized flux across the system for different values of r/ac,ZF. Fig. 2 illustrates the flux behavior in this fissile system. The maximum value of the flux is naturally at the axis of the cylinder (r = 0), and decreases as we move towards the surface. This behavior is expected from the symmetry of the system. Finally, the flux across a finite cylinder is given by

   2k ! 1 X ð1Þk 2:40483 pz ; /ðr; zÞ ¼ C r cos k a Hc c 4 k!k! k¼0

ð34Þ

Cassell, J.S., Williams, M.M.R., 2004. A solution of the neutron diffusion equation for a hemisphere with mixed boundary conditions. Annals of Nuclear Energy 31, 1987–2004. Chowdhury, M.S.H., Hashim, I., 2008. Analytical solutions to heat transfer equations by homotopy-perturbation method revisited. Physics Letters A 372, 1240–1243. Cveticanin, L., 2009. Application of homotopy-perturbation to non-linear partial differential equations. Chaos, Solitons & Fractals 40, 221–228. Duderstadt, J.J., Hamilton, L.J., 1976. Nuclear Reactor Analysis. John Wiley & Sons. Ganji, D.D., Afrouzi, G.A., Talarposhti, R.A., 2007. Application of variational iteration method and homotopy-perturbation method for nonlinear heat diffusion and heat transfer equations. Physics Letters A 368, 450–457. He, J.H., 1999. Homotopy perturbation technique. Computer Methods in Applied Mechanics and Engineering 178, 257–262. He, J.H., 2000. A coupling method of a homotopy technique and perturbation technique for non-linear problems. International Journal of Non-Linear Mechanics 35, 37–43. He, J.H., 2003. Homotopy perturbation method: a new nonlinear analytical technique. Applied Mathematics and Computation 135, 73–79. He, J.H., 2006a. Some asymptotic methods for strongly nonlinear equations. International Journal of Modern Physics B 20, 1141–1199. He, J.H., 2006b. Homotopy perturbation method for solving boundary value problems. Physics Letters A 350, 87–88. Khasawneh, K., Dababneh, S., Odibat, Z., 2009. A solution of the neutron diffusion equation in hemispherical symmetry using the homotopy perturbation method. Annals of Nuclear Energy 36, 1711–1717. Marinca, V., Heriçanu, N., 2008. Application of optimal homotopy asymptotic method for solving nonlinear equations arising in heat transfer. International Communications in Heat and Mass Transfer 35, 710–715. Odibat, Z.M., 2007. A new modification of the homotopy perturbation method for linear and nonlinear operators. Applied Mathematics and Computation 189, 746–753.