source: branches/rate/sample/estimation/BatchReactor.mso @ 502

using "types";

#*
* Batch Reactor model based on gPROMS sample
*#
Model BatchReactor
        PARAMETERS
        Density(4) as dens_mol(Brief="Component Densities");
        Ac(4) as Real(Brief="Specific enthalpy coeficient A", Unit='J/mol/K');
        Bc(4) as Real(Brief="Specific enthalpy coeficient B", Unit='J/mol/K/K');
        
        k0 as Real(Brief="Arrenius constant", Unit='m^3/s/mol');
        Ea as Real(Brief="Activation energy", Unit='J/mol');
        Hr as Real(Brief="Reaction enthalpy", Unit='J/mol');
        
        Nu(4) as Real(Brief="Stoichiometric coefficients");
        R as Real(Brief="Ideal gas constant", Default=8.314, Unit='J/K/mol');
        Tref as temperature(Brief="Reference temperature", Default=298);

        VARIABLES
        Holdup(4) as mol;
        C(4) as conc_mol;
        TotalVolume as volume;
        ReactionRate as Real(Unit='mol/m^3/s');
        ReactionConstant as Real(Unit='m^3/s/mol');
        TotalH as energy;
        H(4) as energy_mol;
        T as temperature;
        Q as power;

        EQUATIONS
        "Mass balance"
        diff(Holdup) = Nu*ReactionRate*TotalVolume;

        "Energy balance"
        diff(TotalH) = ReactionRate*TotalVolume*(-Hr) + Q;

        "Reaction rate"
        ReactionRate = ReactionConstant*C(1)*C(2);
        
        "Reaction constant"
        ReactionConstant = k0 * exp(-Ea/(R*T));

        "Concentrations"
        Holdup = C * TotalVolume ;
        TotalVolume = sum(Holdup/Density) ;
        "Reactor energy content"
        TotalH = sum(Holdup*H) ;
        H = Ac * ( T - Tref ) + Bc * ( T^2 - Tref^2 ) / 2 ;
end

FlowSheet BatchProcess as BatchReactor
        SET
        Density = [ 10.0E3, 8.0E3, 11.0E3, 11.0E3 ] * 'mol/m^3' ;
        Ac = [ 150.0, 175.0, 200.0, 175.0 ] * 'J/mol/K';
        Bc = [ 0.0, 0.0, 0.0, 0.0 ] * 'J/mol/K/K' ;
        Nu = [ -1, -1, +1, +1 ] ;
        
        k0 = 5E-4 * 'm^3/s/mol';
        Ea = 15000.0 * 'J/mol';
        Hr = 55000.0 * 'J/mol';
        
        SPECIFY
        Q = 505e3 * 'W';
        
        INITIAL
    Holdup = [100.0, 100.0, 1.0, 1.0 ] * 'mol';
    T = 293.0 * 'K';
        
end

Estimation BatchEst as BatchProcess
        ESTIMATE
        # PARAMETER  START     LOWER     UPPER     UNIT
        k0           0.0008    0.0001    0.0010    'm^3/s/mol';
        Ea           12000     10000     20000     'J/mol';
        Hr           50000     35000     75000     'J/mol';

        EXPERIMENTS
        # DATA FILE             WEIGTH
        "BatchReactor.dat"      1;
        
        OPTIONS
        
        NLPSolver(
                MaxIterations = 1000,
                File = "complex"
                #File = "ipopt_emso"
        );
end