source: trunk/sample/controllers/CSTR_noniso_pid.mso @ 1009

#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2008 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 - 2008 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.
*
*--------------------------------------------------------------------
* Sample file for controllers
*----------------------------------------------------------------------
* Author: Argimiro R. Secchi
* $Id: CSTR_noniso_pid.mso 536 2008-06-15 23:18:35Z arge $
*--------------------------------------------------------------------*#

using "controllers/PIDs";

const_valv as positive(Brief = "Valve Constant", Default=1,Lower=0,Upper=100, Unit='m^2.5/h');

Model stream_cstr
        VARIABLES
        Ca     as conc_mol;
        F      as flow_vol;
        T      as temperature;
end

Model CSTR

        PARAMETERS
        ko  as frequency                (DisplayUnit='1/h');
        D       as length;
        A   as area;
        Ea  as energy_mol               (DisplayUnit='kJ/kmol');
        R   as Real                     (Unit='kJ/mol/K');
        ro  as dens_mass                (DisplayUnit='kg/m^3');
        Cp  as cp_mass                  (DisplayUnit='kJ/kg/K');
        U   as heat_trans_coeff (DisplayUnit='kW/m^2/K');
        Hr  as heat_reaction    (DisplayUnit='kJ/kmol');
        pi  as Real                             (Default = 3.141593);
        Cv  as const_valv;

        VARIABLES
        At               as area;       
        T        as temperature;
        Tw       as temperature;
        x                as fraction;
        V        as volume;
        Ca       as conc_mol;
        h        as length;
        tau      as time_h;
        rA           as reaction_mol;
        k        as frequency   (DisplayUnit='1/h');
        q                as heat_rate   (DisplayUnit='kJ/h');
        qr               as heat_rate   (DisplayUnit='kJ/h');
in  Inlet    as stream_cstr;
out Outlet   as stream_cstr;

        SET
        A = pi * D^2 / 4;
        
        EQUATIONS

        "Overall Mass Balance"
        diff(V) = Inlet.F - Outlet.F;
        
        "Component Mass Balance"
        V * diff(Ca) = Inlet.F * (Inlet.Ca - Ca) - (-rA) * V;

        "Average Residence Time"
        tau * Inlet.F = V;

        "Energy Balance"
        ro * V * Cp * diff(T) = Inlet.F * ro * Cp * (Inlet.T - T) + qr - q;

        "Heat Transfer Rate"
        q = U * At * (T - Tw);

        "Reaction Heat Rate"
        qr = (-Hr) * (-rA) * V;
        
        "Reaction Rate"
        -rA = k * Ca;
        
        "Arrhenius Equation"
        k = ko * exp(-Ea/(R*T));
        
        "Geometry"
        A * h = V;
        At = A + pi*D*h;
        
        "Valve Equation"
        Outlet.F = x * Cv * sqrt(h);

        "Perfect Mixture"
        Outlet.Ca = Ca;
        Outlet.T  = T;
end

# Process with controlled CSTR and multiple steady-states
FlowSheet CSTR_controller
        PARAMETERS
        Lmin as length;
        Lmax as length;
        Tmin as temperature;
        Tmax as temperature;

        DEVICES
        FEED as stream_cstr;
        CSTR1 as CSTR;
        PIDL as PID;
        PIDT as PID;
 
        VARIABLES
        Lsp as length;
        Tsp as temperature;
out     TI      as control_signal;
out     LI      as control_signal;
        
        CONNECTIONS
        FEED to CSTR1.Inlet;
        LI       to PIDL.Input;
        TI       to PIDT.Input;

        SET
#       CSTR Parameters
        CSTR1.R   = 8.3144 * 'kJ/kmol/K';
        CSTR1.U   = 300 * 'kJ/h/m^2/K';
        CSTR1.ro  = 1000 * 'kg/m^3';
        CSTR1.Cp  = 4*'kJ/kg/K';
        CSTR1.Hr  = -7000 * 'kJ/kmol';
        CSTR1.Ea  = 6e4 * 'kJ/kmol';
        CSTR1.ko  = 89 * '1/s';
        CSTR1.D   = 3.2 * 'm';
        CSTR1.Cv  = 2.7 * 'm^2.5/h';

#       Operating range for control variables
        Lmax=5*'m';
        Lmin=0*'m';
        Tmax=700*'K';
        Tmin=230*'K';

        PIDL.PID_Select = "Ideal_AWBT";
        PIDT.PID_Select = "Ideal_AWBT";

#   Level control: PID parameters
        PIDL.bias=0;
        PIDL.alpha=0.1;
        PIDL.Action="Direct";
        PIDL.gamma=1;
        PIDL.beta=1;
        PIDL.Clip="Clipped";
        PIDL.Mode="Automatic";
        PIDL.gain=1;
        PIDL.intTime=2.5*'h';
        PIDL.derivTime=0*'s';
        PIDL.tau=1*'s';
        PIDL.tauSet=1*'s';
        PIDL.MinInput=Lmin/'m';
        PIDL.MaxInput=Lmax/'m';

#   Temperature control: PID parameters
        PIDT.bias = 0;
        PIDT.alpha=0.1;
        PIDT.Action="Reverse";
        PIDT.gamma=1;
        PIDT.beta=1;
        PIDT.Clip="Clipped";
        PIDT.Mode="Automatic";
        PIDT.gain=1;
#       PIDT.gain=15; # for a setpoint of 523K or a faster response
        PIDT.intTime=2.5*'h';
        PIDT.derivTime=1*'h';
        PIDT.tau=1*'s';
        PIDT.tauSet=1*'s';
        PIDT.MinInput=Tmin/'K';
        PIDT.MaxInput=Tmax/'K';
        PIDT.MinOutput=Tmin/'K';
        PIDT.MaxOutput=Tmax/'K';
        
        EQUATIONS

        "Dimensionless level to connect PID"
        LI=CSTR1.h/'m';

        "Dimensionless temperature to connect PID"
        TI=CSTR1.T/'K';

        "Manipulated Variables"
        CSTR1.x  = PIDL.Output;
        CSTR1.Tw  = PIDT.Output*'K';

        "Level setpoint"
        PIDL.SetPoint=Lsp/'m';

        "Temperature setpoint"
        PIDT.SetPoint=Tsp/'K';

        "Feed Stream"
        FEED.Ca = 300 * 'kmol/m^3';
        FEED.F = 3.5 * 'm^3/h';

#   Disturbance
        if time < 50 * 'h' then 
                FEED.T = 300 * 'K';
        else
                FEED.T = 285 * 'K'; # change to 350 K to saturate controller
        end

#   Set-point changes
        if time < 100 * 'h' then
          Tsp = 630 * 'K';
        else
          Tsp = 400 * 'K';
#         Tsp = 523 * 'K'; # near unstable steady-state (change also the controller's gain)
        end

        if time < 150 * 'h' then
          Lsp = 1.7 * 'm';
        else
          Lsp = 4 * 'm';
        end

        INITIAL
        CSTR1.Ca = 50 * 'kmol/m^3';
        CSTR1.h = 1.7 * 'm';
        CSTR1.T = 570 * 'K';
        
        OPTIONS
        TimeStep = 1;
        TimeEnd = 250;
        TimeUnit = 'h';
        DAESolver(File = "dassl");
end