In order to study the influence of layered rock mass on tunnel seismic response, considering the oblique incident angles of seismic motion and the nonlinear contact characteristics at the interface, a numerical simulation method for seismic response of hydraulic tunnel in layered rock mass was proposed. First, based on the 3D viscoelastic artificial boundary conditions and the wave field decomposition theory, the input method of an obliquely incident earthquake in layered rock mass was put forward. It can transform the seismic waves into equivalent nodal forces acting on the nodes of artificial boundaries. In view of dynamic interaction characteristics between interlayers in layered rock mass under seismic action, a dynamic contact force algorithm considering the bond-slip characteristics of the interface was presented. Then, the methods were applied to the anti-seismic stability calculation of the hydraulic tunnel of AZAD PATTAN hydropower station in Pakistan. The calculation is divided into three different working conditions, with vertically incident earthquake, with obliquely incident earthquake but no dynamic contact force, with obliquely incident earthquake and dynamic contact force. The results indicate that the stress and displacement of the tunnel structure under seismic action are greatly affected by the angle of incidence. The existence of interlayer shearing and crushing fracture zone exacerbates the seismic response of the tunnel, resulting in the fact that the failure zone near the interface develops further. After considering the contact effect, the stress and displacement response of the haunch is larger than that of the vault, so that the haunch of lining where the cracking damage first occurs is the weak part of lining structure under the action of earthquakes. The damage zone of the lining mainly distributes in the place where the soft rock passes through and interlayer contacts.