We consider the discretization method for solving three-dimensional variable-order (3D-VO) time-fractional partial differential equations. The proposed method is developed based on discrete shifted Hahn polynomials (DSHPs) and their operational matrices. In the process of method implementation, the modified operational matrix (MOM) and complement vector (CV) of integration and pseudooperational matrix (POM) of VO fractional derivative plays an important role in the accuracy of the method. Further, we discuss the error of the approximate solution. At last, the methodology is validated by well test examples in two types of space domains. In order to evaluate the accuracy and applicability of the approach, the results are compared with other methods.