On unsteady flows of pore pressure-activated granular materials

Abstract

We investigate mathematical properties of the system of nonlinear partial differential equations that describe, under certain simplifying assumptions, evolutionary processes in water-saturated granular materials. The unconsolidated solid matrix behaves as an ideal plastic material before the activation takes place and then it starts to flow as a Newtonian or a generalized Newtonian fluid. The plastic yield stress is non-constant and depends on the difference between the given lithostatic pressure and the pressure of the fluid in a pore space. We study unsteady three-dimensional flows in an impermeable container, subject to stick-slip boundary conditions. Under realistic assumptions on the data, we establish long-time and large-data existence theory.

1 Introduction

The purpose of this study is to investigate mathematical properties of a system of nonlinear partial differential equations (PDEs) developed in [2] to describe processes, such as static liquefaction or enhanced oil recovery, in water-saturated (geological) materials. Such materials can be viewed as two component mixtures consisting of a granular solid matrix and a fluid filling the interstitial pore space. More specifically, we investigate the following system of PDEs:

$$\begin{aligned} \mathop {\mathrm {div}}\nolimits \varvec{v}&= 0, \end{aligned}$$

(1.1a)

$$\begin{aligned} \varrho ^\mathrm {m}_\mathrm {s}\left( {\partial _t\varvec{v}} + \mathop {\mathrm {div}}\nolimits {(\varvec{v}{\otimes }\varvec{v})}\right)&= \mathop {\mathrm {div}}\nolimits {\mathbb {S}}-\nabla p + \varrho ^\mathrm {m}_\mathrm {s}\varvec{b}\ , \end{aligned}$$

(1.1b)

$$\begin{aligned} {\partial _tp_\mathrm {f}} + \varvec{v}\cdot \nabla p_\mathrm {f}&= K\Delta p_\mathrm {f}- \mathop {\mathrm {div}}\nolimits (K\varrho ^\mathrm {m}_\mathrm {f}\varvec{b}) + \partial _t p_s + \varvec{v}\cdot \nabla p_s\ , \end{aligned}$$

(1.1c)

$$\begin{aligned} \varvec{v}_\mathrm {f}&= \varvec{v}- \frac{1}{\alpha }{\widehat{\phi }}(p-p_\mathrm {f})\left( \nabla p_\mathrm {f}- \varrho ^\mathrm {m}_\mathrm {f}\varvec{b}\right) , \end{aligned}$$

(1.1d)

where $\mathbb {S}$ and $\mathbb {D}\varvec{v}$ satisfy

$$\begin{aligned} \begin{aligned} \mathbb {D}\varvec{v}=\mathbb {O}&\Rightarrow \ |\mathbb {S}|\le \tau (p_\mathrm {f}),\\ \mathbb {D}\varvec{v}\ne \mathbb {O}&\Rightarrow \ \mathbb {S}=\tau (p_\mathrm {f})\frac{ \mathbb {D}\varvec{v}}{|\mathbb {D}\varvec{v}|}+2 \nu _* \left( |\mathbb {D}\varvec{v}|-\delta _*\right) ^+\frac{\mathbb {D}\varvec{v}}{|\mathbb {D}\varvec{v}|}, \end{aligned} \text{ with } \tau ({p_\mathrm {f}}) := q_*(p_s- p_\mathrm {f})^+.\end{aligned}$$

(1.1e)

The system (1.1) coincides with equations (2.24)–(2.26) stated in [2] (see also [7]) provided that we set $\delta _*=0$ in (1.1e) and we identify the symbols $\varvec{v}$ and $p_\mathrm {f}$ with $\varvec{v}_\mathrm {s}$ and $p_\mathrm {f}^\mathrm {t}$ used in [2]. In [2], Eq. (1.1) are summarized at the end of Sect. 2 as the outcome of derivation starting from the general principles of the theory of interacting continua, also using several well-motivated simplifying assumptions. In (1.1), $\varvec{v}$ represents the velocity of the granular solid matrix, $\varvec{v}_\mathrm {f}$ is the velocity of the interstitial fluid, p stands for the total pressure of the whole mixture and $p_\mathrm {f}$ is the pressure of the fluid in a pore space. The vector and scalar fields $\varvec{v}, \varvec{v}_\mathrm {f}, p$ and $p_\mathrm {f}$ represent the unknowns, the other quantities are given material functions/parameters. More precisely, $\phi ={\widehat{\phi }}(p-p_\mathrm {f})$ is the porosity given as a function of the “effective” pressure $p-p_\mathrm {f}$, $\varrho ^\mathrm {m}_\mathrm {s}$ and $\varrho ^\mathrm {m}_\mathrm {f}$ are the constant material densities of the solid and the fluid, $\varvec{b}$ represents given external forces, $\alpha$ is the drag coefficient, $\nu _*$ the viscosity of the fluid, K is a constant coefficient related to permeability, $\delta _*$ is the non-zero activation parameter and $p_s$ is the given lithostatic pressure. Note that $\varvec{v}_\mathrm {f}$ appears only in (1.1d) and can be always obtained a posteriori once $\varvec{v}, p_\mathrm {f}$ and p are obtained from (1.1a)–(1.1c). Consequently, in what follows, we consider system (1.1) without Eq. (1.1d). It is worth observing that the constitutive relation (1.1e) can be rewritten in a more compact way as an implicit constitutive relation (see Fig. 1):

$$\begin{aligned} \mathbb {S}= \mathbb {Z}+2 \nu _* \left( |\mathbb {D}\varvec{v}|-\delta _*\right) ^+\frac{\mathbb {D}\varvec{v}}{|\mathbb {D}\varvec{v}|} \text{ with } \mathbb {Z} \text{ fulfilling } (|\mathbb {Z}|-\tau (p_\mathrm {f}))^+ + ||\mathbb {D}\varvec{v}|\mathbb {Z}-\tau (p_\mathrm {f})\mathbb {D}\varvec{v}|=0. \end{aligned}$$

(1.2)

We will exploit formulation (1.2) in our analysis. A systematic study of implicit constitutive equations go back to the original works [9] and [10].

We study the system of PDEs (1.1a)–(1.1c) and (1.2) in time-space cylinder $(0, T)\times \Omega$, where $T>0$ and $\Omega \subset {\mathbb {R}}^3$ is a bounded open connected set with Lipschitz boundary $\partial \Omega$. We complete the system by considering the following boundary and initial conditions:

$$\begin{aligned}&\varvec{v}\cdot \varvec{n}=0 \text{ and } \nabla p_\mathrm {f}\cdot \varvec{n}=0 \text{ on } (0,T)\times \partial \Omega , \end{aligned}$$

(1.3a)

$$\begin{aligned}&\varvec{s}= \varvec{z}+\gamma _* \left( |\varvec{v}_\tau |-\beta _*\right) ^+ \frac{\varvec{v}_\tau }{|\varvec{v}_\tau |} \text{ with } \varvec{z} \text{ fulfilling } (|\varvec{z}|-s_*)^+ + ||\varvec{v}_\tau |\varvec{z}-s_*\varvec{v}_\tau |=0 \text{ on } (0,T)\times \partial \Omega , \end{aligned}$$

(1.3b)

$$\begin{aligned}&\varvec{v}(0, \cdot )=\varvec{v}_0 \text{ and } p_\mathrm {f}(0,\cdot )={p}_0 \text{ in } \Omega . \end{aligned}$$

(1.3c)

Here, we used the following notation: $\varvec{n}:\partial \Omega \rightarrow {\mathbb {R}}^3$ stands for the unit outer normal vector, while for any vector $\varvec{z}$ defined on $\partial \Omega$, $\varvec{z}_\tau :=\varvec{z}-(\varvec{z}\cdot \varvec{n})\varvec{n}$ denotes the tangential component of $\varvec{z}$, in particular, $\varvec{s}:=-(\mathbb {S}\varvec{n})_\tau$, and $\gamma _*, \beta _*, s_*$ are non-negative constants. Condition (1.3b) describes the shifted stick-slip (or threshold slip) and it is analogous to that for the stress tensor in the bulk (see (1.2)). It includes as special cases, the stick-slip by taking $\beta _*=0$ while $s_*, \gamma _*>0$, Navier’s slip $\varvec{s}=\gamma _*\varvec{v}_\tau$ by taking $s_*, \beta _*=0$ while $\gamma _*>0$, and perfect slip $\varvec{s}=\varvec{0}$ by setting $s_*=\gamma _*=0$. Note that the no-slip condition is obtained by letting either $s_*\rightarrow +\infty$ or by setting $\beta _*=0$ and letting $\gamma _*\rightarrow +\infty$.

Fig. 1

Representation of the material response described by (1.2)

The main purpose of this study is to establish long-time and large-data theory to the initial- and boundary-value problem described by (1.1a)–(1.1c), (1.2), (1.3), see Theorem 2.1 below. The novelties consist not only in incorporating a more general model with $\delta _*\ge 0$, but more importantly in providing a different proof for more general class of data (particularly for $\varvec{b}$ that is merely $L^2$-integrable). More precisely, we can avoid using $L^\infty$-estimates for $p_\mathrm {f}$ needed in [2]. Consequently, the main tool for taking the limit in the constitutive equations cannot be applied in the form given in [2, Proposition 5.3], but has to be modified in an essential way due to a lower integrability of $p_\mathrm {f}$, but also due to a more complicated material response.

Fig. 2

Representation of the material response described by (1.4), whenever $q>2$

The novel key tool regarding the attainment of the constitutive equations by the limiting objects is proved separately in Proposition 3.1. The key assumptions of this proposition, namely (3.5) and (3.9), call for taking $\varvec{v}^n-\varvec{v}$ as a test function in the weak formulation of balance of linear momentum. However, $\varvec{v}^n-\varvec{v}$ is not an admissible test function in the setting considered here. This difficulty can be overcome by using the $L^\infty$-truncation method which requires to introduce an integrable pressure, as the truncations $(\varvec{v}^n-\varvec{v})_\infty$ are not divergenceless. Following the approach originally developed in [6] (see also [5]), we overcome such difficulty by considering slipping boundary conditions (1.3a)–(1.3b). As pointed out in [3], the analysis for unsteady flows changes remarkably when the no–slip condition is considered.

We use the $L^\infty$-truncation method in proving Theorem 2.1 below. While the truncations $(\varvec{v}^n-\varvec{v})_\infty$ are difficult to make solenoidal, the authors of [4] succeeded to make the Lipschitz approximations $(\varvec{v}^n-\varvec{v})_{1, \infty }$ divergenceless and they thus developed a solenoidal version of the Lipschitz truncation method. This tool allows one to avoid the presence of the pressure in the setting, therefore one may include more general responses as well as boundary conditions. As a matter of fact, we present new results available for systems describing materials that behave after activation $|{\mathbb {D}}\varvec{v}|>\delta _*$, as a power-law fluid, i.e. the constitutive equation (1.2) is replaced by (see Fig. 2)

$$\begin{aligned} \mathbb {S}= \mathbb {Z}+ 2 \nu _*|\mathbb {D}\varvec{v}|^{q-2}\left( |\mathbb {D}\varvec{v}|-\delta _*\right) ^+\frac{\mathbb {D}\varvec{v}}{|\mathbb {D}\varvec{v}|} \text{ with } \mathbb {Z} \text{ fulfilling } (|\mathbb {Z}|-\tau (p_\mathrm {f}))^+ + ||\mathbb {D}\varvec{v}|\mathbb {Z}-\tau (p_\mathrm {f})\mathbb {D}\varvec{v}|=0. \end{aligned}$$

(1.4)

The available results are presented in Theorem 2.2. We are not providing the proof of these results as they can be deduced from the approach used when proving Theorem 2.1 and from the methods used recently, for example, in [3]. Note that the latter results are restricted to models (1.4) with $q>\frac{6}{5}$ (in three dimensions). Recently, another concept of dissipative solution was introduced in [1] and, its long-time and large-data existence is proved independently of what is the value of q (in particular also for $q\in [1, {6}/{5}]$). In fact in the theory developed in [1] the stress tensor can be merely subdifferential of a convex potential depending on ${\mathbb {D}}\varvec{v}$, whose growth is at least linear. There are other approaches to analyze the mathematical properties of Bingham fluids (see e.g. [8] and [11]), but they are usually based on regularity techniques requiring smoother data.

2 Preliminaries and main results

For the sake of simplicity in the given datum in the right-hand side of (1.1c), which has the form $g:=\partial _t p_s -\mathop {\mathrm {div}}\nolimits \varvec{b}$, we omit the effect of $\partial _t p_s$ as it plays the role of a given external force and it can be easily incorporated into the analysis. We also set without loss of any generality $\varrho ^\mathrm {m}_\mathrm {s}=\varrho ^\mathrm {m}_\mathrm {f}=K= 2\nu _* = \gamma _* =q_*=1,$ while we assume $\delta _*, s_*, \beta _*\ge 0$. Finally, to shorten the notation we set $Q:=(0,T)\times \Omega$ and $\Sigma := (0,T)\times \partial \Omega$, where we fix $\Omega$ to be a bounded open set in ${\mathbb {R}}^3$ with either Lipschitz or $C^{1,1}$ boundary $\partial \Omega$; such sets are denoted either $\Omega \in C^{0,1}$ or $\Omega \in C^{1,1}$ respectively.

Before stating the main results, let us summarize the notation. The symbol $\mathbb {D}\varvec{\varphi }$ stands for the symmetric part of the gradient of a vector-valued function $\varvec{\varphi }$, i.e.

$$\begin{aligned} \mathbb {D}\varvec{\varphi }:=\frac{\nabla \varvec{\varphi }+(\nabla \varvec{\varphi })^T}{2}. \end{aligned}$$

The symbols $(L^q(\Omega ), \Vert \cdot \Vert _q)$ and $(W^{1,q}(\Omega ), \Vert \cdot \Vert _{1,q})$ with $q\in [1,\infty ]$, stand respectively for the Lebesgue spaces, the Sobolev spaces with their own norms. If X is a Banach space of scalar functions, then $X^3$ and $X^{3\times 3}$ denote the space of vector-valued functions having three components and the space of tensor-valued functions respectively, with each component belonging to X. For a Banach space X, $L^q(0,T; X)$ denotes a corresponding Bochner space. We make use of the following function spaces

$$\begin{aligned}&L^{q}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits } := \overline{\left\{ \varvec{v}\in C^\infty (\Omega )^3; \, \mathop {\mathrm {div}}\nolimits \varvec{v}=0 \ \text{ in } \ \Omega ; \, \varvec{v}\cdot \varvec{n}=0 \ \text{ on } \ \partial \Omega \right\} }^{\Vert \cdot \Vert _{L^q(\Omega )^3}} \text{ for } q\in [1,\infty ),\\&\quad W_{\varvec{n}}^{1,q} := \{\varvec{v}\in W^{1,q} (\Omega )^3; \varvec{v}\cdot \varvec{n}= 0 \ \text{ on }\ \partial \Omega \},\\&\quad W_{\varvec{n},\mathop {\mathrm {div}}\nolimits }^{1,q} := \{\varvec{v}\in W^{1,q} (\Omega )^3; \, \mathop {\mathrm {div}}\nolimits \varvec{v}=0\ \text{ in }\ \Omega ; \, \varvec{v}\cdot \varvec{n}=0 \ \text{ on } \ \partial \Omega \}, \end{aligned}$$

while ${(W_{\varvec{n}}^{1,q})}^*$, ${(W_{\varvec{n},\mathop {\mathrm {div}}\nolimits }^{1,q})}^*$ are the dual spaces to $W_{\varvec{n}}^{1,q}$ and $W_{\varvec{n},\mathop {\mathrm {div}}\nolimits }^{1,q}$ respectively. In particular, assuming $\Omega \in C^{1,1}$ the following Helmholtz decomposition holds

$$\begin{aligned} W^{1,2}_{\varvec{n}}=W^{1,2}_{\varvec{n},\mathop {\mathrm {div}}\nolimits } \oplus \{\nabla \varphi ; \varphi \in W^{2, 2}(\Omega ), \nabla \varphi \cdot \varvec{n}=0 \text { on } \partial \Omega \}. \end{aligned}$$

Note that such decomposition is not valid for $W^{1,2}_{0}(\Omega )^3$.

We are ready to enunciate the first result, which is the existence of weak solutions to system (1.1a)–(1.1c), (1.2), (1.3), proved in Sect. 5.

Theorem 2.1

For any $\Omega \in C^{1,1}$, $T > 0$ and for any $\varvec{v}_0, p{_0}, \varvec{b}, p_s$ fulfilling

$$\begin{aligned} \varvec{v}_0\in L^{2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }, \ p{_0}\in L^2(\Omega ), \ \varvec{b}\in L^2(Q), \ p_s\in L^5(Q), \end{aligned}$$

there exists a quintuplet $(\varvec{v},p_\mathrm {f}, p, \mathbb {S}, \varvec{s})$:

$$\begin{aligned}&\varvec{v}\in L^\infty (0,T; L^{2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits })\cap L^2(0,T;W^{1,2}_{{\varvec{n}},\mathrm{div}}), \partial _t\varvec{v}\in {(L^2(0,T;W^{1,2}_{\varvec{n}})\cap L^5(Q)^3)}^*,\\&\quad p_\mathrm {f}\in L^\infty (0,T; L^2(\Omega ))\cap L^2(0,T; W^{1,2}(\Omega )), \ \partial _t p_\mathrm {f}\in (L^4(0,T; W^{1,2}(\Omega )))^*,\\&\quad p=p_1+p_2 \text{ where } p_1\in L^2(Q) \text{ and } p_2 \in L^{\frac{5}{4}}(0,T; W^{1, \frac{5}{4}}(\Omega )), \\&\quad \mathbb {S}\in L^{2}(Q)^{3\times 3}, \ \varvec{s}\in L^{\frac{8}{3}}(\Sigma )^3, \end{aligned}$$

satisfying the following weak formulations:

$$\begin{aligned}&\int \limits _0^T\langle \partial _t\varvec{v}, \varvec{w}\rangle + \int \limits _{Q}\mathbb {S}:\mathbb {D}\varvec{w}-\int \limits _{Q} (\varvec{v}\otimes \varvec{v}): \mathbb {D}\varvec{w}+\int \limits _{\Sigma } \varvec{s}\cdot \varvec{w}_{\varvec{\tau }} = \int \limits _{Q}\varvec{b}\cdot \varvec{w}+ \int \limits _{Q} p_1\mathop {\mathrm {div}}\nolimits \varvec{w}- \int \limits _{Q} \nabla p_2\cdot \varvec{w}\nonumber \\&\quad \text{ for } \text{ all } \varvec{w}\in L^2(0,T;W^{1,2}_{\varvec{n}})\cap L^5(Q)^3, \end{aligned}$$

(2.1)

$$\begin{aligned}&\int \limits _0^T \langle \partial _t p_\mathrm {f}, z\rangle - \int \limits _{Q}p_\mathrm {f}\varvec{v}\cdot \nabla z + \int \limits _{Q}\nabla p_\mathrm {f}\cdot \nabla z = \int \limits _{Q}\varvec{b}\cdot \nabla z - \int \limits _{Q} p_s\varvec{v}\cdot \nabla z \ \text{ for } \text{ all } z\in L^4(0,T; W^{1,2}(\Omega )), \end{aligned}$$

(2.2)

and the following constitutive equations:

$$\begin{aligned}&\mathbb {S}= \mathbb {Z}+ \left( |\mathbb {D}\varvec{v}|-\delta _*\right) ^+\frac{\mathbb {D}\varvec{v}}{|\mathbb {D}\varvec{v}|} \text{ with } \mathbb {Z} \text{ fulfilling } (|\mathbb {Z}|-\tau (p_\mathrm {f}))^+ + ||\mathbb {D}\varvec{v}|\mathbb {Z}-\tau (p_\mathrm {f})\mathbb {D}\varvec{v}|=0\nonumber \\&\quad \text{ with } \tau ({p_\mathrm {f}}) = (p_s- p_\mathrm {f})^+ \text{ a.e. } \text{ in } Q, \end{aligned}$$

(2.3)

$$\begin{aligned}&\varvec{s}= \varvec{z}+ \left( |\varvec{v}_\tau |-\beta _*\right) ^+ \frac{\varvec{v}_\tau }{|\varvec{v}_\tau |} \text{ with } \varvec{z} \text{ fulfilling } (|\varvec{z}|-s_*)^+ + ||\varvec{v}_\tau |\varvec{z}-s_*\varvec{v}_\tau |=0 \text{ a.e. } \text{ on } \Sigma , \end{aligned}$$

(2.4)

and attaining the initial conditions in the following sense:

$$\begin{aligned} \lim _{t\rightarrow 0+}\Vert \varvec{v}(t)-\varvec{v}_0\Vert _2=0, \ \ \lim _{t\rightarrow 0+}\Vert p_\mathrm {f}(t)-p_0\Vert _2=0. \end{aligned}$$

(2.5)

The second result concerns system (1.1a)–(1.1c), (1.3), and (1.4).

Theorem 2.2

Let $\Omega \in C^{0,1}$, $T > 0$, and $q>\frac{6}{5}$. Set $m:=\max \{2, q'\}$ and $r:=\max \left\{ q, \frac{5q}{5q-6}\right\} .$ For any $\varvec{v}_0, p{_0}, \varvec{b}, p_s$ fulfilling

$$\begin{aligned} \varvec{v}_0\in L^{2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }, \ p{_0}\in L^2(\Omega ), \ \varvec{b}\in L^{m}(Q), \ p_s\in L^{\frac{10q}{5q-6}}(Q), \end{aligned}$$

there exists a quadruplet $(\varvec{v},p_\mathrm {f}, \mathbb {S}, \varvec{s})$:

$$\begin{aligned}&\varvec{v}\in L^\infty (0,T; L^{2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }) \cap L^q(0,T;W^{1,q}_{\mathbf{n},\mathrm{div}}), \partial _t\varvec{v}\in L^{r'}(0, T; (W^{1,r}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits })^*),\\&\quad p_\mathrm {f}\in L^\infty (0,T; L^2(\Omega ))\cap L^2(0,T; W^{1,2}(\Omega )), \partial _t p_\mathrm {f} \in (L^2(0,T; W^{1, 2}(\Omega ))\cap L^{\frac{10q}{5q-6}}(Q))^*,\\&\quad \mathbb {S}\in L^{q'}(Q)^{3\times 3}, \ \varvec{s}\in L^{2}(\Sigma )^3, \end{aligned}$$

satisfying the initial conditions (2.5), the following weak formulations:

$$\begin{aligned}&\int \limits _0^T\langle \partial _t\varvec{v}, \varvec{\varphi }\rangle + \int \limits _{Q}\mathbb {S}: \mathbb {D}\varvec{\varphi }- \int \limits _{Q}(\varvec{v}\otimes \varvec{v}): \mathbb {D}\varvec{\varphi }+\int \limits _{\Sigma } \varvec{s}\cdot \varvec{\varphi }_{\varvec{\tau }}= \int \limits _{Q}\varvec{b}\cdot \varvec{\varphi }\ \ \ \text { for all } \varvec{\varphi }\in L^r(0,T; W^{1,r}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }), \end{aligned}$$

(2.6)

$$\begin{aligned}&\int \limits _0^T\langle \partial _tp_\mathrm {f}, \varphi \rangle +\int \limits _{Q}\nabla p_\mathrm {f}\cdot \nabla \varphi - \int \limits _{Q} p_\mathrm {f}\varvec{v}\cdot \nabla \varphi =\int \limits _{Q}(\varvec{b}-p_s\varvec{v})\cdot \nabla \varphi \ \ \ \text{ for } \text{ all } \varphi \in L^2(0,T; W^{1, 2}(\Omega ))\cap L^{\frac{10q}{5q-6}}(Q),\nonumber \\ \end{aligned}$$

(2.7)

and the following constitutive equations:

$$\begin{aligned}&\mathbb {S}= \mathbb {Z}+ \left( |\mathbb {D}\varvec{v}|-\delta _*\right) ^+|\mathbb {D}\varvec{v}|^{q-2}\frac{\mathbb {D}\varvec{v}}{|\mathbb {D}\varvec{v}|} \text{ with } \mathbb {Z} \text{ fulfilling } (|\mathbb {Z}|-\tau (p_\mathrm {f}))^+ + ||\mathbb {D}\varvec{v}|\mathbb {Z}-\tau (p_\mathrm {f})\mathbb {D}\varvec{v}|=0\nonumber \\&\quad \text{ with } \tau ({p_\mathrm {f}}) = (p_s- p_\mathrm {f})^+ \text { a.e. in } Q, \end{aligned}$$

(2.8)

$$\begin{aligned}&\varvec{s}= \varvec{z}+ \left( |\varvec{v}_\tau |-\beta _*\right) ^+ \frac{\varvec{v}_\tau }{|\varvec{v}_\tau |} \text{ with } \varvec{z} \text{ fulfilling } (|\varvec{z}|-s_*)^+ + ||\varvec{v}_\tau |\varvec{z}-s_*\varvec{v}_\tau |=0 \text{ a.e. } \text{ on } \Sigma . \end{aligned}$$

(2.9)

This result is stated without the proof here. The proof can be however achieved in the spirit of Theorem 2.1 by employing a solenoidal version of the Lipschitz-truncation method developed in [4], and by using the approximation scheme presented in [3, Theorem 3.3].

3 Attainment of the constitutive equations

In this section, we establish a new scheme how to take the limit in the constitutive equations needed when proving Theorem 2.1.

Proposition 3.1

Let $U\subset Q$ be an arbitrary measurable bounded set and let $\{\mathbb {Z}^n\}_{n=1}^{+\infty }$, $\{\mathbb {D}^n\}_{n=1}^{+\infty }$ and $\{p_\mathrm {f}^n\}_{n=1}^{+\infty }$ be sequences such that

$$\begin{aligned}&\mathbb {Z}^n=\tau (p_\mathrm {f}^n) \frac{\mathbb {D}^n}{|\mathbb {D}^n|+\frac{1}{n}} \text{ with } \tau (p_\mathrm {f}^n)= (p_s-p_\mathrm {f}^n)^+ \text{ a.e. } \text{ in } U, \end{aligned}$$

(3.1)

$$\begin{aligned}&\mathbb {Z}^n \rightharpoonup \mathbb {Z}\ \text {weakly in } L^{2}(U)^{3\times 3}, \end{aligned}$$

(3.2)

$$\begin{aligned}&\mathbb {D}^n \rightharpoonup \mathbb {D}\ \text {weakly in } L^{2}(U)^{3\times 3}, \end{aligned}$$

(3.3)

$$\begin{aligned}&p_\mathrm {f}^n \rightarrow p_\mathrm {f}\text { strongly in } L^2(U) \text{ and } \text{ a.e. } \text{ in } U, \end{aligned}$$

(3.4)

$$\begin{aligned}&\limsup _{n\rightarrow \infty } \int \limits _{U} \mathbb {Z}^n : \mathbb {D}^n \le \int \limits _{U} \mathbb {Z}: \mathbb {D}. \end{aligned}$$

(3.5)

Then

$$\begin{aligned} (|\mathbb {Z}|-\tau (p_\mathrm {f}))^+ + ||\mathbb {D}|\mathbb {Z}-\tau (p_\mathrm {f})\mathbb {D}|=0. \end{aligned}$$

(3.6)

In addition, assume that $\{\mathbb {V}^n\}_{n=1}^{+\infty }$ is a sequence such that

$$\begin{aligned} \mathbb {V}^n=\left( 1-\frac{\delta _*}{|\mathbb {D}^n|}\right) ^+\mathbb {D}^n \ \ \text{ a.e. } \text{ in } U, \end{aligned}$$

(3.7)

fulfilling

$$\begin{aligned}&\mathbb {V}^n \rightharpoonup \mathbb {V}\ \text {weakly in } L^{2}(U)^{3\times 3}, \end{aligned}$$

(3.8)

$$\begin{aligned}&\limsup _{n\rightarrow \infty } \int \limits _{U} \mathbb {V}^n : \mathbb {D}^n \le \int \limits _{U} \mathbb {V}: \mathbb {D}. \end{aligned}$$

(3.9)

Then

$$\begin{aligned} \mathbb {V}= \left( 1-\frac{\delta _*}{|\mathbb {D}|}\right) ^+\mathbb {D}\ \ \text{ a.e. } \text{ in } U. \end{aligned}$$

(3.10)

Proof

First, note that by virtue of (3.4) and the Lipschitz-continuity of $\tau$, it follows that

$$\begin{aligned} \tau (p_\mathrm {f}^n) \rightarrow \tau (p_\mathrm {f}) \text{ strongly } \text{ in } L^2(U). \end{aligned}$$

(3.11)

Now, as for any $\mathbb {A}\in L^2 (U)$, $\mathbb {A}\ne {\mathbb {O}}$, it holds

$$\begin{aligned} \left( \tau (p_\mathrm {f}) \frac{\mathbb {D}^n}{|\mathbb {D}^n|+ \frac{1}{n}} - \tau (p_\mathrm {f}) \frac{\mathbb {A}}{|\mathbb {A}|+\frac{1}{n}}\right) : \left( \mathbb {D}^n -\mathbb {A}\right) \ge 0, \end{aligned}$$

(3.12)

integrating this inequality over U, subtracting and adding $\mathbb {Z}^n$ and using (3.1), we get

$$\begin{aligned}&\int \limits _{U} \left( \tau (p_\mathrm {f}) \frac{\mathbb {D}^n}{|\mathbb {D}^n|+\frac{1}{n}} -\tau (p_\mathrm {f}^n) \frac{\mathbb {D}^n}{|\mathbb {D}^n|+\frac{1}{n}}\right) : \left( \mathbb {D}^n -\mathbb {A}\right) \nonumber \\&\quad + \int \limits _{U}\left( \mathbb {Z}^n - \tau (p_\mathrm {f}) \frac{\mathbb {A}}{|\mathbb {A}|+\frac{1}{n}}\right) :\left( \mathbb {D}^n - \mathbb {A}\right) \ge 0. \end{aligned}$$

(3.13)

Taking limsup as $n\rightarrow \infty$ and employing the facts that the first integral converges to zero due to (3.11), because the term $\frac{\mathbb {D}^n}{|\mathbb {D}^n|+\frac{1}{n}}$ is uniformly bounded in $L^\infty (U)$ and the sequence $\mathbb {D}^n -\mathbb {A}$ is bounded in $L^2 (U)$, we obtain

$$\begin{aligned} \limsup _{n\rightarrow \infty }\int \limits _{U} \left( \mathbb {Z}^n - \tau (p_\mathrm {f}) \frac{\mathbb {A}}{|\mathbb {A}|+\frac{1}{n}}\right) :\left( \mathbb {D}^n - \mathbb {A}\right) \ge 0. \end{aligned}$$

(3.14)

Referring then to the convergences (3.2) and (3.3) and using also (3.5), we conclude that

$$\begin{aligned} \int \limits _{U} \left( \mathbb {Z}- \tau (p_\mathrm {f}) \frac{\mathbb {A}}{|\mathbb {A}|}\right) :\left( \mathbb {D}- \mathbb {A}\right) \ge 0. \end{aligned}$$

(3.15)

Now, for any $\delta >0$, $\varepsilon \in (0, \delta )$ and for arbitrary matrices $\mathbb {C}$ and $\mathbb {B}_1$ bounded in $L^2(U)^{3\times 3}$ and satisfying $|\mathbb {C}|\le 1$ and $\mathbb {B}_1 \ne \mathbb {O}$, consider

$$\begin{aligned} \mathbb {A}:= \mathbb {B}_1 \, \chi _{\{|\mathbb {D}|=0\}} + (\mathbb {D}- \varepsilon \mathbb {C}) \, \chi _{\{|\mathbb {D}|>\delta \}} + \mathbb {D}\, \chi _{\{0<|\mathbb {D}|\le \delta \}}. \end{aligned}$$

Note that such $\mathbb {A}$’s are non-zero in U. Inserting them into (3.15) we obtain

$$\begin{aligned} - \int \limits _{\{|\mathbb {D}|=0\}} \left( \mathbb {Z}- \tau (p_\mathrm {f}) \frac{\mathbb {B}_1}{|\mathbb {B}_1|}\right) :\mathbb {B}_1 + \varepsilon \int \limits _{\{|\mathbb {D}|>\delta \}} \mathbb {C}: \left( \mathbb {Z}- \tau (p_\mathrm {f}) \frac{\mathbb {D}- \varepsilon \mathbb {C}}{|\mathbb {D}- \varepsilon \mathbb {C}|}\right) \ge 0. \end{aligned}$$

(3.16)

Letting first $\varepsilon \rightarrow 0$ in (3.16), we observe that

$$\begin{aligned} \int \limits _{\{|\mathbb {D}|=0\}} \mathbb {Z}: \mathbb {B}_1 \le \int \limits _{\{|\mathbb {D}|=0\}} \tau (p_f) |\mathbb {B}_1| \end{aligned}$$

(3.17)

for any $\mathbb {B}_1\ne {\mathbb {O}}$. Consider, for any $a>0$ and $\omega \subset U$, the matrix $\mathbb {B}_1$ of the form

$$\begin{aligned} \mathbb {B}_1=a\, {\mathbb {I}}\chi _{\{(U\setminus \omega ) \cup \{\mathbb {Z}=0\}\}}+ \frac{\mathbb {Z}}{|\mathbb {Z}|}\chi _{\{\omega \setminus \{\mathbb {Z}=0\}\}}. \end{aligned}$$

It then follows from (3.17) that

$$\begin{aligned} \int \limits _{\{|\mathbb {D}|=0\}\cap \, \omega \,\cap \{\mathbb {Z}\ne 0\}} |\mathbb {Z}| \le \int \limits _{\{|\mathbb {D}|=0\}\cap \, \omega \cap \{\mathbb {Z}\ne 0\}} \tau (p_\mathrm {f}) + a\, C\int \limits _{(U\setminus \omega ) \,\cup \{\mathbb {Z}=0\}}(\tau (p_\mathrm {f})+ |\mathbb {Z}|) \end{aligned}$$

with C positive constant, which implies letting $a\rightarrow 0$

$$\begin{aligned} \int \limits _{\{|\mathbb {D}|=0\}\cap \, \omega \,\cap \{\mathbb {Z}\ne 0\}} |\mathbb {Z}| \le \int \limits _{\{|\mathbb {D}|=0\}\cap \, \omega \cap \{\mathbb {Z}\ne 0\}} \tau (p_\mathrm {f}). \end{aligned}$$

Since $\omega$ is arbitrary, we conclude that

$$\begin{aligned} |\mathbb {Z}|\le \tau (p_\mathrm {f}) \text{ on } \text{ the } \text{ set } \{|\mathbb {D}|=0\}. \end{aligned}$$

(3.18)

Next, letting $|\mathbb {B}_1|\rightarrow 0$ in (3.16), employing (3.17), we get

$$\begin{aligned} \int \limits _{\{|\mathbb {D}|>\delta \}} \mathbb {C}: \left( \mathbb {Z}- \tau (p_\mathrm {f}) \frac{\mathbb {D}- \varepsilon \mathbb {C}}{|\mathbb {D}- \varepsilon \mathbb {C}|}\right) \ge 0, \end{aligned}$$

which, after letting $\varepsilon \rightarrow 0$, leads to

$$\begin{aligned} \int \limits _{\{|\mathbb {D}|>\delta \}} \mathbb {C}: \left( \mathbb {Z}- \tau (p_\mathrm {f}) \frac{\mathbb {D}}{|\mathbb {D}|}\right) \ge 0. \end{aligned}$$

Finally, letting $\delta \rightarrow 0$, we get, for arbitrary $\mathbb {C}$,

$$\begin{aligned} \int \limits _{\{|\mathbb {D}|>0\}} \mathbb {C}: \left( \mathbb {Z}-\tau (p_\mathrm {f}) \frac{\mathbb {D}}{|\mathbb {D}|}\right) \ge 0. \end{aligned}$$

This implies

$$\begin{aligned} \mathbb {Z}=\tau (p_\mathrm {f}) \frac{\mathbb {D}}{|\mathbb {D}|} \text{ when } |\mathbb {D}|\ne 0. \end{aligned}$$

(3.19)

The latter and (3.18) are equivalent to (3.6).

It remains to prove (3.10), which however follows from standard Minty’s argument. Indeed, by the monotonicity, we have

$$\begin{aligned} \limsup _{n\rightarrow \infty }\int \limits _{U}\left( \mathbb {V}^n-\mathbb {A}\left( 1-\frac{\delta _*}{|\mathbb {A}|}\right) ^+\right) :\left( \mathbb {D}^n-\mathbb {A}\right) \ge 0 \end{aligned}$$

for any $\mathbb {A}\in L^2(U)^{3\times 3}$. By virtue of (3.9) and of convergences (3.8) and (3.3) we get

$$\begin{aligned} \int \limits _{U}\left( \mathbb {V}-\mathbb {A}\left( 1-\frac{\delta _*}{|\mathbb {A}|}\right) ^+\right) :\left( \mathbb {D}-\mathbb {A}\right) \ge 0. \end{aligned}$$

Choosing $\mathbb {A}:= \mathbb {D}\pm \varepsilon \mathbb {C}$, with arbitrary $\mathbb {C}\in L^2 (U)^{3\times 3}$ and $\varepsilon >0$, and after the limit as $\varepsilon \rightarrow 0$ we obtain

$$\begin{aligned} \int \limits _U \mathbb {C}: \left( \mathbb {V}-\mathbb {D}\left( 1-\frac{\delta _*}{|\mathbb {D}|}\right) ^+\right) =0 \end{aligned}$$

for any $\mathbb {C}$, which implies (3.10). $\square$

Note that here we provided a proof of (3.6), which is simplified and shorter than the one given in [2, Proposition 5.3].

4 Approximations

In this section, we prepare all the needed tools in order to prove Theorem 2.1. For any $n\in {\mathbb {N}}$, we introduce the following approximating system

$$\begin{aligned} \begin{array}{l} \mathop {\mathrm {div}}\nolimits \varvec{v}= 0 \text { in } Q, \\ \partial _t\varvec{v}+ \mathop {\mathrm {div}}\nolimits (\varvec{v}\otimes \varvec{v}) G_n(|\varvec{v}|^2) -\mathop {\mathrm {div}}\nolimits \mathbb {S}+\nabla p = \varvec{b}\text { in } Q, \\ \partial _t p_\mathrm {f}+ \varvec{v}\cdot \nabla p_\mathrm {f}- \Delta p_\mathrm {f}=-\mathop {\mathrm {div}}\nolimits \varvec{b}+\varvec{v}\cdot \nabla p_s \text { in } Q, \end{array} \end{aligned}$$

(4.1)

where $G_n:{\mathbb {R}}\rightarrow {\mathbb {R}}$ is a smooth function such that $G_n(u)=1$ if $|u|\le n$, $G_n(u)=0$ if $|u|\ge 2n$ and $|G'_n|\le \frac{2}{n}$. Next, we consider the following regularization of the constitutive equations (both in the bulk and on the boundary)

$$\begin{aligned}&\begin{array}{l} \mathbb {S}={\mathcal {S}}_n(p_\mathrm {f}, \mathbb {D}{\varvec{v}})= {\mathcal {Z}}_n(p_\mathrm {f}, \mathbb {D}{\varvec{v}})+\left( 1-\frac{\delta _*}{|\mathbb {D}{\varvec{v}}|}\right) ^+\mathbb {D}{\varvec{v}} \ \ \ \text{ where } {\mathcal {Z}}_n(p_\mathrm {f}, \mathbb {D}{\varvec{v}}) := \tau (p_\mathrm {f})\frac{\mathbb {D}{\varvec{v}}}{|\mathbb {D}{\varvec{v}}|+\frac{1}{n}} \\ \text{ with } \tau (p_\mathrm {f})=(p_s-p_\mathrm {f})^+ \text { in } Q, \end{array} \end{aligned}$$

(4.2)

$$\begin{aligned}&\varvec{s}= {s}_n(\varvec{v}_\tau )=\zeta _n(\varvec{v}_\tau ) +\left( 1-\frac{\beta _*}{|\varvec{v}_\tau |}\right) ^+\varvec{v}_\tau \ \ \ \text{ where } \zeta _n(\varvec{v}_\tau )= s_* \frac{\varvec{v}_{\tau }}{|\varvec{v}_{\tau }|+\frac{1}{n}} \text{ on } \Sigma , \end{aligned}$$

(4.3)

and we complete the problem with boundary and initial conditions

$$\begin{aligned} \varvec{v}\cdot \varvec{n}&=0, \nabla p_\mathrm {f}\cdot \varvec{n}=0 \text{ on } \Sigma , \end{aligned}$$

(4.4)

$$\begin{aligned} \varvec{v}(0)&=\varvec{v}_0, \ \, p_\mathrm {f}(0)={p}_0 \text{ in } \Omega . \end{aligned}$$

(4.5)

Note that both mappings $\mathbb {D}\mapsto {\mathcal {Z}}_n(p_\mathrm {f}, \mathbb {D})$ and $\mathbb {D}\mapsto \left( 1-\frac{\delta _*}{|\mathbb {D}|}\right) ^+\mathbb {D}$ are monotone, i.e.

$$\begin{aligned} (\mathbb {Z}-{\hat{\mathbb {Z}}}):(\mathbb {D}-{\hat{\mathbb {D}}})\ge 0 \text{ for } \text{ any } \mathbb {Z}={\mathcal {Z}}_n(p_\mathrm {f}, \mathbb {D}), {\hat{\mathbb {Z}}}={\mathcal {Z}}_n(p_\mathrm {f}, {\hat{\mathbb {D}}}), \end{aligned}$$

(4.6)

see formula (5.2) in [2], and

$$\begin{aligned} (\mathbb {V}-{\hat{\mathbb {V}}}):(\mathbb {D}-{\hat{\mathbb {D}}})\ge 0 \text{ for } \text{ any } \mathbb {V}=\left( 1-\frac{\delta _*}{|\mathbb {D}|}\right) ^+\mathbb {D}, {\hat{\mathbb {V}}}=\left( 1-\frac{\delta _*}{|{\hat{\mathbb {D}}}|}\right) ^+{\hat{\mathbb {D}}}, \end{aligned}$$

(4.7)

see Lemma B.1 in [3]. Therefore, due to the presence of the truncation in the convective term and the introduced approximations in the constitutive equations, the existence of weak solutions to system (4.1)–(4.5) can be proved through standard techniques of monotone operators, following also the spirit of the proof in [2, Proposition 5.1]. We enunciate the relevant result below and for the reader’s convenience the proof can be found in Appendix.

Proposition 4.1

Let $n\in {\mathbb {N}}$ be fixed. For any

$$\begin{aligned} \varvec{v}_0\in L^{2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }, \ p_0\in L^2(\Omega ),\ \varvec{b}\in L^2(Q) \text{ and } p_s\in L^5(Q), \end{aligned}$$

there exists a weak solution to the problem (4.1)–(4.5). More precisely, for each $n\in {\mathbb {N}}$ there is a quadruplet $(\varvec{v},p_\mathrm {f}, \mathbb {S}, \varvec{s}):= (\varvec{v}^n, p_\mathrm {f}^n, \mathbb {S}^n, \varvec{s}^n)$ such that

$$\begin{aligned}&\varvec{v}\in L^\infty (0,T; L^{2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits })\cap L^2(0,T;W^{1,2}_{{\varvec{n}},\mathrm{div}}), \ \ \partial _t\varvec{v}\in {(L^2(0,T;W^{1,2}_{\varvec{n}}))}^*, \end{aligned}$$

(4.8)

$$\begin{aligned}&p_\mathrm {f}\in L^\infty (0,T; L^2(\Omega ))\cap L^2(0,T; W^{1,2}(\Omega )),\ \ \partial _t p_\mathrm {f}\in (L^4(0,T; W^{1,2}(\Omega )))^*, \end{aligned}$$

(4.9)

$$\begin{aligned}&\mathbb {S}\in L^{2}(Q)^{3\times 3}, \ \ \varvec{s}\in L^\frac{8}{3}(\Sigma )^3, \end{aligned}$$

(4.10)

satisfying

$$\begin{aligned}&\begin{array}{l} \displaystyle \int \limits _0^T \langle \partial _t\varvec{v}, \varvec{w}\rangle +\displaystyle \int \limits _{Q}(\mathbb {S}: \mathbb {D}\varvec{w}+ G_n(|\varvec{v}|^2)\mathop {\mathrm {div}}\nolimits (\varvec{v}\otimes \varvec{v}) \cdot \varvec{w}) + \int \limits _{\Sigma } \varvec{s}\cdot \varvec{w}_{\varvec{\tau }} = \int \limits _{Q} \varvec{b}\cdot \varvec{w}\text { for all } \varvec{w}\in L^2(0, T; W^{1,2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }),\end{array} \end{aligned}$$

(4.11)

$$\begin{aligned}&\int \limits _0^T \langle \partial _t p_\mathrm {f}, z\rangle -\int \limits _{Q}p_\mathrm {f}\varvec{v}\cdot \nabla z +\int \limits _{Q} \nabla p_\mathrm {f}\cdot \nabla z= \int \limits _{Q}(\varvec{b}\cdot \nabla z- p_s\varvec{v}\cdot \nabla z) \ \ \text{ for } \text{ all } z\in L^4(0,T; W^{1,2}(\Omega )), \end{aligned}$$

(4.12)

where

$$\begin{aligned}&\mathbb {S}={\mathcal {S}}_n(p_\mathrm {f}, \mathbb {D}{\varvec{v}}) \text { a.e. in } Q, \end{aligned}$$

(4.13)

$$\begin{aligned}&\varvec{s}= s_n(\varvec{v}_\tau ) \text { a.e. in } \Sigma , \end{aligned}$$

(4.14)

and

$$\begin{aligned} \lim _{t\rightarrow 0^+}\Vert \varvec{v}(t)-\varvec{v}_0\Vert _2=0, \ \lim _{t\rightarrow 0^+}\Vert p_\mathrm {f}(t)-p{_0}\Vert _2=0. \end{aligned}$$

(4.15)

5 Proof of Theorem 2.1

The proof is split in the following steps.

Step 1. Approximations. From Proposition 4.1 and following the reconstruction of the pressure in [2, Theorem 4.1, Step 2 of the proof], we get for each $n\in {\mathbb {N}}$ the existence of $(\varvec{v}^n,p_\mathrm {f}^n, p^n, {\mathbb {S}}^n, \varvec{s}^n)$, with $p^n\in L^2(Q)$, satisfying

$$\begin{aligned}&\begin{array}{l} \displaystyle \int \limits _0^T \langle \partial _t \varvec{v}^n, \varvec{w}\rangle + \int \limits _{Q}(\mathbb {S}^n:\mathbb {D}\varvec{w}+\mathop {\mathrm {div}}\nolimits (\varvec{v}^n\otimes \varvec{v}^n) G_n(|\varvec{v}^n|^2)\cdot \varvec{w}) +\displaystyle \int \limits _{\Sigma }\varvec{s}^n\cdot \varvec{w}_{\varvec{\tau }} = \int \limits _{Q}p^n \mathop {\mathrm {div}}\nolimits \varvec{w}\\ \,\,\,+ \displaystyle \int \limits _{Q} \varvec{b}\cdot \varvec{w}\ \ \text { for all } \varvec{w}\in L^2(0, T; W^{1,2}_{\varvec{n}}), \end{array} \end{aligned}$$

(5.1)

$$\begin{aligned}&\begin{array}{l} \displaystyle \int \limits _0^T\langle \partial _t p_\mathrm {f}^{n}, z\rangle - \displaystyle \int \limits _{Q} (p_\mathrm {f}^{n} \varvec{v}^{n}) \cdot \nabla z +\int \limits _{Q} \nabla p_\mathrm {f}^n\cdot \nabla z =\int \limits _{Q}(\varvec{b}- p_s \varvec{v}^n) \cdot \nabla z\ \ \text{ for } \text{ all } z \in L^4(0, T; W^{1,2}(\Omega )), \end{array} \end{aligned}$$

(5.2)

with $\mathbb {S}^n, \varvec{s}^n$ fulfilling (4.13), (4.14) respectively, and satisfying also (4.15).

Step 2. Uniform estimates with respect to n and limit as $n\rightarrow +\infty$. Setting $\varvec{w}:=\varvec{v}^n$ in (5.1) and $z:=p_\mathrm {f}^n$ in (5.2), following the analogous step as in the proof of [2, Theorem 4.1], we obtain

$$\begin{aligned}&\sup _{n}\left( \Vert \varvec{v}^n\Vert _{L^\infty (0,T; L^2(\Omega )^3)} + \Vert \mathbb {D}\varvec{v}^n\Vert _{2, Q} \right) <+\infty , \end{aligned}$$

(5.3)

$$\begin{aligned}&\sup _n \left( \Vert p_\mathrm {f}^n\Vert _{L^\infty (0,T; L^2(\Omega ))} + \Vert \nabla p_\mathrm {f}^n \Vert _{2, Q}\right) <+\infty , \end{aligned}$$

(5.4)

$$\begin{aligned}&\sup _n\left( \Vert \varvec{v}^n\Vert _{\frac{10}{3}, Q}+\Vert p_\mathrm {f}^n\Vert _{\frac{10}{3}, Q}+\Vert \varvec{v}^n\Vert _{\frac{8}{3}, \Sigma }\right) <+\infty , \end{aligned}$$

(5.5)

$$\begin{aligned}&\sup _n\left( \Vert \mathbb {Z}^n\Vert _{\frac{10}{3}, Q}+\Vert \mathbb {V}^n\Vert _{2, Q}+\Vert \varvec{s}^n\Vert _{\frac{8}{3},\Sigma }\right) <+\infty , \end{aligned}$$

(5.6)

where we set $\mathbb {V}^n:=\left( 1-\frac{\delta _*}{|\mathbb {D}{\varvec{v}^n}|}\right) ^+\mathbb {D}{\varvec{v}^n}$. Consequently, as $\sup _n \Vert G_n(|\varvec{v}^n|^2)\Vert _{L^\infty (Q)}\le 1$ by employing (5.3), (5.5) and Korn’s inequality, it follows that

$$\begin{aligned} \sup _n \Vert G_n(|\varvec{v}^n|^2)\mathop {\mathrm {div}}\nolimits (\varvec{v}^n\otimes \varvec{v}^n)\Vert _{L^{\frac{5}{4}}(Q)}<+\infty . \end{aligned}$$

(5.7)

Now, let us introduce

$$\begin{aligned} p_2^n:= (-\Delta _N )^{-1} \left( G^n(|\varvec{v}^n|^2)\mathop {\mathrm {div}}\nolimits (\varvec{v}^n \otimes \varvec{v}^n) \right) ,\ \ \ p_1^n:= p^n- p_2^n, \end{aligned}$$

then

$$\begin{aligned} \sup _n \left( \Vert p_2^n\Vert _{L^\frac{5}{4} (0,T; W^{1,\frac{5}{4}} (\Omega ))}+ \Vert p_1^n\Vert _{L^2(Q)}\right) < +\infty , \end{aligned}$$

(5.8)

and this implies that

$$\begin{aligned} \sup _n\Vert \partial _t \varvec{v}^n\Vert _{{(L^2(0,T;W^{1,2}_{\varvec{n}})\cap L^5(Q)^3)}^*}<+\infty . \end{aligned}$$

(5.9)

Analogously

$$\begin{aligned} \sup _n\Vert \partial _t p_\mathrm {f}^n\Vert _{{(L^4(0,T; W^{1,2}(\Omega )))}^*}<+\infty . \end{aligned}$$

(5.10)

Then, there exist subsequences of $\{\varvec{v}^n\}$, $\{p_\mathrm {f}^n\}$, $\{\mathbb {Z}^n\}$, $\{\mathbb {V}^n\}$, $\{\varvec{s}^n\}$, $\{p^n_1\}$, $\{p_2^n\}$, which we do not relabel, that converge weakly or *-weakly in the corresponding function spaces. By virtue of the established limits, by the Aubin–Lions compactness lemma and the compact embedding of the Sobolev spaces into the space of traces, we also have

$$\begin{aligned}&p_\mathrm {f}^n \rightarrow p_\mathrm {f}\text { strongly in } L^q(Q) \text{ for } \text{ all } q\in \left[ 1, \frac{10}{3}\right) , \end{aligned}$$

(5.11)

$$\begin{aligned}&\varvec{v}^n \rightarrow \varvec{v}\text { strongly in } L^q(Q)^3 \text{ for } \text{ all } q\in \left[ 1, \frac{10}{3}\right) , \end{aligned}$$

(5.12)

$$\begin{aligned}&\varvec{v}^n_\tau \rightarrow \varvec{v}_\tau \text { strongly in } L^q(\Sigma )^3 \text{ for } \text{ all } q\in \left[ 1, \frac{8}{3}\right) . \end{aligned}$$

(5.13)

Since

$$\begin{aligned} \Vert G_n(|\varvec{v}^n|^2)\Vert _{L^\infty (Q)}\le 1 \text{ and } G_n(|\varvec{v}^n|^2)\rightarrow 1 \text{ strongly } \text{ in } L^q(Q) \text{ for } \text{ all } q\in [1,+\infty ), \end{aligned}$$

it follows from (5.12) that

$$\begin{aligned} G_n(|\varvec{v}^n|^2)\mathop {\mathrm {div}}\nolimits (\varvec{v}^n\otimes \varvec{v}^n)\rightharpoonup \mathop {\mathrm {div}}\nolimits (\varvec{v}\otimes \varvec{v}) \text { weakly in } L^{\frac{5}{4}}(Q)^3. \end{aligned}$$

(5.14)

Finally, with the obtained convergences it is standard to prove that

$$\begin{aligned} \begin{aligned} \int \limits _0^T \langle \partial _t \varvec{v}, \varvec{w}\rangle + \int \limits _{Q}(\mathbb {Z}+\mathbb {V}): \mathbb {D}\varvec{w}+ \int \limits _{Q}\varvec{s}\cdot \varvec{w}_{\varvec{\tau }} -\int \limits _{Q}(\varvec{v}\otimes \varvec{v}): \mathbb {D}\varvec{w}= \int \limits _{Q}p_1 \mathop {\mathrm {div}}\nolimits \varvec{w}\\ -\int \limits _{Q}\nabla p_2\cdot \varvec{w}+\int \limits _{Q}\varvec{b}\cdot \varvec{w}\ \ \text{ for } \text{ all } \varvec{w} \in L^2(0,T;W^{1,2}_{\varvec{n}})\cap L^5(Q)^3 \end{aligned} \end{aligned}$$

(5.15)

and

$$\begin{aligned} \int \limits _0^T\langle \partial _t p_\mathrm {f}, z\rangle - \int \limits _{Q}p_\mathrm {f}\varvec{v}\cdot \nabla z + \int \limits _{Q}\nabla p_\mathrm {f}:\nabla z =\int \limits _{Q}(\varvec{b}-p_s\varvec{v})\cdot \nabla z \ \ \text{ for } \text{ all } z\in L^4(0,T; W^{1,2}(\Omega )). \end{aligned}$$

(5.16)

Step 3. Attainment of the constitutive equations on the boundary. Using that

$$\begin{aligned} \varvec{s}^n\rightharpoonup \varvec{s}\text { weakly in } L^{\frac{8}{3}} (\Sigma ) \end{aligned}$$

and (5.13), it easily follows

$$\begin{aligned} \limsup _{n\rightarrow + \infty }\int \limits _{\Sigma } \varvec{s}^n\cdot \varvec{v}^n_\tau = \int \limits _{\Sigma } \varvec{s}\cdot \varvec{v}_\tau . \end{aligned}$$

Thus, a suitable adjustment of Proposition 3.1 implies that (2.4) is fulfilled.

Step 4. Attainment of the constitutive equations in the bulk. In order to employ Proposition 3.1 we need to prove the limsup property (3.5), but as the solution itself can not be used as test function in (5.15), we follow the strategy as in [2] and perform the $L^\infty$-truncation method. To this aim, we introduce

$$\begin{aligned} \varvec{w}^n:=T_{\lambda _n} (\varvec{v}^n-\varvec{v}):= (\varvec{v}^n-\varvec{v}) \min {\left\{ 1, \frac{\lambda ^n}{|\varvec{v}^n-\varvec{v}|}\right\} } \end{aligned}$$

where $\lambda ^n \in [A,B]$ with $0<A< B<\infty$ will be suitably chosen numbers independent of n, but depending on parameter N tending to $+\infty$, see details below. For the reader’s convenience, we recall all the properties of $\varvec{w}^n$ below,

$$\begin{aligned}&\varvec{w}^n\rightarrow 0 \ \text{ strongly } \text{ in } \ L^s(Q)\ \text{ for } \text{ every } \ s \in [1, +\infty ), \end{aligned}$$

(5.17)

$$\begin{aligned}&\varvec{w}^n\rightarrow 0 \ \text{ strongly } \text{ in } \ L^2(\Sigma ), \end{aligned}$$

(5.18)

$$\begin{aligned}&\varvec{w}^n \rightharpoonup 0\ \text{ weakly } \text{ in } \ L^2 (0,T; W^{1,2}_{\varvec{n}}), \end{aligned}$$

(5.19)

$$\begin{aligned}&|\mathrm{div} \ \varvec{w}^n| \le {\left\{ \begin{array}{ll} 0 &{} \text{ if } \ |\varvec{v}^n - \varvec{v}|\le \lambda ^n\\ \frac{2 \lambda ^n \left( |\nabla \varvec{v}^n| + |\nabla \varvec{v}|\right) }{|\varvec{v}^n - \varvec{v}|} &{} \text{ if } \ |\varvec{v}^n-\varvec{v}|> \lambda ^n, \end{array}\right. } \end{aligned}$$

(5.20)

$$\begin{aligned}&\nabla \varvec{w}^n={\left\{ \begin{array}{ll}\nabla \varvec{v}^n-\nabla \varvec{v}\ \ \ \ \ \text{ if } \ |\varvec{v}^n - \varvec{v}|\le \lambda ^n\\ \begin{array}{l} \frac{\lambda ^n}{|\varvec{v}^n - \varvec{v}|}(\nabla \varvec{v}^n-\nabla \varvec{v})-\lambda ^n(\varvec{v}^n-\varvec{v})\otimes \frac{(\nabla \varvec{v}^n-\nabla \varvec{v}) (\varvec{v}^n-\varvec{v})}{|\varvec{v}^n-\varvec{v}|^3} \\ \text{ if } \ |\varvec{v}^n-\varvec{v}|> \lambda ^n.\end{array} \end{array}\right. } \end{aligned}$$

(5.21)

Inserting $\varvec{w}^n$ as test function in (5.1), using the properties of $\varvec{w}^n$ we get (cfr. [2])

$$\begin{aligned} \limsup _{n\rightarrow \infty }\int \limits _{Q} (\mathbb {Z}^n+\mathbb {V}^n): \mathbb {D}{\varvec{w}^n}\le \limsup _{n\rightarrow \infty } \int \limits _{Q} |p_1^n| |\mathop {\mathrm {div}}\nolimits \varvec{w}^n|. \end{aligned}$$

(5.22)

Now, let us define

$$\begin{aligned} {\overline{\mathbb {Z}}}= {\left\{ \begin{array}{ll}\displaystyle \mathbb {O}\ &{}\text{ if } \ \mathbb {D}\varvec{v}=\mathbb {O},\\ \displaystyle \tau (p_\mathrm {f})\frac{ \mathbb {D}{\varvec{v}}}{|\mathbb {D}{\varvec{v}}|} &{}\text{ if } \ \mathbb {D}\varvec{v}\ne \mathbb {O}. \end{array}\right. } \end{aligned}$$

(5.23)

and

$$\begin{aligned} {\overline{\mathbb {V}}}=\left( 1-\frac{\delta _*}{|\mathbb {D}{\varvec{v}}|}\right) ^+\mathbb {D}{\varvec{v}}. \end{aligned}$$

(5.24)

Employing (5.19) and (5.20), formula (5.22) can be rewritten as

$$\begin{aligned} \limsup _{n\rightarrow \infty } \int \limits _{Q} (\mathbb {Z}^n-{\overline{\mathbb {Z}}}):\mathbb {D}{\varvec{w}^n}+\int \limits _{Q} (\mathbb {V}^n-{\overline{\mathbb {V}}}): \mathbb {D}{\varvec{w}^n}\le \limsup _{n\rightarrow \infty } \int \limits _{\{|\varvec{v}^n-\varvec{v}|>\lambda ^n\}} |p_1^n|\left( |\nabla \varvec{v}^n|+ |\nabla \varvec{v}|\right) \frac{\lambda ^n}{|\varvec{v}^n-\varvec{v}|}.\nonumber \\ \end{aligned}$$

(5.25)

Moving the part of the integral on the left-hand side on the set $\{|\varvec{v}^n-\varvec{v}|>\lambda ^n\}$ to the right, it follows that

$$\begin{aligned} \begin{aligned} \limsup _{n\rightarrow \infty } \int \limits _{\{|\varvec{v}^n-\varvec{v}|\le \lambda ^n\}} \left( {(\mathbb {Z}^n-{\overline{\mathbb {Z}}}): (\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})} + {(\mathbb {V}^n-{\overline{\mathbb {V}}}): (\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})}\right) \\ \le C \limsup _{n\rightarrow \infty } \int \limits _{\{|\varvec{v}^n-\varvec{v}|>\lambda ^n\}} {I^n \frac{\lambda ^n}{|\varvec{v}^n-\varvec{v}|}}\end{aligned}\end{aligned}$$

(5.26)

where

$$\begin{aligned} I^n:= |p_1^n|^2+ |\mathbb {Z}^n|^2 + |{\overline{\mathbb {Z}}}|^2 + |\mathbb {V}^n|^2+|{\overline{\mathbb {V}}}|^2+ |\nabla \varvec{v}^n|^2 + |\nabla \varvec{v}|^2. \end{aligned}$$

Note that it holds (see formula (6.60) in [2])

$$\begin{aligned} \left( (\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}- \mathbb {D}{\varvec{v}})\right) ^- \rightarrow 0\ \text{ strongly } \text{ in } \ L^1(Q), \end{aligned}$$

(5.27)

and analogously the monotonicity implies

$$\begin{aligned} \left( (\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}- \mathbb {D}{\varvec{v}})\right) ^- \rightarrow 0\ \text{ strongly } \text{ in } \ L^1(Q).\end{aligned}$$

(5.28)

Let $N\in {\mathbb {N}}$ and fix $A=N$, $B=N^{N+1}$. For $i\in \{1,\ldots ,N\}$ let us define

$$\begin{aligned} Q_i^n:= \{N^i\le |\varvec{v}^n-\varvec{v}|\le N^{i+1}\}. \end{aligned}$$

Since

$$\begin{aligned} \sum _{i=1}^N \int \limits _{Q_i^n} I^n \le C^*, \end{aligned}$$

for every n there exists $i_n\in \{1,\ldots ,N\}$ such that

$$\begin{aligned} \int \limits _{Q_{i_n}^n} I^n \le \frac{C^*}{N}. \end{aligned}$$

Set $\lambda ^n:= N^{i_n}$, then it holds

$$\begin{aligned} \int \limits _{\{|\varvec{v}^n-\varvec{v}|> \lambda ^n\}} I^n \frac{\lambda ^n}{|\varvec{v}^n-\varvec{v}|}= \int \limits _{Q_{i_n}^n} I^n \frac{N^{i_n}}{|\varvec{v}^n-\varvec{v}|} + \int \limits _{\{|\varvec{v}^n-\varvec{v}|> N^{i_n+1}\}} I^n \frac{N^{i_n}}{|\varvec{v}^n-\varvec{v}|} \le \frac{C^*}{N} \end{aligned}$$

where we keep the symbol $C^*$ for a different constant. The latter relation, (5.26), (5.27) and (5.28) give

$$\begin{aligned} \begin{aligned} \limsup _{n\rightarrow +\infty } \left( \int \limits _{\{|\varvec{v}^n-\varvec{v}|\le \lambda ^n\}} |(\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})| + \int \limits _{\{|\varvec{v}^n-\varvec{v}|\le \lambda ^n\}} |(\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})| \right) \le \frac{C^*}{N}. \end{aligned}\end{aligned}$$

(5.29)

Using that

$$\begin{aligned} \begin{aligned} \int \limits _{Q} \sqrt{|(\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|} = \int \limits _{\{|\varvec{v}^n-\varvec{v}|\le N\}} \sqrt{|(\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|} \\ + \int \limits _{\{|\varvec{v}^n-\varvec{v}|> N\}} \sqrt{|(\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|} \end{aligned}\end{aligned}$$

(5.30)

and

$$\begin{aligned} \begin{aligned} \int \limits _{ Q} \sqrt{|(\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|} = \int \limits _{\{|\varvec{v}^n-\varvec{v}|\le N\}} \sqrt{|(\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|} \\ + \int \limits _{\{|\varvec{v}^n-\varvec{v}|> N\}} \sqrt{|(\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|}, \end{aligned}\end{aligned}$$

(5.31)

by Hölder’s and Chebyshev’s inequalities we obtain

$$\begin{aligned} \limsup _{n\rightarrow +\infty }\int \limits _{Q} \sqrt{|(\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|} \le \frac{2{\overline{C}}}{\sqrt{N}}, \end{aligned}$$

(5.32)

and

$$\begin{aligned} \limsup _{n\rightarrow +\infty }\int \limits _{ Q} \sqrt{|(\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})|}\le \frac{2{\overline{C}}}{\sqrt{N}}, \end{aligned}$$

(5.33)

which means, by letting $N\rightarrow \infty$, that

$$\begin{aligned}&(\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})\rightarrow 0 \ \text{ a.e. } \text{ in } \ Q, \\&(\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})\rightarrow 0 \ \text{ a.e. } \text{ in } \ Q. \end{aligned}$$

Egoroff’s theorem then gives that for all $\varepsilon >0$ there exists $U\subset Q$, $|Q\setminus U|\le \varepsilon$ such that

$$\begin{aligned}&\int \limits _U {(\mathbb {Z}^n-{\overline{\mathbb {Z}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})}\rightarrow 0 , \end{aligned}$$

(5.34)

$$\begin{aligned}&\int \limits _U {(\mathbb {V}^n-{\overline{\mathbb {V}}}):(\mathbb {D}{\varvec{v}^n}-\mathbb {D}{\varvec{v}})}\rightarrow 0. \end{aligned}$$

(5.35)

We conclude, thanks to the weak convergences of $\mathbb {Z}^n, \mathbb {V}^n, \mathbb {D}{\varvec{v}^n}$ respectively to $\mathbb {Z}, \mathbb {V}, \mathbb {D}{\varvec{v}}$ that

$$\begin{aligned} \lim _{n\rightarrow \infty }\int \limits _U \mathbb {Z}^n: \mathbb {D}{\varvec{v}^n} = \lim _{n\rightarrow \infty } \int \limits _U \mathbb {Z}^n: \mathbb {D}{\varvec{v}}= \int \limits _U \mathbb {Z}: \mathbb {D}{\varvec{v}} \end{aligned}$$

and

$$\begin{aligned} \lim _{n\rightarrow \infty }\int \limits _U \mathbb {V}^n: \mathbb {D}{\varvec{v}^n} = \lim _{n\rightarrow \infty } \int \limits _U \mathbb {V}^n: \mathbb {D}{\varvec{v}}= \int \limits _U \mathbb {V}: \mathbb {D}{\varvec{v}}. \end{aligned}$$

Finally, all assumptions of Convergence Lemma 4.1 are fulfilled, thus (1.2) holds a.e. in U. Since $|Q\setminus U|\le \varepsilon$ we can let $\varepsilon \rightarrow 0$ and obtain that (1.2) holds a.e. in Q. Theorem 2.1 is proved.

Appendix

Our goal is to prove Proposition 4.1. Let us recall that in this section we fix $n\in {\mathbb {N}}$, we consider $G_n$ smooth function with the properties stated at the beginning of Sect. 4 and the regularization of the material responses given in (4.2) and (4.3). In what follows, to simplify the notation we drop the indices n.

Proof of Proposition 4.1

The proof is split in the following steps.

Step 1. Approximations. For any $m\in {\mathbb {N}}$, we look for

$$\begin{aligned} \varvec{v}^{m}(t,x):=\sum _{r=1}^m c_r^{m}(t) \varvec{w}^r(x), \ \ \ p_\mathrm {f}^{m}(t, x):=\sum _{r=1}^m d_r^{m}(t) z^r(x) \end{aligned}$$

(6.1)

satisfying

$$\begin{aligned}&\left( \frac{\hbox {d} \varvec{v}^{m}}{\hbox {d}t}, \varvec{w}^r\right) + ({\mathbb {S}}^m, \mathbb {D}\varvec{w}^r ) +( \mathop {\mathrm {div}}\nolimits (\varvec{v}^{m}\otimes \varvec{v}^{m})G(|\varvec{v}^{m}|^2),\varvec{w}^r) \nonumber \\&\quad + ( \varvec{s}^{m}, \varvec{w}^r)_{\partial \Omega } = (\varvec{b},\varvec{w}^r) \ \ r=1,\ldots , m \end{aligned}$$

(6.2)

where

$$\begin{aligned}&\begin{array}{l} {\mathbb {S}}^m:={\mathcal {S}}(p_\mathrm {f}^m, \mathbb {D}\varvec{v}^m)=\tau (p_\mathrm {f}^m)\frac{\mathbb {D}{\varvec{v}^m}}{|\mathbb {D}{\varvec{v}^m}|+\frac{1}{n}}+\mathbb {D}{\varvec{v}^m}\left( 1-\frac{\delta _*}{|\mathbb {D}{\varvec{v}^m}|}\right) ^+ \\ \text{ with } \tau (p_\mathrm {f}^m)=(p_s-p_\mathrm {f}^m)^+, \end{array} \end{aligned}$$

(6.3)

$$\begin{aligned}&\varvec{s}^m:=s(\varvec{v}_\tau ^m)=s_* \frac{\varvec{v}_{\tau }^m}{|\varvec{v}_{\tau }^m|+\frac{1}{n}} +\varvec{v}_\tau ^m\left( 1-\frac{\beta _*}{|\varvec{v}_\tau ^m|}\right) ^+, \end{aligned}$$

(6.4)

and

$$\begin{aligned} \begin{array}{l} \left( \partial _t p_\mathrm {f}^{m}, z^r \right) - ( p_\mathrm {f}^{m}\varvec{v}^{m},\nabla z^r ) + ( \nabla p_\mathrm {f}^{m}, \nabla z^r) = (\varvec{b}, \nabla z^r) -(p_s\varvec{v}^{m},\nabla z^r ) \\ \,\,\, r=1,\ldots , m, \end{array}\end{aligned}$$

(6.5)

where $\{\varvec{w}^i\}_{i\in {\mathbb {N}}}$ is an orthogonal basis in $W_{{\mathbf {n}}, \mathop {\mathrm {div}}\nolimits }^{1,2}$ consisting of eigenfunctions of the Stokes operator with boundary conditions $\varvec{w}^i\cdot \varvec{n}=0$ and $[(\mathbb {D}\varvec{w}^i) \varvec{n}]_\tau =\varvec{0}$ on $\partial \Omega$, while $\{z_j\}_{j\in {\mathbb {N}}}$ is an orthogonal basis in $W^{1,2}(\Omega )$ consisting of eigenfunctions of the Laplace operator subject to the Neumann homogeneous boundary conditions. The system is supplemented with the corresponding initial conditions $\varvec{v}^m_0$ and $p_0^m$, obtained by projection $\varvec{v}_0\in L^{2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }$ onto the span of $[ \varvec{w}^1,\ldots ,\varvec{w}^m]$ and respectively $p_0\in L^2(\Omega )$ onto the span of $[z^1,\ldots , z^m]$. Then the local in time existence of $\varvec{v}^{m}$ and $p_\mathrm {f}^{m}$ follows from the Caratheodory theory for systems of ordinary differential equations, whereas the global in time existence is a consequence of the uniform estimates established below.

Step 2. Uniform estimates. Multiplying (6.2) by $c^m_r(t)$ and (6.5) by $d^m_r(t)$ and taking the sum over $r=1, \ldots , m$, we obtain

$$\begin{aligned}&\begin{array}{l} \frac{1}{2} \frac{\hbox {d}}{\hbox {d}t}\Vert \varvec{v}^{m}(t)\Vert _2^2 + \int \limits _{\{|\mathbb {D}\varvec{v}^m|>\delta _*\}}|\mathbb {D}\varvec{v}^m|^2 + \int \limits _{\Omega } \tau (p_\mathrm {f}^m)\frac{|\mathbb {D}{\varvec{v}^m}|^2}{|\mathbb {D}{\varvec{v}^m}|+\frac{1}{n}} \\ \,\,\,+\int \limits _{\partial \Omega } \varvec{s}^{m} \cdot \varvec{v}_{ \tau }^{m}=( \varvec{b}, \varvec{v}^{m})+ \int \limits _{\{|\mathbb {D}\varvec{v}^m|>\delta _*\}} \delta _*|\mathbb {D}\varvec{v}^m|, \end{array} \end{aligned}$$

(6.6)

$$\begin{aligned}&\frac{1}{2} \frac{\hbox {d}}{\hbox {d}t} \Vert p_\mathrm {f}^{m}(t)\Vert _2^2 + \Vert \nabla p_{f}^{m}(t)\Vert _2^2 = (\varvec{b}, \nabla p_\mathrm {f}^m)-(p_s\varvec{v}^m, \nabla p_\mathrm {f}^m), \end{aligned}$$

(6.7)

Adding $\int \limits _{\{|\mathbb {D}\varvec{v}^m|\le \delta _*\}}|\mathbb {D}\varvec{v}^m|^2$ to both sides of (6.6)

$$\begin{aligned} \frac{1}{2} \frac{\hbox {d}}{\hbox {d}t}\Vert \varvec{v}^{m}(t)\Vert _2^2 + \int \limits _{\Omega }|\mathbb {D}\varvec{v}^m|^2 + \int \limits _{\Omega } \tau (p_\mathrm {f}^m)\frac{|\mathbb {D}{\varvec{v}^m}|^2}{|\mathbb {D}{\varvec{v}^m}|+\frac{1}{n}} +\int \limits _{\partial \Omega } s_n(\varvec{v}_\tau ^{m}) \cdot \varvec{v}_{ \tau }^{m}\le (\varvec{b}, \varvec{v}^{m}) + \int \limits _{\Omega }\delta _*|\mathbb {D}\varvec{v}^m|+ \delta _*^2|\Omega | \end{aligned}$$

and then by Young’s inequality, we get

$$\begin{aligned} \frac{1}{2} \frac{\hbox {d}}{\hbox {d}t}\Vert \varvec{v}^{m}(t)\Vert _2^2 + \frac{1}{2} \int \limits _{\Omega }|\mathbb {D}\varvec{v}^m|^2 + \int \limits _{\Omega } \tau (p_\mathrm {f}^m)\frac{|\mathbb {D}{\varvec{v}^m}|^2}{|\mathbb {D}{\varvec{v}^m}|+\frac{1}{n}} +\int \limits _{\partial \Omega } s_n(\varvec{v}_\tau ^{m}) \cdot \varvec{v}_{ \tau }^{m}\le (\varvec{b}, \varvec{v}^{m}) +\frac{3}{2} \delta _*^2|\Omega |. \end{aligned}$$

(6.8)

Integrating in time, by Korn’s and Young’s inequalities, using also the fact that the last two terms on the left-hand side of (6.8) are non-negative, one concludes that

$$\begin{aligned} \sup _{t\in [0,T]} \Vert \varvec{v}^{m} (t)\Vert _2^2 + \int \limits _{Q} |\mathbb {D}{\varvec{v}^{m}}|^2 \le C(\varvec{b}, \varvec{v}_0, \delta _*, |Q|). \end{aligned}$$

(6.9)

By the interpolation inequality

$$\begin{aligned}&\Vert z\Vert _{\frac{10}{3}}\le \Vert z\Vert _2^{\frac{2}{5}}\Vert z\Vert _6^{\frac{3}{5}}\le C\Vert z\Vert _2^{\frac{2}{5}}\Vert z\Vert _{1,2}^{\frac{3}{5}} \end{aligned}$$

(6.10)

and the trace inequalities (see [6, Lemma 1.11]), we obtain

$$\begin{aligned} \sup _m \left( \Vert \varvec{v}^m\Vert _{\frac{10}{3}, Q} + \Vert \varvec{v}^m_\tau \Vert _{\frac{8}{3}, \Sigma }\right) <+\infty . \end{aligned}$$

(6.11)

As a consequence integrating in time (6.7) we deduce that

$$\begin{aligned} \sup _{t\in [0,T]}\Vert p_\mathrm {f}^{m}(t)\Vert _2^2+ \int \limits _{Q}|\nabla p_\mathrm {f}^{m}|^2 \le C\Vert \varvec{b}\Vert _{2,Q}^2+ C\Vert p_s\Vert _{5,Q}^2 \Vert \varvec{v}^m\Vert _{\frac{10}{3},Q}^2+ \Vert {p}_0 \Vert _2^2 \end{aligned}$$

(6.12)

and thus

$$\begin{aligned} \sup _{t\in (0,T)}\Vert p_\mathrm {f}^m (t)\Vert _2 + \Vert \nabla p_\mathrm {f}^m\Vert _{L^2(Q)}\le C(\varvec{b}, p_s, p_0). \end{aligned}$$

(6.13)

Again (6.10) gives

$$\begin{aligned} \sup _m \Vert p_\mathrm {f}\Vert _{\frac{10}{3}, Q} < +\infty . \end{aligned}$$

(6.14)

Recalling the explicit formulas for $\mathbb {S}^m$ and $\varvec{s}^m$ it then follows

$$\begin{aligned} \sup _m\left( \Vert \mathbb {S}^m\Vert _{2, Q} +\Vert \varvec{s}^m\Vert _{\frac{8}{3}, \Sigma }\right) <+\infty . \end{aligned}$$

(6.15)

Employing the inequality

$$\begin{aligned} \Vert z\Vert _{4}\le \Vert z\Vert _2^{\frac{1}{4}}\Vert z\Vert _6^{\frac{3}{4}}\le C\Vert z\Vert _2^{\frac{1}{4}}\Vert z\Vert _{1,2}^{\frac{3}{4}} \end{aligned}$$

we deduce corresponding uniform estimates for $\varvec{v}^m$ and $p_\mathrm {f}^m$ respectively in $L^4(Q)^3$ and $L^4(Q)$, then by virtue of them it results

$$\begin{aligned} \sup _m \Vert \partial _tp_\mathrm {f}^m\Vert _{(L^4(0,T; W^{1,2}(\Omega )))^*}<+\infty . \end{aligned}$$

(6.16)

Analogously and by virtue of the truncation in the convective term, we also get

$$\begin{aligned} \sup _m \Vert \partial _t\varvec{v}^m\Vert _{{(L^2(0,T; W^{1,2}_{\varvec{n}, \mathop {\mathrm {div}}\nolimits }))}^*}<+\infty . \end{aligned}$$

(6.17)

Step 3. Limit. By virtue of uniform estimates established above there exist subsequences of $\{\varvec{v}^m\}, \{p_\mathrm {f}^m\}$, $\{\mathbb {S}^m\}$ and $\{\varvec{s}^m\}$, converging respectively weakly (or *-weakly) to $\varvec{v}, p_\mathrm {f}, \mathbb {S}$ and $\varvec{s}$ in the corresponding function spaces. Furthermore, Aubin-Lions compactness Lemma and its variant including the trace theorem imply the following strong convergences:

$$\begin{aligned}&\varvec{v}^m\rightarrow \varvec{v} \text{ a.e. } \text{ in } Q \text{ and } \text{ strongly } \text{ in } L^q(Q)^3 \text{ for } \text{ any } q\in \left[ 1, \frac{10}{3}\right) , \end{aligned}$$

(6.18)

$$\begin{aligned}&p_\mathrm {f}^m\rightarrow p_\mathrm {f} \text{ a.e. } \text{ in } Q \text{ and } \text{ strongly } \text{ in } L^q(Q)^3 \text{ for } \text{ any } q\in \left[ 1, \frac{10}{3}\right) , \end{aligned}$$

(6.19)

$$\begin{aligned}&\varvec{v}^m_\tau \rightarrow \varvec{v}_\tau \text{ a.e. } \text{ in } \Sigma \text{ and } \text{ strongly } \text{ in } L^q(\Sigma )^3 \text{ for } \text{ any } q\in \left[ 1, \frac{8}{3}\right) . \end{aligned}$$

(6.20)

As a consequence $\varvec{v}, p_\mathrm {f}, \mathbb {S}$ and $\varvec{s}$ fulfill the weak formulations stated in Proposition 4.1.

Step 4. Attainment of the constitutive equations. The convergence

$$\begin{aligned} \varvec{s}^m \rightharpoonup \varvec{s} \text{ weakly } \text{ in } L^{\frac{8}{3}}(\Sigma ) \end{aligned}$$

together with (6.20) ensures that

$$\begin{aligned} \lim _{m\rightarrow +\infty } \int \limits _{\Sigma }\varvec{s}^m\cdot \varvec{v}_\tau ^m=\int \limits _{\Sigma }\varvec{s}\cdot \varvec{v}_\tau . \end{aligned}$$

(6.21)

Then, thanks to the monotonicity it is standard to prove that

$$\begin{aligned} \varvec{s}=s(\varvec{v}_\tau ). \end{aligned}$$

Next, it follows from the monotonicity that

$$\begin{aligned} 0\le \int \limits _{Q}\left( \mathbb {S}^m- {\mathcal {S}}(p_\mathrm {f}^m, \mathbb {A})\right) :\left( \mathbb {D}\varvec{v}^m-\mathbb {A}\right) \ \ \text{ for } \text{ all } \mathbb {A}\in L^2(Q). \end{aligned}$$

(6.22)

Now, note that by (6.19),

$$\begin{aligned}&{\mathcal {S}}(p_\mathrm {f}^m, \mathbb {A}):=(p_s-p_\mathrm {f}^m)^+\frac{\mathbb {A}}{|\mathbb {A}|+\frac{1}{n}}+\mathbb {A}\left( 1-\frac{\delta _*}{|\mathbb {A}|}\right) ^+\nonumber \\&\quad \rightarrow (p_s-p_\mathrm {f})^+\frac{\mathbb {A}}{|\mathbb {A}|+\frac{1}{n}}+\mathbb {A}\left( 1-\frac{\delta _*}{|\mathbb {A}|}\right) ^+=:{\mathcal {S}}(p_\mathrm {f}, \mathbb {A}) \text{ strongly } \text{ in } L^2(Q) \end{aligned}$$

(6.23)

while, as $\varvec{v}$ can play the role of a test function in the established weak formulation, it is standard to obtain

$$\begin{aligned} \limsup _{m\rightarrow +\infty } \int \limits _{Q}\mathbb {S}^m:\mathbb {D}\varvec{v}^m\le \int \limits _{Q}\mathbb {S}:\mathbb {D}\varvec{v}. \end{aligned}$$

(6.24)

Finally, thanks to the convergences

$$\begin{aligned}&\mathbb {D}\varvec{v}^m\rightharpoonup \mathbb {D}\varvec{v} \text{ weakly } \text{ in } L^2(Q), \\&\mathbb {S}^m\rightharpoonup \mathbb {S} \text{ weakly } \text{ in } L^2(Q), \end{aligned}$$

(6.23) and (6.24), the limit as $m\rightarrow + \infty$ in (6.22) gives

$$\begin{aligned} 0\le \int \limits _{Q}\left( \mathbb {S}- {\mathcal {S}}(p_\mathrm {f}, \mathbb {A})\right) :\left( \mathbb {D}\varvec{v}-\mathbb {A}\right) \, \hbox { for all } \mathbb {A}\in L^2(Q). \end{aligned}$$

(6.25)

At this point, it is standard to choose $\mathbb {A}=\mathbb {D}\varvec{v}\pm \varepsilon \mathbb {B}$ for arbitrary $\mathbb {B}\in L^2(Q)$ and $\varepsilon >0$ and arrive at

$$\begin{aligned} 0= \int \limits _{Q}\mathbb {B}:\left( \mathbb {S}- {\mathcal {S}}(p_\mathrm {f},\mathbb {D}\varvec{v})\right) \ \ \text{ for } \text{ all } \mathbb {B}\in L^2(Q), \end{aligned}$$

which implies $\mathbb {S}={\mathcal {S}}(p_\mathrm {f},\mathbb {D}\varvec{v})$ a.e. in Q. The proof of Proposition 4.1 is complete. $\square$