Question 23.CS.5: The Roman natural philosopher, Pliny the Elder, purportedly ...
The Roman natural philosopher, Pliny the Elder, purportedly had an intermittent fountain in his garden. As in Fig. 23.12, water enters a cylindrical tank at a constant flow rate Q_{in} and fills until the water reaches y_{high}. At this point, water siphons out of the tank through a circular discharge pipe, producing a fountain at the pipe’s exit. The fountain runs until the water level decreases to y_{low}, whereupon the siphon fills with air and the fountain stops. The cycle then repeats as the tank fills until the water reaches y_{high}, and the fountain flows again.
When the siphon is running, the outflow Q_{out} can be computed with the following formula based on Torricelli’s law:
Q_{\text {out }}=C \sqrt{2 g y \pi r^2} (23.28)
Neglecting the volume of water in the pipe, compute and plot the level of the water in the tank as a function of time over 100 seconds. Assume an initial condition of an empty tank y(0) = 0, and employ the following parameters for your computation:
\begin{array}{lll} R_T=0.05 \mathrm{~m} & r=0.007 \mathrm{~m} & y_{\text {low }}=0.025 \mathrm{~m} \\ y_{\text {high }}=0.1 \mathrm{~m} & C=0.6 & g=9.81 \mathrm{~m} / \mathrm{s}^2 \\ Q_{\mathrm{in}}=50 \times 10^{-6} \mathrm{~m}^3 / \mathrm{s} && \\ \end{array}
When the fountain is running, the rate of change in the tank’s volume V (m³) is determined by a simple balance of inflow minus the outflow:
\frac{d V}{d t}=Q_{\text {in }}-Q_{\text {out }} (23.29)
where V = volume (m³). Because the tank is cylindrical, V = π R_t^2 y. Substituting this relationship along with Eq. (23.28) into Eq. (23.29) gives
\frac{d y}{d t}=\frac{Q_{\mathrm{in}}-C \sqrt{2 g y} \pi r^2}{\pi R_t^2} (23.30)
When the fountain is not running, the second term in the numerator goes to zero. We can incorporate this mechanism in the model by introducing a new dimensionless variable siphon that equals zero when the fountain is off and equals one when it is flowing:
\frac{d y}{d t}=\frac{Q_{\text {in }}-\text { siphon } \times C \sqrt{2 g y} \pi r^2}{\pi R_t^2} (23.31)
In the present context, siphon can be thought of as a switch that turns the fountain off and on. Such two-state variables are called Boolean or logical variables, where zero is equivalent to false and one is equivalent to true.
Next we must relate siphon to the dependent variable y. First, siphon is set to zero whenever the level falls below y_{low}. Conversely, siphon is set to one whenever the level rises above y_{high}. The following M-file function follows this logic in computing the derivative:
function dy = Plinyode(t,y)
global siphon
Rt = 0.05; r = 0.007; yhi = 0.1; ylo = 0.025;
C = 0.6; g = 9.81; Qin = 0.00005;
if y(1) <= ylo
siphon = 0;
elseif y(1) >= yhi
siphon = 1;
end
Qout = siphon * C * sqrt(2 * g * y(1)) * pi * r ^ 2;
dy = (Qin − Qout) / (pi * Rt ^ 2);
Notice that because its value must be maintained between function calls, siphon is declared as a global variable. Although the use of global variables is not encouraged (particularly in larger programs), it is useful in the present context.
The following script employs the built-in ode45 function to integrate Plinyode and generate a plot of the solution:
global siphon
siphon = 0;
tspan = [0 100]; y0 = 0;
[tp,yp] = ode45(@ Plinyode,tspan,y0);
plot(tp,yp)
xlabel('time, (s)')
ylabel('water level in tank, (m)')
As shown in Fig. 23.13, the result is clearly incorrect. Except for the original filling period, the level seems to start emptying prior to reaching y_{high}. Similarly, when it is draining, the siphon shuts off well before the level drops to y_{low}.
At this point, suspecting that the problem demands more firepower than the trusty ode45 routine, you might be tempted to use one of the other MATLAB ODE solvers such as ode23s or ode23tb. But if you did, you would discover that although these routines yield somewhat different results, they would still generate incorrect solutions.
The difficulty arises because the ODE is discontinuous at the point that the siphon switches on or off. For example, as the tank is filling, the derivative is dependent only on the constant inflow and for the present parameters has a constant value of 6.366 × 10^{−3} m/s.
However, as soon as the level reaches y_{high}, the outflow kicks in and the derivative abruptly drops to −1.013 × 10^{−2} m/s. Although the adaptive step-size routines used by MATLAB work marvelously for many problems, they often get heartburn when dealing with such discontinuities.
Because they infer the behavior of the solution by comparing the results of different steps, a discontinuity represents something akin to stepping into a deep pothole on a dark street.
At this point, your first inclination might be to just give up. After all, if it’s too hard for MATLAB, no reasonable person could expect you to come up with a solution. Because professional engineers and scientists rarely get away with such excuses, your only recourse is to develop a remedy based on your knowledge of numerical methods.
Because the problem results from adaptively stepping across a discontinuity, you might revert to a simpler approach and use a constant, small step size. If you think about it, that’s precisely the approach you would take if you were traversing a dark, pothole-filled street. We can implement this solution strategy by merely replacing ode45 with the constant- step rk4sys function from Chap. 22 (Fig. 22.8). For the script outlined above, the fourth line would be formulated as
[tp,yp] = rk4sys(@ Plinyode,tspan,y0,0.0625);
As in Fig. 23.14, the solution now evolves as expected. The tank fills to y_{high} and then empties until it reaches y_{low}, when the cycle repeats.
There are a two take-home messages that can be gleaned from this case study. First, although it’s human nature to think the opposite, simpler is sometimes better. After all, to paraphrase Einstein, “Everything should be as simple as possible, but no simpler.” Second, you should never blindly believe every result generated by the computer. You’ve probably heard the old chestnut, “garbage in, garbage out” in reference to the impact of data quality on the validity of computer output. Unfortunately, some individuals think that regardless of what went in (the data) and what’s going on inside (the algorithm), it’s always “gospel out.” Situations like the one depicted in Fig. 23.13 are particularly dangerous—that is, although the out put is incorrect, it’s not obviously wrong. That is, the simulation does not go unstable or yield negative levels. In fact, the solution moves up and down in the manner of an intermittent fountain, albeit incorrectly.
Hopefully, this case study illustrates that even a great piece of software such as MATLAB is not foolproof. Hence, sophisticated engineers and scientists always examine numerical output with a healthy skepticism based on their considerable experience and
knowledge of the problems they are solving.
