In this work we introduce the development of a three-phase incompressible Navier–Stokes/Cahn–Hilliard numerical method to simulate three-phase flows, present in many industrial operations. The numerical method is then applied to successfully solve oil transport problems, such as those found in the oil and gas industry. The three-phase model adopted in this work is a Cahn–Hilliard diffuse interface model, which was derived by Boyer and Lapuerta (2006). The Cahn–Hilliard model is coupled to the kinetic-energy stable incompressible Navier–Stokes equations model derived by Manzanero et al. (2020). The spatial discretization uses a high-order discontinuous Galerkin spectral element method which yields highly accurate results in arbitrary geometries. An implicit–explicit (IMEX) method is adopted as temporal scheme for the Cahn–Hilliard equation, while Runge–Kutta 3 (RK3) is used for the Navier–Stokes equations. The developed numerical tool is validated with a manufactured solution test case and used to simulate multiphase flows in pipes, including and a three-phase T-shaped pipe intersection.