MATLAB program to solve otto cycle and plot PV diagram.

Aim: Program to solve otto cycle and plot PV diagram and finding thermal efficiency.

Assumptions: There are four process in otto cycle:

1) Isentropic Compression

2) Isochoric heat addition

3) Isentropic expansion

4) Isochoric heat rejection.

There are two constant volume process so we can say:

v1=v4

v2=v3

input parameters are p1 and t1 and output parameter is t3, by above 3 parameters we can calculate v1,v2,p2,p3,p4,t2,t4.

Output: PV diagram

https://drive.google.com/file/d/1i_aVC-av1pGmqdFkha-AaStYZKAJFWoc/view?usp=drivesdk

Thermal_Efficiency of a Engine is = 39%

Program:

close all
clear all
clc


% otto_cycle_pv

% inputs

p1 = 101325
t1 = 560
t3 = 2400
gamma = 1.2

% engine_parameters

bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
CR = 12;

v_swept = (pi/4)*bore^2*stroke;
v_clearence = v_swept/(CR-1);

v1 = v_swept+v_clearence;
v2 = v_clearence;


% p1v1^gamma = p2v2^gamma | p2 = p1*(v1/v2)^gamma | CR = v1/v2

p2 = p1*(CR)^gamma

%p1v1/t1 = p2v2/t2

t2 = (p2*v2*t1)/(p1*v1)

constant_c = p1*v1^gamma
V_compression = otto_cycle_pv_kinematics(0.1,0.1,0.15,12,180,0);
P_compression = constant_c./V_compression.^gamma;

% p2v2/t2 = p3v3/t3 | p3 = p2v2t3/v3t2
 
 p3 = p2*(t3/t2)
 
 % p3v3^gamma = p4v4^gamma | v2=v3 | v4=v1
 
 v3 = v2;
 v4 = v1;
 p4 = p3*(v2/v1)^gamma
 
 constant_c = p3*v3^gamma
V_expansion = otto_cycle_pv_kinematics(0.1,0.1,0.15,12,180,0);
P_expansion = constant_c./V_expansion.^gamma;

 % p3*v3/t3 = p4*v4/t4
 
 t4 = (p4*v4*t3)/(p3*v3);
 
 % therma_efficiency
 
 T_E = (1-(t4-t1)/(t3-t2))*100;
 
%plotting

figure(1)
hold on
plot(v1,p1,'*','color','r');
plot(V_compression,P_compression)
plot(v2,p2,'*','color','r');
plot([v2 v3],[p2 p3],'linewidth',1)
plot(v3,p3,'*','color','r');
plot(V_expansion,P_expansion)
plot(v4,p4,'*','color','r');
plot([v4 v1],[p4 p1],'linewidth',1)


xlabel('volume');
ylabel('pressure');





function [V] = otto_cycle_pv_kinematics(bore,stroke,con_rod,CR,start_crank,end_crank)

a = stroke/2;
R = con_rod/a;

v_s = (pi/4)*bore^2*stroke;
v_c = (v_s)/(CR-1);

theta = linspace(start_crank,end_crank,100);

term1 = 0.5*(CR-1);
term2 = R+1-cosd(theta);
term3 = (R^2-sind(theta).^2).^0.5;

V = (1+term1*(term2-term3)).*v_c

end