6.24 LAB: Vapor Pressure Calculation using the Peng-Robinson Equation of State Introduction The purpose of this lab is to implement the calculation of vapor pressure of a pure liquid substance using a real-gas equation of state. The vapor pressure is the pressure at which a pure substance boils at a specific temperature (conversely, the boiling point temperature is the temperature at which a pure liquid substance boils at a specific pressure). We envision here a function whose inputs are (1) the identity of some pure substance (for example, water, benzene, hexane, ammonia, whatever), and (2) the temperature, and whose output is the vapor pressure of that substance at that temperature. Vapor-pressure calculations routinely performed as part of the design of chemical process units, such a distillation columns and flash separators Relevant Equations The basis of this lab is the fact that vapor pressure is implicit in real-gas equations of state. Before we consider how our program to uncover vapor pressure will eventually work, let's first consider what we mean by real-gas equation of state. The ideal gas law RT PV= nRT or P=_ where V = Tn is an example of an equation of state, but it is not realistic enough to describe vapor-liquid equilibrium (e.g., boiling). The actual interrelationships among pressure, volume, moles, and temperature are much more complicated than the ideal gas law for real substances One industry-standard "real-gas" equation of state is the Peng-Robinson (PR) equation (D.-Y. Peng and D. B. Robinson, Ind. Eng. Chem. Fundam. 1976:15:59) RT a(T) Here, b is a species-specific parameter representing volume-per-molecule and a(T) is a species-specific function of temperature representing intermolecular interactions. For most useful engineering calculations, the PR equation (and other real-gas equations of state - everybody has their own favorite) is cast in the form of a cubic equation in the compressibility factor Z 0 Z3c2Z2 +cZ +co, where PV RT K 0.37464 15422w-0.26992w" A- 0.45724e B 0.07780 P (T. : and Te P T P This may look complicated, but it really only has four parameters: the critical temperature Tc, the critical pressure Pc, the acentricity factor (Omega"), and the gas constant R. The temperature T and pressure P are considered variables, not parameters. In this cubic form, it is clear that the PR equation can have three roots. For certain values of the coefficients co, c1, and c2, there might exist three real roots, while for other values, there might exist only one real root (the other two are complex conjugates of each other). In the case that only one real root is found, there is only one phase (either liquid or vapor) at the specified T and P. However, in the case that three real roots are found then the lowest one, say ZL is associated with a liquid, and the highest one, say Zv, is associated with a vapor. (Here, the middle one is ignored.) Note that we use the following units: P-] bar, T [-] K. The existence of both a liquid root and a vapor root implies liquid-vapor coexistence, but it does not guarantee that P is the vapor pressure at T. To decide whether the liquid and vapor phases signified by the two roots of P-R are in fact in equilibrium with each other, one has to confirm that their fugacities are equal. Fugacity is a thermodynamic property derived by engineers to make phase-equilibrium calculations easy. For substances obeying the P-R equation of state, one can compute the fugacity fat a specific Z using this equation f(T, P) = P exp 2-1-ln(Z-B) 2V2B Z+(1- V2)B Here, the parameters A and B are the same ones used in PR above. So, if we are lucky enough to find real values for ZL and Zy at the requested T and our guess for P we just plug them each into the fugacity equation independently to get two fugacity values, call them f and fv, respectively. We then ask whether or not they have the same value. If so, then our guess P is PVaP(T). If not, we have to have an algorithm to adjust P until we get pvap, signified by f fv The Algorithm The convergent algorithm for computing Pvep at T using P-R is as follows 1. Choose Tc, Pc, and for your substance; 2 Specify temperature T 3. Guess a value for P (1 bar is good a lot of the time) 4. Compute the coefficients Co, C1, and C2 using the equations for A and B that depend on T, P Tc, Pc, and 5. Find the roots of the associated cubic using numpy, roots ( [10, c 2, c 1, caj ) (see below). The vapor-phase root Zv is the first element returned by numpy.roots) (if it is real), and the liquid-phase root Zu is the third element returned by numpy.roots () (if it is real). 6. If both Zv and ZL are real numbers (not complex), then compute the two fugacities using Eqn. 2: fvaf(Zy) and fL-fL). If not, go to step 3. and choose a different initial guess for P (if Z is really small, use a lower guess, if Z is close to 1, use a higher guess) 7. If fv and fi are equal, then P is Pvap, and we're done. If not, adjust the value of P by multiplying it by fu/fy and go to step 4 Using numpy.roots( to Find Roots of a Cubic Part of your solution to this lab will require finding the roots of a cubic equation. The roots() function of the Numpy module makes this very easy. To use it, you should first import the NumPy module. (You may need to install on your own machine it first; don't worry -the python in zyBooks has it already, and if you install Anaconda on your own laptop as recommended, you already have it.) Then, you can use the numpy.roots () function to find the root of any polynomial.roots( takes single array argument containing the coefficients of the polynomial, in a specific order. The first element is the coefficient on the highest power of x, the second is the coefficient on the second highest, and so forth. An array of n elements implies an n-1-order polynomial. Consider the interactive Python session below: >>>from numpy import roots # return an arrayf -1, 1 as the r ots t x^2-1 array ([-1,1.]) # compute r ts f x^3-1 and store in array r >>> print (r[0]) (-0.5+0.866025403784j) >>> print (r [1]) (-0.5-0.866025403784j) >>> print (r [1].real) 0.5 >>> print (r [1].imag) 0.86602540784 >>> print (r [2]) # compute r ts f x^3-6 x^2 + 11 x - 6 = (x-3) (x-2) (x-1); n te order f elements in array r >>>r- roots ([1,-6,11,-6]) >>> print (r[0]) >>> print (r [1]) 2 >>> print (r [2]) Peng-Robinson vapor pressure calculation Written by: YOUR NAME def kappa '" argument (s) go here''': your code here def AB(""argument (s) go here''') : your code here def pr fug (arguments (s) go here'''): your code here def Pvap "arguments (s) go here'): crit0 numiter0 P - Pinit while abs (1.0-crit) > ep and numiter>>from numpy import roots # return an arrayf -1, 1 as the r ots t x^2-1 array ([-1,1.]) # compute r ts f x^3-1 and store in array r >>> print (r[0]) (-0.5+0.866025403784j) >>> print (r [1]) (-0.5-0.866025403784j) >>> print (r [1].real) 0.5 >>> print (r [1].imag) 0.86602540784 >>> print (r [2]) # compute r ts f x^3-6 x^2 + 11 x - 6 = (x-3) (x-2) (x-1); n te order f elements in array r >>>r- roots ([1,-6,11,-6]) >>> print (r[0]) >>> print (r [1]) 2 >>> print (r [2]) Peng-Robinson vapor pressure calculation Written by: YOUR NAME def kappa '" argument (s) go here''': your code here def AB(""argument (s) go here''') : your code here def pr fug (arguments (s) go here'''): your code here def Pvap "arguments (s) go here'): crit0 numiter0 P - Pinit while abs (1.0-crit) > ep and numiter