Numerical solution of yield viscoplastic fluid flow is hindered by the singularity inherent to the Herschel-Bulkley model. A finite difference method over the boundary-fitted orthogonal coordinate system is utilized to investigate numerically the fully developed steady flow of non-Newtonian yield viscoplastic fluid through concentric and eccentric annuli. The fluid rheology is described with the Herschel-Bulkley model. The numerical simulation based on a continuous viscoplastic approach to the Herschel-Bulkley model is found in poor accordance with the experimental data on volumetric flow rate of a bentonite suspension. A strict mathematical model for Herschel-Bulkley fluid flow is established and the corresponding numerical procedures are proposed. However, only the case of flow of a Herschel-Bulkley fluid in a concentric annulus is resolved based on the presumed flow structure by using the common optimization technique. Possible flow structures in an eccentric annulus are presumed, and further challenges in numerical simulation of the Herschel-Bulkley fluid flow are suggested.