Nonclassical Particle Transport in the 1-D Diffusive Limit

0.78 MB·

by  Richard Vasques

Nonclassical Particle Transport in the 1-D Diffusive Limit
Richard Vasquesa,∗, Rachel Slaybaugha, Kai Kryckib
aDept. of Nuclear Engineering, University of California, Berkeley
bAachen Institute for Nuclear Training GmbH, Aachen, Germany
1. Introduction
6
1
A nonclassical linear Boltzmann equation has been recently proposed to address
0
2 nonexponential attenuation of the particle flux in certain inhomogeneous random

n media applications [1]. In particular, this effect arises in Pebble Bed reactor cores, in
a
which the locations of the pebbles are spatially correlated [2].
J
In this paper, we investigate nonclassical particle transport taking place in a 1-D
5
random periodic diffusive system. We provide computational results that validate

]
h the theoretical predictions, demonstrating for the first time that the solution of the
-t nonclassical particle transport equation is well-approximated by the solution of the
l
c nonclassical diffusion equation.
u
For simplicity, we consider the following assumptions: (i) transport occurs in rod
n
[ geometry, in which particles can only move in the directions µ = ±1; (ii) transport

1 is monoenergetic; and (iii) scattering is isotropic. In this case, the 1-D nonclassical
v linear Boltzmann equation is written as
5
9 ∂ψ± ∂ψ±
4 (x,s)± (x,s)+Σ (s)ψ±(x,s) (1.1)
t
2 ∂s ∂x
0 δ(s) (cid:20) (cid:90) ∞ (cid:21)
(cid:2) (cid:3)
. = c Σ (s(cid:48)) ψ±(x,s(cid:48))+ψ∓(x,s(cid:48)) ds(cid:48) +Q(x) ,
1 2 t
0 0
6 where x = position, s = the path-length traveled by the particle since its previous
1
interaction (birth or scattering), ψ± is the nonclassical angular flux in the directions
:
v
±1, c is the scattering ratio (probability of scattering), and Q(x) is an isotropic
i
X
internal source. The function Σ (s) represents the collision probability (ensemble-
t
r
averaged over all possible physical realizations of the system), such that
a
the probability that a particle, scattered or
Σ (s)ds = born at any point x, will experience a
t
collision between x+s and x+(s+ds).
In this situation, the probability density function for a particle’s distance-to-collision
is given by
p(s) = Σt(s)e−(cid:82)0sΣt(s(cid:48))ds(cid:48), (1.2)
∗Corresponding author: [email protected]
Dept. of Nuclear Engineering, University of California, Berkeley, 4103 Etcheverry Hall, MC 1730,
Berkeley, CA 94720-1730
Preprint submitted to ArXiv.org January 12, 2016
such that its mth moment is defined as
(cid:90) ∞
(cid:10) (cid:11)
sm = smp(s)ds.
0
If Σ (s) = Σ = constant (classical total cross section), we obtain the exponential
t t
p(s) = Σ e−Σts, (1.3)
t
and Eq. (1.1) reduces to the classical linear Boltzmann equation
∂Ψ± Σ Q(x)
± (x)+Σ Ψ±(x) = s (cid:2)Ψ±(x)+Ψ∓(x)(cid:3)+
t
∂x 2 2
for the classical angular flux
(cid:90) ∞
Ψ±(x) = ψ±(x,s)ds.
0
2. Asymptotic Analysis
Following [1], we scale Σ = O(1), 1−c = O(ε2), Q = O(ε2), ∂ψ/∂s = O(1), and
t
∂ψ/∂x = O(ε), with ε (cid:28) 1. In this scaling, Eq. (1.1) yields
∂ψ± ∂ψ±
(x,s)±ε (x,s)+Σ (s)ψ±(x,s) (2.1)
t
∂s ∂x
1−ε2(1−c) (cid:90) ∞
(cid:2) (cid:3)
= δ(s) Σ (s(cid:48)) ψ±(x,s(cid:48))+ψ∓(x,s(cid:48)) ds(cid:48)
t
2
0
Q(x)
+ε2δ(s) .
2
Let us define ψˆ±(x,s) such that
ψ±(x,s) ≡ ψˆ±(x,s)e−(cid:82)0sΣt(s(cid:48))ds(cid:48) .
(cid:10) (cid:11)
s
Then, using Eq. (1.2), Eq. (2.1) becomes the following equation for ψˆ±(x,s):
∂ψˆ± ∂ψˆ±
(x,s)±ε (x,s)
∂s ∂x
1−ε2(1−c) (cid:90) ∞(cid:104) (cid:105)
= δ(s) ψˆ±(x,s(cid:48))+ψˆ∓(x,s(cid:48)) p(s(cid:48))ds(cid:48)
2
0
Q(x)
(cid:10) (cid:11)
+ε2δ(s) s .
2
This equation is mathematically equivalent to:
∂ψˆ± ∂ψˆ±
(x,s)±ε (x,s) = 0 s > 0, (2.2a)
∂s ∂x
2
and
ψˆ±(x,0) = (2.2b)
1−ε2(1−c) (cid:90) ∞ (cid:104) (cid:105)
= p(s(cid:48)) ψˆ±(x,s(cid:48))+ψˆ∓(x,s(cid:48)) ds(cid:48)
2
0
Q(x)
(cid:10) (cid:11)
+ε2 s ,
2
where ψˆ±(x,0) = ψˆ±(x,0+). Integrating Eq. (2.2a) over 0 < s(cid:48) < s, we obtain:
∂ (cid:90) s
ψˆ±(x,s) = ψˆ±(x,0)±ε ψˆ±(x,s(cid:48))ds(cid:48)
∂x
0
1−ε2(1−c) (cid:90) ∞ (cid:104) (cid:105)
= p(s(cid:48)) ψˆ±(x,s(cid:48))+ψˆ∓(x,s(cid:48)) ds(cid:48)
2
0
Q(x) ∂ (cid:90) s
+ε2(cid:10)s(cid:11) ∓ε ψˆ±(x,s(cid:48))ds(cid:48).
2 ∂x
0
Introducing into this equation the ansatz
∞
(cid:88)
ψˆ±(x,s) = εnψˆ±(x,s)
n
n=0
and equating the coefficients of different powers of ε, we obtain for n ≥ 0:
1 (cid:90) ∞ (cid:104) (cid:105)
ψˆ±(x,s) = p(s(cid:48)) ψˆ±(x,s(cid:48))+ψˆ∓(x,s(cid:48)) ds(cid:48) (2.3)
n 2 n n
0
∂ (cid:90) s
∓ ψˆ± (x,s(cid:48))ds(cid:48)
∂x n−1
0
1−c (cid:90) ∞ (cid:104) (cid:105)
− p(s(cid:48)) ψˆ± (x,s(cid:48))+ψˆ∓ (x,s(cid:48)) ds(cid:48)
2 n−2 n−2
0
Q(x)
(cid:10) (cid:11)
+δ s ,
n,2
2
with ψˆ± = ψˆ± = 0.
−1 −2
Equation (2.3) with n = 0 has the general solution
ˆ
φ (x)
ψˆ±(x,s) = 0 ,
0 2
ˆ
where φ (x) is undetermined at this point. For n = 1, Eq. (2.3) has a particular
0
solution of the form:
ˆ
sdφ
ψˆ± (x,s) = ∓ 0(x),
part 2 dx
and its general solution is given by
(cid:34) (cid:35)
ˆ
1 dφ
ψˆ±(x,s) = φˆ±(x)∓s 0(x) ,
1 2 1 dx
3
ˆ
where φ (x) is undetermined.
1
Equation (2.3) with n = 2 has a solvability condition, which is obtained by adding
the equations for ψ+ and ψ− and operating on them by (cid:82)∞p(s)(·)ds; the solvability
2 2 0
condition yields
(cid:10)s2(cid:11)d2φˆ
0 ˆ (cid:10) (cid:11)
0 = (x)−(1−c)φ (x)+ s Q(x).
2 dx2 0
We can rewrite this equations as
(cid:10)s2(cid:11)d2φˆ 1−c
0 ˆ
− (x)+ φ (x) = Q(x), (2.4)
2(cid:10)s(cid:11) dx2 (cid:10)s(cid:11) 0
which is the nonclassical diffusion equation for Eq. (1.1).
Therefore, the solution ψ±(x,s) of Eq. (2.1) satisfies
ψ±(x,s) = φˆ0(x)e−(cid:82)0sΣt(s(cid:48))ds(cid:48) +O(ε), (2.5)
(cid:10) (cid:11)
2 s
ˆ
where φ (x) satisfies Eq. (2.4). The classical angular flux can be obtained to leading
0
order by integrating Eq. (2.5) over 0 < s < ∞:
(cid:90) ∞ φˆ (x)
Ψ±(x) = ψ±(x,s)ds = 0 +O(ε).
2
0
(cid:10) (cid:11) (cid:10) (cid:11)
As expected, if p(s) is given by Eq. (1.3), s = 1/Σ , s2 = 2/Σ2, and it is easy
t t
to verify that Eq. (2.4) reduces to the classical diffusion equation
1 d2φˆ
0 ˆ
− (x)+Σ φ (x) = Q(x).
Σ dx2 a 0
t
3. The 1-D Random Periodic System
Let us consider an infinite rod consisting of periodically arranged solid and void
layers of equal width, given by (cid:96) = 1. We are interested in a finite random periodic
system with total width given by 2X = 2(cid:96)M, where the integer M (the total length
of each material in the system) satisfies M = ε−1. Random realizations of this system
can be obtained by randomly selecting any continuous 2X segment of the infinite rod
described above. Vacuum boundary conditions are assigned at x = ±X.
In this finite random periodic system, the cross sections and source are stochastic
functions of space. For the numerical results provided in this paper, we define the
parameters at each spatial point x in the system by
(cid:26)
1, if x is in solid
Σ (x) = ,
t 0, if x is in void
(cid:26)
2M−2, if x is in solid
Q(x) = ,
0, if x is in void
4
and the absorption ratio by
1−c = 1M−2.
These parameters are in agreement with the assumptions of our asymptotic analysis:
Σ is O(1), (1−c) and Q are O(ε2), and the system is optically thick with 2X being
t
O(1/ε). As M increases, ε decreases, and the 1-D system approaches the diffusive
limit.
4. Numerical Results
An analytical expression for p(s) can be obtained for this type of 1-D system
[3]. Taking into account the parameters considered in this paper, the path length
distribution function is given by
(cid:26)
(2n+1−s)e−(s−n), if 2n ≤ s ≤ 2n+1
p(s) =
(s−2n−1)e−(s−n−1), if 2n+1 ≤ s ≤ 2(n+1)
for n = 0,1,2,... (see Figure 1); this yields the moments
(cid:10) (cid:11)
s = 2,
2e
(cid:10) (cid:11)
s2 = 5+ ≈ 8.1640.
e−1
The ensemble-averaged collision probability can be written as [1]
p(s)ds
Σ (s)ds = .
t (cid:82)s
1− p(s(cid:48))ds(cid:48)
0
Finally, since the volume fraction of solid and void materials are the same in any
realization of the system, the ensemble-averaged source is simply given by
Q(x)
(cid:10) (cid:11)
Q = = 1M−2.
2
Following the procedure described in [3], we rewrite Eq. (1.1) in its initial value
form as
∂ψ± ∂ψ±
(x,s)± (x,s)+Σ (s)ψ±(x,s) = 0, (4.1a)
t
∂s ∂x
c (cid:90) ∞
ψ±(x,0) = Σ (s(cid:48))[ψ+(x,s(cid:48))+ψ−(x,s(cid:48))]ds(cid:48) (4.1b)
t
2
0
(cid:10) (cid:11)
Q
+ .
2
We then adapt the HLL finite volume method introduced in [4]; this method is of
first order in the pseudo-time variable s and in the spatial variable x. In our calcula-
tions, we have cut off the integration at s = 45, and have chosen the trapezoidal
max
5
rule. We have used the mesh interval (cid:52)x = 1/400 (= smallest ε2 considered in this
paper), and a CFL number 0.5 (that is, (cid:52)s = 1/800).
This system is solved in a source-iteration manner, where we iterate between the
discretized forms of Eq. (4.1b) and Eq. (4.1a). The ensemble-averaged nonclassical
scalar flux is given by Φ(x) = (cid:82)45[ψ+(x,s)+ψ−(x,s)]ds.
0
We consider three different problems, with M = 10,15, and 20. As M increases,
we expect the solution of Eq. (1.1) to converge to the solution of the nonclassical
diffusion formulation given by Eq. (2.4). It was shown in [4] that the contraction
rate for the source iteration is given by the scattering ratio c. To converge the
solution, 299 source iterations were needed for c = 0.99 (M = 10); 579 iterations for
c = 224/225 ≈ 0.9956 (M = 15); and 936 iterations for c = 0.9975 (M = 20).
Since our asymptotic analysis does not include boundary conditions, we solve Eq.
(2.4) using the extrapolated endpoint boundary conditions, taking the extrapolation
distance to be the diffusion coefficient:
(cid:32) (cid:10) (cid:11)(cid:33) (cid:32) (cid:10) (cid:11)(cid:33)
s2 s2
ˆ ˆ
φ X + = φ −X − = 0.
0 (cid:10) (cid:11) 0 (cid:10) (cid:11)
2 s 2 s
As anticipated, the solutions plotted in Figures 2, 3, and 4 confirm our claim: the
nonclassical scalar flux and the nonclassical diffusion scalar flux increasingly agree as
M increases and the system approaches the diffusive limit.
5. Conclusions
This paper provides numerical results that demonstrate the validity of the non-
classical diffusion approximation to the nonclassical transport equation in certain 1-D
diffusive systems. To our knowledge, this is the first time computational results that
validate the predictions of the asymptotic theory for the nonclassical model are pre-
sented. This result provides a more solid foundation in which to improve this theory
for relevant nuclear applications.
Acknowledgments
This paper was prepared by Richard Vasques and Rachel Slaybaugh under award
numberNRC-HQ-84-14-G-0052fromtheNuclearRegulatoryCommission. Thestate-
ments, findings, conclusions, and recommendations are those of the authors and do
not necessarily reflect the view of the U.S. Nuclear Regulatory Commission.
References
[1] Vasques, R., Larsen, E.W. Non-classical particle transport with angular-dependent
path-length distributions. I: Theory. Annals of Nuclear Energy 2014; 70: 292-300.
[2] Vasques, R., Larsen, E.W. Non-classical particle transport with angular-dependent
path-length distributions. II: Application to pebble bed reactor cores. Annals of Nu-
clear Energy 2014; 70: 301-311.
6
[3] Vasques, R., Krycki, K. On the accuracy of the non-classical transport equation in 1-D
random periodic media. In: Proc. Joint International Conference on Mathematics and
Computation, Supercomputing in Nuclear Applications and the Monte Carlo Method;
Nashville, TN, April 19-23, 2015.
[4] Krycki, K., Berthon, C., Frank, M., Turpault, R. Asymptotic preserving numerical
schemes for a non-classical radiation transport model for atmospheric clouds. Mathe-
matical Methods in the Applied Sciences 2013; 36: 2101-2116.
7
Figure 1: Probability density function for distance-to-collision.
0.35
Transport
Diffusion
0.3
0.25
x
u
Fl
r 0.2
a
al
c
S
0.15
0.1
0.05
-10 -5 0 5 10
x
Figure 2: Scalar fluxes obtained by solving the nonclassical transport and the non-
classical difusion equations with M = 10.
8
0.35
Transport
Diffusion
0.3
0.25
x
u 0.2
Fl

r
a
al 0.15
c
S
0.1
0.05
0
-15 -10 -5 0 5 10 15
x
Figure 3: Scalar fluxes obtained by solving the nonclassical transport and the non-
classical difusion equations with M = 15.
0.35
Transport
Diffusion
0.3
0.25
x
u 0.2
Fl

r
a
al 0.15
c
S
0.1
0.05
0
-20 -15 -10 -5 0 5 10 15 20
x
Figure 4: Scalar fluxes obtained by solving the nonclassical transport and the non-
classical difusion equations with M = 20.
9