@Migorski2013; @Goeleven2003]. In this paper, we focus on boundary hemivariational inequality problems, which originate from the mathematical model of elastic contact system. To solve this kind of hemivariational inequalities, a finite element method had been implemented [@Haslinger1999] while the thorough numerical analysis has not been established until recent. In [@Han2017], Han et al. derived a Céa’s inequality in an abstract framework, and figure out the influence of solution’s regularity to the numerical computation. To our knowledge, proper assumptions will balance the solvability and generality of mathematical models, and this is extremely important in nonlinear problems. Hence, we adopt the assumptions in [@Han2017] to prove our main results.

To our knowledge the study on hemivariational inequalities with the coefficients setting in small periodic configurations is few. The homogenization result could be found in [@Liu2008]. However, the result or H-convergence property does not provide a priori convergence rate which is pivotal in numerical analysis. It explains why we need to build an $O(\epsilon^{1/2})$ estimation. In the following sections, we set model problem in a scalar form merely for the simplicity of symbols, and the extension to elastic system will be straight.

The rest of paper is organized as follows. In \[sec:preliminaries\] we 