source: trunk/eml/stage_separators/tank.mso @ 992

#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*----------------------------------------------------------------------
* Author: Paula B. Staudt
* $Id: tank.mso 981 2016-05-23 12:46:53Z arge $
*--------------------------------------------------------------------*#

using "streams";

Model VesselVolume

ATTRIBUTES
        Pallete         = false;
        Brief           = "Model to calculate vessel total volume and vessel filled volume from 
different geometries and orientations.";
        Info            =
"== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;
";
        
PARAMETERS
        Pi                      as positive             (Brief="Pi value", Default=3.141593,Hidden=true, Symbol="\pi");
        Gconst                  as acceleration (Brief="Gravity Acceleration",Default=9.81,Hidden=true);

        Orientation     as Switcher     (Valid=["vertical","horizontal"],Default="vertical");
        Heads                   as Switcher     (Valid=["elliptical","hemispherical","flat"],Default="flat");
        Diameter                as length               (Brief="Vessel diameter", Default= 2, Symbol="D_{i}");
        Lenght                  as length               (Brief="Side length of the cylinder shell", Default= 6,  Symbol="L_{vessel}");
        
        Vhead_elliptical                as volume               (Brief="Elliptical Head Total Volume",Hidden=true, Symbol="V_{head}^{elliptical}");
        Vhead_hemispherical     as volume               (Brief="Hemispherical Head Total Volume",Hidden=true, Symbol="V_{head}^{hemispherical}");
        Vcylinder                               as volume               (Brief="Cylinder Total Volume",Hidden=true, Symbol="V_{cylinder}");
        radius                                  as length               (Brief="Vessel radius",Hidden=true, Symbol="R_{cylinder}");

SET

        Gconst = 9.81 * 'm/(s^2)';
        Vhead_elliptical        = (Pi*Diameter^3)/12;
        Vhead_hemispherical = (Pi*Diameter^3)/6;
        Vcylinder = 0.25*(Pi*Diameter^2)*Lenght;
        radius = 0.5*Diameter;
        
VARIABLES

        Vtotal          as volume               (Brief="Vessel total volume",Protected=true, Symbol="V_{total}");
        Vfilled         as volume               (Brief="Vessel volume content",Protected=true, Symbol="V_{filled}");
        Level           as length               (Brief="liquid height", Protected=true);
        Across          as area                 (Brief="Vessel cylinder shell Cross section area", Hidden=true, Symbol="A_{cross}");
        
EQUATIONS

switch Orientation

case "vertical":

"Vessel Cross Section Area"
        Across = 0.25*(Pi*Diameter^2);

switch Heads

case "elliptical":

"Vessel Total Volume"
        Vtotal = Vhead_elliptical +     Vcylinder;

if Level < 0.25*Diameter then

"Vessel Filled Volume"
        Vfilled = 0.25*Pi*(((Diameter*Level)/(0.25*Diameter))^2)*(0.25*Diameter-Level/3);

else

"Vessel Filled Volume"
        Vfilled = 0.25*Pi*(Diameter^2)*(Level - 0.25*Diameter/3);

end

case "hemispherical":

"Vessel Total Volume"
        Vtotal = Vhead_hemispherical + Vcylinder;

if Level < 0.5*Diameter then

"Vessel Filled Volume"
        Vfilled = 0.25*Pi*(Level^2)*(2*Diameter-4*Level/3);

else

"Vessel Filled Volume"
        Vfilled = 0.25*Pi*((2/3)*((0.5*Diameter)^3) - (0.25*(Diameter)^3) + Level*Diameter^2);

end

case "flat":

"Vessel Total Volume"
        Vtotal = Vcylinder;

"Vessel Filled Volume"
        Vfilled = Across*Level;

end

case "horizontal":

"Vessel Cross Section Area"
        Across = (radius^2)*acos((radius-Level)/radius)/'rad'-(radius-Level)*sqrt((2*radius*Level-Level^2));

switch Heads

case "elliptical":

"Vessel Total Volume"
        Vtotal = Vhead_elliptical +     Vcylinder;

"Vessel Filled Volume"
        Vfilled = 0.5236*Level^2*(1.5*Diameter-Level) + Across*Lenght;

case "hemispherical":

"Vessel Total Volume"
        Vtotal = Vhead_hemispherical + Vcylinder;

"Vessel Filled Volume"
        Vfilled = 1.0472*Level^2*(1.5*Diameter-Level) + Across*Lenght;

case "flat":

"Vessel Total Volume"
        Vtotal = Vcylinder;

"Vessel Filled Volume"
        Vfilled = Across*Lenght;

end

end

end

Model TankVL

ATTRIBUTES
        Pallete         = true;
        Icon            = "icon/TankVL"; 
        Brief           = "Model of a Tank With Thermodynamic Equilibrium.";
        Info            =
"== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium.

== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;
        
== SPECIFY ==
* the Inlet stream;
* the outlet flows: OutletVapour.F and OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the TankVL model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the Tank temperature (OutletLiquid.T);
* Levelpercent_Initial : the Tank liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
";
        
PARAMETERS
outer PP                as Plugin       (Brief = "External Physical Properties", Type="PP");
outer NComp     as Integer      (Brief = "Number of components", Lower = 1); 

        Mw(NComp)               as molweight    (Brief="Mol Weight", Hidden=true);

        Gconst          as acceleration (Brief="Gravity Acceleration",Default=9.81,Hidden=true);
        low_flow        as flow_mol     (Brief = "Low Flow",Default = 1E-6, Hidden=true);
        zero_flow       as flow_mol     (Brief = "No Flow",Default = 0, Hidden=true);
        KfConst         as Real                 (Brief="Constant for K factor pressure drop", Unit= 'mol/(s*(Pa^0.5))',Default = 1, Hidden=true);
        Kfactor as positive (Brief="K factor for pressure drop", Lower = 1E-8, Default = 2);
        
        NormalFlow      as Switcher     (Brief="Normal Flow", Valid = ["on", "off"], Default = "on",Hidden=true);

        Levelpercent_Initial                    as positive     (Brief="Initial liquid height in Percent", Default = 0.70);
        Temperature_Initial                             as temperature  (Brief="Initial Liquid Temperature", Default = 330);
        Composition_Initial(NComp)              as fraction             (Brief="Initial Composition", Default = 0.10);

SET

        Mw=PP.MolecularWeight();
        Gconst = 9.81 * 'm/(s^2)';
        low_flow = 1E-6 * 'kmol/h';
        zero_flow = 0 * 'kmol/h';
        KfConst = 1*'mol/(s*(Pa^0.5))';

VARIABLES

        Geometry                as VesselVolume (Brief="Vessel Geometry", Symbol=" ");

in      Inlet                   as stream                       (Brief="Feed Stream", PosX=0.22, PosY=0, Symbol="_{in}");
out     OutletLiquid    as liquid_stream        (Brief="Liquid outlet stream", PosX=0.43, PosY=1, Symbol="_{out}^{Liquid}");
out     OutletVapour    as vapour_stream        (Brief="Vapour outlet stream", PosX=0.68, PosY=0, Symbol="_{out}^{Vapour}");
in      InletQ                  as power                        (Brief="Heat Duty", PosX=0.735, PosY=1, Protected =true,Symbol="Q_{in}");

        TotalHoldup(NComp)              as mol  (Brief="Molar Holdup in the Vessel", Protected=true);
        LiquidHoldup                    as mol  (Brief="Molar liquid holdup", Protected=true);
        VapourHoldup                    as mol  (Brief="Molar vapour holdup", Protected=true);
        
        E                       as energy               (Brief="Total Energy Holdup in the Vessel", Protected=true);
        vL                      as volume_mol   (Brief="Liquid Molar Volume", Protected=true);
        vV                      as volume_mol   (Brief="Vapour Molar volume", Protected=true);
        Pdrop           as press_delta  (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P", Protected=true);
        Peq                     as pressure             (Brief="Equilibrium pressure on the liquid surface", Protected=true, Symbol="\Delta P_{eq}");
        Pstatic         as pressure             (Brief="Static head at the bottom of the tank", Protected = true, Symbol="P_{static}^{Liquid}");

out     TI as control_signal    (Brief="Temperature Indicator", PosX=0.525, PosY=0, Protected=true);
out     PI as control_signal    (Brief="Pressure Indicator", PosX=0.368, PosY=0, Protected=true);
out     LI as control_signal    (Brief="Level Indicator", PosX=1, PosY=0.6, Protected=true);

INITIAL

"Initial level Percent"
        LI = Levelpercent_Initial;
        
"Initial Outlet Liquid Temperature"
        OutletLiquid.T = Temperature_Initial;
        
"Initial Outlet Liquid Composition Normalized"
        OutletLiquid.z(1:NComp - 1) = Composition_Initial(1:NComp - 1)/sum(Composition_Initial);

EQUATIONS

switch NormalFlow

case "on":
        Inlet.F = Kfactor *sqrt(Pdrop)*KfConst;

        when Inlet.F < low_flow switchto "off";

case "off":
        Inlet.F = zero_flow;

        when Inlet.P > OutletVapour.P switchto "on";

end

"Component Molar Balance"
        diff(TotalHoldup)=Inlet.F*Inlet.z - OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z;
        
"Energy Balance"
        diff(E) = Inlet.F*Inlet.h - OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h + InletQ;
        
"Molar Holdup"
        TotalHoldup = LiquidHoldup*OutletLiquid.z + VapourHoldup*OutletVapour.z; 
        
"Energy Holdup"
        E = LiquidHoldup*OutletLiquid.h + VapourHoldup*OutletVapour.h - OutletLiquid.P*Geometry.Vtotal;
        
"Mol fraction normalisation"
        sum(OutletLiquid.z)=1.0;

"Mol fraction normalisation"
        sum(OutletLiquid.z)=sum(OutletVapour.z);

"Liquid Volume"
        vL = PP.LiquidVolume(OutletLiquid.T, Peq, OutletLiquid.z);

"Vapour Volume"
        vV = PP.VapourVolume(OutletVapour.T, Peq, OutletVapour.z);
        
"Chemical Equilibrium"
        PP.LiquidFugacityCoefficient(OutletLiquid.T, Peq, OutletLiquid.z)*OutletLiquid.z = 
                PP.VapourFugacityCoefficient(OutletVapour.T, Peq, OutletVapour.z)*OutletVapour.z;
        
"Thermal Equilibrium"
        OutletVapour.T = OutletLiquid.T;
        
"Mechanical Equilibrium"
        OutletVapour.P = Peq;

"Static Head"   
        Pstatic = PP.LiquidDensity(OutletLiquid.T, Peq, OutletLiquid.z) * Gconst * Geometry.Level;

"Mechanical Equilibrium for the Liquid Phase"
        OutletLiquid.P = Peq + Pstatic; 

"Pressure Drop"
        Pdrop = Inlet.P - OutletVapour.P;

"Geometry Constraint"
        Geometry.Vtotal = LiquidHoldup * vL + VapourHoldup * vV;

"Temperature indicator"
        TI * 'K' = OutletLiquid.T;

"Pressure indicator"
        PI * 'atm' = OutletLiquid.P;

"Level indicator"
        LI*Geometry.Vtotal= Geometry.Vfilled;

"Liquid Level"
        LiquidHoldup * vL = Geometry.Vfilled;

end

Model TankL

ATTRIBUTES
        Pallete         = true;
        Icon            = "icon/TankL"; 
        Brief           = "Model of a Tank.";
        Info            =
"== ASSUMPTIONS ==
* liquid phase only;

== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;
        
== SPECIFY ==
* the Inlet stream;
* the OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the TankL model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the Tank temperature (OutletLiquid.T);
* Levelpercent_Initial : the Tank liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
";
        
PARAMETERS
outer PP                as Plugin       (Brief = "External Physical Properties", Type="PP");
outer NComp     as Integer      (Brief = "Number of components", Lower = 1); 

        Gconst                  as acceleration (Brief="Gravity Acceleration",Default=9.81,Hidden=true);
        low_flow        as flow_mol     (Brief = "Low Flow",Default = 1E-6, Hidden=true);
        zero_flow       as flow_mol     (Brief = "No Flow",Default = 0, Hidden=true);
        KfConst         as Real                 (Brief="Constant for K factor pressure drop", Unit= 'mol/(s*(Pa^0.5))',Default = 1, Hidden=true);
        Kfactor as positive (Brief="K factor for pressure drop", Lower = 1E-8, Default = 2);
        
        NormalFlow      as Switcher     (Brief="Normal Flow", Valid = ["on", "off"], Default = "on",Hidden=true);


        Levelpercent_Initial                    as positive     (Brief="Initial liquid height in Percent", Default = 0.70);
        Temperature_Initial                             as temperature  (Brief="Initial Liquid Temperature", Default = 330);
        Composition_Initial(NComp)              as fraction             (Brief="Initial Composition", Default = 0.10);

SET

        Gconst = 9.81 * 'm/(s^2)';
        low_flow = 1E-6 * 'kmol/h';
        zero_flow = 0 * 'kmol/h';
        KfConst = 1*'mol/(s*(Pa^0.5))';

VARIABLES

        Geometry                as VesselVolume (Brief="Vessel Geometry", Symbol=" ");

in      Inlet                   as stream                       (Brief="Feed Stream", PosX=0.22, PosY=0, Symbol="_{in}");
out     OutletLiquid    as liquid_stream        (Brief="Liquid outlet stream", PosX=0.43, PosY=1, Symbol="_{out}^{Liquid}");
in      InletQ                  as power                        (Brief="Heat Duty", PosX=0.735, PosY=1, Protected =true,Symbol="Q_{in}");

        TotalHoldup(NComp)              as mol  (Brief="Molar Holdup in the Vessel", Protected=true);
        
        E                       as energy               (Brief="Total Energy Holdup in the Vessel", Protected=true);
        vL                      as volume_mol   (Brief="Liquid Molar Volume", Protected=true);
        Pstatic         as pressure             (Brief="Static head at the bottom of the tank", Protected = true, Symbol="P_{static}^{Liquid}");

out     TI as control_signal    (Brief="Temperature Indicator", PosX=0.525, PosY=0, Protected=true);
out     PI as control_signal    (Brief="Pressure Indicator", PosX=0.368, PosY=0, Protected=true);
out     LI as control_signal    (Brief="Level Indicator", PosX=1, PosY=0.6, Protected=true);

INITIAL

"Initial level Percent"
        LI = Levelpercent_Initial;
        
"Initial Outlet Liquid Temperature"
        OutletLiquid.T = Temperature_Initial;
        
"Initial Outlet Liquid Composition Normalized"
        OutletLiquid.z(1:NComp - 1) = Composition_Initial(1:NComp - 1)/sum(Composition_Initial);

EQUATIONS

#*
switch NormalFlow

case "on":
        Inlet.F = Kfactor *sqrt(Pstatic)*KfConst;

        when Inlet.F < low_flow switchto "off";

case "off":
        Inlet.F = zero_flow;

        when Inlet.P > OutletLiquid.P switchto "on";

end
*#
"Component Molar Balance"
        diff(TotalHoldup)=Inlet.F*Inlet.z - OutletLiquid.F*OutletLiquid.z;
        
"Energy Balance"
        diff(E) = Inlet.F*Inlet.h - OutletLiquid.F*OutletLiquid.h + InletQ;

"Energy Holdup"
        E = sum(TotalHoldup)*OutletLiquid.h;

"Static Head"   
        Pstatic = PP.LiquidDensity(OutletLiquid.T, Inlet.P, OutletLiquid.z) * Gconst * Geometry.Level;

"Mechanical Equilibrium"
        Inlet.P + Pstatic = OutletLiquid.P;

"Liquid Volume"
        vL = PP.LiquidVolume(OutletLiquid.T, OutletLiquid.P, OutletLiquid.z);

"Molar Holdup"
        TotalHoldup = OutletLiquid.z*sum(TotalHoldup); 
        
"Liquid Level"
        Geometry.Vfilled = sum(TotalHoldup) * vL;
        
"Temperature indicator"
        TI * 'K' = OutletLiquid.T;

"Pressure indicator"
        PI * 'atm' = OutletLiquid.P;

"Level indicator"
        LI*Geometry.Vtotal= Geometry.Vfilled;

end

Model SumpTank

ATTRIBUTES
        Pallete         = true;
        Icon            = "icon/SumpTank"; 
        Brief           = "Model of a tank with 2 material inlet streams and with thermodynamic equilibrium.";
        Info            =
"== ASSUMPTIONS ==
* perfect mixing of both phases;
* thermodynamics equilibrium.

== SET ==
*Orientation: vessel position - vertical or horizontal;
*Heads (bottom and top heads are identical)
**elliptical: 2:1 elliptical heads (25% of vessel diameter);
**hemispherical: hemispherical heads (50% of vessel diameter);
**flat: flat heads (0% of vessel diameter);
*Diameter: Vessel diameter;
*Lenght: Side length of the cylinder shell;
        
== SPECIFY ==
* the Inlet stream;
* the outlet flows: OutletVapour.F and OutletLiquid.F;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model).

== OPTIONAL ==
* the SumpTank model has three control ports
** TI OutletLiquid Temperature Indicator;
** PI OutletLiquid Pressure Indicator;
** LI Level Indicator;

== INITIAL CONDITIONS ==
* Initial_Temperature :  the Tank temperature (OutletLiquid.T);
* Levelpercent_Initial : the Tank liquid level in percent (LI);
* Initial_Composition : (NoComps) OutletLiquid compositions.
";
        
PARAMETERS
outer PP                as Plugin       (Brief = "External Physical Properties", Type="PP");
outer NComp     as Integer      (Brief = "Number of components", Lower = 1); 

        Mw(NComp)               as molweight    (Brief="Mol Weight", Hidden=true);
        pi                      as positive             (Brief="Pi value", Default=3.141593,Hidden=true, Symbol="\pi");
        g                               as acceleration         (Brief="Gravity Acceleration",Default=9.81,Hidden=true);
        
        Heads                   as Switcher     (Valid=["elliptical","hemispherical","flat"],Default="flat");
        Diameter                as length               (Brief="Vessel diameter", Symbol="D_{i}");
        Lenght                  as length               (Brief="Side length of the cylinder shell", Symbol="L_{vessel}");
        
        Vhead_elliptical                as volume               (Brief="Elliptical Head Total Volume",Hidden=true, Symbol="V_{head}^{elliptical}");
        Vhead_hemispherical     as volume               (Brief="Hemispherical Head Total Volume",Hidden=true, Symbol="V_{head}^{hemispherical}");
        Vcylinder                               as volume               (Brief="Cylinder Total Volume",Hidden=true, Symbol="V_{cylinder}");
        radius                                  as length               (Brief="Vessel radius",Hidden=true, Symbol="R_{cylinder}");
        
        Levelpercent_Initial                    as positive     (Brief="Initial liquid height in Percent", Default = 0.70);
        Temperature_Initial                             as temperature  (Brief="Initial Liquid Temperature", Default = 330);
        Composition_Initial(NComp)              as fraction             (Brief="Initial Composition", Default = 0.10);

SET

        Mw=PP.MolecularWeight();

        g = 9.81 * 'm/(s^2)';
        Vhead_elliptical        = (pi*Diameter^3)/12;
        Vhead_hemispherical = (pi*Diameter^3)/6;
        Vcylinder = 0.25*(pi*Diameter^2)*Lenght;
        radius = 0.5*Diameter;

VARIABLES

in      InletLiquid     as stream                       (Brief="Feed Stream", PosX=0.22, PosY=0, Symbol="_{in}");
out     OutletLiquid    as liquid_stream        (Brief="Liquid outlet stream", PosX=0.50, PosY=1, Symbol="_{out}^{Liquid}");
in      InletVapour     as stream                       (Brief="Vapour outlet stream", PosX=1, PosY=0.20, Symbol="_{out}^{Vapour}");
out     OutletVapour    as vapour_stream        (Brief="Vapour outlet stream", PosX=0.68, PosY=0, Symbol="_{out}^{Vapour}");
        InletQ                  as power                        (Brief="Heat Duty", Protected =false,Symbol="Q_{in}");

        Vtotal                  as volume                       (Brief="Vessel total volume",Protected=true, Symbol="V_{total}");
        Vfilled                 as volume                       (Brief="Vessel volume content",Protected=true, Symbol="V_{filled}");

        TotalHoldup(NComp)              as mol  (Brief="Molar Holdup in the Vessel", Protected=true);
        LiquidHoldup                    as mol  (Brief="Molar liquid holdup", Protected=true);
        VapourHoldup                    as mol  (Brief="Molar vapour holdup", Protected=true);
        
        E                       as energy               (Brief="Total Energy Holdup in the Vessel", Protected=true);
        vL                      as volume_mol   (Brief="Liquid Molar Volume", Protected=true);
        vV                      as volume_mol   (Brief="Vapour Molar volume", Protected=true);
        Level           as length               (Brief="liquid height", Protected=true);
        Across          as area                 (Brief="Vessel cylinder shell Cross section area", Hidden=true, Symbol="A_{cross}");
        #Pdrop          as press_delta  (Brief = "Pressure Drop", DisplayUnit = 'kPa', Symbol ="\Delta P", Protected=true);

        Peq             as pressure             (Brief="Equilibrium pressure on the liquid surface", Protected=true, Symbol="P_{eq}");
        Pstatic         as pressure             (Brief="Static head at the bottom of the tank", Protected = true, Symbol="P_{static}^{Liquid}");

out     LI as control_signal    (Brief="Level Indicator", PosX=1, PosY=0.73, Protected=true);
out     TI as control_signal    (Brief="Temperature Indicator", PosX=1, PosY=0.40, Protected=true);

INITIAL

"Initial level Percent"
        LI = Levelpercent_Initial;
        
"Initial Outlet Liquid Temperature"
        OutletLiquid.T = Temperature_Initial;
        
"Initial Outlet Liquid Composition Normalized"
        OutletLiquid.z(1:NComp - 1) = Composition_Initial(1:NComp - 1)/sum(Composition_Initial);

EQUATIONS

"Vessel Cross Section Area"
        Across = 0.25*(pi*Diameter^2);

switch Heads

case "elliptical":

"Vessel Total Volume"
        Vtotal = Vhead_elliptical +     Vcylinder;

if Level < 0.25*Diameter then

"Vessel Filled Volume"
        Vfilled = 0.25*pi*(((Diameter*Level)/(0.25*Diameter))^2)*(0.25*Diameter-Level/3);

else

"Vessel Filled Volume"
        Vfilled = 0.25*pi*(Diameter^2)*(Level - 0.25*Diameter/3);

end

case "hemispherical":

"Vessel Total Volume"
        Vtotal = Vhead_hemispherical + Vcylinder;

if Level < 0.5*Diameter then

"Vessel Filled Volume"
        Vfilled = 0.25*pi*(Level^2)*(2*Diameter-4*Level/3);

else

"Vessel Filled Volume"
        Vfilled = 0.25*pi*((2/3)*((0.5*Diameter)^3) - (0.25*(Diameter)^3) + Level*Diameter^2);

end

case "flat":

"Vessel Total Volume"
        Vtotal = Vcylinder;

"Vessel Filled Volume"
        Vfilled = Across*Level;

end

"Component Molar Balance"
        diff(TotalHoldup) = InletLiquid.F*InletLiquid.z + InletVapour.F*InletVapour.z- OutletLiquid.F*OutletLiquid.z - OutletVapour.F*OutletVapour.z;
        
"Energy Balance"
        diff(E) = InletLiquid.F*InletLiquid.h + InletVapour.F*InletVapour.h - OutletLiquid.F*OutletLiquid.h - OutletVapour.F*OutletVapour.h + InletQ;
        
"Molar Holdup"
        TotalHoldup = LiquidHoldup*OutletLiquid.z + VapourHoldup*OutletVapour.z; 
        
"Energy Holdup"
        E = LiquidHoldup*OutletLiquid.h + VapourHoldup*OutletVapour.h - OutletLiquid.P*Vtotal;
        
"Mol fraction normalisation"
        sum(OutletLiquid.z)=1.0;

"Mol fraction normalisation"
        sum(OutletLiquid.z)=sum(OutletVapour.z);

"Liquid Volume"
        vL = PP.LiquidVolume(OutletLiquid.T, Peq, OutletLiquid.z);

"Vapour Volume"
        vV = PP.VapourVolume(OutletVapour.T, Peq, OutletVapour.z);
        
"Chemical Equilibrium"
        PP.LiquidFugacityCoefficient(OutletLiquid.T, Peq, OutletLiquid.z)*OutletLiquid.z = 
                PP.VapourFugacityCoefficient(OutletVapour.T, Peq, OutletVapour.z)*OutletVapour.z;
        
"Thermal Equilibrium"
        OutletVapour.T = OutletLiquid.T;
        
"Mechanical Equilibrium for the Vapour Phase"
        OutletVapour.P = Peq;
        
"Static Head"   
        Pstatic = PP.LiquidDensity(OutletLiquid.T, Peq, OutletLiquid.z) * g * Level;

"Mechanical Equilibrium for the Liquid Phase"
        OutletLiquid.P = Peq + Pstatic; 

#*"Pressure Drop"
        Pdrop = min([InletLiquid.P,InletVapour.P]) - OutletLiquid.P;
        #OutletLiquid.P  = InletLiquid.P - Pdrop;
*#

"Geometry Constraint"
        Vtotal = LiquidHoldup * vL + VapourHoldup * vV;

"Level indicator"
        LI*Vtotal= Vfilled;

"Temperature indicator"
        TI * 'K' = OutletVapour.T;

"Liquid Level"
        LiquidHoldup * vL = Vfilled;

end

Model tank_simplified
        ATTRIBUTES
        Pallete         = true;
        Icon            = "icon/Tank"; 
        Brief           = "Model of a simplified tank.";
        Info            =
"== Specify ==
* the Inlet flow rate;

== Initial Conditions ==
* the tank initial level (Level);
";

        PARAMETERS
        k as Real (Brief="Valve Constant", Unit = 'm^2.5/h', Default=4);
        A as area (Brief="Tank area", Default=2);

        VARIABLES
        Level as length(Brief="Tank level");
in      Fin  as flow_vol(Brief="Input flow", PosX=0.3037, PosY=0);
out     Fout as flow_vol(Brief="Output flow", PosX=1, PosY=1);

        EQUATIONS
        "Mass balance"
        diff(A*Level) = Fin - Fout;

        "Valve equation"
        Fout = k*sqrt(Level);           
end