Boundary element analysis of 2D elasticity problems by a new expanding element interpolation method is proposed in this paper. The expanding element is made up based on a traditional discontinuous element by adding virtual nodes along the perimeter of the element. Its shape functions constructed on both source nodes and virtual nodes are referred as fine shape functions and boundary variables are interpolated by the fine shape functions. The internal nodes of the original discontinuous element are referred as source nodes and its shape function as raw shape function. The raw shape functions are used to provide additional constraint equations between variables on virtual nodes and source nodes. The expanding element inherits the advantages of both the continuous and discontinuous elements while overcomes their disadvantages. With the expanding element, both continuous and discontinuous fields on the domain boundary can be accurately approximated, and the interpolation accuracy increases by two orders compared with the original discontinuous element. While the boundary integral equations are collocated at source nodes, the size of the final system of linear equations has not change. In addition, the expanding elements take full advantages of the characteristic that the trial function of the boundary integral equation can be discontinuous. At last, a few numerical examples are presented to verify the accuracy and convergence of the proposed method.