Abstract: Galerkin meshfree methods with arbitrary order smooth shape functions exhibit superior accuracy advantages in structural analysis. However, the smooth meshfree shape functions generally do not have the interpolatory property and thus the enforcement of essential boundary conditions in Galerkin meshfree methods is not trivial. The variationally consistent Nitsche’s method shows very good performance regarding convergence and stability and is widely used to impose essential boundary conditions. In this work, a consistent and efficient method is proposed to impose meshfree essential boundary conditions. The proposed method is based upon the Hellinger-Reissner (HR) principle, where the displacements are represented by the conventional meshfree shape functions and the stresses are approximated by reproducing kernel smoothed gradients defined in each background integration cells. The resulting meshfree discrete equations share almost identical forms with those derived from Nitsche’s method. It is shown that the stabilized term in Nitsche’s method is a natural outcome from the HR variational principle, but there is absolutely no need to use any artificial parameter to maintain the coercivity of stiffness matrix. Moreover, under the reproducing kernel gradient smoothing framework, the costly derivatives of conventional meshfree shape functions are completely avoided and the integration constraint is automatically fulfilled. Numerical results demonstrate that the proposed approach and Nitsche’s method yield comparable solution accuracy, nonetheless, much higher efficiency is observed for the proposed methodology that imposes the essential boundary conditions for Galerkin meshfree formulation via the HR variational principle.