diff --git a/docs/source/methods/energy_deposition.rst b/docs/source/methods/energy_deposition.rst index 4675aee3fb..0babd71246 100644 --- a/docs/source/methods/energy_deposition.rst +++ b/docs/source/methods/energy_deposition.rst @@ -25,19 +25,36 @@ KERMA (Kinetic Energy Release in Materials) [Mack97]_ coefficients for reaction :math:`\times` cross-section (e.g., eV-barn) and can be used much like a reaction cross section for the purpose of tallying energy deposition. -KERMA coefficients can be computed using the energy-balance method with -a nuclear data processing code like NJOY, which performs the following -iteration over all reactions :math:`r` for all isotopes :math:`i` -requested +KERMA coefficients can be computed using the energy-balance method with a +nuclear data processing code like NJOY, which estimates the KERMA coefficients +using the following equation: .. math:: - k_{i, r}(E) = \left(E + Q_{i, r} - \bar{E}_{i, r, n} + k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, x} + \right)\sigma_{i, r}(E), + +where the summation is over each secondary particle type :math:`x`. This +equation states that the energy deposited is equal to the energy of the incident +particle plus the reaction :math:`Q` value less the energy of secondary +particles that are transported away from the reaction site. For neutron +interactions, the energy-balance KERMA coefficient is + +.. math:: + + k_{i, r}(E) = \left(E + Q_{i, r} - \sum\limits_x \bar{E}_{i, r, n} - \bar{E}_{i, r, \gamma}\right)\sigma_{i, r}(E), -removing the energy of neutral particles (neutrons and photons) that are -transported away from the reaction site :math:`\bar{E}`, and the reaction -:math:`Q` value. +where :math:`\bar{E}_{i, r, n}` is the average energy of secondary neutrons and +:math:`\bar{E}_{i, r, \gamma}` is the average energy of secondary photons. For +photon and charged particle interactions, the :math:`Q` value is zero and thus +the KERMA coefficient is + +.. math:: + :label: energy-balance-photon + + k_{i, r}(E) = \left(E - \sum\limits_x \bar{E}_{i, r, x} + \right)\sigma_{i, r}(E). ------- Fission @@ -120,7 +137,7 @@ run with :math:`N918` reflecting fission heating computed from NJOY. This modified heating data is stored as the MT=901 reaction and will be scored if ``heating-local`` is included in :attr:`openmc.Tally.scores`. -Coupled neutron-photon transport +Coupled Neutron-Photon Transport -------------------------------- Here, OpenMC instructs ``heatr`` to assume that energy from photons is not @@ -138,6 +155,50 @@ Let :math:`N301` represent the total heating number returned from this This modified heating data is stored as the MT=301 reaction and will be scored if ``heating`` is included in :attr:`openmc.Tally.scores`. +Photons and Charged Particles +----------------------------- + +In OpenMC, energy deposition from photons or charged particles is scored using +the energy balance method based on Equation :eq:`energy-balance-photon`. Special +consideration is given to electrons and positrons as described below. + ++++++++++++++++++ +Charged Particles ++++++++++++++++++ + +OpenMC tracks photons interaction by interaction so the energy deposited in each +collision is easily attributed back to the nuclide and reaction for which the +photon interacted with. Charged particles (electrons and photons) aren't tracked +in the same way. For charged particles, OpenMC assumes that all their energy +(less the energy of bremsstrahlung radiation) is deposited in the material in +which they were born. In this way it is harder to trace how much energy should +be attributed in each nuclide. + +According to the CSDA approximation (see :ref:`ttb`) the energy deposited by a +charged particle with kinetic energy :math:`T` in the :math:`i`-th element can +be calculated as: + +.. math:: + + E_{i} = \int_{0}^{R(T)} w_{i}S_{\text{col,i}} dx + +where :math:`R(T)` is the CSDA range of the charged particle, +:math:`S_{\text{col},i}` is the collision stopping power of the charged particle +in the :math:`i`-th element and :math:`w_i` is the mass fraction of the +:math:`i`-th element. According to the Bethe formula the collision stopping +power of the :math:`i`-th element is proportional to :math:`Z_i/A_i`, so the +fractional collision stopping power from the :math:`i`-th element is: + +.. math:: + + \frac{w_{i}S_{\text{col},i}(T)}{S_{\text{col}}(T)} = + \frac{\frac{w_{i}Z_{i}}{A_{i}}}{\sum_{i}\frac{w_{i}Z_{i}}{A_{i}}} = + \frac{\gamma_i Z_{i}}{\sum_{i}\gamma_i Z_{i}}. + +where :math:`\gamma_i` is the atomic fraction of the :math:`i`-th element. +Therefore, the energy deposited by charged particles should be attributed to +a given element according to its fractional charge density. + ---------- References ---------- diff --git a/include/openmc/material.h b/include/openmc/material.h index 31fab2ae27..fe587a86f9 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -107,6 +107,10 @@ public: //! \return Density in [g/cm^3] double density_gpcc() const { return density_gpcc_; } + //! Get charge density in [e/b-cm] + //! \return Charge density in [e/b-cm] + double charge_density() const { return charge_density_; }; + //! Get name //! \return Material name const std::string& name() const { return name_; } @@ -177,6 +181,7 @@ public: xt::xtensor atom_density_; //!< Nuclide atom density in [atom/b-cm] double density_; //!< Total atom density in [atom/b-cm] double density_gpcc_; //!< Total atom density in [g/cm^3] + double charge_density_; //!< Total charge density in [e/b-cm] double volume_ {-1.0}; //!< Volume in [cm^3] vector p0_; //!< Indicate which nuclides are to be treated with //!< iso-in-lab scattering diff --git a/src/material.cpp b/src/material.cpp index b8846d7cbf..32384ddf12 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -454,12 +454,15 @@ void Material::normalize_density() // Calculate nuclide atom densities atom_density_ *= density_; - // Calculate density in g/cm^3. + // Calculate density in [g/cm^3] and charge density in [e/b-cm] density_gpcc_ = 0.0; + charge_density_ = 0.0; for (int i = 0; i < nuclide_.size(); ++i) { int i_nuc = nuclide_[i]; double awr = settings::run_CE ? data::nuclides[i_nuc]->awr_ : 1.0; + int z = settings::run_CE ? data::nuclides[i_nuc]->Z_ : 0.0; density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO; + charge_density_ += atom_density_(i) * z; } } @@ -982,12 +985,15 @@ void Material::set_density(double density, const std::string& units) // Recalculate nuclide atom densities based on given density atom_density_ *= density; - // Calculate density in g/cm^3. + // Calculate density in g/cm^3 and charge density in [e/b-cm] density_gpcc_ = 0.0; + charge_density_ = 0.0; for (int i = 0; i < nuclide_.size(); ++i) { int i_nuc = nuclide_[i]; double awr = data::nuclides[i_nuc]->awr_; + int z = settings::run_CE ? data::nuclides[i_nuc]->Z_ : 0.0; density_gpcc_ += atom_density_(i) * awr * MASS_NEUTRON / N_AVOGADRO; + charge_density_ += atom_density_(i) * z; } } else if (units == "g/cm3" || units == "g/cc") { // Determine factor by which to change densities @@ -998,6 +1004,7 @@ void Material::set_density(double density, const std::string& units) density_gpcc_ = density; density_ *= f; atom_density_ *= f; + charge_density_ *= f; } else { throw std::invalid_argument { "Invalid units '" + std::string(units.data()) + "' specified."}; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 5c6386fb42..62b0029074 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -325,6 +325,39 @@ double score_neutron_heating(const Particle& p, const Tally& tally, double flux, return score; } +//! Helper function to obtain particle heating [eV] + +double score_particle_heating(const Particle& p, const Tally& tally, + double flux, int rxn_bin, int i_nuclide, double atom_density) +{ + if (p.type() == ParticleType::neutron) + return score_neutron_heating( + p, tally, flux, rxn_bin, i_nuclide, atom_density); + if (i_nuclide == -1 || i_nuclide == p.event_nuclide() || + p.event_nuclide() == -1) { + // Get the pre-collision energy of the particle. + auto E = p.E_last(); + // The energy deposited is the difference between the pre-collision + // and post-collision energy... + double score = E - p.E(); + // ...less the energy of any secondary particles since they will be + // transported individually later + score -= p.bank_second_E(); + score *= p.wgt_last(); + + // if no event_nuclide (charged particle) scale energy deposition by + // fractional charge density + if (i_nuclide != -1 && p.event_nuclide() == -1) { + const auto& mat {model::materials[p.material()]}; + int z = data::nuclides[i_nuclide]->Z_; + auto i = mat->mat_nuclide_index_[i_nuclide]; + score *= (z * mat->atom_density_[i] / mat->charge_density()); + } + return score; + } + return 0.0; +} + //! Helper function for nu-fission tallies with energyout filters. // //! In this case, we may need to score to multiple bins if there were multiple @@ -1007,23 +1040,8 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, break; case HEATING: - if (p.type() == Type::neutron) { - score = score_neutron_heating( - p, tally, flux, HEATING, i_nuclide, atom_density); - } else { - if (i_nuclide == -1 || i_nuclide == p.event_nuclide()) { - // The energy deposited is the difference between the pre-collision - // and post-collision energy... - score = E - p.E(); - // ...less the energy of any secondary particles since they will be - // transported individually later - score -= p.bank_second_E(); - - score *= p.wgt_last(); - } else { - score = 0.0; - } - } + score = score_particle_heating( + p, tally, flux, HEATING, i_nuclide, atom_density); break; default: @@ -1539,19 +1557,8 @@ void score_general_ce_analog(Particle& p, int i_tally, int start_index, break; case HEATING: - if (p.type() == Type::neutron) { - score = score_neutron_heating( - p, tally, flux, HEATING, i_nuclide, atom_density); - } else { - // The energy deposited is the difference between the pre-collision and - // post-collision energy... - score = E - p.E(); - // ...less the energy of any secondary particles since they will be - // transported individually later - score -= p.bank_second_E(); - - score *= p.wgt_last(); - } + score = score_particle_heating( + p, tally, flux, HEATING, i_nuclide, atom_density); break; default: diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 874e32b85f..c16524689e 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/tests/regression_tests/photon_production/results_true.dat @@ -67,8 +67,8 @@ tally 3: 0.000000E+00 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 +1.774484E+05 +3.148794E+10 0.000000E+00 0.000000E+00 0.000000E+00 @@ -79,8 +79,8 @@ tally 3: 0.000000E+00 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 +7.692488E+03 +5.917437E+07 0.000000E+00 0.000000E+00 0.000000E+00 diff --git a/tests/unit_tests/test_data_multipole.py b/tests/unit_tests/test_data_multipole.py index 4c2ba96d62..bc7136a97c 100644 --- a/tests/unit_tests/test_data_multipole.py +++ b/tests/unit_tests/test_data_multipole.py @@ -52,7 +52,7 @@ def test_from_endf(): pytest.importorskip('vectfit') endf_data = os.environ['OPENMC_ENDF_DATA'] endf_file = os.path.join(endf_data, 'neutrons', 'n-001_H_001.endf') - return openmc.data.WindowedMultipole.from_endf( + assert openmc.data.WindowedMultipole.from_endf( endf_file, log=True, wmp_options={"n_win": 400, "n_cf": 3}) @@ -60,5 +60,5 @@ def test_from_endf_search(): pytest.importorskip('vectfit') endf_data = os.environ['OPENMC_ENDF_DATA'] endf_file = os.path.join(endf_data, 'neutrons', 'n-095_Am_244.endf') - return openmc.data.WindowedMultipole.from_endf( + assert openmc.data.WindowedMultipole.from_endf( endf_file, log=True, wmp_options={"search": True, 'rtol':1e-2}) diff --git a/tests/unit_tests/test_nuclide_heating.py b/tests/unit_tests/test_nuclide_heating.py new file mode 100644 index 0000000000..e00d57b4ff --- /dev/null +++ b/tests/unit_tests/test_nuclide_heating.py @@ -0,0 +1,39 @@ +import openmc +from pytest import approx + + +def test_nuclide_heating(run_in_tmpdir): + mat = openmc.Material() + mat.add_nuclide("Li6", 0.5) + mat.add_nuclide("Li7", 0.5) + mat.set_density("g/cm3", 1.0) + + sphere = openmc.Sphere(r=20, boundary_type="reflective") + inside_sphere = openmc.Cell(fill=mat, region=-sphere) + model = openmc.Model() + model.geometry = openmc.Geometry([inside_sphere]) + + model.settings.particles = 1000 + model.settings.batches = 1 + model.settings.photon_transport = True + model.settings.electron_treatment = "ttb" + model.settings.cutoff = {"energy_photon": 1000} + model.settings.run_mode = "fixed source" + model.settings.source = openmc.IndependentSource( + energy=openmc.stats.delta_function(10.0e6), + particle="photon" + ) + + # Create two tallies, one with heating by nuclide and one with total heating + tally1 = openmc.Tally() + tally1.scores = ["heating"] + tally1.nuclides = mat.get_nuclides() + tally2 = openmc.Tally() + tally2.scores = ["heating"] + model.tallies = [tally1, tally2] + + # Run the model + model.run(apply_tally_results=True) + + # Make sure the heating results are consistent + assert tally1.mean.sum() == approx(tally2.mean.sum())