How to Perform a Thermodynamic Analysis of an Otto Cycle Engine Using Matlab

Understanding the Air‑Standard Otto Cycle Model

The Otto cycle is the thermodynamic idealization of the four‑stroke spark‑ignition engine. It replaces the complex intake, compression, combustion, and exhaust processes with four reversible steps that isolate the essential energy conversion mechanisms. In the air‑standard model, the working fluid is air (treated as an ideal gas with constant specific heats), combustion is approximated as instantaneous heat addition at constant volume, and exhaust is replaced by constant‑volume heat rejection. This simplification allows direct closed‑form calculation of temperatures, pressures, work output, and thermal efficiency, making it an indispensable first step for understanding engine performance scaling.

The four processes that form the cycle are:

Isentropic compression (1→2): The piston compresses the charge without heat transfer, raising temperature and pressure.
Constant‑volume heat addition (2→3): An instantaneous release of energy (representing combustion) at top dead centre raises the temperature and pressure to their peak values.
Isentropic expansion (3→4): The high‑pressure gases push the piston down, delivering net work to the crankshaft.
Constant‑volume heat rejection (4→1): The cylinder contents lose heat back to the initial state, closing the loop.

Because the model uses a constant specific heat ratio and ignores gas exchange, it provides a quick benchmark that can be systematically refined. Implementing this cycle in MATLAB gives you a flexible, interactive tool for exploring how compression ratio, heat addition, and fluid properties affect performance. The air‑standard assumption also provides a convenient upper bound on efficiency; real engines always fall short due to friction, heat losses, and incomplete combustion.

Setting Up a Structured MATLAB Environment

Begin by creating a new script or function file. A clean workspace makes debugging and future modifications straightforward. Use a descriptive header that documents the purpose, input parameters, and outputs. Pre‑allocate arrays for state variables when sweeping over parameters to avoid dynamic resizing overhead. Adopt meaningful variable names such as compression_ratio rather than just r (though r is standard in thermodynamics texts) and keep units consistent throughout (SI units recommended). A typical header block might read:

% otto_cycle_analysis.m
% Air‑standard Otto cycle analysis
% Inputs:
%   r        – compression ratio (dimensionless)
%   T1, P1   – intake state (K, Pa)
%   gamma    – specific heat ratio (cp/cv)
%   Q_in     – specific heat addition (J/kg)
% Outputs:
%   state_table – array of T, P, v at four corners
%   W_net, eta  – specific net work (J/kg) and thermal efficiency

Using MATLAB’s built‑in table or struct arrays to store the four states improves readability when you later plot or export results. For instance, create a structure array state with fields T, P, v indexed from 1 to 4. This makes referencing state properties intuitive: state(3).T gives the temperature after combustion.

Defining the Critical Engine Parameters

The accuracy of the simulation depends on the input data. At minimum, specify the compression ratio, intake pressure and temperature, the specific heat ratio, and the energy added per unit mass. Typical values for a gasoline engine are:

Compression ratio r: 8–12 for modern spark‑ignition engines; higher ratios improve efficiency but increase knock risk. In research engines, ratios as high as 14 are sometimes explored with advanced fuels.
Intake temperature T₁: 300 K (ambient) or 350 K if accounting for intake air heating from the manifold and residual gases.
Intake pressure P₁: 101325 Pa (1 atm) for naturally aspirated; up to 200 kPa or more for turbocharged engines.
Specific heat ratio γ: 1.4 for dry air; 1.3–1.35 for fuel‑air mixtures and combustion products due to the presence of tri‑atomic molecules like CO₂ and H₂O.
Specific heat addition Q_in: approximately 2.5–3 MJ/kg for gasoline based on lower heating value and stoichiometric air‑fuel ratio. For a typical gasoline, the lower heating value is about 44 MJ/kg, but mixed with air (roughly 14.7:1 air‑fuel ratio), the mixture energy density becomes about 2.9 MJ/kg of mixture.

Define these as MATLAB variables at the top of the script so they can be easily changed for sensitivity studies.

r      = 9.0;          % compression ratio
T1     = 300;          % K
P1     = 101325;       % Pa
gamma  = 1.35;         % approximate for fuel‑air mixture
Q_in   = 2500e3;       % J/kg
R_air  = 287;          % J/(kg·K)

You may also include the molar mass of air (28.97 g/mol) if you later switch to mole‑based calculations. Keeping the gas constant explicit avoids unit errors.

Computing the Four Corner States

With the parameters set, the air‑standard relations give the state points directly. The analysis is per unit mass, so we work with specific volume. Setting the clearance volume to unity (v₂ = 1) simplifies the relative scaling; then v₁ = r because the compression ratio v₁/v₂ = r. This normalization is convenient because all volumes are expressed relative to the clearance volume.

State 1 (Intake, Bottom Dead Centre)

v1 = r;                % relative to clearance volume
% T1 and P1 are inputs

State 2 (Isentropic Compression to Top Dead Centre)

For an isentropic process with constant γ:

v2 = 1;
T2 = T1 * r^(gamma-1);
P2 = P1 * r^gamma;

These relations derive from the ideal gas law and the polytropic relation P v^γ = constant. The temperature rise across compression is significant: for r=9 and γ=1.35, T₂ ≈ 300 × 9^(0.35) ≈ 690 K.

State 3 (Constant‑Volume Heat Addition)

First compute the specific heat at constant volume:

cv = R_air / (gamma-1);
T3 = T2 + Q_in / cv;
P3 = P2 * (T3 / T2);
v3 = v2;   % no volume change

The temperature at state 3 can exceed 2500 K for typical gasoline combustion. Peak pressure P₃ often reaches 60–80 bar, which sets structural requirements for the engine block and head.

State 4 (Isentropic Expansion to Original Volume)

The expansion ratio from TDC to BDC is r, but in terms of specific volumes v₃ = v₂ = 1 and v₄ = v₁ = r, so the volume ratio is 1/r.

v4 = v1;
T4 = T3 * (1/r)^(gamma-1);
P4 = P3 * (1/r)^gamma;

These four points fully define the cycle in terms of temperature, pressure, and specific volume. Store them in arrays for later plotting and calculation. A useful validation is to check that P₄ is slightly above P₁; the difference drives the exhaust blowdown.

Calculating Net Work, Thermal Efficiency, and Mean Effective Pressure

The net specific work is the difference between heat added and heat rejected. For the air‑standard model, the heat rejected in process 4→1 is:

Q_out = cv * (T4 - T1);
W_net = Q_in - Q_out;

The thermal efficiency is:

eta  = W_net / Q_in;            % from energy balance
eta_direct = 1 - 1/r^(gamma-1); % direct formula
% Both should match within numerical round-off

Mean effective pressure (MEP) is a performance metric that compares net work per cycle to the displacement volume. For the air‑standard cycle, the specific displacement volume is v₁ – v₂ = r – 1. Therefore:

mep = W_net / (v1 - v2);        % Pa

MEP is independent of engine size and is a convenient indicator of engine design quality. Typical values for naturally aspirated gasoline engines are 8–12 bar (0.8–1.2 MPa). Output MEP alongside efficiency to provide a more complete picture. You may also convert to bar for readability: mep_bar = mep / 1e5.

Including Variable Specific Heats for Improved Accuracy

The constant‑γ assumption can overestimate efficiency because real exhaust gases have lower specific heat ratios at high temperatures. To account for this, replace the constant γ with a temperature‑dependent correlation. A robust method is to use the NASA polynomial coefficients for c_p of major species (O₂, N₂, CO₂, H₂O). The NIST Chemistry WebBook provides these coefficients for pure substances. For the air‑standard model, you can approximate the working fluid as either dry air or a “combustion product” mixture with a fixed composition.

Implementing variable γ requires an iterative solution because the isentropic exponent itself depends on the unknown temperature. Use MATLAB’s fzero or a simple loop. For example, during isentropic compression, the relation becomes:

% Use an anonymous function to find T2 such that
% integral of cp dT from T1 to T2 = R * ln(v1/v2)
% This replaces the constant‑gamma formula.
% For air, cp(T) can be approximated by:
% cp = a0 + a1*T + a2*T^2 + a3*T^3 + a4*T^4
% Using NASA 7-coefficient polynomials.
fun = @(T2) integral(@(T) polyval(cp_coeffs, T), T1, T2) - R_air * log(r);
T2 = fzero(fun, T1 * r^0.4);

The NASA polynomial resource offers downloadable coefficients for air and combustion species. Modify your script to load these coefficients and call them during state calculations. You will observe that variable specific heats reduce the predicted efficiency by 2–5 percentage points at typical compression ratios, which aligns closely with real engine data.

For a quick approximation without full polynomial integration, use a mean specific heat ratio that varies linearly with temperature: γ(T) = 1.4 – 0.00005*(T – 300). This simple correlation captures the trend and can be implemented in a while loop until convergence.

Building a Modular Function‑Driven Analysis

Rather than a monolithic script, encapsulate the cycle logic in a MATLAB function. This makes it reusable for parameter sweeps and reduces code duplication. A function skeleton is shown below:

function [states, performance] = ottoCycle(r, T1, P1, Q_in, gamma, R_air)
    % states: structure with T, P, v for states 1-4
    % performance: structure with W_net, eta, mep
    % (calculation code from previous section goes here)
end

For variable specific heats, add a fluid flag:

function [states, performance] = ottoCycle(r, T1, P1, Q_in, fluid_type)
    % fluid_type: 'constant_gamma' or 'variable_cp'
    % (logic branches based on flag)
end

With the function in place, a driver script can sweep through compression ratios:

r_values = 6:1:14;
for i = 1:length(r_values)
    [~, perf] = ottoCycle(r_values(i), 300, 101325, 2500e3, 'constant_gamma');
    eta(i) = perf.eta;
    mep(i) = perf.mep;
end
plot(r_values, eta*100, 'b-', r_values, mep/1e5, 'r--');

This modular approach immediately shows the diminishing returns of increasing compression ratio beyond 12. The plot highlights a key trade‑off: efficiency continues to rise but at a decreasing rate, while MEP often peaks around r=10 for typical heat inputs.

Generating Diagnostic P‑V and T‑S Diagrams

Visualisation helps verify the model and communicate results. MATLAB’s plotting functions produce publication‑quality diagrams.

Pressure‑Volume (P‑V) Diagram

The compression and expansion curves follow P v^γ = constant. Generate a series of intermediate points for smooth curves:

n = 50;
v_comp = linspace(v2, v1, n);   % reverse direction for plotting
P_comp = P2 * (v_comp ./ v2).^(-gamma);
v_exp  = linspace(v3, v4, n);
P_exp  = P3 * (v_exp ./ v3).^(-gamma);
figure;
hold on;
plot(v_comp, P_comp, 'b', 'LineWidth', 1.5);
plot(v_exp, P_exp, 'b', 'LineWidth', 1.5);
plot([v2 v3], [P2 P3], 'r--', 'LineWidth', 1.2);
plot([v4 v1], [P4 P1], 'r--', 'LineWidth', 1.2);
xlabel('Specific volume (relative)'); ylabel('Pressure (Pa)');
title('Air‑Standard Otto Cycle P–V Diagram');
legend('Compression', 'Expansion', 'Heat addition', 'Heat rejection');
grid on;

The enclosed area represents the net work. You can numerically integrate the loop using trapz to verify the computed W_net. For example, integrate the upper curve (compression + expansion) and subtract the lower curve (heat rejection line). Refer to the MATLAB documentation on numerical integration for details.

Temperature‑Entropy (T‑S) Diagram

For an ideal gas with constant specific heats, the entropy change relative to state 1 is:

% s = cp * log(T / T1) - R * log(P / P1)
cp = gamma * R_air / (gamma - 1);
s1 = 0;
s2 = cp * log(T2/T1) - R_air * log(P2/P1);
s3 = cp * log(T3/T1) - R_air * log(P3/P1);
s4 = cp * log(T4/T1) - R_air * log(P4/P1);

Plot the four processes with appropriate line styles. For variable γ, use numerical integration of dQ/T to compute entropy. The T‑S diagram reveals the isentropic nature of processes 1‑2 and 3‑4 (vertical lines) and the constant‑volume heating/rejection (curved paths). The area under the heating curve represents the heat added, while the area under the cooling curve represents the heat rejected; their difference is the net work.

To improve the T‑S plot, overlay lines of constant pressure (isobars) from the ideal gas relation. This helps contextualise the cycle in a broader thermodynamic framework.

Performing Sensitivity Studies

Once the baseline model works, systematically vary the input parameters to understand their influence. Common studies include:

Compression ratio sweep: Plot efficiency and MEP versus r. Efficiency increases rapidly at low r and then plateaus; MEP often peaks at a moderate r because of the trade‑off between pressure rise and heat addition. Use subplots to show both metrics on the same figure.
Heat addition variation: Changing Q_in models different fuel energies or air‑fuel ratios. Increasing Q_in raises peak pressure and work but may exceed material limits; include a constraint that P₃ < 150 bar (15 MPa) to reflect typical engine strength.
Specific heat ratio effect: Lower γ reduces efficiency. Overlay constant‑γ and variable‑γ results to show the correction. This is especially important when analysing alternative fuels like natural gas or hydrogen.

For a two‑parameter sweep, use nested loops and store results in a matrix. Then create a contour plot or colour map to visualise the combined effect:

r_vec = 6:1:14;
Q_vec = 2e6:0.2e6:3e6;
for i = 1:length(r_vec)
    for j = 1:length(Q_vec)
        [~, perf] = ottoCycle(r_vec(i), 300, 101325, Q_vec(j), 1.35, 287);
        eta_mat(i,j) = perf.eta;
    end
end
contourf(r_vec, Q_vec/1e6, eta_mat', 20);
xlabel('Compression ratio'); ylabel('Q_in (MJ/kg)'); 
colorbar; title('Thermal Efficiency (%)');

Such sensitivity maps are invaluable for trade‑off analysis in engine design. They help identify operating regions where efficiency is high while peak pressure remains acceptable.

Validating Your Results

Before using the model for engineering decisions, check against known benchmarks. The air‑standard efficiency formula η = 1 − 1/r^(γ−1) provides a direct validation for the constant‑γ case. Additionally, compare calculated peak pressure P₃ against data from published engine maps or textbook examples (e.g., Stone, Introduction to Internal Combustion Engines). Discrepancies often arise from the assumed γ or Q_in. For a thorough validation, download a public cylinder‑pressure trace from resources like the Sandia National Laboratories engine database and overlay your P–V loop. This exercise reveals where the air‑standard model diverges from reality, motivating more complex modelling.

Another validation technique is to check the energy balance: the net work plus rejected heat should equal the input heat within floating‑point precision. Also verify that the compression and expansion processes satisfy the isentropic relation T v^(γ‑1) = constant to detect coding errors.

Extending the Model for Realistic Effects

The basic framework can be expanded to include real‑world phenomena:

Finite combustion duration: Replace the instantaneous heat addition with a Wiebe function that distributes energy release over a crank‑angle interval. Use MATLAB’s ode45 to integrate the energy equation together with the piston motion. The Wiebe function is defined as x_b(θ) = 1 – exp(–a * ((θ‑θ₀)/Δθ)^(m+1)).
Heat transfer losses: Incorporate a Woschni‑type correlation for convective heat transfer to the cylinder walls. This requires crank‑angle resolved temperature and pressure. The correlation gives the heat transfer coefficient as a function of bore, piston speed, temperature, and pressure.
Exhaust gas recirculation (EGR): Adjust the specific heat ratio and gas composition to reflect dilution. The MATLAB Simscape Engine Cylinder block provides a ready‑made component for such studies, but you can also implement a simplified dilution model by linearly combining the specific heats of air and inert gas.

For a quick extension, start with the Wiebe function. Example approximations can be found in the MATLAB File Exchange. This simple addition captures the effect of spark timing on efficiency and peak pressure.

Integrating the Analysis into a Larger Workflow

In an engineering project, the cycle analysis must feed into vehicle‑level simulations or report generation. Use MATLAB’s writetable to export state tables and performance metrics to CSV for further processing in Excel. Combine plots into a PDF using the exportgraphics function. Structuring the code with functions and a main driver script from the start will make this integration seamless. For example:

results = table(r_values', eta', mep', ...
    'VariableNames', {'CompressionRatio','Efficiency','MEP_bar'});
writetable(results, 'otto_sweep.csv');

You can also call MATLAB from Python or Simulink to couple the thermodynamic model with other subsystems. For a hybrid vehicle study, the Otto cycle efficiency map can be imported as a lookup table in Simulink.

Conclusion

Building an air‑standard Otto cycle model in MATLAB provides a powerful yet simple platform for exploring spark‑ignition engine thermodynamics. By systematically defining parameters, computing state points, and calculating performance metrics such as net work, thermal efficiency, and mean effective pressure, you gain immediate insight into how design choices affect output. The modular approach presented here allows you to extend the model to variable specific heats, finite combustion, and heat transfer without rewriting core logic. The resulting tool is not just for academic exercises—it delivers actionable answers to the “what‑if” questions that drive engine development. With the next steps of sensitivity studies and validation against real data, your MATLAB script becomes a reliable foundation for more advanced engine simulation.