We consider unsteady flows of incompressible fluids with a general implicit constitutive equation relating the deviatoric part of the Cauchy stress S and the symmetric part of the velocity gradient D in such a way that it leads to a maximal monotone (possibly multivalued) graph and the rate of dissipation is characterized by the sum of a Young function depending on D and its conjugate being a function of S. Such a framework is very robust and includes, among others, classical power-law fluids, stress power-law fluids, fluids with activation criteria of Bingham or Herschel-Bulkley type, and shear rate-dependent fluids with discontinuous viscosities as special cases. The appearance of S and D in all the assumptions characterizing the implicit relationship G(D, S) = 0 is fully symmetric. We establish long-time and large-data existence of weak solution to such a system completed by the initial and the Navier slip boundary conditions in both the subcritical and supercritical cases. We use tools such as Orlicz functions, properties of spatially dependent maximal monotone operators, and Lipschitz approximations of Bochner functions taking values in Orlicz-Sobolev spaces.