Ada语言实现-水和水蒸气热力性质IAPWS-IF1997
对外接口文件：IF97-Interfaces.ads
-- File: IF97.ads
with Ada.Numerics.Generic_Elementary_Functions;
package IF97 is
subtype Double is Long_Float;
package Double_Math is new Ada.Numerics.Generic_Elementary_Functions(Double);
use Double_Math;
-- 普通水的气体常数R=0.461526kJ/(kg*K)
R:constant:=0.461_526;
-- 临界温度
Tc:constant:=647.096;
-- 临界压力
Pc:constant:=22.064;
-- 临界密度ρc=322
Rho_c:constant:=322.0;
function "**"(a,b:double) return Double;
pragma import(c,"**","pow");
end IF97;
-----------------------------------------------------------------------------------------
-- File: IF97-Region1.ads
package IF97.Region1 is
--
-- 1区为常规水,T in 273.15..623.15K,p in Region4.Eq.(30) Ps(T)..100.0MPa
--
--
-- 下列方程参数压力单位为MPa,温度单位为K
--
-- Gibbs自由能
function g(p,T:Double) return Double;
-- γ(p,T)=g(p,T)/RT
function gamma(p,T:Double) return Double;
-- 比容
function v(p,T:Double) return Double;
-- 比内能
function u(p,T:Double) return Double;
-- 比熵
function s(p,T:Double) return Double;
-- 比焓
function h(p,T:Double) return Double;
-- 比等压热容
function cp(p,T:Double) return Double;
-- 等容热容
function cv(p,T:Double) return Double;
-- 声速
function w(p,T:Double) return Double;
--
-- 导出方程
--
-- 根据压力与焓求温度
function T_ph(p,h:Double) return Double;
-- 根据压力与熵求温度
function T_ps(p,s:Double) return Double;
end IF97.Region1;
-----------------------------------------------------------------------------------
-- File: IF97-Region1.adb
package body IF97.Region1 is
package Table1 is
I:constant array(1..34) of Double:=(1..8=>0.0,
9..14=>1.0,
15..19=>2.0,
20..22=>3.0,
23..25=>4.0,
26=>5.0,
27..28=>8.0,
29=>21.0,
30=>23.0,
31=>29.0,
32=>30.0,
33=>31.0,
34=>32.0);
J:constant array(1..34) of Double:=(1=>-2.0,
2=>-1.0,
3=>0.0,
4=>1.0,
5=>2.0,
6=>3.0,
7=>4.0,
8=>5.0,
9=>-9.0,
10=>-7.0,
11=>-1.0,
12=>0.0,
13=>1.0,
14=>3.0,
15=>-3.0,
16=>0.0,
17=>1.0,
18=>3.0,
19=>17.0,
20=>-4.0,
21=>0.0,
22=>6.0,
23=>-5.0,
24=>-2.0,
25=>10.0,
26=>-8.0,
27=>-11.0,
28=>-6.0,
29=>-29.0,
30=>-31.0,
31=>-38.0,
32=>-39.0,
33=>-40.0,
34=>-41.0);
n:constant array(1..34) of Double:=(0.146_329_712_131_67,
-0.845_481_871_691_14,
-0.375_636_036_720_40e1,
0.338_551_691_683_85e1,
-0.957_919_633_878_72,
0.157_720_385_132_28,
-0.166_164_171_995_01e-1,
0.812_146_299_835_68e-3,
0.283_190_801_238_04e-3,
-0.607_063_015_658_71e-3,
-0.189_900_682_184_19e-1,
-0.325_297_487_705_05e-1,
-0.218_417_171_754_14e-1,
-0.528_383_579_699_30e-4,
-0.471_843_210_732_67e-3,
-0.300_017_807_930_26e-3,
0.476_613_939_069_87e-4,
-0.441_418_453_308_46e-5,
-0.726_949_962_975_94e-15,
-0.316_796_448_450_54e-4,
-0.282_707_979_853_12e-5,
-0.852_051_281_201_03e-9,
-0.224_252_819_080_00e-5,
-0.651_712_228_956_01e-6,
-0.143_417_299_379_24e-12,
-0.405_169_968_601_17e-6,
-0.127_343_017_416_41e-8,
-0.174_248_712_306_34e-9,
-0.687_621_312_955_31e-18,
0.144_783_078_285_21e-19,
0.263_357_816_627_95e-22,
-0.119_476_226_400_71e-22,
0.182_280_945_814_04e-23,
-0.935_370_872_924_58e-25
);
end Table1;
package Table2 is
I:constant array(1..20) of Double:=(1..6=>0.0,
7..13=>1.0,
14..15=>2.0,
16..17=>3.0,
18=>4.0,
19=>5.0,
20=>6.0
);
J:constant array(1..20) of Double:=(0.0,1.0,2.0,6.0,22.0,32.0,
0.0,1.0,2.0,3.0,4.0,10.0,
32.0,10.0,32.0,10.0,
32.0,32.0,32.0,32.0
);
n:constant array(1..20) of Double:=(-0.238_724_899_245_21e3,
+0.404_211_886_379_45e3,
+0.113_497_468_817_18e3,
-0.584_576_160_480_39e1,
-0.152_854_824_131_40e-3,
-0.108_667_076_953_77e-5,
-0.133_917_448_726_02e2,
+0.432_110_391_835_59e2,
-0.540_100_671_705_06e2,
+0.305_358_922_039_16e2,
-0.659_647_494_236_38e1,
+0.939_654_008_783_63e-2,
+0.115_736_475_053_40e-6,
-0.258_586_412_820_73e-4,
-0.406_443_630_847_99e-8,
+0.664_561_861_916_35e-7,
+0.806_707_341_030_27e-10,
-0.934_777_712_139_47e-12,
+0.582_654_420_206_01e-14,
-0.150_201_859_535_03e-16
);
end Table2;
package Table3 is
I:constant array(1..20) of Double:=(1..6=>0.0,
7..12=>1.0,
13..17=>2.0,
18..19=>3.0,
20=>4.0
);
J:constant array(1..20) of Double:=(0.0,1.0,2.0,3.0,11.0,31.0,
0.0,1.0,2.0,3.0,12.0,31.0,
0.0,1.0,2.0,9.0,31.0,10.0,
32.0,32.0
);
n:constant array(1..20) of Double:=(+0.174_782_680_583_07e3,
+0.348_069_308_928_73e2,
+0.652_925_849_784_55e1,
+0.330_399_817_754_89,
-0.192_813_829_231_96e-6,
-0.249_091_972_445_73e-22,
-0.261_076_364_893_32,
+0.225_929_659_815_86,
-0.642_564_633_952_26e-1,
+0.788_762_892_705_26e-2,
+0.356_721_106_073_66e-9,
+0.173_324_969_948_95e-23,
+0.566_089_006_548_37e-3,
-0.326_354_831_397_17e-3,
+0.447_782_866_906_32e-4,
-0.513_221_569_085_07e-9,
-0.425_226_570_422_07e-25,
+0.264_004_413_606_89e-12,
+0.781_246_004_597_23e-28,
-0.307_321_999_036_68e-30
);
end Table3;
-- Gibbs自由能
function g(p:Double;T:Double) return Double is
use Table1;
Pi:Double:=p/16.53;
tau:Double:=T/1386.0;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum+n(L)*((7.1-Pi)**I(L))*((tau-1.222)**J(L));
end loop;
return sum*R*T;
end g;
function gamma(p:Double;T:Double) return Double is
use Table1;
Pi:Double:=p/16.53;
tau:Double:=1386.0/T;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum+n(L)*((7.1-Pi)**I(L))*((tau-1.222)**J(L));
end loop;
return sum;
end gamma;
function Gamma_Pi(p,T:Double) return Double is
use Table1;
Pi:Double:=p/16.53;
tau:Double:=1386.0/T;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum-(n(L))*I(L)*((7.1-Pi)**(I(L)-1.0))*((tau-1.222)**J(L));
end loop;
return sum;
end Gamma_Pi;
function Gamma_Pi_Pi(p,T:Double) return Double is
use Table1;
Pi:Double:=p/16.53;
tau:Double:=1386.0/T;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum+n(L)*I(L)*(I(L)-1.0)*((7.1-Pi)**(I(L)-2.0))*((tau-1.222)**J(L));
end loop;
return sum;
end Gamma_Pi_Pi;
function Gamma_Pi_Tau(p,T:Double) return Double is
use Table1;
Pi:Double:=p/16.53;
tau:Double:=1386.0/T;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum-n(L)*I(L)*((7.1-Pi)**(I(L)-1.0))*J(L)*((tau-1.222)**(J(L)-1.0));
end loop;
return sum;
end Gamma_Pi_Tau;
function Gamma_tau(p,T:Double) return Double is
use Table1;
Pi:Double:=p/16.53;
tau:Double:=1386.0/T;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum+n(L)*((7.1-Pi)**I(L))*J(L)*((tau-1.222)**(J(L)-1.0));
end loop;
return sum;
end Gamma_tau;
function Gamma_tau_tau(p,T:Double) return Double is
use Table1;
Pi:Double:=p/16.53;
tau:Double:=1386.0/T;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum+n(L)*((7.1-Pi)**I(L))*J(L)*(J(L)-1.0)*((tau-1.222)**(J(L)-2.0));
end loop;
return sum;
end Gamma_tau_tau;
function v(p,T:Double) return Double is
Gp:Double:=Gamma_Pi(p,T);
begin
return ((p/16.53)*Gp*R*1000.0*T/(p*1.0e6));
end v;
function u(p,T:Double) return Double is
Gp:Double:=Gamma_Pi(p,T);
Gt:Double:=Gamma_tau(p,T);
begin
return ((1386.0/T)*Gt-(p/16.53)*Gp)*R*T;
end u;
function s(p,T:Double) return Double is
G:Double:=gamma(p,T);
Gt:Double:=Gamma_tau(p,T);
begin
return ((1386.0/T)*Gt-G)*R;
end s;
function h(p,T:Double) return Double is
Gt:Double:=Gamma_tau(p,T);
begin
return ((1386.0/T)*Gt)*R*T;
end h;
function cp(p,T:Double) return Double is
Gtt:Double:=Gamma_tau_tau(p,T);
begin
return (-(1386.0/T)**2.0)*Gtt*R;
end cp;
function cv(p,T:Double) return Double is
Gp:Double:=Gamma_Pi(p,T);
Gpp:Double:=Gamma_Pi_Pi(p,T);
Gpt:Double:=Gamma_Pi_Tau(P,T);
Gtt:Double:=Gamma_tau_tau(p,T);
tau:Double:=1386.0/T;
begin
return (-(tau**2.0)*Gtt+(Gp-tau*Gpt)**2.0/Gpp)*R;
end cv;
function w(p,T:Double) return Double is
Gp:Double:=Gamma_Pi(p,T);
Gpp:Double:=Gamma_Pi_Pi(p,T);
Gpt:Double:=Gamma_Pi_Tau(P,T);
Gtt:Double:=Gamma_tau_tau(p,T);
tau:Double:=1386.0/T;
begin
return (R*1000.0*T*Gp**2.0/(((Gp-tau*Gpt)**2.0/(tau**2.0*Gtt)-Gpp)))**0.5;
end w;
function T_ph(p,h:Double) return Double is
use Table2;
pi:Double:=p/1.0;
eta:Double:=h/2500.0;
sum:Double:=0.0;
begin
for L in 1..20 loop
sum:=sum+n(L)*(pi**I(L))*((eta+1.0)**J(L));
end loop;
return sum*1.0;
end T_ph;
function T_ps(p,s:Double) return Double is
use Table3;
pi:Double:=p/1.0;
sigma:Double:=s/1.0;
sum:Double:=0.0;
begin
for L in 1..20 loop
sum:=sum+n(L)*(pi**I(L))*((sigma+2.0)**J(L));
end loop;
return sum*1.0;
end T_ps;
end IF97.Region1;
----------------------------------------------------------------
-- File: IF97-B23.ads
-- IF97 2区与3区边界
package IF97.B23 is
-- 饱和蒸汽与过热蒸汽分界
-- 区2与区3边界为简单的压力与温度关系
-- 压力单位1MPa,温度单位K
-- Eq.(5) 根据给定的温度求区2与区3边界上的压力
function P(T:double) return double with Pre=>(T>=623.15);
-- Eq.(6) 根据给定的压力求区2与区3边界上的温度
function T(P:double) return double with Pre=>(P>=16.5291642526);
end IF97.B23;
---------------------------------------------------------------------
-- File: IF97-B23.adb
package body IF97.B23 is
n1:constant:=+0.348_051_856_289_69e+3;
n2:constant:=-0.116_718_598_799_75e+1;
n3:constant:=+0.101_929_700_393_26e-2;
n4:constant:=+0.572_544_598_627_46e+3;
n5:constant:=+0.139_188_397_788_70e+2;
function P(T:double) return double is
begin
return n1+n2*T+n3*T**2;
end P;
function T(P:double) return double is
use Double_Math;
begin
return n4+((P-n5)/n3)**0.5;
end T;
end IF97.B23;
-------------------------------------------------------------------
-- File: IF97-Region4.ads
package IF97.Region4 is
-- Eq.(30) 根据饱和温度求饱和压力
function Ps(T:Double) return Double with Pre=>T<=647.096;
-- Eq.(31) 根据饱和压力求饱和温度
function Ts(p:Double) return Double with Pre=>p<=22.064;
end IF97.Region4;
-------------------------------------------------------------------
-- File: IF97-Region4.adb
package body IF97.Region4 is
package Table34 is
n:constant array(1..10) of Double:=(+0.116_705_214_527_67e4,
-0.724_213_167_032_06e6,
-0.170_738_469_400_92e2,
+0.120_208_247_024_70e5,
-0.323_255_503_223_33e7,
+0.149_151_086_135_30e2,
-0.482_326_573_615_91e4,
+0.405_113_405_420_57e6,
-0.238_555_575_678_49,
+0.650_175_348_447_98e3
);
end Table34;
function Ps(T:Double) return Double is
use Table34;
Thita:Double:=T+n(9)/(T-n(10));
A:Double:=Thita**2+n(1)*Thita+n(2);
B:Double:=n(3)*Thita**2+n(4)*Thita+n(5);
C:Double:=n(6)*Thita**2+n(7)*Thita+n(8);
begin
return (2.0*C/(-B+(B**2-4.0*A*C)**0.5))**4;
end Ps;
-- 根据饱和压力求饱和温度
function Ts(p:Double) return Double is
use Table34;
beta:Double:=p**0.25;
E:Double:=beta**2+n(3)*beta+n(6);
F:Double:=n(1)*beta**2+n(4)*beta+n(7);
G:Double:=n(2)*beta**2+n(5)*beta+n(8);
D:Double:=2.0*G/(-F-(F**2-4.0*E*G)**0.5);
begin
return (n(10)+D-((n(10)+D)**2-4.0*(n(9)+n(10)*D))**0.5)/2.0;
end Ts;
end IF97.Region4;
------------------------------------------------------------------------
-- File: IF97-Region2.ads
package IF97.Region2 is
--
-- 以下的基本方程均考虑了单相稳定区与亚稳定区
--
-- Gibbs自由能
function g(p,T:Double) return Double;
function gamma(p,T:Double) return Double;
-- 比体积
function v(p,T:Double) return Double;
-- 比内能
function u(p,T:Double) return Double;
-- 比熵
function s(p,T:Double) return Double;
-- 比焓
function h(p,T:Double) return Double;
-- 比等压热容
function cp(p,T:Double) return Double;
-- 比等体积热容
function cv(p,T:Double) return Double;
-- 声速
function w(p,T:Double) return Double;
--
-- 导出方程
--
-- 根据压力与焓求温度
function T_ph(p,h:Double) return Double;
-- 根据压力与熵求温度
function T_ps(p,s:Double) return Double;
end IF97.Region2;
--------------------------------------------------------------------
-- File: IF97-Region2.adb
with IF97.B23; use IF97.B23;
with IF97.Region4; use IF97.Region4;
package body IF97.Region2 is
package Table10 is
J:constant array(1..9) of Double:=(0.0,1.0,-5.0,-4.0,-3.0,-2.0,-1.0,2.0,3.0);
n:constant array(1..9) of Double:=(-0.969_276_865_002_17e1,
+0.100_866_559_680_18e2,
-0.560_879_112_830_20e-2,
+0.714_527_380_814_55e-1,
-0.407_104_982_239_28,
+0.142_408_191_714_44e1,
-0.438_395_113_194_50e1,
-0.284_086_324_607_72,
+0.212_684_637_533_07e-1
);
n2:constant array(1..9) of Double:=(-0.969_372_683_930_49e1,
+0.100_872_759_700_06e2,
-0.560_879_112_830_20e-2,
+0.714_527_380_814_55e-1,
-0.407_104_982_239_28,
+0.142_408_191_714_44e1,
-0.438_395_113_194_50e1,
-0.284_086_324_607_72,
+0.212_684_637_533_07e-1
);
end Table10;
package Table11 is
I:constant array(1..43) of Double:=(1..5=>1.0,
6..10=>2.0,
11..15=>3.0,
16..18=>4.0,
19=>5.0,
20..22=>6.0,
23..25=>7.0,
26..27=>8.0,
28=>9.0,
29..31=>10.0,
32..33=>16.0,
34=>18.0,
35..37=>20.0,
38=>21.0,
39=>22.0,
40=>23.0,
41..43=>24.0
);
J:constant array(1..43) of Double:=(0.0,1.0,2.0,3.0,6.0,
1.0,2.0,4.0,7.0,36.0,
0.0,1.0,3.0,6.0,35.0,
1.0,2.0,3.0,7.0,3.0,
16.0,35.0,0.0,11.0,25.0,
8.0,36.0,13.0,4.0,10.0,
14.0,29.0,50.0,57.0,20.0,
38.0,48.0,21.0,53.0,39.0,
26.0,40.0,58.0
);
n:constant array(1..43) of Double:=(-0.177_317_424_732_13e-2,
-0.178_348_622_923_58e-1,
-0.459_960_136_963_65e-1,
-0.575_812_590_834_32e-1,
-0.503_252_787_279_30e-1,
-0.330_326_416_702_03e-4,
-0.189_489_875_163_15e-3,
-0.393_927_772_433_55e-2,
-0.437_972_956_505_73e-1,
-0.266_745_479_140_87e-4,
+0.204_817_376_923_09e-7,
+0.438_706_672_844_35e-6,
-0.322_776_772_385_70e-4,
-0.150_339_245_421_48e-2,
-0.406_682_535_626_49e-1,
-0.788_473_095_593_67e-9,
+0.127_907_178_522_85e-7,
+0.482_253_727_185_07e-6,
+0.229_220_763_376_61e-5,
-0.167_147_664_510_61e-10,
-0.211_714_723_213_55e-2,
-0.238_957_419_341_04e2,
-0.590_595_643_242_70e-17,
-0.126_218_088_991_01e-5,
-0.389_468_424_357_39e-1,
+0.112_562_113_604_59e-110,
-0.823_113_408_979_98e1,
+0.198_097_128_020_88e-7,
+0.104_069_652_101_74e-18,
-0.102_347_470_959_29e-12,
-0.100_181_793_795_11e-8,
-0.808_829_086_469_85e-10,
+0.106_930_318_794_09,
-0.336_622_505_741_71,
+0.891_858_453_554_21e-24,
+0.306_293_168_762_32e-12,
-0.420_024_676_982_08e-5,
-0.590_560_296_856_39e-25,
+0.378_269_476_134_57e-5,
-0.127_686_089_346_81e-14,
+0.730_876_105_950_61e-28,
+0.554_147_153_507_78e-16,
-0.943_697_072_412_10e-6
);
end Table11;
package Table16 is
I:constant array(1..13) of Double:=(1..4=>1.0,
5..7=>2.0,
8..9=>3.0,
10..11=>4.0,
12..13=>5.0
);
J:constant array(1..13) of Double:=(0.0,2.0,5.0,11.0,
1.0,7.0,16.0,4.0,
16.0,7.0,10.0,9.0,
10.0
);
n:constant array(1..13) of Double:=(-0.733_622_601_865_06e-2,
-0.882_238_319_431_46e-1,
-0.723_345_552_132_45e-1,
-0.408_131_785_344_55e-2,
+0.200_978_033_802_07e-2,
-0.530_459_218_986_42e-1,
-0.761_904_090_869_70e-2,
-0.634_980_376_573_13e-2,
-0.860_430_930_285_88e-1,
+0.753_215_815_227_70e-2,
-0.792_383_754_461_39e-2,
-0.228_881_607_784_47e-3,
-0.264_565_014_828_10e-2
);
end Table16;
package Table19 is
n:constant array(1..5) of Double:=(+0.905_842_785_147_23e3,
-0.679_557_863_992_41,
+0.128_090_027_301_36e-3,
+0.265_265_719_084_28e4,
+0.452_575_789_059_48e1
);
end Table19;
-- Table20用于子区域2a
package Table20 is
I:constant array(1..34) of Double:=(1..6=>0.0,
7..15=>1.0,
16..23=>2.0,
24..25=>3.0,
26..28=>4.0,
29..31=>5.0,
32..33=>6.0,
34=>7.0
);
J:constant array(1..34) of Double:=(0.0,1.0,2.0,3.0,7.0,20.0,
0.0,1.0,2.0,3.0,7.0,9.0,
11.0,18.0,44.0,0.0,2.0,
7.0,36.0,38.0,40.0,42.0,
44.0,24.0,44.0,12.0,32.0,
44.0,32.0,36.0,42.0,34.0,44.0,
28.0
);
n:constant array(1..34) of Double:=(+0.108_989_523_182_88e4,
+0.849_516_544_955_35e3,
-0.107_817_480_918_26e3,
+0.331_536_548_012_63e2,
-0.742_320_167_902_48e1,
+0.117_650_487_243_56e2,
+0.184_457_493_557_90e1,
-0.417_927_005_496_24e1,
+0.624_781_969_358_12e1,
-0.173_445_631_081_14e2,
-0.200_581_768_620_96e3,
+0.271_960_654_737_96e3,
-0.455_113_182_858_18e3,
+0.309_196_886_047_55e4,
+0.252_266_403_578_72e6,
-0.617_074_228_683_39e-2,
-0.310_780_466_295_83,
+0.116_708_730_771_07e2,
+0.128_127_984_040_46e9,
-0.985_549_096_232_76e9,
+0.282_245_469_730_02e10,
-0.359_489_714_107_03e10,
+0.172_273_499_131_97e10,
-0.135_513_342_407_75e5,
+0.128_487_346_646_50e8,
+0.138_657_242_832_26e1,
+0.235_988_325_565_14e6,
-0.131_052_365_450_54e8,
+0.739_998_354_747_66e4,
-0.551_966_970_300_60e6,
+0.371_540_859_962_33e7,
+0.191_277_292_396_60e5,
-0.415_351_648_356_34e6,
-0.624_598_551_925_07e2
);
end Table20;
package Table21 is
I:constant array(1..38) of Double:=(1..8=>0.0,
9..16=>1.0,
17..20=>2.0,
21..24=>3.0,
25..30=>4.0,
31..33=>5.0,
34=>6.0,
35..36=>7.0,
37..38=>9.0
);
J:constant array(1..38) of Double:=(0.0,1.0,2.0,12.0,18.0,24.0,28.0,40.0,
0.0,2.0,6.0,12.0,18.0,24.0,28.0,40.0,
2.0,8.0,18.0,40.0,
1.0,2.0,12.0,24.0,
2.0,12.0,18.0,24.0,28.0,40.0,
18.0,24.0,40.0,28.0,2.0,28.0,1.0,40.0
);
n:constant array(1..38) of Double:=(0.148_950_410_795_16e4,
0.743_077_983_140_34e3,
-0.977_083_187_978_37e2,
0.247_424_647_056_74e1,
-0.632_813_200_160_26,
0.113_859_521_296_58e1,
-0.478_118_636_486_25,
0.852_081_234_315_44e-2,
0.937_471_473_779_32,
0.335_931_186_049_16e1,
0.338_093_556_014_54e1,
0.168_445_396_719_04,
0.738_757_452_366_95,
-0.471_287_374_361_86,
0.150_202_731_397_07,
-0.217_641_142_197_50e-2,
-0.218_107_553_247_61e-1,
-0.108_297_844_036_77,
-0.463_333_246_358_12e-1,
0.712_803_519_595_51e-4,
0.110_328_317_899_99e-3,
0.189_552_483_879_02e-3,
0.308_915_411_605_37e-2,
0.135_555_045_549_49e-2,
0.286_402_374_774_56e-6,
-0.107_798_573_575_12e-4,
-0.764_627_124_548_14e-4,
0.140_523_928_183_16e-4,
-0.310_838_143_314_34e-4,
-0.103_027_382_121_03e-5,
0.282_172_816_350_40e-6,
0.127_049_022_719_45e-5,
0.738_033_534_682_92e-7,
-0.110_301_392_389_09e-7,
-0.814_563_652_078_33e-13,
-0.251_805_456_829_62e-10,
-0.175_652_339_694_07e-17,
0.869_341_563_441_63e-14
);
end Table21;
package Table22 is
I:constant array(1..23) of Double:=(-7.0,-7.0,-6.0,-6.0,-5.0,-5.0,
-2.0,-2.0,-1.0,-1.0,0.0,0.0,
1.0,1.0,2.0,6.0,6.0,6.0,6.0,
6.0,6.0,6.0,6.0
);
J:constant array(1..23) of Double:=(0.0,4.0,0.0,2.0,0.0,2.0,0.0,1.0,
0.0,2.0,0.0,1.0,4.0,8.0,4.0,0.0,
1.0,4.0,10.0,12.0,16.0,20.0,22.0
);
n:constant array(1..23) of Double:=(-0.323_683_985_552_42e13,
0.732_633_509_021_81e13,
0.358_250_899_454_47e12,
-0.583_401_318_515_90e12,
-0.107_830_682_174_70e11,
0.208_255_445_631_71e11,
0.610_747_835_645_16e6,
0.859_777_225_355_80e6,
-0.257_457_236_041_70e5,
0.310_810_884_227_14e5,
0.120_823_158_659_36e4,
0.482_197_551_092_55e3,
0.379_660_012_724_86e1,
-0.108_429_848_800_77e2,
-0.453_641_726_766_60e-1,
0.145_591_156_586_98e-12,
0.112_615_974_072_30e-11,
-0.178_049_822_406_86e-10,
0.123_245_796_908_32e-6,
-0.116_069_211_309_84e-5,
0.278_463_670_885_54e-4,
-0.592_700_384_741_76e-3,
0.129_185_829_918_78e-2
);
end Table22;
package Table25 is
I:constant array(1..46) of Double:=(1..6=>-1.5,
7..9=>-1.25,
10..15=>-1.0,
16..17=>-0.75,
18..21=>-0.5,
22..25=>-0.25,
26..29=>0.25,
30..36=>0.5,
37..40=>0.75,
41..42=>1.0,
43..44=>1.25,
45..46=>1.5
);
J:constant array(1..46) of Double:=(-24.0,-23.0,-19.0,-13.0,-11.0,-10.0,
-19.0,-15.0,-6.0,-26.0,-21.0,-17.0,
-16.0,-9.0,-8.0,-15.0,-14.0,-26.0,-13.0,
-9.0,-7.0,-27.0,-25.0,-11.0,-6.0,
1.0,4.0,8.0,11.0,0.0,1.0,5.0,6.0,
10.0,14.0,16.0,0.0,4.0,9.0,17.0,7.0,
18.0,3.0,15.0,5.0,18.0
);
n:constant array(1..46) of Double:=(-0.392_359_838_619_84e6,
0.515_265_738_272_70e6,
0.404_824_431_610_48e5,
-0.321_937_909_239_02e3,
0.969_614_242_186_94e2,
-0.228_678_463_717_73e2,
-0.449_429_141_243_57e6,
-0.501_183_360_201_66e4,
0.356_844_635_600_15,
0.442_353_358_481_90e5,
-0.136_733_888_117_08e5,
0.421_632_602_078_64e6,
0.225_169_258_374_75e5,
0.474_421_448_656_46e3,
-0.149_311_307_976_47e3,
-0.197_811_263_204_52e6,
-0.235_543_994_707_60e5,
-0.190_706_163_020_76e5,
0.553_756_698_831_64e5,
0.382_936_914_373_63e4,
-0.603_918_605_805_67e3,
0.193_631_026_203_31e4,
0.426_606_436_986_10e4,
-0.597_806_388_727_18e4,
-0.704_014_639_268_62e3,
0.338_367_841_075_53e3,
0.208_627_866_351_87e2,
0.338_341_726_561_96e-1,
-0.431_244_284_148_93e-4,
0.166_537_913_564_12e3,
-0.139_862_920_558_98e3,
-0.788_495_479_998_72,
0.721_324_117_538_72e-1,
-0.597_548_393_982_83e-2,
-0.121_413_589_539_04e-4,
0.232_270_967_338_71e-6,
-0.105_384_635_661_94e2,
0.207_189_254_965_02e1,
-0.721_931_552_604_27e-1,
0.207_498_870_811_20e-6,
-0.183_406_579_113_79e-1,
0.290_362_723_486_96e-6,
0.210_375_278_936_19,
0.256_812_397_299_99e-3,
-0.127_990_029_337_81e-1,
-0.821_981_026_520_18e-5
);
end Table25;
package Table26 is
I:constant array(1..44) of Double:=(1..2=>-6.0,
3..4=>-5.0,
5..7=>-4.0,
8..11=>-3.0,
12..15=>-2.0,
16..20=>-1.0,
21..27=>0.0,
28..33=>1.0,
34..36=>2.0,
37..39=>3.0,
40..41=>4.0,
42..44=>5.0
);
J:constant array(1..44) of Double:=(0.0,11.0,0.0,11.0,0.0,1.0,11.0,
0.0,1.0,11.0,12.0,0.0,1.0,6.0,
10.0,0.0,1.0,5.0,8.0,9.0,0.0,1.0,
2.0,4.0,5.0,6.0,9.0,0.0,1.0,2.0,3.0,
7.0,8.0,0.0,1.0,5.0,0.0,1.0,3.0,0.0,
1.0,0.0,1.0,2.0
);
n:constant array(1..44) of Double:=(0.316_876_650_834_97e6,
0.208_641_758_818_58e2,
-0.398_593_998_035_99e6,
-0.218_160_585_188_77e2,
0.223_697_851_942_42e6,
-0.278_417_034_458_17e4,
0.992_074_360_714_80e1,
-0.751_975_122_991_57e5,
0.297_086_059_511_58e4,
-0.344_068_785_485_26e1,
0.388_155_642_491_15,
0.175_112_950_857_50e5,
-0.142_371_128_544_49e4,
0.109_438_033_641_67e1,
0.899_716_193_084_95,
-0.337_597_400_989_58e4,
0.471_628_858_183_55e3,
-0.191_882_419_936_79e1,
0.410_785_804_921_96,
-0.334_653_781_720_97,
0.138_700_347_775_05e4,
-0.406_633_261_958_38e3,
0.417_273_471_596_10e2,
0.219_325_494_345_32e1,
-0.103_200_500_090_77e1,
0.358_829_435_167_03,
0.525_114_537_260_66e-2,
0.128_389_164_507_05e2,
-0.286_424_372_193_81e1,
0.569_126_836_648_55,
-0.999_629_545_849_31e-1,
-0.326_320_377_784_59e-2,
0.233_209_225_767_23e-3,
-0.153_348_098_574_50,
0.290_722_882_399_02e-1,
0.375_347_027_411_67e-3,
0.172_966_917_024_11e-2,
-0.385_560_508_445_04e-3,
-0.350_177_122_926_08e-4,
-0.145_663_936_314_92e-4,
0.564_208_572_672_69e-5,
0.412_861_500_746_05e-7,
-0.206_846_711_188_24e-7,
0.164_093_936_747_25e-8
);
end Table26;
package Table27 is
I:constant array(1..30) of Double:=(1..2=>-2.0,
3=>-1.0,
4..7=>0.0,
8..11=>1.0,
12..14=>2.0,
15..17=>3.0,
18..20=>4.0,
21..23=>5.0,
24..25=>6.0,
26..30=>7.0
);
J:constant array(1..30) of Double:=(0.0,1.0,0.0,0.0,1.0,2.0,3.0,
0.0,1.0,3.0,4.0,0.0,1.0,2.0,0.0,
1.0,5.0,0.0,1.0,4.0,0.0,1.0,2.0,
0.0,1.0,0.0,1.0,3.0,4.0,5.0
);
n:constant array(1..30) of Double:=(0.909_685_010_053_65e3,
0.240_456_670_884_20e4,
-0.591_623_263_871_30e3,
0.541_454_041_280_74e3,
-0.270_983_084_111_92e3,
0.979_765_250_979_26e3,
-0.469_667_729_594_35e3,
0.143_992_746_047_23e2,
-0.191_042_042_304_29e2,
0.532_991_671_11971e1,
-0.212_529_753_759_34e2,
-0.311_473_344_137_60,
0.603_348_408_946_23,
-0.427_648_397_025_09e-1,
0.581_855_972_552_59e-2,
-0.145_970_082_847_53e-1,
0.566_311_756_310_27e-2,
-0.761_558_645_845_77e-4,
0.224_403_429_193_32e-3,
-0.125_610_950_134_13e-4,
0.633_231_326_609_34e-6,
-0.205_419_896_753_75e-5,
0.364_053_703_900_82e-7,
-0.297_598_977_892_15e-8,
0.101_366_185_297_63e-7,
0.599_257_196_923_51e-11,
-0.206_778_701_051_64e-10,
-0.208_742_781_818_86e-10,
0.101_621_668_250_89e-9,
-0.164_298_282_813_47e-9
);
end Table27;
function bValidEq15(p,T:Double) return Boolean is
begin
if T in 273.15..623.15 and (p>0.0 and p<=Ps(T)) then
return True;
elsif T in 623.15..863.15 and (p>0.0 and p<=B23.P(T)) then
return True;
elsif T in 863.15..1073.15 and (p>0.0 and p<=100.0) then
return True;
end if;
return false;
end bValidEq15;
function bValidEq18(p,T:Double) return Boolean is
begin
if T in 273.15..584.1495 and p in 0.0..10.0 then
return True;
end if;
return False;
end bValidEq18;
-- 考虑单相稳定区与亚稳定区
function g(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
sum:=sum+log(pi);
declare
use Table10;
begin
if bValidEq15(p,T) then
for L in 1..9 loop
sum:=sum+n(L)*tau**J(L);
end loop;
elsif bValidEq18(p,T) then
for L in 1..9 loop
sum:=sum+n2(L)*tau**J(L);
end loop;
end if;
end;
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*(pi**I(L))*((tau-0.5)**J(L));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*(pi**I(L))*((tau-0.5)**J(L));
end loop;
end;
end if;
return sum*R*T;
end g;
function gamma(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
sum:=sum+log(pi);
declare
use Table10;
begin
if bValidEq15(p,T) then
for L in 1..9 loop
sum:=sum+n(L)*tau**J(L);
end loop;
elsif bValidEq18(p,T) then
for L in 1..9 loop
sum:=sum+n2(L)*tau**J(L);
end loop;
end if;
end;
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*(pi**I(L))*((tau-0.5)**J(L));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*(pi**I(L))*((tau-0.5)**J(L));
end loop;
end;
end if;
return sum;
end gamma;
-- 稳定单相区与亚稳定区都包含
function Go(p,T:Double) return Double is
use Table10;
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=log(pi);
begin
if bValidEq15(p,T) then
for L in 1..9 loop
sum:=sum+n(L)*tau**J(L);
end loop;
elsif bValidEq18(p,T) then
for L in 1..9 loop
sum:=sum+n2(L)*tau**J(L);
end loop;
end if;
return sum;
end Go;
function Gop(p,T:Double) return Double is
pi:Double:=p/1.0;
begin
return 1.0/pi+0.0;
end Gop;
function Gopp(p,T:Double) return Double is
pi:Double:=p/1.0;
begin
return -1.0/(pi**2)+0.0;
end Gopp;
function Got(p,T:Double) return Double is
use Table10;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
for L in 1..9 loop
sum:=sum+n(L)*J(L)*(tau**(J(L)-1.0));
end loop;
elsif bValidEq18(p,T) then
for L in 1..9 loop
sum:=sum+n2(L)*J(L)*(tau**(J(L)-1.0));
end loop;
end if;
return sum;
end Got;
function Gott(p,T:Double) return Double is
use Table10;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
for L in 1..9 loop
sum:=sum+n(L)*J(L)*(J(L)-1.0)*(tau**(J(L)-2.0));
end loop;
elsif bValidEq18(p,T) then
for L in 1..9 loop
sum:=sum+n2(L)*J(L)*(J(L)-1.0)*(tau**(J(L)-2.0));
end loop;
end if;
return sum;
end Gott;
function Gopt(p,T:Double) return Double is
begin
return 0.0;
end Gopt;
function Gr(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*(pi**I(L))*((tau-0.5)**J(L));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*(pi**I(L))*((tau-0.5)**J(L));
end loop;
end;
end if;
return sum;
end Gr;
function Grp(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*I(L)*(pi**(I(L)-1.0))*((tau-0.5)**J(L));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*I(L)*(pi**(I(L)-1.0))*((tau-0.5)**J(L));
end loop;
end;
end if;
return sum;
end Grp;
function Grt(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*(pi**I(L))*J(L)*((tau-0.5)**(J(L)-1.0));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*(pi**I(L))*J(L)*((tau-0.5)**(J(L)-1.0));
end loop;
end;
end if;
return sum;
end Grt;
function Grpp(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*I(L)*(I(L)-1.0)*(pi**(I(L)-2.0))*((tau-0.5)**J(L));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*I(L)*(I(L)-1.0)*(pi**(I(L)-2.0))*((tau-0.5)**J(L));
end loop;
end;
end if;
return sum;
end Grpp;
function Grtt(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*(pi**I(L))*J(L)*(J(L)-1.0)*((tau-0.5)**(J(L)-2.0));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*(pi**I(L))*J(L)*(J(L)-1.0)*((tau-0.5)**(J(L)-2.0));
end loop;
end;
end if;
return sum;
end Grtt;
function Grpt(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
sum:Double:=0.0;
begin
if bValidEq15(p,T) then
declare
use Table11;
begin
for L in 1..43 loop
sum:=sum+n(L)*I(L)*(pi**(I(L)-1.0))*J(L)*((tau-0.5)**(J(L)-1.0));
end loop;
end;
elsif bValidEq18(p,T) then
declare
use Table16;
begin
for L in 1..13 loop
sum:=sum+n(L)*I(L)*(pi**(I(L)-1.0))*J(L)*((tau-0.5)**(J(L)-1.0));
end loop;
end;
end if;
return sum;
end Grpt;
-- 比体积
function v(p,T:Double) return Double is
pi:Double:=p/1.0;
begin
return pi*(Gop(p,T)+Grp(p,T))*R*1.0e3*T/(p*1.0e6);
end v;
-- 比内能
function u(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
begin
return R*T*(tau*(Got(p,T)+Grt(p,T))-pi*(Gop(p,T)+Grp(p,T)));
end u;
-- 比熵
function s(p,T:Double) return Double is
tau:Double:=540.0/T;
begin
return R*(tau*(Got(p,T)+Grt(p,T))-(Go(p,T)+Gr(p,T)));
end s;
-- 比焓
function h(p,T:Double) return Double is
tau:Double:=540.0/T;
begin
return R*T*(tau*(Got(p,T)+Grt(p,T)));
end h;
-- 比等压热容
function cp(p,T:Double) return Double is
tau:Double:=540.0/T;
begin
return R*(-(tau**2)*(Gott(p,T)+Grtt(p,T)));
end cp;
-- 比等体积热容
function cv(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
begin
return R*(-(tau**2)*(Gott(p,T)+Grtt(p,T))-(1.0+pi*Grp(p,T)-tau*pi*Grpt(p,T))**2/(1.0-pi**2*Grpp(p,T)));
end cv;
-- 声速
function w(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=540.0/T;
begin
return (R*1000.0*T*(1.0+2.0*pi*Grp(p,T)+pi**2*Grp(p,T)**2)/((1.0-pi**2*Grpp(p,T))+(1.0+pi*Grp(p,T)-tau*pi*Grpt(p,T))**2/(tau**2*(Gott(p,T)+Grtt(p,T)))))**0.5;
end w;
--
-- 导出方程
--
function P_B2bc(h:Double) return Double is
use Table19;
begin
return n(1)+n(2)*h+n(3)*(h**2.0);
end P_B2bc;
function h_B2bc(p:Double) return Double is
use Table19;
begin
return n(4)+((p-n(5))/n(3))**0.5;
end h_B2bc;
function T2a_ph(p,h:Double) return Double is
use Table20;
pi:Double:=p/1.0;
eta:Double:=h/2000.0;
sum:Double:=0.0;
begin
for L in 1..34 loop
sum:=sum+n(L)*(pi**I(L))*((eta-2.1)**J(L));
end loop;
return sum;
end T2a_ph;
function T2b_ph(p,h:Double) return Double is
use Table21;
pi:Double:=p/1.0;
eta:Double:=h/2000.0;
sum:Double:=0.0;
begin
for L in 1..38 loop
sum:=sum+n(L)*(pi-2.0)**I(L)*(eta-2.6)**J(L);
end loop;
return sum;
end T2b_ph;
function T2c_ph(p,h:Double) return Double is
use Table22;
pi:Double:=p/1.0;
eta:Double:=h/2000.0;
sum:Double:=0.0;
begin
for L in 1..23 loop
sum:=sum+n(L)*(pi+25.0)**I(L)*(eta-1.8)**J(L);
end loop;
return sum;
end T2c_ph;
function T2a_ps(p,s:Double) return Double is
use Table25;
pi:Double:=p/1.0;
sigma:Double:=s/2.0;
sum:Double:=0.0;
begin
for L in 1..46 loop
sum:=sum+n(L)*(pi**I(L))*((sigma-2.0)**J(L));
end loop;
return sum;
end T2a_ps;
function T2b_ps(p,s:Double) return Double is
use Table26;
pi:Double:=p/1.0;
sigma:Double:=s/0.7853;
sum:Double:=0.0;
begin
for L in 1..44 loop
sum:=sum+n(L)*(pi**I(L))*((10.0-sigma)**J(L));
end loop;
return sum;
end T2b_ps;
function T2c_ps(p,s:Double) return Double is
use Table27;
pi:Double:=p/1.0;
sigma:Double:=s/2.9251;
sum:Double:=0.0;
begin
for L in 1..30 loop
sum:=sum+n(L)*(pi**I(L))*((2.0-sigma)**J(L));
end loop;
return sum;
end T2c_ps;
function T_ph(p,h:Double) return Double is
begin
if p>0.0 and p<4.0 then
return T2a_ph(p,h);
elsif (p>=4.0 and p<=100.0) and p<P_B2bc(h) then
return T2b_ph(p,h);
elsif (p>=4.0 and p<=100.0) and p>=P_B2bc(h) then
return T2c_ph(p,h);
end if;
return 0.0;
end T_ph;
function T_ps(p,s:Double) return Double is
begin
if p>0.0 and p<4.0 and s>=5.85 then
return T2a_ps(p,s);
elsif p>=4.0 and p<=100.0 and s>=5.85 then
return T2b_ps(p,s);
elsif p>=4.0 and p<=100.0 and s<5.85 then
return T2c_ps(p,s);
end if;
return 0.0;
end T_ps;
end IF97.Region2;
------------------------------------------------------------------
-- File: IF97-Region3.ads
package IF97.Region3 is
-- 霍尓姆兹自由能
function f(rho,T:Double) return Double;
function p(rho,T:Double) return Double; -- return MPa
function u(rho,T:Double) return Double; -- return kJ*kg-1
function s(rho,T:Double) return Double; -- return kJ*kg-1*K-1
function h(rho,T:Double) return Double;
function cv(rho,T:Double) return Double;
function cp(rho,T:Double) return Double;
function w(rho,T:Double) return Double;
-- 3区根据压力温度压力求密度,unit参数表示精度
function Rho3(p,T:Double;unit:Double:=1.0e-10) return Double;
end IF97.Region3;
-----------------------------------------------------------------------
-- File: IF97-Region3.adb
with IF97.Region4;
package body IF97.Region3 is
package Table30 is
I:constant array(2..40) of Double:=(2..8=>0.0,
9..12=>1.0,
13..18=>2.0,
19..23=>3.0,
24..27=>4.0,
28..30=>5.0,
31..33=>6.0,
34=>7.0,
35=>8.0,
36..37=>9.0,
38..39=>10.0,
40=>11.0
);
J:constant array(2..40) of Double:=(0.0,1.0,2.0,7.0,10.0,12.0,
23.0,2.0,6.0,15.0,17.0,0.0,
2.0,6.0,7.0,22.0,26.0,0.0,
2.0,4.0,16.0,26.0,0.0,2.0,
4.0,26.0,1.0,3.0,26.0,0.0,
2.0,26.0,2.0,26.0,2.0,26.0,0.0,
1.0,26.0
);
n:constant array(1..40) of Double:=(0.106_580_700_285_13e1,
-0.157_328_452_902_39e2,
0.209_443_969_743_07e2,
-0.768_677_078_787_16e1,
0.261_859_477_879_54e1,
-0.280_807_811_486_20e1,
0.120_533_696_965_17e1,
-0.845_668_128_125_02e-2,
-0.126_543_154_777_14e1,
-0.115_244_078_066_81e1,
0.885_210_439_843_18,
-0.642_077_651_816_07,
0.384_934_601_866_71,
-0.852_147_088_242_06,
0.489_722_815_418_77e1,
-0.305_026_172_569_65e1,
0.394_205_368_791_54e-1,
0.125_584_084_243_08,
-0.279_993_296_987_10,
0.138_997_995_694_60e1,
-0.201_899_150_235_70e1,
-0.821_476_371_739_63e-2,
-0.475_960_357_349_23,
0.439_840_744_735_00e-1,
-0.444_764_354_287_39,
0.905_720_707_197_33,
0.705_224_500_879_67,
0.107_705_126_263_32,
-0.329_136_232_589_54,
-0.508_710_620_411_58,
-0.221_754_008_730_96e-1,
0.942_607_516_650_92e-1,
0.164_362_784_479_61,
-0.135_033_722_413_48e-1,
-0.148_343_453_524_72e-1,
0.579_229_536_280_84e-3,
0.323_089_047_037_11e-2,
0.809_648_029_962_15e-4,
-0.165_576_797_950_37e-3,
-0.449_238_990_618_15e-4
);
end Table30;
function f(rho,T:Double) return Double is
use Table30;
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
sum:Double:=n(1)*log(delta0);
begin
for L in 2..40 loop
sum:=sum+n(L)*(delta0**I(L))*(tau**J(L));
end loop;
return sum*R*T;
end f;
function phi(rho,T:Double) return Double is
use Table30;
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
sum:Double:=n(1)*log(delta0);
begin
for L in 2..40 loop
sum:=sum+n(L)*(delta0**I(L))*(tau**J(L));
end loop;
return sum;
end phi;
function Phi_d(rho,T:Double) return Double is
use Table30;
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
sum:Double:=n(1)/delta0;
begin
for L in 2..40 loop
sum:=sum+n(L)*I(L)*(delta0**(I(L)-1.0))*(tau**J(L));
end loop;
return sum;
end phi_d;
function Phi_dd(rho,T:Double) return Double is
use Table30;
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
sum:Double:=-n(1)/(delta0**2);
begin
for L in 2..40 loop
sum:=sum+n(L)*I(L)*(I(L)-1.0)*(delta0**(I(L)-2.0))*(tau**J(L));
end loop;
return sum;
end phi_dd;
function Phi_t(rho,T:Double) return Double is
use Table30;
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
sum:Double:=0.0;
begin
for L in 2..40 loop
sum:=sum+n(L)*J(L)*(delta0**I(L))*(tau**(J(L)-1.0));
end loop;
return sum;
end phi_t;
function Phi_tt(rho,T:Double) return Double is
use Table30;
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
sum:Double:=0.0;
begin
for L in 2..40 loop
sum:=sum+n(L)*J(L)*(J(L)-1.0)*(delta0**I(L))*(tau**(J(L)-2.0));
end loop;
return sum;
end phi_tt;
function Phi_dt(rho,T:Double) return Double is
use Table30;
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
sum:Double:=0.0;
begin
for L in 2..40 loop
sum:=sum+n(L)*I(L)*J(L)*(delta0**(I(L)-1.0))*(tau**(J(L)-1.0));
end loop;
return sum;
end phi_dt;
function p(rho,T:Double) return Double is
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
begin
return rho*R*1000.0*T*(delta0*Phi_d(rho,T))/1.0e6;
end p;
function u(rho,T:Double) return Double is
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
begin
return R*T*(tau*Phi_t(rho,T));
end u;
function s(rho,T:Double) return Double is
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
begin
return R*(tau*phi_t(rho,T)-phi(rho,T));
end s;
function h(rho,T:Double) return Double is
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
begin
return R*T*(tau*phi_t(rho,T)+delta0*phi_d(rho,T));
end h;
function cv(rho,T:Double) return Double is
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
begin
return R*(-tau**2*phi_tt(rho,T));
end cv;
function cp(rho,T:Double) return Double is
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
begin
return R*(-tau**2*phi_tt(rho,T)+((delta0*Phi_d(rho,T)-delta0*tau*phi_dt(rho,T))**2.0)/(2.0*delta0*Phi_d(rho,T)+delta0**2*Phi_dd(rho,T)));
end cp;
function w(rho,T:Double) return Double is
delta0:Double:=rho/322.0;
tau:Double:=647.096/T;
begin
return (R*1000.0*T*(2.0*delta0*phi_d(rho,T)+delta0**2*Phi_dd(rho,T)-(delta0*Phi_d(rho,T)-delta0*tau*phi_dt(rho,T))**2/(tau**2*phi_tt(rho,T))))**0.5;
end w;
function Rho3(p,T:Double;unit:Double:=1.0e-10) return Double is
use IF97.Region4;
Rho_old,Rho_new,Rho_p,dlt:double;
begin
if T in 623.15..647.096 and p<Ps(T) then
Rho_old:=100.0;
else
Rho_old:=600.0;
end if;
loop
dlt:=Rho_old/322.0;
Rho_p:=1000.0*R*T*(2.0*Rho_old*Phi_d(Rho_old,T)+Rho_old**2*phi_dd(Rho_old,T)/322.0);
Rho_new:=Rho_old+(1.0e6*p-1000.0*R*T*Rho_old**2*phi_d(Rho_old,T)/322.0)/Rho_p;
if abs(Rho_new-Rho_old)<unit then
return Rho_new;
end if;
Rho_old:=Rho_new;
end loop;
end Rho3;
end IF97.Region3;
------------------------------------------------------------------------------------------------------------
-- File: IF97-Region5.ads
package IF97.Region5 is
function g(p,T:Double) return Double;
function v(p,T:Double) return Double;
function u(p,T:Double) return Double;
function s(p,T:Double) return Double;
function h(p,T:Double) return Double;
function cp(p,T:Double) return Double;
function cv(p,T:Double) return Double;
function w(p,T:Double) return Double;
end IF97.Region5;
-------------------------------------------------------------------------------------------------------------
-- File: IF97-Region5.adb
package body IF97.Region5 is
package Table37 is
J:constant array(1..6) of Double:=(0.0,1.0,-3.0,-2.0,-1.0,2.0);
n:constant array(1..6) of Double:=(-0.131_799_836_742_01e2,
0.685_408_416_344_34e1,
-0.248_051_489_334_66e-1,
0.369_015_349_803_33,
-0.311_613_182_139_25e1,
-0.329_616_265_389_17
);
end Table37;
package Table38 is
I:constant array(1..6) of Double:=(1.0,1.0,1.0,2.0,2.0,3.0);
J:constant array(1..6) of Double:=(1.0,2.0,3.0,3.0,9.0,7.0);
n:constant array(1..6) of Double:=(0.157_364_048_552_59e-2,
0.901_537_616_739_44e-3,
-0.502_700_776_776_48e-2,
0.224_400_374_094_85e-5,
-0.411_632_754_534_71e-5,
0.379_194_548_229_55e-7
);
end Table38;
function Go(p,T:Double) return Double is
use Table37;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=Double_Math.log(pi);
begin
for L in 1..6 loop
sum:=sum+n(L)*tau**J(L);
end loop;
return sum;
end Go;
function Gop(p,T:Double) return Double is
begin
return 1.0/(p/1.0);
end Gop;
function Gopp(p,T:Double) return Double is
begin
return -1.0/((p/1.0)**2);
end Gopp;
function Gopt(p,T:Double) return Double is
begin
return 0.0;
end Gopt;
function Got(p,T:Double) return Double is
use Table37;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*J(L)*(tau**(J(L)-1.0));
end loop;
return sum;
end Got;
function Gott(p,T:Double) return Double is
use Table37;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*J(L)*(J(L)-1.0)*(tau**(J(L)-2.0));
end loop;
return sum;
end Gott;
function Gr(p,T:Double) return Double is
use Table38;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*(pi**I(L))*(tau**J(L));
end loop;
return sum;
end Gr;
function Grp(p,T:Double) return Double is
use Table38;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*I(L)*(pi**(I(L)-1.0))*(tau**J(L));
end loop;
return sum;
end Grp;
function Grpp(p,T:Double) return Double is
use Table38;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*I(L)*(I(L)-1.0)*(pi**(I(L)-2.0))*(tau**J(L));
end loop;
return sum;
end Grpp;
function Grt(p,T:Double) return Double is
use Table38;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*J(L)*(pi**I(L))*(tau**(J(L)-1.0));
end loop;
return sum;
end Grt;
function Grtt(p,T:Double) return Double is
use Table38;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*J(L)*(J(L)-1.0)*(pi**I(L))*(tau**(J(L)-2.0));
end loop;
return sum;
end Grtt;
function Grpt(p,T:Double) return Double is
use Table38;
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
sum:Double:=0.0;
begin
for L in 1..6 loop
sum:=sum+n(L)*I(L)*J(L)*(pi**(I(L)-1.0))*(tau**(J(L)-1.0));
end loop;
return sum;
end Grpt;
function g(p,T:Double) return Double is
begin
return R*T*(Go(p,T)+Gr(p,T));
end g;
function v(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
begin
return R*1000.0*T*pi*(gop(p,T)+grp(p,T))/(p*1.0e6);
end v;
function u(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
begin
return R*T*(tau*(got(p,T)+grt(p,T))-pi*(gop(p,T)+grp(p,T)));
end u;
function s(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
begin
return R*(tau*(got(p,T)+grt(p,T))-(go(p,T)+gr(p,T)));
end s;
function h(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
begin
return R*T*(tau*(got(p,T)+grt(P,T)));
end h;
function cp(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
begin
return R*(-(tau**2)*(gott(p,T)+grtt(p,T)));
end cp;
function cv(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
begin
return R*(-tau**2*(gott(p,T)+grtt(p,T))-(1.0+pi*grp(p,T)-tau*pi*grpt(p,T))**2/(1.0-pi**2*grpp(p,T)));
end cv;
function w(p,T:Double) return Double is
pi:Double:=p/1.0;
tau:Double:=1000.0/T;
begin
return (R*1000.0*T*(1.0+2.0*pi*Grp(p,T)+(pi**2)*(grp(p,T)**2))/((1.0-pi**2*grpp(p,T))+(1.0+pi*grp(p,T)-tau*pi*grpt(p,T))**2/(tau**2*(gott(p,T)+grtt(p,T)))))**0.5;
end w;
end IF97.Region5;
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------
-- File: IF97-Interfaces.ads,IF97的对外接口
package IF97.Interfaces is
-- 压力:MPa;温度:K;比容:m3/(kg*K);密度:kg/m3;
-- 已知压力温度求比容
function v(p,T:Double) return Double;
-- 已知压力温度求密度
function r(p,T:Double) return Double;
-- 已知压力温度求内能
function u(p,T:Double) return Double;
-- 已知压力温度求熵
function s(p,T:Double) return Double;
-- 已知压力温度求焓
function h(p,T:Double) return Double;
-- 已知压力温度求等压热容
function cp(p,T:Double) return Double;
-- 已知压力温度求等容热容
function cv(p,T:Double) return Double;
-- 已知压力温度求声速
function w(p,T:Double) return Double;
-- 根据温度求饱和压力(4-2,2-3),即饱和蒸汽
function Pss(T:Double) return Double;
-- 根据压力求饱和温度(4-2,2-3),即饱和蒸汽
function Tss(P:Double) return Double;
-- 根据温度求饱和水压力(4-2,1-3)
function Psw(T:Double) return Double;
-- 根据压力求饱和水温度
function Tsw(P:Double) return Double;
end IF97.Interfaces;
---------------------------------------------------
-- File: IF97-Interfaces.adb,IF97的对外接口
with IF97.B23;
with IF97.Region1;
with IF97.Region2;
with IF97.Region3;
with IF97.Region4;
with IF97.Region5;
package body IF97.Interfaces is
function Region(p,T:Double) return Integer is
begin
-- 这里重叠了1区与4区<=623.15的部分
if T in 273.15..623.15 and p>=IF97.Region4.Ps(T) and p<=100.0 then
return 1;
elsif T in 273.15..623.15 and p>=0.0 and p<IF97.Region4.Ps(T) then
return 2;
elsif T in 623.15..863.15 and p>=0.0 and p<B23.P(T) then
return 2;
elsif T in 863.15..1073.15 and p>=0.0 and p<=100.0 then
return 2;
elsif T in 623.15..B23.T(P) and p in B23.P(T)..100.0 then
return 3;
elsif T in 273.15..647.096 and p in 0.0..22.064 then
return 4;
elsif T in 1073.15..2273.15 and p in 0.0..50.0 then
return 5;
end if;
return 0;
end Region;
function v(p,T:Double) return Double is
R:Integer:=Region(p,T);
begin
if R=1 then
return Region1.v(p,T);
elsif R=2 then
return Region2.v(p,T);
elsif R=3 then
return 1.0/Region3.Rho3(p,T);
elsif R=4 then
if T in 273.15..623.15 then
return Region1.v(p,T);
else
return 1.0/Region3.Rho3(p,T);
end if;
elsif R=5 then
return Region5.v(p,T);
end if;
return -1.0;
end v;
function r(p,T:Double) return Double is
begin
return 1.0/v(p,T);
end r;
function u(p,T:Double) return Double is
R:Integer:=Region(p,T);
begin
if R=1 then
return Region1.u(p,T);
elsif R=2 then
return Region2.u(p,T);
elsif R=3 then
return Region3.u(Region3.Rho3(p,T),T);
elsif R=4 then
if T in 273.15..623.15 then
return Region1.u(p,T);
else
return Region3.u(Region3.Rho3(p,T),T);
end if;
elsif R=5 then
return Region5.u(p,T);
end if;
return -1.0;
end u;
function s(p,T:Double) return Double is
R:Integer:=Region(p,T);
begin
if R=1 then
return Region1.s(p,T);
elsif R=2 then
return Region2.s(p,T);
elsif R=3 then
return Region3.s(Region3.Rho3(p,T),T);
elsif R=4 then
if T in 273.15..623.15 then
return Region1.s(p,T);
else
return Region3.s(Region3.Rho3(p,T),T);
end if;
elsif R=5 then
return Region5.s(p,T);
end if;
return -1.0;
end s;
function h(p,T:Double) return Double is
R:Integer:=Region(p,T);
begin
if R=1 then
return Region1.h(p,T);
elsif R=2 then
return Region2.h(p,T);
elsif R=3 then
return Region3.h(Region3.Rho3(p,T),T);
elsif R=4 then
if T in 273.15..623.15 then
return Region1.h(p,T);
else
return Region3.h(Region3.Rho3(p,T),T);
end if;
elsif R=5 then
return Region5.h(p,T);
end if;
return -1.0;
end h;
function cp(p,T:Double) return Double is
R:Integer:=Region(p,T);
begin
if R=1 then
return Region1.cp(p,T);
elsif R=2 then
return Region2.cp(p,T);
elsif R=3 then
return Region3.cp(Region3.Rho3(p,T),T);
elsif R=4 then
if T in 273.15..623.15 then
return Region1.cp(p,T);
else
return Region3.cp(Region3.Rho3(p,T),T);
end if;
elsif R=5 then
return Region5.cp(p,T);
end if;
return -1.0;
end cp;
function cv(p,T:Double) return Double is
R:Integer:=Region(p,T);
begin
if R=1 then
return Region1.cv(p,T);
elsif R=2 then
return Region2.cv(p,T);
elsif R=3 then
return Region3.cv(Region3.Rho3(p,T),T);
elsif R=4 then
if T in 273.15..623.15 then
return Region1.cv(p,T);
else
return Region3.cv(Region3.Rho3(p,T),T);
end if;
elsif R=5 then
return Region5.cv(p,T);
end if;
return -1.0;
end cv;
function w(p,T:Double) return Double is
R:Integer:=Region(p,T);
begin
if R=1 then
return Region1.w(p,T);
elsif R=2 then
return Region2.w(p,T);
elsif R=3 then
return Region3.w(Region3.Rho3(p,T),T);
elsif R=4 then
if T in 273.15..623.15 then
return Region1.w(p,T);
else
return Region3.w(Region3.Rho3(p,T),T);
end if;
elsif R=5 then
return Region5.w(p,T);
end if;
return -1.0;
end w;
-- 根据温度求饱和压力
function Pss(T:Double) return Double is
begin
if T >= 623.15 then
return B23.P(T);
else
return Region4.Ps(T);
end if;
end Pss;
-- 根据压力求饱和温度
function Tss(P:Double) return Double is
begin
if P>=16.5291642526 then
return B23.T(P);
else
return Region4.Ts(P);
end if;
end Tss;
function Psw(T:Double) return Double is
begin
if T > 623.15 then
return Region3.P(Region3.Rho3(B23.P(T),623.15),623.15);
else
return Region4.Ps(T);
end if;
end Psw;
function Tsw(P:Double) return Double is
begin
if p >= 16.5291642526 then
declare
Ts:double:=B23.T(P);
begin
return (if Ts<623.15 then Ts else 623.15);
end;
else
return Region4.Ts(P);
end if;
end Tsw;
end IF97.Interfaces;
1032
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
