Figure 7.9.1 shows a spherical tank for storing water. The tank is filled through a hole in the top and drained through a hole in the bottom. The following model for the liquid height h is developed in the Problems 7.38 and 7.39.
\pi(2R h-h^{2})\frac{d h}{d t}=-C_{d}A_{o}\sqrt{2g h} (1)
For water, {{C}}_{d} = 0.6 is a common value.
Use MATLAB to solve this equation to determine how long it will take for the tank to empty if the initial height is 9 ft. The tank has a radius of R = 5 ft and has a 1-in.-diameter hole in the bottom. Use g = 32.2 ft/sec². Discuss how to check the solution.
With C_{d}=0.6,R=5,g=32.2,\mathrm{and}\,A_{o}=\pi(1/24)^{2}, equation (1) becomes
{\frac{d h}{d t}}=-{\frac{0.0334{\sqrt{h}}}{10h-h^{2}}} (2)
We can use our physical insight to guard against grossly incorrect results. We can also check the preceding expression for dh/dt for singularities. The denominator does not become zero unless h = 0 or h = 10, which correspond to a completely empty and a completely full tank. So we will avoid singularities if 0 < h(0) < 10.
We can use the following approximation to estimate the time to empty. Replace h on the right side of equation (2) with its average value, namely, (9 − 0)/2 = 4.5 feet. This gives dh/dt = −0.00286, whose solution is h(t) = h(0) − 0.00286t = 9 − 0.00286t. According to this equation, if h(0) = 9, the tank will be empty at t = 9/0.00286 = 3147 sec, or 52 min. We will use this value as a “reality check” on our answer.
The function file based on this equation is
function hdot = height(t,h)
hdot = -(0.0334*sqrt(h))/(10*h-h^2);
The file is called as follows, using the ode45 solver.
[t, h] = ode45(@height, [0, 2475], 9);
plot(t,h),xlabel('Time (sec)'),ylabel('Height (ft)')
The resulting plot is shown in Figure 7.9.2. Note how the height changes more rapidly when the tank is nearly full or nearly empty. This is to be expected because of the effects of the tank’s shape. The tank empties in 2475 sec, or 41 min. This value is not grossly different from our rough estimate of 52 min, so we should feel comfortable accepting the numerical results.
The value of the final time of 2475 sec was found by increasing the final time until the plot showed that the height became zero. You could use a while loop to do this, by increasing the final time in the loop while calling ode45 repeatedly.