### A Pluto.jl notebook ### # v0.17.5

using Markdown using InteractiveUtils

# c110b070-742d-11ec-20d4-d3aa4debfe54 using CairoMakie

# f0c4e102-16b3-476d-a398-70df362aaf0f set_theme!(theme_light()); update_theme!(fontsize=20, linewidth=3)

# 6b27e69f-5394-4e8b-bbd9-ec1ea169900b md" ## a conical tank emptying of liquid

derive (on pencil and paper) a dynamic model of the liquid level, $h=h(t)$ [m], in a conical tank as liquid autonomously flows out of it, through a pipe at the bottom with a small valve providing a narrow constriction for the outflow.


**tank geometry** is an inverted, right circular cone (see [here](https://en.wikipedia.org/wiki/Cone)). * right cone $\implies$ axis passes through the center of the base and is orthogonal to the base * circular cone $\implies$ base is a circle * inverted cone $\implies$ base is on the top * height of the cone is $H$ [m] * the radius of the circle forming the base is $R$ [m]

**initial condition**: * the tank is initially full to the very brim

**other assumptions**: * the liquid is of constant density $ ho$ [kg/m$^3$] * flow rate out of the tank is driven by hydrostatic pressure and is related to the liquid level as $c\sqrt{h}$ [m$^3$/s]

your dynamic model should be an ODE in $h=h(t)$ and involve *only* the variables $H$, $R$, $c$, and $h$.

!!! hint write a mass balance in differential form. look for two similar triangles.

my model: ```math \begin{equation} \frac{dh}{dt}=f(h) \end{equation} ``` where ```math \begin{equation} f(h):= \text{FILL IN} \end{equation} ``` "

# 87faaf32-c035-4578-9474-400f37943bb1 md" ## finite difference methods

**goal**: write code to find the numerical solution to your dynamic model, via a Forward Euler finite difference method, using the parameter settings below. "

# d7bc2e84-2c23-4497-995f-39ebde143803 begin H = 2.0 # m R = 0.4 # m c = 0.00175 # m ^ (5/2) / s. end

# a5b69862-2448-4225-bf4e-55e4756de6a9 md" declare the initial liquid level as a variable `h`. (the tank is initially completely full of liquid.) "

# f3c9aa61-fd22-482b-9808-5f423f23dbc9

# 626e5ea3-ff66-4455-bc05-aafa1725ba74 md" code up the function `f(h)`, defined in relation to your dynamic model for $h=h(t)$ via:

$\dfrac{dh}{dt} = f(h)$

the function $f(h)$ characterizes the dynamics of $h$. "

# b7599381-2f2c-4cdf-898b-214a91436597

# 9760b691-be94-4d8c-a5d0-56568b9810bd md" to set up our finite difference approximation, define the time step `t` to be 0.05 s. "

# 43e1b7c5-b71c-45c9-85c0-80eb3b9b837a

# ed761a64-1708-4b85-8cc6-cf83d51342dc md" we wish to simulate the model in the time interval $t \in [0, 3]$ min. the variable `t` is defined to be the final time **in seconds**, after we take the `n` time steps to arrive at it."

# 7993b7e5-3e02-4863-a610-3471008afcb8

# 3295e2fe-baea-4766-9210-f21f1ec435c1 md" compute and define `n`---the number of time steps we wish to take---in terms of `t` and `t` so that we can change the time step later and have `n` automatically update, accordingly.

!!! note it's important to make `n` an integer by casting it as an `Int` so we can later write a `for` loop like `for i = 1:n`. you can achieve this by rounding `n` that you computed to the nearest integer. e.g., see what `round(Int, 4.01)` does. "

# f27010b6-6f09-4217-9df9-3feb462cb488

# ee8abcf8-0565-48ff-8d86-9666961907ac md" discretize time by constructing an array `t` with the time points $t_i$. i.e. element $i$ of the array `t` should be equal to $t_i= \Delta t (i-1)$ for $i \in \{1, 2, ..., n+1\}$. we go up to $n+1$ because the first element of `t` should be zero. the last entry of the `t` array should be `t`.

!!! hint preallocate an array of zeros, loop through each entry and overwrite entries with the appropriate values. "

# bfdc82ed-b34d-412e-afe7-3529aaeaecd7

# 84960114-4700-483b-82e1-92711acf8e74 md" we will store the approximations to the solution, $h_i \approx h(t_i)$'s, in an array `h`. pre-allocate this array with zeros, of the appropriate size to correspond to `t`, so $n+1$ elements. we'll overwrite these zeros later.

i.e. entry $i$ of the array `h` will hold our approximation to $h(t_{i-1})$ for $i \in \{1, 2, ..., n+1\}$.

the first element of `h` should contain the initial condition, `h`. assign entry `1` of `h` to be the initial condition `h`. "

# e57b612d-605d-4cbd-bcf8-e1a89ec9259a

# d284fba1-56f7-42ec-a615-48ddeb4ff903 md" implement the forward Euler method to approximate the solution to this ODE and fill in the rest of `h`. march ahead in time for `n` steps, via a `for` loop. use all of the following that you defined earlier: * the `h` array * the `t` array * `t` * `f(h)`

!!! note you should see a `DomainError`. the domain error comes from passing a negative value of $h$ into a square root function contained in $f(h)$. this happens at the time point when the tank is fully emptied and the dynamics are \"over\", but there is a little bit of error emanating from the finite difference approximation. to fix it, *when the tank is finally empty* inside the for loop, stop the time stepping since we know $h(t)=0$ for this time point and beyond. "

# b98527cd-2424-417a-9751-cedfafb0fc14

# fe939caa-a5fa-41c6-9a27-3dd0c82825cb md" plot your numerical solution $h(t)$, via plotting `h` vs `t`.

include an `xlabel` and `ylabel`. be sure to indicate the units. this is proper data visualization practice.

!!! hint if you see a little blip in your plot at the point where the tank empties of liquid, you need to assign `h[i]` at that time point to be zero during your forward Euler routine. "

# 59e9b8ca-a72a-4929-8d5f-1402eff23140

# 10e851c4-0d3c-472f-a585-4be551d2871c md" does the liquid level drop faster early in the emptying process or later? why is that? think about *both* the tank geometry and the hydrostatic pressure at the bottom of the tank that is driving flow out of it. they are _competing_ effect, with one dominating!

[... your answer here ...] "

# bb4755b6-a642-4db1-aa52-798bbeeeea3d md" derive an analytical solution to the differential equation for $h(t)$ by separating variables and integrating.

```math \begin{equation} h(t) = ??? \end{equation} ``` "

# 7e05bec1-3c0b-4d51-81cc-f1b64f50b525 md" use the analytical solution to determine the time $t_e$ at which the tank empties, i.e. the time that the liquid level first reaches zero: ```math \begin{equation} t_e= \arg\min_{t : \hspace{0.4em} h(t)=0} h(t) \end{equation} ```

compute $t_e$ and assign it as a variable `t` below. is it consistent with your plot above (it should be!)? "

# 0d5304e7-74b4-4ff8-82d7-8d6f67419ecf

# 4b02b8d8-84cb-464c-8619-48409f90047e md" plot the analytical solution on top of the numerical solution to ensure they match (they should). since two curves are on the same plot panel, (i) include a legend and (ii) make one of the curves dashed.

!!! note the analytical solution has a $\sqrt(\cdot)$ function in it and is valid only for $tt_e$, $h(t)=0$. thus when you write a function `h_true(t)` to plot the analytical solution, you must write an `if` statement to check if the `t` is less than or equal to the time at which the tank empties. "

Julia problem