In this paper, a viscoelasticity-plastic damage constitutive equation for naturally fractured shale is deduced, coupling nonlinear tensile-shear mixed fracture mode. Dynamic perforation-erosion on fluid re-distribution among multi-clusters are considered as well. DFN-FEM (discrete fracture network combined with finite element method) was developed to simulate the multi-cluster complex fractures propagation within temporary plugging fracturing (TPF). Numerical results are matched with field injection and micro-seismic monitoring data. Based on geomechanical characteristics of Weiyuan deep shale gas reservoir in Sichuan Basin, SW China, a multi-cluster complex fractures propagation model is built for TPF. To study complex fractures propagation and the permeability-enhanced region evolution, intersecting and competition mechanisms between the fractures before and after TPF treatment are revealed. Simulation results show that: fracture from middle cluster is restricted by the fractures from side-clusters, and side-clusters plugging is benefit for multi fractures propagation in uniformity; optimized TPF timing should be delayed within a higher density or strike of natural fractures; Within a reservoir-featured natural fractures distribution, optimized TPF timing for most clustered method is 2/3 of total fluid injection time as the optimal plugging time under different clustering modes.