This section describes the coupled MG and PW techniques used to solve
the neutron transport equation.

Energy/lethargy ranges for MG and PW calculations

The combined MG/PW CENTRM calculation is performed over the energy range
spanned by the group structure in the input MG library. The energy
boundaries for the “IGM” neutron groups specified on the MG library
divide the entire energy range into energy intervals. The lowest energy
group contained in the UMR is defined to be “MGHI”; while the highest
energy group in the LMR is designated “MGLO.” The boundary between the
PW and UMR energy intervals is set by the energy value “DEMAX,” while
“DEMIN” is the boundary between the PW and LMR. The default values of
0.001 eV and 20 keV for DEMIN and DEMAX, respectively, can be modified
by user input, but the input values are altered by the code to
correspond to the closest group boundaries. Hence, DEMAX is always equal
to the lower energy boundary of group MGHI and DEMIN the upper energy
boundary of MGLO. The PW calculation is performed in terms of lethargy
(u), rather than energy (E). The origin (u=0) of the lethargy coordinate
corresponds to the energy E=DEMAX, which is the top of the PW range. See
Fig. 207.

The highest energy group of the thermal range is defined by the
parameter “IFTG,” obtained from the MG library. If DEMIN is less than
the upper energy boundary of IFTG, the PW range extends into thermal. In
this case, scattering in the PW region of the thermal range is based on
the PW scattering kernel data; and the LMR calculation uses 2D transfer
matrices for incoherent and coherent scattering on the MG library.
Coupling between the MG and PW thermal calculations is treated, and
outer iterations are required to address effects of upscattering.


With the exception of hydrogen moderation, elastic downscattering
coupling the UMR and PW ranges, occurs only within a limited subrange
of the UMR called the “transition region”. The highest energy group in
the transition region is designated “MGTOP.” The precise definition of
the transition region is given in Scattering sources for the PW range.

Energy boundaries of the group structure on the input MG library
correspond to the IGM+1 values, { G_{1}, G_{2,} …
G_{g,} G_{g+1}, …, G_{IGM+1}}. It is convenient to
designate the number of groups in the UMR, PW, and LMR ranges equal to
NG_{U}, NG_{P}, and NG_{L}, respectively, so that IGM
= NG_{U} + NG_{P} + NG_{L}; or in terms of the
parameters MGHI and MGLO introduced previously:

NG_{U} = MGHI; NG_{P} = MGLO − MGHI − 1; NG_{L} = IGM
− MGLO + 1.


The flux per unit lethargy is calculated for a discrete energy (or
lethargy) mesh spanning the MG structure. Groups in the UMR and LMR each
contain a single energy mesh point, while groups in the PW range
generally contain several points. The number of mesh points in the
UMR, PW, and LMR is equal respectively to NG_{U}, N_{P}, and
NG_{L}; and the total number of points in the entire energy mesh
is designated as “N_{T},” which is equal to NG_{U} +
N_{P} + NG_{L}. Thus the lethargy (u) mesh consists of the
set of points: {u:sub:1,….u_{NGU,}
u_{NGU+1,}….u_{NGU+NP},
u_{NGU+NP+1},…u_{NT}}. Based on the lethargy origin at
E=DEMAX, the lethargy “u_{n}” associated with any energy point
“E_{n}” is equal to,

u_{n} = ln(DEMAX/E_{n}).


Lethargy points are arranged in order of increasing value. The lethargy
origin is at point NG_{U}+1, the lower energy boundary of group
MGHI; i.e., u_{NGU+1}=0. Note that the entire UMR
(E>DEMAX) corresponds to negative lethargy values. Lethargy values for
the first NG_{U} and the last NG_{L} points in the mesh are
defined to be the midpoint lethargies of groups in the UMR and LMR
ranges, respectively. For example, for the NG_{U} groups within
the UMR,

u_{1} = 0.5[ln(DEMAX/G_{1}) + ln(DEMAX/G_{2})];

u_{NGU} = 0.5[ln(DEMAX/G_{MGHI}) + ln(DEMAX/G_{MGHI + 1})];


and similarly for the NG_{L} groups in the LMR,

u_{NGU + NP + 1} = 0.5[ln(DEMAX/G_{MGLO}) + ln(DEMAX/G_{MGLO + 1})]

u_{NT} = 0.5[ln(DEMAX/G_{IGM}) + ln(DEMAX/G_{IGM + 1})]


The remaining N_{P} points in the mesh (i.e., values
u_{NGU+1} to u_{NGU+NP}) are contained within the
NG_{P} groups that span the PW range. By definition the first
point in the PW range is the lower energy boundary of group MGHI. The
other mesh points are computed internally by CENTRM, based on the
behavior of the macroscopic PW total cross sections and other criteria.

The neutron flux, as a function of space and direction, is calculated
for each energy/lethargy point in the mesh by solving the Boltzmann
transport equation. The transport equation at each lethargy point
generally includes a source term representing the production rate due to
elastic and inelastic scatter from other lethargies, which couples the
solutions at different lethargy mesh points. Except in the thermal
range, neutrons can only gain lethargy (lose energy) in a scattering
reaction; thus the PW flux is computed by solving the transport equation
at successive mesh points, sweeping from low to high lethargy values.


The Boltzmann equation for neutron transport

The steady state neutron transport equation shown below represents a
particle balanceper unit phase space, at an arbitrary point ρ in phase
space,

(326)\[\Omega \cdot \nabla \Psi(\rho)+\sum_{t}(\mathrm{r}, \mathrm{u}) \Psi(\rho)=\int_{0}^{\infty} \int_{0}^{4 \pi} \Sigma\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right) \Psi\left(\mathrm{u}^{\prime}, \Omega^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}+\mathrm{Q}_{\mathrm{ext}}(\rho)\]

where:

ψ(p) = angular flux (per lethargy) at phase space coordinate ρ;

ρ = (r,u,Ω) = phase space point defined by the six independent
variables;

r = (x_{1},x_{2},x_{3}) = space coordinates;

u = ln(E_{ref}/E) = lethargy at energy E, relative to an origin
(u=0) at E_{ref};

Ω = (μ,ζ) = neutron direction defined by polar cosine μ and azimuthal
angle ζ;

Σ_{t}(r,u) = macroscopic total cross section;

Σ(u′→u;μ_{0}) = double differential scatter cross section;

μ_{0} = cosine of scatter angle, measured in laboratory coordinate
system;

Q_{ext}(ρ) = external source term, including fission source;


The left and right sides of (326) respectively, are equal to the neutron
loss and production rates, per unit volumedirectionlethargy. In CENTRM
the spatial distribution of the fission source is input as a component
of the external source Q; hence, a fixed source rather than an
eigenvalue calculation is required for the transport solution.

The angular dependence of the doubledifferential macroscopic scatter
cross section of an arbitrary nuclide “j” is represented by a finite
Legendre expansion of arbitrary order L:

(327)\[\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right)=\sum_{=0}^{\mathrm{L}} \frac{2+1}{2} \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \quad \mathrm{P}\left(\mu_{0}\right)\]

where P_{ℓ}(μ:sub:0) = Legendre polynomial evaluated at the
laboratory scattering cosine μ_{0}; and

\(\Sigma^{(j)}\left(u^{\prime} \rightarrow u\right)\) = cross section moments of nuclide j, defined by the expression


(328)\[\Sigma^{(j)}\left(u^{\prime} \rightarrow u\right)=\int_{1}^{1} \Sigma^{(j)}\left(u^{\prime} \rightarrow u ; \mu_{0}\right) P\left(\mu_{0}\right) d \mu_{0}\]

After substitution of the above Legendre expansions for the scattering
data of each nuclide, and applying the spherical harmonic addition
theorem in the usual manner, the scattering source on the right side of
(326) becomes [BG70]:

(329)\[\mathrm{S}(\mathrm{r}, \mathrm{u}, \Omega) \equiv \int_{0}^{\infty} \int_{0}^{4 \pi} \Sigma\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right) \Psi\left(\mathrm{u}^{\prime}, \Omega^{\prime}\right) \mathrm{d} \Omega^{\prime} \mathrm{d} \mathrm{u}^{\prime}=\sum_{\mathrm{k}=1}^{\mathrm{LK}} \frac{2+1}{2} \mathrm{Y}_{\mathrm{k}}(\Omega) \mathrm{S}_{\mathrm{k}}(\mathrm{r}, \mathrm{u})\]

wherein,

\(\mathrm{Y}_{\mathrm{k}}(\Omega)=\mathrm{Y}_{\mathrm{k}}(\mu, \zeta)\) = the spherical harmonic function evaluated at direction Ω

S_{k} = spherical harmonic moments of the scatter source, per unit letharagy.


The summation index “ℓk” indicates a double sum over ℓ and k indices; in
the most general case it is defined as:

\[\sum_{\mathrm{k}=1}^{\mathrm{LK}}=\sum_{=0}^{\mathrm{L}} \sum_{\mathrm{k}=0}\]

where “L” is the input value for the maximum order of scatter (input
parameter “ISCT”).

Due to symmetry conditions, some of the source moments may be zero. The
parameter LK is defined to be the total number of nonzero moments
(including scalar flux) for the particular geometry of interest, and is
equal to,

LK = L + 1 for 1D slabs and spheres;

LK = L*(L+4)/4+1 for 1D cylinders, and

LK = L*(L+3)/2+1 for 2D MoC cells


More details concerning the 1D Boltzmann equation can be found in the
XSDRNPM chapter of the SCALE manual.

The S_{ℓk} moments in (329) correspond to expansion coefficients in
a spherical harmonic expansion of the scatter source. These can be
expressed in terms of the cross section and flux moments by

(330)\[\mathrm{S}_{\mathrm{k}}(\mathrm{u})=\sum_{\mathrm{j}} \int_{\mathrm{u}^{\prime}} \mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{d} \mathrm{u}^{\prime}=\sum_{\mathrm{j}} \int_{\mathrm{u}^{\prime}} \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}\]

where ψ_{ℓk}(u) = spherical harmonic moments of the angular
flux;

(331)\[= \int_{0}^{4 \pi} \mathrm{Y}_{\mathrm{k}}(\Omega) \Psi(\Omega) \mathrm{d} \Omega\]

and S_{ℓk}^{(j)}(u′→u) = moments of the differential
scatter rate from lethargy u′ to u, for nuclide “j”;

(332)\[= \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right)\]

The ψ_{ℓk} flux moments are the well known coefficients appearing
in a spherical harmonic expansion of the angular flux. These usually are
the desired output from the transport calculation. In particular, the
ℓ=0, k=0 moment corresponds to the scalar flux [indicated here as
Φ(r,u)],

(333)\[\Psi_{0,0}(\mathrm{r}, \mathrm{u})=\Phi(\mathrm{r}, \mathrm{u})=\int_{0}^{4 \pi} \Psi(\mathrm{r}, \mathrm{u}, \Omega) d \Omega\]

In general the epithermal component of the scatter source in (329)
contains contributions from both elastic and inelastic scatter
reactions; however, inelastic scatter is only possible above the
threshold energy corresponding to the lowest inelastic level. The
inelastic Q values for most materials are typically above 40 keV;
therefore, elastic scatter is most important for slowing down
calculations in the resolved resonance range of most absorber materials
of interest. For example, the inelastic Q values of ^{238}U, iron,
and oxygen are approximately 45 keV, 846 keV, and 6 MeV, respectively;
while the upper energy of the ^{238}U resolved resonance range is
20 keV in ENDF/BVII. The inelastic thresholds of some fissile materials
like ^{235}U and ^{239}Pu are on the order of 10 keV;
however, with the exception of highly enriched fast systems, these
inelastic reactions usually contribute a negligible amount to the
overall scattering source. CENTRM assumes that continuum inelastic
scatter is isotropic in the laboratory system, while discrete level
inelastic scatter is isotropic in the center of mass (CM) coordinate
system.

Over a broad energy range, elastic scatter from most moderators can
usually be assumed isotropic
(swave) in the neutronnucleus CM coordinate system. In the case of
hydrogen, this is true up to approximately 13 MeV; for carbon up to
2 MeV; and for oxygen up to 100 keV. However, it is well known that
isotropic CM scatter does not result in isotropic scattering in the
laboratory system. For s‑wave elastic scatter the average
scattercosine in the laboratory system is given by: \(\bar{\mu}_{0}=0.667 / \mathrm{~A};^{3}\) where
“A” is the mass number (in neutron mass units) of the scattering
material. This relation indicates that swave, elastic scattering
from low A materials tends to be more anisotropic in the laboratory,
and that the laboratory scattering distribution approaches isotropic
\(\left(\bar{\mu}_{0}=0 ; \theta_{0}=90\right)\) as A becomes large. For example, the \(\bar{\mu}_{0}\)
of hydrogen
is 0.667 (48.2°); while it is about 0.042 (87.6°) for oxygen. Because
s‑wave scattering from heavy materials is nearly isotropic in the
laboratory system, the differential scattering cross section (and thus
the scattering source) can usually be expressed accurately by a low
order Legendre expansion. On the other hand light moderators like
hydrogen may require more terms—depending on the flux anisotropy—to
accurately represent the elastic scatter source in the laboratory
system. The default settings in CENTRM are to use P0 (isotropic lab
scatter) for mass numbers greater than A=100, and P1 for lighter
masses.

An analytical expression can be derived for the crosssection moments in
the case of twobody reactions, such as elastic and discretelevel
inelastic scattering from “stationary” nuclei. Stationary here implies
that the effect of nuclear motion on neutron scattering kinematics is
neglected.


Note

The stationary nucleus approximation for treating
scattering kinematics does not imply that the effect of nuclear motion
on Doppler broadening of resonance cross sections is ignored, since this
effect is included in the PW crosssection data.


In CENTRM the stationary nucleus approximation is applied above the
thermal cutoff, typically around 35 eV, but is not valid for low energy
neutrons. CENTRM has the capability to perform a PW transport
calculation in the thermal energy range using tabulated thermal
scattering law data for bound molecules, combined with the analytical
freegas kernel for other materials. In this case the crosssection
moments appearing in (328) include upscattering effects. The expressions
used in CENTRM to compute the PW scatter source moments in the thermal
range are given in The Boltzmann equation within the PW range.

The following two sections discuss the evaluation of the scatter source
moments for epithermal elastic and inelastic reactions, respectively.


Epithermal Elastic Scatter

Consider a neutron with energy E′, traveling in a direction Ω′, that
scatters elastically from an arbitrary material “j,” having a
mass A^{(j)} in neutronmass units. Conservation of kinetic energy
and momentum requires that there be a unique relation between the angle
that the neutron scatters (relative to the initial direction) and its
final energy E after the collision. If the nucleus is assumed to be
stationary in the laboratory coordinate system, then the
cosine (μ_{0}) of the scatter angle (θ_{0}) measured in the
laboratory system, as a function of the initial and final energies, is
found to be

(334)\[\mu_{0} \equiv \Omega^{\prime} \cdot \Omega=\mathrm{G}^{(\mathrm{j})}\left(\mathrm{E}^{\prime}, \mathrm{E}\right) ,\]

where the kinematic correlation function G relating E′, E, and
μ_{0} for elastic scatter is equal to

(335)\[\begin{split} \begin{array}{l}
\mathrm{G}^{(\mathrm{j})}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\mathrm{a}_{1}^{(\mathrm{j})} \times\left[\mathrm{E} / \mathrm{E}^{\prime}\right]^{\frac{1}{2}}\mathrm{a}_{2}^{(\mathrm{j})} \times\left[\mathrm{E}^{\prime} / \mathrm{E}\right]^{\frac{1}{2}} \\
\text { and } \mathrm{a}_{1}^{(\mathrm{j})}=\left(\mathrm{A}^{(\mathrm{j})}+1\right) / 2 \quad ; \quad \mathrm{a}_{2}^{(\mathrm{j})}=\left(\mathrm{A}^{(\mathrm{j})}1\right) / 2
\end{array} .\end{split}\]

The final energy E of an elastically scattered neutron is restricted to
the range,

(336)\[\alpha^{(j)} E^{\prime} \leq E \leq E^{\prime}\]

where α^{(j)} = maximum fractional energy lost by elastic scatter

(337)\[= \left[\mathrm{a}_{2}^{(\mathrm{j})} / \mathrm{a}_{1}^{(\mathrm{j})}\right]^{2}\]

The corresponding range of scattered neutrons in terms of lethargy is equal to

(338)\[\mathrm{u}^{\prime} \leq \mathrm{u} \leq \mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})}\]

where

(339)\[\begin{split} \begin{aligned}
\mathrm{u}, \mathrm{u}^{\prime} &=\mathrm{u}(\mathrm{E}), \mathrm{u}^{\prime}\left(\mathrm{E}^{\prime}\right)=\text { lethargies corresponding to } \mathrm{E} \text { and } \mathrm{E}^{\prime}, \text { respectively; and } \\
\varepsilon^{(\mathrm{j})} &=\text { maximum increase in lethargy, per elastic scatter }=\ln \left[1 / \alpha^{(j)}\right]
\end{aligned}\end{split}\]

The doubledifferential scatter kernel of nuclide j (per unit lethargy
and solid angle) for swave elastic scatter of neutrons from
stationary nuclei, is equal to

(340)\[\begin{split} \begin{aligned}
\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u} ; \mu_{0}\right) &=\frac{\mathrm{E} \Sigma^{(\mathrm{i})}\left(\mathrm{u}^{\prime}\right)}{\mathrm{E}^{\prime}\left(1\alpha^{(\mathrm{j})}\right)} \delta\left[\mu_{0}\mathrm{G}^{(\mathrm{j})}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)\right], \text { for } \mathrm{u}^{\prime} \leq \mathrm{u} \leq \mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})} \\
&=0 \quad \mathrm{u}<\mathrm{u}^{\prime} \text { or } \mathrm{u}>\mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})}
\end{aligned}\end{split}\]

The presence of the Dirac delta function completely correlates the angle
of scatter and the values of the initial and final energies.
Substituting the double differential crosssection expression from (340)
into (328) gives the singledifferential Legendre moments of the
cross section, per final lethargy:

(341)\[\begin{split} \begin{aligned}
\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) &=\frac{\mathrm{E} \mathrm{P}\left[\mathrm{G}^{(\mathrm{j})}\right] \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right)}{\mathrm{E}^{\prime}\left(1\alpha^{(\mathrm{j})}\right)}, \text { for } \mathrm{u}^{\prime} \leq \mathrm{u} \leq \mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})} \\
&=0 \quad \mathrm{u}^{\prime} \text { or } \mathrm{u}>\mathrm{u}^{\prime}+\varepsilon^{(\mathrm{j})}
\end{aligned}\end{split}\]

where P_{ℓ} = Legendre polynomial evaluated at argument
G^{(j)} equal to the scatter cosine.

When the above expressions are used in (330) , the following is obtained
for the ℓk moment of the epithermal elastic scattering source at
lethargy u:

(342)\[\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}(\mathrm{u})=\sum_{\mathrm{j}} \int_{\left.\mathrm{u}\varepsilon^{(\mathrm{i}}\right)}^{\mathrm{u}} \mathrm{S}_{\mathrm{k}}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{du}^{\prime}=\sum_{\mathrm{j}} \int_{\mathrm{u}\varepsilon^{(j)}}^{\mathrm{u}} \frac{\mathrm{E} \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{P}\left[\mathrm{G}^{(\mathrm{j})}\right]}{\mathrm{E}^{\prime}\left(1\alpha^{(\mathrm{j})}\right)} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{du}^{\prime} .\]



Epithermal Inelastic Scatter

If the input value of DEMAX is set above the inelastic threshold of some
materials in the problem, then inelastic scattering can occur in the
PW range. The pointwise transport calculation may optionally include
discretelevel and continuum inelastic reactions in computing the
PW scatter source moments. The multigroup calculations always consider
inelastic reactions.

Discretelevel inelastic reactions are twobody interactions, so that
kinematic relations can be derived relating the initial and final
energies and the angle of scatter. It can be shown that the kinematic
correlation function for discretelevel inelastic scatter can be written
in a form identical to that for elastic scatter by redefining the
parameter a_{1} in (335) to be the energy dependent function [Wil00],

(343)\[\mathrm{a}_{1}^{(\mathrm{~m}, \mathrm{j})} = \frac{\left(\mathrm{A}^{(\mathrm{j})}+1\right)}{2}+\frac{\left(\mathrm{Q}^{(\mathrm{m}, \mathrm{j})}\right) \mathrm{A}^{(\mathrm{j})}}{2 \mathrm{E}}\]

The parameter Q^{(m, j)} is the Qvalue for the m_{th} level
of nuclide “j”. The Q value is negative for inelastic scattering, while
it is zero for elastic scatter. The threshold energy in the laboratory
coordinate system is proportional to the Qvalue of the inelastic level,
and is given by:

\[\mathrm{E}_{\mathrm{T}}^{(\mathrm{m}, \mathrm{j})}=\frac{\mathrm{A}^{(\mathrm{j})}+1}{\mathrm{~A}^{(\mathrm{j})}} \times\left(\mathrm{Q}^{(\mathrm{m}, \mathrm{j})}\right)\]

The range of energies that can contribute to the scatter source at E,
due to inelastic scatter from the m_{th} level of nuclide j is
defined to be [E:sub:L , E:sub:H ], where
E_{H} >E:sub:L >E:sub:T . This energy range has a
corresponding lethargy range of [u:sub:LO , u:sub:HI ] which is
equal to,

\[\begin{split}\begin{array}{l}
\mathrm{u}_{\mathrm{LO}}=\mathrm{u}\ln \left(\frac{1}{\alpha_{1}^{(\mathrm{j})}\left(\mathrm{E}_{\mathrm{H}}\right)}\right)=\mathrm{u}\varepsilon_{1}^{(\mathrm{j})} \\
\mathrm{u}_{\mathrm{HI}}=\mathrm{u}\ln \left(\frac{1}{\alpha_{2}^{(\mathrm{j})}\left(\mathrm{E}_{\mathrm{L}}\right)}\right)=\mathrm{u}\varepsilon_{2}^{(\mathrm{j})}
\end{array}\end{split}\]

The energydependent alpha parameters in the above expressions are defined as,

\[\begin{split}\begin{array}{l}
\alpha_{1}(\mathrm{E})=\left(\frac{\mathrm{A} \Delta^{(\mathrm{m}, \mathrm{j})}(\mathrm{E})1}{\mathrm{~A}+1}\right)^{2} \\
\alpha_{2}(\mathrm{E})=\left(\frac{\mathrm{A} \Delta^{(\mathrm{m}, \mathrm{j})}(\mathrm{E})+1}{\mathrm{~A}+1}\right)^{2}
\end{array}\end{split}\]

where

(344)\[\Delta^{(\mathrm{m}, \mathrm{j})}(\mathrm{E})=\sqrt{1\frac{\mathrm{E}_{\mathrm{T}}^{(\mathrm{m}, \mathrm{j})}}{\mathrm{E}}}\]

Modifying the epithermal elastic scatter source in (342) to include
discretelevel inelastic scatter gives the following expression

(345)\[\mathrm{S}_{\mathrm{k}}(\mathrm{u})=\sum_{m, j} \int_{u_{L O}}^{u_{H I}} \frac{\mathrm{E}}{\mathrm{E}^{\prime}} \frac{\Sigma^{(\mathrm{m}, \mathrm{j})}\left(\mathrm{E}^{\prime}\right) \mathrm{P}\left[\mathrm{G}^{(\mathrm{m}, \mathrm{j})}\right]}{\left(1\alpha^{j}\right) \Delta^{(\mathrm{m}, \mathrm{j})}\left(\mathrm{E}^{\prime}\right)} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{d} u^{\prime}\]

Detailed expressions for the lethargy limits are given in [Wil00]. Since
Δ^{(m,j)} is equal to unity for elastic scatter, the above
equation reduces to (339) if there is no discretelevel inelastic
contribution.

At high energies, the inelastic levels of the nucleus become a
continuum. In this case CENTRM represents the energy distribution of the
scattered neutrons by an evaporation spectrum with an isotropic angular
distribution in the lab system; thus, only the P_{0} moment
appears in the continuum inelastic scattering source. Including
continuum inelastic reactions in the PW calculation usually has a small
impact on the spectrum used for resonance selfshielding, and may
adversely impact the computer memory requirements and execution time.
Therefore, by default, CENTRM does not include continuum inelastic
reactions in the pointwise solution; however, it is always included in
the UMR solution.



Thermal Scatter

Since thermal neutrons have energies comparable to the mean kinetic
energy of molecules in thermal equilibrium, the scattering kernels must
account for molecular motion. The scatter moments include both
downscatter as well as upscatter contributions; hence, the integration
limits appearing in (342) must be extended from the lowest to the highest
energy of the thermal range. Furthermore the crosssection moments
correspond to the Legendre expansion coefficients of the thermal scatter
kernel, which has a substantially different form than the epithermal
kernel discussed in the previous two sections. In general the
ℓ_{th} Legendre moment of the thermal scattering kernel at
temperature T, describing scattering from E to E′, is given by

\[\sigma \quad\left(\mathrm{E}^{\prime} \rightarrow \mathrm{E} ; \mathrm{T}\right)=\frac{\sigma_{b}}{T} \sqrt{\frac{\mathrm{E}}{\mathrm{E}^{\prime}}} e^{\frac{\beta\left(\mathrm{E}^{\prime} \rightarrow \mathrm{E}\right)}{2}} \int \mathrm{P}\left(\mu_{0}\right) \mathrm{S}[\alpha, \beta ; T] \quad d \mu_{0}\]

where β(E′→E) and α(E′,E,μ_{0}) are dimensionless variables
(functions of temperature) defining the energy and momentum exchange,
respectively, of the collision [BG70]; σ_{b} is the rigidly
bound scatter cross section, which is proportional to the free atom
cross section; and S(α, β; T) describes the temperaturedependent
thermal scattering law.

If atomic bonding effects are neglected, the atoms of a material behave
like a gas in thermal equilibrium at the temperature of the medium. In
this case S(α, β) can be expressed by an analytical function. CENTRM
uses the free gas model for all nuclides except those materials that
have thermal scattering laws available in the ENDF/B nuclear data files.
The ENDF/B scattering law data account for the effects of molecular
bonding and possible polyatomic crystalline structure. While freegas
kernels are computed internally in CENTRM, the kernel moments describing
bound thermal scatterers are stored in a data file that can be accessed
by CENTRM.



Bound thermal kernels

Thermal scattering from bound atoms is classified either as an
“inelastic reaction,” in which the neutron energy may change, or an
“elastic reaction,” in which the neutron changes direction, but does not
change energy. In ENDF/B the former reactions are treated as incoherent
inelastic scattering with a doubly differential kernel describing the
secondary neutron energy and angle distribution. The latter reactions
are usual treated as coherent elastic scatter characterized by
diffractive interference of the scattered deBroglie waves, although a
few materials are modeled by the incoherent elastic approximation.
Legendre moments for thermal elastic kernels describe the secondary
angular distribution with no energy exchange, at a given neutron energy.
Bound scatter kernels have been processed by the AMPX code system for
most of the ~25 compounds with thermal scatter laws in ENDF/B, and are
stored in individual kinematics files distributed with the SCALE code
system. These include materials such as: H in water, H and C in
polyethylene, H and Zr in ZrH, C in graphite, deuterium in heavy water,
Be metal, Be in BeO, etc. The CRAWDAD module processes scattering kernel
data for individual nuclides into a combined library used in CENTRM, and
also interpolates the kernels to the appropariate temperatures.

The bound scatter kernels are tabulated at different energy points from
the flux solution mesh; therefore it is necessary to map the data onto
the desired energy mesh in the CENTRM calculation. Because thermal
elastic scattering results in no energy loss, the elastic moments only
appear in the withinpoint term of the scattering source in the CENTRM
thermal calculation. Thus the coherent elastic data is easily
interpolated since it only involves a single energy index and
temperature. However, the incoherent inelastic moments are 2D arrays in
terms of the initial and final energies, so that a 2D interpolation
must be done for each temperature. CENTRM uses a simple type of
“unitbase transform” method to interpolate incoherent inelastic kernels
onto the flux solution mesh. The method attempts to preserve the
absolute peak of the secondary energy distribution, at given initial
energy. For waterbound hydrogen and several other moderators, this is
quite adequate, since the kernel generally has only a single maximum.
However, if more than one local extrema is present, such as for
graphite, the other local peaks are not explicitly preserved in the
interpolation method. For this reason it is necessary to include a
fairly dense set of initial energies in the tabulated kernels of
graphite and similar materials, to avoid gross changes in the kernel
shape at adjacent initial energy panels.


Free gas thermal kernels

CENTRM computes freegas kernels using the approach proposed by
Robinson [Rob81] as a modification to the original
FLANGE [HF71] methodology.
Legendre moments of the freegas scatter kernel per unit lethargy are
expressed as,

(346)\[\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right)=\mathrm{A}^{(\mathrm{j})} \Sigma_{\mathrm{free}}^{(\mathrm{j})} \frac{\mathrm{E}}{\mathrm{E}^{\prime}} \mathrm{e}^{\beta / 2} \sum_{\mathrm{n}=0} \mathrm{W}_{\mathrm{n}} \mathrm{H}_{\mathrm{n}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)\]

where W_{ℓn} are constant coefficients associated with the
Legendre polynomial of order ℓ; Σ_{free} is the constant freeatom
cross section for the material; and H_{n} are the αmoments of the
freegas scatter law, given as

\[\mathrm{H}_{\mathrm{n}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\frac{1}{\sqrt{\pi}} \int_{\alpha_{\mathrm{L}}}^{\alpha_{\mathrm{H}}} \alpha^{\mathrm{n}} \times\left(\frac{\mathrm{e}^{\frac{\alpha^{2}+\beta^{2}}{4 \alpha}}}{2 \sqrt{\alpha}}\right) \mathrm{d} \alpha\]

The limits on the above integral correspond to:

\[\alpha_{\mathrm{L}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\alpha\left(\mathrm{E}^{\prime}, \mathrm{E}, \mu_{0}=1\right) \quad \text { and } \quad \alpha_{\mathrm{H}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=\alpha\left(\mathrm{E}^{\prime}, \mathrm{E}, \mu_{0}=1\right) .\]

The alpha moments for n > 0 can be evaluated very efficiently using a
recursive relation [Wil00]:

\[\mathrm{H}_{\mathrm{n}}\left(\mathrm{E}^{\prime}, \mathrm{E}\right)=2(2 \mathrm{n}1) \mathrm{H}_{\mathrm{n}1}+\beta^{2} \mathrm{H}_{\mathrm{n}2}\left[\mathrm{F}_{\mathrm{n}}\left(\sqrt{\alpha_{\mathrm{H}}}, \beta\right)\mathrm{F}_{\mathrm{n}}\left(\sqrt{\alpha_{\mathrm{L}}}, \beta\right)\right]\]

where F_{n} is the function,

\[\mathrm{F}_{\mathrm{n}}(\mathrm{t}, \beta)=\frac{\mathrm{t}^{2 \mathrm{n}1} \mathrm{e}^{\frac{1}{4}\left(\frac{\beta^{2}}{\mathrm{t}^{2}}+\mathrm{t}^{2}\right)}}{\sqrt{\pi} / 2}\]

Analytical expressions for the initial two moments,
H_{0} and H:sub:` −1`, are given in [Rob81].

The standard freegas kernel is based on the assumption of a constant
free atom cross section. When averaged over the molecular velocity
distribution, this gives a 1/v variation in the effective freegas
cross section at low energies. To approximately account for nuclear
structure effects on the energy dependence of the thermal cross section
(e.g., low energy resonances), the freegas moments are multiplied by
the ratio σ_{s}(E)/σ_{FG}(E), where σ_{s} is the
Doppler broadened scatter cross section processed from ENDF/B data; and
σ_{FG} is the effective freegas cross section,

\[\sigma_{\mathrm{FG}}\left(\mathrm{E}^{\prime}\right)=\frac{\sigma_{\mathrm{free}}}{\mathrm{y}^{2}}\left[\left(\mathrm{y}^{2}+1 / 2\right) \operatorname{erf}(\mathrm{y})+\frac{\mathrm{y} \mathrm{e}^{\mathrm{y}^{2}}}{\sqrt{\pi}}\right]\]

where \(\mathrm{y}^{2}=\mathrm{~A} \frac{\mathrm{E}}{\mathrm{kT}}\).




Submoment expansion of the epithermal scattering source

One difficulty in computing the epithermal scatter source moments is
that the Legendre polynomial in the integrand of (342) and (345) is a function
of both the initial and final lethargy (or energy) of the scattered
neutrons, due to the correlation function G^{(j)}(E,E′). At each
lethargy u this requires that the u′integral be recomputed over all
lower lethargies, for every nuclide and moment. A more efficient
algorithm would be possible if the differential scattering moments
appearing in the integrand could be factored into a product of a
function of u multiplied by a function of u′ such as

(347)\[\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right)=\mathrm{F}^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \quad \mathrm{H}^{(\mathrm{j})}(\mathrm{u})\]

where F and H are the two factors (to be specified later).

If this is done, the ufunction can be factored from the scatter source
integrals, leaving only integrals over the u′function as shown below:

\[\mathrm{S}^{(\mathrm{j})}(\mathrm{u})=\int_{\mathrm{u}^{\prime}} \mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{du}^{\prime}=\mathrm{H}^{(\mathrm{j})}(\mathrm{u}) \int_{\mathrm{u}^{\prime}} \mathrm{F}^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{du}^{\prime}\]

Because the factored integrand does not depend on the variable u, a
running summation over all u′ points can be accumulated and saved as the
calculation sweeps from low to high lethargy. For example, note that the
ℓ = 0 moment in (342) is already separable into a product of u times u′
because P_{0} is equal to one at all values of G. Thus the
isotropic component of the elastic differential scatter rate (per unit
lethargy) from u′ to u is proportional to E/E′, where

\[\mathrm{E}=\mathrm{E}(\mathrm{u})=\mathrm{E}_{\mathrm{ref}} \mathrm{e}^{\mathrm{u}}, \quad \text { and } \quad \mathrm{E}^{\prime}=\mathrm{E}^{\prime}\left(\mathrm{u}^{\prime}\right)=\mathrm{E}_{\mathrm{ref}} \mathrm{e}^{\mathrm{u}^{\prime}}\]

Therefore, the two separable factors in the lowest moment,
S^{(j)}_{0.0}(u′→u), are identified as,

\[\begin{split}\begin{array}{l}
\mathrm{H}(\mathrm{u})=\mathrm{E} /\left(1\alpha^{(\mathrm{j})}\right), \quad \text { and } \\
\mathrm{F}\left(\mathrm{u}^{\prime}\right)=\Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \Psi_{00}\left(\mathrm{u}^{\prime}\right) / \mathrm{E}^{\prime}
\end{array}\end{split}\]

However, the higher order Legendre moments contain the term
P_{ℓ}(G) in the integrand; and the expression for G(E′,E) is a
difference of two terms that depend on both E and E′. A new method
called a “sub‑moment expansion” has been developed for CENTRM that
allows the Legendre polynomials appearing in the differential scatter
moments to be factored into the desired separable form. Each spherical
harmonic moment of the scattering source appears expanded in a series of
factored “sub‑moments.”

The Legendre polynomial of order ℓ is a polynomial containing terms up
to the ℓ^{th} power. Applying the binomial expansion theorem and
some algebraic manipulation, the standard expression for P_{ℓ}
evaluated at “G” can be expressed as

(348)\[\mathrm{P}_{\ell}(\mathrm{G})=\frac{\mathrm{E}^{\prime}}{\mathrm{E}} \times \mathrm{a}_{1}^{\ell} \sum_{\mathrm{K}=\ell}^{\ell} \tilde{\mathrm{g}}_{\ell, \mathrm{K}}(\mathrm{E}) \quad \mathrm{h}_{\mathrm{K}}(\mathrm{E}) \mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}^{\prime}\right)\]

where h_{k}(E)=E^{1+K/2}; and the expansion coefficients
\(\tilde{\mathrm{g}}_{\ell, \mathrm{k}}\) are equal to,
\(\tilde{\mathrm{g}}_{\ell \mathrm{K}}=\frac{\mathrm{g}_{\ell \mathrm{K}}}{N_{\ell} \times \alpha_{1}^{\ell}}\)
where the g_{ℓ,K} (no tilde)
parameters were defined in [WA95] to be:

(349)\[\begin{split}\begin{aligned}
&\mathrm{g}_{, \mathrm{K}}=\frac{\left(1+(1)^{+\mathrm{K}}\right)}{2} \sum_{K^{\prime}=0}^{\frac{K}{2}}(1)^{K^{\prime}} b_{2 K^{\prime}+K},\left(\begin{array}{r}
2 K^{\prime}+K \\
K^{\prime}
\end{array}\right) \quad a_{1}^{K+K^{\prime}} a_{2}^{K^{\prime}} ; \quad \text { for } \quad \mathrm{K} \geq 0\\
&\text { and }\\
&=\left(a_{2} / a_{1}\right)^{K} \quad g_{,K} \quad ; \quad \text { for } \quad K<0
\end{aligned}\end{split}\]

In (348) –(349), the constants b_{m,ℓ} and N_{ℓ} are the standard
Legendre constants and normalization factors, respectively, which are
tabulated in Table 269 for orders through P_{7}; and \(\left(\begin{array}{c}
\mathrm{m} \\
\mathrm{i}
\end{array}\right)=\)
the binomial expansion coefficients^{(20)} \(= \frac{\mathrm{m} !}{(\mathrm{m}\mathrm{i}) ! \quad \mathrm{i} !}\)








Table 269 Constants appearing in Legendre polynomials.
 

The explicit dependence of the constants a_{1} and a_{2} on
the nuclide index j [see (335)] has been suppressed to simplify notation.
For discretelevel inelastic scatter the parameter a_{1} is an
energy dependent function given by (343), but for elastic scatter this
reduces to the constant in (335). Note that the g_{ℓ,K} value is zero unless ℓ and K are both
even or both odd, respectively, so that about half the terms appearing
in the summation of (348) vanish. Table 270 through Table 272 give values
for the submoment expansion coefficients for
several nuclides.

The submoment expansion of the scattering source, including both
elastic and discretelevel inelastic reactions, is obtained by
substituting the expansion of the Legendre polynomial from (348) into
(346), giving

(350)\[\mathrm{S}_{\mathrm{k}}(\mathrm{u})=\sum_{m, j} \sum_{K=} \mathrm{Z}_{, \mathrm{K}}^{(\mathrm{m}, \mathrm{j})}(\mathrm{E}) \mathrm{h}_{\mathrm{K}}(\mathrm{E}) \int_{u_{L O}^{(m, j)}}^{u_{H}^{(m, j)}} \psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \Sigma^{(\mathrm{m}, \mathrm{j})}\left(\mathrm{E}^{\prime}\right) \frac{\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}^{\prime}\right)}{\Delta^{(\mathrm{m}, \mathrm{j})}\left(\mathrm{E}^{\prime}\right)} \mathrm{du}^{\prime}\]

where \(Z_{\ell \mathrm{K}}^{(\mathrm{m}, \mathrm{j})}(E)=a_{1}^{\ell}(E) \frac{\tilde{g}_{\ell, K}^{(m, j)}(E)}{\left(1\alpha^{(j)}\right)}\).
For elastic scatter, the Z coefficients are independent
of energy.

With this approach the scatter source moments in (351) have been further
expanded into a summation of “submoments” identified by index K
(although some of these terms are equal to zero, due to the behavior of
the g_{ℓ,K `coefficients). Each term has the desired factored
form expressed in :eq:`eq7422}; i.e., separable in terms of the variables u and
u′ with

(351)\[\mathrm{H}_{, \mathrm{K}}^{(\mathrm{j})}(\mathrm{u})=\mathrm{Z}_{, \mathrm{K}}^{(\mathrm{j})}(\mathrm{E}) \mathrm{h}_{\mathrm{K}}(\mathrm{E}), \quad \text { and } \quad \mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{)})}\left(\mathrm{u}^{\prime}\right)=\frac{\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}^{\prime}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right)}{\Delta^{(\mathrm{m}, \mathrm{j})}\left(\mathrm{E}^{\prime}\right)}\]

so that the lk_{th} moment of the scatter source can be written
as

(352)\[\mathrm{S}_{\mathrm{k}}(\mathrm{u})=\sum_{m, j} \sum_{K=1} \mathrm{H}_{, \mathrm{K}}^{(\mathrm{j})}(\mathrm{u}) \int_{\left.\mathrm{u}_{\mathrm{LO}}^{(\mathrm{m}, \mathrm{j}}\right)}^{\mathrm{u}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{m}, \mathrm{j})}} \mathrm{F}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}\]


Characteristics and Properties of the SubMoment Expansion

The expansion in (348) becomes numerically unstable for heavy nuclides
(large A), with high Legendre orders. Using double precision arithmetic,
it was found that the accuracy of the expansion begins to break down for
heavy nuclides (A100) if the order of scatter exceeds P_{5};
although the expansion for lighter nuclides (viz, moderators) is very
accurate even for scattering orders as high as P_{7} or more. For
this reason CENTRM has an option to restrict the Legendre expansion to
lower orders for heavy masses, while using the input value of “ISCT” for
lighter nuclides. The restricted Legendre order and mass cutoff value
can be controlled by user input, but the default is P_{0} (i.e.,
isotropic lab scattering) for A>100. Table 273 shows the maximum error
observed in the series representation of Legendre polynomials up to
P_{5}, for selected mass numbers. These values were obtained by
evaluating the series expansion for Pℓ(G(x)) in (348), and
comparing to the exact value computed at 11 equally spaced values for
E/E′. The observed error in the P_{5} polynomial expansion is < 1%
even for heavy materials such as ^{238}U, while nuclides whose
mass is < 100 are computed nearly exactly by the expansion.

Although the accuracy of the submoment expansion is good through
P_{7} scattering in moderators, Legendre expansions above P3 are
not recommended because the number of terms in the expansion increases
rapidly with increasing scattering order, especially for 2D MoC and 1D
cylindrical systems. The number of spherical harmonic moments appearing
in the scattering source depends on the order (L=ISCAT) of the Legendre
expansion used to represent the differential scatter cross section, as
well as on the type of geometry (slab, spherical, cylindrical, or 2D
MoC) used in the transport calculation. The submoment method further
expands each source moment. Table 274 shows the number of moments in
the crosssection expansion, and the number of moments and submoments in
the scatter source expansion, as a function of scatter order and
geometry type. Although the use of cumulative integrals discussed below
allows the submoments to be evaluated rapidly, the large number of
terms becomes prohibative for high scattering orders. Fortunately a
P_{1} Legendre order is sufficient for most selfshielding
calculations, and orders beyond P_{2} should seldom be required
for reactor physics and criticality applications.








Table 270 Coefficients in expansion of Pℓ[G(x)]* for hydrogen (A = 1).
 








Table 271 Coefficients in expansion on Pℓ[G(x)]* for oxygen (A = 16).
 








Table 272 Coefficients in expansion of Pℓ[G(x)]* for U238 (A = 236).
 








Table 273 Fractional error ^{(1)} in series expansion [Eq. (F18.2.25)] of Legendre polynomials.
 








Table 274 Number of moments and submoments as function of scattering order.
 



Scattering moments expressed with cumulative integral operator

It will be convenient to express the scatter source moments in terms of
an integral operator \(\mathbb{C}\), designated here as the “cumulative integral.” The
domain of this operator is the vector space of all integrable lethargy
functions. The operator is defined for an arbitrary domain element
f(u’), at an arbitrary lethargy limit U, to be:

(353)\[\mathbb{C}(\mathrm{f} ; \mathrm{U})=\int_{\mathrm{u}_{0}}^{\mathrm{U}} \mathrm{f}\left(\mathrm{u}^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}\]

where u_{0} is an arbitrary reference point. In implementing this
method in CENTRM, it is convenient to set u_{0}=u_{L};
i.e., the negative lethargy value corresponding to highest energy of the
transition range.

The cumulative integral at some lethargy mesh point u_{n} is
related to the value at the previous lethargy mesh point u_{n−1}
by the expression

(354)\[f\left(\text{f} ; \text{u}_{n}\right) = f\left(\text{f} ; \text{u}_{n1}\right)+\int_{u_{n1}}^{u_{n}} f\left(u^{\prime}\right) d u^{\prime}\]

where u_{n} > u_{n−1}.

Note that only a single panel of integration over the interval
[u_{n−1}, u_{n}] must be performed to update the cumulative
integrals.

The submoment expansion of the scatter source in (350) can be expressed
in terms of the cumulative integral operator as follows:

(355)\[\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}(\mathrm{u})=\sum_{j} \sum_{K=} \mathrm{H}_{, \mathrm{K}}^{(\mathrm{j})}(\mathrm{u}) \times\left[\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; u_{H I}^{(m, j)}\right)\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; u_{L O}^{(m, j)}\right)\right]\]

For elastic scatter the value of \(\mathrm{u}_{\mathrm{LO}}^{(\mathrm{m}, \mathrm{j})}\) is equal to (u−ε^{(j)}),
and \(\mathrm{u}_{\mathrm{HI}}^{(\mathrm{m}, \mathrm{j})}\) is equal to u.



The Boltzmann equation within the PW range

In contrast to the “pseudopointwise” fluxes obtained from the MG
transport calculation, a true PW solution is performed for the
N_{P} lethargy points between DEMAX and DEMIN. The PW solution is
performed within a loop over energy groups: i.e., for each of the
NG_{P} groups in the PW range there is an additional loop over all
lethargy mesh points contained inside the group. This approach
facilitates coupling of the scatter source from the UMR to the PW range
and from the PW and LMR.

Evaluating (326) at each of the N_{P} energy meshpoints in the
PW range gives a system of integrodifferential equations that can be
solved to obtain the PW flux moments, per lethargy, for the
N_{P} energy mesh points in the range DEMAX to DEMIN—which
correspond to the lethargy points, \(\mathrm{U}_{\mathrm{NGU}+1}, \ldots \mathrm{U}_{\mathrm{NGU}+\mathrm{NP}}\).
Again linear variation of
the flux between lethargy points is assumed, to obtain a continuous
spectrum. Substituting (329) into (326), the PW transport equation at mesh
point n is found to be,

(357)\[\Omega \bullet \nabla \Psi_{\mathrm{n}}(\mathrm{r}, \Omega)+\Sigma_{\mathrm{t}, \mathrm{n}}(\mathrm{r}) \Psi_{\mathrm{n}}(\mathrm{r}, \Omega)=\sum_{\mathrm{k}} \frac{2+1}{2} \mathrm{Y}_{\mathrm{k}}(\Omega) \quad \mathrm{S}_{\mathrm{k}, \mathrm{n}}(\mathrm{r})+\mathrm{Q}_{\mathrm{n}}(\mathrm{r}, \Omega)\]

for \(\mathrm{n}=\left(\mathrm{NG}_{\mathrm{U}}+1\right), \ldots .,\left(\mathrm{NG}_{\mathrm{U}}+\mathrm{N}_{\mathrm{P}}\right)\)

where

\[\begin{split}\begin{aligned}
\sum_{\mathrm{t}, \mathrm{n}}(\mathbf{r}) &=\sum_{\mathrm{t}}\left(\mathbf{r}, \mathrm{u}_{\mathrm{n}}\right) \\
\Psi_{\mathrm{n}}(\mathbf{r}, \Omega) &=\Psi\left(\mathbf{r}, \Omega, \mathrm{u}_{\mathrm{n}}\right) \\
\mathrm{S}_{\mathrm{k}, \mathrm{n}}(\mathbf{r}) &=\mathrm{S}_{\mathrm{k}}\left(\mathbf{r}, \mathrm{u}_{\mathrm{n}}\right)
\end{aligned}\end{split}\]

Aside from the definition of the crosssection data, the above equation
appears identical in form to the MG transport equation, and can be
solved with virtually the same algorithm as the MG solution, once the
scatter source moments are determined. The same computer routines in
CENTRM calculate both the MG and PW fluxes. However, a major conceptual
difference between the PW and MG transport equations is that the PW
equation describes a differential neutron balance per unit lethargy at
an energy point, while the MG equation represents an integral balance
over an interval of lethargy points. Although this type of point
solution is not inherently conservative over the intervals defined by
the energy mesh, the particle balance for each interval has been found
to be very good. It should also be noted that exact particle
conservation is not a strict requirement for this type of application
where flux spectra rather than particle balances are primarily of
interest.

In the PW range the scatter source is composed of (a) MGtoPW scatter
from the UMR and possibly upscatter from the LMR if the PW range extends
into thermal, and (b) PWtoPW scatter from points in the PW range. The
submoment expansion method described previously is used in CENTRM to
provide an efficient method of evaluating the PWtoPW downscatter
source for the epithermal range, which includes most of the resolved
resonances.


Scattering sources for the PW range

In the case of elastic scatter from nuclide “j,” only the lethargy
interval below u_{n} −ε^{(j)} can scatter to a lethargy
point u_{n} in the PW range. If u_{n} −ε^{(j)} is
negative, then some portion of the source at u_{n} is due to
UMRtoPW from energies above DEMAX, since zerolethargy is equal to the
top energy of the PW range. Otherwise, the elastic source is entirely
PWtoPW.

For any given nuclide j, the lowest lethargy in the UMR range that
contributes to the elastic scatter source in the PW range is equal to
−ε^{(j)}. Let “jL” represent the lightest nonhydrogen nuclide
(i.e., having the smallest A value greater than unity) in the system.
The associated fractional energy loss for this material is indicated as
α_{L}, so that the highest energy neutron in the UMR range that
can scatter into the PW range from an elastic collision with any
nonhydrogenous moderator will have an energy equal to
DEMAX/α_{L}. The corresponding lethargy is equal to be the
negative value −ε^{(L)}, or −ln(1/α_{L}). The value of
−ε^{(L)} is actually adjusted in CENTRM to coincide with the
immediately preceding multigroup boundary, which has a lethargy value
designated as u_{L}. The interval of negative lethargy in the UMR
between u_{L} and 0 has been defined previously to be transition
range, because the elastic slowingdown source from this interval
provides a transition between the UMR and PW solutions, respectively.
The transition range always contains an integer number of groups,
corresponding to MGTOP to MGHI. The total downscatter source from the
UMR to lethargy u_{n} is composed of elastic and inelastic
contributions from the transition range between [u_{L},0]; and
contributions from the “high” energy range from lethargies below
u_{L}. The high contribution comes from inelastic and hydrogen
elastic reactions in the energy interval above the transition range.

The downscatter source at u_{n} in the PW range can thus be
expressed as the sum of three distinct contributions — S_{HI},
S_{Tr}, and S_{PW} —, that correspond to scatter from the
high region of the UMR, the transition region of the UMR, and the PW
ranges, respectively. The source moments appearing in (357) can thus be
expressed as:

(358)\[\begin{split} \begin{aligned}
\mathrm{S}_{\mathrm{k}, \mathrm{n}}(\mathrm{r}) &=\mathrm{S}_{\mathrm{k}, \mathrm{HI}}\left(\mathrm{r}, \mathrm{u}_{\mathrm{n}}\right)+\mathrm{S}_{\mathrm{k}, \mathrm{Tr}}\left(\mathrm{r}, \mathrm{u}_{\mathrm{n}}\right)+\mathrm{S}_{\mathrm{k}, \mathrm{PR}}\left(\mathrm{r}, \mathrm{u}_{\mathrm{n}}\right) \\
&=\int_{\infty}^{\mathrm{u}_{\mathrm{L}}} \mathrm{S}_{\mathrm{k}}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{d} \mathrm{u}^{\prime}+\int_{\mathrm{u}_{\mathrm{L}}}^{0} \mathrm{S}_{\mathrm{k}}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{du}^{\prime}+\int_{0}^{\mathrm{u}_{\mathrm{n}}} \mathrm{S}_{\mathrm{k}}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{du}^{\prime}
\end{aligned}\end{split}\]



Downscatter source from high region of the UMR to the PW range (SHI)

The high region of the UMR corresponds to groups 1 through MGTOP1. The
MGtoPW scattering source (S_{HI}) from high energy region
originates in the energy range above DEMAX/α_{L}; i.e., lethargies
below u_{L} (see Fig. 207). In this region, inelastic
reactions may scatter neutrons to the PW range; but due to the
definition of u_{L}, the only elastic reactions that scatter to
the PW range are due to hydrogen. Therefore in general, the MG matrices
describing scatter from groups in high region to groups in the PW range
correspond to discrete and continuum inelastic reactions, and elastic
scatter from hydrogen. If g′ is an arbitrary group in the UMR range
above the transition interval and g is a fixed group interval in the
PW range, then the rate that neutrons scatter from all groups g′ in the
high region to all energy points in g, for a given direction Ω, is
obtained from the usual expression for MGtoMG transfers, and is equal
to

\[\mathrm{S}_{\mathrm{g}}(\mathrm{r}, \Omega)=\sum_{\mathrm{k}} \frac{2+1}{2} \quad \mathrm{Y}_{\mathrm{k}}(\Omega) \mathrm{S}_{\mathrm{k}, \mathrm{g}}\]

where MGLO > g > MGHI, and the MG source moments are,

(359)\[\mathrm{S}_{\mathrm{k}, \mathrm{g}}=\sum_{\mathrm{g}^{\prime}=1}^{\text {MGTOP1 }} \Sigma_{, \mathrm{g}^{\prime} \rightarrow \mathrm{g}} \Psi_{\mathrm{k}, \mathrm{g}^{\prime}}\]

while (359) gives the moments of the overall scatter rate from all groups
in the high range into the entire PW group g, it is necessary to
determine how the group source should be distributed over the PW energy
mesh contained within the group; i.e., it is desired to extract the PW
source moments, from the group moments by applying some “intragroup”
distribution H_{ℓk,g}(E) such that,

(360)\[\mathrm{S}_{\mathrm{k}, \mathrm{HI}}(\mathrm{u})=\mathrm{S}_{\mathrm{k}, \mathrm{g}} \quad \mathrm{H}_{\mathrm{k}, \mathrm{g}}(\mathrm{E}), \quad \text { for } \mathrm{u}(\mathrm{E}) \varepsilon \operatorname{group} \mathrm{g}\]

The intragroup distribution has units of “per unit lethargy,” and its
integral over the group is normalized to unity. This form of the scatter
source preserves the MG moments S_{ℓk,g}, whenever
S_{ℓk,HI}(u) is integrated over group g, insuring that the
correct number of neutrons (as determined from the UMR calculation) will
always be transferred from the high range into the PW group. Only the
distribution within the group is approximate.

Recall that the scatter source of concern here is due only to elastic
scatter from hydrogen and inelastic scatter from all other materials. In
the case of swave elastic scatter from hydrogen, the P_{0} and
P_{1} moments per unit lethargy, respectively, can be rigorously
expressed in the form of (360) with

(361)\[ \begin{array}{lllll}
\mathrm{H}_{0} & \propto \mathrm{E} & , & \text { and } & \mathrm{H}_{1} & \propto \mathrm{E}^{3 / 2}
\end{array}\]

These expressions can be inferred directly from the moments of the
scatter kernel in (341) . The higher order scatter moments for hydrogen
have a somewhat more complicated form containing sums of energy
functions; but since these moments are usually less important than the
first two moments, a less rigorous treatment of their intragroup
distribution is used. The intragroup distribution due to inelastic
scatter depends on the Q values for the individual levels, and these are
not available on the multigroup libraries. Fortunately, the scatter
source in the PW range is not very sensitive to the assumed intragroup
distribution for inelastic scatter, as long as the total inelastic
source for the group is computed correctly. As a reasonable tradeoff
between rigor and complexity, the high energy component of the UMRtoPW
scatter source is approximated using H_{0} for the intragroup
distribution of all P_{0} moments, and H_{1} for all higher
order moments. This approximation produces the correct intragroup
variation for the lowest two moments of the hydrogen scatter source, but
the higher order moments of hydrogen and the inelastic scatter source
are not distributed exactly throughout the group. However, the
integrated source moments are correct in all cases. Again, it should be
emphasized that the approximations discussed here only apply to the
UMRtoPW component designated as S_{HI}, which comes from
reactions above the transition range (energies above
E_{HI}/α_{L}). This is often a small contribution to the
overall PW source term.



Scattering sources from UMR transition region and epithermal PW range

Most coupling between the UMR and the PW range is due usually to elastic
scatter from energies immediately above DEMAX. The contribution to the
PW source due to downscatter source from this transition range has been
designated S_{Tr}(u:sub:n). The other component of the PW
source, S_{PW}(u_{n}), accounts for the scattering source
coming from all lethargies lower than u_{n} in the PW range. It is
convenient to combine the two sources together as the PW epithermal
source called “S_{Ep},” which has an lk_{th} moment given by
(347),

(362)\[\begin{split} \begin{aligned}
S_{k, E p} &=\int_{\mathrm{u}_{\mathrm{L}}}^{\mathrm{u}_{\mathrm{n}}} \mathrm{S}_{\mathrm{k}}\left(\mathrm{r}, \mathrm{u}^{\prime} \rightarrow \mathrm{u}\right) \mathrm{du}^{\prime} \\
&=\sum_{\mathrm{j}} \sum_{\mathrm{K}=} \mathrm{Z}_{, \mathrm{k}}^{(\mathrm{j})} \mathrm{h}_{\mathrm{K}}(\mathrm{E}) \int_{\mathrm{u}_{\mathrm{n}}\varepsilon^{(\mathrm{j})}}^{\mathrm{u}_{\mathrm{k}}} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \quad \mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}^{\prime}\right) \mathrm{du}^{\prime}
\end{aligned}\end{split}\]

This is done because CENTRM uses the submoment expansion technique to
compute both the PWtoPW epithermal source from the PW range as well as
the MGtoPW source from the transition range of the UMR. Note that
elastic scattering from the transition range only impacts the PW scatter
source at the initial mesh points in the PW range; i.e., those contained
in the interval 0 < u_{n} <ε:sup:(j), for nuclide j. Beyond
these mesh points the elastic scatter source is due only to PWtoPW
scatter, as illustrated in Fig. 207.

The epithermal elastic source at u_{n}, coming from the range
u_{L} to u_{n}, is expressed as an integral over the
immediately preceding lethargy mesh interval from u_{n−1} to
u_{n} plus the integral over the remaining lethargy interval, as
illustrated in Fig. 208. The former integral is designated as
I(u_{n−1},u_{n}) and the latter as
I(u_{L},u_{n−1}), so that

\[\mathrm{S}_{\mathrm{k}, \mathrm{E}}\left(\mathrm{u}_{\mathrm{n}}\right)=\mathrm{I}\left(\mathrm{u}_{\mathrm{n}1}, \mathrm{u}_{\mathrm{n}}\right)+\mathrm{I}\left(\mathrm{u}_{\mathrm{L}}, \mathrm{u}_{\mathrm{n}1}\right)\]


The lethargy mesh in CENTRM is constrained such that the maximum
lethargy gain in an elastic reaction (ε^{(j)}) is always greater
than the maximum mesh interval size, which insures that
I(u_{n−1},u_{n}) always includes the full panel from
u_{n−1} to u_{n}. In the above and subsequent equations the
explicit dependence of S_{Ep} on independent variables other than
lethargy is not shown for notational convenience. The integral
I(u_{n−1},u_{n}) is evaluated approximately by applying the
trapezoidal rule, which leads to,

(363)\[\mathrm{I}\left(\mathrm{u}_{\mathrm{n}}, \mathrm{u}_{\mathrm{n}1}\right)=\int_{\mathrm{u}_{\mathrm{n}1}}^{\mathrm{u}_{\mathrm{n}}} \mathrm{S}_{\mathrm{k}}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}_{\mathrm{n}}\right) \mathrm{du}^{\prime} \sim \frac{\left[\mathrm{S}_{\mathrm{k}}\left(\mathrm{u}_{\mathrm{n}} \rightarrow \mathrm{u}_{\mathrm{n}}\right)+\mathrm{S}_{\mathrm{k}}\left(\mathrm{u}_{\mathrm{n}1} \rightarrow \mathrm{u}_{\mathrm{n}}\right)\right]}{2} \times\left(\mathrm{u}_{\mathrm{n}}\mathrm{u}_{\mathrm{n}1}\right)\]

Using the submoment expansion from (350), (363) can be written for elastic scatter as

(364)\[\mathrm{I}\left(\mathrm{u}_{\mathrm{n}1}, \mathrm{u}_{\mathrm{n}}\right)=\Sigma_{\mathrm{n} \rightarrow \mathrm{n}} \Psi_{\mathrm{k}, \mathrm{n}}+\sum_{\mathrm{K}} Z_{\mathrm{K}}^{(\mathrm{j})} \mathrm{h}_{\mathrm{K}}\left(\mathrm{E}_{\mathrm{n}}\right) \mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{n}1}\right) \Sigma^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{n}1}\right) \Psi_{\mathrm{k}, \mathrm{n}1} \frac{\Delta \mathrm{u}_{\mathrm{n}1}}{2} .\]

The first term on the right side of (364) corresponds to the
“withinpoint” component of elastic scatter from u_{n} to
u_{n}, which only occurs for straight ahead scatter
(μ_{0}=1). The withinpoint cross section is defined as,

(365)\[\Sigma_{\mathrm{n} \rightarrow \mathrm{n}}=\frac{\Delta \mathrm{u}_{\mathrm{n}1}}{2} \sum_{\mathrm{j}} \frac{\Sigma_{\mathrm{n}}^{(\mathrm{j})}}{\left(1\alpha^{(\mathrm{j})}\right)} .\]

In deriving this term the following relation has been used for each nuclide:

(366)\[\sum_{\mathrm{K}} \mathrm{Z}_{\mathrm{K}}=\frac{1}{1\alpha} .\]

The I(u_{L},u_{n−1}) portion of the integral in (362) is
equal to

(367)\[\mathrm{I}\left(\mathrm{u}_{\mathrm{L}}, \mathrm{u}_{\mathrm{n}1}\right)=\sum_{j} \sum_{K=} \mathrm{Z}_{, \mathrm{K}}^{(\mathrm{j})} \mathrm{h}_{\mathrm{K}}\left(\mathrm{E}_{\mathrm{n}}\right) \int_{\mathrm{U}_{\mathrm{n}}\varepsilon^{(j)}}^{\mathrm{u}_{\mathrm{n}1}} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \quad \mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime} .\]

Note that the lower lethargy limit of the integral is restricted to
u_{n} − ε^{(j)}, since this is the maximum limit of lethargy
that can scatter to u_{n} in an elastic reaction. In terms of the
cumulative integral operator, the integral in (367) over the interval
[u_{n} − ε^{(j)}, u_{n−1}] is equal to

(368)\[\int_{\mathrm{u}_{\mathrm{n}}\varepsilon(\mathrm{j})}^{\mathrm{u}_{\mathrm{n}}1} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{h}_{\mathrm{K}}\left(\mathrm{E}^{\prime}\right)^{1} \mathrm{du}^{\prime}=\left[\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}1}\right)\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}}\varepsilon^{(\mathrm{j})}\right)\right]\]

where F_{ℓk,K} has been defined in (351). In order to evaluate (368)
it is necessary to determine the cumulative integral values at
u_{n−1} and u_{n} − ε^{(j)}. The lethargy u_{n−1}
will always correspond to a mesh point value, but in general
u_{n} − ε^{(j)} can fall between mesh points. Evaluation of
the cumulative integrals at an arbitrary limit such as
u_{n} − ε^{(j)} is performed in CENTRM by interpolation of
previously calculated values stored for all the mesh points below
u_{n} during the transport calculation at lower lethargies. The
interpolated value of the cumulative integral at
u_{n} − ε^{(j)} that is subtracted in (368) is called the
“excess integral” in CENTRM. At each lethargy point, excess
integrals must be found as a function of space, nuclide, moment, and
submoment. Also note that for some initial mesh points
(i.e., u_{n} < ε^{(j)}) the value u_{n} − ε^{(j)} can
be negative, indicating that a portion of the PW scatter source at
u_{n} is due to elastic scattering from the negative lethargy
range above DEMAX. This means that cumulative integrals must be known
for mesh intervals in the transition as well as in the PW range. Values
of the cumulative integrals at all points within the transition range
are first computed from the results from the UMR calculation, prior to
the PW transport calculation (but after the UMR calculation). Additional
cumulative integrals are then calculated successively during the
PW transport solution at all mesh points and are stored as the
calculation proceeds from low to high lethargy. Thus in evaluating
S_{ℓk,Ep}(u_{n}), the cumulative integrals at every space
interval already will have been stored at all energy points up to (n−1),
in an array called CUM^{(j)}_{ℓk,K}, for each nuclide j,
moment ℓk, and submoment K:

(369)\[\mathrm{CUM}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})}=\left\{\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}^{\prime}}\right), \quad \mathrm{n}^{\prime}=1, \mathrm{n}1\right\}\]

so that the excess integral values can be interpolated from the above
array. The first N_{Tr} elements of the array
CUM^{(j)}_{ℓk,K} correspond to lethargy points in the
transition range, and the remainder are in the PW range, where




















N_{Tr}  =  G_{U} −
g_{Tr} + 1; 

g_{Tr}  =  MGTOP, the highest
energy group in the
transition range;
(i.e., the group
whose high energy
boundary corresponds
to u_{L}); 
G_{U}  =  Lowest energy group
in the transition
range. 


Elastic cumulative integrals contained in array
CUM^{(j)}_{ℓk,K} are calculated at each lethargy point
u_{n} with the expression:

(370)\[\begin{split} \begin{array}{l}
f\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}}\right)=\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}1}\right)+\int_{\mathrm{u}_{\mathrm{n}1}}^{\mathrm{u}_{\mathrm{n}}} \mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{d} u^{\prime} \\
=\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}1}\right)+\int_{\mathrm{u}_{\mathrm{n}}1}^{\mathrm{u}_{\mathrm{n}}} \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \quad \mathrm{h}_{\mathrm{K}}\left(\mathrm{E}^{\prime}\right)^{1} \mathrm{du}^{\prime}
\end{array}\end{split}\]

After completing the calculation of PW angular fluxes and moments at
u_{n} the integral over the most current lethargy panel
[u_{n−1},u_{n}] is evaluated with the trapezoidal
approximation, resulting in an updated cumulative integral array
containing the value at lethargy u_{n}:

(371)\[\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}}\right) ; \quad \mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}1}\right)+\Delta \mathrm{u}_{\mathrm{n}1} \frac{\left[\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{n}1}\right) \Sigma_{\mathrm{n}1}^{(\mathrm{j})} \Psi_{\mathrm{k}, \mathrm{n}1}+\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{n}}\right) \Sigma_{\mathrm{n}}^{(\mathrm{j})} \Psi_{\mathrm{k}, \mathrm{n}}\right]}{2} ,\]

where the cumulative integrals at the preceding mesh point are known
from the previous calculation, and the flux moments
Ψ_{ℓk}(u_{n}) are determined from the transport calculation
at the current lethargy point. Only a single panel of integration is
required to update the cumulative integrals, significantly reducing the
amount of computation compared to recomputing the entire summation again
at each new energy point. The integration is performed rapidly with the
trapezoidal approximation, which should be accurate since the energy
mesh is defined to reproduce the macroscopic cross sections linearly
between mesh points. In order to avoid loss of numerical significance,
the set of stored cumulative integrals is periodically “renormalized,”
by translating to a new reference lethargy point (recall that only the
differences of cumulative integrals is needed).

Elastic cumulative integrals for the transition range are calculated
with a slightly different expression, using MG flux moments obtained in
the UMR calculation. Because the transition interval is part of the UMR,
it is convenient to evaluate cumulative integrals at lethargy values
corresponding to group boundaries. This requires approximating the
energy distribution of the flux spectrum within each group in the
transition range. To evaluate the cumulative interval in the transition
range of some nuclide j, the scalar flux per energy (at a given space
location) within a transition group is approximated as:
Φ(E) = M^{(j)}/E, where M^{(j)} is a normalization constant
defined so that the MG outscatter rate (i.e., slowingdown density) from
the group is preserved. It can be shown that this normalization
condition requires that

(372)\[\mathbf{M}^{(j)}=\frac{\left[\Sigma_{t, g^{\prime}}^{(j)}\Sigma_{a, g^{\prime}}^{(j)}\Sigma_{g^{\prime}}^{(j)},\right] \Delta u_{g^{\prime}}}{\xi^{(j)} \Sigma_{s, g^{\prime}}^{(j)}} \times\left[\phi_{g^{\prime}} / \Delta u_{g^{\prime}}\right]\]

where ξ is the average lethargy gain in an elastic reaction and
Σ_{g′g′}, is the withingroup MG scatter cross section. Thus the
scalar flux per unit lethargy used to evaluate cumulative integrals of
nuclide j is:

\[\phi\left(u^{\prime}\right)=M^{(j)} ; \quad \text { for } u^{\prime} \varepsilon g^{\prime}, \text { and } g^{\prime} \varepsilon \text { Transition region of } U M R\]

Withingroup energy spectra for the higher order fluxmoments could be
approximated in similar manner by preserving the higher order Legendre
moments of the slowingdown density, but CENTRM simply uses the same
form in (372) for all flux moments, so that in general the withingroup
energy distribution for any ℓk_{th} moment in the transition range
is approximated as,

(373)\[\Psi_{k}\left(u^{\prime}\right)=\frac{\left[\sum_{t, g^{\prime}}^{(j)}\Sigma_{a, g^{\prime}}^{(j)}\Sigma_{g^{\prime}}^{(j)}\right] \Delta u_{g^{\prime}}}{\xi^{(j)} \sum_{s, g^{\prime}}^{(j)}} \times\left[\Psi_{k, g^{\prime}} / \Delta u_{g^{\prime}}\right]\]

for u′εg′, and g′ε transition region of UMR. Therefore the following
integrals can be evaluated:

(374)\[\int_{\mathrm{ug}^{\prime}}^{\mathrm{u}_{\mathrm{g}^{\prime}+1}} \mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}^{\prime}\right) \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{du}^{\prime}=\frac{\Sigma_{\mathrm{r}, \mathrm{g}^{\prime}}^{(\mathrm{j})} \Delta \mathrm{u}_{\mathrm{g}^{\prime}}}{\xi^{(\mathrm{j})}} \frac{\Psi_{\mathrm{k}, \mathrm{g}^{\prime}}}{\Delta \mathrm{u}_{\mathrm{g}^{\prime}}} \int_{\mathrm{ug}^{\prime}}^{\mathrm{u}_{\mathrm{g}^{\prime}+1}} \mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}^{\prime}\right) \mathrm{du}^{\prime} .\]

Integration of the h_{k}^{−1} function is performed
analytically to give the cumulative integral at any group boundary
u_{g} in the transition range:

(375)\[\begin{split} \begin{array}{l}
f\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{g}}\right) ; \sum_{\mathrm{g}^{\prime}=\mathrm{g}}^{\mathrm{G}_{\mathrm{U}}} \frac{\Sigma_{\mathrm{r}, \mathrm{g}^{\prime}}^{(\mathrm{j})}{\xi^{(\mathrm{j})}} \mathrm{g}_{\mathrm{g}^{\prime}}}{\frac{\Psi_{\mathrm{k}, \mathrm{g}^{\prime}}}{\Delta \mathrm{u}_{\mathrm{g}^{\prime}}} \times\left[\frac{2}{\mathrm{~K}+2}\right]}\left[\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{g}^{\prime}+1}\right)\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{g}^{\prime}}\right)\right] \\
\mathrm{g}=\mathrm{g}_{\mathrm{Tr}}, \quad \mathrm{g}_{\mathrm{Tr}+1}, \quad \mathrm{G}_{\mathrm{U}}
\end{array}\end{split}\]

(375) is used to obtain the initial N_{Tr} values of the
cumulative integrals, corresponding to the transition range. If the
lower limit of the integral in (368) is negative, then the cumulative
integral at u_{n} − ε^{(j)} is interpolated from among the set
of N_{Tr} tabulated values generated by (375); otherwise it is
interpolated from the values that were computed with (371). The following
algorithm is used to interpolate cumulative integrals for negative
lethargy arguments (i.e., in the transition range ):

(376)\[\begin{split} \begin{array}{l}
f\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}\right)=\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{g}}\right)+\frac{\mathrm{h}_{\mathrm{K}}^{1}(\mathrm{E})\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{g}}\right)}{\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{g}+1}\right)\mathrm{h}_{\mathrm{K}}^{1}\left(\mathrm{E}_{\mathrm{g}}\right)} \\
\times\left[\mathrm{f}\left(F_{k, K}^{(j)} ; u_{g+1}\right)\mathrm{f}\left(F_{k, K}^{(j)} ; u_{g}\right)\right]
\end{array}\end{split}\]

.or u(E) ε g; and g ε transition range of UMR.

Because the energy mesh in the PW range is very fine, simple linear
interpolation of the cumulative integrals is used for positive lethargy
arguments.

The complete epithermal elastic scatter source S(r,Ω,u_{n})
appearing in (357) at any mesh point u_{n} corresponds to a
spherical harmonic expansion using the previously derived moments of
S_{HI} and S_{Ep}. This angular scatter source is equal
to,

(377)\[\begin{split} \begin{array}{l}
\mathrm{S}\left(\mathbf{r}, \Omega, \mathrm{u}_{\mathrm{n}}\right)=\Sigma_{\mathrm{n} \rightarrow \mathrm{n}} \Psi_{\mathrm{n}}(\mathbf{r}, \Omega) \\
+\sum_{\mathrm{k}} \frac{2+1}{2} \mathrm{Y}_{\mathrm{k}}(\Omega)\left\{\mathrm{H}\left(\mathrm{E}_{\mathrm{n}}\right) \sum_{\mathrm{g}^{\prime}=1}^{\mathrm{g}_{\mathrm{Tr}}1} \sum_{, \mathrm{g}^{\prime} \rightarrow \mathrm{g}} \psi_{\mathrm{k}, \mathrm{g}^{\prime}}\right. \\
\left.+\sum_{j} \sum_{\mathrm{K}} Z_{\mathrm{K}}^{(\mathrm{j})} \mathrm{h}_{\mathrm{K}}\left(\mathrm{E}_{\mathrm{n}}\right)\left[0.5 \Delta \mathrm{u}_{\mathrm{n}1} \mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{n}1}\right)+\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}1}\right)\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{n}}\varepsilon^{(\mathrm{j})}\right)\right]\right\}
\end{array}\end{split}\]

The above expression was written explicitly for the case of elastic
scatter; however, the discrete level inelastic PW source can be
incorporated with little modification. The only changes are that
additional cumulative integral terms corresponding to each inelastic
level will appear in (377); the cumulative integrals for the inelastic
levels must be computed by integrating the more general expression in
(351); and the lethargy arguments for the inelastic cumulative integrals
are the generalized lethargy limits u_{LO} and u_{HI}
defined in Submoment expansion of the epithermal scattering source and [Wil00].

Note that (377) contains the term Σ_{n→n} Ψ_{n}(r,Ω) which
can be subtracted from both sides of the transport equation in Eq. to
give a slightly altered form of the PW transport equation that contains
a modified scatter source and a modified total cross section. The
modified source component is identical to the expression in (377) with
the withinpoint term Σ_{n→n} Ψ_{n}(r,Ω) removed. The
modified total cross section, represented by \(\Sigma_{\mathrm{t}, \mathrm{n}}\) has the
appearance of a “transport‑corrected” cross section given below:

(378)\[\Sigma_{\mathrm{t}, \mathrm{n}}=\Sigma_{\mathrm{t}, \mathrm{n}}\Sigma_{\mathrm{n} \rightarrow \mathrm{n}}\]

An interesting and significant consequence of this operation is that the
right side of Eq. no longer contains the unknown flux
Ψ_{n}(r,Ω) since the withinpoint term is eliminated. The
resulting modified transport equation has the same form as a purely
absorbing medium with a known source term; and thus can be solved
without requiring scattersource iterations in the epithermal range.
However, iterations may still be required for cell cases with two
reflected or albedo boundary conditions.



PW thermal scatter source

There are significant differences in the CENTRM epithermal and thermal
PW transport solutions. In the epithermal range neutrons can only lose
energy in scattering reactions, so that a single sweep from high to low
energy (i.e., low to high lethargy) is required in the solution. On the
other hand, since low energy neutrons may gain as well as lose energy in
scattering reactions, outer iterations are required to converge the
thermal scattering source. Furthermore, the PW scatter kernels
Σ_{ℓ}(u′→u) in the epithermal range represent twobody
interactions (such as elastic and discretelevel inelastic reactions)
between a neutron and a stationary nucleus. The simple kinematic
relations for these cases allow the efficient submoment expansion
method to be utilized in computing scattering source moments. Thermal
scattering reactions are not two body reactions, but rather represent an
effective average over the molecular velocity distribution; thus, there
is no simple kinematic relationship between neutron energy loss and the
angle of scatter relative to its initial direction. In solving the
transport equation for thermal neutrons, the scatter source at lethargy
u_{n} is approximated as a summation over the “N” mesh points in
the thermal range,

(379)\[\int_{\text {thermal }} \Sigma^{(\mathrm{j})}\left(\mathrm{u}^{\prime} \rightarrow \mathrm{u}_{\mathrm{n}}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}=\sum_{\mathrm{m}=1}^{\mathrm{N}} \mathrm{W}_{\mathrm{m}} \Sigma^{(j)}\left(\mathrm{u}_{\mathrm{m}} \rightarrow \mathrm{u}_{\mathrm{n}}\right) \Psi_{\mathrm{k}}\left(\mathrm{u}_{\mathrm{m}}\right)\]

where

m = 1 is the thermal/epithermal boundary point;

m = N is the lowest thermal energy point; and

W_{m} are standard quadrature weights for trapezoidal integration
with N1 lethargy panels:


\[\begin{split}\begin{aligned}
\mathrm{W}_{\mathrm{m}}=& 0.5 \times\left(\Delta \mathrm{u}_{\mathrm{m}}+\Delta \mathrm{u}_{\mathrm{m}+1}\right) \quad ; \quad \text { for } \mathrm{m}=2,3, \ldots \mathrm{N}1 \\
& 0.5 \times \Delta \mathrm{u}_{\mathrm{m}} \quad ; \quad \text { for } \mathrm{m}=1 \quad \text { or } \quad \mathrm{N}
\end{aligned}\end{split}\]

Pointtopoint crosssection moments in the thermal range are computed
from the freegas or bound kernels evaluated at the desired initial
(u_{m}) and final (u_{n}) lethargy mesh points. For a given
outer iteration, the summation in (379) is evaluated using the most
recently computed flux moments. In many instances the main purpose of
the CENTRM calculation will be to obtain a PW spectrum for resonance
selfshielding calculations. In these cases the thermal flux does not
have to be converged very tightly to obtain a reasonable thermal
spectrum for selfshielding low energy resonances, so that only a few
outer iterations are typically employed.

An additional complication in the thermal calculation is that inner
iterations are necessary to converge the “withinpoint” (no energy loss)
contribution of the thermal scattering source, due to the presence of
PW flux moments at lethargy point m = n. No inner iterations are
required to converge the withinpoint elastic scatter term in the
epithermal PW calculation because there can be no change in the neutron
direction if there is no energy loss, unlike the thermal range.

A spacedependent rebalance calculation for the entire thermal energy
band is performed between outer iterations in order to speed up
convergence of the solution. Reaction rates and leakage values appearing
in the thermalband rebalance equation are obtained by integrating
PW values over the thermal range. Other acceleration techniques, such as
overrelaxation, extrapolation, and renormalization, are also employed.



Downscatter source from the epithermal PW range to the LMR

MG transport calculations performed in the energy range below DEMIN,
which includes the thermal energy range, are coupled to the epithermal
PW range transport calculations by the slowing down source. The
epithermal PWtoLMR scatter source represents the contribution to the
multigroup source in some fixed group g contained in the LMR, from
scatter reactions in the epithermal range above DEMIN. The lethargy
value corresponding to the energy DEMIN (i.e., the bottom energy of the
PW range) will be indicated as u_{PW}, thus
u_{PW} = ln(DEMAX/DEMIN); while the lethargy corresponding to the
thermal energy boundary will be designated as u_{TH}. The cutoff
lethargy for the epithermal PW range will correspond to:
u_{cut} = min(u_{PW},u_{TH}). If there is no PW thermal
calculation in CENTRM, then u_{cut} = u_{PW}; otherwise,
u_{cut} = u:sub:TH. For a given nuclide j, the lowest lethargy
in the epithermal PW range from which a neutron can scatter elastically
into the LMR is equal to (u:sub:cut − ε^{(j)}). If the value of
(u − ε^{(j)}) is greater than u_{cut}, then an elastic
collision with nuclide j cannot moderate an epithermal neutron from the
PW range to u. Therefore in general for a given material zone, only a
limited number of nuclides (possibly none) and a limited portion of the
epithermal PW energy range may be able to scatter neutrons elastically
to any particular group in the LMR. Utilizing the elastic scatter kernel
and applying a submoment expansion to the resulting expression, the
source moment describing scatter from the PW epithermal range to a
lethargy u in the LMR is found to be

\[\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}(\mathrm{u})=\sum_{\mathrm{K}} Z_{\mathrm{K}}^{(\mathrm{j})} \quad \mathrm{h}_{\mathrm{K}}(\mathrm{E}) \int_{\mathrm{U}\varepsilon^{6}}^{\mathrm{u}_{\mathrm{Cut}}} \mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})}\left(\mathrm{u}^{\prime}\right) \mathrm{d} \mathrm{u}^{\prime}\]

where u in group g; and g ε LMR.

The integral in the above expression can be evaluated from cumulative
integrals stored during the epithermal PW transport calculation. Thus
the source moment per unit lethargy at u in the LMR range, due to
epithermal scattering from nuclide j, can be written as,

(380)\[\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}(\mathrm{u})=\sum_{\mathrm{K}} \mathrm{Z}_{\mathrm{K}}^{(\mathrm{j})} \quad \mathrm{h}_{\mathrm{K}}(\mathrm{E}) \quad\left[\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}_{\mathrm{cut}}\right)\mathrm{f}\left(\mathrm{F}_{\mathrm{k}, \mathrm{K}}^{(\mathrm{j})} ; \mathrm{u}\varepsilon^{(\mathrm{j})}\right)\right] ,\]

for uε group g and u−ε(j) < u_{PW}.

The source per unit lethargy in Eq. is integrated over the “sink group”
g in the LMR to determine the desired MG scatter source moment due to
reactions in the epithermal PW range. The actual integral over group g
is performed numerically by introducing a threepoint (two panel)
integration mesh within the group, as follows:

\[\begin{split}\begin{aligned}
\mathrm{u}_{\mathrm{I}} &=\text { initial integration point in group } \mathrm{g}=\text { lethargy at top energy of group } \mathrm{g}=\mathrm{u}_{\mathrm{g}} \\
\mathrm{u}_{\mathrm{F}}^{(\mathrm{j})} &=\text { final integration point in group } \mathrm{g} \\
&=\operatorname{MIN}\left\{\mathrm{u}_{\mathrm{g}+1} ; \mathrm{u}_{\mathrm{cut}}+\varepsilon^{(j)}\right\}, \quad \text { where } \mathrm{u}_{\mathrm{g}+1}=\text { lethargy at bottom of energy of group } \mathrm{g} \\
\mathrm{u}_{\mathrm{A}}^{(\mathrm{j})} &=\text { middle integration point in group } \mathrm{g}=0.5\left(\mathrm{u}_{\mathrm{F}}^{(\mathrm{j})}+\mathrm{u}_{\mathrm{I}}\right)
\end{aligned}\end{split}\]

Note that the final and middle points of integration
(i.e., u_{F}^{(j)} and u:sub:A^{(j)}) may be nuclide
dependent; and if (u_{I} − ε^{(j)}) > u_{cut}, then
nuclide j does not contribute to the pointwisetoLMR scatter source in
g. Applying the twopanel Simpson’s approximation for integration over
group g results in

(381)\[\mathrm{S}_{\mathrm{k}, \mathrm{g}}^{(\mathrm{)}}=\Delta^{(\mathrm{j})} / 3\left[\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{I}}\right)+4 \mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{A}}^{(\mathrm{j})}\right)+\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{F}}^{(\mathrm{j})}\right)\right]\]

where Δ^{(j)} = 0.5(u_{F}^{(j)} − u_{I}).

The values for S_{ℓk}^{(j)}(u_{I}),
S_{ℓk}^{(j)}(u_{A}^{(j)}), and
S_{ℓk}^{(j)}(u_{F}^{(j)}) in (381) are obtained
by evaluating (380) at the lethargy values u_{I},
u_{A}^{(j)}, and u_{F}^{(j)}, respectively. Use
of more than two panels for the group integration was found to have an
insignificant impact.

The complete epithermal PWtoLMR source in group g is finally obtained
by summing (380) over all nuclides and then substituting the spherical
harmonic moments into the Legendre expansion of the MG scatter source,
resulting in

(382)\[\mathrm{S}_{\mathrm{PW} \rightarrow \mathrm{g}}=\sum_{\mathrm{k}} \frac{2+1}{2} \mathrm{Y}_{\mathrm{k}}(\Omega) \sum_{\mathrm{j}} \Delta^{(\mathrm{j})} / 3\left[\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{I}}\right)+4 \mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{A}}^{(\mathrm{j})}\right)+\mathrm{S}_{\mathrm{k}}^{(\mathrm{j})}\left(\mathrm{u}_{\mathrm{F}}^{(\mathrm{j})}\right)\right]\]



Thermal scatter sources from LMR and PW range

If the value of DEMIN is specified to be below the thermal energy
boundary, the portion of the PW range between DEMIN and the thermal
cutoff, as well as the entire LMR, will be contained in the thermal
range. In this case thermal neutrons will downscatter from the thermal
PW range to the LMR, and upscatter from the LMR to the thermal PW range.

The latter thermal source (LMRtoPW) is computed in exactly the same
manner as used to compute the UMRtoPW source S_{HI}, described
in Downscatter source from high region of the UMR to the PW range (SHI). On the other hand, the scatter source from the
thermal PW to the LMR is computed with a similar approach as given in
the previous section for epithermal PWtoLMR scatter. In this case (382)
is used as before, except the source moments are not obtained from the
submoment expansion in (380), but rather by evaluating the PW thermal
scatter expression in (379).

In performing the transport calculation for any group g in the LMR
range, the PWtoMG source component in (382) is added to the MGtoMG
scattering into g from all groups in the UMR and LMR ranges,
respectively, to obtain the total scatter source.

