This work concerns the modeling and simulation of a mixture of $N$ ($N\geqslant 2$) immiscible incompressible fluids, with possibly very different physical properties such as densities, viscosities and pairwise surface tensions. We call such problems N-phase flows, where N refers to the number of fluid components in the system. N-phase flows involve multiple types of fluid interfaces, multiple pairwise surface tensions and three-phase lines, as well as the interactions among multiple fluid components. The applications of such flows are numerous. A dramatic example is the BP oil spill in the Gulf of Mexico in 2010 (involving at least three fluid components: oil, water, gas/air). These problems are also encountered in our everyday environment such as showers, non-sticking surfaces, sprays, rain drops, air bubbles and ocean waves.

Simulation of N-phase flows presents two essential challenges: (1) how to model and formulate the N-phase system in a physically consistent manner; (2) how to solve the system in a numerically efficient manner.

We present in this work an N-phase model and a physical formulation within the phase field framework that is thermodynamically consistent. The model/formulation satisfies the conservations of mass and momentum, the second law of thermodynamics, and the Galilean invariance principle. We also present a numerical algorithm for solving the N-phase system, which has several features that make it computationally very efficient. In particular, when solving the system of ($N-1$) strongly-coupled fourth-order phase field equations within each time step, our algorithm involves only the solution of $2(N-1)$ {\em de-coupled} Helmholtz-type equations. We present numerical experiments for several problems involving multiple fluid phases, large density contrasts and large viscosity contrasts. In particular, we compare with the Langmuir-de Gennes theory to demonstrate that the presented method produces physically accurate results for multiple fluid phases.