Based on the model for calculating the critical cooling rate of bulk metallic glasses under non-isothermal condition, a new model for calculating the critical cooling rate of alloy systems based on viscosity calculation is put forward in this paper. The critical cooling rates of Zr-Ni-Al-Cu quaternary alloys were calculated and the glass forming range in Zr66.67(NixAlyCuz)33.33 quaternary alloy system is predicted. The results are in good agreement with the experimental values. The calculation results show that there is high glass forming ability in Zr66.67(NixAlyCuz)33.33 quaternary alloy system. The glass forming range is near eutectic point in which critical cooling rates are less than 100 K/s. There is not crystallization in the temperature range higher than 1000K, the fraction of crystals increases rapidly in the temperature range from 980 K to 870 K, and it is not changed in the temperature range below 870 K. In addition, the effects of Al,Cu,Ni mole fraction on the critical cooling rate are analyzed in this paper.