We explore the confinement-deconfinement phase transition (PT) of the first order (FO) arising in SU(N) pure Yang-Mills theory, based on Polyakov loop models (PLMs), in light of the induced gravitational wave (GW) spectra. We demonstrate that the PLMs with the Haar measure term, involving models successful in QCD with N = 3, are potentially incompatible with the large N scaling for the thermodynamic quantities and the latent heat at around the criticality of the FOPT reported from the lattice simulations. We then propose a couple of models of polynomial form, which we call the 4-6 PLM (with four- and six-point interactions among the basic PL fields which have center charge 1) and 4-8 PLM (with four- and eight-point interactions), and discuss how such models can naturally arise in the presence of a heavy PL with charge 2. We show that those models give the consistent thermodynamic and large N properties at around the criticality. The predicted GW spectra are shown to have high enough sensitivity to be probed in the future prospected interferometers such as BBO and DECIGO.