- Email: [email protected]
A solution of the neutron diffusion equation in hemispherical symmetry using the homotopy perturbation method
Annals of Nuclear Energy 36 (2009) 1711–1717
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A solution of the neutron diffusion equation in hemispherical symmetry using the homotopy perturbation method Kafa Khasawneh a, Saed Dababneh a,b,*, Zaid Odibat a a b
Prince Abdullah Bin Ghazi Faculty of Science and IT, Al-Balqa Applied University, P.O. Box 7051, 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 June 2009 Received in revised form 8 September 2009 Accepted 9 September 2009 Available online 7 October 2009
a b s t r a c t The homotopy perturbation method is used to formulate a new analytic solution of the neutron diffusion equation both for a sphere and a hemisphere of fissile material. Different boundary conditions are investigated; including zero flux on boundary, zero flux on extrapolated boundary, and radiation boundary condition. The interaction between two hemispheres with opposite flat faces is also presented. Numerical results are provided for one-speed fast neutrons in 235U. A comparison with Bessel function based solutions demonstrates that the homotopy perturbation method can exactly reproduce the results. The computational implementation of the analytic solutions was found to improve the numeric results when compared to finite element calculations. Ó 2009 Published by Elsevier Ltd.
1. Introduction
2. Theory
There has been variety of methods used over the last decades to solve the one-group neutron diffusion equation. Up to our knowledge, the nonlinear analytical homotopy perturbation method (HPM) has never been tested as a means to solve reactor physics problems. The merit of this contribution is that it provides a description for the application of HPM to solve time-independent neutron diffusion equations, for different geometries at different boundary conditions. These boundary conditions are the zero flux at the boundary ZF, the zero flux at the extrapolated boundary EBC, and the radiation boundary condition RBC. The numerical results obtained in this work are compared to the outcome of other Bessel function based solutions, both obtained in this work and reported by Cassell and Williams (2004), in addition to transport theory calculations reported in the same reference. The aim of this contribution is to improve, or at least widen the methodologies used to investigate the flux behavior inside the reactor core. The general theory including the essential HPM derivations, and its application to the sphere and hemisphere are given in Section 2. The RBC boundary condition is also detailed in order to set the grounds for the calculations. Numerical results are provided for critical radius calculations, as well as for flux distribution in the critical hemisphere (Section 3). The accompanying technical issues faced during the computations are described together with their corresponding remedies.
2.1. The homotopy perturbation method (HPM)
* Corresponding author. Address: Prince Abdullah Bin Ghazi Faculty of Science and IT, Al-Balqa Applied University, P.O. Box 7051, Salt 19117, Jordan. Tel.: +962 7 95606613; fax: +962 6 5341608. E-mail addresses: [email protected], [email protected] (S. Dababneh). 0306-4549/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.anucene.2009.09.001
The HPM, proposed first by He (1999, 2000) for solving linear and nonlinear differential and integral equations, has been the subject of extensive analytical and numerical studies. The method, which is a coupling of a homotopy technique and a perturbation technique, deforms continuously to a simple problem which is easily solved. This method, which does not require a small parameter in an equation, nor it incorporates unjustified assumptions, 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, 2006). Analysis of the method and different illustrative examples (He, 1999, 2000, 2003, 2006; Odibat, 2007) present HPM as a promissing means to solve general differential equations. The HPM is applied to neutron diffusion in a bare sphere and hemisphere, in Sections 2.2 and 2.3, respectively. 2.2. Bare sphere We will consider first the rather simple example of a bare spherical reactor of radius a. The spherical symmetry of the system implies that the flux is a function of r only. The time-independent diffusion equation can thus be written as
1712
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
r2 /ðrÞ þ B2 /ðrÞ ¼ 0;
ð1Þ
where
B2 ¼
mRf Ra
ð2Þ
D
is the buckling of the reactor. Substituting the Laplacian in spherical coordinates and using x ¼ Br, we get 2
x2
d /ðxÞ 2
dx
þ 2x
d/ðxÞ þ x2 /ðxÞ ¼ 0: dx
ð3Þ
In order to solve Eq. (3) using HPM, we construct the following homotopy:
Hðu; pÞ ¼ x2 u00 ðxÞ þ 2xu0 ðxÞ þ px2 uðxÞ;
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 x2 u00 ðxÞ þ 2xu0 ðxÞ ¼ 0, to uðxÞ; 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, 2
uðxÞ ¼ u0 ðxÞ þ pu1 ðxÞ þ p u2 ðxÞ þ
ð5Þ
Therefore, the solution of Eq. (3) can be readily obtained as
uðxÞ ¼ lim½u0 ðxÞ þ pu1 ðxÞ þ p2 u2 ðxÞ þ : p!1
u0 ð0Þ : finite; u1 ð0Þ : 0; u2 ð0Þ : 0;
ð8Þ
where A and C are constants. The solution should be finite at r = zero; thus u0 ðxÞ reduces to
u0 ðxÞ ¼ C:
Writing the Laplacian in spherical coordinates and using
@ 2 /ðr; lÞ 2 @/ðr; lÞ 1 @ 2 @/ðr; lÞ þ ð1 l Þ þ @r 2 r @r r2 @ l @l þ B2 /ðr; lÞ ¼ 0;
where the buckling B2 is again given by Eq. (2). Applying the separation of variables,
ð9Þ
2
ð15Þ ð16Þ
In order to solve the radial part using HPM, we first consider x ¼ Br and rewrite Eq. (15) as 2
x2
d RðxÞ 2
dx
þ 2x
dRðxÞ þ ðx2 nðn þ 1ÞÞRðxÞ ¼ 0: dx
ð17Þ
Then we construct the homotopy 2
Hðu; pÞ ¼ x2
d uðxÞ 2
dx
þ 2x
duðxÞ þ ðpx2 nðn þ 1ÞÞuðxÞ; dx
p 2 ½0; 1: ð18Þ
ð10Þ
Following the same manner as in Section 2.2, if we substitute Eq. (5) into the homotopy (18), and solve the resulting equations, we get the following solution:
Consequently, the final solution of Eq. (1), in a closed form, is given by 1 X ð1Þk C /ðBrÞ ¼ ðBrÞ2k : ð2k þ 1Þ! k¼0
ð14Þ
r 2 d RðrÞ 2r dRðrÞ þ þ B2 r 2 ¼ nðn þ 1Þ; RðrÞ dr RðrÞ dr 2 1 d dHðlÞ ð1 l2 Þ ¼ nðn þ 1Þ: dl HðlÞ dl
Similarly, by solving the other equations in (7), we obtain
8 u1 ðxÞ ¼ C6 x2 ; > > > > C > x4 ; < u2 ðxÞ ¼ 120 > ... > > > > : u ðxÞ ¼ ð1Þk C x2k : k ð2kþ1Þ!
ð13Þ
we get the eigenvalue equations
ð7Þ
The solution of the first equation in (7) is
A x
ð12Þ
l ¼ cosðhÞ, Eq. (12) reduces to
/ðr; lÞ ¼ RðrÞHðlÞ;
uk ð0Þ : 0:
u0 ðxÞ ¼ þ C;
Dr2 /ðr; hÞ þ ðmRf Ra Þ/ðr; hÞ ¼ 0:
ð6Þ
Substituting Eq. (5) into Eq. (4), and considering the terms with identical powers of p, a set of equations can be obtained of the form
8 0 2 00 p : x u ðxÞ þ 2xu00 ðxÞ ¼ 0; > > > 1 2 000 > > p : x u ðxÞ þ 2xu01 ðxÞ ¼ x2 u0 ðxÞ; > > < 2 2 100 p : x u2 ðxÞ þ 2xu02 ðxÞ ¼ x2 u1 ðxÞ; > > .. > > > . > > : k 2 00 p : x uk ðxÞ þ 2xu0k ðxÞ ¼ x2 uk1 ðxÞ;
Fig. 1. A bare hemisphere of radius a.
Rn ðBrÞ ¼
1 X k¼0
ð11Þ
This flux must clearly satisfy the relevant boundary conditions. 2.3. Bare hemisphere After considering the flux in the rather simple system consisting of a multiplying sphere, which is a function of r only, the bare hemisphere of radius a (Fig. 1) is now tackled. Obviously, the flux in this reactor is a function of both r and h. The time-independent diffusion equation is
An
ð1Þk Cð2n þ 2ÞCðn þ k þ 1Þ ðBrÞnþ2k : Cðn þ 1ÞCðk þ 1ÞCð2n þ 2k þ 2Þ
ð19Þ
Investigating the homotopy solution given in Eq. (19), it is straightforward to derive the recurrence relations
R2mþ1 ðBrÞ BR2mþ2 ðBrÞ ¼ þ BR2m ðBrÞ; r ð4m þ 5Þð4m þ 3Þ
ð20Þ
and
mþ1 R0 ðBrÞ Rm ðBrÞ þ m ¼ ð2m þ 1ÞRm1 ðBrÞ: Br B
ð21Þ
While the solution of the angular part (Eq. (16)) is just the Legendre polynomials
1713
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
Hn ðlÞ ¼ Pn ðlÞ ¼
n X 1 k¼0
k
d ðl2 1Þk : k d 2 k! lk
ð22Þ
with
f ðxÞ ¼
The final solution of Eq. (13) is therefore
/ðr; lÞ ¼
1 X
An Rn ðBrÞP n ðlÞ:
ð23Þ
n¼0
2.4. The radiation boundary condition Applying the two simple boundary conditions (namely ZF and EBC) to small systems yields inaccurate results (Cassell and Williams, 2004). A mixed boundary condition, namely RBC, gives more accurate results for the critical radius. In addition, it takes into consideration the interaction between fissile systems. Since the boundary conditions for the curved and flat surfaces are different, this mixed boundary condition proved to be useful as noticed by Cassell and Williams (2004) and references therein. Fig. 2 depicts two hemispheres each with radius a with the distance between their flat faces equals 2b. The RBC can be written in the form (Cassell and Williams, 2004)
n r/ ¼ gðrÞ/;
ð24Þ
where gðrÞ vary over the surface. The condition on the curved surface is
@/ðr; lÞ ¼ g 1 /ða; lÞ; @r r¼a and on the flat surface where h ¼
1 @/ðr; lÞ ¼ g 2 /ðr; 0Þ: r @ l l¼0
ð25Þ p 2
the boundary condition is
ð26Þ
We will use Marshak’s P1 boundary condition (Cassell and Williams, 2004)
1 g1 ¼ ; 2D gðb=aÞ ; g2 ¼ 2D
ð27Þ ð28Þ
where gðxÞ as given by Williams (2002) is
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R1 1 2 0 dllexp f ðxÞ 1 l2 =l gðxÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; R1 1 þ 3 0 dll2 exp f ðxÞ 1 l2 =l
x 2 pffiffiffi þ 2x : 1þx p
ð30Þ
Since this work is concerned with flux variations in a multiplying bare sphere and hemisphere, the interpretation of gðxÞ and hence g 2 is of importance. For a sphere we have b ¼ 0 and hence x ¼ 0, causing g and g 2 to vanish (Eqs. (30), (29) and (28)). For such a spherical system, the flux should not vary with the angle h (compare Eq. (26)). For a hemisphere, on the other hand, b and x go to infinity, and thus g becomes equal to unity. This corresponds to nonvanishing neutron current in both radial and angular directions (Eqs. (26) and (25)). In order to solve for /ðr; lÞ as given by Eq. (23), we now apply the radiation boundary condition. On the flat surface we have
X 1X An Rn ðBrÞP0n ð0Þ ¼ g 2 An Rn ðBrÞPn ð0Þ; r
ð31Þ
and from the properties of Legendre polynomials,
P0n ð0Þ ¼ 0;
n : even
P n ð0Þ ¼ 0;
n : odd
1n Cðn þ 1=2Þ pffiffiffiffi ; n ¼ 0; 1; 2; . . . P 2n ð0Þ ¼ pn! n 2ð1Þ Cðn þ 3=2Þ pffiffiffiffi P02nþ1 ð0Þ ¼ ; n ¼ 0; 1; 2; . . . pn! Eq. (31) reduces to,
1X 2ð1Þn Cðn þ 3=2Þ pffiffiffiffi A2nþ1 R2nþ1 ðBrÞ r pn! X 1n Cðn þ 1=2Þ pffiffiffiffi : A2n R2n ðBrÞ ¼ g2 pn!
ð32Þ
Using the recurrence relation (20), Eq. (32) can be rewritten as
ð1Þn C n þ 12 1 pffiffiffiffi 2B n þ A2nþ1 g 2 A2n 2 pn! n¼0
1 X ð1Þn1 C n þ 12 R2n pffiffiffiffi A2n1 ð2nBÞ þ ð4n þ 1Þð4n 1Þ pn! n¼0
1 X
R2n
¼ 0:
ð33Þ
The amplitudes An are related to each other by
ð29Þ
Bð2n þ 1ÞA2nþ1
2nB A2n1 ¼ g 2 A2n ; ð4n þ 1Þð4n 1Þ
ð34Þ
which results from the simplification of Eq. (33). From Eq. (34) we can write the odd order amplitudes A2nþ1 in terms of the even order ones A2n as
A2nþ1 ¼
n X
cnk A2k ;
ð35Þ
k¼0
where cnk
cnk ¼
ð2nÞ!! ð2k 1Þ!!ð4k þ 1Þ!! g 2 : ð2kÞ!! ð2n þ 1Þ!!ð4n þ 1Þ!! B
ð36Þ
Similarly, and by applying the RBC on the curved surface in order to solve for /ðr; lÞ, we get 1 X
An R0n ðBaÞPn ðlÞ ¼ g 1
n¼0
1 X
An Rn ðBaÞP n ðlÞ;
ð37Þ
n¼0
or 1 X
Fig. 2. Two hemispheres with separation distance 2b.
n¼0
An fn Pn ðlÞ ¼ 0;
ð38Þ
1714
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
where
fn ¼
R0n ðBaÞ
þ g 1 Rn ðBaÞ:
ð39Þ
Using the recurrence relations (20) and (21), the parameters fn as given by Eq. (39) can be rewritten in the form,
fn ¼
n B þ g 1 Rn ðBaÞ Rnþ1 ðBaÞ: a ð2n þ 3Þ
ð40Þ
In order to make the computational implementation of these procedures more convenient, it is favorable to reshape Eq. (38). This can be achieved by introducing the so called shifted Legendre polynomials Pn ðlÞ, where
Pn ðlÞ ¼ Pn ð2l 1Þ;
ð41Þ
which form a complete orthogonal set over 0 < l < 1, while classical Legendre polynomials Pn ðlÞ form an orthogonal set over 1 < l < 1. Then n X
Pn ðlÞ ¼
bnk Pk ðlÞ;
ð42Þ
k¼0
where
bnk ¼ ð2k þ 1Þ
Z
1
dlPm ðlÞPk ð2l 1Þ:
ð43Þ
0
functions to the results of Cassell and Williams (2004), who used the same nuclear data with the application of Bessel functions only, and with a different computational approach. 3.1. Critical radius and fuel calculations For the rather simple case of the bare sphere, a computer code based on Eq. (11) has been written to calculate, by iteration, the critical value of Ba which was found to be 3.1416. The flux converges to zero at this value of Bac . If the ZF boundary condition is applied, the critical radius ac;ZF is 10.528 cm. On the other hand, and since the flux is assumed to converge to zero at ac þ 2D using the EBC, then ac;EBC is equal to 8.656 cm. The hemisphere example, on the other hand, incorporates interesting features. The general solution given by Eq. (23) can be simplified when applying the ZF boundary condition. For the flux to vanish on the flat surface ðl ¼ 0Þ, the properties of the Legendre polynomials imply that the even amplitudes must be equal to zero. In addition, the ZF boundary condition on the curved surface implies the need to find the first zero of any single homotopy-based solution Rn ðBrÞ of the radial part (Eq. (19)), or a combination of these eigenfunctions. By investigating the properties of these functions, the amplitudes corresponding to all eigenvalues n > 1 must also vanish. This reduces the summation in Eq. (23) to the simple form
Therefore, Eq. (38) can be rewritten as 1 X
An fn
n¼0
n X
/ðr; lÞ ¼ A1 R1 ðBrÞP 1 ðlÞ ¼ bnk P k ð
lÞ ¼ 0:
Interchanging the order of summation, and noting that bnk ¼ 0 for n < k, we finally find 1 1 X X k¼0
! bnk fn An Pk ðlÞ ¼ 0:
ð45Þ
n¼k
From Eq. (45), we conclude that 1 X
bnk fn An ¼ 0;
k ¼ 0; 1; 2; 3; . . .
ð46Þ
n¼0
The even and odd terms can be separated, yielding 1 X
½b2n;k f2n A2n þ b2nþ1;k f2nþ1 A2nþ1 ¼ 0:
ð47Þ
n¼0
The odd amplitudes A2nþ1 can be written in terms of the even amplitudes A2n by making use of Eq. (35), this leads to 1 X
b nk A2n ¼ 0; M
k ¼ 0; 1; 2; . . . ;
ð48Þ
n¼0
A1
ð1Þk ð6Þðk þ 1Þ ðBrÞ1þ2k l: Cð2k þ 4Þ
ð50Þ
Therefore, for 1 MeV neutrons diffusing in a hemisphere of pure U, and using the same multiplying medium parameters used for the bare sphere, another iteration code based on Eq. (50) yields Bac ¼ 4:4934; ac;ZF ¼ 15:085 cm and ac;EBC ¼ 13:185 cm. The results obtained above for the critical radius, in the case of a sphere and hemisphere, and for both ZF and EBC boundary conditions, are listed in Table 1. Cassell and Williams (2004) reported similar results using Bessel functions instead of the homotopy constructed in this work. Since it is expected that the ZF and EBC boundary conditions do not result in accurate calculation of the critical radius, especially for small systems, the RBC has been also applied. In order to get numerical results, the infinite summation in Eqs. (48) and (49), and correspondingly the infinite dimensional matrix, will be limited by a cut-off value N, thus
235
N X
b nk A2n ¼ 0; M
k ¼ 0; 1; 2; . . . ; N;
ð51Þ
n¼0
where
b nk ¼ b f2n þ M 2n;k
where
b nk ¼ b f2n þ M 2n;k
k¼0
ð44Þ
k¼0
1 X
N X
b2lþ1;k f2lþ1 cl;n :
ð52Þ
l¼n 1 X
b2lþ1;k f2lþ1 cl;n :
ð49Þ
l¼n
Consequently, and in view of Eq. (35), the cut-off at N implies that the series in the general solution (Eq. (23)) will be also limited to 2N þ 2 terms, thus
3. Results In order to provide numerical results, 1 MeV neutrons diffusing in pure 235U are considered. The following data will correspondingly be used:
m ¼ 2:42; rf ¼ 1:336b; rs ¼ 5:959b; rc ¼ 0:153b;
Table 1 Critical radius ac (in cm) for a sphere and hemisphere using different boundary conditions.
Nc ¼ 0:0478 1024 atoms cm3 : The diffusion coefficient D ¼ 1=ð3Nc rtotal Þ is then equal to 0.9363 cm. We adopt these data in what follows in order to compare the results obtained in this work using both HPM and Bessel
a
BC
Sphere
% Errora
Hemisphere
% Errora
ZF EBC RBC
10.528 8.656 8.441
24.7 2.5 –
15.085 13.185 11.804
27.8 11.7 –
Compared with RBC.
1715
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
/ðr; lÞ ¼
2Nþ1 X
An Rn ðBrÞPn ðlÞ:
ð53Þ
n¼0
In order to satisfy Eq. (51) for any vector of amplitudes A2n , the b nk should be determinant of the ðN þ 1Þ ðN þ 1Þ square matrix M equal to zero. Therefore, and in order to calculate the critical radius using the RBC, a third computer code has been especially written b nk for an for this work. The code calculates the determinant of M increasing radius a. The value of a affects the determinant through the parameters fn (Eq. (40)). The values of the determinant are either positive or negative, with a reflection point corresponding to the critical radius ac . We will refer to this computational methodology as the ‘‘reflection point technique” in what follows. On the other hand, the calculated ac depends on the value of the interaction parameter g through cnk (Eqs. (36) and (28)). Hence, the results for the sphere and the infinitely separated hemispheres are obtained by using g ¼ 0 and g ¼ 1, respectively. We aim at comparing the results obtained in this work with those reported by Cassell and Williams (2004). Therefore, we chose to use N ¼ 4, i.e. 5 5 matrix and 10 terms in the expansion of Eq. (53). The critical radius for different systems (different values of the interaction parameter g) are listed in Table 2. Tabulated also are the values of the fuel mass per hemisphere. It is evident that ac and the corresponding mass of the fuel increase with the interaction parameter g, to compensate for the increasing leakage. The results in Table 2 are similar to those reported by Cassell and Williams (2004), indicating that HPM can very well reproduce the corresponding values obtained by using Bessel functions. The results, on the other hand, depend not only on the choice between HPM and Bessel functions, but on the implemented computational method as well. Therefore, the critical radius of the hemisphere ðg ¼ 1Þ, accurate to seven significant figures, is calculated using variable number of terms and listed in Table 3. It is worth mentioning here, that mainly for purposes of flux calculations, Cassell and Williams (2004) proposed an alternative for-
mula for bnk in order to ensure that the boundary condition is satisfied at discrete points. This formula is
bnk ¼ Pn ðlk Þ;
ð54Þ
with
lk ¼ cosðkp=2NÞ:
ð55Þ
This alternative bnk (Method 2 as referred by Cassell and Williams (2004) is used in Eq. (52) instead of Eq. (43) (Method 1). Since both forms of bnk are independent of either Bessel functions or HPM, it has neutral effect on the overall comparison between the results of this work and those reported by Cassell and Williams (2004). As described earlier in this subsection, the critical radius is calculated in this work by the reflection point technique. We
Table 4 Angular distribution over the curved surface using Method 2 and 62 terms (N = 30). Theta
HPM
Bessel functions
0 10 20 30 40 50 60 70 80 90
0.6943497 0.6904001 0.6688981 0.6336942 0.5846480 0.5219508 0.4459712 0.3564137 0.2520653 0.1282844
0.6943497 0.6904001 0.6688981 0.6336942 0.5846480 0.5219508 0.4459712 0.3564137 0.2520653 0.1282844
Table 2 Variation of the critical radius ac and critical mass M c per hemisphere with the interaction parameter g. g
ac ðcmÞ
Mc ðkgÞ
0 (sphere) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
8.44123 9.14444 9.73615 10.21877 10.60685 10.91835 11.16983 11.37485 11.54391 11.68495 11.80396
23.431 29.788 35.953 41.569 46.487 50.704 54.289 57.334 59.928 62.152 64.070
Fig. 3. Angular distribution over the curved surface.
Table 3 Critical radius for the hemisphere using two different methods. Terms
4 6 8 10 12 22 32 42 52
Method 1
Method 2
HPM this work
Bessel this work
Bessel Cassell and Williams (2004)
HPM this work
Bessel this work
Bessel Cassell and Williams (2004)
11.80012 11.80389 11.80395 11.80396 11.80396 11.80396 11.80396 11.80396 –
11.80012 11.80389 11.80395 11.80396 11.80396 11.80396 11.80396 – –
11.80008 11.80385 11.80391 11.80392 11.80392 11.80392 11.80392 11.80392 –
12.04073 11.79314 11.79867 11.80112 11.80234 11.80359 11.80380 11.80387 –
12.04073 11.79314 11.79867 11.80112 11.80234 11.80359 11.80380 – –
12.04069 11.79310 11.79863 11.80108 11.80230 11.80355 11.80377 11.80383 11.80386
1716
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
calculated the critical radius as a function of the number of terms, using both Bessel functions as well as HPM. The calculations are repeated twice, adopting Method 1 and Method 2 described above. The results are listed in Table 3. The calculated critical radius does converge using Method 1, while Method 2, which is proposed mainly for purposes of flux calculations, does not show any similar convergence, whatever number of terms are used. Method 1 in Table 3 reveals the fact that the calculation of the critical radius of the hemisphere, to accuracy of seven significant figures, converges with 10 terms using the reflection point technique, exactly similar to that reported by Cassell and Williams (2004). This convergence is independent of whether we are using HPM or Bessel functions with the reflection point technique. The small difference in the convergence value obtained in this work compared to Cassell and Williams (2004), however, is related to the difference in the implemented computational method.
3.2. Flux distribution It is worthwhile now to adopt Method 2 and examine the flux distribution inside the hemisphere. For this purpose, and in order to establish the computational grounds, we will write the flux given by Eq. (53) as
/ðr; lÞ ¼
N X ½A2n R2n ðBrÞP2n ðlÞ þ A2nþ1 R2nþ1 ðBrÞP2nþ1 ðlÞ:
ð56Þ
n¼0
Furthermore, and using Eq. (35), we finally arrive at the form that has been used in the code to compute the flux in the bare hemisphere, that is
/ðr; lÞ ¼
N X
" A2n R2n ðBrÞP2n ðlÞ þ
n¼0
N X
#
cln R2lþ1 ðBrÞP2lþ1 ðlÞ :
ð57Þ
l¼n
In order to calculate the amplitudes A2n , Eq. (51) has been rearranged as N X
b nk A2n ¼ M b 0k A0 ; M
ð58Þ
n¼1
and A1 has been normalized to unity. Hence, A0 ¼ c1 . The subse00 quent flux calculations can thus be scaled to absolute numbers by adjustment of the amplitudes. In addition, all the calculated flux values are normalized to the volume averaged flux /
¼ /
Ra R1 0
r 2 /ðr; lÞdldr 0 ; Ra R1 r 2 dldr 0 0
ð59Þ
the azimuth angle being neutral. Hence
¼ 3 / a3
Fig. 4. Flux distribution across the core of a bare hemisphere at different angles.
Z
a
dr r 2
0
Z
1
dl/ðr; lÞ:
ð60Þ
0
The value of a in Eq. (60) is 11.80396 cm, the convergence value of the critical radius (Table 3). For example, using 46 terms ðN ¼ 22Þ provides values for the volume averaged flux of
Table 5 Flux variation across the critical hemisphere at h ¼ p=4. r=a
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
4 Terms
12 Terms Bessel functions
HPM
Bessel functions
HPM
Bessel functions
0.887119 1.25636 1.55550 1.76260 1.86340 1.85253 1.73389 1.52015 1.23146 0.893463 0.534903
0.887119 1.25636 1.55550 1.76260 1.86340 1.85253 1.73389 1.52015 1.23146 0.893463 0.534903
0.881442 1.24814 1.54501 1.75055 1.85094 1.84118 1.72544 1.51647 1.23433 0.904402 0.554956
0.881442 1.24814 1.54501 1.75055 1.85094 1.84118 1.72544 1.51647 1.23433 0.904402 0.554956
0.881631 1.24841 1.54535 1.75095 1.85137 1.84163 1.72586 1.51684 1.23462 0.90457 0.554941
0.881631 1.24841 1.54535 1.75095 1.85137 1.84163 1.72586 1.51684 1.23462 0.90457 0.554941
46 Terms 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
22 Terms
HPM
0.881673 1.24847 1.54543 1.75104 1.85147 1.84172 1.72595 1.51692 1.23469 0.904621 0.554972
62 Terms 0.881673 1.24847 1.54543 1.75104 1.85147 1.84172 1.72595 1.51692 1.23469 0.904621 0.554972
0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72595 1.51692 1.23469 0.904622 0.554973
86 Terms 0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72595 1.51692 1.23469 0.904622 0.554973
0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72596 1.51692 1.23469 0.904622 0.554974
0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72596 1.51692 1.23469 0.904622 0.554974
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
1717
ism developed in this work and illustrated by Eqs. (57) and (58), it is worthwhile to compare our results to those obtained using deterministic finite element approximation to the neutron transport equation, which is implemented in the EVENT code, and reported in Cassell and Williams (2004). Fig. 5 illustrates that, though HPM can exactly reproduce the results obtained using Bessel functions, the computational implementation illustrated by Eqs. (57) and (58) yields results closer to those obtained using EVENT, especially when smaller number of terms in the expansion are used. 4. Conclusions The homotopy perturbation method, applied to the neutron diffusion problem in hemispherical symmetry, can exactly reproduce the results obtained using Bessel functions. On the other hand, the computational implementation of the analytical solutions obtained in this work does affect the results, and improve them when compared to those obtained using finite element codes like EVENT. Fig. 5. Percent error in the normalized flux values calculated relative to EVENT. Closed and open symbols correspond to results reported by Cassell and Williams (2004) and in this work, respectively.
0.092076 and 0.633771, using Bessel functions and HPM, respectively. This clearly indicates that scaling to absolute flux values differs accordingly. The flux behavior on the curved surface of the critical hemisphere is inferred from Table 4, for both HPM and Bessel functions with 62 terms in Eq. (56) ðN ¼ 30Þ. The remarkable agreement between both calculations is evident, taking into consideration the fact that the angular dependence itself is determined by the Legendre polynomials. Fig. 3 illustrates the fast convergence of the flux calculations on the curved surface with the number of terms. The critical flux across the core has been also calculated using HPM, at different angles, and is shown in Fig. 4. It is evident that at the flat surface ðh ¼ 90 Þ the flux decreases as we move towards the curved surface. Unlike the sphere, the maximum flux (namely 2.29786) is not at r ¼ 0 but at ðr; hÞ ¼ ð5:4 cm;0 Þ. The value of the flux at r ¼ 0 is 0.881671, which is significant compared with the maximum value, confirming the fact that the ZF boundary condition / ¼ 0 at the flat surface is not valid for such a small hemisphere. In order to compare flux calculations using HPM to those using Bessel functions (both calculated in this work), data in Table 5 provide normalized flux values across the hemisphere for h ¼ 45 using Method 2. It is again evident that HPM can perfectly reproduce the normalized flux values calculated using Bessel functions, establishing HPM as a new technique that can be reliably used for similar calculations. Furthermore, recalling that the flux has been computed, both using HPM and Bessel functions, using the formal-
Acknowledgment The financial support provided by King Abdullah II Fund for Development KAFD is gratefully acknowledged by K. Kh. References 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. 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., 2006. Homotopy perturbation method for solving boundary value problems. Physics Letters A 350, 87–88. He, J.H., 2006. Some asymptotic methods for strongly nonlinear equations. International Journal of Modern Physics B 20, 1141–1199. 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. Williams, M.M.R., 2002. Neutron transport in two interacting slabs. II: the influence of transverse leakage. Annals of Nuclear Energy 29, 1551–1596.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A solution of the neutron diffusion equation in hemispherical symmetry using the homotopy perturbation method Kafa Khasawneh a, Saed Dababneh a,b,*, Zaid Odibat a a b
Prince Abdullah Bin Ghazi Faculty of Science and IT, Al-Balqa Applied University, P.O. Box 7051, 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 June 2009 Received in revised form 8 September 2009 Accepted 9 September 2009 Available online 7 October 2009
a b s t r a c t The homotopy perturbation method is used to formulate a new analytic solution of the neutron diffusion equation both for a sphere and a hemisphere of fissile material. Different boundary conditions are investigated; including zero flux on boundary, zero flux on extrapolated boundary, and radiation boundary condition. The interaction between two hemispheres with opposite flat faces is also presented. Numerical results are provided for one-speed fast neutrons in 235U. A comparison with Bessel function based solutions demonstrates that the homotopy perturbation method can exactly reproduce the results. The computational implementation of the analytic solutions was found to improve the numeric results when compared to finite element calculations. Ó 2009 Published by Elsevier Ltd.
1. Introduction
2. Theory
There has been variety of methods used over the last decades to solve the one-group neutron diffusion equation. Up to our knowledge, the nonlinear analytical homotopy perturbation method (HPM) has never been tested as a means to solve reactor physics problems. The merit of this contribution is that it provides a description for the application of HPM to solve time-independent neutron diffusion equations, for different geometries at different boundary conditions. These boundary conditions are the zero flux at the boundary ZF, the zero flux at the extrapolated boundary EBC, and the radiation boundary condition RBC. The numerical results obtained in this work are compared to the outcome of other Bessel function based solutions, both obtained in this work and reported by Cassell and Williams (2004), in addition to transport theory calculations reported in the same reference. The aim of this contribution is to improve, or at least widen the methodologies used to investigate the flux behavior inside the reactor core. The general theory including the essential HPM derivations, and its application to the sphere and hemisphere are given in Section 2. The RBC boundary condition is also detailed in order to set the grounds for the calculations. Numerical results are provided for critical radius calculations, as well as for flux distribution in the critical hemisphere (Section 3). The accompanying technical issues faced during the computations are described together with their corresponding remedies.
2.1. The homotopy perturbation method (HPM)
* Corresponding author. Address: Prince Abdullah Bin Ghazi Faculty of Science and IT, Al-Balqa Applied University, P.O. Box 7051, Salt 19117, Jordan. Tel.: +962 7 95606613; fax: +962 6 5341608. E-mail addresses: [email protected], [email protected] (S. Dababneh). 0306-4549/$ - see front matter Ó 2009 Published by Elsevier Ltd. doi:10.1016/j.anucene.2009.09.001
The HPM, proposed first by He (1999, 2000) for solving linear and nonlinear differential and integral equations, has been the subject of extensive analytical and numerical studies. The method, which is a coupling of a homotopy technique and a perturbation technique, deforms continuously to a simple problem which is easily solved. This method, which does not require a small parameter in an equation, nor it incorporates unjustified assumptions, 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, 2006). Analysis of the method and different illustrative examples (He, 1999, 2000, 2003, 2006; Odibat, 2007) present HPM as a promissing means to solve general differential equations. The HPM is applied to neutron diffusion in a bare sphere and hemisphere, in Sections 2.2 and 2.3, respectively. 2.2. Bare sphere We will consider first the rather simple example of a bare spherical reactor of radius a. The spherical symmetry of the system implies that the flux is a function of r only. The time-independent diffusion equation can thus be written as
1712
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
r2 /ðrÞ þ B2 /ðrÞ ¼ 0;
ð1Þ
where
B2 ¼
mRf Ra
ð2Þ
D
is the buckling of the reactor. Substituting the Laplacian in spherical coordinates and using x ¼ Br, we get 2
x2
d /ðxÞ 2
dx
þ 2x
d/ðxÞ þ x2 /ðxÞ ¼ 0: dx
ð3Þ
In order to solve Eq. (3) using HPM, we construct the following homotopy:
Hðu; pÞ ¼ x2 u00 ðxÞ þ 2xu0 ðxÞ þ px2 uðxÞ;
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 x2 u00 ðxÞ þ 2xu0 ðxÞ ¼ 0, to uðxÞ; 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, 2
uðxÞ ¼ u0 ðxÞ þ pu1 ðxÞ þ p u2 ðxÞ þ
ð5Þ
Therefore, the solution of Eq. (3) can be readily obtained as
uðxÞ ¼ lim½u0 ðxÞ þ pu1 ðxÞ þ p2 u2 ðxÞ þ : p!1
u0 ð0Þ : finite; u1 ð0Þ : 0; u2 ð0Þ : 0;
ð8Þ
where A and C are constants. The solution should be finite at r = zero; thus u0 ðxÞ reduces to
u0 ðxÞ ¼ C:
Writing the Laplacian in spherical coordinates and using
@ 2 /ðr; lÞ 2 @/ðr; lÞ 1 @ 2 @/ðr; lÞ þ ð1 l Þ þ @r 2 r @r r2 @ l @l þ B2 /ðr; lÞ ¼ 0;
where the buckling B2 is again given by Eq. (2). Applying the separation of variables,
ð9Þ
2
ð15Þ ð16Þ
In order to solve the radial part using HPM, we first consider x ¼ Br and rewrite Eq. (15) as 2
x2
d RðxÞ 2
dx
þ 2x
dRðxÞ þ ðx2 nðn þ 1ÞÞRðxÞ ¼ 0: dx
ð17Þ
Then we construct the homotopy 2
Hðu; pÞ ¼ x2
d uðxÞ 2
dx
þ 2x
duðxÞ þ ðpx2 nðn þ 1ÞÞuðxÞ; dx
p 2 ½0; 1: ð18Þ
ð10Þ
Following the same manner as in Section 2.2, if we substitute Eq. (5) into the homotopy (18), and solve the resulting equations, we get the following solution:
Consequently, the final solution of Eq. (1), in a closed form, is given by 1 X ð1Þk C /ðBrÞ ¼ ðBrÞ2k : ð2k þ 1Þ! k¼0
ð14Þ
r 2 d RðrÞ 2r dRðrÞ þ þ B2 r 2 ¼ nðn þ 1Þ; RðrÞ dr RðrÞ dr 2 1 d dHðlÞ ð1 l2 Þ ¼ nðn þ 1Þ: dl HðlÞ dl
Similarly, by solving the other equations in (7), we obtain
8 u1 ðxÞ ¼ C6 x2 ; > > > > C > x4 ; < u2 ðxÞ ¼ 120 > ... > > > > : u ðxÞ ¼ ð1Þk C x2k : k ð2kþ1Þ!
ð13Þ
we get the eigenvalue equations
ð7Þ
The solution of the first equation in (7) is
A x
ð12Þ
l ¼ cosðhÞ, Eq. (12) reduces to
/ðr; lÞ ¼ RðrÞHðlÞ;
uk ð0Þ : 0:
u0 ðxÞ ¼ þ C;
Dr2 /ðr; hÞ þ ðmRf Ra Þ/ðr; hÞ ¼ 0:
ð6Þ
Substituting Eq. (5) into Eq. (4), and considering the terms with identical powers of p, a set of equations can be obtained of the form
8 0 2 00 p : x u ðxÞ þ 2xu00 ðxÞ ¼ 0; > > > 1 2 000 > > p : x u ðxÞ þ 2xu01 ðxÞ ¼ x2 u0 ðxÞ; > > < 2 2 100 p : x u2 ðxÞ þ 2xu02 ðxÞ ¼ x2 u1 ðxÞ; > > .. > > > . > > : k 2 00 p : x uk ðxÞ þ 2xu0k ðxÞ ¼ x2 uk1 ðxÞ;
Fig. 1. A bare hemisphere of radius a.
Rn ðBrÞ ¼
1 X k¼0
ð11Þ
This flux must clearly satisfy the relevant boundary conditions. 2.3. Bare hemisphere After considering the flux in the rather simple system consisting of a multiplying sphere, which is a function of r only, the bare hemisphere of radius a (Fig. 1) is now tackled. Obviously, the flux in this reactor is a function of both r and h. The time-independent diffusion equation is
An
ð1Þk Cð2n þ 2ÞCðn þ k þ 1Þ ðBrÞnþ2k : Cðn þ 1ÞCðk þ 1ÞCð2n þ 2k þ 2Þ
ð19Þ
Investigating the homotopy solution given in Eq. (19), it is straightforward to derive the recurrence relations
R2mþ1 ðBrÞ BR2mþ2 ðBrÞ ¼ þ BR2m ðBrÞ; r ð4m þ 5Þð4m þ 3Þ
ð20Þ
and
mþ1 R0 ðBrÞ Rm ðBrÞ þ m ¼ ð2m þ 1ÞRm1 ðBrÞ: Br B
ð21Þ
While the solution of the angular part (Eq. (16)) is just the Legendre polynomials
1713
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
Hn ðlÞ ¼ Pn ðlÞ ¼
n X 1 k¼0
k
d ðl2 1Þk : k d 2 k! lk
ð22Þ
with
f ðxÞ ¼
The final solution of Eq. (13) is therefore
/ðr; lÞ ¼
1 X
An Rn ðBrÞP n ðlÞ:
ð23Þ
n¼0
2.4. The radiation boundary condition Applying the two simple boundary conditions (namely ZF and EBC) to small systems yields inaccurate results (Cassell and Williams, 2004). A mixed boundary condition, namely RBC, gives more accurate results for the critical radius. In addition, it takes into consideration the interaction between fissile systems. Since the boundary conditions for the curved and flat surfaces are different, this mixed boundary condition proved to be useful as noticed by Cassell and Williams (2004) and references therein. Fig. 2 depicts two hemispheres each with radius a with the distance between their flat faces equals 2b. The RBC can be written in the form (Cassell and Williams, 2004)
n r/ ¼ gðrÞ/;
ð24Þ
where gðrÞ vary over the surface. The condition on the curved surface is
@/ðr; lÞ ¼ g 1 /ða; lÞ; @r r¼a and on the flat surface where h ¼
1 @/ðr; lÞ ¼ g 2 /ðr; 0Þ: r @ l l¼0
ð25Þ p 2
the boundary condition is
ð26Þ
We will use Marshak’s P1 boundary condition (Cassell and Williams, 2004)
1 g1 ¼ ; 2D gðb=aÞ ; g2 ¼ 2D
ð27Þ ð28Þ
where gðxÞ as given by Williams (2002) is
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R1 1 2 0 dllexp f ðxÞ 1 l2 =l gðxÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; R1 1 þ 3 0 dll2 exp f ðxÞ 1 l2 =l
x 2 pffiffiffi þ 2x : 1þx p
ð30Þ
Since this work is concerned with flux variations in a multiplying bare sphere and hemisphere, the interpretation of gðxÞ and hence g 2 is of importance. For a sphere we have b ¼ 0 and hence x ¼ 0, causing g and g 2 to vanish (Eqs. (30), (29) and (28)). For such a spherical system, the flux should not vary with the angle h (compare Eq. (26)). For a hemisphere, on the other hand, b and x go to infinity, and thus g becomes equal to unity. This corresponds to nonvanishing neutron current in both radial and angular directions (Eqs. (26) and (25)). In order to solve for /ðr; lÞ as given by Eq. (23), we now apply the radiation boundary condition. On the flat surface we have
X 1X An Rn ðBrÞP0n ð0Þ ¼ g 2 An Rn ðBrÞPn ð0Þ; r
ð31Þ
and from the properties of Legendre polynomials,
P0n ð0Þ ¼ 0;
n : even
P n ð0Þ ¼ 0;
n : odd
1n Cðn þ 1=2Þ pffiffiffiffi ; n ¼ 0; 1; 2; . . . P 2n ð0Þ ¼ pn! n 2ð1Þ Cðn þ 3=2Þ pffiffiffiffi P02nþ1 ð0Þ ¼ ; n ¼ 0; 1; 2; . . . pn! Eq. (31) reduces to,
1X 2ð1Þn Cðn þ 3=2Þ pffiffiffiffi A2nþ1 R2nþ1 ðBrÞ r pn! X 1n Cðn þ 1=2Þ pffiffiffiffi : A2n R2n ðBrÞ ¼ g2 pn!
ð32Þ
Using the recurrence relation (20), Eq. (32) can be rewritten as
ð1Þn C n þ 12 1 pffiffiffiffi 2B n þ A2nþ1 g 2 A2n 2 pn! n¼0
1 X ð1Þn1 C n þ 12 R2n pffiffiffiffi A2n1 ð2nBÞ þ ð4n þ 1Þð4n 1Þ pn! n¼0
1 X
R2n
¼ 0:
ð33Þ
The amplitudes An are related to each other by
ð29Þ
Bð2n þ 1ÞA2nþ1
2nB A2n1 ¼ g 2 A2n ; ð4n þ 1Þð4n 1Þ
ð34Þ
which results from the simplification of Eq. (33). From Eq. (34) we can write the odd order amplitudes A2nþ1 in terms of the even order ones A2n as
A2nþ1 ¼
n X
cnk A2k ;
ð35Þ
k¼0
where cnk
cnk ¼
ð2nÞ!! ð2k 1Þ!!ð4k þ 1Þ!! g 2 : ð2kÞ!! ð2n þ 1Þ!!ð4n þ 1Þ!! B
ð36Þ
Similarly, and by applying the RBC on the curved surface in order to solve for /ðr; lÞ, we get 1 X
An R0n ðBaÞPn ðlÞ ¼ g 1
n¼0
1 X
An Rn ðBaÞP n ðlÞ;
ð37Þ
n¼0
or 1 X
Fig. 2. Two hemispheres with separation distance 2b.
n¼0
An fn Pn ðlÞ ¼ 0;
ð38Þ
1714
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
where
fn ¼
R0n ðBaÞ
þ g 1 Rn ðBaÞ:
ð39Þ
Using the recurrence relations (20) and (21), the parameters fn as given by Eq. (39) can be rewritten in the form,
fn ¼
n B þ g 1 Rn ðBaÞ Rnþ1 ðBaÞ: a ð2n þ 3Þ
ð40Þ
In order to make the computational implementation of these procedures more convenient, it is favorable to reshape Eq. (38). This can be achieved by introducing the so called shifted Legendre polynomials Pn ðlÞ, where
Pn ðlÞ ¼ Pn ð2l 1Þ;
ð41Þ
which form a complete orthogonal set over 0 < l < 1, while classical Legendre polynomials Pn ðlÞ form an orthogonal set over 1 < l < 1. Then n X
Pn ðlÞ ¼
bnk Pk ðlÞ;
ð42Þ
k¼0
where
bnk ¼ ð2k þ 1Þ
Z
1
dlPm ðlÞPk ð2l 1Þ:
ð43Þ
0
functions to the results of Cassell and Williams (2004), who used the same nuclear data with the application of Bessel functions only, and with a different computational approach. 3.1. Critical radius and fuel calculations For the rather simple case of the bare sphere, a computer code based on Eq. (11) has been written to calculate, by iteration, the critical value of Ba which was found to be 3.1416. The flux converges to zero at this value of Bac . If the ZF boundary condition is applied, the critical radius ac;ZF is 10.528 cm. On the other hand, and since the flux is assumed to converge to zero at ac þ 2D using the EBC, then ac;EBC is equal to 8.656 cm. The hemisphere example, on the other hand, incorporates interesting features. The general solution given by Eq. (23) can be simplified when applying the ZF boundary condition. For the flux to vanish on the flat surface ðl ¼ 0Þ, the properties of the Legendre polynomials imply that the even amplitudes must be equal to zero. In addition, the ZF boundary condition on the curved surface implies the need to find the first zero of any single homotopy-based solution Rn ðBrÞ of the radial part (Eq. (19)), or a combination of these eigenfunctions. By investigating the properties of these functions, the amplitudes corresponding to all eigenvalues n > 1 must also vanish. This reduces the summation in Eq. (23) to the simple form
Therefore, Eq. (38) can be rewritten as 1 X
An fn
n¼0
n X
/ðr; lÞ ¼ A1 R1 ðBrÞP 1 ðlÞ ¼ bnk P k ð
lÞ ¼ 0:
Interchanging the order of summation, and noting that bnk ¼ 0 for n < k, we finally find 1 1 X X k¼0
! bnk fn An Pk ðlÞ ¼ 0:
ð45Þ
n¼k
From Eq. (45), we conclude that 1 X
bnk fn An ¼ 0;
k ¼ 0; 1; 2; 3; . . .
ð46Þ
n¼0
The even and odd terms can be separated, yielding 1 X
½b2n;k f2n A2n þ b2nþ1;k f2nþ1 A2nþ1 ¼ 0:
ð47Þ
n¼0
The odd amplitudes A2nþ1 can be written in terms of the even amplitudes A2n by making use of Eq. (35), this leads to 1 X
b nk A2n ¼ 0; M
k ¼ 0; 1; 2; . . . ;
ð48Þ
n¼0
A1
ð1Þk ð6Þðk þ 1Þ ðBrÞ1þ2k l: Cð2k þ 4Þ
ð50Þ
Therefore, for 1 MeV neutrons diffusing in a hemisphere of pure U, and using the same multiplying medium parameters used for the bare sphere, another iteration code based on Eq. (50) yields Bac ¼ 4:4934; ac;ZF ¼ 15:085 cm and ac;EBC ¼ 13:185 cm. The results obtained above for the critical radius, in the case of a sphere and hemisphere, and for both ZF and EBC boundary conditions, are listed in Table 1. Cassell and Williams (2004) reported similar results using Bessel functions instead of the homotopy constructed in this work. Since it is expected that the ZF and EBC boundary conditions do not result in accurate calculation of the critical radius, especially for small systems, the RBC has been also applied. In order to get numerical results, the infinite summation in Eqs. (48) and (49), and correspondingly the infinite dimensional matrix, will be limited by a cut-off value N, thus
235
N X
b nk A2n ¼ 0; M
k ¼ 0; 1; 2; . . . ; N;
ð51Þ
n¼0
where
b nk ¼ b f2n þ M 2n;k
where
b nk ¼ b f2n þ M 2n;k
k¼0
ð44Þ
k¼0
1 X
N X
b2lþ1;k f2lþ1 cl;n :
ð52Þ
l¼n 1 X
b2lþ1;k f2lþ1 cl;n :
ð49Þ
l¼n
Consequently, and in view of Eq. (35), the cut-off at N implies that the series in the general solution (Eq. (23)) will be also limited to 2N þ 2 terms, thus
3. Results In order to provide numerical results, 1 MeV neutrons diffusing in pure 235U are considered. The following data will correspondingly be used:
m ¼ 2:42; rf ¼ 1:336b; rs ¼ 5:959b; rc ¼ 0:153b;
Table 1 Critical radius ac (in cm) for a sphere and hemisphere using different boundary conditions.
Nc ¼ 0:0478 1024 atoms cm3 : The diffusion coefficient D ¼ 1=ð3Nc rtotal Þ is then equal to 0.9363 cm. We adopt these data in what follows in order to compare the results obtained in this work using both HPM and Bessel
a
BC
Sphere
% Errora
Hemisphere
% Errora
ZF EBC RBC
10.528 8.656 8.441
24.7 2.5 –
15.085 13.185 11.804
27.8 11.7 –
Compared with RBC.
1715
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
/ðr; lÞ ¼
2Nþ1 X
An Rn ðBrÞPn ðlÞ:
ð53Þ
n¼0
In order to satisfy Eq. (51) for any vector of amplitudes A2n , the b nk should be determinant of the ðN þ 1Þ ðN þ 1Þ square matrix M equal to zero. Therefore, and in order to calculate the critical radius using the RBC, a third computer code has been especially written b nk for an for this work. The code calculates the determinant of M increasing radius a. The value of a affects the determinant through the parameters fn (Eq. (40)). The values of the determinant are either positive or negative, with a reflection point corresponding to the critical radius ac . We will refer to this computational methodology as the ‘‘reflection point technique” in what follows. On the other hand, the calculated ac depends on the value of the interaction parameter g through cnk (Eqs. (36) and (28)). Hence, the results for the sphere and the infinitely separated hemispheres are obtained by using g ¼ 0 and g ¼ 1, respectively. We aim at comparing the results obtained in this work with those reported by Cassell and Williams (2004). Therefore, we chose to use N ¼ 4, i.e. 5 5 matrix and 10 terms in the expansion of Eq. (53). The critical radius for different systems (different values of the interaction parameter g) are listed in Table 2. Tabulated also are the values of the fuel mass per hemisphere. It is evident that ac and the corresponding mass of the fuel increase with the interaction parameter g, to compensate for the increasing leakage. The results in Table 2 are similar to those reported by Cassell and Williams (2004), indicating that HPM can very well reproduce the corresponding values obtained by using Bessel functions. The results, on the other hand, depend not only on the choice between HPM and Bessel functions, but on the implemented computational method as well. Therefore, the critical radius of the hemisphere ðg ¼ 1Þ, accurate to seven significant figures, is calculated using variable number of terms and listed in Table 3. It is worth mentioning here, that mainly for purposes of flux calculations, Cassell and Williams (2004) proposed an alternative for-
mula for bnk in order to ensure that the boundary condition is satisfied at discrete points. This formula is
bnk ¼ Pn ðlk Þ;
ð54Þ
with
lk ¼ cosðkp=2NÞ:
ð55Þ
This alternative bnk (Method 2 as referred by Cassell and Williams (2004) is used in Eq. (52) instead of Eq. (43) (Method 1). Since both forms of bnk are independent of either Bessel functions or HPM, it has neutral effect on the overall comparison between the results of this work and those reported by Cassell and Williams (2004). As described earlier in this subsection, the critical radius is calculated in this work by the reflection point technique. We
Table 4 Angular distribution over the curved surface using Method 2 and 62 terms (N = 30). Theta
HPM
Bessel functions
0 10 20 30 40 50 60 70 80 90
0.6943497 0.6904001 0.6688981 0.6336942 0.5846480 0.5219508 0.4459712 0.3564137 0.2520653 0.1282844
0.6943497 0.6904001 0.6688981 0.6336942 0.5846480 0.5219508 0.4459712 0.3564137 0.2520653 0.1282844
Table 2 Variation of the critical radius ac and critical mass M c per hemisphere with the interaction parameter g. g
ac ðcmÞ
Mc ðkgÞ
0 (sphere) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
8.44123 9.14444 9.73615 10.21877 10.60685 10.91835 11.16983 11.37485 11.54391 11.68495 11.80396
23.431 29.788 35.953 41.569 46.487 50.704 54.289 57.334 59.928 62.152 64.070
Fig. 3. Angular distribution over the curved surface.
Table 3 Critical radius for the hemisphere using two different methods. Terms
4 6 8 10 12 22 32 42 52
Method 1
Method 2
HPM this work
Bessel this work
Bessel Cassell and Williams (2004)
HPM this work
Bessel this work
Bessel Cassell and Williams (2004)
11.80012 11.80389 11.80395 11.80396 11.80396 11.80396 11.80396 11.80396 –
11.80012 11.80389 11.80395 11.80396 11.80396 11.80396 11.80396 – –
11.80008 11.80385 11.80391 11.80392 11.80392 11.80392 11.80392 11.80392 –
12.04073 11.79314 11.79867 11.80112 11.80234 11.80359 11.80380 11.80387 –
12.04073 11.79314 11.79867 11.80112 11.80234 11.80359 11.80380 – –
12.04069 11.79310 11.79863 11.80108 11.80230 11.80355 11.80377 11.80383 11.80386
1716
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
calculated the critical radius as a function of the number of terms, using both Bessel functions as well as HPM. The calculations are repeated twice, adopting Method 1 and Method 2 described above. The results are listed in Table 3. The calculated critical radius does converge using Method 1, while Method 2, which is proposed mainly for purposes of flux calculations, does not show any similar convergence, whatever number of terms are used. Method 1 in Table 3 reveals the fact that the calculation of the critical radius of the hemisphere, to accuracy of seven significant figures, converges with 10 terms using the reflection point technique, exactly similar to that reported by Cassell and Williams (2004). This convergence is independent of whether we are using HPM or Bessel functions with the reflection point technique. The small difference in the convergence value obtained in this work compared to Cassell and Williams (2004), however, is related to the difference in the implemented computational method.
3.2. Flux distribution It is worthwhile now to adopt Method 2 and examine the flux distribution inside the hemisphere. For this purpose, and in order to establish the computational grounds, we will write the flux given by Eq. (53) as
/ðr; lÞ ¼
N X ½A2n R2n ðBrÞP2n ðlÞ þ A2nþ1 R2nþ1 ðBrÞP2nþ1 ðlÞ:
ð56Þ
n¼0
Furthermore, and using Eq. (35), we finally arrive at the form that has been used in the code to compute the flux in the bare hemisphere, that is
/ðr; lÞ ¼
N X
" A2n R2n ðBrÞP2n ðlÞ þ
n¼0
N X
#
cln R2lþ1 ðBrÞP2lþ1 ðlÞ :
ð57Þ
l¼n
In order to calculate the amplitudes A2n , Eq. (51) has been rearranged as N X
b nk A2n ¼ M b 0k A0 ; M
ð58Þ
n¼1
and A1 has been normalized to unity. Hence, A0 ¼ c1 . The subse00 quent flux calculations can thus be scaled to absolute numbers by adjustment of the amplitudes. In addition, all the calculated flux values are normalized to the volume averaged flux /
¼ /
Ra R1 0
r 2 /ðr; lÞdldr 0 ; Ra R1 r 2 dldr 0 0
ð59Þ
the azimuth angle being neutral. Hence
¼ 3 / a3
Fig. 4. Flux distribution across the core of a bare hemisphere at different angles.
Z
a
dr r 2
0
Z
1
dl/ðr; lÞ:
ð60Þ
0
The value of a in Eq. (60) is 11.80396 cm, the convergence value of the critical radius (Table 3). For example, using 46 terms ðN ¼ 22Þ provides values for the volume averaged flux of
Table 5 Flux variation across the critical hemisphere at h ¼ p=4. r=a
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
4 Terms
12 Terms Bessel functions
HPM
Bessel functions
HPM
Bessel functions
0.887119 1.25636 1.55550 1.76260 1.86340 1.85253 1.73389 1.52015 1.23146 0.893463 0.534903
0.887119 1.25636 1.55550 1.76260 1.86340 1.85253 1.73389 1.52015 1.23146 0.893463 0.534903
0.881442 1.24814 1.54501 1.75055 1.85094 1.84118 1.72544 1.51647 1.23433 0.904402 0.554956
0.881442 1.24814 1.54501 1.75055 1.85094 1.84118 1.72544 1.51647 1.23433 0.904402 0.554956
0.881631 1.24841 1.54535 1.75095 1.85137 1.84163 1.72586 1.51684 1.23462 0.90457 0.554941
0.881631 1.24841 1.54535 1.75095 1.85137 1.84163 1.72586 1.51684 1.23462 0.90457 0.554941
46 Terms 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
22 Terms
HPM
0.881673 1.24847 1.54543 1.75104 1.85147 1.84172 1.72595 1.51692 1.23469 0.904621 0.554972
62 Terms 0.881673 1.24847 1.54543 1.75104 1.85147 1.84172 1.72595 1.51692 1.23469 0.904621 0.554972
0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72595 1.51692 1.23469 0.904622 0.554973
86 Terms 0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72595 1.51692 1.23469 0.904622 0.554973
0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72596 1.51692 1.23469 0.904622 0.554974
0.881673 1.24847 1.54543 1.75104 1.85147 1.84173 1.72596 1.51692 1.23469 0.904622 0.554974
K. Khasawneh et al. / Annals of Nuclear Energy 36 (2009) 1711–1717
1717
ism developed in this work and illustrated by Eqs. (57) and (58), it is worthwhile to compare our results to those obtained using deterministic finite element approximation to the neutron transport equation, which is implemented in the EVENT code, and reported in Cassell and Williams (2004). Fig. 5 illustrates that, though HPM can exactly reproduce the results obtained using Bessel functions, the computational implementation illustrated by Eqs. (57) and (58) yields results closer to those obtained using EVENT, especially when smaller number of terms in the expansion are used. 4. Conclusions The homotopy perturbation method, applied to the neutron diffusion problem in hemispherical symmetry, can exactly reproduce the results obtained using Bessel functions. On the other hand, the computational implementation of the analytical solutions obtained in this work does affect the results, and improve them when compared to those obtained using finite element codes like EVENT. Fig. 5. Percent error in the normalized flux values calculated relative to EVENT. Closed and open symbols correspond to results reported by Cassell and Williams (2004) and in this work, respectively.
0.092076 and 0.633771, using Bessel functions and HPM, respectively. This clearly indicates that scaling to absolute flux values differs accordingly. The flux behavior on the curved surface of the critical hemisphere is inferred from Table 4, for both HPM and Bessel functions with 62 terms in Eq. (56) ðN ¼ 30Þ. The remarkable agreement between both calculations is evident, taking into consideration the fact that the angular dependence itself is determined by the Legendre polynomials. Fig. 3 illustrates the fast convergence of the flux calculations on the curved surface with the number of terms. The critical flux across the core has been also calculated using HPM, at different angles, and is shown in Fig. 4. It is evident that at the flat surface ðh ¼ 90 Þ the flux decreases as we move towards the curved surface. Unlike the sphere, the maximum flux (namely 2.29786) is not at r ¼ 0 but at ðr; hÞ ¼ ð5:4 cm;0 Þ. The value of the flux at r ¼ 0 is 0.881671, which is significant compared with the maximum value, confirming the fact that the ZF boundary condition / ¼ 0 at the flat surface is not valid for such a small hemisphere. In order to compare flux calculations using HPM to those using Bessel functions (both calculated in this work), data in Table 5 provide normalized flux values across the hemisphere for h ¼ 45 using Method 2. It is again evident that HPM can perfectly reproduce the normalized flux values calculated using Bessel functions, establishing HPM as a new technique that can be reliably used for similar calculations. Furthermore, recalling that the flux has been computed, both using HPM and Bessel functions, using the formal-
Acknowledgment The financial support provided by King Abdullah II Fund for Development KAFD is gratefully acknowledged by K. Kh. References 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. 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., 2006. Homotopy perturbation method for solving boundary value problems. Physics Letters A 350, 87–88. He, J.H., 2006. Some asymptotic methods for strongly nonlinear equations. International Journal of Modern Physics B 20, 1141–1199. 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. Williams, M.M.R., 2002. Neutron transport in two interacting slabs. II: the influence of transverse leakage. Annals of Nuclear Energy 29, 1551–1596.