  We consider a system of nonlinear partial differential equations modelling
the steady motion of an incompressible non-Newtonian fluid, which is chemically
reacting. The governing system consists of a steady convection-diffusion
equation for the concentration and the generalized steady Navier-Stokes
equations, where the viscosity coefficient is a power-law type function of the
shear-rate, and the coupling between the equations results from the
concentration-dependence of the power-law index. This system of nonlinear
partial differential equations arises in mathematical models of the synovial
fluid found in the cavities of moving joints. We construct a finite element
approximation of the model and perform the mathematical analysis of the
numerical method in the case of two space dimensions. Key technical tools
include discrete counterparts of the Bogovski\u{\i} operator, De Giorgi's
regularity theorem in two dimensions, and the Acerbi-Fusco Lipschitz truncation
of Sobolev functions, in function spaces with variable integrability exponents.
