An important class of non-Newtonian fluids presents a yield stress limit which must be exceeded before significant deformation can occur – the so-called viscoplastic fluids or Bingham plastics. In order to model the stress-strain relation in these fluids, some fitting have been proposed such as the linear Bingham equation and the non-linear Herschel-Bulkley and Casson models. Analytical solutions exist for such models in simple flows. For general flow fields, it is necessary to develop numerical techniques to track down yielded/unyielded regions. This can be avoided by introducing into the models a continuation parameter, which facilitates the solution process and produces virtually the same results as the ideal models by the right choice of its value. Viscoplastic materials like slurries, pastes, and suspension materials have a yield stress, i.e. a critical value of stress below which they do not flow are also called Bingham plastics, after Bingham. Viscoplastic materials can be well approximated uniformly at all levels of stress as liquids that exhibit infinitely high viscosity in the limit of low shear rates followed by a continuous transition to a viscous liquid. This approximation could be made more and more accurate at even vanishingly small shear rates by means of a material parameter that controls the exponential growth of stress. Thus, a new impetus was given in 1987 with the publication by Papanastasiou of such a modification of the Bingham model with an exponential stress-growth term. The new model basically rendered the original discontinuous Bingham viscoplastic model as a purely viscous one, which was easy to implement and solve and was valid for all rates of deformation. The early efforts by Papanastasiou and his co-workers were taken up by the author and his coworkers, who in a series of papers solved many benchmark problems and presented useful solutions always providing the yielded/unyielded regions in flow fields of interest. Since the early 1990s, other workers in the field also used the Papanastasiou model for many different problems. Papanastasiou Papanastasiou in 1987, who took into account earlier works in the early 1960s as well as a well-accepted practice in the modelling of soft solids and the sigmoidal modelling behaviour of density changes across interfaces. He introduced a continuous regularization for the viscosity function which has been largely used in numerical simulations of viscoplastic fluid flows, thanks to its easy computational implementation. As a weakness, its dependence on a non-rheological (numerical) parameter, which controls