In this thesis, an analysis of the single-relaxation-time (SRT) and multi-relaxationtime (MRT) lattice Boltzmann model for 2-D simple flow problems is presented. This analysis provides another method to determine coefficients of equilibrium moments in the MRT model. Derived result indicates that commonly adopted boundary condition scheme for SRT model can be directly used on MRT model without any modification. In addition, when analyzing the equation with the forcing term, a discrete forcing term is discovered from the derived equations as well. This discrete forcing term will generate some artificial velocities and contaminates the final fluid field, if appropriate treatment is not applied. It shows that the discrete force effect can be eliminated by choosing a specific value of the relaxation parameter in both SRT and MRT model. However, the MRT model provides a much flexible option for users to select relaxation time τ and accurate results are expected to be obtained. The D2Q9 MRT lattice Boltzmann model is then utilized to simulate 2D deep cavity flows to locate the critical Reynolds number for cavity with different width depths. Simulation results indicate that for square cavity flow, the critical Reynolds number, where the transition phenomena occur, is located in the region between 8325 and 8350. For cavity flows with aspect ratio 2 and 3, the critical Reynolds number resides at 6025-6050, and 5725-5750, respectively. MRT lattice Boltzmann simulation with GPU technology implementation is also presented in this thesis. 2D square cavity flow and 3D cubic cavity flow problems are selected as the test to investigate the parallel performance. Acceleration up to 20 times speedup can be achieved for both 2D and 3D simulation with current GPU program. However, since only global memory is used in present code, efforts for optimization are still required to improve the overall parallel efficiency.