This paper proposes a soil constitutive model suitable for geotechnical earthquake simulation. The Hardin curve of soil dynamic stress-strain relationship and its Pyke’s modification for irregular cyclic loading are used. The boundary surface is established on the deviatoric stress plane according to von Mises criterion. After connecting reverse loading point and the current stress point for boundary surface projection, a plastic hardening modulus is derived regarding projection ratio as hardening parameter. The specific increment formulation is given. This constitutive model is implemented by customizing the secondary development of materials in finite element method software ADINA and verified with dynamic triaxial test. The comparison of numerical simulation and test results illustrates that this constitutive model can truthfully reflect the stress-strain relation of soil. For seismic responses of a site, a 2D numerical simulation is executed using the boundary surface constitutive model in ADINA. The results are compared with those of SHAKE91 to illustrate the rationality of adopting this constitutive model to geotechnical earthquake engineering problems.