We design and numerically validate a recovery based linear finite element method for solving the biharmonic equation.The main idea is to replace the gradient operator ▽ on linear finite element space by G(▽) in the weak formulation of the biharmonic equation,where G is the recovery operator which recovers the piecewise constant function into the linear finite element space.By operator G,Laplace operator △ is replaced by ▽ · G(▽).Furthermore,the boundary condition on normal derivative ▽u·n is treated by the boundary penalty method.The explicit matrix expression of the proposed method is also introduced.Numerical examples on the uniform and adaptive meshes are presented to illustrate the correctness and effectiveness of the proposed method.