From 33ff33dc3a8a54b9fe20e531c6849d27021ed139 Mon Sep 17 00:00:00 2001 From: tuopuu Date: Fri, 8 Jul 2016 16:31:11 +0300 Subject: [PATCH] Reduced equation of state tools. --- tools/maxwell_equal_area_rule.py | 258 ++++++++++++++++++++++++++ tools/reos.py | 304 +++++++++++++++++++++++++++++++ tools/tests/test_eos_tools.py | 146 +++++++++++++++ 3 files changed, 708 insertions(+) create mode 100644 tools/maxwell_equal_area_rule.py create mode 100644 tools/reos.py create mode 100644 tools/tests/test_eos_tools.py diff --git a/tools/maxwell_equal_area_rule.py b/tools/maxwell_equal_area_rule.py new file mode 100644 index 0000000..3c68c87 --- /dev/null +++ b/tools/maxwell_equal_area_rule.py @@ -0,0 +1,258 @@ +"""Computation of liquid-vapour densities. + +The densities are computed by utilizing the Maxwell's equal area rule +(aka Maxwell construction). +""" +import numpy as np +from scipy.integrate import quad +from scipy.optimize import brentq, fminbound +from scipy.signal import argrelmin, argrelmax + + +def _foff(arg1, arg2, f, C): + return f(arg1, arg2) - C + + +def _finters(a, b, f, arg2, C): + """Find an intersection of a 2 argument function f and a constant C. + + Parameters + ---------- + a : float + Lower limit for the search interval + b : float + Upper limit for the search interval + f : function of 2 arguments + arg2 : float + Second argument of f + C : float + Constant + + Returns + ------- + The intersection : float + + Raises + ------ + RuntimeError : If an intersection is not found. + + """ + # brentq from scipy + x0, r = brentq(_foff, a, b, (arg2, f, C), full_output=True) + + if not r.converged: + raise RuntimeError("An intersection not found: abort!") + + return x0 + + +def _area_calc(Pr, Tr, reosf, head_ab, body_ab, tail_ab): + """Compute the Maxwell's areas. + + That is, compute the areas demarcated + by the isotherm and a reduced pressure. + + Parameters + ---------- + Pr : float + Reduced pressure, aka constant pressure + Tr : float + Reduced temperature + reosf : function + Reduced equation of state + head_ab : iterable of float (at least two elements) + Head part of the isotherm (xstart and xend), aka head + body_ab : iterable of float (at least two elements) + Body part of the isotherm (xstart and xend), aka body + tail_ab : iterable of float (at least two elements) + Tail part of the isotherm (xstart and xend), aka tail + + Returns + ------- + A tuple with two elements: + 1. A tuple including the intersections (floats) between the constant + pressure and the head (1.1.), body (1.2.), and tail (1.3.) part. + + 2. A tuple including the area + 2.1. below the constant pressure (negative float), + 2.2. above the constant pressure (positive float). + + """ + h0 = _finters(head_ab[0], head_ab[1], reosf, Tr, Pr) + b0 = _finters(body_ab[0], body_ab[1], reosf, Tr, Pr) + t0 = _finters(tail_ab[0], tail_ab[1], reosf, Tr, Pr) + + # quad from scipy + a1 = quad(reosf, h0, b0, (Tr,))[0] - (b0 - h0)*Pr + a2 = quad(reosf, b0, t0, (Tr,))[0] - (t0 - b0)*Pr + + return ((h0, b0, t0), (a1, a2)) + + +def _area_diff(Pr, Tr, reosf, head_ab, body_ab, tail_ab): + """Compute the difference between Maxwell's areas. + + That is, compute the difference between areas + demarcated by the isotherm and a reduced pressure. + + Parameters + ---------- + Pr : float + Reduced pressure, aka constant pressure + Tr : float + Reduced temperature + reosf : function + Reduced equation of state + head_ab : iterable of float (at least two elements) + Head part of the isotherm (xstart and xend), aka head + body_ab : iterable of float (at least two elements) + Body part of the isotherm (xstart and xend), aka body + tail_ab : iterable of float (at least two elements) + Tail part of the isotherm (xstart and xend), aka tail + + Returns + ------- + The difference between Maxwell's areas (absolute value). + + """ + (intersects, areas) = _area_calc(Pr, Tr, reosf, head_ab, body_ab, tail_ab) + return np.fabs(areas[0] + areas[1]) + + +def liquid_vapour_density(Trs, reos, iSearchMax=40.0, RhorMin=1e-12): + """Compute liquid-vapour densities for the given reduced temperatures. + + The densities are computed by utilizing the Maxwell's equal area rule. + + Parameters + ---------- + Trs : iterable of float + Reduced temperatures + reos : ReducedEquationOfState + iSearchMax : float + Upper limit for the reduced molar volume, Vrmol, in the + search of the primary extrema (default value 40.0) + RhorMin : float + Lower limit for the reduced (vapour) density; limits the range + of solutions, i.e. a constraint for the equal area optimization + problem (default value 1e-12; VrmolMax := 1/RhorMin) + + Returns + ------- + List of tuples, where each tuple has 6 elements (floats): + 1. Reduced temperature, and associated + 2. Reduced equilibrium pressure, + 3. Reduced molar volume for the vapour phase, + 4. Reduced molar volume for the liquid phase, + 5. Reduced density for the vapour phase, and + 6. Reduced density for the liquid phase. + + Raises + ------ + ValueError : + If a reduced temperature parameter is >= 1. + RuntimeError : + If a primary minimum or maximum of an isotherm is not found. + RuntimeError : + If a reduced equilibrium pressure is not found (i.e. solution + to the equal area optimization problem is not converged). + + Algorithm + --------- + 1. Find the primary minimum and maximum of an isotherm. + - The two extrema are referred to as priMinPr and priMaxPr. + - Locations of the extrema are referred to as priMinVr and priMaxVr. + - The search of extrema locations is executed in two steps: + A. Delimit the locations to finite intervals; initial search + from the range ]VrmolMin, iSearchMax (parameter)], where + VrmolMin is defined by the reduced equation of state. + B. Pinpoint the locations in the intervals (refined search). + + 2. Decompose the isotherm into three parts. + - Head: ]VrmolMin, priMinVr]. + - Body: ]priMinVr, priMaxVr]. + - Tail: ]priMaxVr, VrmolMax := 1/RhorMin (parameter)]. + + 3. Find the reduced equilibrium pressure, PrEq, from the body part. + - Utilize the Maxwell's equal area rule (aka Maxwell construction). + - Search PrEq from the range + ]max(priMinPr, 0, Pr(VrmolMax)), priMaxPr[. + + 4. Find the liquid-vapour densities associated with PrEq. + - The intersection of PrEq and head defines the + reduced molar volume (and density) for the liquid phase. + - The intersection of PrEq and tail defines the + reduced molar volume (and density) for the vapour phase. + """ + # Exclude the lower limit for the reduced molar volume, as defined by + # the reduced equation of state, in order to avoid division by zero + iSearchRes = 1e-3 + VrmolMin, VrmolMax = reos.VrmolMin + iSearchRes, 1.0/RhorMin + Vrmols = np.arange(VrmolMin, iSearchMax, iSearchRes) + + vcnt = Vrmols.shape[0] + reosf = reos.Pr + results = [] + + for Tr in Trs: + if Tr >= 1.0: + raise ValueError('Reduced temperature paramater must be < 1.') + + Prs = np.fromfunction(lambda i: reosf(Vrmols[i], Tr), + (vcnt,), dtype=int) + + # 1.A Initial search (argrelmin and argrelmax from scipy) + priMinArr = argrelmin(Prs)[0] + priMaxArr = argrelmax(Prs)[0] + + if len(priMinArr) == 0 or len(priMaxArr) == 0: + msg = 'Initial search of a primary minimum or maximum failed.\n' + msg += 'A possible remedy: try to extend the search interval ' + msg += 'for finding primary extrema,\ni.e. increase the ' + msg += 'value of the iSearchMax parameter.' + raise RuntimeError(msg) + + mina, minb = Vrmols[priMinArr[0]-1], Vrmols[priMinArr[0]+1] + maxa, maxb = Vrmols[priMaxArr[0]-1], Vrmols[priMaxArr[0]+1] + + # 1.B refined search (fminbound from scipy) + priMinVr, priMinPr, e1, nf2 = fminbound(reosf, mina, minb, (Tr,), + xtol=1e-10, full_output=True) + priMaxVr, priMaxPr, e2, nf2 = fminbound(reosf, maxa, maxb, (Tr,), + xtol=1e-10, full_output=True) + + if e1 != 0 or e2 != 0: + msg = 'Refined search of a primary minimum or maximum failed ' + msg += 'due to unknown reason.' + raise RuntimeError(msg) + + # 2. Decompose the isotherm into three parts. + head_ab = [VrmolMin, priMinVr] + body_ab = [priMinVr, priMaxVr] + tail_ab = [priMaxVr, VrmolMax] + + # 3. Find the reduced equilibrium pressure from the isotherm body + PrVrmolMax = reosf(VrmolMax, Tr) + priMinPr = max(priMinPr, 0, PrVrmolMax) + + if PrVrmolMax < priMinPr: + tail_ab[1] = _finters(priMaxVr, VrmolMax, reosf, Tr, priMinPr) + + PrEq, fv, e3, nf3 = fminbound(_area_diff, priMinPr, priMaxPr, + (Tr, reosf, head_ab, body_ab, tail_ab), + xtol=1e-10, full_output=True) + + if e3 != 0: + msg = 'Reduced equilibrium pressure not found, i.e. solution\n' + msg += 'to the Maxwell\'s equal area optimization problem ' + msg += 'not converged.' + raise RuntimeError(msg) + + # 4. Find the liquid-vapour densities associated with PrEq. + VrmolLiq = _finters(head_ab[0], head_ab[1], reosf, Tr, PrEq) + VrmolVap = _finters(tail_ab[0], tail_ab[1], reosf, Tr, PrEq) + RhorVap, RhorLiq = 1.0/VrmolVap, 1.0/VrmolLiq + + results.append((Tr, PrEq, VrmolVap, VrmolLiq, RhorVap, RhorLiq)) + + return results diff --git a/tools/reos.py b/tools/reos.py new file mode 100644 index 0000000..d3bab6a --- /dev/null +++ b/tools/reos.py @@ -0,0 +1,304 @@ +"""Reduced equations of state for fluids. + +Reduced equations of state are general in the sense that they +do not involve material parameters (e.g. physico-chemical constants) +which characterize a particular, real fluid. +""" +import numpy as np +from abc import ABCMeta, abstractproperty, abstractmethod + + +class ReducedEquationOfState(object): + + """Abstract base class for the reduced equations of state.""" + + __metaclass__ = ABCMeta + + @abstractproperty + def Zc(self): + """Compressibility factor at the critical point.""" + + @abstractproperty + def VrmolMin(self): + """Lower limit for the reduced molar volume.""" + + @abstractmethod + def Pr(self, Vrmol, Tr): + """Equation of state for the reduced pressure. + + Parameters + ---------- + Vrmol : float + Reduced molar volume > VrmolMin; + Note, Vrmol = 1/reduced density + Tr : float + Reduced temperature > 0 + + Returns + ------- + Reduced pressure : float + """ + + +class ReducedIdealGas(ReducedEquationOfState): + + """Reduced ideal gas equation of state.""" + + @property + def Zc(self): + """Compressibility factor at the critical point.""" + return 1.0 + + @property + def VrmolMin(self): + """Lower limit for the reduced molar volume.""" + return 0.0 + + def Pr(self, Vrmol, Tr): + """Equation of state for the reduced pressure. + + Parameters + ---------- + Vrmol : float + Reduced molar volume > VrmolMin; + Note, Vrmol = 1/reduced density + Tr : float + Reduced temperature > 0 + + Returns + ------- + Reduced pressure : float + """ + return Tr/Vrmol + + +class ReducedVanDerWaals(ReducedEquationOfState): + + """Reduced Van der Waals equation of state.""" + + @property + def Zc(self): + """Compressibility factor at the critical point.""" + return 3.0/8.0 + + @property + def VrmolMin(self): + """Lower limit for the reduced molar volume.""" + return 1.0/3.0 + + def Pr(self, Vrmol, Tr): + """Equation of state for the reduced pressure. + + Parameters + ---------- + Vrmol : float + Reduced molar volume > VrmolMin; + Note, Vrmol = 1/reduced density + Tr : float + Reduced temperature > 0 + + Returns + ------- + Reduced pressure : float + """ + rep = 8*Tr/(3*Vrmol - 1) + att = 3/(Vrmol*Vrmol) + + return rep - att + + +class ReducedRedlichKwong(ReducedEquationOfState): + + """Reduced Redlich-Kwong equation of state.""" + + def __init__(self): + """Precompute coefficients of the equation of state.""" + self._cff = np.power(2, 1.0/3.0) - 1 + + @property + def Zc(self): + """Compressibility factor at the critical point.""" + return 1.0/3.0 + + @property + def VrmolMin(self): + """Lower limit for the reduced molar volume.""" + return self._cff + + def Pr(self, Vrmol, Tr): + """Equation of state for the reduced pressure. + + Parameters + ---------- + Vrmol : float + Reduced molar volume > VrmolMin; + Note, Vrmol = 1/reduced density + Tr : float + Reduced temperature > 0 + + Returns + ------- + Reduced pressure : float + """ + rep = 3*Tr/(Vrmol - self._cff) + att = 1/(self._cff*np.sqrt(Tr)*Vrmol*(Vrmol + self._cff)) + return rep - att + + +class ReducedSoave(ReducedEquationOfState): + + """Reduced Soave equation of state. + + Attributes + ---------- + Omega: float + Acentric factor (default value for water) + """ + + def __init__(self, Omega=0.344): + """Precompute coefficients of the equation of state.""" + self._Omega = Omega + self._cff = np.power(2, 1.0/3.0) - 1 + self._alpha_cff = 0.48508 + 1.55171*Omega - 0.15613*Omega*Omega + + @property + def Zc(self): + """Compressibility factor at the critical point.""" + return 1.0/3.0 + + @property + def Omega(self): + """Acentric factor""" + return self._Omega + + @property + def VrmolMin(self): + """Lower limit for the reduced molar volume.""" + return self._cff + + def Pr(self, Vrmol, Tr): + """Equation of state for the reduced pressure. + + Parameters + ---------- + Vrmol : float + Reduced molar volume > VrmolMin; + Note, Vrmol = 1/reduced density + Tr : float + Reduced temperature > 0 + + Returns + ------- + Reduced pressure : float + """ + rep = 3*Tr/(Vrmol - self._cff) + + alpha1 = 1 + self._alpha_cff*(1 - np.sqrt(Tr)) + alpha = alpha1*alpha1 + + att = alpha/(self._cff*Vrmol*(Vrmol + self._cff)) + return rep - att + + +class ReducedPengRobinson(ReducedEquationOfState): + + """Reduced Peng-Robinson equation of state. + + Attributes + ---------- + Omega : float + Acentric factor (default value for water) + """ + + def __init__(self, Omega=0.344): + """Precompute coefficients of the equation of state.""" + self._Omega = Omega + self._cff = 0.2534 + self._alpha_cff = 0.37464 + 1.54226*Omega - 0.26992*Omega*Omega + + @property + def Zc(self): + """Compressibility factor at the critical point.""" + return 0.307 + + @property + def Omega(self): + """Acentric factor""" + return self._Omega + + @property + def VrmolMin(self): + """Lower limit for the reduced molar volume.""" + return self._cff + + def Pr(self, Vrmol, Tr): + """Equation of state for the reduced pressure. + + Parameters + ---------- + Vrmol : float + Reduced molar volume > VrmolMin; + Note, Vrmol = 1/reduced density + Tr : float + Reduced temperature > 0 + + Returns + ------- + Reduced pressure : float + """ + rep = 3.2573*Tr/(Vrmol - self._cff) + + alpha1 = 1 + self._alpha_cff*(1 - np.sqrt(Tr)) + alpha = 4.8514*alpha1*alpha1 + + att = alpha/(Vrmol*Vrmol + 2*self._cff*Vrmol - self._cff*self._cff) + return rep - att + + +class ReducedCarnahanStarling(ReducedEquationOfState): + + """Reduced Carnahan-Starling equation of state.""" + + def __init__(self): + """Precompute coefficients of the equation of state.""" + self._cff1 = 2.785855166 + self._cff2 = 0.1304438842 + self._cff3 = 3.852462257 + + @property + def Zc(self): + """Compressibility factor at the critical point.""" + return 0.3589562 + + @property + def VrmolMin(self): + """Lower limit for the reduced molar volume.""" + return self._cff2 + + def Pr(self, Vrmol, Tr): + """Equation of state for the reduced pressure. + + Parameters + ---------- + Vrmol : float + Reduced molar volume > VrmolMin; + Note, Vrmol = 1/reduced density + Tr : float + Reduced temperature > 0 + + Returns + ------- + Reduced pressure : float + """ + aux1 = self._cff2/Vrmol + aux2 = aux1*aux1 + + aux3 = 1 - aux1 + aux4 = aux3*aux3*aux3 + + rep1 = self._cff1*Tr/Vrmol + rep2 = (1 + aux1 + aux2 - aux1*aux2)/aux4 + rep = rep1*rep2 + + att = self._cff3/(Vrmol*Vrmol) + return rep - att diff --git a/tools/tests/test_eos_tools.py b/tools/tests/test_eos_tools.py new file mode 100644 index 0000000..a23f1ca --- /dev/null +++ b/tools/tests/test_eos_tools.py @@ -0,0 +1,146 @@ +"""Test reduced equation of state tools.""" +import unittest +import numpy as np +from simphony.tools import reos +from simphony.tools.maxwell_equal_area_rule import ( + liquid_vapour_density) + + +class MaxwellTestCase(unittest.TestCase): + + """Test application of the Maxwell's equal area rule.""" + + def test_maxwell(self): + """Density computation with several reduced equations of state.""" + ig = reos.ReducedIdealGas() + vdw = reos.ReducedVanDerWaals() + rk = reos.ReducedRedlichKwong() + so = reos.ReducedSoave() + pr = reos.ReducedPengRobinson() + cs = reos.ReducedCarnahanStarling() + + # Note, sorted + Trs = [0.35, 0.55, 0.75, 0.95] + Trs_len = len(Trs) + + vdw_results = liquid_vapour_density(Trs, vdw) + rk_results = liquid_vapour_density(Trs, rk) + so_results = liquid_vapour_density(Trs, so) + pr_results = liquid_vapour_density(Trs, pr) + cs_results = liquid_vapour_density(Trs, cs) + + self.assertEqual(Trs_len, len(vdw_results)) + self.assertEqual(Trs_len, len(rk_results)) + self.assertEqual(Trs_len, len(so_results)) + self.assertEqual(Trs_len, len(pr_results)) + self.assertEqual(Trs_len, len(cs_results)) + + for r in np.arange(0, Trs_len-1): + self.assertEqual(Trs[r], vdw_results[r][0]) + self.assertEqual(Trs[r], rk_results[r][0]) + self.assertEqual(Trs[r], so_results[r][0]) + self.assertEqual(Trs[r], pr_results[r][0]) + self.assertEqual(Trs[r], cs_results[r][0]) + + # Reduced molar volume must be + # greater for the vapour phase than for the liquiq phase + self.assertGreater(vdw_results[r][2], vdw_results[r][3]) + self.assertGreater(rk_results[r][2], rk_results[r][3]) + self.assertGreater(so_results[r][2], so_results[r][3]) + self.assertGreater(pr_results[r][2], pr_results[r][3]) + self.assertGreater(cs_results[r][2], cs_results[r][3]) + + # Reduced density must be + # less for the vapour phase than for the liquiq phase + self.assertLess(vdw_results[r][4], vdw_results[r][5]) + self.assertLess(rk_results[r][4], rk_results[r][5]) + self.assertLess(so_results[r][4], so_results[r][5]) + self.assertLess(pr_results[r][4], pr_results[r][5]) + self.assertLess(cs_results[r][4], cs_results[r][5]) + + # Reduced equilibrium pressure must be smaller + # for lower reduced temperature (note that Trs must be sorted) + self.assertLess(vdw_results[r][1], vdw_results[r+1][1]) + self.assertLess(rk_results[r][1], rk_results[r+1][1]) + self.assertLess(so_results[r][1], so_results[r+1][1]) + self.assertLess(pr_results[r][1], pr_results[r+1][1]) + self.assertLess(cs_results[r][1], cs_results[r+1][1]) + + # Reduced molar volume for the vapour phase must be greater + # for lower reduced temperature (note that Trs must be sorted) + self.assertGreater(vdw_results[r][2], vdw_results[r+1][2]) + self.assertGreater(rk_results[r][2], rk_results[r+1][2]) + self.assertGreater(so_results[r][2], so_results[r+1][2]) + self.assertGreater(pr_results[r][2], pr_results[r+1][2]) + self.assertGreater(cs_results[r][2], cs_results[r+1][2]) + + # Reduced molar volume for the liquid phase must be smaller + # for lower reduced temperature (note that Trs must be sorted) + self.assertLess(vdw_results[r][3], vdw_results[r+1][3]) + self.assertLess(rk_results[r][3], rk_results[r+1][3]) + self.assertLess(so_results[r][3], so_results[r+1][3]) + self.assertLess(pr_results[r][3], pr_results[r+1][3]) + self.assertLess(cs_results[r][3], cs_results[r+1][3]) + + # Reduced density for the vapour phase must be smaller + # for lower reduced temperature (note that Trs must be sorted) + self.assertLess(vdw_results[r][4], vdw_results[r+1][4]) + self.assertLess(rk_results[r][4], rk_results[r+1][4]) + self.assertLess(so_results[r][4], so_results[r+1][4]) + self.assertLess(pr_results[r][4], pr_results[r+1][4]) + self.assertLess(cs_results[r][4], cs_results[r+1][4]) + + # Reduced density for the liquiq phase must be greater + # for lower reduced temperature (note that Trs must be sorted) + self.assertGreater(vdw_results[r][5], vdw_results[r+1][5]) + self.assertGreater(rk_results[r][5], rk_results[r+1][5]) + self.assertGreater(so_results[r][5], so_results[r+1][5]) + self.assertGreater(pr_results[r][5], pr_results[r+1][5]) + self.assertGreater(cs_results[r][5], cs_results[r+1][5]) + + try: + liquid_vapour_density(Trs, ig) + except RuntimeError: + pass + else: + msg = 'Computation of densities must fail with Ideal-gas reos.' + raise AssertionError(msg) + + try: + liquid_vapour_density([1.01], vdw) + except ValueError: + pass + else: + raise AssertionError('Tr > 1 must fail with VdW reos.') + + try: + liquid_vapour_density([1.01], rk) + except ValueError: + pass + else: + raise AssertionError('Tr > 1 must fail with R-K reos.') + + try: + liquid_vapour_density([1.01], so) + except ValueError: + pass + else: + raise AssertionError('Tr > 1 must fail with Soave reos.') + + try: + liquid_vapour_density([1.01], pr) + except ValueError: + pass + else: + raise AssertionError('Tr > 1 must fail with P-R reos.') + + try: + liquid_vapour_density([1.01], cs) + except ValueError: + pass + else: + raise AssertionError('Tr > 1 must fail with C-S reos.') + + +if __name__ == '__main__': + unittest.main()