Source Code for GEN

Contents - Index

Source Code for GEN_EOS.DLL

The following source code implements the Compress, EnthDep, EntrDep, and FugCoef EES external functions in file GEN_EOS.DLL. This code is written in DELPHI 10.

library GEN_EOS;
{This library uses the Redlich-Kwong-Soave equation of state in reduced coordinates
to provide a general equation of state.}

uses
SysUtils;

const
doExample = -1;

{***********************************************************************}

type
ParamRecPtr = ^ParamRec;
ParamRec =
    record
      Value : Extended;
      Next : ParamRecPtr;
    end;

{***********************************************************************}
{There are 4 functions; names are separated with commas}
procedure DLFNames(Names : PChar); export; stdcall;
begin
StrCopy(Names,'COMPRESS,ENTHDEP,ENTRDEP,FUGCOEF');
end;

{***********************************************************************}
procedure DLPNames(Names : PChar); export; stdcall;
{no DLP procedures}
begin
StrCopy(Names,'');
end;

{***********************************************************************}
{no FDL procedures}
procedure FDLNames(Names : PChar); export; stdcall;
begin
StrCopy(Names,'');
end;

{***************** General functions used by all routines**********************}
function CountValues (P: ParamRecPtr): Integer;
var
    N: Integer;
begin
    N := 0;
    while (P <> nil) do begin
      N := N + 1;
      P := P^.Next
    end;
    Result := N;
end;

procedure GetCheckInputs(FuncName:shortString; ParamPtr:ParamRecPtr; var PString:shortString; var A,B,m,Tr:extended);
    var Ptr: ParamRecPtr;
        Pr,w:extended;
        NArgs:integer;
    begin
       PString:='';
       NArgs:=CountValues(ParamPtr);
       if (NArgs <> 2) and (NArgs <> 3) then begin
         PString := 'Wrong number of arguments for '+FuncName+' function.';
         exit;
       end;
       A:=-999; B:=-999; m:=-999;
       Ptr := ParamPtr;
       Tr := Ptr^.value;
       if (Tr<0) then begin
         PString:='Temperature input to '+FuncName+' is negative.';
         exit;
       end;
       Ptr := Ptr^.next;
       Pr := Ptr^.value;
       if (Pr<0) then begin
         PString:='Pressure input to '+FuncName+' is negative.';
         exit;
       end;
       Ptr := Ptr^.next;
       if (Ptr <> nil) then
         w := Ptr^.Value
       else
         w := 0.1; {acentric factor}
       m := 0.480 + 1.574 * w - 0.176 * w * w;
       A := 0.42747 * sqr(1 + m * (1 - sqrt(Tr))) * Pr / sqr(Tr);
       B := 0.08664 * Pr / Tr;
    end;

function Compressibility(A,B:extended): extended;
    const
      NMax = 50; {maximum number of iterations}
    var
      NIter, NCut: integer;
      F1, F2, Fnew, Z1, Z2, Znew, DeltaZ: extended;

    function F (Z: extended): extended;
      begin {Redlich Kwong Soave eqn of state in reduced form}
        F := Z * Z * Z - Z * Z + Z * (A - B - B * B) - A * B;
      end; {F}

     begin {Compressibility}
       Z1 := 1;{ideal gas}
       F1 := F(Z1);
       Z2 := 0.9 * Z1;
       F2 := F(Z2);
       NIter := 0;
       while (NIter < NMax) and (abs(F2) > 1E-6) do begin
         NIter := NIter + 1;
         DeltaZ := F2 / (F2 - F1) * (Z2 - Z1);
         NCut := 0;
         repeat
           Znew := Z2 - DeltaZ;
           Fnew := F(Znew);
           DeltaZ := DeltaZ / 2;
           NCut := NCut + 1;
         until (NCut > 5) or (abs(Fnew) < abs(F2));
         Z1 := Z2;
         F1 := F2;
         Z2 := Znew;
         F2 := Fnew;
       end; {while}
       if (NIter < NMax) then
         Compressibility := Z2
       else
         Compressibility := -1; {function did not converge.}
    end;
{***********************************************************************}

function COMPRESS (var PString: ShortString; Mode: integer;
       Inputs: ParamRecPtr): extended; export; stdCall;
var
   A,B,Z,m,Tr:extended;
begin
   Compress:=-999;
   if (Mode = doExample) then
     PString := 'Compress(Tr, Pr [,w])'
   else begin
     GetCheckInputs('Compress',Inputs,PString,A,B,m,Tr);
     if (length(PString)>0) then exit; {error}
     Z:=Compressibility(A,B);
     if (Z<0) then
       PString:='Compress function did not converge.'
     else
       Compress:=Z;
    end;
end; {Compress}

{***********************************************************************}
function ENTHDEP (var PString: string; Mode: integer;
    Inputs: ParamRecPtr): extended; export; stdCall;
var
     Z, A, B, m, Tr,temp: extended;
begin {EnthDep}
    EnthDep:=-999;
    if (Mode = doExample) then
      PString := 'EnthDep(Tr, Pr [,w])'
    else begin
      GetCheckInputs('EnthDep',Inputs,PString,A,B,m, Tr);
      if (length(PString)>0) then exit; {error}
      temp := (sqr(1 + m) - m * (1 + m) * sqrt(Tr));
      Z := Compressibility(A,B);
      if (Z > 0) then
        EnthDep := -Tr * (Z - 1) + 0.42747 / 0.08664 * temp * ln((Z + B) / Z)
      else
        PString:='EnthDep function did not converge';
    end; {else}
end; {EnthDep}

function ENTRDEP (var PString: string; Mode: integer;
    Inputs: ParamRecPtr): extended; export; stdCall;
{Returns the entropy departure divided by R using RKS equation of state}
var
    Tr, Z, A, B, m, EnthalpyDep, lnFugacityCoef, temp: extended;
begin
    EntrDep:=-999;
    if (Mode = doExample) then
      PString := 'EntrDep(Tr, Pr [,w])'
    else begin
      GetCheckInputs('EntrDep',Inputs,PString,A,B,m,Tr);
      if (length(PString)>0) then exit; {error}
      Z := Compressibility(A,B);
      if (Z > 0) then begin
        temp := (sqr(1 + m) - m * (1 + m) * sqrt(Tr));
        EnthalpyDep := (Z - 1) - 0.42747 / (0.08664 * Tr) * temp * ln((Z + B) / Z);
                  {EnthalpyDep is non-dimensionalized with RT}
        lnFugacityCoef := Z - 1 - ln(Z - B) - A / B * ln((Z + B) / Z);
        EntrDep := -(EnthalpyDep - lnFugacityCoef);
      end else
        PString:='EntrDep function did not converge';
    end;
end;

function FUGCOEF (var PString: Shortstring; Mode: integer;
    Inputs: ParamRecPtr): extended; export; stdCall;
var
    Tr, Z, A, B, m: extended;
begin {fugCoef}
    FugCoef:=-999;
    if (Mode = doExample) then
       PString := 'FugCoef(Tr, Pr [,w])'
    else begin
      GetCheckInputs('FugCoef',Inputs,PString,A,B,m,Tr);
      if (length(PString)>0) then exit; {error}
      Z := Compressibility(A,B);
      if (Z > 0) then
        fugCoef := exp(z - 1 - ln(Z - B) - A / B * ln((Z + B) / Z))
      else
        PString:='FugCoef function did not converge.';
    end;{else}
end; {fugCoef}

exports
DLFNames,
DLPNames,
FDLNames,
COMPRESS,
ENTHDEP,
ENTRDEP,
FUGCOEF;

begin

end.