The underlying physical mechanism of the residual-based large eddy simulation (LES) based on the variational multiscale (VMS) method is clarified. Resolved large-scale energy transportation equation is mathematically derived for turbulent kinetic energy budget analysis. Firstly, statistical results of benchmark turbulent channel flow at Re_τ=180 obtained using a coarse mesh are compared with the results obtained by the classical LES with the Smagorinsky and dynamic subgrid stress (SGS) model. The present LES shows an advantage in predicting the statistical results of the incompressible turbulent flows. Secondly, the contributions of the unresolved small-scale presentation terms (Term I-IV in Eq. (10)) to the turbulent kinetic dissipation are analysed for the VMS method. The results show that the turbulent kinetic dissipation provided by the numerical diffusion in the VMS method is smaller in the inner layer, larger in the outer layer of the channel flow than those by the Smagorinsky and dynamic SGS model. The turbulent kinetic dissipation in the VMS method is mainly given by the numerical diffusion provided by one of the “cross-stress” terms (Term I, same as the stabilization term in the SUPG method) and LSIC term (Term IV). The other one of the “cross-stress” terms (Term II) gives rise to the positive turbulent kinetic energy budget, and does not dissipate the turbulent kinetic energy. The so-called “Reynolds stress” term (Term III) dissipates the turbulent energy but provides a very small numerical diffusion. Finally, on the basis of the turbulent kinetic energy dissipation analysis, a new residual-based stabilized finite element formulation is proposed by modifying the large-scale equation in the VMS method. Numerical experiments of 2D lid-driven cavity flow and 3D incompressible turbulent channel flow are tested to validate the proposed formulation. It is shown that all the stabilization terms in the proposed formulation produce additional numerical diffusions and physically increase the total turbulent kinetic dissipation. Consequently, an apparent improvement in both the first-order and second-order statistical quantities are pursued by the new stabilized finite element formulation.