Inspiration provided by the high strength-to-weight ratios and energy absorption capability of biological materials such as bee honeycomb, feathers, cork, and trabecular bone, a synthetic composite---cellular composite---has been developed and it increasingly attracts significant attention in the scientific community. To computationally investigate the plastic behavior of cellular composites we adopted an elastoplastic constitutive models [1], that captures material anisotropy, tension-compression asymmetry, and nonlinear isotropic-kinematic hardening/softening, to simulate the mechanical behavior of the matrix and inclusion in the composite. In order to compute the stress response exactly for every material point, we further explored the internal symmetries of the model, which relies on Lie algebra and Lie group theory, and develop its return-free integration, which updates stress responses on the yield surface automatically during the plastic phase. The proposed model and the return-free integration method were implemented in the user material subroutine (UMAT) of the finite element software ABAQUS. On the other hand, the representative block approach were adopted to simulate the effective properties of the cellular composites and a proper yield point determination were selected to detect the yield point of cellular composites [2]. In addition, probing paths and pre-loading paths were designed. Using this finite element framework, we conducted the simulation to detect the yield point of cellular composites under different pre-loading paths and then to investigate their yield surface evolution. The simulation results show the isotropic hardening, kinematic hardening, and distortional hardening of cellular composites under different pre-loading paths. Furthermore, according to the simulation results, we extracted the plastic flows of the cellular composites. The influence of microstructure on the yield surface evolution and the plastic flow rule were addressed in this work.