Nonlinear time series analysis was applied for the first time to time series obtained from large eddy simulations (LES) of an internal combustion engine (ICE). The aim of the study was to obtain more information about the cycle-to-cycle variation (CCV) in the studied simplified ICE geometry than what is available from standard methods. Phase space reconstructions were created from the time series and then estimates for the largest Lyapunov exponent were calculated. The time delays used in the phase space reconstructions were determined using average mutual information while the proper embedding dimensions were chosen according to the method of false nearest neighbours. Quantitative information on the behaviour of the flow and the CCV was acquired from three-dimensional phase space reconstructions. Introduced modifications to the flow were clearly visible in the phase space reconstructions of energy and dissipation, indicating that these quantities are appropriate for monitoring and analysing the state of the system. The estimates for the largest Lyapunov exponents were positive for almost all time series, indicating chaotic dynamics. The permutation spectrum test was used to confirm the chaoticity of the CCV. The present results indicate that the used methods offer a promising new framework for characterising the CCV from the viewpoint of nonlinear time series analysis.