This work explores novel ideas to improve the accuracy of integral approximation to differential operators (divergence, gradient and Laplacian) in the simulation of 3D thermal buoyant flows with meshless Lagrangian Blobs methods. Basically, we investigate and develop an integral discretization of the differential operators of the field equations, by using convolutions of truncated 3D-Taylor series expansions with a kernel function defined on a compact support around the blob centre of a given particle. This allows to overtake: the irregular distribution of cells in the compact support around the given blob, the deficiency of cells in the compact support due to the presence of a boundary cutting the compact support of nearby blobs. The accuracy and the order of approximation of such discretizations are determined in regular and randomly distorted grids of various sizes, and compared with the widely used particle strength exchange formulations. The analysis of the effects of using the new formulations to solve problems at realistic values of the Grashof number demonstrates the validity and the benefits of the novel findings. (C) 2008 Elsevier Ltd. All rights reserved.