\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2017 (2017), No. 155, pp. 1--28.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2017 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2017/155\hfil Non-Newtonian fluids] {On non-Newtonian fluids with convective effects} \author[S. Herr\'on, E. J. Villamizar-Roa \hfil EJDE-2017/155\hfilneg] {Sigifredo Herr\'on, \'Elder J. Villamizar-Roa } \address{Sigifredo Herr\'on \newline Universidad Nacional de Colombia - Sede Medell\'in, Escuela de Matem\'aticas, A.A. 3840, Medell\'in, Colombia} \email{sherron@unal.edu.co} \address{\'Elder J. Villamizar-Roa (corresponding author)\newline Universidad Industrial de Santander, Escuela de Matem\'aticas, A.A. 678, Bucaramanga, Colombia} \email{jvillami@uis.edu.co} \thanks{Submitted January 12, 2017. Published June 28, 2017.} \subjclass[2010]{35Q35, 76D03, 76D05, 35D30, 35D35} \keywords{Non-Newtonian fluids; shear-dependent viscosity; weak solutions; \hfill\break\indent strong solutions; uniqueness} \begin{abstract} We study a system of partial differential equations describing a steady thermoconvective flow of a non-Newtonian fluid. We assume that the stress tensor and the heat flux depend on temperature and satisfy the conditions of $p,q$-coercivity with $p>\frac{2n}{n+2}$, $q>\frac{np}{p(n+1)-n}$, respectively. Considering Dirichlet boundary conditions for the velocity and a mixed and nonlinear boundary condition for the temperature, we prove the existence of weak solutions. We also analyze the existence and uniqueness of strong solutions for small and suitably regular data. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction} This article analyzes a system of partial differential equations describing a steady thermoconvective flow of a non-Newtonian fluid in a bounded domain $\Omega$ of $\mathbb{R}^n$, $n=2,3$, with smooth enough boundary $\partial\Omega$. The model is given by the system of PDEs \begin{equation}\label{intro1} \begin{gathered} -\operatorname{div}\big(\mu(\cdot,\theta)\mathbf{T}(D(\mathbf{u}))\big) +\operatorname{div}(\mathbf{u}\otimes\mathbf{u})+\nabla\pi =\theta\mathbf{f}\quad \text{in } \Omega,\\ \operatorname{div}\ \mathbf{u}=0 \quad \text{in } \Omega,\\ -\operatorname{div}(\kappa(\cdot,\theta)\mathbf{a}(\nabla \theta)) +\mathbf{u}\cdot\nabla\theta =g \quad \text{in } \Omega, \end{gathered} \end{equation} where the unknowns are $\mathbf{u}:\Omega\to\mathbb{R}^n$, $\theta:\Omega\to\mathbb{R}$ and $\pi:\Omega\to\mathbb{R}$ denoting the velocity, the temperature and the pressure of the fluid, respectively. The field $\mathbf{f}$ denotes the given external body forces and $g$ represents the heat source. The symbol $\mathbf{T}:\mathbb{M}^{n\times n}_{\rm sym}\to \mathbb{M}^{n\times n}_{\rm sym}$ denotes the extra stress tensor and $\mathbf{a}$ indicates the constitutive law for diffusivity. The symbol $D(\mathbf{u})$ represents the symmetric part of the velocity gradient $\nabla\mathbf{u}$, that is, $D(\mathbf{u})=\frac12(\nabla\mathbf{u}+\nabla^T\mathbf{u})$; the functions $\mu(\cdot,\theta)>0$, $\kappa(\cdot,\theta)>0$ denote the kinematic viscosity and thermal conductivity, respectively. Equations \eqref{intro1}$_1$ and \eqref{intro1}$_3$ correspond to the momentum and heat equations respectively; the second equation in \eqref{intro1} corresponds to the incompressibility condition. We assume that the functions $\boldsymbol{\eta}\to \mathbf{T}(\boldsymbol{\eta})$ and $\boldsymbol{\chi}\to \mathbf{a}(\boldsymbol{\chi})$ are continuous in $\mathbb{M}^{n\times n}_{\rm sym}$ and $\mathbb{R}^n$ respectively, and satisfy the following conditions for some $p,q>1$ (see notation in Section 2): \begin{itemize} \item[i)] (Coercivity) There exist $\tau_1,\alpha_1>0$ such that \begin{equation}\label{coercivity} \begin{gathered} \mathbf{T}(\boldsymbol{\eta}):\boldsymbol{\eta} \geq \tau_1|\mathbf{\boldsymbol{\eta}}|^p,\\ \mathbf{a}(\boldsymbol{\chi})\cdot\boldsymbol{\chi} \geq \alpha_1|\boldsymbol{\chi}|^q, \end{gathered} \end{equation} for all $\boldsymbol{\eta}\in \mathbb{M}^{n\times n}_{\rm sym}$, $\boldsymbol{\chi}\in\mathbb{R}^n$. \item[(ii)] (Polynomial growth) There exist $\tau_2,\alpha_2>0$ such that \begin{equation}\label{polinomial-growth} \begin{gathered} |\mathbf{T}(\boldsymbol{\eta})|\leq \tau_2(1+|\boldsymbol{\eta}|)^{p-1},\\ |\mathbf{a}(\boldsymbol{\chi})|\leq \alpha_2|\boldsymbol{\chi}|^{q-1}, \end{gathered} \end{equation} for all $\boldsymbol{\eta}\in \mathbb{M}^{n\times n}_{\rm sym}$, $\boldsymbol{\chi}\in\mathbb{R}^n$. \item[(iii)] (Strict monotonicity) \begin{equation}\label{monotonicity} \begin{gathered} \left(\mathbf{T}(\boldsymbol{\eta})-\mathbf{T}(\boldsymbol{\xi})\right): (\boldsymbol{\eta}-\boldsymbol{\xi})>0,\quad \forall \boldsymbol{\eta}, \boldsymbol{\xi}\in \mathbb{M}^{n\times n}_{\rm sym}, \boldsymbol{\eta}\neq\boldsymbol{\xi},\\ \left(\mathbf{a}(\boldsymbol{\varsigma})-\mathbf{a}(\boldsymbol{\chi})\right) \cdot (\boldsymbol{\varsigma}-\boldsymbol{\chi})>0,\quad \forall \boldsymbol{\varsigma},\boldsymbol{\chi}\in \mathbb{R}^n,\; \boldsymbol{\varsigma} \neq\mathbf{x}. \end{gathered} \end{equation} \end{itemize} The general non-linear tensor function $\mathbf{T}$ and constitutive law for the heat flux $\mathbf{a}$ allow to consider a large class of non-Newtonian fluids subjected to heat effects, which have physical motivations as described in \cite{Consiglieri,Rajagopal,Rodrigues} and references therein. Typical prototypes of extra stress tensors used in applications are $\mathbf{T}_1(\boldsymbol{\eta})=2\mu(1+|\boldsymbol{\eta}|^2)^{(p-2)/2} \boldsymbol{\eta}$ and $\mathbf{T}_2(\boldsymbol{\eta})=2\mu(1+|\boldsymbol{\eta}|)^{p-2}\boldsymbol{\eta}$ with $p>1$. In these cases, if $p=2$ and $\mathbf{a}$ is the identity, we obtain the classical Boussinesq equation (see \cite{Boldfrini-Lorca,Ferreira0,Ferreira1,Ferreira2,Villamizar1,Villamizar2}). We also consider the following hypotheses on the viscosity and the thermal conductivity functions $\mu,\kappa$. It is assumed that $\mu,\kappa:\Omega\times \mathbb{R}\to\mathbb{R}$ are Carath\'eodory functions (i.e., for each fixed $\theta$ the functions $x\mapsto \mu(x,\theta), \ x\mapsto \kappa(x,\theta)$ are (Lebesgue) measurable in $\Omega$ and, the functions $\theta\mapsto \mu(x,\theta)$, $\theta\mapsto \kappa(x,\theta)$ are continuous for almost every $x\in\Omega$) such that \begin{equation}\label{kap} \begin{gathered} 0<\mu_1\leq \mu(x,\theta)\leq \mu_2 \quad\text{a.e. } x\in\Omega, \; \forall \theta\in\mathbb{R},\\ 0<\kappa_1\leq \kappa(x,\theta)\leq \kappa_2 \quad\text{a.e. } x\in\Omega, \; \forall \theta\in\mathbb{R}. \end{gathered} \end{equation} System \eqref{intro1} is complemented with the mixed boundary conditions \begin{equation}\label{intro2} \begin{gathered} \mathbf{u}=0 \quad \text{on} \ \partial\Omega,\\ \theta=0 \text{ on } \Gamma_0,\quad \kappa(\cdot,\theta)\mathbf{a}(\nabla\theta)\cdot\mathbf{n}+\gamma\theta=h \text{ on } \Gamma:=\partial\Omega\setminus\overline{\Gamma_0}, \end{gathered} \end{equation} where $\gamma$ is a non-negative constant, $\mathbf{n}$ denotes the unit outward normal on the boundary $\partial\Omega$, and $\Gamma_0$ is a open subset of $\partial\Omega$. Boundary conditions \eqref{intro2}$_2$ include several physical boundary conditions like those appearing in several natural convection problems \cite{Ferreira1,Villamizar1}. The existence of weak solutions in the case of Navier-Stokes equations for flows with shear-dependent viscosity is known in $W^{1,p}(\Omega)$ for $p\geq 2n/(n+2)$. For the case $p\geq 3n/(n+2)$, the existence of weak solutions was obtained by Lions \cite{Lions} and Ladyzhenskaya \cite{Ladyz} by using monotone operators theory. In \cite {Ruzicka}, using the $L^\infty$-truncation method, the authors obtained the existence of weak solutions for $p\geq 2n/(n+1)$. This method is based on the construction of a special class of test functions, and a characterization of the pressure, which permit the almost everywhere convergence of $D(\mathbf{u}^m)$ to $D(\mathbf{u})$, where $\mathbf{u}^m$ corresponds to a sequence of approximated solutions $\mathbf{u}^m$ of the original problem. However, this method only works for $p\geq 2n/(n+1)$ because of the required $L^1$-integrability of the nonlinear term $(\mathbf{u}\cdot\nabla)\mathbf{u}$. To consider the case $p\geq 2n/(n+2)$, in \cite{Frehse} the Lipschitz truncation method was applied, which permits controling the nonlinear term $(\mathbf{u}\cdot\nabla)\mathbf{u}$ using a test function class smoother than the test functions used in the $L^\infty$-truncation method. On the other hand, focusing on the boundary-value problem \eqref{intro1}-\eqref{intro2}, the existence of weak solutions for $p> 2n/(n+1)$ and $q>np/(p(n+1)-n)$ was obtained in \cite{Consiglieri}. Motivated by this facts, in the first part of this paper, we extend the results of \cite{Consiglieri} to the case $p> 2n/(n+2)$ and $q>np/(p(n+1)-n)$. The second part of this article concentrates on the existence of regular solutions to the boundary value problem \eqref{intro1}-\eqref{intro2}. In the case of Navier-Stokes equations for flows with shear-dependent viscosity, there are few works concerning the regularity of weak solutions (cf. \cite{Arada,daVeiga,Crispo,Kaplicky} and some references therein). The most recent results for the steady Navier-Stokes equations for flows with shear-dependent viscosity are due to Arada \cite{Arada}. In \cite{Arada}, the author assumed that $\mathbf{T}$ is a classical power law stress tensor of the form $\mathbf{T}(\boldsymbol{\eta})=\mathbf{T}_1(\boldsymbol{\eta}) :=2\mu(1+| \boldsymbol{\eta}|^2)^{\frac{p-2}{2}}\boldsymbol{\eta}$ or $\mathbf{T}(\boldsymbol{\eta})=\mathbf{T}_2(\boldsymbol{\eta}) :=2\mu(1+| \boldsymbol{\eta}|)^{{p-2}}\boldsymbol{\eta}$, where $\mu>0$ is a viscosity coefficient and $p>1$. He proved the existence of strong solutions $\mathbf{u}\in \mathbf{W}^{2,q}(\Omega)$, $q>n$, by assuming that $\| \mathbf{f}\|_q/\mu$ is small enough. Some uniqueness results were also established. However, to the best of our knowledge, there are no results of existence of strong solutions for the steady problem \eqref{intro1}-\eqref{intro2}. In the second part of this paper, we will study the existence of a strong solution for small and suitably regular data by taking $\mathbf{T}=\mathbf{T}_1$ or $\mathbf{T}=\mathbf{T}_2$. To ease the exposition, we also simplify the boundary conditions on temperature $\theta$; however, a similar analysis can be adapted for other types of boundary data. Our approach is based on regularity results for the Stokes problem and the Laplace equation, as well as a fixed-point argument. Observe that $\mathbf{T}_1$ depends on the differentiable term $| D(\mathbf{u})|^2$ while $\mathbf{T}_2$ depends merely on the Lipschitz continuous term $| D(\mathbf{u})|$; thus, in the case $\mathbf{T}=\mathbf{T}_1$ we can use the classical regularity results for the Stokes system to solve the velocity equation for a fixed temperature. However, in the case $\mathbf{T}=\mathbf{T}_2$, to overcome the difficulty caused by the lack of regularity of $\mathbf{T}_2$, we first introduce a family of penalized problems, then, we establish the existence of penalized strong solutions and finally, we carry out the pass to the limit in the sequence of penalized problems, as the penalization term goes to zero. This article is organized as follows. In Section 2, we introduce the notation. Section 3 is devoted to the existence of weak solutions. In Section 4, we analyze the existence of strong solutions in both cases: with the differentiable stress tensor $\mathbf{T}_1$, and with the Lipschitz continuous stress tensor $\mathbf{T}_2$. In Section 4, we also give conditions on the data which ensure that the obtained strong solution agrees with weak solutions. \section{Notation} In this section, we establish some general notation to be used throughout this article. As usual, $C_0^\infty(\Omega)$ denotes the set of all $C^\infty$-functions with compact support in $\Omega$, while $\mathbf {C}_{0,\sigma}^\infty(\Omega)$ consists of functions $\boldsymbol{\Phi}\in C_0^\infty(\Omega)$ such that $\operatorname{div}\boldsymbol{\Phi}=0$. For $p,q>1$ we set \begin{gather*} \mathbf{H}_q:=\overline{\mathbf{C}_{0,\sigma}^\infty(\Omega)}^{\|\cdot\|_q} =\{\mathbf{u}\in \mathbf{L}^q(\Omega): \operatorname{div}\mathbf{u}=0, \; \mathbf{u}\cdot\mathbf{n}=0\ \text{on } \partial\Omega\},\\ \mathbf{V}_p:=\overline{\mathbf{C}_{0,\sigma}^\infty(\Omega)}^{\|\nabla \cdot\|_p}=\{\mathbf{u}\in \mathbf{W}_0^{1,p}(\Omega): \operatorname{div}\mathbf{u}=0\},\\ X_{q}:=\{\theta\in W^{1,q}(\Omega)\cap L^{2}(\Gamma): \theta=0 \text{ on } \Gamma_0\}. \end{gather*} Here $\mathbf{V}_p$ and $X_{q}$ are Banach spaces with the norms $\| D(\mathbf{u})\|_p$ and $\|\theta\|_{X_{q}}=\|\nabla\theta\|_q +\|\theta\|_{2,\Gamma}$. As usual, $\| \cdot\| _p$ denotes the $L^p$-norm. Notice that, due to the trace theorem, $X_{q}=\{\theta\in W^{1,q}(\Omega): \theta=0 \text{ on } \Gamma_0\}$ if $q\geq2n/(n+1)$. For $x, y\in \mathbb{R}$ we denote $(x, y)^+=\max\{x, y\}$, $x^+=\max \{x, 0\}, \ S_p=(|p-2|, 2)^+$. Frequently, we will use the notation $\langle \cdot,\cdot\rangle_{X'}$ (or simply $\langle \cdot,\cdot\rangle$ if there is no ambiguity) to represent the duality product between $X'$ and $X$, for the Banach space $X$. We also introduce the constants \[ 2r_p=1+(p-3)^+-(p-4)^+,\quad \gamma_p=\frac{[(p,3)^+-2]^{(p,3)^{+}-2}}{[(p,3)^+-1]^{(p,3)^{+}-1}}\,. \] For $m\in\mathbb{N}$ and $1n$ and $\delta>0$, let us denote by $B_\delta$ the convex set defined by \begin{equation}\label{ball} B_\delta=\{[\boldsymbol{\xi},\omega]\in \mathbf{V}_{2,q} \times W^{2,r}(\Omega): C_E\| \nabla\boldsymbol{\xi}\| _{{1,q}} \leq \delta, \; C_{\tilde E}\| \nabla\omega\| _{{1,r}}\leq \delta\}, \end{equation} where $C_E$ is the norm of the embedding of $W^{1,q}(\Omega)$ into $L^{\infty}(\Omega)$ and $C_{\tilde E}$ is the norm of the embedding of $W^{1,r}(\Omega)$ into $L^{\infty}(\Omega)$. Also, we consider the space $\mathbf{V}_{2,q}\times( W^{2,r}(\Omega)\cap W^{1,r}_0(\Omega))$ endowed with the norm \[ \| [\boldsymbol{\xi},\omega]\|_{1,q,r}:=\max\{\| \nabla\boldsymbol{\xi}\| _{1,q}, \| \nabla\omega\| _{1,r}\}. \] Throughout the paper, $\mathbb{M}^{n\times n}$ denotes the space of all real $n\times n$ matrices and $\mathbb{M}^{n\times n}_{\rm sym}$ its subspace of all symmetric $n\times n$ matrices. We use the following summation convention on repeated indices: $\boldsymbol{\eta}:\boldsymbol{\xi}:=\boldsymbol{\eta}_{ij}\boldsymbol{\xi}_{ij}$ for $\boldsymbol{\eta}:\boldsymbol{\xi}\in \mathbb{M}^{n\times n}$, $(\mathbf{u}\otimes\mathbf{v})_{ij}:= u^iv^j$ for $\mathbf{u},\mathbf{v}\in\mathbb{R}^n$ and $\mathbf{u}\cdot\mathbf{v}:= u^iv^i$. Also we set $| \mathbf{u}|:= (\mathbf{u}\cdot \mathbf{u})^{1/2}$ and $| \boldsymbol{\eta}|:= (\boldsymbol{\eta}:\boldsymbol{\eta})^{1/2}$ for $\mathbf{u}\in \mathbb{R}^n$, $\boldsymbol{\eta}\in \mathbb{M}^{n\times n}$. Finally, the letter $C$ stands for several positive constants that may change line by line; also $C_P=C_P(n,s,\Omega)$ denotes the Poincar\'e constant corresponding to the general Poincar\'e inequality $\|\cdot\|_s\leq C_P\|\nabla(\cdot)\|_s$. \section{Weak solutions} The aim of this section is to prove the existence of weak solutions to problem \eqref{intro1}-\eqref{intro2} for the case $\frac{2n}{n+2}\frac{np}{p(n+1)-n}$. The existence of weak solutions for $p>\frac{2n}{n+1}$ was analyzed in \cite{Consiglieri}. We assume that $\mathbf{f}\in \mathbf{L}^\infty(\Omega), g\in (W^{1,q}(\Omega))', h\in L^{2}(\Gamma)$. First we establish the notion of weak solution to \eqref{intro1}-\eqref{intro2}. \begin{definition} \rm We say that a pair $[\mathbf{u}, \theta]\in \mathbf{V}_p\times X_{q}$ is a weak solution to problem \eqref{intro1}-\eqref{intro2} if \begin{gather*} \int_\Omega \mu(\cdot,\theta)\mathbf{T}\left(D(\mathbf{u})\right): D(\mathbf{\Phi})\,dx -\int_\Omega (\mathbf{u}\otimes\mathbf{u}):D(\mathbf{\Phi})\,dx = \int_\Omega\theta\mathbf{f}\cdot\mathbf{\Phi}\,dx,\\ \forall\mathbf{\Phi}\in{\mathbf{C}_{0,\sigma}^\infty(\Omega)}, \\ \int_\Omega\kappa(\cdot,\theta)\mathbf{a}(\nabla\theta)\cdot\nabla\phi\,dx +\int_\Omega\phi\mathbf{u}\cdot\nabla\theta\,dx +\gamma\int_\Gamma\theta\phi\,d\Gamma =\langle g,\phi\rangle_{(W^{1,q}(\Omega))'} +\int_\Gamma h\phi\,d\Gamma,\\\ \forall\phi\in {C_{0}^\infty(\Omega)}. \end{gather*} \end{definition} The purpose of this section is to prove the following theorem on existence of weak solutions. \begin{theorem}\label{main-theorem} Let $p>\frac{2n}{n+2}$, $q>\frac{np}{p(n+1)-n}$, $\mathbf{f}\in \mathbf{L}^\infty(\Omega)$, $g\in (W^{1,q}(\Omega))', h\in L^{2}(\Gamma)$. There exists a weak solution $[\mathbf{u}, \theta]\in \mathbf{V}_p\times X_{q}$ to problem \eqref{intro1}-\eqref{intro2}. \end{theorem} To prove Theorem \ref{main-theorem}, we first consider a suitable sequence of approximate problems (see \eqref{penal1}-\eqref{intro2aprox} below); we establish the existence of approximate solutions, as well as some a priori estimates. In a second step, we describe the passing to the limit of the sequence of approximate solutions. Finally, we analyze the almost everywhere convergence in $\Omega$ of $[D(\mathbf{u}^m),\nabla \theta^m]\to [D(\mathbf{u}),\nabla \theta]$ through the Lipschitz truncation method. \subsection{Approximate solutions} For $m\in\mathbb{N}$ and $t>\max\{\frac{2p}{p-1}, \frac{nq}{q(n+1)-2n}\}$, we define the approximated problem: Find a weak solution $[\mathbf{u}^{m}, \theta^{m}]$ of the system \begin{equation}\label{penal1} \begin{gathered} -\operatorname{div}\left(\mu(\cdot,\theta^m)\mathbf{T}(D(\mathbf{u}^m))\right) +\operatorname{div}(\mathbf{u}^m\otimes\mathbf{u}^m)+\frac1{m} | \mathbf{u}^m|^{t-2}\mathbf{u}^m+\nabla\pi=\theta^m\mathbf{f}\quad \text{in } \Omega,\\ \operatorname{div}\ \mathbf{u}^m =0 \quad \text{in } \Omega,\\ -\operatorname{div}(\kappa(\cdot,\theta^m)\mathbf{a}(\nabla \theta^m)) +\mathbf{u}^m\cdot\nabla\theta^m =g \quad \text{in } \Omega, \end{gathered} \end{equation} with the boundary conditions \begin{equation}\label{intro2aprox} \begin{gathered} \mathbf{u}^m =0 \quad \text{on } \partial\Omega,\\ \theta^m =0 \text{ on } \Gamma_0,\quad \kappa(\cdot,\theta^m)\mathbf{a}(\nabla\theta^m)\cdot\mathbf{n} +\gamma\theta^m=h \text{ on } \Gamma:=\partial\Omega\setminus\overline{\Gamma_0}. \end{gathered} \end{equation} Following the ideas presented in \cite{Frehse,Frehse97}, we obtain the existence of a weak solution $[\mathbf{u}, \theta]$ of \eqref{intro1}-\eqref{intro2} as the limit of a sequence of weak solutions $[\mathbf{u}^m, \theta^m]$ of \eqref{penal1}-\eqref{intro2aprox}. A weak solution of the system \eqref{penal1}-\eqref{intro2aprox} is a pair $[\mathbf{u}^m,\theta^m]\in \mathbf{V}_p\times X_{q}$ satisfying \begin{gather} \label{prob-aprox} \begin{aligned} &\int_\Omega \mu(\cdot,\theta^m)\mathbf{T}\left(D(\mathbf{u}^m)\right): D(\mathbf{\Phi})\,dx -\int_\Omega (\mathbf{u}^m\otimes\mathbf{u}^m):D(\mathbf{\Phi})\,dx\\ &+\frac1{m}\int_\Omega | \mathbf{u}^m|^{t-2}\mathbf{u}^m\cdot\mathbf{\Phi}\,dx\\ &= \int_\Omega\theta^m\mathbf{f}\cdot\mathbf{\Phi}\,dx, \end{aligned}\\ \label{prob-theta} \begin{aligned} &\int_\Omega \kappa(\cdot,\theta^m)\mathbf{a}\left(\nabla\theta^m\right) \cdot\nabla\phi\,dx -\int_\Omega \theta^m\mathbf{u}^m\cdot\nabla\phi\,dx +\gamma\int_\Gamma\theta^m\phi\,d\Gamma \\ &=\langle g,\phi\rangle_{(W^{1,q}(\Omega))'} +\int_\Gamma h\phi\,d\Gamma, \end{aligned} \end{gather} for all $\mathbf{\Phi}\in \mathbf{C}_{0,\sigma}^\infty(\Omega), \phi\in C_{0}^\infty(\Omega)$. The following lemma provides the existence of a weak solution to \eqref{penal1}-\eqref{intro2aprox}. \begin{lemma}\label{lemma-estimate-aprox} Let $p>2n/(n+2)$, $t\geq 2p'$, $q>\frac{np}{p(n+1)-n}$. Assume that $\mathbf{f}\in \mathbf{L}^\infty(\Omega)$, $g\in (W^{1,q}(\Omega))'$, $h\in L^{2}(\Gamma)$. Then, there exists a unique weak solution $[\mathbf{u}^m,\theta^m]\in (\mathbf{V}_p\cap \mathbf{H}_t)\times X_{q} $ of \eqref{prob-aprox}-\eqref{prob-theta}. Moreover, the following uniform estimates hold \begin{gather} \frac{\tau_1\mu_1}{2}\| \mathbf{u}^m\|_{1,p}^p +\frac{1}{m}\| \mathbf{u}^m\|_t^t \leq C_1\|\mathbf{f}\|^{p'}_\infty\left(\| g\|_{(W^{1,q})'}^{q'} +\| h\|_{2,\Gamma}^{2}\right)^{p'/q},\label{estimate-aprox} \\ \frac{\alpha_1\kappa_1}{2}\|\nabla \theta ^m\|_q^q +\frac{\gamma}{2}\|\theta ^m\|_{2,\Gamma}^{2} \leq C_2\left(\| g\|_{(W^{1,q})'}^{q'} +\| h\|_{2,\Gamma}^{2}\right), \label{estimate-aprox2} \end{gather} for some constants $C_1,C_2>0$ independent on $m$. \end{lemma} \begin{proof} The proof follows by standard arguments of the monotone operator theory (cf. \cite{Frehse,Frehse97}). The uniform estimates \eqref{estimate-aprox}-\eqref{estimate-aprox2} follow by taking $\mathbf{\Phi}=\mathbf{u}^m$ and $\phi=\theta^m$ in \eqref{prob-aprox} and \eqref{prob-theta}, respectively, and using the assumptions on $\mathbf{T},\mathbf{a},\mathbf{f},g$ and $h$. \end{proof} \subsection{Existence of weak solutions} The existence of a weak solution to the problem \eqref{intro1}-\eqref{intro2} will be obtained as the limit, as $m$ goes to infinity, in the sequence of solutions $[\mathbf{u}^m,\theta^m]$ of \eqref{prob-aprox}-\eqref{prob-theta}. We use the Lipschitz truncation method used previously in \cite{Frehse} in the context of incompressible fluids with shear-dependent viscosity (without heat effects). Following \cite{Frehse}, we introduce the sequence of approximate pressures $\pi^m$, observing that in \eqref{prob-aprox} we can consider test functions $\Phi$ from $\mathbf{V}_p\cap \mathbf{V}_r=\mathbf{V}_r$ with $r=np/[(n+2)p-2n]$. Notice that for this value of $r$ and $2n/(n+2)2n/(n+2)$, and writing $\mathbf{u}^m=(\mathbf{u}^m-\mathbf{u})+\mathbf{u}$, for $\mathbf{\Phi}\in{\mathbf{C}_{0}^\infty(\Omega)}$ and letting $m\to\infty$ we obtain \[ \int_\Omega (\mathbf{u}^m\otimes\mathbf{u}^m):D(\mathbf{\Phi})\,dx \to \int_\Omega (\mathbf{u}\otimes\mathbf{u}):D(\mathbf{\Phi})\,dx. \] Also, since $\theta^m\to \theta$ in $L^1(\Omega)$ and a.e. in $\Omega$, and $\mu$ is a Carath\'eodory function, then $\mu(\cdot,\theta^m)\to \mu(\cdot,\theta)$ a.e. in $\Omega$. Then, collecting the last convergences, we have \begin{equation} \int_\Omega \mu(\cdot,\theta)\boldsymbol{\chi}: D(\mathbf{\Phi})\,dx +\int_\Omega (\mathbf{u}\otimes\mathbf{\Phi}): D(\mathbf{u})\,dx= \int_\Omega\theta\mathbf{f}\cdot\mathbf{\Phi}\,dx +\int_\Omega \pi\operatorname{div}\mathbf{\Phi}\,dx, \label{formulation-pressure} \end{equation} for all $\mathbf{\Phi}\in \mathbf{C}^\infty_0(\Omega)$ and consequently for all $\mathbf{\Phi}\in \mathbf{W}^{1,r}_0(\Omega)$. As before, since $\theta^m\to \theta$ in $L^1(\Omega)$ and a.e. in $\Omega$, and $\kappa$ is a Carath\'eodory function, then $\kappa(\cdot,\theta^m)\to \kappa(\cdot,\theta)$ a.e. in $\Omega$. Then, from the uniform estimates \eqref{est20}, \eqref{rev1} and \eqref{rev2} we also get \begin{equation} \int_\Omega \Big(\kappa(\cdot,\theta)\boldsymbol{\chi}_1 -\theta\mathbf{u}\Big)\cdot\nabla\phi\,dx +\gamma\int_\Gamma\theta\phi\,d\Gamma=\langle g,\phi\rangle_{(W^{1,q}(\Omega))'} +\int_\Gamma h\phi\,d\Gamma,\label{lim21} \end{equation} for all $\phi\in C^\infty_0(\Omega)$ and consequently for all $\phi\in X_{q}$. It remains to prove that $\boldsymbol{\chi}=\mathbf{T}\left(D(\mathbf{u}^m)\right)$ and $\boldsymbol{\chi}_1=\mathbf{a}\left(\nabla\theta^m\right)$. To this end, it is sufficient to prove that \begin{equation} [D(\mathbf{u}^m),\nabla \theta^m]\to [D(\mathbf{u}),\nabla\theta]\quad \text{in measure on } \Omega,\label{convergence-aprox} \end{equation} or almost everywhere convergence on compact subsets of $\Omega$. Having proved \eqref{convergence-aprox}, through a diagonal procedure, we can find a subsequence of $([\mathbf{u}^m,\theta^m])_{m\in\mathbb{N}}$, still denoted by $([\mathbf{u}^m,\theta^m])_{m\in\mathbb{N}}$, such that \begin{equation} [D(\mathbf{u}^m),\nabla \theta^m]\to [D(\mathbf{u}),\nabla \theta] \quad \text{almost everywhere in } \Omega.\label{conve1} \end{equation} Thus, by using Vitali's theorem we obtain \begin{gather} \int_\Omega \mu(\cdot,\theta^m)\mathbf{T}\left(D(\mathbf{u}^m)\right): D(\mathbf{\Phi})\,dx \to \int_\Omega \mu(\cdot,\theta )\mathbf{T}\left(D(\mathbf{u})\right): D(\mathbf{\Phi})\,dx,\label{conve2}\\ \int_\Omega \kappa(\cdot,\theta^m)\mathbf{a} \left(\nabla\theta^m\right)\cdot\nabla\phi\,dx \to \int_\Omega \kappa(\cdot,\theta)\mathbf{a} \left(\nabla\theta\right)\cdot\nabla\phi\,dx.\label{conve3} \end{gather} Once we have \eqref{conve2} and \eqref{conve3} we conclude the proof of Theorem \ref{main-theorem}. In Subsections 3.3 and 3.4, we analyze the convergence of $[D(\mathbf{u}^m),\nabla(\theta^m)]$ to $[D(\mathbf{u}),\theta]$ almost everywhere in $\Omega$. This part is closely related to \cite[Sections 3 and 4]{Frehse}; however, we expose it with some details for the reader's convenience. \subsection{Almost everywhere convergence of $D(\mathbf{u}^m)$ to $D(\mathbf{u})$} To prove the convergence of $D(\mathbf{u}^m)$ to $D(\mathbf{u})$ almost everywhere in $\Omega$, we prove that for an arbitrary $\eta_1>0$, there exists a subsequence of $(\mathbf{u}^m)_{m\in\mathbb{N}}$, still denoted by $(\mathbf{u}^m)_{m\in\mathbb{N}}$, such that for some $\rho_1\in(0,1)$, it holds that \begin{equation}\label{conve4} \lim_{m\to \infty}\int_\Omega[(\mathbf{T}\left(D(\mathbf{u}^m)\right) -\mathbf{T}\left(D(\mathbf{u})\right)):D(\mathbf{u}^m-\mathbf{u})]^{\rho_1} dx \leq \eta_1. \end{equation} Following \cite[Section 3]{Frehse}, we first consider a decomposition of the pressure. Consider the Stokes problems \begin{equation}\label{Stokes} \begin{gathered} -\Delta \mathbf{u}^{I_m}+\nabla\pi^{I_m}=\mathbf{H}^{I_m}\quad \text{in } \Omega,\; I=1,2,3,4,5,\\ \operatorname{div}\ \mathbf{u}^{I_m}=0 \quad \text{in } \Omega,\\ \mathbf{u}^{I_m}=\mathbf{0} \quad \text{on } \partial\Omega, \end{gathered} \end{equation} where \begin{gather*} \mathbf{H}^{1_m}=-\operatorname{div}\left(\mu(\cdot,\theta^m)\mathbf{T} (D(\mathbf{u}^m))\right)\in \mathbf{W}^{-1,p'}(\Omega),\\ \mathbf{H}^{2_m}=\operatorname{div}\left( \mathbf{u}^m\otimes(\mathbf{u}^m -\mathbf{u})\right)\in \mathbf{W}^{1,r'}(\Omega),\\ \mathbf{H}^{3_m}=\operatorname{div}\left((\mathbf{u}^m-\mathbf{u})\otimes \mathbf{u}\right)\in \mathbf{W}^{1,r'}(\Omega),\quad \mathbf{H}^{4_m}=\frac{1}{m}| \mathbf{u}^m|^{t-2}\mathbf{u}^m \in \mathbf{L}^{t'}(\Omega),\\ \mathbf{H}^{5_m}=-\theta^m\mathbf{f}\in \mathbf{L}^{q}(\Omega). \end{gather*} It is well known that there exists a weak solution $[\mathbf{u}^{I_m},\pi^{I_m}]$ of \eqref{Stokes}, for $I=1,2,3,4,5$; that is, there exist \begin{gather*} \begin{aligned} &[\mathbf{u}^{1_m},\mathbf{u}^{2_m},\mathbf{u}^{3_m},\mathbf{u}^{4_m}, \mathbf{u}^{5_m}] \\ &\in \mathbf{W}_0^{1,p'}(\Omega)\times \mathbf{W}^{1,r'}_0(\Omega)\times \mathbf{W}^{1,r'}_0(\Omega) \times \mathbf{W}^{2,t'}(\Omega)\times \mathbf{W}^{2,q'}(\Omega), \end{aligned} \\ [\pi^{1_m},\pi^{2_m},\pi^{3_m},\pi^{4_m},\pi^{5_m}]\in L^{p'} (\Omega)\times L^{r'}(\Omega)\times L^{r'}(\Omega) \times W^{1,t'}(\Omega)\times W^{1,q'}(\Omega), \end{gather*} satisfying \begin{equation}\label{weak_stokes} \int_\Omega\nabla \mathbf{u}^{I_m}:\nabla\mathbf{\Phi}dx -\int_\Omega \pi^{I_m}\operatorname{div}\mathbf{\Phi}dx =\langle \mathbf{H}^{I_m},\mathbf{\Phi}\rangle,\quad \forall \mathbf{\Phi}\in\mathbf{C}^\infty_0(\Omega),\; I=1,2,3,4,5. \end{equation} Moreover, the following estimates hold: \begin{gather}\label{stokes3} \| \pi^{1_m}\|_{p'}\leq C\| \mathbf{H}^{1_m}\|_{{-1,p'}} \leq C\mu_2\| \mathbf{T}(D(\mathbf{u}^m))\|_{p'},\\ \| \pi^{2_m}\|_{r'}\leq C\| \mathbf{H}^{2_m}\|_{{-1,r'}} \leq C\| \mathbf{u}^m\otimes(\mathbf{u}^m-\mathbf{u})\|_{r'} \leq C\| \mathbf{u}^m\|_{2r'}\| \mathbf{u}^m-\mathbf{u}\|_{2r'},\label{stokes4}\\ \| \pi^{3_m}\|_{r'}\leq C\| \mathbf{H}^{3_m}\|_{{-1,r'}} \leq C\| \mathbf{u}^m\otimes(\mathbf{u}^m-\mathbf{u})\|_{r'} \leq C\| \mathbf{u}^m\|_{2r'}\| \mathbf{u}^m-\mathbf{u}\|_{2r'},\label{stokes5}\\ \| \nabla\pi^{4_m}\|_{t'}\leq C\| \mathbf{H}^{4_m}\|_{t'} \leq \frac{C}{m^{1/t}}\big( \frac{1}{m^{1/t}} \| \mathbf{u}^m\|_{t}\big)^{t-1},\label{stokes6}\\ \| \nabla\pi^{5_m}\|_{q}\leq C\| \mathbf{H}^{5_m}\|_{q} \leq \| \theta^m\|_q\| \mathbf{f}\|_\infty. \end{gather} Since $2r'=np/(n-p)$, from \eqref{estt1}, \eqref{stokes4} and \eqref{stokes5}, as $m$ goes to $\infty$, we obtain \begin{equation}\label{conve6} [\pi^{2_m},\pi^{3_m}]\to [0,0]\ \text{in}\ L^{s}(\Omega) \times L^{s}(\Omega)\quad \text{for all } s\in[1,r'). \end{equation} Furthermore, by using \eqref{estimate-aprox}, as $m$ goes to $\infty$, it holds that \begin{equation}\label{convet} \nabla\pi^{4_m}\to \mathbf{0}\quad \text{in } \mathbf{L}^{t'}(\Omega). \end{equation} Adding the weak formulations \eqref{weak_stokes} and using \eqref{formulation-pressure21} we obtain \begin{equation}\label{dma1} \begin{aligned} &\sum_{I=1}^5[\int_\Omega\nabla \mathbf{u}^{I_m}:\nabla\mathbf{\Phi}dx -\int_\Omega \pi^{I_m}\operatorname{div}\mathbf{\Phi}dx] \\ &=\int_\Omega (\mathbf{u}\otimes\mathbf{u}):D(\mathbf{\Phi})\,dx +\int_\Omega \pi^m\operatorname{div}\mathbf{\Phi}\,dx,\quad \forall\mathbf{\Phi}\in \mathbf{W}^{1,r}_0(\Omega). \end{aligned} \end{equation} Taking $\mathbf{\Phi}\in\mathbf{V}_r$ in \eqref{dma1} we obtain \begin{equation}\label{dma2} \sum_{I=1}^5\int_\Omega\nabla \mathbf{u}^{I_m}:\nabla\mathbf{\Phi}dx =\int_\Omega (\mathbf{u}\otimes\mathbf{u}):D(\mathbf{\Phi})\,dx\quad \forall\mathbf{\Phi}\in \mathbf{V}_r. \end{equation} From \eqref{dma2} and using that $\sum_{I=1}^5\mathbf{u}^{I_m}\in \mathbf{W}_0^{1,r'}(\Omega)$ we obtain \begin{equation}\label{dma3} \sum_{I=1}^5\mathbf{u}^{I_m}=\mathbf{U}\in \mathbf{W}_0^{1,r'}(\Omega),\quad \forall m\in\mathbb{N}. \end{equation} Finally, taking the term $\int_\Omega \pi^m\operatorname{div}\mathbf{\Phi}\,dx$ in \eqref{dma1}, replacing it in \eqref{formulation-pressure21} and using \eqref{dma3} we obtain \begin{equation}\label{formulation-new2} \begin{aligned} &\int_\Omega \mu(\cdot,\theta^m)\mathbf{T}\left(D(\mathbf{u}^m)\right): D(\mathbf{\Phi})\,dx +\frac1{m}\int_\Omega | \mathbf{u}^m|^{t-2}\mathbf{u}^m\cdot\mathbf{\Phi}\,dx \\ &= \int_\Omega (\mathbf{u}^m\otimes\mathbf{u}^m):D(\mathbf{\Phi})\,dx -\int_\Omega (\mathbf{u}\otimes\mathbf{u}):D(\mathbf{\Phi})\,dx \\ &\quad +\int_\Omega\nabla \mathbf{U}:\nabla\mathbf{\Phi}dx -\sum_{I=1}^5\int_\Omega\pi^{I_m}\operatorname{div}\mathbf{\Phi}\,dx +\int_\Omega\theta^m\mathbf{f}\cdot\mathbf{\Phi}dx, \end{aligned} \end{equation} for all $\mathbf{\Phi}\in \mathbf{W}_0^{1,r}(\Omega)$. Now we are in a position to prove \eqref{conve4}. Let us define \begin{equation}\label{conve8} \mathcal{X}^m:= C(1+| D(\mathbf{u}^m)|^p+| D(\mathbf{u})|^p+| \pi^{1_m}|^{p'}). \end{equation} Then, from \eqref{est20} and \eqref{stokes3} we have \begin{equation}\label{dma30} \int_\Omega\mathcal{X}^m dx\leq K_1, \end{equation} for some positive constant $K_1$ independent on $m$. Fixed $p\in(\frac{2n}{n+2},\frac{2n}{n+1}]$, let $\varepsilon_1>0$ be small enough to be chosen below (see \eqref{epsilon1}). Then, from \cite[Proposition 4.1]{Frehse} there exists a subsequence of $(\mathbf{u}^m)_{m\in\mathbb{N}}$, still denoted by $(\mathbf{u}^m)_{m\in\mathbb{N}}$, and $\lambda_1\geq \frac{1}{\varepsilon_1}$ (independent on $m$), such that \begin{equation}\label{conve9} \int_{\mathcal{B}^m_{\lambda_1}} \mathcal{X}^m dx\leq \varepsilon_1,\quad \mathcal{B}^m_{\lambda_1}:=\{x\in\Omega: {\lambda_1}\lambda_1^2 \}, \end{gather*} it holds \begin{gather}\label{conve12} | \mathcal{A}^m_{\lambda_1}| \leq | \mathcal{B}^m_{\lambda_1} | +| \mathcal{C}^m_\lambda |, \\ | \mathcal{A}^m_{\lambda_1}|+| \mathcal{B}^m_{\lambda_1} | \leq \frac{C}{\lambda_1^p}\| \nabla (\mathbf{u}^m-\mathbf{u})\|^p_{p}, \label{conve13}\\ | \mathcal{C}^m_{\lambda_1} | \leq \frac{C}{\lambda_1^{2p}}\| \nabla (\mathbf{u}^m-\mathbf{u})\|^p_{p}, \label{conve14}\\ \| \nabla(\mathbf{u}^m-\mathbf{u})_{\lambda_1}\|^p_p \leq C \| \nabla(\mathbf{u}^m-\mathbf{u})\|^p_p\leq K_1.\label{conve21} \end{gather} Now, we consider $(\mathbf{u}^m-\mathbf{u})_{\lambda_1}$ as a test function in \eqref{formulation-new2} and add in both sides of the obtained equation the term \begin{equation} -\int_\Omega \mu(\cdot,\theta^m)\mathbf{T}\left(D(\mathbf{u})\right): D((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx, \end{equation} to obtain \begin{equation}\label{formulation-new4} \begin{aligned} &\int_\Omega \mu(\cdot,\theta^m)[\mathbf{T}\left(D(\mathbf{u}^m)\right) -\mathbf{T}\left(D(\mathbf{u})\right)]: D((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx\\ & +\frac{1}{m}\int_\Omega | \mathbf{u}^m|^{t-2}\mathbf{u}^m\cdot(\mathbf{u}^m -\mathbf{u})_{\lambda_1}\,dx \\ &= \int_\Omega [(\mathbf{u}^m\otimes\mathbf{u}^m)-(\mathbf{u}\otimes\mathbf{u})]: D((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx\\ &\quad +\int_\Omega\nabla \mathbf{U}:\nabla(\mathbf{u}^m-\mathbf{u})_{\lambda_1} dx \\ &\quad -\sum_{I=1}^5\int_\Omega\pi^{I_m}\operatorname{div} ((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx +\int_\Omega\theta^m\mathbf{f}\cdot(\mathbf{u}^m-\mathbf{u})_{\lambda_1} dx. \end{aligned} \end{equation} Notice that $\mathbf{u}^m-\mathbf{u}=(\mathbf{u}^m-\mathbf{u})_{\lambda_1}$ on $\Omega\setminus \mathcal{A}^m_{\lambda_1}$, and then, $\operatorname{div}(\mathbf{u}^m-\mathbf{u})_{\lambda_1}=0$ almost everywhere on $\Omega\setminus \mathcal{A}^m_{\lambda_1}$. Therefore, from \eqref{formulation-new4} we obtain \begin{align*} Z^m&:=\int_{\Omega\setminus \mathcal{A}^m_{\lambda_1}} \mu(\cdot,\theta^m)[\mathbf{T}\left(D(\mathbf{u}^m)\right) -\mathbf{T}\left(D(\mathbf{u})\right)]: D((\mathbf{u}^m-\mathbf{u}))\,dx \\ &= -\int_{\mathcal{A}^m_{\lambda_1}}\mu(\cdot,\theta^m)[\mathbf{T} \left(D(\mathbf{u}^m)\right)-\mathbf{T}\left(D(\mathbf{u})\right)]: D((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx \\ &\quad -\int_{\mathcal{A}^m_{\lambda_1}}\pi^{1_m} \operatorname{div}((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx \\ &\quad +\int_\Omega [(\mathbf{u}^m\otimes(\mathbf{u}^m-\mathbf{u}) +(\mathbf{u}^m-\mathbf{u})\otimes \mathbf{u}]: D((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx \\ &\quad +\int_\Omega[\nabla \mathbf{U}-\mu(\cdot,\theta^m)\mathbf{T} \left(D(\mathbf{u})\right)]:\nabla(\mathbf{u}^m-\mathbf{u})_{\lambda_1} dx \\ &\quad +\int_\Omega [ \nabla \pi^{4_m}-\frac{1}{m} \int_\Omega | \mathbf{u}^m|^{t-2}\mathbf{u}^m]\cdot (\mathbf{u}^m-\mathbf{u})_{\lambda_1} dx \\ &\quad -\int_{\mathcal{A}^m_{\lambda_1}}(\pi^{2_m}+\pi^{3_m}+\pi^{5_m}) \operatorname{div}((\mathbf{u}^m-\mathbf{u})_{\lambda_1})\,dx +\int_\Omega\theta^m\mathbf{f}\cdot(\mathbf{u}^m-\mathbf{u})_{\lambda_1} dx\\ &:=\sum^7_{i=1}Z^m_i. \end{align*} From \eqref{est20}, \eqref{estt1}, \eqref{estt2}, \eqref{conve6}, \eqref{convet}, \eqref{estt3}, \eqref{dma5} we obtain \begin{equation}\label{limf1} \lim_{m\to \infty} (Z^m_3+Z^m_5+Z^m_6+Z^m_7)=0. \end{equation} Moreover, from \eqref{dma5} and since $\nabla \mathbf{U}\in \mathbf{L}^{r'}(\Omega)$, $\mu_1\leq \mu(x,\theta)\leq \mu_2,\ \text{a.e.}\ x\in\Omega$ and $\mathbf{T}(D(\mathbf{u}))\in \mathbf{L}^{p'}(\Omega)$, we obtain that $\lim_{m\to \infty} Z^m_4=0$. Now we deal with $\lim_{m\to \infty} Z^m_1+Z^m_2$. From the H\"older inequality, \eqref{estt3} and \eqref{conve12}-\eqref{conve21} it holds that \begin{equation} \begin{aligned} &| Z^m_1+Z^m_2| \\ &\leq \Big| \int_{\mathcal{B}^m_{\lambda_1}\cup {\mathcal{C}^m_{\lambda_1}}} ( \mu(\cdot,\theta^m)[\mathbf{T}\left(D(\mathbf{u}^m)\right) -\mathbf{T}\left(D(\mathbf{u})\right)]:D((\mathbf{u}^m -\mathbf{u})_{\lambda_1}) \\ &\quad -\pi^{1_m}\operatorname{div}((\mathbf{u}^m -\mathbf{u})_{\lambda_1})\,)dx\Big| \\ &\leq \mu_2\tau_2C\Big(\int_{\mathcal{B}^m_{\lambda_1}} \mathcal{X}^mdx\Big)^{1/p'}\| \nabla(\mathbf{u}^m -\mathbf{u})_{\lambda_1}\|_{p,\mathcal{B}^m_{\lambda_1}} \\ &\quad +\mu_2\tau_2C{\lambda_1}\Big(\int_{\mathcal{C}^m_{\lambda_1}} \mathcal{X}^mdx\Big)^{1/p'}| \mathcal{C}^m_{\lambda_1}|^{1/p} \\ &\leq C\mu_2\tau_2(\varepsilon_1^{1/p'}K_1^{1/p} +C{\lambda_1} K_1^{1/p'}(C{\lambda_1}^{-2p}K_1)^{1/p}) \\ &\leq C\mu_2\tau_2(\varepsilon_1^{1/p'}K_1^{1/p} +\frac{K_1}{{\lambda_1}}) \\ &\leq C\mu_2\tau_2(\varepsilon_1^{1/p'}K_1^{1/p} +\varepsilon_1K_1). \end{aligned} \label{dma11} \end{equation} In \eqref{dma11}, $\| \nabla(\mathbf{u}^m-\mathbf{u})_{\lambda_1}\|_{p,\mathcal{B}^m_{\lambda_1}}$ denotes the $L^p(\mathcal{B}^m_{\lambda_1})$-norm of $\nabla(\mathbf{u}^m-\mathbf{u})_{\lambda_1}$. Since $\lim_{m\to \infty} Z^m_4=0$, from \eqref{limf1} and \eqref{dma11} we obtain \begin{equation}\label{dma23} \lim_{m\to \infty}Z^m\leq C\mu_2\tau_2(\varepsilon_1^{1/p'}K_1^{1/p} +\varepsilon_1K_1). \end{equation} Therefore, fixed $\rho_1\in(0,1)$, by using the H\"older inequality and \eqref{dma30} we obtain \begin{equation}\label{conve41} \begin{aligned} S^m&:= \mu_1^{\rho_1}\int_\Omega[(\mathbf{T}\left(D(\mathbf{u}^m)\right) -\mathbf{T}\left(D(\mathbf{u}))\right):D(\mathbf{u}^m-\mathbf{u})]^{\rho_1} dx \\ &\leq \int_{\Omega\setminus \mathcal{A}^m_{\lambda_1}} [\mu(\cdot,\theta^m)(\mathbf{T}\left(D(\mathbf{u}^m)\right) -\mathbf{T}\left(D(\mathbf{u})\right)):D(\mathbf{u}^m-\mathbf{u})]^{\rho_1} dx \\ &\quad+\int_{\mathcal{A}^m_{\lambda_1}}[\mu(\cdot,\theta^m)(\mathbf{T} \left(D(\mathbf{u}^m)\right)-\mathbf{T}\left(D(\mathbf{u})\right)): D(\mathbf{u}^m-\mathbf{u})]^{\rho_1} dx \\ &\leq (Z^m)^{\rho_1}| \Omega\setminus{\mathcal{A}^m_{\lambda_1}}|^{1-\rho_1} +C(\mu_2\tau_2K_1)^{\rho_1}| {\mathcal{A}^m_{\lambda_1}}|^{1-\rho_1}. \end{aligned} %\label{dma20} \end{equation} Then, taking $\varepsilon_1>0$ small enough such that \begin{equation}\label{epsilon1} C(\mu_2\tau_2)^{\rho_1}| \Omega|^{1-\rho_1}(\varepsilon_1^{1/p'}K_1^{1/p} +\varepsilon_1K_1)^{\rho_1}+C(\mu_2\tau_2)^{\rho_1} K_1\varepsilon_1^{p(1-\rho_1)}<\rho_1, \end{equation} from \eqref{dma23}-\eqref{epsilon1} we have \begin{equation}\label{dma21} \begin{aligned} &\lim_{m\to \infty}S^m \\ &\leq | \Omega|^{1-\rho_1}C(\mu_2\tau_2)^{\rho_1}(\varepsilon_1^{1/p'}K_1^{1/p} +\varepsilon_1K_1)^{\rho_1}+C(\mu_2\tau_2K_1)^{\rho_1} (C\lambda_1^{-p}K_1)^{1-\rho_1} \\ &\leq C(\mu_2\tau_2)^{\rho_1}| \Omega|^{1-\rho_1}(\varepsilon_1^{1/p'}K_1^{1/p} +\varepsilon_1K_1)^{\rho_1}+C(\mu_2\tau_2)^{\rho_1}K_1\varepsilon_1^{p(1-\rho_1)} <\rho_1. \end{aligned} \end{equation} Thus, we conclude \eqref{conve4} and therefore the convergence of $D(\mathbf{u}^m)$ to $D(\mathbf{u})$ almost everywhere in $\Omega$. \subsection{Almost everywhere convergence of $\nabla \theta^m$ to $\nabla \theta$} Let be a fixed value $p\in(2n/(n+2), 2n/(n+1)]$ and $q>\frac{2np}{p(n+2)-n}$. To prove the convergence of $\nabla \theta^m$ to $\nabla \theta$ almost everywhere in $\Omega$, we proceed in the same spirit of Subsection 3.3. We prove that for an arbitrary $\eta_2>0$, there exists a subsequence of $(\theta^m)_{m\in\mathbb{N}}$, still denoted by $(\theta^m)_{m\in\mathbb{N}}$, such that for some $\rho_2\in(0,1)$, it holds that \begin{equation}\label{conve4b} \lim_{m\to \infty}\int_\Omega[(\mathbf{a}\left(\nabla\theta^m\right) -\mathbf{a}\left(\nabla\theta\right))\cdot\nabla(\theta^m-\theta)]^{\rho_2} dx \leq \eta_2. \end{equation} Let us define \begin{equation}\label{conve8b} \mathcal{E}^m:= C(1+| \nabla\theta^m|^q+| \nabla\theta|^q). \end{equation} Then, from \eqref{est20} we have \begin{equation}\label{dma30b} \int_\Omega\mathcal{E}^m dx\leq K_2, \end{equation} for some positive constant $K_2$ independent on $m$. Let $\varepsilon_2>0$ small enough to be chosen below (see \eqref{epsi2}). Reasoning as in Subsection 3.3 (see also \cite[Proposition 4.1]{Frehse}), there exists a subsequence of $(\theta^m)_{m\in\mathbb{N}}$, still denoted by $(\theta^m)_{m\in\mathbb{N}}$, and $\lambda_2\geq \frac{1}{\varepsilon_2}$ (independent on $m$), such that \begin{equation}\label{conve9b} \int_{\mathcal{D}^m_{\lambda_2}} \mathcal{E}^m dx \leq \varepsilon_2,\quad \mathcal{D}^m_{\lambda_2} :=\{x\in\Omega: \lambda_2\lambda_2^2 \}, \end{gather*} it holds \begin{gather}\label{conve12b} | \mathcal{F}^m_{\lambda_2}| \leq | \mathcal{D}^m_{\lambda_2} | +| \mathcal{G}^m_{\lambda_2} |,\\ | \mathcal{D}^m_{\lambda_2}|+| \mathcal{F}^m_{\lambda_2} | \leq \frac{C}{\lambda_2^q}\| \nabla (\theta^m-\theta)\|^q_{q},\label{conve13b}\\ | \mathcal{G}^m_{\lambda_2} | \leq \frac{C}{\lambda_2^{2q}}\| \nabla (\theta^m-\theta)\|^q_{q},\label{conve14b}\\ \| \nabla(\theta^m-\theta)_{\lambda_2}\|^q_q \leq C \| \nabla(\theta^m-\theta)\|^q_q\leq K_2.\label{conve21b} \end{gather} Now we consider $(\theta^m-\theta)_{\lambda_2}$ as a test function in \eqref{prob-theta}, and add in both sides of the obtained equation the term \begin{equation} -\int_\Omega \kappa(\cdot,\theta^m)\mathbf{a}\left(\nabla\theta\right)\cdot \nabla((\theta^m-\theta)_{\lambda_2})\,dx, \end{equation} this gives \begin{equation}\label{formulation-new4b} \begin{aligned} &\int_\Omega \kappa(\cdot,\theta^m)[\mathbf{a}\left(\nabla\theta^m\right) -\mathbf{a}\left(\nabla\theta\right)]\cdot \nabla((\theta^m-\theta)_{\lambda_2})\,dx \\ &= \int_\Omega (\theta^m\mathbf{u}^m)\cdot\nabla((\theta^m-\theta)_{\lambda_2})\,dx -\int_\Omega \kappa(\cdot,\theta^m)\mathbf{a}\left(\nabla\theta\right)\cdot \nabla((\theta^m-\theta)_{\lambda_2})\,dx \\ &\quad -\gamma\int_\Gamma\theta^m(\theta^m-\theta)_{\lambda_2}\,d\Gamma +\langle g,(\theta^m-\theta)_{\lambda_2}\rangle_{(W^{1,q}(\Omega))'}\\ &\quad +\int_\Gamma h(\theta^m-\theta)_{\lambda_2}\,d\Gamma. \end{aligned} \end{equation} Notice that $\theta^m-\theta=(\theta^m-\theta)_{\lambda_2}$ on $\Omega\setminus \mathcal{F}^m_{\lambda_2}$. Therefore, from \eqref{formulation-new4b} we obtain \begin{equation} \begin{aligned} Y^m &:= \int_{\Omega\setminus \mathcal{F}^m_{\lambda_2}}\kappa(\cdot,\theta^m) [\mathbf{a}(\nabla\theta^m)-\mathbf{a}(\nabla\theta)]\cdot \nabla(\theta^m-\theta)\,dx \\ &= -\int_{\mathcal{F}^m_{\lambda_2}}\kappa(\cdot,\theta^m) [\mathbf{a}(\nabla\theta^m)-\mathbf{a}(\nabla\theta)]\cdot \nabla((\theta^m-\theta)_{\lambda_2})\,dx \\ &\quad +\int_\Omega (\theta^m\mathbf{u}^m)\cdot\nabla((\theta^m-\theta)_{\lambda_2}) \,dx \\ &\quad -\int_\Omega \kappa(\cdot,\theta^m)\mathbf{a}\left(\nabla\theta\right)\cdot \nabla((\theta^m-\theta)_{\lambda_2})\,dx -\gamma\int_\Gamma\theta^m(\theta^m-\theta)_{\lambda_2}\,d\Gamma \\ &\quad +\langle g,(\theta^m-\theta)_{\lambda_2}\rangle_{(W^{1,q}(\Omega))'} +\int_\Gamma h(\theta^m-\theta)_{\lambda_2}\,d\Gamma:=\sum^6_{i=1}Y^m_i. \end{aligned} \end{equation} Using \eqref{est20} we obtain \begin{equation}\label{limf1b} \lim_{m\to \infty} (Y^m_2+Y^m_4+Y^m_5+Y^m_6)=0. \end{equation} Moreover, from \eqref{dma5b} and since $\kappa_1\leq \kappa(x,\theta)\leq \kappa_2$, a.e. $x\in\Omega$ and $\mathbf{a}(\nabla\theta)\in \mathbf{L}^{q'}(\Omega)$, we obtain that $\lim_{m\to \infty} Y^m_3=0$. Now we deal with $\lim_{m\to \infty} Y^m_1$. From the properties of $(\theta^m-\theta)_{\lambda_2}$, the H\"older inequality and \eqref{conve21b} it holds \begin{equation} \begin{aligned} | Y^m_1| &\leq \big| \int_{\mathcal{D}^m_{\lambda_2} \cup {\mathcal{G}^m_{\lambda_2}}}( \kappa(\cdot,\theta^m) [\mathbf{a}(\nabla\theta^m)-\mathbf{a}(\nabla\theta)] \cdot\nabla((\theta^m-\theta)_{\lambda_2}))dx\big| \\ &\leq \kappa_2\alpha_2C \Big(\int_{\mathcal{D}^m_\lambda}\mathcal{E}^mdx\Big)^{1/q'}\| \nabla(\theta^m-\theta)_{\lambda_2}\|_{q,\mathcal{D}^m_{\lambda_2}} \\ &\quad +\kappa_2\alpha_2C{\lambda_2}\Big(\int_{\mathcal{G}^m_{\lambda_2}} \mathcal{E}^mdx\Big)^{1/q'}| \mathcal{G}^m_\lambda|^{1/q} \\ &\leq C\kappa_2\alpha_2(\varepsilon_2^{1/q'}K_2^{1/q} +{\lambda_2} K_2^{1/q'}(C{\lambda_2}^{-2q}K_2)^{1/q}) \\ &\leq C\kappa_2\alpha_2(\varepsilon_2^{1/q'}K_2^{1/q} +\frac{K_2}{\lambda_2})\leq C\kappa_2\alpha_2(\varepsilon_2^{1/q'}K_2^{1/q} +\varepsilon_2K_2). \end{aligned}\label{dma11b} \end{equation} Thus, since $\lim_{m\to \infty} Y^m_3=0$, from \eqref{limf1b} and \eqref{dma11b} we obtain $\lim_{m\to \infty}Y^m\leq C\kappa_2\alpha_2(\varepsilon_2^{1/q'}K_2^{1/q} +\varepsilon_2K_2)$. Therefore, fixed $\rho_2\in(0,1)$ we obtain \begin{equation}\label{conve41b} \begin{aligned} L^m&:= \kappa_1^{\rho_2}\int_\Omega[(\mathbf{a} \left(\nabla\theta^m\right)-\mathbf{a}\left(\nabla\theta)\right) \cdot\nabla(\theta^m-\theta)]^{\rho_2} dx \\ &\leq \int_{\Omega\setminus \mathcal{F}^m_{\lambda_2}} [\kappa(\cdot,\theta^m)(\mathbf{a}\left(\nabla\theta^m\right) -\mathbf{a}\left(\nabla\theta)\right)\cdot\nabla(\theta^m-\theta)]^{\rho_2} dx \\ &\quad +\int_{\mathcal{F}^m_{\lambda_2}}[\kappa(\cdot,\theta^m)(\mathbf{a} \left(\nabla\theta^m\right)-\mathbf{a}\left(\nabla\theta)\right) \cdot\nabla(\theta^m-\theta)]^{\rho_2} dx \\ &\leq (Y^m)^{\rho_2}| \Omega\setminus{\mathcal{F}^m_{\lambda_2}}|^{1-\rho_2} +C(\kappa_2\alpha_2K_2)^{\rho_2}| {\mathcal{F}^m_{\lambda_2}}|^{1-\rho_2}. \end{aligned}%\label{dma20b} \end{equation} Then, taking $\varepsilon_2>0$ small enough such that \begin{equation}\label{epsi2} C| \Omega|^{1-\rho_2}(\kappa_2\alpha_2)^{\rho_2}(\varepsilon_2^{1/q'}K_2^{1/q} +\varepsilon_2K_2)^{\rho_2}+C(\kappa_2\alpha_2)^{\rho_2} K_2\varepsilon_2^{q(1-\rho_2)}<\rho_2, \end{equation} we have \begin{equation}\label{dma21b} \begin{aligned} &\lim_{m\to \infty}Y^m \\ &\leq C| \Omega|^{1-\rho_2}(\alpha_2\kappa_2)^{\rho_2}(\varepsilon_2^{1/q'} K_2^{1/q}+\varepsilon_2K_2)^{\rho_2}+C(\kappa_2\alpha_2 K_2)^{\rho_2}(C{\lambda_2}^{-q})^{1-\rho_2} \\ &\leq C| \Omega|^{1-\rho_2}(\kappa_2\alpha_2)^{\rho_2}(\varepsilon_2^{1/q'} K_2^{1/q}+\varepsilon_2K_2)^{\rho_2}+C(\kappa_2\alpha_2)^{\rho_2} K_2\varepsilon_2^{q(1-\rho_2)} <\rho_2. \end{aligned} \end{equation} Thus, we conclude \eqref{conve4b} and therefore the convergence of $\nabla\theta^m$ to $\nabla\theta$ almost everywhere in $\Omega$. \section{Strong solutions} In this section we analyze the existence of a strong solution considering the tensor stress $\mathbf{T}(\boldsymbol{\eta})=\mathbf{T}_1(\boldsymbol{\eta}) :=2\mu(1+| \boldsymbol{\eta}|^2)^{\frac{p-2}{2}}\boldsymbol{\eta}$ or $\mathbf{T}(\boldsymbol{\eta})=\mathbf{T}_2(\boldsymbol{\eta}) :=2\mu(1+| \boldsymbol{\eta}|)^{{p-2}}\boldsymbol{\eta}$, with $p>1$. We also simplify the boundary conditions on the temperature $\theta$. In fact, we assume Dirichlet boundary condition for the temperature; however, our approach can be adapted in order to analyze other boundary conditions. Indeed, we want to study the existence of strong solution for the problem \begin{equation}\label{eqnwithT} \begin{gathered} -\operatorname{div}(\mathbf{T}(D\mathbf{u})) +\operatorname{div}(\mathbf{u}\otimes\mathbf{u})+\nabla\pi =\theta\mathbf{f}\quad \text{in } \Omega,\\ \operatorname{div} \mathbf{u}=0 \quad \text{in } \Omega,\\ -\operatorname{div}(\kappa(\cdot,\theta)\nabla \theta) +\mathbf{u}\cdot\nabla\theta =g \quad \text{in } \Omega,\\ \mathbf{u}=0\quad \text{on } \partial\Omega,\\ \theta=0\quad \text{on } \partial\Omega, \end{gathered} \end{equation} with $\mathbf{T}$ defined as above. Also, throughout this section we assume that $\kappa:\Omega\times \mathbb{R}\to \mathbb{R}$ is a $C^1$-function such that $0<\kappa_1\leq \kappa(x,\theta)\leq \kappa_2$ a.e. $x\in\Omega$ and for all $\theta\in\mathbb{R}$ and, it satisfies $|\kappa'(\cdot,a)-\kappa'(\cdot,b)|\leq \lambda'| a-b|$, for all $a,b\in\mathbb{R}$ and $\kappa'(\cdot, 0)=0$, with $\kappa_1,\kappa_2$ and $\lambda'$ are positive constants. Under mild conditions on the data $\mathbf{f}\in \mathbf{L}^q,g\in L^r(\Omega)$, we obtain the existence of strong solution $[\mathbf{u},\theta]\in \mathbf{W}^{2,q} (\Omega)\times W^{2,r}(\Omega)$, for $q,r>n$. Our approach is based on regularity results for the Stokes problem and Laplace equation, and a fixed-point argument. Observe that $\mathbf{T}_1$ depends on the differentiable term $| D(\mathbf{u})|^2$ while $\mathbf{T}_2$ depends on the merely Lipschitz continuous term $| D(\mathbf{u})|$; thus, in case $\mathbf{T}=\mathbf{T}_1$ we can use the classical regularity results for the Stokes system to solve the velocity equation for a fixed temperature. However, in the case $\mathbf{T}=\mathbf{T}_2$, in order to overcome the difficulty caused by the lack of regularity of $\mathbf{T}_2$, we first introduce a family of penalized problems, then, we establish existence of penalized strong solutions and finally, we pass to the limit in the sequence of penalized problems, as the penalization term goes to zero. Next, we recall a classical result concerning the existence and uniqueness of solutions to the Stokes system, as well as some technical results. \begin{lemma}[{\cite[Theorem 6.1]{Galdi}}] \label{lemma A1} Let $m\geq -1$ be an integer and let $\Omega$ be a bounded domain in $\mathbb{R}^n$ $(n=2, 3)$ with boundary $\partial\Omega$ of class $C^k$ with $k=(m+2,2)^+$. Then for any $\boldsymbol{\tau}\in \mathbf{W}^{m, \rho}(\Omega)$, the following system \begin{gather*} -\Delta \mathbf{u}+\nabla\pi =\boldsymbol{\tau} \quad \text{in } \Omega,\\ \operatorname{div} \mathbf{u} =0 \quad \text{in } \Omega,\\ \mathbf{u}=0 \quad \text{on } \partial\Omega, \end{gather*} admits a unique solution $[\mathbf{u}, \pi]\in \mathbf{W}^{m+2,\rho}(\Omega)\times W^{m+1,\rho}(\Omega)$. Moreover, the following estimate holds \begin{equation*}\label{A1} \| \nabla \mathbf{u}\| _{{m+1,\rho}}+\| \pi\| _{{m+1,\rho}/\mathbb{R}} \leq C_{m}\| \boldsymbol{\tau}\| _{{m,\rho}}, \end{equation*} where $C_m\equiv C_m(n, \rho, \Omega)$ is a positive constant. \end{lemma} \begin{proposition}[{\cite[Proposition A.4]{Arada}}] \label{A4} Let $\gamma_p=\frac{[(p,3)^+-2]^{(p,3)^{+}-2}}{[(p,3)^+-1]^{(p,3)^{+}-1}}$ and let $F:\mathbb{R}^+\longrightarrow\mathbb{R}$ be defined by \[ F(\delta)=A\delta^2-\delta+E\delta\mathcal{F}(\delta)+D, \] where $A, E, D$ are positive constants and $\mathcal{F}(x)=x^{2r_p}(1+x)^{(p-4)^+}$. Thus, if the following assertion holds \[ AD+ED^{2r_p}(1+D)^{(p-4)^+}\leq\gamma_p, \] then $F$ possesses at least one root $\delta_0$. Moreover, $\delta_0>D$ and for every $\beta\in [1, 2]$ the following estimate holds \[ \frac{\beta-1}{\beta}\delta_0+\frac{2-\beta}{\beta}A\delta_0^2 +\frac{2r_p+1-\beta}{\beta}E\delta_0\mathcal{F}(\delta_0) +\frac{E(p-4)^{+}}{\beta}\delta_0^{2r_p+2}(1+\delta_0)^{(p-4)^{+}-1}\leq D. \] \end{proposition} \begin{proposition}[{\cite[ Proposition A.5]{Arada}}] \label{A5} Let $\gamma_p=\frac{[(p,3)^+-2]^{(p,3)^{+}-2}}{[(p,3)^+-1]^{(p,3)^{+}-1}}$ and let $L:\mathbb{R}^+\longrightarrow\mathbb{R}$ be defined by \[ L(\rho)=A\rho^2-\rho+E\rho\mathcal{G}(\rho)+D, \] where $A, E, D$ are positive constants and $\mathcal{G}(x)=x(1+x)^{(p-3)^+}$. Hence, if the following assertion holds \[ AD+ED(1+D)^{(p-3)^+}\leq\gamma_p, \] then $L$ possesses at least one root $\rho_1$. Moreover, $\rho_1>D$ and for every $\beta\in [1, 2]$ the following estimate holds \[ \frac{\beta-1}{\beta}\rho_1+\frac{2-\beta}{\beta}A\rho_1^2 +\frac{2-\beta}{\beta}E\rho_1\mathcal{G}(\rho_1) +\frac{E(p-3)^{+}}{\beta}\rho_1^3(1+\rho_1)^{(p-3)^{+}-1}\leq D. \] \end{proposition} \begin{theorem}[\cite{tpf}] \label{Banach} Let $X$ and $Y$ be Banach spaces such that $X$ is reflexive and $X\hookrightarrow Y$. Let $B$ be a non-empty, closed, convex and bounded subset of $X$ and let $\mathcal{A}:B\longrightarrow B$ be a mapping such that \[ \|\mathcal{A}(u)-\mathcal{A}(v) \|_Y\leq K \| u-v \|_Y\quad \forall u,v\in B\ (0n$ and $\mathbf{T}=\mathbf{T}_1$, $p>1, \ \mu>0$. There exist positive constants $\overline {C}=\overline {C}(\lambda',\kappa_1,C_{-1},C_0,C_E,C_{\tilde E},C_P)$ and $m_2=m_2(\lambda',c_2,C_P,C_{\tilde E})$ such that if $\| g\|_r/\kappa_1^2 1$, $\mu>0$, $\mathbf{f}\in \mathbf{L}^q(\Omega),g\in L^r(\Omega)$, $q,r>n$. There exist positive constants $M_1=M_1(C_0,C_E,C_P)$ and $m_2=m_2(\lambda',c_2,C_P,C_{\tilde E})$ such that if $\| g\|_r/\kappa_1^2 0$. Here $B_{\delta_0}$ is the closed ball defined in \eqref{ball}. \end{proposition} \begin{proof} Let $[\boldsymbol{\xi}, \omega]\in B_\delta$. From Lemma \ref{lemma A1}, $\mathbf{u}\in \mathbf{V}_{2,q}$ and it satisfies \begin{equation}\label{gradient of u1} \| \nabla \mathbf{u}\| _{1, q} \leq \frac{C_0}{2\mu}\left(\| \omega \mathbf{f}\| _q +\| \boldsymbol{\xi}\cdot\nabla\boldsymbol{\xi}\| _q +\| \operatorname{div}\,(2\mu\sigma(|D\boldsymbol{\xi}|^2) D\boldsymbol{\xi})\| _q\right). \end{equation} Notice that \begin{equation} \begin{aligned} \| \omega \mathbf{f}\| _q &\leq \| \omega\| _\infty\| \mathbf{f}\| _q \leq C_{\tilde E}(C_P+1)\| \nabla \omega\| _r\| \mathbf{f}\| _q \\ &\leq\delta (C_P+1)\| \mathbf{f}\| _q\leq \frac{(C_P+1)^2\delta^2}{2} +\frac{\| \mathbf{f}\|^2 _q}{2}. \end{aligned}\label{estim50} \end{equation} On the other hand, reasoning as in \cite[Proposition 3.1]{Arada} (also see \cite{Arada2}), we obtain \begin{equation}\label{estim51} \| \boldsymbol{\xi}\cdot\nabla\boldsymbol{\xi}\| _q +\| \operatorname{div}\,(2\mu\sigma(|D\boldsymbol{\xi}|^2)D\boldsymbol{\xi})\| _q \leq \frac{4\mu S_p}{C_E}\delta\mathcal{F}(\delta)+\frac{C_P}{C_E}\delta^2, \end{equation} where $\mathcal{F}(x)=x^{2r_p}(1+x)^{(p-4)^+}$. Thus, if $M_1=\frac{C_0}{2}\max\{\frac{4}{C_E}, \frac{C_P}{C_E}+\frac{(C_P+1)^2}{2},\frac{1}{2} \}$, from \eqref{gradient of u1}-\eqref{estim51} we obtain \begin{equation*}\label{delta-estimate1} \| \nabla \mathbf{u}\| _{1, q}\leq \frac{M_1}{\mu} \left(\| \mathbf{f}\| _q^2+\mu S_p\delta\mathcal{F}(\delta)+\delta^2\right). \end{equation*} On the other hand, from classical elliptic regularity (see \cite{Gilbarg}), there exists a constant $c_2>0$ such that \begin{equation}\label{delta-estimate2} \begin{aligned} \| \nabla\theta\| _{1, r} &\leq \frac{c_2}{\kappa_1}\|\kappa'(\cdot,\omega)| \nabla \omega|^2\|_r +\frac{c_2}{\kappa_1}\| \boldsymbol{\xi}\cdot\nabla\omega\| _r +\frac{c_2}{\kappa_1}\| g\|_r \\ &\leq \frac{c_2}{\kappa_1}\lambda'\| \omega\|_\infty\| \nabla \omega\|^2_{2r} + \frac{c_2}{\kappa_1}\| \boldsymbol{\xi}\|_\infty\|\nabla\omega\| _r +\frac{c_2}{\kappa_1}\| g\|_r \\ &\leq \frac{c_2C_{\tilde{E}}(C_P+1)\lambda'}{\kappa_1}\|\nabla\omega\|_r\| \nabla \omega\|^2_{2r}+\frac{c_2}{\kappa_1}(C_P+1) C_E\| \nabla\boldsymbol{\xi}\|_q\|\nabla\omega\|_r \\ &\quad +\frac{c_2}{\kappa_1}\| g\|_r \\ &\leq \frac{c_2C(C_P+1)\lambda'}{\kappa_1C^2_{\tilde{E}}}\delta^3 +\frac{c_2\delta^2}{\kappa_1C_{\tilde{E}}}(C_P+1)+\frac{c_2}{\kappa_1}\| g\|_r. \end{aligned} \end{equation} It can be assumed that $\delta\leq 1$. Thus, in order to ensure that $\mathcal{A}(B_{\delta})\subseteq B_{\delta}$ it is enough to observe that \begin{equation}\label{delta-estimate15} \begin{gathered} \| \nabla \mathbf{u}\| _{1, q} \leq \frac{M_1}{\mu}\left(\| \mathbf{f}\| _q^2 +\mu S_p\delta\mathcal{F}(\delta)+\delta^2\right)\leq \delta,\\ \| \nabla\theta\| _{1, r}\leq \frac{c_2(C_P+1)}{\kappa_1C_{\tilde{E}}} \Big(\frac{C\lambda'}{C_{\tilde E}}+1\Big)\delta^2+\frac{c_2}{\kappa_1}\| g\|_r \leq \delta. \end{gathered} \end{equation} Using Proposition \ref{A4} with $A=\frac {M_1}{\mu}, E=M_1S_p$ and $D=\frac{M_1\| \mathbf{f}\|^2 _q}{\mu}$, there exists $\delta_1>\frac{M_1\| \mathbf{f}\|^2_q}{\mu}$ such that \begin{equation*}\label{delta-estimate16} \frac{M_1}{\mu}\left(\| \mathbf{f}\|^2 _q+\mu S_p\delta_1\mathcal{F}(\delta_1) +\delta_1^2\right)\leq \delta_1, \end{equation*} provided that \[ M_1^2\frac{\| \mathbf{f}\| _q^2}{\mu^2}+M_1S_p \Big(M_1\frac{\| \mathbf{f}\|_q^2}{\mu} \Big)^{2r_p} \Big(1+M_1\frac{\| \mathbf{f}\|_q^2}{\mu}\Big)^{(p-4)^+}\leq\gamma_p, \] which holds from the hypothesis \eqref{condition bola en bola}. Also, it holds ($\beta=2$ in Proposition \ref{A4}) that \[ %\label{estim59} \delta_1\leq \frac{2M_1\| \mathbf{f}\|^2_q}{\mu}. \] On the other hand, we will consider $\| g\|_r$ such that $\frac{\| g\|_r}{\kappa_1^2}<\frac{C^2_{\tilde{E}}}{4c_2^2(C_P+1) (C\lambda'+C_{\tilde{E}})}\equiv m_2$ and $\delta^-3$, $\gamma_p>1/4^{(p-2,1)^{+}}$. Therefore, setting $\overline{C}=(M_1, \overline{C}_2)^+$ and because of \eqref{condicion4.5} implies \eqref{condition bola en bola} and \eqref{estim67}, we see that the proof of Theorem \ref{mainthm1} is a consequence of Propositions \ref{bola en la bola}, \ref{contraction} and Theorem \ref{Banach}. To apply Theorem \ref{Banach} we consider the spaces $X=\mathbf{V}_{2,q}\times (W^{2,r}(\Omega)\cap W_0^{1,r}(\Omega))$ and $Y=\mathbf{W}_0^{1,q}(\Omega)\times W_0^{1,r}(\Omega)$. \end{proof} \subsection{Power law stress for $\mathbf{T}=\mathbf{T}_2$} In this subsection we prove the existence of strong solutions for the boundary-value problem \eqref{eqnwithT} in the case $\mathbf{T}(\boldsymbol{\eta})=\mathbf{T}_2(\boldsymbol{\eta}) =2\mu(1+|\boldsymbol{\eta}|)^{p-2}\boldsymbol{\eta}$. The purpose of this subsection is prove the following theorem. \begin{theorem}\label{exis2} Let $\mathbf{f}\in \mathbf{L}^q(\Omega), g\in L^r(\Omega)$ with $q, r>n$ and $\mathbf{T}=\mathbf{T}_2$, $p>1$, $\mu>0$. There exist positive constants $\overline {\lambda}=\overline {\lambda}(\lambda', \kappa_1,C_0,C_{-1},C_E,C_P,C_{\tilde E})$ and $m_2=m_2(\lambda',c_2,C_P,C_{\tilde E})$ such that if $\| g\|_r/\kappa_1^2 0$. This is the content of the next theorem. \begin{theorem}\label{penalizedthm} Let $\mathbf{f}\in \mathbf{L}^q(\Omega), g\in L^r(\Omega)$ with $q, r>n$, $p>1,\mu>0$, and $0<\varepsilon<1$. There exist positive constants $\overline {\lambda}=\overline {\lambda}(\lambda',\kappa_1, C_0,C_{-1},C_E,C_P,C_{\tilde E})$ and $m_2=m_2(\lambda',c_2,C_P,C_{\tilde E})$, such that if $\| g\|_r/\kappa_1^2 n, p>1$ and $\mu>0$. There exist constants $\overline {\lambda}_1=\overline {\lambda}_1(C_0,C_P, C_E)>0$ and $m_2=m_2(\lambda',c_2,C_P,C_{\tilde E})>0$ such that if $\| g\|_r/\kappa_1^2 0$. Here, $B_{\rho}$ is the closed ball defined in \eqref{ball}. \end{proposition} \begin{proof} From Lemma \ref{lemma A1} and reasoning as in the proof of Proposition \ref{bola en la bola} and \cite[Proposition 4.2]{Arada} we obtain that $\mathbf{u_\varepsilon}\in\mathbf{V}_{2,q}$ and it satisfies \begin{equation}\label{gradient of u} \begin{aligned} \| \nabla \mathbf{u_\varepsilon}\| _{1, q} &\leq \frac{C_0}{2(1+\varepsilon)^{p-2}\mu}\left(\| \omega \mathbf{f}\| _q +2\mu\| \operatorname{div}(\sigma_\varepsilon(|D\boldsymbol{\xi}|^2) D\boldsymbol{\xi})\| _q+\| \boldsymbol{\xi}\cdot\nabla\boldsymbol{\xi}\| _q\right)\\ &\leq \frac{C_0}{\mu}\big(\| \omega\| _\infty \| \mathbf{f}\| _q +2\mu\| \operatorname{div}\,(\sigma_\varepsilon(|D\boldsymbol{\xi}|^2) D\boldsymbol{\xi})\| _q+\| \boldsymbol{\xi}\cdot\nabla\boldsymbol{\xi}\| _q\big)\\ &\leq \frac{C_0}{\mu}\Big(\rho(C_P+1) \| \mathbf{f}\| _q +2\mu\| \operatorname{div}\,(\sigma_\varepsilon(|D\boldsymbol{\xi}|^2) D\boldsymbol{\xi})\| _q +\| \boldsymbol{\xi}\cdot\nabla\boldsymbol{\xi}\| _q\Big)\\ &\leq \frac{C_0}{\mu}\Big(\frac{\rho^2(C_P+1)^2}{2} +\frac{\| \mathbf{f}\| _q^2}{2}+8\mu\overline{S}_p\mathcal{G}(\| D\xi\|_\infty) \| \nabla\xi\|_{1,q}+\frac{C_P}{C_E}\rho^2\Big)\\ &\leq \frac{C_0}{\mu}\Big(\big[\frac{(C_P+1)^2}{2}+\frac{C_P}{C_E}\big] \rho^2+\frac{\| \mathbf{f}\|^2 _q}{2}+8\mu\rho\overline{S}_p\mathcal{G}(\rho) /C_E\Big)\\ &\leq \frac{\overline{\lambda}_1}{\mu} \left(\| \mathbf{f}\|^2 _q+\mu\rho\overline{S}_p\mathcal{G}(\rho)+\rho^2\right), \end{aligned} \end{equation} where $\overline{\lambda}_1=C_0\max\{\frac{1}{2}, \frac{8}{C_E}, \frac{(C_P+1)^2}{2}+\frac{C_P}{C_E} \}, \mathcal{G}(x)=x(1+x)^{(p-3)^+}$ and $\overline{S}_p=(| p-2|,1)^+\,2^{(p-3)^+}$. As in \eqref{delta-estimate2}, we obtain \begin{equation}\label{unif1} \| \nabla\theta_\varepsilon\| _{1, r} \leq \frac{c_2C(C_P+1)\lambda'}{\kappa_1C^2_{\tilde{E}}}\rho^3 +\frac{c_2\rho^2}{\kappa_1C_{\tilde{E}}}(C_P+1)+\frac{c_2}{\kappa_1}\| g\|_r\,. \end{equation} As in the proof of Proposition \ref{bola en la bola}, we can assume $\rho\leq 1$. Thus, in order to have $\mathcal{A}(B_{\rho})\subseteq B_{\rho}$ it is sufficient to notice that \begin{equation*} \frac{\overline{\lambda}_1}{\mu}\left(\| \mathbf{f}\|^2 _q+\mu\rho\overline{S}_p \mathcal{G}(\rho)+\rho^2\right)\leq \rho,\quad \text{and}\quad \frac{c_2(C_P+1)}{\kappa_1C_{\tilde{E}}}\Big(\frac{C\lambda'}{C_{\tilde E}}+1\Big) \rho^2+\frac{c_2}{\kappa_1}\| g\|_r\leq \rho. \end{equation*} By the hypothesis, from Proposition \ref{A5} with $A=\frac {\overline{\lambda}_1}{\mu}, E=\overline{\lambda}_1\overline{S}_p$ and $D=\frac{\overline{\lambda}_1\| \mathbf{f}\|^2 _q}{\mu}$, there exists $\rho_1>\frac{\overline{\lambda}_1\| \mathbf{f}\|^2_q}{\mu}$ such that \begin{equation*}\label{delta-estimate6} \frac{\overline{\lambda}_1}{\mu}\left(\| \mathbf{f}\|^2 _q +\mu\rho_1\overline{S}_p\mathcal{G}(\rho_1)+\rho_1^2\right)\leq \rho_1. \end{equation*} Moreover, \begin{equation}\label{estim9} \rho_1\leq \frac{2\overline{\lambda}_1\| \mathbf{f}\|^2_q}{\mu}\cdot \end{equation} The proof follows in the same way as in the end of the proof the Proposition \ref{bola en la bola}. Namely, we consider $\| g\|_r$ and $\kappa_1$ such that $\frac{\| g\|_r}{\kappa_1^2}3, \gamma_p>1/4^{(p-2,1)^{+}}$. The proof of Theorem \ref{penalizedthm} is a consequence of Propositions \ref{bola en la bola2}, \ref{contraction2} and Theorem \ref{Banach} when taking $\overline{\lambda}=(\overline{\lambda}_1, \overline{\lambda}_0)^+$ and keeping in mind that \eqref{condicion4.9} implies \eqref{condicion4.10} and \eqref{estim67-epsilon}. \end{proof} \begin{proof}[Proof of Theorem \ref{exis2}] The existence of a strong solution $[\mathbf{u}, \theta]\in\mathbf{V}_{2,q}\times (W^{2,r}(\Omega) \cap W^{1,r}_0(\Omega))$ is obtained as the limit of a subsequence of the penalized solutions $[\mathbf{u}_\varepsilon, \theta_\varepsilon]$ provided by Theorem \ref{penalizedthm}. Notice that for each $\varepsilon>0$, $[\mathbf{u}_\varepsilon, \theta_\varepsilon]$ satisfies the following weak formulation \begin{gather}\label{pe1} \begin{aligned} &\int_\Omega \left(2\mu(1+\sqrt{\varepsilon^2 +| D\mathbf{u}_\varepsilon|^2})^{p-2}D\mathbf{u}_\varepsilon\right): D(\mathbf{\Phi})\,dx-\int_\Omega (\mathbf{u}_\varepsilon \otimes\mathbf{u}_\varepsilon):D(\mathbf{\Phi})\,dx \\ &= \int_\Omega\theta_\varepsilon\mathbf{f}\cdot\mathbf{\Phi}\,dx,\quad \forall\mathbf{\Phi}\in \mathbf{V}_p, \end{aligned}\\ \label{pe2} \int_\Omega\kappa(x,\theta_\varepsilon)\nabla\theta_\varepsilon\cdot \nabla\phi\,dx +\int_\Omega\phi\mathbf{u}_\varepsilon\cdot\nabla\theta_\varepsilon\,dx =\int_\Omega g\phi\,dx,\quad\forall\phi\in W^{1,q}_0(\Omega). \end{gather} From \eqref{gradient of u}, \eqref{unif1} and \eqref{estim9} we have that $([\mathbf{u}_\varepsilon, \theta_\varepsilon])_{\varepsilon}$ is uniformly bounded in $\mathbf{V}_{2,q}\times (W^{2,r}(\Omega)\cap W^{1,r}_0(\Omega))$. Then, there exists a subsequence of $([\mathbf{u}_\varepsilon, \theta_\varepsilon])_{\varepsilon}$, still denoted by $([\mathbf{u}_\varepsilon, \theta_\varepsilon])_{\varepsilon}$, and $[\mathbf{u}, \theta]$ such that \begin{gather*} [\mathbf{u}_\varepsilon, \theta_\varepsilon]\rightharpoonup [\mathbf{u}, \theta] \quad \text{weakly in } \mathbf{V}_{2,q}\times (W^{2,r}(\Omega) \cap W^{1,r}_0(\Omega)),\\ [\mathbf{u}_\varepsilon, \theta_\varepsilon]\to [\mathbf{u}, \theta]\quad \text{strongly in } \mathbf{C}^{1,\alpha_1}(\overline{\Omega}) \times C^{1,\alpha_2}(\overline{\Omega}), \alpha_1<1-\frac{n}{q}, \alpha_2<1-\frac{n}{r}. \end{gather*} Thus, recalling that $\kappa:\Omega\times \mathbb{R}\to \mathbb{R}$ is a $C^{1}$-function and passing to the limit as $\varepsilon$ tends to zero in \eqref{pe1}-\eqref{pe2}, we obtain \begin{gather}\label{pe1b} \begin{aligned} &\int_\Omega \left(2\mu(1+| D\mathbf{u}|)^{p-2}D\mathbf{u}\right): D(\mathbf{\Phi})\,dx -\int_\Omega (\mathbf{u}\otimes\mathbf{u}):D(\mathbf{\Phi})\,dx\\ &= \int_\Omega\theta\mathbf{f}\cdot\mathbf{\Phi}\,dx,\quad \forall\mathbf{\Phi}\in \mathbf{V}_p, \end{aligned} \\ \label{pe2b} \int_\Omega\kappa(x,\theta)\nabla\theta\cdot\nabla\phi\,dx +\int_\Omega\phi\mathbf{u}\cdot\nabla\theta\,dx =\int_\Omega g\phi\,dx,\quad\forall\phi\in W^{1,q}_0(\Omega). \end{gather} The regularity of $[\mathbf{u}, \theta]$ follows from \eqref{gradient of u} and \eqref{unif1}. This completes the proof of Theorem \ref{exis2}. \end{proof} \subsection{Uniqueness} We finish this section with the following uniqueness result which gives conditions on the data to ensure that the obtained strong solution agrees with the weak solution. \begin{theorem}\label{uniqueness1} Let $p\geq2$ and consider $[\mathbf{u}_1,\theta_1]$ a weak solution of \eqref{eqnwithT} with $\mathbf{T}=\mathbf{T}_1,\mathbf{T}_2$, and let $[\mathbf{u}_2,\theta_2]$ be a strong solution of \eqref{eqnwithT} provided by Theorem \ref{mainthm1} or Theorem \ref{exis2}. If \begin{align*} 1-\Big(\frac{C_P^4C_k^3}{4\mu^2\kappa_1} \| \mathbf{f}\|_2\| g\|_2+\frac{C^2_PC_k}{2\mu}\| \mathbf{f}\|_2 +\frac{C^3_PC_k}{\kappa^2_1}\| g\|_2 +\frac{2\lambda'C^5_P}{\kappa^3_1}\| g\|_2\Big)>0, \end{align*} then $[\mathbf{u}_1,\theta_1]=[\mathbf{u}_2,\theta_2]$. Here $C_P$ denotes a general Poincar\'e constant and $C_k$ denotes the Korn constant. \end{theorem} \begin{proof} First of all, $[\mathbf{u}_1,\theta_1]$ being a weak solution of \eqref{eqnwithT} implies that $[\mathbf{u}_1,\theta_1] \in \mathbf{V}_p\times W_0^{1,q}(\Omega)$ and it satisfies \begin{gather*} \int_\Omega \mathbf{T}\left(D(\mathbf{u})\right): D(\mathbf{\Phi})\,dx-\int_\Omega (\mathbf{u}\otimes\mathbf{u}): D(\mathbf{\Phi})\,dx= \int_\Omega\theta\mathbf{f}\cdot\mathbf{\Phi}\,dx,\quad \forall\mathbf{\Phi}\in \mathbf{V}_p, \\ \int_\Omega\kappa(x,\theta)\nabla\theta\cdot\nabla\phi\,dx +\int_\Omega\phi\mathbf{u}\cdot\nabla\theta\,dx =\langle g,\phi\rangle_{(W_0^{1,q}(\Omega))'},\quad\forall\phi\in W^{1,q}_0(\Omega). \end{gather*} Considering the difference between the weak formulations of $[\mathbf{u}_1,\theta_1],[\mathbf{u}_2,\theta_2]$, we can obtain \begin{gather*} \begin{aligned} &\int_\Omega (\mathbf{T}\left(D(\mathbf{u}_1))-\mathbf{T}(D(\mathbf{u}_2))\right): D(\mathbf{u}_1-\mathbf{u}_2)dx \\ &=\int_\Omega (\mathbf{u}_1-\mathbf{u}_2) \nabla \mathbf{u}_2(\mathbf{u}_1-\mathbf{u}_2)dx + \int_\Omega(\theta_1-\theta_2)\mathbf{f}\cdot(\mathbf{u}_1-\mathbf{u}_2)dx, \end{aligned} \\ \int_\Omega (\kappa(x,\theta_1)\nabla\theta_1-\kappa(x,\theta_2) \nabla\theta_2)\cdot\nabla(\theta_1-\theta_2)dx =\int_\Omega(\mathbf{u}_1-\mathbf{u}_2)\cdot\nabla\theta_2(\theta_1-\theta_2)dx. \end{gather*} Notice that since $p\geq 2$, we have the strict monotonicity condition $\left(\mathbf{T}(\boldsymbol{\eta})-\mathbf{T}(\boldsymbol{\xi})\right): (\boldsymbol{\eta}-\boldsymbol{\xi})\geq 2\mu| \boldsymbol{\eta}-\boldsymbol{\xi}|^2$. Then, using the H\"older, Poincar\'e and Korn inequalities we obtain \begin{equation}\label{j1} \begin{aligned} &2\mu\| D(\mathbf{u}_1)-D(\mathbf{u}_2)\|_2^2 \\ &\leq \int_\Omega (\mathbf{T}\left(D(\mathbf{u}_1)) -\mathbf{T}(D(\mathbf{u}_2))\right): D(\mathbf{u}_1-\mathbf{u}_2)dx \\ &\leq \| \mathbf{u}_1-\mathbf{u}_2\|_4^2\| \nabla\mathbf{u}_2\|_2 +\| \theta_1-\theta_2\|_4\| \mathbf{u}_1-\mathbf{u}_4\|_2\|\mathbf{f}\|_2 \\ &\leq C_P^2(\| \nabla(\mathbf{u}_1-\mathbf{u}_2)\|_2^2\| \nabla\mathbf{u}_2\|_2 +\| \nabla(\theta_1-\theta_2)\|_2\| \nabla(\mathbf{u}_1 -\mathbf{u}_2)\|_4\|\mathbf{f}\|_2) \\ &\leq C_P^2C_k^2\| D(\mathbf{u}_1-\mathbf{u}_2)\|_2^2\| D\mathbf{u}_2\|_2 \\ &\quad+C^2_PC_k\| \nabla(\theta_1-\theta_2)\|_2\| D(\mathbf{u}_1-\mathbf{u}_2)\|_2 \|\mathbf{f}\|_2. \end{aligned} \end{equation} On the other hand, since \begin{align*} \kappa(x,\theta_1)\nabla\theta_1-\kappa(x,\theta_2)\nabla\theta_2 =\kappa(x,\theta_1)\nabla(\theta_1-\theta_2)+(\kappa(x,\theta_1)-\kappa(x,\theta_2)) \nabla\theta_2, \end{align*} and using the assumptions on the boundedness and regularity of $\kappa$, as well as the H\"older, Poincar\'e and Korn inequalities we obtain \begin{equation} \begin{aligned} &\kappa_1\| \nabla(\theta_1-\theta_2)\|_2^2 \\ &\leq \| \mathbf{u}_1-\mathbf{u}_2\|_4\| \nabla\theta_2\|_2\| \theta_1-\theta_2\|_4 +\lambda'\| \nabla\theta_2\|_2\| \theta_1-\theta_2\|^2_6(\| \theta_1\|_6 +\| \theta_2\|_6) \\ &\leq C^2_PC_k \| D(\mathbf{u}_1-\mathbf{u}_2)\|_2\| \nabla\theta_2\|_2 \| \nabla(\theta_1-\theta_2)\|_2 \\ &\quad +\lambda'C^3_P\| \nabla\theta_2\|_2\| \nabla(\theta_1-\theta_2)\|^2_2 (\| \nabla\theta_1\|_2+\| \nabla\theta_2\|_2). \end{aligned} \end{equation} Note that $\| D\mathbf{u}_2\|_2$ satisfies the estimate \begin{align*} 2\mu\| D\mathbf{u}_2\|_2^2\leq \| \mathbf{u}_2\|_4 \| \theta_2\|_4\| \mathbf{f}\|_2 \leq C_PC_k\| D\mathbf{u}_2\|_2\| \nabla\theta_2\|_2\| \mathbf{f}\|_2, \end{align*} which implies that \begin{equation}\label{j2} \| D\mathbf{u}_2\|_2\leq \frac{C_PC_k}{2\mu}\| \nabla\theta_2\|_2\| \mathbf{f}\|_2. \end{equation} Moreover, \begin{align*} \kappa_1\| \nabla\theta_i\|^2_2 \leq\int_\Omega\kappa(x,\theta_i)|\nabla\theta_i|^2 \leq C_P\| g\|_2\| \nabla\theta_i\|_2,\quad i=1,2, \end{align*} which yields \begin{equation}\label{j3} \kappa_1(\| \nabla\theta_1\|_2+\| \nabla\theta_1\|_2)\leq 2C_P\| g\|_2. \end{equation} Thus, from \eqref{j1}-\eqref{j3} we obtain \begin{equation} \label{j4} \begin{aligned} 2\mu\| D(\mathbf{u}_1-\mathbf{u}_2)\|_2^2 &\leq \frac{C_P^3C_k^3}{2\mu}\| D(\mathbf{u}_1 -\mathbf{u}_2)\|_2^2\| \nabla\theta_2\|_2\|\mathbf{f}\|_2 \\ &\quad +C^2_PC_k\| \nabla(\theta_1-\theta_2)\|_2\| D(\mathbf{u}_1-\mathbf{u}_2)\|_2 \|\mathbf{f}\|_2 \\ &\leq \frac{C_P^4C_k^3}{2\mu\kappa_1}\| D(\mathbf{u}_1-\mathbf{u}_2)\|_2^2 \| g\|_2\|\mathbf{f}\|_2 \\ &\quad +C^2_PC_k\| \nabla(\theta_1-\theta_2)\|_2\| D(\mathbf{u}_1-\mathbf{u}_2)\|_2 \|\mathbf{f}\|_2, \end{aligned} %\label{j5} \end{equation} and \begin{equation} \begin{aligned} &\kappa_1\| \nabla(\theta_1-\theta_2)\|_2^2 \\ &\leq \frac{C^3_PC_k}{\kappa_1} \| D(\mathbf{u}_1-\mathbf{u}_2)\|_2 \| \nabla(\theta_1-\theta_2)\|_2\| g\|_2 +\frac{2\lambda'C^5_P}{\kappa^2_1}\| \nabla(\theta_1-\theta_2)\|^2_2 \| g\|^2_2. \end{aligned}\label{j6} \end{equation} Then, taking $\| \mathbf{f}\|_2,\| g\|_2,\mu,\kappa$ such that \[ 1-\Big(\frac{C_P^4C_k^3}{4\mu^2\kappa_1}\| \mathbf{f}\|_2\| g\|_2 +\frac{C^2_PC_k}{2\mu}\| \mathbf{f}\|_2 +\frac{C^3_PC_k}{\kappa^2_1}\| g\|_2 +\frac{2\lambda'C^5_P}{\kappa^3_1}\| g\|_2\Big)>0, \] from \eqref{j4}, \eqref{j6} we obtain that $[\mathbf{u}_1,\theta_1]=[\mathbf{u}_2,\theta_2]$. \end{proof} \subsection*{Acknowledgments} The authors are grateful to the anonymous referees for their helpful comments. The first author was partially supported by Fondo Nacional de Financiamiento para la Ciencia, la Tecnolog\'ia y la Innovaci\'on Francisco Jos\'e de Caldas, contrato Colciencias FP44842-087-2015. The second author was partially supported by Fondo Nacional de Financiamiento para la Ciencia, la Tecnolog\'ia y la Innovaci\'on Francisco Jos\'e de Caldas, contrato Colciencias FP 44842-157-2016. \begin{thebibliography}{99} \bibitem {Amrouche} C. Amrouche, V. Girault; \emph{Decomposition of vector spaces and application to the Stokes problem in arbitrary dimension}, Czechoslovak Math. J., 44 (1994), 109-140. \bibitem{Arada} N. Arada; \emph{A note on the regularity of flows with shear-dependent viscosity,} Nonlinear Analysis, 75 (2012), 5401-5415. \bibitem{Arada2} N. Arada, A. Sequeira; \emph{Steady flows of shear-dependent Oldroyd-B fluids around an obstacle,} J. math. fluid mech., 7 (2005), 451-483. \bibitem {daVeiga} H. Beirão da Veiga; \emph{On the global regularity of shear thinning flows in smooth domains,} J. Math. Anal. Appl., 349 (2009), 335-360. \bibitem{Boldfrini-Lorca} J. L. Boldrini, S. Lorca; \emph{Stationary Solutions for Generalized Boussinesq Models,} J. Diff. Eqns., 124 (1996), 389-406. \bibitem {Consiglieri} L. Consiglieri; \emph{Steady-state flows of thermal viscous incompressible fluids with convective-radiation effects,} Math. Models Methods Appl. Sci., 16 (2006), 2013-2027. \bibitem {Crispo} F. Crispo, C. R. Grisanti; \emph{On the existence, uniqueness and $C^{1,\gamma}(\Omega)\cap W^{2,2}(\Omega)$ regularity for a class of shear-thinning fluids}, J. Math. Fluid Mech., 10 (2008), 455-487. \bibitem {Ferreira0} L. C. F. Ferreira, E. J. Villamizar-Roa; \emph{Well-posedness and asymptotic behaviour for the convection problem in $\mathbb{R}^n$,} Nonlinearity, 19 (2006), 2169-2191. \bibitem {Ferreira1} L. C. F. Ferreira, E. J. Villamizar-Roa; \emph{On the stability problem for the Boussinesq equations in weak-Lp spaces,} Commun. Pure Appl. Anal., 9, No. 3 (2010), 667-684. \bibitem {Ferreira2} L. C. F. Ferreira, E. J. Villamizar-Roa; \emph{Existence of solutions to the convection problem in a pseudomeasure-type space,} Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 464, No. 2096 (2008), 1983-1999. \bibitem {Frehse} J. Frehse, J. M\'alek, M. Steinhauer; \emph{On analysis of steady flows of fluids with shear-dependent viscosity based on the Lipschitz truncation method}, SIAM J. Math. Anal., 34, No. 5 (2003), 1064-1083. \bibitem {Frehse97} J. Frehse, J. M\'alek, M. Steinhauer; \emph{An existence result for fluids with shear dependent viscosity-steady flows}, Nonlinear Anal., 30 (1997), 3041-3049. \bibitem {Galdi} G. P. Galdi; \emph{An introduction to the mathematical theory of the Navier Stokes Equations}, Vol I, Springer, 1994. \bibitem{Gilbarg} D. Gilbarg, N. S. Trudinger; \emph{Elliptic partial Differential Equations of Second Order}. Springer-Verlag, 2001. \bibitem {Kaplicky} P. Kaplick\'y, J. M\'alek, J. Star\'a; \emph{$C^{1,\gamma}$-solutions to a class of nonlinear fluids in two dimensions-stationary Dirichlet problem}, Zap. Nauchn. Sem. S.- Peterburg. Otdel. Mat. Inst. Steklov. (POMI), 259 (1999), 89-121. \bibitem{tpf} O. Kreml, M. Pokorn\'y; \emph{On the local strong solutions for the FENE dumbbell model}, Discrete Contin. Dyn. Syst. Ser. S, 3, No. 2 (2010), 311-324. \bibitem{Ladyz} O. A. Ladyzenskaya; \emph{Mathematical Problems in the Dynamics of a Viscous Incompressible Fluid}, 2nd rev. edn. (Nauka, 1970) [English transl. of 1st edn., The Mathematical Theory of Viscous Incompressible Flow (Gordon and Breach, 1969)]. \bibitem{Lions} J. L. Lions; \emph{Quelques m\'ethodes de r\'esolution des probl\`emes aux limites non lin\'eaires,} (Dunod et Gauthier-Villars, 1969). \bibitem{Rajagopal} K. R. Rajagopal; \emph{Mechanics of non-Newtonians fluids, in Recent Developments in Theoretical Fluid Mechanics}, eds. G. P. Galdi and J. Ne$\check{c}$as, Pitman Res. Notes in Math., Vol. 291 (Pitman, 1993), 129-162. \bibitem{Rodrigues} J. F. Rodr\'igues; \emph{Weak solutions for thermoconvective flows of Boussinesq-Stefan type}, in Mathematical Topics in Fluid Mechanics, eds. J. F. Rodrigues and A. Sequeira (Longman, 1992), 93-116. \bibitem{Ruzicka} M. R\v{z}i\v{c}ka; \emph{A note on steady flow of fluids with shear dependent viscosity}, Nonlinear Anal., 30 (1997), 3029-3039. \bibitem{Villamizar1} E. J. Villamizar-Roa, M. A. Rodr\'iguez-Bellido, M. A. Rojas-Medar; \emph{The Boussinesq system with mixed nonsmooth boundary data}, C. R. Math. Acad. Sci. Paris 343, No. 3 (2006), 191-196. \bibitem{Villamizar2} E. J. Villamizar-Roa, E. Ortega-Torres; \emph{On a generalized Boussinesq model around a rotating obstacle: existence of strong solutions.} Discrete Contin. Dyn. Syst. Ser. B, 15, No. 3 (2011), 825-847. \end{thebibliography} \end{document}