Changeset 243


Ignore:
Timestamp:
Apr 16, 2007, 2:19:39 PM (16 years ago)
Author:
Paula Bettio Staudt
Message:

Updated reactive distillation models andsample

Location:
trunk
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/eml/controllers/PIDs.mso

    r176 r243  
    9797        if (Parameters.tau equal 0) then
    9898                "Input first order filter"
    99                 (Parameters.tau + 1e-3*"s")*diff(Internal.inputFilt)= Ports.input - Internal.inputFilt;
     99                (Parameters.tau + 1e-3*'s')*diff(Internal.inputFilt)= Ports.input - Internal.inputFilt;
    100100        else
    101101                "Input first order filter"
    … …  
    105105        if (Parameters.tauSet equal 0) then
    106106                "setPoint first order filter"
    107                 (Parameters.tauSet + 1e-3*"s")*diff(Internal.setPointFilt)= Ports.setPoint - Internal.setPointFilt;
     107                (Parameters.tauSet + 1e-3*'s')*diff(Internal.setPointFilt)= Ports.setPoint - Internal.setPointFilt;
    108108        else
    109109                "setPoint first order filter"
    … …  
    132132        if (Parameters.derivTime equal 0) then
    133133                "Derivative term filter"       
    134                 Parameters.alpha*(Parameters.derivTime + 1e-3*"s")*diff(Internal.dFilt) = Internal.errorD - Internal.dFilt;
     134                Parameters.alpha*(Parameters.derivTime + 1e-3*'s')*diff(Internal.dFilt) = Internal.errorD - Internal.dFilt;
    135135        else
    136136                "Derivative term filter"       

  • trunk/eml/stage_separators/column.mso

    r210 r243  
    930930Model ReactiveDistillation
    931931        PARAMETERS
    932         outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
    933         outer NComp as Integer;
    934         NTrays as Integer(Brief="Number of trays", Default=2);
     932        outer PP as Plugin(Type="PP");
     933        outer NComp as Integer;
     934        NTrays as Integer(Brief="Number of trays", Default=2);
     935        topdown as Integer(Brief="Trays counting (1=top-down, -1=bottom-up)", Default=1);
     936        top as Integer(Brief="Number of top tray");
     937        bot as Integer(Brief="Number of bottom tray");
     938        alfacond as Real;
     939
     940        VapourFlow as Switcher(Valid = ["on", "off"], Default = "off");
     941       
     942        SET
     943        top = (NTrays-1)*(1-topdown)/2+1;
     944        bot = NTrays/top;
    935945       
    936946        VARIABLES
    … …  
    942952       
    943953        EQUATIONS
    944         if ( reb.OutletV.P > 1 * 'atm' ) then
    945                 "Pressure Drop through the tray"
    946                 reb.OutletV.F = trays(1).Ah/reb.vV * sqrt((reb.OutletV.P - 1*'atm') / (0.15*reb.rhoV) );
    947         else
     954       
     955        switch VapourFlow
     956                case "on":
     957                "Pressure Drop through the condenser"
     958                cond.InletV.F*trays(top).vV / 'm^2' =
     959                        sqrt((trays(top).OutletV.P - cond.OutletL.P + 1e-8 * 'atm')/(trays(top).rhoV*alfacond));
     960                when trays(top).OutletV.P < cond.OutletL.P switchto "off";
     961               
     962                case "off":
    948963                "Prato selado"
    949                 reb.OutletV.F = 0.0 * 'mol/s';
     964                cond.InletV.F = 0.0 * 'mol/s';
     965                when trays(top).OutletV.P > cond.OutletL.P + 1e-3 * 'atm' switchto "on";
    950966        end
    951        
    952         CONNECTIONS
    953         #vapor
    954         reb.OutletV to trays([NTrays]).InletV;
    955         trays([2:NTrays]).OutletV to trays([1:NTrays-1]).InletV;
    956         trays(1).OutletV to cond.InletV;
    957 
     967
     968        CONNECTIONS
     969        #vapor
     970        reb.OutletV to trays(bot).InletV;
     971        trays([top+topdown:topdown:bot]).OutletV to trays([top:topdown:bot-topdown]).InletV;
     972        trays(top).OutletV to cond.InletV;
     973       
    958974        #liquid
    959975        cond.OutletL to sp.Inlet;       
    960976        sp.Outlet2 to p.Inlet;
    961         p.Outlet to trays(1).InletL;
    962         trays([1:NTrays-1]).OutletL to trays([2:NTrays]).InletL;
    963         trays(NTrays).OutletL to reb.InletL;
    964 end
     977        p.Outlet to trays(top).InletL;
     978        trays([top:topdown:bot-topdown]).OutletL to trays([top+topdown:topdown:bot]).InletL;
     979        trays(bot).OutletL to reb.InletL;
     980       
     981end

  • trunk/eml/stage_separators/condenser.mso

    r210 r243  
    135135Model condenserReact
    136136        PARAMETERS
    137         outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
     137        outer PP as Plugin(Type="PP");
    138138        outer NComp as Integer;
    139139        V as volume (Brief="Condenser total volume");
    … …  
    158158        Q as heat_rate (Brief="Heat supplied");
    159159        Vol as volume;
    160         r as reaction_mol (Brief = "Specific reaction rate");
     160        r3 as reaction_mol (Brief = "Reaction resulting ethyl acetate", DisplayUnit = 'mol/l/s');
    161161        C(NComp) as conc_mol (Brief = "Molar concentration", Lower = -1);
    162162
    … …  
    165165        OutletL.z = vL * C;
    166166       
     167        "Reaction"
     168        r3 = exp(-7150*'K'/OutletL.T)*(4.85e4*C(1)*C(2) - 1.23e4*C(3)*C(4)) * 'l/mol/s';
     169       
    167170        "Component Molar Balance"
    168171        diff(M) = InletV.F*InletV.z - OutletL.F*OutletL.z
    169                                 - OutletV.F*OutletV.z + stoic*r*ML*vL;
     172                                - OutletV.F*OutletV.z + stoic*r3*ML*vL;
    170173
    171174        "Energy Balance"
    172175        diff(E) = InletV.F*InletV.h - OutletL.F*OutletL.h
    173                                 - OutletV.F*OutletV.h + Q + Hr * r * ML*vL;
     176                                - OutletV.F*OutletV.h + Q + Hr * r3 * ML*vL;
    174177
    175178        "Molar Holdup"
    … …  
    206209
    207210        sum(OutletL.z)=sum(OutletV.z);
     211
    208212end

  • trunk/eml/stage_separators/reboiler.mso

    r210 r243  
    191191Model reboilerReact
    192192        PARAMETERS
    193         outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
     193        outer PP as Plugin(Type="PP");
    194194        outer NComp as Integer;
    195195        Across as area (Brief="Cross Section Area of reboiler");
    … …  
    217217        startup as Real;
    218218        rhoV as dens_mass;
    219         r as reaction_mol (Brief = "Specific reaction rate");
     219        r3 as reaction_mol (Brief = "Reaction resulting ethyl acetate", DisplayUnit = 'mol/l/s');
    220220        C(NComp) as conc_mol (Brief = "Molar concentration", Lower = -1);
    221221
    … …  
    224224        OutletL.z = vL * C;
    225225       
     226        "Reaction"
     227        r3 = exp(-7150*'K'/OutletL.T)*(4.85e4*C(1)*C(2) - 1.23e4*C(3)*C(4)) * 'l/mol/s';
     228
    226229        "Component Molar Balance"
    227230        diff(M)= Inlet.F*Inlet.z + InletL.F*InletL.z
    228                 - OutletL.F*OutletL.z - OutletV.F*OutletV.z + stoic*r*ML*vL;
     231                - OutletL.F*OutletL.z - OutletV.F*OutletV.z + stoic*r3*ML*vL;
    229232       
    230233        "Energy Balance"
    231234        diff(E) = Inlet.F*Inlet.h + InletL.F*InletL.h
    232                 - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q + Hr * r * vL*ML;
     235                - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q + Hr * r3 * vL*ML;
    233236       
    234237        "Molar Holdup"
    … …  
    260263       
    261264        "Geometry Constraint"
    262         V = ML*vL + MV*vV;
     265        V = ML*vL + MV*vV;             
    263266
    264267        "Chemical Equilibrium"

  • trunk/eml/stage_separators/tray.mso

    r210 r243  
    159159
    160160        PARAMETERS
    161         outer PP as Plugin(Brief = "External Physical Properties", Type="PP");
     161        outer PP as Plugin(Type="PP");
    162162        outer NComp as Integer;
    163163        V as volume(Brief="Total Volume of the tray");
    … …  
    175175        Hr as energy_mol;
    176176        Pstartup as pressure;
     177       
     178        VapourFlow as Switcher(Valid = ["on", "off"], Default = "off");
     179        LiquidFlow as Switcher(Valid = ["on", "off"], Default = "off");
    177180       
    178181        VARIABLES
    … …  
    197200        rhoL as dens_mass;
    198201        rhoV as dens_mass;
    199         r as reaction_mol (Brief = "Specific reaction rate");
     202        r3 as reaction_mol (Brief = "Reaction resulting ethyl acetate", DisplayUnit = 'mol/l/s');
    200203        C(NComp) as conc_mol (Brief = "Molar concentration", Lower = -1); #, Unit = "mol/l");
    201204       
    … …  
    204207        OutletL.z = vL * C;
    205208       
     209        "Reaction"
     210        r3 = exp(-7150*'K'/OutletL.T)*(4.85e4*C(1)*C(2) - 1.23e4*C(3)*C(4))*'l/mol/s';
     211       
    206212        "Component Molar Balance"
    207213        diff(M)=Inlet.F*Inlet.z + InletL.F*InletL.z + InletV.F*InletV.z
    208                 - OutletL.F*OutletL.z - OutletV.F*OutletV.z + stoic*r*ML*vL;
     214                - OutletL.F*OutletL.z - OutletV.F*OutletV.z + stoic*r3*ML*vL;
    209215       
    210216        "Energy Balance"
    211217        diff(E) = ( Inlet.F*Inlet.h + InletL.F*InletL.h + InletV.F*InletV.h
    212                 - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q ) + Hr * r * vL*ML;
     218                - OutletL.F*OutletL.h - OutletV.F*OutletV.h + Q ) + Hr * r3 * vL*ML;
    213219       
    214220        "Molar Holdup"
    … …  
    242248        rhoV = PP.VapourDensity(InletV.T, InletV.P, InletV.z);
    243249
    244         if Level > (beta * hw) then
     250        switch LiquidFlow
     251                case "on":
    245252                "Francis Equation"
    246                 OutletL.F = (1.84/'s'*lw*((Level-(beta*hw))/(beta))^2/vL);
    247         else
     253                OutletL.F*vL = 1.84*'1/s'*lw*((Level-(beta*hw)+1e-6*'m')/(beta))^2;
     254                when Level < (beta * hw) switchto "off";
     255               
     256                case "off":
    248257                "Low level"
    249258                OutletL.F = 0 * 'mol/h';
     259                when Level > (beta * hw) + 1e-6*'m' switchto "on";
    250260        end
    251261
     262        switch VapourFlow
     263                case "on":
     264                #InletV.P = OutletV.P + Level*g*rhoL + rhoV*alfa*(InletV.F*vV/Ah)^2;
     265                InletV.F*vV = sqrt((InletV.P - OutletV.P - Level*g*rhoL + 1e-8 * 'atm')/(rhoV*alfa))*Ah;
     266                when InletV.P < OutletV.P + Level*g*rhoL switchto "off";
    252267               
    253         "Pressure Drop through the tray"
    254         OutletV.F = (1 + tanh(1 * (OutletV.P - Pstartup)/'Pa'))/2 *
    255                 Ah/vV * sqrt(2*(OutletV.P - InletL.P + 1e-8 * 'atm') / (alfa*rhoV) );
     268                case "off":
     269                InletV.F = 0 * 'mol/s';
     270                when InletV.P > OutletV.P + Level*g*rhoL + 3e-2 * 'atm' switchto "on";
     271                #when InletV.P > OutletV.P + Level*beta*g*rhoL + 1e-2 * 'atm' switchto "on";
     272        end
    256273
    257274        "Chemical Equilibrium"

  • trunk/sample/stage_separators/sample_columnReact.mso

    r213 r243  
    2828FlowSheet Startup_ReactiveDistillation
    2929        PARAMETERS
    30         PP              as Plugin(Brief="Physical Properties",
    31                 Type="PP",
     30        PP      as Plugin(Brief="Physical Properties", Type="PP",
    3231                Components = [ "acetic acid", "ethanol",  "ethyl acetate", "water"],
    3332                LiquidModel = "UNIFAC",
    34                 VapourModel = "Ideal",
    35                 Derivatives = 1
     33                VapourModel = "Ideal"
    3634        );
    3735        NComp   as Integer;
    … …  
    5048
    5149        VARIABLES
     50        Qcmin as heat_rate (Brief="Condenser Heat supplied");
     51        Qcmax as heat_rate (Brief="Condenser Heat supplied");
     52        Qrmin as heat_rate (Brief="Reboiler Heat supplied");
     53        Qrmax as heat_rate (Brief="Reboiler Heat supplied");
     54        Fmin as flow_mol;
     55        Fmax as flow_mol;
     56        Frmin as flow_mol;
     57        Frmax as flow_mol;
     58       
    5259        Lreb_ad as Real;
    5360        Lrebmin as length;
    … …  
    6976        zero to col.trays([6:col.NTrays]).Inlet;
    7077       
     78       
    7179        EQUATIONS
    72         "Equações do PID para controle de nível"
    73         Lreb_ad*(Lrebmax-Lrebmin)=col.reb.Level-Lrebmin;
    74         PIDLreb.Ports.input=Lreb_ad;
    75         Lcond_ad*(Lcondmax-Lcondmin)=col.cond.Level-Lcondmin;
    76         PIDLcond.Ports.input=Lcond_ad;
    77         "Equações do PID para controle de temperatura"
    78         Treb_ad*(Trebmax-Trebmin)=col.reb.OutletL.T-Trebmin;
    79         PIDTreb.Ports.input=Treb_ad;
    80         Tcond_ad*(Tcondmax-Tcondmin)=col.cond.OutletL.T-Tcondmin;
    81         PIDTcond.Ports.input=Tcond_ad;
    82        
    83         col.sp.frac = 0.1;
    84 
    85         col.sp.Outlet1.F = 0.5 * PIDLcond.Ports.output * 'mol/s';
    86         col.reb.OutletL.F = 1.736 * PIDLreb.Ports.output * 'mol/s';
    87 
    88         PIDTreb.Ports.setPoint= (367 * 'K' - Trebmin)/(Trebmax-Trebmin);
    89         PIDTcond.Ports.setPoint= (344 * 'K' - Trebmin)/(Trebmax-Trebmin);
    90         col.cond.Q = 10*'J/s' - (5 * PIDTcond.Ports.output) * 'kJ/s';
    91 
    92 #verificando a partida do refervedor
    93         if time < 200 * 's' then
     80        col.sp.frac = 0.09;
     81
     82        #verificando a partida do refervedor
     83        if time < 400 * 's' then
    9484                col.reb.startup = 1;
    9585        else
    9686                col.reb.startup = 0;
    9787        end
    98 
     88 
    9989        if col.reb.startup then
     90                col.cond.Q = 0 * PIDTcond.Ports.output * 'kJ/s';
    10091                col.reb.Q = 0 * PIDTreb.Ports.output * 'kJ/s';
     92               
     93                PIDTreb.Ports.input = PIDTreb.Ports.setPoint;
    10194        else
    102                 col.reb.Q = 1e1 * PIDTreb.Ports.output * 'kJ/s';
     95                col.cond.Q = Qcmin+(Qcmax-Qcmin)*PIDTcond.Ports.output;
     96                col.reb.Q = Qrmin+(Qrmax-Qrmin)*PIDTreb.Ports.output;
     97               
     98                PIDTreb.Ports.input=Treb_ad;
    10399        end
    104        
    105         "Reaction - Trays"
    106         col.trays.r = exp(-7150*'K'/col.trays.OutletL.T)*
    107         (4.85e4*col.trays.C(1)*col.trays.C(2) - 1.23e4*col.trays.C(3)*col.trays.C(4))*'l/mol/s';
    108         "Reaction - Reboiler"
    109         col.reb.r = exp(-7150*'K'/col.reb.OutletL.T)*(4.85e4*col.reb.C(1)*col.reb.C(2)
    110                         - 1.23e4*col.reb.C(3)*col.reb.C(4)) * 'l/mol/s';
    111         "Reaction - Condenser"
    112         col.cond.r = exp(-7150*'K'/col.cond.OutletL.T)*(4.85e4*col.cond.C(1)*col.cond.C(2)
    113                         - 1.23e4*col.cond.C(3)*col.cond.C(4)) * 'l/mol/s';
    114 
    115100
    116101        SPECIFY
    … …  
    138123        PIDLreb.Parameters.beta=1;
    139124        PIDLreb.Parameters.gain=1;
    140         PIDLreb.Parameters.intTime=100*'s';
     125        PIDLreb.Parameters.intTime=10*'s';
    141126        PIDLreb.Parameters.derivTime=1*'s';
    142127        PIDLreb.Options.action=-1;
    … …  
    144129        PIDLreb.Options.autoMan=0;
    145130        PIDLreb.Ports.setPoint=(0.5 * 'm' - Lrebmin)/(Lrebmax-Lrebmin);
     131        Lreb_ad*(Lrebmax-Lrebmin)=col.reb.Level-Lrebmin;
     132        PIDLreb.Ports.input=Lreb_ad;
     133        col.reb.OutletL.F = Frmin + (Frmax-Frmin) * PIDLreb.Ports.output;
    146134
    147135        PIDLcond.Parameters.tau = 1*'s';       
    148         PIDLcond.Parameters.tauSet=1*'s';       
    149         PIDLcond.Parameters.bias = 0;
     136        PIDLcond.Parameters.tauSet=1*'s';
     137        PIDLcond.Parameters.bias = 0.5;
    150138        PIDLcond.Parameters.alpha=1;
    151139        PIDLcond.Parameters.gamma=1;
    … …  
    158146        PIDLcond.Options.autoMan=0;
    159147        PIDLcond.Ports.setPoint=(0.5 * 'm' - Lcondmin)/(Lcondmax-Lcondmin);
     148        Lcond_ad*(Lcondmax-Lcondmin)=col.cond.Level-Lcondmin;
     149        PIDLcond.Ports.input=Lcond_ad;
     150        col.sp.Outlet1.F = Fmin + (Fmax-Fmin) * PIDLcond.Ports.output;
    160151
    161152        PIDTreb.Parameters.tau = 1*'s';
    162         PIDTreb.Parameters.tauSet=1*'s';       
     153        PIDTreb.Parameters.tauSet=1*'s';
    163154        PIDTreb.Parameters.bias = 0.2;
    164155        PIDTreb.Parameters.alpha=0.2;
    … …  
    166157        PIDTreb.Parameters.beta=1;
    167158        PIDTreb.Parameters.gain=0.9;
    168         PIDTreb.Parameters.intTime=10*'s';
     159        PIDTreb.Parameters.intTime=100*'s';
    169160        PIDTreb.Parameters.derivTime=1*'s';
    170161        PIDTreb.Options.action=1;
    171162        PIDTreb.Options.clip=1;
    172163        PIDTreb.Options.autoMan=0;
    173 
    174         PIDTcond.Parameters.tau = 1*'s';
    175         PIDTcond.Parameters.tauSet=1*'s';       
    176         PIDTcond.Parameters.bias = 0.2;
     164        PIDTreb.Ports.setPoint= (366 * 'K' - Trebmin)/(Trebmax-Trebmin);
     165        Treb_ad*(Trebmax-Trebmin)=col.reb.OutletL.T-Trebmin;
     166
     167        PIDTcond.Parameters.tau = 1*'s';       
     168        PIDTcond.Parameters.tauSet=1*'s';
     169        PIDTcond.Parameters.bias = 0.5;
    177170        PIDTcond.Parameters.alpha=0.2;
    178171        PIDTcond.Parameters.gamma=1;
    … …  
    181174        PIDTcond.Parameters.intTime=10*'s';
    182175        PIDTcond.Parameters.derivTime=1*'s';
    183         PIDTcond.Options.action=-1;
     176        PIDTcond.Options.action=1;
    184177        PIDTcond.Options.clip=1;
    185178        PIDTcond.Options.autoMan=0;
     179        PIDTcond.Ports.setPoint= (346 * 'K' - Tcondmin)/(Tcondmax-Tcondmin);
     180        Tcond_ad*(Tcondmax-Tcondmin)=col.cond.OutletL.T-Tcondmin;
     181        PIDTcond.Ports.input=Tcond_ad;
    186182       
    187183        "Valores limites para normalizações"
    … …  
    194190        Tcondmax=380*'K';
    195191        Tcondmin=250*'K';
     192        Qcmin = -100 * 'kJ/s';
     193        Qcmax = 0 * 'kJ/s';
     194        Qrmin = 0 * 'kJ/s';
     195        Qrmax = 150 * 'kJ/s';
     196        Fmin = 0 * 'kmol/h';
     197        Fmax = 2 * 'kmol/h';
     198        Frmin = 0 * 'kmol/h';
     199        Frmax = 5 * 'kmol/h';
    196200       
    197201        col.cond.OutletV.F = 0 * 'kmol/h';
    … …  
    206210        col.cond.Across = 6 * 'l/m';
    207211       
    208         col.trays.V =  0.0961 * 'm^3'; # 0.34* (3.14*0.6*0.6/4)
    209         col.trays.Ah = 0.04 * 'm^2';#0.2 * 'm^2';
     212        col.trays.V =  0.0961 * 'm^3';
     213        col.trays.Ah = 0.04 * 'm^2';
    210214        col.trays.lw = 0.457 * 'm';
    211215        col.trays.hw = 0.05 * 'm';
    212216        col.trays.Q = 0 * 'kW';
    213217        col.trays.beta = 0.8;
    214         col.trays.alfa = 5;
    215         col.trays.Ap = 0.07 * 'm^2'; #0.24 * 'm^2'; # 3.14*0.6*0.6/4 - 15%
     218        col.trays.alfa = 30;
     219        col.alfacond = 100000;
     220        col.trays.Ap = 0.07 * 'm^2';
    216221
    217222        col.trays.Hr = 0 * 'kJ/mol';
    … …  
    232237
    233238        # reboiler
    234         col.reb.OutletL.T = 300 *'K';
     239        col.reb.OutletL.T = 300 * 'K';
    235240        col.reb.Level = 0.1 * 'm';
    236241        col.reb.OutletL.z([1:3]) = [0.4962, 0.4808, 0];
    … …  
    242247
    243248        OPTIONS
     249        TimeStep = 100;
     250        TimeEnd = 50000;
     251        NLASolver = "sundials";
    244252        DAESolver = "dassl";
    245         RelativeAccuracy = 1e-2;
    246         TimeEnd = 20000;
    247         TimeStep = 100;
    248         TimeUnit = 's';
    249 #       time = [0:200, 210:10:2000, 2100:100:20000];
    250253end
Note: See TracChangeset for help on using the changeset viewer.