To characterize fractured rock masses, self-developed programs are utilized to generate fracture networks. A nonlinear flow model considering the shear dilatancy effect is established, and a numerical solution method for modelling nonlinear flow in fractured rock masses during shear is proposed on the basis of extended finite element analysis. The contour plots reveal distinct patterns in the water pressure and flow distributions within fractures. The reduction in the lateral pressure coefficient and increase in the shear stiffness of the joints facilitate a more homogeneous distribution of the water pressure gradient. Under the same vertical stress, increasing the lateral pressure coefficient or decreasing the shear stiffness leads to a more pronounced shear dilatancy effect on fractures. Consequently, an increase in fracture aperture and permeability occurs, and the flow of the fractured rock mass is enhanced. With the same vertical stress, an increase in horizontal stress and a decrease in shear stiffness lead to a gradual reduction in the linear and nonlinear coefficients of Forchheimer’s law. Specifically, the influence of the lateral pressure coefficient on the linear coefficient is greater than that on the nonlinear coefficient.