1 Introduction to Chemical Engineering computing : Extension to Python Bruce A. Finlayson Copyright, 2017 This pdf illustrates how to use the programming language Python to solve the problems posed in the book Introduction to Chemical Engineering computing , Bruce A. Finlayson, Wiley (2006-2014). The material mirrors the use of MATLAB in the book, and solves the examples in Chapters 2, 3, 4, and 8. In addition, a new Appendix F gives summaries of many Python commands and examples that can be used independently of any Chemical Engineering applications. Together with the book, Chemical Engineering applications are illustrated using Microsoft Excel , MATLAB , Aspen Plus , Comsol Multiphysics , and Python. This pdf is intended to be used in conjunction with the book, and the treatment using Python mirrors that using MATLAB. Appendix F includes information for setting up Python on your computer. This pdf is free for personal use and is available at , Supplement Using Python.
2 Introduction to Chemical Engineering computing , Chapter 2 with Python Copyright, Bruce A. Finlayson, 2017 2 CHAPTER 2. EQUATIONS OF STATE USING PYTHON SOLVING EQUATIONS OF STATE (SINGLE EQUATION IN ONE UNKNOWN) Nonlinear algebraic equations can be solved using Python, too. First, you have to define the problem to solve by defining a function; then, you check it; finally, you issue a command to solve it. See Appendix F for additional details. Step 1 Import the fsolve program from SciPy and define the function (note the indentation after the declaration def). from scipy . optimize import fsolve def f(x): return x ** 2 - 2 * x - 8 Step 2 Check the function. Issue the command: print(f( )) to get the result: You can easily calculate Eq. ( ) to see that for x = , the function value is Now you know the function f is correct. Note that we used a value for x that meant that every term in the function was important.
3 If we had used x = 1e-5, then the x*x term would be negligible unless we computed it to ten significant figures; hence we wouldn t have checked the entire function. The value x = is not a good choice either, since an incorrect function x- 2*x-8 would return the same value as x*x-2*x-8, hence the error would not be discovered. This is a trivial example, and it is more important for more complicated problems. Step 3 To find the value of x that makes f(x) = 0 in Python, use the fsolve function. In the command window, issue the following command. x = fsolve (f, 0) # one root is at x = print (' The root is % ' % x) This command solves the following problem for x: starting from an initial guess of 0. The answer is You can test the result by saying: print(f(x))) which gives [ ]. Sometimes the function will have more than one solution, and that can be determined only by using the command with a different initial guess.
4 To summarize the steps, step 1 defined the problem you wished to solve and evaluated it for some x, step 2 checked your programming, and step 3 instructed Python to solve the problem. It is tempting to skip the second step checking your programming but remember: if the programming is wrong, you will solve the wrong problem. The last command gives a further check that the zero of the function has been found. When examining the command x = fsolve (f, x0), the f defines which problem to solve, the x0 is your best guess of the solution and fsolve tells Python to vary x, starting from x0 until the f is zero. In all the commands, the f can be replaced by other things, say prob1. The answer can f(x)=0 Using Python in Chapter 2 Copyright, Bruce A. Finlayson, 2017 3 also be put into another variable name: z = x. z = x print(z) In the last example the result is put into the variable z.
5 The options vector allows you to set certain quantities, like the tolerance. Add to the fsolve command: xtol= For the example used above, you can find the other root by running the program with x0 = 3. Multiple roots can be found only if you search for them starting with different guesses. Example of a Chemical Engineering Problem Solved Using Python Find the specific volume of n-butane at 500 K and 18 atm using the Redlich-Kwong equation of state. Step 1 First, you need to define the function that will calculate the f(x), here specvol(v), given the temperature, pressure, and thermodynamic properties. The file is shown below. def specvol(v): # in K, atm, l/gmol # for n-butane Tc = pc = T = p = R = aRK = * (R * Tc) ** 2 / pc aRK = aRK * (Tc / T) ** bRK = * (R * Tc / pc) return p * v ** 3 - R * T * v ** 2 + (aRK - p * bRK ** 2 - R * T * bRK) * v - aRK * bRK This function, called specvol, defines the problem you wish to solve.
6 Step 2 To test the function specvol you issue the command: print(specvol(2)) and get , which is correct. The specvol function causes Python to compute the value of the function specvol when v = 2. You should check these results line by line, especially the calculation of aRK, bRK, and y (just copy the code except for the return statement and calculate aRK and bRK with a Alternatively, you can use the spreadsheet you developed, put in v = and see what f(v) is; it should be the same as in MATLAB since the cubic function and parameters are the same. Step 3 Next you issue the command: v = fsolve(specvol,2) Introduction to Chemical Engineering computing , Chapter 2 with Python Copyright, Bruce A. Finlayson, 2017 4 print(v) and get In specvol the 2 is an initial guess of the answer. To check, you might evaluate the function to find how close to zero f(v) is. print (specvol(v)) and get Of course you expect this to be zero (or very close to zero) because you expect Python to work properly.)
7 If Python can t find a solution, it will tell you. If you use an initial guess of , you might get the specific volume of the liquid rather than the gas. Python gives Another Example of a Chemical Engineering Problem Solved Using Python Next rearrange the Python code to compute the compressibility factor for a number of pressure values. The compressibility factor is defined in Eq. ( ). ( ) For low pressures, where the gas is ideal, the compressibility factor will be close to As the pressure increases, it will change. Thus, the results will indicate the pressure region where the ideal gas is no longer a good assumption. The following code solves for the Redlich-Kwong, Redlich-Kwong-Soave, and Peng-Robinson equations of state and plots the compressibility factor versus pressure as in Figure from scipy . optimize import fsolve import numpy as np import pylab as plt # n-butane Redlich-Kwong, Eq.
8 ( ) def specvolRK(v, p): # in K, atm, l/gmol # for n-butane Tc = pc = T = 500 R = aRK = * (R * Tc) ** 2 / pc aRK = aRK * (Tc / T) ** bRK = * (R * Tc / pc) return p * v ** 3 - R * T * v ** 2 + (aRK - p * bRK ** 2 - R * T * bRK) * v - aRK * bRK Z=pvRTUsing Python in Chapter 2 Copyright, Bruce A. Finlayson, 2017 5 # n-butane Redlich-Kwong-Soave, Eq. ( ) def specvolRKS(v, p): # in K, atm, l/gmol # for n-butane Tc = pc = T = 500 R = acentric = mRKS = + ( - *acentric)*acentric alphaRKS = (1 + mRKS *(1-(T/Tc)** )) ** 2 aRKS = * alphaRKS * (R * Tc) ** 2 / pc bRKS = * (R * Tc / pc) return p * v ** 3 - R * T * v ** 2 + (aRKS - p * bRKS ** 2 - R * T * bRKS) * v - aRKS * bRKS # n-butane Peng-Robinson, Eq. ( ) def specvolPR(v, p): # in K, atm, l/gmol # for n-butane Tc = pc = T = 500 R = acentric = mPR = + ( - *acentric)*acentric alphaPR = (1 + mPR *(1-(T/Tc)** )) ** 2 aPR = * alphaPR * (R * Tc) ** 2 / pc bPR = * (R * Tc / pc) return p*v**3+(bPR*p - R*T)*v**2+(aPR- *p*bPR**2- *R*T*bPR)*v + (p*bPR**3 + R*T*bPR**2-aPR*bPR) T = 500 R = pressure = (1, 27, 5) print(pressure) print(pressure[0]) print(pressure[5]) zcompRK = (6,dtype=float) zcompRKS = (6,dtype=float) zcompPR = (6,dtype=float) print(zcompRK) for i in range(0, 6, 1): p = pressure[i] guess = R*T/p v = fsolve(specvolRK, guess, p) z = p * v / (R * T) zcompRK[i] = z Introduction to Chemical Engineering computing , Chapter 2 with Python Copyright, Bruce A.
9 Finlayson, 2017 6 v = fsolve(specvolRKS,v,p) z = p * v / (R * T) zcompRKS[i] = z v = fsolve(specvolPR,v,p) z = p * v / (R * T) zcompPR[i] = z print(zcompRK) print(zcompRKS) print(zcompPR) plt . plot (pressure,zcompRK,'o-g',label='Redlich-K wong') plt . plot (pressure,zcompRKS,'x-b',label='Redlich- Kwong-Soave') plt . plot (pressure,zcompPR,'s-r',label='Peng-Robi nson') plt . legend(loc='best)') plt . xlabel('Pressure (atm)') plt . ylabel('Z') plt . title ('n-Butane') plt . show() Figure Compressibility factor for n-butane, using Python The first three commands bring in the needed routines fsolve, numpy, and pylab (for plotting. Then there are three definitions of functions that define the equation governing the specific volume, Eq. ( ) and ( ). The main program sets the temperature, provides a vector of 6 pressures, equidistant from 1 to 27; pressure = [1, 6, 11, 16, 21, 26]. The index goes from 0 to 5.
10 The vectors of compressibilites are also set with 6 values for each equation of state, starting with 0, to be filled as the calculations proceed. Then a loop calculation is made for i from 0 to 5. Using Python in Chapter 2 Copyright, Bruce A. Finlayson, 2017 7 For each pressure in turn, the guess for the Redlich-Kwong equation of state is the result from the ideal gas law. The result from the Redlich-Kwong equation of state is used for the guess when solving the Redlich-Kwong-Soave equation of state, and that solution is used as the gues for the Peng-Robinson equation of state. For each i, after the compressibility is found it is put in a vector for that equation of state. Finally the results are plotted, with three curves on one plot. Using Python in Chapter 3 Copyright, Bruce A. Finlayson, 2017 8 CHAPTER 3. VAPOR-LIQUID EQUILIBRIA USING PYTHON Example using Python Before working this example, explore Appendix F so that you have all the program parts needed here.
Show more