A theoretical model and numerical methods were developed for testing different bone internal remodelling stimuli. The keystone of the study was the formulation of a stimulus based on the mechanical invariants of the stress tensor, which took into account bone non-homogeneity and anisotropy. A non-site specific remodelling rate equation was then used for the apparent density whereas anisotropy was fixed and evaluated from anatomic observations. An node-based semi-implicit algorithm with adaptive stepsize was implemented for solving the evolution equation. To preclude numerical artifacts (non-convergence, instability), a phase space description was proposed. As an illustration, the evolution of apparent density distribution surrounding the femoral stem after a Total Hip Replacement was simulated. Three stimuli were tested: the strain energy density stimulus, the octahedral shear stress stimulus, and an anisotropic plastic yield stress stimulus.