Numerical Solution of Coupled Differential Equations

1
Numerical Solution of Coupled Differential
Equations

• Packed Bed Micro-Reactor Case

2
PBR Micro-Reactor Schematic Representation
Flow per channel
FTO/10 6 x 10-8 mol/s
Po 1.5 atm
Dia 400 mm
(yCH3OH)0 1
FTO 1 x 10-7 mol/s
3
Overall Objective Compute Pressure and
Conversion Profile
P ? X?
P, X
z
Po Xo
Wcat or Z
4
Coupled Differential Equations
1. GMBE Differential Form
2. Pressure Drop in Packed Bed Differential
Form
5
Integration by ode23 and ode45 Matlab Command

• w, y ode45 (pbr, w0,wf, y0)
• where
• pbr is a string variable containing the name
  of the m-file for the derivatives.
• w0 is the inlet catalyst mass
• wf is the total catalyst mass from inlet to
  outlet
• y0 is the initial condition vector for the state
  variables
• w a (column) vector of catalyst mass
• y an array of state variables as a function of
  catalyst mass

6
Matlab script file syntax and other details
7
Purpose of function files
function outputfunction_name (input1, input2)
As indicated above, the function file generates
the value of outputs every time it called upon
with certain sets of inputs of dependent and
independent variables For instance the pbr.m
file generates the value of output (dy), every
iteration it is called upon with inputs of
independent variable catalyst weight (w) and
dependent variables (y) NOTE For pbr.m file,
the output dy is actually dX/dW and dP/dW
function dypbr (w, y)
8
The function pbr returns the numerical values of
dy or OUTPUTS for corresponding sets of INPUTS y
defined at given value of y
function dypbr(W,y) global dtube operating
conditions To300273 Inlet temperature Po1.51
01.325 Inlet pressure in kPa
In our case, OUTPUTS (dy) INPUTS d(P)/dW y(2,)
P d(X)/dW y(1,) X
9
Function file name
function dypbr(W,y) global dtube operating
conditions To300273 Inlet temperature Po1.51
01.325 Inlet pressure in kPa
10
Semi-colon suppresses the SCREEN PRINTING of the
statement To300273
function dypbr(W,y) global dtube operating
conditions To300273 Inlet temperature Po1.51
01.325 Inlet pressure in kPa
11
defined variables Py(1,) Xy(2,)
For convenience sake, I have defined the variable
y(1,) as Pn and variable y(2,) as X
12
defined variables Py(1,) Xy(2,)
COLON indicates that the y is a variable. That
is y(1,) varies with time
13
momentum balance equation or pressure drop
equation dy(1,)-(lamda/dw_to_dz_conv)(1/gasden)
mole balance equation in terms of conversion
changing with catalyst weight or
dX/dw dy(2,)(kfpCH3OH)/Fao
These set of equations calculate the values of
the OUPUTS dy based on previously defined/known
constants and/or INPUTS
14
RUN FILE
15
This file contains executable commands only
global dtube dtube 4e-4 microreactor diameter
400 micron Po1.5101325 y0 Po 0 initial
condition for P and X. W00 Wmax3.1416(dt/2)2
(5e-2)(1-0.3)1400 wpir2L(1-porosity)cat
_density WspanW0 Wmax Defines integration
range W,yode15s('pbr', Wspan, y0, ) hold
on plot(W,y(,1)/Po,'b-.', W,y(,2),'r-')
16
This statement defines the initial condition of
the dependent variable
global dtube dtube 4e-4 microreactor diameter
400 micron Po1.5101325 y0 Po 0 initial
condition for P and X. W00 Wmax3.1416(dt/2)2
(5e-2)(1-0.3)1400 wpir2L(1-porosity)cat
_density WspanW0 Wmax Defines integration
range W,yode15s('pbr', Wspan, y0, ) hold
on plot(W,y(,1)/Po,'b-.', W,y(,2),'r-')
17
This statement defines the Limits of Integration
over which the Differential Equation must be
INTEGRATED over.
global dtube dtube 4e-4 microreactor diameter
400 micron Po1.5101325 y0 Po 0 initial
condition for P and X. W00 Wmax3.1416(dt/2)2
(5e-2)(1-0.3)1400 wpir2L(1-porosity)cat
_density WspanW0 Wmax Defines integration
range W,yode15s('pbr', Wspan, y0, ) hold
on plot(W,y(,1)/Po,'b-.', W,y(,2),'r-')
18
This is the Main COMMAND. An output in form of a
matrix of W and ys will be generated by solving
the differential equations defined/calculated by
pbr The ODE solver ode15s is used to
integrate the differential equations over the
integration limit defined in Wspan as WW0 to
WWmax given initial conditions for ys as y0
W00 Wmax3.1416(dtube/2)2(5e-2)(1-0.3)1400
wpir2L(1-porosity)cat_density WspanW0
Wmax Defines integration range W,yode15s('p
br', Wspan, y0, ) hold on plot(W,y(,1)/Po,'b
-.', W,y(,2),'r-')
19
Another useful command is
wk1write(filename', t, y)
This command will write the output t, y that
is, all ys as a function of t obtained from
the ODE solver in a file called filename
20

• Getting started with Matlab

21
Requirements

• You need Matlab on the computer
• Computers in the Dupuis Cluster and some in ILC
  have Matlab pre-loaded

22
Matlab Environment

• Double click on MATLAB icon on desktop
• Search for Matlab if icon does not appear on
  desktop
• A Matlab window should appear with the following
  prompt
• EDUgtgt

23
Current Directory

• Matlab will open up in default directory
• Change the Current directory to the location
  where the function and run files are stored.

Note You need create or have a copy of the
function files before you can run the files
24
Running the files

• Type the runfile name at the Matlab command, i.e.
• EDUgtgt pbrrun