Merge pull request #967 from smharper/boric_acid

Add convenience function for boric acid Materials

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

docs/source/pythonapi/model.rst

@@ -5,16 +5,12 @@
 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:
    :template: myfunction.rst
 
+   openmc.model.borated_water
    openmc.model.get_hexagonal_prism
    openmc.model.get_rectangular_prism
    openmc.model.subdivide

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,108 @@ 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.  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
+    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 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]
+    # = 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).
 

openmc/model/funcs.py

@@ -5,6 +5,99 @@ from numbers import Real
 
 from openmc import XPlane, YPlane, Plane, ZCylinder
 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.),

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^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)

tests/unit_tests/test_material.py

@@ -139,3 +139,29 @@ 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_borated_water():
+    # Test against reference values from the BEAVRS benchmark.
+    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'
+    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
+
+    # Test the Celsius conversion.
+    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.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.borated_water(975, 566.5, 15.51, density=0.9)
+    assert m.density == pytest.approx(0.9, 1e-3)