  Implicit constitutive theory provides a very general framework for fluid flow
models, including both Newtonian and generalized Newtonian fluids, where the
Cauchy stress tensor and the rate of strain tensor are assumed to be related by
an implicit relation associated with a maximal monotone graph. For
incompressible unsteady flows of such fluids, subject to a homogeneous
Dirichlet boundary condition on a Lipschitz polytopal domain $\Omega \subset
\mathbb{R}^d$, $d \in \{2,3\}$, we investigate a fully-discrete approximation
scheme, using a spatial mixed finite element approximation combined with
backward Euler time-stepping. We show convergence of a subsequence of
approximate solutions, when the velocity field belongs to the space of
solenoidal functions contained in $L^\infty(0,T;L^2(\Omega)^d)\cap
L^q(0,T;W^{1,q}_0(\Omega)^d)$, provided that $q\in
\big(\frac{2d}{d+2},\infty\big)$, which is the maximal range for $q$ with
respect to existence of weak solutions. This is achieved by a technique based
on splitting and regularizing, the use of a solenoidal parabolic Lipschitz
truncation method, a local Minty-type monotonicity result, and various weak
compactness results.
