%% Block Parameters - Compression/Expansion Thermodynamics & Combustion
% Two separate polytropic processes are used to model combustion. One for the process before combustion (subscripted with c indicating compression), the other for the process after combustion (subscripted with e indicating expansion). The formulation is largely based on a literature by Lars Eriksson and Lars Nielsen, "Modeling and Control of Engines and Drivelines" (Wiley 2014).
%% Parameters for plotting
% Cylinder pressures at intake valve close (IVC):

clear p_ivc_kPa
p_ivc_kPa(1) = 30;
p_ivc_kPa(2) = 80;

% Cylinder temperature at IVC:

T_ivc_K = 273.15 + 25;

% Spark advance:

clear spa_degBTDC
spa_degBTDC(1) = 30;
spa_degBTDC(2) = 10;
spa_rad = (360 - spa_degBTDC)./180*pi;

% Air-to-Fuel ratio:

afr = 14.7;

% Stoichiometric number
AFR_s = 14.7;  % gasoline

%% Parameter values in Block Parameters window
% Make sure the Block Parameters window is selected before pushing Get button.

 
blk = gcb;

% Bore:

B_mm_blk = eval(get_param(blk, 'B'))

% Stroke:

L_stroke_mm_blk = eval(get_param(blk, 'L_stroke'))

% Compression ratio:

r_c_blk = eval(get_param(blk, 'r_c'))

% Intake valve close angle:

IVC_deg_blk = eval(get_param(blk, 'IVC'))

% Exhaust valve close angle:

EVO_deg_blk = eval(get_param(blk, 'EVO'))

% Lower heating value for gasoline:

q_LHV_kJ_perkg_blk = eval(get_param(blk, 'q_LHV'))

% Constant volume specific heat for compression model:

Cv_c_J_perkg_perK_blk = eval(get_param(blk, 'Cv_c'))

% Specific heat ratio for compression model:

k_c_blk = eval(get_param(blk, 'k_c'))

% Specific heat ratio for expansion model:

k_e_blk = eval(get_param(blk, 'k_e'))

% Burn duration:

comb_bd_degCA_blk = eval(get_param(blk, 'comb_bd'))
comb_bd_rad = comb_bd_degCA_blk/180*pi;

% Wiebe parameters:

comb_a_blk = eval(get_param(blk, 'comb_a'))
comb_m_blk = eval(get_param(blk, 'comb_m'))

%
%% Derived parameters
% Crank radius:

R_crank_mm = L_stroke_mm_blk/2

% Displacement volume:

V_displ_mm3 = pi*(B_mm_blk/2)^2*L_stroke_mm_blk

% Clearance volume:

V_c_mm3 = V_displ_mm3/(r_c_blk - 1)

% Temperature increase due to Otto combustion:

dTcomb_K = q_LHV_kJ_perkg_blk * min(1, afr/AFR_s) ...
    / ((afr + 1)*Cv_c_J_perkg_perK_blk*1e-3)
% Cylinder volume at IVC
IVC_rad = IVC_deg_blk/180*pi;
V_ivc_mm3 = sdl_carem_cylinder_volume_mm3_matlab(...
    IVC_rad, V_c_mm3, R_crank_mm, L_stroke_mm_blk, B_mm_blk);

% Cylinder volume at TDC
V_tdc_mm3 = sdl_carem_cylinder_volume_mm3_matlab(...
    0, V_c_mm3, R_crank_mm, L_stroke_mm_blk, B_mm_blk);

% Cylinder volume at spark angle (Otto-cycle's 2-3 expansion)
% for 2 spark angles
V_23_mm3_sa1 = sdl_carem_cylinder_volume_mm3_matlab(...
    spa_rad(1), V_c_mm3, R_crank_mm, L_stroke_mm_blk, B_mm_blk);
V_23_mm3_sa2 = sdl_carem_cylinder_volume_mm3_matlab(...
    spa_rad(2), V_c_mm3, R_crank_mm, L_stroke_mm_blk, B_mm_blk);

% Compute pressures and temperatures as a function of crank position

% X axis - crank position between IVC and EVO
x_degCA = transpose(linspace(IVC_deg_blk, EVO_deg_blk, 100));
x_rad = x_degCA/180*pi;

% Cylinder volume profile
V_cyl_mm3 = sdl_carem_cylinder_volume_mm3_matlab(...
    x_rad, V_c_mm3, R_crank_mm, L_stroke_mm_blk, B_mm_blk);

% Cylinder pressure right before Otto combustion
% for 2 p_IVC and 2 spark angles:
p_2_kPa_pivc1_sa1 = p_ivc_kPa(1).*(V_ivc_mm3./V_23_mm3_sa1).^k_c_blk;
p_2_kPa_pivc1_sa2 = p_ivc_kPa(1).*(V_ivc_mm3./V_23_mm3_sa2).^k_c_blk;
p_2_kPa_pivc2_sa1 = p_ivc_kPa(2).*(V_ivc_mm3./V_23_mm3_sa1).^k_c_blk;
p_2_kPa_pivc2_sa2 = p_ivc_kPa(2).*(V_ivc_mm3./V_23_mm3_sa2).^k_c_blk;

% Cylinder temperature right before Otto combustion
% for 2 spark angles:
T_2_K_sa1 = T_ivc_K.*(V_ivc_mm3./V_23_mm3_sa1).^(k_c_blk - 1);
T_2_K_sa2 = T_ivc_K.*(V_ivc_mm3./V_23_mm3_sa2).^(k_c_blk - 1);

% Cylinder temperature right after Otto combustion
% for 2 spark angles:
T_3_K_sa1 = T_ivc_K.*(V_ivc_mm3./V_23_mm3_sa1).^(k_c_blk - 1) + dTcomb_K;
T_3_K_sa2 = T_ivc_K.*(V_ivc_mm3./V_23_mm3_sa2).^(k_c_blk - 1) + dTcomb_K;

% Cylinder pressure right after Otto combustion
% for 2 p_IVC and 2 spark angles:
p_3_kPa_pivc1_sa1 = p_2_kPa_pivc1_sa1.*T_3_K_sa1./T_2_K_sa1;
p_3_kPa_pivc1_sa2 = p_2_kPa_pivc1_sa2.*T_3_K_sa2./T_2_K_sa2;
p_3_kPa_pivc2_sa1 = p_2_kPa_pivc2_sa1.*T_3_K_sa1./T_2_K_sa1;
p_3_kPa_pivc2_sa2 = p_2_kPa_pivc2_sa2.*T_3_K_sa2./T_2_K_sa2;

% Compression model cylinder pressure profile
% for 2 p_IVC:
p_c_kPa_pivc1 = p_ivc_kPa(1).*( V_ivc_mm3 ./ V_cyl_mm3 ).^k_c_blk;
p_c_kPa_pivc2 = p_ivc_kPa(2).*( V_ivc_mm3 ./ V_cyl_mm3 ).^k_c_blk;

% Compression model cylinder temperature profile
T_c_K = T_ivc_K .*( V_ivc_mm3 ./ V_cyl_mm3 ).^(k_c_blk - 1);

p_e_kPa_pivc1_sa1 = p_3_kPa_pivc1_sa1.*( V_tdc_mm3 ./ V_cyl_mm3 ).^k_e_blk.*( 1 - (0.3/90).*spa_degBTDC(1) );
p_e_kPa_pivc1_sa2 = p_3_kPa_pivc1_sa2.*( V_tdc_mm3 ./ V_cyl_mm3 ).^k_e_blk.*( 1 - (0.3/90).*spa_degBTDC(2) );
p_e_kPa_pivc2_sa1 = p_3_kPa_pivc2_sa1.*( V_tdc_mm3 ./ V_cyl_mm3 ).^k_e_blk.*( 1 - (0.3/90).*spa_degBTDC(1) );
p_e_kPa_pivc2_sa2 = p_3_kPa_pivc2_sa2.*( V_tdc_mm3 ./ V_cyl_mm3 ).^k_e_blk.*( 1 - (0.3/90).*spa_degBTDC(2) );

T_e_K_sa1 = T_3_K_sa1.*( V_tdc_mm3 ./ V_cyl_mm3 ).^(k_e_blk - 1).*(0.5*(1+0.8)-(1-0.8)/100.*spa_degBTDC(1) );
T_e_K_sa2 = T_3_K_sa2.*( V_tdc_mm3 ./ V_cyl_mm3 ).^(k_e_blk - 1).*(0.5*(1+0.8)-(1-0.8)/100.*spa_degBTDC(2) );

PR_sa1 = 1 - exp( -comb_a_blk.*( (x_rad - spa_rad(1))./comb_bd_rad).^(comb_m_blk + 1) );
PR_sa1(x_rad < spa_rad(1)) = 0;

PR_sa2 = 1 - exp( -comb_a_blk.*( (x_rad - spa_rad(2))./comb_bd_rad).^(comb_m_blk + 1) );
PR_sa2(x_rad < spa_rad(2)) = 0;

p_cyl_pivc1_sa1 = (1 - PR_sa1).*p_c_kPa_pivc1 + PR_sa1.*p_e_kPa_pivc1_sa1;
p_cyl_pivc1_sa2 = (1 - PR_sa2).*p_c_kPa_pivc1 + PR_sa2.*p_e_kPa_pivc1_sa2;
p_cyl_pivc2_sa1 = (1 - PR_sa1).*p_c_kPa_pivc2 + PR_sa1.*p_e_kPa_pivc2_sa1;
p_cyl_pivc2_sa2 = (1 - PR_sa2).*p_c_kPa_pivc2 + PR_sa2.*p_e_kPa_pivc2_sa2;

T_cyl_sa1 = (1 - PR_sa1).*T_c_K + PR_sa1.*T_e_K_sa1;
T_cyl_sa2 = (1 - PR_sa2).*T_c_K + PR_sa2.*T_e_K_sa2;

%
%% Plot
% Compression model: p_c \left(\theta \;\right)=p_3 {\left(\frac{V_{\textrm{ivc}} }{V\left(\theta \;\right)}\right)}^{k_c }, T_c \left(\theta \;\right)=T_3 {\left(\frac{V_{\textrm{ivc}} }{V\left(\theta \;\right)}\right)}^{k_c -1}
% Expansion model*: p_e \left(\theta \;\right)=B_p \;p_3 {\left(\frac{V_{\textrm{ivc}} }{V\left(\theta \;\right)}\right)}^{k_e }, T_e \left(\theta \;\right)=B_T \;T_3 {\left(\frac{V_{\textrm{ivc}} }{V\left(\theta \;\right)}\right)}^{k_e -1}
% Combustion model: \textrm{PR}\left(\theta \;\right)=1-\exp \left(-a{\left(\frac{\theta -\alpha {\;}_{\textrm{spark}} }{\Delta \;\;}\right)}^{m+1} \right)
% Combined cylinder pressure: p_{\textrm{cyl}} \left(\theta \;\right)=\left(1-\textrm{PR}\left(\theta \;\right)\right)p_c \left(\theta \;\right)+\textrm{PR}\left(\theta \;\right)p_e \left(\theta \;\right)
% Combined cylinder temperature: T_{\textrm{cyl}} \left(\theta \;\right)=\left(1-\textrm{PR}\left(\theta \;\right)\right)T_c \left(\theta \;\right)+\textrm{PR}\left(\theta \;\right)T_e \left(\theta \;\right)
% *B_p and B_T in expansion model are introduced to break the symmetry of pressure and temperature profile shapes.
%% Cylinder pressure

subplot(3,1,[1 2])
hold on
p = plot(x_degCA, p_cyl_pivc1_sa1/1e3);
p.LineWidth = 2;
plot(x_degCA, p_e_kPa_pivc1_sa1/1e3)
plot(x_degCA, p_c_kPa_pivc1/1e3)  % plot in MPa, instead of kPa
xlim([IVC_deg_blk, EVO_deg_blk])
title('Cylinder pressure (MPa)')
legend({'Combined' 'Expansion model' 'Compression model'},...
    'Location','northwest', 'NumColumns',1);
grid on

subplot(3,1,3)
hold on
plot(x_degCA, PR_sa1)
title('Combustion model, PR')
xlabel('Crank angle, \theta (deg)')
xlim([IVC_deg_blk, EVO_deg_blk])
grid on

figure;  hold on;  grid on
plot(x_degCA, p_cyl_pivc1_sa1/1e+3)
plot(x_degCA, p_cyl_pivc1_sa2/1e+3)
plot(x_degCA, p_cyl_pivc2_sa1/1e+3)
plot(x_degCA, p_cyl_pivc2_sa2/1e+3)
xlim([IVC_deg_blk, EVO_deg_blk])
title('Cylinder pressure (MPa)')
xlabel('Crank angle, \theta (deg)')
leg = legend({...
    [num2str(p_ivc_kPa(1)) ' kPa, ' num2str(spa_degBTDC(1)) ' deg'],...
    [num2str(p_ivc_kPa(1)) ' kPa, ' num2str(spa_degBTDC(2)) ' deg'],...
    [num2str(p_ivc_kPa(2)) ' kPa, ' num2str(spa_degBTDC(1)) ' deg'],...
    [num2str(p_ivc_kPa(2)) ' kPa, ' num2str(spa_degBTDC(2)) ' deg']},...
    'Location','northwest', 'NumColumns',1);
title(leg, {'Pres. at IVC (kPa)'; 'Spark Angle (degBTDC)'})
hold off

%% Cylinder temperature

figure;  hold on;  grid on
plot(x_degCA, T_cyl_sa1 - 273.15)
plot(x_degCA, T_cyl_sa2 - 273.15)
xlim([IVC_deg_blk, EVO_deg_blk])
title('Cylinder temperature (degC)')
xlabel('Crank angle, \theta (deg)')
leg = legend({...
    [num2str(spa_degBTDC(1)) ' deg'],...
    [num2str(spa_degBTDC(2)) ' deg']},...
    'Location','northwest', 'NumColumns',1);
title(leg, 'Spark Angle (degBTDC)')
hold off

%