10/15/24 6:25 PM /Users/eleonoradisi.../newtry_lab1_AS.
m 1 of 3
clc; clear; close all;
%% Defining all the known variables of the system, using SI
v_fs = 90E-3; % Forward stroke velocity (VFS)
[m/s]
p_m = 20E6; % Main line pressure (pm) [Pa]
p_t = 0.22E6; % Tank pressure (pT) [Pa]
Lp = 28; % Distance pump-actuator (pipe
length, Lp) [m]
Dp = 15E-3; % Pipe diameter (Dp) [m]
R = 32E-3; % Bend radius (R) [m]
ep = 0.0035E-3; % Pipe roughness (ε) [m]
b = 0.15; % Actuatorbs arm (b) [m]
theta_max = deg2rad(28); % Maximum angular excursion
(θmax) [rad]
theta_nom = deg2rad(22); % Nominal angular excursion
(θnom) [rad]
a_Ma = 300; % Momentum equation coefficient
(aMa) [N*m]
b_Ma = 68000; % Momentum equation slope (bMa)
[N*m/rad]
g = 9.81; % Gravity acceleration [m/s^2]
K_dcv = 3.9; % Constant of direction control
valve
K_nrv = 5.7; % Constant of non-return valve
K_br = 0.18; % Constant of bend radius
K_ex = 1; % Constant of exit
gamma = 1.3; % Polytropic tranformation
constant
M = 50; % Equivalent mass of the
piston-flap assembly [kg]
rho = 900; % Fluid density [kg/m^3]
nu = 1.6E-5; % Kinematic viscosity [m^2/s]
E = 2.71828; % Euler's number
%% TASK 1 : Define the actuator area
Ma_max = a_Ma + b_Ma*(theta_max); % Maximum momentum resulting from the arm's
rotation.
F_max = Ma_max/b; % Maximum load
% Then, setting the force balance equation = 0 (static condition) and using
% the main line pressure :
A_a = F_max/p_m;
disp(['The actuator area is: ', num2str(A_a)])
%% TASK 2 : Define the nominal working pressure
% Educational assumption : neglect friction forces
% The nominal momentum is computed using the nominal angle
Ma_nom = a_Ma + b_Ma*(theta_nom); % Nominal momentum resulting from the arm's
rotation
F_nom = Ma_nom/b; % Nominal load
% As the actuator is required to move the nominal load at constant velocity, the
static
% force equilibrium equation can be applied as well.
p_a = F_nom/A_a;
10/15/24 6:25 PM /Users/eleonoradisi.../newtry_lab1_AS.m 2 of 3
disp(['The nominal working pressure is: ', num2str(p_a)])
%% TASK 3 : Define the value of Darcy coefficient, solving the Colebrook formula with
an iterative method
Q = A_a*v_fs; % volumetric flow rate [m^3/s]
A_p = pi*(Dp/2)^2; % Cross-sectional area of the pipes [m^2]
v_p = Q/A_p; % Velocity of the fluid through the pipes [m/s]
Re_d = (v_p*Dp)/nu; % Reynold's number across the pipe
f(1) = 0.01; % Estimated initial guess
tol = exp(-6); % Convergence tolerance
max_iter = 10; % Maximum number of iterations
% Define the function and its derivative (Colebrook equation)
F = @(f) 1/sqrt(f) + 2*log10(ep/(3.7*Dp) + 2.51/(Re_d*sqrt(f)));
F_prime = @(f) -0.5/f^(3/2) + 2*(1/log(10))*((ep/Dp)/(3.7) + 2.51/(Re_d*sqrt(f))) *
(-2.51/(2*Re_d*sqrt(f)^3));
% Iterative method for finding Darcy friction factor
for i = 1:max_iter
% Update friction factor using the Newton-Raphson method
f(i+1) = f(i) - F(f(i)) / F_prime(f(i));
% Check for convergence
if abs(f(i+1) - f(i)) < tol
fprintf('Converged: The Darcy friction factor is %.5f\n', f(i+1));
break; % Exit the loop if converged
end
% If we reach the max number of iterations without convergence
if i == max_iter
fprintf('Did not converge in %d iterations. Last f = %.5f\n', max_iter, f
(i+1));
break;
end
end
%% TASK 4 : Estimate the friction head losses in the pipe
% The friction head losses coefficent has to be computed first.
h_f = f(end) * (Lp / (2 * g)) * ((16 * Q^2) / ((pi^2) * (Dp^5)));
h_f_check = (f(end) * Lp * (v_p^2)) / (2 * Dp * g);
% Now, the friction pressure head losses (deltap_pipe) can be determined as follow :
deltap_pipe = rho * g * h_f;
fprintf("Friction head losses in the pipe : %.2e Pa\n", deltap_pipe);
%% TASK 5: Estimate the concentrated head losses in the direction control valve
deltap_DCV = K_dcv*rho*(v_p^2/2);
fprintf("The concentrated head loss in the Direction Control Valve : %.2d Pa",
deltap_DCV);
%% TASK 6: Estimate the concentrated head losses in all non-return valves
deltap_NRV_one = K_nrv*rho*(v_p^2)/2; % In one non-return valve
deltap_NRV = 2*K_nrv*rho*(v_p^2)/2; % In all non-return valves
fprintf("Concentrated head loss in all non-return valves: %.2d Pa", deltap_NRV);
%% TASK 7 : Estimate the concentrated head losses in all bend radii
10/15/24 6:25 PM /Users/eleonoradisi.../newtry_lab1_AS.m 3 of 3
deltap_BR_one = K_br*rho*(v_p^2)/2;
deltap_BR = 2*K_br*rho*(v_p^2)/2;
fprintf("Concentrated head loss in all bend radii: %.2d Pa\n", deltap_BR);
%% TASK 8 : Estimate the concentrated head losses in
deltap_EX = K_ex*rho*(v_p^2)/2;
fprintf("Concentrated head loss at the exit: %.2d Pa\n", deltap_EX);
%% TASK 9 : Define the concentrate head loss in the flow control valve
deltap_FCV = p_m-p_a-deltap_BR-deltap_pipe-deltap_DCV-deltap_NRV-deltap_EX;
fprintf("Concentrated head loss in the flow control valve: %.2d Pa\n", deltap_FCV);
%% TASK 10 : Define the rquired area ratio of the flow control valve
eqn1 = @(A_t) ((2.3 * rho * 0.5 * (Q^2)) / ((1 + (A_t / A_p)^1.8) * (A_t^2))) -
deltap_FCV;
A_t0 = A_p;
options = optimoptions('fsolve', 'Display', 'off'); % Using optimoptions for newer
MATLAB versions
A_t = fsolve(eqn1, A_t0, options);
AR = A_t / A_p;
fprintf("Task 10: Required aspect ratio of the flow control valve: %.2f\n", AR);
