Question

Using the Peng-Robinson EOS, do the following: Construct a P-V diagram for methane. Use Isotherms in 10-degree increments starting at 75 degrees Kelvin until you are approximately 40 degrees above the critical point. Calculate the vapor pressure of methane as a function of temperature. Compare your computed results with the literature and describe the similarities and differences. On the P-V diagram, use the limits of stability derived in class to construct the metastable regime in your computer code. Use the requirement that fugacity's are equal in the liquid and vapor states to construct the coexistence lines. Compare the envelope of coexistence curves with the result in part (d). ln f^v (T, P)/P = 1/RT integral^V = Z^v RT/P _V = infinity (RT/V - P) DV - ln Z^v + (Z^V - 1) = (Z^V - 1) - ln (Z^v - B) - A/2 Squareroot B ln [Z^v + (1 + Squareroot 2) B/Z^v + (1 - Squareroot 2) B] ln f^L (T, P)/P = 1/RT integral^V = Z^L RT/P _V = infinity (RT/V - P) DV - ln Z^L + (Z^L - 1) = (Z^L - 1) - ln (Z^L - B) - A/2 Squareroot 2B ln [Z^L + (1 + Squareroot 2) B/Z^L + (1 - Squareroot 2) B]
MATLAB code needed! Can garentee lots of up votes. Any questions just comment

Answer 1

%
% function result = PengRobinson(T,P,Tc,Pc,w,MW,Liquido)
% Parameters: T,P,w,Tc,Pc,w,MW,Liquido
% T: Temperature [=] K                                          
% P: Presure [=] Pa                                             
% Tc: critical temperature [=] K                               
% Pc: critical presure [=] Pa                                   
% w: accentic factor
% MW: molar weigth [=] kg/mol
% Liquido:  if Liquido == 1, then calculates liquid fugacity;  
%           if Liquido == 0 then calculates vapor fugacity
% Example:

% [Z fhi density] = PengRobinson(273,2*1.013*1e5,304.21,7.382*1e6,0.225,0.044,1)

function [Z,fhi,density] = PengRobinson(T,P,Tc,Pc,w,MW,Liquido)

R = 8.314; % gas constant [=] J/(mol K)

% Reduced variables
Tr = T/Tc ;
Pr = P/Pc ;

% Parameters of the EOS for a pure component
m = 0.37464 + 1.54226*w - 0.26992*w^2;
alfa = (1 + m*(1 - sqrt(Tr)))^2;
a = 0.45724*(R*Tc)^2/Pc*alfa;
b = 0.0778*R*Tc/Pc;
A = a*P/(R*T)^2;
B = b*P/(R*T);

% Compressibility factor
Z = roots([1 -(1-B) (A-3*B^2-2*B) -(A*B-B^2-B^3)]);

ZR = [];
for i = 1:3
   if isreal(Z(i))
        ZR = [ZR Z(i)];   
   end
end

if Liquido == 1
    Z = min(ZR);   
else
    Z = max(ZR);
end

% Fugacity coefficient
fhi = exp(Z - 1 - log(Z-B) - A/(2*B*sqrt(2))*log((Z+(1+sqrt(2))*B)/(Z+(1-sqrt(2))*B)));
if isreal(fhi)
    density=P*MW/(Z*R*T);
    result = [Z fhi density];
else
    'No real solution for "fhi" is available in this phase'
    result=['N/A' 'N/A' 'N/A'];
end