How to solve a Bernoulli Equation
17 views (last 30 days)
Show older comments
I have to solve this equation:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/869540/image.png)
It has to start from known initial state and simulating forward to predetermined end point displaying output of all flow stages.
I have translated it into matlab code:
heightA+Id*DeltaL+turbulence*v^2/2*gravityG == heightA+turbulence*v^2/2*gravityG+Ifsr*DeltaL
Now I wonder should I embed it into ode45 or fsolve function?
12 Comments
William Rose
on 23 Jan 2022
Thank you for the picture of the dam with a pipe passing through it.
I think we are having difficulty due to translation issues. I am sorry that I only speak English.
I do not understand which values you already know, and which values you are trying to calculate, and with what equations.
You used the words Segment and Section. What is a Segment and whtat is a Section, in this problem?
Is the pipe circular cross section, with radius R and length L? Or length
?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870915/image.png)
You have introduced a new variable, n, in the equation for
. What is n? What are the units for n? The reason I am asking is that I want to understand the units for
. What do you think the units are for
?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870920/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870925/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870930/image.png)
You say energy
, but the units on the right hand side and length, not energy.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870935/image.png)
I keep asking about units, because if the units do not match then the equaito cannot be correct.
In the diagram of the dam, I woudl say p1=0, v1=0 (pressure and velocity at the surface).
I also belive that P2=0 asince it appears to be at the outflow of the pipe. The pressure at the inlet is obviously important. Let us cal the inlet pressure P3. Then
, where
and ρ = density of the fluid.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870940/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870945/image.png)
How is R computed?
The table include values for P. Does that refer to P1 or P2 or some other pressure? How was P computed?
You have given equaitons for
and
, but your equation has
on the left hand side. What is the equation for
?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870950/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870955/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870960/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/870965/image.png)
Do you have an equation for pressur drop per unit legth of pipe? If you do, what is that equation?
Accepted Answer
William Rose
on 24 Jan 2022
We have the following equation for tubulent flow in a rectangular channel, or canal:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872485/image.jpeg)
I assume that h1 and h2 are the depth of water in the channel, and v1 and v2 are the mean velocities, at two points. I assume that Id, Ifsr, Delta L, g, and alpha are known constants.
If we know h1 and v1 , then there will be an infinite number of combinations of h2,v2 which will satisfy the equation. (One equation, two unknowns.) Conservation of flow provides the other equation needed to obtain a unique solution for h2,v2. If the channel width is B at both points, then the flows are
and
. We assume the flow is the same at both points. Therefore
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872490/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872495/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872500/image.png)
Now we solve:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872505/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872510/image.png)
where
and
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872515/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872520/image.png)
We solve Eq.3 for h2. Then we use eq. 2 to solve for v2:
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872525/image.png)
Example:
Use values from the comments above: Id=0.0029, Ifsr=0.0092, alpha=1.1, Delta L=13.95 m, g=9.81 m/s^2. Choose h1=1 m, which is less than 75% of h_n=1.94, computed earlier. Choose v1= 10 m/s, so that flow Q1=B*h1*v1=20 m^3/s (assuming B=2, as given), therefore Q1 is less than the maximim Q=20.96 in the spreadsheet. Let the initial guess for h2 be h1.
clear; %clear variables from previous activity
Id=0.0029; Ifsr=0.0092; alpha=1.1; DL=13.95; g=9.81; h1=1; v1=10;
A1=h1+(Id-Ifsr)*DL+alpha*v1^2;
A2=alpha*(h1*v1)^2;
h20=h1; %initial guess for h2
h2=fsolve(@(h2)A1-h2-A2/h2^2,h20);
v2=v1*h1/h2;
fprintf('h1=%.4f m, v1=%.3f m/s, h2=%.4f m, v2=%.3f m/s\n',h1,v1,h2,v2);
Let us use a for loop to vary h1 from 0.6 to 1.4 m in step of 0.2 m:
clear; %clear previous values from memory
Id=0.0029; Ifsr=0.0092; alpha=1.1; DL=13.95; g=9.81; v1=10;
for h1=0.6:0.2:1.4 %depth at position 1 (m)
A1=h1+(Id-Ifsr)*DL+alpha*v1^2;
A2=alpha*(h1*v1)^2;
h20=h1; %initial guess for h2
options=optimoptions('fsolve','Display','off'); %option: turn off fsolve display
h2=fsolve(@(h2)A1-h2-A2/h2^2,h20,options);
v2=v1*h1/h2;
fprintf('h1=%.4f m, v1=%.3f m/s, h2=%.4f m, v2=%.3f m/s\n',h1,v1,h2,v2);
end
You can try varying other quantities.
The next example is like the previous one, but h1, h2, and v2 are stored in one-dimensional arrays. The resuts are displayed graphically.
clear; %clear previous values from memory
Id=0.0029; Ifsr=0.0092; alpha=1.1; DL=13.95; g=9.81; v1=10;
h1=0.6:0.2:1.4; %vector h1
h2=zeros(size(h1)); %allocate memory for vector h2, with the same size as h1
v2=zeros(size(h1)); %allocate memory for vector v2
for i=1:length(h1)
A1=h1(i)+(Id-Ifsr)*DL+alpha*v1^2;
A2=alpha*(h1(i)*v1)^2;
h20=h1(i); %initial guess for h2
options=optimoptions('fsolve','Display','off'); %option: turn off fsolve display
h2(i)=fsolve(@(h2)A1-h2-A2/h2^2,h20,options);
v2(i)=v1*h1(i)/h2(i);
end
figure; %create figure for plotting
subplot(2,1,1); %divide figure into 2 rows by 1 column, and plot in row 1
plot(h1,h2,'-rx'); %plot h1 versus h2, with red "x" symbols and connecting lines
xlabel('h1 (m)'); ylabel('h2 (m)'); title('h1 versus h2'); grid on
subplot(2,1,2); %divide figure into 2 rows by 1 column, and plot in row 2
plot(h1,v2,'-rx'); %plot h1 versus v2, with red "x" symbols and connecting lines
xlabel('h1 (m)'); ylabel('v2 (m/s)'); title('h1 versus v2'); grid on
We made two plots in one figure by using the subplot() command.
2 Comments
More Answers (1)
William Rose
on 23 Jan 2022
The pipe under the dam appears to be horizontal. I assume that R is the pipe diameter and that
is its length and that it has a circular cross section. Then the pressure loss along the pipe is Pin-Pout.
.
since it is the pressure at the top of the surface above the dam.
since it is the pipe outflow. Therefore
. is the pressure drop, or head loss, from the beginning to the end of the pipe. The Darcy Weisbach equation says
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871340/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871345/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871350/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871355/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871360/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871365/image.png)
where v=average velocity, ρ=density, and
=Darcy friction factor, which is dimensionless. The Darcy fricton factor complicated. If the flow is laminar, then
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871370/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871375/image.png)
where μ=dynamic viscosity, and the other quantities are defined above.
If the flow is turbulent, then the Darcy friction factor
depends on the velocity and the pipe diameter and pipe roughness. Is this all part of your study? Are you supposed to know abotu this to solve the problem?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871380/image.png)
In your table, A does not equal
, therefore A is not the cross secitnal area of the pipe. IN your table,
, which seems strange to me, In your tble, P is nt proportional to h. Therefore P is not the pressure at the inlet to the pipe. This also surprises me.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871385/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/871390/image.png)
4 Comments
William Rose
on 24 Jan 2022
You can compute hn exactly:
satisfies the equation
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872155/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872160/image.jpeg)
In the equation above, Q, B, n, and
are all known constants. Therefore we can solve for
in Matlab or Excel. In my version of Excel, I have the Analysis Tool Pack installed. Therefore I have the "Solver" utility available in Excel. I create a cell in Excel that has the formula above. I tell the Solver utility to make the formula be equal to zero, by adjusting the value in the cell containing
. The initial value of
is 2.0. Excel immediately changes the value in the cell containing
to 1.9366. The image below shows my Excel workbook before I run the solver (left) and after I run the solver (right).
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872165/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872170/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872175/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872180/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872185/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872190/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872195/image.jpeg)
In Matlab, you can solve for
as follows:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872200/image.png)
B=2; Id=.01; n=.014; Q=20.96; hn0=2;
hn=fsolve(@(hn)Q*n/Id^.5-(B*hn)^(5/3)/(B+2*hn)^(2/3),hn0);
fprintf('hn=%.4f\n',hn);
The small difference in the value of
found by Excel and Matlab is because they have different default tolerances for finding a solution.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872205/image.png)
___________
The Excel sreadsheet which you posted has Q=20,96, but in another post, you said "Q = flow capacity (calculated earlier to be 7,860 m3/s)".
The Excel spreadsheet which you posted has Id=0,01, but in another post, you said "Id = slope gradient of pipe (calculated to be 0,029%)", and an even earlier post , you showed an image which included Id=0.0029.
When you solve for
, use the correct values for Q and Id.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872210/image.png)
Now that we have a way to compute
, does that help you solve your problem?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/872215/image.png)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!