In this work cyclic variability is examined in a methane-fueled spark-ignition (SI) engine. A new numerical approach is developed and integrated into the ignition model of an in-house computational fluid dynamics (CFD) code, with the aim to estimate the effect of this highly-complex phenomenon in engines. This sub-model focuses on small-scale turbulence, with its main feature concerning the prediction of the flame propagation once the relative position of the initial flamelet and the local turbulent eddy is defined. This is accomplished by using a blending of the laminar and turbulent flame speeds at the spark cell, according to the flame surface that lies outside the local eddy, smoothing the transition from the laminar to turbulent combustion regimes. A methodology for estimating the coefficient of variation (COV) of indicated mean effective pressure (IMEP) is then proposed based on the simulation of a large number of closed engine cycles. A validated CFD code including this sub-model is applied for the cyclic variability estimation of a SI engine fueled with methane, for which a complete set of measured data is available. The numerical results are then compared with the measured data, focusing on the minimum and maximum peak pressure and IMEP. The final result is the calculation of the COV of IMEP, once the COV value of the multi-cycle simulations converges according to the criteria selected. Overall, the numerical results are close to the measured data with a relative difference of COV of IMEP about 25%, showing that the developed code can be used to estimate, at least qualitatively, the main parameters of cyclic variability in engines, with the intention to introduce in a next version more mechanisms that contribute to this phenomenon. The ultimate goal is then to direct attention to the effects of different fuels on the cyclic variations and the associated pollutant emissions.