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

30 次查看(过去 30 天)
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
  2 个评论
Walter Roberson
Walter Roberson 2018-3-1
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
Seanice Thompson 2018-3-1
编辑:Seanice Thompson 2018-3-1
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.

请先登录,再进行评论。

回答(2 个)

Walter Roberson
Walter Roberson 2018-3-1
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

Nikolaj Maack Bielefeld
编辑:Nikolaj Maack Bielefeld 2020-3-28
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')
  2 个评论
Walter Roberson
Walter Roberson 2020-3-28
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)

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Line Plots 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by