Clustering of flagellated microswimmers such as sperm is often mediated by hydrodynamic interactions between them. To better understand the interaction of microswimmers in viscoelastic fluids, we perform two-dimensional simulations of two swimming sheets, using a viscoelastic version of the smoothed dissipative particle dynamics method that implements the Oldroyd-B fluid model. Elasticity of sheets (stiff versus soft) defines two qualitatively different regimes of clustering, where stiff sheets exhibit a much more robust clustering than soft sheets. A formed doublet of soft sheets generally swims faster than a single swimmer, while a pair of two stiff sheets normally shows no speed enhancement after clustering. A pair of two identical swimmers is stable for most conditions, while differences in the beating amplitudes and/or frequencies between the two sheets can destroy the doublet stability. Clustering of two distinct swimmers is most stable at Deborah numbers of De = ≈ 1 ( is the relaxation time of a viscoelastic fluid and is the beating frequency), in agreement with experimental observations. Therefore, the clustering of two swimmers depends non-monotonically on De. Our results suggest that the cluster stability is likely a dominant factor which determines the cluster size of collectively moving flagellated swimmers.