Numerical Solution of Coupled Differential Equations

1. Numerical Solution of Coupled Differential Equations Packed Bed Micro-Reactor Case

2. 5 cm or 0.05 m 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 wfis 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 output=function_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 dy=pbr (w, y)

8. The function pbr returns the numerical values of dyor OUTPUTS for corresponding sets of INPUTSydefined at given value ofy function dy=pbr(W,y) global dtube % operating conditions To=300+273; %Inlet temperature Po=1.5*101.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 name should match file name Function file name function dy=pbr(W,y) global dtube % operating conditions To=300+273; %Inlet temperature Po=1.5*101.325; %Inlet pressure in kPa

10. Semi-colon suppresses the SCREEN PRINTING of the statement “To=300+273” function dy=pbr(W,y) global dtube %operating conditions To=300+273; %Inlet temperature Po=1.5*101.325; %Inlet pressure in kPa

11. %defined variables P=y(1,:); X=y(2,:); For convenience sake, I have defined the variable y(1,:) as Pn and variable y(2,:) as X

12. %defined variables P=y(1,:); X=y(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,:)=(kf*pCH3OH)/Fao; These set of equations calculate the values of the OUPUTSdy 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 Po=1.5*101325; y0 = [Po 0]; %initial condition for P and X. W0=0; Wmax=3.1416*(dt/2)^2*(5e-2)*(1-0.3)*1400;% w=pi*r^2*L*(1-porosity)*cat_density Wspan=[W0 Wmax];% Defines integration range [W,y]=ode15s('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 Po=1.5*101325; y0 = [Po 0]; %initial condition for P and X. W0=0; Wmax=3.1416*(dt/2)^2*(5e-2)*(1-0.3)*1400;% w=pi*r^2*L*(1-porosity)*cat_density Wspan=[W0 Wmax];% Defines integration range [W,y]=ode15s('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 Po=1.5*101325; y0 = [Po 0]; %initial condition for P and X. W0=0; Wmax=3.1416*(dt/2)^2*(5e-2)*(1-0.3)*1400;% w=pi*r^2*L*(1-porosity)*cat_density Wspan=[W0 Wmax];% Defines integration range [W,y]=ode15s('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 y’s 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 W=W0 to W=Wmax given initial conditions for y’s as y0 W0=0; Wmax=3.1416*(dtube/2)^2*(5e-2)*(1-0.3)*1400;% w=pi*r^2*L*(1-porosity)*cat_density Wspan=[W0 Wmax];% Defines integration range [W,y]=ode15s('pbr', 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 y’s 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 • EDU>>

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. EDU>> pbrrun