PALLAS-2DCY-FC, a calculational method and radiation transport code in two-dimensional (R, Z) geometry

ARCHIEF

PAPERS

OF

SHIP RESEARCH INSTITUTE

PALLAS-2DCY-FC, A Calculational Method and Radiation

Transport Code in TwoPimensonal (R, Z) Geometry

By

Kiyoshi TAKEUCHI

July 1979 Ship Research Institute

Tokyo, Japan

Lab. y.

Scheepsbouwkwide

Technische Hogeschool

Deli t

Papers of Ship Research Institute, No. 57 (July 1979)

PALLAS-2DCY-FC, A CALCULATIONAL METHOD AND

RADIATION TRANSPORT CODE IN TWO-DIMENSIONAL (R, Z) GEOMETRY*

By Kiyoshi TAKEUCHI

ABSTRACT

A numerical method for solving the steady state Boltzmann transport equation is presented for practical radiation shielding calculations. The

method is based on a direct integration of the equation along flight path of radiation in the direction of motion at each discrete-ordinate direction.

Instead of using a Legendre expansion approximation for treatment of

anisotropic scattering, the differential scattering cross section is directly used, which results in always positive solutions. In addition, an

iteration-convergence technique is not used in calculation for within-group scattering

source. As a consequence shortening of computer time and getting always reasonable solutions are attained.

A computer code, PALLAS-2DCY-FC has been written in

FORTRAN-IV on FACOM 230/75 and CDC 6600, which calculates radiation transport in two-dimensional (r, z) geometry. It includes a first collision source option for various source geometries such as point, line, disk and cylindri-cal sources. Informations on input specification are given for users with detailed data notes and sample problems.

Contents Abstract

Program abstract

Introduction Theory

Evaluation of the source term Q(r, 12, E)

1. 1 Calculation of the elastic scattering term

1. 2 Calculation of neutron inelastic scattering integration

1. 3 Calculation of photon scattering integration Solution of the integral transport equation, Eq. (2)

Calculation of the within group scattering radiations

without using an iteration technique

Expression in two-dimensional (r. z) geometry

First collision source

PALLAS group constants

Angular quadrature

Slowing down calculation into thermal group and

con-vergence

Received on January 8, 1979 Tokai Branch

Aldetmg Scheepsbuw- en Scheepvaartkuode

Te.hnk' e _{'oc;chroI, Deift}

DOCUMENTAT1E

I :1&-37

1. Input specification

2. Detailed data flotes

3. External and internal data files 4. Flow diagram with subroutine

5. Program mnemonic and program variable 6. Sample problems

References

I. PROGRAM ABSTRACT

Name of Program; PALLAS-2DCY-FC.

Computer for which designed or operable; FACOM 230/75, CDC6600.

Nature of physical problem solved; PALLAS-2DCY-FC solves the steady state Boltzmann transport equation in two-dimensional (r, z)

geometry. Only fixed-source problem can be solved, so that its

application is restricted to radiation shielding problems.

Method of solution; The method of direct integration of the

trans-port equation is used, and the integration of the equation is performed

along the flight path of particle in the direction of motion at each

discrete-ordinate direction. Neither a Legendre expansion

approxi-mation nor an iteration and convergence technique are used for treatment of anisotropic scattering and getting the solution of the

equation. Instead, the differential scatterring cross section data are directly used for neutron, and the Klein-Nishina formula is used for

photon. The use of fine energy meshes should be remommended for computing the equation.

.5. Restrictions; 254 K words are required for the total core storage for

(55x55) spatial and 28 fixed directional meshes.

(3. Typical running time; About 0.0205 sec/spatial mesh/group is required

as C.P.U. time on FACOM 230/75.

7. Usual features; Group constants for neutron shielding calculation

are read in from the PALLAS library tape. Linear attenuation co-efficients and pair production coco-efficients for photon are read in from

cards. Special features which make PALLAS-2DCY-FC well adapted

to shielding problems include the capability to output fluxes at

ex-terior or inex-terior boundaries for use in bootstrap problems, and an

analytic first-collision source calculation for various source configu-rations such as point, line, disk and cylindrical sources that mitigates ray-streaming effects in certain problems.

8. Machine requirements; Card input, printed and the scratch data

sets may be located on any external storage device.

9. Languages; FORTRAN IV.

10. Material available; Input description, source deck for FORTRAN,

where R=i'

, = and

"=iR'Q, and .Y,X,(7", E)=const.

II. Introduction

The code PALLAS-2DCY-FC is based on the method of direct

inte-gration of the steady state transport equation for neutrons and photons to describe the radiation transport in (r, z) two-dimensional geometry.

As this method is fundamentally based on the discrete ordinates method, the Boltzmann transport equation is approximated by a set of equations describing relation between solution and variables on each mesh point of

network in a phase space of energy, space and angles. In the network,

an energy mesh approximation is used and the angular variables are meshed into a set of discrete directions with weights. The spatial

vari-able is partitioned into (r, z) mesh points. Particle multiplication is not

allowed in the current code. Only fixed source problem can be treated

with the vaccum boundary condition at the outermost boundaries.

Inner iterations are not necessary for getting the solution. As a result, one can get always a solution without any fear of no convergence.

Outer iterations are not allowed, so that its application is restricted to

only radiation shielding problems. A Legendre expansion approximation is not used for treatment of anisotropy-calculation. Instead, the differ-ential scattering cross section is used, which is formed by using the full

data of Legendre coefficients for neutron from the ENDF/B file or by

using the Klein-Nishina formula for photon. As a result, always positive

reasonable scattering source is calculated, so that no negative fix-up

technique is necessary. To mitigate ray-streaming effects an analytic

first-collision source technique is utilized for various source configurations

such as point, line, disk and cylindrical sources. III. Theory

The time-independent transport equation is written with usual nota-tions as

V()+X,(, E)«, Q, E)

=$5

dE'd?Y, Q', E')I,(i, E'E, .?'

.

+S(, [2, E).

₍₁₎

Above equation can be rewritten in the following integral form

as-suming that the total cross section is constant within any spatial interval (i', ).

Further,

Q(, ?.),

E)=JJQ, E'-+E,

û', E')d?.?'dE'+S(i, û, E).

(3)

The angular variable Q is represented by discrete-ordinate directional

points on a unit sphere, in which !7,=û(O», Çpq). The O is the polar angle and the ç is the azimuthal angle. In PALLAS calculation, wcos O is used instead of O (see Figs. i and 2).

z

Fig. 1. Moving direction Q of radiation in PALLAS

II

22

. Q33

6145

1-6U 32 31 »43

Ç42 iI

53

5l

62

_6i

Fig. 2. Angular mesh points and their order in the current PALLAS code.

(9)

1.

Evaluation of the source term Q(,

, E)

in order to perform numerical integration of the first term on the right hand side of Eq. (3), the flux is transformed into the energy flux

I(, !2, E);

I(, Q, E)=E(, Q, E).

Further, the differential scattering cross section is represented by

I,(, E'E, ')=I,(, E')f,.,(E', 4ô(cos &a)

+1)2

2AE'

for neutron elastic scattering,

=.Y (, E')

47r

for neutron inelastic scattering,

K(E'E)/E

=n(r) i

2ir

'E'

for photon scattering,

where E) and E) are respectively the total elastic and ine lastic scattering cross sections, A is the mass number of nuclide and ¡

is the cosine of the polar angle of elastic scattering in the center of mass

system, while is the cosine of the same angle (C-)) in the laboratory system. The azimuthal angle of elastic scattering is represented by fr.

Further, f,(E, u) is the scattering distribution function and is normalized to 1.0 as

2 f f,,(E, p)d1i= 1,

(7)

and f,,(E'-±E) is the slowing down probability assuming that the inelastic scattering is isotropic in the laboratory system.

For photon scattering the compton scattering cross section is exactly

represented by the Klein-Nishina formula, so that the scattering kernel

K(E', E) is derived from the Klein-Nishina formula and is given2

3

/2'

K(2', A) =

+

(í)

+ 2(2' - A)+(A'A)2},

(8)

where 2 is the compton wavelength;

2=

m0c20.511

E(MeV) E

Substituting the expression (4) into the integral term of Eq. (3), we can write the integral term of Eq. (3) for neutron with the energy flux

notation as

QQ, ,E)=JJJ.Y,1(, E')f,(E', ¡t)ô(cos

9)

(A

X I(i, ', E')dedfrdE'.

From the relation between ¡t and E of

(A+1)2 (

E \

2A \

E')'

2A E'2

dE

=---(AH-1)2 E

the above expression (10) is rewritten as

Q; _=f5X1(, E')f,,(E', p)I(, ', E')dtdifr. (13)

Numerical integration on ¡t is performed by using Gaussian quadrature.

For this purpose ¡t is represented by discrete points ¡t,,, and f,1(E', ¡t) is

also represented by f;'(E') at discrete meshes ¡j,,,. On the other hand the

numerical integration on the azimuthal angle fr is performed by using

the formula of scattering angle of the following relation (14).

cos O' = cos O cos e + /(i - cos2 8i - cos2 ) cos /,. (14)

The relation (14) may be rewritten as

w'=ox+ sf

2)(1_2jcos /r. (15)

The value of can be determined by using the following expression (16)

representing the relation of scattering angles between in the center of

mass system and in the laboratory system;

Apm+1

-- 2At,. + 1 for ¡t=í,,,. (16)

If the w' in the relation (15) is represented by an angular mesh w. for w=w,, and a=«,, then the value of fr is calculated for a set of discrete

values (w,,, w, a,,,) by transforming the relation (15) as

f

Wpm

\ (1-1a)

(17)

(10)

Fig. 3. Relation between azimuthal angle of scattering 'fr and directional cosine w,, in scattering.

Now it is possible to calculate numerically the Q, at and E for a

discrete angle !2=L?j,q(S2(w,, ç,,q)). To carry out the numerical calculation,

we represent the ifr obtained from the expression (17) by a discrete mesh

1/1. in accordance with the angular mesh w,, as shown in Fig. 3. The

weight w for the -/r,, to be used in the numerical integration over fr is

determined as

= fr(w;,(iJ, am) - 1/r(w,, a,,,), (1.8)

where the w and w, are the values of upper and lower bounds of w,, range

respectively. Thus the value of Q) at i and E for is calculated by numerical integrations as

Q, Q, E) =

W,. E,,,)f,,",'(E,,,)I(r, Q,,,, E,,,), (19)

where W,,, is the weight for in the Gaussian quadrature, and Q,,,

can be determined by using the alternate expression of the relation of

scattering angle;

a = ww' + /(i w2)(1_ w') cos (20)

To determine ç' the relation (20) is transformed and written in the form

Fig. 4. Variation of azimuthal angle of flight direction in before and after scattering.

ç'; =ç1±cos(-

(21)

If the two values of ç'.. belong to the regions of Jç. and 4çb,,,,,

respec-tively, then ç5 can be represented by discrete meshes çí,,, and çl, (see

Fig. 4). The remaining unknown variable E,. is determined for E=E!

and 1o=i,,, by using the relation between energy and scattering angle in the center of mass;

(A+1) _E

A2+2Ap+1

For simplicity in practical calculations, the energy E,. is determined in

the following simple manner, which is adopted in the current PALLAS

code: In the PALLAS program the energy mesh E, is chosen, in princi-ple, at a midpoint of the E energy region. An exception is in the case

of the energy region of energy mesh E prior to the energy mesh E at

which the calculation is to be carried out. The energy region of E_1 is

neccessarily chosen from just above the E up to the midpoint energy

between E and E , on account of avoiding the use of an inner-iteration

technique for the calculation of within-group scattering, which will be

described in detail in the succeeding section 3. Then if the E,,, calculated

(25)

by expression (22) belongs to the energy region of k-th energy mesh E,

the values of Ei), f7(Ek.) and I(, ?,,,, Ek) may be chosen as those of

X,(r, E,,), f(E,,,) and E,,,).

1. 2 Calculation of neutron inelastic scattering integration

Based on the assumption that the inelastic scattering is isotropic in the laboratory system, the integral term of the expression (3) is simply written in the form

Q,(r, 17, E)=$ X1,(, E')f''

')E10(,

E')dE', (23)

47rE

where

I,(, E')=$I(, 17, E')dQ. (24)

The above expression (23) is readily calculated numerically for the energy

mesh E;

Q(r, 17, E1) =

_ki

E,,) f(F, E,,-+E,)E I(f, Eh.)JEk 47rEk

X,,(f, E,,-+E)E,

L(, E,)zJE,,

k=1 4irE1

where X1,,(, E-+E,)=I1(, Ek)ff,(, E1-E,). The inelastic slowing down cross sections X1,,(, E-E,) for each nuclide are prepared in the PALLAS

library.

1. 3 Calculation of photon scattering integration

Using the expression (8) for the photon scattering cross section, we can write the photon scattering integral term of the expression (3) as

QG( _{Q, 2)= JJ n(? I(, 17', A')d17'dÄ'.}

(26)

In the compton scattering the relation between the scattering angle and

the change of the wavelength is given in the following form

p=1±A'-A. (27)

Then d.=dA, and the expression (26) is rewritten as

Q6Q, 17, 2)= if n() K(A', 2) 1Q, 17', 2')d/rd2'. (28)

neutron elastic scattering integral calculation. In the photon calculation

a,. = , and ¡t,,, = i + 2,, 2. Then

¡ WWpj4

i/tm afr(w,,, w,,, p,,) = cos

/(i

-,)

and

w,, = ifr(w, (I) p,/,,,) - ifr((t), w,,, pm)I. (30)

On the other hand the numerical integration on 2' is carried out by

re-presenting the wavelength by discrete mesh A. Thus

Q(;(,

pq,2,) = w,n(i) K(Am, A.J( Ä,,,, A,)J2,,,, (31)

where Q,,, is determined in the same way as in the case of neutron elastic scattering calculation and

Ç' =q5pq±CO5

(

P!)

- w,)(1 - w,) (32)

The mesh assignment for q5' may be made in the same way as in the

neutron case and ç5'=q5,, and q5'=q5,,,1.

From the calculations described above, Q(, ?, E) in the expression (3) is calculated numerically for a set of discrete values

and E=E,) by the relation:

Q(,, Lpq, E) = Qfi(i, E,) + pq E,)

± S(1, _,, E.,) (33)

for neutron,

Q6(,, ?pq, E,)+S(4, b,,,,, E,) (34)

for photon.

2.

Solution of the integral transport equation, Eq. (2)

The spatial variable is partitioned into r,(i=i, 2, , MRR) so that the total cross section may be constant in any interval (.,, ). To per-form the integration over R' in Eq. (2), the function Q(ì', , E) must be approximated by a rather simple function with respect to r.

In the

original PALLAS code3, this function is approximated by a linear

func-tion with respect to . Recently N. Sasamoto has introduced an

expo-nential approximation4 for this function to the current PALLAS code,

which has resulted in much more accurate radiation transport calculations in attenuating media in shields.

The former approximation brings the following expression with respect

to i:

R Q(",Q, E) Q(") = a(?) + b(1)R', (35) then j'R (a(',) + b(1)R') exp (- X,R')dR' [Q(1){X,R+exp(X,R)-1} + Q(1){1 (1 + ,R) exp (XR)}J. (36)

On the other hand, the latter approximation results in _{an expression4} Q(i",i, E)mQ(?")==a(r,) exp {b(r1)R'}, (37)

then

a(i,) exp { b(,)R'} exp(

-{Q(r1_1)exp(X,R)Q(r)}R

(38)

in

'R

Since the first term on the right hand side of Eq. (2) is an analytical

expression, the numerical error in a solution of the transport equation,

Eq. (2), depends on that of the numerical integration of the second term.

In the current PALLAS-2DCY-FC code, the choice of the expression (36)

or (38) _{is determined depending on the ratio of Q(.j/Q(i1) except for} Q(4)=O; In the case of l/2<Q(1)/Q(1)<2, the expression (36) is auto-matically chosen in the code, while in the other values of the ratio the

expression (38) is automatically chosen. In the case of Q(i,)=O, the ex-pression (36) is chosen in the code.

From the integral calculation described above, Eq. (2) is solved and

written with the notation of the energy flux as

Qpq, E)=I(4, Q, E3) exp {X(Y, E)R1}

+ F{Q(r,- ), Q(), Ei), R.}, (39)

where

E) = E1fJ(F1, Q0, E,),

1Expression (36) for QQ,)=O, or Q(P,)+O and

F(

)_J

1f2<Q(i)/Q(4)<2,

Expression (38) for

3.

New technique for treatment of the within-group scattering

radiations instead of using an iteration technique

Those radiations scattered within a certain small angle may be thought

to be unscattered ones based on the idea of the removal theory.

Con-sequently PALLAS program adopts this idea and the radiations scattered

within the small region of 4.?7pq from their flight direction of , are treated as unscattered. For simplicity in the practical calculation of PALLAS code, this small region of scattering angle is assumed to be constant for all the angular meshes such as the size of weight W, of the

first mesh , in the scattering distribution function f,(E,, ,) for neutron

or the weight JA,(=(A,-2,)/2) in the scattering

kernel K(2, 2f,). Then the scattered radiations to be treated as unscattered are calculated as follows:

Q,,,,, E,) = 27rW,X,(, E)f,(E, ¡i,)I(, Q,,,,, E,) (40)

for neutron,

pq, A)=27rn(ì)

K(2,, _{I(, !,,,,, E1)4A (41)}

2r

for photon.

Subtracting the expression (40) or (41) from the both sides of Eq. (1), we obtain the following equation with the notation of the energy flux.

P'(QI)+I,I(i, Qpg, E)A(r, E1)I(i, Qpq, E)

=

ft dE'd'(E,/E')I(, i)', E')X,(, E'E5,

- A(, E)I(F, Q,,,,, E,) +S(J2,,,,, E),

where

(42)

E)f,(E, p,)

A(i, E,)= for neutron, (43)

n()K(2, 2)4À

for photon.

The above equation (42) can be rewritten as follows:

(I) +

_{dE' J} dY(E,/E')IX. + S, (44)

where

E>Ej

X(i, E)=X,Q, E)A(i-, E). (45)

and

I(,, ?pq, E,)=I(, Q, E,) exp (XR,)

+ F{Q'(, ), Q'(,), E,), R,}, (46)

where

= Q'' + Q + S for neutron,

=Q'°+S

for photon,

Q'

WmWX,j(, E,,,)f,i(Em, /im)I(, n,,,, E,,), (47)

= " wnQ) K(Arn, A,) _{I(i, (fl3 Am)JA,,.}

(48)

= i ,t _2,r

4. _{Expression in two-dimensional (r, z) geometry}

In this section the final equation, Eq. (46), is expressed in

two-dimen-sional (r, z) geometry. _{For this purpose the spatial meshes} , and r,

are expressed, respectively, by (r,, z,) and (r,,, z') or (r',Zk±1) as shown

in Fig. 5. Then,

I(,

E,) =I(r,, 2k Q,,,, E,),

Fig. 5. Illustration of spatial mesh for evaluation of 1(r,_i) and Q'(r,1).

where Q'(, .fpq, E,)= Q'(r1, Z, ?.pq, E?), for i=1,2,

,MRR

k=1,2,,MZZ.

z,,Z,,1I (O

I(±, z',

, E,)

I( _,?, E,) = I(', z,,.1, , E,)

z', ?, E,),

Q'(r,.,, z', , E,)

, 2, E) = Q'(r', z,,+, .), E,)

Q'(r,, z', ., E,),

r - <r' < r

} (+r, +z) direction, ((a) in Fig. 5) r<r'<r1+1 } z,, - 1<z' < z,, z,,<z'<z,,1 } (+r, --z) direction, ((c) in Fig. 5) z (r, +z) direction, ((b) in Fig. 5) r

Fig. 6. Azimuthal angle of flight direction varies with its motion in (r, z) geometry.

r1<r'<r,+1

} (r, z) direction. ((d) in Fig. 5)

To evaluate I and _Q' in expressions (50) and (51), an interpolation must

be used wìth respect to

r or z meshes. Furthermore, the direction

Q(r2(1,, _{)) included in expressions (50) and (51) must be determined since}

is not equal to Qpq(J2(w,,,ç5)) as shown in Fig. 6. If the value of çl is

in the region between ç,,,, and çi,,, , then the functions I and

Q' in

ex-pressions (50) and (51) may be determined by an interpolation with

re-spect to (Ø,

_,, ). After all the functions I and Q' are determined by

applying double interpolations with respect to spatial and directional meshes.

First collision source

For a point source on a certain point z on the z-axis, the uncollided flux is calculated quite easily as

ESe

I

(rz)=

4tR2

where

R= /r2+(z_z,)1.

In addition to a point source, Y. Kanai has proposed

a calculational

technique5 based on the point kernel method for various sources such as

cylindrical, disk and line sources and has added it to the

PALLAS-2DCY-FC code. The expression is written as

=J dr,. exp (-j'

E)ds) ô(Q . Q 1)

I(

? ₍₅₄₎

r-r,I

where I" (,, ) is the uncollided energy flux at t, and flight direction

, in shields, and r, and , are respectively, the position and flight

di-rection on the surface of a cylindrical volume source.

The use of the FC (First collision) option is very much effective for

mitigating the ray effect, though it results in

_{an increase in running} time.

PALLAS group constants

The group constants X,, X, and X or reaction constant X,,, for

neu-tron are generated by using the SUPERTOG code included in

RADHEAT-V3' code system. _{For this purpose the upper and lower bounds of the}

energy mesh E, must be specified. _{In the PALLAS calculation neutron}

slowing down calculation is usually made in the equal lethargy interval

Table 1. Energy group structures applied to PALLAS library. 0.1, 0.2. 0.4 and 0.8 lethargy width

structures are included.7

Group

0.1 Lethagry width structure

(eV)

E0 (eV) EL (eV)

i i.4208E+07 i.4018E+07 1.3499E+07

2 1.2856E+07 1.3499E+07 1.2214E+07

3 1.1633E+07 1.2214E+07 1.1052E+07

4 1.0526E+07 i.1052E+07 1.0000E+07

5 9.5242E+06 1.0000E+07 9.0484E+06

6 8.6178E +06 9.0484E +06 8.1873E +06

7 7.7977E+06 8.1873E+06 7.4082E+06

8 7.0557E+06 7.4082E+06 6.7032E+06

9 6.3842E+06 6.7032E+06 6.0653E+06

10 5.7767E-f-06 6.0653E+06 5.4881E+06

il

5.2270E+06 5.4881E+06 4.9658E+06

12 4.7296E+06 4. 9658E +06 4. 4933E +06

13 4.2795E+06 4.4933E+06 4.0657E-f-06

14 3.8722E+06 4.0657E+06 3.6788E+06

15 3.5037E+06 3.6788E+06 3.3287E+06

16 3.1703E+06 3.3287E+06 3.0119E+06

17 2.8686E+06 3.0119E+06 2.7253E+06

18 2.5956E+06 2.7253E-i-06 2.4660E+06

19 2.3488E+06 2.4660E+06 2.2313E+06

20 2.1251E+06 2.2313E+06 2.0190E+06

21 1.9229E-f-06 2.0190E+06 1.8268E+06

22 1.7399E+06 1.8268E+06 1.6530E+06

23 1.5743E+06 1.6530E+06 1.4957E+06

24 1.4245E+06 1.4957E+06 1.3534E+06

25 1.2890E+06 i.3534E+06 1.2246E+06

26 1.1663E+06 i.2246E+06 1.1080E+06

27 1.0553E+06 1.1080E+06 1.0026E+06

28 9.5488E+05 1.0026E+06 9.0718E+05

29 8. 6401E +05 9.0718E +05 8.2085E +05

30 7.8179E+05 8.2085E+05 7.4273E+05

31 7.0739E+05 7.4273E+05 6. 7205E +05 32 6.4008E+05 6.7205E+05 6.0810E-f-05 33 5.7917E+05 6.0810E+05 5.5023E-f-05

34 5.2405E+05 5.5023E+05 4.9787E+05

35 4.7418E-f-05 4.9787E+05 4.5049E+05 36 4.2906E+05 4.5049E +05 4.0762E +05

37 3.8823E+05 4.0762E+05 3.6883E+05

38 3.5128E+05 3.6883E+05 3.3373E+05

39 3.1785E-f-05 3.3373E+05 3.0197E+05

40 2.8761E+05 3.0197E+05 2.7324E+05

41 2.6024E +05 2.7324E +05 2.4723E +05 42 2.3547E+05 2.4723E+05 2.2371E-f-05

43 2.1306E+05 2.2371E+05 2.0242E+05

44 1.9279E+05 2.0242E+05 1.8316E+05

45 1.7444E+05 1.8316E+05 1.6573E+05

46 1.5784E+05 1.6573E-f-05 i.4996E+05

47 1.4282E+05 1.4996E+05 1.3569E+05

48 1.2923E+05 1.3569E+05 1.2277E+05

49 1.1693E+05 1.2277E+05 1.1109E+05

Table 1. (continued)

E1 (eV)

Group

E5 (eV)

0.2 Lethagry width structure

E, (eV)

i 1.4208E+07 1.5625E + 07 1.2792E +07

2 1.1633E+07 1.2792E + 07 1. 0473E +07

3 9.5242E+06 1. 0473 E + 07 8.5749E +06

4 7.7977E+06 8.5749E +06 7.0206E +06

5 6.3842E+06 7.0206E +06 5.7479E +06

6 5.2270E+06 5. 7479E +06 4.7060E +06

7 4.2795E+06 4. 7060E +06 3. 8530E +06

8 3.5037E+06 3. 8530E +06 3.1545E +06

9 2.8686E+06 3.1545E +06 2.5827E +06

10 2.3486E+06 2. 5827E +06 2. 1145 E + 06

11 l.9229E+06 2. 1145E +06 1. 7312E +06

12 1.5743E+06 1. 7312E +06 1. 4174E +06

13 l.2890E+06 1. 4174 E +06 1.1605E +06

14 l.0553E+06 1.1605E +06 9. 5013E +05

15 8.6401E+05 9.5013E +05 7.7790E +05

16 7.0739E+05 7.7790E +05 6.3689E +05

17 5.7917E +05 6.3689E +05 5.2144E +05

18 4.7418E+05 5. 2144E +05 4. 2692E +05

19 3. 8823E +05 4. 2692E +05 3. 4953 E + 05

20 3.1785E+05 3. 4953E +05 2. 8617E +05

21 2.6024E+05 2. 8617E +05 2.3430E +05

22 2.1306E+05 2.3430E +05 1. 9183E +05

23 1.7444E +05 19183E +05 1.5706E +05

24 1.4282E+05 1.5706E +05 1.2859E + 05

25 1.1693E+05 1. 2859E +05 1.0528E +05

26 9.5735E+04 1. 0528E +05 8.6194E +04

27 7.8382E+04 8. 6194E +04 7.0569E +04

28 6.4173E+04 7. 0569 E + 04 5.7777E +04

29 5.2541E+04 5.7777E +04 4. 7304E + 04

30 4.3017E+04 4. 7304E +04 3.8729E +04

31 3.5219E+04 3. 8729 E ± 04 3.1709E +04

32 2.8835E+04 3.1709E +04 2. 5961E +04

33 2.3608E+04 2. 5961E + 04 2.1255E +04

34 1.9329E+04 2. 1255 E + 04 1. 7402E +04

35 1.5825E+04 1. 7402E +04 1.4202E +04

36 1.2956E+04 1.4248E +04 1.1665E +04

37 1.0608E+04 1.1665E +04 9. 5505 E + 03

38 8.6849E+03 9. 5505E +03 7.8193E +03

39 7.1106E+03 7.8193E +03 6.4019E +03

40 5.8217E+03 6.4019E +03 5. 2414 E + 03

41 4.7664E+03 5. 2414E +03 4. 2913 E + 03

42 3.9024E+03 4. 2913E + 03 3. 5134 E + 03 43 3. 1950E +03 3. 5134E +03 2. 8766 E + 03

44 2.6158E+03 2.8766E +03 2. 3551E +03

45 2. 1417E +03 2. 3551E + 03 1. 9282E +03

46 1. 7535E +03 1. 9282E +03 1.5787E +03

47 1.4356E+03 1. 5787E +03 1.2925E +03

48 1.1754E+03 1.2925E +03 1. 0582E +03

49 9.6232E+02 1. 0582E +03 8.6640E +02

so that the upper bound of E may be chosen as (U1 + U)/2 and the lower

bound may be chosen as (U1-F U1)/2. In which U is lethargy of E1 and

U, = lfl(Emx/Ej). The weighting function for averaging cross sections is the fission+1/E spectrum for the present PALLAS library7.

It should be

noted that the use of the weighting function of fission/71(E)+1/X1(E)E

spectrum is much better. We will generate next PALLAS library with

this better weighting spectrum. Anisotropy in elastic scattering is treated

Table 1. (continued)

0.4 Lethargy width structure

Eq (eV) (eV) E1 (eV)

Group

1 1.4208E+07 1.7013E+07 1.1404E+07

2 9.5242E+06 1.1404E+07 7.6443E+06

3 6.3842E+06 7.6443E+06 5.1242E+06

4 4.2795E+06 5.1242E+06 3.4348E+06

5 2.8686E+06 3.4348E+06 2.3024E+06

6 1.9229E +06 2.3024E +06 1.5434E +06 7 1.2980E-u06 1.5434E+06 1.0345E+06 8.6401E+05 1.0345E+06 6.9348E+05 5.7917E+05 6. 9348E +05 i 4.6485E+05 lU 3.8823E+05 4.6485E+05 3.1160E-f-05

11 2.6024E+05 3.1160E+05 2.0887E+05

12 1.7444E+05 2.0887E+05 1.4001E+05

13 1.1693E±05 1.4001E+05 9.3852E+04

14 7.8382E+04 9.3582E+04 6.2911E+04

15 5.2541E+04 6.2911E+04 4.2170E+04

16 3.5219E+04 4.2170E+04 2.8268E+04

17 2.3608E+04 2.8268E+04 1.8948E+04

18 1.5825E+04 1.8948E+04 1.2702E+04

19 1.0608E+04 1.2702E+04 8.5141E+03

20 7.1106E+03 8.5141E+03 5.7072E+03

21 4.7664E+03 5.7072E+03 3.8256E+03

22 3.1950E+03 3.8256E+03 2.5644E+03

23 2. 1417E +03 2.5644E +03 1.7190E +03

24 1.4356E+03 1.7190E+03 1.1523E+03

25 9.6232E+02 1.1523E+03 7.7238E+02

26 6.4506E+02 7.7238E+02 5.1774E+02

27 4.3240E+02 5.1774E+02 3.4705E+02

28 2.8984E+02 3.4705E+02 2.3264E+02

29 1.9429E+02 2.3264E+02 1.5594E+02

30 1.3024E+02 1.5594E+02 1.0453E+02

31 8. 7299E +01 i 1.0453E+02 7.0069E+01

32 5.8519E+01 7.0069E±01 4.6968E+01

33 3.9226E+01 4.6968E+01 3.1484E+01

34 2.6294E+01 3.1484E+01 2.1104E+O1

35 1.7625E+01 2.1104E+01 1.4147E+01

36 1.1815E+01 1.4147E+01 9.4828E+00

37 7.9196E+00 9.4828E+00 6.3565E+00

38 5. 3087E +00 6. 3565E +00 4.2609E+00

39 3.5585E+00 4.2609E+00 2.8562E+00

40 2.3853E-f-00 2.8562E+00 1.9145E+00

41 1.5989E+00 1.9145E+00 1.2834E+00

42 1.0718E-j-00 1.2834E+00 8.6026E-01

43 7.1845E-01 8.6026E-01 5.7665E-01

44 4.8159E-01 5.7665E-01 3.8654E-01

Group

Table 1. (continued)

0.8 Lethargy width structure

E, (eV) (eV) E1 (eV)

1 1.4208E+07 1.9607E+07 _8.8099E+06

2 6.3842E+06 8.8099E+06 3.9586E+06

3 2.8686E+06 3.9588E+06 1.7787E+06

4 1. 2890E +06 1.7787E+06 7.9922E+05

5 5.7917E+05 7.9922E+05 3.5911E+05

6 2.6024E+05 3.5911E+05 1.6136E+05

7 1.1693E+05 1.6136E+05 7.2504E+04

8 5.2541E+04 7.2504E+04 3.2578E+04

9 2.3608E+04 3.2578E+04 1.4638E+04

10 1.0608E+04 1.4638E+04 6.5774E+03

11 4.7664E+03 6.5774E+03 2.9554E+03

12 2.1417E+03 2.9554E+03 1.3279E+03

13 9.6232E+02 1.3279E+03 5.9669E+02

14 4.3240E+02 5.9669E+02 2.6811E+02

15 1.9429E+02 2.6811E+02 1.2047E+02

16 8.7299E+01 1.2047E+02 5.4130E+01

17 3.9226E+01 5.4130E+01 2.4322E+01

18 1.7629E+01 2.4322E+01 1.0929E+01

19 7.9196E+00 1.0929E+01 4.9106E+00

20 3.5585E+00 4.9106E+00 2.2065E+00

21 1.5989E+00 2.2065E+00 9.9143E-01

22 7.1845E-01 9.9143E-01 4.4548E-01

23 3.2281E-01 4.4548E-01 2.0017E-01

precisely by using the angular distribution ¡(E, ¡) for each nuclide, which is produced at all the energy meshes E based on the Legendre coefficients

taken from the data in ENDF/B-IV. Inelastic slowing down cross sections

f(E-8E) for each nuclide are generated by using SUPERTOG routine

included in RADHEAT-V3 and compiled in PALLAS library.

On the other hand, cross sections for photon are directly read in

from cards. The linear attenuation coefficients _¡ may be prepared for

energy meshes E,, which are obtained by an interpolation of the compiled

data of photon cross sections8. In addition, the linear pair production

cross sections 1»(cm') may be also prepared by the same way as in the

production of u. The Klein-Nishina formula is directly applied in the

PALLAS program for obtaining compton scattering cross sections.

Sec-ondary photon production cross section data for PALLAS can be prepared using the data compiled in POPUP-49 in the RADHEAT-V3.

The neutron cross section data compiled in the PALLAS library are

generated for four energy-mesh structures; energy meshes divided by an

equal lethargy width of 0.1, 0.2, 0.4 or 0.8. These energy meshes are given together with their upper and lower energy bounds in Table 1 in

order of lethargy widths of 0.1, 0.2, 0.4 and 0.8. PALLAS executes the

neutron energy group calculation in a series of these lethargy-width

structures, for instance a series of structures of 0.2, 0.4 and 0.8 lethargy widths for fast, intermediate and low energy ranges (see Data note 1).

PALLAS code calculates the transport equation without using

itera-tive calculations for within group scattering radiations, which brings

shortening computer time. Then multi-group radiation calculations in

fine energy mesh structures have become possible even in two-dimensional

transport calculations.

At least 35 or 20 energy meshes are used in

two-dimensional calculations respectively for neutron or photon. Much

finer energy meshes can be chosen within a reasonable computer time.

In order to treat precisely the anisotropy in elastic scattering of

neutrons, Legendre polynomial expansion approximation is not applied

in the PALLAS calculation. Instead, the scattering distribution function

[(E, t) is directly used, which is formed by using the data of Legendre

coefficient taken from the ENDF/B-IV. The maximum number of terms

Table 2. Maximum number of terms of Legendre polynomials and the energy

group number and the corresponding energy above which the anisotropy of elastic scattering is taken into account.7

maximum number of terms of Legendre polynomials.

the energy group and the corresponding energy above which the anisotropy

of elastic scattering is taken into account.

1.42+7 reads 1.42X107.

Jn=0. i 4u=0.2 Au=O.4 .iu=0.8

Nuclide LmY

.

Grp. E ¿**)_{(eV) (e\)} Grp. _(eV) Grp. (eV)E

il-1 i O 1.42+7n 0 1.42+7 0 1.42+7 0 1.42+7 Li-6 9 50 1.06+4 37 1.06+4 19 1.06+4 10 1. 06+4 Li-7 10 43 2.13+5 22 2.13+5 11 2.60+5 6 2.60+5 B-10 9 46 1.58+5 23 1.74+5 12 1.74+5 6 2.60+5 B-11 9 46 1.58+5 23 1.74+5 12 1.74+5 6 2.60+5 C-12 7 50 1.06+5 37 1.06+4 19 1.06+4 10 1.06+4 N-14 11 26 1.17+6 13 1. 29+6 7 1. 29+6 4 1.29+6 O-16 11 50 1.06+5 25 1. 17+5 13 1. 17+5 7 1.17+5 Na 15 50 1.06+5 31 3.52+4 16 3.52+4 8 5.25+4 Mg 16 50 1.06+5 31 3. 52+4 16 3. 52+4 8 5.25+4 Al-27 11 43 2.13+5 22 2.13+5 11 2.60+5 6 2.60+5 Si 11 50 1.06+5 25 1.17+5 13 1.17+5 7 1. 17+5 Ca 11 50 1.06+5 25 1.17+5 13 1.17+5 7 1. 17+5 Cr 16 50 1.06+5 48 1.17+3 24 1.43+3 12 2. 14+3 Mn-55 16 32 6. 40+5 16 7.07+5 8 8.64+5 4 1.29+6 Fe 12 39 3. 18+5 20 3.18+5 10 3.88+5 5 5.79+5 Ni 13 50 1.06+5 28 6.41+4 14 7.83+4 7 1.17+5 Zr 16 50 1.06+5 31 3.52+4 16 3.52+4 8 5.25+4 Mo 15 50 1.06+5 25 1.17+5 13 1.17+5 7 1.17+5 Pb 15 41 2.60+5 21 2.60+5 11 2.60+5 6 2.60+5 U-235 16 43 2. 13+5 22 2. 13+5 ii 2.60+5 6 2.60+5 U-238 16 50 1.06+5 35 1.58+4 18 1.58+4 9 2.36+4

Table 3. Identification numbers of PALLAS library7)

of Legendre coefficients included in the ENDF/B-IV data library is

uti-lized for producing the distribution function f(E, p), which is given in

Table 2 for each nuclide compiled in the PALLAS library.

The numbers composed of 4 digits used for identifying nuclides are given in Table 3, which are utilized when neutron cross section data are

selected from the PALLAS library. The number of nuclides included in

the library is restricted for use in shielding calculations. The increase of the number of nuclides is encouraged for wide use of PALLAS code.

Since PALLAS calculation of neutron scattering is made precisely

for each nuclide, the number of nuclides composed of material must be specified in each material region by an input, and identification numbers and nuclear densities of these nuclides are also specified in each material region by input data.

On the other hand for photon cross section data, linear attenuation and pair production coefficients are prepared for each material region, and in addition, the compton scattering cross sections are calculated

Nuclide Ju=0.1 group width Ju=0.2 No. ID No in lethargy 4uO.4 Ju=0.8

No. ID No No. ID No No. ID No

H-1 i 1011 23 2011 45 4011 67 8011 Li-6 2 1036 24 2036 46 4036 68 8036 Li-7 3 1037 25 2037 47 4037 69 8037 B-10 4 1050 26 2050 48 4050 70 8650 B-11 5 1051 27 2051 49 4051 71 8051 C-12 6 1062 28 2062 50 4062 72 8062 N-14 7 1074 29 2074 51 4074 73 8074 O-16 8 1086 30 2086 52 4086 74 8086 Na 9 1113 31 2113 53 4113 75 8113 Mg 10 1120 32 2120 54 4120 76 8120 Al-27 11 1137 33 2137 55 4137 77 8137 Si 12 1140 34 2140 56 4140 78 8140 Ca 13 1200 35 2200 57 4200 79 8200 Cr 14 1240 36 2240 58 4240 80 8240 Mn-55 15 1255 37 2255 59 4255 81 8255 Fe 16 1260 38 2260 60 4260 82 8260 Ni 17 1280 39 2280 61 4280 83 8280 Zr 18 1400 40 2400 62 4400 84 8400 Mo 19 1420 41 2420 63 4420 85 8420 Pb 20 1820 42 2820 64 4820 86 8820 U-235 21 1925 43 2925 65 4925 87 8925 U-238 22 1928 44 2928 66 4928 88 8928

analytically based on the Klein-Nishina formula for each material region.

Angular quadrature

An angular quadrature set is fixed in the PALLAS code as shown in

Fig. 2. The number of angular mesh points and their distribution on a

unit hemisphere are approximately equivalent to the ordinary S,

quadra-ture in S,, code except for angular mesh points on both ends of z axis

designated as Q,,, Q,2, Q, and Q82 in the figure. These mesh points are

added especially for accurate calculation of radiation streaming through a duct whose axis coincides with z axis.

There are no constraint conditions on production of the PALLAS

angular quadrature because of no Legendre expansion of the angular flux.

Although the angular quadrature mesh points shown in Fig. 2 are

ap-proximately symmetric, any asymmetrical quadrature sets can be utilized.

The weights for the angular quadrature mesh points are equal to the

differential areas on a unit sphere. From this condition the dis-crete level of w,, may be chosen at the midpoint of the range Jw,,, while

the mesh point of ç,,,, may be also chosen at the midpoint of the range

Jpq.

Slowing down calculation into thermal group and convergence

Slowing down calculation from higher energies into thermal group is carried out through two ways: The first is a calculation of scattering due to hydrogen, in which the calculation is made in the same way as the usual

calculation described in the section 1.1 of III. While the other is a

cal-culation of scattering due to all nuclides other than hydrogen, in which

only the slowing down calculation from the energy region just above the

thermal group is taken into account based on the assumption that the

elastic scattering is isotropic in the laboratory system and inelastic scat-tering component is disregard for simplicity as follows;

Q(r, z, Q, E,,,) X,(E,h_l)IO(r, z, E,,,_I)JE,h_I

4U

4ir

e

(55)

where E,,, is the thermal energy, and ¿lU and ¿ are respectively the

leth-argy width and the average logarithmic energy decrement per collision

given by

(A-1)2

ln-1

2A A+1 (56)

The iterations are continued until the scalar flux converges according

Max{ çb';() }<EPSRN

çY(r)

for all spatial meshes.

IV. INFORMATIONS FOR USERS

1. Input Specification

Card i Title card PROBLEM (20A4) Card 2 Control Integers for a PROBLEM

KNDG, KIN, NORF, KTST, MONOE, MNODRE (613)

KNDG= 1, only neutron calculation. =4, only photon calculation.

=3, neutron and photon calculations.

KIN==0, no effect.

=n, n coupled calculations (4U=0.2, 0.4 and 0.8) (data note 1). NORF=0, reflection boundary condition at the bottom boundary.

>0, no reflection at the bottom boundary.

KTST=0, no effect.

= 1, check of input data read in and no transport calculation.

=3, test calculation up to KTST-th energy mesh. MONOE=0, no effect.

=10, monoenergy problem (data note 2).

MNODRE=0, no effect.

=n, monodirectional problem (data note 2).

Repeat n times from card 3 if KIN=nrO.

Card 3 JJ, IR, IZ, IUNCL, IFIS, TNPTP (613)

JJ= number of energy meshes.

IR=number of r-regions.

IZ=number of z-regions.

IUNCL=0, no effect.

=1, uncollided flux is separately computed with or without a First Collision Source (data note 3).

IFIS=0, source energy spectrum S(E) is read in from cards.

=1, fission spectrum is computed in the program. Do not

enter S(E).

INPTP=_{0, nuclear data for neutron are taken from PALLAS library.} =_{1, nuclear data for neutron are read in from cards.}

CARD 4 NBND, IBZ1, IBZ2, IBR, LTAP, JOAK (613) NBND=0, volume source.

= 1, point Source.

=10, top and/or bottom boundary flux problem.

=30, the same as NBND=20 and NORF=0 (data note 4).

=40, the same as NBND=20 in the case of utilizing tapes

in which specially prepared data are stored (data note 4). IBZ1=z mesh at which top boundary is defined.

IBZ2 = z mesh at which bottom boundary is defined.

IBR=r mesh at which cylindrical surface boundary is defined. LTAP=0, no effect.

>0, boundary fluxes are read in from Tape 22 in which a

previous calculation is stored.

JOAK = 0, no effect.

>0, print boundary angular fluxes up to JOAK-th energy mesh for check.

Card 5 EMAX, HH, SNORM, RDST (4E10.3) EMAX=maximum energy in MeV.

HH=lethargy interval for neutron, for photon HH=0.0.

SNORM = source normalization for source volume. f S(r,

RDST=0, no effect.

>0, first radial distance in cm if one wants to define a certain value (data note 5).

Card 6 MER(n), n=1, IR (2013)

Number of meshes in n-th radial region (data note 6).

Card 7 RR(n), n=1, IR (8E10.3)

Thickness in cm in n-th radial region.

Card 8 MEZ(n), n=1, IZ (2013)

Number of meshes in n-th axial region (data note 6).

Card 9 ZZ(n), n=l, IZ (8E10.3)

Thickness in cm in n-th axial region.

Card 10 E(j), j=1, JJ (8E10.3)

Photon energy at j-th energy mesh in MeV. Do not enter for neu-tron because neuneu-tron energy is specified by lethargy interval.

The next card is specified in accordance with the value of NBND.

For NBND=0, volume source option is chosen;

S(r, z. E)=S(r)*S(z)*S(E)n/cm4 sec MeV.

Card 11-A MSR, MSZ (213)

MSR = r-mesh up to which radial source distribution is read in. MSZ=z-mesh up to which axial source distribution is read in.

Card 11-B RZR, ZZR1, ZZR2 (3E10.3)

RZR=radial distance of source in cm.

ZZR1=top axial distance of cylindrical source. ZZR2=bottom axial distance of cylindrical source.

Card 11-C SR('n), n=1, MSR (8E10.3)

Radial distribution of source intensity at each energy mesh.

Card 11-D SZ(n), n=1, MSZ (8E10.3)

Axial distribution of source intensity at each energy mesh.

Card 11-E SE(j), j=1, JJ (8E10.3)

Source energy spectrum (/MeV) if IFIS=0.

Do not enter if IFIS>0.

For NBND=1, point source option is chosen.

Card 11-A MZSC, NSC, LIAG, MZIR (413)

MZSC=z-mesh at which a point source is set.

NSC=z region in which a point source is set.

LIAG= limit of angular meshes at which radiation is emitted (data

note 7).

MZLR=z-mesh of boundary between air and ground for skyshine

calculation. Card 11-B ZSC (E10.3)

Z distance at which a point source is set. Card 11-C SE(n), n= 1, JJ (8E10.3)

Source energy spectrum (/MeV).

For NBND=10, top and/or bottom boundary flux problem is chosen.

BOUN (f7, R)= S(R, Q)*S(E)n/cm1 sr sec MeV

for Q>0 with respect to z axis.

BNMZ (Q, R)= S(R, Q)*S(E)n/cm2 sr sec MeV

for Q<0 with respect to z axis.

LTAP=0

Card 11-A LRL, ISOC, ICONT, IXR, IYR (513)

LRL=r mesh up to which boundary fluxes are read in.

ISOC=0, angular distribution of boundary flux is read in. >0, isotropic angular distribution.

Do not enter angular distribution.

ICONT=0, radial distribution of boundary flux is read in. >0, constant radial distribution.

Do not enter radial distribution.

IXR and JYR; boundary fluxes at r meshes from IXR through IYR are printed for check.

Card 11-B RZR, ZZR1, ZZR2 (3E10.3)

If one wants to define the volume _{source used in a previous}

cal-culation, one may define its position. RZR=ZZR1=ZZR2=0, no

effect.

Card 11-C SE(j), j=1, JJ (8E 10.3) Source energy spectrum (/MeV).

Card 11-D-1 SN(m, ip), ip=ipl, ip2/m=1, LRL (8E10.3)

srsecMeV if ISOC=ICONT=0.

ipl=1 and ip2=IHIC for IBZ1>0, ipl=IHIC±1 and ip2=IQT for

IBZ2>0, where IHIC angular meshes belong to Q>0 with respect to z axis and IQT is the total number of angular meshes.

Card 11-D-2 SN(ip), ip=ipl, ip2 (8E10.3)

Angular flux with constant radial distribution

if ISOC=0 and

ICONT>0.

Card 11-D-3 SN(M), M=1, LRL (8E10.3)

Radial distribution with isotropic angular distribution if ISOC>0

and ICONT=0.

Note if ISOC>0 and ICONT>0 do not enter Card 11-D.

BOUN (Q, R) and/or BNMZ (Q, R) are calculated in the program as

B(Q, R)=1/47r.

For NBND =10 and LTAP >0,

Card 11-A LRR, LZZ, LZZ1, LZZ2, LRL, IXR, IYR (713)

LRR=total r mesh in a previous calculation. LZZ=total z mesh in a previous calculation.

LZZ1=z mesh in a previous calculation at which top boundary

condition is defined for a new calculation.

LZZ2=z mesh in a previous calculation at which bottom boundary condition is defined for a new calculation.

LRL= r mesh up to which top boundary values are defined. Card 11-B ROLD (i), i=1, LRR (8E10.3)

R distances in cm in a previous calculation.

Card 11-C RZR, ZZR1, ZZR2 (8E10.3)

RZR=ZZR1=ZZR2=0, no effect.

For NBND = 20, cylindrical surface flux problem is chosen.

Card 11-A LZL1, LZL2, IXZ, IYZ (413)

LZL1 and LZL2; boundary condition is defined at z meshes from

LZL2 through IZL1 (LZL1 LZL2).

IXZ and IYZ; boundary condition is printed at z meshes from IXZ

through IYZ.

For LTAP=0,

Card 11-B BOUNR (ip, M), ip=1, IQT (8E10.3)

Repeat M=LZL2, LZL1.

Angular fluxes at all the angular meshes at each z mesh in units

of n/cm sr sec MeV.

Repeat from the first through last energy mesh. For LTAP>0,

Card 11-B LRR, LZZ, LRRL (313)

LRR=total r meshes in a previous calculation. LZZ=total z meshes in a previous calculation.

condition is defined for a new calculation.

Card 11-C ZOLD(n), n=1, LZZ (2013)

Z distances in cm in a previous calculation.

If IBZ1>0 and/or IBZ2>0, top and/or bottom boundary _{fluxes are}

specified in the option of NBND=l0.

For NBND==30, LTAP>0 and NORF=0, top and/or bottom_boundary

fluxes and cylindrical surface boundary fluxes are taken from a previous

calculation, which is stored in Tape 22, with condition of NORF=0.

BOUN (Q, R) and/or BNMZ (Q, R) and BOUNR (Q, Z) are specified at

each energy mesh.

Card 11-A _{LRL, LZL1, LZL2, IXR, IYR, IXZ, IYZ (713)}

All these integers are the same control words as those used

previ-ously.

Card 11-B RZR, ZZR1, ZZR2 (3E10.3)

RZR=ZZR1=ZZR2=0, no effect.

Card 11-C _{LRR, LZZ, LRRL, LZZ1, LZZ2 (513)} All these data are the same as used previously.

Card 11-D _{NRLD(n), n=1, LRL (2013)}

R-mesh numbers in a previous calculation, which

_{are used for}

boundary mesh numbers in a new calculation.

Card 11-E _{NZLD(n), n=1, (IBZ1IBZ2+1) (2013)}

Z-mesh numbers in a previous calculation, which

_{are used for}

boundary mesh numbers in a new calculation.

Card 11-F _{NZND(n), n=1, (IBZ1IBZ2+1) (2013)}

Z-mesh numbers of boundary in a new calculation.

Card 12 LTHAL, LCIJT (213)

LTHAL=0, no effect.

>0, execute thermal group calculation (data note _8).

Lcut=0, no effect.

>0, termination of iterative thermal group calculations. Do not enter this card if KNDG=4.

Card 13 EPSRN (E10.3)

Convergence criterion. Do not enter this card if KNDG=4.

Card 14 NOEL (i, j), j=1, IR/i=1, IZ (2013)

Number of elements in each material region (data _{note 9).}

Card 15 NEK (i, i), i= 1, JR/i = 1, IZ (2013)

Numbers of material identification (data note 9). Card 16 ISOR (i,j), j=1, IR/i=1, IZ (2013)

Numbers of source identification (data note 10).

if JUNCL>0. Do not enter them if JUNCL=0.

JSOR(i, j) = 1, source region.

ISOR(í, j)=0, no source region.

Numbers indicating the use of a first collision source in uncollided

flux calculation in regions of JUNG (i, j)= 1. If IUNC(i, j)=0, no use

of a first collision source in (i, j) regions (data note 10). Do not enter them if IUNCL=0.

Card 18 NPHI (13)

Top and bottom faces of a cylindrical source are divided into NPHI

sections by equal azimuthal angle for calculation of uncollided fluxes above and below the cylindrical source by a first collision

source (data note 10). If IUNCL=0, do not enter it.

Input data for nuclear data.

Card 19 MATERIAL (6A4)

Name of material.

For neutron and in the case of INPTP=0,

Card 20 NUCLIDE (6A4)

Symbol of nuclide.

Card 21 MATNO, AMAS, AN (NUC) (15, 2E10.3)

MATNO=rmaterial number (4 digits given in Table 3).

AMAS = mass number of nuclide.

AN (NUC) = nuclear density (X 1024). For neutron and in the case of INPTP>0,

Card 20 NUCLIDE (6A4)

Card 21 AMAS, AN (NUC) (2E10.3) Card 22 SIGT(j), j=1, JJ (8E10.3)

a(E,), microscopic total cross section (barn).

Card 23 SIGMA(j), j=1, JJ (8E10.3)

o,(E), microscopic elastic scattering cross section (barn).

Card 24 LL, JLL (213)

LL=number of terms of Legendre expansion coefficient.

JLL= energy mesh up to which Legendre expansion coefficients are read in.

Card 25 FMU(j, i), 1= 1, LL/j= 1, JLL (8E10.3)

Legendre expansion coefficients for =1, LL at E,. Card 26 INEL (13)

INEL=0, inelastic scattering data are not read in. >0, inelastic scattering data are read in.

Card 27 JIN, J2N (213)

Inelastic scattering cross sections and/or (n, 2n) reaction cross sec-tions are not zeros up to JTIN-th and/or J2N-th energy mesh.

Card 28 CN(j, k), k= 1, J2N/j= 1, JJ (8E10.3)

Slowing down matrix from E, to E, energy mesh due to (n, 2n)

reaction.

Card 29 CIB(k,j), k=1, JTN/j=1, JJ (8E10.3)

scattering.

For photon,

Card 20 AN (NUC) (E10.3)

Electron density (X 1024).

Card 21 CRT(j), j=1, JJ (8E10.3)

Linear attenuation coefficients (cm') for material. Card 22 SIGMP (j). 1=1, JJ (8E10.3)

Linear pair production coefficients (cm-) for material. Input data specifying external data and output

Card-0 i ITP2O, ITP21, KANK, ISKIP, ITP14, ITP19 (613) 1TP20=0, no effect.

= J, calculation starts using the data calculated previously and stored in Tape 20.

ITP21=0, no effect.

= 1, calculated angular fluxes are stored in Tape 21. KANK=0, no effect.

=1 or 2, if ITP2O=1 (data note 11). ISIP=0, no effect.

=k, skip k groups from the first group when reading the data

from Tape 20.

ITP14=0, no effect.

= 1, calculated reaction rates (or dose rates) are stored in Tape

14.

ITP19=0, no effect.

= 1, calculated scalar fluxes are stored in Tape 19.

Card-0 2 _{MRK, MZK, MDS, MZDS, JEF, NOTR (613)}

MRK=number of r meshes at which angular fluxes are printed. MZK=number of z meshes at which angular fluxes are printed.

MDS=0, no effect.

=n, n reaction rates are calculated.

MZDS=0, no effect.

>0, currents at MZDS-th (z) mesh for +z direction and

z

direction are printed at all radial meshes.

IEF=0, angular and scalar fluxes are printed in units of n/cmsec (sr) MeV.

=1, energy X angular and scalar fluxes are printed in units of MeV/cm2 sec (sr) MeV.

NOTR=0, no effect.

>0, responses (reaction rates or dose rates) are printed at

each energy mesh. Usually they are printed after all

com-putations.

Card-0 3 KR(n), n=1, MRK (2013)

If MRK=0, no angular fluxes are printed.

Card-0 4 KZ(n), n = 1, MZK (2013)

Axial mesh numbers at which angular fluxes are printed if MZK>0.

Card-0 5 DOSE(j, n), j=1, JJ/n=1, MDS (8E10.3)

Reaction constants (or dose rate conversion factors) at energy mesh

E.

If MDS=0, do not enter them.

2.

Detailed data notes

(1) Data note i

The PALLAS neutron transport calculations are usually made based

on the nuclear data from the PALLAS library, in which the energy mesh

structures are defined depending upon the lethargy width as shown in

Table 1. The structure of 0.1 lethargy width is generally used for

neu-tron transport calculations by one-dimensional PALLAS code'. For energy

structures used in PALLAS two-dimensional neutron calculations, the

structure of 0.2 lethargy width is recommended for the energy range of

MeV because of a steep gradient in the shape of energy spectrum

com-pared with i/E. Then the structure of 0.4 lethargy width is used for

calculation in the intermediate energy range and the 0.8 lethargy-width

structure is used in the low energy range including thermal group. Other

combinations of structures can be utilized, for instance the 0.4+0.8

lethargy-width structures or only the 0.8 lethargy-width structure. It

should be desirable from the viewpoint of accuracy of calculated result to use the combination of 0.2+0.4+0.8 lethargy-width structures. As an

example we present the accuracy of a neutron transport calculation in

water. The neutron fluxes in the intermediate and low energy regions

including thermal group calculated with the 0.8 lethargy-width structure

of 23 energy meshes are underestimated by a factor of approximately 2

compared with those by the combination of the 0.2+0.4+0.8 lethargy-width

structures of total 35 energy meshes.

Group 1 2 3 4 5 6 7 8 9 0.2 cthargy-width structure

No. i i i i i i t i I (KIN = 2) ¿ ' 't' 0.4 lethargy-Group 1 2 3 4 5 6 7 6 9 10 'i width structure No. J I t I i i L i i i i (KIN = 1) Group 1 2 3 4 5 6 7 8 No . I I I i I I i 0.8 lethargy-width 20 21 22 23 structure (KIN = 0)

Fig. 7. Examples of series calculations by combinations of energy structures of 0.2, 0.4 and 0.8 lethargy widths.

O z A r-'Reflection boundary candi tion Previous calculation O A New calculation

Fig. 8. Illustration of calculation with NBND=30.

r

When one uses the series calculation of lethargy-width _structures,

one should input an integer of i or 2 into the KIN; in the

_{case of}

KIN=2, the structures of 0.2+0.4+0.8 lethargy widths are used, while

in the case of KIN=1, the structures of

_{0.2+0.4 or 0.4+0.8 lethargy} widths are used. _{Concrete examples of combinations in energy-mesh}

structure are shown in Fig.

7.

As seen in the figure the calculated

angular fluxes at energy meshes of odd numbers _{are picked up for use}

in the next calculation in a series calculation.

Data note 2

Since the PALLAS calculation is made based on the continuous source

energy distribution for the purpose of avoiding iterative calculations, it

is inevitable to use a special treatment in the calculation of

discrete-energy source problems. Thus in the case of a mono-energy source prob-lem one should input as MONOE>0.

In the case of a problem of monodirectional perpendicular _{or slant}

incidence on the plane at the first mesh of z axis, one may use the input

option of MNODRE=n, where n is the mesh number of w (the discrete

value of cosine of incident polar angle). _{The monodirectional uncollided}

flux is calculated analytically for a direction of w,,.

Data note 3

The input IUNCL=1 may be used when

one wants to calculate

analytically the uncollided flux for mitigating _{ray-streaming effects in}

such a problem as a deep penetration _{or a transport problem in a very}

large air medium. Data note 4

One may use the input option of NBND=30 when one wants to

per-form a calculation using the top, bottom and cylindrical surface boundary

fluxes taken from the data calculated previously with the reflection boundary condition at the bottom as shown in Fig. 8. On the other

hand for input option of NBND = 40, the top, bottom and cylindrical

surface boundary fluxes have been taken from a previous calculation and

stored in Tape 14, 18 and 17, respectively. Then a new calculation can

be made by using these tapes. Data note 5

When RDST=0, the first radial distance is automatically defined by

dividing equally the thickness of the first radial region. If one wants to

> r

Fig. 10. Angular mesh numbers of polar angle. 1 2 3&4 5 C 748 9 10 C 1/

Fig. 9. Spatial mesh assignment.

14 H 11 10 b -3

NE K= z

Fig. il. An example of assignment of identification numbers of material regions in the case of IR=1O.

r

define a certain value in the first radial distance, one may input the value

into RDST.

Data note 6

As shown in Fig. 9, two meshes are assigned at inner boundaries.

However, no mesh assignment is at the r=O cm as seen in the figure. Data note 7

There are 8 discrete levels on polar angle in the current PALLAS code, so that LIAG=8 means the isotropic point source and LIAG<8

means that radiations are emmitted into the first LIAG directions shown

in Fig. lo.

Data note 8

The last energy meshes in the structures of 0.4 and 0.8 lethargy widths

may be assumed to be thermal group. Then thermal neutrons can be

calculated with the conventional iteration-convergence technique applied in group calculation. For this thermal group calculation, one may specify

the inputs as LTHAL>O, LCUT>O and a certain value of convergence

criterion into EPSRN. Data note 9

Mixture cross sections for neutron are not utilized in the current PALLAS code. Instead, cross sections for each nuclide are read in and

used for calculation in each material region. _{For this purpose the}

num-ber of elements in each material region NOEL(i, j) must be specified by

WATER 3 WtTER 3 CONCRETE 33 IRON 4 WATER 3 WATER 3 CONCRETE 33 IRON 4 REF. 2 REF. 2 WATER 3 IRON 4 CORE 1 REFLECTOR 2 WATER 3 IRON 4

input. Although the number of all elements included in a material region should be specified, it should be recommended from the viewpoint of

ef-fective shortening the computer time that those nuclides which have

small values of macroscopic cross sections (nuclear densities X

micro-scopic total cross sections) are omitted and their nuclear densities are

added to those of other nuclides.

The identification numbers of material regions, NEK(i,j), are assumed

at first to be integers calculated by j+(il)XIR. When the same

ma-terials are repeatedly specified in material regions, their identification

numbers are replaced by those of the previous same materials as shown in Fig. 11.

(10) Data note 10

In order to define the (r, z)-positions of surface of cylindrical source,

integer of 1 or O is specified in each material region.

In the case of a

disk source, ISOR(i,j)=l should be specified in the material region whose

bottom boundary coincides with the disk source. Further in the case of

a line source, ISOR(i,j)=l should be specified in the material region whose first radial mesh (which may be considered as cylindrical axis)

coincides with the line source.

The uncollided angular flux can be analytically calculated at spatial

mesh points in all material regions outside source regions by using a

first collision source. However, restriction should be considered on the

material regions in which the first collision source is applied, because its

use results in very much increase in computer time. From the viewpoint

of practical calculation it is recommended that its use is restricted to those

material regions which are positioned at large distances from source regions.

Since the angular flux is only one value with respect to radial mesh

z

Fig. 12. Division of top and bottom surfaces by an equal

Table 4. PALLAS-2DCY-FC file requirements Contents

PALLAS library

Scrach

Calculated data for RS, zS, PSYS

Input

Output

Calculated data WT and

GZI

9 _{Calculated data for}

FGM

Angular fluxes FN

Top boundary fluxes,

BOUN

Calculated data for CIA

WT and GZI for only hydrogen

Reaction rates or dose

rates IJncollided angular fluxes Cylindrical surface boundary fluxes, BOUNR

Bottom boundary fluxes,

BNMZ

Scalai fluxes

20 Angular fluxes calculated previously

21 _{Calculated angular}

fluxes for all groups

Remarks

Neutron cross sections are read

from PALLAS library. For working.

Data for RS, ZS, PSYS are

calcu-lated in Subroutine MISIMA and

used in KYOTO.

Data for WT and GZI are calcu-lated in Sub. ATAMI and used in

NAGOYA.

Data for FGM are calculated in

Sub. YOKHAM and used in

NA-GOYA.

FN for each energy mesh are stored

temporally in this unit. Large storage is necessary, for instance (55Y55 28).

Data are stored in Sub. TOKYO

and used in MAIN.

Inelastic scattering matrices CIA are calculated in Sub. YOKHAM

and used in NAGOYA.

Data are calculated in Sub. ATAMI and used in NAGOYA.

Data are stored in Sub. OSAKA.

Data are stored in Sub. TOKYO

and used in MAIN.

Data are stored in Sub. TOKYO

and used in MAIN.

Data are stored in Sub. OSAKA. This file is used only when input ITP2O>0 (data note 11).

This file is used only when input

ITP21=l.

at z=z, in (r, z) geometry, then top and bottom surfaces of _{a cylindrical}

source or a surface of disk source must be divided azimuthally into _NPHI

sections for calculation of the uncollided angular fluxes at upper and/or lower positions of source as shown in Fig. 12. The values of NPHI are

recommended to be 10 for _{nearer positions to the top or bottom surface}

of cylindrical source and 6-8 for distant positions.

(li) Data note 11

In the case of ITP2O=J, KANK=l and ISKIP=O, the angular fluxes

are read from Tape 20 one by one for all energy _{groups up to J-th group}

lo 11 12 13 14 15 17 18 19 Logical unit i 2 3 5 6 8

MAIN

Input data for

Title and

con-trol integers for a problem >

Input data for specifying external data and output Only KTSTO Re ad (20) Only KNDG J= 1 IUNCL ) O JHII< O -- KNDG < 3 KNDG 3 IUNCL = O or IUNCL 70 and IUNC (N1, N2) =0 Only JHIK 70 and IUNCL >0 Store uncllided flux - 15 JHIK = O KNDG 3 and LTHAL 7 0 and no convergence - or MAIS < LCUT J <JJ J=J+l ) KIN > O KIN = KIN - i TOKYO

Read in input data for program and calculate or prepare boundary

data of BOIJN,BNMZ,BOUNR and energy

mesh and spatial mesh assignments.

1'

Y 0KM AM

Read in input data for nucear data and calculate the data of CRT, CIA,

SIGMA and FON.

ATAMI MISIMA

Il'

READ

Read neutron cross sections from PALLAS

library.

ZS and PSYS are calculated. GIFU

Determination of number of

meshes of _pq

PSt.)

_,Tection of spatial meshes belong to (N1,N2) region.

NAGOYA

Calculation of neutron scattering source, SN (rizkOpq)

HIMEJI

Calculation of photon scattering source,

KYOTO

Radiation transport calculation, FN(rizkOpq)

UN CLF X

Uncollided angular fluxes are calculated by a first collision source.

OSAKA

sPrint anoular fluxes and scalar fluxes

and calculate reaction rates. If J = JJ,

ITP14 >0 reaction rates are stored in Tape

14 or 24; ITpl9 O scaler fluxes are stored in Tape 19 or 29.

KIN = O and

ITP21 >0

Calculated angular fluxes are stored in Tape 21.

Fig. 13. Simplified logical flow diagram with subroutine.

Program mnemonic

WP (IP) WWP (IP) WBP (IP)

WPQ (IPQ)

PSY (IP, IPQ) PSYB (IP, IPQ)

AMTJ (M) WAMU (M) DR (I) DZ (I) CRT (N,, N2, J) RD (M) FN (MR, MZ, IPQ SN (MR, MZ, IPQ) ALPH (M) ALHY (M) FGM (M, J) KMAX (M) KK (M, NONII) SIGMA (J, NUC) WT (N, IP, J) GZI (N, IP, J) RO (NTJC) RS (M, IP, IPQ) ZS (M, IP, IPQ) PSYS (M, IP, IPQ)

A (N1, N2) ION >0 <0 MUON >0 <0 MRR MZZ CIA (JMK, NONH) GZAI (NUC) WAVE (J) DWAVE (J) GMU ZD (K) ORGN MRB (N2) MZB (N1) RB (N2) ZB (N1) EDN (N1, N2) IQT

Program variable and remark

Wj)

4w,,; weight for Wp Boundary value of 4w,,

JQ,,,,; weight for angular mesh point Boundary value of

Il,,,

Ji,,, ; weight for p,,

4v; radial interval in I-th region 4z; axial interval in I-tb region

X,(r, z, E,); macroscopic total cross section

in N1z and N8r region

r,; radial distance (cm) ¡(r,, Zk, Qpq)

Q'(r,, Zk, Qpq) Ir,,,

er,,, for hydrogen

.Y,(r, Ei)f(E, ) Maximum E,.

E5(1i,,,) for nuclide nimber NONH ,,(r, E,) Wp, n.,,,)ç6(w', wi,, a,,,) ( a,,,wPw,, COS p=l/A rj.1 Z Q=Q(w, çl) A(6,E) (±r, +z) direction (±r, z) direction (r, ±z) direction (+r, ±z) direction

Total number of r meshes Total number of z meshes

.,,,(r, Ek'Ej) .4 U 4x ç A1 4A [1 = i + A' A Zk; axial distance (cm) (ZZR1 + ZZR2)/2

Boundary mesh numbers in

N2r

region

Boundary mesh numbers in N1z region Radial boundary distance

Axial boundary distance

nC); electron density in (N1z, N2r) region

Total angular meshes

NEUTRON CALCULATION SAMPLE PROBLEM 1 02

00

0.0

42111011

22111011

22111011

_11111111

12343038

22343080

22333038

(Input data for external data and output)

o o o

o ii

1 1 1

0 00

31 1

1.500E-01 1.490E-01 1.400E-01 1.480E-01 1.470E-01 1.460E-01 1.455E-01 1.450E-01

1.450E-01

Fig. 14. List of input data of sample problem 1.

loo

981

000

14.20 10 3 2 36.20 1217 4 44.0 lo 12 36.20 2.967E 02 2.967E 02 1.00 1.00

200

010

000

0.20 6 2 2 7.0 9 60.7 44.0 2.967E 2.967E 1.00 1.00 5 02 02 0.0 2 1.250 10.2 0.0 2.967E 1.00 1.00 02 0.0 16.50 31.8 2.967E 1.00 1.00 02 1.25 2.967E 1.00 02 2.60 2.967E 1.00 02 8.90 2.967E 1.00 02 14.0 2.967E 1.00 3 3 3 3 CORE HY 3 3 3 3

(Input data for nuclear data) 2011 1.000E 00 4.716E-02 Oxy 2086 1. 600E 01 2.977E-02 FE 2260 5.585E 01 1.159E-02 U-238 2928 2.380E 02 3.110E-03 WATER HY 2011 1.000E 00 6.680E-02 Oxy 2086 1. 600E 01 3.340E-02 IRON FE 2260 5. 585E 01 8.479E-02 LEAD PB 2820 2.080E 02 3.297E-02 AIR N-14 2074 1.400E 01 4.000E-05

GAMMA-RAY CALCULATION SAMPLE PROBLEM 2

(Input data for external data and output)

000011

121 000

33

13

3.800E-03 3.400E-03 3. 100E-03 2.740E-03 2.320E-03 1.780E-03 Fig. 15. List of input data of sample problem 2.

i i i i 1

0 il

i i i i i

0 11

1 1 i 1 i

0 11

1 1 1 1 1

0 ii

1 ii 111 ii

1 2 3 4 3

0 38

2 2 3 4 3

0 38

2 2 3 4 3

0 38

2 2 3 3 3

0 38

3 3 3 3 3

3 33

(Input data for nuclear data)

CORE 0.887

0. 1252 0. 1306 0. 1379 0.1482 0.1645 0.2013

2.435E-02 2.000E-02 1.352E-02 7.971E-03 2.236E-03 0.0 WATER

0. 3345

3.818E-02 4. 100E-02 4.496E-02 4.958E-02 5.577E-02 6.608E-02

0.0 0.0 0.0 0.0 0.0 0.0

IRON

2.203

0.2787 0.2930 0.3116 0.3364 0.3726 0.4365

3.566E-02 2.800E-02 1.783E-02 9.530E-03 2.450E-03 0.0 LEAD 2.693 04736 0.4800 0.4926 0.5182 0.5660 0.7122 0. 1484 0. 1180 0.0833 0.05001 0.0136 AIR 3.345E-04

3.818E-05 4.100E-05 4.496E-05 4.958E-05 5.577E-05 4.608E-05

0.0 0.0 0.0 0.0 0.0 0.0

400

685

000

3.250 8 3 2 36.20 12 710 44.00 3.250 8 12 36.20 1.00 1.00 1.00 5.187E 04

200

000

000

0.0 9 2 2 7.00 4 9 18.0 2.80 44.0 1.00 1.00 1.00 1.741E 0.0 52 1.250 42.696 2.40 0.0 1.00 1.00 1.00 06 2.408E 07 0.0 16.50 10.20 2.00 1.00 1.00 1.00 1.907E 08 1.250 31.80 1.60 1.00 1.00 4.790E 08 2.60 1.150 1.00 1.00 1.520E 09 8.90 1.00 1.00 14.0 1.00 1.00