diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 8644f8e7c..c6202117c 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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') diff --git a/src/photon_header.F90 b/src/photon_header.F90 index 15203e4ff..162b4d526 100644 --- a/src/photon_header.F90 +++ b/src/photon_header.F90 @@ -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) diff --git a/src/photon_physics.F90 b/src/photon_physics.F90 index aea6a6588..4b81e2598 100644 --- a/src/photon_physics.F90 +++ b/src/photon_physics.F90 @@ -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 !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index 982584844..204158f05 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 !===============================================================================