In order to effectively simulate the process of random propagation of hydraulic fractures in fractured shale reservoirs, a new method of random propagation of hydraulic fractures based on the finite element mesh embedded with zero-thickness cohesive elements is proposed. This new method is based on the topological data structure of element nodes and the splitting mode of mesh nodes. The accuracy and effectiveness of the random propagation method are verified by comparing with the analytical solution of KGD model and two kinds of laboratory experiments. Meanwhile, the influences of horizontal in-situ stress difference and reservoir heterogeneity on the process of random propagation of hydraulic fractures are evaluated by running numerical examples. The results show that: (1) the new method makes up for the deficiency that the cohesive element built-in ABAQUS platform can not effectively simulate the random propagation of hydraulic fractures; (2) under a higher horizontal in-situ stress difference condition, the stronger the heterogeneity of a shale reservoir is, the easier it is to reopen a high-angle natural fracture intersecting with hydraulic fractures. The proposed method can accurately describe the random propagation behavior of complex hydraulic fractures, and thus provide a novel means for numerical simulation of naturally fractured shale reservoirs.