%% TASK 3 - MODEL REFINEMENT AND VALIDATION % Ford Sigma 1.6 L spark-ignition engine % %% 1. USER SETTINGS filename = 'Sigma_1.6L.xls'; sheetName = '6500'; pressureColumn = 3; % 3 = MEAN(B); 4:7 = individual cylinders 1:4 exportFigures = true; % Model constants retained from Task 1 Bore = 0.079; % m Stroke = 0.0814; % m Rod = 0.137; % m Clearance = 39.27e-6; % m^3 LHV = 44e6; % J/kg Rgas = 287; % J/(kg K) rhoAir = 1.167; % kg/m^3 volEff = 0.90; combEff = 0.90; AFR = 14.7; gamma = 1.25; % Original Task 1 combustion settings (baseline comparison) baseIgnition = -12; % deg ATDC (negative = BTDC) baseIgnDelay = 5; % deg baseSOC = baseIgnition + baseIgnDelay; baseDuration = 30; % deg CA baseA = 5; baseM = 2; baseNcomp = 1.33; baseNexp = 1.27; %% 2. IMPORT MEASURED DATA % Columns: A crank angle, B volume, C mean pressure, D:G cylinders 1:4. raw = readmatrix(filename, 'Sheet', sheetName, 'Range', 'A13:G732'); raw = raw(all(isfinite(raw(:,1:3)),2),:); CA720 = raw(:,1); % crank angle, deg VfileCC = raw(:,2); % supplied volume, cm^3 Pmeas720 = raw(:,pressureColumn); % pressure, bar % Select one closed cycle centred on firing TDC. closed = CA720 >= -180 & CA720 <= 180; CA = CA720(closed); Pmeas = Pmeas720(closed); Vmeas = VfileCC(closed)*1e-6; % m^3 % Ensure column vectors and increasing crank angle. [CA, order] = sort(CA(:)); Pmeas = Pmeas(order); Vmeas = Vmeas(order); % Remove duplicate crank-angle samples, if present. [CA, uniqueIndex] = unique(CA, 'stable'); Pmeas = Pmeas(uniqueIndex); Vmeas = Vmeas(uniqueIndex); % Model grid is one degree from -180 to +180. CAmodel = (-180:1:180)'; Pmeas = interp1(CA, Pmeas, CAmodel, 'linear', 'extrap'); Vmeas = interp1(CA, Vmeas, CAmodel, 'linear', 'extrap'); CA = CAmodel; % Light smoothing is used only for derivative/timing calculations. Psm = movmean(Pmeas, 5); %% 3. ENGINE GEOMETRY A = pi*Bore^2/4; Vd = A*Stroke; Vc = Clearance; CR = (Vd + Vc)/Vc; r = Stroke/2; theta = deg2rad(CA); pistonDisplacement = r*(1-cos(theta)) + ... (Rod - sqrt(Rod^2 - (r*sin(theta)).^2)); V = Vc + A*pistonDisplacement; % The geometry calculated here should closely match the workbook volume. volumeDifference = max(abs(V - Vmeas))*1e6; % cm^3 %% 4. MEASURED POLYTROPIC INDICES % Windows avoid valve events and the main combustion period. compWindow = CA >= -140 & CA <= -35 & Psm > 0; expWindow = CA >= 50 & CA <= 140 & Psm > 0; fitComp = polyfit(log(V(compWindow)), log(Psm(compWindow)), 1); fitExp = polyfit(log(V(expWindow)), log(Psm(expWindow)), 1); nCompMeasured = -fitComp(1); nExpMeasured = -fitExp(1); % Apply physically reasonable fallback values if the fit is corrupted. if ~isfinite(nCompMeasured) || nCompMeasured < 1.05 || nCompMeasured > 1.60 warning('Compression index fit was %.3f; using Task 1 value %.2f.', ... nCompMeasured, baseNcomp); nCompMeasured = baseNcomp; end if ~isfinite(nExpMeasured) || nExpMeasured < 1.05 || nExpMeasured > 1.60 warning('Expansion index fit was %.3f; using Task 1 value %.2f.', ... nExpMeasured, baseNexp); nExpMeasured = baseNexp; end %% 5. COMBUSTION TIMING FROM APPARENT HEAT RELEASE % Closed-system first-law estimate (net apparent heat release): % dQ/dtheta = gamma/(gamma-1)*P*dV/dtheta % + 1/(gamma-1)*V*dP/dtheta Ppa = Psm*1e5; dVdCA = gradient(V, CA); dPdCA = gradient(Ppa, CA); dQnet = gamma/(gamma-1).*Ppa.*dVdCA + ... 1/(gamma-1).*V.*dPdCA; % J/deg CA % Use only the normal SI combustion window and suppress negative increments % when constructing cumulative burn fraction. combWindow = CA >= -40 & CA <= 80; dQpositive = zeros(size(dQnet)); dQpositive(combWindow) = max(dQnet(combWindow), 0); Qcum = cumtrapz(CA, dQpositive); if max(Qcum) <= 0 error('No positive apparent heat release was detected. Check the selected sheet and pressure column.'); end xbMeasured = Qcum/max(Qcum); CA10 = crossingAngle(CA, xbMeasured, 0.10); CA50 = crossingAngle(CA, xbMeasured, 0.50); CA90 = crossingAngle(CA, xbMeasured, 0.90); % Initial values for calibration. CA10 is used as a robust measured timing % indicator; the optimizer subsequently adjusts SOC and duration. SOC0 = CA10 - 4; duration0 = max(15, min(70, (CA90 - CA10)/0.75)); %% 6. TASK 1 BASELINE MODEL Pinitial = Pmeas(1); % anchors simulation to measured pressure at -180 deg Qfuel = rhoAir*Vd*volEff/AFR * LHV * combEff; Pbaseline = simulatePressure(CA, V, Pinitial, baseNcomp, baseNexp, ... gamma, Qfuel, baseSOC, baseDuration, baseA, baseM, 1.0); %% 7. REFINED MODEL CALIBRATION % Parameters: [SOC, duration, Wiebe m, heat-release scale]. % fminsearch is used so no Optimization Toolbox is required. p0 = [SOC0, duration0, 2.0, 1.0]; fitWindow = CA >= -60 & CA <= 140; objective = @(p) modelObjective(p, CA, V, Pmeas, Pinitial, ... nCompMeasured, nExpMeasured, gamma, Qfuel, fitWindow, baseA); options = optimset('Display','final', 'MaxIter',500, 'MaxFunEvals',1500, ... 'TolX',1e-5, 'TolFun',1e-5); pBest = fminsearch(objective, p0, options); SOCrefined = pBest(1); durationRefined = pBest(2); mRefined = pBest(3); heatScale = pBest(4); Prefined = simulatePressure(CA, V, Pinitial, nCompMeasured, nExpMeasured, ... gamma, Qfuel, SOCrefined, durationRefined, baseA, mRefined, heatScale); xbRefined = wiebeFraction(CA, SOCrefined, durationRefined, baseA, mRefined); dxbRefined = [diff(xbRefined); 0]; dQrefined = Qfuel*heatScale*dxbRefined; %% 8. VALIDATION METRICS baselineMetrics = calculateMetrics(CA, V, Pmeas, Pbaseline, Vd, fitWindow); refinedMetrics = calculateMetrics(CA, V, Pmeas, Prefined, Vd, fitWindow); improvementRMSE = 100*(baselineMetrics.RMSE_bar-refinedMetrics.RMSE_bar) ... /baselineMetrics.RMSE_bar; measuredPeak = max(Pmeas); measuredPeakAngle = CA(Pmeas == measuredPeak); measuredPeakAngle = measuredPeakAngle(1); measuredWork = trapz(V, Pmeas*1e5); measuredIMEP = measuredWork/Vd/1e5; % Include measured values and explicit model errors in the exported table. results = table(... [NaN; baselineMetrics.RMSE_bar; refinedMetrics.RMSE_bar], ... [NaN; baselineMetrics.NRMSE_percent; refinedMetrics.NRMSE_percent], ... [NaN; baselineMetrics.R2; refinedMetrics.R2], ... [measuredPeak; baselineMetrics.PeakPressure_bar; refinedMetrics.PeakPressure_bar], ... [measuredPeakAngle; baselineMetrics.PeakAngle_deg; refinedMetrics.PeakAngle_deg], ... [measuredIMEP; baselineMetrics.IMEP_bar; refinedMetrics.IMEP_bar], ... [0; baselineMetrics.PeakPressure_bar-measuredPeak; refinedMetrics.PeakPressure_bar-measuredPeak], ... [0; baselineMetrics.PeakAngle_deg-measuredPeakAngle; refinedMetrics.PeakAngle_deg-measuredPeakAngle], ... [0; baselineMetrics.IMEP_bar-measuredIMEP; refinedMetrics.IMEP_bar-measuredIMEP], ... 'VariableNames', {'RMSE_bar','NRMSE_percent','R2','PeakPressure_bar', ... 'PeakAngle_deg','IMEP_bar','PeakPressureError_bar', ... 'PeakAngleError_deg','IMEPError_bar'}, ... 'RowNames', {'Measured','Task1Baseline','RefinedModel'}); fprintf('\n============================================================\n'); fprintf('TASK 3: MODEL REFINEMENT AND VALIDATION - %s rpm\n', sheetName); fprintf('============================================================\n'); fprintf('Compression ratio : %.3f\n', CR); fprintf('Max geometry/workbook difference : %.4f cm^3\n', volumeDifference); fprintf('Measured compression index : %.3f\n', nCompMeasured); fprintf('Measured expansion index : %.3f\n', nExpMeasured); fprintf('Measured CA10 / CA50 / CA90 : %.1f / %.1f / %.1f deg ATDC\n', ... CA10, CA50, CA90); fprintf('Refined Wiebe SOC : %.2f deg ATDC\n', SOCrefined); fprintf('Refined burn duration : %.2f deg CA\n', durationRefined); fprintf('Refined Wiebe shape factor m : %.3f\n', mRefined); fprintf('Refined heat-release scale : %.3f\n', heatScale); fprintf('Measured peak pressure : %.2f bar at %.1f deg ATDC\n', ... measuredPeak, measuredPeakAngle); fprintf('Measured indicated work / IMEP : %.2f J / %.2f bar\n', ... measuredWork, measuredIMEP); fprintf('RMSE improvement vs Task 1 : %.1f %%\n\n', improvementRMSE); disp(results) writetable(results, 'Task3_validation_metrics.csv', 'WriteRowNames', true); %% 9. FIGURES FOR THE REPORT % Figure 1: measured data against baseline and refined models figure('Name','Pressure validation','Color','w'); plot(CA, Pmeas, 'k', 'LineWidth',1.8); hold on; plot(CA, Pbaseline, '--', 'LineWidth',1.3); plot(CA, Prefined, 'LineWidth',1.6); xline(0, ':', 'TDC'); xlabel('Crank angle (deg ATDC)'); ylabel('Cylinder pressure (bar)'); title(sprintf('Measured and simulated cylinder pressure at %s rpm',sheetName)); legend('Measured','Task 1 baseline','Refined model','Location','best'); grid on; xlim([-180 180]); if exportFigures, exportgraphics(gcf,'Task3_pressure_validation.png','Resolution',300); end % Figure 2: focused combustion comparison figure('Name','Combustion validation','Color','w'); plot(CA, Pmeas, 'k', 'LineWidth',1.8); hold on; plot(CA, Pbaseline, '--', 'LineWidth',1.3); plot(CA, Prefined, 'LineWidth',1.6); xline(0, ':', 'TDC'); xline(CA10, ':', 'CA10'); xline(CA50, ':', 'CA50'); xline(CA90, ':', 'CA90'); xlabel('Crank angle (deg ATDC)'); ylabel('Cylinder pressure (bar)'); title('Combustion-region pressure comparison'); legend('Measured','Task 1 baseline','Refined model','Location','best'); grid on; xlim([-50 100]); if exportFigures, exportgraphics(gcf,'Task3_combustion_pressure.png','Resolution',300); end % Figure 3: measured apparent heat release and refined Wiebe behaviour figure('Name','Combustion behaviour','Color','w'); yyaxis left plot(CA, dQnet, 'LineWidth',1.2); hold on; plot(CA, dQrefined, 'LineWidth',1.4); ylabel('Heat release rate (J/deg CA)'); yyaxis right plot(CA, xbMeasured, '--', 'LineWidth',1.4); plot(CA, xbRefined, 'LineWidth',1.5); ylabel('Cumulative burned fraction'); xlabel('Crank angle (deg ATDC)'); title('Measured apparent heat release and refined Wiebe model'); legend('Measured apparent heat release','Refined model heat release', ... 'Measured cumulative fraction','Refined Wiebe fraction','Location','best'); grid on; xlim([-50 100]); if exportFigures, exportgraphics(gcf,'Task3_heat_release.png','Resolution',300); end % Figure 4: pressure-volume validation figure('Name','PV validation','Color','w'); plot(V*1e6, Pmeas, 'k', 'LineWidth',1.8); hold on; plot(V*1e6, Pbaseline, '--', 'LineWidth',1.3); plot(V*1e6, Prefined, 'LineWidth',1.6); xlabel('Cylinder volume (cm^3)'); ylabel('Cylinder pressure (bar)'); title('Measured and simulated pressure-volume diagrams'); legend('Measured','Task 1 baseline','Refined model','Location','best'); grid on; if exportFigures, exportgraphics(gcf,'Task3_PV_validation.png','Resolution',300); end % Figure 5: residual error figure('Name','Pressure residuals','Color','w'); plot(CA, Pbaseline-Pmeas, '--', 'LineWidth',1.2); hold on; plot(CA, Prefined-Pmeas, 'LineWidth',1.4); yline(0, ':'); xlabel('Crank angle (deg ATDC)'); ylabel('Model - measured pressure (bar)'); title('Cylinder-pressure residuals'); legend('Task 1 baseline error','Refined-model error','Location','best'); grid on; xlim([-180 180]); if exportFigures, exportgraphics(gcf,'Task3_pressure_residuals.png','Resolution',300); end %% LOCAL FUNCTIONS function xb = wiebeFraction(CA, SOC, duration, a, m) duration = max(duration, 1e-6); z = (CA-SOC)/duration; xb = zeros(size(CA)); active = z >= 0 & z <= 1; xb(active) = (1-exp(-a*z(active).^(m+1))) / (1-exp(-a)); xb(z > 1) = 1; end function P = simulatePressure(CA, V, P0, nComp, nExp, gamma, Qfuel, ... SOC, duration, a, m, heatScale) xb = wiebeFraction(CA, SOC, duration, a, m); dxb = diff(xb); P = zeros(size(CA)); P(1) = P0; for k = 1:numel(CA)-1 if CA(k) <= 0 n = nComp; else n = nExp; end dQ = Qfuel*heatScale*max(dxb(k),0); P(k+1) = P(k)*(V(k)/V(k+1))^n + ... (gamma-1)*dQ/(V(k+1)*1e5); end end function J = modelObjective(p, CA, V, Pmeas, P0, nComp, nExp, ... gamma, Qfuel, fitWindow, a) SOC = p(1); duration = p(2); m = p(3); heatScale = p(4); % Soft constraints keep the search physically meaningful. penalty = 0; penalty = penalty + 1e4*max(0,-40-SOC)^2 + 1e4*max(0,SOC-20)^2; penalty = penalty + 1e3*max(0,10-duration)^2 + 1e3*max(0,duration-90)^2; penalty = penalty + 1e3*max(0,0.2-m)^2 + 1e3*max(0,m-6)^2; penalty = penalty + 1e3*max(0,0.2-heatScale)^2 + 1e3*max(0,heatScale-2.5)^2; Pmodel = simulatePressure(CA,V,P0,nComp,nExp,gamma,Qfuel, ... SOC,duration,a,m,heatScale); residual = Pmodel(fitWindow)-Pmeas(fitWindow); J = mean(residual.^2) + penalty; end function metrics = calculateMetrics(CA, V, Pmeas, Pmodel, Vd, evalWindow) e = Pmodel(evalWindow)-Pmeas(evalWindow); y = Pmeas(evalWindow); metrics.RMSE_bar = sqrt(mean(e.^2)); metrics.NRMSE_percent = 100*metrics.RMSE_bar/(max(y)-min(y)); metrics.R2 = 1-sum(e.^2)/sum((y-mean(y)).^2); [metrics.PeakPressure_bar, idx] = max(Pmodel); metrics.PeakAngle_deg = CA(idx); work = trapz(V, Pmodel*1e5); metrics.IMEP_bar = work/Vd/1e5; end function angle = crossingAngle(CA, y, target) idx = find(y >= target, 1, 'first'); if isempty(idx) angle = NaN; elseif idx == 1 angle = CA(1); else angle = interp1(y(idx-1:idx), CA(idx-1:idx), target, 'linear'); end end