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

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');

Answer 1

[Shadaab Siddiqie]

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.

[Nana Kwaku Okai-Mensah]

Yes it does. Thanks very much