Mathematical models and numerical simulations can provide an essential insight into the mechanisms through which local cell-cell interactions affect tissue-level cell morphology. Among such morphological phenomena, cellular patterns observed in developing sensory epithelia have gained keen attention of researchers in recent years, because they are thought to be of utmost importance for accurate sensory functions. However, most of current computational approaches to cellular rearrangements lack solid mathematical background and involve experimentally unreachable parameters, whereby only weak and ambiguous conclusions can be made based on simulation results. Here we present a simple mathematical model for tissue morphogenesis together with a level set-based numerical scheme for its solution as a tool to rigorously investigate evolving cellular patterns. This combined framework of a model and a numerical method features minimum possible number of physical parameters and guarantees reliability of simulation results, including correct handling of topology changes, such as cell intercalations. In this framework, we adopt the viewpoint of free energy minimization principle, and take cellular rearrangement as a gradient flow of a weighted surface energy associated with cell membrane, where the weights are related to physical parameters of the cells, for example, cell-cell adhesion and cell contractility. We present the applicability of this model to a wide range of tissue morphological phenomena, such as cell sorting, engulfment or internalization. In particular, we stress that this method is the first one to be successful in computationally reproducing the experimentally observed development of cellular mosaic patterns in sensory epithelia. Thanks to its simplicity and reliability, the model is able to capture the essence of biological phenomena, and may give a strong helping hand in deciphering unsolved questions of morphology.