A possibility for the simulation of a biomolecular system in an aqueous solvent is to use a continuum model for the solvent. The evaluation of so-called solvation energy coming from the electrostatic interaction between the solute and the water molecules surrounding it is important for a Monte Carlo simulation with this model. In these simulations, we have to deal with a potential problem with jumping coefficient and with a known boundary condition at infinity. One of the advanced ways to solve the problem is to use a multigrid method on locally refined grids around the solute molecule. In this paper, we focus on the error analysis of the numerical solution and the numerical solvation energy obtained for the locally refined grids. Based on a rigorous error analysis via discrete approximation of the Green function, we show the proper way to construct a composite grid, to discretize the discontinuity in the diffusion coefficient and to interpolate the solutions at the interface between the fine and coarse grids. The error analysis developed is confirmed by the numerical experimental results.