In this study, laminar natural convection heat transfer to Bingham plastic fluids from two differentially heated isothermal cylinders confined in a square enclosure (with isothermal walls) has been investigated numerically. The governing partial differential equations have been solved over the ranges of the dimensionless parameters, namely, Rayleigh number, 102to 106, Prandtl number, 10 to 100, and Bingham number, 0.01 to 100, for seven locations of inner cylinders as ±0.25, ±0.2, ±0.1 and 0. These values correspond to the range of Grashof number varying from 10 to 105. The detailed flow and temperature fields are visualized in terms of the streamlines and isotherm contours. Further insights are developed by examining the iso-shear rate contours and the yield surfaces delineating the fluid-like and solid-like regions. The corresponding heat transfer results are analyzed in terms of the distribution of the local Nusselt number along the cylinder surface together with its surface averaged value as functions of the Rayleigh number, Prandtl number, Bingham number, and positions of the cylinders. It is found that the average Nusselt number increases with the increasing values of the Rayleigh number and decreases with the increasing Bingham number. For sufficiently large values of the Bingham number, the average Nusselt number reaches its asymptotic value wherein heat transfer takes place solely by conduction. Based on the present numerical results, simple correlations for the prediction of the average Nusselt number and the limiting Bingham number have been developed. Also, a dimensionless criterion denoting the cessation of convection regime is outlined for this configuration. [ABSTRACT FROM PUBLISHER]