water.pas

delphi water math Programning with delphi.

        NAC := NA[I3];
         For IA1 := 1 To NAC do
          begin
            zwt := 1 + NZ[I3, IA1] * BKA * TR - BKA * TR * H4B / HIL1;
            H4A := H4A + BG[I3, IA1] * exp(NZ[I3, IA1]*ln(X)) * zwt;            
          end;
        H4 := H4 + (H4A / HIL1);
       end;
      H5A := 0;
      For I4 := 1 To 7 do
	begin
          H5A := H5A + (1 + TR * (10 * DPRL / PRL + (I4 - 1) * BKA)) *
          B9[I4] * exp((I4 - 1)*ln(X));
        end;

      if PR/PRL > 0 then
         zwt := exp(10*ln(PR / PRL))
      else         
         zwt := exp(10*ln(-PR / PRL));

      H5 := PR * zwt * H5A;
      HR := H1 - H2 - H3 - H4 + H5;
      HSUB2 := HR * Hc;
end;


function BL(TR:real):real;

Var
L:array[1..3] of real;
begin
L[1] := 15.74373327;
L[2] := -34.17061978;
L[3] := 19.31380707;
BL := L[1] + L[2] * TR + L[3]* exp(2*ln(TR));
end;

function DBL(TR:real):real;
{functie voor intern gebruik t.b.v HSUB2}
Var
L:array[1..2] of real;
begin
L[1] := -34.17061978;
L[2] := 19.31380707 ;
DBL := L[1] + L[2] * TR;
end;

Function VSUB2(T, P:real):real;

const
      Tc = 647.3;
      Pc = 221.2;
      Vc = 0.00317;

var
NA:array[1..8]of integer;
NZ:array[1..8,1..3]of integer;
NL:array[1..3]of integer;
NX:array[1..3,1..2]of integer;
BG:array[1..8,1..3]of real;
BK:array[1..3,1..2]of real;
B9:array[1..7]of real;
BKA, RI1:real;
TR,PR,X,PRL,DPRL,VR:real;
S1,S2,S3,S4,S2A,S3A:real;
I,I1,I2,I3,I4:integer;
IA,IB,IC,NAC,NAB,LAM:integer;
S3B,S3HI,S3NEN,S3ZAE,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;
      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;
      RI1 := 4.260321148;
      BKA := 0.7633333333;


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

X := Exp(BKA * (1 - TR));
PRL := BL(TR);
DPRL := DBL(TR);
S1 := RI1 * TR / PR;
S2 := 0;
For I := 1 To 5 do
  begin
    S2A := 0;
    NAB := NA[I];
      For IA := 1 To NAB do
        begin
          S2A := S2A + BG[I, IA] * exp(NZ[I, IA]*ln(X));
        end;
      S2 := S2 + I * exp((I - 1)*ln(PR)) * S2A;
  end;
 S3 := 0;
 For I1 := 6 To 8 do
   begin
    LAM := I1 - 5;
    S3A := 0;
    NAB := NA[I1];
    For IA := 1 To NAB do
     begin
       S3A := S3A + BG[I1, IA] * exp(NZ[I1, IA]*ln(X));
     end;
    S3ZAE := (I1 - 2) * exp((1 - I1)*ln(PR)) * S3A;
    S3B := 0;
    NAC := NL[LAM];
    For IB := 1 To NAC do
      begin
        S3B := S3B + BK[LAM, IB] * exp(NX[LAM, IB]*ln(X));
      end;
    S3HI := exp((2 - I1)*ln(PR)) + S3B;
    S3NEN := S3HI * S3HI;
    S3 := S3 + S3ZAE / S3NEN;
 end;

S4 := 0;
For I2 := 1 To 7 do
 begin
   S4 := S4 + B9[I2]* exp((I2 - 1)*ln(X));
 end;

if PR/PRL > 0 then
    zwt := exp(10*ln(PR / PRL))
else         
    zwt := exp(10*ln(-PR / PRL));

S4 := 11 * zwt * S4;
VR := S1 - S2 - S3 + S4;
VSUB2 := VR * Vc;

end;

Function VSUB1(T, P:real):real;
const
      Tc = 647.3;
      Pc = 221.2;
      Vc = 0.00317;
Var
AK:array[1..12]of real;
AG:array[1..12]of real;
TR,PR,VR:real;
Y,Y1,Z:real;
S1,S2,S3,S4,S5,S6,zwt:real;
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;


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);
S1 := AG[1] * AK[5] / exp((5 / 17)*ln(Z));

if AK[6] > TR then
         zwt := exp(10*ln(AK[6] - TR))
      else         
         zwt := exp(10*ln(TR - AK[6]));


S2 := AG[2] + AG[3] * TR + AG[4] * TR * TR + AG[5] * zwt + AG[6] / (AK[7] + exp(19*ln(TR)));
S3 := (AG[7] + 2 * AG[8] * PR + 3 * AG[9] * PR * PR) / (AK[8] + exp(11*ln(TR)));
S4 := AG[10] *exp(18*ln(TR)) * (AK[9] + TR * TR) * ((-3) / exp(4*ln(AK[10] + PR)) + AK[11]);
S5 := 3 * AG[11] * (AK[12] - TR) * PR * PR;
S6 := 4 * AG[12] * exp(3*ln(PR)) /exp(20*ln(TR));
VR := S1 + S2 - S3 - S4 + S5 + S6;
VSUB1 := VR * Vc;
end;

Function TCON(T, V: real):real;
Var
AK: array[1..12] of real;
BK: array[1..4] of real;
BG: array[1..12] of real;
DK: array[1..12] of real;
CG: array[1..12] of real;
RHO, TR,RHOR,Sum,ALAM0,ALAMQ,DTR,DL1,DL2,DL3,DL4,DL41,DL5,DLAM,HI:real;
Q,R,S:real;
I:integer;
begin
      AK[1] := 0.0102811;
      AK[2] := 0.0299621;
      AK[3] := 0.0156146;
      AK[4] := -0.00422464;
      BK[1] := -0.39707;
      BK[2] := 0.400302;
      BK[3] := 1.06;
      BG[1] := -0.171587;
      BG[2] := 2.39219;
      DK[1] := 0.0701309;
      DK[2] := 0.011852;
      DK[3] := 0.00169937;
      DK[4] := -1.02;
      CG[1] := 0.642857;
      CG[2] := -4.11717;
      CG[3] := -6.17937;
      CG[4] := 0.00308976;
      CG[5] := 0.0822994;
      CG[6] := 10.0932;
 RHO := 1 / V;
 TR := T / 647.3;
 RHOR := RHO / 317.7;

 Sum := 0;
 For I := 1 To 4 do
   begin
     Sum := Sum + AK[I] * exp((I - 1)*ln(TR));
   end;
   ALAM0 := Sqrt(TR) * Sum;

  HI := RHOR + BG[2];
  ALAMQ := BK[1] + BK[2] * RHOR + BK[3] * Exp(BG[1] * HI * HI);

  DTR := Abs(TR - 1) + CG[4];
  Q := 2 + CG[5] * exp((-0.6)*ln(DTR));
  R := Q * TR + 1;
  If TR < 1 Then S := CG[6] * exp((-0.6)*ln(DTR));
  If TR >= 1 Then S := 1 / DTR;

  DL1 := (DK[1] * exp(10*ln(1/ TR)) + DK[2]) * exp(1.8*ln(RHOR));
  DL2 := CG[1] * (1 - exp(2.8*ln(RHOR)));
  DL3 := DK[3] * S * exp(Q*ln(RHOR));
  DL4 := Q / R * (1 - exp(R*ln(RHOR)));
  DL41 := CG[2] * exp(1.5*ln(TR)) + CG[3] * exp(5*ln(1 / RHOR));
  If Abs(DL41) <= 700 Then DL5 := DK[4] * Exp(DL41);
  If Abs(DL41) > 700 Then DL5 := 0;
  DLAM := DL1 * Exp(DL2) + DL3 * Exp(DL4) + DL5;
  TCON := ALAM0 + ALAMQ + DLAM;



end;

Function ETAF(T, V: real):real;
Var
FK:array[1..4]of real;
GK:array[1..6,1..5]of real;
zwt1,zwt2:real;
TUM,VUM,Sur,ETA0:real;

IT,JT:integer;

begin

      FK[1] := 0.0181583;
      FK[2] := 0.0177624;
      FK[3] := 0.0105287;
      FK[4] := -0.0036744;
      GK[1, 1] := 0.501938;
      GK[2, 1] := 0.162888;
      GK[3, 1] := -0.130356;
      GK[4, 1] := 0.907919;
      GK[5, 1] := -0.551119;
      GK[6, 1] := 0.146543;
      GK[1, 2] := 0.235622;
      GK[2, 2] := 0.789393;
      GK[3, 2] := 0.673665;
      GK[4, 2] := 1.207552;
      GK[5, 2] := 0.0670665;
      GK[6, 2] := -0.084337;
      GK[1, 3] := -0.274637;
      GK[2, 3] := -0.743539;
      GK[3, 3] := -0.959456;
      GK[4, 3] := -0.687343;
      GK[5, 3] := -0.497089;
      GK[6, 3] := 0.195286;
      GK[1, 4] := 0.145831;
      GK[2, 4] := 0.263129;
      GK[3, 4] := 0.347247;
      GK[4, 4] := 0.213486;
      GK[5, 4] := 0.100754;
      GK[6, 4] := -0.032932;
      GK[1, 5] := -0.0270448;
      GK[2, 5] := -0.0253093;
      GK[3, 5] := -0.0267758;
      GK[4, 5] := -0.0822904;
      GK[5, 5] := 0.0602253;
      GK[6, 5] := -0.0202595;
TUM := 647.27 / T;
VUM := 0.003147 / V;
Sur := 0;
For IT := 1 To 4 do
begin
Sur := Sur + FK[IT] * exp((IT - 1)*ln(TUM));
end;
ETA0 := Sqrt(1 / TUM) / Sur;
Sur := 0;
For IT := 1 To 6 do
  begin
  if TUM > 1 then
     zwt1 := exp((IT - 1)*ln(TUM - 1))
  else
     begin
        zwt1 := ln(1 - TUM);
         if odd(IT-1) then
           zwt1 := -exp((IT - 1)*zwt1)
         else
           zwt1 := exp((IT - 1)*zwt1);
     end;

  For JT := 1 To 5 do
    begin
    if VUM > 1 then
       zwt2 := exp((JT - 1)*ln(VUM - 1))
    else
       begin
          zwt2 := ln(1 - VUM);
          if odd(JT-1) then
             zwt2 := -exp((JT - 1)*zwt2)
          else
             zwt2 := exp((JT - 1)*zwt2);
       end;
    Sur := Sur + GK[IT, JT] * zwt1 * zwt2;
    end;
  end;
ETAF := ETA0 * Exp(VUM * Sur);

end;

Function HCAP(T,P:real;IAG:integer):real;
Const
     Tc = 647.3;
     Pc = 221.2;
     Vc = 0.00317;
Var
TS1,TS2,HS1,HS2:real;
begin
TS1:= T - 0.5;
TS2:= T + 0.5;

if IAG = 2 then
    begin
      HS1:= HSUB2(TS1,P);
      HS2:= HSUB2(TS2,P);
    end
else
    begin
      HS1:= HSUB1(TS1,P);
      HS2:= HSUB1(TS2,P);
    end;

HCAP := (HS2 - HS1);

end;

end.