Appendix B, Problem 15 - Perform Muskat-McDowell flash calculation on three component system.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 15\n",
    "Estimate K values with modified Wilson correlation and perform flash calculation of 3 component mixture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from tabulate import tabulate\n",
    "R = 10.73146 # Universal gas constant, psia-ft3/R-lbm mol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Per Table B-18 and B-19\n",
    "comps = ['C1', 'C3', 'nC5']\n",
    "z = [0.20, 0.32, 0.48]\n",
    "pc = [667.8, 616.3, 488.6]\n",
    "tc = [343.0, 665.7, 845.4]\n",
    "w = [0.0115, 0.1454, 0.251]\n",
    "\n",
    "# As specified by problem scenario\n",
    "pk = 2000 # Convergence pressure psia \n",
    "p = 500 # psia\n",
    "t = 160 # deg F\n",
    "A2 = 0.7 # problem characterizes this as A0, but I assume this is a misprint as A0 doesn't factor into these equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup a function we can call for the modified Wilson correlation\n",
    "def k_mod_wilson(w, tc, pc, t, p, pk, A2):\n",
    "    A1 = 1-((p-14.7)/(pk-14.7))**A2 # Eq 3.160 - assuming pressures given in psia, Eq 3.160 is in terms of psig\n",
    "    tr = t/tc\n",
    "    pr = p/pc\n",
    "    return (pc/pk)**(A1-1)*math.exp(5.37*A1*(1+w)*(1-1/tr))/pr # Eq 3.159"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Table B-19 Component Properties (Problem 15)\n",
      "\n",
      "Component      pci    tci      wi     Ki\n",
      "-----------  -----  -----  ------  -----\n",
      "C1           667.8  343    0.0115  9.208\n",
      "C3           616.3  665.7  0.1454  1.439\n",
      "nC5          488.6  845.4  0.251   0.358\n"
     ]
    }
   ],
   "source": [
    "# Calculate all the Ki's\n",
    "ncomps = len(comps)\n",
    "K = [k_mod_wilson(w[i], tc[i], pc[i], t+460, p, pk, A2) for i in range(ncomps)]\n",
    "\n",
    "# And print out summary\n",
    "print('Table B-19 Component Properties (Problem 15)\\n')\n",
    "header=['Component', 'pci', 'tci', 'wi', 'Ki']\n",
    "table = []\n",
    "for i in range(ncomps):\n",
    "    table.append([comps[i], pc[i], tc[i], w[i], round(K[i],3)])\n",
    "print(tabulate(table,headers=header))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Muscat-McDowell two-phase split\n",
    "def MM(Fv, *args):\n",
    "    z, K = args\n",
    "    return sum(z[i]/(Fv+(1/(K[i]-1))) for i in range(len(z))) # Eq 4.39"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vapor fraction associated with first guess Ki's: 0.48242\n"
     ]
    }
   ],
   "source": [
    "from scipy.optimize import brentq\n",
    "# Calculate vapor fraction to solve the Muscat-McDowell equation\n",
    "data = (z, K)\n",
    "Fv = brentq(MM,0,1, args = data) # Use SciPy Python solver brentq instead of Newton iterations (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.brentq.html)\n",
    "\n",
    "# And print result\n",
    "print(\"Vapor fraction associated with first guess Ki's:\", round(Fv,5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8FdX9//HXJwsESFizEnYIEAKyIyIgYVFUBNz3rS7VVmtrrUu1aq1tbau1X3/1q1+t+wJWVEBEEZAgCLLviwiyQ9hBArKf3x93otcYyM12J8l9Px+P++DOzJmZz0nC/dxzZuYcc84hIiISiii/AxARkcpDSUNEREKmpCEiIiFT0hARkZApaYiISMiUNEREJGRKGiIViJk5M2vlvX/VzB4v4XGaeceKOUWZdmY2t6SxeseobmYrzSy5NMeRykNJQwAws3VmdsTMEgusX+h9+DTzKa4iPzjN7FEvxl8VWP9rb/2jIZwnx8wOmdl+M/vWzOaZ2f1mVr2UVajI/gQ8mb/g/Q18Z2Z5Qa+GpzqAc+4w8DJwX3kGGpQEg2NbVJ7nlMIpaUiwtcCV+Qtm1gGo4V84xbIKuL7Auuu89aG6wzmXAKQBvwWuAMabmZVNiBWHmaUB2cDoApsucM7FB722hHC4t4Hrw5Rg6wbF1jEM55MClDQk2BsEPmjzXQ+8HlzA+0Z+c9DyDWY23XtvZva0mW03s31mttjM2nvbapjZU2a23ts23cxqeNveNbNcb/3nZpblrb8VuBq41/tm+eEpYp8D1AzaN4tAwptTIP5hXuvpWzNbY2aDCx7IOXfAOZcDDAXOAM739q1uZv8ysy3e61/5H5RmNtXMLvbe9/a+FZ/nLQ80s4VBMfzMzFaY2R4zm2BmTU9Rr9K62sw2mNlOM3swaP0gYL5z7lBRBzCzT8zsjgLrFpnZRQDOuU3AHqBnIfs29Fov9YPWdfbiiTWzVt7Pbp+37p3iVM77nezN/zvz1iV551SXWTlQ0pBgXwK1zSzTzKKBy4E3i7H/2UBfoDVQ19t/l7ftSaAr0AuoD9wLnPC2fQxkAMnAfOAtAOfcC977v3vfLC8o4vzBSa+whNfDW/c7L76+wLqTHcw5twGYC/TxVj1I4IOxE9AR6AE85G2bCvTz3vcFvgHOClqe6sUwHPg9cBGQBEwDRhRRr9LoDbQBBgAPm1mmt74D8FWIx3ibH7dA2wFNgY+Cyqwg8DP5Ea+lMhO4OGj1VcAo59xRAl1knwL1gEbA/wsxpvzjHwbeD44PuAyY6pzbXpxjSWiUNKSg/A/eQcBKYHMx9j0KJABtAXPOrXDObTWzKOBnwF3Ouc3OuePOuRnef3iccy875/Z7y48CHc2sTglifxO40sxiCXQtFUx4NwEvO+cmOudOeLGsLOKYWwgkOQi0eh5zzm13zu0A/ghc622byo+TxF+Dls/ytgP8HPir97M5BvwF6FSOrY0/Oue+c84tAhbxwwd7XWB/IeVHe9/c95pZftfVBwVivBp4P//359nvHbMw3ycdr6vvCm8dBP5mmgINnXOHnHPTi6jPzqD47il4fM9VQceXMqakIQW9QeA/3Q0U+KZeFOfcZ8C/gWeBbWb2gpnVBhKBOGBNwX3MLNrMnvC6ir7lh2/+iQXLeuWvDroQ+nGB828AVhP4IP7aObexwO6NC4uhCOnAbu99Q2B90Lb13joIfJtubWYpBFoirwONLXBjQQ/gc69cU+B/8j/4vGObd56QmdmyoJ9Dn1MUzQ16fxCI997vIZDgCxrunKvrvYYDOOf2E2hVXOGVuQKvNRgkAdh7khhGAWd4F9X7Ao5ACwsCLU4DZnt1+tkp6gKQGBRf/kX8z4AaZna6l9g6EUh0Ug6UNORHnHPrCVwQP49As7+gA0DNoOXUAvs/45zrCmQR6Kb6HbATOAS0LOR4VwHDgIFAHaCZtz7/4vOPhmF2zr0VdCH03EKO9zqBi9iFJbyNJ4mhUGbWmECXWv4H3BYCH/r5mnjrcM4dBOYBdwFLnXNHgBnA3cAa59zOoBh+HvTBV9c5V8M5NyPUuLzzZQX9HKYVvcdPLCbw+wnVCAKtuDMIXCuaUmB7JoGWTGGx7iXQBXUZgd/3COcNr+2cy3XO3eKca0igFfa/5t1yHCrn3AngvwRaG1cB47xEJ+VASUMKcxPQ3zl3oJBtC4GLzKym95/7pvwNZtbd+7YXSyC5HAKOe/+pXwb+6V0YjTazM7yLyAnAYQLXPmoSaCUE2wa0KEbs7xC4tvLfQra9BNxoZgPMLMrM0s2sbcFCXt3OAsYAs4Hx3qYRwEPehdZE4GF+3AU2FbiDH7qicgosAzwPPBB0wb6OmV1ajPqVlYlAFzOLC7H8eAIJ8zHgHe93CoCZpRPowvvyFPu/TaDb82KCuo7M7FIza+Qt7iHwJeF4qJUocPzLCXSdqWuqHClpyE8459Y450720NfTwBECH+av8eNuitrAiwT+868nkAjyuxDuAZYQuJtpN/A3An9/r3tlNwPL+ekHz0tAuwJ97KeK/Tvn3CTn3HeFbJsN3OjVYR+BD/PglsO/zWy/V7d/Ae8Bg4M+IB8ncGF8sVeX+d66fFMJJMHPT7KMc+4Dr+4jve64pUBhLaZy5ZzbRqBbZ1iI5fMvOA/kpx/KVwGvFbjGUdBYAjc7bPOur+TrDswyszyvzF3OubWh1eJH8c0i8EWlIYEbK6ScmCZhEolM3l1QrwE9XAk/CLzW4iKgr+5WigxKGiIiEjJ1T4mISMiUNEREJGRKGiIiErKTDptcWSUmJrpmzZqVeP8DBw5Qq1atsguoEoi0OkdafUF1jhSlqfO8efN2OueSiipX5ZJGs2bNmDu35FME5OTk0K9fv7ILqBKItDpHWn1BdY4Upamzma0vupS6p0REpBiUNEREJGRKGiIiEjJfk4aZvWyBCXuWnmR7P29yloXe6+FwxygiIj/w+0L4qwSG0j7VENzTnHNDwhOOiIiciq8tDefc5/wwV4GIiFRwvo89ZWbNCIx/376Qbf0IjDS6icC8Bfc455YVUu5W4FaAlJSUriNHjixxPHl5ecTHxxddsAqJtDpHWn1BdY4Upalzdnb2POdct6LK+d09VZT5QFPnXJ6ZnQeMJjC88o94c0m/ANCtWzdXkvuU9x48wmsz1lP3xAaG6N7uKi3S6guqc6QIR50r9N1TzrlvnXN53vvxQKw3+U2Zi4oynvnsa+bmlmT+FxGRyFChk4aZpXoT0WNmPQjEu6s8zlU7LpaOjeqwbJeShojIyfh9y+0IYCbQxsw2mdlNZnabmd3mFbkEWGpmi4BngCtKOllMKHpnJLF23wn2HTxaXqcQEanUfL2m4Zy7sojt/yZwS25Y9MlI5JnJXzPzm50Mbp8WrtOKiFQaFbp7Ktw6Na5LXDRM+3qn36GIiFRIShpBYqOjaFM/mumrlTRERAqjpFFA+wbRrN91kI27D/odiohIhaOkUUBWYjSAWhsiIoVQ0iggrZaRWjuO6bquISLyE0oaBZgZZ7ZK5Is1Ozl+wt8hVkREKholjUL0yUhk78GjLNuyz+9QREQqFCWNQpzZKjBSiW69FRH5MSWNQiQlVKdtaoKua4iIFKCkcRJ9MhKZt34P3x3RWFQiIvmUNE7izFaJHDl+gtnrNEeUiEg+JY2TOL15A6pFRzH96x1+hyIiUmEoaZxEjWrRdG1aTxfDRUSCKGmcQu+MRFbm7mfH/sN+hyIiUiEoaZxCn4zArbcz1qi1ISICShqnlNWwDnVqxKqLSkTEo6RxCtFRxpmtGjD9652U44SBIiKVhpJGEXq3SiL320Os2ZHndygiIr5T0ihC/nUNdVGJiChpFKlx/Zo0T6zF5BXb/Q5FRMR3ShohGHJaGjPW7GT7t4f8DkVExFdKGiEY1imdEw7GLtridygiIr5S0ghBq+R42qfXZsxCJQ0RiWxKGiEa3imdJZv36S4qEYloShohuqBjQ6IMxizY7HcoIiK+UdIIUUrtOHq1TGT0wi160E9EIpaSRjEM69SQDbsPMn/DXr9DERHxhZJGMQxun0r1mCjGLFQXlYhEJiWNYkiIi2VgZgrjFm/l6PETfocjIhJ2ShrFNKxTQ3YfOMJ0DSsiIhFISaOY+rVJpk6NWEari0pEIpCSRjFVi4nivA5pfLpsGwcOH/M7HBGRsFLSKIELO6fz3dHjfLo81+9QRETCSkmjBLo1rUd63RqMXqBhRUQksihplEBUlDG0U0Omr97JzrzDfocjIhI2SholNLxTOsdPOMZp5FsRiSBKGiXUJjWBtqkJvL9gs4YVEZGIoaRRCled3oTFm/Yxd/0ev0MREQkLJY1SuKRrI+rWjOWFz7/xOxQRkbDwNWmY2ctmtt3Mlp5ku5nZM2a22swWm1mXcMd4KjWrxXBtz6ZMWrFN82yISETwu6XxKjD4FNvPBTK8163Ac2GIqViuO6MZsdFRvDR9rd+hiIiUO1+ThnPuc2D3KYoMA153AV8Cdc0sLTzRhSYpoToXd0nnvXmbdPutiFR55vedP2bWDBjnnGtfyLZxwBPOuene8mTgPufc3ALlbiXQEiElJaXryJEjSxxPXl4e8fHxxdpnS94Jfj/9O4a1jOXCjGolPrdfSlLnyizS6guqc6QoTZ2zs7PnOee6FVUupkRHDx8rZN1Pspxz7gXgBYBu3bq5fv36lfiEOTk5lGT/z3bN4fP1e3ji+j7UqBZd4vP7oaR1rqwirb6gOkeKcNTZ72saRdkENA5abgRUyKfpbu3bkj0HjzJq/ia/QxERKTcVPWmMBa7z7qLqCexzzm31O6jCdG9Wj46N6/LStG84fkIP+4lI1eT3LbcjgJlAGzPbZGY3mdltZnabV2Q88A2wGngR+IVPoRbJzLi1TwvW7TrIxOXb/A5HRKRc+HpNwzl3ZRHbHfDLMIVTaudkpdC4fg1enPYNg9un+h2OiEiZq+jdU5VKTHQUN/duwbz1e5i3/lR3EouIVE5KGmXs0m6NqFNDQ4uISNWkpFHG8ocW+XT5NpZv+dbvcEREypSSRjm4pU8L6tSI5c/jl2vYdBGpUpQ0ykGdmrH8ekAGX6zexZSvtvsdjohImVHSKCdX92xKi8Ra/PmjFRw9fsLvcEREyoSSRjmJjY7i9+dlsmbHAUbM3uB3OCIiZUJJoxwNyEymV8sGPD1xFfu+O+p3OCIipaakUY7MjAfPz2Tvd0d5dspqv8MRESk1JY1yltWwDpd0acSrX6xj/a4DfocjIlIqShphcM85bYiOMv72yUq/QxERKRUljTBIqR3HbWe1ZPySXOas0/AiIlJ5KWmEyS19m5NaO47HP1rBCQ2dLiKVlJJGmNSsFsPvzmnDoo17GTlno9/hiIiUiJJGGF3YOZ0zWzXgzx8tZ+Pug36HIyJSbEoaYRQVZfz9ko6YGfe8u0jdVCJS6ShphFl63Ro8PKQds9bu5tUZ6/wOR0SkWJQ0fHBpt0b0b5vM3z5ZyZodeX6HIyISMiUNH5gZT1zUgbjYaO55dxHHNKChiFQSSho+Sa4dx5+Gt2fBhr28ME2z/IlI5aCk4aMLTkvjvA6pPD1xFStzNcufiFR8Sho+MjP+NKw9dWrEcvc7izhyTN1UIlKxKWn4rEF8df58YQeWb/2Wf0zQ2FQiUrEpaVQA52Slcm3Pprw4bS1jFm72OxwRkZNS0qgg/jCkHT2a1ee+9xazdPM+v8MRESmUkkYFUS0mimev7kK9mtX4+Rvz2JV32O+QRER+QkmjAklKqM7/XduVHXmHuePtBRzV8xsiUsEoaVQwpzWqy18v7MDMb3bxl/Er/A5HRORHYvwOQH7q4q6NWLblW17+Ym1gutiujfwOSUQEUEujwvr9eW3p1bIBv/9gCYs27vU7HBERQEmjwoqJjuLfV3UhOaE6N702h280sKGIVABKGhVY/VrVeO1nPXAOrvnPLDbv/c7vkEQkwilpVHAtk+J57Wc92H/4GNf+ZxY7dSuuiPhISaMSaJ9eh1du6M6Wfd9x3Uuz2ffdUb9DEpEIpaRRSXRrVp//u7YbX2/fz02vzuHgkWN+hyQiEUhJoxI5q3US/3NFZ+Zv2MPP35jH4WPH/Q5JRCKMkkYlc16HNJ646DSmfb2Tu0Ys1HDqIhJWShqV0GXdG/PIBe34ZFkut7w+V11VIhI2ISUNMzvDzJ41s8VmtsPMNpjZeDP7pZnVKenJzWywmX1lZqvN7P5Ctt/gnW+h97q5pOeqam48szl/u7gD077ewbUvzWbfQV0cF5HyV2TSMLOPgZuBCcBgIA1oBzwExAFjzGxocU9sZtHAs8C53vGuNLN2hRR9xznXyXv9p7jnqcou796EZ6/qwpJN+7j8hZls33/I75BEpIoLZeypa51zOwusywPme6+nzCyxBOfuAax2zn0DYGYjgWHA8hIcK2Kd2yGN+LgYfv7GPC59fiZv3nQ6jevX9DssEamiimxp5CcMM/uNmRU6cl4hSSUU6cDGoOVN3rqCLva6xUaZWeMSnKfK65ORxJs3n87eg0e55PkZrNq23++QRKSKMudcaAXNHgEuA3YDI4FRzrltJT6x2aXAOc65m73la4Eezrk7g8o0APKcc4fN7DbgMudc/0KOdStwK0BKSkrXkSNHljQs8vLyiI+PL/H+ftq0/wRPzj3E0ROOOzrFkdkgOqT9KnOdSyLS6guqc6QoTZ2zs7PnOee6FVnQOVesF3Aa8GdgJTCpuPsHHecMYELQ8gPAA6coHw3sK+q4Xbt2daUxZcqUUu3vt/U7D7gBT+W4Fg985F6fsdadOHGiyH0qe52LK9Lq65zqHClKU2dgrgvhs7skt9xuB3KBXUByCfbPNwfIMLPmZlYNuAIYG1zAzNKCFocCmpWoCE0a1OSDX/TirNZJ/GHMMh4cvVTPcohImQk5aZjZ7WaWA0wGEoFbnHOnlfTEzrljwB0E7spaAfzXObfMzB4LuhvrV2a2zMwWAb8Cbijp+SJJQlwsL17Xjdv7teTtWRu45qVZmnNcRMpEcWbuawr82jm3sKxO7pwbD4wvsO7hoPcPEOi2kmKKjjLuG9yWtqkJ3DtqMUP//QUvXteNdg1r+x2aiFRioTyncRGAc+5+YH25RyRlalindN697QyOn3Bc/NwMPliwye+QRKQSC6V76qGg95PLKxApP6c1qsvYO86kfXptfvPOIn7730UcOKyhR0Sk+EJJGnaS91KJJNeOY8QtPfnVgAw+WLCJC/7fdJZu3ud3WCJSyYSSNGqYWWcz6wrEee+75L/KO0ApOzHRUdw9qDVv39KTg0eOc9H/zuDl6Wvzb2kWESlSKBfCc4F/FvIewAE/edhOKraeLRow/q4+3DtqEY+NW07HpGhO636YBvHV/Q5NRCq4IpOGc65fGOKQMKtfqxovXteN12as4/Fxyzn76c95bFh7zj8treidRSRihXL3VO8ittc2s/ZlF5KEi5lxw5nNebRXDdLr1eCXb8/n9jfnsWO/nukQkcKFck3jYjObYWYPm9n5ZtbDzPqa2c/M7A1gHFCjnOOUctQoIYr3b+/FfYPbMnnldgY9PZXRCzbrWoeI/EQo3VO/MbN6wCXApUAq8B2Bp7ifd859Ub4hSjjEREdxe7+WDGqXzO9GLebX7yxk3OItPD68A6l14vwOT0QqiJCGEXHO7QFeByYC04CFwGFgQPmFJn5olZzAqNt68dD5mUxfvZMBT+Xwn2nfcPS4xq8SkeLNET4GuAA4SmASpjzgQHkEJf6KjjJu7tOCCb/uS4/m9Xn8oxWc/8w0Zq7Z5XdoIuKz4ow91cg5N7jcIpEKp2mDWrx8Q3cmrdjOHz9cxpUvfsnQjg158PxMUmqry0okEhWnpTHDzDqUWyRSIZkZg9qlMOnus/jVgAw+WZZL/ydz+L+pazh87Ljf4YlImIVyy+0SM1sM9Abmm9lX3vSr+eslAsTFRnP3oNZM/E1ferZowF8/XsmAp6YyZuFmTpzQXVYikSKU7qkh5R6FVBpNG9TipRu6M/3rnfxl/AruGrmQ/0xbywPntaVXy0S/wxORchbKLbcaDl1+ondGIuPu7M3ohZt5csJXXPXiLPq3Teb+c9vSOiXB7/BEpJyUZLpXEQCiooyLujTis3v6cd/gtsxZu5vB//qcu99ZyNqdurFOpCpS0pBSi4uN5vZ+LZl6bzY39W7O+KVbGfjPqfzu3UVs2HXQ7/BEpAwpaUiZqV+rGg+e347P783m+jOaMXbRFvo/lcP97y1m0x4lD5GqQElDylxyQhwPXxBIHtf0bMr78zeT/WQO941azDc78vwOT0RKQUlDyk1K7TgeHZrF1Hv7cWWPJoxeuJkB/5zKL96ap1kDRSqp4jwRLlIiaXVq8Niw9tzZP4NXvljLGzPXM35JLn0yEvlFv1b0bFEfM80kLFIZqKUhYZOUUJ17B7fliwf6c9/gtqzYup8rX/ySYc9+wZiFmzUookgloKQhYVc7Lpbb+7Vk+n3ZPD68PXmHjnHXyIX0/fsUnstZw76DR/0OUUROQklDfBMXG801PZsy6e6zePmGbrRIqsXfPllJz79O5g+jl+qiuUgFpGsa4ruoKKN/2xT6t01h+ZZvefmLtbwzZyNvfLmePhmJXNuzKQMyU4iO0nUPEb8paUiF0q5hbZ68tCP3DW7LyNkbeHv2Bm59Yx7pdWtw1elNuLx7YxLjq/sdpkjEUveUVEhJCdW5c0AG0+7N5vlrutAssSb/mPAVZ/x1MneNXMDMNbs0h7mID9TSkAotJjqKwe3TGNw+jdXb83jzy/W8N38TYxZuoXliLS7v3phLujZS60MkTNTSkEqjVXI8jw7NYvbvB/LUpR1JjK/GEx+vpOdfJnP7m/OYumoHxzW3h0i5UktDKp0a1aK5uGsjLu7aiNXb9zNy9kbem7+Jj5fmklYnjou6pHNJ18Y0T6zld6giVY6ShlRqrZITeGhIO343uA2Tlm9n1LyNPJezhmenrKFb03pc0rUR55+WRkJcrN+hilQJShpSJVSPieb809I4/7Q0tn17iA8WbObduRu5//0lPPrhMs7JSmV453T6tNLsgiKloaQhVU5K7ThuO6slP+/bgoUb9zJq3ibGLd7KmIVbSIyvRucGJ6jfai8d0utozCuRYlLSkCrLzOjcpB6dm9TjkQuyyPlqO6MXbubTZblM/PcXtEiqxbCO6Qzt1FDXP0RCpKQhEaFaTBRnZ6VydlYqH02cwv46LflgwWb+NXkVT09aRYf0Ogzt2JAhHdNIq1PD73BFKiwlDYk4tWKN83s04YoeTdi67zs+WryVsYu28OfxK/jLxyvo3qw+F3RsyOCsVJIS9PyHSDAlDYloaXVqcHOfFtzcpwVrdx7gw0VbGLtoC38YvZRHxiylZ4sGnH9aGoOzUmmgBwhFlDRE8jVPrMWvBmRwZ/9WrNqWx0eLtzBuyVYe/GApfxi9lF4tEzmvQxpnZ6XoCXSJWL4mDTMbDPwPEA38xzn3RIHt1YHXga7ALuBy59y6cMcpkcXMaJOaQJvUNvxmUGtW5u5n/JKtjFu8ld9/sISHRi+hR/P6nNs+jcHtU0mpHed3yCJh41vSMLNo4FlgELAJmGNmY51zy4OK3QTscc61MrMrgL8Bl4c/WolUZkZmWm0y02pz96DWrNi6n0+WbuXjpbk8MnYZj4xdRtem9Ti3fSrnZKXSuH5Nv0MWKVd+tjR6AKudc98AmNlIYBgQnDSGAY9670cB/zYzcxreVHxgZrRrWJt2DWtz99ltWL19Px8vyWX80lwe/2gFj3+0gqyGtRmclco57VPJSI7XcyBS5fiZNNKBjUHLm4DTT1bGOXfMzPYBDYCdYYlQ5BRaJSdw54AE7hyQwfpdB5iwLJdPluby1MRVPDVxFS0Sa3F2VirnZKXQsVFdojSJlFQB5teXdjO7FDjHOXezt3wt0MM5d2dQmWVemU3e8hqvzK4Cx7oVuBUgJSWl68iRI0scV15eHvHx8SXevzKKtDqXd333HjrB/O3HmbftGCt3n+C4g7rVjU7J0XRJjiazQTSxYU4gkfY7BtW5uLKzs+c557oVVc7PlsYmoHHQciNgy0nKbDKzGKAOsLvggZxzLwAvAHTr1s3169evxEHl5ORQmv0ro0irczjqO9z7d9/Bo0z5ajufLs8l56sd5Gw8THz1GM5qk8TZ7VLo1yaZOjXKfzDFSPsdg+pcXvxMGnOADDNrDmwGrgCuKlBmLHA9MBO4BPhM1zOkMqlTM5bhndMZ3jmdQ0ePM3PNLj5dnsvE5dv4aPFWYqKM01vUZ1BmCgPbpdConi6kS8XmW9LwrlHcAUwgcMvty865ZWb2GDDXOTcWeAl4w8xWE2hhXOFXvCKlFRcbTXbbZLLbJvP4cMfCjXuZuHwbk1Zs49EPl/Poh8vJTKvNoMxkBrZLoX3DOroOIhWOr89pOOfGA+MLrHs46P0h4NJwxyVS3qKjjK5N69G1aT3uP7cta3ceYNLybUxcvo1/T1nNM5+tJjmhOgMyUxiYmcyZrRKJi432O2wRPREuUhE0T6zFLX1bcEvfFuw5cIQpX21n0optfLhoCyNmbyAuNorerZIYmJlM/7bJJOuBQvGJkoZIBVOvVjUu6tKIi7o04vCx48z6ZjeTV2xj0opAIgHo2KgOAzJTGJCZTLu02noeRMJGSUOkAqseE03f1kn0bZ3Eo0MdK3P3f59Anp60in9OXEXDOnFkt01mYGYKZ7RsoG4sKVdKGiKVRPCQJnf0z2DH/sNMWRlofXywYDNvzdpAjdhozmyVqG4sKTdKGiKVVFJCdS7r3pjLujfm0NHjzFob6MaaHNSNdVqjOgxom0KdA8dxzqkbS0pNSUOkCoiLjeas1kmc1TqJPw51fLVtP5NXbGfyim38a/IqnIPnlk2mf9sUBrQN3I1Vo5q6saT4lDREqhgzo21qbdqm1uaX2a3YlXeY58Z8zhbqfX83VvWYKM5slUj/tskMyEzWFLcSMiUNkSquQXx1eqfH0q9fV44cO8GstbsCrZCV2/hs5XYeGg1ZDWszoG0yAzJT6JCuhwrl5JQ0RCJItZgo+mQk0ScjiUcuaMfq7XlMWrGdz1b+8FBhUkL17xM908znAAAQUUlEQVRIb3VjSQFKGiIRyszISEkgIyWB2/u1ZPeBI+R8tZ3JK7YzbvFWRs7Z+H031kDvmRDNUihKGiICQP2ghwqPHDvB7LW7mbQiMDbWZyu3wweBu7EGZqYwMDOFzLQE3Y0VgZQ0ROQnqsVE0Tsjkd4ZiTxyQTtWbcv7PoEEP1Q4sF0ggfRs0YBqMVF+hy1hoKQhIqdkZrRJTaBNagK/zG71/UOFE1ds4925m3h95npf5ggRfyhpiEixFHyocPrXO71WyPbv5wjp0bw+Z7dLYVBWKul1dTtvVaKkISIlFhcbHeiiapfCiROOBRv3MmnFNj5dlvv9HCHt02szKDOVs7NSaJuq6yCVnZKGiJSJqKA5Qu4b3JY1O/KY6M0R8q/Jq3h60ioa16/B2e1SOScrla5N6xGt50EqHSUNESkXLZPiaXlWPLed1ZLt+w8xecV2JizL5Y2Z63lp+loS46sxqF0K52Sl0qtloi6kVxJKGiJS7pIT4riyRxOu7NGE/YeOMuWrHUxYlsvYhVsYMXsjCdVj6J+ZzOCsVPq1SdYDhRWYkoaIhFVCXCxDOzZkaMeGHDp6nBlrdvLJ0lwmLt/GmIVbqBEbTb82SQxun0r/tskkxOlOrIpESUNEfBMXG03/tin0b5vCseOBBwo/XprLJ8ty+XhpLtWio+iTkch5HdIY2C5Ft/JWAEoaIlIhxERH0atVIr1aJfLHoVnM27CHj5fk8snSrUxeuZ3YaKNvRpISiM+UNESkwomKMro3q0/3ZvV56PxMFm7ay/jFW/l4ae73CaRPRhJDTktjULsUdWGFkZKGiFRoUVFGlyb16NKkHg+en8nCjXsZv2Qr45fk8tnK7VSLiSK7TRJDTmvIgMxkalbTx1p50k9XRCoNM6Nzk3p0blKP35+XyfwNexm3eAsfLd7KhGXbqBEbzYDMZIZ2bIidcH6HWyUpaYhIpWT2w8OED53fjjnrdvPhoi18vDSXcYu3UjMGhu5ZzLBO6ZzevL4mliojShoiUulFRxk9WzSgZ4sGPDo0i+mrd/LChPmMXbSFkXM2klo7jgs6pjGsUzpZDWtrKJNSUNIQkSolNjqK7DbJ2NY4Tu/Vh0krtjFm4WZe+WIdL05bS+uUeC7s3IjhnRtqbvQSUNIQkSqrRrVoLujYkAs6NmTPgSN8tGQr78/fxN8+WcnfJ6ykV8sGXNi5EYPbpxJfXR+HodBPSUQiQr1a1bimZ1Ou6dmUdTsP8MGCzXywYDP3vLuIP4xeyrkdUrm0a2Nd/yiCkoaIRJxmibX4zaDW/HpgBvPW7+G9+Zv4cNFW3p+/mcb1a3BJl8Zc1CWdxvVr+h1qhaOkISIRy8zo1qw+3ZrV5+EhWUxYlsu78zZ+P5R7r5YNuLx7Y87JSiUuVoMogpKGiAgQuP4xvHM6wzuns2nPQd6fv5l3523krpELqVszlgs7p3Nljya0TknwO1RfKWmIiBTQqF5NfjUggzuyWzFjzS5GzNnAm1+u55Uv1tGlSV2u6N6EIR3TIvLp88irsYhIiKKijN4ZifTOSGT3gSO8P38TI2Zv4N73FvOnj5ZzcZdGXNOzCa2SI6f1oaQhIhKC+rWqcXOfFtzUuzlz1u3hrVnreXvWBl6dsY7Tm9fnmp5NOScrtcrPQKikISJSDGZGj+b16dG8Pg8POcx/527i7dnruXPEAhLjq3NF98Zc07MpqXXi/A61XChpiIiUUIP46tzeryU/79uCqV/v4K0v1/Nszmqen7qGwe1TufHMZnRpUq9KDVuipCEiUkpRUUZ2m2Sy2ySzYddBXp+5jnfmbmTc4q10SK/DDb2aMaRjGtVjKv9tu1W7801EJMyaNKjJQ0Pa8eUDA/jT8PYcPHKM3767iDOf+IxnJn/N7gNH/A6xVHxJGmZW38wmmtnX3r/1TlLuuJkt9F5jwx2niEhJ1aoew7U9mzLp7rN446YetE+vwz8nrqLXE5N5aPQS1u484HeIJeJX99T9wGTn3BNmdr+3fF8h5b5zznUKb2giImXHLDA1bZ+MJFZt289/pn3Df+ds4q1ZGxiUmcKtfVvQrVl9v8MMmV/dU8OA17z3rwHDfYpDRCRsWqck8PdLOjL9/mzuyG7F7HW7ueT5mVz2/EymfLUd5yr+bIN+JY0U59xWAO/f5JOUizOzuWb2pZkpsYhIlZCcEMdvz27DjPv78+gF7di05yA3vjKH85+ZzrjFWzhegaeqtfLKbGY2CUgtZNODwGvOubpBZfc4535yXcPMGjrntphZC+AzYIBzbk0h5W4FbgVISUnpOnLkyBLHnZeXR3x8fIn3r4wirc6RVl9QnSu6YyccM7cc46O1R8k94EipaZzXIpYzG8YQU4xh2ktT5+zs7HnOuW5FlSu3pHHKk5p9BfRzzm01szQgxznXpoh9XgXGOedGnapct27d3Ny5c0scW05ODv369Svx/pVRpNU50uoLqnNlcfyE49NluTybs5qlm78lvW4N7ujfiou7NArpSfPS1NnMQkoafnVPjQWu995fD4wpWMDM6plZde99InAmsDxsEYqIhFl0lHFuhzQ+vKM3r9zYncSE6jzw/hL6P5XDyNkbOHr8hN8h+pY0ngAGmdnXwCBvGTPrZmb/8cpkAnPNbBEwBXjCOaekISJVnlngYcHRv+jFKzd2p0F8de5/fwnZT/qfPHy55dY5twsYUMj6ucDN3vsZQIcwhyYiUmHkJ49+rZPIWbWDf01cxf3vL+G5qWv47dltGNIhLexT0+qJcBGRCu77lscvz+Sl67tRIzaaX41YwJD/N52cMN+qq7GnREQqCTNjQGYK2W2S+XDxFp76dBU3vDKH05vX597BbcMSg1oaIiKVTFSUMaxTOpPuPovHhmWxZscBLn5uBs8uPFTurQ61NEREKqlqMVFcd0YzLu7SiFe+WMvK1WvLfRh2JQ0RkUquVvUY7uifQU7U5nI/l7qnREQkZEoaIiISMiUNEREJmZKGiIiETElDRERCpqQhIiIhU9IQEZGQKWmIiEjIfJmEqTyZ2Q5gfSkOkQjsLKNwKotIq3Ok1RdU50hRmjo3dc4lFVWoyiWN0jKzuaHMXlWVRFqdI62+oDpHinDUWd1TIiISMiUNEREJmZLGT73gdwA+iLQ6R1p9QXWOFOVeZ13TEBGRkKmlISIiIVPSEBGRkEVk0jCzwWb2lZmtNrP7C9le3cze8bbPMrNm4Y+ybIVQ57vNbLmZLTazyWbW1I84y1JRdQ4qd4mZOTOr9LdnhlJnM7vM+10vM7O3wx1jWQvhb7uJmU0xswXe3/d5fsRZVszsZTPbbmZLT7LdzOwZ7+ex2My6lGkAzrmIegHRwBqgBVANWAS0K1DmF8Dz3vsrgHf8jjsMdc4Ganrvb4+EOnvlEoDPgS+Bbn7HHYbfcwawAKjnLSf7HXcY6vwCcLv3vh2wzu+4S1nnvkAXYOlJtp8HfAwY0BOYVZbnj8SWRg9gtXPuG+fcEWAkMKxAmWHAa977UcAAK++Jd8tXkXV2zk1xzh30Fr8EGoU5xrIWyu8Z4E/A34FD4QyunIRS51uAZ51zewCcc9vDHGNZC6XODqjtva8DbAljfGXOOfc5sPsURYYBr7uAL4G6ZpZWVuePxKSRDmwMWt7krSu0jHPuGLAPaBCW6MpHKHUOdhOBbyqVWZF1NrPOQGPn3LhwBlaOQvk9twZam9kXZvalmQ0OW3TlI5Q6PwpcY2abgPHAneEJzTfF/f9eLDFldaBKpLAWQ8H7jkMpU5mEXB8zuwboBpxVrhGVv1PW2cyigKeBG8IVUBiE8nuOIdBF1Y9Aa3KambV3zu0t59jKSyh1vhJ41Tn3lJmdAbzh1flE+Yfni3L9/IrElsYmoHHQciN+2lz9voyZxRBo0p6qOVjRhVJnzGwg8CAw1Dl3OEyxlZei6pwAtAdyzGwdgb7fsZX8Yniof9tjnHNHnXNrga8IJJHKKpQ63wT8F8A5NxOIIzCwX1UV0v/3korEpDEHyDCz5mZWjcCF7rEFyowFrvfeXwJ85rwrTJVUkXX2umr+j0DCqOz93FBEnZ1z+5xzic65Zs65ZgSu4wx1zs31J9wyEcrf9mgCNz1gZokEuqu+CWuUZSuUOm8ABgCYWSaBpLEjrFGG11jgOu8uqp7APufc1rI6eMR1TznnjpnZHcAEAndevOycW2ZmjwFznXNjgZcINGFXE2hhXOFfxKUXYp3/AcQD73rX/Dc454b6FnQphVjnKiXEOk8Azjaz5cBx4HfOuV3+RV06Idb5t8CLZvYbAt00N1TmL4FmNoJA92Kid53mESAWwDn3PIHrNucBq4GDwI1lev5K/LMTEZEwi8TuKRERKSElDRERCZmShoiIhExJQ0REQqakISIiIVPSkErNzHLM7JwC635tZv8bxhgeNbPNZrbQez1RBsesa2a/CFpuaGajSntc71g53qiw+fFeUhbHlcigW26lUjOznwM9nXM3Bq37ksDzB9PK6ZzRzrnjQcuPAnnOuSdD3SeEczQDxjnn2pci1JMdOwe4p5I/yCg+UUtDKrtRwBAzqw7ff9g2BKabWbw3N8h8M1tiZsPyy5jZSjN7zZtvYJSZ1fS2DfDmXVjizVuQf9x1ZvawmU0HLg0lsIL7mNktZjbHzBaZ2XtB50wxsw+89YvMrBfwBNDSawn8w4t5qVc+zsxe8WJcYGb5T3jfYGbvm9knZva1mf091B+imf2tQMvmUTP7baj7S+RQ0pBKzXuaeTaQP1pr/vwnjsBw5xc657oQGDrjKbPvh7hvA7zgnDsN+Bb4hZnFAa8ClzvnOhAYMeH2oNMdcs71ds6NLCSU3wR195xzkn3ed851d851BFYQGBMJ4Blgqre+C7AMuB9Y45zr5Jz7XYFz/dKrewcCg/G95sUO0Am4HOgAXG5mjSncW0HxNiAwpPjlQdsvA949yb4SwZQ0pCoYwQ9DvVzhLUNgtM+/mNliYBKB4aFTvG0bnXNfeO/fBHoTSCRrnXOrvPWvEZjwJt87p4jhae8DvpNzbsJJ9mlvZtPMbAlwNZDlre8PPAfgnDvunNtXRH17A2945VcC6wmMIQUw2RtX6xCwHDjZDIxXB8W7yzm3AEj2rp10BPY45zYUEYdEoIgbe0qqpNHAPy0wrWUN59x8b/3VQBLQ1Tl31AKj2eZ/Iy94Mc9R+JDSwQ6UILbgfV4FhjvnFpnZDQTGDyqJU8UZPDrxcYr3f3wUgQE6Uwm0PER+Qi0NqfScc3lADvAyP7QyIDCk/XYvYWTz42/dTSwwtwIEunimAyuBZmbWylt/LTC1DENNALaaWSyBhJZvMl43mJlFm1ltYL9XvjCf5+9vZq2BJgSGOC+tkQRaapcQSCAiP6GkIVXFCKAjP/6G/BbQzczmEviQXRm0bQVwvdd1VR94zuvSuZHASL9LgBPA82UY4x+AWcDEArHcBWR755wHZHnXar4ws6Vm9o8Cx/lfINor/w6BUVtLPf+Jc24ZgUS1uSyH0paqRbfcSsQpz9tZRao6tTRERCRkammIiEjI1NIQEZGQKWmIiEjIlDRERCRkShoiIhIyJQ0REQnZ/wehCts0BO1vEQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Just for laughs, plot the associated h(Fv) plot\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "Fv_plot = np.linspace(0,1)\n",
    "plt.plot(Fv_plot,[MM(Fv_plot, *data) for Fv_plot in Fv_plot])\n",
    "plt.xlabel('Vapor Fraction Fv')\n",
    "plt.ylabel('h(Fv)')\n",
    "plt.title(\"Muscat-McDowell - h(Fv) vs Fv\")\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Table B-21 Final Flash Calculation Results (Problem 15)\n",
      "\n",
      "Component      zi  Ki      ci         zi/(Fv+ci)    zi/(Fv+ci)**2         xi        yi\n",
      "-----------  ----  ------  -------  ------------  ---------------  ---------  --------\n",
      "C1           0.2   9.2082  0.1218     0.330988          0.547766   0.0403243  0.371313\n",
      "C3           0.32  1.439   2.2777     0.115936          0.0420036  0.26407    0.380006\n",
      "nC5          0.48  0.3575  -1.5564   -0.446924          0.416128   0.695606   0.248681\n",
      "Total        1                       -7.4496e-14        1.0059     1          1\n"
     ]
    }
   ],
   "source": [
    "# And print out summary\n",
    "c = [1/(K[i]-1) for i in range(ncomps)]\n",
    "x = [z[i]/(Fv*(K[i]-1)+1) for i in range(ncomps)] # Eq 4.41\n",
    "y = [K[i]*x[i] for i in range(ncomps)]\n",
    "print('Table B-21 Final Flash Calculation Results (Problem 15)\\n')\n",
    "header=['Component', 'zi', 'Ki', 'ci', 'zi/(Fv+ci)', 'zi/(Fv+ci)**2', 'xi', 'yi']\n",
    "table = []\n",
    "for i in range(ncomps):\n",
    "    table.append([comps[i], z[i], round(K[i],4), round(c[i],4), z[i]/(Fv+c[i]), z[i]/(Fv+c[i])**2, x[i], y[i] ])\n",
    "table.append(['Total', sum(z), ' ', ' ', sum(z[i]/(Fv+c[i]) for i in range(ncomps)), sum(z[i]/(Fv+c[i])**2  for i in range(ncomps)), sum(x), sum(y) ])\n",
    "print(tabulate(table,headers=header))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}