Hydraulic fracturing is a widely used production stimulation technology for conventional and unconventional reservoirs. The cohesive element is used to explain the tip fracture process. In this paper, the cohesive zone model was used to simulate hydraulic fracture initiation and propagation at the same time rock deformation and fluid exchange. A numerical model for fracture propagation in poro-viscoelastic formation is considered. In this numerical model, we incorporate the pore-pressure effect by coupling fluid diffusion with shale matrix viscoelasticity. The numerical procedure for hydraulically driven fracture propagation uses a poro-viscoelasticity theory to describe the fluid diffusion and matrix creep in the solid skeleton, in conjunction with pore-pressure cohesive zone model and ABAQUS was used as a platform for the numerical simulation. The simulation results are compared with the available solutions in the literature. The higher the approaching angle, the higher the differential stress, tensile stress difference, injection rate, and injection fluid viscosity, and it will be easier for hydraulic fracture crossing natural fracture. These results could provide theoretical guidance for predicting the generation of fracture network and gain a better understanding of deformational behavior of shale when fracturing.Energies 2019, 12, 1254 2 of 12 side, it may cause additional resources to leak. Therefore, understanding of how the hydraulic fractures interact with natural fractures is necessary and is beneficial to our future simulations and improvements. Overall, continuous research on hydraulic fracturing is necessary and worthwhile. Hydraulic fracture intersection with natural fracture was investigated through both laboratory and mine-back experimental approaches [4][5][6][7]. When hydraulic fractures intersect with natural fractures, there are three different developments. Hydraulic fractures cross over natural fractures without changing the volume or shape of natural fractures; hydraulic fractures divert into the natural fracture and expand the volume of the flow path; hydraulic fractures enter natural fractures and form a complex fracture network structure.Biot first proposed viscoelastic theory to study the linear viscoelasticity and anisotropy of porous media [8]. The influence of the permeability on the deformation and the settlement problem of the loaded column was discussed. A micromechanical method for explaining the theory of pore viscoelasticity of Biot is proposed by Abousleiman et al. [9]. The problem of drilling under plane-strain deformation is discussed and the effects of combined poro-and visco-elasticity are investigated. Simakin and Ghassemi [10] put forward another poro-viscoelastic model by considering the relaxation of the deviatoric stress and the symmetric effective stress. These constitutive models can be used to simulate the coupling process between fluid flow and creep deformation of matrix rocks, but they cannot be directly applied to simulate hydraulic fracturing, especially the ...