Help with Thermodynamics equations

[Ángel Martin]

To all resident Thermodynamics experts out there....

Recently I've been documenting a bunch of programs I wrote back at engineering school; you may have noticed the contributions in the Program Library section. The latest one is based on the principle of corresponding states to approximate the PVT properties of a real gas - with more or less accuracy depending on the region.

There's another program that builds on this one to also calculate the enthalpy and entropy values, making the calculations also applicable to saturated vapor and liquid regions. However the results it yields don't really correlate with the values from the substance charts. I suspect the issue lies in the subroutines that evaluate the "discrepancies" of those two magnitudes (i.e. the value differences between the real gas and the ideal gas for each enthalpy and entropy). In particular, the formulas used don't seem to make much sense, as some factors are logarithms of dimensional magnitudes.

Since I reversed-engineered those formulas from the program code, I can't be sure whether the error is at the conceptual stage (i.e. bad formula) or due to errors in the code. Needless to say the little I knew about thermodynamics has faded away since then, about 30 years ago and without any exposure to these subjects.

So here are are the suspects:-

The Equation of State used is a variation on Martin's Cubic in volume:
Pr = Tr / [Zc.Vr – B ] – A / { Tr^n [ Zc.Vr – B + 1/8 ]^2 }

with {Zc, A, B, n} non-dimensional constants specific for each substance.
Volumes are in specific terms, i.e. "per mass" (like m^3/kg)
All reduced magnitudes are non-dimensional (Vr=V/Vc, Tr=T/Tc, Pr=P/Pc).

The Discrepancy of Entropy (DS) is obtained with the formula below:

DS(T,P) = (R/M) { (n+1) + n.A.(M.Z"c/Zc) / [Zc.Vr + 1/8 - B].Tr^(n+1) – Ln(M.Z") + Ln[Vr / (Vr - B/Zc)^(M.Z"c/Zc)] }

where M is the molecular weight; Z"= P.V/R.T ; and Z"c = Pc.Vc/R.Tc

The Discrepancy of Enthalpy (DH) is calculated using the Hemholz's Free Energy Discrepancy (DF), as per the relationship:

DH = DF - U - T.DS
where U ~= P.V - T.R/M is the internal energy.

and the formula for DF shown below:
DF(T,P) = (R.T/M).Ln(M.Z") + Pc.Vc {A / [Zc.Tr^n.(Zc.Vr + 1/8 - B)^2] – Tr.[Ln(Vr)/M.Z"c - Ln(Vr-B/Zc)/Zc] }


Do these equations look ok to you?

I can't really tell where they come from, obviously must be derived from the equation of state but how? Besides I'm not sure they're dimensionally correct...

any feedback is welcome:-

[SlideRule]

(12-31-2016 09:27 AM)Ángel Martin Wrote:  To all resident Thermodynamics experts out there....
any feedback is welcome:-

I graduated from University in '85 w undergraduate degrees in Civil Engineering & Physical Science .. but practiced military engineering for the next 28 years .. so it's been a while. I've followed your 41-program posts without comment but with interest .. I'll bone up (so to speak) on my Physical (or should it be Chemical) Thermodynamics and give it a go!

SlideRule

ps: is this a proper interpretation of your modified Martin's cubic Equation of State?
[attachment=4329]

I'm perusing the Article "Development of an Equation of State for Gases" by Joseph J. Martin & Yu-Ghun Hou pages 142-151 in Vol.1 No.2 of A.I.Ch.E. Journal and am having a bit of difficulty pairing your equations to one of the 60 equations in the article. I guess I'm also far too rusty, it's not clicking yet.
edit:1Jan'17 - replaced P w/ Pr in bmp

[Ángel Martin]

(12-31-2016 01:58 PM)SlideRule Wrote:  ps: is this a proper interpretation of your modified Martin's cubic Equation of State?

I'm perusing the Article "Development of an Equation of State for Gases" by Joseph J. Martin & Yu-Ghun Hou pages 142-151 in Vol.1 No.2 of A.I.Ch.E. Journal and am having a bit of difficulty pairing your equations to one of the 60 equations in the article. I guess I'm also far too rusty, it's not clicking yet.

Yes I think this is the same one - although the pressure in the left term should be Pr reduced, (i.e. Pr = P/Pc) instead of the nominal pressure. You're definitely on the right track - I don't know why that EOS was the selected one for the study back then, but it seemed to strike the teacher's fancy on a number of points. I have a reprint of another Joseph Martin's articles ("Cubic Equations of State: Which?") from I&EC Fundamentals, Vol. 16 page 81 - May 1976. It must be an earlier one because that equation is not covered - the exponent in the reduced temperature must have been a later addition...

The PCS method is very intriguing, even if I understand the accuracy isn't going to be its forte. Still the results should be within some acceptable range, which right now they aren't. There are multiple points where the model can be a fluke, for instance:

1.- it assumes a reference state ("dead state" for entropy) at zero degrees Celsius (273 K), where it assigns 100 kcal/kg for enthalpy and 1 kcal/kg.K for entropy... does this make sense?? Obviously this has a huge impact in the final values, although it would be irrelevant when calculating magnitude changes between two {P,T} conditions.

2.- it uses the Pitzer correlation to calculate the enthalpy of vaporization - but the formula is applied with To = 273 K instead of the boiling temperature (Tb) at the given pressure. It justifies it saying the error introduced by this is small, and the advantage is that there's no need to know the Tb. The formula used is:

Hvap = R.(Tc/M) {10.95 w [(1- Tr)^0.456 + 7.08*(1- Tr)^0.354 ]
with w the acentric factor of the substance.

3.- the internal energy U is expressed as U = P.V- R.T/M . I don't remember seeing this anywhere else so can't tell it's a proper expression. Dimensionally it is correct, though.

4.- and then there is the issue with the equations used to calculate the discrepancies of enthalpy and entropy. Are those correct? What's the definition used to obtained them from the EOS? I can identify some partial derivatives in the expressions, like dPr/dTr(v=cte) embedded into the entropy discrepancy equation...

Attached is the MOD files (includes also the unit conversion MOD that should be present) in case you want to check the results on V41. The program is labeled "FREON". There are five routines for common refrigerants included in the ROM so you can enter the name at the initial prompt to load all the constants automatically. These are "R11", R12", "R13", "R22", and "NH3".

Then enter T and P - with their units - and the program will do its thing, determining the region (superheated, saturated, or liquid) and calculating the results. Turbo mode will help, the execution it's relatively slow. A manual is already posted at TOS (first entry in the 'What's New" section), with the equations listed there as well.

Thanks for looking into this one!

[SlideRule]

Angel

I am perusing your MATLAB™ programs for the calculation of phase equilibrium and other thermodynamic properties using different equations of state (classical cubic equations, cubic equations with excess Gibbs energy mixing rules, group contribution equations and SAFT equations). All programs are open-source and have been designed to be easily reusable using an object-oriented programming methodology. We also provide a detailed documentation that describes the application of the programs. Some of our experiences in the application of these programs in education and research are described in:
Ángel Martín, María Dolores Bermejo, Fidel A. Mato, María José Cocero. Teaching Advanced Equations of State in Applied Thermodynamics Courses Using Matlab. Education for Chemical Engineers 6 (2011) e114-e121. at url:Advanced Equations of State - Advanced Thermodynamics
as well as the numerous references & sources (but specifically your contributions) at url:Useful web sites about thermodynamics
with some anticipation for further personal enlightenment: not there yet. This is an interesting eclectic pursuit in a rigorous field of previous study but current neglect. Thanks for the bump!

BEST!
SlideRule

[Ángel Martin]

Same name, but not the same person! Regardless it's a good source for references and maybe also for the equations I'm looking for, thanks for the pointer!

ÁM

[Ángel Martin]

at last some progress - I've figured out the method used to calculate the discrepancies, which it really is nothing more than applying the definitions in a form that is compatible with the equation of state available:

Dis(F) = R.T.Ln(Z) + ITG{(P – R.T/V).dV} at constant T ; and
Dis(S) = - d[Dis(F)]/dT ; at constant V

The first equation involves a couple of integrals of the type ITG {dx/(ax+b)} and ITG { dx/(ax+b)^2}. Integration limits are the quid, as they define the origins used for the different properties. The final expression is very close to the reversed-engineered formula, except a term of the equation where an error was coded.

The second equation is a bear to handle, laborious but conceptually not complicated - all you need is your theory of derivatives and pay close attention to the details. There seems to be at least a couple of ways to write the final expression (that I have come up with) but neither can be reconciled with the formula I reversed-engineered from the code... so I'm assuming there were errors in the original program.

I also got rid of the "forced" origins for enthalpy and entropy. In summary, the results are now very close to the substance charts within the range of applicability - next step is to figure out the valid ranges!

I removed the MOD files from this thread, will post an update in the software section when ready.

If only I had documented all this when i wrote the program 30 years ago...

Cheers,
ÁM