It is known that velocity fields computed by using an immersed boundary-lattice Boltzmann method (IB-LBM) with a single-relaxation time (SRT) show unphysical distortion when the relaxation time, τ, is high. The authors proposed an immersed boundary-finite difference lattice Boltzmann method (IB-FDLBM) using SRT to predict liquid-solid flows. In simulations with IB-FDLBM, numerical errors in the velocity fields appear as in IBLBMs when τ is high. A two-relaxation time (TRT) collision operator is therefore implemented into IBFDLBM in this study to reduce numerical errors at high τ. Simulations of circular Couette flows show that the proposed method gives accurate predictions at high τ, provided that the magic parameter, which is a function of the relaxation times, is less than unity. In addition, predicted drag coefficients of a circular cylinder and a sphere at low Reynolds numbers show reasonable agreements with theoretical solutions and measured data. NOMENCLATURE A negative viscosity term c lattice speed, c = 1 c discrete particle velocity, c = (cx, cy) CD drag coefficient D particle diameter f particle velocity distribution function f eq equilibrium distribution function f + symmetric part of distribution function f anti-symmetric part of distribution function f eq,+ symmetric part of equilibrium distribution function f eq,anti-symmetric part of equilibrium distribution function F external force, F = (Fx, Fy) FD drag force Gi direct forcing term for fi l norm Nm number of Lagrangian points at immersed boundary p pressure Q number of discrete particle velocities Re Reynolds number Rin radius of inner cylinder Rout radius of outer cylinder t time u fluid velocity, u = (u, v) uθ azimuthal velocity component uP velocity of solid body uT analytical azimuthal velocity component U0 free stream velocity Uθ rotation velocity of circular cylinder W weighting function x Eulerian coordinates, x = (x, y) XL coordinates of Lagrangian point, XL = (XL, YL) Δ domain of the smoothed delta function δ smoothed-delta function δh one-dimensional smoothed-delta function ΔS area segment of solid body Δt time step size ΔV computational cell volume Δx lattice width in the x direction Δy lattice width in the y direction Δz lattice width in the z direction γ Euler constant e Knudsen number Λ magic parameter μ viscosity ν kinematic viscosity ρ density τ+ relaxation time for f + in two-time relaxation model τ− relaxation time for f in two-time relaxation model Ω collision operator Subscript i direction of discrete particle velocity i direction opposite to i I,J,K indexes of lattice point L index of Lagrangian node Superscript n discrete time INTRODUCTION The lattice Boltzmann method is now regarded as one of the promising methods for simulating fluid flows. Due to its simplicity and suitability for parallel computation, it has been widely used for predicting various flows such as turbulent flows (Martinez et al., 1994) and two-phase flows (Shan and Chen, 1993). In particular, the combination of the immersed boundary method and lattice Boltzmann method (IB-LBM) using a single relaxation time (SRT) has been adopted to reasonably predict liquidsolid two-phase flows (Feng and Michaelides, 2004; Feng and Michaelides, 2005; Feng and Michaelides, 2009; Dupuis et al., 2008). However, Le & Zhang (2009) carried out simulations of circular Couette flows and pointed out that large non-physical velocity distortion in the vicinity of immersed boundaries is caused when the relaxation time, τ, is high.