“HumidAirProps” violates the first law of thermodynamic for isentropic process

[thnttt]

For an isentropic compression process, the relation dh=vdp can be found in any text book. When I calculate the compression work for high pressure humid air, the enthalpy obtained with function “HumidAirProps” seems to produce incorrect results. Below are the codes I used in Mathematica 9.0.1, Win7 x64.

{P1, T1, [Phi]1} = {5000, 303.15, 1}; P2 = 10000; {W, V1, H1, S1} = {HumidAirProps["W", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["Vha", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["Hha", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["S", "P", P1, "T", T1, "R", [Phi]1]}; {T2, V2, H2, S2} = {HumidAirProps["T", "P", P2, "S", S1, "W", W], HumidAirProps["Vha", "P", P2, "S", S1, "W", W], HumidAirProps["Hha", "P", P2, "S", S1, "W", W], S1};

Ws1 = NIntegrate[ HumidAirProps["Vha", "P", p, "S", S1, "W", W], {p, P1, P2}]; Ws2 = H2 - H1;

In[415]:= {Ws1,Ws2} Out[415]= {67.7215,70.9963}

"Ws1" is calculated by numerical integration along the isentropic path, while "Ws2" is the enthalpy difference between the starting and ending point of the compression process. "Ws1" should equal "Ws2". However, this isn't true for any value of relative humidity [Phi]1. Besides, "Ws2" does not vary continuously at relative humidity [Phi]1=0. The large difference between these methods confuses me a lot, and it will lead to large error in the prediction of compressor efficiency.

I tested the function "Props" for dry air, and got the right result "Ws1"="Ws2". So I wonder If there is a bug in "HumidAirProps" or the method described in ASHRAE RP-1485.

Thanks!

[ibell]

Do you end up with the same entropies at the ends of the compression processes when you integrate? If the process is truly isentropic, I suppose you should. Or your numerical integration isn't very accurate. I can't test Mathematica, so if you redid this simple example in python or MATLAB I could help you.

On Thu, Apr 24, 2014 at 4:09 PM, thnttt notifications@github.com wrote:

For an isentropic compression process, the relation dh=vdp can be found in any text book. When I calculate the compression work for high pressure humid air, the enthalpy obtained with function “HumidAirProps” seems to produce incorrect results. Below are the codes I used in Mathematica 9.0.1, Win7 x64.

{P1, T1, [Phi]1} = {5000, 303.15, 1}; P2 = 10000; {W, V1, H1, S1} = {HumidAirProps["W", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["Vha", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["Hha", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["S", "P", P1, "T", T1, "R", [Phi]1]}; {T2, V2, H2, S2} = {HumidAirProps["T", "P", P2, "S", S1, "W", W], HumidAirProps["Vha", "P", P2, "S", S1, "W", W], HumidAirProps["Hha", "P", P2, "S", S1, "W", W], S1};

Ws1 = NIntegrate[ HumidAirProps["Vha", "P", p, "S", S1, "W", W], {p, P1, P2}]; Ws2 = H2 - H1;

In[415]:= {Ws1,Ws2} Out[415]= {67.7215,70.9963}

"Ws1" is calculated by numerical integration along the isentropic path, while "Ws2" is the enthalpy difference between the starting and ending point of the compression process. "Ws1" should equal "Ws2". However, this isn't true for any value of relative humidity [Phi]1. Besides, "Ws2" does not vary continuously at relative humidity [Phi]1=0. The large difference between these methods confuses me a lot, and it will lead to large error in the prediction of compressor efficiency.

I tested the function "Props" for dry air, and got the right result "Ws1"="Ws2". So I wonder If there is a bug in "HumidAirProps" or the method described in ASHRAE RP-1485.

Thanks!

— Reply to this email directly or view it on GitHubhttps://github.com/ibell/coolprop/issues/211 .

[ibell]

Another thought, the enthalpies and entropies are per kg of dry air, perhaps that is throwing things off too

On Thu, Apr 24, 2014 at 4:14 PM, Ian Bell ian.h.bell@gmail.com wrote:

Do you end up with the same entropies at the ends of the compression processes when you integrate? If the process is truly isentropic, I suppose you should. Or your numerical integration isn't very accurate. I can't test Mathematica, so if you redid this simple example in python or MATLAB I could help you.

On Thu, Apr 24, 2014 at 4:09 PM, thnttt notifications@github.com wrote:

For an isentropic compression process, the relation dh=vdp can be found in any text book. When I calculate the compression work for high pressure humid air, the enthalpy obtained with function “HumidAirProps” seems to produce incorrect results. Below are the codes I used in Mathematica 9.0.1, Win7 x64.

{P1, T1, [Phi]1} = {5000, 303.15, 1}; P2 = 10000; {W, V1, H1, S1} = {HumidAirProps["W", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["Vha", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["Hha", "P", P1, "T", T1, "R", [Phi]1], HumidAirProps["S", "P", P1, "T", T1, "R", [Phi]1]}; {T2, V2, H2, S2} = {HumidAirProps["T", "P", P2, "S", S1, "W", W], HumidAirProps["Vha", "P", P2, "S", S1, "W", W], HumidAirProps["Hha", "P", P2, "S", S1, "W", W], S1};

Ws1 = NIntegrate[ HumidAirProps["Vha", "P", p, "S", S1, "W", W], {p, P1, P2}]; Ws2 = H2 - H1;

In[415]:= {Ws1,Ws2} Out[415]= {67.7215,70.9963}

"Ws1" is calculated by numerical integration along the isentropic path, while "Ws2" is the enthalpy difference between the starting and ending point of the compression process. "Ws1" should equal "Ws2". However, this isn't true for any value of relative humidity [Phi]1. Besides, "Ws2" does not vary continuously at relative humidity [Phi]1=0. The large difference between these methods confuses me a lot, and it will lead to large error in the prediction of compressor efficiency.

I tested the function "Props" for dry air, and got the right result "Ws1"="Ws2". So I wonder If there is a bug in "HumidAirProps" or the method described in ASHRAE RP-1485.

Thanks!

— Reply to this email directly or view it on GitHubhttps://github.com/ibell/coolprop/issues/211 .

[thnttt]

Thank you for your fast reply.

In my codes above, the compression process has the same entropy, and the parameters “Hha” and “Vha” are used to get the enthalpy and volume at per kg of humid air.

It seems “Sha” for entropy is invalid, so I have to use “S” per kg of dry air. I don’t think this is the cause of my problem.

Sorry I don’t have Python or MATLAB installed. Instead, I can explain my calculation procedure: At the compressor inlet surface, pressure “P1” and temperature “T1” of humid air with relative humidity “[Phi]1” are known; so humidity ratio “W”, entropy “S1” and enthalpy “Hha1” can be obtained by “HumidAirProps”; the compressor outlet conditions are calculated with a given pressure “P2”, and the same humidity ratio “W” and entropy “S1” as those of inlet.

Let's discuss the accuracy of the numerical integration for compressor work later. Here is another interesting phenomenon I found, and maybe it should be solved first.

Below are the compressor parameters at two different relative humidity “[Phi]1”(0 and 0.0001). Notice the severe enthalpy jump at compressor outlet (Hha, from 85.4 to 90.97). Why?

[ibell]

Lets start off with an even smaller W (say 1e-10) - though I am not sure that is the cause. You are adding such a small amount of water, it shouldn't be a problem... It looks like maybe it is calculating the wrong drybulb temperature, I'll have to see why. If you then use the drybulb temperature you obtain + p + W, do you get the right S again?

Also, Ill add Sha as an input for consistency's sake.

On Fri, Apr 25, 2014 at 4:15 AM, thnttt notifications@github.com wrote:

Thank you for your fast reply.

In my codes above, the compression process has the same entropy, and the parameters “Hha” and “Vha” are used to get the enthalpy and volume at per kg of humid air.

It seems “Sha” for entropy is invalid, so I have to use “S” per kg of dry air. I don’t think this is the cause of my problem.

Sorry I don’t have Python or MATLAB installed. Instead, I can explain my calculation procedure: At the compressor inlet surface, pressure “P1” and temperature “T1” of humid air with relative humidity “[Phi]1” are known; so humidity ratio “W”, entropy “S1” and enthalpy “Hha1” can be obtained by “HumidAirProps”; the compressor outlet conditions are calculated with a given pressure “P2”, and the same humidity ratio “W” and entropy “S1” as those of inlet.

Let's discuss the accuracy of the numerical integration for compressor work later. Here is another interesting phenomenon I found, and maybe it should be solved first.

Below are the compressor parameters at two different relative humidity “[Phi]1”(0 and 0.0001). Notice the severe enthalpy jump at compressor outlet (Hha, from 85.4 to 90.97). Why?

[image: 11]https://cloud.githubusercontent.com/assets/7394920/2796962/cfed3f02-cc1d-11e3-9cbc-db64d215e3a1.jpg [image: 22]https://cloud.githubusercontent.com/assets/7394920/2796963/cfeeb9d6-cc1d-11e3-9477-c51b5482cbfd.jpg

— Reply to this email directly or view it on GitHubhttps://github.com/ibell/coolprop/issues/211#issuecomment-41353132 .

[thnttt]

I think the problem I encountered in the isentropic compression process may be caused by the strange behavior of the Entropy subroutine in "HumidAirProps".

For the dry air at any certain condition (W=0), even adding a very small amount of water vapor (W=1e-10), "HumidAirProps" will show a severe drop of entropy (especially at high pressure and low temperature conditions).

Another interesting thing happens when W is a negative value (W = -1e-10):

In[139]:= HumidAirProps["S","P",100,"T",380,"W",-10.^-10] Out[139]= Indeterminate

In[140]:= HumidAirProps["H","P",100,"T",380,"W",-10.^-10] Out[140]= 107.706

The function of enthalpy and volume are continuous at W=0 and have real values when W<0. I know W<0 is meaningless. I just wonder why the Entropy function is so special.

[ibell]

Interesting. I'll have to take a look and see what is going on. Thanks for the very clear description of the problem.

On Fri, Apr 25, 2014 at 1:27 PM, thnttt notifications@github.com wrote:

I think the problem I encountered in the isentropic compression process may be caused by the strange behavior of the Entropy subroutine in "HumidAirProps".

For the dry air at any certain condition (W=0), even adding a very small amount of water vapor (W=1e-10), "HumidAirProps" will show a severe drop of entropy (especially at high pressure and low temperature conditions).

[image: 33]https://cloud.githubusercontent.com/assets/7394920/2799983/80af98d4-cc6a-11e3-9ac6-264487759db5.jpg [image: 44]https://cloud.githubusercontent.com/assets/7394920/2799984/80b06bf6-cc6a-11e3-902f-a12505ec6b59.jpg

Another interesting thing happens when W is a negative value (W = -1e-10):

In[139]:= HumidAirProps["S","P",100,"T",380,"W",-10.^-10] Out[139]= Indeterminate

In[140]:= HumidAirProps["H","P",100,"T",380,"W",-10.^-10] Out[140]= 107.706

The function of enthalpy and volume are continuous at W=0 and have real values when W<0. I know W<0 is meaningless. I just wonder why the Entropy function is so special.

— Reply to this email directly or view it on GitHubhttps://github.com/ibell/coolprop/issues/211#issuecomment-41382728 .

[thnttt]

Is there any good news about this problem?

I still need to calculate the isentropic compression work, but the strange results yielded by "HumidAirProps" worries me a lot. REFPROP doesn't support humid air, and I don't think the ideal gas assumption is valid for humid air above 5MPa.

So please help me! Thank you very much!

[ibell]

Sorry at the moment, no, I am working hard on a new version of CoolProp and haven't had the chance to take a look yet. Can you do a plot for fixed T,P with varied W? I'm curious to see whether it is truly discontinuous or just a very sharp change.

On Wed, May 21, 2014 at 10:43 AM, thnttt notifications@github.com wrote:

Is there any good news about this problem?

I still need to calculate the isentropic compression work, but the strange results yielded by "HumidAirProps" worries me a lot. REFPROP doesn't support humid air, and I don't think the ideal gas assumption is valid for humid air above 5MPa.

So please help me! Thank you very much!

— Reply to this email directly or view it on GitHubhttps://github.com/ibell/coolprop/issues/211#issuecomment-43728981 .

[thnttt]

Thank you for your attention! Here are the entropy and enthalpy varying with W.

When using "Props" for dry air, the two ways to calculate the isentropic compression work (described above) produce consistent results. However, "HumidAirProps" gives incorrect result even when W=0.

[ibell]

Interesting. The HAProps definition of the enthalpy is not the same as for the pure fluid - you can see that by reading http://rp.ashrae.biz/page/ASHRAE-D-RP-1485-20091216.pdf . It is not entirely surprising behavior really I guess. Perhaps there is still an error here that needs to be remedied though

On Wed, May 21, 2014 at 1:27 PM, thnttt notifications@github.com wrote:

Thank you for your attention! Here are the entropy and enthalpy varying with W.

[image: 1]https://cloud.githubusercontent.com/assets/7394920/3039133/4b0f13c6-e0d8-11e3-8204-bc89a5536fdf.jpg [image: 2]https://cloud.githubusercontent.com/assets/7394920/3039134/4b348ef8-e0d8-11e3-95d7-7db9ce7347d7.jpg

When using "Props" for dry air, the two ways to calculate the isentropic compression work (described above) produce consistent results. However, "HumidAirProps" gives incorrect result even when W=0.

— Reply to this email directly or view it on GitHubhttps://github.com/ibell/coolprop/issues/211#issuecomment-43742782 .

[ibell]

The basic problem is that if you read that document, they have terms that are log(psi_w), which is -infinity at psi_w = 0, psi_w being the mole fraction of water vapor. We have a conditional in our code ( https://github.com/CoolProp/CoolProp/blob/master/src/HumidAirProp.cpp#L774-L778) to avoid this problem. I have emailed the authors before but we didn't come to a nice solution.

On Wed, May 21, 2014 at 1:39 PM, Ian Bell ian.h.bell@gmail.com wrote:

Interesting. The HAProps definition of the enthalpy is not the same as for the pure fluid - you can see that by reading http://rp.ashrae.biz/page/ASHRAE-D-RP-1485-20091216.pdf . It is not entirely surprising behavior really I guess. Perhaps there is still an error here that needs to be remedied though

On Wed, May 21, 2014 at 1:27 PM, thnttt notifications@github.com wrote:

Thank you for your attention! Here are the entropy and enthalpy varying with W.

[image: 1]https://cloud.githubusercontent.com/assets/7394920/3039133/4b0f13c6-e0d8-11e3-8204-bc89a5536fdf.jpg [image: 2]https://cloud.githubusercontent.com/assets/7394920/3039134/4b348ef8-e0d8-11e3-95d7-7db9ce7347d7.jpg

When using "Props" for dry air, the two ways to calculate the isentropic compression work (described above) produce consistent results. However, "HumidAirProps" gives incorrect result even when W=0.

— Reply to this email directly or view it on GitHubhttps://github.com/ibell/coolprop/issues/211#issuecomment-43742782 .

[thnttt]

Thanks. Now I know the discontinuity at W=0 is an inherent character of the Entropy model.

Here is the performance data of a compressor stage calculated with HAProps :

Notice the abnormal isentropic efficiency, which should be less than the polytropic efficiency. Would you give me some advice on how to obtain a reasonable isentropic efficiency for high pressure humid air. Thank you!