Transcript Boundary Value Problems and Partial Differential Equations (PDEs)

Boundary Value Problems and
Partial Differential Equations (PDEs)
 d 2 y (t , z ) dy (t , z )

dy (t , z )
 f
,
, y, t , z 
2
dt
dz
 dz

Daniel Baur
ETH Zurich, Institut für Chemie- und Bioingenieurwissenschaften
ETH Hönggerberg / HCI F128 – Zürich
E-Mail: [email protected]
http://www.morbidelli-group.ethz.ch/education/index
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
1
Boundary Value Problems (BVP) for ODEs
 Problem definition:
 Find a solution for a system of ODEs
dy (t )
 y  f (t , y )
dt
 Subject to the boundary conditions (BCs):
g ( y (t0 ), t0 )  0
h( y (t f ), t f )  0
 The total number of BCs has to be equal to the number of
equations!
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
2
Shooting Method
 A first approach is to transform the BVP into an initial value
problem (IVP), by guessing the missing initial conditions
and using the BC to refine the guess, until convergence is
reached
Too high: reduce initial velocity!
Too low: increase initial velocity!
Target
 This way, the same algorithms as for IVPs can be used,
but the convergence can be very problematic
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
3
Collocation Method
 A more sophisticated approach is the collocation method;
it is based on approximating the unknown function with a
sum of polynomials multiplied with unknown coefficients
N 1
y (t )   ai Pi (t )
i 1
 The coefficients are determined by forcing the
approximated solution to satisfy the ODE at a number of
points equal to the number of coefficients
 Matlab has a built-in function bvp4c which implements this
method; it can also solve singular value problems
dy
y
 S  f ( y, t )
dt
t
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
4
Example of a BVP
 Consider a tubular reactor
cin
cout
nA  B
kn
0
L
 We can model it as a plug flow reactor (PFR) with backmixing by using the following partial differential equation
c A (t , x)
 2cA
c A
n
 Dax

v

k
c
n A
t
x 2
x
Dax is the (effective) axial dispersion coefficient [m2/s],
v is the linear flow velocity [m/s] and n is the reaction order
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
5
Tubular Reactor: Dimensionless Form
 Let us cast the model in dimensionless form by defining
t v
 

L
cA
u
cin
t
x
z
L

u
1  2u u
n



Da

u
 Pe z 2 z
The numerical solution of a problem is usually much simpler if it
is dimensionless (most variables will range from 0 to 1).
 Where Pe is the Peclet and Da is the Damköhler number
Lv
Pe 
,
Dax
L  kn   cin 
Da 
v
n 1
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
6
Steady State Assumption, Boundary Conditions
 By assuming steady state, the time variable vanishes and
we get an ODE
1 d 2 u du
n
0


Da

u
Pe dz 2 dz
 This equation is subject to the Danckwerts BCs (mass
balance over inlet, continuous profile at the outlet)
u ( z  0)  1 
du
dz
1 du

Pe dz
z 0
0
z 1
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
7
Transformation into a first order ODE
 bvp4c solves first order ODEs, so if we remember the
«trick» and transform our ODE, we get
y1  u
y2
dy1 du

dz dz
dy1
 y2
dz
dy2 d 2u
 du

 2  Pe    Da  u n   Pe   y2  Da  y1n 
dz dz
 dz

 And the boundary conditions
1
y1 ( z  0)  1 
 y2 ( z  0)
Pe
y2 ( z  1)  0
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
8
Partial Differential Equations
 Problem definition:
In a partial differential equation (PDE), the solution
depends on more than one independent variable, e.g.
space and time
 The function is usually subject to both inital conditions and
boundary conditions
 In our example
u (, z )
1  2u u
n



Da

u

Pe z 2 z
u (  0)  0
z
plus the Danckwerts BCs which apply at all times
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
9
Characterization of Second Order PDEs
 Second order PDEs take the general form
Au xx  2 Bu xy  Cu yy 
0
where A, B and C are coefficients that may depend on x
and y
 These PDEs fall in one of the following categories
1. B2 – AC < 0: Elliptic PDE
2. B2 – AC = 0: Parabolic PDE
3. B2 – AC > 0: Hyperbolic PDE
 There are specialized solvers for some types of PDEs,
hence knowing its category can be useful for solving a PDE
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
10
Numerical Solution of PDEs
 In general, it can be very difficult to solve PDEs numerically
 One approach is to discretize all but one dimension of the
solution; this way a system of ODEs is obtained that can be
solved more easily
 Note that these ODE systems are usually very stiff
 There are different ways of discretizing a dimension, for
example the collocation method we saw earlier
 Sophisticated algorithms refine the discretization in places
where the solution is still inaccurate
 Matlab has a built-in solver for parabolic and elliptic PDEs
in two dimensions, pdepe
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
11
Unsteady State Tubular Reactor
 Let us consider the start-up of a tubular reactor, i.e.
u (, z )
1  2u u
n



Da

u

Pe z 2 z
u (  0)  0
z
1 du
u ( z  0)  1 

Pe dz
du
dz
0

z 0

z 1
 We can easily see that this is always a parabolic PDE
(B = C = 0), hence the Matlab solver is applicable
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
12
Method of Finite Differences
 We can apply the so-called finite differences method, if we
remember numerical differentiation
df ( x)
dx
f ( x0  h)  f ( x0  h)
2h
 Also, we can easily derive a similar expression for the
second order derivative
d 2 f ( x)
dx 2
f ( x0  h)  2 f ( x0 )  f ( x0  h)
h2
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
13
Method of Finite Differences (Continued)
 The PDE for the start-up of the tubular reactor reads
u (, z )
1  2u u
n



Da

u

Pe z 2 z
 Applying the method of finite differences, we get
du z
1 u z z  2u z  u z z u z z  u z z
n



Da

u
z
d Pe
z 2
2z
where Δz = 1/N is the discretization step, N being the
number of grid points
 If we number the grid points with i = 1...N, we get
dui
1 ui 1  2ui  ui 1 ui 1  ui 1
n



Da

u
i
d Pe
1/ N 2
2/ N
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
14
Method of Finite Differences (Continued)
 What happens at the boundaries u1 and uN?
 One possibility is to invent pseudo grid-points u0 and uN+1
that fulfill the boundary conditions
 In our case
u0
u0
1 u
 1
Pe z
z 0
1  u1  u0 
 1
Pe z
1  N / Pe  u1 


1  N / Pe
u N 1  u N
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
15
Method of Finite Differences (Continued)
 Rearranging the equations gives us
 N2 N 
 N2
 N2 N 
dui
n 1 
 ui 1 
   ui  2
 Da  ui   ui 1 
 
d
 Pe 2 
 Pe

 Pe 2 
 With the initial conditions and «boundary conditions»
ui (  0)  0,
u0
i 1
N
1  N / Pe  u1 


1  N / Pe
u N 1  u N
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
16
Assignment 1
1. Solve the dimensionless tubular reactor using bvp4c for
50 different values of the Peclet number between 0.01 and
100 and for a reaction of first and second order
 In both cases use Da = 1
2. Plot the conversion at the end of the reactor (1-cA/cin) vs.
Peclet for both reaction orders; Also plot the ratio between
the conversions of the first order and second order
reaction
 What is better for these reactions, a lot of back-mixing (Pe small,
CSTR) or ideal plug flow (Pe large, PFR)?
 What influence does the reaction order have overall and at low or
high Peclet numbers?
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
17
Usage of bvp4c
 bvp4c uses a call of the form
 sol = bvp4c(@ode_fun, @bc_fun, solinit, options);
 ode_fun is a function taking as inputs a scalar t and a vector y,
returning as an output dy / dt
 bc_fun is a function taking as inputs vectors where the boundary
conditions are evaluated, returning as output the residual at the
boundary
 solinit initializes the solution by using
solinit = bvpinit(range, @initfun)
 options is an options structure resulting from
options = bvpset('FJacobian', @jac_fun)
 sol is a struct containing the solution and other parameters
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
18
Usage of bvp4c (Continued)
 In our case use the following
 dy1
 y2
dy  dz

dz  dy2
 Pe   y2  Da  y1n 
 dz
1

 ya (2)
 ya (1)  1 
residual  
Pe
 yb (2)
0

J
n 1
 Pe  Da  n  y1
1
Pe 
function dy = ode_fun(t,y,...)
( 0)
( 0)
function res = bc_fun(ya,yb,...)
function J = jac_fun(t,y,...)
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
19
Assignment 2
1. Use pdepe to solve the startup of the tubular reactor.
 Consider only the first order reaction with Pe = 100 and Da = 1.
 Plot the conversion at the end of the reactor vs. dimensionless time.
 At what time does the solution reach steady state, i.e. how many
reactor volumes of solvent will you need? Compare the solution to
what you have found in assignment 1, if the difference is smaller
than 0.1%, assume that steady state has been reached.
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
20
Usage of pdepe
 pdepe uses the following syntax
 sol =
pdepe(m,@pde_fun,@ic_fun,@bc_fun,xmesh,tspan);
 m is a parameter that describes the symmetry of the problem (slab =
0, cylindrical = 1, spherical = 2); in our case slab, so m = 0
 @pde_fun is a function that describes the PDE in this form:
u  u
 

c  x, t , u , 
 xm  xm f
x  t
x 

u   
u 

x
,
t
,
u
,

s
x
,
t
,
u
,




x   
x 

 @ic_fun is a function that takes as input a vector x and returns the
initial conditions at t = 0
 @bc_fun is a function that describes the boundary conditions, taking
as an input xl, ul, xr, ur and t, returning the BCs in a form:
u 

p  x, t , u   q  x, t  f  x , t , u ,   0
x 

that is, pl, ql, pr and qr
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
21
Usage of pdepe (Continued)
 In our case, use the following
u  u
 

c  x, t , u , 
 xm  xm f
x  t
x 

u   
u 

 x, t , u ,    s  x, t , u , 
x   
x 

u (, z ) 1  2u u

  Da  u n
2

Pe z
z
u 

p  x, t , u   q  x, t  f  x , t , u ,   0
x 

1 u
1

 u ( z  0)  0
Pe z z 0
du
dz
4 Variable Inputs
5 Variable Inputs
0
z 1
Daniel Baur / Numerical Methods for Chemical Engineers / BVP and PDE
22