Document Sample

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 |

OTHER DOCS BY mikeholy

How are you planning on using Docstoc?
BUSINESS
PERSONAL

Feel free to Contact Us with any questions you might have.