Unformatted text preview: Page 1 of 6 MATLAB Review + Simulation CHE 361 A vertical cylindrical tank is constructed with holes in its side such that the flow rate of liquid
out of the tank is proportional to the height of liquid in the tank. The height of the tank is 2.00
meters and the circular cross section (A) is 1 m2. When the liquid inlet flow (Fi) is constant at
2.00 m3/min, the tank is initially half full (at a steady state for t < 0). If the inlet flow is suddenly
changed to 5.00 m3/min, how much time will elapse, to the nearest second, before the liquid
overflows from the top of the tank?
A dh(t )
= Fi (t ) − Fo (t ) = Fi (t ) − Kh(t )
dt or simplifying for the derivative only, with h(t = 0) = h0 initial condition
dh Fi − Kh
=
= the RHS of the ODE
dt
A Using Laplace transforms, we solved for the general case with Fi (t ) = Fi for all t ≥ 0
K
Kt
− Kt
− Kt
⎞ −A
⎞
Fi ⎛
F i ⎞ −A t F i ⎛ F i
Fi ⎛
A
h(t ) =
+ ⎜ h0 − ⎟ e =
−⎜
− h0 ⎟ e = h0 e + ⎜ 1 − e A ⎟
K⎝
K⎠
K ⎝K
K⎝
⎠
⎠ With A = 1.00 m, h0 = 1.00 m 2 , Fi = 5.00 m3 / min and K = 2.00 m 2 / min
h(t ) = 2.5 − 1.5e −2t for h(t ) ≤ H = the height of the tank = 2.00 m Here is a simple MATLAB program h_of_t.m to show how the ODE for the height of liquid
in the tank can be solved numerically using the builtin MATLAB function ode45 and then use
MATLAB commands to make a plot of h(t). as a solid line with circles at calculated points. Page 2 of 6 Contents of the program file named = h_of_t.m % h_of_t.m
% An example MATLAB program for Check Problem 3 simulation.
% Dr. K.Levien Chemical Engineering 361 Jan. 13, 201X
% Oregon State University
%
format short, format compact
global A Fi K
% variables shared with functions
%
A = 1.00;
% cylindrical tank cross sectional area (m^2)
Fi = 5.00;
% constant feedflow at t >= 0 (m^3/min)
K = 2.00;
% proportional constant for flow out of tank = Kh (m^2/min)
%
h0 = 1.00;
% initial height of liquid in tank (m)
tend = 1.00;
% end point of simulation (min)
tspan = [0 tend]; % solve ODE from time = 0 to time = tend
%
[t_min,h] = ode45(@dhdt,tspan,h0);
%
t_sec = t_min*60; % time in the ode was in min, convert to s for plot
%
figure
% open a new Figure Window
plot(t_sec,h,'bo','LineWidth',2) % use thicker line for copy and paste
xlabel('Time (s)'), ylabel('Height (m)'), grid
title('Liquid Height in Tank, Benny Beaver CHE 361 1/13/1X') Contents of the function file named = dhdt.m function [RHS] = dhdt(t,h)
%dhdt.m = RHS function to calculate derivative of Prob3 tank
% Dr. K.Levien Chemical Engineering 361 Jan. 13, 201X
% Oregon State University
%
global A Fi K
% variables shared with MAIN program
%
RHS = (Fi  K*h)/A; Often MATLAB is used to compare experimental results versus a theoretical model.
A second example P3tank.m showns how this can be done and then creates a second plot
where the numerical solution of the ODE simulation is compared to the analytical solution
found using the Laplace transform by plotting both results as curves of h(t) .
In the Command Window we simply type the name of the program file (without the .m
extention). Each “input” command allows us to specify a number at “run time”. >> P3tank
Enter constant feedflow at t >= 0 (m^3)/min (5.0) : 5
Enter duration of simulation in minutes (1.0) : 1
>> Page 3 of 6 Here are two plots corresponding to two common situations:
Figure 1 = Plot a model curve versus points of data.
Figure 2 = Plot simulations based on two models. In this example there is only one model, the
linear ordinary differential equation, but it is solved 2 ways = an analytical solution using
Laplace transforms and a numerical way using the MATLAB builtin function “ode45”.
We use a MAIN program file with name “P3tank.m” and a function file with name “dhdt.m”.
Window copy and paste: Page 4 of 6 Figure copy and paste after adjusting font size and line size using FilePreferencesFigure Copy
Template Using the Preferences you can copy Figures and paste into reports without including the Page 5 of 6 Microsoft Window “look”. Sample Prob. 3 Model + Data, Benny Beaver CHE 361 1/13/12
2.5 Height (m) 2 1.5 Model
Data
1 0 10 20 30
Time (s) 40 50 60 Sample Prob. 3 Analytical + Numerical, Benny Beaver CHE 361 1/13/12
2.5 Height (m) 2 1.5 Analytical
Numerical
1 0 10 20 30
Time (s) 40 50 60 Page 6 of 6 First line of MAIN program file = P3tank.m
% P3tank.m
% An example MATLAB program for Check Problem 3 simulation.
% An overflowing tank problem.
% Dr. K.Levien Chemical Engineering 361 Jan. 13, 2012
% Oregon State University
%
format short, format compact
global A Fi K
% variables shared with function(s)
%
A = 1.00; % cylindrical tank cross sectional area (m^2)
h0 = 1.00; % initial height of liquid in tank (m)
K = 2.00; % proportional constant for flow out of tank = Kh (m^2/min)
Fi = input('Enter constant feedflow at t >= 0 (m^3)/min (5.0) : ');
tend = input('Enter duration of simulation in minutes (1.0) : ');
%
tdata = [0 5 10 20 25]';
% time in experiment in seconds
hdata = [1 1.25 1.4 1.7 1.85]'; % height in experiment in meters
%
t = linspace(0,tend,201)'; % create column vector of time values
h = Fi/K + (h0Fi/K)*exp(K*t/A); % calculate corresponding heights
tsec = t*60; % convert time in minutes into time in seconds
%
figure % open a Figure Window
plot(tsec,h,'b'), grid
hold on
plot(tdata,hdata,'or')
xlabel('Time (s)'), ylabel('Height (m)')
title('Sample Prob. 3 Model + Data, Benny Beaver CHE 361 1/13/12')
legend('Model','Data','Location','Best')
hold off
%
tspan = [0 tend]; % specifies beginning and ending time values
[t_ode,h_ode] = ode45(@dhdt,tspan,h0); % uses ode45 to solve the ODE
%
figure
plot(tsec,h,'b','LineWidth',2), grid
hold on
plot(t_ode*60,h_ode,'ro') % t_ode was in minutes, so convert to s to plot
xlabel('Time (s)'), ylabel('Height (m)')
title('Sample Prob. 3 Analytical + Numerical, Benny Beaver CHE 361 1/13/12')
legend('Analytical','Numerical','Location','Best')
hold off Note first line of function file = dhdt.m
function [RHS] = dhdt(t,h)
%dhdt.m = RHS function to calculate derivative of Prob3 tank
% Dr. K.Levien Chemical Engineering 361 Jan. 13, 2012
% Oregon State University
%
global A Fi K
% variables shared with MAIN program
%
RHS = (Fi  K*h)/A; ...
View
Full
Document
This note was uploaded on 03/01/2012 for the course CHE 361 taught by Professor Staff during the Winter '08 term at Oregon State.
 Winter '08
 Staff

Click to edit the document details