To simulate Air Standard Cycle(Otto-cycle) by creating PV diagram and to calculate output thermal efficency and work done by using MATLAB

Objective- To simulate Air Standard Cycle(Otto-cycle) by creating PV diagram and to calculate output thermal efficency and work done by using MATLAB.

Theory-

Process of Otto cycle

process 1-2 isentropic compression

process 2-3 constant volume heat addition

process 3-4 isentropic expension

process 4-1 constant volume heat rejection

Ideal gas equation

PV/T = constant

Isentropic process

PV^gamma = constant

Geometric parameters

Bore of Cylinder(bore_dia) = 0.1

Stroke of Cylinder(stroke) = 0.1

Connecting rod length(c_n) = 0.2

Compression ratio(cr) = 12

Swept volume(v_s)=pie/4*bore_dia^2*stroke

Clearance volume(v_c)=v_s/(cr-1)

Crank radius(a)=stroke/2

Ratio of Connecting rod length and Crank Radius(R).

Equation for calculating volume at any crank angle

V=(1+0.5(cr-1)*[R+1-cos(theta)-((R^2-sin^2(theta))^0.5))]*v_c

Inputs

Initial pressure(p1)=101325 Pa

initial temperatur(t1)=500 K

At stage-3 temperature(t3)=2300 K

STEP1:- Create a function.

function [v]=Engine_kinematic(bore,stroke,c_n,cr,v_c,start_crank,end_crank)

a=stroke/2
R=c_n/a

theta=linspace(start_crank,end_crank,200)
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

STEP2:- Creat the code for calculation and plot

% poogram for Otto Cycle
close all
clear all
clc

%Input geometrical parameter
stroke=0.1
bore=0.1
c_n=0.2
cr=12
gamma=1.4
%Swept_volume
v_s=(pi/4)*bore^2*stroke
%Clearence_volume
v_c=v_s/(cr-1)

%Point 1
p1=101325
t1=300
v1=v_s+v_c

%point 2
v2=v_c 
p2=(p1*v1^gamma)/v2^gamma  %p1v1^gamma=p2v2^gamma
k=(p1*v1)/t1               %p1v1/t1=p2v2/t2=constant(k)
t2=(p2*v2)/k

%Compression Process
v_compression=Engine_kinematic(bore,stroke,c_n,cr,v_c,180,0)
constant=p1*v1^gamma
p_compression=constant./v_compression.^gamma

%point 3
v3=v2
t3=2300
p3=(k*t3)/v3

%Expension Process
v_expension=Engine_kinematic(bore,stroke,c_n,cr,v_c,0,180)
constant=p3*v3^gamma
p_expension=constant./v_expension.^gamma


%point 4
v4=v1
p4=(p3*v3^gamma)/v4^gamma
t4=(p4*v4)/k

Volume=[v1 v_compression v2 v3 v_expension v4 v1]
Pressure=[p1 p_compression p2 p3 p_expension p4 p1]

Work_done=trapz(Volume,Pressure)

Thermal_efficency=1-1/(cr^(gamma-1))

%Plotting
figure(1)
hold on
plot(v1,p1,'*','color','r')
plot(v2,p2,'*','color','r')
plot(v3,p3,'*','color','r')
plot(v4,p4,'*','color','r')
plot([v1 v_compression v2 v3 v_expension v4 v1],[p1 p_compression p2 p3 p_expension p4 p1])
xlabel("Volume(M^3)")
ylabel('Pressure(Pa)')
hold off
saveas(1,"Otto_cycle.png")

Result:-

WORK_DONE = 678.82 Joul

THERMAL_EFFICIENCY = 62.99%