High-order and conservative phase space direct solvers that preserve the Euler asymptotic limit of the Boltzmann-BGK equation for modelling rarefied gas flows are explored and studied. The approach is based on the conservative discrete ordinate method for velocity space by using Gauss Hermite or Simpsons quadrature rule and conservation of macroscopic properties are enforced on the BGK collision operator. High-order asymptotic-preserving time integration is adopted and the spatial evolution is performed by high-order schemes including a finite difference weighted essentially non-oscillatory method and correction procedure via reconstruction schemes. An artificial viscosity dissipative model is introduced into the Boltzmann-BGK equation when the correction procedure via reconstruction scheme is used. The effects of the discrete velocity conservative property and accuracy of high-order formulations of kinetic schemes based on BGK model methods are provided. Extensive comparative tests with one-dimensional and two-dimensional problems in rarefied gas flows have been carried out to validate and illustrate the schemes presented. Potentially advantageous schemes in terms of stable large time step allowed and higher-order of accuracy are suggested.