CU Mechanical EngineeringFluid Mechanics Matlab tutorials |
Writing Matlab Functions: solving for the head loss through a pipe
In this example, we will write a function that computes the amount of pressure lost when a fluid is pushed through a pipe.
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 that 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 exercise, the name of the function will be 'pipeloss' and the first line of the function is as follows:
function [out1,out2,out3] = pipeloss(Q,L,A,Dh,e,KL,Am,nu)
Within the parenthesis are the defining parameters, which are: Q, which represents the flow rate through a system (in cubic meters per second); L, the vector of pipe lengths; A, the vector of cross-sectional areas of ducts (A(1) is area of a pipe with length L(1) and hydraulic diameter Dh(1)); Dh, the vector of pipe diameters; e, the vector of pipe roughness; KL, the vector of minor loss coefficients; Am, the vector of areas associated with minor loss coefficients (for a flow rate, Q, through minor loss element 1, the area, Am(1) gives the appropriate velocity from V = Q/Am(1). Thus, different minor loss elements will have different velocities for the same flow rate. For a sudden expansion the characteristic velocity is the upstream velocity, so Am for that element is the area of the upstream duct); and nu, the kinematic viscosity of the fluid.
It is necessary to establish the gravitational force on the system. To do this, input:
g = 9.81;
The first sequence of calculations computes the losses for any straight sections of pipes. The first step determines if the vector of pipe lengths, 'L' is empty. If this vector is nonempty, the losses are calculated using the moody function developed in the previous excersize to find the friction factor. The final term, 'hv' represents the viscous losses in the straight sections of the pipe.
if isempty(L)
hv = 0; f = [];
else
V = Q./A;
f = zeros(size(L));
for k=1:length(f)
f(k) = moody(e(k)/Dh(k),V(k)*Dh(k)/nu);
end
hv = f.*(L./Dh).*(V.^2)/(2*g);
end
Next, we compute the minor losses due to disturbances such as pipe bends, sharp corners or expansion / contraction orifices. Once again, the first test is to determine if the vector of minor losses is empty. If not, the minor losses are computed. It should be noted that the period before the operation indicates that the operation is element wise. This means that given a vector 'a', the expression 'a.^2' will square each element of the vector as opposed to performing an inner product of the vector with itself. Computing the minor losses is as follows:
if isempty(KL) hm = 0; % no minor lossess else hm = KL.*((Q./Am).^2)/(2*g); % minor lossess end
Finally, we are ready to process the information and assign values to the output parameters. The value 'nargout' is a scalar representing the number of arguments the function will return. If one argument is specified, the function will return only the total head loss. If two arguments are specified, vectors of the viscous and minor losses will be outputted. If three arguments are called for, the vectors of visous and minor losses will be returned in addition to the corresponding friction factors. Should more than three arguments be requested, the script terminates with an error. This is accomplished by the following lines of code:
if nargout==1
out1 = sum(hv) + sum(hm); % return hL = total head loss
elseif nargout==2
out1 = hv; out2 = hm; % return vectors of viscous and minor losses
elseif nargout==3
out1 = hv; out2 = hm; % return vectors of viscous and minor losses
out3 = f; % and friction factors
else
error('Only 1, 2 or 3 return arguments are allowed');
end
The m-file can now be saved as 'pipeloss.m' The following are several sample outputs from the code illustrating the ways to call the function with different types of return arguments.
>> [out1] = pipeloss(5,3,0.04,0.1,0.001, ,0.1,1e-5)
??? [out1] = pipeloss(5,3,0.04,0.1,0.001, ,0.1,1e-5)
|
Error: Expression or statement is incorrect--possibly unbalanced (, {, or [.
>> [out1] = pipeloss(5,3,0.04,0.1,0.001,1,0.1,1e-5)
out1 =
1.0342e+003
>> [out1 out2] = pipeloss(5,3,0.04,0.1,0.001,1,0.1,1e-5)
out1 =
906.7754
out2 =
127.4210
>> [out1 out2 out3] = pipeloss(5,3,0.04,0.1,0.001,1,0.1,1e-5)
out1 =
906.7754
out2 =
127.4210
out3 =
0.0380
The entire matlab script can be seen below or downloaded here (right-click to save):
function [out1,out2,out3] = pipeloss(Q,L,A,Dh,e,KL,Am,nu)
g = 9.81;
if isempty(L)
hv = 0; f = [];
else
V = Q./A;
f = zeros(size(L));
for k=1:length(f)
f(k) = moody(e(k)/Dh(k),V(k)*Dh(k)/nu);
end
hv = f.*(L./Dh).*(V.^2)/(2*g);
end
if isempty(KL)
hm = 0; % no minor lossess
else
hm = KL.*((Q./Am).^2)/(2*g); % minor lossess
end
if nargout==1
out1 = sum(hv) + sum(hm); % return hL = total head loss
elseif nargout==2
out1 = hv; out2 = hm; % return vectors of viscous and minor losses
elseif nargout==3
out1 = hv; out2 = hm; % return vectors of viscous and minor losses
out3 = f; % and friction factors
else
error('Only 1, 2 or 3 return arguments are allowed');
end