In this thesis, we use Lattice Boltzmann method to simulate microflows. According to the previous work, the cubic form of the equilibrium distribution function, feq, is deemed to be able to improve the velocity prediction in microflow. Hence, we choose quadratic D2Q9 model and three cubic models, D2Q13, D2Q17 and D2Q21 as our bases to analyze effects of higher-order term in feq on velocity prediction. Moreover, modification is also applied to these models, such as Stops’ wall function (SWF). SWF can not only lower the slip velocity but also predict a nonlinear behavior in near-wall region. Here, wall function is applied to the modification of relaxation time. In order to predict the slip velocity, we use several different boundary conditions to simulate microflow, including bounceback boundary condition, diffuse-scattering boundary condition, and β-weighted diffusive-bounceback boundary condition. First of all, we use a constant force in the streamwise direction with periodic boundary condition at the inlet and outlet. The results show that when we use β-weighted diffusive-bounceback boundary condition with SWF, the slip velocity at the wall can be captured correctly. These results are compared with linearized Boltzmann solution data and DSMC results. Finally, Knudsen minimum effect is exhibited for these models in flow rate simulation. Second, we utilize extrapolated formulas for pressure boundary condition at the inlet and outlet. We test two different Knudsen number using LBM with those three different boundary conditions and SWF. These results are compared with Direct Simulation Monte Carlo (DSMC) data. However, both of the pressure distribution and the exit velocity profile simulated by the LBM model deviate from the DSMC data.