Many practical flow problems are defined with the circular boundary. Fluid flows within a circular boundary are however susceptible to a singularity problem when the cylindrical coordinates are employed. To remove this singularity a method has been developed in this study which uses the local harmonic functions in discretization of derivatives as well as interpolation. This paper describes the basic reason for introducing the harmonic functions and the overall numerical methods. The numerical methods are evaluated in terms of the accuracy and the stability. The Lamb-dipole flow is selected as a test flow. We will see that the harmonic-function method indeed gives more accurate solutions than the conventional methods in which the polynomial functions are utilized.