CU Mechanical EngineeringFluid Mechanics Matlab tutorials |
Writing Matlab Functions: solving for the friction coefficient
Th ability to write custom functions within Matlab lends an extrodinary amount of power to the user. In this example, we will write a function that will compute the friction factor as a function of the Reynolds number and the relative roughness: epsilon / diameter.
To begin, all functions within matlab must begin by declaring that m-file is a function. This is accomplished by having the first line of the m-file having the word 'function' followed by the name of the variable which will be returned. This is consequently set equal to the name of the name of the file with the quantities that are passed to the function in parenthesis. For the purpose of this excersize, the name of the function will be 'moody' and the first line of the function is as follows:
function f = moody(ed,Re)
Within the parenthese is 'ed' which is the relative roughness and 'Re' is the Reynolds number.
The next step is to use a series of conditional tests to determine if the Reynolds number is in an appropriate range for the application of the moody chart. The first conditional statement tests to see if the Reynolds number is positive and if not displays a warning. The 'error' command will cause the script to quit without finishing if the condition proves to be true.
if Re<0
error(sprintf('Reynolds number = %f cannot be negative',Re));
Immediately following this line is another conditional statement to determine if the flow is laminar. If so, the corresponding friction factor us returned.
elseif Re<2300 f = 64/Re; return % laminar flow end
The 'return' command causes the script to terminate and since the value of 'f' is assigned, this is the value returned by the function. The conditional statments are terminated by the use use of the 'emd' command.
The function only executes beyond this point of the Reynolds number is greater than 2300. The next conditional statement determines if the relative roughness is within the range of the moody chart.
if ed>0.05
warning(sprintf('epsilon/diameter ratio = %f is not on Moody chart',ed));
end
After testing for low Reynolds numbers and high relative roughness, it is time for the meat of the function. This begins by defining an inline function with 'f', 'ed' and 'Re' being the parameters.
findf = inline('1.0/sqrt(f) + 2.0*log10( ed/3.7 + 2.51/( Re*sqrt(f)) )','f','ed','Re');
Since the process involves finding the zero of a nonseparable equation, the process is iterative and so in the next line we start with an initial guess.
fi = 1/(1.8*log10(6.9/Re + (ed/3.7)^1.11))^2; % initial guess at f
The last step is to define the tolerance needed to converge on a solution and then the 'fzero' command is used to find the zero of the inline function defined above.
dfTol = 5e-6;
f = fzero(findf,fi,optimset('TolX',dfTol,'Display','off'),ed,Re);
The function can now be saved as 'moody.m' and can be executed from the Matlab command prompt. Several sample inputs and corresponding outputs follow:
>> moody(0.001,5000)
ans =
0.0385
>> moody(0.1,5000)
Warning: epsilon/diameter ratio = 0.100000 is not on Moody chart
> In moody at 19
ans =
0.1049
The entire matlab script can be seen below or downloaded here (right-click to save):
function f = moody(ed,Re,verbose)
if Re<0
error(sprintf('Reynolds number = %f cannot be negative',Re));
elseif Re<2000
f = 64/Re; return % laminar flow
end
if ed>0.05
warning(sprintf('epsilon/diameter ratio = %f is not on Moody chart',ed));
end
findf = inline('1.0/sqrt(f) + 2.0*log10( ed/3.7 + 2.51/( Re*sqrt(f)) )','f','ed','Re');
fi = 1/(1.8*log10(6.9/Re + (ed/3.7)^1.11))^2; % initial guess at f
dfTol = 5e-6;
f = fzero(findf,fi,optimset('TolX',dfTol,'Display','off'),ed,Re);