In recent years, a new ultrasonic nondestructive testing that utilizes nonlinear ultrasonic behavior generated in cracks or at material joints has emerged. However, the mechanism of generating subharmonics and higher harmonics used in testing is still not well understood theoretically. In this research, in order to simulate the higher harmonics, dynamic contact problems of a solid-solid interface with nonlinear boundary conditions are investigated by convolution quadrature time-domain boundary element method (CQ-BEM), which can stably analyze wave propagation with small time increments. Numerical results show that higher harmonics are excited by the P-wave incidence to the interface defect, but that only odd order harmonics are generated in case of S-wave incidence.