I need to predict isobaric and Isotropic heat capacity vs. temperature and density from equation 49 and 50 in the attached paper. I don't know how to write these equations.

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;
% z is alph^r (equation 42)
%F=p/(rho*temp*R)-Z;
p = Z*rho*temp*R;
end

Answer 1

[William Rose]

@Hanadi,

[edit: fix typos, and add a question about units, at the end, and change units for Cv in the comments of my code below]

Equations 49 and 50, for Cp/R and Cv/R, involve second derivatives of the dimensionless Helmholtz energy, α.

You have two options:

A. Calculate the analytic second derivatives with respect to τ of equations 41 and 42, and use the formulas you get to write functions for Cv (eq.49) and Cp (eq.50). The right hand sides of eqs.41,42 indicate that the analytic approach will be tedious but not too bad. You just have to be careful.

or

B. Write functions for and , and compute their second derivatives numerically. The form for and could be as follows:

function a0=alpha0(delta,tau)
% alpha0=component of dimensionless Helmholtz energy
% See eq.41 of Lemmon & Jacobsen, JPhysChemRefData 34(1), 2005
% Inputs
%   delta=dimensionless density
%   tau=dimensionless temperature
% Globals used
%   t=vector t(k)
%   a=vector a(k)
%   b=vector b(k)
% Output=a0 [dimensionless]
a0=1; % code for equation 41 goes here
end

and

function ar=alphar(delta,tau,d,t,L,m)
% alpha0=component of dimensionless Helmholtz energy
% See eq.42 of Lemmon & Jacobsen, JPhysChemRefData 34(1), 2005
% Inputs
%   delta=dimensionless density
%   tau=dimensionless temperature
% Globals used
%   d=vector d(k)
%   t=vector t(k)
%   L=vector L(k)
%   m=vector m(k)
% Output=ar [dimensionless]
ar=1; % code for equation 42 goes here
end

Then implement equations 49 and 50 by computing the second derivative of these functions numerically. This could be done with something like the following code.

function Cv=heatCapIsochoric(delta,tau)
% Cv=isochoric heat capacity
% See eq.49 of Lemmon & Jacobsen, JPhysChemRefData 34(1), 2005
% Inputs
%   delta=dimensionless density
%   tau=dimensionless temperature
% Output=Cv [J/degK] or [J/(mol-degK)]
deltaTau=1;  % choose a small delta tau to estimate the second derivative numerically
d2a0=alpha0(tau+deltaTau)-2*alpha0(tau)+alpha0(tau-deltaTau);
d2ar=alphar(tau+deltaTau)-2*alphar(tau)+alphar(tau-deltaTau);
Cv=-R*tau^2*(d2a0/(deltaTau^2)+d2ar/(deltaTau^2)); % equation 49, times R
end

and similar for Cp, eq.50.

The equations above involve the dimensionless temperature τ and dimensionless density δ, so you will have to do the necessary conversions to get Cp and Cv at real world temperature T and density ρ.

Check all my equations and code carefully, since there could be mistakes. I am just trying to give a general idea.

By the way, is there a mistake in the equation for dimensionless temperature, which appears just after eq.40 of Lemmon and Jacbosen? From the published paper:

I would have expected . But it's not my field.

Also, I have a bit of a question about the units. The right hand sides of eq.49,50 should be dimensionless. Therefore the left hand sides also should be dimensionless. The left hand sides are and . Cp and Cv usually have units of [J/degK], but R has units of [J/(mole-degK)]. This would mean the left side has units of moles. Is this a problem? Maybe Cp and Cv are molar heat capacities.

Good luck!

[William Rose]

@Hanadi, I think the answer to my question about the units of equations 49, 50 in Lemmon & Jacobsen (2005) is that Cv and Cp in those equations are the molar heat capacities, and therefore Cp and Cv, in eqs.49,50, have units of [J/(mol-degK)]. Therefore I edited a comment in my code above, where I state the units for Cv.

[Hanadi]

@William Rose Thank you. How can I plot heat capacity Vs. Temperature ?

[William Rose]

@Hanadi,

You could do something like this:

Tc=339.173;                 % critical temp (°K)
rhoc=4.777;                 % citical density (mol/dm^3)
global ak bk dk Lk mk tk;   % global constants (vectors)
% add lines to define the globals
T=Tmin:Tmax;                % temp range (°K)
tau=Tc./T;                  % dimensionless temperature (vector)
                            % Lemmon&Jacobsen(2005); should it be tau=T/Tc?
N=length(T);
Cv=zeros(1,N); Cp=zeros(1:N); % allocate arrays
rho=4;              % density for Cv, Cp calculations (mol/dm^3) 
delta=rho/rhoc;     % dimensionless density
for i=1:N
    Cv(i)=heatCapIsochoric(delta,tau(i));   % compute Cv
    Cp(i)=heatCapIsobaric(delta,tau(i));    % compute Cp
end

Plot results:

figure; plot(T,Cv,'-r.',T,Cp,'-b.');
title('Heat Capacity'); grid on; legend('Cv','Cp');
xlabel('Temperature (K)'); ylabel('C_v, C_p')

Of course you have to write funcitons alpha0() and alphar() which are called by heatCapIsochoric() and by heatCapIsobaric().

Done.

[Hanadi]

I tried, but each time i got error.

I am tring to use these temperature values

223.160 233.159 233.159 243.146 243.156 253.151 253.152 253.154 253.154 263.147 263.149 263.152 263.152 273.146 273.150 273.150 273.151 283.144 283.144 293.141 293.141 293.184 298.144 298.183
303.139 303.139 303.139 303.181 303.181 313.136 313.136 313.136 313.136 313.136 313.136 313.136 313.136 317.148 318.138 319.638 323.133 323.133 323.133 323.133 323.133 323.133 323.250 328.138 328.176 333.130 333.130 333.130 333.130 333.130 333.251 334.252 335.134 336.252 337.133 337.133 338.252

[William Rose]

@Hanadi, please attach the current version of the code (.m files) you are using that produced the results above.

[William Rose]

@Hanadi,

Of the 18 temperature values you cited, three values are repeated twice , one is repeated three times, and the last is repeated four times. Two other values are within one or two millidegrees of other values. If you round to the nearest hundredth of a degree, there are only 7 different values in your list of 18 values. Is there some reason that you are trying to use the same temperatures more than once? Same temperature, different densities?

[Torsten]

@William Rose

I attach the document from which temperature/pressure data are taken. It's Table 1 therein.

[William Rose]

@Torsten,

Thank you. That make sense now: when one does an experiment, it is often useful to collect data multiple times under the same, or nearly the same conditions. Of course for this simulation, which has no stochastic component, the result will be exactly the same when the conditions are the same. Thanks again.

[Hanadi]

Thank you for your help, this my first time I use MATLAB I saw some videos in Youtube, but they were simple. I need to do some plot for my project. So, I have short time as well as so many basic questions. According to your suggestion I need to solve 3 functions, isn't it? the first one is

function a0=alpha0(delta,tau)

end
%% what is the value of delta? 
%% tau is written in code that I'd sent ...
%% then I will save this file. After that what should i write in the next code in order to solve for a0?

[Torsten]

In the meantime, you could have spent 2 hours of your time to pass the online introductory MATLAB course to learn the basics of the language:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

If you have a specific question about code you have written, you can come here and ask.

But you cannot expect that we do what you are supposed to do.