Mesh parameterization is the process of finding one-to-one mapping between an input mesh and a parametric domain. It has been considered as a fundamental tool for digital geometric processing which is required to develop several applications of digital geometries. In this paper, we propose a novel 3D volume parameterization by means that a harmonic mapping is established between a 3D volume mesh and a unit solid cube. To do that, we firstly partition the boundary of the given 3D volume mesh into the six different rectangular patches whose adjacencies are topologically identical to those of a surface cube. Based on the partitioning result, we compute the boundary condition as a precondition for computing a volume mesh parameterization. Finally, the volume mesh parameterization with a low-distortion can be accomplished by performing a harmonic mapping, which minimizes the harmonic energy, with satisfying the boundary condition. Experimental results show that our method is efficient enough to compute 3D volume mesh parameterization for several models, each of whose topology is identical to a solid sphere.