From 8afa5ded101cb4fa2268ee7e127f150936f5db3b Mon Sep 17 00:00:00 2001 From: Ethan Peterson Date: Sun, 12 Jan 2020 12:08:28 -0500 Subject: [PATCH] Adding mix_materials method to Material class --- openmc/material.py | 95 ++++++++++++++++++++++++++++++- tests/unit_tests/test_material.py | 24 ++++++++ 2 files changed, 118 insertions(+), 1 deletion(-) diff --git a/openmc/material.py b/openmc/material.py index 375e1dca1e..996d6a3327 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1,4 +1,5 @@ -from collections import OrderedDict +from collections import OrderedDict, defaultdict +from collections.abc import Iterable from copy import deepcopy from numbers import Real, Integral from pathlib import Path @@ -914,6 +915,98 @@ class Material(IDManagerMixin): return element + @classmethod + def mix_materials(cls, materials, fracs, percent_type='ao', name=None): + """Mix materials together based on atom, weight, or volume fractions + + Parameters + ---------- + materials : Iterable of openmc.Material + Materials to combine + fracs : Iterable of float + Fractions of each material to be combined + percent_type : {'ao', 'wo', 'vo'} + Type of percentage, must be one of 'ao', 'wo', or 'vo', to signify atom + percent (molar percent), weight percent, or volume percent, + optional. Defaults to 'ao' + name : str + The name for the new material, optional. Defaults to concatenated + names of input materials with percentages indicated inside + parentheses. + + Returns + ------- + openmc.Material + Mixture of the materials + + """ + + cv.check_type('materials', materials, Iterable, Material) + cv.check_type('fracs', fracs, Iterable, Real) + cv.check_value('percent type', percent_type, {'ao', 'wo', 'vo'}) + + fracs = np.asarray(fracs) + void_frac = 1. - np.sum(fracs) + + # Warn that fractions don't add to 1, set remainder to void, or raise + # an error if percent_type isn't 'vo' + if not np.isclose(void_frac, 0.): + if percent_type in ('ao', 'wo'): + msg = ('A non-zero void fraction is not acceptable for ' + 'percent_type: {}'.format(percent_type)) + raise ValueError(msg) + else: + msg = ('Warning: sum of fractions do not add to 1, void ' + 'fraction set to {}'.format(void_frac)) + warnings.warn(msg) + + # Calculate appropriate weights which are how many cc's of each + # material are found in 1cc of the composite material + amms = np.asarray([mat.average_molar_mass for mat in materials]) + mass_dens = np.asarray([mat.get_mass_density() for mat in materials]) + if percent_type == 'ao': + wgts = fracs * amms / mass_dens + wgts /= np.sum(wgts) + elif percent_type == 'wo': + wgts = fracs / mass_dens + wgts /= np.sum(wgts) + elif percent_type == 'vo': + wgts = fracs + + # If any of the involved materials contain S(a,b) tables raise an error + sab_names = set(sab[0] for mat in materials for sab in mat._sab) + if sab_names: + msg = ('Currently we do not support mixing materials containing ' + 'S(a,b) tables') + raise NotImplementedError(msg) + + # Add nuclide densities weighted by appropriate fractions + nuclides_per_cc = defaultdict(float) + mass_per_cc = defaultdict(float) + for mat, wgt in zip(materials, wgts): + for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().values(): + nuc_per_cc = wgt*1.e24*atoms_per_bcm + nuclides_per_cc[nuc] += nuc_per_cc + mass_per_cc[nuc] += nuc_per_cc*openmc.data.atomic_mass(nuc) / \ + openmc.data.AVOGADRO + + # Create the new material with the desired name + if name is None: + name = '-'.join(['{}({})'.format(m.name, f) for m, f in + zip(materials, fracs)]) + new_mat = openmc.Material(name=name) + + # Compute atom fractions of nuclides and add them to the new material + tot_nuclides_per_cc = np.sum([dens for dens in nuclides_per_cc.values()]) + for nuc, atom_dens in nuclides_per_cc.items(): + new_mat.add_nuclide(nuc, atom_dens/tot_nuclides_per_cc, 'ao') + + # Compute mass density for the new material and set it + new_density = np.sum([dens for dens in mass_per_cc.values()]) + new_mat.set_density('g/cm3', new_density) + + return new_mat + @classmethod def from_xml_element(cls, elem): """Generate material from an XML element diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 62fff7cdf1..cbff16e403 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -218,3 +218,27 @@ def test_from_xml(run_in_tmpdir): assert m2.density == 10.0 assert m2.density_units == 'kg/m3' assert mats[2].density_units == 'sum' + + +def test_mix_materials(): + m1 = openmc.Material() + m1.add_nuclide('U235', 1.) + m1dens = 10.0 + m1amm = m1.average_molar_mass + m1.set_density('g/cm3', m1dens) + m2 = openmc.Material() + m2.add_nuclide('Zr90', 1.) + m2dens = 2.0 + m2amm = m2.average_molar_mass + m2.set_density('g/cm3', m2dens) + f0, f1 = 0.6, 0.4 + dens3 = (f0*m1amm + f1*m2amm) / (f0*m1amm/m1dens + f1*m2amm/m2dens) + dens4 = 1. / (f0 / m1dens + f1 / m2dens) + dens5 = f0*m1dens + f1*m2dens + m3 = openmc.Material.mix_materials([m1, m2], [f0, f1], percent_type='ao') + m4 = openmc.Material.mix_materials([m1, m2], [f0, f1], percent_type='wo') + m5 = openmc.Material.mix_materials([m1, m2], [f0, f1], percent_type='vo') + assert m3.density == pytest.approx(dens3) + assert m4.density == pytest.approx(dens4) + assert m5.density == pytest.approx(dens5) +