Creating for loop for a matrix of data
Show older comments
Hi, I'm trying to create a plot that shows system head vs flowrate.
The flowrate (Q1x) goes from 0 to 0.05 in 100 intervals. But I'm struggling to use fzero to calculate the friction factors corresponding to the flowrates.
A1 = 0.0254;
rho = 1024;
D1 = 0.23;
mu = 0.0025;
e = 4.60e-05;
L1 = 77.6;
V1 = 0.5254;
g = 9.81;
h_J = 16.3347;
h1 = 31.9276;
%Colebrook-White equation
ceq = @(f,edD,Re) 1/sqrt(f)+2*log10(edD/3.7+2.51/(Re*sqrt(f)));
Q1x = linspace(0,0.05,100);
n = length(Q1x);
for i = 1:n
V1x = Q1x/A1;
Re_1x = (rho*V1x*D1)/mu;
%Friction factor:
format shortE; %Display 5 s.f.
edD_1x(1:1,1:100) = e/D1; %e/D
f_i1x(1:1,1:100) = 1e-1;
ceq_set1x = @(f) ceq(f,edD_1x,Re_1x);
options1x(1:1,1:100) = optimset('Display','iter');
options = '';
fD_1x(i) = fzero(ceq_set1x,f_i1x,options1x);
f1x = fD_1x/4;
h1_frictionx = (2*f_i1x*L1*V1^2)/(D1*g);
h_system(i) = abs(h1+h1_frictionx-h_J);
end
%Plotting system head vs flowrate
figure(1)
plot(Q1x,h_system(i),'b-')
xlable('Flowrate, Q(m^3/s)')
ylable('System head, hS (m)')
title('Graph of system head (m) vs flowrate (m^3/s)')
Everytime when I run the script, there's an error for the fD_1x(i) = fzero(ceq_set1x,f_i1x,options1x) part that says 'Second argument must be a scalar or vector of length 2'.
How can I fix it?
Answers (1)
Walter Roberson
on 9 May 2023
f_i1x(1:1,1:100) = 1e-1;
That is a vector of length 100. And it is constant, not changed anywhere in the loop, so it is not clear why it was not set before the loop if it is going to be used as a vector.
options1x(1:1,1:100) = optimset('Display','iter');
A vector of options structures is... unexpected.
fD_1x(i) = fzero(ceq_set1x,f_i1x,options1x);
You are passing in a vector of length 100 for the second parameter, and a vector of options of length 100 for the third parameter.
I get the impression that you are thinking that you can pass in a vector of problems to work on, and that fzero() will handle them all. But that is not the case: fzero() can only handle one problem per call.
1 Comment
Starting Q1x below 1e-4 can result in NaN or in fzero attempting to use complex locations.
I had to modify the h1_frictionx calculation. You were using f_1x which is independent of the results calculated by fzero, so you were getting constant output. I modified it to use f1x that was just calculated.
format shortE; %Display 5 s.f.
A1 = 0.0254;
rho = 1024;
D1 = 0.23;
mu = 0.0025;
e = 4.60e-05;
L1 = 77.6;
V1 = 0.5254;
g = 9.81;
h_J = 16.3347;
h1 = 31.9276;
%Colebrook-White equation
ceq = @(f,edD,Re) 1/sqrt(f)+2*log10(edD/3.7+2.51/(Re*sqrt(f)));
Q1x = linspace(1e-4,0.05,100);
n = length(Q1x);
edD_1x = e/D1; %e/D
f_i1x = 1e-1;
options1x = optimset('Display','final');
for i = 1:n
V1x = Q1x(i)/A1;
Re_1x = (rho*V1x*D1)/mu;
%Friction factor:
ceq_set1x = @(f) ceq(f,edD_1x,Re_1x);
ceq_set1x(f_i1x)
fD_1x(i) = fzero(ceq_set1x,f_i1x,options1x);
f1x = fD_1x(i)/4;
h1_frictionx = (2*f1x*L1*V1^2)/(D1*g);
h_system(i) = abs(h1+h1_frictionx-h_J);
end
%Plotting system head vs flowrate
figure(1)
plot(Q1x,h_system,'b-')
xlabel('Flowrate, Q(m^3/s)')
ylabel('System head, hS (m)')
title('Graph of system head (m) vs flowrate (m^3/s)')
Categories
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
on 9 May 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
