Abstract: As site response analysis work in the field of earthquake engineering increasingly advances towards large-scale, complex and integrated applications, the accuracy and efficiency optimization problems of the algorithm have regained significant focus. This paper comprehensively investigates the computational performance of finite-element mesh models constructed by irregular elements in site response analysis through numerical experiments and theoretical analyses. The main findings are as follows: (1) The commonly used “one-tenth wavelength” mesh size exhibits certain inaccuracies; a mesh size of one-fifteenth to one-twentieth wavelength can achieve higher precision. (2) Triangular and quadrilateral elements demonstrate equivalent accuracy at identical mesh sizes. Similarly, when comparing asymmetric, elongated elements with symmetric, well-proportioned elements of equal long-edge dimensions, the former exhibits higher accuracy along its short-edge direction. (3) The stable time step for explicit time integration schemes correlates with the minimum characteristic dimension of elements—such as the short side of rectangles, the minimum height of triangular elements, or the perpendicular distance between the nearest two edges of irregular quadrilaterals. The time step for irregular finite element meshes is governed by the smallest mesh size in the whole discrete model. (4) Asymmetric, elongated elements and localized fine-scale meshes in irregular models severely constrain the limit time step of computational stability, significantly reducing computational efficiency. Therefore, adopting symmetric, well-proportioned element shapes and improving mesh size uniformity within each region of identical wave velocity can substantially enhance computational efficiency. Considering that site model data in practical seismic response analyses often lack high precision, we recommend utilizing the principles outlined in this study to establish highly efficient discrete models.