Initial and subsequent yield surfaces of cellular composites are computationally explored and investigated in this study. The representative block approach was selected and finite element models of cellular composites which consists of one kind of material embedded in hexagonal cell well of another material were established for yield surface detection. After the mesh convergence, the boundary effect and the size effect analysis, proper arrangement of finite element models were determined and the yield surfaces of cellular composites were detected by using the equivalent plastic strain increment approach of yield point determination. Under different types of pre-loading including proportional and rectangular displacement controlled paths, yield points on the subsequent yield surfaces of cellular composites were detected and then corresponding yield functions were estimated according to the three-stage identification of yield surface. According to our computational yield surface detection, the yield surface evolution of cellular composites is explored and phenomena including the Bauschinger effect and the hardening behavior (isotropic, kinematic, rotation, distortional) are observed. Furthermore, the influence of matrix-inclusion material property and microstructure on the yield surface evolution of cellular composites were investigated in this study.