A two dimensional model is developed to simulate the process of keyhole formation on aluminum slab during laser drilling by millisecond pulsed laser. In this model, phase change, including melt and vaporization, is considered by not only introducing the equivalent specific heat capacity but also explicitly adding the source terms of gas dynamics in the thermal-hydrodynamic equations and level-set equation firstly. Besides, in order to trace the free surface, a modified level-set method is developed. All possible effects which can impact the dynamic behavior of the keyhole are taken into account, containing gravity, recoil pressure of the metallic vapor, surface tension, and Marangoni effect. This simulation is based on the same experiment condition where 1064nm Nd:YAG single pulsed laser with 3ms pulse width, 18J peak energy and 300μm spot radius is used. The results of the keyhole morphology and the melt ejection trajectory angle are in good agreement with experiment results, indicating the model is reasonable and feasible.