How to solve a Bernoulli Equation

31 次查看（过去 30 天）
KarolN
KarolN 2022-1-21
评论： William Rose 2022-1-25
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 个评论
显示 10更早的评论
William Rose
William Rose 2022-1-23
@KarolN,
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 ?
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 ?
You say energy , but the units on the right hand side and length, not energy.
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.
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 ?
Do you have an equation for pressur drop per unit legth of pipe? If you do, what is that equation?
KarolN
KarolN 2022-1-23
编辑：KarolN 2022-1-23
@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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

William Rose
William Rose 2022-1-24
Ran in:
@KarolN,
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 个评论
KarolN
KarolN 2022-1-25
@William Rose PERFECT! Thank you very much for your devotion to this difficult issue! It was a fascinating discussion!
William Rose
William Rose 2022-1-25
@KarolN, you're welcome. Good luck with your civil engineering.

请先登录，再进行评论。

更多回答（1 个）

William Rose
William Rose 2022-1-23
@KarolN,
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 个评论
显示 2更早的评论
KarolN
KarolN 2022-1-24
@William Rose nevertheless I wholeheartedly thank you for your interest.
What I try to solve in this thread: I need some coding example, how to make a loop with this equation running through a matrix h =[] which contains subsequent water levels, and then displays output for all elements of the equation ex. Ifsr, deltaL etc.
William Rose
William Rose 2022-1-24
Ran in:
@KarolN,
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?

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Structural Mechanics 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by