In nonlinear ultrasonic testing (NLUT) based on contact acoustic nonlinearity (CAN), higher- and sub-harmonic waves are measured and analyzed for detection of cracks. Theoretical explanation of the harmonic generation is desirable for accurate NLUT. In particular, the steady-state motion of crack faces is supposed to significantly affect behavior of higher- and sub-harmonic wave generation. However, the detail of wave scattering by cracks with CAN has not been clarified yet. Therefore, in order to investigate the detail of steady-state vibration accompanied with contact of crack faces, the authors coupled the boundary element method (BEM) and harmonic balance method. In the proposed method, a concept of the steady-state and nonlinear vibration of crack faces is introduced as asymptotic behavior of vibration that can be recognized after sufficiently elapsed time. The scattering model is for the antiplane shear wave field that causes the crack face friction. Several numerical calculations were performed to examine and verify the accuracy of the proposed formulation. The present results showed almost complete agreement with those obtained by the conventional time-domain BEM. It was also found that the nonlinear effects for the present model were mainly due to the third-order higher-harmonic waves.