Case Study 1 - solving nonlinear equations (example van der Waals equation)

Case 1 - van der Waals v1_1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Constants  \n",
    "R is gas constant (atm.litre/g-mol.K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 0.08206"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Critical parameters for ammonia  \n",
    "Tc is critical temperature (K)   \n",
    "Pc is critical pressure (atm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tc = 405.5\n",
    "Pc = 111.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### van der Waals equation of state:  \n",
    "    V is molar volume (liters/g-mol)  \n",
    "    P is pressure (atm)  \n",
    "    T is temperature (K)   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def vdw_calc(V):\n",
    "    a = 27.0 / 64.0 *(R**2 * Tc**2 / Pc)\n",
    "    b = R*Tc / (8.0 * Pc)\n",
    "    converg_error = (P + a/V**2)*(V-b) - R*T\n",
    "    return converg_error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Compresibility factor  \n",
    "Z is compressibility factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def comp_factor():\n",
    "    Z = P*V / (R*T)\n",
    "    return Z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Problem part a)\n",
    "Calculate the molar volume and compressibility factor for gaseous ammonia at a pressure P = 56 atm and a temperature T = 450 K using the van der Waals equation of state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "molar volume = 0.5749 liters/g-mol\n",
      "compressibility factor = 0.8718\n"
     ]
    }
   ],
   "source": [
    "T = 450.0 \n",
    "P = 56.0\n",
    "\n",
    "V = fsolve(vdw_calc, 1.0)[0]\n",
    "Z = comp_factor()\n",
    "\n",
    "print('molar volume = %.4f liters/g-mol' %V)\n",
    "print('compressibility factor = %.4f' %Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Problem part b)\n",
    "Repeat the calculations for the following reduced pressures: Pr = 1, 2, 4, 10, and 20. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "reduced pressure = 1.0\n",
      "  molar volume = 0.2335 liters/g-mol\n",
      "  compressibility factor = 0.7038\n",
      "\n",
      "reduced pressure = 2.0\n",
      "  molar volume = 0.0773 liters/g-mol\n",
      "  compressibility factor = 0.4658\n",
      "\n",
      "reduced pressure = 4.0\n",
      "  molar volume = 0.0607 liters/g-mol\n",
      "  compressibility factor = 0.7313\n",
      "\n",
      "reduced pressure = 10.0\n",
      "  molar volume = 0.0509 liters/g-mol\n",
      "  compressibility factor = 1.5334\n",
      "\n",
      "reduced pressure = 20.0\n",
      "  molar volume = 0.0462 liters/g-mol\n",
      "  compressibility factor = 2.7835\n",
      "\n"
     ]
    }
   ],
   "source": [
    "T = 450.0 \n",
    "Pr_list = [1,2,4,10,20]\n",
    "\n",
    "for Pr in Pr_list:\n",
    "    P = Pr * Pc\n",
    "    V = fsolve(vdw_calc, 1.0)[0]\n",
    "    Z = comp_factor()\n",
    "    \n",
    "    print('reduced pressure = %.1f' %Pr)\n",
    "    print('  molar volume = %.4f liters/g-mol' %V)\n",
    "    print('  compressibility factor = %.4f' %Z)\n",
    "    print('')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Problem part c)  \n",
    "How does the compressibility factor vary as a function of Pr?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4nHd57//3rc22rN1abcvxGsdbYltOnIQQYhKyGJJAEiCcwClwTlNOoZwutD96SkvL9euvFz0051AoO6GkUJxCoGFJYhvHTkJWLzixZce7402rZe3WMpr798c8FhNFy1jWLJI+r+vSpZl5nmfmo0ejufVs99fcHREREYC0ZAcQEZHUoaIgIiL9VBRERKSfioKIiPRTURARkX4qCiIi0k9FQURE+qkoiIhIPxUFERHpl5HsABeruLjY586de9HLdXR0MH369LEPdIlSNRekbrZUzQWpm025Ll6qZhttrp07dza6e8mIM7r7uPqqqqry0di6deuolou3VM3lnrrZUjWXe+pmU66Ll6rZRpsL2OExfMZq95GIiPRTURARkX4qCiIi0k9FQURE+qkoiIhIPxUFERHpp6IgIiL9VBRERFJcOOy8cqyJutauuL/WuLuiWURkMmnq6GFjdS21LV30hIri/noqCiIiKcjd2XWimRcON5KRnsb6FRUsLs9l2+n4vq6KgohIimnu7GHTvjpOnzvP/JLp3LKkjOlTEvNxraIgIpIi3J3XTrXwm8ONANy6rIylFXmYWcIyqCiIiKSA1q5eNlfXcaKpk8tmZHPL0jLypmYmPIeKgohIErk7+2pa2XagAYCbl5SyYlZ+QrcOoqkoiIgkSXt3iC376zja0MGswmnctrSc/OzEbx1EU1EQEUkwd+dgXTtPv15PqC/MOxaXsKqyIGlbB9FUFEREEqizJ8TTr9dzqK6divyp3LqsnKLpWcmO1U9FQUQkQQ7Xt7Nlfx3doTA3LCqmak4haWnJ3zqIpqIgIhJnXb19bDtQz/6aNkpyp3DP6nJKcqckO9agVBREROLoeGMHm/fV0dnTx9r5RaydN4P0FNs6iKaiICISB92hPp472Mie0y3MyMnirpUzKcubmuxYI1JREBEZYyebOtm0r462rl7WzC3kuvkzyEgfH02pVRRERMZIb1+Y3xxuZPeJZgqyM/nAmkpmFkxLdqyLoqIgIjIGzjSfZ1N1Lec6e1k5p4C3LSgmK2N8bB1EU1EQEbkEob4wLx49y843zpE7NZP7qmZTWZSd7FijpqIgIjJKda1dbKyu5Wx7Dytm5fP2y4uZkpGe7FiXREVBROQi9YWdl4+dZfuxc2RnpfPeVbOYVzw92bHGhIqCiMhFaGzvZmN1LfWt3SypyOOmxSVMzRzfWwfR4lYUzKwSeAQoB8LAt9z9ywPmuQl4HDgWPPRTd/9CvDKJiIxWOOzsOnGOF46cZUpGGndeVcHC0txkxxpz8dxSCAF/5u67zCwX2Glmm91934D5nnP398Qxh4jIJWnu7GFTdR2nm8+zsDSHm5eUkp01MXe0xO2ncvcaoCa43WZm+4FZwMCiICKSki4Mj/ncoQbS0ozbl5dzRXluSrS4jhdz9/i/iNlc4Flgubu3Rj1+E/AYcAo4A3zG3asHWf5B4EGAsrKyqg0bNlx0hvb2dnJyckaRPr5SNRekbrZUzQWpm025Ll5DczsHO6ZQ3+mUZhuryzKYlpH8YjDadbZu3bqd7r5mxBndPa5fQA6wE7hnkGl5QE5wez1waKTnq6qq8tHYunXrqJaLt1TN5Z662VI1l3vqZlOu2IXDYd93psX/7Lsb/StbDvruE+c8HA4nO1a/0a4zYIfH8Jkd151iZpZJZEvgh+7+00EKUmvU7SfM7GtmVuzujfHMJSIymM6eEFv213O4vp3cLOPD115GQXbqDICTCPE8+8iA7wL73f2hIeYpB+rc3c3sGiANOBuvTCIiQ4keAOfti4ppzTg96QoCxPfso7cBHwH2mNnu4LH/BcwBcPdvAPcB/8PMQsB54P5gM0dEJCG6evt45mAD+860vmkAnG3Hk3/8IBniefbRb4Bh16q7fxX4arwyiIgM52RTJxura2nvDrF2XhFr56f2ADiJMDFPtBURGUZ0i+vC7Ew+eHUlFfnjq8V1vKgoiMikUtNyno17f9fi+oaFxWSOkwFwEkFFQUQmhb6w89LRs2w/3kTOlIxx3+I6XlQURGTCa2iLNLFraOtm2cw8brx8YjWxG0sqCiIyYYXDzs4T53jxyFmmZqZx18qZLChJzSuoU4WKgohMSOc6eti0r5YzzV0sKsvh5ivKmJalrYORqCiIyITi7rx6qoXfBE3s7lhRzuKyid3EbiypKIjIhNHa1cvm6jpONHUytzibW5aUkTs1M9mxxhUVBREZ99yd/TVtbDtYjzvcvKSUFbPytXUwCioKIjKudfaE+PX+eo7UtzOrYBq3LiublD2LxoqKgoiMW4fr29iyv57uUJgbLy9mVWUhaZO8TcWlUlEQkXGnq7ePbQca2F/TSmneFO5dVk5xzpRkx5oQVBREZFx542wHm/fV0dHdx9r5RaydpyZ2Y0lFQUTGhZ5QmOcPN7L7ZDNF07P44NUzKc+fmuxYE46KgoikvDPN59lYXUvL+V5WX1bI9QtmqIldnAxbFMwsHfi0u/+fBOUREekX3cQud2om965WE7t4G7YouHufmd0NqCiISEI1tkea2NW3RprYvWNxCVMy1KYi3mLZffS8mX0VeBTouPCgu++KWyoRmbTcnV0nmnnhcCNZGWncedVMFpaqiV2ixFIUrg++fyHqMQfeOfZxRGQyaznfy6bqWk6dO8+C0hxuWVJKdpYOfSbSiGvb3dclIoiITF7uzr6aVrYdaADgXUvLWDYzT20qkmDEomBm+cDngRuDh54BvuDuLfEMJiKTQ2dPiC376zlc386swmnctqyc/GlqYpcssWyXPQzsBT4Q3P8I8D3gnniFEpHJobYjzA9eeoOuXrWpSBWxFIUF7n5v1P2/M7Pd8QokIhNfd6iPZw828uKZEFevyOB9q8opyVWbilQQS1E4b2Y3uPtvAMzsbcD5+MYSkYnqdPN5Nu6tpbWrl8sL07n/6koydCFayoilKHwCeCQ4tgBwDvi9+EUSkYko1BfmpaNN7Hijibypmbx/TSWHdp9RQUgxsRSFVne/yszyANy91czmxTmXiEwgDW3dPFVdS2NbNytm5fP2y4uZkpHOoWQHk7eIpSg8Bqx299aox34CVMUnkohMFOGws+vEOV44cpapmWncvXIm80t0IVoqG7IomNkVwDIg38yizzTKA9SaUESG1dLZy8Z9tZw+d56FpTncrAvRxoXhfkOLgfcABcCdUY+3Ab8fz1AiMn65O9VnWnnmYORCtNuWlbOkIlcXoo0TQxYFd38ceNzMrnP3FxOYSUTGqejxkmcXTuNWXYg27sR09pGZ7Xf3ZgAzKwT+yd0/Ht9oIjKeHK5vZ8v+OnpCYW68vITVcwq0dTAOxVIUrrxQEADc/ZyZrYpjJhEZR7pDfTxzoIHqM5Hxkm/TeMnjWixFIc3MCt39HICZFcW4nIhMcKfOdbKxuo62rl7Wziti7XyNlzzexfLh/k/AC2b2k+D++4G/H2khM6sEHgHKgTDwLXf/8oB5DPgysB7oBD6qcRpEUl+oL8wLR86y68Q58qdl8oE1lcwsmJbsWDIGYmmd/YiZ7QTWAQbc4+77YnjuEPBn7r7LzHKBnWa2ecCydwCLgq+1wNeD7yKSourbuti4t5bG9h6unJ3P2xeVkJWhq5Iniph2A7l7tZk1EFyfYGZz3P3ECMvUADXB7TYz2w/MAqKLwt3AI+7uwEtmVmBmFcGyIpJCwmFn54lzvBhciPbeVbOYVzw92bFkjMUynsJdRHYhzQTqgcuA/UQubIuJmc0FVgEvD5g0CzgZdf9U8JiKgkgKaensZWN1Laebz3N5WS7vvKKUaVkaL3kissg/6cPMYPYqkaE3f+3uq8xsHfAhd38wphcwyyEyMM/fu/tPB0z7FfAPUR1YtwB/4e47B8z3IPAgQFlZWdWGDRti+uGitbe3k5OTepfXp2ouSN1sqZoLUjfbaHO5O2+0htnT2IcZXFWSwewcG7NTTVN1fUHqZhttrnXr1u109zUjzujuw34BO4LvrwJpwe1XRloumC8T2Aj86RDTv0mkwFy4fwCoGO45q6qqfDS2bt06quXiLVVzuadutlTN5Z662UaTq6O71x/ffdof2nTAf7zjpLec70mJXImSqtlGm+vCZ/lIX7EcU2gO/tt/FvihmdUTOYg8rODMou8C+939oSFm+znwKTPbQOQAc4vreIJI0h1r7GDzvtpgRDRdiDaZxFIU7ga6gD8BHgDygS/EsNzbiAzduSdqpLb/BcwBcPdvAE8QOR31MJFTUj92MeFFZGz1hMI8d6iB1061UJw7RSOiTULDdUm91t1fcveOqIe/H+sTe+Q4wbD/WgSbNJ+M9TlFJH5qW7p4am8Nzed7qbqskOsXzNAAOJPQcFsKXwNWA5jZi+5+XWIiiUgihcPOK8ebePloE9OnpHPv6tlUFmUnO5YkyXBFIfq/fI2fIDIBNXf2sLG6ljPNXVxRnsu6K0qZmqlTTSez4YpCWtARNS3qdn+hcPemeIcTkfjwqDEPzGD9igoWl+cmO5akgOGKQj6wk98VguieRA7Mj1coEYmf6DEPKouyuXVZGXlTNeaBRAw3yM7cBOYQkQQ42tDOr/fX6VRTGZJaYItMAqGws2V/nU41lRGpKIhMcLUtXWw9GaIk3KJTTWVEKgoiE1T0qaZhR6eaSkxi6ZL6JeB77l6dgDwiMgaaO3t4am8tNS2RU02XWIYKgsQklm3I14FvmdnLZvYJM8uPdygRGR13Z+/pFn748gmaOntYv6KCO1ZUkJWug8kSm1hGXvsO8B0zW0ykN9FrZvY88G133xrvgCISm86eEJv31XG0oUOnmsqoxXRMwczSgSuCr0YibbT/1Mz+wN3vj2M+EYnB0YZ2Nu+royekU03l0sRyTOEh4E7gaeD/c/dXgklfNLMD8QwnIsMb2NX03qpyinN0qqmMXixbCnuBz7l75yDTrhnjPCISI3U1lXiI5R30wMCCEAybibu3xCWViAwpHHZeOnqWR7efJBR27l09mxsvL1FBkDEx3HgKU4FsoHhAM7w8YGYCsonIANGnmi6pyOWmxepqKmNruN1HfwD8MZECEN0MrxX4l3iGEpE3c3f21bSy7YC6mkp8DdcQ78vAl83sj9z9KwnMJCJRunr72LK/noN1bcwunMZty8t1qqnEzXC7j97p7k8Dp83snoHT3f2ncU0mIpxs6mRjdS0d3X3csKiYqjmFpKXpVFOJn+F2H72DyGmodw4yzQEVBZE46Qs7LxxpZOcb5yiYlsn911RSlqcBECX+htt99Png+8cSF0dEzrZ381R1LfWt3Vw5O5+3LyohK0NnFkliDLf76E+HW9DdHxr7OCKTl7vz2qkWnjvUQEZ6GndeNZOFpTnJjiWTzHC7j3Rqg0iCRPctmluczbuWlpMzRZ3tJfGG2330d4kMIjJZHWvsYPO+Wrp7w9y0uISVlepbJMkz3O6jv3D3fzSzrxA5sPwm7v7puCYTmeB6+8L85nAju080a4hMSRnDbZ/uD77vSEQQkcmkvq2LjXtraWzvYdWcAm5YWKw2FZIShtt99Ivg+/cBzCwvctfbEpRNZMJxd3adaOb5w41MzUzjntWzuGzG9GTHEukXS+vsNcD3iBx4NjNrBj7u7jvjHU5kImnvDrFxby0nmjpZUJrDu5aUMS1LfYsktcRyesPDwB+6+3MAZnYDkSJxZTyDiUwkh+vb2Lyvnr5wmFuWlLF8Vp4OJktKiqUotF0oCADu/hsz0y4kkRj0hMI8c7CBvadbKMubyu3LyymanpXsWCJDGu7so9XBzVfM7JvAj4ichfRBYFv8o4mMb9GD4Fwzr4hr588gXX2LJMUNt6XwTwPufz7q9ltOURWRiHDY2fHGOV48cpbpU9K5d/VsKouykx1LJCbDnX20LpFBRCaCjl7nJ7tOcfrceRaX5/LOKzQIjowvw+0++rC7/2CoHkjqfSTyZq/XtrL1ZC/zMrq5bVk5SypydTBZxp3hdh9dOHl6VD2QzOxh4D1AvbsvH2T6TcDjwLHgoZ+6+xdG81oiydQd6mPr6w3sr2klN8v48NrLyM/WIDgyPg23++ibwffR9kD6V+CrwCPDzPOcu79nlM8vknS1LV08saeG1q5erp0/g/MZp1UQZFwb8bp6M/tHM8szs0wz22JmjWb24ZGWc/dngaYxSSmSYsJh55VjTTy6/SRhd96/ppLrFswgTbuLZJwz9+FPJDKz3e6+0szeB7wX+BNgq7tfNeKTm80FfjnM7qPHgFPAGeAz7l49xPM8CDwIUFZWVrVhw4aRXvot2tvbyclJvd70qZoLUjdbsnOdDzk760I0dDqzctNYWZJOVrqlRLahKNfFS9Vso821bt26ne6+ZsQZ3X3YL6A6+P5t4Pbg9qsjLRfMNxfYO8S0PCAnuL0eOBTLc1ZVVflobN26dVTLxVuq5nJP3WzJzHWors2/tvWwf/XpQ773dLOHw+E3Tdc6uzipmss9dbONNheww2P4jI3liuZfmNnrwHngD82sBOi6qBI1eDFqjbr9hJl9zcyK3b3xUp9bZKz19oV59mADr52KXJl8x/JyCnVlskxAIxYFd/+smX0RaHX3PjPrAO6+1Bc2s3Kgzt3dzK4hcnzj7KU+r8hYq2/r4qm9tZxt72HN3EKuX1CsK5NlwoqlS+r7gaeCgvA5YDXw/wK1Iyz3I+AmoNjMThG5IjoTwN2/AdwH/A8zCxHZCrk/2MQRSQnuzm9PNvObQ41My4xcmTxnhq5Mloktlt1Hf+3uPw66o94GfAn4OrB2uIXc/UMjTP8qkVNWRVJOZ0+ITdV1HGvsYH7JdN61tIzsLI2ZLBNfLO/yvuD7u4Gvu/vjZva38YskklzHGzvYFIyZvO6KUq6ana8rk2XSiKUonA66pN4CfNHMphDD9Q0i402oL8zzR86y641zFOdkcc/q2RTnaMxkmVxiKQofAG4HvuTuzWZWAfx5fGOJJFZTRw9P7Kmhoa2blZUF3LComEyNmSyTUCxnH3WaWT1wA3AICAXfRcY9d6f6TCvbDtSTkZ7GXStnsqAk9S5YEkmUWM4++jywBlhMZBjOTOAHwNviG00kvrp6+/j1/joO1bUzpyib25aXkzNFB5NlcovlL+B9wCpgF4C7nzGzUXVOFUkVJ5s62VhdS0d3H29fVEzVZYU6mCxCbEWhJ7jAzAHMbPpIC4ikqr6w8/LRs7xyvIn8aZncf00lZXlTkx1LJGXEUhT+Izj7qMDMfh/4OJE+SCLjSktnL0/uraGmpYtlM/O4aXEpWRk6mCwSLZYDzV8ys3cBrUSOK/yNu2+OezKRMfR6bStb9tdjButXVLC4XHtARQYzbFEws3Rgo7vfAqgQyLjTEwqz9UA9+860MrNgKrcvryB/mgbBERnKsEUh6HfUaWb57t6SqFAiY6G+rYsn99RyrrOHtfOKuHb+DNLUyE5kWLEcU+gC9pjZZqDjwoPu/um4pRK5BO7O7pPNPBfVyK6ySI3sRGIRS1H4VfAlkvLO9/SxaV8tRxvUyE5kNGI50Px9M8sCrgAcOODuPXFPJnKRTp3r5Km9tXT29PGOxSWsqizQtQciFymWK5rXA98EjgAGzDOzP3D3J+MdTiQW4bDz8rEmXj52lvxpmXzwal17IDJasWxXPwSsc/fDAGa2gMjuJBUFSbq2rl6e3FvL6XPnWVKRx7orSpiSkZ7sWCLjVixFof5CQQgcBerjlEckZkca2tlUXUfYnduWlbN0Zl6yI4mMe7EUhWozewL4DyLHFN4PbDezewDc/adxzCfyFn1hZ+uBenafaKY0bwrrl1dQOD0r2bFEJoRYisJUoA54R3C/ASgC7iRSJFQUJGGaOnp45lSIwr5mVs0p4IaFxWRo3AORMRPL2UcfS0QQkeG4O/tqWtl2oIGukHP3ypnM17gHImMulrOP5gF/BMyNnt/d74pfLJHf6Q718fT+el6vbWN24TTmz8lUQRCJk1h2H/0n8F3gF0A4vnFE3qyutYsn9tTQcr6X6xbM4Jq5RTz77JFkxxKZsGJqc+Hu/xz3JCJR3J1dJ5p5/nAj2VnpvH9NJbMKpiU7lsiEF0tR+HIwJOcmoPvCg+6+K26pZFLr7AmxqbqOY40dLCjN4dalZUzN1LUHIokQS1FYAXwEeCe/233kwX2RMXWyKdKqoqu3j3deUcqVs/PVqkIkgWIdo3m++h1JPEW3qijMzuLuVTMpzVWrCpFEi6UovAoUoKuYJU7au0M8tbeWk02dLKnI451XaJhMkWSJpSiUAa+b2XbefExBp6TKJXvjbAdP7a2lty/MrcvKWDYzP9mRRCa1WIrC5+OeQiadcNh58ehZth9vYsb0LNavmM2MnCnJjiUy6cVyRfMzZlYGXB089Iq7a1eSjFpbVy9P7qnldPN5ls/K56bFJWSqVYVISojliuYPAP8b2EZkPIWvmNmfu/tP4pxNJqBjjR1srK6lL+zcsaKcK8rV2VQklcSy++ivgKsvbB2YWQnwa0BFQWLWF3ZeONLIjuPnKMmdwrtXqLOpSCqKpSikDdhddBbQtr7ErOV8L0/uqaGmpYurKvO5cVGJOpuKpKhYisJTZrYR+FFw/4PEMOqamT0MvIfIID3LB5luwJeB9UAn8FFdJT3xHK5vZ9O+Wtzh3VdWcHlZbrIjicgwYjnQ/OfBgDo3EDmm8C13/1kMz/2vwFeBR4aYfgewKPhaC3w9+C4TQF/Yee5QA7890UxZ3lTWryinIFu7i0RS3ZBFwcwWAmXu/nwwutpPg8dvNLMF7j5sq0p3f9bM5g4zy93AI+7uwEtmVmBmFe5ec9E/haSUls5efrWnhrrWLg2EIzLODPeX+n+BtkEe7wymXapZwMmo+6eCx2QcO1TXxg9efoPm8z3cedVMblpcqoIgMo5Y5B/1QSaY7R3sWEAwbY+7rxjxySNbCr8c4pjCr4B/cPffBPe3AH/h7jsHmfdB4EGAsrKyqg0bNoz00m/R3t7O1OzpZKSlVnO19vZ2cnJSc8CYi8nWF3b2NPZxrCVM0VRjTXkG0zPjs64nyjpLJOW6eKmabbS51q1bt9Pd14w4o7sP+gUcHs20AfPNBfYOMe2bwIei7h8AKkZ6zqqqKh+NH/5ii39j22Fvau8e1fLxsnXr1mRHGFKs2Zrau/3fXjzuD2064M8cqPdQXzglciVDqmZTrouXqtlGmwvY4TF8bg+3Xb/dzH5/4INm9t+At/w3Pwo/B/6rRVwLtHgcjycUTDHCDk/ujVw4JWPjQG0b//7KCdq7Q9y9ciY3Xl5CeoptjYlI7IY7++iPgZ+Z2QP8rgisAbKItNMelpn9CLgJKDazU0R6KGUCuPs3gCeInI56mMhxio+N7keITXam8a6lZfzi1TM8f7iRGy8viefLTXihvjDPHWpk98lmZhZM5Y4VFeRNzUx2LBG5REMWBXevA643s3XAhWMCv3L3p2N5Ynf/0AjTHfhkrEHHwsLSHK6qzGfnG+eoLMpmXvH0RL78hNFyvpcn9tRQ29JF1WWFvG1hsbYORCaIWK5T2ApsTUCWhHj7ohJON3exqbqWD197GdOnxHL9nlxwtKGdjdV1hN2586oKFpbqYjSRiWTSnSuYmZ7G+uXl9ITCbHm9/sJBbhlBOOw8f7iRx3efIW9aBg+snaOCIDIBTcp/k2fkTOG6BTN47lAjB+ra1KlzBO3dIZ7cU8Opc+dZEbS61rUHIhPTpCwKAKvnFHK4vp2trzdQWZit3UhDONnUyZN7a+gJhbltWTlLZ6qAikxkk/bfvbQ049Zl5YT6wjyt3Uhv4e4caOrjsV2nmJKRzv3XzFFBEJkEJm1RACiansV1C2ZwuL6dw/XtyY6TMs739PH47jPsO9vH5WW5fOiaORRrqEyRSWFSFwWI7EYqyZ3CMwcb6AmFkx0n6Wpbuvjhy29woqmTq0rTuWN5OVkZk/5tIjJpTPq/9rQ0451XlNLWFeKVY03JjpM07s7uk838x46TmBkfvLqS+fnpRIa9EJHJQkdXgZkF01g6M49dJ86xdGYeRZNsmMjuUB+/3lfPwbo25pdM57Zl5UzNTGd/soOJSMJN+i2FCyI9/41tBybXQeeGtm5+9PIJDte3c8OiYu66aiZTM9OTHUtEkkRFITB9SgbXzZ/BG2c7OdIwOQ46V59p4dHtJ+jpC3Nv1Syunluk3UUik5yKQpSrZhdQND2LF46cJTyBO6n29oXZvK+OTdV1lOdP44G1lzG7MDvZsUQkBagoRElLM65fMIOz7T28XjvYoHPj37mOHh7dfpK9p1tYO6+Ie1bN0oV7ItJPnwYDLCzNoTRvCi8dPcvi8twJ1f3zUF0bm/bVkWbGe1fNUpdYEXkLbSkMYGZcv6CYlvO97D3dkuw4Y6Iv7Gw7UM8vX6thxvQsHrh2jgqCiAxKWwqDmDsjm1kF03jlWBNLZ+aROY6bv7V29fLEazXUtHSxck4BNy7SyGgiMrTx+2kXR2bG9Qtn0N4d4rVT43dr4XhjB//+8gnOdvTw7isrWLe4VAVBRIalLYUhzC7MZnbhNH574hwrKwvG1YdpOOy8dOwsrxxrYkbOFN6zooLCSXZBnoiMjrYUhrFmbhFtXSEOjKMzkTp7Qvzst6d5+WgTSyvyuP/qShUEEYmZthSGMXdGNsU5Wex8o4klFbkpf2HXqXOdPLmnlq7ePt61tIzls/KTHUlExhltKQzDzFh9WSGN7T2cOnc+2XGG5O7sON7EYztPk5lufPCaShUEERkVFYURXF6Wy9TMdF491ZzsKIPq6u3j56+e4blDjSwonc6H1s6hNHdqsmOJyDil3UcjyExPY9nMPH57opm2rl5yp2YmO1K/utYufvVaDW1dIW5aXMLKyoKU38UlIqlNWwoxuHJ2Po6z93RrsqMAkd1Fr55s5tHtJwm784GrZ7NqTqEKgohcMm0pxKAgO4vYJEGSAAANoklEQVTKwmz217Ry7fzkdhLtCYXZsr+O12vbmFucze3LKpiWpVbXIjI2tKUQo6Uz82g538vp5uQdcD7b3s2G7Sc4UNfG9Qtm8N6Vs1QQRGRMaUshRgtKcsjKSGPfmdaktJneX9PKlv11ZKance/q2VQWqdW1iIw9FYUYZWWksbA0h0P17ay7Ipywfki9fWGeOdDAntMtzCqcxvoVFeSo1bWIxIl2H12EJeV59ITCvHG2IyGvd66jhw3bT7LndAvXzCvivtWzVRBEJK70CXMRZhdOIzsrnYN17SwszY3rax2obePX++tIT9PYByKSOCoKFyEtzVhYmsP+mlZ6++KzCynUF+bZQw28erKFmQVTuWNFBXkpdG2EiExs2n10kS4vy6W3zzneOPa7kJo7e3h0x0lePdnCmrmF3FdVqYIgIgmlLYWLNKvgd7uQFpWN3S6k6KEy71o5kwUlOWP23CIisVJRuEhjvQupL+xsfb2e3SebKc+fyvoVFeRP09aBiCRHXHcfmdntZnbAzA6b2WcHmf5RM2sws93B13+PZ56xMla7kBrautl2KsTuk82smlPAB9ZUqiCISFLFbUvBzNKBfwHeBZwCtpvZz91934BZH3X3T8UrRzxc6i4kd+e3J5t5/lAjPX3o7CIRSRnx3H10DXDY3Y8CmNkG4G5gYFEYd9LSjEVlOew700pPKExWRuwbXO3dITbvq+V4YyfzS6ZzuWeoIIhIyjB3j88Tm90H3O7u/z24/xFgbfRWgZl9FPgHoAE4CPyJu58c5LkeBB4EKCsrq9qwYcNF52lvbycnZ+wO3jZ1hXnmZIgVxeksLBy5/5C7c6ItzJ7GPsIOy4vTmZeXRkdHx5jmGktjvc7GSqrmgtTNplwXL1WzjTbXunXrdrr7mhFndPe4fAHvB74Tdf8jwFcGzDMDmBLc/gTw9EjPW1VV5aOxdevWUS03nB/vOOnfeuaI94b6hp2vuaPHH9t50h/adMAffeWEn23vjmuusZKq2VI1l3vqZlOui5eq2UabC9jhMXx2x/NA8ymgMur+bODMgIJ01t27g7vfBqrimGfMrZ1XRHt3iD2nWwad3tkTYtuBer7/4nFqWrpYd0Up718zm6LpWYkNKiISo3geU9gOLDKzecBp4H7gv0TPYGYV7l4T3L0L2B/HPGNuduE0KouyefZgI1kZaSytyAPgXGcv+8608uqpZnr7wiybmc+184tSatQ2EZHBxK0ouHvIzD4FbATSgYfdvdrMvkBkM+bnwKfN7C4gBDQBH41XnngwM+68qoJfvlrDpuo6th1owAy6e8OYRdptX79gBjNypiQ7qohITOJ68Zq7PwE8MeCxv4m6/ZfAX8YzQ7xNyUjnvatm8XptK3WtXbhDcc4UFpTmqKOpiIw7+tQaA+lpxrKZ+SybmZ/sKCIil0QN8UREpJ+KgoiI9FNREBGRfioKIiLST0VBRET6qSiIiEg/FQUREemnoiAiIv3i1jo7XsysAXhjFIsWA41jHGcspGouSN1sqZoLUjebcl28VM022lyXuXvJSDONu6IwWma2w2PpJZ5gqZoLUjdbquaC1M2mXBcvVbPFO5d2H4mISD8VBRER6TeZisK3kh1gCKmaC1I3W6rmgtTNplwXL1WzxTXXpDmmICIiI5tMWwoiIjKCCVcUzOx2MztgZofN7LODTJ9iZo8G0182s7kJyFRpZlvNbL+ZVZvZ/xxknpvMrMXMdgdffzPYc8Up33Ez2xO87o5BppuZ/XOwzl4zs9UJyLQ4al3sNrNWM/vjAfMkbJ2Z2cNmVm9me6MeKzKzzWZ2KPheOMSyvxfMc8jMfi8Buf63mb0e/K5+ZmYFQyw77O89Drn+1sxOR/2+1g+x7LB/w3HK9mhUruNmtnuIZeO5zgb9nEj4+8zdJ8wXkWE/jwDzgSzgVWDpgHn+EPhGcPt+4NEE5KoAVge3c4GDg+S6CfhlktbbcaB4mOnrgScBA64FXk7C77WWyHnWSVlnwI3AamBv1GP/CHw2uP1Z4IuDLFcEHA2+Fwa3C+Oc61YgI7j9xcFyxfJ7j0OuvwU+E8Pveti/4XhkGzD9n4C/ScI6G/RzItHvs4m2pXANcNjdj7p7D7ABuHvAPHcD3w9u/wS42cwsnqHcvcbddwW324D9wKx4vuYYuxt4xCNeAgrMrCKBr38zcMTdR3PR4phw92eJjCMeLfq99H3gvYMsehuw2d2b3P0csBm4PZ653H2Tu4eCuy8Bs8fq9S4lV4xi+RuOW7bgs+ADwI/G8jVjMcznRELfZxOtKMwCTkbdP8VbP3z75wn+cFqAGQlJBwS7q1YBLw8y+Toze9XMnjSzZYnKBDiwycx2mtmDg0yPZb3G0/0M/UearHUGUObuNRD5gwZKB5kn2evu40S28gYz0u89Hj4V7NZ6eIjdIMleX28H6tz90BDTE7LOBnxOJPR9NtGKwmD/8Q88vSqWeeLCzHKAx4A/dvfWAZN3Edk9chXwFeA/E5Ep8DZ3Xw3cAXzSzG4cMD2Z6ywLuAv48SCTk7nOYpXMdfdXQAj44RCzjPR7H2tfBxYAK4EaIrtpBkra+gp8iOG3EuK+zkb4nBhysUEeG9V6m2hF4RRQGXV/NnBmqHnMLAPIZ3SbuRfFzDKJ/KJ/6O4/HTjd3VvdvT24/QSQaWbF8c4VvN6Z4Hs98DMim/DRYlmv8XIHsMvd6wZOSOY6C9Rd2I0WfK8fZJ6krLvgQON7gAc82Ok8UAy/9zHl7nXu3ufuYeDbQ7xe0t5rwefBPcCjQ80T73U2xOdEQt9nE60obAcWmdm84D/M+4GfD5jn58CFI/P3AU8P9UczVoL9lN8F9rv7Q0PMU37h2IaZXUPkd3M2nrmC15puZrkXbhM5SLl3wGw/B/6rRVwLtFzYnE2AIf9zS9Y6ixL9Xvo94PFB5tkI3GpmhcHukluDx+LGzG4H/h/gLnfvHGKeWH7vY50r+jjU+4Z4vVj+huPlFuB1dz812MR4r7NhPicS+z6Lx1H0ZH4ROVPmIJEzGP4qeOwLRP5AAKYS2RVxGHgFmJ+ATDcQ2ZR7DdgdfK0HPgF8IpjnU0A1kbMtXgKuT9D6mh+85qvB619YZ9HZDPiXYJ3uAdYkKFs2kQ/5/KjHkrLOiBSmGqCXyH9l/43IsagtwKHge1Ew7xrgO1HLfjx4vx0GPpaAXIeJ7F++8F67cLbdTOCJ4X7vcc71b8H75zUiH3QVA3MF99/yNxzvbMHj/3rhvRU1byLX2VCfEwl9n+mKZhER6TfRdh+JiMglUFEQEZF+KgoiItJPRUFERPqpKIiISD8VBUlZZtYXdKPca2a/GKrb5zDL/62ZfSaO+eZGd9oc8Pj5IPs+M/uGmelvTcYFvVEllZ1395XuvpzIVeefTHagi3DE3VcCVxLpdPmmJmZmlp6oIMFFh/pbl5jojSLjxYtENfgysz83s+1Bc7W/i3r8ryzSi//XwOKox7eZ2ZrgdrGZHQ9up5vZlyzSI/81M/uj4PEqM3smaHy2MarNQFXQgO9FYihSHmm6+AKw0CLjP2w1s38nchEXZvZhM3sl2Kr4ZpAn3cz+NdhC2mNmfxLM++lgy+M1M9sQPPamraFgmbnB134z+xqRHlGVZnarmb1oZrvM7MdBjx2RN1FRkJQX/Fd9M0G7AzO7FVhEpO/MSqDKzG40syoibRFWEelhc3UMT/8gMA9Y5e5XAj8M+s98BbjP3auAh4G/D+b/HvBpd78uxuzZQfY9wUPXELkSdqmZLQE+SKTJ2kqgD3gg+Jlmuftyd18RvCZEeulfyPmJGF5+MZGW56uADuBzwC0eaei2A/jTWH4GmVwykh1AZBjTLDIC1lxgJ5Ee8RDp63Ir8Nvgfg6RIpEL/MyDfj9mFkvPnFuItIEIAbh7k5ktB5YDm4PWSulAjZnlAwXu/kyw7L8Radg3mAVBdgced/cnzewm4BV3PxbMczNQBWwPXmcakWZnvwDmm9lXgF8Bm4L5XyNStP6T2DrCvuGR8S8gMjjSUuD54LWyiGx9ibyJioKksvPuvjL4MP4lkd01/0ykF9M/uPs3o2e2yHCdQ/VtCfG7LeOp0YsNsowB1QO3BoID3bH2hblwTGGgjgGv8313/8uBM5nZVUQGTvkkkUFfPg68m8ioYXcBf22R8SOify5488828LU2u/uHYswvk5R2H0nKc/cW4NPAZ4JdOxuBj1/YJ25ms8ysFHgWeJ+ZTQu6Wd4Z9TTHifxXDpHuuBdsAj5hkbbJmFkRcAAoMbPrgscyzWyZuzcDLWZ2Q7DsA5f4o20B7guyXxiL9zKLtP9Oc/fHgL8GVgcHiivdfSvwF0ABkS2k40SGlsQiY2fPG+K1XgLeZmYLg3mzzezyS8wvE5C2FGRccPffmtmrwP3u/m/B/vgXg10h7cCH3X2XmT1KpLvkG8BzUU/xJeA/zOwjwNNRj38HuBx4zcx6gW+7+1fN7D7gn4OtlAzg/xLpjPkx4GEz6+QSW2C7+z4z+xyRkbzSiHTt/CRwHvhe1BlDf0lkF9YPgjwG/B93bzazx4i0Nd9NpO30wSFeq8HMPgr8yMymBA9/bqj5ZfJSl1QREemn3UciItJPRUFERPqpKIiISD8VBRER6aeiICIi/VQURESkn4qCiIj0U1EQEZF+/z9kkBmgd0KFtwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "T = 450.0 \n",
    "Pr_list = np.arange(0.1, 20, 0.1)\n",
    "\n",
    "Z_list = []\n",
    "for Pr in Pr_list:\n",
    "    P = Pr * Pc\n",
    "    V = fsolve(vdw_calc, 1.0)[0]\n",
    "    Z_list.append(comp_factor())\n",
    "    \n",
    "plt.plot(Pr_list, Z_list, alpha=0.5)\n",
    "plt.xlabel('Reduced Pressure')\n",
    "plt.ylabel('Compressibility Factor')\n",
    "plt.grid(b=True)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.0"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": false,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}