jjgomera / iapws

python libray for IAPWS standard calculation of water and steam properties
GNU General Public License v3.0
170 stars 64 forks source link

IAPWS97: is there any way to get multiple solutions? #22

Closed sodynamic closed 8 years ago

sodynamic commented 8 years ago

Hi, For some input values, there may be multiple solutions in the whole IAPWS 97 data. For instance, steam=IAPWS97(T=300, h=120) the steam may be wet steam or unsanturated water. The program will feed no warning when I set this status, and it seems to feed me AttributeError no matter any attribute I want to get. Is there any possible way that I can get every solution of the steam at these input values? Or what should I do to limit the solution at a certain phase? Thanks.

jjgomera commented 8 years ago

Hi, See bug #21 Thanks

sodynamic commented 8 years ago

Hi, Thank you for your reply. Though I believe this should be a different problem. The following code is what I try at my anaconda. ` from iapws import IAPWS97 s1=IAPWS97(P=0.1,h=2600,) print(s1.P,s1.h) 0.1 2600.0 s2=IAPWS97(T=300,h=120) print(s2.T,s1.h) Traceback (most recent call last): File "", line 1, in AttributeError: 'IAPWS97' object has no attribute 'T' print(s1.Liquid.nu,s1.Vapor.nu) 2.949538851800334e-07 2.0698362356268167e-05 print(s2.Liquid.nu,s2.Vapor.nu) Traceback (most recent call last): File "", line 1, in AttributeError: 'IAPWS97' object has no attribute 'Liquid'

`

As you can see above, this problem is not about 2 phases.

And I have tested these two cases at a different software (Chinese Language). For P=0..1 h=2600 case, The result is: Region: Region4; 湿蒸汽或饱和线

压力 P = 0.10000000 MPa 温度 T = 99.61 ℃ 比焓 H = 2600.00 kJ/kg 比熵 S = 7.1577 kJ/(kg.℃) 比容 V = 1.6378152 m^3/kg 干度 X = 0.9668

But for the T=300, h=120 case, the results are: Region: Region4; 湿蒸汽或饱和线

温度 T = 26.40 ℃ 比焓 H = 120.00 kJ/kg 压力 P = 0.00344411 MPa 比熵 S = 0.4179 kJ/(kg.℃) 比容 V = 0.1539467 m^3/kg 干度 X = 0.0038

另外的一个高压点为:

Region:Region1; 低于350℃的未饱和水 温度 T = 26.40 ℃ 比焓 H = 120.00 kJ/kg 压力 P = 10.13572057 MPa 比熵 S = 0.3841 kJ/(kg.℃) 比容 V = 0.0009989 m^3/kg

which is quite different with bug #21. This problem is about 2 solutions, not 2 phases.

sodynamic commented 8 years ago

I have another code that deals with IAPWS 97, which is very old and not compatible with new compiler, and it has three input arguments to set up a steam. The first two inputs, which like yours, are the properties of the steam, such as pressure, temperature, enthalpy and so on. The third input indicates the region of IAPWS that the solution may locate. If the third input is 0, it means find solution in whole IAPWS 97 data. If the third input is 4, it means find the soluntion from and only from the region 4 of IAPWS 97.

I don't know if your code has the same ability. Say, an input argument has a default value 0, which means for mose cases the code finds solution from all regions of IAPWS 97 data. If the user input the specific region number, such as 2 or 3, then this the code will find the solution in the user-defined region.

jjgomera commented 8 years ago

Ok, you're right, the problem it's my code don't implement the pair T,h as a input parameter, so if you used the state isn't defined:

In [1]: from iapws import IAPWS97
In [2]: st = IAPWS97(T=300, h=120)
In [3]: st.status
Out[3]: 0

That pair T,h don't define well the region, if you see a T-h plot you can see region very small with great state change near saturated liquid line, that's a problem for a good convergence in the mathematical resolution. For that reason, the iapws association don't define a backward equation for T-h. It's a bad definition pair.

iT-h plot

That bad convergence can give erroneous multiple solutions, as I suppose happen with the Chinese software you talk about. By the way, 300K is 26.85ºC, not 26.4ºC as that software return. That point (T=300, h=120) I'm sure its in two phase region, and only with a solution.

In [1]: from iapws import IAPWS97
In [2]: st_liquid = IAPWS97(T=300, x=0)
In [3]: st_vapor = IAPWS97(T=300, x=1)
In [4]: st_liquid.h, st_vapor.h
Out[4]: (112.57499081240724, 2549.8930083073265)
In [5]: from iapws.iapws97 import _PSat_T
In [6]: st = IAPWS97(P=_PSat_T(300), h=120)
In [7]: st.h
Out[7]: 120.0
In [8]: st.T
Out[8]: 299.9999999999999
In [9]: st.P
Out[9]: 0.0035365894130130106
In [10]: st.x
Out[10]: 0.003046385056975116

In [11]: st2 = IAPWS97(T=300, x=st.x)
In [12]: st.h
Out[12]: 120.0
In [13]: st.T
Out[13]: 299.9999999999999
In [14]: st.P
Out[14]: 0.0035365894130130106
In [15]: st.x
Out[15]: 0.003046385056975116

If you want to calculate a state in region 4, I recommend you use the T-x or P-x input pair.

sodynamic commented 8 years ago

OK, I think I can use this alternative way. Thanks!