Abstract Multidimensional numerical models are useful tools for understanding the heat transfer mechanisms and performance optimization of thermoelectric generators (TEGs). In this study, two three-dimensional numerical models are developed for TEGs based on different formulations, but with similar abilities for heat and electricity transfer analysis and performance prediction. Model 1 solves the conservation equations of the Seebeck potential and the Ohmic potential separately, and the total built-in potential can be obtained based on the solved Seebeck and Ohmic potentials. Model 2 solves the conservation equation of the total built-in potential directly, and the conservation equation for the Ohmic potential is also solved. The comparison between Model 1 and Model 2 shows that Model 2 is slightly more precise for power output prediction. The detailed formulations of these two models are described, and the difference among the present and previous models is also discussed. Some important modeling aspects are elucidated for the TEG models, such as the conservation equations and boundary conditions. Parametric studies are carried out based on various thermal boundary conditions. The influence of the TEG semiconductor shape on performance is investigated in details. It is found that for the nearly same volume of semiconductor materials, changing the shape from normal cuboid (constant cross-sectional area) to hexahedrons (variable cross-sectional area) could increase the power output significantly. The reason is that the temperature gradient could be enhanced when proper variable cross-sectional areas are used.