The lattice Bhatnagar–Gross–Krook (LBGK) model has become the most popular one in the lattice Boltzmann method for simulating the convection heat transfer in porous media. However, the LBGK model generally suffers from numerical instability at low fluid viscosity and effective thermal diffusivity. In this paper, a modified LBGK model, which incorporates the shear rate and temperature gradient in the equilibrium distribution functions, is developed for incompressible thermal flows in porous media at the representative elementary volume scale. With two additional parameters, the relaxation times in the collision process can be fixed at a proper value invariable to the viscosity and the effective thermal diffusivity. In addition, by constructing a modified equilibrium distribution function and a source term in the evolution equation of temperature field, the present model can recover the macroscopic equations correctly through the Chapman–Enskog analysis, which is another key point different from previous LBGK models. Several benchmark problems are simulated to validate the present model with the local computing scheme for the shear rate and temperature gradient, and the numerical results agree well with analytical solutions and/or those well-documented data in previous studies. It is also shown that the present model and the computational schemes for the gradient operators have a second-order accuracy in space, and better numerical stability of the present modified LBGK model than previous LBGK models is demonstrated.
