How to combine the functions to get this graph - down below

조회 수: 1 (최근 30일)
The graph at the top is what I'm trying to achieve using the above equations. So I wrote a code using the using the equations and the necessary constants from the exact paper. However, my graph seems to be decreasing as compared to the paper's increasing graph.
Please, help me out.
%% Knowns
Xo = 0.0254; % (m) the clearance distance (the compression possible during the compression stroke)
R = 0.305; % (m) the radius of the flywheel
A = 0.001883; % (m^2) the piston area.
C = 0.0113; % (kgm^2) is load constant.
I = 3.171; % (kgm^2) the moment of inertia.
n = 1.3; % the polytropic index
Change_t = 0.001; % (Change in time/s) The step for the Euler approximation
% IVP
P1 = 0.1e6; % (MPa) The pressure at compression
P3 = 10.3e6; % (MPa) The pressure expansion
w = 50; % (rad/s) initial that starts the Otto cycle
theta1 = 0; % Crank angle that starts the Otto cycle
theta3 = pi;
%% Finding the angular velocity (w) and crank angle (theta)
wnew = zeros(1,1000); % Store all angular velocity values
theta_vals = zeros(1,1000); % Store all crank angle values (theta)
wnew(1) = w;
theta = 0;
theta_vals(1) = theta;
i = 1;
iter = 1; % number of of iterations
while iter <= 1000
if i > 0 && i <= pi
diff_w = P3*((R-R*cos(theta1) + Xo)/(R-R*cos(theta)+Xo))^n * ((A*R*sin(theta))/I) - ((C*w^2)/I); % First equation Ɵ'' function in paper
% Applying Euler's method of approximation
w = w + Change_t*diff_w;
wnew(iter+1) = w;
theta = theta + Change_t*w;
theta_vals(iter+1) = theta;
else (i > pi && i <= 2*pi)
diff_w = P1*((R-R*cos(theta3) + Xo)/(R-R*cos(theta)+Xo))^n * ((A*R*sin(theta))/I) - ((C*w^2)/I); % Second equation from Ɵ'' function in paper
w = w + Change_t*diff_w;
wnew(iter+1) = w;
theta = theta + Change_t*w;
theta_vals(iter+1) = theta;
end
i = i + 0.001;
iter = iter + 1;
end
% Ploting the result
wnew(wnew==0) = [];
time = ones(1,length(wnew)-1);
time = [0 time];
time = cumsum(time);
time = 0.001 * time; % Step for time
plot(time,wnew,'-b'); grid on;
xlabel('Time (seconds)'); ylabel('Angular Velocity (rad/sec)');
legend('Angular velocity over time','Location','northeast');

채택된 답변

Shadaab Siddiqie
Shadaab Siddiqie 2021년 4월 27일
From my understanding you are getting wrong result from your experiment. This might be because C value in the paper might be negative.
Thus changing lines
diff_w = P3*((R-R*cos(theta1) + Xo)/(R-R*cos(theta)+Xo))^n * ((A*R*sin(theta))/I) - ((C*w^2)/I); % First equation Ɵ'' function in paper
diff_w = P1*((R-R*cos(theta3) + Xo)/(R-R*cos(theta)+Xo))^n * ((A*R*sin(theta))/I) - ((C*w^2)/I); % Second equation from Ɵ'' function in paper
to
diff_w = P3*((R-R*cos(theta1) + Xo)/(R-R*cos(theta)+Xo))^n * ((A*R*sin(theta))/I) + ((C*w^2)/I); % First equation Ɵ'' function in paper
diff_w = P1*((R-R*cos(theta3) + Xo)/(R-R*cos(theta)+Xo))^n * ((A*R*sin(theta))/I) + ((C*w^2)/I); % Second equation from Ɵ'' function in paper
Would result in
Hope this helps you to figure out the issue.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Directed Graphs에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by