Fluid dynamics problem: Equalize the volume of water in tanks at differing heights

[Juanda]

I was wrong about this, it gets messy.

We need to be using "Unsteady -Bernoulli"¹ which says in general:

$$ P_1 + \rho g z_1 + \frac{1}{2}\rho v_1^2 = \int_1^2 \rho \frac{\partial v_s}{\partial t} ds + P_2 + \rho g z_2 + \frac{1}{2}\rho v_2^2 \tag{1}$$

Messy is an understatement for me. I had never seen this version of Bernoulli. Fluids are hard... Anyways, this just proves this problem is a chance to learn some new things for me.

After reading the PDF you shared here's my attempt to show the equations before moving into coding where traceability will be much harder. I want to ensure we agree on the methodology before crunching numbers and plotting graphs. I'll express everything in terms of flow instead of velocities which I believe is more closely related to what we're interested in.
I'll need to redraw the diagram and rename a few things. I hate doing it at this stage but I need to do it to keep track of the relevant points and variables assigned to them.

Unsteady Bernoulli $(1)$ from the PDF you shared.

Using $(1)$ between the points $P$ and $D$ gives $(2)$.
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds +P_D + \frac{1}{2}\rho v_D^2+\rho g (H+z_D) - P_P-\frac{1}{2}\rho v_P^2 - \rho g z_P = 0 \tag{2}$$

Taking into account the height of $P$ is zero and writing the velocities in terms of flow yields $(3)$
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds +P_D + \frac{1}{2}\rho \left ( \frac{Q_2}{S_T} \right ) ^2+\rho g (H+z_D) - P_P-\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2 = 0 \tag{3}$$

Considering ideal gas law and constant temperature then the pressure of the air inside tank 2 can be expressed in relation to its water level. That it'll be equation $4$.
$$P_DV_2=m_2RT \ \left \{ V_2=S_T(H_T-z_D) \right \} \rightarrow P_D= \frac{m_2RT}{S_T(H_T-z_D)} \tag{4}$$

The water level $z_D$ is related to the water flow $(5)$.
$$z_D = z_{D_0}+\frac{Q_2}{S_T}t \tag{5}$$

Combining $(4)$ and $(5)$, $(6)$ is obtained.
$$P_D= \frac{m_2RT}{S_T(H_T-\left (z_{D_0}+\frac{Q_2}{S_T}t \right ))} \tag{6}$$

The equation $(5)$ can also be used on $(3)$. The results from $(6)$ will be implemented too to give $(7)$.
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds +\left (\frac{m_2RT}{S_T(H_T-\left (z_{D_0}+\frac{Q_2}{S_T}t \right ))} \right ) + \frac{1}{2}\rho \left ( \frac{Q_2}{S_T} \right ) ^2+\rho g (H+\left (z_{D_0}+\frac{Q_2}{S_T}t \right )) - P_P-\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2 = 0 \tag{7}$$

Pressure at P ($P_P$) is an input value for the problem so I'll consider it known. Its expression is $(8)$.
$$P_P = P_{atm}+\rho g h_{pump}(Q)\rightarrow P_P = P_{atm}+ \rho g \left ( \alpha - \beta Q^2 \right ) \tag{8}$$

Using the information from $(8)$ into $(7)$ gives $(9)$.
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds +\left (\frac{m_2RT}{S_T(H_T-\left (z_{D_0}+\frac{Q_2}{S_T}t \right ))} \right ) + \frac{1}{2}\rho \left ( \frac{Q_2}{S_T} \right ) ^2+\rho g (H+\left (z_{D_0}+\frac{Q_2}{S_T}t \right )) - \left (P_{atm}+ \rho g \left ( \alpha - \beta Q^2 \right ) \right )-\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2 = 0 \tag{9}$$

Finally, it is necessary to evaluate the path integral in $(9)$ that goes from $P$ to $D$. Before doing it I want to clarify the change of variable to have it in terms of the flow instead of the velocity $(10)$.
$$\rho\frac{\partial v}{\partial t}ds \rightarrow \left \{ v=\frac{Q}{S} \right \} \rightarrow \frac{\rho}{S}\frac{\partial Q}{\partial t}ds \tag{10}$$

Here is something that mathematically might not be expressed correctly because I'll go from a general variable $Q$ between the points $P$ and $D$ to the specific flows $Q_2$ and $Q$ in the pipes. I'll write the "generic" flow with $Q_g$ and the generic section with $S_g$ to be able to tell the difference. This will be the equation $(11)$.
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds = \int_{P}^{D} \frac{\rho}{S_g}\frac{\partial Q_g}{\partial t}ds \tag{11}$$

According to the document from MIT, that integral from $P$ to $D$ can be simplified to the integral from $P$ to $F$ because $S_T>>S_P$ which only considers what's going on in the pipes. In this case, it is necessary to note the distinction in the flow from $P$ to $A$ and from $B$ to $F$ because the velocities/flow rates will not be the same. Equation $(12)$.
$$\int_{P}^{D} \frac{\rho}{S_g}\frac{\partial Q_g}{\partial t}ds = \int_{P}^{A} \frac{\rho}{S_p}\frac{\partial Q}{\partial t}ds+\int_{B}^{F} \frac{\rho}{S_p}\frac{\partial Q_2}{\partial t}ds=\frac{\rho}{S_p}\left ( \frac{\partial Q}{\partial t}L_{PA} + \frac{\partial Q_2}{\partial t}(H+L_{JF})\right ) \tag{12}$$

The expression from $(12)$ can be plugged in $(9)$ to obtain the final expression of Unsteady Bernoulli from $(2)$ for Tank 2. Equation $(13)$.
$$\frac{\rho}{S_p}\left ( \frac{\partial Q}{\partial t}L_{PA} + \frac{\partial Q_2}{\partial t}(H+L_{JF})\right ) +\left (\frac{m_2RT}{S_T(H_T-\left (z_{D_0}+\frac{Q_2}{S_T}t \right ))} \right ) + \frac{1}{2}\rho \left ( \frac{Q_2}{S_T} \right ) ^2+\rho g (H+\left (z_{D_0}+\frac{Q_2}{S_T}t \right )) - \left (P_{atm}+ \rho g \left ( \alpha - \beta Q^2 \right ) \right )-\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2 = 0 \tag{13}$$

Following the same logic it is possible to obtain the expression for Unsteady Bernoulli applied to Tank 1.
Finally, with the conservation of mass $Q = Q_1+Q_2$, I feel there is enough information to solve the system and obtain the functions $Q(t)$, $Q_1(t)$ and $Q_2(t)$ although I'd love to have confirmation about it if you know it for sure.

If we are on the same page about the procedure then I'd be ready to put it into Python to try solving it which will probably be a challenge in itself because I don't have a lot of experience with it.

 

[erobz]

Messy is an understatement for me. I had never seen this version of Bernoulli. Fluids are hard... Anyways, this just proves this problem is a chance to learn some new things for me.

After reading the PDF you shared here's my attempt to show the equations before moving into coding where traceability will be much harder. I want to ensure we agree on the methodology before crunching numbers and plotting graphs. I'll express everything in terms of flow instead of velocities which I believe is more closely related to what we're interested in.
I'll need to redraw the diagram and rename a few things. I hate doing it at this stage but I need to do it to keep track of the relevant points and variables assigned to them.
View attachment 330727Unsteady Bernoulli $(1)$ from the PDF you shared.
View attachment 330728Using $(1)$ between the points $P$ and $D$ gives $(2)$.
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds +P_D + \frac{1}{2}\rho v_D^2+\rho g (H+z_D) - P_P-\frac{1}{2}\rho v_P^2 - \rho g z_P = 0 \tag{2}$$

Taking into account the height of $P$ is zero and writing the velocities in terms of flow yields $(3)$
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds +P_D + \frac{1}{2}\rho \left ( \frac{Q_2}{S_T} \right ) ^2+\rho g (H+z_D) - P_P-\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2 = 0 \tag{3}$$

I believe the last term on the right should be:

$$\frac{1}{2}\rho \left ( \frac{Q_2}{S_p} \right ) ^2 $$

We are summing along a single streamline, what you have written is the kinetic energy along both streamlines.

The water level $z_D$ is related to the water flow $(5)$.
$$z_D = z_{D_0}+\frac{Q_2}{S_T}t \tag{5}$$

That is only if $Q,Q_2,Q_1$ are constant in time, I think it should be:

$$z_D = z_{D_0}+ \int \frac{Q_2}{S_T} dt \tag{5}$$

That propagates into a few other equations.

The equation $(5)$ can also be used on $(3)$. The results from $(6)$ will be implemented too to give $(7)$.
$$\int_{P}^{D}\rho\frac{\partial v}{\partial t}ds +\left (\frac{m_2RT}{S_T(H_T-\left (z_{D_0}+\frac{Q_2}{S_T}t \right ))} \right ) + \frac{1}{2}\rho \left ( \frac{Q_2}{S_T} \right ) ^2+\rho g (H+\left (z_{D_0}+\frac{Q_2}{S_T}t \right )) - P_P-\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2 = 0 \tag{7}$$

Pressure at P ($P_P$) is an input value for the problem so I'll consider it known. Its expression is $(8)$.
$$P_P = P_{atm}+\rho g h_{pump}(Q)\rightarrow P_P = P_{atm}+ \rho g \left ( \alpha - \beta Q^2 \right ) \tag{8}$$

You could just pretend $P_p$ is a constant for now for slightly more simple expressions, but that is up to you. It's going to be a beast to solve either way.

 

[Juanda]

I believe the last term on the right should be:

$$\frac{1}{2}\rho \left ( \frac{Q_2}{S_p} \right ) ^2 $$

We are summing along a single streamline, what you have written is the kinetic energy along both streamlines.

I'm integrating along a single streamline in this case. However, along that streamline, the velocity is not constant because the flow $Q$ must divide into $Q_1$ and $Q_2$. Since the cross-section of the pipes is constant, the velocity of the fluid after the separation must be smaller. At the junction, there's a discontinuity I'm choosing to ignore which is why I created the points $A$, $B$, and $C$ around it and "infinitely" close to it. That discontinuity would of course disappear with a more detailed analysis such as CFD but this already is complex enough the way it is.
At the point $P$, we have $v_P$ which in terms of the flow it is $Q/S_p$ ($S_p$ is the cross-section of the pipe) and at the point $D$ we have $v_D$ which in terms of the flow is $Q_2/S_T$ ($S_T$ is the cross-section of the tank).
Therefore, I think the equation $(3)$ is correct.

That is only if $Q,Q_2,Q_1$ are constant in time, I think it should be:

$$z_D = z_{D_0}+ \int \frac{Q_2}{S_T} dt \tag{5}$$

That propagates into a few other equations.

You're absolutely right here. I'll need to update that when I try to solve it.

You could just pretend $P_P$ is a constant for now for slightly more simple expressions, but that is up to you. It's going to be a beast to solve either way.

It's true I could. However, I'll have to discretize the rest of the problem anyways so maybe there's not much extra effort needed to describe the pump as a function of the flow. I just need to be careful when choosing the values of $\alpha$ and $\beta$ to make sure there'll be no backflow. I think the equations should predict that scenario without needing modifications (assuming the pipes are always full) but I'm interested in seeing first the "simple" case.

 

[erobz]

I'm integrating along a single streamline in this case. However, along that streamline, the velocity is not constant because the flow $Q$ must divide into $Q_1$ and $Q_2$. Since the cross-section of the pipes is constant, the velocity of the fluid after the separation must be smaller. At the junction, there's a discontinuity I'm choosing to ignore which is why I created the points $A$, $B$, and $C$ around it and "infinitely" close to it. That discontinuity would of course disappear with a more detailed analysis such as CFD but this already is complex enough the way it is.
At the point $P$, we have $v_P$ which in terms of the flow it is $Q/S_p$ ($S_p$ is the cross-section of the pipe) and at the point $D$ we have $v_D$ which in terms of the flow is $Q_2/S_T$ ($S_T$ is the cross-section of the tank).
Therefore, I think the equation $(3)$ is correct.

I don't think that a streamline begins until after we are outside the infinitesimal ball we're calling the junction. You should be going from $B$ to $D$, and $C$ to $E$ IMO.

Imagine a cylinder is rolling down a hill of mass $2M$, and it splits. On the path where the object is becomes the object, to where that particular object is going, I'm not going to uses kinetic energies involving $2M$ to describe a particular portion of its mechanical energy just after it became an object of mass $M$. I think the same is true here.

But I'm not going to argue further. You're modeling it, do whatever you want.

 

[Juanda]

I don't think that a streamline begins until after we are outside the infinitesimal ball we're calling the junction. You should be going from $B$ to $D$, and $C$ to $E$ IMO.

Why would you think that though? The flow starts at $P$ and goes to both $D$ and $E$ according to the diagram.

It'd be similar to something like this. Would you have trouble applying Bernoulli in this situation to find the state at A, B, C and D? With Bernoulli and mass conservation is possible. The only difference when compared with the problem we're discussing is that the flow isn't reuniting again after the junction but accumulating in the tanks.

Besides, in the problem with the tanks, it is necessary to use Bernoulli P→E and P→D so that the equations are related and can be solved simultaneosly in a system.

I'm not trying to argue for the sake of arguing but I'd prefer we both agree on how to model the problem before trying to implement it in the simulation.

 

[erobz]

If we go directly from A to D in the side problem you propose we learn nothing about the division of flow. We must visit the junctions in our equations.

We go from $P$ ( Which is just outside the pump) to the junction, then from the junction to each of the endpoints independently.

 

[erobz]

If we go directly from A to D in the side problem you propose we learn nothing about the division of flow. We must visit the junctions in our equations.

We go from $P$ ( Which is just outside the pump) to the junction, then from the junction to each of the endpoints independently.

Ok, after a quick example mimicking the problem at hand a bit more( one were the streamlines don't reunite), I convinced myself that we don't always have to visit nodes (the problem you chose as an instructive example, we do have to visit the nodes if we are to learn anything of the flow division - and we also would need to consider head loss to solve it). It could be that I just like visiting nodes out of habit, and may not be the most "streamlined" approach. I'm not 100% sure at the moment if what we are discussing isn't just going to give the same results in the end (but I believe they will)...have to think on that a bit. But I see now that (at least) your equations are good.

 

[erobz]

This is my take on it. on the incoming line where it splits you put a flow divider. there whole purpose is to split the flow to near 50/50.

We are basically just pontificating on the possibility of this at this point to try to learn something.

 

[Joe591]

We are basically just pontificating on the possibility of this at this point.

I'm not sure what you mean

 

[Juanda]

Ok, after a quick example mimicking the problem at hand a bit more( one were the streamlines don't reunite), I convinced myself that we don't always have to visit nodes (the problem you chose as an instructive example, we do have to visit the nodes if we are to learn anything of the flow division - and we also would need to consider head loss to solve it). It could be that I just like visiting nodes out of habit, and may not be the most "streamlined" approach. I'm not 100% sure at the moment if what we are discussing isn't just going to give the same results in the end (but I believe they will)...have to think on that a bit. But I see now that (at least) your equations are good.

I'm glad we agree. Knowing information about all the nodes within a streamline would be possible by applying Bernoulli to them but I chose the points I considered most relevant (exiting the pump and at the top of the tanks).
I'm then ready to try to write all this in code to solve it. Let's see if I can actually do it and get some physically possible results that make any sense.

 

[erobz]

@Joe591 Some say the flow just splits 50/50 without intervention, other suggest control ( like yourself). @Juanda and I are trying to tease something out of underlying physics for "fun".

 

[Joe591]

Some say the flow just splits 50/50 without intervention, other suggest control ( like yourself). @Juanda and I are trying to tease something out of underlying physics for "fun".

oh I see. yeah, I'm reasonably convinced it will require intervention. I tend to solve my problems with technology rather than math though, lol...

I must admit I don't have the faintest idea about the math that your using, haha

 

[erobz]

I'm glad we agree. Knowing information about all the nodes within a streamline would be possible by applying Bernoulli to them but I chose the points I considered most relevant (exiting the pump and at the top of the tanks).
I'm then ready to try to write all this in code to solve it. Let's see if I can actually do it and get some physically possible results that make any sense.

Can you write up the system of equations you are going to solve in final reduced form if you get a chance and talk about the numerical approach? I think this is going to be challenging system to solve (even programmatically).

 

[Juanda]

Can you write up the system of equations you are going to solve in final reduced form if you get a chance and talk about the numerical approach? I think this is going to be challenging system to solve (even programmatically).

I still owe you this. It's not that I have forgotten it but I know it's going to take at the very least a couple of hours being focussed to be able to do it and I can't get to it now. I'll try next week though.

 

[Juanda]

@erobz I think I got the system of equations you were talking about in post #48.
Equations will be referred to this diagram.

View attachment 330727

The dynamic state when filling it up can be defined by the equations $(1)$, $(2)$ and $(3)$.
Equation $(1)$ is unsteady Bernoullin from $P$ to $D$.
$$
\overbrace{\frac{\rho}{S_p}\left ( \frac{\partial Q}{\partial t}L_{PA} + \frac{\partial Q_2}{\partial t}(H+L_{JF})\right )}^{\int_{P}^{D} \frac{\rho}{S_g}\frac{\partial Q_g}{\partial t}ds}
+
\overbrace{\left (\frac{m_2RT}{S_T(H_T-\left (z_{D_0}+ \int \frac{Q_2}{S_T} dt \right ))} \right )}^{P_D}
+
\overbrace{\frac{1}{2}\rho \left ( \frac{Q_2}{S_T} \right ) ^2}^{\frac{1}{2}\rho v_D^2}
+
\overbrace{\rho g (H+\left (z_{D_0}+ \int \frac{Q_2}{S_T} dt \right ))}^{\rho g (H+z_D)}
-
\overbrace{\left (P_{atm}+ \rho g \left ( \alpha - \beta Q^2 \right ) \right )}^{P_P}
-
\overbrace{\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2}^{\frac{1}{2}\rho v_P^2}
=
0 \tag{1}
$$

Equation $(2)$ is unsteady Bernoulli from $P$ to $E$.
$$
\overbrace{\frac{\rho}{S_p}\left ( \frac{\partial Q}{\partial t}L_{PA} + \frac{\partial Q_1}{\partial t}L_{CG}\right )}^{\int_{P}^{E} \frac{\rho}{S_g}\frac{\partial Q_g}{\partial t}ds}
+
\overbrace{\left (\frac{m_1RT}{S_T(H_T-\left (z_{E_0}+ \int \frac{Q_1}{S_T} dt \right ))} \right )}^{P_E}
+
\overbrace{\frac{1}{2}\rho \left ( \frac{Q_1}{S_T} \right ) ^2}^{\frac{1}{2}\rho v_E^2}
+
\overbrace{\rho g \left (z_{E_0}+ \int \frac{Q_1}{S_T} dt \right )}^{\rho g z_E}
-
\overbrace{\left (P_{atm}+ \rho g \left ( \alpha - \beta Q^2 \right ) \right )}^{P_P}
-
\overbrace{\frac{1}{2}\rho \left ( \frac{Q}{S_p} \right ) ^2}^{\frac{1}{2}\rho v_P^2}
=
0 \tag{2}
$$

Equation $(3)$ is mass conservation.
$$Q=Q_1+Q_2 \tag{3}$$

From those 3 equations, I'm hoping it's possible to obtain the functions $Q(t)$, $Q_1(t)$ and $Q_3(t)$ from which it would also be possible to obtain the functions $z_E(t)$ and $z_D(t)$.

It'd be necessary to define the geometry of the problem (pipe lengths, tank's cross-section, tank's height, tanks' height difference), the amount of air in each tank, the system's temperature, pump's properties, gravity, fluid's density, and the initial conditions (initial heights and flows).

To set up the initial conditions $(t=0)$ I'd prefer to start the system from an equilibrium point. I feel like it's not necessary but I think it's convenient because it's more realistic and the oscillatory nature of the result will be suppressed (or so I think).

To find the initial conditions there are different ways of doing it. Previously in the thread (post #27), we discussed using Bernoulli from $E$ to $D$ and the conservation of mass. I'll change the method to using Bernoulli from $P$ to $E$ and Bernoulli from $P$ to $D$.
$$
\overbrace{\frac{m_2RT}{S_T(H_T-z_{D_0})}}^{P_{D_0}}
+
\rho g(H+z_{D_0})
=
P_{P_0} \tag{4}
$$
$$\overbrace{\frac{m_1RT}{S_T(H_T-z_{E_0})}}^{P_{E_0}}
+
\rho gz_{E_0}
=
P_{P_0} \tag{5}
$$

The change is because of two main reasons. First, the initially proposed method required to know the initial water in the system which is not realistic. It's easier to install a manometer at $P$ and read the pressure when the valve is closed. Secondly, doing it this way is more consistent with the previous set of equations $(1)$, $(2)$ and $(3)$ and the equations $(6)$ and $(7)$ that will be written later and describe the new equilibrium at $t \rightarrow \infty$ when the amount of mass that entered the system is unknown but the pressure at $P$ is known.

So, for the initial conditions, we have $(4)$ and $(5)$.
During the fill-up process the equations are $(1)$, $(2)$ and $(3)$.
Finally, for $t \rightarrow \infty$ we know the flow tends to 0 so we can use "normal" Bernoulli again just like it was done to set up the equations $(4)$ and $(5)$. Also, there will be no difference between having the inlet valve open or closed. That defines the equations $(6)$ and $(7)$
$$\overbrace{\frac{m_2RT}{S_T(H_T-z_{D_f})}}^{P_{D_f}}
+
\rho g(H+z_{D_f})
=
P_{atm}+\rho g \alpha \tag{6}
$$
$$\overbrace{\frac{m_1RT}{S_T(H_T-z_{E_f})}}^{P_{E_f}}
+
\rho gz_{E_f}
=
P_{atm}+\rho g \alpha \tag{7}
$$

It's useful to define some "checkpoints" to make sure the results are correct. Some of them are:

• Functions for flow and height should not oscillate since we started from an equilibrium point with zero velocity.
• We know that the evolution described by $(1)$, $(2)$ and $(3)$ when starting from $(4)$ and $(5)$ with zero velocity must tend towards the result from $(6)$ and $(7)$ with zero velocity as time goes on.
• The functions $z_E$ and $z_D$ can never be negative. It'd imply the tanks have no water and from the beginning, we said the pipes are always full so the very minimum is that the height is zero. From the previous points, it can also be understood that their slopes can't be negative either.
• Since I'll be working with absolute pressure. No pressures can be negative either.

I think from there the last point would be to try to code it to find a numerical solution. I'm not too experienced in this but by looking at the presented math I feel it should be possible.

 

[erobz]

To set up the initial conditions $(t=0)$ I'd prefer to start the system from an equilibrium point. I feel like it's not necessary but I think it's convenient because it's more realistic and the oscillatory nature of the result will be suppressed (or so I think).

I think starting from a state of non-equilibrium would pose some issues. The initial conditions would be no longer be $Q_1, Q_2,Q = 0$ nor would the initial flow accelerations be $\dot Q_1, \dot Q_2 = 0$