Fix superancillary implementation, update check notebook

402164ad · Ian Bell · 98072740 · 402164ad · 402164ad
Commit 402164ad authored 3 years ago by Ian Bell
--- a/include/teqp/models/cubics.hpp
+++ b/include/teqp/models/cubics.hpp
@@ -77,8 +77,8 @@ public:
        auto b = get_b(T, z);
        auto Ttilde = R(z)*T*b/get_a(T,z);
        return std::make_tuple(
-            CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOL_CODE, Ttilde)*b,
-            CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOV_CODE, Ttilde)*b
+            CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOL_CODE, Ttilde)/b,
+            CubicSuperAncillary::supercubic(superanc_index, CubicSuperAncillary::RHOV_CODE, Ttilde)/b
        );
    }


--- a/notebooks/CheckCubics.ipynb
+++ b/notebooks/CheckCubics.ipynb
@@ -37,7 +37,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
@@ -53,12 +53,13 @@
     "data": {
      "text/plain": [
       "{'Delta1': 2.414213562373095,\n",
+       " 'Delta2': -0.41421356237309515,\n",
       " 'OmegaA': 0.4572355289213822,\n",
       " 'OmegaB': 0.07779607390388846,\n",
       " 'kind': 'Peng-Robinson'}"
      ]
     },
-     "execution_count": 7,
+     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
@@ -76,7 +77,7 @@
  },
  {
   "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
@@ -92,12 +93,13 @@
     "data": {
      "text/plain": [
       "{'Delta1': 1.0,\n",
+       " 'Delta2': 0.0,\n",
       " 'OmegaA': 0.4274802335403417,\n",
       " 'OmegaB': 0.08664034996495766,\n",
       " 'kind': 'Soave-Redlich-Kwong'}"
      ]
     },
-     "execution_count": 6,
+     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
@@ -115,10 +117,36 @@
  },
  {
   "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
   "metadata": {},
-   "outputs": [],
-   "source": []
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "{'Delta1': 2.414213562373095, 'Delta2': -0.41421356237309515, 'OmegaA': 0.4572355289213822, 'OmegaB': 0.07779607390388846, 'kind': 'Peng-Robinson'}\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "([24221.746187758014, 1030.361868803444],\n",
+       " [24222.536725733073, 1029.779329419418])"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "model = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]])\n",
+    "print(model.get_meta())\n",
+    "[rhoL, rhoV ] = model.superanc_rhoLV(150)\n",
+    "\n",
+    "import CoolProp.CoolProp as CP\n",
+    "[rhoL, rhoV], CP.PropsSI('Dmolar','T',150,'Q',[0,1],'PR::METHANE').tolist()"
+   ]
  }
 ],
 "metadata": {

 %% Cell type:markdown id: tags:

 # Check Cubic EOS

 The values will not be precisely the same because Bell and Jager, JRN, 2016 used the conventional truncated values for $\Omega_a$ and $\Omega_b$ while teqp uses the exact values from Bell and Deiters, IECR, 2021

 %% Cell type:code id: tags:

 ``` python
 %matplotlib inline
 import matplotlib.pyplot as plt, numpy as np
 import matplotlib as mpl
 mpl.rcParams['figure.dpi']= 100
 import teqp
 ```

 %% Cell type:code id: tags:

 ``` python
 # Values taken from http://dx.doi.org/10.6028/jres.121.011
 Tc_K = [190.564, 154.581, 150.687]
 pc_Pa = [4599200, 5042800, 4863000]
 acentric = [0.011, 0.022, -0.002]
 z = np.array([0.5, 0.3, 0.2]) # mole fractions
 ```

 %% Cell type:code id: tags:

 ``` python
 model = teqp.canonical_PR(Tc_K, pc_Pa, acentric)

 T = 800 # K
 rho = 5000 # mol/m^3
 print(model.get_Ar01(T, rho, z), 9.6019116018396e-002)
 print(model.get_Ar02n(T, rho, z)[2], 2.3611667278971e-002)
 print(model.get_Ar03n(T, rho, z)[3], 1.7433108161805e-003)
 model.get_meta()
 ```

 %% Output

    0.0960116460031759 0.096019116018396
    0.02360946208541689 0.023611667278971
    0.0017427456394558435 0.0017433108161805

    {'Delta1': 2.414213562373095,
+     'Delta2': -0.41421356237309515,
     'OmegaA': 0.4572355289213822,
     'OmegaB': 0.07779607390388846,
     'kind': 'Peng-Robinson'}

 %% Cell type:code id: tags:

 ``` python
 model = teqp.canonical_SRK(Tc_K, pc_Pa, acentric)

 T = 800 # K
 rho = 5000 # mol/m^3
 print(model.get_Ar01(T, rho, z), 1.2741566551477e-001)
 print(model.get_Ar02n(T, rho, z)[2], 2.4895937945147e-002)
 print(model.get_Ar03n(T, rho, z)[3], 6.0986538256190e-003)
 model.get_meta()
 ```

 %% Output

    0.12741589269633935 0.12741566551477
    0.024896233252879204 0.024895937945147
    0.006098736292510815 0.006098653825619

    {'Delta1': 1.0,
+     'Delta2': 0.0,
     'OmegaA': 0.4274802335403417,
     'OmegaB': 0.08664034996495766,
     'kind': 'Soave-Redlich-Kwong'}

 %% Cell type:code id: tags:

 ``` python
+model = teqp.canonical_PR([Tc_K[0]], [pc_Pa[0]], [acentric[0]])
+print(model.get_meta())
+[rhoL, rhoV ] = model.superanc_rhoLV(150)
+
+import CoolProp.CoolProp as CP
+[rhoL, rhoV], CP.PropsSI('Dmolar','T',150,'Q',[0,1],'PR::METHANE').tolist()
 ```
+
+%% Output
+
+    {'Delta1': 2.414213562373095, 'Delta2': -0.41421356237309515, 'OmegaA': 0.4572355289213822, 'OmegaB': 0.07779607390388846, 'kind': 'Peng-Robinson'}
+
+    ([24221.746187758014, 1030.361868803444],
+     [24222.536725733073, 1029.779329419418])