tf = 72; % h
h0 = 2; % m
% Solve the DAE
opts = odeset('Reltol',1e-5);
[t,h] = ode45(@tankf1,[0 tf],h0,opts);
% Plot the graph
plot(t,h)
title('Tank model solved by Method 1')
xlabel('Time (h)')
ylabel('Level (m)')
grid on
function dh_dt = tankf1(t,h)
% Tank model parameters and constants
A = 12.5; % m^2
CV1 = 3.41; % (m^3/h)/kPa^0.5
CV2 = 3.41; % (m^3/h)/kPa^0.5
P0 = 100; % kPa
P3 = 100; % kPa
rho = 1000; % kg/m^3
g = 9.81; % m/s^2
P1 = P1function(t); % kPa
dh_dt = (1/A)*((CV1*sqrt(P1-P0-rho*g*h/1000))-(CV2*sqrt(P0+rho*g*h/1000-P3)));
% Set extra parameters
end
function P1 = P1function(t)
% P1FUNCTION: Upstream pressure (kPa) as a function of time (h)
if t<=10
P1 = 139.5;
elseif t<=11
P1 = 139.5 + 60.5*(t-10);
elseif t<=40
P1 = 200;
elseif t<=55
P1 = 200 - 30*(t-40)/15;
else % t>55
P1 = 170;
end
end
