Numerical simulations are undertaken for unsteady flows of an ideal Bingham fluid in a circular Couette viscometer. The main difficulties in such simulations are caused by the non-differentiability of the constitutive equation and the need to determine the position and shape of the yield surface separating the yielded zones from the unyielded ones. In this work, these difficulties are overcome by using a numerical method based on variational inequalities, i.e. the augmented Lagrangian/Uzawa method. The start-up and cessation of circular Couette flows of a Bingham fluid are solved numerically assuming that only one of the cylinders is rotating. An improved theoretical upper bound for the stopping time in the case of cessation is derived. The numerical estimates for the stopping time compare well with the theoretical bounds. Moreover, with the adopted method the evolution of the velocity profiles and the locations of yielded/unyielded surfaces are accurately calculated. In flow cessation, we observe an interesting effect, namely the appearance of a small unyielded region adjoined to the outer cylinder shortly before cessation.