Iterative method to find flowrate in pipe in series

Hello,

I've been trying to finish this code for a few days now. It's quite frustrating since I know exactly what to do but I can't write it in Matlab. I have to find lambda and the Reynolds numbers for the 2 pipes so that I can then calculate the flowrate Q. I need to use an iterative method (using the previous values Matlab calculated to use them in the next loop) but I can't seem to make mine work. My idea is to use a sequence and matlab to loop it until the difference between the previous value of lambdas and the value calculated in the loop is smaller than the tolerance number. I would appreciate it if someone could put me in the right path on how to properly use the iterative method in this case.

The code doesn't work yet because I can't get the lambda values to be close enough yet.

Thank you

clear;
clc;
% The physical properties
rho = 1000;
mu = 1.31e-3;
g = 9.81;
% The system properties
d1 = 0.1;%input ('Diameter of pipe 1=');
l1 = 2000;%input ('length of pipe 1=');
d2 = 0.05;%input ('Diameter of pipe 2=');
l2 = 100;%input ('length of pipe 2=');
Hc=5; %headloss (45-40m)
A1= (pi*(d1/2)^2/4);
A2= (pi*(d2/2)^2/4);
tol=10e-6;%tolerance
i=1;
% Initial guess of the Reynolds Number
Re1(1)= 2.8e5;

for i=1:100 %loop to repeat 100x
    
lam1(1)= 64/Re1(1); %lambda in pipe 1
V1= sqrt(2*g*d1*Hc/(lam1(1)*l1));% Calculating the velocity in pipe 1 from the Darcy Eq
V2= V1*(A1/A2); %using continuity eq
Re2(1) = rho*(V1*(A1/A2))*d2/mu; %Re in the 2nd pipe
lam2(1)= 64/ Re2(1); %lam2 in 2nd pipe
if lam1(1)-lam2(1)>tol %i'd like to use a while loop that would stop when the 2 lambda values are close enough (tolerance)
Re1(i+1)= rho*V1*d1/mu;
Re2(i+1)= rho*V2*d2/mu;
lam1= 64/Re1(i+1);
lam2= 64/Re2(i+1);
else
    break
end
    end
end
disp ('end')
Q1 = V1*pi*d1^2/4;
Q2 = V2*pi*d2^2/4;

Answer 1

[Alan Stevens]

Why would you assume the friction factor in the second pipe is the same as that in the first?

You ony apply the head loss to the first pipe. In this case you know everything needed to calculate the volumetric flowrate without iteration:

% g*Hc = lam*(L/D)*(1/2)*V^2           (1)  Darcy  L = length D = diameter
% Re = rho*V*D/mu;                     (2)  Reynolds number
% lam = 64/Re;                         (3)  laminar friction factor
% From (2) and (3):  lam = 64*mu/(rho*V*D);  (4)
% Use (4) in (1):    g*Hc = 64*mu/(rho*V*D)*(L/D)*(1/2)*V^2;   (5)
% Rearrange (5):     V = g*Hc*rho*D^2/(32*mu*L)    (6)
% Hence volumetric fowrate:  Q = V*(pi*D^2/4);

[Michael]

I forgot to point out that both pipes have an absolute roughness of 0.01mm. I was told that we had to use the sum of darcy's equation (in both pipes) as part of our approach.

I will try your idea later and see if it works out for me

[Alan Stevens]

You were using a laminar friction factor, so pipe roughness is irrelevant. If your head loss is total, across both pipes, then you need to use both the fact that V1A1 = V2A2 and that H1 + H2 = Hc. You can then express everything in terms of, say, H1, and will need to solve for that

[Michael]

This is what I wrote. I don't think my code is right, shouldn't the lambdas come out as the same? I get the same flowrate but the velocity seems a bit unrealistic for a pipe (96m/s for pipe1). I was wondering in your code, how does matlab know what V is at the beginning when calculating lam for example (since it's an unknown) when V is actually calculated later in the loop? I always get an error message if I put V and it's an unkown

I hope i'm clear!

best

clear;
clc;
% The physical properties
rho = 1000;
mu = 1.31e-3;
g = 9.81;
D1 = 0.1;%input ('Diameter of pipe 1=');
L1 = 2000;%input ('length of pipe 1=');
D2 = 0.05;%input ('Diameter of pipe 2=');
L2 = 100;%input ('length of pipe 2=');
Hc=5; %headloss   
A1= (pi*(D1/2)^2/4);
A2= (pi*(D2/2)^2/4);
% Initial guess of the Reynolds Number
Re= 2.8e5;
for i= 1:200
%laminar friction factor
lam = 64/Re;   
%Reynolds number
V1 = sqrt(g*Hc*rho*D1^2/(32*mu*L1));
Re = rho*V1*D1/mu;
Hc = lam*(L1/D1)*(1/2)*V1^2*g;
V2 = V1*(A1/A2);
Re2= (rho*V2*D2)/mu;
lam2=64/Re2;
end
Q = V1*(pi*D1^2/4);

[Alan Stevens]

Apologies - the square root for V1 shouldn't be there!

The following code

% Pipe flow.  Two pipes in series.  Assuming
% 1. Given head loss is applied to first pipe only
% 2. The flow regime is laminar.

% g*Hc = lam*(L/D)*(1/2)*V^2           (1)  Darcy  L = length D = diameter
% Re = rho*V*D/mu;                     (2)  Reynolds number
% lam = 64/Re;                         (3)  laminar friction factor
% From (2) and (3):  lam = 64*mu/(rho*V*D);  (4)
% Use (4) in (1):    g*Hc = 64*mu/(rho*V*D)*(L/D)*(1/2)*V^2;   (5)
% Rearrange (5):     V = g*Hc*rho*D^2/(32*mu*L)    (6)
% Hence volumetric fowrate:  Q = V*(pi*D^2/4);

% The physical properties
rho = 1000;
mu = 1.31e-3;
g = 9.81;
D1 = 0.1;%input ('Diameter of pipe 1=');
L1 = 2000;%input ('length of pipe 1=');
D2 = 0.05;%input ('Diameter of pipe 2=');
L2 = 100;%input ('length of pipe 2=');
Hc=5; %headloss   
A1= (pi*(D1/2)^2/4);
A2= (pi*(D2/2)^2/4);

V1 = g*Hc*rho*D1^2/(32*mu*L1);
Re1 = rho*V1*D1/mu;
Q = V1*A1;
V2 = Q/A2;
Re2 = rho*V2*D2/mu;
disp([V1 V2])
disp([Re1 Re2])
disp(Q)

Produces the following for velocities (I'm assuming your units are all SI): V1 = 5.85 m/s V2 = 23.4 m/s

For Reynolds numbers: Re1 = 4.466*10^5 Re2 = 8.932*10^5

For volumetric flow rate Q = 0.0115 m^3/s

These Reynolds numbers are too high for laminar flow, so either (1) the head loss is meant to apply to the total across both pipes or (2) you need to assume turbulent flow using, say, the Colebrook-White friction factor correlation, or (3) both asumptions need to change if change (1) still results in Reynolds numbers too high for laminar flow.

Note that in turbulent flow, my expression for V1 does not hold.

[Michael]

Here is my updated code. It gives me the same lambdas to 0.001 difference. I just wrote Q1 Q2 to make sure i got the same ones. It seems right to me since the values seem correct V1= 0.36m/s and V2= 1.43 m/s for those pipes. I used the swamee jain eq which is an approximation of the colebrook white one.

clc;clear;
%The physical properties
rho = 1000;
mu = 1.31e-3;
g = 9.81;
D1 = 0.1;%input ('Diameter of pipe 1=');
L1 = 2000;%input ('length of pipe 1=');
D2 = 0.05;%input ('Diameter of pipe 2=');
L2 = 100;%input ('length of pipe 2=');
Hc=5; %headloss   
A1= (pi*(D1/2)^2/4);
A2= (pi*(D2/2)^2/4);
V1 = g*Hc*rho*D1^2/(32*mu*L1);
Re1 = rho*V1*D1/mu;
Q = V1*A1;
V2 = Q/A2;
Re2 = rho*V2*D2/mu;
% Re1 and Re2 > 4000 so use Swamee Jain eq
pR = 0.01; %pipe roughness in mm (epsilon/D)
f1= 0.25/(log10((pR/3.7)+5.74/Re1^0.9))^2;
f2= 0.25/(log10((pR/3.7)+5.74/Re2^0.9))^2;
%lambda = f
V1= sqrt (2*g*D1*Hc/(f1*L1));
Re1 = rho*V1*D1/mu;
V2= V1*(A1/A2);
Re2 = rho*V2*D2/mu;
Q1 = V1*A1;
Q2= V2*A2;
disp([V1 V2])
disp([Re1 Re2])
disp([Q1 Q2])
disp ([f1 f2])

[Alan Stevens]

You have used the absolute roughness (pr = 0.01mm) in the friction factor correlation, where you should have the relative roughness pr/D. Also, be careful with units; you give pr in mm but your D1 and D2 are in m I think.

However, with those corrections the two friction factors are still almost the same!