We present a modeling and computational approach to study fusion of multicellular aggregates during tissue and organ fabrication, which forms the foundation for the scaffold-less biofabrication of tissues and organs known as bioprinting. It is known as the phase field method, where multicellular aggregates are modeled as mixtures of multiphase complex fluids whose phase mixing or separation is governed by interphase force interactions, mimicking the cell-cell interaction in the multicellular aggregates, and intermediate range interaction mediated by the surrounding hydrogel. The material transport in the mixture is dictated by hydrodynamics as well as forces due to the interphase interactions. In a multicellular aggregate system with fixed number of cells and fixed amount of the hydrogel medium, the effect of cell differentiation, proliferation, and death are neglected in the current model, which can be readily included in the model, and the interaction between different components is dictated by the interaction energy between cell and cell as well as between cell and medium particles, respectively. The modeling approach is applicable to transient simulations of fusion of cellular aggregate systems at the time and length scale appropriate to biofabrication. Numerical experiments are presented to demonstrate fusion and cell sorting during tissue and organ maturation processes in biofabrication.