Show how to get critical point at fixed composition

Closes #141
usnistgov · Jul 19, 2024 · 848fe3f · 848fe3f
1 parent 658be3c
commit 848fe3f
Showing 1 changed file with 58 additions and 0 deletions.
diff --git a/doc/source/algorithms/critical_curves.ipynb b/doc/source/algorithms/critical_curves.ipynb
@@ -292,6 +292,64 @@
     "plt.ylim(20, 300)\n",
     "plt.legend(loc='best');"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "933a88ef",
+   "metadata": {},
+   "source": [
+    "## Critical points\n",
+    "\n",
+    "If you have a relatively simple critical curve for a binary mixture and you would like to obtain the critical point at a given composition, you can obtain it reliably in a multi-step process:\n",
+    "\n",
+    "1. Trace the binary curve \n",
+    "2. Interpolate to find the approximate value for critical temperature and density\n",
+    "3. Use Newton's method to solve the criticality conditions from estimated critical point\n",
+    "4. Calculate whatever else you want"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4618df2a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = get_critical_curve([name0, othername], ipure)\n",
+    "model = teqp.build_multifluid_model([name0, othername], teqp.get_datapath())\n",
+    "line, = plt.plot(df['T / K'], df['p / Pa']/1e6, '-')\n",
+    "rhotot = df['rho0 / mol/m^3']+df['rho1 / mol/m^3']\n",
+    "\n",
+    "# Build cubic interpolators between the variables\n",
+    "tinterp = scipy.interpolate.interp1d(df['z0 / mole frac.'], df['t'], kind='cubic') # interpolator from mole fraction to tracing parameter\n",
+    "Tinterp = scipy.interpolate.interp1d(df['t'], df['T / K'], kind='cubic')\n",
+    "rhointerp = scipy.interpolate.interp1d(df['t'], rhotot, kind='cubic')\n",
+    "\n",
+    "# Iterate over composition\n",
+    "for x0 in np.arange(0.1, 1, 0.05):\n",
+    "\n",
+    "    molefracs = np.array([x0, 1-x0])\n",
+    "        \n",
+    "    # Initial guess for critical point from interpolation\n",
+    "    t = tinterp(x0)\n",
+    "    T0 = Tinterp(t)\n",
+    "    rho0 = rhointerp(t)\n",
+    "    p0 = rho0*model.get_R(molefracs)*T0*(1+model.get_Ar01(T0, rho0, molefracs))\n",
+    "    plt.plot(T0, p0/1e6, 'o')\n",
+    "    \n",
+    "    # Newton iteration for the correct critical point\n",
+    "    def resid(x):\n",
+    "        return model.get_criticality_conditions(T=x[0], rhovec=x[1]*molefracs)\n",
+    "    soln = scipy.optimize.fsolve(resid, x0=[T0, rho0])\n",
+    "    T, rho = soln\n",
+    "    p = rho*model.get_R(molefracs)*T*(1+model.get_Ar01(T, rho, molefracs))\n",
+    "    plt.plot(T, p/1e6, '+')\n",
+    "#     print(x0, T0, rho0, p0, resid([T0, rho0]), T, rho, resid([T, rho]))\n",
+    "\n",
+    "elap = timeit.default_timer()-tic\n",
+    "plt.gca().set(xlabel='$T$ / K', ylabel='$p$ / MPa')\n",
+    "plt.tight_layout(pad=0.2)"
+   ]
   }
  ],
  "metadata": {