The hydrodynamic (HD) model of semiconductor devices is solved numerically in three-dimensions (3-D) for the MOSFET device. The numerical instabilities of the HD model are analyzed to develop a stable discretization. The formulation is stabilized by using a new, higher-order discretization for the relaxation-time approximation (RTA) term of the energy-balance (EB) equation. The developed formulation is used to model the MOSFET.