The nonlinear evolution of hydrodynamically unstable flames is studied numerically within the context of a hydrodynamic model, where the flame is confined to a surface separating the fresh mixture from the hot combustion products. The numerical scheme uses a variable-density Navier–Stokes solver in conjunction with a level-set front-capturing technique for the numerical treatment of the propagating front. Unlike most previous studies that were limited to the weakly nonlinear Michelson–Sivashinsky equation valid for small density changes, the present work places no restriction on the density contrast and thus elucidates the effect of thermal expansion on flame dynamics. It is shown that the nonlinear development leads to corrugated flames with a transverse dimension that is significantly larger than the wavelength corresponding to the most amplified disturbance predicted by the linear theory, and which is determined by the overall size of the system. The flame structure consists of wide troughs and relatively narrow cusp-like crests, and propagates ‘steadily’ at a constant speed, larger than the speed of a planar flame. The propagation speed increases as the cells widen, but eventually reaches a constant value that remains independent of the mixture's composition and of the transverse length. The dependence of the incremental increase in speed on thermal expansion is found to be nearly linear; for realistic values of thermal expansion it may be as large as 15% to 20%. In sufficiently large domains the dynamics is found to be extremely sensitive to background noise that may result, for example, from weak turbulence. Small-scale wrinkles appear sporadically on the flame surface and travel along its surface, causing a significant increase in the overall speed of propagation, up to twice the laminar flame speed.