A numerical procedure about nonlinear flow in three-dimension discrete fracture networks was developed by solving the Reynolds equation and Forchheimer equation using the Galerkin method. The numerical simulation results of three intersecting models were compared with experimental data and a good agreement was observed. Based on further applications of the model to fracture networks with the consideration of Barton’s and Chen’s empirical model, it was found that nonlinear fluid flows in fracture networks could be appropriately described by the Forchheimer equation. The smoother the fracture surfaces and the greater the connectivity of the fracture network are, the lower the hydraulic gradient for the onset of nonlinear flows is.