How do I complete my code to plot the Moody Chart?

function [f] = frictionFactor(Re, ed) %Re = Reynolds Number, ed = eps/d, relative roughness

colebrook = @(f) 1/sqrt(f)+2*log10((ed/3.7)+(2.51)/(Re*sqrt(f)));

if Re > 4000 %turbulent

      f = fzero(colebrook, [0.008, 0.1]);

elseif Re < 2000 %laminar

      f = 64/Re;

else %transitional

      f = (((Re-2000)/(4000-2000))*(0.1-0.008))+0.008;

end

array = linspace(0.000001, 0.05, 21);

for i = array

     colebrook(i)
  end

end

[Walter Roberson]

It is not clear to me why you calculate f and then ignore it when you for i = array ?

I somehow suspect that the values in array are intended to represent different Re values that you want to evaluate colebrook with after figuring out what f value you want to use ?

[Seanice Thompson]

I was told to create an array of roughnesses and plug them into the colebrook equation. Then, plot the moody. I'm just not sure how to go about that. The for loop portion is what I need help with.

Answer 1

[Walter Roberson]

function [f, rough] = frictionFactor(Re, ed) %Re = Reynolds Number, ed = eps/d, relative roughness

colebrook_fed = @(f, ed) 1/sqrt(f)+2*log10((ed/3.7)+(2.51)/(Re*sqrt(f)));

if Re > 4000 %turbulent

        f = fzero(@(f) colebrook_fed(f, ed), [0.008, 0.1]);

elseif Re < 2000 %laminar

        f = 64/Re;

else %transitional

        f = (((Re-2000)/(4000-2000))*(0.1-0.008))+0.008;

end

array = linspace(0.000001, 0.05, 21);

narray = length(array);
rough = zeros(1, narray);
for K = 1 : narray
  i = array(K);
  rough(K) = colebrook_fed(f, i);
end

Answer 2

[Nikolaj Maack Bielefeld]

You could also have used the symbolic math implicit plot command:

close all; clear; clc

% symbolic math: y = f, x = Re
syms x y

% relative roughness
relrough = [0 2e-7 1e-6 5e-6 10e-6 50e-6 100e-6 200e-6 400e-6 600e-6 ...
    800e-6 1e-3 2e-3 4e-3 6e-3 8e-3 10e-3 15e-3 20e-3 30e-3 40e-3 50e-3];

% Colebrook equation
eqn = 1/sqrt(y) == -2*log10(relrough/3.7+2.51/(x*sqrt(y)));

% implicit plot with x- and y-limits
fimplicit(eqn,[2000 10e8 0.006 0.1]) 

set(gca, 'XScale', 'log') % logarithmic x-axis
set(gca, 'YScale', 'log') % logarithmic y-axis
title('Moody Chart')
xlabel('Reynolds Number')
ylabel('Friction Factor')

[Walter Roberson]

Because of the three different ranges of values, your lower bound on x for the fimplicit should be 4000. You would need to also hold on and plot the other two parts (you could plot both together if you used piecewise)

[Nikolaj Maack Bielefeld]

Yes, I agree Walter.

I just posted such a plot in this thread:

https://www.mathworks.com/matlabcentral/answers/110994-labelling-three-axes-of-moody-plot