Hydraulic fracturing in combination with horizontal well has been widely used to stimulate shale oil and gas from reservoirs. Numerical simulation of hydraulic fracturing is a great challenge due to the geological complexity of reservoirs. In this paper, a fully coupled fluid flow, plastic deformation and fracture propagation model is presented based on cohesive zone method. An implicit algorithm to discretize and integrate the fully coupled differential equations is presented. Plastic deformation is observed in the vicinity of hydraulic fracture during the stimulation process. The effect of plastic and stimulation parameters on fracture propagation behaviors has been investigated. The energy dissipation induced by the plastic deformation increases the fracture width. The fracture propagation behaviors both in viscosity and toughness dominated regimes are discussed. The evolution formula of non-dimensional fracture length at different regimes is presented by fitting numerical results.