Example - T=380%K Tc=370%K Pc=42 bar R=8314 for i=1:50 if i==1 P(i)=30 else P(i)=P(i-1 1 end Pr(i)= P(i/Pc Tr(i)=T/Tc k(i)=0.37464

Info iconThis preview shows pages 1–3. Sign up to view the full content.

View Full Document Right Arrow Icon
Example No 13 Develop a Matlab code for predicting the state property of a fluid based on Peng-Robinson equation of state. You may refer to your note (pg 28-39) or any thermodynamics textbooks for documentation on Peng-Robinson equation of state. To test your Matlab code, determine the PV relation for propane at 360 K, 370 K and 380K. Plot the P-V diagram for propane for the above temperatures at a range of pressure between 30 – 80 bar Table 1. Data Material Propane Molecular Weight, MW 44.09 kg/kgmol Critical Pressure, Pc 42 bar Critical Temperature, Tc 370 K Acentric factor, ϖ 0.153 Solution Peng-Robinson Equation of State C r T T T = c r P P P = 2 26992 . 0 54226 . 1 37464 . 0 ϖ - ϖ + = κ ( 29 [ ] 2 2 1 r 2 r r T 1 1 T P 45724 . 0 A - κ + = r r T P 07780 . 0 B = ( 29 ( 29 ( 29 0 AB B B Z B 2 B 3 A Z 1 B Z 2 3 2 2 3 = - + + - - + - + specific volume P ZRT = ν Density ν = ρ MW
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Matlab Code: % Equation of State: Peng-Robinson % M=44.09; %kg/kgmol w=0.153; %P=120; % bar
Background image of page 2
Background image of page 3
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: T=380; %K Tc=370; %K Pc=42; % bar R=8314; % for i=1:50 if i==1 P(i)=30; else P(i)=P(i-1)+1; end Pr(i)= P(i)/Pc; Tr(i)=T/Tc; k(i)=0.37464+1.54226*w-0.2692*w^2; A(i)=0.45724*(Pr(i)/Tr(i)^2)*(1+k(i)*(1-Tr(i)^0.5))^2; B(i)=0.07780*(Pr(i)/Tr(i)); a(i)=1; b(i)=(B(i)-1); c(i)=(A(i)-3*B(i)^2-2*B(i)); d(i)=(B(i)^3+B(i)^2-A(i)*B(i)); % % Newton-Raphson xr=0.01; es=0.001; maxit=5000; iter = 0; while (1) xrold = xr; xr = xr - (a(i)*xr^3+b(i)*xr^2+c(i)*xr+d(i))/(3*a(i)*xr^2+2*b(i)*xr+c(i)); iter = iter + 1; if xr ~= 0, ea = abs((xr - xrold)/xr) * 100; end ; if ea <= es | iter >= maxit, break , end ; end Z(i) = xr; % v(i)=(Z(i)*R*T/(P(i)*1e5)); rho(i)=M/v(i); end % plot(v,P) xlabel( 'V(m3/kg)' ) ylabel( 'P(bar)' ) % 20 30 40 50 60 70 80 0.2 0.4 0.6 0.8 1 P(bar) V(m3/kg) Figure xx: PV diagram...
View Full Document

This note was uploaded on 03/21/2010 for the course ENGINEERIN Advance Pr taught by Professor - during the Spring '10 term at Universiti Teknologi Malaysia.

Page1 / 3

Example - T=380%K Tc=370%K Pc=42 bar R=8314 for i=1:50 if i==1 P(i)=30 else P(i)=P(i-1 1 end Pr(i)= P(i/Pc Tr(i)=T/Tc k(i)=0.37464

This preview shows document pages 1 - 3. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online