The cohesive zone model is widely used in fracture mechanics. When the fracture process zone (FPZ) in front of the crack tip is too large to be neglected, the nonlinear behavior must be considered. That is to say, in this circumstance the linear fracture mechanics is no longer valid. In order to take into account the nonlinear behavior in FPZ, many fracture models have been proposed, among which, the cohesive zone model (CZM) might be one of the simplest and has been widely used. However, there still remain some problems in the existing numerical methods; for instance, length of the fracture process zone cannot be obtained accurately; dense meshes are required, etc. In order to get over these difficulties, a new analytical singular element is proposed in the present study and further extended into the cohesive zone model for crack propagation problems. In this singular element, the cohesive traction is approximately expressed in the form of polynomial expanding though Lagrange interpolation. The special solution corresponding to each expanding term is specified analytically. Each special solution strictly satisfies the requirements of both differential equations of interior domain and the corresponding traction expanding terms. The real cohesive traction acting on the cohesive crack surface is thus expressed in a natural and strict way. Then the special solution can be transformed into nodal forces of the present singular element. Assembling the stiffness matrix and nodal force into the global FEM system, the cohesive crack problem can be analyzed. An efficient iteration procedure is also proposed to solve the nonlinear problem. Finally, the cohesive crack propagation under arbitrary external loading can be simulated, and the length of FPZ, crack tip opening displacement (CTOD) and other parameters in the cohesive crack problem can be obtained simultaneously. The validity of the present method is illustrated by numerical examples.