From 7536de564480e592678de7173ea1496353acbcb6 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 10 Apr 2019 13:46:44 -0400 Subject: [PATCH] Account for fluorescent photons release in heating cross section computing --- openmc/data/photon.py | 76 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 61 insertions(+), 15 deletions(-) diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 974b9323cf..aa01eac9b4 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -161,6 +161,7 @@ class AtomicRelaxation(EqualityMixin): self.binding_energy = binding_energy self.num_electrons = num_electrons self.transitions = transitions + self._e_fluorescence = {} @property def binding_energy(self): @@ -390,6 +391,41 @@ class AtomicRelaxation(EqualityMixin): _SUBSHELLS, range(len(_SUBSHELLS))) group.create_dataset('transitions', data=df.values.astype(float)) + def energy_fluorescence(self, shell): + """Compute expected energy of fluorescent photons for the shell + """ + + if shell not in self.binding_energy: + raise KeyError('Invalid shell {}.'.format(shell)) + + if shell in self._e_fluorescence: + # Already computed + return self._e_fluorescence[shell] + else: + e = 0.0 + if shell not in self.transitions or self.transitions[shell].empty: + e = self.binding_energy[shell] + else: + df = self.transitions[shell] + for index, row in df.iterrows(): + e_row = 0.0 + primary = row['secondary'] + secondary = row['tertiary'] + if secondary is None: + # Fluorescent photon release in radiative transition + e_row += row['energy (eV)'] + else: + # Fill pole left by auger eletron + e_row += self.energy_fluorescence(secondary) + + # Fill photonelectron pole + e_row += self.energy_fluorescence(primary) + + # Expected fluorescent photon energy + e += e_row * row['probability'] + + self._e_fluorescence[shell] = e + return e class IncidentPhoton(EqualityMixin): r"""Photon interaction data. @@ -965,20 +1001,20 @@ class IncidentPhoton(EqualityMixin): """Compute heating cross sections (KERMA) Photon energy is deposited as energy loss in three reactions: - incoherent scattering, photoelectric effect and pair production. + incoherent scattering, pair production and photoelectric effect. The point-wise heating cross section is calculated as: .. math:: - \sigma_{Hx} &= (E - \overline{E}_x(E)) \times \sigma_x(E), x \in \left \{ I, PE, PP \right\} + \sigma_{Hx} &= (E - \overline{E}_x(E)) \times \sigma_x(E), x \in \left \{ I, PP, PE \right\} \overline{E}_I (E) &= \frac {\int E' \sigma_I (E,E',\mu) d\mu} {\int \sigma_I (E,E',\mu) d\mu} - \overline{E}_{PE} &= 0 - \overline{E}_{PP} &= 2 m_e c^2 = 1.022 \times 10^6 eV - The differential cross section for incoherent scattering can be - found in the theory manual. + \overline{E}_{PE} &= E_{fluorescent photons} + + The differential cross section representation for incoherent + scattering can be found in the theory manual. """ @@ -987,41 +1023,51 @@ class IncidentPhoton(EqualityMixin): for mt in (504, 515, 517, 522): if mt not in self: warn('Cross sections for MT={} do not exist. The total heating ' - 'cross sections may be wrong'.format(mt)) + 'to be calculated is probably incorrect'.format(mt)) else: energy = np.union1d(energy, self[mt].xs.x) heating_xs = np.zeros_like(energy) + # Incoherent scattering if 504 in self: rx = self[504] - def dsigma_dmu(mu, E): + def _dsigma_dmu(mu, E): k = E / MASS_ELECTRON_EV kout = k / (1.0 + k * (1.0 - mu)) x = E * np.sqrt(0.5 * (1.0 - mu)) / PLANCK_C return np.pi * R0**2 * (kout/k)**2 * (kout/k + k/kout + mu*mu - 1.0) * rx.scattering_factor(x) - def dsigma_dmu_times_eout(mu, E): + def _eout_dsigma_dmu(mu, E): Eout = E / (1.0 + E / MASS_ELECTRON_EV * (1 - mu)) - return dsigma_dmu(mu, E) * Eout - def average_eout(E): - integral_sigma = quad(dsigma_dmu, -1.0, 1.0, + return Eout * _dsigma_dmu(mu, E) + def _eout_average(E): + #TODO: the integrals are currently time consuming + integral_sigma = quad(_dsigma_dmu, -1.0, 1.0, args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] - integral_sigma_e = quad(dsigma_dmu_times_eout, -1.0, 1.0, + integral_sigma_e = quad(_eout_dsigma_dmu, -1.0, 1.0, args=(E,), epsabs=0.0, epsrel=1.0e-3)[0] return integral_sigma_e / integral_sigma - e_out = np.vectorize(lambda x: average_eout(x))(energy) + e_out = np.vectorize(lambda x: _eout_average(x))(energy) heating_xs += (energy - e_out) * rx.xs(energy) + # Pair production, electron field if 515 in self: heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[515].xs(energy) + # Pair production, nuclear field if 517 in self: heating_xs += (energy - 2*MASS_ELECTRON_EV)*self[517].xs(energy) + # Photoelectric effect if 522 in self: - heating_xs += (energy - 0.)*self[522].xs(energy) + # Account for fluorescent photons + for mt, rx in self.reactions.items(): + if mt >= 534 and mt <= 572: + shell = _REACTION_NAME[mt][1] + e_f = self.atomic_relaxation.energy_fluorescence(shell) + heating_xs += (energy - e_f) * rx.xs(energy) heat_rx = PhotonReaction(525) heat_rx.xs = Tabulated1D(energy, heating_xs, [energy.size], [5])