This article studies grid generation for the calculation of airfoil ice accretion; in particular, it analyzes how to deal with the effect brought by geometrically complex and time-dependent leading edge accretion. As the first step of icing calculation, a grid is generated by solving the elliptic partial differential equations(PDEs). The forcing terms are automatically chosen in the manner of Hilgenstock method such that orthogonality and spacing control on the icing smoothed boundary are achieved. In order to improve the efficiency, quality and reliability of grids, two partition ideas and grid reconstruction methods are introduced and validated simultaneously. Finally, validation of the usefulness of the method developed is completed on two and three dimensional airfoils using the NUAA-ICE3D software, which has integrated the grid modules.