Abstract. We consider in this paper numerical approximations of phase-field models for twophase complex fluids. We first reformulate the phase-field models derived from an energetic variational formulation into a form which is suitable for numerical approximation and establish their energy laws. Then, we construct two classes, stabilized and convex-splitting, of decoupled time discretization schemes for the coupled nonlinear systems. These schemes are unconditionally energy stable and lead to decoupled elliptic equations to solve at each time step. Furthermore, these elliptic equations are linear for the stabilized version. Stability analysis and ample numerical simulations are presented.Key words. phase-field, liquid crystals, multiphase flows, Navier-Stokes, Allen-Cahn, CahnHilliard, stability AMS subject classifications. 65M12, 65M70, 65Z05 DOI. 10.1137/130921593 1. Introduction. Complex fluids have great practical utility since the microstructure can be manipulated to produce useful mechanical, optical, or thermal properties [4,23]. Numerical studies of multiphase complex fluids is a rich area for research with important applications that can be addressed. Some numerical difficulties are (i) the moving interfaces between various components; (ii) open issues regarding wellposedness of many multiphase models leading to confusion as to whether numerical pathologies are due to the model or the methods; and most important, (iii) nonlinear coupling (hydrodynamics, interface, and microstructure) makes it very difficult to design easy-to-implement, energy stable numerical schemes.For the first two issues above, the diffuse-interface/phase-field models, whose origin can be traced back to [38] and [46], have been proved efficient with much success. A particular advantage of the phase-field approach is that models can often be derived from an energy-based variational formalism (energetic variational approaches), leading to well-posed nonlinear coupled systems that satisfy thermodynamics-consistent energy dissipation laws. This makes it possible to carry out mathematical analysis and design numerical schemes which satisfy a corresponding discrete energy dissipation law. In recent years, the phase-field method has become one of the major tools to study a variety of interfacial phenomena (cf. [17,2,36,19,32,49] and recent review papers [44,22] and the references therein).For phase-field models which satisfy a dissipative energy law, it is especially desirable to design numerical schemes that preserve the energy dissipation law at the