# ODE Implicit BVP ETH Zrich Eidgenssische Technische pellet

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