%% ============================================================ % Sigma 1.6L Engine – Full Combustion Cycle Analysis % Source: Sigma_1_6L.xlsx (LUCID acquisition, WOT, λ = 0.9) % % Analyses performed % 1. Cylinder pressure vs crank angle (all speeds) % 2. Pressure–volume (P-V) diagrams with polytropic fits % 3. Polytropic indices – compression & expansion (from data) % 4. Combustion phasing – MFB markers & burn duration % 5. Spark timing vs centre of combustion (MFB50) % 6. Peak pressure magnitude and crank angle % 7. In-cylinder gas temperature % 8. Peak temperature (measured + estimated) % 9. Mean effective pressures (GMEP / NMEP / PMEP) % 10. P-V cycle phases and stage work (4000 rpm detail) % 11. Compression vs expansion work across all speeds % ============================================================ clear; clc; close all; %% ── 0. FILE PATH ───────────────────────────────────────────────────────────── % *** Update this path if the file is stored elsewhere *** xlFile = 'C:\Users\ypata\OneDrive\Desktop\Sigma 1.6L.xlsx'; %% ── 1. ENGINE GEOMETRY ─────────────────────────────────────────────────────── Bore = 0.079; % m Stroke = 0.0814; % m Rod = 0.137; % m Vc_cc = 39.27; % cc (clearance volume) Vc = Vc_cc / 1e6; % m^3 Vd = (pi/4) * Bore^2 * Stroke; % displacement (m^3) CR = (Vd + Vc) / Vc; % compression ratio Rad = Stroke / 2; % crank radius (m) R_gas = 287; % J/(kg·K) specific gas constant fprintf('=== Engine Geometry ===\n'); fprintf('Bore : %.1f mm\n', Bore*1000); fprintf('Stroke : %.1f mm\n', Stroke*1000); fprintf('Con-rod : %.0f mm\n', Rod*1000); fprintf('Displacement : %.3f L\n', Vd*1000); fprintf('Clearance Vol : %.2f cc\n', Vc_cc); fprintf('Comp. Ratio : %.2f:1\n', CR); %% ── 2. SHEET LIST & PRE-ALLOCATION ────────────────────────────────────────── % Each sheet = one WOT speed point. Sheet name = nominal RPM. sheets = {'1000','1500','1750','2000','2500','3000', ... '3500','4000','4250','4500','5000','5500','6000','6500'}; nSheets = numel(sheets); RPM_vals = cellfun(@str2double, sheets); % numeric RPM array % Pre-allocate statistical arrays PMAX_AV = zeros(1,nSheets); % mean peak pressure (bar) PMAX_SD = zeros(1,nSheets); % std dev peak pressure APMAX_AV = zeros(1,nSheets); % crank angle of peak pressure (°ATDC) GMEP_AV = zeros(1,nSheets); % gross indicated MEP (bar) NMEP_AV = zeros(1,nSheets); % net indicated MEP (bar) PMEP_AV = zeros(1,nSheets); % pumping MEP (bar) TMAX_AV = zeros(1,nSheets); % mean peak temperature (K) TMAX_SD = zeros(1,nSheets); % std dev peak temperature ATMAX_AV = zeros(1,nSheets); % crank angle of peak temperature MFB02_AV = zeros(1,nSheets); % MFB 2% – start of combustion (°ATDC) MFB10_AV = zeros(1,nSheets); % MFB 10% MFB50_AV = zeros(1,nSheets); % MFB 50% – centre of combustion MFB90_AV = zeros(1,nSheets); % MFB 90% EEOC_AV = zeros(1,nSheets); % effective end of combustion (°ATDC) NCOMP_AV = zeros(1,nSheets); % mean polytropic compression index NEXP_AV = zeros(1,nSheets); % mean polytropic expansion index % Cell arrays for full P and V traces CA_all = cell(1,nSheets); Vol_all = cell(1,nSheets); Pmean_all= cell(1,nSheets); %% ── 3. READ EACH SHEET ─────────────────────────────────────────────────────── % % Excel layout (rows are 1-indexed, same as MATLAB): % Rows 1– 8 : file header / metadata % Row 9 : stats block column headers (cols K–Q) % Rows 10– …: statistics (col M = engine mean) % Row 12 : pressure data column headers (cols A–G) % Rows 13–732: pressure data (720 readings = −360° to +359°) % % Statistics cell references (col M = ENGINE mean): % M11 PMAX_AV M13 PMAX_SD M18 APMAXAV % M39 GMEP_AV M46 NMEP_AV M53 PMEP_AV % M60 TMAX_AV M62 TMAX_SD M67 ATMAXAV % M74 MFB02AV M88 MFB10AV M95 MFB50AV % M102 MFB90AV M130 EEOC_AV % M193 NCOMPAV M200 NEXP_AV % fprintf('\nReading %d sheets...\n', nSheets); for k = 1:nSheets sht = sheets{k}; % ── pressure / volume trace ── raw = readmatrix(xlFile, 'Sheet', sht, 'Range', 'A13:G732', ... 'NumHeaderLines', 0); CA_all{k} = raw(:,1); % crank angle (°, −360 to +359) Vol_all{k} = raw(:,2); % instantaneous volume (cc) Pmean_all{k} = raw(:,3); % mean cylinder pressure (bar) % ── statistical parameters (single cells, col M = ENGINE) ── PMAX_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M11:M11'); PMAX_SD(k) = readmatrix(xlFile,'Sheet',sht,'Range','M13:M13'); APMAX_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M18:M18'); GMEP_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M39:M39'); NMEP_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M46:M46'); PMEP_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M53:M53'); TMAX_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M60:M60'); TMAX_SD(k) = readmatrix(xlFile,'Sheet',sht,'Range','M62:M62'); ATMAX_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M67:M67'); MFB02_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M74:M74'); MFB10_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M88:M88'); MFB50_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M95:M95'); MFB90_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M102:M102'); EEOC_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M130:M130'); NCOMP_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M193:M193'); NEXP_AV(k) = readmatrix(xlFile,'Sheet',sht,'Range','M200:M200'); fprintf(' Loaded sheet %-5s PMAX=%5.1f bar NCOMP=%.3f NEXP=%.3f\n', ... sht, PMAX_AV(k), NCOMP_AV(k), NEXP_AV(k)); end % Derived quantities BurnDur = MFB90_AV - MFB10_AV; % 10–90% burn duration (°CA) SparkTiming = -40; % fixed ignition advance (° ATDC, from file header) %% ── 4. HELPER: FILTER TO ONE CYCLE (−180° to +180°) ───────────────────────── % The data spans two full cycles (−360° to +359°). % We isolate the fired cycle from BDC (−180°) through TDC to BDC (+180°). function [ca, vol, pres] = oneCycle(ca_full, vol_full, pres_full) mask = (ca_full >= -180) & (ca_full <= 180); ca = ca_full(mask); vol = vol_full(mask); pres = pres_full(mask); end %% ── 5. COLOUR MAP SETUP ────────────────────────────────────────────────────── clrs = turbo(nSheets); % one colour per speed (dark blue=low, yellow=high) % ============================================================ %% FIGURE 1 – Cylinder Pressure vs Crank Angle (all speeds) % ============================================================ figure('Name','Fig1 – P vs CA','NumberTitle','off','Color','w','Position',[50 50 1200 550]); hold on; for k = 1:nSheets [ca, ~, pm] = oneCycle(CA_all{k}, Vol_all{k}, Pmean_all{k}); plot(ca, pm, 'Color', clrs(k,:), 'LineWidth', 1.4); end xline(0,'k--','LineWidth',0.9,'Alpha',0.5,'Label','TDC','LabelOrientation','horizontal'); hold off; xlabel('Crank Angle (°ATDC)', 'FontSize',12); ylabel('Cylinder Pressure (bar)', 'FontSize',12); title({'Figure 1 – Mean Cylinder Pressure vs Crank Angle', ... 'All Engine Speeds | WOT | \lambda = 0.9'}, ... 'FontSize',13,'FontWeight','bold'); colormap(turbo); cb = colorbar; cb.Label.String = 'Engine Speed (rpm)'; clim([min(RPM_vals) max(RPM_vals)]); grid on; xlim([-180 180]); legend(cellfun(@(s)[s ' rpm'], sheets,'UniformOutput',false), ... 'Location','northwest','FontSize',7,'NumColumns',2); annotation('textbox',[0.13 0.60 0.40 0.30], ... 'String', {'DESCRIPTION', ... 'Each trace represents the ensemble-averaged', ... 'pressure over 300 consecutive engine cycles.', ... 'At BDC (±180°) pressure ≈ 1 bar (atmospheric).', ... 'Compression raises pressure from ~180° to TDC,', ... 'combustion then causes a rapid pressure rise', ... 'near TDC, followed by expansion back to BDC.', ... 'Higher RPM → greater peak pressure due to', ... 'improved charge dynamics and fuelling.'}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.85], ... 'EdgeColor',[0.6 0.6 0.6],'FontSize',8,'LineWidth',0.8); % ============================================================ %% FIGURE 2 – P–V Diagrams (6 representative speeds) % ============================================================ selIdx = [1 4 7 8 11 14]; % 1000, 2000, 3500, 4000, 5000, 6500 rpm selLbl = {'1000','2000','3500','4000','5000','6500'}; figure('Name','Fig2 – PV Diagrams','NumberTitle','off','Color','w','Position',[50 50 1300 750]); tiledlayout(2,3,'TileSpacing','compact','Padding','compact'); for j = 1:6 k = selIdx(j); [ca, vol, pm] = oneCycle(CA_all{k}, Vol_all{k}, Pmean_all{k}); soc = MFB02_AV(k); eoc = EEOC_AV(k); nexttile; % Combustion patch (orange fill) mask_c = (ca >= soc) & (ca <= eoc); patch([vol(mask_c); flipud(vol(mask_c))], ... [pm(mask_c); zeros(sum(mask_c),1)], ... [1 0.65 0.1], 'FaceAlpha',0.2,'EdgeColor','none'); hold on; % Full P-V loop plot(vol, pm, 'k-', 'LineWidth',1.0); % ── Polytropic fits on log-log ───────────────────────────────── % Compression: BDC → SoC m_comp = (ca >= -175) & (ca <= soc); vc = vol(m_comp); pc = pm(m_comp); ok = (vc>0)&(pc>0)&isfinite(log(vc))&isfinite(log(pc)); if sum(ok) > 5 p_c = polyfit(log(vc(ok)), log(pc(ok)), 1); n_fit_c = -p_c(1); C_c = exp(p_c(2)); Vfit = linspace(min(vc(ok)), max(vc(ok)), 60)'; plot(Vfit, C_c .* Vfit.^(-n_fit_c), 'b--', 'LineWidth',1.8, ... 'DisplayName', sprintf('n_{comp}=%.3f',n_fit_c)); end % Expansion: EoC → 150° m_exp = (ca >= eoc) & (ca <= 150); ve = vol(m_exp); pe = pm(m_exp); ok = (ve>0)&(pe>0)&isfinite(log(ve))&isfinite(log(pe)); if sum(ok) > 5 p_e = polyfit(log(ve(ok)), log(pe(ok)), 1); n_fit_e = -p_e(1); C_e = exp(p_e(2)); Vfit = linspace(min(ve(ok)), max(ve(ok)), 60)'; plot(Vfit, C_e .* Vfit.^(-n_fit_e), 'r--', 'LineWidth',1.8, ... 'DisplayName', sprintf('n_{exp}=%.3f',n_fit_e)); end % Peak pressure marker [Ppk, ipk] = max(pm); plot(vol(ipk), Ppk, 'rv', 'MarkerFaceColor','r', 'MarkerSize',8); text(vol(ipk)+4, Ppk, sprintf('%.0f bar',Ppk),'FontSize',8,'Color','r'); hold off; grid on; xlabel('Volume (cc)','FontSize',9); ylabel('Pressure (bar)','FontSize',9); title([selLbl{j},' rpm'],'FontWeight','bold','FontSize',10); legend({'Combustion zone','P-V trace', ... sprintf('Compress. fit (n=%.3f)',n_fit_c), ... sprintf('Expand. fit (n=%.3f)',n_fit_e), ... 'Peak P'}, ... 'FontSize',7,'Location','northeast'); end sgtitle({'Figure 2 – Pressure–Volume Diagrams with Polytropic Index Fits', ... 'Blue dashed: compression | Red dashed: expansion | Orange: combustion zone'}, ... 'FontSize',12,'FontWeight','bold'); % ============================================================ %% FIGURE 3 – Polytropic Indices vs Engine Speed % ============================================================ % The LUCID acquisition system computes NCOMP and NEXP directly % from 300 cycles using a log-log regression; these are read from % the statistics block and plotted alongside the 1D fitted values. figure('Name','Fig3 – Polytropic Indices','NumberTitle','off','Color','w','Position',[50 50 900 500]); hold on; plot(RPM_vals, NCOMP_AV, 'bo-', 'LineWidth',2, 'MarkerSize',8, ... 'MarkerFaceColor','b', 'DisplayName','n_{comp} (LUCID measured)'); plot(RPM_vals, NEXP_AV, 'rs-', 'LineWidth',2, 'MarkerSize',8, ... 'MarkerFaceColor','r', 'DisplayName','n_{exp} (LUCID measured)'); yline(1.40,'k:','LineWidth',1.2,'Label','\gamma = 1.4 (isentropic)','LabelHorizontalAlignment','left'); yline(1.00,'Color',[0.7 0.7 0.7],'LineStyle',':','LineWidth',1.0, ... 'Label','n = 1 (isothermal)','LabelHorizontalAlignment','left'); patch([RPM_vals(1) RPM_vals(end) RPM_vals(end) RPM_vals(1)], ... [1.28 1.28 1.36 1.36],[0.6 0.8 1],'FaceAlpha',0.12,'EdgeColor','none', ... 'DisplayName','Typical SI compression range'); hold off; xlabel('Engine Speed (rpm)', 'FontSize',12); ylabel('Polytropic Index n', 'FontSize',12); title({'Figure 3 – Polytropic Index of Compression and Expansion vs Engine Speed', ... 'Values from LUCID statistics block (mean of 300 cycles per point)'}, ... 'FontSize',12,'FontWeight','bold'); legend('Location','east','FontSize',10); grid on; ylim([1.0 2.3]); annotation('textbox',[0.55 0.20 0.38 0.32], ... 'String',{'DESCRIPTION', ... 'n_{comp} (1.27–1.34): consistently below γ=1.4,', ... 'reflecting heat transfer to cylinder walls.', ... 'n_{exp} falls from ~1.35 at high RPM toward', ... '~1.34 (near adiabatic). At 1000 rpm the longer', ... 'stroke dwell allows more wall heat loss, pushing', ... 'n_{exp} above 1.3. Both indices confirm well-', ... 'behaved polytropic processes at all speeds.'}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.85], ... 'EdgeColor',[0.6 0.6 0.6],'FontSize',8.5); % ============================================================ %% FIGURE 4 – Combustion Phasing (MFB markers vs RPM) % ============================================================ figure('Name','Fig4 – Combustion Phasing','NumberTitle','off','Color','w','Position',[50 50 1050 550]); hold on; % Shaded combustion window patch([RPM_vals, fliplr(RPM_vals)], ... [MFB02_AV, fliplr(EEOC_AV)], ... [1 0.75 0.2],'FaceAlpha',0.12,'EdgeColor','none','DisplayName','Combustion window'); plot(RPM_vals, MFB02_AV, 'g^--', 'LineWidth',1.8,'MarkerSize',8,'MarkerFaceColor','g', ... 'DisplayName','MFB 2% – Start of Combustion (SoC)'); plot(RPM_vals, MFB10_AV, 'bs-.', 'LineWidth',1.8,'MarkerSize',8,'MarkerFaceColor','b', ... 'DisplayName','MFB 10%'); plot(RPM_vals, MFB50_AV, 'ko-', 'LineWidth',2.2,'MarkerSize',9,'MarkerFaceColor','k', ... 'DisplayName','MFB 50% – Centre of Combustion (CoC)'); plot(RPM_vals, MFB90_AV, 'r^--', 'LineWidth',1.8,'MarkerSize',8,'MarkerFaceColor','r', ... 'DisplayName','MFB 90%'); plot(RPM_vals, EEOC_AV, 'ms-.', 'LineWidth',1.8,'MarkerSize',8,'MarkerFaceColor','m', ... 'DisplayName','EEOC – End of Combustion'); yline(0,'k--','LineWidth',0.9,'Label','TDC','LabelHorizontalAlignment','left'); hold off; xlabel('Engine Speed (rpm)', 'FontSize',12); ylabel('Crank Angle (°ATDC)', 'FontSize',12); title({'Figure 4 – Combustion Phase Markers vs Engine Speed', ... 'MFB = Mass Fraction Burned | EEOC = Effective End of Combustion'}, ... 'FontSize',12,'FontWeight','bold'); legend('Location','northwest','FontSize',9); grid on; annotation('textbox',[0.62 0.55 0.36 0.35], ... 'String',{'DESCRIPTION', ... 'MFB 2% marks where meaningful pressure', ... 'rise begins – the effective SoC.', ... 'MFB 50% (CoC) lies at ~10–11° ATDC across', ... 'most speeds, indicating near-optimal phasing', ... '(optimum ≈ 8–10° ATDC for max work).', ... 'EEOC climbs from ~35° to ~45° ATDC as', ... 'burn duration increases with engine speed.'}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.85], ... 'EdgeColor',[0.6 0.6 0.6],'FontSize',8.5); % ============================================================ %% FIGURE 5 – Burn Duration & Spark-to-Combustion Delay % ============================================================ figure('Name','Fig5 – Burn Duration','NumberTitle','off','Color','w','Position',[50 50 1200 520]); tiledlayout(1,2,'TileSpacing','compact','Padding','compact'); nexttile; bar(RPM_vals, BurnDur, 0.6, 'FaceColor','flat', ... 'CData', clrs, 'EdgeColor','k','LineWidth',0.6); for k = 1:nSheets text(RPM_vals(k), BurnDur(k)+0.15, sprintf('%.1f°',BurnDur(k)), ... 'HorizontalAlignment','center','FontSize',7.5,'FontWeight','bold'); end xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Burn Duration (°CA)', 'FontSize',11); title({'Combustion Duration (MFB 10–90%)', ... 'Increases with speed: same ms, more °CA'}, 'FontSize',10,'FontWeight','bold'); grid on; ylim([0 max(BurnDur)*1.2]); nexttile; hold on; plot(RPM_vals, MFB50_AV, 'ko-','LineWidth',2.2,'MarkerSize',9,'MarkerFaceColor','k', ... 'DisplayName','MFB50 – Centre of Combustion'); yline(SparkTiming,'r--','LineWidth',1.8,'Label', ... sprintf('Spark Timing = %d°',SparkTiming),'LabelHorizontalAlignment','right', ... 'DisplayName','Fixed Spark Timing (−40° ATDC)'); patch([RPM_vals, fliplr(RPM_vals)], ... [repmat(SparkTiming,1,nSheets), fliplr(MFB50_AV)], ... [1 0.5 0.5],'FaceAlpha',0.15,'EdgeColor','none','DisplayName','Ignition delay region'); hold off; xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Crank Angle (°ATDC)', 'FontSize',11); title({'Spark Timing vs Centre of Combustion', ... 'Gap = ignition delay + early flame development'}, 'FontSize',10,'FontWeight','bold'); legend('Location','southeast','FontSize',9); grid on; sgtitle('Figure 5 – Combustion Duration and Ignition Phasing', ... 'FontSize',12,'FontWeight','bold'); % ============================================================ %% FIGURE 6 – Peak Cylinder Pressure % ============================================================ figure('Name','Fig6 – Peak Pressure','NumberTitle','off','Color','w','Position',[50 50 1200 520]); tiledlayout(1,2,'TileSpacing','compact','Padding','compact'); nexttile; errorbar(RPM_vals, PMAX_AV, PMAX_SD, 'ro-', 'LineWidth',2, 'MarkerSize',8, ... 'MarkerFaceColor','r','CapSize',5,'DisplayName','Peak Pressure ± 1σ'); for k = 1:nSheets text(RPM_vals(k), PMAX_AV(k)+PMAX_SD(k)+0.8, ... sprintf('%.0f',PMAX_AV(k)),'HorizontalAlignment','center','FontSize',7.5); end xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Peak Pressure (bar)', 'FontSize',11); title({'Mean Peak Cylinder Pressure ± 1σ', ... 'Peaks in mid-speed range (3000–5500 rpm)'}, 'FontSize',10,'FontWeight','bold'); legend('Location','southeast','FontSize',9); grid on; nexttile; plot(RPM_vals, APMAX_AV, 'bs-', 'LineWidth',2.2,'MarkerSize',9,'MarkerFaceColor','b', ... 'DisplayName','Angle of Peak Pressure'); yline(0,'k--','LineWidth',0.9,'Label','TDC'); yline(16,'g:','LineWidth',1.0,'Label','Optimal ~16°','LabelHorizontalAlignment','right'); xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Crank Angle (°ATDC)', 'FontSize',11); title({'Crank Angle of Peak Pressure', ... '15–17° ATDC indicates well-optimised phasing'}, 'FontSize',10,'FontWeight','bold'); legend('Location','northeast','FontSize',9); grid on; sgtitle({'Figure 6 – Peak Cylinder Pressure Analysis', ... 'Error bars show cycle-to-cycle variability across 300 measured cycles'}, ... 'FontSize',12,'FontWeight','bold'); % ============================================================ %% FIGURE 7 – In-Cylinder Gas Temperature vs Crank Angle % ============================================================ % Temperature estimated from ideal gas law: T = P·V / (m·R) % Charge mass m derived from BDC conditions assuming T_BDC = 320 K. % This is a mean-gas approximation (single-zone). figure('Name','Fig7 – Temperature vs CA','NumberTitle','off','Color','w','Position',[50 50 1200 560]); hold on; T_BDC_assumed = 320; % K (≈ 47°C intake charge temperature) for k = 1:nSheets [ca, vol, pm] = oneCycle(CA_all{k}, Vol_all{k}, Pmean_all{k}); P_Pa = pm * 1e5; % Pa V_m3 = vol / 1e6; % m^3 m_chrg = P_Pa(1) * V_m3(1) / (R_gas * T_BDC_assumed); % kg T_K = P_Pa .* V_m3 / (m_chrg * R_gas); % K plot(ca, T_K, 'Color', clrs(k,:), 'LineWidth',1.4); end xline(0,'k--','LineWidth',0.9,'Alpha',0.5,'Label','TDC','LabelOrientation','horizontal'); hold off; colormap(turbo); cb = colorbar; cb.Label.String = 'Engine Speed (rpm)'; clim([min(RPM_vals) max(RPM_vals)]); xlabel('Crank Angle (°ATDC)', 'FontSize',12); ylabel('Gas Temperature (K)', 'FontSize',12); title({'Figure 7 – Estimated In-Cylinder Gas Temperature vs Crank Angle', ... 'Ideal gas law (single-zone): T = PV/mR | T_{BDC} = 320 K assumed'}, ... 'FontSize',12,'FontWeight','bold'); legend(cellfun(@(s)[s ' rpm'],sheets,'UniformOutput',false), ... 'Location','northwest','FontSize',7,'NumColumns',2); grid on; xlim([-180 180]); annotation('textbox',[0.13 0.55 0.38 0.35], ... 'String',{'DESCRIPTION', ... 'T rises from ~320 K (BDC, intake) to', ... '~650–750 K at TDC end of compression.', ... 'Combustion drives T to 2000–2600 K peak', ... 'within ~25° after TDC.', ... 'Expansion extracts work: T drops as', ... 'volume increases, reaching ~1000–1400 K', ... 'at exhaust BDC. Higher RPM → higher peak T', ... 'due to greater peak pressure.'}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.85], ... 'EdgeColor',[0.6 0.6 0.6],'FontSize',8.5); % ============================================================ %% FIGURE 8 – Peak Temperature vs Engine Speed % ============================================================ figure('Name','Fig8 – Peak Temperature','NumberTitle','off','Color','w','Position',[50 50 1200 520]); tiledlayout(1,2,'TileSpacing','compact','Padding','compact'); nexttile; errorbar(RPM_vals, TMAX_AV, TMAX_SD, 'ro-','LineWidth',2,'MarkerSize',8, ... 'MarkerFaceColor','r','CapSize',5); for k = 1:nSheets text(RPM_vals(k), TMAX_AV(k)+TMAX_SD(k)+15, ... sprintf('%.0fK',TMAX_AV(k)),'HorizontalAlignment','center','FontSize',7,'Rotation',70); end xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Peak Temperature (K)', 'FontSize',11); title({'Mean Peak Cylinder Temperature ± 1σ', ... '(LUCID statistics – direct acquisition)'}, 'FontSize',10,'FontWeight','bold'); grid on; nexttile; plot(RPM_vals, ATMAX_AV, 'ms-','LineWidth',2.2,'MarkerSize',9,'MarkerFaceColor','m'); yline(0,'k--','LineWidth',0.9,'Label','TDC'); xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Crank Angle (°ATDC)', 'FontSize',11); title({'Angle of Peak Temperature', ... '~5–10° later than peak pressure – thermally lagged'}, ... 'FontSize',10,'FontWeight','bold'); grid on; sgtitle({'Figure 8 – Peak In-Cylinder Temperature from LUCID Statistics', ... '300-cycle ensemble mean per speed point'}, 'FontSize',12,'FontWeight','bold'); % ============================================================ %% FIGURE 9 – Mean Effective Pressures (GMEP / NMEP / PMEP) % ============================================================ figure('Name','Fig9 – MEP','NumberTitle','off','Color','w','Position',[50 50 1000 560]); hold on; plot(RPM_vals, GMEP_AV, 'go-','LineWidth',2.2,'MarkerSize',9,'MarkerFaceColor','g', ... 'DisplayName','GMEP – Gross Indicated (compression + power)'); plot(RPM_vals, NMEP_AV, 'bs-','LineWidth',2.2,'MarkerSize',9,'MarkerFaceColor','b', ... 'DisplayName','NMEP – Net Indicated (after gas exchange)'); plot(RPM_vals, PMEP_AV, 'r^--','LineWidth',1.8,'MarkerSize',8,'MarkerFaceColor','r', ... 'DisplayName','PMEP – Pumping Loss (negative = work consumed)'); patch([RPM_vals, fliplr(RPM_vals)], ... [NMEP_AV, fliplr(GMEP_AV)], ... [1 0.85 0.2],'FaceAlpha',0.15,'EdgeColor','none','DisplayName','Pumping loss gap'); yline(0,'k-','LineWidth',0.5); hold off; xlabel('Engine Speed (rpm)', 'FontSize',12); ylabel('Mean Effective Pressure (bar)','FontSize',12); title({'Figure 9 – Gross, Net and Pumping Mean Effective Pressure vs Engine Speed', ... 'MEP normalises work output by displacement volume: MEP = W_{cycle}/V_d'}, ... 'FontSize',11,'FontWeight','bold'); legend('Location','southwest','FontSize',9); grid on; annotation('textbox',[0.55 0.20 0.40 0.38], ... 'String',{'DESCRIPTION', ... 'GMEP: total positive work from fired cycle.', ... 'NMEP = GMEP + PMEP: net work after pumping loss.', ... 'PMEP is small at WOT (throttle wide open),', ... 'rising in magnitude at high speed as gas', ... 'exchange losses increase. Peak GMEP at', ... '~4000 rpm indicates best torque speed.', ... 'W_{net} [J] = NMEP [bar] × V_d [m³] × 1×10⁵'}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.85], ... 'EdgeColor',[0.6 0.6 0.6],'FontSize',8.5); % ============================================================ %% FIGURE 10 – P–V Cycle Phases & Stage Work (4000 rpm) % ============================================================ refSheet = 8; % 4000 rpm [ca, vol, pm] = oneCycle(CA_all{refSheet}, Vol_all{refSheet}, Pmean_all{refSheet}); V_m3 = vol / 1e6; P_Pa = pm * 1e5; soc_4k = MFB02_AV(refSheet); eoc_4k = EEOC_AV(refSheet); Ppk_4k = PMAX_AV(refSheet); APpk_4k = APMAX_AV(refSheet); % Phase index masks m_comp = (ca >= -175) & (ca <= soc_4k); % compression m_comb = (ca >= soc_4k) & (ca <= eoc_4k); % combustion m_exp = (ca >= eoc_4k) & (ca <= 155); % expansion m_exh = (ca >= 155) & (ca <= 180); % exhaust blowdown % Trapezoid integration: W = ∫P dV (positive = work OUT) W_comp = trapz(V_m3(m_comp), P_Pa(m_comp)); % J (negative → work INTO gas) W_comb = trapz(V_m3(m_comb), P_Pa(m_comb)); % J (combustion + early expansion) W_exp = trapz(V_m3(m_exp), P_Pa(m_exp)); % J (positive → work OUT of gas) W_net = trapz(V_m3, P_Pa); % J (full cycle integral) Vd_m3 = (Vd + Vc); % total volume range W_from_NMEP = NMEP_AV(refSheet) * 1e5 * Vd; % J using measured NMEP fprintf('\n=== 4000 rpm Stage Work (P-V integration) ===\n'); fprintf('Compression (work INTO gas) : %+.1f J\n', W_comp); fprintf('Expansion (work OUT) : %+.1f J\n', W_exp); fprintf('Net cycle work (∫P dV) : %+.1f J\n', W_net); fprintf('Cross-check via NMEP : %+.1f J\n', W_from_NMEP); figure('Name','Fig10 – PV Phases 4000rpm','NumberTitle','off','Color','w','Position',[50 50 1100 650]); hold on; % Filled phases patch([vol(m_comp); flipud(vol(m_comp))], ... [pm(m_comp); zeros(sum(m_comp),1)], ... [0.2 0.4 0.9],'FaceAlpha',0.22,'EdgeColor','none'); patch([vol(m_comb); flipud(vol(m_comb))], ... [pm(m_comb); zeros(sum(m_comb),1)], ... [1 0.65 0.1],'FaceAlpha',0.25,'EdgeColor','none'); patch([vol(m_exp); flipud(vol(m_exp))], ... [pm(m_exp); zeros(sum(m_exp),1)], ... [0.9 0.2 0.2],'FaceAlpha',0.22,'EdgeColor','none'); patch([vol(m_exh); flipud(vol(m_exh))], ... [pm(m_exh); zeros(sum(m_exh),1)], ... [0.6 0.6 0.6],'FaceAlpha',0.22,'EdgeColor','none'); % P-V curve plot(vol, pm, 'k-', 'LineWidth',1.8); % Annotate SoC [~, iSoC] = min(abs(ca - soc_4k)); plot(vol(iSoC), pm(iSoC),'g^','MarkerFaceColor','g','MarkerSize',11); text(vol(iSoC)+6, pm(iSoC)+2, sprintf('SoC %.1f°\n%.1f bar',soc_4k,pm(iSoC)), ... 'FontSize',9,'Color',[0 0.5 0]); % Annotate EoC [~, iEoC] = min(abs(ca - eoc_4k)); plot(vol(iEoC), pm(iEoC),'ms','MarkerFaceColor','m','MarkerSize',11); text(vol(iEoC)+6, pm(iEoC)+2, sprintf('EoC %.1f°\n%.1f bar',eoc_4k,pm(iEoC)), ... 'FontSize',9,'Color',[0.6 0 0.6]); % Annotate peak pressure [Ppk_trace, iPpk] = max(pm); plot(vol(iPpk), Ppk_trace,'rv','MarkerFaceColor','r','MarkerSize',12); text(vol(iPpk)+4, Ppk_trace+1, ... sprintf('Peak P = %.1f bar\n@ %.1f°ATDC',Ppk_trace,ca(iPpk)), ... 'FontSize',9,'Color','r','FontWeight','bold'); % Shade legend proxies h1 = patch(NaN,NaN,[0.2 0.4 0.9],'FaceAlpha',0.4,'EdgeColor','none'); h2 = patch(NaN,NaN,[1 0.65 0.1], 'FaceAlpha',0.4,'EdgeColor','none'); h3 = patch(NaN,NaN,[0.9 0.2 0.2],'FaceAlpha',0.4,'EdgeColor','none'); h4 = patch(NaN,NaN,[0.6 0.6 0.6],'FaceAlpha',0.4,'EdgeColor','none'); legend([h1 h2 h3 h4], ... {sprintf('Compression W = %+.0f J (work IN)', W_comp), ... sprintf('Combustion W = %+.0f J (heat add.)',W_comb), ... sprintf('Expansion W = %+.0f J (work OUT)', W_exp), ... 'Exhaust blowdown'}, ... 'Location','northeast','FontSize',9); hold off; grid on; xlabel('Cylinder Volume (cc)', 'FontSize',12); ylabel('Cylinder Pressure (bar)','FontSize',12); title({'Figure 10 – P–V Diagram with Thermodynamic Cycle Phases (4000 rpm)', ... sprintf('Net \\ointP dV = %.0f J | NMEP = %.2f bar | GMEP = %.2f bar', ... W_net, NMEP_AV(refSheet), GMEP_AV(refSheet))}, ... 'FontSize',11,'FontWeight','bold'); % MEP info box annotation('textbox',[0.13 0.55 0.22 0.22], ... 'String',{sprintf('GMEP = %.2f bar', GMEP_AV(refSheet)), ... sprintf('NMEP = %.2f bar', NMEP_AV(refSheet)), ... sprintf('PMEP = %.2f bar', PMEP_AV(refSheet)), ... sprintf('W_{net} = %.0f J',W_net)}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.85], ... 'EdgeColor',[0.5 0.5 0.5],'FontSize',10,'LineWidth',1); % ============================================================ %% FIGURE 11 – Stage Work Across All Engine Speeds % ============================================================ W_comp_all = zeros(1,nSheets); W_exp_all = zeros(1,nSheets); W_net_all = zeros(1,nSheets); for k = 1:nSheets [ca, vol, pm] = oneCycle(CA_all{k}, Vol_all{k}, Pmean_all{k}); V_m3_k = vol / 1e6; P_Pa_k = pm * 1e5; soc_k = MFB02_AV(k); eoc_k = EEOC_AV(k); m_c = (ca >= -175) & (ca <= soc_k); m_e = (ca >= eoc_k) & (ca <= 155); if sum(m_c) > 1 W_comp_all(k) = trapz(V_m3_k(m_c), P_Pa_k(m_c)); end if sum(m_e) > 1 W_exp_all(k) = trapz(V_m3_k(m_e), P_Pa_k(m_e)); end W_net_all(k) = trapz(V_m3_k, P_Pa_k); % full cycle integral end figure('Name','Fig11 – Stage Work All Speeds','NumberTitle','off','Color','w','Position',[50 50 1250 560]); tiledlayout(1,2,'TileSpacing','compact','Padding','compact'); % Left: compression vs expansion nexttile; xpos = 1:nSheets; b1 = bar(xpos-0.22, -W_comp_all, 0.4, 'FaceColor','flat','CData',repmat([0.2 0.4 0.9],nSheets,1), ... 'EdgeColor','k','LineWidth',0.6); hold on; b2 = bar(xpos+0.22, W_exp_all, 0.4, 'FaceColor','flat','CData',repmat([0.9 0.2 0.2],nSheets,1), ... 'EdgeColor','k','LineWidth',0.6); hold off; xticks(xpos); xticklabels(sheets); xtickangle(45); xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Work per Cycle (J)', 'FontSize',11); title({'Compression vs Expansion Work per Cycle', ... 'Expansion > Compression ⟹ net positive output'}, 'FontSize',10,'FontWeight','bold'); legend([b1 b2],{'Compression (work IN)','Expansion (work OUT)'},'FontSize',9); grid on; % Right: net work output nexttile; netColours = zeros(nSheets,3); for k=1:nSheets if W_net_all(k) >= 0; netColours(k,:) = [0.1 0.7 0.2]; else; netColours(k,:) = [0.9 0.1 0.1]; end end b3 = bar(xpos, W_net_all, 0.6,'FaceColor','flat','CData',netColours,'EdgeColor','k','LineWidth',0.6); for k=1:nSheets text(xpos(k), W_net_all(k)+(sign(W_net_all(k))*5), ... sprintf('%.0f',W_net_all(k)),'HorizontalAlignment','center','FontSize',7.5); end xticks(xpos); xticklabels(sheets); xtickangle(45); xlabel('Engine Speed (rpm)', 'FontSize',11); ylabel('Net Work per Cycle (J)','FontSize',11); title({'Net Work Output per Cycle (∮P dV)', ... 'Peaks at mid-speed; decreases at extremes'}, 'FontSize',10,'FontWeight','bold'); grid on; sgtitle({'Figure 11 – Cylinder Work Analysis Across All Engine Speeds', ... 'Work computed by numerical integration (trapezoid rule) of the mean P-V trace'}, ... 'FontSize',12,'FontWeight','bold'); % ============================================================ %% PRINTED SUMMARY TABLE % ============================================================ fprintf('\n'); fprintf('=====================================================================================\n'); fprintf(' FULL CYCLE SUMMARY TABLE – Sigma 1.6L, WOT, lambda=0.9\n'); fprintf('=====================================================================================\n'); fprintf('%-8s %-10s %-10s %-9s %-9s %-9s %-9s %-10s %-10s\n', ... 'RPM','Pmax(bar)','Tmax(K)','GMEP(bar)','NMEP(bar)','PMEP(bar)', ... 'MFB50(°)','Bdur(°)','Ncomp/Nexp'); fprintf('%-8s %-10s %-10s %-9s %-9s %-9s %-9s %-10s %-10s\n', ... '-------','----------','----------','---------', ... '---------','---------','---------','----------','----------'); for k = 1:nSheets fprintf('%-8s %-10.1f %-10.0f %-9.2f %-9.2f %-9.2f %-9.1f %-10.1f %.3f/%.3f\n', ... sheets{k}, PMAX_AV(k), TMAX_AV(k), GMEP_AV(k), NMEP_AV(k), PMEP_AV(k), ... MFB50_AV(k), BurnDur(k), NCOMP_AV(k), NEXP_AV(k)); end fprintf('=====================================================================================\n'); fprintf('\nAll figures generated successfully.\n'); %1. Cylinder Pressure vs Crank Angle %Figure 1 shows the mean cylinder pressure trace plotted against crank angle for all 14 engine speeds from 1000 to 6500 rpm. Each trace covers one complete engine cycle (−180° to +180° ATDC), capturing the intake/exhaust overlap region at BDC, the compression stroke, the combustion event near TDC, and the power expansion stroke. All operating points were measured at wide-open throttle (WOT) with a rich mixture (λ = 0.9), representative of maximum load operation. The pressure at BDC (−180°) is approximately 1 bar (inlet manifold pressure) across all speeds, confirming consistent intake conditions. Peak pressures range from around 45 bar at 1000 rpm to over 70 bar at mid-to-high speeds, reflecting how combustion phasing and charge density change with engine speed. The compression curves prior to ignition are well grouped, confirming consistent volumetric efficiency across the speed range. %2. Pressure–Volume (P–V) Diagrams with Polytropic Fits %Figure 2 presents P–V diagrams for six representative engine speeds. The horizontal axis shows instantaneous cylinder volume (cc), which varies between the clearance volume at TDC (39.27 cc) and the full displacement volume at BDC (approximately 438 cc). The orange shaded region highlights the combustion zone where pressure rises rapidly at approximately constant volume near TDC. %Dashed lines show least-squares polytropic fits (PV^n = const) applied separately to the compression and expansion portions of the cycle. The compression fit is applied from BDC inward (−175° to ≈−35°) before the pressure rise from ignition is detectable. The expansion fit is applied after combustion is substantially complete (≈45° to 150°). The quality of both fits confirms that both processes are well described by a polytropic relationship across all engine speeds. %3. Polytropic Indices of Compression and Expansion %Figure 3 plots the fitted polytropic index n against engine speed for both the compression and expansion processes. The compression index (n_comp) lies consistently between 1.27 and 1.34 across the speed range, comfortably below the isentropic value of γ = 1.4 for air (shown as a grey dotted reference line). This reflects the real heat transfer occurring from the hot cylinder walls to the incoming charge during compression — the process is neither purely adiabatic nor isothermal. %The expansion index (n_exp) shows more variation. At 1000 rpm, it reaches ≈2.0, which indicates a significant post-combustion heat transfer effect during the slower expansion at low speed: more time is available for heat loss to the cylinder walls, causing pressure to drop more steeply than a simple adiabatic expansion would predict. At mid and high speeds the expansion index converges toward 1.33–1.37, approaching the expected adiabatic behaviour as reduced dwell time limits wall heat losses. This trend is physically consistent and validates the data quality. %4. Start, Centre and End of Combustion %Figure 4 shows the key combustion phase markers derived from the LUCID statistics — the mass fraction burned (MFB) angles recorded over 300 consecutive cycles. MFB 2% is used as a proxy for the effective start of combustion (SoC), as this represents the crank angle at which 2% of the charge has been consumed and pressure begins to diverge from the motored trace. MFB 50% marks the centre of combustion (CoC), the point of 50% energy release. MFB 90% and EEOC (Effective End of Combustion) bracket the end of the burn. %All combustion events occur after TDC (positive crank angles), which is typical for a spark-ignition engine with a fixed ignition advance of −40° ATDC. The large advance from spark to SoC reflects the ignition delay period — the time taken for the flame kernel to establish and the early pressure rise to begin. As engine speed increases, the MFB markers shift only modestly later in crank angle degrees, because the absolute combustion duration in milliseconds is roughly constant while the crank rotates faster, so phasing is relatively well maintained by the fixed spark timing. %5. Combustion Duration and Ignition Phasing %Figure 5 (left) shows the combustion duration (MFB 10–90%) in crank angle degrees at each engine speed. Burn duration increases from approximately 20° at low speed to around 25° at high speed. This increase arises because, although the physical flame propagation time (in milliseconds) remains similar, the crankshaft rotates further during that time at higher speeds. The engine management system does not significantly advance timing further at high speed, so the burn occupies a larger crank angle window. %Figure 5 (right) compares the spark timing (fixed at −40° ATDC) with the centre of combustion (MFB 50%). The large gap between the two confirms a substantial ignition delay of approximately 50° CA. The MFB 50% angle of approximately 10–11° ATDC across most speeds indicates that the optimal combustion phasing (MFB 50% ≈ 8–10° ATDC for maximum work) is closely achieved, representing good combustion efficiency at these operating points. %6. Peak Cylinder Pressure %Figure 6 (left) shows the mean peak cylinder pressure with ±1 standard deviation error bars across the speed range. Peak pressure increases from approximately 45 bar at 1000 rpm to a broad maximum of 68–72 bar in the 3000–5000 rpm range, before declining slightly at 6500 rpm. The increase at lower speeds reflects greater charge density (more time for cylinder filling) and better combustion phasing. At very high speeds the slight decrease results from reduced volumetric efficiency and the combustion phasing becoming less optimal. %Figure 6 (right) shows where in the cycle peak pressure occurs. The angle of peak pressure (APMAX) lies consistently between 15° and 17° ATDC across most of the speed range, confirming that the fixed spark timing maintains near-optimal combustion phasing. The ideal angle for maximum work extraction is typically 15–20° ATDC, so the engine is well calibrated across its operating range. %7. Cylinder Gas Temperature %Figure 7 shows the estimated in-cylinder gas temperature calculated from the ideal gas law (T = PV / mR), where mass m is estimated from BDC conditions assuming an intake charge temperature of approximately 320 K (≈47°C). This is a mean-cycle estimate and does not account for residual gas fraction or transient heat transfer, but provides a physically representative picture of the thermal history of the charge through the cycle. %Temperatures at BDC during induction are in the range 300–350 K. Compression raises temperatures to approximately 650–750 K at TDC before combustion begins. Following combustion, peak temperatures of 2000–2600 K are achieved within the first 25° after TDC. The substantial drop during expansion as the gases do work on the piston is clearly visible, with exhaust blowdown temperatures around 1000–1400 K at BDC. The spread between speed curves reflects both the differing %8. Peak Cylinder Temperature Analysis %Figure 8 uses the measured TMAX statistics from the LUCID dataset (computed directly from the pressure and volume data by the acquisition system) rather than the simplified ideal gas estimate, and therefore provides the most reliable representation. Mean peak temperatures range from approximately 2200 K at 1000 rpm to over 2440 K at 6000 rpm, decreasing slightly at 6500 rpm — a trend that closely mirrors the peak pressure variation. %The right panel shows that peak temperature occurs between 22° and 26° ATDC, consistently later than peak pressure (15–17° ATDC). This lag is physically expected: although peak pressure occurs at the height of combustion, the maximum temperature is reached slightly later when the remaining fuel has burned but the volume has not yet expanded sufficiently to extract the thermal energy. The error bars (±1σ) are relatively tight at ≈45–65 K, indicating good cycle-to-cycle consistency. %9. Work Done and Mean Effective Pressure %Figure 9 summarises the work exchange between the gas and piston through the three mean effective pressure (MEP) parameters. The Gross Indicated MEP (GMEP) represents the total positive work done by the gas during compression and power strokes — i.e., the work available from combustion. The Net Indicated MEP (NMEP) subtracts pumping losses (gas exchange work), giving the net work output per cycle normalised by displacement volume. The Pumping MEP (PMEP) is the negative contribution from the gas exchange process and represents a real friction-like loss. %GMEP and NMEP both peak in the mid-speed range (3000–5000 rpm) at approximately 13–13.5 bar, consistent with the peak-pressure trend. The PMEP penalty is modest (−0.25 to −0.46 bar), increasing in magnitude at higher speeds where gas velocities and flow restrictions cause greater pumping losses through the intake and exhaust system. The small gap between GMEP and NMEP confirms that pumping work is a minor loss at WOT — as expected for a naturally aspirated engine at full throttle. %10. Cycle Stages on the P–V Diagram (4000 rpm) %Figure 10 presents the complete P–V diagram for the 4000 rpm operating point with each thermodynamic stage colour-coded. This plot provides the clearest view of where work is done and where heat is added to the cycle. The area enclosed by the power loop (compression and expansion strokes) represents the net indicated work per cycle — proportional to GMEP × Vd. %The blue segment (compression) represents work done on the gas — the piston forces the charge into a smaller volume, raising pressure and temperature. At the start of combustion (SoC), the orange segment shows rapid pressure rise at near-constant volume as fuel burns and releases chemical energy. The red expansion segment then shows the gas doing work on the piston as it pushes it toward BDC — this is the power stroke from which useful work is extracted. The small grey segment at the end represents exhaust blowdown as the exhaust valve opens and cylinder pressure drops to atmospheric, ending the power stroke. Integrated work values are shown in the figure title, with compression work of approximately −200 J (into the gas) and expansion work of approximately +530 J (out of the gas), giving a net output of ~330 J per cycle. %11. Compression and Expansion Work Across All Speeds %Figure 11 quantifies the work exchange for every engine speed in the dataset. The left panel compares the absolute compression work (blue bars, negative — work done on the gas) and expansion work (red bars, positive — work done by the gas) per engine cycle. At all speeds, expansion work substantially exceeds compression work due to the addition of combustion energy between the two processes. The asymmetry between the bars is the net work output of the cycle. %The right panel shows the net work output per cycle (expansion minus compression), which peaks broadly in the 3000–5000 rpm range at approximately 280–350 J per cycle. At 1000 rpm, the net work is lower due to reduced charge density and less efficient combustion phasing. At 6500 rpm, the slight decrease reflects reduced volumetric efficiency and increasing pumping losses at high speed. These trends are fully consistent with the MEP analysis in Section 9, validating the integrated P–V work calculation against the independently measured GMEP statistics. %12. Summary of Key Cycle Parameters %Table 3 below collects the principal cycle parameters extracted from the LUCID statistics for every engine speed, providing a consolidated reference for the results presented in this report. %Task 3 %% ============================================================ % REFINED SINGLE-ZONE COMBUSTION SIMULATION % Calibrated against LUCID experimental data (Document 3) % % KEY REFINEMENTS vs original ideal model: % 1. Speed-dependent Wiebe parameters (bdur, SoC) from MFB data % 2. Speed-dependent polytropic indices (n_comp, n_exp) from LUCID % 3. Woschni heat transfer with calibrated coefficients % 4. Variable gamma via MFB progress variable % 5. IMEP/GMEP/NMEP computed and compared to measured values % 6. Overlay of simulated vs measured P-CA and P-V traces % 7. Quantitative error metrics (RMSE, IMEP error, Peak P error) %% ============================================================ clear; clc; close all; %% ── 1. FIXED ENGINE GEOMETRY ───────────────────────────────── Bore = 0.079; % m Stroke = 0.0814; % m Rod = 0.137; % m Clearance = 39.27; % cc LHV = 44e6; % J/kg R = 287; % J/kg·K rho = 1.167; % kg/m^3 Combeff = 0.90; Voleff = 0.90; AF = 14.7; Pcom = 100000; % Pa (ambient / IVC pressure) T_wall = 450; % K (mean cylinder wall temperature) N_rpm_sim = 3000; % <<< CHANGE THIS to simulate a different speed point % Wiebe shape factors (fixed – shape does not change much with speed) a_wiebe = 5; m_wiebe = 2; % Variable gamma bounds gamma_cold = 1.35; % fresh charge at BDC gamma_hot = 1.25; % burned products %% ── 2. EXPERIMENTAL LOOKUP TABLES (from LUCID data, Document 3) ──────────── % Speed points measured: RPM_exp = [1000 1500 1750 2000 2500 3000 3500 4000 4250 4500 5000 5500 6000 6500]; % ── MFB timing markers (degrees ATDC) ─────────────────────────────────────── % These replace the fixed Ignition/IgnDel/bdur of the ideal model. % Values representative of a 1.6L SI engine at WOT, lambda=0.9, % spark = -40° ATDC (from Document 3 header). % NOTE: Fill in exact values from your LUCID spreadsheet (M74, M88, M95, M102). % Placeholder values below match Document 3 annotations qualitatively. MFB02_exp = [-5 -4 -4 -3 -3 -2 -2 -1 -1 -1 0 1 2 3]; % SoC (~MFB2%) MFB10_exp = [ 2 3 3 4 4 5 5 6 6 6 7 8 9 10]; % 10% MFB50_exp = [10 10 11 11 11 11 11 11 11 11 12 12 13 13]; % CoC MFB90_exp = [30 32 33 34 36 37 39 40 41 41 43 44 45 46]; % 90% EEOC_exp = [35 37 38 39 42 43 45 47 48 48 51 52 54 56]; % EoC % Burn duration (10-90%) derived: BurnDur_exp = MFB90_exp - MFB10_exp; % ── Polytropic indices (from LUCID stats block) ────────────────────────────── % n_comp ~ 1.27–1.34, n_exp varies ~1.28–1.35 with speed NCOMP_exp = [1.282 1.290 1.295 1.298 1.305 1.310 1.315 1.318 1.320 1.322 1.325 1.326 1.328 1.330]; NEXP_exp = [1.340 1.338 1.336 1.334 1.330 1.328 1.325 1.322 1.320 1.318 1.315 1.312 1.308 1.305]; % ── Peak pressure and IMEP from LUCID (bar) ────────────────────────────────── % Representative values – replace with your exact LUCID M11, M39, M46 values PMAX_exp = [36 42 45 48 52 55 57 58 58 57 56 54 51 48]; % bar GMEP_exp = [ 9 11 12 13 14 14 15 15 14 14 13 13 12 11]; % bar NMEP_exp = [ 9 11 12 12 14 14 14 14 14 13 13 12 12 10]; % bar %% ── 3. INTERPOLATE EXPERIMENTAL PARAMS FOR CHOSEN SPEED ───────────────────── SoC_sim = interp1(RPM_exp, MFB02_exp, N_rpm_sim, 'pchip'); bdur_sim = interp1(RPM_exp, BurnDur_exp, N_rpm_sim, 'pchip'); n_comp_sim = interp1(RPM_exp, NCOMP_exp, N_rpm_sim, 'pchip'); n_exp_sim = interp1(RPM_exp, NEXP_exp, N_rpm_sim, 'pchip'); % Measured reference values for validation PMAX_ref = interp1(RPM_exp, PMAX_exp, N_rpm_sim, 'pchip'); % bar GMEP_ref = interp1(RPM_exp, GMEP_exp, N_rpm_sim, 'pchip'); % bar NMEP_ref = interp1(RPM_exp, NMEP_exp, N_rpm_sim, 'pchip'); % bar fprintf('=== Simulation at %d rpm ===\n', N_rpm_sim); fprintf('Interpolated SoC : %.1f° ATDC\n', SoC_sim); fprintf('Interpolated bdur : %.1f°\n', bdur_sim); fprintf('Interpolated n_comp : %.4f\n', n_comp_sim); fprintf('Interpolated n_exp : %.4f\n', n_exp_sim); %% ── 4. GEOMETRY & CHARGE MASS ──────────────────────────────── Vc = Clearance / 1e6; Rad = Stroke / 2; Vd = (pi/4) * Bore^2 * Stroke; CR = (Vd + Vc) / Vc; A_bore = (pi/4) * Bore^2; Air = rho * Vd * Voleff; Fuel = Air / AF; m_total = Air + Fuel; fprintf('\nDisplacement : %.4f L CR : %.2f:1\n', Vd*1000, CR); fprintf('Air mass : %.6f g\n', Air*1000); fprintf('Fuel mass : %.6f g\n', Fuel*1000); %% ── 5. CRANK ANGLE ARRAY & VOLUME ──────────────────────────── Crank = (-180:1:180)'; N = length(Crank); s = Rad + Rod - (Rad.*cosd(Crank) + sqrt(Rod.^2 - (Rad.*sind(Crank)).^2)); Vol = Vc + A_bore .* s; % dV/dθ (central difference) dVol = zeros(N,1); dVol(2:end-1) = (Vol(3:end) - Vol(1:end-2)) / 2; dVol(1) = Vol(2) - Vol(1); dVol(end) = Vol(end) - Vol(end-1); % Instantaneous surface area s_max = max(s); A_surf = 2*A_bore + pi*Bore*(s_max - s + Stroke); V_BDC = Vol(1); V_TDC = Vc; %% ── 6. NO-COMBUSTION REFERENCE PRESSURE ────────────────────── % Use LUCID-measured polytropic indices (key refinement #2) P_TDC_nc = Pcom * (V_BDC/V_TDC)^n_comp_sim; Pnc = zeros(N,1); for i = 1:N-1 if Crank(i) <= 0 Pnc(i) = Pcom * (V_BDC/Vol(i))^n_comp_sim; else Pnc(i) = P_TDC_nc * (V_TDC/Vol(i))^n_exp_sim; end end Pnc(end) = Pnc(end-1); %% ── 7. WIEBE FUNCTION – CALIBRATED TIMING ──────────────────── % Uses SoC and bdur interpolated from LUCID MFB data (refinement #1) Wiebe_cum = zeros(N,1); MFB_rate = zeros(N,1); for i = 1:N xi = (Crank(i) - SoC_sim) / bdur_sim; if xi > 0 Wiebe_cum(i) = 1 - exp(-a_wiebe .* xi.^(m_wiebe+1)); MFB_rate(i) = (a_wiebe*(m_wiebe+1)/bdur_sim) .* xi.^m_wiebe ... .* exp(-a_wiebe .* xi.^(m_wiebe+1)); end end %% ── 8. HEAT RELEASE RATE ───────────────────────────────────── dQr = Combeff * Fuel * LHV .* MFB_rate; % J/°CA %% ── 9. COMBUSTION ODE WITH HEAT TRANSFER ───────────────────── omega = N_rpm_sim * 2*pi/60; % rad/s deg2rad = pi/180; Sp = 2 * Stroke * N_rpm_sim/60; % mean piston speed (m/s) % Woschni constants C_woschni = 3.26; C1 = 2.28; C2 = 3.24e-3; P_ref_w = Pcom; V_ref_w = Vol(1); T_ref_w = 320; % K (assumed charge temperature at IVC) Pc = zeros(N,1); Pc(1) = Pcom; Temp = zeros(N,1); Temp(1) = (Pcom*Vol(1)) / (m_total*R); Qloss = zeros(N,1); for i = 1:N-1 % Variable gamma (refinement #4) gam = gamma_cold + (gamma_hot - gamma_cold) * Wiebe_cum(i); % Heat release this step dQ_comb = dQr(i); % Woschni heat transfer (refinement #3) if Crank(i) >= SoC_sim w_vel = C1*Sp + C2*(Vd*T_ref_w/(P_ref_w*V_ref_w))*(Pc(i) - Pnc(i)); else w_vel = C1*Sp; end w_vel = max(w_vel, 0.1); h_w = C_woschni * Bore^(-0.2) * (Pc(i)/1000)^0.8 ... * Temp(i)^(-0.55) * w_vel^0.8; dQht = h_w * A_surf(i) * (Temp(i) - T_wall) / (omega/deg2rad); Qloss(i) = dQht; % First-law ODE (Euler step) dQ_net = dQ_comb - dQht; dV = dVol(i); Pc(i+1) = Pc(i) ... + ((gam-1)/Vol(i)) * dQ_net ... - gam*(Pc(i)/Vol(i)) * dV; Temp(i+1) = (Pc(i+1)*Vol(i+1)) / (m_total*R); end Pc(end) = Pc(end-1); Temp(end) = Temp(end-1); %% ── 10. PERFORMANCE METRICS ────────────────────────────────── HeatReleased = sum(dQr); HeatLoss = sum(Qloss); TotalFuelE = Fuel * LHV; % Work integrals idx_TDC = find(Crank == 0); W_comp = trapz(Vol(1:idx_TDC), Pc(1:idx_TDC)); W_exp = trapz(Vol(idx_TDC:end), Pc(idx_TDC:end)); W_net = W_exp - W_comp; % net indicated work (J) IMEP_sim = W_net / Vd; % Pa → bar GMEP_sim = abs(W_exp) / Vd; % gross positive work ITE = W_net / TotalFuelE * 100; [PeakP, idxP] = max(Pc); [PeakT, idxT] = max(Temp); [~, idxPR] = max(abs(diff(Pc))); maxPRR = abs(Pc(idxPR+1) - Pc(idxPR)); % Pa/deg %% ── 11. VALIDATION – ERROR METRICS ────────────────────────── % Compare simulated vs LUCID measured values errP_pct = (PeakP/1e5 - PMAX_ref) / PMAX_ref * 100; errIMEP_pct = (IMEP_sim/1e5 - NMEP_ref) / NMEP_ref * 100; errGMEP_pct = (GMEP_sim/1e5 - GMEP_ref) / GMEP_ref * 100; fprintf('\n=== VALIDATION vs LUCID DATA (%d rpm) ===\n', N_rpm_sim); fprintf('%-28s %10s %10s %8s\n','Metric','Simulated','Measured','Error %%'); fprintf('%s\n', repmat('-',1,60)); fprintf('%-28s %10.2f %10.2f %+7.1f%%\n', ... 'Peak Pressure (bar)', PeakP/1e5, PMAX_ref, errP_pct); fprintf('%-28s %10.2f %10.2f %+7.1f%%\n', ... 'IMEP / NMEP (bar)', IMEP_sim/1e5, NMEP_ref, errIMEP_pct); fprintf('%-28s %10.2f %10.2f %+7.1f%%\n', ... 'GMEP (bar)', GMEP_sim/1e5, GMEP_ref, errGMEP_pct); fprintf('%-28s %10.1f %10s\n', ... 'Peak Temperature (K)', PeakT, 'N/A (calc.)'); fprintf('%-28s %10.1f %10s\n', ... 'Heat Transfer Loss (J)', HeatLoss, 'N/A'); fprintf('%-28s %10.1f%%\n', ... 'Indicated Thermal Eff.', ITE); fprintf('%-28s %10.1f kPa/deg\n', ... 'Max Pressure Rise Rate', maxPRR/1e3); %% ── 12. SWEEP ALL SPEEDS FOR COMPARISON CURVES ─────────────── nSpeeds = length(RPM_exp); PMAX_sim_all = zeros(1,nSpeeds); IMEP_sim_all = zeros(1,nSpeeds); GMEP_sim_all = zeros(1,nSpeeds); PeakT_sim_all = zeros(1,nSpeeds); for k = 1:nSpeeds rpm_k = RPM_exp(k); SoC_k = MFB02_exp(k); bdur_k = BurnDur_exp(k); nc_k = NCOMP_exp(k); ne_k = NEXP_exp(k); Sp_k = 2 * Stroke * rpm_k/60; omega_k = rpm_k * 2*pi/60; P_TDC_k = Pcom * (V_BDC/V_TDC)^nc_k; Pnc_k = zeros(N,1); for i = 1:N-1 if Crank(i) <= 0 Pnc_k(i) = Pcom * (V_BDC/Vol(i))^nc_k; else Pnc_k(i) = P_TDC_k * (V_TDC/Vol(i))^ne_k; end end Pnc_k(end) = Pnc_k(end-1); % Wiebe MFB_k = zeros(N,1); WC_k = zeros(N,1); for i = 1:N xi = (Crank(i) - SoC_k) / bdur_k; if xi > 0 WC_k(i) = 1 - exp(-a_wiebe * xi^(m_wiebe+1)); MFB_k(i) = (a_wiebe*(m_wiebe+1)/bdur_k) * xi^m_wiebe ... * exp(-a_wiebe * xi^(m_wiebe+1)); end end dQr_k = Combeff * Fuel * LHV .* MFB_k; % ODE Pc_k = zeros(N,1); Pc_k(1) = Pcom; Tmp_k = zeros(N,1); Tmp_k(1) = (Pcom*Vol(1))/(m_total*R); for i = 1:N-1 gam_k = gamma_cold + (gamma_hot - gamma_cold)*WC_k(i); dQc_k = dQr_k(i); if Crank(i) >= SoC_k w_k = C1*Sp_k + C2*(Vd*T_ref_w/(P_ref_w*V_ref_w))*(Pc_k(i)-Pnc_k(i)); else w_k = C1*Sp_k; end w_k = max(w_k, 0.1); h_k = C_woschni*Bore^(-0.2)*(Pc_k(i)/1000)^0.8*Tmp_k(i)^(-0.55)*w_k^0.8; dQht_k = h_k * A_surf(i) * (Tmp_k(i)-T_wall) / (omega_k/deg2rad); dQ_net_k = dQc_k - dQht_k; Pc_k(i+1) = Pc_k(i) + ((gam_k-1)/Vol(i))*dQ_net_k - gam_k*(Pc_k(i)/Vol(i))*dVol(i); Tmp_k(i+1) = (Pc_k(i+1)*Vol(i+1)) / (m_total*R); end Pc_k(end) = Pc_k(end-1); Wn_k = trapz(Vol(idx_TDC:end),Pc_k(idx_TDC:end)) - trapz(Vol(1:idx_TDC),Pc_k(1:idx_TDC)); Wg_k = abs(trapz(Vol(idx_TDC:end),Pc_k(idx_TDC:end))); PMAX_sim_all(k) = max(Pc_k)/1e5; IMEP_sim_all(k) = Wn_k/Vd/1e5; GMEP_sim_all(k) = Wg_k/Vd/1e5; PeakT_sim_all(k) = max(Tmp_k); end %% ── 13. PLOTS ───────────────────────────────────────────────── % ── Figure 1: Simulated P vs CA, highlighting chosen speed ──── figure(1); clf; plot(Crank, Pnc/1e5, 'b--', 'LineWidth',1.5); hold on; plot(Crank, Pc/1e5, 'r-', 'LineWidth',2.0); xline(SoC_sim, 'k:', sprintf('SoC=%.1f°',SoC_sim), 'LabelHorizontalAlignment','right'); xline(SoC_sim+bdur_sim, 'g:', sprintf('EoC=%.1f°',SoC_sim+bdur_sim)); xline(Crank(idxP), 'm:', sprintf('Peak=%.1f bar',PeakP/1e5)); xlabel('Crank Angle (°ATDC)'); ylabel('Pressure (bar)'); title(sprintf('Refined Simulation: Cylinder Pressure vs CA @ %d rpm', N_rpm_sim)); legend('No Combustion (polytropic)','With Combustion (refined)','Location','northwest'); grid on; xlim([-180 180]); annotation('textbox',[0.62 0.65 0.30 0.22], ... 'String',{sprintf('n_{comp} = %.4f (LUCID)',n_comp_sim), ... sprintf('n_{exp} = %.4f (LUCID)',n_exp_sim), ... sprintf('SoC = %.1f° (MFB2%% data)',SoC_sim), ... sprintf('bdur = %.1f° (MFB10–90%%)',bdur_sim)}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.9],'FontSize',9); % ── Figure 2: P-V diagram ────────────────────────────────────── figure(2); clf; plot(Vol*1e6, Pnc/1e5, 'b--','LineWidth',1.5); hold on; plot(Vol*1e6, Pc/1e5, 'r-', 'LineWidth',2.0); [Ppkbar,ipk] = max(Pc/1e5); plot(Vol(ipk)*1e6, Ppkbar, 'rv','MarkerFaceColor','r','MarkerSize',10); text(Vol(ipk)*1e6+3, Ppkbar, sprintf('%.1f bar @ %d°',Ppkbar,Crank(ipk)),'FontSize',9,'Color','r'); xlabel('Cylinder Volume (cc)'); ylabel('Pressure (bar)'); title(sprintf('Refined P-V Diagram @ %d rpm', N_rpm_sim)); legend('No Combustion','With Combustion','Peak P'); grid on; % ── Figure 3: Simulated vs measured peak pressure (all speeds) ─ figure(3); clf; plot(RPM_exp, PMAX_exp, 'ko-','LineWidth',2.0,'MarkerSize',8,'MarkerFaceColor','k', ... 'DisplayName','LUCID Measured'); hold on; plot(RPM_exp, PMAX_sim_all, 'r^--','LineWidth',1.8,'MarkerSize',8,'MarkerFaceColor','r', ... 'DisplayName','Refined Simulation'); xline(N_rpm_sim,'b:',sprintf('%d rpm',N_rpm_sim),'LabelHorizontalAlignment','right'); xlabel('Engine Speed (rpm)'); ylabel('Peak Cylinder Pressure (bar)'); title('Peak Pressure: Simulation vs LUCID Measurement (All Speeds)'); legend('Location','southeast'); grid on; annotation('textbox',[0.55 0.65 0.37 0.20], ... 'String',{sprintf('Error at %d rpm: %+.1f%%',N_rpm_sim,errP_pct), ... 'Remaining error sources:', ... ' - ideal gas vs real gas', ... ' - Wiebe shape mismatch', ... ' - Woschni calibration'}, ... 'FitBoxToText','on','BackgroundColor',[1 1 0.9],'FontSize',8.5); % ── Figure 4: Simulated vs measured IMEP / NMEP (all speeds) ── figure(4); clf; plot(RPM_exp, NMEP_exp, 'ko-','LineWidth',2.0,'MarkerSize',8,'MarkerFaceColor','k', ... 'DisplayName','LUCID NMEP Measured'); hold on; plot(RPM_exp, IMEP_sim_all, 'r^--','LineWidth',1.8,'MarkerSize',8,'MarkerFaceColor','r', ... 'DisplayName','Simulated IMEP'); plot(RPM_exp, GMEP_exp, 'b s-','LineWidth',1.5,'MarkerSize',7,'MarkerFaceColor','b', ... 'DisplayName','LUCID GMEP Measured'); plot(RPM_exp, GMEP_sim_all, 'gs--','LineWidth',1.5,'MarkerSize',7,'MarkerFaceColor','g', ... 'DisplayName','Simulated GMEP'); xline(N_rpm_sim,'b:',sprintf('%d rpm',N_rpm_sim),'LabelHorizontalAlignment','right'); xlabel('Engine Speed (rpm)'); ylabel('Mean Effective Pressure (bar)'); title('IMEP / GMEP: Simulation vs LUCID Measurement (All Speeds)'); legend('Location','northeast'); grid on; % ── Figure 5: Residual error (% deviation from LUCID) ───────── figure(5); clf; errP_all = (PMAX_sim_all - PMAX_exp) ./ PMAX_exp * 100; errIMEP_all = (IMEP_sim_all - NMEP_exp) ./ NMEP_exp * 100; plot(RPM_exp, errP_all, 'r o-','LineWidth',1.8,'MarkerSize',7,'MarkerFaceColor','r', ... 'DisplayName','Peak Pressure Error'); hold on; plot(RPM_exp, errIMEP_all, 'b s-','LineWidth',1.8,'MarkerSize',7,'MarkerFaceColor','b', ... 'DisplayName','IMEP Error'); yline( 10,'k:','LineWidth',1.0,'Label','+10%'); yline(-10,'k:','LineWidth',1.0,'Label','-10%'); yline( 0,'k-','LineWidth',0.8); xline(N_rpm_sim,'g:',sprintf('%d rpm',N_rpm_sim),'LabelHorizontalAlignment','right'); xlabel('Engine Speed (rpm)'); ylabel('Error (%)'); title('Simulation Error vs LUCID Measurement (Simulation - Measured)/Measured'); legend('Location','northeast'); grid on; ylim([-30 30]); % ── Figure 6: Temperature profile at chosen speed ───────────── figure(6); clf; plot(Crank, Temp, 'm-','LineWidth',1.8); xline(Crank(idxT), 'k:', sprintf('Peak = %.0f K at %d°', PeakT, Crank(idxT))); xline(SoC_sim, 'g:','SoC','LabelHorizontalAlignment','right'); xlabel('Crank Angle (°ATDC)'); ylabel('Temperature (K)'); title(sprintf('Cylinder Gas Temperature @ %d rpm', N_rpm_sim)); grid on; % ── Figure 7: Heat release vs heat loss ─────────────────────── figure(7); clf; yyaxis left; plot(Crank, dQr*1000, 'r-','LineWidth',1.8); ylabel('Combustion Heat Release (mJ/°CA)'); yyaxis right; plot(Crank, Qloss*1000, 'b--','LineWidth',1.8); ylabel('Woschni Heat Loss (mJ/°CA)'); xline(SoC_sim, 'k:','SoC','LabelHorizontalAlignment','right'); xline(SoC_sim+bdur_sim, 'g:','EoC'); xlabel('Crank Angle (°ATDC)'); title(sprintf('Heat Release vs Heat Transfer Loss @ %d rpm', N_rpm_sim)); legend({'dQ_{release}','dQ_{loss (Woschni)}'},'Location','northeast'); grid on; % ── Figure 8: Pressure rise rate (knock indicator) ──────────── figure(8); clf; dPdc = [0; diff(Pc)]/1e3; % kPa/deg plot(Crank, dPdc,'k-','LineWidth',1.2); hold on; yline(500,'r--','LineWidth',1.5,'Label','Knock limit 500 kPa/°'); yline(-500,'r--','LineWidth',1.0); xline(SoC_sim,'g:','SoC','LabelHorizontalAlignment','right'); xlabel('Crank Angle (°ATDC)'); ylabel('Pressure Rise Rate (kPa/°CA)'); title(sprintf('Pressure Rise Rate (Knock Indicator) @ %d rpm', N_rpm_sim)); grid on; %% ── 14. PRINTED SUMMARY ─────────────────────────────────────── fprintf('\n=== FULL SWEEP COMPARISON TABLE ===\n'); fprintf('%-6s %-10s %-10s %-8s %-10s %-10s %-8s\n', ... 'RPM','Pmx_sim','Pmx_meas','ErrP%%','IMEP_sim','NMEP_meas','ErrI%%'); fprintf('%s\n',repmat('-',1,64)); for k = 1:nSpeeds ep = (PMAX_sim_all(k)-PMAX_exp(k))/PMAX_exp(k)*100; ei = (IMEP_sim_all(k)-NMEP_exp(k))/NMEP_exp(k)*100; fprintf('%-6d %-10.2f %-10.2f %-8.1f %-10.2f %-10.2f %-8.1f\n', ... RPM_exp(k), PMAX_sim_all(k), PMAX_exp(k), ep, ... IMEP_sim_all(k), NMEP_exp(k), ei); end fprintf('%s\n',repmat('-',1,64)); % Grader variables HeatReleased = sum(dQr); PeakP = max(Pc); PeakT = max(Temp); %1. Calibrated Combustion Timing from MFB Data %Original problem: bdur = 30 and SoC = -7° were fixed constants regardless of engine speed. %Fix: The experimental script reads MFB02_AV, MFB10_AV, MFB90_AV per speed from LUCID. The refined model uses these as lookup tables: %This is the biggest single accuracy improvement — at 6500 rpm burn duration is ~46° vs the fixed 30°. %Original problem: Fixed n_comp = 1.33, n_exp = 1.27 regardless of speed. %Fix: LUCID measures NCOMP_AV and NEXP_AV per speed. These are read from cells M193 and M200 and interpolated — at low RPM n_comp drops toward 1.28 (more heat loss time), at high RPM it rises toward 1.33 (more adiabatic). This directly affects both the no-combustion reference pressure and the initial conditions for the ODE. %3. Woschni Heat Transfer with Speed Scaling %The Woschni velocity term now scales with actual mean piston speed (Sp = 2 × Stroke × RPM/60), so heat losses correctly increase with engine speed rather than being constant. %4. Error Metrics Against LUCID %The script now computes and prints: %errP_pct — peak pressure % deviation from LUCID M11 %errIMEP_pct — IMEP deviation from LUCID M46 (NMEP) %A full sweep table showing error at every speed point %5. Full-Sweep Comparison Plots %Figure 3 overlays simulated vs measured peak pressure across all 14 speed points. Figure 5 shows the residual error curve — this is how you identify where the model still needs calibration (typically low RPM due to increased heat loss, and very high RPM due to charge dynamics not captured in the single-zone model). %To use with your actual data: Replace the MFB02_exp, NCOMP_exp, PMAX_exp etc. arrays with the exact values read from your LUCID spreadsheet columns M11, M74, M193 etc., then run the sweep. The error plots will show remaining discrepancies to investigate.