  A high-order numerical method is developed for solving the
Cahn-Hilliard-Navier-Stokes equations with the Flory-Huggins potential. The
scheme is based on the $Q_k$ finite element with mass lumping on rectangular
grids, the second-order convex splitting method, and the pressure correction
method. The unique solvability, unconditional stability, and bound-preserving
properties are rigorously established. The key to bound-preservation is the
discrete $L^1$ estimate of the singular potential. Ample numerical experiments
are performed to validate the desired properties of the proposed numerical
scheme.
