© 2017 Elsevier Ltd Recently, the ultrasonic testing (UT) based on the contact acoustic nonlinearity (CAN) has attracted notice as a new technique for detection of closed cracks which cannot be detected by using the linear UT (LUT). In the nonlinear UT (NLUT), detection and size measurement of flaws are conducted via the frequency spectrum analysis of nonlinear ultrasonic waves which consist of higher- and sub-harmonic waves. Although the mechanism of higher-harmonic generation due to CAN has been understood mostly from theoretical and experimental points of view, that of sub-harmonic generation has not been revealed yet. In addition, there are few numerical and theoretical studies on the phenomenon of sub-harmonic generation due to the CAN. In this paper, the boundary integral equation for two-dimensional (2-D) elastic wave scattering by cracks is formulated and numerically solved to investigate the behavior of the sub-harmonic waves. The clapping motion and dynamic friction on crack faces are modeled by considering contact boundary conditions, and their interaction with the characteristic of frequency response in the corresponding linear system is investigated. In this study, time-domain numerical simulations are implemented for both an interior crack in an unbounded elastic solid and a surface breaking crack in an elastic half-space. Some frequency response analyses for a linear system corresponding to the time-domain nonlinear simulations are also performed by using the frequency-domain boundary element method. From obtained numerical results, the causes of sub-harmonic generation are discussed.