crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Address #1023 comments

Browse source
This commit is contained in:
amandalund 2018-07-09 21:33:48 -05:00
parent 899c1107db
commit d29923f6bf
9 changed files with 69 additions and 61 deletions
Download patch file Download diff file Expand all files Collapse all files

4
docs/source/io_formats/nuclear_data.rst
View file

@ -181,7 +181,7 @@ Incident Photon Data
:Datasets: - **I** (*double*) -- Mean excitation energy in [eV]
- **energy** (*double[]*) -- Energies in [eV]
- **s_collision** (*double[]*) -- Collisiong stopping power in [eV-cm\ :sup:`2`\ /g]
- **s_collision** (*double[]*) -- Collision stopping power in [eV-cm\ :sup:`2`\ /g]
- **s_radiative** (*double[]*) -- Radiative stopping power in [eV-cm\ :sup:`2`\ /g]
**/<element>/subshells/**
@ -199,7 +199,7 @@ Incident Photon Data
:Attributes:
- **threshold_idx** (*int*) -- Index on the energy
grid that the reaction threshold
grid of the reaction threshold
-------------------------------
Thermal Neutron Scattering Data

62
docs/source/methods/photon_physics.rst
View file

@ -5,7 +5,7 @@ Photon Physics
==============
Photons, being neutral particles, behave much in the same manner as neutrons,
traveling in straight lines and experiencing occassional collisions which change
traveling in straight lines and experiencing occasional collisions which change
their energy and direction. Photons undergo four basic interactions as they pass
through matter: coherent (Rayleigh) scattering, incoherent (Compton) scattering,
photoelectric effect, and pair/triplet production. Photons with energy in the
@ -31,11 +31,11 @@ cross section is
.. math::
:label: thomson
\frac{d\sigma}{d\mu} = \pi r_0^2 ( 1 + \mu^2 )
\frac{d\sigma}{d\mu} = \pi r_e^2 ( 1 + \mu^2 )
where :math:`\mu` is the cosine of the scattering angle and :math:`r_0` is the
classical radius of the electron. Thomson scattering can generally occur when
the photon energy is much less than rest mass energy of the particle.
where :math:`\mu` is the cosine of the scattering angle and :math:`r_e` is the
classical electron radius. Thomson scattering can generally occur when the
photon energy is much less than the rest mass energy of the particle.
In practice, most elastic scattering of photons off electrons happens not with
free electrons but those bound in atoms. This process is known as Rayleigh
@ -50,30 +50,33 @@ The differential cross section for Rayleigh scattering is given by
.. math::
:label: coherent-xs
\frac{d\sigma(E,E',\mu)}{d\mu} = \pi r_0^2 ( 1 + \mu^2 ) \left [ ( F(x, Z) +
F'(E) )^2 + F''(E)^2 \right ]
\frac{d\sigma(E,E',\mu)}{d\mu} &= \pi r_e^2 ( 1 + \mu^2 )~\left| F(x,Z)
+ F' + iF'' \right|^2 \\
&= \pi r_e^2 ( 1 + \mu^2 ) \left [ ( F(x,Z)
+ F'(E) )^2 + F''(E)^2 \right ]
where :math:`F(x,Z)` is a form factor as a function of the momentum transfer
:math:`x` and the atomic number :math:`Z` and :math:`F' + iF''` is a factor that
accounts for `anomalous scattering`_ which can occur near absorption edges. In a
Monte Carlo simulation, when coherent scattering occurs, we only need to sample
the scattering angle using the differential cross section in :eq:`coherent-xs`
since the energy of the photon does not change. In OpenMC, anomalous scattering
is ignored such that differential cross section becomes
:math:`x` and the atomic number :math:`Z` and the term :math:`F' + iF''`
accounts for `anomalous scattering`_ which can occur near absorption edges. In
a Monte Carlo simulation, when coherent scattering occurs, we only need to
sample the scattering angle using the differential cross section in
:eq:`coherent-xs` since the energy of the photon does not change. In OpenMC,
anomalous scattering is ignored such that the differential cross section
becomes
.. math::
:label: coherent-xs-openmc
\frac{d\sigma(E,E',\mu)}{d\mu} = \pi r_0^2 ( 1 + \mu^2 ) F(x, Z)^2
\frac{d\sigma(E,E',\mu)}{d\mu} = \pi r_e^2 ( 1 + \mu^2 ) F(x, Z)^2
To construct a proper probability density, we need to normalize the differential
cross section in :eq:`coherent-xs-openmc` by the integrated coherent scattering
cross section:
To construct a proper probability density, we need to normalize the
differential cross section in :eq:`coherent-xs-openmc` by the integrated
coherent scattering cross section:
.. math::
:label: coherent-pdf-1
p(\mu) d\mu = \frac{\pi r_0^2}{\sigma(E)} ( 1 + \mu^2 ) F(x, Z)^2 d\mu.
p(\mu) d\mu = \frac{\pi r_e^2}{\sigma(E)} ( 1 + \mu^2 ) F(x, Z)^2 d\mu.
Since the form factor is given in terms of the momentum transfer, it is more
convenient to change variables of the probability density to :math:`x^2`. The
@ -84,23 +87,24 @@ momentum transfer is traditionally expressed as
x = \kappa \alpha \sqrt{1 - \mu}
where the coefficient :math:`\kappa` can be shown to be
where :math:`\alpha` is the ratio of the photon energy to the electron rest
mass, and the coefficient :math:`\kappa` can be shown to be
.. math::
:label: kappa
\kappa = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329,
:math:`m_e` is the mass of the electron, :math:`c` is the speed of light
where :math:`m_e` is the mass of the electron, :math:`c` is the speed of light
in a vacuum, and :math:`h` is Planck's constant. Using :eq:`momentum-transfer`,
we have that :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 =
we have :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 =
-1/(\kappa\alpha)^2`. The probability density in :math:`x^2` is
.. math::
:label: coherent-pdf-x2
p(x^2) dx^2 = p(\mu) \left | \frac{d\mu}{dx^2} \right | dx^2 = \frac{2\pi
r_0^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left (
r_e^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left (
\frac{1 + \mu^2}{2} \right ) \left ( \frac{F(x, Z)^2}{A(\bar{x}^2, Z)} \right ) dx^2
where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for
@ -116,7 +120,7 @@ and :math:`A(x^2, Z)` is the integral of the square of the form factor:
.. math::
:label: coherent-int-ff
A(x^2, Z) = \int_0^{x^2} F(\chi, Z)^2 d\chi^2.
A(x^2, Z) = \int_0^{x^2} F(x,Z)^2 dx^2.
As you see, we have multiplied and divided the probability density by the
integral of the squared form factor so that the density in :eq:`coherent-pdf-x2`
@ -127,7 +131,7 @@ run-time to do a table search on the cumulative distribution function:
.. math::
:label: coherent-form-factor-cdf
\frac{\int_0^{x^2} F(\chi,Z)^2 d\chi^2}{\int_0^{\bar{x}^2} F(x,Z)^2 dx^2}
\frac{\int_0^{x^2} F(x,Z)^2 dx^2}{\int_0^{\bar{x}^2} F(x,Z)^2 dx^2}
Once a trial :math:`x^2` value has been selected, we can calculate :math:`\mu`
and perform rejection sampling using the Thomson scattering differential cross
@ -162,7 +166,7 @@ the two authors who discovered it:
.. math::
:label: klein-nishina
\frac{d\sigma_{KN}}{d\mu} = \pi r_0^2 \left ( \frac{\alpha'}{\alpha} \right
\frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right
) \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1 \right
]
@ -188,7 +192,7 @@ differential cross section for incoherent scattering is given by
.. math::
:label: incoherent-xs
\frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_0^2 \left (
\frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_e^2 \left (
\frac{\alpha'}{\alpha} \right )^2 \left [ \frac{\alpha'}{\alpha} +
\frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] S(x,Z)
@ -213,6 +217,10 @@ Doppler Energy Broadening
LA-UR-04-0487_ and LA-UR-04-0488_
Compton Electrons
+++++++++++++++++
Photoelectric Effect
--------------------
@ -253,7 +261,7 @@ Thick-Target Bremsstrahlung Approximation
+++++++++++++++++++++++++++++++++++++++++
.. _Koblinger: http://www.tandfonline.com/doi/abs/10.13182/NSE75-A26646
.. _Koblinger: https://doi.org/10.13182/NSE75-A26663
.. _anomalous scattering: http://pd.chem.ucl.ac.uk/pdnn/diff1/anomscat.htm

11
openmc/data/photon.py
View file

@ -110,10 +110,11 @@ class AtomicRelaxation(EqualityMixin):
"""Atomic relaxation data.
This class stores the binding energy, number of electrons, and electron
transitions possible from ioniziation for each subshell with an atom. All of
the data originates from an ENDF-6 atomic relaxation sub-library
(NSUB=6). Instances of this class are not normally instantiated directly but
rather created using the factory method :math:`AtomicRelaxation.from_endf`.
transitions possible from ioniziation for each electron subshell of an
atom. All of the data originates from an ENDF-6 atomic relaxation
sub-library (NSUB=6). Instances of this class are not normally instantiated
directly but rather created using the factory method
:math:`AtomicRelaxation.from_endf`.
Parameters
----------
@ -353,7 +354,7 @@ class IncidentPhoton(EqualityMixin):
Dictionary of bremsstrahlung DCS data with keys 'electron_energy'
(incident electron kinetic energy values in [eV]), 'photon_energy'
(ratio of the energy of the emitted photon to the incident electron
kinetic energy), and 'dcs' (cross sectin values in [b]). The cross
kinetic energy), and 'dcs' (cross section values in [b]). The cross
sections are in scaled form: :math:`(\beta^2/Z^2) E_k (d\sigma/dE_k)`,
where :math:`E_k` is the energy of the emitted photon.
compton_profiles : dict

2
src/constants.F90
View file

@ -78,7 +78,7 @@ module constants
MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu
MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2
MASS_PROTON = 1.007276466879_8, & ! mass of a proton in amu
MASS_ELECTRON = 0.5109989461e6_8, & ! electron mass energy equivalent in eV
MASS_ELECTRON_EV = 0.5109989461e6_8, & ! electron mass energy equivalent in eV/c^2
FINE_STRUCTURE = 137.035999139_8, & ! inverse fine structure constant
PLANCK_C = 1.2398419739062977e4_8,& ! Planck's constant times c in eV-Angstroms
AMU = 1.660539040e-27_8, & ! 1 amu in kg

2
src/input_xml.F90
View file

@ -2610,7 +2610,7 @@ contains
do l = 1, filt % n_bins
if (filt % particles(l) /= NEUTRON) then
call warning("Particle filter other than NEUTRON used with &
&photon transport turn off. All tallies for particle &
&photon transport turned off. All tallies for particle &
&type " // trim(to_str(filt % particles(l))) // " will have no scores")
end if
end do

4
src/material_header.F90
View file

@ -809,7 +809,7 @@ contains
! Issy-les-Moulineaux, France (2011).
if (positron_) then
do i = 1, n_e
t = log(ONE + 1.0e6_8*ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON))
t = log(ONE + 1.0e6_8*ttb_e_grid(i)/(Z_eq_sq*MASS_ELECTRON_EV))
r = ONE - exp(-1.2359e-1_8*t + 6.1274e-2_8*t**2 - 3.1516e-2_8*t**3 + &
7.7446e-3_8*t**4 - 1.0595e-3_8*t**5 + 7.0568e-5_8*t**6 - &
1.808e-6_8*t**7)
@ -846,7 +846,7 @@ contains
x = x_l + (k - k_l) * (x_r - x_l) / (k_r - k_l)
! Ratio of the velocity of the charged particle to the speed of light
beta = sqrt(e*(e + TWO*MASS_ELECTRON)) / (e + MASS_ELECTRON)
beta = sqrt(e*(e + TWO*MASS_ELECTRON_EV)) / (e + MASS_ELECTRON_EV)
! Compute the integrand of the PDF
f(j) = x / (beta**2 * stopping_power(j) * w)

28
src/photon_physics.F90
View file

@ -107,20 +107,20 @@ contains
! Note that the parameter used here does not correspond exactly to the
! momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the
! parameter as defined by Hubbell, where the actual data comes from
x = MASS_ELECTRON/PLANCK_C*alpha*sqrt(HALF*(ONE - mu))
x = MASS_ELECTRON_EV/PLANCK_C*alpha*sqrt(HALF*(ONE - mu))
! Calculate S(x, Z) and S(x_max, Z)
form_factor_x = el % incoherent_form_factor % evaluate(x)
if (form_factor_xmax == ZERO) then
form_factor_xmax = el % incoherent_form_factor % evaluate(&
MASS_ELECTRON/PLANCK_C*alpha)
MASS_ELECTRON_EV/PLANCK_C*alpha)
end if
! Perform rejection on form factor
if (prn() < form_factor_x / form_factor_xmax) then
if (use_doppler_) then
call compton_doppler(el, alpha, mu, e_out, i_shell)
alpha_out = e_out/MASS_ELECTRON
alpha_out = e_out/MASS_ELECTRON_EV
else
i_shell = 0
end if
@ -168,16 +168,16 @@ contains
e_b = el % binding_energy(i_shell)
! Determine p_z,max
e = alpha*MASS_ELECTRON
e = alpha*MASS_ELECTRON_EV
if (e < e_b) then
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV
exit
end if
pz_max = -FINE_STRUCTURE*(e_b - (e - e_b)*alpha*(ONE - mu)) / &
sqrt(TWO*e*(e - e_b)*(ONE - mu) + e_b**2)
if (pz_max < ZERO) then
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV
exit
end if
@ -230,7 +230,7 @@ contains
quad = b**2 - FOUR*a*c
if (quad < 0) then
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON
e_out = alpha/(1 + alpha*(1 - mu))*MASS_ELECTRON_EV
exit
end if
quad = sqrt(quad)
@ -280,7 +280,7 @@ contains
do
! Determine maximum value of x^2
x2_max = (MASS_ELECTRON/PLANCK_C*alpha)**2
x2_max = (MASS_ELECTRON_EV/PLANCK_C*alpha)**2
! Determine F(x^2_max, Z)
F_max = el % coherent_int_form_factor % evaluate(x2_max)
@ -517,20 +517,20 @@ contains
end do
! Compute the kinetic energy of the electron and the positron
E_electron = (alpha*e - ONE)*MASS_ELECTRON
E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON
E_electron = (alpha*e - ONE)*MASS_ELECTRON_EV
E_positron = (alpha*(ONE - e) - ONE)*MASS_ELECTRON_EV
! Sample the scattering angle of the electron. The cosine of the polar
! angle of the direction relative to the incident photon is sampled from
! p(mu) = C/(1 - beta*mu)^2 using the inverse transform method.
beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON)) &
/ (E_electron + MASS_ELECTRON)
beta = sqrt(E_electron*(E_electron + TWO*MASS_ELECTRON_EV)) &
/ (E_electron + MASS_ELECTRON_EV)
rn = TWO*prn() - ONE
mu_electron = (rn + beta)/(rn*beta + ONE)
! Sample the scattering angle of the positron
beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON)) &
/ (E_positron + MASS_ELECTRON)
beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON_EV)) &
/ (E_positron + MASS_ELECTRON_EV)
rn = TWO*prn() - ONE
mu_positron = (rn + beta)/(rn*beta + ONE)

16
src/physics.F90
View file

@ -195,7 +195,7 @@ contains
p % event_nuclide = i_element
! Calculate photon energy over electron rest mass equivalent
alpha = p % E/MASS_ELECTRON
alpha = p % E/MASS_ELECTRON_EV
! For tallying purposes, this routine might be called directly. In that
! case, we need to sample a reaction via the cutoff variable
@ -226,7 +226,7 @@ contains
end if
! Create Compton electron
E_electron = (alpha - alpha_out)*MASS_ELECTRON - e_b
E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b
mu_electron = (alpha - alpha_out*mu) &
/ sqrt(alpha**2 + alpha_out**2 - TWO*alpha*alpha_out*mu)
phi = TWO*PI*prn()
@ -241,7 +241,7 @@ contains
end if
phi = phi + PI
p % E = alpha_out*MASS_ELECTRON
p % E = alpha_out*MASS_ELECTRON_EV
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu, phi)
p % event_MT = INCOHERENT
return
@ -274,8 +274,8 @@ contains
SAMPLE_MU: do
r = prn()
if (FOUR * (ONE - r) * r >= prn()) then
rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON))&
/ (E_electron + MASS_ELECTRON)
rel_vel = sqrt(E_electron * (E_electron + TWO * MASS_ELECTRON_EV))&
/ (E_electron + MASS_ELECTRON_EV)
mu = (TWO * r + rel_vel - ONE) / &
(TWO * rel_vel * r - rel_vel + ONE)
exit SAMPLE_MU
@ -354,7 +354,7 @@ contains
! energy locally (electron_treatment = ELECTRON_LED) or creates secondary
! bremsstrahlung photons from electron deflections with charged particles
! (electron_treatment = ELECTRON_TTB). Two annihilation photons of energy
! MASS_ELECTRON (0.511 MeV) are created and travel in opposite directions.
! MASS_ELECTRON_EV (0.511 MeV) are created and travel in opposite directions.
!===============================================================================
subroutine sample_positron_reaction(p)
@ -380,8 +380,8 @@ contains
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
! Create annihilation photon pair traveling in opposite directions
call p % create_secondary( uvw, MASS_ELECTRON, PHOTON, .true.)
call p % create_secondary(-uvw, MASS_ELECTRON, PHOTON, .true.)
call p % create_secondary( uvw, MASS_ELECTRON_EV, PHOTON, .true.)
call p % create_secondary(-uvw, MASS_ELECTRON_EV, PHOTON, .true.)
p % E = ZERO
p % alive = .false.

1
src/tallies/tally.F90
View file

@ -298,7 +298,6 @@ contains
if (t % find_filter(FILTER_ENERGYOUT) > 0) then
! Normally, we only need to make contributions to one scoring
! bin. However, in the case of fission, since multiple fission
! neutrons were emitted with different energies, multiple
! outgoing energy bins may have been scored to. The following
! logic treats this special case and results to multiple bins