In the past three decades, lattice Boltzmann modeling, which originated from the lattice gas automata method [9], has been developed into an efficient numerical approach to simulate and study different complex fluid flow problems ranging from turbulent [10,11] and multiphase flows [12-16] to thermal [17,18] and particulate flows [19,20]. For multiphase flows, there are in general four main categories of lattice Boltzmann models (LBMs): the Shan-Chen pseudopotential model [21,22], the free energy model [13,23,24], the colorgradient model [25], and the phase-field model [26,27]. Due to its simplicity and versatility, we apply the pseudopotential model which presents the intermolecular interactions with a density-dependent pseudopotential. The phase separation in this model is achieved by imposing a neighbor attraction between different phases; thus no interface tracking or reconstruction is needed. The original pseudopotential model suffers from the numerical instability when applied to large