Flow past an isolated circular cylinder and two cylinders in tandem is numerically simulated, under the influence of buoyancy aiding and opposing the flow. A modified velocity correction method is employed, which has second order accuracy in both space and time. The influence of buoyancy on the temporal fluid flow patterns is investigated, with respect to streamlines, isotherms and streaklines. Comparisons are made with respect to mean center line velocities, drag coefficients, Strouhal number and streakline patterns. Degeneration of naturally occurring Karman vortex street into a twin eddy pattern is noticed in the Reynolds number (Re) range of 41-200, under buoyancy aided convection. On the contrary, buoyancy opposed convection could trigger vortex shedding even at a low Re range of 20-40, where only twin eddies are found in the natural wake. Temporal evolution of unsteady eddy patterns is visualized by means of numerical particle release (NPR). Zones of vortex shedding and twin vortices are demarcated on a plot of Richardson number against Strouhal number. Root mean square (RMS) lift coefficients (CL,RMS) and average drag coefficient are obtained as a function of Richardson number (Ri).