Previous studies of the response of the Moon to gravitational waves have been carried out using analytical or semi-analytical models assuming ideal lunar structures. Such models are advantageous for their high-speed calculation but fail to account for the extremely heterogeneous subsurface and/or interior structures of the Moon. Numerical calculations are needed, but it is challenging to model the topography and lateral heterogeneity of the Moon. In addition, the computational cost is great especially when performing the GW simulation for a long time. As a first step towards overcoming the above difficulties, we employ a two-dimensional finite element method to numerically simulate the lunar response to gravitational waves. We verify our method by comparing our numerical results with those semi-analytical solutions. Based on such comparison, we also analyze the limitation of the two-dimensional simulation. Our work breaks a new way towards the precise simulation of realistic lunar response to gravitational waves in the future and lays down a solid foundation for three-dimensional numerical simulations.