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.

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;

Answer 1

[Torsten]

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

[Hanadi]

How can I have a plot?

[Torsten]

You want to plot negative pressures ?

Ok, give me a list of densities and temperatures of the same length and I'll compute the corresponding "pressures" with your code.

[Hanadi]

I can not garantie the same length. for example i need to make like graph 2 in the attached file ... I hope you can see the file.

[Torsten]

Why can't you guarantee the same length ? You need density and temperature to compute pressure, thus the three arrays must have the same length. Of course, you can repeat temperatures for varying densities, e.g.

rho = 0.01:0.01:14;
temp = 180*ones(size(rho));
p = arrayfun(@(rho,temp)EOS(rho,temp),rho,temp)

p = 1×1400

14.7560 29.0917 43.0071 56.5022 69.5770 82.2318 94.4669 106.2825 117.6792 128.6574 139.2177 149.3609 159.0877 168.3989 177.2954 185.7783 193.8485 201.5072 208.7556 215.5949 222.0265 228.0518 233.6723 238.8894 243.7047 248.1200 252.1368 255.7570 258.9824 261.8149

plot(rho,p)

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

[Hanadi]

Thank you

[Image Analyst]

If this Answer solves your original question, then could you please click the "Accept this answer" link to award the answerer with "reputation points" for their efforts in helping you? They'd appreciate it. Thanks in advance. 🙂 Note: you can only accept one answer (so pick the best one) but you can click the "Vote" icon for as many Answers as you want. Voting for an answer will also award reputation points.

[Hanadi]

Can you run the same code using (for-loop) ?

[Torsten]

rho = 0.01:0.01:14;
temp = 180;
for i = 1:numel(rho)
    p(i) = EOS(rho(i),temp);
end
plot(rho,p)