In order to investigate the hydraylic fracture development of the specimens and simulate the cracks drived by fluid flow in rocks, a flow-stress-damage (FSD) model, implemented with parallel Rock Failure Process Analysis code (parallel -RFPA3D), is presented. The numerical code is based on linear elastic damage mechanics on mesoscopic scale and FEM. For simulating the complete progressive 3D failure and macroscopic mechanical behaviors of rock materials, rock properties such as elastic constants, peak strength, and poisson ratio are randomly distributed to reflect the initial random distributed weakness in mesoscopic scale. The FSD model is used to represent the permeability variation at the two stages, that is stress-dependent permeability for pre-failure and deformation-dependent permeability for post-peak stage of rock at the elemental scale. The results of the simulation with 680,000-element cylindrical rock specimen coincide well with reported experimental results and the process of crack from initiation to the instability extensions is captured vividly. The results and the process indicate that the FSD model works well and parallel-RFPA3D incorporated with FSD model is a valid tool of understanding the physical essence of the evolution of fracture with large-scale elements and fluid flows in rocks.