In this paper, we investigate a low-order robust numerical method for the linear elasticity problem. The method is based on a Bernardi–Raugel-like \(\varvec{H}(\textrm{div})\) -conforming method proposed first for the Stokes flows in [Li and Rui, IMA J. Numer. Anal. 42 (2022) 3711–3734]. Therein, the lowest-order \(\varvec{H}(\textrm{div})\) -conforming Raviart–Thomas space ( \(\varvec{RT}_0\) ) was added to the classical conforming \(\varvec{P}_1\times P_0\) pair to meet the inf-sup condition, while preserving the divergence constraint and some important features of conforming methods. Due to the inf-sup stability of the \(\varvec{P}_1\oplus \varvec{RT}_0\times P_0\) pair, a locking-free elasticity discretization with respect to the Lamé constant \(\lambda \) can be naturally obtained. Moreover, our scheme is gradient-robust for the pure and homogeneous displacement boundary problem, that is, the discrete \(\varvec{H}^1\) -norm of the displacement is \(\mathcal {O}(\lambda ^{-1})\) when the external body force is a gradient field. We also consider the mixed displacement and stress boundary problem, whose \(\varvec{P}_1\oplus \varvec{RT}_0\) discretization should be carefully designed due to a consistency error arising from the \(\varvec{RT}_0\) part. We propose both symmetric and nonsymmetric schemes to approximate the mixed boundary case. The optimal error estimates are derived for the energy norm and/or \(\varvec{L}^2\) -norm. Numerical experiments demonstrate the accuracy and robustness of our schemes.
