Graziani F. (editor) Computational Methods in Transport

Подождите немного. Документ загружается.

522 E.W. Larsen

In spite of these random aspects of individual particle histories, (1) contain

no randomness. That is, if these equations are solved, one obtains a result

that seems to carry no information about what a single neutron might do.

How can such a result be meaningful?

The answer is that (1) describe the mean neutron ﬂux, averaged over

an inﬁnite number of neutron histories. Typical neutron transport problems

often involve well over 10

neutrons, so the random ﬂuctuations from the

mean – which do exist – are usually small. If one were to simulate a very large

number of individual neutron histories and average the results, an estimate

of the neutron ﬂux would be obtained that agrees with the solution of (1).

The principle underlying Monte Carlo simulations is to ignore the mathe-

matical description of the problem (1), and to directly simulate a large number

of particle histories, just as discussed above. The goal is to simulate and av-

erage a suﬃciently large number of particle histories that a useful estimate

of the ﬂux is obtained.

The principle underlying deterministic simulations is to ignore the ran-

dom aspects of individual particle histories and solve (1), by (i) discretizing

these equations with respect to each of its variables – thereby converting the

equations into a (typically) large system of algebraic equations – and (ii)

solving this algebraic system of equations.

Thus, although Monte Carlo and deterministic methods solve the same

physical problem, they are based on fundamentally diﬀerent principles of

numerical computation.

Next, we discuss in somewhat greater detail (i) the ways in which Monte

Carlo and deterministic methods are implemented, and (ii) the comparative

advantages and disadvantages of these two appoaches. We begin with a dis-

cussion of Monte Carlo methods.

4 Stochastic (Monte Carlo) Methods

Neutron and photon Monte Carlo methods are based directly on the physics

of neutron and photon transport, as experienced by individual neutrons and

photons. In the Monte Carlo process, individual particle histories are sim-

ulated (and averaged together); one does not need to work with the linear

Boltzmann (1a). The details of Monte Carlo particle transport simulations

are well-described in documents that are devoted to the subject ([4, 5, 9, 12],

and [15]; see also [18] for shielding applications), and we will only touch on

a few main points here.

The ﬁrst important issue to discuss is the fact that, in principle, Monte

Carlo simulations include no errors in the modeling of the physical system or

the treatment of the transport physics. If one perfectly knows the geometry

of the system and the the cross sections at each point within the system, then

the analog Monte Carlo process (in which particle histories are simulated in

a manner that is completely faithful to the actual physics) will, in the limit

Neutron Transport Problems and Simulation Techniques 523

of an inﬁnite number of Monte Carlo particles, produce an exact estimate of

the ﬂux. This fact is one of the most compelling reasons for employing Monte

Carlo methods.

However, actual Monte Carlo simulations do not produce exact solutions

because no present or future computer will ever be able to process an inﬁ-

nite number of particle histories. Realistic Monte Carlo calculations require

that the simulation terminate after a ﬁnite number of histories have been

processed. Thus, all Monte Carlo simulations will inherently contain statis-

tical errors, due to the fact that the simulation ends after N histories have

been processed. Can one say anything about the magnitude of these statisti-

cal errors?

Fortunately, the answer to this question is “Yes.” To explain, let us con-

sider a typical problem in which a system V , a neutron source within V ,and

a detector region R within V are all prescribed (see Fig. 5).

Source

Detector

Fig. 5. Source-Detector Problem

We wish to calculate:

P = the probability that a neutron, born in the

source region, will be captured within R. (14)

Analog Monte Carlo simulations of this problem will generate and process

the histories of a speciﬁed number N of Monte Carlo particles. These particles

are all born within the source region, in a manner that is consistent with the

speciﬁed source. Each particle will stream through V , undergoing collisions,

etc., just as described above. Each history will end with the particle either

being captured at some point within V , or leaking out of V . When the n

particle history ends, we assign a tally:

1 if the particle is captured in R,

0 if the particle is not captured in R.

(15)

524 E.W. Larsen

After the N histories are completed, the probability P is estimated by the

average:



n=1

. (16)

The Central Limit Theorem provides a relationship between the exact

probability P and the estimate P

. This relationship holds for any Monte

Carlo simulation in which tallys are averaged, as they are in (16). The Central

Limit Theorem states that for each problem, there exists a constant σ (called

the standard deviation) such that, for each a>0,

lim

N→∞

prob



|P − P

σ/

√





1/2



−t

dt . (17)

This says that the statistical errors in Monte Carlo estimates will on the

average decrease as the number of particle histories increase, and that the

errors will be on the order of σ/

√

|P − P

| = O



√



. (18)

This result is both good and bad. It is good because it shows that, indeed,

as N →∞, Monte Carlo estimate P

of P will converge to the correct result.

It is bad because it shows that the convergence of P

to P is (i) slow, and (ii)

not guaranteed to be (and almost never is) monotonic. For instance, if one

has an estimate P

and one wants to generate a more accurate estimate with

an error which is smaller by a factor of 10, one should rerun the Monte Carlo

simulation with 100 times the number of particles as the ﬁrst calculation.

Unfortunately, the Central Limit Theorem does not guarantee that for a

single simulation the error will reduce by precisely a factor of 10; it only

says that if one were to run a very large number of these more expensive

simulations, the average error will reduce by a factor of 10 over the average

error in the less expensive simulation.

The O(1/

√

N) behavior of the error is a fundamental feature of Monte

Carlo simulations. Undeniably, this is a drawback to Monte Carlo simulations;

they converge slowly.

Another important aspect of Monte Carlo simulations is that the standard

deviation σ in (18) can be large. This is the case if the detector region in

Fig. 5 is distant from the source region, in which case P is very small (a

source neutron ending its history by being captured in the detector region

is a rare event). In this case, if the physical source emits 10

neutrons and

only 10

of these end their history by being absorbed in the detector, then

P =10

−7

, and on the average, only one out of a million Monte Carlo particles

will score in R. For simulations of rare events, it is almost never advisable

to use analog Monte Carlo, in which the physics is reproduces faithfully; the

resulting simulation is much too ineﬃcient.

Neutron Transport Problems and Simulation Techniques 525

For simulations of rare events – which occur inevitably in deep shielding

calculations – it is necessary to employ nonanalog techniques. In doing this,

the “rules of the game” are modiﬁed in such a way that the same value of P

is estimated, but the value of σ is reduced.

The most basic technique for doing this is absorption weighting.Here

particles are assigned an extra variable: w =particle weight. Particles are

generally born with w = 1. The same rules concerning the calculation of

distance-to-collision apply to these particles. However, when a particle un-

dergoes a collision, the rules change. If p

is the probability that the neutron

will be captured:

+ Σ

, (19)

then when a collision occurs, (i) we do not allow the particle to be captured

(we do permit all other scattering and ﬁssion processes to occur with their

relative probabilities), and (ii) we reduce the particle’s weight by the multi-

plicative factor p

. (By doing this, we greatly extend the lifetime of Monte

Carlo particles, ensuring that many more of them will eventually travel into

the detector.) If a Monte Carlo particle with weight w enters the detector

region R,weﬁxw and use analog Monte Carlo until the particle is captured

or leaks out of R. If the particle does leak out of R, we turn the nonanalog

process back on. Thus, Monte Carlo particles will end their histories in only

two ways: they can leak out of the system, or they can be captured within

In this procedure, when a particle ends its history, we assign:

ﬁnal value of w if the particle is captured in R,

0 if the particle leaks out of V .

(20)

After N particles are generated and their histories are processed, P

again calculated by (16). It is clear that the calculational cost per particle is

greater with absorption weighting than with analog Monte Carlo, because the

absorption-weighted Monte Carlo particles will have longer histories. How-

ever, many more absorption-weighted particles will score in the detector re-

gion, so the absorption-weighted statistical errors should be smaller. Given

the same amount of computing time, which method (analog or absorption-

weighting) will generally produce the more accurate answer?

For rare-event problems, absorption-weighted Monte Carlo simulations

are usually more eﬃcient than analog simulations. Also, in the limit N →

∞, the absorption-weighted estimate of P

will, like the analog estimate,

converge to the exact value of P . (Although absorption-weighted Monte Carlo

is not a faithful representation of the physics, it is a fair game.)

Many more sophisticated nonanalog Monte Carlo techniques have been

developed, all of which make use of the notion of particle weight, and all

of which are designed to make the Monte Carlo simulation run more eﬃ-

ciently. (In diﬃcult problems, Monte Carlo with absorption weighting is still

526 E.W. Larsen

unacceptably ineﬃcient.) One of the most widely-used of these nonanalog

techniques employs the use of weight windows. This procedure is based on

the observation that it can be helpful to control the range of particle weights

at each point in V . The weight window technique operates as follows:

1. The system V is divided into a subregions, each no more than a few mean

free paths across.

2. In each (j

) subregion, a lower weight w

−

and an upper weight w

are

assigned, with w

−

. (Usually, w

≈ 10w

−

3. The Monte Carlo process with absorption weighting is started. In addition

to the rules of absorption weighting described above, the following rules

are also imposed:

4. If a particle has weight w which is greater than the upper weight window

at that point, the particle is split into n particles, each with weight w/n

lying as close as possible to the center of the weight window. Each of these

n particles must then be tracked, and their histories run to completion.

5. If a particle has weight w which is lower than the lower weight window at

that point, the particle is subjected to Russian roulette:ifw

is the “cen-

ter” of the weight window, then with probability p = w/w

the particle

survives with increased weight w

, and with probability 1 −p the particle

is killed.

The eﬀect of the weight window is to place rough bounds on the numbers

of Monte Carlo particles at diﬀerent locations across V .Ithasbeenshown

that if the upper and lower weight window weights are chosen optimally, the

resulting Monte Carlo simulation becomes much more eﬃcient than absorp-

tion weighting. Other nonanalog strategies, such as geometric splitting,are

based on related ideas.

The use of nonanalog Monte Carlo techniques, such as weight windows,

for practical deep penetration problems in nuclear engineering, places a sig-

niﬁcant burden on the code user. The user must specify all the geometric

details of the physical system V and the necessary neutron sources. In ad-

dition, the user must supply all the biasing parameters – for example, the

numerical values of each of the weight windows (assuming that the weight

window technique is used). If V is large, with many subregions, and if the

weight windows also depend on energy, then there can be many hundreds of

these parameters that must be chosen. Determining the biasing parameters

in a way that leads to near-optimal eﬃciency of the Monte Carlo simulation

can be a diﬃcult and time-consuming task, relying on the experience, skill,

and luck of the code user.

Furthermore, the biasing parameters depend strongly on the location and

nature of the detector; if the detector changes, the biasing parameters also

change. This reveals another fact about Monte Carlo simulations for diﬃcult

nuclear engineering problems requiring nonanalog biasing: they are not user-

friendly.

Neutron Transport Problems and Simulation Techniques 527

To conclude this section, we brieﬂy summarize the advantages and dis-

advantages of Monte Carlo methods for neutron transport problems. First,

Monte Carlo methods for neutron and photon transport are, at their root,

quite simple – they directly apply the physics of particle transport to the

simulation of individual particle histories. This can be done in a way that

requires no approximations – each particle can be simulated in a manner

that perfectly adheres to the underlying physics. Therefore, Monte Carlo

simulations are capable of producing solutions which are free of any type of

truncation error.

However, Monte Carlo simulations of diﬃcult problems are often very

costly to set up and run. To make the Monte Carlo code run with accept-

able eﬃciency, the code user must specify a large number of biasing para-

meters, which are specialized to each diﬀerent problem. Determining these

parameters can be diﬃcult and time-consuming. Also, even when the bias-

ing parameters are well-chosen, Monte Carlo simulations converge slowly and

non-monotonically with increasing run time. Thus, while Monte Carlo solu-

tions are free of truncation errors, they are certainly not free of statistical

errors, and it can be challenging to obtain Monte Carlo solutions with suﬃ-

ciently small statistical errors, and with acceptable cost. Finally, the nonana-

log techniques that have been developed for making Monte Carlo simulations

acceptably eﬃcient are useful for source-detector problems – in which a de-

tector response in a small portion of phase space is desired – but are not

useful for obtaining eﬃcient global solutions, over all of phase space. Gener-

ally, Monte Carlo solutions work best when very limited information about

the ﬂux (e.g. a single detector response) is desired in a given simulation.

5 Deterministic Methods

As we discussed earlier, deterministic methods for solving neutron transport

problems are based on solving (1). The ﬁrst step is to discretize these equa-

tions with respect to all of the independent variables. We now brieﬂy discuss

each of these discretizations.

First, we discretize the energy variable E. To do this, we conﬁne the

energies of the particles to a ﬁnite interval E

min

<E<E

max

and divide this

interval into G energy groups as depicted in Fig. 6.

max

G-1

g-1

energy

group

min

Fig. 6. Multigroup Energy Grid

528 E.W. Larsen

Next, we approximate all the cross sections in (1a) by histograms in energy,

each histogram having one value within each energy group. (There are many

prescriptions for doing this, but however it is done, it is done carefully, in order

to minimize the total number of groups G.) Finally, the resulting equation is

integrated over each group. This yields the multigroup transport equations:

Ω

·∇ψ

(r,Ω)+Σ

t,g

(r)ψ

(r,Ω)





4π

s,g



→g

(r,Ω



· Ω)ψ



(r,Ω



)dΩ



(r,E)

4π





4π

νΣ

f,g



(r)ψ



(r,Ω



)dΩ



(r)

4π

∈ V,Ω∈ 4π, 1 <g<G, (21)

where

(r,Ω)=



g−1

ψ(r,Ω,E)dE (22)

are the group ﬂuxes. This multigroup approximation is used almost univer-

sally to treat the energy variable.

The angular variable is usually discretized in one of two ways. The more

common S

(discrete ordinates) approximation assumes that, instead of being

able to travel in an arbitrary direction of ﬂight Ω

on the unit sphere, neutrons

are only able to travel in a ﬁnite number M

of speciﬁed directions, Ω

,for

1 ≤ n ≤ M

. Furthermore, each discrete ordinate Ω

is associated with a

cone on the unit sphere of area w

. The set {(Ω

)|1 ≤ n ≤ M

} is called

an angular quadrature set of order N.(In1-D,M

= N , but for historical

reasons, in 2-D and 3-D, M

>N.) With this approximation, (21) becomes:

Ω

·∇ψ

g,n

(r)+Σ

t,g

(r)ψ

g,n

(r)



s,g



→g

(r,Ω



· Ω

)ψ



(r)w



(r,E)

4π



νΣ

f,g



(r)ψ



(r)w



(r)

4π

∈ V,1 ≤ n ≤ M

, 1 <g<G. (23)

An alternate P

(spherical harmonics) method is available, in which the

group angular ﬂuxes ψ

(r,Ω) are expanded in a ﬁnite number of spherical

harmonic functions

(r,Ω) ≈



n=0



m=−n

g,n,m

(r)Y

n,m

(Ω) , (24)

Neutron Transport Problems and Simulation Techniques 529

and a coupled system of equations for the expansion coeﬃcients ψ

g,n,m

(r)are

obtained by multiplying (21) by the spherical harmonic functions Y

n,m

(Ω)

and integrating Ω

over the unit sphere. The S

approach is more widely-used

because the system of S

equations closely resembles the original Boltzmann

equation; the comparable system of P

equations has a very diﬀerent and

more complicated mathematical structure.

Finally, to discretize the spatial variable, one ﬁrst takes the physical sys-

tem and imposes on it a spatial grid. Then, one discretizes each of (23) on

this grid. Many diﬀerent spatial discretization techniques have been devel-

oped, among them ﬁnite diﬀerence, ﬁnite element,andnodal methods. It is

not possible to discuss these methods here in any detail. Generally, all the

methods work well for problems on which the spatial cells are optically thin.

However, as the spatial cells become more coarse, the simpler ﬁnite-diﬀerence

methods will generally exhibit undesirable behavior – they can become nega-

tive, or oscillatory. While the coarsening of the spatial mesh will degrade any

discretization method, some methods are much more sensitive than others.

Certain ﬁnite element methods are quite insensitive to the use of coarse spa-

tial grids for optically thick problems dominated by scattering. Much eﬀort

has gone into the development of spatial discretization methods that are as

accurate as possible on coarse grids [16].

This concludes the ﬁrst part of a deterministic method – the discretization

of the linear Boltzmann equation. One now usually has a very large system

of algebraic equations to be solved, of the form:

AΨ = B, (25)

where A is a matrix representing the linear Boltzmann operator, Ψ is the

vector of unknowns, and B is a vector of interior and boundary sources. How

large is this system? In a typical 3-D shielding calculation, G = 150 energy

groups is often needed to achieve suﬃcient resolution of the energy variable,

= 24 discrete ordinates (3 ordinates per octant) is the minimum number

needed to resolve the angular variable, and a spatial grid of 100 × 100 × 100

cells is not unreasonable (and is often a low estimate). This implies 3.6 ×10

unknowns – a huge number, which is often much smaller than is seen in

some problems. (The large number of unknowns is a direct and inevitable

consequence of the high-dimensionality of phase space.) Generally, there is

little hope to solve (25) by directly inverting the matrix A. Thus, the second

part of a deterministic method – developing an eﬃcient solution strategy for

solving (25) – comes to the forefront.

Since it is generally impractical to directly invert A, solution strategies for

solving (25) have focused on iterative methods. A great deal of work has been

done on this topic and is described in a recent comprehensive review article

[22]. It is impossible to discuss the ﬁne points of this work here. However, the

general situation is as follows:

530 E.W. Larsen

1. A straightforward iteration method, called Source Iteration (SI),converges

in one iteration if the problem has no scattering or ﬁssion. For problems

with scattering and ﬁssion, SI converges rapidly if the processes of capture

and leakage are not weak. In this case, neutron histories are generally

short; it is unlikely that a neutron will scatter a very large number of times

before being captured or leaking out of V . However, if a problem has a

diﬀusive subregion, which is optically thick and is dominated by scattering,

then SI will generally converge very slowly – because the particles in this

subregion will have long histories.

2. Numerous acceleration methods have been developed to increase the rate

of convergence of the SI technique. Most notably, the diﬀusion synthetic

acceleration (DSA) method can, under ideal conditions, produce spectac-

ular speedups over SI. However, DSA has stability issues associated with

the use of advanced discretization schemes and unstructured grids [22].

3. The use of Krylov subspace methods, with DSA as a preconditioner, has

been a major topic of research in recent years. It now appears that for

most problems of practical interest, the proper combination of DSA +

Krylov will yield an eﬃcient solution technique [25].

We conclude this section with a brief summary of the advantages and

disadvantages of deterministic methods. First, by their nature, determinis-

tic methods automatically produce a global solution on the chosen grid in

phase space. However, these solutions will contain truncation errors associ-

ated with the discretizations of the energy, direction, and spatial variables.

While truncation errors generally decrease as the phase space grid is reﬁned,

some errors can be notoriously diﬃcult to suppress, and because of the high

dimensionality of phase space, it may be diﬃcult (due to computer memory

limitations) to employ a suﬃciently ﬁne grid. Also, deterministic methods are

generally formulated on spatial grids that have a limited capacity to model

curved surfaces. Thus, the geometric modeling of the physical system may

itself contain “truncation” errors. Finally, iterative techniques for solving the

ﬁnal discrete set of equations are not always as eﬃcient as desired. As one

might expect, all of these deﬁciencies are topics of current research.

Detailed information on deterministic neutron transport methods can be

found in [2,3,10,11,13,14,16,22], and [25]. The dates of publication of these

references indicate how up-to-date they are; this is relevant because the ﬁeld

of deterministic transport simulations has changed immensely in the past 20

years. Two recent detailed reviews are [22] and [25].

6 Automatic Variance Reduction (Hybrid) Methods

We have already discussed a serious deﬁciency of classic Monte Carlo codes:

often, they are not user-friendly. (In addition to specifying the geometry and

Neutron Transport Problems and Simulation Techniques 531

material structure of the system, the Monte Carlo code user must also specify

biasing parameters. No analogous task is required of an S

code user.)

In the past 10 years, it has been recognized that the biasing parameters –

which (to repeat) have traditionally been chosen by the code user – can in fact

be calculated by a relatively inexpensive deterministic adjoint calculation. In

principle, and in fact, it is possible to: (i) ﬁrst, run a suitable (inexpensive)

deterministic adjoint problem, (ii) use the results of that calculation to deter-

mine biasing parameters, (iii) implement these biasing parameters in a Monte

Carlo code, and (iv) run the Monte Carlo code with the machine-generated

biasing parameters. This general process has been described as an Automatic

Variance Reduction method, or a Hybrid Monte Carlo-Deterministic method.

The commercial hybrid neutron transport code MCBEND empoys de-

terministic diﬀusion with Monte Carlo transport [17, 20]. A similar eﬀort

combining deterministic S

and Monte Carlo codes is described in [19]. This

and other related work is discussed in a recent review [23], and case studies

have been reported [24], showing the eﬃcacy of the approach. The computer-

generated biasing parameters are usually obtained by the computer much

more quickly than by human (trial and error) eﬀort, and are usually more ef-

ﬁcient than the ones obtained by the code user (the subsequent Monte Carlo

simulation runs much more eﬃciently). The principle drawback to the user is

that presently, two codes must be used: a deterministic code for the prelim-

inary adjoint calculation, and a subsequent Monte Carlo code. The problem

must be set up on both codes, and the results of the deterministic code must

be shipped to the Monte Carlo code. Thus, in current implementations, hy-

brid methods are not as convenient to use as they could be, if they were

implemented in one relatively easy-to-use code.

Due to their relative unfamiliarity and the awkwardness of using two codes

to implement them, the use of hybrid methods is currently very limited. How-

ever, hybrid methods oﬀer computational advantages that, sooner or later,

will motivate code developers to implement them in user-friendly ways. When

this happens, hybrid methods should become much more widely-adopted.

7 Discussion

The simulation of neutron and photon transport processes for nuclear re-

actors is generally diﬃcult and costly. During the past 50-odd years, two

distinct computational techniques have been developed to implement these

simulations. These techniques (Monte Carlo and deterministic) are comple-

mentary, in ways that have been discussed above, and that have led to their

wide use in diﬀerent types of problems. Current research in Monte Carlo and

deterministic methods is aimed at making these methods run more eﬃciently

and accurately. Perhaps the best sources of literature on these topics are

(i) journal articles published in Nuclear Science and Engineering (the pre-

mier research journal of the American Nuclear Society), (ii) proceedings of