%CAC THONG SO BAN DAU
R = 0.508; % dm
D = 0.829; %dm
S = 1.016; %dm
l = 1.44; %dm
Tsn= 10;
lamda = R/l;
mtt = 16.58; % g/cm2
mK = 16.58; % g/cm2
mnp = 12.435; % g/cm2
mA = 0.3*mtt; % g/cm2
mB = 0.7*mtt; % g/cm2
mj = mA + mnp; % g/cm2
mr = mB + mK; % g/cm2
Vd = 0.549; %lit
Sp = (pi*D^2)/4; % dien tich piston
Pa = 0.09; %MN/m
Pr = 0.115; %MN/m
Pz = 8.555; %MN/m
Pc = 2.241; %MN/m
Vc = Vd/(Tsn -1); %dm^3
Vz = 0.9*Vc;
n1 = 1.373;
n2 = 1.231;
n = 6600; %vong/phut
w = (pi*n)/30; %rad/s
Xa = R*(1-cosd(180) + (lamda/4).*(1-cosd(2.*180)));
Va = Xa*Sp + Vc; %(lit)
%Góc xupap nap mo som: 16 do Truoc DCT
%Góc xupap nap dong muon: 50 do Sau DCD
%Góc xupap xa mo som: 46 do Truoc DCD
%Góc xupap xa dong muon: 13 do Sau DCT
%Góc phun nhiên lieu som: 17 do
%HIEU CHINH THAI – NAP
ahc1 = [0 17 46];
phc1 = [Pr 0.1 Pa];
a1 = linspace(0,46,1000);
X1 = R*(1-cosd(a1) + (lamda/4).*(1-cosd(2.*a1)));
V1 = X1*Sp + Vc;
Pl = interp1(ahc1,phc1,a1, 'spline');
J1 = R*(w^2).*(cosd(a1)+lamda.*cosd(2.*a1));
Vt1 = R*w*(sind(a1) + (lamda/2).*sind(2*a1));
%QUA TRINH NAP
a2 = linspace(46,180,1000);
X2 = R*(1-cosd(a2)+(lamda/4)*(1-cosd(2*a2)));
V2 = X2*Sp+Vc;
P2 = linspace(Pa,Pa,1000);%MN/m2
J2 = R*(w^2).*(cosd(a2)+lamda.*cosd(2.*a2));
Vt2 = R*w*(sind(a2) + (lamda/2).*sind(2*a2));
% QUA TRINH NEN
a3 = linspace(180,343,1000);
X3 = R*(1-cosd(a3)+(lamda/4)*(1-cosd(2*a3)));
V3 = X3*Sp+Vc;
P3 = Pa*((Va./V3).^n1);%MN/m2
J3 = R*(w^2).*(cosd(a3)+lamda.*cosd(2.*a3));
Vt3 = R *w*(sind(a3) + (lamda/2).*sind(2*a3));
% HIEU CHINH NEN - CHAY
Pcc =(Pz-Pc)/3 + Pc;
ahc4 = [343 353 363];
phc4 = [max(P3) 2.35 Pcc];
a4 = linspace(343,363,1000);
X4 = R*(1-cosd(a4)+(lamda/4)*(1-cosd(2*a4)));
V4 = X4*Sp+Vc;
P4 = interp1(ahc4,phc4,a4, 'spline');
J4 = R*(w^2).*(cosd(a4)+lamda.*cosd(2.*a4));
Vt4 = R*w*(sind(a4) + (lamda/2).*sind(2*a4));
% Diem z”
X41 = R*(1-cosd(373) + (lamda/4).*(1-cosd(2.*373)));
V41 = X41*Sp + Vc;
% QUA TRINH GIAN NO
a6 = linspace(373,494,1000);
X6 = R*(1-cosd(a6) + (lamda/4).*(1-cosd(2.*a6)));
V6 = X6.*Sp + Vc;
P6 = Pz.*(Vz./V6).^n2;
J6 = R*(w^2).*(cosd(a6)+lamda.*cosd(2.*a6));
Vt6 = R*w*(sind(a6) + (lamda/2).*sind(2*a6));
% HIEU CHINH CHAY - GIAN NO
ahc5=[363 368 373];
phc5=[Pcc 5.91 max(P6)];
a5=linspace(363,373,1000);
X5 = R*(1-cosd(a5) + (lamda/4).*(1-cosd(2.*a5)));
V5 = X5.*Sp + Vc;
P5 = interp1(ahc5,phc5,a5,'spline');
J5 = R*(w^2).*(cosd(a5)+lamda.*cosd(2.*a5));
Vt5 = R*w*(sind(a5) + (lamda/2).*sind(2*a5));
% HIEU CHINH GIAN NO - THAI
% Diem b’
Pb = min(P6);
% Diem b”
Vb = Va;
Pbb =(Pb-Pr)/2+Pr;
ahc7 = [494 540 586];
phc7 = [Pb Pbb Pr];
a7 = linspace(494,586,1000);
X7 = R*(1-cosd(a7) + (lamda/4).*(1-cosd(2.*a7)));
V7 = X7.*Sp + Vc;
P7 = interp1(ahc7,phc7,a7,'spline');
J7 = R*w^2.*(cosd(a7)+lamda.*cosd(2.*a7));
Vt7 = R*w*(sind(a7) + (lamda/2).*sind(2*a7));
% QUA TRINH THAI
a8 = linspace(586,726,1000);
X8 = R*(1-cosd(a8) + (lamda/4).*(1-cosd(2.*a8)));
V8 = X8*Sp + Vc;
P8 = linspace(Pr,Pr,1000);
J8 = R*w^2.*(cosd(a8)+lamda.*cosd(2.*a8));
Vt8 = R*w*(sind(a8) + (lamda/2).*sind(2*a8));
V = [V1 V2 V3 V4 V5 V6 V7 V8];
P = [Pl P2 P3 P4 P5 P6 P7 P8];
a0 = [a1 a2 a3 a4];
ap= [a1 a2 a3 a4 a5 a6 a7 a8];
X = [X1 X2 X3 X4 X5 X6 X7 X8];
Vt = [Vt1 Vt2 Vt3 Vt4 Vt5 Vt6 Vt7 Vt8];
J = [J1 J2 J3 J4 J5 J6 J7 J8];
Pj = -0.00017961*R/10*(pi*(640/3))^2.*(cosd(ap)+lamda.*cosd(2.*ap));
Pkt =(P-0.1);
P1 = Pkt+Pj;
% DO THI P-V
figure(1)
plot(V,P);
title('DO THI P-V');
xlabel('The tich V (dm3)');
ylabel('Ap suat P (MN/dm2)');
grid on
hold on
% DO THI CHUYEN VI PISTON THEO GOC
figure(2)
X0 = [X1 X2 X3 X4];
plot(a0,X0)
title('DO THI CHUYEN VI');
xlabel('Goc quay truc khuyu (do)');
ylabel('Chuyen vi x (dm)') ;
grid on
hold on
% DO THI VAN TOC
figure(3)
Vt0 = [Vt1 Vt2 Vt3 Vt4];
plot(a0,Vt0)
grid on
hold on
title('DO THI VAN TOC PISTON');
xlabel('Goc quay truc khuyu (do)');
ylabel('Van toc v (dm/s) ');
% DO THI GIA TOC
figure(4)
J0 = [J1 J2 J3 J4];
plot(a0,J0)
grid on
hold on
title('DO THI GIA TOC PISTON');
xlabel('Goc quay truc khuyu (do)');
ylabel('Gia toc J (dm/s2) ')
% DO THI Pkt - Pj - P1
figure(5)
plot(ap,Pkt,'r');
hold on
plot(ap,Pj,'g');
hold on
plot(ap,P1,'b');
hold on
title('DO THI Pkt Pj P1');
xlabel('Goc quay truc khuyu (Do)');
ylabel('Pkt (MN/m2) Pj (MN/m2) P1 (MN/m2)');
legend('Pkt','Pj','P1');
grid on
%% xuat bang excel
% Xuất bảng áp suất MCCT
Vxn = linspace(Va, Vc, 50);
Pxn = Pa .* (Va ./ Vxn).^n1;
Vxg = linspace(Vz, Vb, 50);
Pxg = Pz .* (Vz ./ Vxg).^n2;
% Xuất dữ liệu vào file Excel
writematrix([Vxn(:), Pxn(:)], 'qtnen.xlsx');
writematrix([Vxg(:), Pxg(:)], 'qtgianno.xlsx');
% Lập bảng giá trị kết quả tính toán động lực học
writematrix([ap(:), Pkt(:), Pj(:), P1(:)], 'Bangdgluchoc.xlsx');
% Lập bảng giá trị kết quả tính toán động học
writematrix([a0(:), X0(:), Vt0(:), J0(:)], 'donghoc.xlsx');
