In this paper, a systematic approach is proposed to calculate the torsional rigidity of a circular bar containing multiple circular inclusions. To fully capture the circular geometries, the kernel function is expanded to the degenerate form and the boundary density is expressed into Fourier series. The approach is seen as a semi-analytical manner since error purely attributes to the truncation of Fourier series. By collocating the null-field point exactly on the real boundary and matching the boundary condition, a linear algebraic system is obtained. After obtaining the unknown Fourier coefficients, the solution can be obtained by using the integral representation. Finally, torsion problems are revisited to check the validity of our method. Torsional rigidities for a circular bar with an eccentric inclusion are compared well with the exact solution, BEM data and the Tang's results. Convergence study shows that only a few number of Fourier series terms can yield acceptable results. The torsional rigidities of two limiting case of cavity and rigid inclusion are also obtained using the present approach. Five gains of well-posed model, singularity free, free of boundary-layer effect, exponential convergence and mesh-free approach are achieved. A general-purpose program was developed to determine the torsional rigidity for a circular bar with arbitrary number, radii, positions and shear moduli of circular inclusions.