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

Accurately simulate pair production

Browse source
This commit is contained in:
amandalund 2018-04-26 18:44:39 -05:00
parent e9bd7a203e
commit e12449754e
4 changed files with 230 additions and 36 deletions
Download patch file Download diff file Expand all files Collapse all files

59
openmc/data/photon.py
View file

@ -106,6 +106,22 @@ _STOPPING_POWERS = {}
# for each element are in a 2D array with shape (n, k) stored on the key 'Z'.
_BREMSSTRAHLUNG = {}
# Reduced screening radii for Z = 1-99 from F. Salvat, J. M. Fernández-Varea,
# and J. Sempau, "PENELOPE-2011: A Code System for Monte Carlo Simulation of
# Electron and Photon Transport," OECD-NEA, Issy-les-Moulineaux, France (2011).
_REDUCED_SCREENING_RADIUS = [
122.81, 73.167, 69.228, 67.301, 64.696, 61.228, 57.524, 54.033, 50.787,
47.851, 46.373, 45.401, 44.503, 43.815, 43.074, 42.321, 41.586, 40.953,
40.524, 40.256, 39.756, 39.144, 38.462, 37.778, 37.174, 36.663, 35.986,
35.317, 34.688, 34.197, 33.786, 33.422, 33.068, 32.740, 32.438, 32.143,
31.884, 31.622, 31.438, 31.142, 30.950, 30.758, 30.561, 30.285, 30.097,
29.832, 29.581, 29.411, 29.247, 29.085, 28.930, 28.721, 28.580, 28.442,
28.312, 28.139, 27.973, 27.819, 27.675, 27.496, 27.285, 27.093, 26.911,
26.705, 26.516, 26.304, 26.108, 25.929, 25.730, 25.577, 25.403, 25.245,
25.100, 24.941, 24.790, 24.655, 24.506, 24.391, 24.262, 24.145, 24.039,
23.922, 23.813, 23.712, 23.621, 23.523, 23.430, 23.331, 23.238, 23.139,
23.048, 22.967, 22.833, 22.694, 22.624, 22.545, 22.446, 22.358, 22.264
]
class AtomicRelaxation(EqualityMixin):
"""Atomic relaxation data.
@ -351,19 +367,6 @@ class IncidentPhoton(EqualityMixin):
Number of protons in the target nucleus
atomic_relaxation : openmc.data.AtomicRelaxation or None
Atomic relaxation data
compton_profiles : dict
Dictionary of Compton profile data with keys 'num_electrons' (number of
electrons in each subshell), 'binding_energy' (ionization potential of
each subshell), and 'J' (Hartree-Fock Compton profile as a function of
the projection of the electron momentum on the scattering vector,
:math:`p_z` for each subshell). Note that subshell occupancies may not
match the atomic relaxation data.
stopping_powers : dict
Dictionary of stopping power data with keys 'energy' (in eV), 'density'
(mass density in g/cm:sup:`3`), 'I' (mean excitation energy),
's_collision' (collision stopping power in eV cm:sup:`2`/g),
's_radiative' (radiative stopping power in eV cm:sup:`2`/g), and
'density_effect' (density effect parameter).
bremsstrahlung : dict
Dictionary of bremsstrahlung DCS data with keys 'electron_energy'
(incident electron kinetic energy values in eV), 'photon_energy'
@ -371,9 +374,27 @@ class IncidentPhoton(EqualityMixin):
kinetic energy), and 'dcs' (cross sectin values in mb). 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
Dictionary of Compton profile data with keys 'num_electrons' (number of
electrons in each subshell), 'binding_energy' (ionization potential of
each subshell), and 'J' (Hartree-Fock Compton profile as a function of
the projection of the electron momentum on the scattering vector,
:math:`p_z` for each subshell). Note that subshell occupancies may not
match the atomic relaxation data.
reactions : collections.OrderedDict
Contains the cross sections for each photon reaction. The keys are MT
values and the values are instances of :class:`PhotonReaction`.
reduced_screening_radius : float
Reduced screening radius :math:`R m_e c/\hbar`, where R is the screening
radius for an atom of atomic number Z under the assumption that the
Coulomb field of the nucleus is exponentially screened by atomic electrons.
:math:`\hbar/m_e c` is the Compton wavelength of the electron.
stopping_powers : dict
Dictionary of stopping power data with keys 'energy' (in eV), 'density'
(mass density in g/cm:sup:`3`), 'I' (mean excitation energy),
's_collision' (collision stopping power in eV cm:sup:`2`/g),
's_radiative' (radiative stopping power in eV cm:sup:`2`/g), and
'density_effect' (density effect parameter).
summed_reactions : collections.OrderedDict
Contains summed cross sections. The keys are MT values and the values
are instances of :class:`PhotonReaction`.
@ -418,6 +439,14 @@ class IncidentPhoton(EqualityMixin):
def name(self):
return ATOMIC_SYMBOL[self.atomic_number]
@property
def reduced_screening_radius(self):
if self.atomic_number < 100:
return _REDUCED_SCREENING_RADIUS[self.atomic_number - 1]
else:
raise IndexError('No reduced screening radius for '
'Z={}.'.format(self.atomic_number))
@atomic_number.setter
def atomic_number(self, atomic_number):
cv.check_type('atomic number', atomic_number, Integral)
@ -755,6 +784,10 @@ class IncidentPhoton(EqualityMixin):
shell_group.attrs['designators'] = np.array(designators, dtype='S')
# Write reduced screening radius
if Z < 100:
group.attrs['reduced_screening_radius'] = self.reduced_screening_radius
# Write Compton profiles
if self.compton_profiles:
compton_group = group.create_group('compton_profiles')

22
src/photon_header.F90
View file

@ -53,6 +53,11 @@ module photon_header
! dictionary gives an index in shells(:)
type(ElectronSubshell), allocatable :: shells(:)
! Pair production data
real(8) :: reduced_screening_radius
real(8) :: coulomb_correction
real(8) :: correction_factor_coeffs(4)
! Compton profile data
real(8), allocatable :: profile_pdf(:,:)
real(8), allocatable :: profile_cdf(:,:)
@ -124,6 +129,7 @@ contains
integer :: n_k
integer :: n_e
character(3), allocatable :: designators(:)
real(8) :: a
real(8) :: c
real(8), allocatable :: matrix(:,:)
@ -274,6 +280,22 @@ contains
call read_dataset(this % binding_energy, rgroup, 'binding_energy')
this % electron_pdf(:) = this % electron_pdf / sum(this % electron_pdf)
! Get reduced screening radius
call read_attribute(this % reduced_screening_radius, group_id, &
'reduced_screening_radius')
! Compute the high-energy Coulomb correction
a = this % Z / FINE_STRUCTURE
this % coulomb_correction = a**2*(ONE/(ONE + a**2) + 0.202059_8 &
- 0.03693_8*a**2 + 0.00835_8*a**4 - 0.00201_8*a**6 + 0.00049_8*a**8 &
- 0.00012_8*a**10 + 0.00003_8*a**12)
! Compute the coefficients of the correction factor
this % correction_factor_coeffs(1) = -0.1774_8 - 12.10_8*a + 11.18_8*a**2
this % correction_factor_coeffs(2) = 8.523_8 + 73.26_8*a - 44.41_8*a**2
this % correction_factor_coeffs(3) = -13.52_8 - 121.1_8*a + 96.41_8*a**2
this % correction_factor_coeffs(4) = 8.946_8 + 62.05_8*a - 63.41_8*a**2
! Read Compton profiles
dset_id = open_dataset(rgroup, 'J')
call get_shape(dset_id, dims2)

137
src/photon_physics.F90
View file

@ -378,6 +378,143 @@ contains
end subroutine atomic_relaxation
!===============================================================================
! PAIR_PRODUCTION samples the kinetic energy and direction of the electron and
! positron created when a photon is absorbed near an atomic nucleus. The
! simulation procedure follows the semiempirical model outlined in F. Salvat, J.
! M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: A Code System for Monte
! Carlo Simulation of Electron and Photon Transport," OECD-NEA,
! Issy-les-Moulineaux, France (2011).
!===============================================================================
subroutine pair_production(elm, alpha, E_electron, E_positron, uvw_electron, &
uvw_positron)
type(PhotonInteraction), intent(in) :: elm
real(8), intent(in) :: alpha
real(8), intent(out) :: E_electron
real(8), intent(out) :: E_positron
real(8), intent(out) :: uvw_electron(3)
real(8), intent(out) :: uvw_positron(3)
integer :: i
real(8) :: f
real(8) :: a
real(8) :: b
real(8) :: r
real(8) :: rn
real(8) :: beta
real(8) :: mu
real(8) :: phi
real(8) :: e, e_min, e_max
real(8) :: t1, t2, t3, t4
real(8) :: u1, u2
real(8) :: phi1, phi2
real(8) :: phi1_max, phi2_max
real(8) :: c(4)
! Compute the minimum and maximum values of the electron reduced energy,
! i.e. the fraction of the photon energy that is given to the electron
e_min = ONE/alpha
e_max = ONE - ONE/alpha
! The reduced screening radius r is the ratio of the screening radius to
! the Compton wavelength of the electron, where the screening radius is
! obtained under the assumption that the Coulomb field of the nucleus is
! exponentially screened by atomic electrons. This allows us to use a
! simplified atomic form factor and analytical approximations of the
! screening functions in the pair production DCS instead of computing the
! screening functions numerically.
r = elm % reduced_screening_radius
! The analytical approximation of the DCS underestimates the cross section
! at low energies. The correction factor f compensates for this.
a = sqrt(TWO/alpha)
c = elm % correction_factor_coeffs
f = c(1)*a + c(2)*a**2 + c(3)*a**3 + c(4)*a**4
! Calculate phi_1(1/2) and phi_2(1/2). The unnormalized PDF for the reduced
! energy is given by p = 2*(1/2 - e)^2*phi_1(e) + phi_2(e), where phi_1 and
! phi_2 are non-negative and maximum at e = 1/2.
b = TWO*r/alpha
t1 = TWO*log(ONE + b**2)
t2 = b*atan(ONE/b)
t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2))
t4 = FOUR*log(r) - FOUR*elm % coulomb_correction + f
phi1_max = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4
phi2_max = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4
! To aid sampling, the unnormalized PDF can be expressed as
! p = u_1*U_1(e)*pi_1(e) + u_2*U_2(e)*pi_2(e), where pi_1 and pi_2 are
! normalized PDFs on the interval (e_min, e_max) from which values of e can
! be sampled using the inverse transform method, and
! U_1 = phi_1(e)/phi_1(1/2) and U_2 = phi_2(e)/phi_2(1/2) are valid
! rejection functions. The reduced energy can now be sampled using a
! combination of the composition and rejection methods.
u1 = TWO/THREE*(HALF - ONE/alpha)**2*phi1_max
u2 = phi2_max
do
rn = prn()
! Sample the index i in (1, 2) using the point probabilities
! p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2)
if (prn() < u1/(u1 + u2)) then
i = 1
! Sample e from pi_1 using the inverse transform method
if (rn >= HALF) then
e = HALF + (HALF - ONE/alpha)*(TWO*rn - ONE)**(ONE/THREE)
else
e = HALF - (HALF - ONE/alpha)*(ONE - TWO*rn)**(ONE/THREE)
end if
else
i = 2
! Sample e from pi_2 using the inverse transform method
e = ONE/alpha + (HALF - ONE/alpha)*TWO*rn
end if
! Calculate phi_i(e) and deliver e if rn <= U_i(e)
b = r/(TWO*alpha*e*(ONE - e))
t1 = TWO*log(ONE + b**2)
t2 = b*atan(ONE/b)
t3 = b**2*(FOUR - FOUR*t2 - THREE*log(ONE + ONE/b**2))
if (i == 1) then
phi1 = 7.0_8/THREE - t1 - 6.0_8*t2 - t3 + t4
if (prn() <= phi1/phi1_max) exit
else
phi2 = 11.0_8/6.0_8 - t1 - THREE*t2 + HALF*t3 + t4
if (prn() <= phi2/phi2_max) exit
end if
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
! Sample the direction 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)
rn = TWO*prn() - ONE
mu = (rn + beta)/(rn*beta + ONE)
phi = TWO*PI*prn()
uvw_electron(1) = mu
uvw_electron(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw_electron(3) = sqrt(ONE - mu*mu)*sin(phi)
! Sample the direction of the positron
beta = sqrt(E_positron*(E_positron + TWO*MASS_ELECTRON)) &
/ (E_positron + MASS_ELECTRON)
rn = TWO*prn() - ONE
mu = (rn + beta)/(rn*beta + ONE)
phi = TWO*PI*prn()
uvw_positron(1) = mu
uvw_positron(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw_positron(3) = sqrt(ONE - mu*mu)*sin(phi)
end subroutine pair_production
!===============================================================================
! THICK_TARGET_BREMSSTRAHLUNG
!===============================================================================

48
src/physics.F90
View file

@ -12,7 +12,7 @@ module physics
use particle_header, only: Particle
use photon_header
use photon_physics, only: rayleigh_scatter, compton_scatter, &
atomic_relaxation, &
atomic_relaxation, pair_production, &
thick_target_bremsstrahlung
use physics_common
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
@ -172,10 +172,15 @@ contains
real(8) :: alpha ! photon energy divided by electron rest mass
real(8) :: alpha_out ! outgoing photon energy over electron rest mass
real(8) :: mu ! scattering cosine
real(8) :: mu_electron ! electron scattering cosine
real(8) :: mu_positron ! positron scattering cosine
real(8) :: phi ! azimuthal angle
real(8) :: E_electron ! electron energy
real(8) :: uvw(3) ! new direction
real(8) :: rel_vel ! relative velocity of electron
real(8) :: E_electron ! electron energy
real(8) :: E_positron ! positron energy
real(8) :: uvw_electron(3) ! new electron direction
real(8) :: uvw_positron(3) ! new positron direction
! Kill photon if below energy cutoff -- an extra check is made here because
! photons with energy below the cutoff may have been produced by neutrons
@ -273,29 +278,26 @@ contains
end do
end if
prob = prob_after
! Pair production
prob = prob + micro_photon_xs(i_element) % pair_production
if (prob > cutoff) then
call pair_production(elm, alpha, E_electron, E_positron, uvw_electron, &
uvw_positron)
! Create secondary electron
call p % create_secondary(uvw_electron, E_electron, ELECTRON, .true.)
! Create secondary positron
call p % create_secondary(uvw_positron, E_positron, POSITRON, .true.)
p % event_MT = PAIR_PROD
p % alive = .false.
p % E = ZERO
end if
end associate
! Pair production
prob = prob + micro_photon_xs(i_element) % pair_production
if (prob > cutoff) then
! Sample angle isotropically
mu = TWO*prn() - ONE
phi = TWO*PI*prn()
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu)*cos(phi)
uvw(3) = sqrt(ONE - mu*mu)*sin(phi)
! Compute the kinetic energy of each particle
E_electron = HALF * (p % E - 2 * MASS_ELECTRON)
! Create electron-positron pair traveling in opposite directions
call p % create_secondary( uvw, E_electron, ELECTRON, .true.)
call p % create_secondary(-uvw, E_electron, POSITRON, .true.)
p % event_MT = PAIR_PROD
p % alive = .false.
p % E = ZERO
end if
end subroutine sample_photon_reaction
!===============================================================================