  In this paper, we present an efficient numerical algorithm for solving the
time-dependent Cahn--Hilliard--Navier--Stokes equations that model the flow of
two phases with different densities. The pressure-correction step in the
projection method consists of a Poisson problem with a modified right-hand
side. Spatial discretization is based on discontinuous Galerkin methods with
piecewise linear or piecewise quadratic polynomials. Flux and slope limiting
techniques successfully eliminate the bulk shift, overshoot and undershoot in
the order parameter, which is shown to be bound-preserving. Several numerical
results demonstrate that the proposed numerical algorithm is effective and
robust for modeling two-component immiscible flows in porous structures and
digital rocks.
