Docstoc

ODE Implicit BVP ETH Zrich Eidgenssische Technische pellet

Document Sample
ODE Implicit BVP ETH Zrich Eidgenssische Technische pellet Powered By Docstoc
					     Integration of Ordinary
     Differential Equations
          dy/dt = f(t,y)

                Marco Lattuada
    Swiss Federal Institute of Technology - ETH
Institut für Chemie und Bioingenieurwissenschaften
ETH Hönggerberg/ HCI F135 – Zürich (Switzerland)

          E-mail: lattuada@chem.ethz.ch
http://www.morbidelli-group.ethz.ch/education/index
                Exercise: Simulation of a Batch Reactor
                                                     • Hp: T = const & V = const

                                                     • Variables:            CA, CB , CC

                                                     • Equations:

                                                     R1: A → 2B with R1 = k1A
                                                     R2: B → C with R2 = k2B




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 2
                                               Analytical Solution
 Suppose to have the following system:



 where A is a (N x N) matrix of constant coefficients.
 The analytical solution is:



 where: B is a matrix of constant coefficients
        l the vector of the eigenvalues of A

 The coefficients of the matrix B can be found by imposing that the solution satisfies
 the initial differential equation and the initial values of y




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 3
                                               Analytical Solution
 The eigenvalues of J are:




 an the solution of the system of differential equations is (C A = 1 and CB = 0 at t = 0):




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 4
                                                      Forward Euler
 Forward Euler recursive formula:



 If we substitute the expression for our problem:



 we obtain:



 According to the largest eigenvalue of J, the maximum step size is:




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 5
                                 Solution with Forward Euler




                                                                             1+hl = 1+0.2*(-10) = -1


Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 6
                                                         Stiff Systems
 The system of ODE is stiff when:




 In this case:
 • the solution of one component requires very small steps to be stable
 • the precision on the other components is much higher than the required precision

 Physical interpretation:
 • There is one mechanism (e.g. of reaction) which is much faster than the others




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 7
                                                    Backward Euler
 Backward Euler recursive formula:



 If we substitute the expression for our problem:



 we obtain:



 If we consider a single ODE:

                             A-stable algorithm: k is smaller than zero for every hl < 0
                                                      for h > 0 (l < 0)
                             Strongly A-stable algorithm: k tends to zero for hl → -∞



Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 8
                                                  Trapezoidal Rule
 Trapezoidal Rule recursive formula:



 If we substitute the expression for our problem:



 we obtain:




 If we consider a single ODE:




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 9
                                Solution with Trapezoid Rule




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 10
                                                       Exercise No. 1
 Implicit Methods
 1. Use the following values of the parameters: k 1 = 1, k2 = 100, and the following
    initial values: CA(t=0) = 1 and CB(t=0) = 0.
 2. Integrate the system of ODEs from t = 0 to t = 5 with the Forward and Backward
    Euler algorithms. Which is the maximum step size for both algorithms needed to
    obtain a stable solution?
 3. Compare the numerical solution obtained with the two algorithms at each time
    step with the analytical solution on both C A and CB. For each time step compute
    the error of the numerical solution and find the maximum of the relative error for
    both CA and CB. Which is the maximum step size (h) for the two algorithms
    which allows to have a maximum relative error on the value of C B of 0.1%?
    Which is the corresponding maximum relative error on C A?
 4. Repeat the integration using the Matlab algorithms ode113 (multi-step solver for
    non-stiff problems) and ode15s (multi-step solver for stiff problem). Compare
    the two results and the number of steps needed by the two solvers. Note that the
    call of the two algorithm is identical.
Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 11
                          Boundary Value Problems (BVP)
 Definition of the problem:
 Integration of a system of ordinary differential equations (ODEs):




 Subject to the Boundary Conditions (B.C.):


                     The total number of B.C. has to be equal to the
                     number of equations!


Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 12
                                                   Shooting method
Idea: transform the BVP in an initial value problem (IVP), by
guessing some of the initial conditions and using the B.C. to refine
the guess, until convergence is reached



                                                                                          Target

                                                     Too high: reduce the initial velocity!


                 Convergence can be problematic
                                                       Too low: increase the initial velocity!
                    Use the same algorithms used for IVP


Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 13
                                               Collocation method
 In order to solve a boundary value problem, a common technique is the collocation
 method.
 It is based on approximating the unknown function with the sum of polynomials
 Pi(t), multiplied by unknown coefficients:




 The coefficient are determined by forcing the approximated solution to satisfy the
 differential equations at a number of points equal to the number of unknown
 coefficients.
 Matlab has a built-in routine, “bvp4c” which implements the collocation method,
 which can also solve singular value problems in the form:

                                                                             S is a constant matrix

Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 14
                    Example of Boundary Value Problem
 Diffusion and reaction of a gaseous species inside a catalyst pellet:
                                                         Calculate the concentration profile of a
                                                         gaseous species diffusing and reacting inside
                     Cs
                                                         a spherical porous catalyst pellet.
                                                         Determine the efficiency of the pellet

                           0                   Rp Data:
                                                  Cs= concentration at the surface (outside)
                                                  ks= reaction rate constant per unit surface
                                                  Sv = surface to volume ratio in the pellet
                                                  Rp= pellet radius
                                                  D= diffusion coefficient of the gaseous
                                                  component

Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 15
                                                     Mass Balance
 The difference between the diffusive flux of gas entering the shell and the diffusing
   flux of gas leaving the shell equal the rate of consumption of gas.

                                                                 Diffusive Flux=

                               r
               dr

                                                                                   Reaction rate=

 ACCUMULATION=IN-OUT+PRODUCTION




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 16
         Diffusion and reaction inside a catalyst pellet
 Mass Balance equation (in spherical coordinates, steady state):


                                                                                            r
                                                                                          r+dr
       Mass Transfer                             Chemical
       by Diffusion                            Reaction Rate

 B.C.                                                         No discontinuity of the profile at
                                                              the center of the pellet

                                                               Concentration at the outer
                                                               surface known

Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 17
                      Dimensionless form of the equations
 By defining the following dimensionless variables:




        Equation in dimensionless form: when the equation is
        written in dimensionless form, the effect of all physical
                              f = Thiele Modulus: reaction
        parameters is included in f. The solution is numerically
                                    rate/diffusion rate
        much simples if variables are dimensionless (most of them
        are in the B.C. of be rewritten
 the equation and range canvalues [0-1]). as follows:




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 18
                              Efficiency of the catalyst pellet
 The efficiency of the catalyst pellet is defined as follows:
                                                                        It is the ratio between the cumulative
                                                                          reaction rate in the pellet over the
                                                                               reaction rate in the case the
                                                                        concentration was everywhere equal
                                                                              to the surface concentration.

 By manipulating the expression in the definition given above, and
 using the dimensionless variables, it can be rewritten as:

                                                             The efficiency depends only on f and on
                                                            the slope of the profile at the surface.


Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 19
                                              Non-isothermal case

 In the case the catalyst pellet cannot be considered isothermal, one
 needs to integrate also the energy balance (with B.C.):




                                                                    kT= thermal conductivity
                                                                    -DHR= heat of reaction
                                                                    Ts= surface temperature




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 20
                         How to deal with the temperature
 In this case there is no need to solve the thermal balance, since it
 can be proved that:



 Therefore, it suffices to solve the following equation:




                                   Dimensionless activation energy (sensitivity of the
                                   reaction to the change of T


                                                         Prater Number: >0, exothermic, <0 endothermic


Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 21
                                                       Exercise No. 2
 Consider the dimensionless equation for the material balance of the gaseous
    component inside the catalyst pellet
 1. Solve the material balance equation for the isothermal case (Eqs. on page 17)
    using “bvp4c”, and calculate the efficiency of the catalyst pellet as a function of
    the Thiele modulus, in the range between 10-2 and 102 for a first and second
    order reaction
 2. Solve the material balance equation for the non-isothermal case (Eqs. on page
    20) using “bvp4c”, and calculate the efficiency of the catalyst pellet as a function
    of the Thiele modulus, in the range between 10-2 and 102 for a first order
    reaction, for Prater number values b=-0.4, and b=+0.4 and g=10.




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 22
                                                  How to use bvp4c
Example of how to call bvp4c:




function dydx = fun_ode(x,y)
 dydx = [ y(2); y(1)];

function res = fun_bc(ya,yb)
 res = [ ya(1)-1; yb(1)];

S=[0,0;0,-1];

options = bvpset('SingularTerm',S,’ 'FJacobian',@fun_Jac);


Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 23
                                             How to use bvp4c - 2
 solinit = bvpinit(linspace(0,1,30),@initfun);

 initfun(x): function that calculates the initial guess of the solution

 function yinit=initfun (x)
 yinit = [ exp(-x)
        -exp(-x) ];

 function dydx = fun_Jac(x,y)
  dydx = [ 0,1; 1,0];

 sol = bvp4c(@fun_ode,@fun_bc,solinit,options);




Marco Lattuada – Statistical and Numerical Methods f or Chemical Engineers
Implicit ODE Solvers – Page # 24

				
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
views:12
posted:5/16/2011
language:German
pages:24