The starting premise of any soft discrete element method simulation, widely used in granular physics and granular mechanics, is the modelling of grain-grain contact force. Most of models often used in the literature including the famous ones by Hertz-Mindlin and Luding, do not present the algorigthy of total elastic potential, or the rate of dissipation which is mainly due to the partially frictional character of the forces. This renders the question of thermodynamic consistency unsettled. A model that possesses explicit expressions for both is proposed here. It is conceptually closely related to the continuum-mechanical theory of granular solid hydrodynarmics (GSH). This theory contains expressions for the total elastic potential and the thermal energy, it accounts for energy conservation and the positivity of entropy production, and it clarifies the equilibrium properties of granular media. All these are lacking (or hidden) in the contact models widely used in the literature. A preliminary calculation shows that the restitution coefficient varies with the impact velocity, which is an added bonus, and demonstrates the model's increased realism. For simplicity, the equations presented in this work are limited to the 2D-case and neglect granular rotations. Nevertheless, the generalization to the 3D-case and the inclusion of granular rotations are carefully discussed, clarifying how to treat rolling and the torsional forces in a thermodynamically consistent fashion. A key point of the present approach, and the major difference to other force models, is the fact that, starting from the characteristic thermodynamic potential, we employ the Onsager reciprocity relation to set up the transport coefficients. The contact forces (usually postulated) are then derived from them. This difference is both conceptually and methodologically relevant. We discussed in detail off-diagonal transport coefficients, especially the so called gear ratio that is particular to granular matter. It reflects the difference between the elastic and the total strain, and is closely related to the slip movement of contact surface, which occur during shear, rolling and torsional deformations. It is relevant to both the macroscopic GSH scales, and the mesoscopic granular scale.