Changeset 497 for branches/packed

Timestamp:

Apr 11, 2008, 2:03:21 PM (14 years ago)

Author:

Paula Bettio Staudt

Message:

Problem solved: model of a packed columnstage and its examples working

Location:

branches/packed

Files:

• eml/stage_separators/PackedStage.mso (modified) (3 diffs)
• eml/stage_separators/column.mso (modified) (4 diffs)
• eml/stage_separators/tray.mso (modified) (3 diffs)
• eml/streams.mso (modified) (1 diff)
• sample/stage_separators/sample_column.mso (modified) (6 diffs)
• sample/stage_separators/sample_tray.mso (modified) (2 diffs)

branches/packed/eml/stage_separators/PackedStage.mso

@@ -50 +50 @@
         uV as velocity (Brief="volume flow rate of vapor, m^3/m^2/s", Lower=-10, Upper=100);
         dp as length (Brief="Particle diameter", Default=1e-3, Lower=0, Upper=10);
-        invK as positive (Brief="Wall factor");
-        Rev as Real (Brief="Reynolds number of the vapor stream", Default=100);
-        Al as area;
-        hl as Real (Brief="Column holdup", Unit='m^3/m^3');
+        invK as positive (Brief="Wall factor", Default=1, Upper=10);
+        Rev as Real (Brief="Reynolds number of the vapor stream", Default=4000);
+        Al as area (Brief="Area occupied by the liquid", Default=0.001, Upper=1);
+        hl as positive (Brief="Column holdup", Unit='m^3/m^3', Default=0.01,Upper=10);
 #       hlS as Real (Brief="Column holdup at loading point", Unit='m^3/m^3');
 #       hlFl as Real (Brief="Column holdup at flooding point", Unit='m^3/m^3');
@@ -60 +60 @@
 #       uLS as velocity;
 #       uLFl as velocity;
-        Qsil as Real (Brief="Resistance coefficient on the liquid load", Lower = 0);
+        Qsil as positive (Brief="Resistance coefficient on the liquid load");
 #       QsiS as Real (Brief="Resistance coefficient at the loading point", Lower = 0);
 #       QsiFl as Real (Brief="Resistance coefficient at the flooding point", Lower = 0);
@@ -232 +232 @@
         #t1.OutletV.F = 180.1 * 'kmol/h';
         t1.deltaP = 0.00189536 * 'atm';
-        t1.Qsil = 1;
+        t1.Qsil = 5;
         
         SET

branches/packed/eml/stage_separators/column.mso

@@ -1281 +1281 @@
         bot = NStages/top;
         stage.hs = H/NStages;   
-        stage.V = stage.hs * stage.Ap;
+        stage.V = stage.hs * stage.d^2*3.14159/4;
         
         CONNECTIONS
@@ -1288 +1288 @@
 end
 
-Model PackedDistillation_kettle_cond
-        PARAMETERS
-        outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
-        outer NComp as Integer;
-        NStages as Integer(Brief="Number of stages", Default=2);
-        topdown as Integer(Brief="Trays counting (1=top-down, -1=bottom-up)", Default=1);
-        top as Integer(Brief="Number of top tray");
-        bot as Integer(Brief="Number of bottom tray");
-        H as length (Brief="Height of packing");
-        K as Real (Brief="Reboiler flow constant", Unit='mol*s^0.5/kg^0.5/m');
-        
-        VARIABLES
-        stage(NStages) as packedStage_BilletSchultes;
+Model PackedDistillation_kettle_cond as Section_Column_Packed
+        PARAMETERS
+        VapourFlow as Switcher(Valid = ["on", "off"], Default = "on");
+        
+        VARIABLES
         cond as condenser;
         reb as reboiler;
@@ -1306 +1298 @@
         pump1 as pump;
         
-        SET
-        top = (NStages-1)*(1-topdown)/2+1;
-        bot = NStages/top;
-        stage.hs = H/NStages;   
-        stage.V = stage.hs * stage.Ap;
-        
         CONNECTIONS
         #vapor
         reb.OutletV to stage(bot).InletV;
-        stage([top+topdown:topdown:bot]).OutletV to stage([top:topdown:bot-topdown]).InletV;
         stage(top).OutletV to cond.InletV;
         
@@ -1322 +1307 @@
         sptop.Outlet2 to pump1.Inlet;
         pump1.Outlet to stage(top).InletL;
-        stage([top:topdown:bot-topdown]).OutletL to stage([top+topdown:topdown:bot]).InletL;
         stage(bot).OutletL to reb.InletL;
-end
+        
+        EQUATIONS
+        switch VapourFlow
+                case "on":
+                stage(bot).InletV.F*stage(bot).vV = sqrt((reb.OutletV.P - stage(bot).OutletV.P)/
+                                        (stage(bot).rhoV*stage(bot).Qsil*20))*(stage(bot).d^2*3.14159/4*stage(bot).e - stage(bot).Al);
+                when stage(bot).InletV.F < 1e-6 * 'kmol/h' switchto "off";
+                
+                case "off":
+                stage(bot).InletV.F = 0 * 'mol/s';
+                when reb.OutletV.P > stage(bot).OutletV.P + 1e-1 * 'atm' switchto "on";
+        end
+
+end

branches/packed/eml/stage_separators/tray.mso

@@ -323 +323 @@
         V = ML* vL + MV*vV;
 end
-        
+
+#*-------------------------------------
+* Model of a packed column stage
+-------------------------------------*#
+Model packedStage
+        PARAMETERS
+outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
+outer NComp as Integer;
+        PPwater as Plugin(Brief="Physical Properties",
+                Type="PP",
+                Components = [ "water" ],
+                LiquidModel = "PR",
+                VapourModel = "PR"
+        );
+
+        V as volume(Brief="Total Volume of the tray");
+        Q as heat_rate (Brief="Rate of heat supply"); 
+        d as length (Brief="Column diameter");  
+
+        a as Real (Brief="surface area per packing volume", Unit='m^2/m^3');
+        g as acceleration;
+        e as Real (Brief="Void fraction of packing, m^3/m^3");
+        Cpo as Real (Brief="Constant for resitance equation"); # Billet and Schultes, 1999.
+        Mw(NComp)       as molweight    (Brief = "Component Mol Weight");
+        hs as length (Brief="Height of the packing stage");
+
+        VARIABLES
+in      Inlet as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
+in      InletL as stream (Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
+in      InletV as stream (Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
+out     OutletL as liquid_stream (Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
+out     OutletV as vapour_stream (Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");
+
+        M(NComp) as mol (Brief="Molar Holdup in the tray", Default=0.01, Lower=0, Upper=100);
+        ML as mol (Brief="Molar liquid holdup", Default=0.01, Lower=0, Upper=100);
+        MV as mol (Brief="Molar vapour holdup", Default=0.01, Lower=0, Upper=100);
+        E as energy (Brief="Total Energy Holdup on tray", Default=-500);
+        vL as volume_mol (Brief="Liquid Molar Volume");
+        vV as volume_mol (Brief="Vapour Molar volume");
+        
+        miL as viscosity (Brief="Liquid dynamic viscosity", DisplayUnit='kg/m/s');
+        miV as viscosity (Brief="Vapor dynamic viscosity", DisplayUnit='kg/m/s');
+        rhoL as dens_mass;
+        rhoV as dens_mass;
+        
+        deltaP as pressure;
+        
+        uL as velocity (Brief="volume flow rate of liquid, m^3/m^2/s", Lower=-10, Upper=100);
+        uV as velocity (Brief="volume flow rate of vapor, m^3/m^2/s", Lower=-10, Upper=100);
+        dp as length (Brief="Particle diameter", Default=1e-3, Lower=0, Upper=10);
+        invK as positive (Brief="Wall factor", Default=1, Upper=10);
+        Rev as Real (Brief="Reynolds number of the vapor stream", Default=4000);
+        Al as area (Brief="Area occupied by the liquid", Default=0.001, Upper=1);
+        hl as positive (Brief="Column holdup", Unit='m^3/m^3', Default=0.01,Upper=10);
+        Qsil as positive (Brief="Resistance coefficient on the liquid load");
+
+        SET
+        Mw = PP.MolecularWeight();
+
+        EQUATIONS
+        "Component Molar Balance"
+        diff(M)=Inlet.F*Inlet.z + InletL.F*InletL.z + InletV.F*InletV.z
+                - OutletL.F*OutletL.z - OutletV.F*OutletV.z;
+
+        "Energy Balance"
+        diff(E) = ( Inlet.F*Inlet.h + InletL.F*InletL.h + InletV.F*InletV.h
+                - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q );
+        
+        "Molar Holdup"
+        M = ML*OutletL.z + MV*OutletV.z;
+        
+        "Energy Holdup"
+        E = ML*OutletL.h + MV*OutletV.h - OutletL.P*V;
+        
+        "Mol fraction normalisation"
+        sum(OutletL.z)= 1.0;
+        
+        "Liquid Volume"
+        vL = PP.LiquidVolume(OutletL.T, OutletL.P, OutletL.z);
+        "Vapour Volume"
+        vV = PP.VapourVolume(OutletV.T, OutletV.P, OutletV.z);
+        
+        "Chemical Equilibrium"
+        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z = 
+                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
+        
+        "Thermal Equilibrium"
+        OutletV.T = OutletL.T;
+        
+        "Mechanical Equilibrium"
+        OutletL.P = OutletV.P;
+        
+        "Geometry Constraint"
+        V*e = ML*vL + MV*vV;
+        
+        "Liquid Density"
+        rhoL = PP.LiquidDensity(OutletL.T, OutletL.P, OutletL.z);
+        "Vapour Density"
+        rhoV = PP.VapourDensity(InletV.T, InletV.P, InletV.z);
+        "Liquid viscosity"
+        miL = PP.LiquidViscosity(OutletL.T, OutletL.P, OutletL.z);
+        "Vapour viscosity"
+        miV = PP.VapourViscosity(InletV.T, InletV.P, InletV.z);
+
+        "Area occupied by the liquid"
+        Al = ML*vL/hs;
+
+        "Volume flow rate of liquid, m^3/m^2/s"
+        uL * Al = OutletL.F * vL;
+        "Volume flow rate of vapor, m^3/m^2/s"
+        uV * ((d^2*3.14159/4)*e - Al) = OutletV.F * vV;
+        
+        "Liquid holdup"
+        hl = ML*vL/V/e;
+        
+        "Particle diameter"
+        dp = 6 * (1-e)/a;
+        
+        "Wall Factor"
+        invK = (1 + (2*dp/(3*d*(1-e))));
+        
+        "Reynolds number of the vapor stream"
+        Rev*invK = dp*uV*rhoV / (miV*(1-e));
+        
+        deltaP = InletV.P - OutletV.P;
+        
+        "Pressure drop and Vapor flow"
+        deltaP/hs  = Qsil*a*uV^2*rhoV*invK / (2*(e-hl)^3);
+
+        "Liquid holdup"
+        hl = (12*miL*a^2*uL/rhoL/g)^1/3;
+end
+
+
+#Models not working!!!!
+
 Model packedStage_Navaes as trayBasic
         PARAMETERS
@@ -571 +706 @@
 end
 
-Model packedStage
-        PARAMETERS
-outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
-outer NComp as Integer;
-        PPwater as Plugin(Brief="Physical Properties",
-                Type="PP",
-                Components = [ "water" ],
-                LiquidModel = "PR",
-                VapourModel = "PR"
-        );
-
-        V as volume(Brief="Total Volume of the tray");
-        Q as heat_rate (Brief="Rate of heat supply"); 
-        Ap as area (Brief="Plate area = Atray - Adowncomer");
-
-        a as Real (Brief="surface area per packing volume", Unit='m^2/m^3');
-        g as acceleration;
-        e as Real (Brief="Void fraction of packing, m^3/m^3");
-        ds as length (Brief="Column diameter");
-        Cpo as Real (Brief="Constant for resitance equation"); # Billet and Schultes, 1999.
-        Mw(NComp)       as molweight    (Brief = "Component Mol Weight");
-        hs as length (Brief="Height of the packing stage");
-
-        VARIABLES
-in      Inlet as stream (Brief="Feed stream", PosX=0, PosY=0.4932, Symbol="_{in}");
-in      InletL as stream (Brief="Inlet liquid stream", PosX=0.5195, PosY=0, Symbol="_{inL}");
-in      InletV as stream (Brief="Inlet vapour stream", PosX=0.4994, PosY=1, Symbol="_{inV}");
-out     OutletL as liquid_stream (Brief="Outlet liquid stream", PosX=0.8277, PosY=1, Symbol="_{outL}");
-out     OutletV as vapour_stream (Brief="Outlet vapour stream", PosX=0.8043, PosY=0, Symbol="_{outV}");
-
-        M(NComp) as mol (Brief="Molar Holdup in the tray");
-        ML as mol (Brief="Molar liquid holdup");
-        MV as mol (Brief="Molar vapour holdup");
-        E as energy (Brief="Total Energy Holdup on tray");
-        vL as volume_mol (Brief="Liquid Molar Volume");
-        vV as volume_mol (Brief="Vapour Molar volume");
-        
-        miL as viscosity (Brief="Liquid dynamic viscosity", DisplayUnit='kg/m/s');
-        miV as viscosity (Brief="Vapor dynamic viscosity", DisplayUnit='kg/m/s');
-        rhoL as dens_mass;
-        rhoV as dens_mass;
-        
-        deltaP as pressure;
-        
-        uL as velocity (Brief="volume flow rate of liquid, m^3/m^2/s", Lower = -10, Default = 0.007);
-        uV as velocity (Brief="volume flow rate of vapor, m^3/m^2/s", Lower = -10, Default = 1.14);
-        dp as length (Brief="Particle diameter", Default=1e-3);
-        invK as Real (Brief="Wall factor");
-        Rev as Real (Brief="Reynolds number of the vapor stream", Lower = 0, Default=100);
-        Al as area;
-        hl as Real (Brief="Column holdup", Unit='m^3/m^3');
-        Qsil as Real (Brief="Resistance coefficient on the liquid load", Lower = 0);
-
-        SET
-        Mw = PP.MolecularWeight();
-
-        EQUATIONS
-        "Component Molar Balance"
-        diff(M)=Inlet.F*Inlet.z + InletL.F*InletL.z + InletV.F*InletV.z
-                - OutletL.F*OutletL.z - OutletV.F*OutletV.z;
-
-        "Energy Balance"
-        diff(E) = ( Inlet.F*Inlet.h + InletL.F*InletL.h + InletV.F*InletV.h
-                - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q );
-        
-        "Molar Holdup"
-        M = ML*OutletL.z + MV*OutletV.z;
-        
-        "Energy Holdup"
-        E = ML*OutletL.h + MV*OutletV.h - OutletL.P*V;
-        
-        "Mol fraction normalisation"
-        sum(OutletL.z)= 1.0;
-        
-        "Liquid Volume"
-        vL = PP.LiquidVolume(OutletL.T, OutletL.P, OutletL.z);
-        "Vapour Volume"
-        vV = PP.VapourVolume(OutletV.T, OutletV.P, OutletV.z);
-        
-        "Chemical Equilibrium"
-        PP.LiquidFugacityCoefficient(OutletL.T, OutletL.P, OutletL.z)*OutletL.z = 
-                PP.VapourFugacityCoefficient(OutletV.T, OutletV.P, OutletV.z)*OutletV.z;
-        
-        "Thermal Equilibrium"
-        OutletV.T = OutletL.T;
-        
-        "Mechanical Equilibrium"
-        OutletL.P = OutletV.P;
-        
-        "Geometry Constraint"
-        V*e = ML*vL + MV*vV;
-        
-        "Liquid Density"
-        rhoL = PP.LiquidDensity(OutletL.T, OutletL.P, OutletL.z);
-        "Vapour Density"
-        rhoV = PP.VapourDensity(InletV.T, InletV.P, InletV.z);
-        "Liquid viscosity"
-        miL = PP.LiquidViscosity(OutletL.T, OutletL.P, OutletL.z);
-        "Vapour viscosity"
-        miV = PP.VapourViscosity(InletV.T, InletV.P, InletV.z);
-
-        "Area occupied by the liquid"
-        Al = ML*vL/hs;
-
-        "Volume flow rate of liquid, m^3/m^2/s"
-        uL * Al = OutletL.F * vL;
-        "Volume flow rate of vapor, m^3/m^2/s"
-        uV * (Ap*e - Al) = OutletV.F * vV;
-        
-        "Liquid holdup"
-        hl = ML*vL/V/e;
-        
-        "Particle diameter"
-        dp = 6 * (1-e)/a;
-        
-        "Wall Factor"
-        invK = (1 + (2*dp/(3*ds*(1-e))));
-        
-        "Reynolds number of the vapor stream"
-        Rev*invK = dp*uV*rhoV / (miV*(1-e));
-        
-        deltaP = InletV.P - OutletV.P;
-        
-        "Pressure drop and Vapor flow"
-        deltaP/hs  = Qsil*a*uV^2*rhoV*invK / (2*(e-hl)^3);
-
-        "Liquid holdup"
-        hl = (12*miL*a^2*uL/rhoL/g)^1/3;
-end
-
 Model packedStage_old
         PARAMETERS
@@ -980 +985 @@
         (InletV.P - OutletV.P)/hs  = Qsio*a*uV^2*rhoV*invK / (2*(e-hl)^3);
 end
-
-FlowSheet test
-        PARAMETERS
-        a as Real (Brief="surface area per packing volume", Unit='m^2/m^3');
-#       N as Real (Brief="Number of elements per volume",Unit='1/m^3');
-        niL as Real (Brief="Liquid dynamic viscosity", Unit='kg/m/s');
-        niV as Real (Brief="Vapor dynamic viscosity", Unit='kg/m/s');
-        g as acceleration;
-        rhoL as dens_mass;
-        rhoV as dens_mass;
-        e as Real (Brief="Void fraction of packing, m^3/m^3");
-        V as volume;
-        Across as area;
-        ds as length (Brief="Column diameter");
-        d as length (Brief="size of an element of packing");
-        h as length (Brief="Height of packing");
-        C as Real (Brief="Constant for resitance factor equation");
-        Cp as Real (Brief="Constant for resitance at loading point factor equation");
-        
-        vL as volume_mol;
-        vV as volume_mol;
-        ML as mol;
-        Mw as molweight;
-        dP as pressure (DisplayUnit='atm');
-        
-        VARIABLES
-        hL as Real (Brief="Liquid holdup", Default = 0.424);
-        VL as volume (Brief="Liquid volume", Default = 0.025);
-        uL as velocity (Brief="volume flow rate of liquid, m^3/m^2/s", Default = 0.007);
-        uV as velocity (Brief="volume flow rate of vapor, m^3/m^2/s", Default = 1.14);
-#       n as Real;
-        FV as flow_mol(Default = 149);
-        FL as flow_mol(Default = 222);
-        ksi as Real (Brief="Coefficient of Resistance", Default = 0.784);
-        ksil as Real (Brief="Coefficient of Resistance", Default = 0.032);
-        Rev as Real(Default = 0.966);
-        hLs as Real(Default = 0.037);
-        
-        EQUATIONS
-        VL = vL * ML;
-        hL = VL/V;
-        uL * Across = FL * vL;
-        uV * Across = FV * vV;
-        
-        ksi * C^2 * (uL/uV * sqrt(rhoV/rhoL) * (niL/niV)^5.8)^3 = g/1*'s^2/m';
-        
-        a^2 * niL * uL = hL^1 *(g*rhoL/3 - ksi*a*rhoV*uV^2/(4*hL*(e-hL)^2));
-        
-        dP/h = ksil *(a/2 + 2/ds)*(uV^2*rhoV/(e-hL)^3);
-        
-        ksil = Cp * (exp(uL*rhoL/a/niL/200)*(hL/hLs)^0.3) * (64/Rev+(1.8/Rev)) *
-                                ((e-hL)/e);#1.5
-        
-        Rev = uV * (d-2*hL/a) * rhoV/ niV;
-        
-        hLs = (12*a^2*niL*uL/g/rhoL)^0.333;
-
-        SPECIFY
-#       FV = 147.1 * 'kmol/h';
-#       FL = 229.5 * 'kmol/h';
-#       ksi = 0.809623;
-        
-        SET
-        Mw = 75 * 'g/mol';
-        vL = 9.5e-5 * 'm^3/mol';
-        vV = 0.022 * 'm^3/mol';
-        niL = 0.00032 * 'kg/m/s';
-        niV = 8.2e-5 * 'kg/m/s';
-        rhoL = 809 * 'kg/m^3';
-        rhoV  = 4.63 * 'kg/m^3';
-        ML = 0.268 * 'kmol';
-        dP = 0.1984 * 0.001 * 'atm';#0.1984 * 'atm';
-        
-        V = 0.06 * 'm^3';#0.06 * 'm^3';
-        Across = 0.8 * 'm^2';
-        h = 0.075 * 'm';
-        ds = 1.009 * 'm';
-        d = 50 * 'mm';
-
-#       ksi = 0.8;
-        C = 2.37;
-        Cp = 0.662;
-        e = 0.78;
-        a = 120 * 'm^2/m^3';
-#       N = 6400 * '1/m^3';
-
-end
-
-#*
-        "Reynolds number of the liquid stream"
-        Rel = uL*rhoL / (a*miL);
-
-        if Rel < 5 then
-                if Rel < 1e-4 then
-                        ah = 0 * '1/m';
-                else
-                        "Hydraulic surface area"
-                        ah = a * Ch * Rel^0.15 * (uL^2*a/g)^0.1; #(Rel^2*a^3*miL^2/rhoL^2/g)^0.1; #
-                end
-        else
-                ah = a * Ch * 0.85 * Rel^0.25 * (uL^2*a/g)^0.1; #(Rel^2*a^3*miL^2/rhoL^2/g)^0.1; #
-        end
-        
-
-*#      

branches/packed/eml/streams.mso

@@ -35 +35 @@
 
         VARIABLES
-        F as flow_mol                   (Brief = "Stream Molar Flow Rate", Lower=-1e2);
+        F as flow_mol                   (Brief = "Stream Molar Flow Rate");
         T as temperature                (Brief = "Stream Temperature");
         P as pressure                   (Brief = "Stream Pressure");

branches/packed/sample/stage_separators/sample_column.mso

@@ -557 +557 @@
         sec.stage.Qsil = 0.1;
         #sec.stage(1).OutletV.F = 150 * 'kmol/h';
-        sec.stage.deltaP = 0.0001 * 'atm';
+        sec.stage.deltaP = 0.001 * 'atm';
 
         SET
@@ -563 +563 @@
         sec.NStages = 8;
         sec.stage.Q = 0 * 'kW';
-        sec.stage.Ap = 3.94 * 'ft^2';
-        sec.stage.ds = 1.009 * 'm';
+        sec.stage.d = 1.009 * 'm';
         sec.stage.Cpo = 0.763;
         sec.stage.e = 0.951;
@@ -571 +570 @@
         INITIAL
         sec.stage.OutletL.T =[283:(325-283)/(sec.NStages-1):325] *'K';
-        sec.stage.ML = 0.1 * 'kmol';
-        sec.stage.OutletL.z([1:4]) = [0.2, 0.2, 0.2, 0.2];
+        sec.stage.ML = 0.01 * 'kmol';
+        sec.stage.OutletL.z([1:4]) = [0.3, 0.2, 0.002, 0.1];
 
         OPTIONS
@@ -582 +581 @@
 end
 
-FlowSheet Packed_kettle_cond_Test_1
+FlowSheet Packed_kettle_cond_Test
         PARAMETERS
         PP      as Plugin(Brief="Physical Properties",
                 Type="PP",
-                Components = [ "methanol", "water"],
-                LiquidModel = "IdealLiquid",
-                VapourModel = "Ideal"
+                Components = [ "isobutane", "n-pentane", "propylene",
+                "benzene", "isobutene" ],
+                LiquidModel = "PR",
+                VapourModel = "PR"
         );
         NComp   as Integer;
@@ -605 +605 @@
         
         CONNECTIONS
-        feed.Outlet to col.reb.Inlet;
-        zero to col.stage([1:col.NStages]).Inlet;
+        feed.Outlet to col.stage(5).Inlet;
+        zero to col.reb.Inlet;
+        zero to col.stage([1:4]).Inlet;
+        zero to col.stage([6:col.NStages]).Inlet;
         Qc.OutletQ to col.cond.InletQ;
         Qr.OutletQ to col.reb.InletQ;
         
-        VARIABLES
-        deltaP(col.NStages) as Real (Unit = 'atm/m');
-        
         SPECIFY
-        feed.Outlet.F = 0.98 * 'mol/min';
-        feed.Outlet.T = (50+273.15) * 'K';
-        feed.Outlet.P = 0.91 * 'kPa';
-        feed.Outlet.z = [0.16, 0.84];
+        feed.Outlet.F = 113.4 * 'kmol/h';
+        feed.Outlet.T = 291 * 'K';
+        feed.Outlet.P = 168.3 * 'kPa';
+        feed.Outlet.z = 1/NComp;
         
         zero.F = 0 * 'kmol/h';
@@ -626 +625 @@
         zero.h = 0 * 'J/mol';
         
-        col.sptop.frac = 0.4;
-        col.sptop.Outlet1.F = 0.10 * 'mol/min';
-        col.reb.OutletL.F = 0.88 * 'mol/min';
+        col.sptop.Outlet2.F = 85 * 'kmol/h';
+        col.reb.OutletL.F = 28.4 * 'kmol/h';
+        col.sptop.frac = 0.444445;
         col.cond.OutletV.F = 0 * 'kmol/h';
-        Qr.OutletQ.Q = 700 * 'cal/min';
-        Qc.OutletQ.Q = -700 * 'cal/min';
-        col.pump1.dP = 0.1 * 'atm';
-        
-        EQUATIONS
-        col.reb.OutletV.F = 0.2 * 'mol*min^0.5/kg^0.5/m' * sqrt(Qr.OutletQ.Q);
-        deltaP = (col.stage.InletV.P - col.stage.OutletV.P)/col.stage.hs;
-        
-        SET
-        col.H = 1 * 'm';
-        col.NStages = 2;
-
-        col.cond.V = 1 * 'l';
-        col.cond.Across = 100 * 'cm^2';
-#       col.K = 0.2 * 'mol*min^0.5/kg^0.5/m';
-        
-        col.reb.V = 2 * 'l';
-        col.reb.Across = 200 * 'cm^2';
-
+        Qr.OutletQ.Q = 3.7743e6 * 'kJ/h';
+        Qc.OutletQ.Q = -3.71e6 * 'kJ/h';
+        col.pump1.dP = 16 * 'kPa';
+
+        col.stage.Qsil = 2;
+        col.stage.deltaP = 0.003 * 'atm';
+
+        SET
+        col.H = 3.5 * 'm';
+        col.NStages = 8;
         col.stage.Q = 0 * 'kW';
-        col.stage.Ap = 56.74 * 'cm^2';
-        col.stage.ds = 85 * 'mm';
-        col.stage.e = 0.662;
-        col.stage.a = 185.4 * 'm^2/m^3';
+        col.stage.d = 2.24 * 'ft';
         col.stage.Cpo = 0.763;
-        col.stage.Qsio = 1;
+        col.stage.e = 0.951;
+        col.stage.a = 112.6 * 'm^2/m^3';
+
+        col.cond.V = 2 * 'm^3';
+        col.cond.Across = 1 * 'm^2';
+        col.reb.V = 2 * 'm^3';
+        col.reb.Across = 1 * 'm^2';
         
         INITIAL
         # condenser
-        col.cond.OutletL.T = (63.5+273.15) *'K';
-        col.cond.Level = 2 * 'cm';
-        col.cond.OutletL.z(1) = 0.8; #0.16;
+        col.cond.OutletL.T = 260 *'K';
+        col.cond.Level = 1 * 'm';
+        col.cond.OutletL.z([1:4]) = [0.65, 0.05, 0.01, 0.01];
 
         # reboiler
-        col.reb.OutletL.T = (82+273.15) *'K';
-        col.reb.Level = 2 * 'cm';
-        col.reb.OutletL.z(1) = 0.7;
-
-        # column stages
-        col.stage.OutletL.T = [(62.5+273.15):((83+273.15)-(62.5+273.15))/(col.NStages-1):(83+273.15)] * 'K';
-        #col.stage.Level = 2 * 'cm';
-        col.stage.ML = 0.2 * 'mol';
-        col.stage.OutletL.z(1) = 0.5; #[0.79, 0.65, 0.5, 0.25, 0.16];
+        col.reb.OutletL.T = 330 *'K';
+        col.reb.Level = 1 * 'm';
+        col.reb.OutletL.z([1:4]) = [0.1, 0.7, 0.01, 0.01];
+
+        # column trays
+        col.stage.OutletL.T = [290:(330-290)/(col.NStages-1):330] * 'K';
+        col.stage.ML = 0.1 * 'kmol';
+        col.stage.OutletL.z([1:4]) = [0.15, 0.5, 0.001, 0.1];
 
         OPTIONS
-        DAESolver(File="dassl");
+        TimeStep = 10;
+        TimeEnd = 500;
+end
+
+FlowSheet PackedColumn_ctrl
+        PARAMETERS
+        PP      as Plugin(Brief="Physical Properties",
+                Type="PP",
+                Components = [ "isobutane", "n-pentane", "propylene",
+                        "benzene", "isobutene" ],
+                LiquidModel = "PR",
+                VapourModel = "PR"
+        );
+        NComp   as Integer;
+        
+        Qcmin as heat_rate (Brief="Minimum Condenser Heat supplied");
+        Qcmax as heat_rate (Brief="Maximum Condenser Heat supplied");
+        Qrmin as heat_rate (Brief="Minimum Reboiler Heat supplied");
+        Qrmax as heat_rate (Brief="Maximum Reboiler Heat supplied");
+        Frmin as flow_mol (Brief="Minimum bottom flow rate");
+        Frmax as flow_mol (Brief="Maximum bottom flow rate");
+        Fcmin as flow_mol (Brief="Minimum reflux flow rate");
+        Fcmax as flow_mol (Brief="Maximum reflux flow rate");
+        Hmint as length (Brief="Minimum liquid level in top tank");
+        Hmaxt as length (Brief="Maximum liquid level in top tank");
+    Hminb as length (Brief="Minimum liquid level in reboiler");
+        Hmaxb as length (Brief="Maximum liquid level in reboiler");
+        Pmax as pressure (Brief="Maximum column pressure");
+        Pmin as pressure (Brief="Minimum column pressure");
+        Tmax as temperature (Brief="Maximum column temperature");
+        Tmin as temperature (Brief="Minimum column temperature");
+
+        VARIABLES
+        Qc as energy_source (Brief="Heat rate removed from condenser");
+        Qr as energy_source (Brief="Heat rate supplied to reboiler");
+        Had_top as Real (Brief="Dimensionless condenser level");
+        Had_bot as Real (Brief="Dimensionless reboiler level");
+        Pad as Real (Brief="Dimensionless pressure");
+        Tad as Real (Brief="Dimensionless temperature");
+        RR      as positive (Brief="Reflux ratio");
+
+        SET
+        NComp = PP.NumberOfComponents;
+
+        DEVICES
+        col as PackedDistillation_kettle_cond;
+        feed as source;
+        zero as stream;
+        TCcond as PIDIncr;
+        LCtop as PIDIncr;
+        LCbot as PIDIncr;
+        PC as PIDIncr;
+
+        CONNECTIONS
+        feed.Outlet to col.stage(5).Inlet;
+        zero to col.reb.Inlet;
+        zero to col.stage([1:4]).Inlet;
+        zero to col.stage([6:col.NStages]).Inlet;
+        Qc.OutletQ to col.cond.InletQ;
+        Qr.OutletQ to col.reb.InletQ;
+
+        EQUATIONS
+   "Temperature Controller"
+        TCcond.Parameters.tau = 0*'s';  
+        TCcond.Parameters.tauSet = 0*'s';       
+        TCcond.Parameters.alpha = 0.3;
+        TCcond.Parameters.bias = 0.5;   
+        TCcond.Parameters.gamma = 1;
+        TCcond.Parameters.beta = 1;
+        TCcond.Options.action = 1;
+        TCcond.Options.clip = 1;
+        TCcond.Options.autoMan = 0;
+        TCcond.Parameters.intTime = 60*'s';
+        TCcond.Parameters.gain = 0.6;
+        TCcond.Parameters.derivTime = 1*'s';
+        TCcond.Ports.setPoint = ((15+273.15) * 'K' - Tmin)/(Tmax-Tmin);
+        TCcond.Ports.input = Tad;
+        Tad = (col.cond.OutletL.T-Tmin)/(Tmax-Tmin);
+        Qc.OutletQ.Q = Qcmin+(Qcmax-Qcmin)*TCcond.Ports.output; 
+
+        "Pressure Controller"
+        PC.Parameters.tau = 0*'s';      
+        PC.Parameters.tauSet = 0*'s';   
+        PC.Parameters.alpha = 0.3;
+        PC.Parameters.bias = 0; 
+        PC.Parameters.gamma = 1;
+        PC.Parameters.beta = 1;
+        PC.Options.action = -1;
+        PC.Options.clip = 1;
+        PC.Options.autoMan = 0;
+        PC.Parameters.intTime = 50*'s';
+        PC.Parameters.gain = 0.5;
+        PC.Parameters.derivTime = 1*'s';
+        PC.Ports.setPoint = (4.0*'bar'-Pmin)/(Pmax-Pmin);
+        PC.Ports.input = Pad;
+        Pad = (col.cond.OutletV.P-Pmin)/(Pmax-Pmin);
+        col.cond.OutletV.F = (Fcmin+(Fcmax-Fcmin)*PC.Ports.output);     
+        
+        "Ttop Level Controller"
+        LCtop.Parameters.tau = 0*'s';   
+        LCtop.Parameters.tauSet = 0*'s';        
+        LCtop.Parameters.alpha = 0.3;
+        LCtop.Parameters.bias = 0.5;    
+        LCtop.Parameters.gamma = 1;
+        LCtop.Parameters.beta = 1;
+        LCtop.Options.action = -1;
+        LCtop.Options.clip = 1;
+        LCtop.Options.autoMan = 0;
+        LCtop.Parameters.intTime = 10*'s';
+        LCtop.Parameters.gain = 1;
+        LCtop.Parameters.derivTime = 0*'s';
+        LCtop.Ports.setPoint = (1.0 * 'm' - Hmint)/(Hmaxt-Hmint);
+        LCtop.Ports.input = Had_top;
+        Had_top = (col.cond.Level-Hmint)/(Hmaxt-Hmint);
+        col.sptop.Outlet1.F = Fcmin + (Fcmax-Fcmin) * LCtop.Ports.output;
+
+        "Tbottom Level Controller"
+        LCbot.Parameters.tau = 0*'s';   
+        LCbot.Parameters.tauSet = 0*'s';        
+        LCbot.Parameters.alpha = 0.3;
+        LCbot.Parameters.bias = 0.5;    
+        LCbot.Parameters.gamma = 1;
+        LCbot.Parameters.beta = 1;
+        LCbot.Options.action = -1;
+        LCbot.Options.clip = 1;
+        LCbot.Options.autoMan = 0;
+        LCbot.Parameters.intTime = 100*'s';
+        LCbot.Parameters.gain = 1;
+        LCbot.Parameters.derivTime = 0*'s';
+        LCbot.Ports.setPoint = (1.0 * 'm' - Hminb)/(Hmaxb-Hminb);
+        LCbot.Ports.input = Had_bot;
+        Had_bot = (col.reb.Level-Hminb)/(Hmaxb-Hminb);
+        col.reb.OutletL.F = Frmin + (Frmax-Frmin) * LCbot.Ports.output;
+
+        RR * (col.cond.OutletV.F + col.sptop.Outlet1.F) = col.sptop.Outlet2.F;
+        
+        if time < 1 * 'h' then
+                col.sptop.Outlet2.F = 75 * 'kmol/h'; # reflux
+        else
+                col.sptop.Outlet2.F = 85 * 'kmol/h'; # reflux
+        end
+
+        SPECIFY
+        feed.Outlet.F = 113.4 * 'kmol/h';
+        feed.Outlet.T = 291 * 'K';
+        feed.Outlet.P = 5 * 'bar';
+        feed.Outlet.z = 1/NComp;
+        
+        zero.F = 0 * 'kmol/h';
+        zero.T = 300 * 'K';
+        zero.P = 1 * 'atm';
+        zero.z = 1/NComp;
+        zero.v = 0;
+        zero.h = 0 * 'J/mol';
+        
+        Qr.OutletQ.Q = 4e6 * 'kJ/h';
+        col.pump1.dP = 16 * 'kPa';
+        
+        col.stage.Qsil = 2;
+        col.stage.deltaP = 0.003 * 'atm';
+
+        SET
+        col.H = 3.5 * 'm';
+        col.NStages = 8;
+        col.stage.Q = 0 * 'kW';
+        col.stage.d = 2.24 * 'ft';
+        col.stage.Cpo = 0.763;
+        col.stage.e = 0.951;
+        col.stage.a = 112.6 * 'm^2/m^3';
+
+        col.cond.V = 2 * 'm^3';
+        col.cond.Across = 1 * 'm^2';
+        col.reb.V = 2 * 'm^3';
+        col.reb.Across = 1 * 'm^2';
+
+        # Controllers type
+        TCcond.PID_Select = "Ideal_AW";
+        PC.PID_Select = "Ideal_AW";
+        LCtop.PID_Select = "Ideal_AW";
+        LCbot.PID_Select = "Ideal_AW";
+
+        Qrmax = 5e6 * 'kJ/h';
+        Qrmin = 1e6 * 'kJ/h';
+        Frmin = 0 * 'kmol/h';
+        Frmax = 60 * 'kmol/h';
+        Fcmin = 0 * 'kmol/h';
+        Fcmax = 120 * 'kmol/h';
+        Hmint = 0 * 'm';
+        Hmaxt = 2 * 'm';
+        Hminb = 0 * 'm';
+        Hmaxb = 2 * 'm';
+        Pmin = 0.5 * 'bar';
+        Pmax = 6 * 'bar';
+        Qcmax = -5e5 * 'kJ/h';
+        Qcmin = -5e6 * 'kJ/h';
+        Tmax = (30+273.15) * 'K';
+        Tmin = (-20+273.15) * 'K';
+        
+        INITIAL
+        # condenser
+        col.cond.OutletL.T = 260 *'K';
+        col.cond.Level = 1 * 'm';
+        col.cond.OutletL.z([1:4]) = [0.2, 0.2, 0.4, 0.05];
+
+        # reboiler
+        col.reb.OutletL.T = 350 *'K';
+        col.reb.Level = 1 * 'm';
+        col.reb.OutletL.z([1:4]) = [0.1, 0.4, 0.1, 0.3];
+
+        # column trays
+        col.stage.OutletL.T = [290:(330-290)/(col.NStages-1):330] * 'K';
+        col.stage.ML = 0.1 * 'kmol';
+        col.stage.OutletL.z([1:4]) = [0.15, 0.3, 0.25, 0.2];
+
+        OPTIONS
         TimeStep = 0.1;
-        TimeEnd = 50;
+        TimeEnd = 5;
+        TimeUnit = 'h';
+        #InitialFile = "Column_ctrl.rlt";
 end

branches/packed/sample/stage_separators/sample_tray.mso

@@ -369 +369 @@
         inV.z = [0.0584, 0.9416];#[0.5, 0.5];#
         
-        #t1.OutletV.P = 145 * 'kPa';
-        t1.OutletV.F = 190 * 'kmol/h';
+        t1.OutletV.P = 145 * 'kPa';
+        #t1.OutletV.F = 190 * 'kmol/h';
         t1.Qsil = 10;
         
@@ -445 +445 @@
         inV.z = [0.265, 0.233, 0.150, 0.014, 0.338];
         
-        #t1.deltaP = 0.01 * 'atm';
-        t1.OutletV.F = 165 * 'kmol/h';
+        t1.deltaP = 0.01 * 'atm';
+        #t1.OutletV.F = 165 * 'kmol/h';
         t1.Qsil = 10;