clear; clc; close all;
P_D = 3000 + 5 * 4 * 10^6; %[Pa] Discharge pressure
P_S = 3000; %[Pa] Suction pressure
N = 9; % Number of pistons
beta = 1.45 * 10^9; % [Pa] Bulk modulus
w = 183.3; % [rad/s] Angular speed
V_p0 = 1.05 * 10^-6; % [m^3] Cylinder volume for y_p=0
A_p2 = 82.23 * 10^-6; % [m^2] Actual piston area
A_p = 83.48 * 10^-6; % [m^2] Effective piston area
R = 22.38 * 10^-3; % [m] Piston pitch radius
a = 1.45 * 10^-3; % [m] Axial distance from piston-slipper joint to yoke pivot
K_0 = 0.238 * 10^-3; % [m^2] Overlap area constant
C_d0 = 0.75; % Discharge coefficient for overlap area
n_0 = 1.489; % Overlap area constant
C_dn = 0.75; % Discharge coefficient for the relief notch
gama = 0.384; % [rad]
row = 912; % [kg/m^3] Density
K_l = 95.8 * 10^-15; % [m^3/s.Pa] Leakage constant for a piston
alfa = 0.33;
theta1 = 77 * 2 * pi / 360;
theta2 = 81 * 2 * pi / 360;
theta3 = 99 * 2 * pi / 360;
theta4 = 258 * 2 * pi / 360;
theta5 = 261 * 2 * pi / 360;
theta6 = 279 * 2 * pi / 360;
theta_span = linspace(0, pi, 1000); % Reduce points for efficiency
[th, Y] = ode45(@(th, Y) kavanagh_fcn(th, Y), theta_span, [P_D, 0, -25]); % Initialize with alfa_dot=-25
% Extract Variables
P_values = Y(:, 1);
T_values = Y(:, 2);
alfa_dot_values = Y(:, 3); % Extract stored alfa_dot values
% Plot Results
figure;
subplot(3, 1, 1);
plot(th, P_values, 'b', 'LineWidth', 1.5);
xlabel('Theta [rad]');
ylabel('Pressure [Pa]');
title('Pressure vs. Theta');
grid on;
subplot(3, 1, 2);
plot(th, T_values, 'r', 'LineWidth', 1.5);
xlabel('Theta [rad]');
ylabel('Torque');
title('Torque vs. Theta');
grid on;
subplot(3, 1, 3);
plot(th, alfa_dot_values, 'g', 'LineWidth', 1.5);
xlabel('Theta [rad]');
ylabel('Alfa Dot');
title('Alfa Dot vs. Theta');
grid on;
% ODE Function
function dYdth = kavanagh_fcn(th, Y)
P = Y(1);
T = Y(2);
alfa_dot = -25 + 50 * (th / pi); % Compute alfa_dot
% Constants
P_D = 3000 + 4 * 10^6;
P_S = 3000;
beta = 1.45 * 10^9;
w = 183.3;
V_p0 = 1.05 * 10^-6;
A_p2 = 82.23 * 10^-6;
K_0 = 0.238 * 10^-3;
C_d0 = 0.75;
n_0 = 1.489;
C_dn = 0.75;
gama = 0.384;
row = 912;
K_l = 95.8 * 10^-15;
N = 9;
A_p = 83.48 * 10^-6;
R = 22.38 * 10^-3;
a = 1.45 * 10^-3;
alfa = 0.33;
theta1 = 77 * 2 * pi / 360;
theta2 = 81 * 2 * pi / 360;
theta3 = 99 * 2 * pi / 360;
theta4 = 258 * 2 * pi / 360;
theta5 = 261 * 2 * pi / 360;
theta6 = 279 * 2 * pi / 360;
% Compute dp/dth based on theta region
dpdth = 0;
if th >= theta1 && th < theta2
dpdth = beta / (w * (V_p0 - A_p2 * ((R * sin(th) * sin(alfa) - a) / cos(alfa)))) * ...
(A_p2 * (((R * sin(th) - a * sin(alfa)) / (cos(alfa)^2)) * alfa_dot) + ...
w * R * cos(th) * tan(alfa) + ...
K_0 * C_d0 * ((theta2 - th)^n_0) * sqrt(2 / row * abs(P_D - P)) - ...
R^2 * C_dn * (tan(gama)^2) * ((th - theta1)^2) * sqrt(2 / row * abs(P - P_S)) - K_l * P);
elseif th >= theta2 && th < theta3
dpdth = beta / (w * (V_p0 - A_p2 * ((R * sin(th) * sin(alfa) - a) / cos(alfa)))) * ...
(A_p2 * (((R * sin(th) - a * sin(alfa)) / (cos(alfa)^2)) * alfa_dot) + ...
w * R * cos(th) * tan(alfa) - ...
R^2 * C_dn * (tan(gama)^2) * ((th - theta1)^2) * sqrt(2 / row * abs(P - P_S)) - K_l * P);
elseif th >= theta4 && th < theta5
dpdth = beta / (w * (V_p0 - A_p2 * ((R * sin(th) * sin(alfa) - a) / cos(alfa)))) * ...
(A_p2 * (((R * sin(th) - a * sin(alfa)) / (cos(alfa)^2)) * alfa_dot) + ...
w * R * cos(th) * tan(alfa) - ...
K_0 * C_d0 * ((theta5 - th)^n_0) * sqrt(2 / row * abs(P - P_S)) + ...
R^2 * C_dn * (tan(gama)^2) * ((th - theta4)^2) * sqrt(2 / row * abs(P_D - P)) - K_l * P);
elseif th >= theta5 && th < theta6
dpdth = beta / (w * (V_p0 - A_p2 * ((R * sin(th) * sin(alfa) - a) / cos(alfa)))) * ...
(A_p2 * (((R * sin(th) - a * sin(alfa)) / (cos(alfa)^2)) * alfa_dot) + ...
w * R * cos(th) * tan(alfa) - ...
R^2 * C_dn * (tan(gama)^2) * ((th - theta4)^2) * sqrt(2 / row * abs(P_D - P)) - K_l * P);
end
% Compute dT/dth
dTdth = N / (2 * pi) * (-P * A_p * (R * sin(th) - a * sin(alfa)) / (cos(alfa)^2));
% Store alfa_dot
dalfa_dot_dth = 50 / pi;
dYdth = [dpdth; dTdth; dalfa_dot_dth];
end
