In this article we propose and numerically implement a mathematical model for the simulation of three-dimensional semiconductor devices characterized by an heterogeneous material structure. The model consists of a system of nonlinearly coupled timedependent diffusion-reaction partial differential equations with convection terms describing the principal electrical, thermal and chemical phenomena that determine the macroscopic electrical response of the device under the action of externally applied electrical and thermal forces. The system is supplied with suitable initial, boundary and interface conditions that account for the interaction occurring among the various regions of the device with the surrounding environment. Temporal semi-discretization of the problem is carried out with the Backward Euler Method while a fixed-point iteration of Gummel type is used for system decoupling. Numerical approximation of the linearized subproblems is carried out using an exponentially fitted stabilized Finite Element Method on unstructured tetrahedral grids. Several computational experiments are included to validate the physical accuracy of the proposed computational algorithm in the study of realistic device structures.