Numerical solution of Darcy equation for headloss and velocity 1

[chrisjc]

Hi

Im having problems with a numerical solution to these equations, and i think it relates more to the process of iterative solutions rather than to the fluid theory. For some initial values a convergent solution is found, but often the Darcy frictional headloss equation returns a value greater than the initial head meaning that the next iteration results in an error since the solution to velocity in terms of head is complex.

Given: head, diameter of the pipe, length of pipe, and roughness i employ the following procedure in an attempt to derive frictional headloss and velocity:

' HEAD LOSS
Function HeadLoss(h, d, l, k) As Double
Dim hf, ff, v, vg
vg = 0.001 ' Inital guess for velocity
Do
re = Reynolds(vg, d) ' Calculate Reynolds number
ff = ff_haaland(k, d, re) ' Calculate friction factor
hf = hf_darcy(ff, l, vg, d) ' Calculate frictional head loss
v = Sqr(2 * 9.81 * (h - hf)) ' Calculate velocity *** ERROR if hf >= h ***

If Abs(v - vg) > 0.00001 Then ' If calculated v is different from vg recalculate with v
vg = v
Else
HeadLoss = hf
Exit Do
End If
Loop
End Function

I guess incrementing the velocity guess by a small amount on each iteration until the calculated velocity equals the guessed velocity but i thought the above method should solve it in less iterations. I'm not sure if i'm going about this in the right way any ideas?

Thanks
Chris

 

[BigInch]

Maybe you made a very bad guess for friction factor. I can't tell with what you posted, or rather didn't post.

And the iteration is to find f for a given velocity, f = 1/f^0.5

Given V
calculate R

Do
f
Exit Do

Once you have f, do not recalculate velocity, find the head loss.
Start with f=0.02

I attached my iterative solution for Colebrook-White ff
Open my attachment then insert a value of 0.001 in the orange cell where I have 0.02. If you start with a bad guess, you arrive at the wrong solution for f.


Download this,
On page 676 you find Colebrook-White ff

http://www.chem.mtu.edu/~fmorriso/cm310/MoodyLFpaper1944.pdf

"People will work for you with blood and sweat and tears if they work for what they believe in......" - Simon Sinek

 

[BigInch]

If you start with a good guess, you only need 3 or 4 iterations to get close enough where any further values don't contribute to the accuracy of the solution of head loss very much at all. I typically only carry out to 4 iterations.

"People will work for you with blood and sweat and tears if they work for what they believe in......" - Simon Sinek

 

[chrisjc]

Hi thanks for the link,

Sorry, for obscure post!

My objective is to calculate head-loss for a given pressure head and set of pipe dimensions so velocity, reynolds number, fiction factor and head loss are all unknown.

From Bernoulli i have velocity = sqrt(2 g (h - headloss_friction))
But headloss_friction depends on reyolds number, friction_factor, and ultimately, velocity.

The VB code starts with a low velocity estimate, then calculates reynolds number, friction factor, headloss, and finally, velocity.
If the calculated velocity is different from the estimate, then it is used for the estimate on the next iteration.

If the pipe diameter is large then the solution tends to converge quickly, however if it is small - but not too small to prevent flow under the given head - then one of the early iterations results in a head loss greater than the head driving the flow so the calculation fails due to the complex root.

I have got round this issue by just incrementing the velocity guess by a small amount on each iteration, but it's not very efficient and it'd be good to know how to get around this kind of numerical problem!

Thanks

 

[BigInch]

Pressure head does not enter into the problem. That's a boundry condition.
You are confusing pressure head with differential head!

You can get a velocity from pressure head, but not when you consider friction loss. Friction loss depends on the differential head, pressure head at inlet - pressure head at outlet. You have only assumed an inlet pressure head, therefore you are in effect trying to solve for the differential head using 0 as outlet head. As soon as you say outlet head = 0, then you have already solved your head loss problem, haven't you, so why calculate it again considering friction factor. It PROBABLY WON'T MATCH unless you have such a high flow (high velocity) that pressure at the outlet is 0.... or less if that's possible. For all flowrates, velocities that give you a pressure drop less than inlet head - 0, your answer will be wrong. For all flowrates that give a pressure drop greater than inlet head, your answer will be wrong. It will only be correct for the flowrate and velocity that gives you a differential head loss of H_in - 0



"People will work for you with blood and sweat and tears if they work for what they believe in......" - Simon Sinek

 

[IRstuff]

Are these equations so intractable that they can't be solved in Excel or Mathcad or SMath Studio?

TTFN
faq731-376

 

[BigInch]

No point in solving equations based on incorrect understanding of the problem.

"People will work for you with blood and sweat and tears if they work for what they believe in......" - Simon Sinek

 

[prex]

crisjc, the reason for the 'instability' of your calculation is that, when the initial guess for v is too far from the actual value, the solution does not converge.
You should start with a better guess: the simplest way is to start with the velocity corresponding to zero friction losses. To make also more stable the calculation you can limit the change of the calculated v[sub]new[/sub] at each iteration, e.g.
[tt]if (v-v[sub]new[/sub])/v>0.5 then v[sub]new[/sub]=0.5v
else if (v[sub]new[/sub]-v)/v>0.5 then v[sub]new[/sub]=1.5v[/tt]

prex
[URL unfurl="true"]http://www.xcalcs.com[/url] : Online engineering calculations
[URL unfurl="true"]http://www.megamag.it[/url] : Magnetic brakes and launchers for fun rides
[URL unfurl="true"]http://www.levitans.com[/url] : Air bearing pads

 

[chrisjc]

Hi BigItch

Thanks for the help, it would be good to know if ive got the theory backward and inside out!

Perhaps it is my poor use of terminology leading to confusion though, i think i meant hydrostatic head,potential head,differential head or whichever one it is!

The way i understand the problem is:
The potential energy (hydrostatic head) of elevated fluid is converted to kinetic energy by a nozzle discharging at a lower potential energy, and both the input fluid and discharged fluid are at atmospheric pressure. The rate of discharge depends on the size of the nozzle and the gauge pressure just before the fluid enters the nozzle. The complication arises when the system enters its dynamic, steady flow, state because the friction in the pipe causes a loss of pressure just before the nozzle and I need a numerical solution for velocity and frictional head loss because each depends on the other.

I think this model ignores the loss of gauge pressure just behind the nozzle caused by the conversion of potenital energy into kinetic energy of the fluid in the pipe because the pipe is large compared to the nozzle, but this bit is a little vague in my head at the the moment! Im trying to work out the optimal pipe size, rotor diameter, etc for a pelton turbine BTW.

thanks again

 

[chrisjc]

Hi Prex,

Yep absolutely, i did try this for starting values but unfortunately the head-loss is quite significant for some of the various pipe diameters that my model must consider... I am looking at small capacity rural hydroelectric system that will have relatively long penstock pipework which will be the major cost of the system, so in order to meet the budget I may have to sacrifice quite a bit of potential head.

I will try limiting the values as suggested.

Cheers

 

[chrisjc]

Hi IRstuff,


No MathCad's solve block black-artiness solves the lot in one go. I am using excel because I must solve for many different values and I couldn't work out how to parametrize the solve block in mathcad. It works fine when assigning the Find(x,y,z,...) results to an explicit matrix, but for some reason it wont let me parameterize it due to multiple datatypes (as stated in the help page for find)... any ideas how to get round this?

I'll have do some reading maybe excel has some functions that can help solve systems of equations quickly!

 

[katmar]

chrisjc, your terminology is probably what is causing you and us the most confusion.

The differential head is the difference between the inlet and discharge conditions. The head at the inlet can be made up of pressure (for example from a pump) plus the static head difference between inlet and outlet. It looks like your discharge conditions are atmospheric so that makes it easier. For example, if you have a pump delivering 20 ft of head, and the discharge is 5 ft above the inlet then your differential head is 15 ft.

Now you have to guess (calculate) a flowrate that gives a pressure loss exactly equal to the available differential. The pressure loss is made up of friction in the pipe and fittings, plus any entrance and exit losses. All of these can be included in the Darcy-Weisbach equation by determinig resistance coefficients (K values) for each of them. It looks to me as though you are confusing the issue by adding in the Bernoulli term for velocity head separately. The velocity head is the same thing as the exit loss and should all be included in your calculation for head loss. You have to keep iterating until your head loss is exactly equal to your available head.

With liquids the velocity head is generally a minor part of the overall loss, and if you are using this to calculate the velocity you are using the tail to wag the dog. Any small error will be magnified and the solution will be overly sensitive.

For your second estimate of the velocity (after your initial guess) make use of the fact that the losses are roughly proportional to the square of the flow. ie
(Q1/Q2)^2 = DP1/DP2. If your last flow estimate was Q1 and it resulted in head loss DP1, and you eventually want the head loss to be equal to the differential head available (DP2) then Q2 = Q1 x sqrt(DP2/DP1). This will be a lot more stable. Good initial guesses help too. You don't give any units but if your guess for velocity was 0.001 ft/s or 0.001 m/s then that is very slow and would mean that somewhere along the line you will have to cross over from laminar to turbulent flow. This always complicates the calculation and makes it less stable.

Katmar Software - AioFlo Pipe Hydraulics

http://katmarsoftware.com

"An undefined problem has an infinite number of solutions"

 

[chrisjc]

Many thanks Katmar,

Your post has greatly clarified the problem.

The velocity head and pressure drop terms were separated purely because that is how I dug the equations out of my textbook. Joining the terms into one energy balance equation makes a lot more sense. As it happens what I'm really getting at is what diameter pipe is required for a specific flow-rate, given the available static head, or rather, what combination of pipe and nozzle diameters give the most cost-effective design at the required flow-rate. But i am on the right track... I'll try to be a bit more sensible with the estimates and guess increments!

Cheers

 

[BigInch]

You seemed to be saying that Pressure head produces some velocity = V^2/2/g, but that is without considering nozzle outlet coefficient, for one, and the friction effect of any downstream pipe line, or other flow element, and without assuming any backpressure on the outlet side. Those will tend to reduce the velocity you had just assumed. Flow and the resulting velocity through anything, or through some pipe diameter, or hole in a tank, a pipe, or flow through a length of pipe is always controlled by the inlet head - outlet head, differential pressure. Flow through a significant length of pipe, has friction as well. Since you only had assumed the inlet head, outlet head by default, became zero. The only valid answer you could have was the flowrate & resulting velocity in the pipe that you needed to make that true; Head Loss = P_in - P_out = P_in - zero = P_in. That's only going to be true for the right combinations of P_in, pipe diameter, viscosity and length of pipe. Changing any one of those parameters will change all the others except for one, [Δ]H, which was kind of (losley) defined as = P_in -0.


"People will work for you with blood and sweat and tears if they work for what they believe in......" - Simon Sinek

 

[BigInch]

Make the spreadsheet better and add a laminar flow friction calculation.
if Reynolds number, Re < 2300, then f = 64/ Re

"People will work for you with blood and sweat and tears if they work for what they believe in......" - Simon Sinek

 

[BigInch]

gremlin attack
EVEN BETTER STILL, GRAB THIS CORRECTED VERSION.


"People will work for you with blood and sweat and tears if they work for what they believe in......" - Simon Sinek

 

[CARF]

Dear all,

Here are "all" equations, including the iterative algorithm of
Colebrook-White.

Link

A freeware fluid flow calculator based on these equations can be downloaded on bottom of the webpage.

Or a direct web application:

Link

Hope this helps to solve your problem,

CARF