Atomization spray of non-Newtonian liquid plays a pivotal role in various engineering applications, especially for the energy utilization. To operate spray systems efficiently and well understand the effects of liquid rheological properties on the whole spray process, a comprehensive model using Euler-Lagrangian approaches was established to simulate the evolution of the atomization spray for viscoelastic liquid. Based on the Oldroyd model, the viscoelastic linear dispersion relation was introduced into the primary atomization; an extended viscoelastic version of Taylor analogy breakup (TAB) model was proposed; and the coalescence criteria was modified by rheological parameters, such as the relaxation time, the retardation time and the zero shear viscosity. The predicted results are validated with experimental data varying air-liquid mass flow ratio (ALR). Then, numerical calculations are conducted to investigate the characteristics of viscoelastic liquid atomization process. Results showed that the evolutionary trend of droplet mean diameter, Weber number and Ohnesorge number of viscoelastic liquids along with axial direction were qualitatively similar to that of Newtonian liquid. However, the mean size of polymer solution increased more gently than that of water at the downstream of the spray, which was beneficial to stable control of the desirable size in the applications. As concerned the effects of liquid physical properties, the surface tension played an important role in the primary atomization, which indicated the benefit of selecting the solvents with lower surface tension for finer atomization effects, while, for the evolution of atomization spray, larger relaxation time and zero shear viscosity increased droplet Sauter mean diameter (SMD) significantly. The zero shear viscosity was effective throughout the jet region, while the effect of relaxation time became weaken at the downstream of the spray field.