path: root/Chemical_Engineering_Thermodynamics/ch10_1.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Chapter 10 : Residual Properties by Equations of State"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.2  Page Number : 334"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Preparation of fugacity and fugacity coefficient\n",
      "\n",
      "import math\n",
      "\n",
      "# Variables\n",
      "T = 40 + 273.15;\t\t\t#[C] - Temperature\n",
      "P_1 = 0;\t\t\t#[bar]\n",
      "P_2 = 10;\t\t\t#[bar]\n",
      "V_liq = 90.45;\t\t\t#[cm**(3)/mol]\n",
      "V_liq = V_liq*10**(-6);\t\t\t#[m**(3)/mol]\n",
      "P_sat = 4.287;\t\t\t#[bar] \n",
      "\n",
      "# For butadiene\n",
      "T_c = 425.0;\t\t\t#[K] - Critical temperature\n",
      "P_c = 43.3;\t\t\t#[bar] - Critical pressure\n",
      "P_c = P_c*10**(5);\t\t\t#[N/m**(2)]\n",
      "w = 0.195;\t\t\t# Acentric factor\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "# Calculations and Results\n",
      "# Let us calculate second virial coefficient at 40 C\n",
      "Tr = T/T_c;\t\t\t# Reduced temperature\n",
      "B_0 = 0.083-(0.422/(Tr)**(1.6));\n",
      "B_1 = 0.139-(0.172/(Tr)**(4.2));\n",
      "#We know,(B*Pc)/(R*Tc) = B_0+(w*B_1)\n",
      "B = ((B_0 + (w*B_1))*(R*T_c))/P_c;\t\t\t#[m**(3)/mol] - Second virial coefficient\n",
      "\n",
      "# math.log(f/P) = (B*P)/(R*T)\n",
      "# f = P*exp((B*P)/(R*T))\n",
      "\n",
      "print \" The table is as follows\"\n",
      "print \" Pbar \\t\\t fbar \\t\\t                 phi\";\n",
      "from numpy import zeros\n",
      "P = [1,2,3,4,4.287,5,6,8,10];\n",
      "f = zeros(9);\n",
      "phi = zeros(9);\n",
      "for i in range(5):\n",
      "    f[i]=P[i]*(math.exp((B*P[i]*10**(5))/(R*T)));\t\t\t#[bar]  # Pressure inside the exponential term has to be in N/m**(2)\n",
      "    phi[i]= (f[i]/P[i]);\n",
      "    print \" %f \\t %f \\t\\t\\t %f\"%(P[i],f[i],phi[i]) \n",
      "\n",
      "f_sat = f[4];\n",
      "\n",
      "# From pressure of 4.287 bar onwards the valid equation to compute fugacity of compressed liquid is given by\n",
      "# f = f_sat*exp[V_liq*(P-P_sat)/(R*T)]\n",
      "\n",
      "for j in range(5,9):\n",
      "    f[j] = f_sat*math.exp((V_liq*(P[j]-P_sat)*10**(5))/(R*T));\t#[bar]  # Pressure inside the exponential term has to be in N/m**(2)\n",
      "    phi[j] = f[j]/P[j]; \n",
      "    print \" %f \\t %f \\t\\t\\t %f\"%(P[j],f[j],phi[j]);\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The table is as follows\n",
        " Pbar \t\t fbar \t\t                 phi\n",
        " 1.000000 \t 0.978336 \t\t\t 0.978336\n",
        " 2.000000 \t 1.914284 \t\t\t 0.957142\n",
        " 3.000000 \t 2.809221 \t\t\t 0.936407\n",
        " 4.000000 \t 3.664484 \t\t\t 0.916121\n",
        " 4.287000 \t 3.902801 \t\t\t 0.910380\n",
        " 5.000000 \t 3.912481 \t\t\t 0.782496\n",
        " 6.000000 \t 3.926097 \t\t\t 0.654349\n",
        " 8.000000 \t 3.953471 \t\t\t 0.494184\n",
        " 10.000000 \t 3.981037 \t\t\t 0.398104\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.3  Page Number : 334"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of enthalpy entropy and internal energy change\n",
      "\n",
      "\n",
      "import math \n",
      "from scipy.integrate import quad \n",
      "\n",
      "# Variables\n",
      "n = 100.;\t\t\t#[mol] - No of moles\n",
      "T_1 = 600.;\t\t\t#[K] - Initial temperature\n",
      "T_2 = 300.;\t\t\t#[K] - Final temperature\n",
      "P_1 = 10.;\t\t\t#[atm] - Initial pressure\n",
      "P_1 = P_1*101325;\t\t\t#[Pa]\n",
      "P_2 = 5.;\t\t\t#[atm] - Final presssure\n",
      "P_2 = P_2*101325;\t\t\t#[Pa]\n",
      "Tc = 369.8;\t\t\t#[K] - Critical temperature\n",
      "Pc = 42.48;\t\t\t#[bar] - Critical pressure\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.152;\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "# Calculations\n",
      "# At 600 K\n",
      "Tr = T_1/Tc;\t\t\t# Reduced temperature\n",
      "B_0 = 0.083-(0.422/(Tr)**(1.6));\n",
      "B_1 = 0.139-(0.172/(Tr)**(4.2));\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B = ((B_0 + (w*B_1))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol] - Second virial coefficient\n",
      "dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);\t\t\t# (dB_0/dT)\n",
      "dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);\t\t\t# (dB_1/dT)\n",
      "dB_dT = ((R*Tc)/(Pc))*(dB0_dT + w*dB1_dT);\t\t\t# dB/dT\n",
      "\n",
      "# Now let us calculate B and dB/dT at 300 K\n",
      "Tr_prime = T_2/Tc;\t\t\t# Reduced temperature\n",
      "B_0_prime = 0.083-(0.422/(Tr_prime)**(1.6));\n",
      "B_1_prime = 0.139-(0.172/(Tr_prime)**(4.2));\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B_prime = ((B_0_prime + (w*B_1_prime))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol] - Second virial coefficient\n",
      "dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);\t\t\t# (dB_0/dT)\n",
      "dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);\t\t\t# (dB_1/dT)\n",
      "dB_dT_prime = ((R*Tc)/(Pc))*(dB0_dT_prime + w*dB1_dT_prime);\t\t\t# dB/dT\n",
      "\n",
      "# The change in enthalpy for ideal gas is given by\n",
      "\n",
      "def f16(T): \n",
      "    return -0.966+7.279*10**(-2)*T-3.755*10**(-5)*T**(2)+7.58*10**(-9)*T**(3)\n",
      "\n",
      "delta_H_ig =  quad(f16,T_1,T_2)[0]\n",
      "\n",
      "delta_H_ig = delta_H_ig*4.184;\t\t\t#[J/mol]\n",
      "\n",
      "# We know that delta_H_ig = delta_U_ig + R*delta_T. Therefore change in internal energy is given by \n",
      "delta_U_ig = delta_H_ig - R*(T_2 - T_1);\t\t\t#[J/mol]\n",
      "\n",
      "# The change in entropy of ideal gas is given by \n",
      "\n",
      "def f17(T): \n",
      "    return Cp_0/T\n",
      "\n",
      "#delta_S_ig =  quad(f17,T_1,T_2) - R*math.log(P_2/P_1)[0]\n",
      "\n",
      "\n",
      "def f18(T): \n",
      "    return (-0.966+7.279*10**(-2)*T-3.755*10**(-5)*T**(2)+7.58*10**(-9)*T**(3))/T\n",
      "\n",
      "delta_S_ig =  quad(f18,T_1,T_2)[0]*4.184 - R*math.log(P_2/P_1)\n",
      "\n",
      "\n",
      "# Now let us calculate the change in enthalpy of gas. We know that\n",
      "# delta_H = delta_H_ig + delta_H_R\n",
      "# delta_H_R = H_2_R - H_1_R\n",
      "H_2_R = B_prime*P_2 - P_2*T_2*dB_dT_prime;\t\t\t# [J/mol]\n",
      "H_1_R = B*P_1 - P_1*T_1*dB_dT;\t\t\t# [J/mol]\n",
      "delta_H_R = H_2_R - H_1_R;\t\t\t# [J/mol]\n",
      "delta_H = delta_H_ig + delta_H_R;\t\t\t# [J/mol]\n",
      "\n",
      "# Let us calculate the residual entropy of gas\n",
      "S_2_R = -P_2*dB_dT_prime;\t\t\t#[J/mol-K]\n",
      "S_1_R = -P_1*dB_dT;\t\t\t#[J/mol-K]\n",
      "delta_S = delta_S_ig + (S_2_R - S_1_R);\t\t\t#[J/mol-K]\n",
      "\n",
      "# Let us calculate the residual internal energy of gas\n",
      "U_2_R = -P_2*T_2*dB_dT_prime;\t\t\t#[J/mol-K]\n",
      "U_1_R = -P_1*T_1*dB_dT;\t\t\t#[J/mol-K]\n",
      "delta_U = delta_U_ig + (U_2_R - U_1_R);\t\t\t#[J/mol-K]\n",
      "\n",
      "# For 100 mol sample,\n",
      "delta_H_ig = delta_H_ig*n*10**(-3);\t\t\t#[kJ/mol]\n",
      "delta_H = delta_H*n*10**(-3);\t\t\t#[kJ/mol]\n",
      "\n",
      "delta_U_ig = delta_U_ig*n*10**(-3);\t\t\t#[kJ/mol]\n",
      "delta_U = delta_U*n*10**(-3);\t\t\t#[kJ/mol]\n",
      "\n",
      "delta_S_ig = delta_S_ig*n*10**(-3);\t\t\t#[kJ/mol]\n",
      "delta_S = delta_S*n*10**(-3);\t\t\t#[kJ/mol]\n",
      "\n",
      "# Results\n",
      "print \" The value of delta_H = %f kJ/mol\"%(delta_H);\n",
      "print \" The value of delta_H_ig ideal gas)= %f J/mol\"%(delta_H_ig);\n",
      "print \" The value of delta_U = %f kJ/mol\"%(delta_U);\n",
      "print \" The value of delta_U_ig ideal gas) = %f kJ/mol\"%(delta_U_ig);\n",
      "print \" The value of delta_S = %f kJ/mol\"%(delta_S);\n",
      "print \" The value of delta_S_ig ideal gas) = %f kJ/mol\"%(delta_S_ig);\n",
      "\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The value of delta_H = -3130.406247 kJ/mol\n",
        " The value of delta_H_ig ideal gas)= -3096.763542 J/mol\n",
        " The value of delta_U = -2867.753839 kJ/mol\n",
        " The value of delta_U_ig ideal gas) = -2847.343542 kJ/mol\n",
        " The value of delta_S = -6.466831 kJ/mol\n",
        " The value of delta_S_ig ideal gas) = -6.358994 kJ/mol\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.4  Page Number : 337"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of molar heat capacity\n",
      "\n",
      "# Variables\n",
      "T = 35 + 273.15;\t\t\t#[K] - Temperature\n",
      "P = 10;\t\t\t#[atm] - Pressure\n",
      "P = P*101325;\t\t\t#[Pa]\n",
      "# Methane obeys the equation of state\n",
      "# Z = 1 + (P*B)/(R*T)\n",
      "\n",
      "# Calculations\n",
      "# At 35 C,\n",
      "B = -50;\t\t\t#[cm**(3)/mol]\n",
      "dB_dT = 1.0;\t\t\t#[cm**(3)/mol-K] - dB/dT\n",
      "dB_dT = dB_dT*10**(-6);\t\t\t#[m**(3)/mol-K]\n",
      "d2B_dT2 = -0.01;\t\t\t#[cm**(3)/mol-K**(2)] - d**2(B)/d(T**2)\n",
      "d2B_dT2 = d2B_dT2*10**(-6);\t\t\t#[m**(3)/mol-K**(2)]\n",
      "\n",
      "# Ideal gas molar heat capacity of methane is given by\n",
      "# Cp_0 = 4.75 + 1.2*10**(-2)*T + 0.303*10**(-5)*T**(2) - 2.63*10**(-9)*T**(3)\n",
      "\n",
      "# The molar heat capacity is given by\n",
      "# Cp = Cp_0 + Cp_R\n",
      "# For virial gas equation of state \n",
      "Cp_R = -P*T*d2B_dT2;\t\t\t#[J/mol-K]\n",
      "\n",
      "# thus heat capacity is given by \n",
      "# Cp = a + b*T + c*T**(2) + d*T**(3) - P*T*d2B_dT2\n",
      "# Putting the values, we get\n",
      "Cp = (4.75 + 1.2*10**(-2)*T + 0.303*10**(-5)*T**(2) - 2.63*10**(-9)*T**(3))*4.184 - P*T*d2B_dT2;\t\t\t#[J/mol-K]\n",
      "\n",
      "# Results\n",
      "print \" The molar heat capacity of methane is %f J/mol-K\"%(Cp);\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The molar heat capacity of methane is 39.349753 J/mol-K\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.5  Page Number : 338"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of final temperature after expansion\n",
      "\n",
      "from scipy.optimize import fsolve \n",
      "\n",
      "# Variables\n",
      "T_1 = 360;\t\t\t#[K] - Initial temperature\n",
      "P_1 = 10;\t\t\t#[bar] - Initial pressure\n",
      "P_1 = P_1*10**(5);\t\t\t#[Pa]\n",
      "Tc = 408.1;\t\t\t#[K] - Critical temperature\n",
      "Pc = 36.48;\t\t\t#[bar] - Critical pressure\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.181;\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "Cv_0 = 106.0;\t\t\t#[J/mol-K]\n",
      "\n",
      "# Calculations\n",
      "# At 360 K\n",
      "Tr = T_1/Tc;\t\t\t# Reduced temperature\n",
      "B_0 = 0.083-(0.422/(Tr)**(1.6));\n",
      "B_1 = 0.139-(0.172/(Tr)**(4.2));\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B = ((B_0 + (w*B_1))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol] - Second virial coefficient\n",
      "dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);\t\t\t# (dB_0/dT)\n",
      "dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);\t\t\t# (dB_1/dT)\n",
      "dB_dT = ((R*Tc)/(Pc))*(dB0_dT + w*dB1_dT);\t\t\t# dB/dT\n",
      "\n",
      "# At state 1\n",
      "V_1 = B + (R*T_1)/P_1;\t\t\t#[m**(3)/mol] - Molar volume\n",
      "# At state 1\n",
      "V_2 = 10*V_1;\t\t\t#[m**(3)/mol] - Molar volume\n",
      "\n",
      "T_2 = -(P_1*T_1*(dB_dT))/Cv_0 + T_1;\t\t\t#[K]\n",
      "\n",
      "T = 350;\n",
      "\n",
      "fault = 10;\n",
      "def f(T): \n",
      "    return  106*(T-T_1)+972.72-((R*T**(2))/(V_2-B_prime))*dB_dT_prime\n",
      "while(fault>0.007):\n",
      "        Tr_prime = T/Tc;\t\t\t# Reduced temperature\n",
      "        B_0_prime = 0.083-(0.422/(Tr_prime)**(1.6));\n",
      "        B_1_prime = 0.139-(0.172/(Tr_prime)**(4.2));\n",
      "            #We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "        B_prime = ((B_0_prime + (w*B_1_prime))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol] - Second virial coefficient\n",
      "        dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);\t\t\t# (dB_0/dT)\n",
      "        dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);\t\t\t# (dB_1/dT)\n",
      "        dB_dT_prime = ((R*Tc)/(Pc))*(dB0_dT_prime + w*dB1_dT_prime);\t\t\t# dB/dT\n",
      "        T_prime = fsolve(f,0.15)\n",
      "        fault=abs(T-T_prime);\n",
      "        T = T + 0.001;\n",
      "\n",
      "# Results\n",
      "print \" The final temperature is %f K\"%(T);\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The final temperature is 351.910000 K\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.6  Page Number : 339"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of fugacity of liquid benzene\n",
      "\n",
      "# Variables\n",
      "T = 220 + 273.15;\t\t\t#[K] - Temperature\n",
      "Tc = 562.2;\t\t\t#[K] - Critical temperature\n",
      "Pc = 48.98;\t\t\t#[bar] - Critical pressure\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.210;\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "P_sat = 1912.86;\t\t\t#[kPa] - Saturation pressure at 220 C\n",
      "P_sat = P_sat*10**(3);\t\t\t#[Pa]\n",
      "Mol_wt = 78.114;\t\t\t#[g/mol] - Molecular weight of benzene\n",
      "\n",
      "# Calculations and Results\n",
      "#(1)\n",
      "# Since liquid and vapour are in equilibrium the fugacity is saturated fugacity (f_sat) and can be calculated using virial gas equation of state\n",
      "# At 220 C\n",
      "Tr = T/Tc;\t\t\t# Reduced temperature\n",
      "B_0 = 0.083-(0.422/(Tr)**(1.6));\n",
      "B_1 = 0.139-(0.172/(Tr)**(4.2));\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B = ((B_0 + (w*B_1))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol] - Second virial coefficient\n",
      "\n",
      "# We know that math.log(f/P) = (B*P)/(R*T)\n",
      "# Thus at saturated conditions\n",
      "# math.log(f_sat/P_sat) = B*P_sat/(R*T)\n",
      "f_sat = P_sat*(math.exp((B*P_sat)/(R*T)));\t\t\t#[Pa]\n",
      "f_sat = f_sat*10**(-3);\t\t\t#[kPa]\n",
      "\n",
      "print \" 1).The fugacity of liquid benzene is %f kPa\"%(f_sat);\n",
      "\n",
      "#(2)\n",
      "P = 2014.7;\t\t\t# [psia] - Total gauge pressure\n",
      "P = 138.94;\t\t\t# [bar]\n",
      "P = P*10**(5);\t\t\t# [Pa]\n",
      "den = 0.63;\t\t\t# [g/cm**(3)] - density of benzene\n",
      "den = den*10**(3);\t\t\t# [kg/m**(3)]\n",
      "\n",
      "# Therefore specific volume is \n",
      "V = 1/den;\t\t\t#[m/**(3)/kg]\n",
      "# Molar volume is given by\n",
      "V = V*Mol_wt*10**(-3);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# Thus fugacity at 220 C and pressure P is given by \n",
      "f = f_sat*(math.exp((V*(P-P_sat))/(R*T)));\n",
      "\n",
      "print \" 2).The fugacity of liquid benzene is %f kPa\"%(f);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The fugacity of liquid benzene is 1551.083216 kPa\n",
        " 2).The fugacity of liquid benzene is 2228.386532 kPa\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.7  Page Number : 341"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of molar enthalpy\n",
      "\n",
      "# Variables\n",
      "# C = -0.067 + 30.7/T\n",
      "# D = 0.0012 - 0.416/T\n",
      "\n",
      "T = 80 + 273.15;\t\t\t#[K]\n",
      "P = 30;\t\t\t#[bar]\n",
      "#P = P;\t\t\t#[N/m**(2)]\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "\n",
      "# Calculations\n",
      "# We have the relation derived in the book\n",
      "# d(G/(R*T)) = (V/(R*T))dP - (H/(R*T**(2)))dT\n",
      "# Writing the same equation for ideal gas and subtracting it from the above equation we get\n",
      "# d(G_R/(R*T)) = (V_R/(R*T))dP - (H_R/(R*T**(2)))dT\n",
      "# Therefore, H_R/(R*T**(2)) = -[del((G_R)/(R*T))/del(T)]_P\n",
      "\n",
      "# Substituting the relation G_R/(R*T) = math.log(f/P), we get\n",
      "# H_R/(R*T**(2)) = -[del(math.log(f/P))/del(T)]_P = -[del(-C*P - D*P**(2))/del(T)]_P\n",
      "# or, H_R/(R*T**(2)) = P*(dC/dT) + P**(2)*dD/dT\n",
      "# Note that in the above equation the partial derivative is replaced by full derivative as C and D are functions of temperature. Therfore we get\n",
      "# H_R/(R*T**(2)) = (30.7*P)/T**(2) + (0.416*P**(2))/T**(2)\n",
      "# H_R/R = - 30.7*P + 0.416*P**(2)\n",
      "\n",
      "# Substituting the given conditions we get\n",
      "H_R = R*(-30.7*P + 0.416*P**(2));\t\t\t#[J/mol]\n",
      "\n",
      "# Results\n",
      "print \" The molar enthalpy of the gas relative to that of the ideal gas at 80 C and 30 bar pressure is H_R = %f J/mol\"%(H_R);\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The molar enthalpy of the gas relative to that of the ideal gas at 80 C and 30 bar pressure is H_R = -4544.432400 J/mol\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.8  Page Number : 341"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Determination of second and third virial coefficients and fugacity\n",
      "\n",
      "# Variables\n",
      "# (1)\n",
      "T = 311;\t\t\t#[K] - Temperature\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "# Pressure in 'bar' is given below\n",
      "P = [0.690,1.380,2.760,5.520,8.280,11.034,13.800,16.550];\n",
      "# Molar volume in 'm**(3)/mol' is given below\n",
      "V = [0.0373,0.0186,0.00923,0.00455,0.00298,0.00220,0.00175,0.00144];\n",
      "\n",
      "# Z = 1 + (B/V) + (C/V**(2))\n",
      "# (Z-1)*V = B + (C/V)\n",
      "\n",
      "from numpy import zeros\n",
      "Z=zeros(8);\n",
      "k=zeros(8);\n",
      "t=zeros(8);\n",
      "\n",
      "# Calculations and Results\n",
      "for i in range(8):\n",
      "    Z[i]=(P[i]*10**(5)*V[i])/(R*T);\n",
      "    k[i]=(Z[i]-1)*V[i];\n",
      "    t[i]=1./V[i];\n",
      "\n",
      "from scipy.stats import linregress\n",
      "from numpy import array\n",
      "t = array(t)\n",
      "k = array(k)\n",
      "C,B,c,d,e = linregress(t.T,k.T)\n",
      "#[C,B,sig]=reglin(t',k');\n",
      "\n",
      "#From the regression, we get intercept = B and slope = C,and thus,\n",
      "\n",
      "print \" 1).The second virial coefficient of CO2 is given by B = %e m**3)/mol\"%(B);\n",
      "print \"     The thied virial coefficient of CO2 is given by C = %e m**6)/mol**2)\"%(C);\n",
      "\n",
      "# (2)\n",
      "P_final = 13.8;\t\t\t#[bar]\n",
      "\n",
      "def f40(P): \n",
      "\t return V-(R*T)/P\n",
      "\n",
      "# We know that R*T*math.log(f/P) =  quad(f40,0,P)[0]\n",
      "\n",
      "# Therefore we have to plot V - (R*T)/P versus P and calculate the area beneath the curve from 0 to 13.8 bar\n",
      "# For this we need the value of the term V - (R*T)/P at P = 0. At low pressure the virial equation becomes\n",
      "# Z = 1 + (B/V)\n",
      "#and V - (R*T)/P = (Z*R*T)/P - (R*T)/P = (1 + (B/V))*((R*T)/P) - (R*T)/P = (B*R*T)/(P*V) = (B/Z)\n",
      "# Thus lim P tending to zero (V - (R*T)/P) = B    ( as P tend to zero, Z tend to 1 ) \n",
      "\n",
      "P_prime = [0.000,0.690,1.380,2.760,5.520,8.280,11.034,13.800];\n",
      "V_prime = [0.000,0.0373,0.0186,0.00923,0.00455,0.00298,0.00220,0.00175];\n",
      "summation = 0;\n",
      "x=zeros(8);\n",
      "y=zeros(8);\n",
      "z=zeros(8);\n",
      "for j in range(2,8):\n",
      "    x[j]=V_prime[j]-(R*T)/(P_prime[j]*10**(5));\t\t\t#[m**(3)/mol]\n",
      "    y[j]=(x[j] + x[j-1])/2;\n",
      "    z[j]=y[j]*((P_prime[j]-P_prime[j-1]))*10**(5);\n",
      "    summation = summation + z[j] ;\t\t\t#[J/mol]\n",
      "\n",
      "summation = summation + 2*z[1] - z[1];\t\t\t# Because in the above calculation,in order to calculate the average a summation of z(2) is not included,only half of it gets added\n",
      "\n",
      "\n",
      "def f41(P): \n",
      "\t return V - (R*T)/P\n",
      "\n",
      "\t\t\t# Now we have, area =  quad(f41,0,13.8*10**(5))[0]\n",
      "\n",
      "\t\t\t# R*T*math.log(f/P) = summation\n",
      "f = P_final*(math.exp(summation/(R*T)));\t\t\t#[bar]\n",
      "\n",
      "\n",
      "print \" 2).The fugacity of steam at 311 K and 13.8 bar pressure is %f bar\"%(f);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The second virial coefficient of CO2 is given by B = -1.515264e-04 m**3)/mol\n",
        "     The thied virial coefficient of CO2 is given by C = 5.618543e-08 m**6)/mol**2)\n",
        " 2).The fugacity of steam at 311 K and 13.8 bar pressure is 12.892780 bar\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.9  Page Number : 344"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Determination of second and third virial coefficients\n",
      "\n",
      "# Variables\n",
      "T = 0 + 273.15;\t\t\t#[K] - Temperature\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\t\t\t# Pressure in 'atm' is given below\n",
      "P = [100,200,300,400,500,600,700,800,900,1000];\n",
      "\t\t\t# The compressibility factor values are\n",
      "Z = [1.069,1.138,1.209,1.283,1.356,1.431,1.504,1.577,1.649,1.720];\n",
      "\n",
      "# Z = 1 + (B/V) + (C/V**(2))\n",
      "# (Z-1)*V = B + (C/V)\n",
      "from numpy import zeros\n",
      "\n",
      "V = zeros(10);\n",
      "k = zeros(10);\n",
      "t = zeros(10);\n",
      "\n",
      "# Calculations and Results\n",
      "for i in range(10):\n",
      "    V[i]=Z[i]*R*T/(P[i]*101325);\t\t\t#[m**(3)/mol]\n",
      "    k[i]=(Z[i]-1)*V[i];\n",
      "    t[i]=1/V[i];\n",
      "\n",
      "from scipy.stats import linregress\n",
      "C,B,c,d,e = linregress(t,k)\n",
      "\t\t\t#[C,B,sig]=reglin(t,k);\n",
      "\n",
      "\t\t\t#From the regression, we get intercept = B and slope = C,and thus,\n",
      "\n",
      "print \" 1).The second virial coefficient of H2 is given by B = %e m**3)/mol\"%(B);\n",
      "print \"     The thied virial coefficient of H2 is given by C = %e m**6)/mol**2)\"%(C);\n",
      "\n",
      "\t\t\t# (2)\n",
      "\t\t\t# We know that, limit P tending to zero (V-(R*T)/P) = B, therfore P = 0,  V-(R*T)/P = B\n",
      "\n",
      "def f42(P): \n",
      "\t return V-(R*T)/P \t\t\t#and determine the integral integrate((V-(R*T)/P)\n",
      "\n",
      "\t\t\t# Now let us tabulate V-(R*T)/P and determine the integral  quad(f42,0,1000)[0]\n",
      "\n",
      "\n",
      "P_prime = [0,100,200,300,400,500,600,700,800,900,1000];\n",
      "Z_prime = [0,1.069,1.138,1.209,1.283,1.356,1.431,1.504,1.577,1.649,1.720];\n",
      "\n",
      "summation = 0;\n",
      "V_prime = zeros(11);\n",
      "x = zeros(11);\n",
      "y = zeros(11);\n",
      "z = zeros(11);\n",
      "for j in range(1,11):\n",
      "    V_prime[j]=Z_prime[j]*R*T/(P_prime[j]*101325);\t\t\t#[m**(3)/mol]\n",
      "    x[j]=V_prime[j]-(R*T)/(P_prime[j]*101325);\n",
      "    y[j]=(x[j] + x[j-1])/2;\n",
      "    z[j]=y[j]*((P_prime[j]-P_prime[j-1]))*101325;\n",
      "    summation = summation + z[j] ;\t\t\t#[J/mol]\n",
      "\n",
      "summation = summation + 2*z[1] - z[1];\t\t\t# Because in the above calculation,in order to calculate the average a summation of z(2) is not included,only half of it gets added\n",
      "\n",
      "\t\t\t# Now we have\n",
      "\t\t\t# R*T*math.log(f/P) = summation\n",
      "P_dash = 1000;\t\t\t#[atm] - Pressure at which fugacity is to be calculated\n",
      "T_dash = 273.15;\t\t\t#[K] - Temperature at which fugacity is to be calculated\n",
      "f = P_dash*math.exp(summation/(R*T_dash));\t\t\t#[atm] \n",
      "\n",
      "print \" 2).The fugacity of H2 at 0 C and 1000 atm pressure is  f = %f atm\"%(f);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The second virial coefficient of H2 is given by B = 1.345242e-05 m**3)/mol\n",
        "     The thied virial coefficient of H2 is given by C = 5.286525e-10 m**6)/mol**2)\n",
        " 2).The fugacity of H2 at 0 C and 1000 atm pressure is  f = 2030.305169 atm\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.10  Page Number : 345"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Determination of work done and the exit temperature\n",
      "\n",
      "from scipy.optimize import fsolve \n",
      "import math\n",
      "from scipy.integrate import quad \n",
      "\n",
      "# Variables\n",
      "P_1 = 1*10.**(6.);\t\t\t#[Pa] - Initial pressure\n",
      "T_1 = 200. + 273.15;\t\t#[K] - Initial temperature\n",
      "P_2 = 8*10.**(6);\t\t\t#[Pa] - Final pressure\n",
      "Tc = 647.1;\t\t\t        #[K] - Critical temperature of water\n",
      "Pc = 220.55;    \t\t\t#[bar] - Critical pressure of water\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.345;\n",
      "R = 8.314;\t\t    \t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "# For the virial gas the following are the relations for residual enthalpy and entropy\n",
      "# H_R = B*P - P*T*(dB/dT) \n",
      "# S_R = -P*(dB/dT)\n",
      "# Where, (dB/dT) = ((R*Tc)/Pc)*((dB_0/dT) + w*(dB_1/dT))\n",
      "# dB0_dT = 0.422*1.6*Tc**(1.6)*T**(-2.6);\t\t\t# (dB_0/dT)\n",
      "# dB1_dT = 0.172*4.2*Tc**(4.2)*T**(-5.2);\t\t\t# (dB_1/dT)\n",
      "\n",
      "# (1)\n",
      "Cp_0 = 29.114;\t\t\t#[J/mol-K] - Specific heat capacity at constant pressure\n",
      "# For the isentropic process entropy change is zero, thus\n",
      "# delta_S = Cp_0*math.log(T_2/T_1) - P_2*(dB/dT)_2 + P_1*(dB/dT)_1 = 0\n",
      "\n",
      "# Calculations\n",
      "# At state 1, \n",
      "Tr_1 = T_1/Tc;\n",
      "B0_1 = 0.083 - 0.422/(Tr_1**(1.6));\n",
      "B1_1 = 0.139 - 0.172/(Tr_1**(4.2));\n",
      "# (B*Pc)/(R*Tc) = B0 + w*B1\n",
      "B_1 = ((B0_1 + w*B1_1)*(R*Tc))/Pc;\t\t\t# [m**(3)/mol] - Second virial coefficient at state 1\n",
      "dB0_dT_1 = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);\t\t\t# (dB_0/dT)\n",
      "dB1_dT_1 = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);\t\t\t# (dB_1/dT)\n",
      "dB_dT_1 = ((R*Tc)/Pc)*((dB0_dT_1) + w*(dB1_dT_1));\t\t\t# (dB/dT)_1\n",
      "\n",
      "# Now let us assume the exit temperature to be 870 K, at this temperature\n",
      "# T_2 = 870;\t\t\t#[K] - \n",
      "# At this temperature\n",
      "# delta_S = Cp_0*math.log(T_2/T_1) - P_2*(dB/dT)_2 + P_1*(dB/dT)_1 = \n",
      "\n",
      "\n",
      "T_2 = 860.;\t\t\t#[K] - Exit temperature\n",
      "# Therefore at state 2, we have\n",
      "Tr_2 = T_2/Tc;\n",
      "B0_2 = 0.083 - 0.422/(Tr_2**(1.6));\n",
      "B1_2 = 0.139 - 0.172/(Tr_2**(4.2));\n",
      "# (B*Pc)/(R*Tc) = B0 + w*B1\n",
      "B_2 = ((B0_2 + w*B1_2)*(R*Tc))/Pc;\t\t\t# [m**(3)/mol] - Second virial coefficient at state 2\n",
      "dB0_dT_2 = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);\t\t\t# (dB_0/dT)\n",
      "dB1_dT_2 = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);\t\t\t# (dB_1/dT)\n",
      "dB_dT_2 = ((R*Tc)/Pc)*((dB0_dT_2) + w*(dB1_dT_2));\t\t\t# (dB/dT)_2\n",
      "\n",
      "delta_H_s = Cp_0*(T_2 - T_1) + B_2*P_2 -P_2*T_2*(dB_dT_2) - B_1*P_1 + P_1*T_1*(dB_dT_1);\t\t\t#[J/mol] - Enthalpy change\n",
      "\n",
      "# As no heat exchange is assumed to take place with the surroundings,work transfer is given by\n",
      "W_1 = - delta_H_s;\t\t\t# [J/mol]\n",
      "\n",
      "# Results\n",
      "print \" 1).The exit temperature is %f K\"%(T_2);\n",
      "print \"     The required amount of work is %f J/mol\"%(W_1);\n",
      "\n",
      "\n",
      "# (2)\n",
      "eff = 0.8;\t\t\t# Adiabatic efficiency\n",
      "delta_H_a = delta_H_s/0.8;\t\t\t# Actual enthalpy change\n",
      "\n",
      "# Now for calculating the value of T_exit\n",
      "# delta_H_a = Cp_0*(T_exit - T_1) + B*P_2 -P_2*T_exit*(dB_dT) - B_1*P_1 + P_1*T_1*(dB_dT_1)\n",
      "# On simplification we get\n",
      "# 29.114*(T_2 - T_1)*B_2*8*10**(6)-8*10**(6)*T_2*(dB/dT)_2 = 12643.77 \n",
      "\n",
      "# Let us assume a temperature of say\n",
      "T = 900.;\t\t\t#[K]\n",
      "fault=10.;\n",
      "\n",
      "def f(T_exit): \n",
      "\t return  delta_H_a - Cp_0*(T_exit - T_1) + B*P_2 -P_2*T_exit*(dB_dT) - B_1*P_1 + P_1*T_1*(dB_dT_1)\n",
      "\t \n",
      "while(fault>0.3):\n",
      "    Tr = T/Tc;\n",
      "    B0 = 0.083 - 0.422/(Tr**(1.6));\n",
      "    B1 = 0.139 - 0.172/(Tr**(4.2));\n",
      "    \t\t\t# (B*Pc)/(R*Tc) = B0 + w*B1\n",
      "    B = ((B0 + w*B1)*(R*Tc))/Pc;\t\t\t# [m**(3)/mol] - Second virial coefficient at state 2\n",
      "    dB0_dT = 0.422*1.6*Tc**(1.6)*T**(-2.6);\t\t\t# (dB_0/dT)\n",
      "    dB1_dT = 0.172*4.2*Tc**(4.2)*T**(-5.2);\t\t\t# (dB_1/dT)\n",
      "    dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT));\t\t\t# (dB/dT)_1\n",
      "\n",
      "    T_exit = fsolve(f,900)\n",
      "    fault=abs(T-T_exit);\n",
      "    T = T + 0.2;\n",
      "\n",
      "Texit = T;\n",
      "\n",
      "\t\t\t# As no heat exchange is assumed to take place with the surroundings,work transfer is given by\n",
      "W_2 = - delta_H_a;\t\t\t# [J/mol]\n",
      "\n",
      "print \" 2).The exit temperature is %f K\"%(Texit);\n",
      "print \"     The required amount of work is %f J/mol\"%(W_2);\n",
      "\n",
      "\n",
      "def f20(T): \n",
      "\t return Cp_0/T\n",
      "\n",
      "T_prime = 700.;\t\t\t#[K]\n",
      "fault1=10.;\n",
      "\n",
      "def f1(T_out): \n",
      "\t return  32.2168*math.log(T_out) + 0.1922*10**(-2)*T_out + 0.5274*10**(-5) \\\n",
      "\t *T_2**(2) - 1.1976*10**(-9)*T_out**(3)-8*10**(6)*dB_dT_prime -216.64\n",
      "\n",
      "while(fault1>0.5):\n",
      "    Tr_prime = T_prime/Tc;\n",
      "    B0_prime = 0.083 - 0.422/(Tr_prime**(1.6));\n",
      "    B1_prime = 0.139 - 0.172/(Tr_prime**(4.2));\n",
      "    B_prime = ((B0_prime + w*B1_prime)*(R*Tc))/Pc;\t\t\t# [m**(3)/mol] - Second virial coefficient at state 2\n",
      "    dB0_dT_prime = 0.422*1.6*Tc**(1.6)*T_prime**(-2.6);\t\t\t# (dB_0/dT)\n",
      "    dB1_dT_prime = 0.172*4.2*Tc**(4.2)*T_prime**(-5.2);\t\t\t# (dB_1/dT)\n",
      "    dB_dT_prime = ((R*Tc)/Pc)*((dB0_dT_prime) + w*(dB1_dT_prime));\t\t\t# (dB/dT)_1\n",
      "    T_out = fsolve(f1,10)\n",
      "    fault1=abs(T_prime-T_out);\n",
      "    T_prime = T_prime + 0.5;\n",
      "\n",
      "T_out = T_prime;\n",
      "\n",
      "# Now we have to calculate enthalpy change as W = -delta_H\n",
      "\n",
      "def f21(T): \n",
      "\t return (7.7 + 0.04594*10**(-2.)*T + 0.2521*10**(-5.)*T**(2.) - 0.8587*10.**(-9.)*T**(3.))*4.184\n",
      "\n",
      "delta_H_3 =  quad(f21,T_1,T_out)[0] + B_prime*P_2 - P_2*T_out*dB_dT_prime - B_1*P_1 + P_1*T_1*dB_dT_1\n",
      "print T_1,T_out\n",
      "\n",
      "\n",
      "W_3 = - delta_H_3;\t\t\t# [J/mol]\n",
      "\n",
      "print \" 3).The exit temperature is %f K\"%(T_out);\n",
      "print \"     The required amount of work is %f J/mol\"%(W_3);\n",
      "\n",
      "#(4)\n",
      "n = 0.8;\t\t\t# Adiabatic efficiency\n",
      "delta_H_a_4 = delta_H_3/n;\t\t\t#[J/mol]\n",
      "W_4 = -delta_H_a_4;\t\t\t#[J/mol]\n",
      "\n",
      "\n",
      "T_prime1 = 700.;\t\t\t#[K]\n",
      "fault2=10.;\n",
      "\n",
      "def f2(T_2): \n",
      "\t return  7.7*4.184*(T_2-T_1) + ((0.04594*4.184*10**(-2))/2)* \\\n",
      "\t (T_2**(2)-T_1**(2)) + ((0.2521*4.184*10**(-5))/3)*(T_2**(3)-T_1**(3)) - \\\n",
      "\t ((0.8587*4.184*10**(-9))/4)*(T_2**(4)-T_1**(4)) + B_prime1*8*10**(6) \\\n",
      "\t - 8*10**(6)*T_2*dB_dT_prime1 + 191.7 + 496.81 - delta_H_a_4\n",
      "\n",
      "while(fault2>0.5):\n",
      "    Tr_prime1 = T_prime1/Tc;\n",
      "    B0_prime1 = 0.083 - 0.422/(Tr_prime1**(1.6));\n",
      "    B1_prime1 = 0.139 - 0.172/(Tr_prime1**(4.2));\n",
      "    \t\t\t# (B*Pc)/(R*Tc) = B0 + w*B1\n",
      "    B_prime1 = ((B0_prime1 + w*B1_prime1)*(R*Tc))/Pc;\t\t\t# [m**(3)/mol] - Second virial coefficient at state 2\n",
      "    dB0_dT_prime1 = 0.422*1.6*Tc**(1.6)*T_prime1**(-2.6);\t\t\t# (dB_0/dT)\n",
      "    dB1_dT_prime1 = 0.172*4.2*Tc**(4.2)*T_prime1**(-5.2);\t\t\t# (dB_1/dT)\n",
      "    dB_dT_prime1 = ((R*Tc)/Pc)*((dB0_dT_prime1) + w*(dB1_dT_prime1));\t\t\t# (dB/dT)_1\n",
      "    T_out1 = fsolve(f2,100)\n",
      "    fault2=abs(T_prime1-T_out1);\n",
      "    T_prime1 = T_prime1 + 0.5;\n",
      "\n",
      "T_out1 = T_prime1;\n",
      "\n",
      "print \" 4).The exit temperature is %f K\"%(T_out1);\n",
      "print \"     The required amount of work is %f J/mol\"%(W_4);\n",
      "\n",
      "\n",
      "# answers are varies because of rounding error. please check it manually."
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The exit temperature is 860.000000 K\n",
        "     The required amount of work is -10667.724150 J/mol\n",
        " 2).The exit temperature is 916.600000 K\n",
        "     The required amount of work is -13334.655188 J/mol\n",
        "473.15"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 755.5\n",
        " 3).The exit temperature is 755.500000 K\n",
        "     The required amount of work is -9282.701079 J/mol\n",
        " 4).The exit temperature is 809.500000 K"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "     The required amount of work is -11603.376349 J/mol\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.11  Page Number : 348"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of temperature and pressure\n",
      "\n",
      "# Variables\n",
      "Vol = 0.15;\t\t\t    #[m**(3)] - Volume of the cylinder\n",
      "P_1 = 100.;\t\t\t    #[bar] - Initial pressure\n",
      "P_1 = P_1*10**(5);\t\t#[Pa]\n",
      "T_1 = 170.;\t\t\t    #[K] - Initial temperature\n",
      "n_withdrawn = 500.;\t\t#[mol] - Withdrawn number of moles\n",
      "R = 8.314;\t\t    \t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "\n",
      "# Calculations and Results\n",
      "#(1)\n",
      "Y = 1.4;\t\t\t# Coefficient of adiabatic expansion\n",
      "n_total = (P_1*Vol)/(R*T_1);\t\t\t#[mol] - Total number of moles\n",
      "n_2 = n_total - n_withdrawn;\t\t\t#[mol] - Left number of moles\n",
      "V_1 = Vol/n_total;\t\t\t#[m**(3)/mol] - Molar volume at initial state.\n",
      "# At final state\n",
      "V_2 = Vol/n_2;\t\t\t#[m**(3)/mol] - Molar volume at final state\n",
      "\n",
      "# During duscharging  P_1*V_1**(Y) = P_2*V_2**(Y), therefore\n",
      "P_2_1 = P_1*((V_1/V_2)**(Y));\t\t\t#[Pa] - Final pressure\n",
      "P_2_1 = P_2_1*10**(-5);\t\t\t#[bar]\n",
      "T_2_1 = ((P_2_1*10**(5))*V_2)/R;\t\t\t#[K] - Final temperature\n",
      "\n",
      "print \" 1).The final temperature %f K\"%(T_2_1);\n",
      "print \"     The final pressure %f bar\"%(P_2_1);\n",
      "\n",
      "\n",
      "def f53(T): \n",
      "\t return Cp_0/T\n",
      "\n",
      "T_prime = 150;\t\t\t#[K]\n",
      "error = 10;\n",
      "while(error>1):\n",
      "    f_T = 18.886*math.log(T_prime) + 4.2*10**(-3)*T_prime - 92.4;\n",
      "    f_T_dash = 18.886/T_prime + 4.2*10**(-3);\n",
      "    T_new = T_prime - (f_T/f_T_dash);\n",
      "    error=abs(T_prime - T_new);\n",
      "    T_prime = T_new;\n",
      "\n",
      "T_2_2 = T_prime;\t\t\t#[K] - Final temperature\n",
      "P_2_2 = ((n_2*R*T_2_2)/Vol)*10**(-5);\t\t\t#[bar] - Final pressure\n",
      "\n",
      "print \" 2).The final temperature %f K\"%(T_2_2);\n",
      "print \"     The final pressure %f bar\"%(P_2_2);\n",
      "\n",
      "#(3)\n",
      "Tc = 126.2;\t\t\t#[K] - Critical temperature of nitrogen\n",
      "Pc = 34.0;\t\t\t#[bar] - Critical pressure of nitrogen\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.038;\t\t\t# Acentric factor\n",
      "\n",
      "\n",
      "dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);\t\t\t# (dB_0/dT) at state 1\n",
      "dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);\t\t\t# (dB_1/dT) at state 1\n",
      "dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT));\t\t\t# (dB/dT) at state 1\n",
      "# The residual entropy at the initial state is given by \n",
      "S_R_1 = -P_1*(dB_dT);\t\t\t#[J/mol-K]\n",
      "\n",
      "# Now let us calculate molar volume at initial state\n",
      "Tr = T_1/Tc;\t\t\t# Reduced temperature\n",
      "B_0 = 0.083-(0.422/(Tr)**(1.6));\n",
      "B_1 = 0.139-(0.172/(Tr)**(4.2));\n",
      "\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B = ((B_0+(w*B_1))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "V_1_3 = B + (R*T_1)/P_1;\t\t\t#[m**(3)/mol]\n",
      "# Therefore number of moles in the initial state is\n",
      "n_1_3 = Vol/V_1_3;\t\t\t#[mol]\n",
      "# Therefore final number of moles is\n",
      "n_2_3 = n_1_3 - n_withdrawn;\n",
      "\n",
      "# Therefore molar volume at final state is\n",
      "V_2_3 = Vol/n_2_3;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# Now let us determine the relation between pressure and temperature in the final state\n",
      "# P_2_3 = (R*T_2_3)/(V_2_3 - B_2)\n",
      "#delta_S = 0, thus delta_S_ig + delta_S_R = 0\n",
      "delta_S_R = - S_R_1;\n",
      "\n",
      "def f54(T): \n",
      "\t return Cp_0/T\n",
      "\n",
      "\n",
      "T_2_prime = 135;\t\t\t#[K]\n",
      "delta = 0.1;\n",
      "error = 10;\n",
      "while(error>0.01):\n",
      "    T_r = T_2_prime/Tc;\t\t\t# Reduced temperature\n",
      "    B_0_3 = 0.083-(0.422/(T_r)**(1.6));\n",
      "    B_1_3 = 0.139-(0.172/(T_r)**(4.2));\n",
      "    B_3 = ((B_0_3+(w*B_1_3))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol]\n",
      "    dB0_dT_3 = 0.422*1.6*Tc**(1.6)*T_2_prime**(-2.6);\t\t\t# (dB_0/dT)\n",
      "    dB1_dT_3 = 0.172*4.2*Tc**(4.2)*T_2_prime**(-5.2);\t\t\t# (dB_1/dT)\n",
      "    dB_dT_3 = ((R*Tc)/Pc)*((dB0_dT_3) + w*(dB1_dT_3));\t\t\t# (dB/dT)\n",
      "    P_2_3 = (R*T_2_prime)/(V_2_3 - B_3);\n",
      "    delta_S = 27.2*(math.log(T_2_prime/T_1)) + 4.2*10**(-3)*(T_2_prime - T_1) - R*(math.log(P_2_3/P_1)) - P_2_3*(dB_dT_3) + delta_S_R;\n",
      "    T_new = T_2_prime + delta;\n",
      "    error=abs(delta_S);\n",
      "    T_2_prime = T_new;\n",
      "\n",
      "T_2_3 = T_2_prime;\t\t\t#[K] - Final temperature\n",
      "\t\t\t# Therefore at T_2_3\n",
      "P_2_3 = P_2_3*10**(-5);\t\t\t#[bar] - Final pressure\n",
      "\n",
      "print \" 3).The final temperature %f K\"%(T_2_3);\n",
      "print \"     The final pressure %f bar\"%(P_2_3);\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The final temperature 131.761821 K\n",
        "     The final pressure 40.991361 bar\n",
        " 2).The final temperature 129.504151 K\n",
        "     The final pressure 40.288995 bar\n",
        " 3).The final temperature 141.300000 K\n",
        "     The final pressure 57.079997 bar\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.12  Page Number : 351"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of change of internal energy enthalpy entropy and exergy\n",
      "\n",
      "# Variables\n",
      "P_1 = 80.;\t\t\t        #[bar] - Initial pressure\n",
      "P_1 = P_1*10.**(5);\t\t\t#[Pa]\n",
      "T_1 = 300. + 273.15;\t\t#[T] - Initial temperature\n",
      "P_2 = 40.;\t\t\t        #[bar] - Final pressure\n",
      "P_2 = P_2*10**(5);\t\t\t#[Pa]\n",
      "T_2 = 300. + 273.15;\t\t#[K] - Final temperature\n",
      "T_0 = 25. + 273.15;\t\t\t#[K] - Surrounding temperature\n",
      "P_0 = 1.;\t\t\t        #[atm] - Surrounding pressure\n",
      "P_0 = P_0*101325;\t\t\t#[Pa]\n",
      "Tc = 647.1;\t\t\t        #[K]\n",
      "Pc = 220.55;\t\t\t    #[bar]\n",
      "Pc = Pc*10.**(5);\t\t\t#[Pa]\n",
      "R = 8.314;\t\t\t        #[J/mol*K] - Universal gas constant\n",
      "\n",
      "\n",
      "# Calculations\n",
      "# For van der Walls equation of state\n",
      "a = (27*R**(2)*Tc**(2))/(64*Pc);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = (R*Tc)/(8*Pc);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "\n",
      "def f(V): \n",
      "\t return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1\n",
      "V_1_1=fsolve(f,0.1)\n",
      "V_1_2=fsolve(f,10)\n",
      "V_1_2=fsolve(f,100)\n",
      "# The largest root is considered because of vapour\n",
      "V_1 = V_1_1;\n",
      "\n",
      "U_R_1 = -a/V_1;\t\t\t#[J/mol] - Internal energy\n",
      "H_R_1 = P_1*V_1 - R*T_1 - a/V_1;\t\t\t#[J/mol] - Enthalpy\n",
      "S_R_1 = R*math.log((P_1*(V_1-b))/(R*T_1));\n",
      "\n",
      "def f1(V): \n",
      "\t return V**(3)-(b+(R*T_2)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2\n",
      "V_2_1 = fsolve(f1,0.1)\n",
      "V_2_2 = fsolve(f1,10)\n",
      "V_2_3 = fsolve(f1,100)\n",
      "V_2 = V_2_1;\n",
      "\n",
      "U_R_2 = -a/V_2;\t\t\t#[J/mol] - Internal energy\n",
      "H_R_2 = P_2*V_2 - R*T_2 - a/V_2;\t\t\t#[J/mol] - Enthalpy\n",
      "S_R_2 = R*math.log((P_2*(V_2-b))/(R*T_2));\n",
      "\n",
      "delta_U_R = U_R_2 - U_R_1;\t\t\t#\n",
      "delta_H_R = H_R_2 - H_R_1;\t\t\t#\n",
      "delta_S_R = S_R_2 - S_R_1;\t\t\t#\n",
      "\n",
      "delta_U_ig = 0;\t\t\t#[J/mol] - As temperature is constant\n",
      "delta_H_ig = 0;\t\t\t#[J/mol] - As temperature is constant\n",
      "delta_S_ig = - R*math.log(P_2/P_1);\t\t\t# [J/mol-K]\n",
      "delta_U = delta_U_R + delta_U_ig;\t\t\t#[J/mol]\n",
      "delta_H = delta_H_R + delta_H_ig;\t\t\t#[J/mol]\n",
      "delta_S = delta_S_R + delta_S_ig;\t\t\t#[J/mol-K]\n",
      "\n",
      "\t\t\t# Change in exergy is given by\n",
      "\t\t\t# delta_phi = phi_1 - phi_2 = U_1 - U_2 + P_0*(V_1 - _V_2) - T_0*(S_1 - S_2)\n",
      "delta_phi = - delta_U + P_0*(V_1 - V_2) - T_0*(-delta_S);\t\t\t#[J/mol]\n",
      "\n",
      "# Results\n",
      "print \" The change in internal energy is %f J/mol\"%(delta_U);\n",
      "print \" The change in enthalpy is %f J/mol\"%(delta_H);\n",
      "print \" The change in entropy is %f J/mol-K\"%(delta_S);\n",
      "print \" The change in exergy is %f J/mol\"%(delta_phi);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The change in internal energy is 615.070900 J/mol\n",
        " The change in enthalpy is 1053.220316 J/mol\n",
        " The change in entropy is 6.930264 J/mol-K\n",
        " The change in exergy is 1389.940735 J/mol\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.13  Page Number : 353"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of change in enthalpy\n",
      "\n",
      "# Variables\n",
      "T_1 = 500.;\t\t\t#[K] - Initial temperature\n",
      "P_1 = 30.;\t\t\t#[atm] - Initial pressure\n",
      "P_1 = P_1*101325;\t\t\t#[Pa]\n",
      "P_2 = 1;\t\t\t#[atm] - Final pressure\n",
      "P_2 = P_2*101325;\t\t\t#[Pa]\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\t\t\t# For chlorine\n",
      "Tc = 417.2;\t\t\t#[K] - Critical temperature\n",
      "Pc = 77.10;\t\t\t#[bar] - Critical pressure\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "\n",
      "#Redlich Kwong equation of state,\n",
      "a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc;\t\t\t# [Pa*m**(6)*K**(1/2)/mol]\n",
      "b = (0.08664*R*Tc)/Pc;\t\t\t# [m**(3)/mol]\n",
      "\n",
      "# Calculations\n",
      "def f1(V): \n",
      "\t return V**(3)-((R*T_1)/P_1)*V**(2)-((b**(2))+((b*R*T_1)/P_1)-(a/(T_1**(1./2)*P_1)))*V-(a*b)/(T_1**(1./2)*P_1)\n",
      "V_1=fsolve(f1,1)\n",
      "V_2=fsolve(f1,10)\n",
      "V_3=fsolve(f1,100)\n",
      "\n",
      "V = V_1;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "\n",
      "Z = (P_1*V_1)/(R*T_1);\t\t\t#compressibility factor\n",
      "\n",
      "# The residual enthalpy at state 1 is given by\n",
      "H_R_1 = (Z-1)*R*T_1 + ((3*a)/(2*b*T_1**(1./2)))*(math.log(V/(V+b)));\t\t\t#[J/mol]\n",
      " \n",
      "# Since chlorine is assumed to behave ideally under the final condition,therefore\n",
      "H_R_2 = 0;\t\t\t# Residual enthalpy at state 2\n",
      "delta_H_R = H_R_2 - H_R_1;\t\t\t#[J/mol] - Residual enthalpy change\n",
      "# and since isothermal conditions are maintained, therfore\n",
      "delta_H_ig = 0;\t\t\t# Enthalpy change under ideal condition\n",
      "delta_H = delta_H_R + delta_H_ig;\t\t\t#[J/mol]\n",
      "\n",
      "# Results\n",
      "print \" The change in enthalpy is given by delta_H = %f J/mol\"%(delta_H);\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The change in enthalpy is given by delta_H = 1053.558471 J/mol\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.14  Page Number : 353"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of final temperature\n",
      "\n",
      "#(1)\n",
      "#This part involves proving a relation in which no mathematical components are involved.\n",
      "#For proof refer to this example 10.14 on page number 354 of the book.\n",
      "\n",
      "\n",
      "#(2)\n",
      "\n",
      "# Variables\n",
      "Vol_1 = 0.1;\t\t\t#[m**(3)] - Initial volume of each compartment\n",
      "n_1 = 400;\t\t\t#[mol] - Initial number of moles in compartment 1\n",
      "V_1 = Vol_1/n_1;\t\t\t#[m**(3)/mol] - Molar volume at state 1\n",
      "T_1 = 294;\t\t\t#[K]\n",
      "Vol_2 = 0.2;\t\t\t#[m**(3)] - Final volume of the compartment after removing the partition.\n",
      "n_2 = n_1;\t\t\t#[mol] - Number of moles remains the same\n",
      "V_2 = Vol_2/n_2;\t\t\t#[m**(3)/mol] - Molar volume at state 2\n",
      "\n",
      "a = 0.1362;\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = 3.215*10**(-5);\t\t\t#[m**(3)/mol]\n",
      "Cv_0 = 12.56;\t\t\t#[J/mol-K] - Heat capacity in ideal gas state\n",
      "\n",
      "# Calculations\n",
      "# For overall system q = 0, and no work is done, therefore delta_U = 0\n",
      "# Therfore from the relation proved in part (1), we have\n",
      "T_2 = T_1 + (a/Cv_0)*(1/V_2 - 1/V_1);\t\t\t#[K]\n",
      "\n",
      "# Results\n",
      "print \" 2).The final temperatutre is %f K\"%(T_2)\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 2).The final temperatutre is 272.312102 K\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.17  Page Number : 356"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Determination of work done and the exit temperature\n",
      "\n",
      "# Variables\n",
      "P_1 = 1*10**(6);\t\t\t#[Pa] - Initial pressure\n",
      "T_1 = 200 + 273.15;\t\t\t#[K] - Initial temperature\n",
      "P_2 = 8*10**(6);\t\t\t#[Pa]\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "Y = 1.4;\t\t\t# Index of expansion\n",
      "Cp_0 = 29.114;\t\t\t#[J/mol-K]\n",
      "\t\n",
      "a = 0.55366;\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = 3.049*10**(-5);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# Calculations\n",
      "V_1 = 3.816*10**(-3);\t\t\t#[m**(3)/mol]\n",
      "Z_1 = (P_1*V_1)/(R*T_1);\n",
      "\n",
      "T_2 = T_1*(P_2/P_1)**((Y-1)/Y);\t\t\t#[K]\n",
      "\n",
      "# At 8 MPa and T_2,\n",
      "# The molar volume of steam following van der Walls equation of state (as reported in the book) is\n",
      "V_2 = 8.41*10**(-4);\t\t\t#[m**(3)/mol]\n",
      "# And the compressibility factor is \n",
      "Z_2 = (P_2*V_2)/(R*T_2);\n",
      "\n",
      "# For van der Walls equation of state we know that\n",
      "# delta_S_R/R = math.log(Z_2/Z_1) + math.log((V_2 - b)/V_2) - math.log((V_1 - b)/V_1)\n",
      "delta_S_R = R*(math.log(Z_2/Z_1) + math.log((V_2 - b)/V_2) - math.log((V_1 - b)/V_1));\t\t\t#[J/mol]\n",
      "\n",
      "# delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1)\n",
      "# The entropy change is therefore\n",
      "# delta_S = delta_S_ig + delta_S_R\n",
      "# But during an isentropic process the total entropy change is zero\n",
      "# Therefore we have to modify the exit temperature so that the entropy change is zero\n",
      "\n",
      "# Let us assume a temperature, say T = 870 K\n",
      "# At 870 K the molar volume of steam following van der Walls equation of state (as reported in the book) is\n",
      "#  V_3 = 8.57*10**(-4);\t\t\t# [m**(3)/mol]\n",
      "# Therefore\n",
      "# Z_3 = (P_2*V_3)/(R*T_2);\n",
      "# At this temperature,\n",
      "# delta_S = Cp_0*math.log(T/T_1) - R*math.log(P_2/P_1) + R*(math.log(Z/Z_1) + R*math.log((V - b)/V) - R*math.log((V_1 - b)/V_1))\n",
      "\n",
      "T = 800;\t\t\t#[K]\n",
      "fault=10;\n",
      "\n",
      "def f1(V): \n",
      "\t return V**(3)-(b+(R*T)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2\n",
      "\t \n",
      "def f2(T): \n",
      "\t return  Cp_0*math.log(T/T_1) - R*math.log(P_2/P_1) + R*(math.log(Z/Z_1) + R*(math.log((V - b)/V)) - R*(math.log((V_1 - b)/V_1)))\n",
      "\n",
      "\n",
      "while(fault>0.3):\n",
      "    \t\t\t# At T and 8 MPa\n",
      "    V = fsolve(f1,1)\n",
      "    Z = (P_2*V)/(R*T);\n",
      "    \n",
      "    T_exit = fsolve(f2,0.1)\n",
      "    fault=abs(T-T_exit);\n",
      "    T = T + 0.5;\n",
      "\n",
      "Texit = T;\n",
      "\n",
      "# Now applying the first law to an adiabatic process we get\n",
      "# W = - delta_H\n",
      "\n",
      "# For van der Walls gas the enthalpy change is given by\n",
      "delta_H_s = Cp_0*(T_exit - T_1) + (Z - 1)*R*T_exit - a/V - (Z_1-1)*R*T_1 + a/V_1;\t\t\t#[J/mol]\n",
      "W = - delta_H_s;\t\t\t#[J/mol]\n",
      "\n",
      "# Results\n",
      "print \" 1).The exit temperature is %f K\"%(Texit);\n",
      "print \"     The work required is given by W = %f J/mol\"%(W);\n",
      "\n",
      "#(2)\n",
      "eff = 0.8;\t\t\t# Adiabatic efficiency\n",
      "delta_H_a = eff*delta_H_s;\t\t\t#[J/mol] - Actual enthalpy change\n",
      "W_2 = - delta_H_a;\n",
      "\n",
      "\t\t\t# Let us assume a temperature, say\n",
      "T_prime= 900;\t\t\t#[K]\n",
      "fault1=10;\n",
      "def f22(V): \n",
      "\t return V**(3)-(b+(R*T_prime)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2\n",
      "\n",
      "def f3(T_prime): \n",
      "\t return  Cp_0*(T_prime - T_1) + (Z_prime - 1)*R*T_prime - a/V_prime - 13230.49\n",
      "\t \n",
      "while(fault1>0.3):\n",
      "    \t\t\t# At T_prime and 8 MPa\n",
      "    V_prime=fsolve(f22,1)\n",
      "    Z_prime = (P_2*V_prime)/(R*T_prime);\n",
      "    T_exit1 = fsolve(f3,100)\n",
      "    fault1=abs(T_prime-T_exit1);\n",
      "    T_prime = T_prime + 0.2;\n",
      "\n",
      "Texit1 = T_prime;\n",
      "\n",
      "print \" 2).The exit temperature is %f K\"%(Texit1);\n",
      "print \"     The work required is given by W = %f J/mol\"%(W_2);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The exit temperature is 916.500000 K\n",
        "     The work required is given by W = -12195.093996 J/mol\n",
        " 2).The exit temperature is 958.400000 K"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "     The work required is given by W = -9756.075197 J/mol\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.19  Page Number : 359"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of molar volume and fugacity\n",
      "\n",
      "# Variables\n",
      "T = 100 + 273.15;\t\t\t#[K] - Temperature\n",
      "Tc = 647.1;\t\t\t#[K] - Critical temperature of water\n",
      "Pc = 220.55;\t\t\t#[bar] - Critical pressure of water\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "# Calculations\n",
      "a = (27*R**(2)*Tc**(2))/(64*Pc);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = (R*Tc)/(8*Pc);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# The cubic form of van der Walls equation of state is given by,\n",
      "# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0\n",
      "\n",
      "# For water vapour at 100 C under saturated conditions pressure is 1 atm, therefore\n",
      "P = 1;\t\t\t#[atm]\n",
      "P = P*101325;\t\t\t#[Pa]\n",
      "\n",
      "# At 100 C and 1 atm \n",
      "def f(V): \n",
      "\t return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P\n",
      "V_1 = fsolve(f,0.1)\n",
      "V_1 = fsolve(f,10)\n",
      "V_1 = fsolve(f,100)\n",
      "V = V_1;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "\n",
      "f = P*(math.exp(math.log((R*T)/(P*(V-b))) + b/(V-b) - (2*a)/(R*T*V)));\t\t\t#[Pa]\n",
      "f = f/101325;\t\t\t#[atm]\n",
      "\n",
      "# Results\n",
      "print \" The molar volume is %f m**3)/mol\"%(V);\n",
      "print \" The fugacity is %f atm\"%(f);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The molar volume is 0.030469 m**3)/mol\n",
        " The fugacity is 0.995168 atm\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.20  Page Number : 359"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of enthalpy and entropy change\n",
      "\n",
      "# Variables\n",
      "P_1 = 6;\t\t\t#[bar] - Initial pressure\n",
      "P_1 = P_1*10**(5);\t\t\t#[Pa]\n",
      "T_1 = 100 + 273.15;\t\t\t#[T] - Initial temperature\n",
      "P_2 = 12;\t\t\t#[bar] - Final pressure\n",
      "P_2 = P_2*10**(5);\t\t\t#[Pa]\n",
      "T_2 = 500 + 273.15;\t\t\t#[K] - Final temperature\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "Y = 1.126;\t\t\t# Index of expansion\n",
      "Cp_0 = (R*Y)/(Y-1);\t\t\t#[J/mol-K]\n",
      "\n",
      "# For propane\n",
      "Tc = 369.8;\t\t\t#[K]\n",
      "Pc = 42.48;\t\t\t#[bar]\n",
      "Pc = Pc*10**(5);\n",
      "w = 0.152;\n",
      "\n",
      "# Calculations and Results\n",
      "#(1)\n",
      "# For van der Walls equation of state\n",
      "a = (27*R**(2)*Tc**(2))/(64*Pc);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = (R*Tc)/(8*Pc);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# The cubic form of van der Walls equation of state is given by,\n",
      "# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0\n",
      "\n",
      "# At state 1 (100 C and 6 bar) \n",
      "def f(V): \n",
      "\t return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1\n",
      "V_1_1 = fsolve(f,1)\n",
      "V_1_2 = fsolve(f,10)\n",
      "V_1_3 = fsolve(f,100)\n",
      "# The largest root is considered because of molar volume of vapour phase is to determined\n",
      "V_1 = V_1_1;\t\t\t#[m**(3)/mol]\n",
      "# Thus compressibility factor is\n",
      "Z_1 = (P_1*V_1)/(R*T_1);\t\t\t#compressibility factor\n",
      "\n",
      "H_R_1 = (Z_1 - 1)*R*T_1 - (a/V_1);\t\t\t# [J/mol]\n",
      "S_R_1 = R*math.log((P_1*(V_1-b))/(R*T_1));\t\t\t# [J/mol-K]\n",
      "\n",
      "# At state 2 (500 C and 12 bar) \n",
      "def f1(V): \n",
      "    return V**(3)-(b+(R*T_2)/P_2)*V**(2)+(a/P_2)*V-(a*b)/P_2\n",
      "V_2_1 = fsolve(f1,1)\n",
      "V_2_2 = fsolve(f1,10)\n",
      "V_2_3 = fsolve(f1,100)\n",
      "# The largest root is considered because of molar volume of vapour phase is to determined\n",
      "V_2 = V_2_1;\t\t\t#[m**(3)/mol]\n",
      "# Thus compressibility factor is\n",
      "Z_2 = (P_2*V_2)/(R*T_2);\t\t\t#compressibility factor\n",
      "\n",
      "H_R_2 = (Z_2 - 1)*R*T_2 - (a/V_2);\t\t\t# [J/mol]\n",
      "S_R_2 = R*math.log((P_2*(V_2-b))/(R*T_2));\t\t\t# [J/mol-K]\n",
      "\n",
      "# Ideal gas entropy change is given by\n",
      "delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1);\t\t\t#[J/mol-K]\n",
      "# Entropy change is given by\n",
      "delta_S = delta_S_ig + (S_R_2 - S_R_1);\t\t\t#[J/mol-k]\n",
      "\n",
      "# Ideal gas enthalpy change is given by\n",
      "delta_H_ig = Cp_0*(T_2 - T_1);\t\t\t#[J/mol]\n",
      "# Enthalpy change is given by\n",
      "delta_H = delta_H_ig + (H_R_2 - H_R_1);\t\t\t#[J/mol]\n",
      "\n",
      "print \"1).The change in enthalpy is %f J/mol\"%(delta_H);\n",
      "print \"    The change in entropy is %f J/mol-K\"%(delta_S);\n",
      "\n",
      "#(2)\n",
      "# Virial equation of state\n",
      "\n",
      "# At state 1 (372.15 K, 6 bar) let us calculate B and dB/dT\n",
      "Tr = T_1/Tc;\t\t\t# Reduced temperature\n",
      "B_0 = 0.083-(0.422/(Tr)**(1.6));\n",
      "B_1 = 0.139-(0.172/(Tr)**(4.2));\n",
      "\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B = ((B_0+(w*B_1))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol]\n",
      "dB0_dT = 0.422*1.6*Tc**(1.6)*T_1**(-2.6);\t\t\t# (dB_0/dT) at state 1\n",
      "dB1_dT = 0.172*4.2*Tc**(4.2)*T_1**(-5.2);\t\t\t# (dB_1/dT) at state 1\n",
      "dB_dT = ((R*Tc)/Pc)*((dB0_dT) + w*(dB1_dT));\t\t\t# (dB/dT) at state 1\n",
      "\n",
      "H_R_1_2 = B*P_1 - P_1*T_1*dB_dT;\t\t\t#[J/mol] - Residual enthalpy at state 1\n",
      "S_R_1_2 = -P_1*(dB_dT);\t\t\t#[J/mol-K] - Residual entropy at state 1\n",
      "\n",
      "# At state 2 (773.15 K, 12 bar)\n",
      "Tr_2 = T_2/Tc;\t\t\t# Reduced temperature\n",
      "B_0_2 = 0.083-(0.422/(Tr_2)**(1.6));\n",
      "B_1_2 = 0.139-(0.172/(Tr_2)**(4.2));\n",
      "\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B_2 = ((B_0_2+(w*B_1_2))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol]\n",
      "dB0_dT_2 = 0.422*1.6*Tc**(1.6)*T_2**(-2.6);\t\t\t# (dB_0/dT) at state 1\n",
      "dB1_dT_2 = 0.172*4.2*Tc**(4.2)*T_2**(-5.2);\t\t\t# (dB_1/dT) at state 1\n",
      "dB_dT_2 = ((R*Tc)/Pc)*((dB0_dT_2) + w*(dB1_dT_2));\t\t\t# (dB/dT) at state 1\n",
      "\n",
      "H_R_2_2 = B_2*P_2 - P_2*T_2*dB_dT_2;\t\t\t#[J/mol] - Residual enthalpy at state 1\n",
      "S_R_2_2 = -P_2*(dB_dT_2);\t\t\t#[J/mol-K] - Residual entropy at state 1\n",
      "\n",
      "delta_H_2 = delta_H_ig + (H_R_2_2 - H_R_1_2);\t\t\t#[J/mol]\n",
      "delta_S_2 = delta_S_ig + (S_R_2_2 - S_R_1_2);\t\t\t#[J/mol]\n",
      "\n",
      "print \"2).The change in enthalpy is %f J/mol\"%(delta_H_2);\n",
      "print \"    The change in entropy is %f J/mol-K\"%(delta_S_2);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1).The change in enthalpy is 29798.211799 J/mol\n",
        "    The change in entropy is 48.649067 J/mol-K\n",
        "2).The change in enthalpy is 29992.807365 J/mol\n",
        "    The change in entropy is 49.021724 J/mol-K\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.21  Page Number : 362"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of fugacity\n",
      "\n",
      "# Variables\n",
      "P = 2.76*10**(6);\t\t\t#[N/m**(2)] - Pressure\n",
      "T = 310.93;\t\t\t#[K] - Temperature\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "# For n-butane\n",
      "Tc = 425.18;\t\t\t#[K] - Critical temperature\n",
      "Pc = 37.97;\t\t\t#[bar] - Critical pressure\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.193;\n",
      "den = 0.61;\t\t\t#[g/cm**(3)]\n",
      "mol_wt = 58;\t\t\t#[g/mol] - Molecular weight of butane\n",
      "\n",
      "# Calculations and Results\n",
      "# math.log(P_sat) = 15.7374 - 2151.63/(T-36.24)\n",
      "P_sat = math.exp(15.7374 - 2151.63/(T-36.24));\t\t\t#[mm Hg]\n",
      "P_sat = (P_sat/760)*101325;\t\t\t#[N/m**(2)]\n",
      "\n",
      "#(1)\n",
      "# Let us determine the second virial coefficient at 310.93 K\n",
      "Tr = T/Tc;\t\t\t# Reduced temperature\n",
      "B_0 = 0.083-(0.422/(Tr)**(1.6));\n",
      "B_1 = 0.139-(0.172/(Tr)**(4.2));\n",
      "#We know,(B*Pc)/(R*Tc) = B_0 + (w*B_1)\n",
      "B = ((B_0+(w*B_1))*(R*Tc))/Pc;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# Fugacity under saturated conditions is given by\n",
      "# math.log(f_sat/P_sat) = (B*P_sat)/(R*T)\n",
      "f_sat = P_sat*(math.exp((B*P_sat)/(R*T)));\t\t\t#[N/m**(2)]\n",
      "\n",
      "# The molar volume is given by\n",
      "V_liq = (1/(den*1000))*(mol_wt/1000);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "f = f_sat*math.exp(V_liq*(P-P_sat)/(R*T));\n",
      "\n",
      "print \" 1).The fugacity of n-butane is %e N/m**2)\"%(f);\n",
      "\n",
      "#(2)\n",
      "# For van der Walls equation of state\n",
      "a = (27*R**(2)*Tc**(2))/(64*Pc);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = (R*Tc)/(8*Pc);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# The cubic form of van der Walls equation of state is given by,\n",
      "# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0\n",
      "\n",
      "# At 100 C and 1 atm \n",
      "def f(V): \n",
      "    return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P\n",
      "V_1 = fsolve(f,0.1)\n",
      "V_1 = fsolve(f,10)\n",
      "V_1 = fsolve(f,100)\n",
      "# The above equation has only 1 real root, other two roots are imaginary\n",
      "V = V_1;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# math.log(f/P) = math.log((R*T)/(P*(V-b))) + b/(V-b) -(2*a)/(R*T*V)\n",
      "f_2 = P*(math.exp(math.log((R*T)/(P*(V-b))) + b/(V-b) -(2*a)/(R*T*V)));\n",
      "\n",
      "print \" 2).The fugacity of n-butane is %e N/m**2)\"%(f_2);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The fugacity of n-butane is 3.293476e+05 N/m**2)\n",
        " 2).The fugacity of n-butane is 8.997352e+05 N/m**2)\n"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.22  Page Number : 363"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of enthalpy change\n",
      "\n",
      "# Variables\n",
      "T = 50+273.15;\t\t\t#[K] - Temperature\n",
      "P = 25.*10**(3);\t\t\t#[Pa] - Pressure\n",
      "y1 = 0.5;\t\t\t#[mol] - mole fraction of equimolar mixture\n",
      "y2 = 0.5;\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "#For component 1 (methyl ethyl ketone)\n",
      "Tc_1  = 535.5;\t\t\t#[K] - Critical temperature\n",
      "Pc_1 = 41.5*10**(5);\t\t\t#[N/m**(2)] - Critical pressure\n",
      "Vc_1 = 267.;\t\t\t#[cm**(3)/mol] - Critical volume\n",
      "Zc_1 = 0.249;\t\t\t# Critical compressibility factor\n",
      "w_1 = 0.323;\t\t\t# acentric factor\n",
      "\n",
      "#For component 2 (toluene)\n",
      "Tc_2 = 591.8;\t\t\t#[K]\n",
      "Pc_2 = 41.06*10**(5);\t\t\t#[N/m**(2)]\n",
      "Vc_2 = 316.;\t\t\t#[cm**(3)/mol]\n",
      "Zc_2 = 0.264;\n",
      "w_2 = 0.262;\n",
      "\n",
      "# Calculations and Results\n",
      "# For equation of state Z = 1 + B/V\n",
      "#For component 1, let us calculate B and dB/dT\n",
      "Tr_1 = T/Tc_1;\t\t\t#Reduced temperature\n",
      "#At reduced temperature\n",
      "B1_0 = 0.083-(0.422/(Tr_1)**(1.6));\n",
      "B1_1 = 0.139-(0.172/(Tr_1)**(4.2));\n",
      "#We know,(B*Pc)/(R*Tc) = B_0+(w*B_1)\n",
      "B_11 = ((B1_0+(w_1*B1_1))*(R*Tc_1))/Pc_1;\t\t\t# [m**(3)/mol-K] \n",
      "dB0_dT_1 = 0.422*1.6*Tc_1**(1.6)*T**(-2.6);\t\t\t# [m**(3)/mol-K] -  (dB_0/dT)\n",
      "dB1_dT_1 = 0.172*4.2*Tc_1**(4.2)*T**(-5.2);\t\t\t# [m**(3)/mol-K] -  (dB_1/dT)\n",
      "dB_dT_1 = ((R*Tc_1)/Pc_1)*((dB0_dT_1) + w_1*(dB1_dT_1));\t\t\t#[m**(3)/mol-K] -  (dB/dT)_\n",
      "\n",
      "#Similarly for component 2\n",
      "Tr_2 = T/Tc_2;\t\t\t#Reduced temperature\n",
      "#At reduced temperature Tr_2,\n",
      "B2_0 = 0.083 - (0.422/(Tr_2)**(1.6));\n",
      "B2_1 = 0.139 - (0.172/(Tr_2)**(4.2));\n",
      "B_22 = ((B2_0+(w_2*B2_1))*(R*Tc_2))/Pc_2;\t\t\t#[m**(3)/mol]\n",
      "dB0_dT_2 = 0.422*1.6*Tc_2**(1.6)*T**(-2.6);\t\t\t# [m**(3)/mol-K] -  (dB_0/dT)\n",
      "dB1_dT_2 = 0.172*4.2*Tc_2**(4.2)*T**(-5.2);\t\t\t# [m**(3)/mol-K] -  (dB_1/dT)\n",
      "dB_dT_2 = ((R*Tc_2)/Pc_2)*((dB0_dT_2) + w_2*(dB1_dT_2));\t\t\t#[m**(3)/mol-K] -  (dB/dT)_\n",
      "\n",
      "#For cross coeffcient, let us calculate B and dB/dT\n",
      "Tc_12 = (Tc_1*Tc_2)**(1./2);\t\t\t#[K]\n",
      "w_12 = (w_1 + w_2)/2.;\n",
      "Zc_12 = (Zc_1 + Zc_2)/2;\n",
      "Vc_12 = (((Vc_1)**(1./3) + (Vc_2)**(1./3))/2)**(3);\t\t\t#[cm**(3)/mol]\n",
      "Vc_12 = Vc_12*10**(-6);\t\t\t#[m**(3)/mol]\n",
      "Pc_12 = (Zc_12*R*Tc_12)/Vc_12;\t\t\t#[N/m**(2)]\n",
      "\n",
      "#Now we have,(B_12*Pc_12)/(R*Tc_12) = B_0+(w_12*B_1)\n",
      "#where B_0 and B_1 are to be evaluated at Tr_12\n",
      "Tr_12 = T/Tc_12;\n",
      "#At reduced temperature Tr_12\n",
      "B_0 = 0.083 - (0.422/(Tr_12)**(1.6));\n",
      "B_1 = 0.139 - (0.172/(Tr_12)**(4.2));\n",
      "B_12 = ((B_0 + (w_12*B_1))*R*Tc_12)/Pc_12;\t\t\t#[m**(3)/mol]\n",
      "dB0_dT_12 = 0.422*1.6*Tc_12**(1.6)*T**(-2.6);\t\t\t# [m**(3)/mol-K] -  (dB_0/dT)\n",
      "dB1_dT_12 = 0.172*4.2*Tc_12**(4.2)*T**(-5.2);\t\t\t# [m**(3)/mol-K] -  (dB_1/dT)\n",
      "dB_dT_12 = ((R*Tc_12)/Pc_12)*((dB0_dT_12) + w_12*(dB1_dT_12));\t\t\t#[m**(3)/mol-K] -  (dB/dT)_12\n",
      "\n",
      "#For the mixture\n",
      "B = y1**(2)*B_11 + 2*y1*y2*B_12 + y2**(2)*B_22;\t\t\t#[m**(3)/moL]\n",
      "\n",
      "# The equation of state can be written as\n",
      "# V**(2) - ((R*T)/P) - (B*R*T)/P = 0\n",
      "# V**(2) - 0.1075*V + 1.737*10**(-4) = 0\n",
      "def f(V): \n",
      "    return V**(2) - 0.1075*V + 1.737*10**(-4)\n",
      "V1 = fsolve(f,0.1)\n",
      "V2 = fsolve(f,1)\n",
      "# We will consider the root which is near to R*T/P\n",
      "V = V1;\n",
      "# dB/dT = y_1**(2)*dB_11/dT + y_2**(2)*dB_22/dT + 2*y_1*y_2*dB_12/dT\n",
      "dB_dT = y1**(2)*dB_dT_1 + y2**(2)*dB_dT_2 + 2*y1*y2*dB_dT_12;\t\t\t#[m**(3)/mol-K]\n",
      "\n",
      "# For equation of state Z = 1 + B/V\n",
      "H_R = (B*R*T)/V - ((R*T**(2))/V)*dB_dT;\t\t\t#[J/mol]\n",
      "\n",
      "print \" 1).The value of H_R for the mixture using virial equation of state is %f J/mol\"%(H_R);\n",
      "\n",
      "#(2)\n",
      "#  For van der Walls equation of state \n",
      "a_11 = (27*R**(2)*Tc_1**(2))/(64*Pc_1);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "a_22 = (27*R**(2)*Tc_2**(2))/(64*Pc_2);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "a_12 = (a_11*a_22)**(1./2);\n",
      "b_1 = (R*Tc_1)/(8*Pc_1);\t\t\t#[m**(3)/mol]\n",
      "b_2 = (R*Tc_2)/(8*Pc_2);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# For the mixture\n",
      "a = y1**(2)*a_11 + y2**(2)*a_22 + 2*y1*y2*a_12;\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = y1*b_1 + y2*b_2;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# From the cubic form of van der Walls equation of state\n",
      "def f1(V): \n",
      "    return V**(3)-(b+(R*T)/P)*V**(2)+(a/P)*V-(a*b)/P\n",
      "V2_1 = fsolve(f1,0.1)\n",
      "V2_2 = fsolve(f1,10)\n",
      "V2_3 = fsolve(f1,100)\n",
      "# The largest root is considered\n",
      "V_2 = V2_1;\n",
      "\n",
      "# The residual enthalpy is given by\n",
      "H_R_2 = P*V_2 - R*T -a/V_2;\t\t\t#[J/mol]\n",
      "\n",
      "print \" 2).The value of H_R for the mixture using van der Walls equation of state is %f J/mol\"%(H_R_2);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The value of H_R for the mixture using virial equation of state is -151.289795 J/mol\n",
        " 2).The value of H_R for the mixture using van der Walls equation of state is -38.476127 J/mol\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.23  Page Number : 366"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of fugacity of water vapour\n",
      "\n",
      "# Variables\n",
      "T = 320 + 273.15;\t\t\t#[K]\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "\n",
      "# For water\n",
      "Tc = 647.1;\t\t\t#[K]\n",
      "Pc = 220.55;\t\t\t#[bar]\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "\n",
      "# The cubic form of Redlich Kwong equation of state is given by,\n",
      "# V**(3) - ((R*T)/P)*V**(2) - ((b_1**(2)) + ((b_1*R*T)/P) - (a/(T**(1/2)*P))*V - (a*b)/(T**(1/2)*P) = 0\n",
      "\n",
      "# At 320 C and 70 bar pressure\n",
      "P_1 = 70;\t\t\t#[bar]\n",
      "P_1 = P_1*10**(5);\t\t\t#[Pa]\n",
      "\n",
      "# Calculations and Results\n",
      "a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc;\t\t\t#[Pa*m**(6)*K**(1/2)/mol]\n",
      "b = (0.08664*R*Tc)/Pc;\t\t\t#[m**(3)/mol]\n",
      "# Solving the cubic equation\n",
      "def f1(V): \n",
      "    return V**(3)-((R*T)/P_1)*V**(2)-((b**(2))+((b*R*T)/P_1)-(a/(T**(1./2)*P_1)))*V-(a*b)/(T**(1./2)*P_1)\n",
      "V1=fsolve(f1,1)\n",
      "V2=fsolve(f1,10)\n",
      "V3=fsolve(f1,100)\n",
      "# The largest root is considered because at 320 C and 70 bar vapour phase exists.\n",
      "V_1 = V1;\t\t\t#[m**(3)/mol]\n",
      "# Thus compressibility factor is\n",
      "Z_1 = (P_1*V_1)/(R*T);\n",
      "\n",
      "# For Redlich-Kwong equation of state\n",
      "# math.log(f/P) = Z - 1 - math.log(V_1/(V_1-b)) + (a/(b*R*(T**(3/2))))*math.log(V/(V+b))\n",
      "f_1 = P_1*(math.exp(Z_1-1-math.log(Z_1)+math.log(V_1/(V_1-b))+(a/(b*R*(T**(3./2))))*math.log(V_1/(V_1+b))));\t\t\t#[Pa]\n",
      "f_1 = f_1*10**(-5);\t\t\t#[bar]\n",
      "\n",
      "print \" The fugacity of water vapour at 320 C and 70 bar pressure is %f bar\"%(f_1);\n",
      "\n",
      "# At 320 C and 170 bar pressure, we have\n",
      "P_2 = 170;\t\t\t#[bar]\n",
      "P_2 = P_2*10**(5);\t\t\t#[Pa]\n",
      "\n",
      "# Solving the cubic equation\n",
      "def f2(V): \n",
      "    return V**(3)-((R*T)/P_2)*V**(2)-((b**(2))+((b*R*T)/P_2)-(a/(T**(1./2)*P_2)))*V-(a*b)/(T**(1./2)*P_2)\n",
      "V4 = fsolve(f2,1)\n",
      "V5 = fsolve(f2,10)\n",
      "V6 = fsolve(f2,100)\n",
      "# The above equation has only 1 real root,other two roots are imaginary. Therefore,\n",
      "V_2 = V6;\t\t\t#[m**(3)/mol]\n",
      "# Thus compressibility factor is\n",
      "Z_2 = (P_2*V_2)/(R*T);\n",
      "\n",
      "# For Redlich-Kwong equation of state\n",
      "# math.log(f/P) = Z - 1 - math.log(V_1/(V_1-b)) + (a/(b*R*(T**(3/2))))*math.log(V/(V+b))\n",
      "f_2 = P_2*(math.exp(Z_2-1-math.log(Z_2)+math.log(V_2/(V_2-b))+(a/(b*R*(T**(3./2))))*math.log(V_2/(V_2+b))));\t\t\t#[Pa]\n",
      "f_2 = f_2*10**(-5);\t\t\t#[bar]\n",
      "\n",
      "print \" The fugacity of water vapour at 320 C and 170 bar pressure is %f bar\"%(f_2);\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The fugacity of water vapour at 320 C and 70 bar pressure is 60.456239 bar\n",
        " The fugacity of water vapour at 320 C and 170 bar pressure is 101.590989 bar\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/usr/lib/python2.7/dist-packages/scipy/optimize/minpack.py:227: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
        "  improvement from the last ten iterations.\n",
        "  warnings.warn(msg, RuntimeWarning)\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.24  Page Number : 367"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Determination of change in internal energy\n",
      "\n",
      "# Variables\n",
      "Vol = 0.057;\t\t\t#[m**(3)] - Volume of car tyre\n",
      "P_1 = 300.;\t\t\t#[kPa] - Initial pressure\n",
      "P_1 = P_1*10**(3);\t\t\t#[Pa]\n",
      "T_1 = 300.;\t\t\t#[K] - Initial temperature\n",
      "P_2 = 330.;\t\t\t#[kPa] - Finnal pressure\n",
      "P_2 = P_2*10**(3);\t\t\t#[Pa]\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "Cv_0 = 21.;\t\t\t#[J/mol-K] - Heat capacity for air \n",
      "\n",
      "# For oxygen\n",
      "Tc_O2 = 154.6;\t\t\t#[K] - Critical temperature\n",
      "Pc_O2 = 50.43;\t\t\t#[bar] - Critical pressure\n",
      "Pc_O2 = Pc_O2*10**(5);\t\t\t#[Pa]\n",
      "y1 = 0.21;\t\t\t# - Mole fraction of oxygen\n",
      "# For nitrogen\n",
      "Tc_N2 = 126.2;\t\t\t#[K] - Critical temperature\n",
      "Pc_N2 = 34.00;\t\t\t#[bar] - Critical pressure\n",
      "Pc_N2 = Pc_N2*10**(5);\t\t\t#[Pa]\n",
      "y2 = 0.79;\t\t\t# - Mole fraction of nitrogen\n",
      "\n",
      "# Calculations and Results\n",
      "# (1)\n",
      "# Assuming ideal gas behaviour. The volume remains the same,therefore,we get\n",
      "# P_1/T_1 = P_2/T_2\n",
      "T_2 = P_2*(T_1/P_1);\t\t\t#[K]\n",
      "\n",
      "n = (P_1*Vol)/(R*T_1);\t\t\t#[mol] - Number of moles\n",
      "delta_U = n*Cv_0*(T_2-T_1);\t\t\t#[J]\n",
      "\n",
      "print \" 1).The change in internal energy for ideal gas behaviour) is %f J\"%(delta_U);\n",
      "\n",
      "#(2)\n",
      "#  For van der Walls equation of state \n",
      "a_O2 = (27*R**(2)*Tc_O2**(2))/(64*Pc_O2);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "a_N2 = (27*R**(2)*Tc_N2**(2))/(64*Pc_N2);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "a_12 = (a_O2*a_N2)**(1./2);\n",
      "b_O2 = (R*Tc_O2)/(8*Pc_O2);\t\t\t#[m**(3)/mol]\n",
      "b_N2 = (R*Tc_N2)/(8*Pc_N2);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# For the mixture\n",
      "a = y1**(2)*a_O2 + y2**(2)*a_N2 + 2*y1*y2*a_12;\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = y1*b_O2 + y2*b_N2;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# From the cubic form of van der Walls equation of state\n",
      "# At 300 K and 300 kPa,\n",
      "def f1(V): \n",
      "    return V**(3)-(b+(R*T_1)/P_1)*V**(2)+(a/P_1)*V-(a*b)/P_1\n",
      "V_1 = fsolve(f1,0.1)\n",
      "V_2 = fsolve(f1,10)\n",
      "V_3 = fsolve(f1,100)\n",
      "# The above equation has only 1 real root, other two roots are imaginary\n",
      "V = V_1;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# Now at P = 330 kPa and at molar volume V\n",
      "# The van der Walls equation of state is\n",
      "# (P + a/V**(2))*(V - b) = R*T\n",
      "T_2_prime = ((P_2 + a/V**(2))*(V - b))/R;\t\t\t#[K]\n",
      "n_prime = Vol/V;\t\t\t#[mol]\n",
      "\n",
      "# Total change in internal energy is given by\n",
      "# delta_U_prime = n_prime*delta_U_ig + n_prime*(U_R_2 - U_R_1)\n",
      "# delta_U_prime = n_prime*Cv_0*(T_2_prime - T_1) + n_prime*a(1/V_2 - 1/V_1);\n",
      "# Since V_1 = V_2 = V, therefore\n",
      "delta_U_prime = n_prime*Cv_0*(T_2_prime - T_1);\n",
      "\n",
      "print \" 2).The change in internal energy for van der Walls equation of state) is %f J\"%(delta_U_prime);\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The change in internal energy for ideal gas behaviour) is 4319.220592 J\n",
        " 2).The change in internal energy for van der Walls equation of state) is 4299.872282 J\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.25  Page Number : 369"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of enthalpy and entropy change\n",
      "\n",
      "# Variables\n",
      "T_1 = 150 + 273.15;\t\t\t#[K] - Initial emperature\n",
      "T_2 = T_1;\t\t\t# Isothermal process\n",
      "P_1 = 100.*10**(3);\t\t\t#[Pa] - Initial pressure\n",
      "P_2 = 450.*10**(3);\t\t\t#[Pa] - Final pressure\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "# For water\n",
      "Tc = 647.1;\t\t\t#[K] - Critical temperature\n",
      "Pc = 220.55;\t\t\t#[bar] - Critical pressure\n",
      "Pc = Pc*10.**(5);\n",
      "w = 0.345;\n",
      "Mol_wt = 18.015;\t\t\t#[g/mol] - Molecular weight of water\n",
      "Cp_0 = 4.18;\t\t\t#[J/mol-K] - Smath.tan(math.radiansard heat capacity of water\n",
      "\n",
      "# Calculations\n",
      "# Both phases are superheated vapour phases because at 150 C the vapour pressure of steam is 4.67 bar and both operating pressures are below saturated pressure.\n",
      "# In Peng-Robinson equation of state\n",
      "m = 0.37464 + 1.54226*w - 0.26992*w**(2);\n",
      "# At T_1 and P_1, we have\n",
      "Tr = T_1/Tc;\n",
      "alpha = (1 + m*(1 - Tr**(1./2)))**(2);\n",
      "a = ((0.45724*(R*Tc)**(2))/Pc)*alpha;\t\t\t#[Pa*m**(6)/mol**(2)]\n",
      "b = (0.07780*R*Tc)/Pc;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# Cubuc form of Peng-Robinson equation of stste is given by\n",
      "# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;\n",
      "# Solving the cubic equation\n",
      "def f(V): \n",
      "    return V**(3)+(b-(R*T_1)/P_1)*V**(2)-((3*b**(2))+((2*R*T_1*b)/P_1)-(a/P_1))*V+b**(3)+((R*T_1*(b**(2)))/P_1)-((a*b)/P_1)\n",
      "V1 = fsolve(f,-1)\n",
      "V2 = fsolve(f,0)\n",
      "V3 = fsolve(f,1)\n",
      "#The largest root is for vapour phase,\n",
      "#The largest root is only considered as the systemis gas\n",
      "V_1 = V3;\t\t\t#[m**(3)/mol]\n",
      "# Thus compressibility factor is\n",
      "Z_1 = (P_1*V_1)/(R*T_1);\n",
      "\n",
      "# At T_2 and P_2, we have\n",
      "# Cubuc form of Peng-Robinson equation of stste is given by\n",
      "# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;\n",
      "# Solving the cubic equation\n",
      "def f1(V): \n",
      "    return V**(3)+(b-(R*T_2)/P_2)*V**(2)-((3*b**(2))+((2*R*T_2*b)/P_2)-(a/P_2))*V+b**(3)+((R*T_2*(b**(2)))/P_2)-((a*b)/P_2)\n",
      "V4 = fsolve(f1,-1)\n",
      "V5 = fsolve(f1,0)\n",
      "V6 = fsolve(f1,1)\n",
      "#The largest root is for vapour phase,\n",
      "#The largest root is only considered as the systemis gas\n",
      "V_2 = V6;\t\t\t#[m**(3)/mol]\n",
      "# Thus compressibility factor is\n",
      "Z_2 = (P_2*V_2)/(R*T_2);\n",
      "\n",
      "# In the Peng-Robinson equation of stste\n",
      "# da/dT = -(a*m)/((alpha*T*Tc)**(1/2))\n",
      "# The residual enthalpy is given by\n",
      "# H_R = R*T*(Z-1) + (((T*(da_dT))-a)/(2*2**(1/2)*b))*math.log((Z+(1+2**(1/2)*((P*b)/(R*T))))/(Z+(1-2**(1/2)*((P*b)/(R*T)))))\n",
      "\n",
      "# At state 1\n",
      "da_dT_1 = -(a*m)/((alpha*T_1*Tc)**(1./2));\t\t\t#[Pa*m**(6)/mol**(2)]\n",
      "H_R_1 = R*T_1*(Z_1-1) + (((T_1*(da_dT_1))-a)/(2*(2**(1./2))*b))* \\\n",
      "math.log((Z_1+(1+2**(1./2))*((P_1*b)/(R*T_1)))/(Z_1+(1-2**(1./2))*((P_1*b)/(R*T_1))));\n",
      "\n",
      "# At state 2\n",
      "da_dT_2 = -(a*m)/((alpha*T_2*Tc)**(1/2.));\t\t\t#[Pa*m**(6)/mol**(2)]\n",
      "H_R_2 = R*T_2*(Z_2-1) + (((T_2*(da_dT_2))-a)/(2*2**(1./2)* \\\n",
      "b))*math.log((Z_2+(1+2**(1./2))*((P_2*b)/(R*T_1)))/(Z_2+(1-2**(1./2))*((P_2*b)/(R*T_1))));\n",
      "\n",
      "# Since the temperature is the same,therefore ideal gas change in enthalpy is zero and thus\n",
      "delta_H = H_R_2 - H_R_1;\t\t\t#[J/mol]\n",
      "delta_H = delta_H/Mol_wt;\t\t\t#[kJ/kg]\n",
      "\n",
      "# The residual entropy relation for a substance following Peng - Robinson equation of state ia\n",
      "# S_R = R*math.log(Z - (P*b)/(R*T)) + (da_dT/(2*2**(1/2)*b))*math.log((Z+(1+2**(1/2))*((P*b)/(R*T)))/(Z+(1-2**(1/2))*((P*b)/(R*T))))\n",
      "\n",
      "# The residual entropy at state 1 is\n",
      "S_R_1 = R*math.log(Z_1 - (P_1*b)/(R*T_1)) + (da_dT_1/(2*2**(1./2)*b))* \\\n",
      "math.log((Z_1+(1+2**(1./2))*((P_1*b)/(R*T_1)))/(Z_1+(1-2**(1./2))*((P_1*b)/(R*T_1))));\n",
      "\n",
      "# The residual entropy at state 2 is\n",
      "S_R_2 = R*math.log(Z_2 - (P_2*b)/(R*T_2)) + (da_dT_2/(2*2**(1./2)*b)) \\\n",
      "*math.log((Z_2+(1+2**(1./2))*((P_2*b)/(R*T_2)))/(Z_2+(1-2**(1./2))*((P_2*b)/(R*T_2))));\n",
      "\n",
      "delta_S_R = S_R_2 - S_R_1;\t\t\t#[J/mol-K]\n",
      "\n",
      "# The ideal gas change in entropy is\n",
      "delta_S_ig = Cp_0*math.log(T_2/T_1) - R*math.log(P_2/P_1);\t\t\t#[J/mol-K]\n",
      "\n",
      "# Therefore\n",
      "delta_S = delta_S_R + delta_S_ig;\t\t\t#[J/mol-K]\n",
      "\n",
      "# Results\n",
      "print \" The enthalpy change is given by delta_H = %f kJ/mol\"%(delta_H);\n",
      "print \" The entropy change is given by  delta_S = %f J/mol-K\"%(delta_S);\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The enthalpy change is given by delta_H = -11.747722 kJ/mol\n",
        " The entropy change is given by  delta_S = -12.826050 J/mol-K\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.26  Page Number : 370"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of final temperature and pressure\n",
      "\n",
      "# Variables\n",
      "Vol = 0.15;\t\t\t#[m**(3)]\n",
      "T_1 = 170;\t\t\t#[K] - Initial emperature\n",
      "P_1 = 100;\t\t\t#[bar] - Initial pressure\n",
      "P_1 = P_1*10**(5);\t\t\t#[Pa]\n",
      "R = 8.314;\t\t\t#[J/mol*K] - Universal gas constant\n",
      "# For nitrogen\n",
      "Tc = 126.2;\t\t\t#[K] - Critical tempeature\n",
      "Pc = 34;\t\t\t#[bar] - Critical pressure\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.038;\n",
      "# Cp_0 = 27.2+4.2*10**(-3)*T\n",
      "\n",
      "# Calculations and Results\n",
      "#(1)\n",
      "# For van der Walls equation of state\n",
      "a = (27*R**(2)*Tc**(2))/(64*Pc);\t\t\t#[Pa-m**(6)/mol**(2)]\n",
      "b = (R*Tc)/(8*Pc);\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# The cubic form of van der Walls equation of state is given by,\n",
      "# V**(3) - (b + (R*T)/P)*V**(2) + (a/P)*V - (a*b)/P = 0\n",
      "# On simplification the equation changes to\n",
      "# V**(3) - 1.799*10**(4)*V**(2) + 1.366*10**(-8)*V - 5.269*10**(-13) = 0\n",
      "\n",
      "# Solving the cubic equation \n",
      "def f(V): \n",
      "    return V**(3)-1.799*10**(-4)*V**(2) + 1.366*10**(-8)*V - 5.269*10**(-13)\n",
      "V1 = fsolve(f,1)\n",
      "V2 = fsolve(f,10)\n",
      "V3 = fsolve(f,100)\n",
      "# The above equation has only 1 real root, other two roots are imagimnary\n",
      "V_1 = V1;\t\t\t#[m**(3)/mol]\n",
      "# Thus total number of moles is given by\n",
      "n_1 = Vol/V_1;\t\t\t#[mol]\n",
      "\n",
      "# After 500 mol are withdrawn, the final number of moles is given by\n",
      "n_2 = n_1 - 500;\t\t\t#[mol]\n",
      "# Thus molar volume at final state is \n",
      "V_2 = Vol/n_2;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# The ideal entropy change is guven by\n",
      "\n",
      "def f24(T): \n",
      "    return 27.2+4.2*10**(-3)*T\n",
      "\n",
      "# delta_S_ig =  quad(f24,T_1,T_2) - R*math.log(P_2/P_1)[0]\n",
      "\n",
      "# The residual entropy change is given by\n",
      "# delta_S_R = R*math.log((P_2*(V_2-b))/(R*T_2)) - R*math.log((P_1*(V_1-b))/(R*T_1)) \n",
      "# delta_S = delta_S_ig = delta_S_R\n",
      "\n",
      "def f25(T): \n",
      "    return 27.2+4.2*10**(-3)*T\n",
      "\n",
      "# delta_S =  quad(f25,T_1,T_2) + R*math.log((V_2-b)/(V_1-b))[0]\n",
      "\n",
      "# During discharging delta_S = 0, thus on simplification we get\n",
      "# 18.886*math.log(T_2) + 4.2*10**(-3)*T_2 - 92.937 = 0\n",
      "# Solving the above equation we get\n",
      "def f1(T_2): \n",
      "    return 18.886*math.log(T_2) + 4.2*10**(-3)*T_2 - 92.937\n",
      "T_2 = fsolve(f1,1)\n",
      "\n",
      "# Thus at T_2, \n",
      "P_2 = (R*T_2)/(V_2-b) - a/V_2**(2);\t\t\t#[N/m**(2)]\n",
      "P_2 = P_2*10**(-5);\t\t\t#[bar]\n",
      "\n",
      "print \" 1).The final temperature is %f K\"%(T_2);\n",
      "print \"     The final pressure is %f bar\"%(P_2);\n",
      "\n",
      "#(2)\n",
      "# In Peng-Robinson equation of state\n",
      "m = 0.37464 + 1.54226*w - 0.26992*w**(2);\n",
      "# At T_1 and P_1, we have\n",
      "Tr = T_1/Tc;\n",
      "alpha = (1 + m*(1 - Tr**(1./2)))**(2);\n",
      "a_2 = ((0.45724*(R*Tc)**(2))/Pc)*alpha;\t\t\t#[Pa*m**(6)/mol**(2)]\n",
      "b_2 = (0.07780*R*Tc)/Pc;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# Cubuc form of Peng-Robinson equation of stste is given by\n",
      "# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;\n",
      "# Solving the cubic equation\n",
      "def f2(V): \n",
      "    return V**(3)+(b_2-(R*T_1)/P_1)*V**(2)-((3*b_2**(2))+((2*R*T_1*b_2)/P_1)-(a_2/P_1))*V+b_2**(3)+((R*T_1*(b_2**(2)))/P_1)-((a_2*b_2)/P_1)\n",
      "V4 = fsolve(f2,-1)\n",
      "V5 = fsolve(f2,0)\n",
      "V6 = fsolve(f2,0.006)\n",
      "#The above equation has only 1 real root,the other two roots are imaginary\n",
      "V_1_2 = V6;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# The number of moles in the initial state is given by\n",
      "n_1_2 = Vol/V_1_2;\t\t\t#[mol]\n",
      "# After 500 mol are withdrawn, the final number of moles is given by\n",
      "n_2_2 = n_1_2 - 500;\t\t\t#[mol]\n",
      "# Thus molar volume at final state is \n",
      "V_2_2 = Vol/n_2_2;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# At the final state the relation between pressure and temperature is\n",
      "# P_2_2 = (R*T_2_2)/(V_2_2-b_2) - a_2/V_2_2**(2)\n",
      "# P_2_2 = 7.23*10**(4)*T_2 - 3.93*10**(7)*a_2\n",
      "\n",
      "# Now let us calculate the residual entropy at initial state\n",
      "Z_1 = (P_1*V_1_2)/(R*T_1);\n",
      "da_dT_1 = -(a*m)/((alpha*T_1*Tc)**(1./2));\t\t\t#[Pa*m**(6)/mol**(2)] - da/dT\n",
      "\n",
      "# The residual entropy change for Peng-Robinson equatiob of state is given by\n",
      "# S_R = R*math.log(Z-(P*b)/(R*T)) + (da_dT/(2*2**(1/2)*b))*math.log((V+(1+2**(1/2))*b))/((V+(1-2**(1/2)*b))));\n",
      "S_R_1 = R*(math.log(Z_1-(P_1*b_2)/(R*T_1))) + (da_dT_1/(2*2**(1.2)*b_2))*(math.log((V_1_2+(1+2**(1./2))*b_2)/(V_1_2+(1-2**(1./2))*b_2)));\n",
      "\n",
      "# The total entropy change is given by\n",
      "# delta_S = delta_S_ig + delta_S_R\n",
      "\n",
      "def f26(T): \n",
      "    return 27.2+4.2*10**(-3)*T\n",
      "\n",
      "# where, delta_S_ig =  quad(f26,T_1,T_2_2) - R*math.log(P_2_2/P_1)[0]\n",
      "\n",
      "# and, P_2_2 = (R*T_2_2)/(V_2_2-b_2) - a_2/V_2_2**(2)\n",
      "# On simplification we get\n",
      "# delta_S = 27.2*math.log(T_2_2-T_1) + 4.2*10**(-3)*(T_2_2-T_1) - R*math.log(P_2_2/P_1) + R*math.log(Z_2-(P_2_2*b)/(R*T_2_2)) + 6226*(da_dT_2) + 9.22\n",
      "\n",
      "# Now we have the determine the value of T_2_2 such that delta_S = 0\n",
      "# Starting with a temperature of 150 K\n",
      "T_prime = 100.;\t\t\t#[K]\n",
      "error = 10.;\n",
      "while(error>0.1):\n",
      "    Tr_prime = T_prime/Tc;\n",
      "    alpha_prime = (1 + m*(1 - Tr_prime**(1./2)))**(2);\n",
      "    a_prime = ((0.45724*(R*Tc)**(2))/Pc)*alpha_prime;\n",
      "    P_prime = 7.23*10**(4)*T_prime - 3.93*10**(7)*a_prime;\n",
      "    Z_prime = (P_prime*V_2_2)/(R*T_prime);\n",
      "    da_dT_prime = -(a_prime*m)/((alpha_prime*T_prime*Tc)**(1./2));\n",
      "    delta_S = 27.2*math.log(T_prime/T_1) + 4.2*10**(-3)*(T_prime-T_1) - R*math.log(P_prime/P_1) + R*math.log(Z_prime-((P_prime*b_2)/(R*T_prime))) + 6226*(da_dT_prime) + 9.22;\n",
      "    error=abs(delta_S);\n",
      "    T_prime = T_prime + 0.3;\n",
      "\n",
      "T_2_2 = T_prime;\t\t\t#[K] - Final temperature\n",
      "P_2_2 = P_prime*10**(-5);\t\t\t#[bar] - Final pressure\n",
      "\n",
      "print \" 2).The final temperature is %f K\"%(T_2_2);\n",
      "print \"     The final pressure is %f bar\"%(P_2_2);\n",
      "\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 1).The final temperature is 133.131844 K\n",
        "     The final pressure is 39.636659 bar\n",
        " 2).The final temperature is 134.200000 K\n",
        "     The final pressure is 40.130674 bar\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.27  Page Number : 374"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculation of vapour pressure\n",
      "\n",
      "# Variables\n",
      "T = 373.15;\t\t\t#[K]\n",
      "Tc = 562.16;\t\t\t#[K]\n",
      "Pc = 48.98;\t\t\t#[bar]\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "R = 8.314;\t\t\t#[J/mol-K] - Universal gas constant\n",
      "\n",
      "# The cubic form of Redlich Kwong equation of state is given by,\n",
      "# V**(3) - ((R*T)/P)*V**(2) - ((b_1**(2)) + ((b_1*R*T)/P) - (a/(T**(1/2)*P))*V - (a*b)/(T**(1/2)*P) = 0\n",
      "\n",
      "# Calculations\n",
      "a = (0.42748*(R**(2))*(Tc**(2.5)))/Pc;\t\t\t#[Pa*m**(6)*K**(1/2)/mol]\n",
      "b = (0.08664*R*Tc)/Pc;\t\t\t#[m**(3)/mol]\n",
      "\n",
      "# At 373.15 K, let us assume the pressure to be 2.5 bar and under these conditions \n",
      "P_1 = 2.5;\t\t\t#[bar]\n",
      "P_1 = P_1*10**(5);\t\t\t#[bar]\n",
      "\n",
      "# Putting the values in Redlich Kwong equation of state, the equation becomes\n",
      "# V**(3) - 0.0124*V**(2) + 8.326*10**(-6)*V - 7.74*10**(-10) = 0\n",
      "# Solving the cubic equation\n",
      "\n",
      "def f(V): \n",
      "    return V**(3) - 0.0124*V**(2) + 8.326*10**(-6)*V - 7.74*10**(-10)\n",
      "V1=fsolve(f,-9)\n",
      "V2=fsolve(f,10)\n",
      "V3=fsolve(f,0.1)\n",
      "# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.\n",
      "V_liq = V1;\t\t\t#[m**(3)/mol] - Molar volume in liquid phase\n",
      "V_vap = V3;\t\t\t#[m**(3)/mol] - Molar volume in vapour phase\n",
      "\n",
      "# Let us calculate the fugacity of vapour phase\n",
      "# math.log(f_vap/P) = b/(V-b) + math.log((R*T)/(P*(V-b))) - (a/(R*T**(1.5)))*(1/(V+b) - (1/b)*math.log(V/(V+b)))\n",
      "f_vap = P_1*math.exp(b/(V_vap-b) + math.log((R*T)/(P_1*(V_vap-b))) - (a/(R*T**(1.5)))*(1/(V_vap+b) - (1/b)*math.log(V_vap/(V_vap+b))));\t\t\t#[Pa]\n",
      "\n",
      "# Let us calculate the fugacity of the liquid phase\n",
      "f_liq = P_1*math.exp(b/(V_liq-b) + math.log((R*T)/(P_1*(V_liq-b))) - (a/(R*T**(1.5)))*(1/(V_liq+b) - (1/b)*math.log(V_liq/(V_liq+b))));\n",
      "\n",
      "\n",
      "# The two fugacities are not same; therefore another pressure is to be assumed. The new pressure is\n",
      "P_new = P_1*(f_liq/f_vap);\t\t\t#[Pa]\n",
      "\n",
      "# At P_new\n",
      "def f1(V): \n",
      "    return V**(3) - ((R*T)/P_new)*V**(2) - (b**(2) + ((b*R*T)/P_new) - a/(T**(1/2)*P_new))*V - (a*b)/(T**(1/2)*P_new)\n",
      "V4=fsolve(f1,-9)\n",
      "V5=fsolve(f1,10)\n",
      "V6=fsolve(f1,0.1)\n",
      "# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.\n",
      "V_liq_2 = V4;\t\t\t#[m**(3)/mol] - Molar volume in liquid phase\n",
      "V_vap_2 = V6;\t\t\t#[m**(3)/mol] - Molar volume in vapour phase\n",
      "\n",
      "f_vap_prime = P_new*math.exp(b/(V_vap_2-b) + math.log((R*T)/(P_new*(V_vap_2-b))) - (a/(R*T**(1.5)))*(1/(V_vap_2+b) - (1/b)*math.log(V_vap_2/(V_vap_2+b))));\t\t\t#[Pa]\n",
      "f_liq_prime = P_new*math.exp(b/(V_liq_2-b) + math.log((R*T)/(P_new*(V_liq_2-b))) - (a/(R*T**(1.5)))*(1/(V_liq_2+b) - (1/b)*math.log(V_liq_2/(V_liq_2+b))));\n",
      "\n",
      "# Since the fugacities of liquid and vapour phasesare almost same the assumed pressure may be taken as vapour pressure at 373.15 K\n",
      "P_new = P_new*10**(-5);\t\t\t#[bar]\n",
      "\n",
      "# Results\n",
      "print \" The vapour pressure of benzene using Redlich Kwong equation of state is %f bar\"%(P_new);\n",
      "\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The vapour pressure of benzene using Redlich Kwong equation of state is 2.666995 bar\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Example 10.28  Page Number : 374"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Determination of vapour pressure\n",
      "\n",
      "# Variables\n",
      "T = 150 + 273.15;\t\t\t#[K]\n",
      "Tc = 647.1;\t\t\t#[K]\n",
      "Pc = 220.55;\t\t\t#[bar]\n",
      "Pc = Pc*10**(5);\t\t\t#[Pa]\n",
      "w = 0.345;\n",
      "R = 8.314;\t\t\t#[J/mol-K] - Universal gas constant\n",
      "\n",
      "# Let us assume a pressure of 100 kPa.\n",
      "P_1 = 100.*10**(3);\t\t\t#[Pa]\n",
      "\n",
      "# Calculations\n",
      "# At 100 kPa and 423.15 K, from Peng-Robinson equation of stste \n",
      "m = 0.37464 + 1.54226*w - 0.26992*w**(2);\n",
      "Tr = T/Tc;\n",
      "alpha = (1 + m*(1 - Tr**(1./2)))**(2);\n",
      "a = ((0.45724*(R*Tc)**(2))/Pc)*alpha;\t\t\t#[Pa*m**(6)/mol**(2)]\n",
      "b = (0.07780*R*Tc)/Pc;\t\t\t#[m**(3)/mol]\n",
      "# Cubic form of Peng-Robinson equation of stste is given by\n",
      "# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;\n",
      "# Solving the cubic equation\n",
      "def f(V): \n",
      "    return V**(3)+(b-(R*T)/P_1)*V**(2)-((3*b**(2))+((2*R*T*b)/P_1)-(a/P_1))*V+b**(3)+((R*T*(b**(2)))/P_1)-((a*b)/P_1)\n",
      "V1 = fsolve(f,-1)\n",
      "V2 = fsolve(f,0)\n",
      "V3 = fsolve(f,1)\n",
      "# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.\n",
      "V_liq = V1;\t\t\t#[m**(3)/mol] - Molar volume in liquid phase\n",
      "V_vap = V3;\t\t\t#[m**(3)/mol] - Molar volume in vapour phase\n",
      "\n",
      "# The compressibility factor is given by\n",
      "Z_vap = (P_1*V_vap)/(R*T);\t\t\t# For liquid phase\n",
      "Z_liq = (P_1*V_liq)/(R*T);\t\t\t# For vapour phase\n",
      "\n",
      "# The math.expression for fugacity of Peng Robinson equation is\n",
      "# math.log(f/P) = (Z-1) - math.log(Z-((P*b)/(R*T))) - (a/(2*2**(1/2)*b*R*T))*math.log((Z+(1+2**(1/2))*((P*b)/(R*T)))/((Z+(1-2**(1/2))*((P*b)/(R*T)))\n",
      "# For vapour phase\n",
      "f_P_vap = math.exp((Z_vap-1) - math.log(Z_vap-((P_1*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_vap+(1+2**(1./2))*((P_1*b)/(R*T)))/(Z_vap+(1-2**(1./2))*((P_1*b)/(R*T)))));\n",
      "# For liquid phase\n",
      "f_P_liq = math.exp((Z_liq-1) - math.log(Z_liq-((P_1*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_liq+(1+2**(1./2))*((P_1*b)/(R*T)))/(Z_liq+(1-2**(1./2))*((P_1*b)/(R*T)))));\n",
      "\n",
      "# Therefore f_liq/f_vap can be calculated as\n",
      "fL_fV = (f_P_liq/f_P_vap);\n",
      "\n",
      "# The two values (f/P)_vap and (f/P)_vap are not same [ (f_P_liq/f_P_vap) >1 ]; therefore another pressure is to be assumed. The new pressure be\n",
      "P_new = P_1*(f_P_liq/f_P_vap);\t\t\t#[Pa]\n",
      "\n",
      "# At P_new and 423.15 K, from Peng-Robinson equation of stste \n",
      "\n",
      "# V**(3)+(b-(R*T)/P)*V**(2)-((3*b**(2))+((2*R*T*b)/P)-(a/P))*V+b**(3)+((R*T*(b**(2))/P)-((a*b)/P)=0;\n",
      "# Solving the cubic equation\n",
      "def f(V): \n",
      "    return V**(3)+(b-(R*T)/P_new)*V**(2)-((3*b**(2))+((2*R*T*b)/P_new)-(a/P_new))*V+b**(3)+((R*T*(b**(2)))/P_new)-((a*b)/P_new)\n",
      "V4 = fsolve(f,-1)\n",
      "V5 = fsolve(f,0)\n",
      "V6 = fsolve(f,1)\n",
      "# The largest root and the smallest root is considered for liquid phase and vapour phase respectively.\n",
      "V_liq_2 = V4;\t\t\t#[m**(3)/mol] - Molar volume in liquid phase\n",
      "V_vap_2 = V6;\t\t\t#[m**(3)/mol] - Molar volume in vapour phase\n",
      "\n",
      "# The compressibility factor is given by\n",
      "Z_vap_2 = (P_new*V_vap_2)/(R*T);\t\t\t# For liquid phase\n",
      "Z_liq_2 = (P_new*V_liq_2)/(R*T);\t\t\t# For vapour phase\n",
      "\n",
      "# For vapour phase\n",
      "f_P_vap_2 = math.exp((Z_vap_2-1) - math.log(Z_vap_2-((P_new*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_vap_2+(1+2**(1./2))*((P_new*b)/(R*T)))/(Z_vap_2+(1-2**(1./2))*((P_new*b)/(R*T)))));\n",
      "# For liquid phase\n",
      "f_P_liq_2 = math.exp((Z_liq_2-1) - math.log(Z_liq_2-((P_new*b)/(R*T))) - (a/(2*2**(1./2)*b*R*T))*math.log((Z_liq_2+(1+2**(1./2))*((P_new*b)/(R*T)))/(Z_liq_2+(1-2**(1./2))*((P_new*b)/(R*T)))));\n",
      "\n",
      "# Therefore f_liq/f_vap can be calculated as\n",
      "fL_fV_2 = (f_P_liq_2/f_P_vap_2);\n",
      "\n",
      "# And new pressure is given by\n",
      "P_new_prime = P_new*(f_P_liq_2/f_P_vap_2);\t\t\t#[Pa]\n",
      "P_new_prime = P_new_prime*10**(-5);\n",
      "\n",
      "# Since the change in pressure is small, so we can take this to be the vapour pressure at 150 C\n",
      "\n",
      "# Results\n",
      "print \" The vapour pressure of water using Peng-Robinson equation of stste  is %f bar\"%(P_new_prime);\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " The vapour pressure of water using Peng-Robinson equation of stste  is 4.675976 bar\n"
       ]
      }
     ],
     "prompt_number": 24
    }
   ],
   "metadata": {}
  }
 ]
}