This study proposes a double-rough-walled fracture model to represent the natural geometries of rough fractures. The rough surface is generated using a modified successive random additions (SRA) algorithm and the aperture distribution during shearing is calculated using a mechanistic model. The shear-flow simulations are performed by directly solving the Navier-Stokes (NS) equations. The results show that the double-rough-walled fracture model can improve the accuracy of fluid flow simulations by approximately 14.99％-19.77％, compared with the commonly used single-rough-walled fracture model. The ratio of flow rate to hydraulic gradient increases by one order of magnitude for fluids in a linear flow regime with increment of shear displacement from 2.2 mm to 2.6 mm. By solving the NS equations, the inertial effect is taken into account and the significant eddies are simulated and numerically visualized, which are not easy to be captured in conventional experiments. The anisotropy of fluid flow in the linear regime during shearing is robustly enhanced as the shearing advances;however, it is either increased or decreased for fluids in the nonlinear flow regime, depending on the geometry of shear-induced void spaces between the two rough walls of the fracture. The present study provides a method to represent the real geometry of fractures during shearing and to simulate fluid flow by directly solving the NS equations, which can be potentially utilized in many applications such as heat and mass transfer, contaminant transport, and coupled hydro-thermo-mechanical processes within rock fractures/fracture networks.