  How to develop efficient numerical schemes while preserving the energy
stability at the discrete level is a challenging issue for the three component
Cahn-Hilliard phase-field model. In this paper, we develop first and second
order temporal approximation schemes based on the "Invariant Energy
Quadratization" approach, where all nonlinear terms are treated
semi-explicitly. Consequently, the resulting numerical schemes lead to a
well-posed linear system with the symmetric positive definite operator to be
solved at each time step. We rigorously prove that the proposed schemes are
unconditionally energy stable. Various 2D and 3D numerical simulations are
presented to demonstrate the stability and the accuracy of the schemes.
