Clear Filters
Clear Filters

IC engine plot graph

4 views (last 30 days)
raziq halim
raziq halim on 11 Dec 2020
Answered: Walter Roberson on 11 Dec 2020
can someone help with this code, I tried to change the variable N speed from 2000:2000:12000 and it seems that my graphs were still constant with 2000 rpm?
% Engine parameters
r = 10.5; % compression ratio
b = 0.077; % bore [m]
s = 0.0858; % stroke [m]
cr = 0.1365; % connecting rod length [m]
ct = s/2; % crank throw [m]
Vd = 0.0016/4;%0.25*pi*s*b^2; % displacement volume [m3]
V2 = Vd./(r-1); % TDC volume [m3]
V1 = r.*V2; % BDC volume [m3]
N = [2000:2000:12000]*(1/60); % engine speed [rev/s]
% Chemical parameters
LHV = 42000000; % lower heating value of gasoline [J/kg]
FA = 14.7; % stoichiometric air/fuel ratio
g = 1.4; % ratio of specific heats - unburned
R = 287; % particular gas constant [J/kgK]
cv = R/(g-1); % constant-volume specific heat at high temps [J/kg-K]
% Cycle parameters
CAD = 0:1:720;
CADrad = CAD*pi/180;
% Displacement and volume
y = cr+ct-(sqrt(cr.^2 - (ct.^2).*sin(CADrad).^2)+ct.*cos(CADrad));
V = V2 + 0.25.*pi.*(b.^2).*y;
% Supercharger calculation
etac = 1; % supercharger component efficiency
cps = 1001; % cp of supercharger [J/kg-K] (at a low temp)
PR = 1; % supercharger pressure ratio
Patm = 101325; % atmospheric pressure [Pa]
Tatm = 300; % atmospheric temperature [K]
% Cycle calculation
CADc = 180:1:540; % crank angles of compression and expansion strokes
Vc = V(179:539); % volumes during compression and expansion strokes [m3]
M = zeros(size(N)); % charge mass [kg]
Mf = zeros(size(N)); % mass of fuel [kg]
Qin = zeros(size(N)); % heat addition [kg]
P = zeros(length(N),length(Vc)); % pressures [Pa]
T = zeros(length(N),length(Vc)); % temperatures [K]
Wcycle = zeros(1,length(N)); % cycle power for 4 cylinders [W]
mdot = zeros(1,length(N)); % mass flow into engine for 4 cylinders [kg]
Ws = zeros(size(N)); % power of the supercharger [W]
Wnet = zeros(1,length(N)); % net power including supercharger for 4 cylinders [W]
eta = zeros(size(N)); % thermal efficiency
for count = 1:length(N)
% supercharger
P1 = PR*Patm; % intake runner pressure [Pa]
T1s = Tatm*(PR)^((g-1)/g); % ideal intake runner temperature [K]
T1 = Tatm - ((Tatm-T1s)/etac); % real intake runner temperature
% intake conditions
P(count,1) = P1; % intake pressure [Pa]
T(count,1) = T1; % intake temperature [K]
M(count) = P(count,1)*Vc(1)/(R*T(count,1)); % charge mass (fuel + air) [kg]
% compression stroke
for count1 = 2:length(CADc(1:181))
P(count,count1) = P(count,count1-1)*(Vc(count1-1)/Vc(count1))^g;
T(count,count1) = P(count,count1)*Vc(count1)/(M(count)*R);
end
% heat addition - happens from 359-360 CAD
Mf(count) = M(count)/(1+FA); % mass of fuel [kg]
Qin(count) = Mf(count)*LHV; % heat in due to fuel, assuming ideal combustion [J]
T(count,182) = T(count,181) + Qin(count)/(M(count)*cv); % temperature from first law [K]
P(count,182) = M(count)*R*T(count,182)/Vc(182); % pressure [Pa]
% explansion stroke
for count2 = 183:length(CADc)
P(count,count2) = P(count,count2-1)*(Vc(count2-1)/Vc(count2))^g;
T(count,count2) = P(count,count2)*Vc(count2)/(M(count)*R);
end
Wcycle(count) = 2.*N(count).*trapz(Vc,P(count,:)); % power from 4 cylinders - 2 power strokes per revolution of the engine[W]
mdot(count) = 2.*M(count).*N(count); % mass flow rate for 4 cylinders - 2 intake strokes per revoluation [kg/s]
Ws(count) = -mdot(count).*cps.*(T1-Tatm); % power from supercharger [W]
Wnet(count) = Wcycle(count) + Ws(count); % net power [W]
eta(count) = Wnet(count)./(2.*Qin(count)*N(count)); % thermal efficiency with two combustion processes per revolution
end
figure
axes('FontSize',14)
hold on
plot(Vc,P,'LineWidth',2)
grid on
legend('No supercharger','Supercharger')
xlabel('Volume [m^3]')
ylabel('Pressure [Pa]')
figure
axes('FontSize',14)
hold on
plot(N,Wnet*2/N,'-o','LineWidth',2)
grid on
xlabel('Supercharger Pressure Ratio')
ylabel('Net Work [J]')
figure
axes('FontSize',14)
hold on
plot(N,eta,'-o','LineWidth',2)
grid on
xlabel('Supercharger Pressure Ratio')
ylabel('Thermal Efficiency')

Answers (1)

Walter Roberson
Walter Roberson on 11 Dec 2020
plot(N,Wnet*2/N,'-o','LineWidth',2)
Wnet and N are both vectors, so you have accidentally invoked matrix division. You need
plot(N,Wnet.*2./N,'-o','LineWidth',2)
However you will still see a straight line. If you look at the way you construct Wnet, you have the sum of two values, with each of the values being a multiple of N, so the N potentially factors out when you do the Wnet ./ N

Categories

Find more on Oil, Gas & Petrochemical in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!