Labelling three axes of moody plot

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.

[Cesur193]

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

Answer 1

[Image Analyst]

Did you try calling zlabel()?

[David]

I did. Adding zlabel('...') under xlabel() and ylabel() doesn't seem to work.

[Image Analyst]

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');

[David]

Still nothing, as if it were never called. The z-axis will also need to be scaled differently to the y-axis. Is there a way that I could possibly restructure my code to accommodate the new axis?

Answer 2

[Nikolaj Maack Bielefeld]

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