  A new numerical scheme for solving incompressible Bingham flows with variable
density, plastic viscosity and yield stress is proposed. The mathematical and
computational difficulties due to the non-differentiable definition of the
stress tensor in the plug regions, i.e. where the strain-rate tensor vanishes,
is overcome by using a projection formulation as in the Uzawa-like method for
viscoplastic flows. This projection definition of the plastic tensor is coupled
with a fractional time-stepping scheme designed for Newtonian incompressible
flows with variable density. The plastic tensor is treated implicitly in the
first sub-step of the fractional time-stepping scheme and a fixed-point
iterative procedure is used for its computation. A pseudo-time relaxation term
is added into the Bingham projection whose effect is to ensure a geometric
convergence of the fixed-point algorithm. This is a key feature of the
bi-projection scheme which provides a fast and accurate computation of the
plastic tensor. Stability and error analyses of the bi-projection scheme are
provided. The use of the discrete divergence-free velocity to convect the
density in the mass conservation equation allows us to derive lower and upper
bounds for the discrete density. The error induced by the pseudo-time
relaxation term is controlled by a prescribed numerical parameter so that a
first-order estimate of the time error is derived for the velocity field and
the density, as well as the dependent parameters that are the plastic viscosity
and the yield stress.
