In this thesis, we use Lattice Boltzmann method with diffuse scattering boundary condition to simulate microflows. According the previous works, the cubic form of the equilibrium distribution function, feq, is said 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 for analyzing effects of higher-order term in feq on velocity prediction. In order to predict the slip velocity, Knudsen layer , and Knudsen minimum effect, other modifications are also applied for these models, such as wall functions, regularizations, and combinations of both. Wall functions are applied for the modification of relaxation time, and regularizations are utilized to modify the nonequilibrium part of distribution function fneq. In the simulations, Couette flow is simulated for Kn=0.25, 0.5, 0.75, 1, and Poiseuille flow is simulated for Kn=0.1, 1, 10. As regularization is applied in LBE, the nonlinear behavior in near-wall region is displayed only for D2Q21 model. One of the wall functions, NWF, is found capable of abating the slip velocity but failing to capture Knudsen layer phenomenon. The other wall function, SWF, can not only lower the slip velocity but also predict a nonlinear behavior in near-wall region. In Couette flow simulation, D2Q21 model with NWF+REG is found to give the most accurate prediction of velocity. In Poiseuille flow simulation, results for D2Q21 model with SWF is in good agreement with DSMC data. Finally, Knudsen minimum effect is exhibited for D2Q21 model in flow rate simulation.