A numerical enthalpy method is proposed for solving the multi-dimensional phase change problem by coupling the enthalpy and the moving boundary. Its basic feature is to sepasate appearent and latent heat terms by transforming the standard enthalpy relation and to couple them into the latent heat term by introducing the idea of phase change rate between the axes and simultaneously to correct the heat flux of the inflowing/outflowing phase change control volume element , so that the accuracy of numerical solution is improved and the moving boundary can be accurately followed.