matlab code
function [Q] = RootFinder_P2(dens, visc, rough, eps)
g = 9.81;
Q = [10 1 1 1 1 2 2.2 6 6 2.3 15]/60;
D = [405 200 200 200 355 355 250 305 305 305 305]*10^-3;
L = [1250 500 400 500 400 600 1100 1100 1250 1000 1000];
Delta = 0.01;
t = 0;
while t == 0
F = flowRate(Q);
t = termCrit(F);
J = jacob(F,Q);
deltaX = (-J\F)’;
Q = Q + deltaX;
end
function [f] = flowRate(Q)
f = [Q(3)+Q(7)+Q(8)-1000/3600
-Q(6)-Q(8)+700/3600
-Q(1)+Q(5)+Q(6)+800/3600
Q(1)+Q(2)-2000/3600
-Q(4)-Q(5)-Q(7)+1300/3600
-node(Q(8),D(8),L(8))+node(Q(6),D(6),L(6))-node(Q(5),D(5),L(5))+node(Q(7),D(7),L(7))
node(Q(5),D(5),L(5))+node(Q(1),D(1),L(1))-node(Q(2),D(2),L(2))-node(Q(4),D(4),L(4))
-node(Q(7),D(7),L(7))+node(Q(4),D(4),L(4))+node(Q(3),D(3),L(3))
];
end
%jacobian
function [J] = jacob(f,Q)
%find number for matrix size
max = size(f);
J = zeros(max);
for i =1:max
for j = 1:max
Q(j) = Q(j) + Delta;
f = flowRate(Q);
J(i,j) = f(i);
Q(j) = Q(j) – Delta;
f = flowRate(Q);
J(i,j) = (J(i,j) – f(i))/Delta;
end
end
end
%test termination criterion
function [t] = termCrit(f)
total = 0;
for i = 1:size(f)
total = total + f(i).^2;
end
Fu = sqrt(total);
if Fu > eps
t = 0;
else
t = 1;
end
end
function [n] = node(Q,D,L)
fric = frictFac(Q,D);
n = fric * (8*L)/(pi^2*g*D^5)*Q*abs(Q);
end
function [f] = frictFac(Q,D)
fe = -1.8*log10(((rough/D)/3.7)^1.11 +6.9/reynolds(Q,D));
f = (1/fe)^2;
end
function [Re] = reynolds(Q,D)
Re = (4*abs(Q)*dens)/(pi*D*visc);
end
end