Exploring the fluid flow mechanism in rock masses is of great significance for preventing water inrush during excavation of underground constructions such as tunnels. Quantitative descriptions of the hydraulic properties of single rock fracture subject to normal stresses and shear displacement are the basis for understanding the coupled hydro-mechanical processes in fractured rock masses; however, the quantitative relationships among stress, deformation, aperture, inertial coefficients have not been developed in previous works. Granite specimens with single fracture were prepared and flow tests with variable water heads were carried out, in which incremental normal stresses were applied at each fixed shear displacement to characterize the evolution of permeability. The surface morphology of fractures was digitalized using a three-dimensional high-resolution scanning system. A self-designed numerical code was employed to calculate the deformation of fractures under normal stresses based on the framework of variational principles in contact mechanics. The fracture deformation and void space variation under different shear displacements and normal stresses were investigated. By extracting the aperture data and solving Navier–Stokes equations using COMSOL software, a series of numerical simulations were performed to investigate the nonlinear flow behavior of fluids within fractures under different shear displacements and normal stresses. The relationships among shear displacement, normal stress, void space distributions and parameters describing the nonlinear flow were quantitatively analyzed. The results showed that the fracture surface damage areas obtained from the experiment agreed well with the numerical simulation results, which verified the reliability of the deformation calculation code. The normal stress and the shear displacement exhibited a decreasing power function and an increasing exponential function with the fracture aperture, respectively. The increase in the shear displacement resulted in the concentration of contact areas. The inertia coefficient B in the Forchheimer equation and the critical hydraulic gradient Jc could quantitatively characterize the nonlinear flow behavior. B and Jc exhibited decreasing power functions with the shear displacement. The increasing rates and ranges of Jc and B decreased gradually as shear advances. When the shear displacement increased from 2 to 8 mm, the range of Jc decreased from 6.10×10–3 to 1.20×10–3 by a rate of 80.32%. The range of B decreased sharply from 2.97×1014 Pa·s2·m–7 to 2.43×1013 Pa·s2·m–7 by a rate of 91.28%. A similar power function relationship existed between RSD~B and between RSD~Jc. Finally, a predictive function was proposed to quantify the onset of nonlinear fluid flow through fractures.