The authors apply the Galenkin variational equation to solve the integral equation with hyper singularity, which can be deduced from the double layer solution for Neumann problem of Laplace equation. The scheme of partial integration in the sense of distributions is introduced to reduce the hyper singularity integral into a weak one with the boundary rotation of unknown function. The numerical implementation with linear boundary elements is presented. The numerical examples illustrate the feasibility and efficiency of the method.