Abstract: In consistence with large and small scales in turbulent flows, shape function space can be divided into resolved and unresolved scale spaces in a frame of finite element method. Introducing the same decomposition of the weighting function space, the variational formulations of Navier-Stokes equations can be divided into two systems of equations: resolved- and unresolved-scale equations. Generally, only the resolved-scale equation is computed, and the unresolved scales are modeled. Based on the unresolved-scale equations, an approximate residual-based unresolved-scale modeling is proposed in the present study. The large-scale equations are then computed by substituting the unresolved-scale modeling. The method is called residual-based large eddy simulation, in which unlike in the classical LES a filtering for Navier-Stokes equations is needed, multiscale decomposition is instead used. Numerical simulations of a turbulent channel flow are implemented with in-house codes of the residual-based large eddy simulation. The results show that, with a low number of elements, the mean streamwise velocity obtained using the present method is in agreement with the DNS data in the inner layer, and it is slightly overpredicted in the outer layer. Underprediction of the Reynolds stress by the present method causes a reduction of turbulence intensity transportation from the streamwise direction to the normal direction. Isosurfaces of the streamwise velocity reveals its capability of capturing the large-eddy structures. Meanwhile, low-speed streaks can be clearly observed in the sublayer near the wall.