The duct ยฏow of Bingham plastic ยฏuids is analysed with the variational inequality-based ยฎnite element method. The problem of tracking the yield surface is solvable through the regularization technique which can be easily incorporated into the existing ยฎnite element code. The existence theorem of this method was established through the theory of variational inequalities. A small positive constant is added to the second shear rate invariant, resulting in an apparent viscosity of ยฎnite magnitude in the unyielding plug zone. This makes the minimization of the non-differential variational integral possible. In order to achieve convergence at small regularization parameter, a zero-order continuation is employed. It is also shown that a ยฎne tessellation of the ยฏow domain is necessary for tracking the yield surfaces unambiguously. Two classes of duct ยฏow, namely axial ยฏows in eccentric annuli and in an L-shaped duct, were investigated. In both cases it was easy to show the presence of the mobile plugs around the duct centres from the axial velocity proยฎles; however, the stagnant plugs at the narrow side in eccentric annuli with large eccentricity and near the apex of right-angled corners in an L-shaped duct could only be identiยฎed from the calculation of the distributions of the second shear rate or shear stress invariant.