water.pas

delphi water math Programning with delphi.

{
     ****************************************************************
     *                                                              *
     * Pascal functions for calculating                             *
     * properties of water and steam                                *
     *                                                              *
     * Authored by Gerard Koops ([e-mail])                      *
     * with some revisions by Bernhard Spang, The Mining Company    *
     * URL:        http://chemengineer.miningco.com                 *
     * email:      [e-mail]                  *
     *                                                              *
     * May be used and modified for free. Redistribution is also    * 
     * allowed if full credit is given to the authors and a link to *
     * http://chemengineer.miningco.com is included.                *
     * Provided "as is" without warranty of any kind.               *
     *                                                              *
     * For documentation see                                        *
     * http://chemengineer.miningco.com/library/weekly/aa083198.htm *
     *                                                              *
     * Reference:                                                   *
     * Properties of Water and Steam in SI-Units,                   *
     * 2nd Revised and Updated Printing, Springer 1979, pp. 175 ff. *
     *                                                              *
     ****************************************************************
}

Unit Water;

{Functions for steam/water properties, ported from visual basic into pascal}
 {porting version 1.1 : 11-07-1998
 Writen By B. Spang}

INTERFACE

function VAPP(T:real):real;
{Saturationpressure out of Temperature
 T in K
 VAPP in bar
 When out of reach, VAPP = -1}

function BOILP(P:real):real;
{Saturationtemperature out of pressure
 0.01 < p < 220 bar
 P in bar
 BOILP in K
 When out of reach, BOILP = -1}

Function HSUB1(T, P:real):real;
{Enthalpy of water out of pressure and temperature.
T in K  273.16 K < T < 623.15 K
P in bar saturationpressure < P < 1000 bar}


Function HSUB2(T, P:real):real;
{Enthalpy of steam out of pressure and temperature.
T in K 273.16 K < T 1073.15 K
P in bar 0 < P < saturationpressure
HSUB2 in KJ/kg}

Function VSUB1(T,P:real):real;
{Specific volume of water out of temperature and pressure
T in K  273.16 K < T < 623.15 K
P in bar saturationpressure < P < 1000 bar
VUB1 in m^3/kg}

Function VSUB2(T,P:real):real;
{Specific volume of steam out of temperature and pressure
T in K 273.16 K < T 1073.15 K
P in bar 0 < P < saturationpressure
HSUB2 in m^3/kg}

Function ETAF(T, V:real):real;
{viscosity of water out of temperature and specific volume
 T upto 800 癈, P upto 1000 bar
 T in K
 v in M^3/kg
 ETAF in E-6 kg/m*s
 WATCH IT!!!!!!!!!!!!!!!!!!}

Function TCON(T, V: real):real;
{Thermal conductivity of water out of temperature and specific volume
up to 1500癈 and 3000 bar
T in K
V in m^3/kg
TCON in W/(K*m)}

Function HCAP (T,P:real;IAG:integer):real;
{Specific heatcapacity, out of temperature and pressure.
 Nummeric solution of differentials step 0.5 K
 T in K
 P in bar
 IAG = 1: fluid    2: vapour
 HCAP in kJ/(kg*K)}


function BL(TR:real):real;
{Hulpfunctie voor HSUB2 en VSUB2}

function DBL(TR:real):real;
{Hulpfunctie voor HSUB2 en VSUB2}

implementation

function VAPP(T:real):real;

const

Tc = 647.3;
Pc = 221.2;
exponent = 0.4;
var

AK: array[1..12] of real;
BK: array[1..4] of real;
Tz: array[1..12] of real;
TR,Som,Pr,U: real;
I:integer;

begin

AK[1] := -4.0596821;
AK[2] := 5.1322555;
AK[3] := -1.1842407;
AK[4] := 0.11779592;
AK[5] := -0.005157642;
AK[6] := -0.0014689537;
AK[7] := 0.00053622818;
AK[8] := 0.00012455399;
AK[9] := -0.000049154288;
AK[10] := 0.000046302565;
AK[11] := 0.000015301334;
AK[12] := -0.00002095453;
BK[1] := 2;
BK[2] := 0.95;
BK[3] := 1.45220717;
BK[4] := -0.84878953;

{afhandeling van overschreiding bereik}
TR := T / Tc;
If (TR < 0.421) or (TR > 1) then 
      begin
            VAPP := -1;
      end
else 
      begin
            U := (BK[1]* exp(0.4*ln(1/TR - BK[2]))  - BK[3])/BK[4];
            Tz[1] := 1;
            Tz[2] := U;
            For I := 3 To 12 do
                Tz[I] := 2*U*Tz[I-1] - Tz[I-2];

            Som := 0;
            For I := 1 to 12 do
                Som := Som + AK[I] * Tz[I];

            PR := Exp(Som);
            VAPP := PR * Pc;
      end;
end;

function BOILP(P:real):real;

Var

Tunt,Tob,Tm,Punt,Pob,Pm,DPunt,DPob,DPm:real;
I : integer;
begin
Tunt := 280.15;
Tob := 646.85;
Tm := 0.5 * (Tunt + Tob);
Punt := VAPP(Tunt);
Pob :=VAPP(Tob);
Pm := VAPP(Tm);
DPunt := Punt - P;
DPob := Pob - P;
DPm := Pm - P;
If (DPunt*Dpob)> 0 then
begin
Tm:=-1;
BOILP:=-1;
end

else begin

REPEAT

      If (DPm * DPunt) < 0 then
        begin
           Tob := Tm;
           Pob := Pm;
           DPob := DPm;
        end
      else
        begin
           Tunt := Tm;
           Punt := Pm;
           DPunt := DPm;
        end;
     Tm := 0.5 *(Tunt+Tob);
     BOILP := Tm;
     Pm := VAPP(Tm);
     DPm := Pm - P;
UNTIL (Tm - Tunt) < 0.00005;
     end;
   end;

Function HSUB1(T, P:real):real;
Const
Tc = 647.3;
Pc = 221.2;
Hc = 70.1204;
L = 17;
m = 19;
Var
AK:array[1..12]of real;
AG:array[1..12]of real;
AGA:array[1..11]of real;
TR,PR:real;
H1,H2,H3,H4,H5,H6,H7,H8,HR:real;
Y,Y1,Z,zwt1:real;
I,I1:integer;
begin


 AK[1] := 0.8438375405;
      AK[2] := 0.0005362162162;
      AK[3] := 1.72;
      AK[4] := 0.07342278489;
      AK[5] := 0.0497585887;
      AK[6] := 0.65371543;
      AK[7] := 0.00000115;
      AK[8] := 0.000015108;
      AK[9] := 0.14188;
      AK[10] := 7.002753165;
      AK[11] := 0.0002995284926;
      AK[12] := 0.204;
      AG[1] := 7.982692717;
      AG[2] := -0.02616571843;
      AG[3] := 0.00152241179;
      AG[4] := 0.02284279054;
      AG[5] := 242.1647003;
      AG[6] := 1.269716088E-10;
      AG[7] := 2.074838328E-07;
      AG[8] := 2.17402035E-08;
      AG[9] := 1.105710498E-09;
      AG[10] := 12.93441934;
      AG[11] := 0.00001308119072;
      AG[12] := 6.047626338E-14;
      AGA[1] := 6824.687741;
      AGA[2] := -542.2063673;
      AGA[3] := -20966.66205;
      AGA[4] := 39412.86787;
      AGA[5] := -67332.77739;
      AGA[6] := 99023.81028;
      AGA[7] := -109391.1774;
      AGA[8] := 85908.41667;
      AGA[9] := -45111.68742;
      AGA[10] := 14181.38926;
      AGA[11] := -2017.271113;

      TR := T / Tc;
      PR := P / Pc;

      Y := 1 - AK[1] * TR * TR - AK[2] / exp(6*ln(TR));
      Y1 := 6 * AK[2] / exp(7*ln(TR)) - 2 * AK[1] * TR;
      Z := Y + Sqrt(AK[3] * Y * Y - 2 * AK[4] * TR + 2 * AK[5] * PR);

      H1 := AGA[1] * TR;
      H2 := 0;
      For I := 1 To 10 do
	begin
        I1 := I + 1;
        H2 := H2 + (I - 2) * AGA[I1] * exp((I - 1)*ln(TR));
        end;

      H3 := AG[1] * (Z * (17 * (Z / 29 - Y / 12) + 5 * TR * Y1 / 12) + AK[4] *
       TR - (AK[3] - 1) * TR * Y * Y1) / exp((5 / 17)*ln(Z));

      if AK[6] > TR then
         zwt1 := exp(9*ln(AK[6] - TR))
      else         
         zwt1 := -exp(9*ln(TR - AK[6]));
 
      H4 := (AG[2] - AG[4]* TR * TR + AG[5] * (9 * TR + AK[6]) *
       zwt1 + AG[6] * (20 * exp(19*ln(TR)) + AK[7]) /
       exp(2*ln(AK[7] + exp(19*ln(TR))))) * PR;
      H5 := exp(2*ln(12 * exp(11*ln(TR)) + AK[8]) / (AK[8] + exp(11*ln(TR)))) *
       (AG[7] * PR + AG[8] * PR * PR + AG[9] * exp(3*ln(PR)));
      H6 := AG[10] * exp(18*ln(TR)) * (L * AK[9] + M * TR * TR) *
      (exp((-3)*ln(AK[10] + PR)) + AK[11] * PR);
      H7 := AG[11] * AK[12] * exp(3*ln(PR));
      H8 := 21 * AG[12] * exp(4*ln(PR)) / exp(20*ln(TR));
      HR := H1 - H2 + H3 + H4 - H5 + H6 + H7 + H8;
      HSUB1 := HR * Hc;



end;

Function HSUB2(T, P:real):real;
Const
 Tc = 647.3;
 Pc = 221.2;
 Hc = 70.1204;

Var
NA:array[1..8]of integer;
NZ:array[1..8,1..3]of real;
NL:array[1..3]of integer;
NX:array[1..3,1..2]of real;
BGA:array[1..6]of real;
BG:array[1..8,1..3]of real;
BK:array[1..3,1..2]of real;
B9:array[1..9]of real;
BKA:real;
H1,H2,H3,H4,H5,H3A,H4A,H4B,H4C,H5A:real;
TR,PR,X,PRL,DPRL:real;
I,I2,I3,I4,IA, IB,IC,IA1,LAM:integer;
NAA,NAB,NAC:integer;
HIl1,HR,zwt:real;
begin
      NA[1] := 2;
      NA[2] := 3;
      NA[3] := 2;
      NA[4] := 2;
      NA[5] := 3;
      NA[6] := 2;
      NA[7] := 2;
      NA[8] := 2;
      NZ[1, 1] := 13;
      NZ[2, 1] := 18;
      NZ[3, 1] := 18;
      NZ[4, 1] := 25;
      NZ[5, 1] := 32;
      NZ[6, 1] := 12;
      NZ[7, 1] := 24;
      NZ[8, 1] := 24;
      NZ[1, 2] := 3;
      NZ[2, 2] := 2;
      NZ[3, 2] := 10;
      NZ[4, 2] := 14;
      NZ[5, 2] := 28;
      NZ[6, 2] := 11;
      NZ[7, 2] := 18;
      NZ[8, 2] := 14;
      NZ[1, 3] := 0;
      NZ[2, 3] := 1;
      NZ[3, 3] := 0;
      NZ[4, 3] := 0;
      NZ[5, 3] := 24;
      NZ[6, 3] := 0;
      NZ[7, 3] := 0;
      NZ[8, 3] := 0;
      NL[1] := 1;
      NL[2] := 1;
      NL[3] := 2;
      NX[1, 1] := 14;
      NX[2, 1] := 19;
      NX[3, 1] := 54;
      NX[1, 2] := 0;
      NX[2, 2] := 0;
      NX[3, 2] := 27;
      BGA[1] := 16.83599274;
      BGA[2] := 28.56067796;
      BGA[3] := -54.38923329;
      BGA[4] := 0.4330662834;
      BGA[5] := -0.6547711697;
      BGA[6] := 0.08565182058;
      BG[1, 1] := 0.06670375918;
      BG[2, 1] := 0.08390104328;
      BG[3, 1] := 0.4520918904;
      BG[4, 1] := -0.5975336707;
      BG[5, 1] := 0.5958051609;
      BG[6, 1] := 0.1190610271;
      BG[7, 1] := 0.1683998803;
      BG[8, 1] := 0.006552390126;
      BG[1, 2] := 1.388983801;
      BG[2, 2] := 0.02614670893;
      BG[3, 2] := 0.1069036614;
      BG[4, 2] := -0.08847535804;
      BG[5, 2] := -0.5159303373;
      BG[6, 2] := -0.09867174132;
      BG[7, 2] := -0.05809438001;
      BG[8, 2] := 0.0005710218649;
      BG[1, 3] := 0;
      BG[2, 3] := -0.03373439453;
      BG[3, 3] := 0;
      BG[4, 3] := 0;
      BG[5, 3] := 0.2075021122;
      BG[6, 3] := 0;
      BG[7, 3] := 0;
      BG[8, 3] := 0;
      BK[1, 1] := 0.4006073948;
      BK[2, 1] := 0.08636081627;
      BK[3, 1] := -0.8532322921;
      BK[1, 2] := 0;
      BK[2, 2] := 0;
      BK[3, 2] := 0.3460208861;
      B9[1] := 193.6587558;
      B9[2] := -1388.522425;
      B9[3] := 4126.607219;
      B9[4] := -6508.211677;
      B9[5] := 5745.984054;
      B9[6] := -2693.088365;
      B9[7] := 523.5718623;
      BKA := 0.7633333333;
TR:= T/Tc;
PR:= P/Pc;

X:=Exp(BKA*(1-TR));
PRL:=BL(TR);
DPRL:=DBL(TR);

H1:=BGA[1] * TR;
H2:=0;
For I := 2 To 6 do
  begin
    zwt := exp((I-2)*ln(TR));
    H2:= H2 + BGA[I] * (I-3) * zwt;
  end;
H3:= 0;
For I2:= 1 To 5 do
  begin
    H3A:=0;
    NAA:= NA[I2];
    For IA:=1 To NAA do
      begin
        H3A := H3A + (BG[I2,IA] * (1 + NZ[I2,IA]* BKA *TR)) * exp(NZ[I2,IA]*ln(X));
      end;
    H3 := H3 + exp(I2*ln(PR))*H3A;
  end;
H4:= 0;
For I3 := 6 To 8 do
  Begin
    H4A := 0;
    H4B := 0;
    LAM := I3 - 5;
        NAB := NL[LAM];
        For IB := 1 To NAB do
          begin
            H4B:= H4B + NX[LAM,IB] * BK[LAM,IB] * exp(NX[LAM, IB]*ln(X));
          end;
        H4C := 0;
        For IC := 1 To NAB do
          begin
            H4C := H4C + BK[LAM, IC] * exp(NX[LAM, IC]*ln(X));
          end;
        HIL1 := exp((2 - I3)*ln(PR)) + H4C;