Alter Engine dimensions, Pr. and Temp to visualize the PV diagram and the corresponding efficiency

Alter Engine dimensions, Pressure and Temperature to visualize the PV diagram and the corresponding efficiency.

clear all
close all
clc

g = 1.4;

% Engine initializing parameters for alteration
%input variables (pressure in Nm^-2) (Temperature in K)
p1= 101325;
t1 = 500;
t3 = 2300;

% Engine Geometeric parameters for alteration
% Engine based parameters are in meters (m)
bore = 0.4;
stroke = 0.2;
con_rod = 0.15;

cr = 12.5;

%***********************COMPUTATION CODE*********************************

% Geometric variable computation

vs = (pi/4)*stroke*bore^2;
vc = vs/(cr-1);

% state variable computation

v1 = vs+vc;
v2 = vc;

p2 = p1*cr^g;

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

%dynamic compression stroke computation
c=p1*v1^g;
V_comp = Engine_kinematics_function(bore, stroke, con_rod, cr, 0,180);
p_comp = c./(V_comp.^g);




%for point 3

v3=v2;

p3 = p2*t3/t2;

%dynamic ignition stroke computation
v_igi = linspace(v2,v3,100);
p_igi = linspace(p2,p3,100);

c=p3*v3^g;
V_exp = Engine_kinematics_function(bore, stroke, con_rod, cr, 180,360);
p_exp = c./(V_exp.^g);

% for state point 4

v4=v1;
p4 = p3*(v3/v4)^g;

%dynamic exhaust stroke computation
v_exhaust = linspace(v4,v1,100);
p_exhaust = linspace(p4,p1,100);

%Thermal efficiency computation
e = 1-1/(cr)^g;



%***********************Display CODE*********************************

fprintf('The thermal efficiency of the Engine is :');
disp(e);


figure(1)
hold on;

plot (v1,p1,'.','color','r')
plot (V_comp,p_comp)
plot (v2,p2,'*','color','r')
plot(v_igi,p_igi)
plot (v3,p3,'*','color','r')
plot(V_exp,p_exp)
plot (v4,p4,'.','color','r')
plot(v_exhaust,p_exhaust)
grid on;
xlabel('Volume (V) : m^3 ')
ylabel('Pressure (P): Pa') 


%************FUNCTION CODE*************************

function [V] = Engine_kinematics_function(bore, stroke, con_rod, cr, startc, endc)

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

vs = (pi/4)*stroke*bore^2;
vc = vs/(cr-1);

theta = linspace(startc,endc,100);

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




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


end