I have a code related to thermodynamic subject, i don't know how to run it in order to get the different plot out of it. for example plot (p,rho). plot of pressure vs density.

3 views (last 30 days)
function F=EOS(rho)
%Based on Lemmon and Jacobsen's HFC-125 EOS paper (2005):
R=8.314472; %(l.kPa)/(mol.K) (universal gas constant)
MM=120.0214; %g/mol (molar mass)
rho_c=4.779; %mol/l (critical density)
temp_c=339.173; %K (critical temperature)
p_c=3617.7; %kPa (critical pressure)
%input pressure and temperature to find the density at:
%p=2.914; %kPa (triple point)
%p=101.325; %kPa (normal point)
%p=16.5*p_c; %kPa (max range of validity)
p=0.2*p_c; %kPa
%temp=172.52; %K (triple point)
%temp=225.061; %K (normal point)
%temp=500; %K (max range of validity)
temp=150; %K
sigma=rho/rho_c; %reduced density
tau=temp_c/temp; %inverse reduced temperature
%The HFC-125 EOS parameters:
%f=[ n_i d_i t_i l_i m_i]
f=[ 5.280760000 1 0.669 0 0.0
-8.676580000 1 1.050 0 0.0
0.750112700 1 2.750 0 0.0
0.759002300 2 0.956 0 0.0
0.014518990 4 1.000 0 0.0
4.777189000 1 2.000 1 0.0
-3.330988000 1 2.750 1 0.0
3.775673000 2 2.380 1 0.0
-2.290919000 2 3.370 1 0.0
0.888826800 3 3.470 1 0.0
-0.623486400 4 2.630 1 0.0
-0.041272630 5 3.450 1 0.0
-0.084553890 1 0.720 2 0.0
-0.130875200 5 4.230 2 0.0
0.008344962 1 0.200 3 0.0
-1.532005000 2 4.500 2 1.7
-0.058836490 3 29.00 3 7.0
0.022966580 5 24.00 3 6.0];
n=f(:,1);
d=f(:,2);
t=f(:,3);
l=f(:,4);
m=f(:,5);
n1=n(1);
n2=n(2);
n3=n(3);
n4=n(4);
n5=n(5);
n6=n(6);
n7=n(7);
n8=n(8);
n9=n(9);
n10=n(10);
n11=n(11);
n12=n(12);
n13=n(13);
n14=n(14);
n15=n(15);
n16=n(16);
n17=n(17);
n18=n(18);
d1=d(1);
d2=d(2);
d3=d(3);
d4=d(4);
d5=d(5);
d6=d(6);
d7=d(7);
d8=d(8);
d9=d(9);
d10=d(10);
d11=d(11);
d12=d(12);
d13=d(13);
d14=d(14);
d15=d(15);
d16=d(16);
d17=d(17);
d18=d(18);
t1=t(1);
t2=t(2);
t3=t(3);
t4=t(4);
t5=t(5);
t6=t(6);
t7=t(7);
t8=t(8);
t9=t(9);
t10=t(10);
t11=t(11);
t12=t(12);
t13=t(13);
t14=t(14);
t15=t(15);
t16=t(16);
t17=t(17);
t18=t(18);
l1=l(1);
l2=l(2);
l3=l(3);
l4=l(4);
l5=l(5);
l6=l(6);
l7=l(7);
l8=l(8);
l9=l(9);
l10=l(10);
l11=l(11);
l12=l(12);
l13=l(13);
l14=l(14);
l15=l(15);
l16=l(16);
l17=l(17);
l18=l(18);
m1=m(1);
m2=m(2);
m3=m(3);
m4=m(4);
m5=m(5);
m6=m(6);
m7=m(7);
m8=m(8);
m9=m(9);
m10=m(10);
m11=m(11);
m12=m(12);
m13=m(13);
m14=m(14);
m15=m(15);
m16=m(16);
m17=m(17);
m18=m(18);
z1=sigma*(d1*n1*sigma^(d1 - 1)*tau^t1 + d2*n2*sigma^(d2 - 1)*tau^t2 + d3*n3*sigma^(d3 - 1)*tau^t3 + d4*n4*sigma^(d4 - 1)*tau^t4 + d5*n5*sigma^(d5 - 1)*tau^t5);
z2=sigma*(d6*n6*sigma^(d6 - 1)*tau^t6*exp(-sigma^l6) + d7*n7*sigma^(d7 - 1)*tau^t7*exp(-sigma^l7) + d8*n8*sigma^(d8 - 1)*tau^t8*exp(-sigma^l8) + d9*n9*sigma^(d9 - 1)*tau^t9*exp(-sigma^l9) + d10*n10*sigma^(d10 - 1)*tau^t10*exp(-sigma^l10) + d11*n11*sigma^(d11 - 1)*tau^t11*exp(-sigma^l11) + d12*n12*sigma^(d12 - 1)*tau^t12*exp(-sigma^l12) + d13*n13*sigma^(d13 - 1)*tau^t13*exp(-sigma^l13) + d14*n14*sigma^(d14 - 1)*tau^t14*exp(-sigma^l14) + d15*n15*sigma^(d15 - 1)*tau^t15*exp(-sigma^l15) - l6*n6*sigma^d6*sigma^(l6 - 1)*tau^t6*exp(-sigma^l6) - l7*n7*sigma^d7*sigma^(l7 - 1)*tau^t7*exp(-sigma^l7) - l8*n8*sigma^d8*sigma^(l8 - 1)*tau^t8*exp(-sigma^l8) - l9*n9*sigma^d9*sigma^(l9 - 1)*tau^t9*exp(-sigma^l9) - l10*n10*sigma^d10*sigma^(l10 - 1)*tau^t10*exp(-sigma^l10) - l11*n11*sigma^d11*sigma^(l11 - 1)*tau^t11*exp(-sigma^l11) - l12*n12*sigma^d12*sigma^(l12 - 1)*tau^t12*exp(-sigma^l12) - l13*n13*sigma^d13*sigma^(l13 - 1)*tau^t13*exp(-sigma^l13) - l14*n14*sigma^d14*sigma^(l14 - 1)*tau^t14*exp(-sigma^l14) - l15*n15*sigma^d15*sigma^(l15 - 1)*tau^t15*exp(-sigma^l15));
z3=sigma*(d16*n16*sigma^(d16 - 1)*tau^t16*exp(-sigma^l16)*exp(-tau^m16) + d17*n17*sigma^(d17 - 1)*tau^t17*exp(-sigma^l17)*exp(-tau^m17) + d18*n18*sigma^(d18 - 1)*tau^t18*exp(-sigma^l18)*exp(-tau^m18) - l16*n16*sigma^d16*sigma^(l16 - 1)*tau^t16*exp(-sigma^l16)*exp(-tau^m16) - l17*n17*sigma^d17*sigma^(l17 - 1)*tau^t17*exp(-sigma^l17)*exp(-tau^m17) - l18*n18*sigma^d18*sigma^(l18 - 1)*tau^t18*exp(-sigma^l18)*exp(-tau^m18));
Z=1+z1+z2+z3;
F=p/(rho*temp*R)-Z;

Accepted Answer

Torsten
Torsten on 3 Dec 2023
Edited: Torsten on 3 Dec 2023
I think there is something wrong in the calculation of the compressibility factor.
p = EOS(2.0,150)
Z = -2.0828
p = -5.1953e+03
function p=EOS(rho,temp)
%Based on Lemmon and Jacobsen's HFC-125 EOS paper (2005):
R=8.314472; %(l.kPa)/(mol.K) (universal gas constant)
MM=120.0214; %g/mol (molar mass)
rho_c=4.779; %mol/l (critical density)
temp_c=339.173; %K (critical temperature)
p_c=3617.7; %kPa (critical pressure)
%input pressure and temperature to find the density at:
%p=2.914; %kPa (triple point)
%p=101.325; %kPa (normal point)
%p=16.5*p_c; %kPa (max range of validity)
%p=0.2*p_c; %kPa
%temp=172.52; %K (triple point)
%temp=225.061; %K (normal point)
%temp=500; %K (max range of validity)
%temp=150; %K
sigma=rho/rho_c; %reduced density
tau=temp_c/temp; %inverse reduced temperature
%The HFC-125 EOS parameters:
%f=[ n_i d_i t_i l_i m_i]
f=[ 5.280760000 1 0.669 0 0.0
-8.676580000 1 1.050 0 0.0
0.750112700 1 2.750 0 0.0
0.759002300 2 0.956 0 0.0
0.014518990 4 1.000 0 0.0
4.777189000 1 2.000 1 0.0
-3.330988000 1 2.750 1 0.0
3.775673000 2 2.380 1 0.0
-2.290919000 2 3.370 1 0.0
0.888826800 3 3.470 1 0.0
-0.623486400 4 2.630 1 0.0
-0.041272630 5 3.450 1 0.0
-0.084553890 1 0.720 2 0.0
-0.130875200 5 4.230 2 0.0
0.008344962 1 0.200 3 0.0
-1.532005000 2 4.500 2 1.7
-0.058836490 3 29.00 3 7.0
0.022966580 5 24.00 3 6.0];
n=f(:,1);
d=f(:,2);
t=f(:,3);
l=f(:,4);
m=f(:,5);
n1=n(1);
n2=n(2);
n3=n(3);
n4=n(4);
n5=n(5);
n6=n(6);
n7=n(7);
n8=n(8);
n9=n(9);
n10=n(10);
n11=n(11);
n12=n(12);
n13=n(13);
n14=n(14);
n15=n(15);
n16=n(16);
n17=n(17);
n18=n(18);
d1=d(1);
d2=d(2);
d3=d(3);
d4=d(4);
d5=d(5);
d6=d(6);
d7=d(7);
d8=d(8);
d9=d(9);
d10=d(10);
d11=d(11);
d12=d(12);
d13=d(13);
d14=d(14);
d15=d(15);
d16=d(16);
d17=d(17);
d18=d(18);
t1=t(1);
t2=t(2);
t3=t(3);
t4=t(4);
t5=t(5);
t6=t(6);
t7=t(7);
t8=t(8);
t9=t(9);
t10=t(10);
t11=t(11);
t12=t(12);
t13=t(13);
t14=t(14);
t15=t(15);
t16=t(16);
t17=t(17);
t18=t(18);
l1=l(1);
l2=l(2);
l3=l(3);
l4=l(4);
l5=l(5);
l6=l(6);
l7=l(7);
l8=l(8);
l9=l(9);
l10=l(10);
l11=l(11);
l12=l(12);
l13=l(13);
l14=l(14);
l15=l(15);
l16=l(16);
l17=l(17);
l18=l(18);
m1=m(1);
m2=m(2);
m3=m(3);
m4=m(4);
m5=m(5);
m6=m(6);
m7=m(7);
m8=m(8);
m9=m(9);
m10=m(10);
m11=m(11);
m12=m(12);
m13=m(13);
m14=m(14);
m15=m(15);
m16=m(16);
m17=m(17);
m18=m(18);
z1=sigma*(d1*n1*sigma^(d1 - 1)*tau^t1 + d2*n2*sigma^(d2 - 1)*tau^t2 + d3*n3*sigma^(d3 - 1)*tau^t3 + d4*n4*sigma^(d4 - 1)*tau^t4 + d5*n5*sigma^(d5 - 1)*tau^t5);
z2=sigma*(d6*n6*sigma^(d6 - 1)*tau^t6*exp(-sigma^l6) + d7*n7*sigma^(d7 - 1)*tau^t7*exp(-sigma^l7) + d8*n8*sigma^(d8 - 1)*tau^t8*exp(-sigma^l8) + d9*n9*sigma^(d9 - 1)*tau^t9*exp(-sigma^l9) + d10*n10*sigma^(d10 - 1)*tau^t10*exp(-sigma^l10) + d11*n11*sigma^(d11 - 1)*tau^t11*exp(-sigma^l11) + d12*n12*sigma^(d12 - 1)*tau^t12*exp(-sigma^l12) + d13*n13*sigma^(d13 - 1)*tau^t13*exp(-sigma^l13) + d14*n14*sigma^(d14 - 1)*tau^t14*exp(-sigma^l14) + d15*n15*sigma^(d15 - 1)*tau^t15*exp(-sigma^l15) - l6*n6*sigma^d6*sigma^(l6 - 1)*tau^t6*exp(-sigma^l6) - l7*n7*sigma^d7*sigma^(l7 - 1)*tau^t7*exp(-sigma^l7) - l8*n8*sigma^d8*sigma^(l8 - 1)*tau^t8*exp(-sigma^l8) - l9*n9*sigma^d9*sigma^(l9 - 1)*tau^t9*exp(-sigma^l9) - l10*n10*sigma^d10*sigma^(l10 - 1)*tau^t10*exp(-sigma^l10) - l11*n11*sigma^d11*sigma^(l11 - 1)*tau^t11*exp(-sigma^l11) - l12*n12*sigma^d12*sigma^(l12 - 1)*tau^t12*exp(-sigma^l12) - l13*n13*sigma^d13*sigma^(l13 - 1)*tau^t13*exp(-sigma^l13) - l14*n14*sigma^d14*sigma^(l14 - 1)*tau^t14*exp(-sigma^l14) - l15*n15*sigma^d15*sigma^(l15 - 1)*tau^t15*exp(-sigma^l15));
z3=sigma*(d16*n16*sigma^(d16 - 1)*tau^t16*exp(-sigma^l16)*exp(-tau^m16) + d17*n17*sigma^(d17 - 1)*tau^t17*exp(-sigma^l17)*exp(-tau^m17) + d18*n18*sigma^(d18 - 1)*tau^t18*exp(-sigma^l18)*exp(-tau^m18) - l16*n16*sigma^d16*sigma^(l16 - 1)*tau^t16*exp(-sigma^l16)*exp(-tau^m16) - l17*n17*sigma^d17*sigma^(l17 - 1)*tau^t17*exp(-sigma^l17)*exp(-tau^m17) - l18*n18*sigma^d18*sigma^(l18 - 1)*tau^t18*exp(-sigma^l18)*exp(-tau^m18));
Z=1+z1+z2+z3
%F=p/(rho*temp*R)-Z;
p = Z*rho*temp*R;
end
  8 Comments

Sign in to comment.

More Answers (0)

Categories

Find more on Mathematics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!