Problem Statement: The following polynomial describes the gas behavior.
0 = P*V4 - R*T*V3 + (R*c/T2-R*T*B0+A0)*V2 + ...
+ (R*c*B0/T2+R*T*B0*b-A0*a)*V - R*c*B0*b/T2
where P = pressure
T = temperature
V = specific volume
A0, a, B0, b, c, and R are constants.
Write a FORTRAN function V(P,T) that returns the
specific volume V as a function of P and T. You can assume that
all the constants in the equation (i.e., A0, a,
B0, b, c, and R) are passed into the function through
a common block.
Solution:
c-----------------------------------------------------------------------
function V(P,T)
c-----------------------------------------------------------------------
c Return the specific volume of an nonideal gas
c calculated from Beattie-Bridgemen Equation of State
c V ... specific volume in m^3/kmole=liter/mole (output)
c P ... pressure in kPa (1atm = 101.325kPa) (input)
c T ... temperature in Kelvin (input)
c All the necessary constants are available in a common block
c-----------------------------------------------------------------------
c Beattie-Bridgeman Equation of State (EOS)
c R*T*(1-e) A
c P = ---------*(V+B) - ----
c V**2 V**2
c where
c A = A0*(1-a/V)
c B = B0*(1-b/V)
c e = c/(V*T**3)
c Multiply both sides of the EOS by V**4 gives a 4th degree polynomial
c P*V**4 = R*T*(V-c/T**3)*(V**2+B0*(V-b)) - V*A0*(V-a)
c Collect terms ...
c 0 = P*V**4 - R*T*V**3 + (R*c/T**2-R*T*B0+A0)*V**2 + ...
c + (R*c*B0/T**2+R*T*B0*b-A0*a)*V - R*c*B0*b/T**2
c-----------------------------------------------------------------------
c Programmin Note:
c Call of ZPORC of IMSL to find the polynomial roots.
c-----------------------------------------------------------------------
c Declare variables ----------------------------------------------------
parameter (n=4)
real coeff(n+1)
complex roots(n)
c All species-dependent constants comes through the common block.
common /cblock/A0, a, B0, b, c, R, Temp
c Make up the coefficient vector ---------------------------------------
coeff(5) = P ! V**4 term
coeff(4) = -R*T ! V**3 term
coeff(3) = R*c/T**2 - R*T*B0+A0 ! V**2 term
coeff(2) = R*c*B0/T**2 + R*T*B0*b - A0*a ! V term
coeff(1) = -R*c*B0*b/T**2 ! constant term
c Call a routine to find the roots -------------------------------------
c CALL ZRHQR(coeff, n, rootr, rooti) ! from Numerical Recipes
CALL ZPORC (n, coeff, roots)
c Return a valid root; pick the root near the ideal gas law. -----------
c The 4th root seems to be the right one
v = roots(4)
end
fugacity óP p*V/(R*T) - 1
ln(--------) = ô ------------- dp
P õ0 p
ln(fugacity/P) = integral (from 0 to P) (p*V/(R*T)-1)/p dp (if you don't see the integral symbol above)
Suggestion: If you have difficulty comprehending the question, first numerically integrate for a uninteresting problem, say V(P)=exp(-P). When you finally get what you are doing, switch to the more practical case of V(P)=vbb(P,T)
For those who have not had thermodynamics: the integrand in the above equation is 0 for the ideal gas law, leading to fugacity=P. Thus, fugacity is a measure of the "real concentration" -- the chemical engineering term is "chemical potential". Note that you do not have to know anything about the meaning of fugacity to do this assignment; just evaluate the given integral. And if you already know the term fugacity, now you know how it is calculated.
Solution:
c-----------------------------------------------------------------------
function fugacity(P,T)
c-----------------------------------------------------------------------
c Return the fugacity based on the Beattie-Bridgemen Equation of State
c fugacity ... fugacity of a gas in kPa (output)
c P ... pressure in kPa (1atm = 101.325kPa) (input)
c T ... temperature in Kelvin (input)
c-----------------------------------------------------------------------
c Fugacity is calculated from the following equation.
c ôP
c ln(fugacity/P) = ³ (p*V/(R*T)-1)/p dp
c õ0
c ln(fugacity/P) = integral (from 0 to P) (p*V/(R*T)-1)/p dp
c-----------------------------------------------------------------------
c Programmin Note:
c Call QDAGS of IMSL to integrate
c-----------------------------------------------------------------------
external fugafcn
common /cblock/A0, a, B0, b, c, R, Temp
c Cannot include the function argument in COMMON; so store it in another variable
Temp = T
c Numerically evaluate the integral first, then find the fugacity ------
c call QDAGS(func, a, b, errabs, errrel, s, errest)
c where func ... name of the function to integrate,
c must be in the form func(x)
c a ... lower integration limit (input)
c b ... upper integration limit (input)
c errabs.. absolute error tolerance (input)
c errrel.. relative error tolerance (input)
c s ... estimate of the integral (output)
c errest ... estimated error (output)
call QDAGS(fugafcn, 0., P, 1.E-4, 0., xintegral, errest)
fugacity = exp(xintegral)*P
end
c-----------------------------------------------------------------------
function fugafcn(P)
c-----------------------------------------------------------------------
c Specify the integrand function in fugacity calculation
c
c Fugacity is calculated from the following equation.
c ln(fugacity/P) = integral (from 0 to P) (p*V/(R*T)-1)/p dp
c-----------------------------------------------------------------------
common /cblock/A0, a, B0, b, c, R, T
c Create an output vector whose size matches the input
fugafcn = ( P*V(P,T) / (R*T) - 1.) / P
end
Combined into a stand-alone program.
Solve the same problem in MATLAB.
Test case: At 500atm and 50C, fugacity(500*101.325, 273.15+50)=51292kPa=506.21atm (The given value is what my MATLAB program, more specifically the MATLAB's "quad8" function, gives. Since I do not have an analytical solution to compare my answer, I do not know the exact answer. Nevertheless, I highly suspect the value because 1) it does not agree with the experimental results and 2) from the class demo, e.g., xqdags.m in 3/18's Handout section, I know MATLAB's quad or quad8 cannot handle more demanding integration problems.)
Solution:
Solution:
Re-work with Mathcad. In finding V(P,T), the specific volume as a function of pressure P and temperature T, either we can provide the polynomial coefficients and find the roots, or we can simply solve the equation of state for V. I prefer tackling the problem with the latter approach because it requires less work.
Solution:
|
