A series of three-dimensional unsteady Reynolds-averaged Navier-Stokes (RANS) simulations are conducted to investigate the formation and oscillation of stall cells over a NACA 0012 aerofoil. The simulations are conducted with various physical and numerical conditions, such as the Reynolds number, angle of attack, chord-to-span ratio and span-wise mesh resolution. Results show a clear relationship between the oscillation of stall cells and the fluctuation of lift (observed between 17 and 19.5° angle of attack at a high chord Reynolds number of one million). This unsteadiness shows some dependency on the span-wise mesh resolution, which is significant with a medium span-wise mesh resolution (corresponding to 10% of the chord) and moderate with finer resolutions (corresponding to 5% and 2.5% of the chord, respectively). In addition, a proper orthogonal decomposition (POD) method is adopted to further analyse the oscillatory characteristics of stall cells. In particular, it is shown that the first few POD modes have clear spatial patterns corresponding to the profile of stall cells and their time coefficients are correlated with the fluctuation of lift, further confirming the correlation between the stall cell oscillations and the lift fluctuation.