Imposing a constant non-zero heat/mass flux at a boundary is of a rather complicated type of problem with severe convergence issues especially in unsteady conditions. The present report, firstly covers the implementation of thermal lattice Boltzmann scheme and the counter-slip thermal energy density approach, to model the constant flux conditions for a typical flow of fluid over heat generating blocks. Analysis of kinetic relations shows that at certain conditions, the inflating density field could result in inconstancy of fluid heat conductivity leading to unacceptable results. Therefore to minimize such effects, density field normalization at the end of each iteration is suggested. Comparing the test case results of the current study with those of macroscopic finite volume approach in the literature, reveals the promising accuracy of the present method. Secondly, the previously unexplored problem of steady and unsteady mixed convection over inline tandem square cylinders with a constant heat flux at the block boundary is investigated. Due to the employment of Neumann thermal boundary condition on the blocks, modified Richardson number is introduced and its consistency with the conventional definition is checked. Detailed results of the study for 10 < Reynolds < 100, blockage ratios of 1/3, 1/4, 1/5 and 1/8 and modified Richardson numbers of 0.0, 0.5 and 1.0 are presented. The intense of the superimposed buoyance force (defined by modified Richardson number) and its impact on the onset of vortex shedding and the global quantities such as lift and drag coefficients and Nusselt number is also carefully explored. (C) 2015 Elsevier Ltd. All rights reserved.