Convexit_concaviteSU.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ATTENTION, work in progress! Ce notebook n'est que partiellement fonctionnel!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tracé des surfaces S(U,V) et U(S,V) pour GP et vdW\n",
    "\n",
    "On a démontré que $S$ est une fonction concave de $U$ et $V$ et que U est une fonction convexe de $S$ et $V$. \n",
    "\n",
    "On utilise les équations d'état du gaz parfait et de van der Waals pour tracer ces surfaces. Pour van der Waals, on voit les zones d'instabilité qui correspondent aux domaines biphasés. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equation d'état de Van der Waals\n",
    "\n",
    "L'équation d'état de Van der Waals est donnée par :\n",
    "\n",
    "$$(p+a\\frac{n^2}{V^2})(V-bn)=nRT$$\n",
    "\n",
    "$a$ et $b$ sont des paramètres du gaz. \n",
    "\n",
    "Au point critique, le gaz est à $T_c$ et $p_c$, et une mole occupe le volume critique molaire $V_c$\n",
    "\n",
    "$$V_c=3b$$\n",
    "\n",
    "$$p_c=\\frac{a}{27b^2}$$\n",
    "\n",
    "$$T_c=\\frac{8a}{27 R b}$$\n",
    "\n",
    "Pour une mole de fluide, l'équation d'état devient:\n",
    "$$p=\\frac{RT}{V-b}-\\frac{a}{V^2}$$\n",
    "\n",
    "En coordonnées normalisées aux valeurs critiques, $p^*=p/p_c$; $T^*=T/T_c$ et $V^*=V/V_c$ et pour une mole\n",
    "\n",
    "$$p^*=\\frac{8T^*}{3V^*-1}-\\frac{3}{V^{*2}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import csv\n",
    "from IPython.display import Image\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "import bokeh\n",
    "from bokeh.io import push_notebook, show, output_notebook\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.models.glyphs import Patch\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Constante des gaz parfaite\n",
    "R=8.314\n",
    "\n",
    "# paramètres a et b pour certains gaz\n",
    "# gaz=(a,b) a en l^2.bar.mol^(-2) et b en l.mol^(-1)\n",
    "\n",
    "co2=(3.64,0.04267)\n",
    "sf6=(7.857,0.08786)\n",
    "h2o=(5.537,0.03049)\n",
    "\n",
    "\n",
    "# valcrit(gaz) =(Tc,Pc,Vc) en unités SI et pour une mole\n",
    "# Calcule les valeurs critiques grâce à a et b\n",
    "\n",
    "def valcrit(gaz):\n",
    "    valcrit=(100*8*gaz[0]/(27*R*gaz[1]),10**5*gaz[0]/(27*gaz[1]**2),3*gaz[1]/10**3)\n",
    "    return valcrit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tracé des surfaces S(U,V) et U(S,V) pour GP et vdW"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pour un gaz parfait \n",
    "$$dU=TdS-pdV=cNRdT$$. \n",
    "$$dS=cNR\\frac{dT}{T}+\\frac{p}{T}dV =\n",
    "cNR\\frac{dT}{T}+NR\\frac{dV}{V}$$ \n",
    "\n",
    "$$S-S_0 =NR\\ln\\Big[\\Big(\\frac{U}{U_0}\\Big)^c\\cdot\\Big(\\frac{V}{V_0}\\Big)\\Big]$$\n",
    "\n",
    "Cette forme est la même que l'équation de Sakur-Tetrode dérivée avec des considérations de phy stat sur un modèle de GP pour un gaz monoatomique (c=3/2). Elle ne fonctionne pas à très basse température. \n",
    "\n",
    "$$\n",
    "\\frac{S}{k_{\\rm B} N} = \\ln\n",
    "\\left[ \\frac VN \\left(\\frac{4\\pi m}{3h^2}\\frac UN\\right)^{3/2}\\right]+\n",
    "{\\frac 52}\n",
    "$$\n",
    "\n",
    "Et donc aussi\n",
    "\n",
    "$$U =U_0\\Big(\\frac{V}{V_0}\\Big)^{-1/c}\\exp\\big({\\frac{S-S_0}{cNR}}\\big)$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# définition de S(U,V) pour un gaz parfait\n",
    "def suvgp(U,V):\n",
    "    c=3/2 #GP monoatomique\n",
    "    S0=0\n",
    "    N=1\n",
    "    R=8.31\n",
    "    fonc = S0+np.log(U**c*V)\n",
    "    return fonc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# définition de U(S,V) pour un gaz parfait\n",
    "def usvgp(S,V):\n",
    "    c=3/2 #GP monoatomique\n",
    "    S0=0\n",
    "    N=1\n",
    "    R=8.31\n",
    "    en_int_rel = V**(-1/c) * np.exp((S-S0)/(c*N*R))\n",
    "    return en_int_rel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Utab = np.linspace(10,80,101)\n",
    "Vtab1 = np.linspace(2,101,91)\n",
    "Vtab = [round(x, 2) for x in Vtab1]\n",
    "SUVGPtab = []\n",
    "    \n",
    "for i in Utab:\n",
    "    SUVGPtab.append(suvgp(i,Vtab).tolist())\n",
    "    \n",
    "df_suvgp = pd.DataFrame(SUVGPtab,columns=Vtab, index=Utab)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Stab1 = np.linspace(70,80,101)\n",
    "Stab = [round(x, 2) for x in Stab1]\n",
    "Vtab = np.linspace(2,11,91)\n",
    "USVGPtab = []\n",
    "for i in Stab:\n",
    "    USVGPtab.append(usvgp(i,Vtab).tolist())\n",
    "    \n",
    "df_usvgp = pd.DataFrame(USVGPtab,columns=Vtab, index=Stab)      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "def plottable_3d_info(df_usvgp: pd.DataFrame):\n",
    "    \"\"\"\n",
    "    Transform Pandas data into a format that's compatible with\n",
    "    Matplotlib's surface and wireframe plotting.\n",
    "    \"\"\"\n",
    "    index = df_usvgp.index\n",
    "    columns = df_usvgp.columns\n",
    "\n",
    "    x, y = np.meshgrid(np.arange(len(columns)), np.arange(len(index)))\n",
    "    z = np.array([[df_usvgp[c][i] for c in columns] for i in index])\n",
    "    \n",
    "    xticks = dict(ticks=np.linspace(0, len(columns) - 1, 6), labels=columns[::len(columns) // 5])\n",
    "    yticks = dict(ticks=np.linspace(0, len(index) - 1, 6), labels=index[::len(index) // 5])\n",
    "    \n",
    "    return x, y, z, xticks, yticks\n",
    "\n",
    "\n",
    "### Transform to Matplotlib friendly format.\n",
    "x, y, z, xticks, yticks = plottable_3d_info(df_usvgp)\n",
    "\n",
    "### Set up axes and put data on the surface.\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(projection='3d')\n",
    "axes.plot_surface(x, y, z, cmap='gist_earth')\n",
    "\n",
    "### Customize labels and ticks (only really necessary with\n",
    "### non-numeric axes).\n",
    "axes.set_xlabel('S (unités arbitraires)')\n",
    "axes.set_ylabel('V (unités arbitraires)')\n",
    "axes.set_zlabel('U (unités arbitraires)')\n",
    "axes.set_zlim3d(bottom=0)\n",
    "axes.view_init(elev=20, azim=-77)\n",
    "plt.xticks(**xticks)\n",
    "plt.yticks(**yticks)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "def plottable_3d_info(df_suvgp: pd.DataFrame):\n",
    "    \"\"\"\n",
    "    Transform Pandas data into a format that's compatible with\n",
    "    Matplotlib's surface and wireframe plotting.\n",
    "    \"\"\"\n",
    "    index = df_suvgp.index\n",
    "    columns = df_suvgp.columns\n",
    "\n",
    "    x, y = np.meshgrid(np.arange(len(columns)), np.arange(len(index)))\n",
    "    z = np.array([[df_suvgp[c][i] for c in columns] for i in index])\n",
    "    \n",
    "    xticks = dict(ticks=np.linspace(0, len(columns) - 1, 7), labels=columns[::len(columns) // 6])\n",
    "    yticks = dict(ticks=np.linspace(0, len(index) - 1, 5), labels=index[::len(index) // 4])\n",
    "    \n",
    "    return x, y, z, xticks, yticks\n",
    "\n",
    "\n",
    "### Transform to Matplotlib friendly format.\n",
    "x, y, z, xticks, yticks = plottable_3d_info(df_suvgp)\n",
    "\n",
    "### Set up axes and put data on the surface.\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(projection='3d')\n",
    "axes.plot_surface(x, y, z, cmap='gist_earth')\n",
    "\n",
    "### Customize labels and ticks (only really necessary with\n",
    "### non-numeric axes).\n",
    "axes.set_xlabel('U (unités arbitraires)')\n",
    "axes.set_ylabel('V (unités arbitraires)')\n",
    "axes.set_zlabel('S (unités arbitraires)')\n",
    "axes.set_zlim3d(bottom=3)\n",
    "axes.view_init(elev=21, azim=-153)\n",
    "plt.xticks(**xticks)\n",
    "plt.yticks(**yticks)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notebook by Cécile Hébert (2023).\n",
    "Except where otherwise noted, the content of this notebook is licensed under MIT licence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}