usnistgov / REFPROP-issues

A repository solely used for reporting issues with NIST REFPROP
26 stars 13 forks source link

Clarification on computation of entropy and PR equation #692

Open csimha opened 3 days ago

csimha commented 3 days ago

Hello

We work in solid mechanics are new to the use of REFPROP; we are using it to develop models of fracture of pipelines conveying CO2 and trying to understand some of the basics of the thermodynamics. The issue is to compute the decompression curve using, say, the PR EOS (without volume translation).

For example, using 298.15 K and 0.10133 MPa (Enthalpy and Entropy set to 0) as the reference state REFPROP yields entropy of 44.3 kJ/kmol-K for T=281.89 K and P= 11.27MPa

We tried to compute this from known formulas 1 Departure function (eq. 4.4-30 of Sandler, Chemical and Engineering Thermodynamics) gives -34.12 kJ/kmol-k

  1. Integral Cp(T)/T from Tr (298.15 to 281.89) - R*ln(P/Pr) gives -41.22 kJ/kmol-k
  2. Entropy of the ideal gas at Tr, Pr is 119.68 kJ/kmol-k (from REFRPROP)

Sum of 1, 2, and 3 gives 44.3 kj/kmol-k which matches the REFPROP result.

Could you please clarify how the entropy of the ideal gas is computed (item 3)?

We then decrement the temperature in steps of 1, and keep the entropy constant (44.3 kJ/kmol-K). For a temperature of 274.89 K and 4.3720 MPa, REFPROP indicates a "Liquid" phase and the Z value of 0.092, which agrees with Maple output below.

However, the PR Z-equation gives 3 real roots (Maple output below) indicating a two-phase state. Please clarify this result.

f1 := z^3 - 0.9489891422 z^2 + 0.2453568192 z - 0.01538344089

f1solve := fsolve(f1, z); f1solve := 0.09267064158, 0.2965571634, 0.5597613372

image

REFPROP Version 10 Windows 10 direct use

Thanks for your help.

ianhbell commented 3 days ago

Is there any particular reason you are using PR instead of the much more accurate multiparameter EOS?

If you are talking decompression, should that be an adiabatic process, so at constant enthalpy? Or you are assuming you are doing perfect work recovery?

To know where your state point is, you need to consider the phase boundary. Recommend to plot it and overlay your point(s) to know for sure. The density solve is only one part of the problem for phase equilibrium

csimha commented 3 days ago

Hello Thanks for the response. We have developed Fortran code for Span-Wagner and are developing GERG2008 code; these are to be integrated as user EOS into the Abaqus finite element solver for modeling rupture of pipelines.

For beginning graduate students to gain an understanding of the process and assumptions PR is simpler - hence, the exercise.

For decompression, we assume that the initial entropy does not change. Also, in the two-phase states, the pressure and temperature are taken to be the same in the phases. Density is a weighted sum of the two densities through the volume fractions.

The idea is to reproduce the PR curve of the sort shown below (from Botros, K. K., et al. "Measurements of decompression wave speed in pure carbon dioxide and comparison with predictions by equation of state." Journal of Pressure Vessel Technology 138.3 (2016): 031302.). They used REFPROP to generate the decompression speed curves.

image

We have been able to reproduce the curve, but wanted to understand the numbers behind it. Will do the overlay plot. Thanks

ianhbell commented 16 hours ago

You assume that the initial entropy doesn't change, or that the entropy is constant throughout the process (my assumption)?

For two-phase states, molar phase fractions are more clearly defined, so I recommend them, in which case $v=\beta v_v+(1-\beta)v_l$ where $\beta$ is the molar phase fraction of the vapor in this case. It makes a nicer analogy to enthalpy, entropy, etc.

Also, for two-phase states, you need to carry out the complete phase equilibrium solution, you are doing that right?

csimha commented 15 hours ago

Yes, we assume that the decompression process is isentropic. The entropy of the initial state is used to obtain new pressure and volume states, with temperature decremented. This is the reason we were inquiring about computation of the ideal gas entropy in the initial state. The algorithm may be found in Figure 4. of https://www.sciencedirect.com/science/article/pii/S0957582018308413?casa_token=kTziMwUvMcMAAAAA:T81t5I3b2xyItA1yJ-efUUrPVD9T0cco7_wGzCjiEI6miLlwCwi2u6_Nxsbv1T-mRj87MEvePTRc

image

In the case of the PR EOS, z=pv/RT, and there is a cubic in z, whose solution indicates the number of phases; when there are 3 real roots, the phase is two state and when there is only one real root, there is only one phase.

ianhbell commented 12 hours ago

You are also doing the iteration to keep entropy constant? It wasn't clear from your description above. Hopefully you are using derivative-based methods do the rootfinding.