Report

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 Lv 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 2z 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 , xm 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 , xm 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