source: trunk/sample/miscellaneous/tenprobs/prob07.mso @ 494

#*---------------------------------------------------------------------
* 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.
*
*----------------------------------------------------------------------
* 7. Diffusion with chemical reaction in a one dimensional slab
*----------------------------------------------------------------------
*
*   Description:
*      This problem is part of a collection of 10 representative
*       problems in Chemical Engineering for solution by numerical methods
*       developed for Cutlip (1998).
*
*   Subject:
*       * Transport Phenomena
*               * Reaction Engineering
*
*       Concepts utilized:
*               Methods for solving second order ODEs with 2 point boundary
*       values typically used in transport phenomena and reaction kinetics.
*
*       Numerical method:
*               * Simultaneous ODEs with split boundary conditions
*               * Resolved by finite difference method
*
*   Reference:
*               * CUTLIP et al. A collection of 10 numerical problems in 
*       chemical engineering solved by various mathematical software 
*       packages. Comp. Appl. in Eng. Education. v. 6, 169-180, 1998. 
*       * More informations and a detailed description of all problems
*       is available online in http://www.polymath-software.com/ASEE
* 
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* GIMSCOP/UFRGS - Group of Integration, Modeling, Simulation, 
*                                       Control, and Optimization of Processes
* $Id$
*--------------------------------------------------------------------*#
using "types";



Model problem
        PARAMETERS
outer N as Integer      (Brief="Number of discrete points", Lower=3);
        
        Co              as conc_mol             (Brief="Constant concentration at the surface");
        D               as diffusivity  (Brief="Binary diffusion coefficient");
        k               as Real                 (Brief="Homogeneous reaction rate constant", Unit='1/s');
        L               as length               (Brief="Bottom surface");
        
        
        VARIABLES
        C(N+2)  as conc_mol             (Brief="Concentration of reactant");
        z(N+2)  as length               (Brief="Distance", Default=1e-3);
        dz              as length_delta (Brief="Distance increment");   


        EQUATIONS
        "Discrete interval"
        dz = (z(N+2) - z(1))/(N+1);
        
        for i in [2:(N+1)]
        "Concentration of reactant"     
                (C(i+1) - 2*C(i)+ C(i-1))/(z(i) - z(i-1))^2 
                        = (k/D)*C(i);
                
        "Discrete length"
                z(i) = z(i-1) + dz;
        end

        # Boundary conditions
        "Initial and boundary condition" # z = 0
        C(1) = Co;
        
        "Upper boundary" # z = L
        (C(N+2) - C(N+1))/(z(N+2) - z(N+1)) = 0*'kmol/m^4';
        
        
        SET
        Co= 0.2*'kmol/m^3';
        D = 1.2e-9*'m^2/s';
        k = 1e-3/'s';
        L = 1e-3*'m';
end



FlowSheet numerical_solution
        PARAMETERS
        N               as Integer;
        
        
        DEVICES
        reac    as problem;
        
        
        SET
        N = 10; # Number of discrete points
        
        
        SPECIFY
        "Lower boundary"
        reac.z(1) = 0*'m';

        "Upper boundary"
        reac.z(N+2) = reac.L;


        OPTIONS
        Dynamic = false;
end



FlowSheet comparative
        PARAMETERS
        N               as Integer;
        
        
        VARIABLES
        C_(N+2) as conc_mol     (Brief="Concentration of reactant by analytical solution");
        
        r_              as Real         (Brief="Pearson product-moment correlation coefficient");
        
        Cm              as conc_mol     (Brief="Arithmetic mean of calculated C");
        C_m             as conc_mol     (Brief="Arithmetic mean of analytical C");
        
        
        DEVICES
        reac    as problem;
        
        
        SET
        N = 10; # Number of discrete points
        
        
        EQUATIONS
        "Analytical solution"
        C_ = reac.Co*cosh(reac.L*sqrt(reac.k/reac.D)*(1 - reac.z/reac.L))
        /cosh(reac.L*sqrt(reac.k/reac.D));

        "Pearson correlation coefficient" # used by softwares like MS Excel
        r_ = (sum((reac.C - Cm)*(C_ - C_m)))/sqrt(sum((reac.C - Cm)^2)*sum((C_ - C_m)^2));
        
        "Arithmetic mean of C"
        Cm = sum(reac.C)/(N+1);
        
        "Arithmetic mean of C_"
        C_m = sum(C_)/(N+1);
        
        
        SPECIFY
        "Lower boundary"
        reac.z(1) = 0*'m';

        "Upper boundary"
        reac.z(N+2) = reac.L;


        OPTIONS
        Dynamic = false;
end


FlowSheet analytical_solution
        PARAMETERS
        Co              as conc_mol;
        L               as length;
        k               as Real(Unit='1/s');
        D               as diffusivity;
        
        
        VARIABLES
        C               as conc_mol             (Default=0.2);
        z               as length               (Default=1e-3);
        
        
        EQUATIONS
        "Change time in z"
        z = time*'m/s';
        
        "Analytical solution"
        C = Co*cosh(L*sqrt(k/D)*(1 - z/L))/cosh(L*sqrt(k/D));
        
        
        SET
        Co= 0.2*'kmol/m^3';
        D = 1.2e-9*'m^2/s';
        k = 1e-3/'s';
        L = 1e-3*'m';
        
        
        OPTIONS
        TimeStart = 0;
        TimeStep = 1e-6;
        TimeEnd = 1e-3;
end