The aim of this paper is to deal with nonsteady behaviour of laminar viscous isochoric fluid by a numerical boundary element method. The problem is composed of thermal, kinetic and kinematic descriptions. Due to the presence of convective therms there is also a need for internal integration cells in the domain in addition to boundary discretization, what in a way reduces the advantages of the procedure. As a case study example buoyancy flow has been dealt with in a square cavity, caused by a temperature difference between the left and right vertical walls. Results are presented for variable Rayleigh number values (Ra—Gr'Pr) during the transient phenomenon.