The multi-messenger observations of binary merger events have provided us unprecedented opportunities to study the state of matter under extreme conditions, such as those inside neutron stars. For instance, depending on the total mass of the merging binary, the post-merger remnant could undergo collapse to black hole at different timescales, producing very different gravitational wave and electromagnetic signals. The post-merger evolution scenarios have been extensively explored for binary neutron star mergers by numerical relativity simulations and the ways of constraining neutron star EoS via multi-messenger observations have been suggested accordingly. In my talk, I will introduce our modifications to the numerical relativity code SACRA-MPI which allows us to perform full GR-hydrodynamic simulations for mergers of binary bare quark stars (i.e., without crusts and with finite surface density) for the first time. Our simulation results reveal significantly different post-merger evolution scenarios for binary quark stars compared with neutron stars. The prospect of distinguishing quark stars from neutron stars with future merger observations according to thoses differences will also be discussed.