The Preisach approach was used to describe the widely observed quadratic field dependence of the direct longitudinal piezoelectric coefficient and its associated hysteresis loop. In this perspective, a four parameters distribution function was put forward and refined experimentally using bias stress variations. It permitted one to fully describe the stress bias dependence of the nonlinear coefficients. Considering the hysteresis loop associated with the proposed distribution function, modeling of the experimental data was made possible by the addition of a viscous term. Moreover, the distribution parameters extracted from the loops were in strict agreement with the values obtained either directly from experiments or from bias stress dependence characterization. This new set of results is considered as a further confirmation of the applicability of the Preisach model to piezoelectricity. (C) 2001 American Institute of Physics.