In this paper, we investigate the properties of binary neutron stars (BNSs) and their mergers by combining population synthesis models for binary stellar evolution (BSE) with cosmological galaxy formation and evolution models. We obtain constraints on BSE model parameters by using the observed Galactic BNSs and local BNS merger rate density (R 0 ) inferred from Gravitational Wave (GW) observations, and consequently estimate the host galaxy distributions of BNS mergers. We find that the Galactic BNS observations imply efficient energy depletion in the common envelope (CE) phase, a bimodal kick velocity distribution, and low mass ejection during the secondary supernova explosion. However, the inferred R 0 does not necessarily require an extremely high CE ejection efficiency and low kick velocities, different from the previous claims, mainly because the latest inferred R 0 is narrowed to a lower value (320 +490 −240 Gpc −3 yr −1 ). The BNS merger rate density resulting from the preferred model can be described by R(z) ∼ R 0 (1 + z) ζ at low redshift (z 0.5), with R 0 ∼ 316-784 Gpc −3 yr −1 and ζ ∼ 1.34-2.03, respectively. Our results also show that R 0 and ζ depend on settings of BSE model parameters, and thus accurate estimates of these parameters by future GW detections will put strong constraints on BSE models. We further estimate that the fractions of BNS mergers hosted in spiral and elliptical galaxies at z ∼ 0 are ∼ 81-84% and ∼ 16-19%, respectively. The BNS merger rate per galaxy can be well determined by the host galaxy stellar mass, star formation rate, and metallicity, which provides a guidance in search for most probable candidates of BNS host galaxies.