Quantitative analysis of flows of fiber suspensions in viscoelastic matrices, which is the general situation for thermoplastic composites, requires constitutive equations which incorporate specific features of the system and its constituents. Matrix viscoelasticity, fiber orientation and fiber/matrix interactions are key parameters to model such systems. In this work, the constituents of the system are represented by two second order symmetric tensors: c(r, t) for the viscoelastic matrix and a(r, t) for the fiber orientation. The time evolution equation for c(r, t) is developed in the generalized Poisson bracket framework with a finitely extensible non-linear elastic (FENE-P) and a Hookean Helmholtz energy functions. Several expressions for the mobility tensor including expressions with fiber matrix interactions are used. The time evolution equation for ar(r, t) is based on the classical Jeffery equation modified to include fiber/fiber interactions in the case of semi-dilute suspensions. The sensitivity of the model to the choice of the mobility tensor together with the effect of fiber volume fraction on the prediction of material Functions in start up and steady shear flows are discussed.