From bb4ac5a7b2365983bca8250d97d151426b7cbedb Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 10 Feb 2018 01:04:45 -0500 Subject: [PATCH 1/5] Add IAPWS water density data --- openmc/data/data.py | 107 +++++++++++++++++++++++++++++ tests/unit_tests/test_data_misc.py | 10 +++ 2 files changed, 117 insertions(+) diff --git a/openmc/data/data.py b/openmc/data/data.py index ded6870188..8c0ed63728 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -1,6 +1,9 @@ import itertools import os import re +from warnings import warn + +from numpy import sqrt # Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions @@ -208,6 +211,110 @@ def atomic_weight(element): return None if weight == 0. else weight +def water_density(temperature, pressure=0.1013): + """Return the density of liquid water at a given temperature and pressure. + + The density is calculated from a polynomial fit using equations and values + from the 2012 version of the IAPWS-IF97 formulation. Only the equations + for region 1 are implemented here. + + Results are invalid for water vapor; pressures above 100 [MPa]; and + temperatures below 273.15 [K], above 623.15 [K], or above saturation. + + Reference: International Association for the Properties of Water and Steam, + "Revised Release on the IAPWS Industrial Formulation 1997 for the + Thermodynamic Properties of Water and Steam", IAPWS R7-97(2012). + + Parameters + ---------- + temperature : float + Water temperature in units of [K] + pressure : float + Water pressure in units of [MPa] + + Returns + ------- + float + Water density in units of [g / cm^3] + + """ + + # Make sure the temperature and pressure are inside the min/max region 1 + # bounds. (Relax the 273.15 bound to 273 in case a user wants 0 deg C data + # but they only use 3 digits for their conversion to K.) + if pressure > 100.0: + warn("Results are not valid for pressures above 100 MPa.") + if pressure < 0.0: + warn("Results are not valid for pressures below zero.") + if temperature < 273: + warn("Results are not valid for temperatures below 273.15 K.") + if temperature > 623.15: + warn("Results are not valid for temperatures above 623.15 K.") + + # IAPWS region 4 parameters + _n4 = [0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2, + 0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2, + -0.48232657361591e4, 0.40511340542057e6, -0.23855557567849, + 0.65017534844798e3] + + # Compute the saturation temperature at the given pressure. + beta = pressure**(0.25) + E = beta**2 + _n4[2] * beta + _n4[5] + F = _n4[0] * beta**2 + _n4[3] * beta + _n4[6] + G = _n4[1] * beta**2 + _n4[4] * beta + _n4[7] + D = 2.0 * G / (-F - sqrt(F**2 - 4 * E * G)) + T_sat = 0.5 * (_n4[9] + D + - sqrt((_n4[9] + D)**2 - 4.0 * (_n4[8] + _n4[9] * D))) + + # Make sure we aren't above saturation. (Relax this bound by .2 degrees + # for deg C to K conversions.) + if temperature > T_sat + 0.2: + warn("Results are not valid for temperatures above saturation " + "(above the boiling point).") + + # IAPWS region 1 parameters + R_GAS_CONSTANT = 0.461526 # kJ / kg / K + _ref_p = 16.53 # MPa + _ref_T = 1386 # K + _n1f = [0.14632971213167, -0.84548187169114, -0.37563603672040e1, + 0.33855169168385e1, -0.95791963387872, 0.15772038513228, + -0.16616417199501e-1, 0.81214629983568e-3, 0.28319080123804e-3, + -0.60706301565874e-3, -0.18990068218419e-1, -0.32529748770505e-1, + -0.21841717175414e-1, -0.52838357969930e-4, -0.47184321073267e-3, + -0.30001780793026e-3, 0.47661393906987e-4, -0.44141845330846e-5, + -0.72694996297594e-15, -0.31679644845054e-4, -0.28270797985312e-5, + -0.85205128120103e-9, -0.22425281908000e-5, -0.65171222895601e-6, + -0.14341729937924e-12, -0.40516996860117e-6, -0.12734301741641e-8, + -0.17424871230634e-9, -0.68762131295531e-18, 0.14478307828521e-19, + 0.26335781662795e-22, -0.11947622640071e-22, 0.18228094581404e-23, + -0.93537087292458e-25] + _I1f = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, + 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32] + _J1f = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, + 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41] + + # Nondimensionalize the pressure and temperature. + pi = pressure / _ref_p + tau = _ref_T / temperature + + # Compute the derivative of gamma (dimensionless Gibbs free energy) with + # respect to pi. + gamma1_pi = 0.0 + for i in range(34): + gamma1_pi -= (_n1f[i] * _I1f[i] * (7.1 - pi)**(_I1f[i] - 1) + * (tau - 1.222)**_J1f[i]) + + # Compute the leading coefficient. This sets the units at + # 1 [MPa] * [kg K / kJ] / [1 / K] + # = 1e6 [N / m^2] * 1e-3 [kg K / N / m] * [1 / K] + # = 1e3 [kg / m^3] + # = 1 [g / cm^3] + coeff = pressure / R_GAS_CONSTANT / temperature + + # Compute and return the density. + return coeff / pi / gamma1_pi + + # Values here are from the Committee on Data for Science and Technology # (CODATA) 2014 recommendation (doi:10.1103/RevModPhys.88.035009). diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py index aeeb04c0fb..b33ef996fe 100644 --- a/tests/unit_tests/test_data_misc.py +++ b/tests/unit_tests/test_data_misc.py @@ -48,3 +48,13 @@ def test_thin(): x_thin, y_thin = openmc.data.thin(x, y) f = openmc.data.Tabulated1D(x_thin, y_thin) assert f(1.0) == pytest.approx(np.sin(1.0), 0.001) + + +def test_water_density(): + dens = openmc.data.water_density + # These test values are from IAPWS R7-97(2012). They are actually specific + # volumes so they need to be inverted. They also need to be divided by 1000 + # to convert from [kg / m^3] to [g / cm^2]. + assert dens(300.0, 3.0) == pytest.approx(1e-3/0.100215168e-2, 1e-6) + assert dens(300.0, 80.0) == pytest.approx(1e-3/0.971180894e-3, 1e-6) + assert dens(500.0, 3.0) == pytest.approx(1e-3/0.120241800e-2, 1e-6) From a06190ee407c14fd7a82cd0edf2fd89ce7e956d1 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 10 Feb 2018 02:23:50 -0500 Subject: [PATCH 2/5] Add convenience function for boric acid Materials --- openmc/material.py | 70 +++++++++++++++++++++++++++++++ tests/unit_tests/test_material.py | 18 ++++++++ 2 files changed, 88 insertions(+) diff --git a/openmc/material.py b/openmc/material.py index e409d6536c..32b7d3a80d 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -973,3 +973,73 @@ class Materials(cv.CheckedList): # Write the XML Tree to the materials.xml file tree = ET.ElementTree(root_element) tree.write(path, xml_declaration=True, encoding='utf-8') + + +def make_boric_acid(boron_ppm, temperature=293., pressure=0.1013, density=None, + **kwargs): + """Return a Material with the composition of boric acid. + + The water density can either be given directly, or it can be determined from + a temperature and pressure. + + Parameters + ---------- + boron_ppm : float + The weight fraction in parts-per-million of elemental boron in the acid. + temperature : float + Water temperature in [K] used to compute water density. + pressure : float + Water pressure in [MPa] used to compute water density. + density : float + Water density in [g / cm^3]. If specified, this value overrides the + temperature and pressure arguments. + **kwargs + All keyword arguments are passed to the created Material object. + + Returns + ------- + openmc.Material + + """ + # Set the density of water, either from an explicitly given density or from + # temperature and pressure. + if density is not None: + water_density = density + else: + water_density = openmc.data.water_density(temperature, pressure) + + # Compute the density of the boric acid. + acid_density = water_density / (1 - boron_ppm * 1e-6) + + # Compute the molar mass of pure water. + hydrogen = openmc.Element('H') + oxygen = openmc.Element('O') + M_H2O = 0.0 + for iso_name, frac, junk in hydrogen.expand(2.0, 'ao'): + M_H2O += frac * openmc.data.atomic_mass(iso_name) + for iso_name, frac, junk in oxygen.expand(1.0, 'ao'): + M_H2O += frac * openmc.data.atomic_mass(iso_name) + + # Compute the molar mass of boron. + boron = openmc.Element('B') + M_B = 0.0 + for iso_name, frac, junk in boron.expand(1.0, 'ao'): + M_B += frac * openmc.data.atomic_mass(iso_name) + + # Compute the number fractions of each element. + frac_H2O = (1 - boron_ppm * 1e-6) / M_H2O + frac_H = 2 * frac_H2O + frac_O = frac_H2O + frac_B = boron_ppm * 1e-6 / M_B + + # Build the material. + if density is None: + out = openmc.Material(temperature=temperature, **kwargs) + else: + out = openmc.Material(**kwargs) + out.add_element('H', frac_H, 'ao') + out.add_element('O', frac_O, 'ao') + out.add_element('B', frac_B, 'ao') + out.set_density('g/cc', acid_density) + out.add_s_alpha_beta('c_H_in_H2O') + return out diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index e92a2ac088..38495cb1ab 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -139,3 +139,21 @@ def test_materials(run_in_tmpdir): mats.cross_sections = '/some/fake/cross_sections.xml' mats.multipole_library = '/some/awesome/mp_lib/' mats.export_to_xml() + + +def test_boric_acid(): + # Test against reference values from the BEAVRS benchmark. + m = openmc.make_boric_acid(975, 566.5, 15.51, material_id=50) + assert m.density == pytest.approx(0.7405, 1e-3) + assert m.temperature == pytest.approx(566.5) + assert m._sab[0][0] == 'c_H_in_H2O' + ref_dens = {'B10':8.0023e-06, 'B11':3.2210e-05, 'H1':4.9458e-02, + 'O16':2.4672e-02} + nuc_dens = m.get_nuclide_atom_densities() + for nuclide in ref_dens: + assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert m.id == 50 + + # Make sure the density override works + m = openmc.make_boric_acid(975, 566.5, 15.51, 0.9) + assert m.density == pytest.approx(0.9, 1e-3) From 9c53c085d6b81125f8cd3835ee40c2fc7f207f63 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 11 Feb 2018 15:01:51 -0500 Subject: [PATCH 3/5] Change boric_acid to borated_water; add more units --- openmc/data/data.py | 7 +-- openmc/material.py | 70 ----------------------- openmc/model/funcs.py | 95 +++++++++++++++++++++++++++++++ tests/unit_tests/test_material.py | 16 ++++-- 4 files changed, 110 insertions(+), 78 deletions(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index 8c0ed63728..45a4406dbe 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -216,10 +216,9 @@ def water_density(temperature, pressure=0.1013): The density is calculated from a polynomial fit using equations and values from the 2012 version of the IAPWS-IF97 formulation. Only the equations - for region 1 are implemented here. - - Results are invalid for water vapor; pressures above 100 [MPa]; and - temperatures below 273.15 [K], above 623.15 [K], or above saturation. + for region 1 are implemented here. Region 1 is limited to liquid water + below 100 [MPa] with a temperature above 273.15 [K], below 623.15 [K], and + below saturation. Reference: International Association for the Properties of Water and Steam, "Revised Release on the IAPWS Industrial Formulation 1997 for the diff --git a/openmc/material.py b/openmc/material.py index 32b7d3a80d..e409d6536c 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -973,73 +973,3 @@ class Materials(cv.CheckedList): # Write the XML Tree to the materials.xml file tree = ET.ElementTree(root_element) tree.write(path, xml_declaration=True, encoding='utf-8') - - -def make_boric_acid(boron_ppm, temperature=293., pressure=0.1013, density=None, - **kwargs): - """Return a Material with the composition of boric acid. - - The water density can either be given directly, or it can be determined from - a temperature and pressure. - - Parameters - ---------- - boron_ppm : float - The weight fraction in parts-per-million of elemental boron in the acid. - temperature : float - Water temperature in [K] used to compute water density. - pressure : float - Water pressure in [MPa] used to compute water density. - density : float - Water density in [g / cm^3]. If specified, this value overrides the - temperature and pressure arguments. - **kwargs - All keyword arguments are passed to the created Material object. - - Returns - ------- - openmc.Material - - """ - # Set the density of water, either from an explicitly given density or from - # temperature and pressure. - if density is not None: - water_density = density - else: - water_density = openmc.data.water_density(temperature, pressure) - - # Compute the density of the boric acid. - acid_density = water_density / (1 - boron_ppm * 1e-6) - - # Compute the molar mass of pure water. - hydrogen = openmc.Element('H') - oxygen = openmc.Element('O') - M_H2O = 0.0 - for iso_name, frac, junk in hydrogen.expand(2.0, 'ao'): - M_H2O += frac * openmc.data.atomic_mass(iso_name) - for iso_name, frac, junk in oxygen.expand(1.0, 'ao'): - M_H2O += frac * openmc.data.atomic_mass(iso_name) - - # Compute the molar mass of boron. - boron = openmc.Element('B') - M_B = 0.0 - for iso_name, frac, junk in boron.expand(1.0, 'ao'): - M_B += frac * openmc.data.atomic_mass(iso_name) - - # Compute the number fractions of each element. - frac_H2O = (1 - boron_ppm * 1e-6) / M_H2O - frac_H = 2 * frac_H2O - frac_O = frac_H2O - frac_B = boron_ppm * 1e-6 / M_B - - # Build the material. - if density is None: - out = openmc.Material(temperature=temperature, **kwargs) - else: - out = openmc.Material(**kwargs) - out.add_element('H', frac_H, 'ao') - out.add_element('O', frac_O, 'ao') - out.add_element('B', frac_B, 'ao') - out.set_density('g/cc', acid_density) - out.add_s_alpha_beta('c_H_in_H2O') - return out diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index 6013f6caee..b8984ada03 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -5,6 +5,7 @@ from numbers import Real from openmc import XPlane, YPlane, Plane, ZCylinder from openmc.checkvalue import check_type, check_value +import openmc.data def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), @@ -252,6 +253,100 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), return prism +def make_borated_water(boron_ppm, temperature=293., pressure=0.1013, + temp_unit='K', press_unit='MPa', density=None, **kwargs): + """Return a Material with the composition of boron dissolved in water. + + The water density can be determined from a temperature and pressure, or it + can be set directly. + + The concentration of boron has no effect on the stoichometric ratio of H + and O---they are fixed at 2-1. + + Parameters + ---------- + boron_ppm : float + The weight fraction in parts-per-million of elemental boron in the + water. + temperature : float + Temperature in [K] used to compute water density. + pressure : float + Pressure in [MPa] used to compute water density. + temp_unit : str + The units used for the `temperature` argument. Valid units are 'K', + 'C', and 'F'. + press_unit : str + The units used for the `pressure` argument. Valid units are 'MPa' and + 'psi'. + density : float + Water density in [g / cm^3]. If specified, this value overrides the + temperature and pressure arguments. + **kwargs + All keyword arguments are passed to the created Material object. + + Returns + ------- + openmc.Material + + """ + # Perform any necessary unit conversions. + check_value('temperature unit', temp_unit, ('K', 'C', 'F')) + if temp_unit == 'K': + T = temperature + elif temp_unit == 'C': + T = temperature + 273.15 + elif temp_unit == 'F': + T = (temperature + 459.67) * 5.0 / 9.0 + check_value('pressure unit', press_unit, ('MPa', 'psi')) + if press_unit == 'MPa': + P = pressure + elif press_unit == 'psi': + P = pressure * 0.006895 + + # Set the density of water, either from an explicitly given density or from + # temperature and pressure. + if density is not None: + water_density = density + else: + water_density = openmc.data.water_density(T, P) + + # Compute the density of the solution. + solution_density = water_density / (1 - boron_ppm * 1e-6) + + # Compute the molar mass of pure water. + hydrogen = openmc.Element('H') + oxygen = openmc.Element('O') + M_H2O = 0.0 + for iso_name, frac, junk in hydrogen.expand(2.0, 'ao'): + M_H2O += frac * openmc.data.atomic_mass(iso_name) + for iso_name, frac, junk in oxygen.expand(1.0, 'ao'): + M_H2O += frac * openmc.data.atomic_mass(iso_name) + + # Compute the molar mass of boron. + boron = openmc.Element('B') + M_B = 0.0 + for iso_name, frac, junk in boron.expand(1.0, 'ao'): + M_B += frac * openmc.data.atomic_mass(iso_name) + + # Compute the number fractions of each element. + frac_H2O = (1 - boron_ppm * 1e-6) / M_H2O + frac_H = 2 * frac_H2O + frac_O = frac_H2O + frac_B = boron_ppm * 1e-6 / M_B + + # Build the material. + if density is None: + out = openmc.Material(temperature=T, **kwargs) + else: + out = openmc.Material(**kwargs) + out.add_element('H', frac_H, 'ao') + out.add_element('O', frac_O, 'ao') + out.add_element('B', frac_B, 'ao') + out.set_density('g/cc', solution_density) + out.add_s_alpha_beta('c_H_in_H2O') + return out + + def subdivide(surfaces): """Create regions separated by a series of surfaces. diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 38495cb1ab..295153a9c1 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -141,9 +141,9 @@ def test_materials(run_in_tmpdir): mats.export_to_xml() -def test_boric_acid(): +def test_borated_water(): # Test against reference values from the BEAVRS benchmark. - m = openmc.make_boric_acid(975, 566.5, 15.51, material_id=50) + m = openmc.model.make_borated_water(975, 566.5, 15.51, material_id=50) assert m.density == pytest.approx(0.7405, 1e-3) assert m.temperature == pytest.approx(566.5) assert m._sab[0][0] == 'c_H_in_H2O' @@ -154,6 +154,14 @@ def test_boric_acid(): assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) assert m.id == 50 - # Make sure the density override works - m = openmc.make_boric_acid(975, 566.5, 15.51, 0.9) + # Test the Celsius conversion. + m = openmc.model.make_borated_water(975, 293.35, 15.51, 'C') + assert m.density == pytest.approx(0.7405, 1e-3) + + # Test Fahrenheit and psi conversions. + m = openmc.model.make_borated_water(975, 560.0, 2250.0, 'F', 'psi') + assert m.density == pytest.approx(0.7405, 1e-3) + + # Test the density override + m = openmc.model.make_borated_water(975, 566.5, 15.51, density=0.9) assert m.density == pytest.approx(0.9, 1e-3) From e0d6640c4fb43edf3213efe62c05c1106575c5d9 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 11 Feb 2018 15:03:48 -0500 Subject: [PATCH 4/5] Add make_borated_water and water_density to docs --- docs/source/pythonapi/data.rst | 1 + docs/source/pythonapi/model.rst | 6 +----- tests/unit_tests/test_data_misc.py | 2 +- 3 files changed, 3 insertions(+), 6 deletions(-) diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 7f75ddaaaf..3c221906db 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -34,6 +34,7 @@ Core Functions openmc.data.atomic_mass openmc.data.linearize openmc.data.thin + openmc.data.water_density openmc.data.write_compact_458_library Angle-Energy Distributions diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 4ba247468a..6aff6d4c25 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -5,11 +5,6 @@ Convenience Functions --------------------- -Several helper functions are available here. Ther first two create rectangular -and hexagonal prisms defined by the intersection of four and six surface -half-spaces, respectively. The last function takes a sequence of surfaces and -returns the regions that separate them. - .. autosummary:: :toctree: generated :nosignatures: @@ -17,6 +12,7 @@ returns the regions that separate them. openmc.model.get_hexagonal_prism openmc.model.get_rectangular_prism + openmc.model.make_borated_water openmc.model.subdivide TRISO Fuel Modeling diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py index b33ef996fe..433d34adb9 100644 --- a/tests/unit_tests/test_data_misc.py +++ b/tests/unit_tests/test_data_misc.py @@ -54,7 +54,7 @@ def test_water_density(): dens = openmc.data.water_density # These test values are from IAPWS R7-97(2012). They are actually specific # volumes so they need to be inverted. They also need to be divided by 1000 - # to convert from [kg / m^3] to [g / cm^2]. + # to convert from [kg / m^3] to [g / cm^3]. assert dens(300.0, 3.0) == pytest.approx(1e-3/0.100215168e-2, 1e-6) assert dens(300.0, 80.0) == pytest.approx(1e-3/0.971180894e-3, 1e-6) assert dens(500.0, 3.0) == pytest.approx(1e-3/0.120241800e-2, 1e-6) From d34bdf9b98bd2c06da5e5282d9516cce041be4cc Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 11 Feb 2018 18:39:06 -0500 Subject: [PATCH 5/5] Address #967 comments --- docs/source/pythonapi/model.rst | 2 +- openmc/data/data.py | 65 +++++------ openmc/model/funcs.py | 186 +++++++++++++++--------------- tests/unit_tests/test_material.py | 8 +- 4 files changed, 129 insertions(+), 132 deletions(-) diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst index 6aff6d4c25..9ed77f3666 100644 --- a/docs/source/pythonapi/model.rst +++ b/docs/source/pythonapi/model.rst @@ -10,9 +10,9 @@ Convenience Functions :nosignatures: :template: myfunction.rst + openmc.model.borated_water openmc.model.get_hexagonal_prism openmc.model.get_rectangular_prism - openmc.model.make_borated_water openmc.model.subdivide TRISO Fuel Modeling diff --git a/openmc/data/data.py b/openmc/data/data.py index 45a4406dbe..a7c0e536f6 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -251,19 +251,19 @@ def water_density(temperature, pressure=0.1013): warn("Results are not valid for temperatures above 623.15 K.") # IAPWS region 4 parameters - _n4 = [0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2, - 0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2, - -0.48232657361591e4, 0.40511340542057e6, -0.23855557567849, - 0.65017534844798e3] + n4 = [0.11670521452767e4, -0.72421316703206e6, -0.17073846940092e2, + 0.12020824702470e5, -0.32325550322333e7, 0.14915108613530e2, + -0.48232657361591e4, 0.40511340542057e6, -0.23855557567849, + 0.65017534844798e3] # Compute the saturation temperature at the given pressure. beta = pressure**(0.25) - E = beta**2 + _n4[2] * beta + _n4[5] - F = _n4[0] * beta**2 + _n4[3] * beta + _n4[6] - G = _n4[1] * beta**2 + _n4[4] * beta + _n4[7] + E = beta**2 + n4[2] * beta + n4[5] + F = n4[0] * beta**2 + n4[3] * beta + n4[6] + G = n4[1] * beta**2 + n4[4] * beta + n4[7] D = 2.0 * G / (-F - sqrt(F**2 - 4 * E * G)) - T_sat = 0.5 * (_n4[9] + D - - sqrt((_n4[9] + D)**2 - 4.0 * (_n4[8] + _n4[9] * D))) + T_sat = 0.5 * (n4[9] + D + - sqrt((n4[9] + D)**2 - 4.0 * (n4[8] + n4[9] * D))) # Make sure we aren't above saturation. (Relax this bound by .2 degrees # for deg C to K conversions.) @@ -273,38 +273,37 @@ def water_density(temperature, pressure=0.1013): # IAPWS region 1 parameters R_GAS_CONSTANT = 0.461526 # kJ / kg / K - _ref_p = 16.53 # MPa - _ref_T = 1386 # K - _n1f = [0.14632971213167, -0.84548187169114, -0.37563603672040e1, - 0.33855169168385e1, -0.95791963387872, 0.15772038513228, - -0.16616417199501e-1, 0.81214629983568e-3, 0.28319080123804e-3, - -0.60706301565874e-3, -0.18990068218419e-1, -0.32529748770505e-1, - -0.21841717175414e-1, -0.52838357969930e-4, -0.47184321073267e-3, - -0.30001780793026e-3, 0.47661393906987e-4, -0.44141845330846e-5, - -0.72694996297594e-15, -0.31679644845054e-4, -0.28270797985312e-5, - -0.85205128120103e-9, -0.22425281908000e-5, -0.65171222895601e-6, - -0.14341729937924e-12, -0.40516996860117e-6, -0.12734301741641e-8, - -0.17424871230634e-9, -0.68762131295531e-18, 0.14478307828521e-19, - 0.26335781662795e-22, -0.11947622640071e-22, 0.18228094581404e-23, - -0.93537087292458e-25] - _I1f = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, - 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32] - _J1f = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, - 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41] + ref_p = 16.53 # MPa + ref_T = 1386 # K + n1f = [0.14632971213167, -0.84548187169114, -0.37563603672040e1, + 0.33855169168385e1, -0.95791963387872, 0.15772038513228, + -0.16616417199501e-1, 0.81214629983568e-3, 0.28319080123804e-3, + -0.60706301565874e-3, -0.18990068218419e-1, -0.32529748770505e-1, + -0.21841717175414e-1, -0.52838357969930e-4, -0.47184321073267e-3, + -0.30001780793026e-3, 0.47661393906987e-4, -0.44141845330846e-5, + -0.72694996297594e-15, -0.31679644845054e-4, -0.28270797985312e-5, + -0.85205128120103e-9, -0.22425281908000e-5, -0.65171222895601e-6, + -0.14341729937924e-12, -0.40516996860117e-6, -0.12734301741641e-8, + -0.17424871230634e-9, -0.68762131295531e-18, 0.14478307828521e-19, + 0.26335781662795e-22, -0.11947622640071e-22, 0.18228094581404e-23, + -0.93537087292458e-25] + I1f = [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 4, + 4, 4, 5, 8, 8, 21, 23, 29, 30, 31, 32] + J1f = [-2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0, 1, 3, -3, 0, 1, 3, 17, -4, + 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41] # Nondimensionalize the pressure and temperature. - pi = pressure / _ref_p - tau = _ref_T / temperature + pi = pressure / ref_p + tau = ref_T / temperature # Compute the derivative of gamma (dimensionless Gibbs free energy) with # respect to pi. gamma1_pi = 0.0 - for i in range(34): - gamma1_pi -= (_n1f[i] * _I1f[i] * (7.1 - pi)**(_I1f[i] - 1) - * (tau - 1.222)**_J1f[i]) + for n, I, J in zip(n1f, I1f, J1f): + gamma1_pi -= n * I * (7.1 - pi)**(I - 1) * (tau - 1.222)**J # Compute the leading coefficient. This sets the units at - # 1 [MPa] * [kg K / kJ] / [1 / K] + # 1 [MPa] * [kg K / kJ] * [1 / K] # = 1e6 [N / m^2] * 1e-3 [kg K / N / m] * [1 / K] # = 1e3 [kg / m^3] # = 1 [g / cm^3] diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index b8984ada03..e974f93b1b 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -8,6 +8,98 @@ from openmc.checkvalue import check_type, check_value import openmc.data +def borated_water(boron_ppm, temperature=293., pressure=0.1013, temp_unit='K', + press_unit='MPa', density=None, **kwargs): + """Return a Material with the composition of boron dissolved in water. + + The water density can be determined from a temperature and pressure, or it + can be set directly. + + The concentration of boron has no effect on the stoichiometric ratio of H + and O---they are fixed at 2-1. + + Parameters + ---------- + boron_ppm : float + The weight fraction in parts-per-million of elemental boron in the + water. + temperature : float + Temperature in [K] used to compute water density. + pressure : float + Pressure in [MPa] used to compute water density. + temp_unit : {'K', 'C', 'F'} + The units used for the `temperature` argument. + press_unit : {'MPa', 'psi'} + The units used for the `pressure` argument. + density : float + Water density in [g / cm^3]. If specified, this value overrides the + temperature and pressure arguments. + **kwargs + All keyword arguments are passed to the created Material object. + + Returns + ------- + openmc.Material + + """ + # Perform any necessary unit conversions. + check_value('temperature unit', temp_unit, ('K', 'C', 'F')) + if temp_unit == 'K': + T = temperature + elif temp_unit == 'C': + T = temperature + 273.15 + elif temp_unit == 'F': + T = (temperature + 459.67) * 5.0 / 9.0 + check_value('pressure unit', press_unit, ('MPa', 'psi')) + if press_unit == 'MPa': + P = pressure + elif press_unit == 'psi': + P = pressure * 0.006895 + + # Set the density of water, either from an explicitly given density or from + # temperature and pressure. + if density is not None: + water_density = density + else: + water_density = openmc.data.water_density(T, P) + + # Compute the density of the solution. + solution_density = water_density / (1 - boron_ppm * 1e-6) + + # Compute the molar mass of pure water. + hydrogen = openmc.Element('H') + oxygen = openmc.Element('O') + M_H2O = 0.0 + for iso_name, frac, junk in hydrogen.expand(2.0, 'ao'): + M_H2O += frac * openmc.data.atomic_mass(iso_name) + for iso_name, frac, junk in oxygen.expand(1.0, 'ao'): + M_H2O += frac * openmc.data.atomic_mass(iso_name) + + # Compute the molar mass of boron. + boron = openmc.Element('B') + M_B = 0.0 + for iso_name, frac, junk in boron.expand(1.0, 'ao'): + M_B += frac * openmc.data.atomic_mass(iso_name) + + # Compute the number fractions of each element. + frac_H2O = (1 - boron_ppm * 1e-6) / M_H2O + frac_H = 2 * frac_H2O + frac_O = frac_H2O + frac_B = boron_ppm * 1e-6 / M_B + + # Build the material. + if density is None: + out = openmc.Material(temperature=T, **kwargs) + else: + out = openmc.Material(**kwargs) + out.add_element('H', frac_H, 'ao') + out.add_element('O', frac_O, 'ao') + out.add_element('B', frac_B, 'ao') + out.set_density('g/cc', solution_density) + out.add_s_alpha_beta('c_H_in_H2O') + return out + + def get_rectangular_prism(width, height, axis='z', origin=(0., 0.), boundary_type='transmission', corner_radius=0.): """Get an infinite rectangular prism from four planar surfaces. @@ -253,100 +345,6 @@ def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.), return prism -def make_borated_water(boron_ppm, temperature=293., pressure=0.1013, - temp_unit='K', press_unit='MPa', density=None, **kwargs): - """Return a Material with the composition of boron dissolved in water. - - The water density can be determined from a temperature and pressure, or it - can be set directly. - - The concentration of boron has no effect on the stoichometric ratio of H - and O---they are fixed at 2-1. - - Parameters - ---------- - boron_ppm : float - The weight fraction in parts-per-million of elemental boron in the - water. - temperature : float - Temperature in [K] used to compute water density. - pressure : float - Pressure in [MPa] used to compute water density. - temp_unit : str - The units used for the `temperature` argument. Valid units are 'K', - 'C', and 'F'. - press_unit : str - The units used for the `pressure` argument. Valid units are 'MPa' and - 'psi'. - density : float - Water density in [g / cm^3]. If specified, this value overrides the - temperature and pressure arguments. - **kwargs - All keyword arguments are passed to the created Material object. - - Returns - ------- - openmc.Material - - """ - # Perform any necessary unit conversions. - check_value('temperature unit', temp_unit, ('K', 'C', 'F')) - if temp_unit == 'K': - T = temperature - elif temp_unit == 'C': - T = temperature + 273.15 - elif temp_unit == 'F': - T = (temperature + 459.67) * 5.0 / 9.0 - check_value('pressure unit', press_unit, ('MPa', 'psi')) - if press_unit == 'MPa': - P = pressure - elif press_unit == 'psi': - P = pressure * 0.006895 - - # Set the density of water, either from an explicitly given density or from - # temperature and pressure. - if density is not None: - water_density = density - else: - water_density = openmc.data.water_density(T, P) - - # Compute the density of the solution. - solution_density = water_density / (1 - boron_ppm * 1e-6) - - # Compute the molar mass of pure water. - hydrogen = openmc.Element('H') - oxygen = openmc.Element('O') - M_H2O = 0.0 - for iso_name, frac, junk in hydrogen.expand(2.0, 'ao'): - M_H2O += frac * openmc.data.atomic_mass(iso_name) - for iso_name, frac, junk in oxygen.expand(1.0, 'ao'): - M_H2O += frac * openmc.data.atomic_mass(iso_name) - - # Compute the molar mass of boron. - boron = openmc.Element('B') - M_B = 0.0 - for iso_name, frac, junk in boron.expand(1.0, 'ao'): - M_B += frac * openmc.data.atomic_mass(iso_name) - - # Compute the number fractions of each element. - frac_H2O = (1 - boron_ppm * 1e-6) / M_H2O - frac_H = 2 * frac_H2O - frac_O = frac_H2O - frac_B = boron_ppm * 1e-6 / M_B - - # Build the material. - if density is None: - out = openmc.Material(temperature=T, **kwargs) - else: - out = openmc.Material(**kwargs) - out.add_element('H', frac_H, 'ao') - out.add_element('O', frac_O, 'ao') - out.add_element('B', frac_B, 'ao') - out.set_density('g/cc', solution_density) - out.add_s_alpha_beta('c_H_in_H2O') - return out - - def subdivide(surfaces): """Create regions separated by a series of surfaces. diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 295153a9c1..2265215417 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -143,7 +143,7 @@ def test_materials(run_in_tmpdir): def test_borated_water(): # Test against reference values from the BEAVRS benchmark. - m = openmc.model.make_borated_water(975, 566.5, 15.51, material_id=50) + m = openmc.model.borated_water(975, 566.5, 15.51, material_id=50) assert m.density == pytest.approx(0.7405, 1e-3) assert m.temperature == pytest.approx(566.5) assert m._sab[0][0] == 'c_H_in_H2O' @@ -155,13 +155,13 @@ def test_borated_water(): assert m.id == 50 # Test the Celsius conversion. - m = openmc.model.make_borated_water(975, 293.35, 15.51, 'C') + m = openmc.model.borated_water(975, 293.35, 15.51, 'C') assert m.density == pytest.approx(0.7405, 1e-3) # Test Fahrenheit and psi conversions. - m = openmc.model.make_borated_water(975, 560.0, 2250.0, 'F', 'psi') + m = openmc.model.borated_water(975, 560.0, 2250.0, 'F', 'psi') assert m.density == pytest.approx(0.7405, 1e-3) # Test the density override - m = openmc.model.make_borated_water(975, 566.5, 15.51, density=0.9) + m = openmc.model.borated_water(975, 566.5, 15.51, density=0.9) assert m.density == pytest.approx(0.9, 1e-3)