  We develop an entropy stable two-phase incompressible
Navier--Stokes/Cahn--Hilliard discontinuous Galerkin (DG) flow solver method.
The model poses the Cahn-Hilliard equation as the phase field method, a
skew-symmetric form of the momentum equation, and an artificial compressibility
method to compute the pressure. We design the model so that it satisfies an
entropy law, including free- and no-slip wall boundary conditions with non-zero
wall contact angle. We then construct a high-order DG approximation of the
model that satisfies the SBP-SAT property. With the help of a discrete
stability analysis, the scheme has two modes: an entropy conserving
approximation with central advective fluxes and the Bassi-Rebay 1 (BR1) method
for diffusion, and an entropy stable approximation with an exact Riemann solver
for advection and interface stabilization added to the BR1 method. The scheme
is applicable to, and the stability proofs hold for, three-dimensional
unstructured meshes with curvilinear hexahedral elements. We test the
convergence of the schemes on a manufactured solution, and their robustness by
solving a flow initialized from random numbers. In the latter, we find that a
similar scheme that does not satisfy an entropy inequality had 30% probability
to fail, while the entropy stable scheme never does. We also solve the static
and rising bubble test problems, and to challenge the solver capabilities we
compute a three-dimensional pipe flow in the annular regime.
