In this study, a thermodynamically consistent constitutive relation is developed for ferroelectric materials. The present model is developed in such a way that an intermediate variable is introduced as a stored variable which accounts for calculating the internal variable. Due to this approach, simulations can be performed without an iterative procedure. Thus, this model improves the computational efficiency. The algorithmic procedure can also handle the evolution of ferroelastic strain as well as the evolution of polarisation simultaneously, without getting into convergence issues. In this work, simulations have been performed to capture various behaviours such as ferroelectric hysteresis with prestress, mechanical depolarisation, ferroelastic hysteresis, and creep behaviour. The simulated results have been compared with the experimental results, and it is observed that the model is good enough to capture various nonlinear hysteresis behaviours. © 2018, Springer-Verlag GmbH Austria, part of Springer Nature.