A new phase method for predicting crack propagation path and mechanical response of fiber reinforced composite laminate is proposed. A new three-dimensional (3D) crack surface density function considering material anisotropy is proposed for composite laminates. The explicit relationship between the model penalty parameter (MPP) in a standard anisotropic phase field method and the material properties is solved analytically. This relationship is then used to construct a new expression of driving force in phase field model to account for longitudinal, transverse normal and transverse shear cracks in composite laminates. The proposed method is implemented into the commercial software ABAQUS through user subroutine user defined element (UEL). Numerical methods on both two-dimensional (2D) lamina and 3D angle ply laminate are provided. The present predictions are in good agreement with the experimental and other numerical results.