To resolve the issue of fracturing design for multi-stage and multi-cluster fracturing in multi-zone reservoirs, a new efficient algorithm for the planar 3D multi-fracture model was proposed. The model considers fluid flow in the wellbore-perforation-fracture system and fluid leak-off into the rock matrix, and uses a 3D boundary integral equation to describe the solid deformation. The solid-fluid coupling system is solved by an explicit integration algorithm, and the fracture front is determined by the uniform tip asymptotic solutions and shortest path algorithm. The accuracy of the algorithm is verified by the analytical solution of radial fracture, results of implicit level set algorithm and results of organic glass fracturing experiment. Compared with the implicit level set algorithm (ILSA), the new algorithm is much higher in computation speed. The numerical case study is conducted based on a horizontal well in shale gas formation of Zhejiang oilfield. The impact of stress heterogeneity among multiple clusters and perforation number distribution on multi-fracture growth and fluid distribution among multiple fractures are analyzed by numerical simulation. The results show that reducing perforation number in each cluster can counteract the effect of stress contrast among perforation clusters. Adjusting perforation number distribution in each cluster can promote uniform flux among clusters, and the perforation number difference should better be 1-2 among clusters. Increasing perforation number in the cluster with high in situ stress is conducive to uniform fluid partitioning. However, uniform fluid partitioning is not equivalent to uniform fracture geometry. The fracture geometry is controlled by the stress interference and horizontal principal stress profile jointly.