Labelling three axes of moody plot
12 views (last 30 days)
Show older comments
I've written the following code, which use another function 'Moody', to plot the moody plot. I've plotted it correctly, I think. With log scales on both the x and y axis. The x and y axes are labelled and scaled. Is it possible for me to label the third axis? This axis denotes the relative roughness and is the axis parallel to the y-axis. Also, is there a way that I could label more of the divisions along the y axis? As it stands there are only two values but I would like to have the values between these labelled also. The code is as follows: Moody:
function f = Moody(relr,Re)
% 'moody' finds friction factor by solving the Colebrook equation (Moody Chart)
% Inputs: e/D, Re.
% Output: f.
% Note: Laminar and turbulent flow are correctly accounted for
if (Re < 0)
error (sprintf('Reynolds number = %f cannot be negative',Re));
elseif (Re < 2000)
f = 64/Re; return %Laminar flow
end
if (relr > 0.05)
warning (sprintf('epsilon/diameter ratio = %f is not on Moody chart',relr));
end
if Re<4000, warning('Re = %f in transition range',Re); end
% --- Use fzero to find f from Colebrook equation.
% coleFun is an anonymous function object to evaluate F(f,e/d,Re)
% fzero returns the value of f such that F(f,e/d/Re) = 0 (approximately)
% fi = initial guess from Haaland equation, see White, equation 6.64a
% Iterations of fzero are terminated when f is known to whithin +/- dfTol
coleFun = @(f,ed,Re) 1.0/sqrt(f) + 2.0*log10( ed/3.7 + 2.51/( Re*sqrt(f)));
fi = 1/(1.8*log10(6.9/Re + (relr/3.7)^1.11))^2; % initial guess at f
dfTol = 5e-6;
f = fzero(coleFun,fi,optimset('TolX',dfTol,'Display','off'),relr,Re);
% --- sanity check:
%if f<0
%error(sprintf('Friction factor = %f, but cannot be negative',f));
%end
Moody Plot:
% myMoody makes a simple Moody chart
% Generates a log-spaced vector of Re values in the range 2500 <= Re < 10^8
Re = logspace(log10(2500),8,50);
relr = [0 0.00001 0.00005 0.0001 0.0002 0.0005 0.001 0.002 0.005 0.01 0.02 0.03 0.04 0.05];
f = zeros(size(Re));
% Plot f(Re) curves for one value of epsilon/D at a time
% Temporarily turn warnings off to avoid lots of messges when 2000 < Re < 4000.
warning('off')
figure('units','normalized','outerposition',[0 0 1 1]);
title('Moody Plot')
hold('on');
% Nested loop to find f values.
for i=1:length(relr)
for j=1:length(Re)
f(j) = Moody(relr(i),Re(j));
end
loglog(Re,f,'k-', 'LineWidth', 1.5)
end
ReLam = [100 2000];
fLam = 64./ReLam;
plot(ReLam,fLam,'r-');
xlabel('Re');
ylabel('f','Rotation',0)
axis([0 1e8 0 2e-1]);
hold('off');
grid on
warning('on');
% MATLAB resets loglog scale? Set it back.
set(gca,'Xscale','log','Yscale','log'); %Formats/Scales axes using handles.
1 Comment
Cesur193
on 7 Jun 2023
Hello. Could you explain me how to add multiple points to your Moody's chart that you created? I am unfortunately a beginner in Matlab programming, it would be great if someone from the community could help me.
thank you very much
Answers (2)
Image Analyst
on 30 Dec 2013
Did you try calling zlabel()?
3 Comments
Image Analyst
on 30 Dec 2013
So there's absolutely no label whatsoever near the z axis? Just like it was never even called? I think sometimes labels don't show up depending on the angular orientation of the coordinate system. Try a different orientation. Or increase the font size:
zlabel('This is the Z axis', 'FontSize', 30, 'FontWeight', 'bold');
Nikolaj Maack Bielefeld
on 29 Mar 2020
Edited: Nikolaj Maack Bielefeld
on 29 Mar 2020
You can do something like that with this script that I created, maybe you can adjust your own with it (plot is below).
I utilize a function made by French mathematician Didier Clamond that computes the friction factor to machine precision with only two iterations. The method is Quartic Iterations. Read the paper here (includes code): https://arxiv.org/pdf/0810.5564.pdf
I also attached the .m file colebrook.m (from the paper) which is neccessary in order to run my script.
"Elapsed time is 0.355480 seconds."
I'm guessing plotting and looping uses most of that time. Quick plotting is actually quite an advantage if you need to plot multiple times to adjust the plot.
You might need to adjust the position of the 'Relative Roughness' and the other texts if you alter the code or stretch the plot since I couldn't figure out how to do it generically. Maybe someone has an idea for that.
close all; clear; clc
tic
% Relative roughness vector
K = [0 1e-6 5e-6 10e-6 50e-6 100e-6 200e-6 400e-6 600e-6 ...
1e-3 2e-3 4e-3 6e-3 8e-3 10e-3 15e-3 20e-3 30e-3 40e-3 50e-3];
% Reynolds numbers vector
R1 = [4e3:1e3:1e4];
R2 = logspace(4,8);
R = [R1 R2] % The Reynolds numbers used
L_K = length(K); % length of the vector
L_R = length(R); % length of the vector
f = zeros(L_K,L_R); % Preallocate friction factor matrix
for ii = 1:length(K) % Loop through all the relative roughnesses
hold on % Plot on top of each plot
% Friction factor matrix from **colebrook.m** see attachment
f(ii,:) = colebrook(R,K(ii));
plot(R,f(ii,:)) % Plot friction factor as function of Re
% Plot loglog
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')
% Labels
xlabel('Reynolds Number')
ylabel('Friction Factor')
% Limits
ylim([0.0055 0.1])
xlim([600 1e8])
% Places the relative roughness text where the line ends
xlim2 = get(gca,'xlim');xmax = xlim2(2);
txt = num2str(K(ii));
text(xmax,f(ii,end),txt)
end
% Plots the laminar region
Re_laminar = 600:100:2300;
f_laminar = 64./Re_laminar;
plot(Re_laminar,f_laminar)
% Places a red box in the transition region
h=fill([Re_laminar(end),4e3,4e3,Re_laminar(end)], ...
[f_laminar(end),f_laminar(end),f(end,1),f(end,1)],'red');
h.FaceAlpha=0.3;
grid on
% Places the relative roughness axis title attempt
xlim3=get(gca,'XLim');
ylim3=get(gca,'YLim');
ht = text(xlim3(1)+0.08*xlim3(2),ylim3(1)+0.8*ylim3(2),'Relative Roughness');
ht2 = text(xlim3(1)+0.008*xlim3(2),ylim3(1)+0.2*ylim3(2),'Turbulent Region');
% Places the text for the laminar region
ht3 = text(0.5*xlim3(1)+0.000008*xlim3(2),0.6*ylim3(1)+0.7*ylim3(2),'Laminar');
set(ht3,'Rotation',-72.5)
% Places the text for transition
txt = 'transition';
text(1300,f_laminar(end)-0.003,txt)
title('MOODY CHART: Friction Factor(Relative Roughness,Reynolds Number)')
hold off
toc
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/280362/image.png)
0 Comments
See Also
Categories
Find more on Line Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!