Gas distribution function plays a crucial role in the description of gas flows at the mesoscopic scale. In the presence of non-equilibrium flow, the distribution function loses its rotational symmetricity, making the mathematical derivation difficult. From both the Chapman-Enskog expansion and the Hermite polynomial expansion (Grad's method), we observe that the non-equilibrium effect is closely related to the peculiar velocity space (C). Based on this recognition, we propose a new methodology to construct the non-equilibrium distribution function from the perspective of polynomial expansion in the peculiar velocity space of molecules. The coefficients involved in the non-equilibrium distribution function can be exactly determined by the compatibility conditions and the moment relationships. This new framework allows constructing non-equilibrium distribution functions at any order of truncation, and the ones at the third and the fourth order have been presented in this paper for illustration purposes. Numerical validations demonstrate that the new method is more accurate than the Grad's method at the same truncation error for describing non-equilibrium effects. Two-dimensional benchmark tests are performed to shed light on the applicability of the new method to practical engineering problems.