A new approach to constitutive modeling for non-Newtonian liquids is presented, with a particular application to a dilute suspension of spheroids (prolate and oblate) subject to Brownian rotations at finite Peclet numbers, but with a general framework that can be applied to other complex fluids, including concentrated suspensions and emulsions. A generalized traceless Oldroyd model is used for the particle contribution to the stress, with five material parameters as functions of one instantaneous flow invariant: the intrinsic energy dissipation rate. All five parameters are found from simultaneously fitting the model to numerical results for two base flows at arbitrary flow intensities: planar extensional flow and simple shear flow. Precise numerical solutions of the Fokker-Planck-Smoluchowski equation for the orientation distribution function are used to prepare the database rheological functions, and also to validate the resulting model in flows with arbitrary kinematics. The present work also verifies the accuracy of various closure models in the literature for dilute suspensions of spheroids in several test flows. Although these closure models give excellent results and are often more accurate than the present approach in the test cases, they are relevant specifically to fiber suspensions, and it is not clear how to generalize the closure approach to other types of non-Newtonian, microstructurally complex liquids. In contrast, the present approach to constitutive modeling is more general and does not hinge on a particular microstructure; it can be applied to various systems with strong hydrodynamical interactions (highly concentrated suspensions, emulsions, etc.). (C) 2016 The Society of Rheology.