ABSTRACT: A new form of DDA method, which replaces traditional linear polynomial series with natural neighbor interpolants, is developed. In its numerical models, some meshfree nodes need be distributed over each block’ domain, then the Sibson natural neighbor interpolation is adopted. Moreover, the function of modeling block materials’ elasto-plastic nonlinearity is developed. In each time step of computation, the blocks’ stress fields are checked and adjusted to satisfy the assumed elasto-plastic constitutive models. To solve possible difficulties of computation convergence, controlling rational time step length or partly releasing rigid contact constraints are supposed.