The voxel-based indirect boundary element method (IBEM) using the Laplace-kernel fast multipole method (FMM) is capable of analysing relatively large-scale problems. Furthermore, the voxel-based IBEM is suitable for acceleration using graphics processing units (GPUs). A typical application of the IBEM is the analysis of electrostatic fields for virtual-human anatomical voxel models such as the model called Duke provided by the IT’IS Foundation. An important property of voxel-version Duke models is that they have different voxel sizes but the same structural feature. This property is useful for examining the O(N) and O(D2) dependencies of the calculation times and the amount of memory required by the FMM-IBEM, where N and D are the number of boundary elements and the reciprocal of the voxel side length, respectively. In this study, the O(N) and O(D2) dependencies of the voxel-based GPU-accelerated FMM-IBEM were confirmed by analysing Duke models with voxel side lengths of 5.0, 2.0, 1.0, and 0.5 mm. The finest model comprised 2.2 billion voxels with 61 million square boundary elements, and a linear equation solver on a personal computer with four GPUs required 1,276 s to obtain a solution. In addition, a technique is proposed to improve the convergence performance of the linear equation solver by considering the non-uniqueness of the electric potential, and its effectiveness is demonstrated.