How to solve a Bernoulli Equation

66 views (last 30 days)
KarolN
KarolN on 21 Jan 2022
Commented: William Rose on 25 Jan 2022
I have to solve this equation:
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
KarolN
KarolN on 23 Jan 2022
@William Rose that's OK I believe difficulties are more due to my ignorance of that issue (I deal with Bernoulli first time in my life).
First I have to disclaim that picture was for general reference only to show how the process looks, it does not neccessarily mirror the set of equations from my instruction.
The other picture:
Shows, how technical drawing based on our calculations, looks like.
The GOAL of calculation is: you have to compute normal depth, and then using Bernoulli equation calculate the flow in such way that it won't sink our canal (it cannot fill the canal more than 75% of normal depth). Afterwards you assume final width, height and total length of canal.
The instruction is titled 'Calculation of Unpressurized Canal', so I guess there's no worthwhile pressure here.
Section, is where where the water level rises. Segment, is length from one section to next section, where some values change on the way. You assume the 'h' for each section as it is the water filling the canal at given section.
n - the coefficient of concrete roughness, constant and = 0,014
Id - slope gradient of pipe, based on average gradient fall in river bed = 0,029%
R - the hydraulic radius computed from:
Where A is cross- sectional area and P is perimeter of dampness computed from:
Where B = pipe width, hn = normal depth computed by trial and error from:
Where Q = flow capacity (calculated earlier to be 7,860 m3/s) Id = slope gradient of pipe (calculated to be 0,029%), A and R described earlier. The left side is constant.
As for your query about units I have to ask my professor. I only got set of equations. I go back to you as soon as I have them.

Sign in to comment.

Accepted Answer

William Rose
William Rose on 24 Jan 2022
We have the following equation for tubulent flow in a rectangular channel, or canal:
Eq.1.
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
Eq. 2
Now we solve:
Eq. 3
where and .
We solve Eq.3 for h2. Then we use eq. 2 to solve for v2: .
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);
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
v2=v1*h1/h2;
fprintf('h1=%.4f m, v1=%.3f m/s, h2=%.4f m, v2=%.3f m/s\n',h1,v1,h2,v2);
h1=1.0000 m, v1=10.000 m/s, h2=1.0004 m, v2=9.996 m/s
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
h1=0.6000 m, v1=10.000 m/s, h2=0.6002 m, v2=9.996 m/s h1=0.8000 m, v1=10.000 m/s, h2=0.8003 m, v2=9.996 m/s h1=1.0000 m, v1=10.000 m/s, h2=1.0004 m, v2=9.996 m/s h1=1.2000 m, v1=10.000 m/s, h2=1.2005 m, v2=9.996 m/s h1=1.4000 m, v1=10.000 m/s, h2=1.4006 m, v2=9.996 m/s
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
William Rose
William Rose on 25 Jan 2022
@KarolN, you're welcome. Good luck with your civil engineering.

Sign in to comment.

More Answers (1)

William Rose
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
where v=average velocity, ρ=density, and =Darcy friction factor, which is dimensionless. The Darcy fricton factor complicated. If the flow is laminar, then
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?
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.
  4 Comments
William Rose
William Rose on 24 Jan 2022
You can compute hn exactly: satisfies the equation
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).
In Matlab, you can solve for as follows:
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);
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
fprintf('hn=%.4f\n',hn);
hn=1.9365
The small difference in the value of found by Excel and Matlab is because they have different default tolerances for finding a solution.
___________
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.
Now that we have a way to compute , does that help you solve your problem?

Sign in to comment.

Products


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by