%
% OILTANK.M Modeling of Flow from a Damaged Oil Tank
%
% This demo models and evaluates the flow of oil from a horizontal tank. This
% simulation requires several numerical procedures including:
% 1. numerical integration of a fairly complicated function (QUADL)
% 2. spline fitting, root finding, and piecewise polynomial interpolation (SPLINE,
% FZERO, and PPVAL)
% 3. numerical solution of a first order differential equation (ODE23)
%
% The cross sectional area of the horizontal cylindrical tank is given by
% A(h) = 2*int[x(y)dy] for 0 < y < h
% where x(y) = sqrt(Dy - y^2)
% This integral is performed numerically for several values of h using QUADL.
% This gives a relationship for A(h) and V(h) = L*A(h).
%
% The volume of oil in the tank, V(h) = L*A(h), is then fit to a set of piecewise
% cubic polynomials using SPLINE. Then, knowing the initial volume of oil, one can
% determine the initial height of oil in the tank, ho. This is done using FZERO
% to find the root of f(h) = V(h) - Vo = 0. PPVAL is also used here to perform the
% piecewise polynomial interpolation - evaluating V(h) for arbitrary h.
%
% With the initial height known, ODE23 is used to numerically integrate the
% differential equation that represents a mass balance on the oil in the tank
% [written in terms of h(t)]. The equation of interest here is
% dh/dt = -(Aout/2L)*sqrt(2g/(D-h)) with h(0) = ho
%
% Finally, now that h(t) is known, PPVAL is used again to determine V(t) = V(h(t))
% and to evaluate the amount of oil that is spilled from the tank, Vs(t) = Vo-V(t).
%
% This main m-file also refers to the following function files:
% oiltanka.m - function file for QUADL routine
% oiltankb.m - function file for FZERO routine
% oiltankc.m - function file for ODE23 routine
%
% File prepared by J. R. White, UMass-Lowell (July 2003)
%
%
% getting started
clear all, close all, nfig = 0;
global PP Vo
%
% specify problem data and convert to consistent units
DIAM = 2.1; % tank diameter (m)
LENGTH = 5.6; % tank length (m)
Vo = 5000; % initial volume of oil (gallons)
cg2m3 = 3.785*0.001; % conversion factor (gallons to m^3)
Vo = Vo*cg2m3; % initial volume of oil (gallons converted to m^3)
Dout = 1.0; % diameter of opening (inches)
cin2m = 2.54/100; % conversion factor (inches to m)
Dout = Dout*cin2m; % diameter of opening (inches converted to m)
AOUT = pi*Dout*Dout/4; % area of opening (m^2)
GACCEL = 9.8; % gravitational accel (m^2/sec)
Tresp = 50; % response time (minutes)
Tresp = Tresp*60; % response time (minutes converted to seconds)
hmin = 0; hmax = DIAM; % min & max value for height of oil (m)
%
% find volume of oil versus height of oil
Nh1 = 11; hd1 = linspace(hmin,hmax,Nh1); Vh1 = zeros(size(hd1));
for i = 2:Nh1
Vh1(i) = 2*LENGTH*quadl('oiltanka',hmin,hd1(i),[],[],DIAM);
end
%
% find piecewise cubic polynomial fit to discrete data
PP = spline(hd1,Vh1);
%
% plot oil volume versus oil height (over fine grid for smooth curve)
Nh2 = 51; hd2 = linspace(hmin,hmax,Nh2); Vh2 = ppval(PP,hd2);
Vhg = Vh2/cg2m3; % oil volume in gallons
nfig = nfig+1; figure(nfig)
subplot(2,1,1),plot(hd2,Vhg,'b-','LineWidth',2),grid
range = axis; range(2)= hmax; axis(range);
title('OilTank: Volume of Oil Versus Oil Height')
xlabel('Oil Height (m)'),ylabel('Oil Volume (gallons)')
%
% knowing the initial oil volume, Vo, determine the initial height, ho
ho = fzero('oiltankb',hmax);
disp('Initial oil height (m) = '), ho
%
% now, with ho known, solve the differential balance equation for h(t)
to = 0; tf = 120*60; tol = 0.00001;
options = odeset('RelTol',tol);
[t,h] = ode23('oiltankc',[to tf],ho,options,AOUT,LENGTH,GACCEL,DIAM);
%
% plot h(t)
subplot(2,2,3),plot(t/60,h,'b-','LineWidth',2),grid
range = axis; range(2)= 120; axis(range);
title('Height of Oil Versus Time')
xlabel('Time (min)'),ylabel('Oil Height (m)')
%
% now find and plot amount of oil spilled
Vs = Vo-ppval(PP,h);
Vsg = Vs/cg2m3; % spilled oil volume in gallons
subplot(2,2,4),plot(t/60,Vsg,'b-','LineWidth',2),grid
range = axis; range(2)= 120; axis(range);
title('Volume of Oil Spilled Versus Time')
xlabel('Time (min)'),ylabel('Oil Spilled (gallons)')
%
range = axis;
xt = range(1) + 0.40*(range(2)-range(1));
yt1 = range(3) + 0.20*(range(4)-range(3));
yt2 = range(3) + 0.08*(range(4)-range(3));
text(xt,yt1,['V_o = ',num2str(Vo/cg2m3),' gallons'])
text(xt,yt2,['D_{out} = ',num2str(Dout/cin2m),' inches'])
%
% end simulation