How to solve a Bernoulli Equation

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?

[Torsten]

No explicit time derivatives appear in the equation. So you can use "fsolve" (or a calculation with pencil and paper because only quadratic and linear terms are present).

[KarolN]

Thanks for answering. Do you have any hints how to set the initial value problem when using fsolve/vpasolve?

[Torsten]

What initial value problem ?

You have a known state (h0,v0). You can calculate the left-hand side of the Bernoulli equation. Then either height h0 or velocity v0 change. If h0 changes to h1, v0 changes to v1 according to Bernoulli equation. If v0 changes to v1, then h0 changes to h1 according to Bernoulli equation.

Maybe you could explain in more detail the application where you try to use Bernoulli.

[William Rose]

@KarolN,

@Torsten is correct, as usual. The equation is not a differential equaiton and it is not an initial value problem. Write the differential equation. Also, make it clear which quantitiess in your equation are constants and which quantities are functions of time, or position, or both.

The units of all terms on both sides appear to be length - assuming and and α are dimensionless.

Your Matlab equation does not match your written equation. You need to add () around the denominator terms:

heightA+Id*DeltaL+turbulence*v^2/(2*gravityG) == heightA+turbulence*v^2/(2*gravityG)+Ifsr*DeltaL

[KarolN]

@Torsten I begin to get it. The application is for unpressurized canal in an earth dam (my semestral project).

@William Rose According to the instruction the left side of the equation is constant.

Id*DeltaL symbolizes lowering of the section

Ifsr*DeltaL symbolizes energy loss on entire length of the section

DeltaL = DeltaE/Id-Ifsr; and DeltaE = E1-E2 is difference of energy on the beginning and end of the section DeltaL

heightA - subsequent water flow into the section

turbulence and gravity are constant values taken from table

All parameters had been calculated already:

All I need to know is how to set the equation to receive output like this:

[KarolN]

Should I run this in a for loop?

[William Rose]

@KarolN, I am travelling but I will look at your modified question later today.

A for loop is a good idea.

Your original equation has h1 and h2. But the table has only "h". What happened to h1 and h2?

[KarolN]

@William Rose Thank you very much for interest in my problem. The 'h' is the height of water flowing through the pipe, and as we go from section to section separated by segments [edited], h1 is initial height and h2 is final height (so h2>h1). Similarly as v is speed, v1 and v2 indicate which section or segment.

The process looks likesimilar to this:

I wrote a rough code that calculates the equation, the problem is, I have no idea how to make it to generate that table professor wants. Obviously you have to code every value separately. And so far only values afflicted by 'h' are displayed.

[William Rose]

@KarolN,

Your original equation has v1 and v2, but the table has only a single v.

Your original equation has h1 and h2, but the table has only a single h.

You said in your comment that , but the Table you provided has different values for on every row

You said in your comment that , but the table has different values for on every row.

You said "Ifsr*DeltaL symbolizes energy loss on entire length of the section". This cannot be correct, because the units do not work: is dimensionless.

You said "DeltaL = DeltaE/Id-Ifsr". Do you mean DeltaL = DeltaE/(Id-Ifsr)? Note that I added parentheses. In any case, this equation cannot be correct, becuase the units do not match, since and are dimensionless.

Please explain these discrepancies.

[KarolN]

@William Rose Right away. I made some mistakes in initial translation of that equations into matlab. DeltaL I got wrong because I ignored the index. Same goes for Ifsr. h1,h2, if1,if2, v1,v2 or E1,E2 are not separate values. They are the same value but with a different index. I made you an infographics with meaning of each formula, due to my errors and disambiguation of what we speak about.

[William Rose]

@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]

@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.

Answer 1

[William Rose]

@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.

[KarolN]

@William Rose PERFECT! Thank you very much for your devotion to this difficult issue! It was a fascinating discussion!

[William Rose]

@KarolN, you're welcome. Good luck with your civil engineering.

Answer 2

[William Rose]

@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.

[KarolN]

@William Rose As fas as I know we do not deal with pressure here, and that is why the canal is assumed to be unpressurized. I guess they simplified the task for us, so we could complete it at all. Fluid mechanics isn't even in civil engineering curriculum.

To answer your question about units: Professor said that: they are standard SI units except 'n' for concrete roughness which is in s/m^1/3.

As for calculating the normal depth 'hn' he gave me a helpful chart in excel:

As you can see, you do it by trial/error substituting values for 'hn' as long as you catch the right value.

The problem with my query is, these are only part of an entire chain of calculations I had made since october, starting with geotechnical and geological data of the region etc. I have amassed so far about 30 pages of them. Therefore it is difficult to explain all that in a single post. and if I sound confusing I apologize for that.

I am interested in mechanism, how to implement it in matlab, so it works. It doesn't matter whether I will get units or some other details wrong, as that will be worked out later during classes. But I need some helping hand with coding all that in some wholesome form: should I use for or while loops, should I run it through a vector or a matrix, vsolve or vpasolve at all etc.

I use matlab for about three months now, and I simply do not have know-how yet for more complicated issues as this one.

[William Rose]

@KarolN,

It is clear that I am not going to understand this problem. It is too complicated for me. Since this is flow in a canal, the equations I gave, which were for flow in a pipe are not applicable.

I would like to help you with the Matlab part, but I still do not know what problem you are trying to solve. I am sorry that I do not have more time to devot to this problem. Good luck with your work.

Sincererly, W. Rose

[KarolN]

@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]

@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?