Abstract: The axisymmetric Mohr Coulomb yield surfaces have edges and corners in the three-dimensional stress space, which leads to difficulties in numerical calculation. Therefore, how to deal with the axisymmetric Mohr-Coulomb criterion efficiently has always been an important research content of the lower bound finite element limit analysis (LB-FELA) method. Firstly, the axisymmetric Mohr-Coulomb criterion is transformed into a set of inequality constraints and three linear equality constraints by introducing the complete plasticity assumption. Then, the inequality constraints are directly transformed into the second-order cone ones, which avoids the smooth approximation of numerical singularities. Finally, the axisymmetric LB-FELA model based on the Mohr Coulomb criterion is transformed into a mathematical optimization one of the second-order cone programming. The linear stress element adopted by the axisymmetric LB-FELA cannot accurately simulate the stress change in the failure region. Therefore, the mesh distribution form has a great impact on the calculation accuracy of the LB-FELA. To solve this problem, an adaptive mesh refinement strategy based on the element yield residual is proposed. By judging the degree of the node stress in the element close to the yield, the elements can be automatically identified and refined in the failure area. By analyzing the stability problems of typical axisymmetric geotechnical projects such as the bearing capacity of circular foundation and the ultimate uplift capacity of vertical anchor plate, it is shown that the proposed method has high calculation efficiency and accuracy, and it has certain theoretical value and application prospect.