From c941d0bc554ec1c973c7f4e6538c8667683c15f5 Mon Sep 17 00:00:00 2001 From: joergbuchwald <joerg.buchwald@ufz.de> Date: Thu, 7 Jul 2022 17:30:47 +0200 Subject: [PATCH] adding library for nonlinear watermodels --- Physics/water/LICENSE | 27 ++++++ Physics/water/README.md | 3 + Physics/water/setup.py | 24 ++++++ Physics/water/water/__init__.py | 0 Physics/water/water/properties/__init__.py | 0 .../water/water/properties/conductivity.py | 47 +++++++++++ Physics/water/water/properties/density.py | 83 +++++++++++++++++++ Physics/water/water/properties/expansivity.py | 50 +++++++++++ .../water/properties/specificheatcapacity.py | 47 +++++++++++ Physics/water/water/properties/template.py | 58 +++++++++++++ Physics/water/water/properties/viscosity.py | 79 ++++++++++++++++++ Physics/water/water/water.py | 32 +++++++ 12 files changed, 450 insertions(+) create mode 100644 Physics/water/LICENSE create mode 100644 Physics/water/README.md create mode 100644 Physics/water/setup.py create mode 100644 Physics/water/water/__init__.py create mode 100644 Physics/water/water/properties/__init__.py create mode 100644 Physics/water/water/properties/conductivity.py create mode 100644 Physics/water/water/properties/density.py create mode 100644 Physics/water/water/properties/expansivity.py create mode 100644 Physics/water/water/properties/specificheatcapacity.py create mode 100644 Physics/water/water/properties/template.py create mode 100644 Physics/water/water/properties/viscosity.py create mode 100644 Physics/water/water/water.py diff --git a/Physics/water/LICENSE b/Physics/water/LICENSE new file mode 100644 index 000000000..7d468793b --- /dev/null +++ b/Physics/water/LICENSE @@ -0,0 +1,27 @@ +Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. +2. Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. +3. Neither the name of the OpenGeoSys Community nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. +IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, +INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT +NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, +OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, +WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY +OF SUCH DAMAGE. + + diff --git a/Physics/water/README.md b/Physics/water/README.md new file mode 100644 index 000000000..bf5ab7ed7 --- /dev/null +++ b/Physics/water/README.md @@ -0,0 +1,3 @@ +# water.py +class for creating nonlinear watermodels to be used in OGS +Most data taken from https://www.engineeringtoolbox.com \ No newline at end of file diff --git a/Physics/water/setup.py b/Physics/water/setup.py new file mode 100644 index 000000000..e5193d94a --- /dev/null +++ b/Physics/water/setup.py @@ -0,0 +1,24 @@ +# -*- coding: utf-8 -*- +from setuptools import setup + +setup(name="water", + version=0.01, + maintainer="Jörg Buchwald", + maintainer_email="joerg_buchwald@ufz.de", + author="Jörg Buchwald", + author_email="joerg.buchwald@ufz.de", + url="https://github.com/joergbuchwald/mini-projects", + classifiers=["Intended Audience :: Science/Research", + "Topic :: Scientific/Engineering :: Visualization", + "Topic :: Scientific/Engineering :: Physics", + "Topic :: Scientific/Engineering :: Mathematics", + "License :: OSI Approved :: MIT License", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.8"], + license="BSD-3 - see LICENSE.txt", + platforms=["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"], + include_package_data=True, + python_requires='>=3.8', + install_requires=["numpy"], + py_modules=["water/water"], + packages=["water/properties"]) diff --git a/Physics/water/water/__init__.py b/Physics/water/water/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/Physics/water/water/properties/__init__.py b/Physics/water/water/properties/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/Physics/water/water/properties/conductivity.py b/Physics/water/water/properties/conductivity.py new file mode 100644 index 000000000..48cabd151 --- /dev/null +++ b/Physics/water/water/properties/conductivity.py @@ -0,0 +1,47 @@ +""" +Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" +# pylint: disable=C0103, R0902, R0914, R0913 +import numpy as np +from water.properties import template +T0 = 273.15 + +class conductivity_1(template.PROPERTY): + def value(self, temperature): + # source: https://www.engineeringtoolbox.com/water-liquid-gas-thermal-conductivity-temperature-pressure-d_2012.html + temp = np.array([0.01, 10., 20., 30., 40., 50., 60., 70., 80., 90., 99.6])+T0 + if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]): + print("The temperature is not within the defined domain") + K = np.array([ + 0.556, 0.579, 0.598, 0.614, 0.629, 0.641, 0.651, 0.660, 0.667, + 0.673, 0.677 + ]) + return np.interp(temperature, temp, K) + def dvalue(self, temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + raise NotImplementedError + def exprtk_dvalue(self): + raise NotImplementedError + +class conductivity_2(template.PROPERTY): + def value(self, temperature): + # use only one-liner to keep it parsable: + return 3.65759470e-08 * temperature**3 - 4.51285141e-05 * temperature**2 + 1.88275537e-02 * temperature - 1.96475421e+00 + def dvalue(self, temperature): + # use only one-liner to keep it parsable: + return 3.65759470e-08 * 3 * temperature**2 - 4.51285141e-05 * 2 * temperature + 1.88275537e-02 + def dvaluenum(self,temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self): + string = self._getcodeasstring("dvalue",None) + string = self._convertpythontoexprtk(string) + return string diff --git a/Physics/water/water/properties/density.py b/Physics/water/water/properties/density.py new file mode 100644 index 000000000..38a01106a --- /dev/null +++ b/Physics/water/water/properties/density.py @@ -0,0 +1,83 @@ +""" +Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" +# pylint: disable=C0103, R0902, R0914, R0913 +import numpy as np +from water.properties import template +T0 = 273.15 + +class density_1(template.PROPERTY): + def value(self, temperature): + # source: https://www.engineeringtoolbox.com/water-density-specific-weight-d_595.html + temp = np.array([ + 0.1, 1.0, 4., 10., 15., 20., 25., 30., 35., 40., 45., 50., 60., + 70., 80., 90., 100., 110., 140., 200., 260., 300., 360. + ]) + T0 + if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]): + print("The temperature is not within the defined domain") + dens = np.array([ + 999.85, 999.9, 999.97, 999.7, 999.1, 998.2, 997., 995.6, 994., + 992.2, 990.2, 988., 983.2, 977.76, 971.8, 965.3, 958.6, 951.0, + 926.1, 865., 783.6, 712.2, 527.6 + ]) + return np.interp(temperature, temp, dens) + def dvalue(self,temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + raise NotImplementedError + def exprtk_dvalue(self): + raise NotImplementedError + +class density_2(template.PROPERTY): + def value(self, temperature, phase_pressure): + # use only one-liner to keep it parsable: + return 1000.1 * (1 - (-6*10**-6 * (temperature-T0)**4 + 0.001667 * (temperature-T0)**3 + -0.197796 * (temperature-T0)**2 + 16.86446 * (temperature-T0) - 64.319951)/10**6 * (temperature-293.15) + 4.65e-10 * (np.maximum(phase_pressure, 0.0) - 1e5)) + def dvalue(self, temperature, phase_pressure, variable="temperature"): + if variable == "temperature": + # use only one-liner to keep it parsable: + return -1000.1 * ((-6*10**-6 * (temperature-T0)**4 + 0.001667 * (temperature-T0)**3 + -0.197796 * (temperature-T0)**2 + 16.86446 * (temperature-T0) - 64.319951)/10**6+(-6*10**-6 * 4 * (temperature-T0)**3 + 0.001667 * 3 * (temperature-T0)**2 + -0.197796 * 2 * (temperature-T0) + 16.86446 )/10**6*(temperature-293.15)) + elif variable == "phase_pressure": + # use only one-liner to keep it parsable: + return np.maximum(np.sign(phase_pressure),0)* 1000.1 * 4.65e-10 + def dvaluenum(self,temperature, phase_pressure, variable="temperature"): + if variable == "temperature": + return np.gradient(self.value(temperature,phase_pressure),temperature) + elif variable == "phase_pressure": + return np.gradient(self.value(temperature,phase_pressure),phase_pressure) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self, variable="temperature"): + string = self._getcodeasstring("dvalue",variable) + string = self._convertpythontoexprtk(string) + return string + +class density_3(template.PROPERTY): + def value(self, temperature, phase_pressure): + # use only one-liner to keep it parsable: + return (5.38784232e-06*temperature**3 - 8.61928859e-03*temperature**2 + 3.44638154e+00*temperature + 5.92596450e+02) * np.exp(4.65e-10 * (np.maximum(phase_pressure, 0.0) - 1e5)) + def dvalue(self, temperature, phase_pressure, variable="temperature"): + if variable == "temperature": + return (3*5.38784232e-06*temperature**2 - 2*8.61928859e-03*temperature + 3.44638154e+00) * np.exp(4.65e-10 * (np.maximum(phase_pressure, 0.0) - 1e5)) + elif variable == "phase_pressure": + return (5.38784232e-06*temperature**3 - 8.61928859e-03*temperature**2 + 3.44638154e+00*temperature + 5.92596450e+02) * np.maximum(np.sign(phase_pressure),0) * 4.65e-10 + def dvaluenum(self,temperature, phase_pressure, variable="temperature"): + if variable == "temperature": + # use only one-liner to keep it parsable: + return np.gradient(self.value(temperature,phase_pressure),temperature) + elif variable == "phase_pressure": + # use only one-liner to keep it parsable: + return np.gradient(self.value(temperature,phase_pressure),phase_pressure) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self, variable="temperature"): + string = self._getcodeasstring("dvalue",variable) + string = self._convertpythontoexprtk(string) + return string diff --git a/Physics/water/water/properties/expansivity.py b/Physics/water/water/properties/expansivity.py new file mode 100644 index 000000000..ae511353f --- /dev/null +++ b/Physics/water/water/properties/expansivity.py @@ -0,0 +1,50 @@ +""" +Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" +# pylint: disable=C0103, R0902, R0914, R0913 +import numpy as np +from water.properties import template +T0 = 273.15 + +class expansivity_1(template.PROPERTY): + def value(self, temperature): + # source: https://www.engineeringtoolbox.com/water-density-specific-weight-d_595.html + temp = np.array([ + 0.0, 4., 10., 20., 30., 40., 50., 60., 70., 80., 90., 140., 200., + 260. + ]) + T0 + if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]): + print("The temperature is not within the defined domain") + beta = np.array([ + -5.e-5, 0.003e-4, 8.8e-5, 2.07e-4, 3.03e-4, 3.85e-4, 4.57e-4, 5.22e-4, + 5.82e-4, 6.40e-4, 6.95e-4, 9.75e-4, 1.59e-3, 2.21e-3 + ]) + return np.interp(temperature, temp, beta) + def dvalue(self, temperature): + return np.gradient(self.value(temperature), temperature) + def exprtk_value(self): + raise NotImplementedError + def exprtk_dvalue(self): + raise NotImplementedError + +class expansivity_2(template.PROPERTY): + def value(self, temperature): + # use only one-liner to keep it parsable: + return -1.46673616e-12 * temperature**4 + 2.46256925e-09 * temperature**3 - 1.51400613e-06 * temperature**2 + 4.11530308e-04 * temperature - 4.15284133e-02 + def dvalue(self, temperature): + # use only one-liner to keep it parsable: + return -1.46673616e-12 * 4 * temperature**3 + 2.46256925e-09 * 3 * temperature**2 - 1.51400613e-06 * 2 * temperature + 4.11530308e-04 + def dvaluenum(self,temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self): + string = self._getcodeasstring("dvalue",None) + string = self._convertpythontoexprtk(string) + return string diff --git a/Physics/water/water/properties/specificheatcapacity.py b/Physics/water/water/properties/specificheatcapacity.py new file mode 100644 index 000000000..ce4ed68d1 --- /dev/null +++ b/Physics/water/water/properties/specificheatcapacity.py @@ -0,0 +1,47 @@ +""" +Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" +# pylint: disable=C0103, R0902, R0914, R0913 +import numpy as np +from water.properties import template +T0 = 273.15 + +class specificheatcapacity_1(template.PROPERTY): + def value(self, temperature): + # source: https://www.engineeringtoolbox.com/specific-heat-capacity-water-d_660.html + temp = np.array([0.01, 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.]) + T0 + if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]): + print("The temperature is not within the defined domain") + cw = np.array([ + 4217., 4191., 4157., 4118., 4074., 4026., 3977., 3925., 3873., + 3820., 3768 + ]) + return np.interp(temperature, temp, cw) + def dvalue(self, temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + raise NotImplementedError + def exprtk_dvalue(self): + raise NotImplementedError + +class specificheatcapacity_2(template.PROPERTY): + def value(self, temperature): + # use only one-liner to keep it parsable: + return 1.55452794e-04 * temperature**3 - 1.64289282e-01 * temperature**2 + 5.25998084e+01 * temperature - 1.06049924e+03 + def dvalue(self, temperature): + # use only one-liner to keep it parsable: + return 1.55452794e-04 * 3 * temperature**2 - 1.64289282e-01 * 2 * temperature + 5.25998084e+01 + def dvaluenum(self,temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self): + string = self._getcodeasstring("dvalue",None) + string = self._convertpythontoexprtk(string) + return string diff --git a/Physics/water/water/properties/template.py b/Physics/water/water/properties/template.py new file mode 100644 index 000000000..8c9e81b00 --- /dev/null +++ b/Physics/water/water/properties/template.py @@ -0,0 +1,58 @@ +import pathlib +class PROPERTY: + def __init__(self): + pass + def value(self): + raise NotImplementedError + def dvalue(self): + raise NotImplementedError + def exprtk_value(self): + raise NotImplementedError + def exprtk_dvalue(self): + raise NotImplementedError + def _getcodeasstring(self, fct, variable): + classname = self.__class__.__name__ + string = "" + currentpath = pathlib.Path(__file__).parent.resolve() + with open(f"{currentpath}/{classname.split('_')[0]}.py") as source: + currentclass = None + currentmethod = None + currentvariable = None + for line in source: + linelist = line.split(" ") + try: + if linelist[0] == "class": + currentclass = linelist[1].split("(")[0] + currentmethod = None + currentvariable = None + except IndexError: + pass + try: + if linelist[4] == "def": + currentmethod = linelist[5].split("(")[0] + except IndexError: + pass + try: + if (linelist[8] == "if") or (linelist[8] == "elif"): + if linelist[9] == "variable": + currentvariable = linelist[11].split('"')[1] + except IndexError: + pass + try: + if classname == currentclass: + if currentmethod == fct: + if currentvariable == variable: + if "return" in linelist: + string = line.split("return")[1] + except IndexError: + pass + return string + + def _convertpythontoexprtk(self, string): + string = string.replace("**", "^") + string = string.replace("np.maximum","max") + string = string.replace("np.minimum","min") + string = string.replace("np.sign","sgn") + string = string.replace("np.","") + string = string.replace("\n","") + return string diff --git a/Physics/water/water/properties/viscosity.py b/Physics/water/water/properties/viscosity.py new file mode 100644 index 000000000..392d3c16f --- /dev/null +++ b/Physics/water/water/properties/viscosity.py @@ -0,0 +1,79 @@ +""" +Copyright (c) 2012-2021, OpenGeoSys Comunity (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" +# pylint: disable=C0103, R0902, R0914, R0913 +import numpy as np +from water.properties import template +T0 = 273.15 + +class viscosity_1(template.PROPERTY): + def value(self, temperature): + # source: https://www.engineeringtoolbox.com/water-dynamic-kinematic-viscosity-d_596.html + temp = np.array([10., 20., 30., 40., 50., 60., 70., 80., 90., 100.]) + T0 + if (np.min(temperature) < temp[0]) or (np.max(temperature) > temp[-1]): + print("The temperature is not within the defined domain") + vis = np.array([ + 0.0013, 0.001, 0.0007978, 0.0006531, 0.0005471, 0.0004658, + 0.0004044, 0.000355, 0.000315, 0.00002822 + ]) + return np.interp(temperature, temp, vis) + def dvalue(self, temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + raise NotImplementedError + def exprtk_dvalue(self): + raise NotImplementedError + +class viscosity_2(template.PROPERTY): + def value(self, temperature): + # use only one-liner to keep it parsable: + return -2.01570959e-09 * temperature**3 + 2.11402803e-06 * temperature**2 - 7.43932150e-04 * temperature + 8.82066680e-02 + def dvalue(self, temperature): + # use only one-liner to keep it parsable: + return -2.01570959e-09 * 3 * temperature**2 + 2.11402803e-06 * 2 * temperature - 7.43932150e-04 + def dvaluenum(self,temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self): + string = self._getcodeasstring("dvalue",None) + string = self._convertpythontoexprtk(string) + return string + +class viscosity_3(template.PROPERTY): + def value(self, temperature): + return 2.414*10**(-5)*10**(247.8/(temperature-140)) + def dvalue(self, temperature): + return -2.414*10**(-5)*247.8*np.log(10) * 10**(247.8/(temperature-140))/(temperature-140)**2 + def dvaluenum(self,temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self): + string = self._getcodeasstring("dvalue",None) + string = self._convertpythontoexprtk(string) + return string + +class viscosity_4(template.PROPERTY): + def value(self, temperature): + return 4.2844e-5+1/(0.157 * (temperature-208.157) * (temperature-208.157)-91.296) + def dvalue(self, temperature): + raise NotImplementedError + def dvaluenum(self,temperature): + return np.gradient(self.value(temperature),temperature) + def exprtk_value(self): + string = self._getcodeasstring("value", None) + string = self._convertpythontoexprtk(string) + return string + def exprtk_dvalue(self): + string = self._getcodeasstring("dvalue", None) + string = self._convertpythontoexprtk(string) + return string diff --git a/Physics/water/water/water.py b/Physics/water/water/water.py new file mode 100644 index 000000000..5f63d5617 --- /dev/null +++ b/Physics/water/water/water.py @@ -0,0 +1,32 @@ +""" +Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" + +# pylint: disable=C0103, R0902, R0914, R0913 +from water.properties import (density, viscosity, conductivity, specificheatcapacity, expansivity) + +class WATER: + """ + class defining water related properties at atmospheric + and arbitrary pressure + """ + def __init__(self, rho=1, mu=0, K=0, c=0, a=0): + self.densities = [density.density_1, + density.density_2, + density.density_3] + self.viscosities = [viscosity.viscosity_1, + viscosity.viscosity_2, + viscosity.viscosity_3, + viscosity.viscosity_4] + self.conductivities = [conductivity.conductivity_1, conductivity.conductivity_2] + self.specificheatcapacities = [specificheatcapacity.specificheatcapacity_1, specificheatcapacity.specificheatcapacity_2] + self.expansivities = [expansivity.expansivity_1, expansivity.expansivity_2] + self.rho = self.densities[rho]() + self.mu = self.viscosities[mu]() + self.K = self.conductivities[K]() + self.c = self.specificheatcapacities[c]() + self.a = self.expansivities[a]() -- GitLab