1 Introduction

In this paper we analyse a model describing the stationary flow of an incompressible non-Newtonian fluid of Bingham type in a bounded domain with the subdifferential boundary conditions. The constitutive law of the Bingham type is one of the simplest models describing the flow behaviour of a fluid with a yield stress. It is characterized by the fact that at low stresses, the medium behavies like a rigid body until a stress function reaches its yield point, while above this limit, the medium already flows like an incompressible viscous fluid. The former and current interest in Bingham type fluids comes and is stimulated by various applications in mathematical hydrodynamics and in a number of technological and industrial processes. The mathematical model was suggested by Bingham [5] in 1916. There is a large number of real substances such as mineral suspensions, adhesives like wall paper paste, wet beach sand, polymer melts and solutions, biological fluids like blood in the capillaries, coal slurries, fire fighting foams, etc. which can be described by the Bingham model.

The system of partial differential equations involving the Bingham constitutive law has been studied in many papers, for instance, in Duvaut and Lions in [14, Chap. VI], where the existence of weak solutions is demonstrated. More recently, the stationary Bingham model has been studied in [8] by using a mixed variational formulation with Lagrange multipliers and in [2] by the variational inequality formulation that uses the abstract result of Lions [23, Theorem 8.5, p. 251].

This paper is inspired by a series of recent papers dealing with mathematical models describing steady-state non-Newtonian fluids with multivalued nonmonotone subdifferential frictional boundary conditions, see [2, 11, 13]. The novelties of the paper are following. The first novelty of the paper is a result, see Theorem 4, on existence of solution to a variational–hemivariational inequality proved by using a surjectivity theorem for pseudomonotone operators. In contrast to [13, 26, 27], in our framework, we remove the assumption on the relaxed monotonicity of the generalized subgradient of a locally Lipschitz potential. In this way, the problem can be described by a wider class of boundary conditions involving nonmonotone law with a nonconvex and nondifferentiable superpotential. As far as we know there does not exist any uniqueness result for this problem.

The second novelty is to deal with the Navier–Stokes problem with the impermeability boundary condition and a nonlinear boundary condition with the Clarke subgradient. The latter represents a generalization of the Navier–Fujita slip condition from [21, 22] in several directions, see Sect. 3. We also generalize some existence results obtained earlier for the Oseen model in [13, 24]. In particular, our existence result in Theorem 4 extends [24, Theorem 4.1] proved for the Oseen model under an additional Clarke regularity of the nonconvex potential and a more restrictive hypothesis on the viscosity function. We observe that the class of monotone operators is too narrow to deal with the sum of the governing term and the convective term in the conservation law. However, the operator equation enters the theory of pseudomontone operators as a sum of a pseudomonotone operator and a compact one. On the other hand, due to nonconvexity of the potential, the problem can not be treated within the theory of variational inequalities, as in [2]. To the best of our knowledge, the Bingham model has not been studied so far by the theory of hemivariational inequaties. The weak formulation of the problem studied in the paper naturally leads to a variational–hemivariational inequality.

The third novelty is to establish continuous dependence of solutions on the second member and demonstrate the existence of solutions to an optimal control problem for the Bingham model. Further, prove a new convergence result which, in particular, shows that if the yield stress tends to zero, then the Bingham fluid tends to behave as a Newtonian one.

Note that the frictional boundary condition is described by a nonlinearity of the form \(k \partial j\). In general, this kind of nonlinearity can not be represented by a purely hemivariational inequality since there is not a potential G with \(\partial G = k \partial j\). This type of nonlinearity leads to a “hemivariational term” on the boundary in the inequality while the Bingham law leads to a “variational part” in the domain. Various models in fluid mechanics have studied within the theory of variational–hemivariational inequalities in [25] for Navier–Stokes equations, in [11, 12] for non-Newtonian and generalized Newtonian fluids, in [13] for stationary and evolutionary Oseen models.

The rest of the paper is structured as follows. After a preliminaries in Sect. 2, in Sect. 3 we deliver both classical and variational formulation of the Bingham model, and state the main existence result of the paper, Theorem 4. Section 4 is devoted to the proof of Theorem 4. In Sect. 5 we provide a result on the continuous dependence of solutions on the second member, analyse optimal control problem for the variational–hemivariational inequality, and study the dependence of the solution on the yield limit. Finally, Sect. 6 provides a short indication to open problems for the future work.

2 Mathematical Background

The study on variational inequalities uses results from nonsmooth analysis and the theory of monotone operators which comprehensive presentation can be found in [7, 9, 10, 26, 28].

Let X be a normed space with norm denoted by \(\Vert \cdot \Vert _X\). Let \(X^*\) stand for its topological dual, and the duality brackets for the pair \((X^*, X)\) is denoted by \(\langle \cdot ,\cdot \rangle _{X^*\times X}\). If no confusion arises, we often skip the subscripts. The weak and norm convergences in X are denoted by “\(\rightharpoonup \)" and “\(\rightarrow \)", respectively. By \({\mathcal {L}}(E, F)\) we denote the space of linear and bounded operators from a normed space E to a normed space F endowed with the operator norm \(\Vert \cdot \Vert _{{\mathcal {L}}(E,F)}\). Given a set \(G \subset X\), we set \(\Vert G \Vert _X = \sup \{ \Vert x \Vert _X \mid x \in G \}\).

Single-valued monotone operators Let X be a Banach space. An operator \(A :X \rightarrow X^*\) is said to be monotone, if for all u, \(v \in X\), it holds \(\langle Au - A v, u-v \rangle \ge 0\). It is called maximal monotone, if it is monotone, and \(\langle Au - w, u - v \rangle \ge 0\) for any \(u \in X\) entails \(w = Av\). Operator A is bounded, if A maps bounded sets of X into bounded sets of \(X^*\). Operator A is called pseudomonotone, if it is bounded and \(u_n \rightharpoonup u\) in X with \(\displaystyle \limsup \langle A u_n, u_n -u \rangle \le 0\) imply \(\displaystyle \langle A u, u - v \rangle \le \liminf \langle A u_n, u_n - v \rangle \) for all \(v \in X\). It is known that if X is a reflexive Banach space, then \(A :X \rightarrow X^*\) is pseudomonotone, if and only if it is bounded and \(u_n \rightharpoonup u\) in X with \(\displaystyle \limsup \,\langle A u_n, u_n - u \rangle \le 0\) imply \(\displaystyle \lim \,\langle A u_n, u_n - u \rangle = 0\) and \(A u_n \rightharpoonup A u\) in \(X^*\).

Multivalued monotone operators For an operator \(T :X \rightarrow 2^{X^*}\), its domain D(T) and graph Gr(T) are defined, respectively, by

$$\begin{aligned} D(T) = \{ x \in X \mid Tx \not = \emptyset \},\quad Gr(T)= \{ (x, x^*) \in X \times X^* \mid x^* \in T x \}. \end{aligned}$$

An operator \(T :X \rightarrow 2^{X^*}\) is called monotone when \(\langle u^* - v^*, u - v \rangle \ge 0\) for all \((u, u^*)\), \((v, v^*) \in Gr(T)\). It is said to be maximal monotone, if it is monotone and maximal in the sense of inclusion of graphs in the family of monotone operators from X to \(2^{X^*}\). Operator T is called coercive, if \(\inf \{\langle u^*, u \rangle \mid u^* \in Tu\} \ge \alpha (\Vert u \Vert _X) \, \Vert u \Vert _X\) for all \(u \in X\), where a function \(\alpha :[0,\infty ) \rightarrow {\mathbb {R}}\) is such that \(\alpha (r) \rightarrow +\infty \) as \(r \rightarrow +\infty \).

Let X be a reflexive Banach space. An operator \(T :X \rightarrow 2^{X^*}\) is pseudomonotone if:

  1. (a)

    for every \(u\in X\), the set \(T u \subset X^*\) is nonempty, closed and convex,

  2. (b)

    T is upper semicontinuous from each finite dimensional subspace of X into \(X^*\) endowed with weak topology,

  3. (c)

    for all sequences \(\{u_n\} \subset X\) and \(\{u_n^*\}\subset X^*\) such that \(u_n \rightharpoonup u\) in X, \(u_n^* \in T u_n\) for all \(n \ge 1\) and \(\limsup \langle u^*_n, u_n - u \rangle \le 0\), we have that for every \(v \in X\), there exists \(u^*(v) \in T u\) such that

    $$\begin{aligned} \langle u^*(v),u-v\rangle \le \liminf \, \langle u_n^*, u_n - v\rangle . \end{aligned}$$

Let X be a reflexive Banach space. An operator \(T :X \rightarrow 2^{X^*}\) is generalized pseudomonotone if for any sequences \(\{ u_n \}\subset X\) and \(\{ u_n^*\} \subset X^*\) such that \(u_n\rightharpoonup u\) in X, \(u_n^*\in T u_n\) for \(n\ge 1\), \(u_n^* \rightharpoonup u^*\) in \(X^*\) and \(\limsup \langle u_n^*,u_n-u \rangle \le 0\), we have \(u^*\in T u\) and \(\lim \, \langle u_n^*, u_n \rangle = \langle u^*, u \rangle \). It is known that if T is pseudomonotone, then it is generalized pseudomonotone. Conversely, if T is a bounded generalized pseudomonotone operator such that for all \(u \in X\), Tu is a nonempty, closed and convex subset of \(X^*\), then T is pseudomonotone.

Generalized subgradient Let X be a Banach space. Let \(\varphi :X \rightarrow {\mathbb {R}}\cup \{ +\infty \}\) be a proper, convex and lower semicontinuous function. The mapping \(\partial \varphi :X \rightarrow 2^{X^*}\) defined by

$$\begin{aligned} \partial \varphi (x) = \left\{ x^*\in X^* \mid \langle x^*, v -x \rangle \le \varphi (v)-\varphi (x) \ \text{ for } \text{ all } \ v \in X\right\} \end{aligned}$$

is called the subdifferential of \(\varphi \). An element \(x^* \in \partial \varphi (x)\) is called a subgradient of \(\varphi \) at x. Let \(h :X \rightarrow {\mathbb {R}}\) be a locally Lipschitz function. The generalized (Clarke) directional derivative of h at the point \(x \in X\) in the direction \(v \in X\) is defined by

$$\begin{aligned} h^{0}(x; v) = \limsup _{y \rightarrow x, \ \lambda \downarrow 0} \frac{h(y + \lambda v) - h(y)}{\lambda }. \end{aligned}$$

The generalized subgradient of h at x is a subset of the dual space \(X^*\) given by

$$\begin{aligned} \partial h (x) = \{\, \zeta \in X^* \mid h^{0}(x; v) \ge {\langle \zeta , v \rangle } \ \text{ for } \text{ all } \ v \in X \, \}. \end{aligned}$$

3 Statement of the Flow Problem

In this section we provide the classical and variational formulations of the flow model of a Bingham fluid, and give a statement of the main existence result.

Let \(\Omega \subset {\mathbb {R}}^d\) be a bounded domain (open and connected set) with Lipschitz boundary \(\Gamma =\partial \Omega \) occupied by the incompressible Bingham fluid for \(d = 2\) and \(d = 3\) (the cases of interest). The boundary \(\Gamma \) consists of smooth, disjoint, and measurable parts \(\Gamma _0\) and \(\Gamma _1\) such that \(\mathrm{meas}(\Gamma _0) > 0\) while the part \(\Gamma _1\) can be empty. The unit outward normal vector on boundary is denoted by \({\varvec{\nu }}\). Moreover, \({\mathbb {M}}^{d}\) denotes the class of symmetric \(d\times d\) matrices.

The classical formulation of the steady-state flow problem reads as follows.

Problem 1

Find a flow velocity \({\varvec{u}}:\Omega \rightarrow {\mathbb {R}}^d\), an extra stress tensor \({\mathbb {S}} :{\mathbb {M}}^d \rightarrow {\mathbb {M}}^d\), and a pressure \(p :\Omega \rightarrow {\mathbb {R}}\) such that

$$\begin{aligned}&- \mathop {\mathrm{Div}}\nolimits {{\mathbb {S}}} + \rho \mathop {\mathrm{Div}}\nolimits ({\varvec{u}}\otimes {\varvec{u}}) + \nabla p = \rho {\varvec{f}}&\text{ in } \&\ \ \Omega , \end{aligned}$$
(1)
$$\begin{aligned}&\left\{ \begin{array}{lll} \displaystyle {\mathbb {S}} = \mu (\Vert {\mathbb {D}}{\varvec{u}}\Vert )\, {\mathbb {D}}{\varvec{u}}+ g \, \frac{{\mathbb {D}}{\varvec{u}}}{\Vert {\mathbb {D}}{\varvec{u}}\Vert } &{}\quad \text{ if } \ \ \ {\mathbb {D}}{\varvec{u}}\not = {\varvec{0}}\\ \Vert {\mathbb {S}} \Vert \le g \qquad \qquad \qquad \qquad &{}\quad \text{ if } \ \ \ {\mathbb {D}} {\varvec{u}}= {\varvec{0}}\end{array}\right.&\mathrm{in} \&\ \ \Omega , \end{aligned}$$
(2)
$$\begin{aligned}&\mathop {\mathrm{div}}\nolimits {\varvec{u}}= 0&\text{ in } \&\ \ \Omega , \end{aligned}$$
(3)
$$\begin{aligned}&{\varvec{u}}=0&\text{ on } \&\ \ \Gamma _0, \end{aligned}$$
(4)
$$\begin{aligned}&\left\{ \begin{array}{lll} u_\nu = 0 \\ -{\varvec{\tau }}_\tau ({\varvec{u}}) \in k({\varvec{u}}_\tau ) \partial j ({\varvec{u}}_\tau ) \end{array}\right.&\mathrm{on} \&\ \ \Gamma _{1}. \end{aligned}$$
(5)

In Problem 1, the flow of an incompressible fluid is governed by the conservation law (1). The tensorial constitutive law (2) is a relation between the extra stress tensor \({\mathbb {S}}\) and the strain velocity tensor \({\mathbb {D}}{\varvec{u}}\) whose components are given by

$$\begin{aligned} {\mathbb {D}}_{ij}({\varvec{u}}) = \frac{1}{2} \Big ( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \Big ) \ \ \text{ for } \ \ {\varvec{u}}= (u_1,\ldots ,u_d). \end{aligned}$$

The law (2) represents the Bingham model which is used to analyze flows of materials for which the imposed stress must exceed a critical yield stress (called plasticity threshold) g to initiate motion, i.e., they behave as rigid bodies when the stress is low but flow as viscous fluids at high stress. Examples of Bingham fluids include cosmetics and personal care products, water suspensions of clay, drilling muds, volcanic lava and magmas, or sewage sludge. For \(g = 0\) and \(\mu (r) = \mu _0\), the constitutive law reduces to the model of Newtonian fluid for the Navier-Stokes equation. Here, \(\mu :[0, \infty ) \rightarrow {\mathbb {R}}\) is a given positive viscosity function, \(g :\Omega \rightarrow [0, +\infty )\) stands for the yield limit, and \(\rho \) is the constant density of the fluid, conveniently put equal to one. The solenoidal (divergence free) condition (3) of the form \(\mathop {\mathrm{div}}\nolimits {\varvec{u}}= \nabla \cdot {\varvec{u}}= \sum _{i=1}^d {{\mathbb {D}}}_{ii}({\varvec{u}}) = 0\) in \(\Omega \) states that the fluid is incompressible. The homogeneous Dirichlet boundary condition (4) means that the fluid adheres to the wall, see [26].

The symbols \(\mathop {\mathrm{Div}}\nolimits \) and \(\mathop {\mathrm{div}}\nolimits \) denote the divergence operators for tensor and vector valued functions \({\mathbb {S}} :\Omega \times {\mathbb {M}}^d \rightarrow {\mathbb {M}}^d\) and \({\varvec{u}}:\Omega \rightarrow {\mathbb {R}}^d\) defined by \(\mathop {\mathrm{Div}}\nolimits {\mathbb {S}} = (S_{ij,j})\) and \(\mathop {\mathrm{div}}\nolimits {\varvec{u}}= (u_{i,i})\), and the index that follows a comma represents the partial derivative with respect to the corresponding component of \({\varvec{x}}\). The tensor product of two vector fields \(\varvec{{v}}\), \({\varvec{w}}\in {\mathbb {R}}^d\) is the second-order tensor \({\varvec{\sigma }}\) defined by

$$\begin{aligned} {\varvec{\sigma }}= \varvec{{v}}\otimes {\varvec{w}}= (v_i w_j)_{1\le i,j\le d}, \ \ \text{ i.e. },\ \ \sigma _{ij} = v_i w_j \ \ \text{ for } \ \ i,j = 1, \ldots , d. \end{aligned}$$

The contracted product of tensors \({\varvec{\sigma }}\) and \({\varvec{\tau }}\) of order 2 is the real number defined by

$$\begin{aligned} {\varvec{\sigma }}: {\varvec{\tau }}= \sum _{1\le i, j \le d} \sigma _{ij} \tau _{ij}. \end{aligned}$$

Alternatively, the contracted tensor product is simply the Euclidean scalar product in the space of the matrices identified with \({\mathbb {R}}^{d\times d}\). For a vector-valued function \(\varvec{{v}}\) on \(\Gamma \), we denote by \(v_\nu \) and \(\varvec{{v}}_\tau \) its normal and tangential components defined by \(v_\nu = \varvec{{v}}\cdot {\varvec{\nu }}\) and \(\varvec{{v}}_\tau = \varvec{{v}}- v_\nu {\varvec{\nu }}\), respectively. Further, the total stress tensor is expressed by

$$\begin{aligned} {\varvec{\sigma }}({\varvec{u}}, p) = -p \, {\mathbb {I}} + {\mathbb {S}} ({\mathbb {D}}({\varvec{u}})) \ \ \ \text{ in } \ \ \Omega , \end{aligned}$$

where \({\mathbb {I}}\) is the identity \(d\times d\) matrix. The traction vector is denoted by \({\varvec{\tau }}({\varvec{u}},p) = {\varvec{\sigma }}({\varvec{u}},p) {\varvec{\nu }}\) on \(\Gamma \), so \(\tau _\nu ({\varvec{u}},p) = {\varvec{\tau }}({\varvec{u}},p) \cdot {\varvec{\nu }}\) and \({\varvec{\tau }}_\tau ({\varvec{u}}) = {\varvec{\tau }}({\varvec{u}},p) - \tau _\nu ({\varvec{u}}, p) {\varvec{\nu }}\) represent normal and tangential components of the traction vector, respectively. Note that \({\varvec{\tau }}\) and \(\tau _\nu \) do depend on the pressure p, and \({\varvec{\tau }}_\tau \) is independent of p. To simplify the notation, we often do not indicate explicitly the dependence of various functions on the spatial variable \({\varvec{x}}\in \Omega \cup \Gamma \). The inner products and norms on \({\mathbb {R}}^d\) and \({\mathbb {M}}^{d}\) are denoted by the intuitive notation.

For the analysis of the problem, we need the following spaces

$$\begin{aligned} {\widetilde{V}}= & {} \{\varvec{{v}}\in {\mathcal {C}}^\infty ({{\overline{\Omega }}};{\mathbb {R}}^d) \mid \mathop {\mathrm{div}}\nolimits \, \varvec{{v}}= 0 \ \text {in}\ \Omega ,\ \varvec{{v}}= 0\ \text {on}\ \Gamma _0, \ v_\nu =0 \ \text {on} \ \Gamma _1\}, \nonumber \\ V= & {} \ \text {closure of}\ {\widetilde{V}}\ \text {in}\ H^{1}(\Omega ;{\mathbb {R}}^d). \end{aligned}$$
(6)

The space V is endowed with the standard norm \(\Vert \varvec{{v}}\Vert = \Vert \varvec{{v}}\Vert _{H^{1}(\Omega ;{\mathbb {R}}^d)}\). We also consider the norm given by \(\Vert \varvec{{v}}\Vert _V = \Vert {{\mathbb {D}}} (\varvec{{v}}) \Vert _{L^2(\Omega ;{\mathbb {M}}^d)}\) for \(\varvec{{v}}\in V\). Using the Korn inequality, see, e.g., [11, Theorem 8],

$$\begin{aligned} c_K \Vert \varvec{{v}}\Vert _{H^{1}(\Omega ;{\mathbb {R}}^d)} \le \Vert {{\mathbb {D}}} (\varvec{{v}}) \Vert _{L^2(\Omega ;{\mathbb {M}}^d)} \ \ \text{ for } \ \ \varvec{{v}}\in V \ \ \text{ with } \ \ c_K > 0, \end{aligned}$$
(7)

it is known that \(\Vert \cdot \Vert _{H^{1}(\Omega ;{\mathbb {R}}^d)}\) and \(\Vert \cdot \Vert _V\) are the equivalent norms on V. By \(\langle \cdot , \cdot \rangle \) we denote the duality brackets for the pair \((V^*, V)\). Let \(H = L^2(\Omega ; {\mathbb {R}}^d)\). From the Rellich-Kondrachov theorem, see [30, Sect. 2.6.1, Theorem 6.1 and Corollary 6.2], we know that the embeddings

$$\begin{aligned} H^1(\Omega ;{\mathbb {R}}^3) \subset L^4(\Omega ;{\mathbb {R}}^3), \ \ \text{ and } \ \ H^1(\Omega ;{\mathbb {R}}^2) \subset L^q(\Omega ;{\mathbb {R}}^2) \ \ \text{ for } \text{ any } \ \ q \ge 1 \end{aligned}$$
(8)

are continuous and compact. By [30, Sect. 2.5.4, Theorems 5.5 and 5.7], the trace operator denoted by

$$\begin{aligned} \gamma :V \subset H^1(\Omega ;{\mathbb {R}}^2) \rightarrow L^2(\Gamma ;{\mathbb {R}}^d) \end{aligned}$$
(9)

is continuous and compact. We denote its norm in the space \({\mathcal {L}}(V,L^2(\Gamma ; {\mathbb {R}}^d))\) by \(\Vert \gamma \Vert \). Moreover, instead of \(\gamma \varvec{{v}}\), for simplicity, we often write \(\varvec{{v}}\).

We need the following hypotheses on the data.

  • \(\underline{H(\mu )}:\)    \(\mu :[0, \infty ) \rightarrow {\mathbb {R}}\) is a function such that

    • (i) \(\mu \) is continuous and \(0 < \mu _0 \le \mu (r) \le \mu _1\) for all \(r \in [0,\infty )\),

    • (ii) \((\mu (\Vert {\varvec{A}}\Vert ) {\varvec{A}}- \mu (\Vert {\varvec{B}}\Vert ) {\varvec{B}}) : ({\varvec{A}}- {\varvec{B}}) \ge 0\) for all \({\varvec{A}}\), \({\varvec{B}}\in {\mathbb {M}}^d\).

  • \(\underline{H(f, g)}:\)    \({\varvec{f}}\in H\), \(g \in L^2(\Omega )\), \(g \ge 0\).

  • \(\underline{H(j)}:\)    \(j :\Gamma _1\times {\mathbb {R}}^d \rightarrow {\mathbb {R}}\) is a function such that

    • (i) \(j(\cdot ,{\varvec{\xi }})\) is measurable on \(\Gamma _1\) for all \({\varvec{\xi }}\in {\mathbb {R}}^d\),

    • (ii) \(j({\varvec{x}},\cdot )\) is locally Lipschitz for a.e. \({\varvec{x}}\in \Gamma _1\),

    • (iii) \(\displaystyle \Vert \partial j({\varvec{x}},{\varvec{\xi }}) \Vert _{{\mathbb {R}}^d} \le c_0 ({\varvec{x}}) + c_1 \, \Vert {\varvec{\xi }}\Vert _{{\mathbb {R}}^d}\)   for all   \({\varvec{\xi }}\in {\mathbb {R}}^d\), a.e. \({\varvec{x}}\in \Gamma _1\) with \(c_0 \in L^2(\Gamma _1)\), \(c_0\), \(c_1 \ge 0\).

  • \(\underline{H(k)}:\)    \(k:\Gamma _1 \times {\mathbb {R}}^d \rightarrow {\mathbb {R}}\) is a function such that

    • (i) \(k(\cdot , {\varvec{\xi }})\) is measurable on \(\Gamma _1\) for all \({\varvec{\xi }}\in {\mathbb {R}}^d\),

    • (ii) \(k({\varvec{x}}, \cdot )\) is continuous on \({\mathbb {R}}^d\) for a.e. \({\varvec{x}}\in \Gamma _1\),

    • (iii) \(0 < k_0 \le k({\varvec{x}}, {\varvec{\xi }}) \le k_1\) for all \({\varvec{\xi }}\in {\mathbb {R}}^d\), a.e. \({\varvec{x}}\in \Gamma _1\).

In hypothesis H(j)(iii), and in what follows, \(\partial j\) denotes the generalized gradient of the function j with respect to its last variable.

Hypothesis \(H(\mu )\) holds for typical models like the Carreau-type and power-law models, see [18]. For instance, the Carreau law is described by

$$\begin{aligned} \mu (t) = \mu _\infty + (\mu _0-\mu _\infty ) (1 + \lambda \, t^2)^{\frac{r-2}{2}} \ \ \text{ for } \ \ t \ge 0, \end{aligned}$$

with \(\lambda > 0\), \(1< r \le 2\) and \(0< \mu _\infty < \mu _0\). This function satisfies \(\mu \in C^1([0,+\infty ))\) and

$$\begin{aligned} \mu _\infty \, (t-s) \le \mu (t) t - \mu (s) s \le \mu _0 \, (t-s) \ \ \text{ for } \text{ all } \ \ t \ge s \ge 0. \end{aligned}$$
(10)

It is also known that if the viscosity \(\mu \) satisfies (10), than \(H(\mu )\) holds with suitable constants \(\mu _0\), \(\mu _1\), \(\mu _2 >0\), see [1, condition (2.3)] and [4, Lemma 2.1]. Clearly the choice \(r=2\) leads to \(\mu (t) = \mu _0\) which is the linear Newtonian constitutive relation. Moreover, hypothesis \(H(\mu )\)(ii) is satisfied when \(\mu \) is a nondecreasing function, for example, \(\mu (r) = \sqrt{r} + 1/2\) for \(r \in [0,4]\), and 5/2 for \(r > 4\), or \(\mu (r) = (\arctan r)^{1/2} + \mu _0\) for \(r \ge 0\), see [2, Remark 3].

We shall explain the boundary condition (5). The first condition in (5) is called the impermeability (no leak) boundary condition, and the second one is called the nonmonotone slip boundary condition.

(i) The slip condition of Navier

$$\begin{aligned} -{\varvec{\tau }}_\tau ({\varvec{u}}) = \kappa \, {\varvec{u}}_{\tau }\ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$
(11)

where \(\kappa :\Gamma _1 \rightarrow (0,\infty )\) is a given function, is a prototype of (5) in which \(k({\varvec{x}}, {\varvec{\xi }}) = \kappa ({\varvec{x}})\) and \(j({\varvec{x}}, \xi ) = \frac{1}{2} \Vert {\varvec{\xi }}\Vert ^2\), see [29]. Condition (11) states that the tangential velocity is proportional to the shear stress. It has been discussed in [22] and the references therein.

(ii) The nonlinear Navier-type slip condition

$$\begin{aligned} -{\varvec{\tau }}_\tau ({\varvec{u}}) = h(\Vert {\varvec{u}}_\tau \Vert ) \, \frac{{\varvec{u}}_\tau }{\Vert {\varvec{u}}_\tau \Vert } \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$
(12)

where \(h:\Gamma _1 \times (0,\infty ) \rightarrow [0,\infty )\) is prescribed, was used in [21] to model the wall slip of non-Newtonian fluids.

(iii) The following threshold slip condition or the Navier–Fujita condition of frictional type

$$\begin{aligned} \left\{ \begin{array}{ll} \Vert {\varvec{\tau }}_\tau ({\varvec{u}})\Vert \le \alpha \ \text{ with } \\ \displaystyle \quad \text{ if } \ \Vert {\varvec{\tau }}_\tau ({\varvec{u}}) \Vert < \alpha , \ \text{ then } \ {\varvec{u}}_\tau = 0, \\[2mm] \quad \text{ if } \ \Vert {\varvec{\tau }}_\tau ({\varvec{u}}) \Vert = \alpha , \ \text{ then } \ {\varvec{u}}_{\tau } \not = 0 \ \text{ and } \ \displaystyle -{\varvec{\tau }}_\tau ({\varvec{u}}) = \alpha \, \frac{{\varvec{u}}_\tau }{\Vert {\varvec{u}}_\tau \Vert } \end{array} \right. \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$
(13)

where \(\alpha :\Gamma _1 \rightarrow (0,\infty )\) is a given prescribed slip threshold, has been treated in [15,16,17, 32, 33]. Condition (13) can be formulated as follows

$$\begin{aligned} \Vert {\varvec{\tau }}_\tau \Vert \le \alpha , \quad {\varvec{\tau }}_\tau ({\varvec{u}}) \cdot {\varvec{u}}_\tau + \alpha \, \Vert {\varvec{u}}_\tau \Vert = 0 \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$

The version of (13) of the form

$$\begin{aligned} \left\{ \begin{array}{ll} \Vert {\varvec{\tau }}_\tau ({\varvec{u}})\Vert \le \alpha \ \text{ with }\\ \displaystyle \quad \text{ if } \ \Vert {\varvec{\tau }}_\tau ({\varvec{u}}) \Vert < \alpha , \ \text{ then } \ {\varvec{u}}_\tau = {\varvec{w}}_\tau , \\[2mm] \quad \text{ if } \ \Vert {\varvec{\tau }}_\tau ({\varvec{u}}) \Vert = \alpha , \ \text{ then } \ {\varvec{u}}_{\tau } \not = {\varvec{w}}_\tau \ \text{ and } \ \displaystyle -{\varvec{\tau }}_\tau ({\varvec{u}}) = \alpha \, \frac{{\varvec{u}}_\tau -{\varvec{w}}_\tau }{\Vert {\varvec{u}}_\tau -{\varvec{w}}_\tau \Vert } \end{array} \right. \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$

where \({\varvec{w}}_\tau \) denotes the tangential velocity of the wall surface \(\Gamma _1\) can be found in [20, 22].

(iv) The nonlinear Navier–Fujita slip condition

$$\begin{aligned} \left\{ \begin{array}{ll} \Vert {\varvec{\tau }}_\tau ({\varvec{u}})\Vert \le \alpha + h(\Vert {\varvec{u}}_\tau - {\varvec{w}}_\tau \Vert )\\ {\varvec{\tau }}_\tau ({\varvec{u}}) \cdot ({\varvec{u}}_\tau - {\varvec{w}}_\tau ) = -(\alpha + h(\Vert {\varvec{u}}_\tau - {\varvec{w}}_\tau \Vert )) \, \Vert {\varvec{u}}_\tau - {\varvec{w}}_\tau \Vert \end{array} \right. \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$
(14)

where \(\alpha :\Gamma _1 \rightarrow (0,\infty )\) and \(h :\Gamma _1 \times [0,\infty ) \rightarrow [0, \infty )\) are given bounded functions such that for a.e. \({\varvec{x}}\in \Gamma _1\), \(h({\varvec{x}}, r) = 0\) if and only if \(r = 0\), has been studied in [21, 22]. This condition is a generalization of Navier’s slip condition and a Coulomb-type friction condition. Condition (14) says that the fluid slips at the boundary points if the magnitude of the tangential traction reaches a critical value \(\alpha \) being a prescribed threshold on the boundary wall, and when the slip occurs, the tangential traction equals to a given nonlinear function of the slip velocity.

We observe that the Navier–Fujita law of the form

$$\begin{aligned} \left\{ \begin{array}{ll} \Vert {\varvec{\tau }}_\tau ({\varvec{u}})\Vert \le F(\Vert {\varvec{u}}_\tau \Vert ),\\ \displaystyle -{\varvec{\tau }}_\tau ({\varvec{u}}) = F(\Vert {\varvec{u}}_\tau \Vert ) \, \frac{{\varvec{u}}_\tau -{\varvec{w}}_\tau }{\Vert {\varvec{u}}_\tau -{\varvec{w}}_\tau \Vert } \ \ \text{ if } \ \ {\varvec{u}}_\tau \not = {\varvec{w}}_\tau , \end{array} \right. \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$
(15)

where \(F({\varvec{x}},r) = \alpha ({\varvec{x}}) + h({\varvec{x}}, r)\) for \({\varvec{x}}\in \Gamma _1\), \(r \in {\mathbb {R}}\) can be formulated in the subdifferential form

$$\begin{aligned} -{\varvec{\tau }}_\tau ({\varvec{u}}) \in F(\Vert {\varvec{u}}_\tau \Vert ) \, \partial j ({\varvec{u}}_\tau ) \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$
(16)

where \(\partial j\) represents the subdifferential of the convex function \(j :{\mathbb {R}}^d \rightarrow {\mathbb {R}}\), \(j({\varvec{\xi }}) = \Vert {\varvec{\xi }}\Vert \) for \({\varvec{\xi }}\in {\mathbb {R}}^d\). In fact, for all \({\varvec{\xi }}\in {\mathbb {R}}^d\), we have

$$\begin{aligned}&-{\varvec{\tau }}_\tau ({\varvec{u}}) \cdot ({\varvec{\xi }}- ({\varvec{u}}_\tau -{\varvec{w}}_\tau )) = F(\Vert {\varvec{u}}_\tau \Vert ) \, \frac{{\varvec{u}}_\tau -{\varvec{w}}_\tau }{\Vert {\varvec{u}}_\tau -{\varvec{w}}_\tau \Vert } \cdot ({\varvec{\xi }}- ({\varvec{u}}_\tau -{\varvec{w}}_\tau ))\\&\quad = F (\Vert {\varvec{u}}_\tau \Vert ) \ \frac{{\varvec{\xi }}\cdot ({\varvec{u}}_\tau -{\varvec{w}}_\tau )}{\Vert {\varvec{u}}_\tau -{\varvec{w}}_\tau \Vert } - F(\Vert {\varvec{u}}_\tau \Vert ) \, \Vert {\varvec{u}}_\tau -{\varvec{w}}_\tau \Vert \le F(\Vert {\varvec{u}}_\tau \Vert )\, (\Vert {\varvec{\xi }}\Vert - \Vert {\varvec{u}}_\tau -{\varvec{w}}_\tau \Vert ) \end{aligned}$$

at the points on \(\Gamma _1\) where \({\varvec{u}}_\tau \not ={\varvec{w}}_\tau \), and

$$\begin{aligned} -{\varvec{\tau }}_\tau ({\varvec{u}}) \cdot ({\varvec{\xi }}- ({\varvec{u}}_\tau -{\varvec{w}}_\tau )) = -{\varvec{\tau }}_\tau ({\varvec{u}}) \cdot {\varvec{\xi }}\le \Vert {\varvec{\tau }}_\tau ({\varvec{u}})\Vert \, \Vert {\varvec{\xi }}\Vert \le F(\Vert {\varvec{u}}_\tau \Vert )\, (\Vert {\varvec{\xi }}\Vert - \Vert {\varvec{u}}_\tau -{\varvec{w}}_\tau \Vert ) \end{aligned}$$

at the points on \(\Gamma _1\) where \({\varvec{u}}_\tau ={\varvec{w}}_\tau \). Since \(F(\cdot ) > 0\), one gets

$$\begin{aligned} -\frac{{\varvec{\tau }}_\tau ({\varvec{u}})}{F (\Vert {\varvec{u}}_\tau \Vert )} \cdot ({\varvec{\xi }}- ({\varvec{u}}_\tau -{\varvec{w}}_\tau )) \le j({\varvec{\xi }}) - j({\varvec{u}}_\tau -{\varvec{w}}_\tau ) \ \ \text{ for } \text{ all } \ \ {\varvec{\xi }}\in {\mathbb {R}}^d \end{aligned}$$

which implies (16).

(v) Note that condition (5) is much more general than the Navier–Fujita slip condition (15) and leads to novel models. It allows to model nonmonotone slip boundary conditions of frictional type governed by nonconvex locally Lipschitz potentials. For example, let \(j :{\mathbb {R}}^d \rightarrow {\mathbb {R}}\) be defined by

$$\begin{aligned} j({\varvec{\xi }}) = \Vert {\varvec{\xi }}\Vert - \theta \, e^{-\Vert {\varvec{\xi }}\Vert } \ \ \text{ for } \ \ {\varvec{\xi }}\in {\mathbb {R}}^d, \end{aligned}$$

where \(\theta \ge 0\) is a given constant. For simplicity, here we take \({\varvec{w}}_\tau = 0\) on \(\Gamma _1\). The function j is nonconvex for \(\theta > 0\) and its generalized gradient is given by

$$\begin{aligned} \partial j({\varvec{\xi }}) = {\left\{ \begin{array}{ll} \displaystyle \ \left( 1 + \theta \, e^{-\Vert {\varvec{\xi }}\Vert } \right) \frac{{\varvec{\xi }}}{\Vert {\varvec{\xi }}\Vert } \ \ \ \ \text{ if } \ \ {\varvec{\xi }}\not = 0,\\ \ (\theta +1)\, \mathbf{B}(0, 1) \ \ \ \ \ \ \ \ \, \text{ if } \ \ {\varvec{\xi }}= 0, \end{array}\right. } \ \ \text{ on } \ \ \Gamma _1, \end{aligned}$$
(17)

where \(\mathbf{B}(0, 1)\) is the closed unit ball in \({\mathbb {R}}^d\). This choice of j leads to the slip condition of the form

$$\begin{aligned} \left\{ \begin{array}{ll} \Vert {\varvec{\tau }}_\tau ({\varvec{u}})\Vert \le (1 + \theta )\, k({\varvec{u}}_\tau ), \\ \displaystyle -{\varvec{\tau }}_\tau ({\varvec{u}}) = k({\varvec{u}}_\tau ) \, \big (1 + \theta \, e^{-\Vert {\varvec{u}}_\tau \Vert }\big )\, \frac{{\varvec{u}}_\tau }{\Vert {\varvec{u}}_\tau \Vert } \ \ \text{ if } \ \ {\varvec{u}}_\tau \not = 0 \end{array} \right. \ \ \text{ on } \ \ \Gamma _1 \end{aligned}$$

This condition illustrates the slip weakening phenomenon in which the tangential traction is a decreasing function of the tangential velocity. It is clear that for \(\theta =0\) we get (15) with \(k({\varvec{\xi }}) = F(\Vert {\varvec{\xi }}\Vert )\) and that j satisfies H(j) with \(c_0(x) = 1 + \theta \) and \(c_1 = 0\), so for this potential the smallness condition (25) below holds trivially. It is obvious that we can deal with a slip law (17) in which \(\theta \in L^\infty (\Gamma _1)\), \(\theta \ge 0\) a.e. on \(\Gamma _1\). This allows for varying the slip condition on different parts of the boundary \(\Gamma _1\).

Example 2

Let \(d=2\) and \({\varvec{\tau }}\) be the (only one) tangential direction on the boundary. Then \({\varvec{\tau }}_\tau = \tau _\tau {\varvec{\tau }}\) and \({\varvec{u}}_\tau = u_\tau {\varvec{\tau }}\), where \(\tau _\tau \) and \(u_\tau \) are real-valued functions. Consider the nonconvex locally Lipschitz function \(j :{\mathbb {R}}\rightarrow {\mathbb {R}}\) given by

$$\begin{aligned} j(r) = {\left\{ \begin{array}{ll} \displaystyle -e^{r} + \big (1-\frac{1}{e}\big ) r +2 &{}{\mathrm{if}} \, r < -1, \\ |r| &{}{\mathrm{if}} \, r \in [-1, 1], \\ \displaystyle -e^{-r} + \big (1-\frac{1}{e} \big ) r + \frac{2}{e} &{}{\mathrm{if}} \, r > 1. \end{array}\right. } \end{aligned}$$

Then the law (5) reduces to the nonmonotone slip condition

$$\begin{aligned} {-\tau }_{\tau } ({\varvec{u}}) = k({u}_{\tau }) \times {\left\{ \begin{array}{ll} \big (e^{-|u_\tau |} + 1 - \frac{1}{e}\big )\, {\hbox {sgn}}(u_\tau ) &{}{\hbox {if}} \, |u_\tau | >1, \\ {[}-1, 1{]} &{}{\hbox {if}} \, u_\tau = 0, \\ {\hbox {sgn}}(u_\tau ) &{}{\hbox {if}} \, u_\tau \in [-1,0) \cup (0, 1], \end{array}\right. } \end{aligned}$$

where \({\mathrm{sgn}}(r) = 1\) if \(r > 0,\) and \(-1\) if \(r < 0\). Thus, on \([-1, 1]\) the slip condition is represented by a monotone graph, and for \(u_\tau < -1\) and \(u_\tau > 1\) the tangential traction is a decreasing function of the tangential velocity.

To provide the variational formulation, we suppose that \({\varvec{u}}\), \({\mathbb {S}}\) and p are sufficiently smooth functions that satisfy Problem 1. Let \(\varvec{{v}}\in V\) be a test function. We multiply (1) by \(\varvec{{v}}- {\varvec{u}}\) and apply the Green-type formula, see [26, Theorems 2.25] to obtain

$$\begin{aligned}&\int _{\Omega } {\mathbb {S}}({{\mathbb {D}}}{\varvec{u}}):{{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx - \int _{\Gamma } ({\mathbb {S}}{\varvec{\nu }})\cdot (\varvec{{v}}- {\varvec{u}})\, d\Gamma - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}(\varvec{{v}}-{\varvec{u}}) \, dx \nonumber \\&\quad + \int _{\Gamma } ({\varvec{u}}\otimes {\varvec{u}}) {\varvec{\nu }}\cdot (\varvec{{v}}-{\varvec{u}}) \, d\Gamma + \int _\Omega \nabla p \, (\varvec{{v}}-{\varvec{u}}) \, dx = \int _\Omega {\varvec{f}}\cdot (\varvec{{v}}-{\varvec{u}}) \, dx . \end{aligned}$$
(18)

Since \(\mathop {\mathrm{div}}\nolimits \varvec{{v}}= \mathop {\mathrm{div}}\nolimits {\varvec{u}}= 0\) in \(\Omega \) and \(v_\nu = u_\nu = 0\) on \(\Gamma \), by [26, Theorem 2.24], we easily get

$$\begin{aligned} \int _\Omega \nabla p \, (\varvec{{v}}-{\varvec{u}}) \, dx = - \int _{\Omega } p \, \mathop {\mathrm{div}}\nolimits (\varvec{{v}}-{\varvec{u}}) \, dx + \int _{\Gamma } p \, (v_\nu - u_\nu ) \, d\Gamma = 0. \end{aligned}$$

Using the above equality and the properties

$$\begin{aligned} {\varvec{u}}= \varvec{{v}}= 0 \ \ \text{ on } \ \ \Gamma _0 \ \ \ \text{ and } \ \ \ ({\varvec{u}}\otimes {\varvec{u}}) {\varvec{\nu }}= 0 \ \ \text{ on } \ \ \Gamma _1 \end{aligned}$$
(19)

in (18), we have

$$\begin{aligned}&\int _{\Omega } {\mathbb {S}}({{\mathbb {D}}}{\varvec{u}}):{{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}(\varvec{{v}}-{\varvec{u}}) \, dx \nonumber \\&\qquad - \int _{\Gamma } ({\mathbb {S}}{\varvec{\nu }})\cdot (\varvec{{v}}- {\varvec{u}})\, d\Gamma = \int _\Omega {\varvec{f}}\cdot (\varvec{{v}}-{\varvec{u}}) \, dx . \end{aligned}$$
(20)

From the condition (5), H(k)(iii) and the definition of the subgradient, it follows that \(-{\varvec{\tau }}_\tau ({\varvec{u}}) = k({\varvec{u}}_\tau ) {\varvec{\xi }}\) and \({\varvec{\xi }}\in \partial j({\varvec{u}}_\tau )\) with \({\varvec{\xi }}\cdot (\varvec{{v}}_\tau -{\varvec{u}}_\tau ) \le j^0 ({\varvec{u}}_\tau ; \varvec{{v}}_\tau -{\varvec{u}}_\tau )\) on \(\Gamma _1\), and consequently

$$\begin{aligned} -{\varvec{\tau }}_\tau ({\varvec{u}}) \cdot (\varvec{{v}}_\tau -{\varvec{u}}_\tau ) \le k({\varvec{u}}_\tau ) j^0 ({\varvec{u}}_\tau ; \varvec{{v}}_\tau -{\varvec{u}}_\tau ) \ \ \text{ on } \ \ \Gamma _1. \end{aligned}$$

Combining the well known decomposition formula, see [26, (6.33)] and the property \({\mathbb {S}}_\tau ({\varvec{u}}) = {\varvec{\tau }}_\tau ({\varvec{u}})\) on \(\Gamma \), we have

$$\begin{aligned}&- \int _{\Gamma } ({\mathbb {S}}{\varvec{\nu }})\cdot (\varvec{{v}}- {\varvec{u}})\, d\Gamma = - \int _{\Gamma _1} ({\mathbb {S}}_\nu (v_\nu -u_\nu ) + {\mathbb {S}}_\tau \cdot (\varvec{{v}}_\tau - {\varvec{u}}_\tau ))\, d\Gamma \nonumber \\&\quad = -\int _{\Gamma _1} {\varvec{\tau }}_\tau ({\varvec{u}}) \cdot (\varvec{{v}}_\tau - {\varvec{u}}_\tau ) \, d\Gamma \le \int _{\Gamma _1} k({\varvec{u}}_\tau ) j^0 ({\varvec{u}}_\tau ; \varvec{{v}}_\tau -{\varvec{u}}_\tau )\, d\Gamma . \end{aligned}$$
(21)

Let \(\Omega _+ = \{ {\varvec{x}}\in \Omega \mid \Vert {\mathbb {D}} {\varvec{u}}({\varvec{x}}) \Vert > 0 \}\) and \(\Omega _0 = \{ {\varvec{x}}\in \Omega \mid \Vert {\mathbb {D}} {\varvec{u}}({\varvec{x}}) \Vert = 0 \}\). Exploiting the constitutive law (2), by the Cauchy–Schwarz inequality, we get

$$\begin{aligned}&\int _{\Omega _+} {\mathbb {S}}({{\mathbb {D}}}{\varvec{u}}):{{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx = \int _{\Omega _+} \left( \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}+ g \, \frac{{\mathbb {D}}{\varvec{u}}}{\Vert {\mathbb {D}}{\varvec{u}}\Vert } \right) : {{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx \nonumber \\&\quad \le \int _{\Omega _+} \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx + \int _{\Omega _+} g \, \frac{{\mathbb {D}}{\varvec{u}}}{\Vert {\mathbb {D}}{\varvec{u}}\Vert } : {{\mathbb {D}}}\varvec{{v}}\, dx - \int _{\Omega _+} g \, \Vert {\mathbb {D}}{\varvec{u}}\Vert \, dx \nonumber \\&\quad \le \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx + \int _{\Omega _+} g \, \Vert {{\mathbb {D}}}\varvec{{v}}\Vert \, dx - \int _{\Omega } g \, \Vert {\mathbb {D}}{\varvec{u}}\Vert \, dx . \end{aligned}$$
(22)

On the other hand, by (2), we have

$$\begin{aligned} \int _{\Omega _0} {\mathbb {S}}({{\mathbb {D}}}{\varvec{u}}):{{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx= & {} \int _{\Omega _0} {\mathbb {S}}({{\mathbb {D}}}{\varvec{u}}):{{\mathbb {D}}}\varvec{{v}}\, dx\nonumber \\\le & {} \int _{\Omega _0} \Vert {\mathbb {S}}({{\mathbb {D}}}{\varvec{u}})\Vert \, \Vert {{\mathbb {D}}}\varvec{{v}}\Vert \, dx \le \int _{\Omega _0} g \, \Vert {{\mathbb {D}}}\varvec{{v}}\Vert \, dx. \end{aligned}$$
(23)

Adding the inequalities (22) and (23), we deduce

$$\begin{aligned} \int _{\Omega } {\mathbb {S}}({{\mathbb {D}}}{\varvec{u}}):{{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx \le \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx + \int _{\Omega } g \, (\Vert {{\mathbb {D}}}\varvec{{v}}\Vert - \Vert {\mathbb {D}}{\varvec{u}}\Vert ) \, dx . \nonumber \\ \end{aligned}$$
(24)

Finally, we use (21) and (24) in (20) and obtain the following variational formulation of Problem 1.

Problem 3

Find a velocity \({\varvec{u}}\in V\) such that

$$\begin{aligned}&\int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} (\varvec{{v}}{-}{\varvec{u}}) \, dx {-} \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}(\varvec{{v}}-{\varvec{u}}) \, dx {+} \int _\Omega g \, (\Vert {\mathbb {D}}\varvec{{v}}\Vert {-} \Vert {\mathbb {D}} {\varvec{u}}\Vert )\, dx \\&\quad +\int _{\Gamma _1} k({\varvec{u}}_\tau ) j^0({\varvec{u}}_\tau ; \varvec{{v}}_\tau -{\varvec{u}}_\tau ) \, d\Gamma \ge \int _\Omega {\varvec{f}}\cdot (\varvec{{v}}-{\varvec{u}}) \, dx \ \ \text{ for } \text{ all } \ \ \varvec{{v}}\in V. \end{aligned}$$

Theorem 4

Under the hypotheses \(H(\mu )\), H(fg), H(j) and H(k), and the smallness condition

$$\begin{aligned} \mu _0 > k_1 c_1 \Vert \gamma \Vert ^2, \end{aligned}$$
(25)

Problem 3 has a solution.

The inequality in Problem 3 is a variational–hemivariational inequality. It combines the convex potential \(\varvec{{v}}\mapsto \int _\Omega g \, \Vert {\mathbb {D}} \varvec{{v}}\Vert \, dx \) and the locally Lipschitz, in general nonconvex, superpotential j. When \(j({\varvec{x}}, \cdot )\) is a convex function for a.e. \({\varvec{x}}\in \Gamma _1\), then from Problem 3 we obtain a purely elliptic variational inequality of second kind: find a velocity \({\varvec{u}}\in V\) such that

$$\begin{aligned}&\int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}(\varvec{{v}}-{\varvec{u}}) \, dx \\&\quad + \, \varrho ({\varvec{u}}, \varvec{{v}}) - \varrho ({\varvec{u}}, {\varvec{u}}) \ge \int _\Omega {\varvec{f}}\cdot (\varvec{{v}}-{\varvec{u}}) \, dx \ \ \text{ for } \text{ all } \ \ \varvec{{v}}\in V \end{aligned}$$

with \(\varrho :V \times V \rightarrow {\mathbb {R}}\) defined by

$$\begin{aligned} \varrho ({\varvec{u}}, \varvec{{v}}) = \int _\Omega g \, \Vert {\mathbb {D}}\varvec{{v}}\Vert \, dx + \int _{\Gamma _1} k({\varvec{u}}_\tau ) j(\varvec{{v}}_\tau ) \, d\Gamma \ \ \text{ for } \ \ {\varvec{u}}, \varvec{{v}}\in V. \end{aligned}$$

Further, for instance, if \(j({\varvec{x}}, {\varvec{\xi }}) = \Vert {\varvec{\xi }}\Vert \), as in the nonlinear Navier–Fujita slip condition (15), the condition H(j) holds with \(c_0 = 1\) and \(c_1 = 0\), so the smallness condition (25) is automatically satisfied and can be removed from hypotheses of Theorem 4. The uniqueness of a solution to Problem 3 under the hypotheses of Theorem 4 is an interesting and challenging open problem.

4 Proof of Theorem 4

The proof is based on an application of a surjectivity result for pseudomonotone and coercive operators, see [10, Theorem 1.3.70]. We divide the proof into several steps.

Step 1 Let \(X = L^2(\Gamma _1;{\mathbb {R}}^d)\). We introduce the following operators and functions

$$\begin{aligned}&A :V \rightarrow V^*, \ \ \langle A{\varvec{u}}, \varvec{{v}}\rangle = \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert )\, {\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} \varvec{{v}}\, dx, \ \ {\varvec{u}}, \varvec{{v}}\in V, \end{aligned}$$
(26)
$$\begin{aligned}&B :V \rightarrow V^*, \ \ \langle B {\varvec{u}}, \varvec{{v}}\rangle = - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}(\varvec{{v}}) \, dx, \ \ {\varvec{u}}, \varvec{{v}}\in V, \end{aligned}$$
(27)
$$\begin{aligned}&J :X \times X \rightarrow {\mathbb {R}}, \ \ J({\varvec{w}}, {\varvec{u}}) = \int _{\Gamma _1} k({\varvec{w}}) j({\varvec{u}})\, d\Gamma , \ \ {\varvec{w}}, {\varvec{u}}\in X, \end{aligned}$$
(28)
$$\begin{aligned}&\varphi :V \rightarrow {\mathbb {R}}, \ \ \varphi (\varvec{{v}}) = \int _\Omega g \, \Vert {\mathbb {D}} \varvec{{v}}\Vert \, dx, \ \ \varvec{{v}}\in V, \end{aligned}$$
(29)
$$\begin{aligned}&{{{\widetilde{{\varvec{f}}}}}} \in V^*, \ \ \langle {{{\widetilde{{\varvec{f}}}}}}, \varvec{{v}}\rangle = \int _\Omega {\varvec{f}}\cdot \varvec{{v}}\, dx, \ \ \varvec{{v}}\in V, \nonumber \\&M :V \rightarrow X, \ \ M \varvec{{v}}= \varvec{{v}}_\tau , \ \ \varvec{{v}}\in V. \end{aligned}$$
(30)

Under the above notation, the variational–hemivariational inequality in Problem 3 can be equivalently written as follows: find \({\varvec{u}}\in V\) such that

$$\begin{aligned} \langle A{\varvec{u}}+ B{\varvec{u}}- {{{\widetilde{{\varvec{f}}}}}}, \varvec{{v}}- {\varvec{u}}\rangle +\int _{\Gamma _1} k({\varvec{u}}_\tau ) j^0({\varvec{u}}_\tau ; \varvec{{v}}_\tau -{\varvec{u}}_\tau ) \, d\Gamma + \varphi (\varvec{{v}}) - \varphi ({\varvec{u}}) \ge 0 \end{aligned}$$
(31)

for all \(\varvec{{v}}\in V\). In order to show existence of solution to (31), we consider the auxiliary inclusion problem.

Problem 5

Find \({\varvec{u}}\in V\) such that

$$\begin{aligned} A{\varvec{u}}+ B{\varvec{u}}+ M^* \partial J(M{\varvec{u}}, M{\varvec{u}}) + \partial \varphi ({\varvec{u}}) \ni {{{\widetilde{{\varvec{f}}}}}}. \end{aligned}$$

In Problem 5, the notation \(\partial \varphi \) denotes the subdifferential in the sense of convex analysis of \(\varphi \), and \(M^* :X^* \rightarrow V^*\) is the adjoint operator to M. The notation \(\partial J\) stands for the generalized subgradient of the function \(J({\varvec{w}}, \cdot )\). Under hypothesis H(j), we know, see [26, Theorem 3.47(i)–(iv)], that J is well defined on \(X \times X\), \(J({\varvec{w}}, \cdot )\) is Lipschitz on bounded subsets of X for all \({\varvec{w}}\in X\) (and hence also locally Lipschitz on X), and the following inequality holds

$$\begin{aligned} J^0({\varvec{w}}, {\varvec{u}}; \varvec{{v}}) \le \int _{\Gamma _1} k({\varvec{w}}) j^0({\varvec{u}}; \varvec{{v}}) \, d\Gamma , \ \ {\varvec{w}}, {\varvec{u}}, \varvec{{v}}\in X. \end{aligned}$$
(32)

We observe that every solution to Problem 5 is also a solution to problem (31). Indeed, let \({\varvec{u}}\in V\) be a solution to Problem 5, which means that

$$\begin{aligned} A {\varvec{u}}+ B {\varvec{u}}+ M^* {\varvec{\xi }}+ {\varvec{\eta }}= {{{\widetilde{{\varvec{f}}}}}}, \end{aligned}$$
(33)

where \({\varvec{\xi }}\in \partial J(M{\varvec{u}}, M{\varvec{u}})\) and \({\varvec{\eta }}\in \partial \varphi ({\varvec{u}})\). Hence, by the definitions of subgradients, we get

$$\begin{aligned} \left\{ \begin{array}{lll} \langle {\varvec{\xi }}, {\varvec{z}}\rangle _{X^* \times X} \le J^0(M{\varvec{u}}, M{\varvec{u}}; {\varvec{z}}) \ \ \text{ for } \text{ all } \ \ {\varvec{z}}\in X, \\[2mm] \langle {\varvec{\eta }}, \varvec{{v}}-{\varvec{u}}\rangle \le \varphi (\varvec{{v}}) - \varphi ({\varvec{u}}) \ \ \text{ for } \text{ all } \ \ \varvec{{v}}\in V. \end{array}\right. \end{aligned}$$
(34)

Let \(\varvec{{v}}\in V\). Then, by (33), we have

$$\begin{aligned} \langle A {\varvec{u}}+ B {\varvec{u}}, \varvec{{v}}-{\varvec{u}}\rangle + \langle {\varvec{\xi }}, M (\varvec{{v}}-{\varvec{u}}) \rangle + \langle {\varvec{\eta }}, \varvec{{v}}-{\varvec{u}}\rangle = \langle {{{\widetilde{{\varvec{f}}}}}}, \varvec{{v}}-{\varvec{u}}\rangle . \end{aligned}$$
(35)

We combine (35) with (34) and (32) to deduce (31). We conclude that \({\varvec{u}}\in V\) is a solution to the inequality (31).

Step 2 We will show that Problem 5 has a solution. To this end, we use the well known result from the theory of monotone operators, see [10, Theorem 1.3.70] which states that if X is a reflexive Banach space and \(T :X \rightarrow 2^{X^*}\) is a multivalued pseudomonotone and coercive operator, then T is surjective.

Let \(T :V \rightarrow 2^{V^*}\) be defined by

$$\begin{aligned} T{\varvec{u}}= A{\varvec{u}}+ B{\varvec{u}}+ M^* \partial J(M{\varvec{u}}, M{\varvec{u}}) + \partial \varphi ({\varvec{u}}) \ \ \text{ for } \ \ {\varvec{u}}\in V. \end{aligned}$$

We will verify the following properties.

\(\bullet \) The operator \(A :V \rightarrow V^*\) given by (26) is bounded, monotone and continuous. We use \(H(\mu )\)(i) and the Hölder inequality to get

$$\begin{aligned} \langle A{\varvec{u}}, \varvec{{v}}\rangle&= \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert )\, {\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} \varvec{{v}}\, dx \le \mu _1 \int _{\Omega } {\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} \varvec{{v}}\, dx \\&\quad \le \mu _1 \Vert {\mathbb {D}}{\varvec{u}}\Vert _{L^2(\Omega , {\mathbb {M}}^d)} \Vert {\mathbb {D}}\varvec{{v}}\Vert _{L^2(\Omega , {\mathbb {M}}^d)} = \mu _1 \Vert {\varvec{u}}\Vert _V \Vert \varvec{{v}}\Vert _V \ \ \text{ for } \text{ all } \ \ {\varvec{u}}, \varvec{{v}}\in V \end{aligned}$$

and

$$\begin{aligned} \Vert A {\varvec{u}}\Vert _{V^*} \le \mu _1 \Vert {\varvec{u}}\Vert _V. \end{aligned}$$

The latter implies that A maps bounded sets in V into the bounded sets in \(V^*\), i.e., A is a bounded operator. Next, from hypothesis \(H(\mu )\)(ii), we directly have

$$\begin{aligned} \langle A {\varvec{u}}- A \varvec{{v}}, {\varvec{u}}-\varvec{{v}}\rangle = \int _{\Omega } \left( \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert )\, {\mathbb {D}}{\varvec{u}}- \mu (\Vert {{\mathbb {D}}}\varvec{{v}}\Vert )\, {\mathbb {D}}\varvec{{v}}\right) : {{\mathbb {D}}} ({\varvec{u}}- \varvec{{v}}) \, dx \ge 0 \end{aligned}$$

for all \({\varvec{u}}, \varvec{{v}}\in V\), which means that A is a monotone operator. We use [10, Theorem 1.5.2] (Krasnoselskii’s theorem for the Nemytskii operators) together with \(H(\mu )\)(ii) to deduce that the operator A is bounded and continuous from V to \(V^*\). Since the operator A is bounded, monotone and hemicontinuous (being continuous), by [26, Theorem 3.69(i)], we infer that A is pseudomonotone.

\(\bullet \) The operator \(B :V \rightarrow V^*\) defined by (27) is bounded and completely continuous. We can calculate that

$$\begin{aligned}&\langle B{\varvec{u}}, \varvec{{v}}\rangle = - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}(\varvec{{v}}) \, dx = \sum _{i,j}^d \int _{\Omega } u_j \frac{\partial u_i}{\partial x_j} v_i \, dx\\&\qquad \qquad \qquad + \, \frac{1}{2} \int _{\Omega } (\mathop {\mathrm{div}}\nolimits {\varvec{u}}) ({\varvec{u}}\cdot \varvec{{v}}) \, dx - \int _{\Gamma } ({\varvec{u}}\cdot \varvec{{v}}) u_\nu \, d\Gamma \ \ \text{ for } \ \ {\varvec{u}}, \varvec{{v}}\in H^1(\Omega ;{\mathbb {R}}^d). \end{aligned}$$

Hence, we have

$$\begin{aligned} \langle B{\varvec{u}}, \varvec{{v}}\rangle = \sum _{i,j}^d \int _{\Omega } u_j \frac{\partial u_i}{\partial x_j} v_i \, dx \ \ \text{ for } \ \ {\varvec{u}}, \varvec{{v}}\in V. \end{aligned}$$

We use the fact that the embedding operator \(V \subset L^4(\Omega ;{\mathbb {R}}^d)\) is continuous and compact for \(d=2\), 3, see (8). This fact and the generalized Hölder inequality imply that the integral in B is well defined as a product of two functions in \(L^4(\Omega ;{\mathbb {R}}^d)\) and a function in \(L^2(\Omega ;{\mathbb {R}}^d)\). The boundedness of the operator B follows from the estimate

$$\begin{aligned} |\langle B{\varvec{u}},\varvec{{v}}\rangle | \le \Vert {\varvec{u}}\Vert _{L^4(\Omega ,{\mathbb {R}}^d)} \Vert \nabla {\varvec{u}}\Vert _{L^2(\Omega ,{\mathbb {M}}^d)} \Vert \varvec{{v}}\Vert _{L^4(\Omega ,{\mathbb {R}}^d)} \le K\Vert {\varvec{u}}\Vert _V^2\Vert \varvec{{v}}\Vert _V \end{aligned}$$

for all \({\varvec{u}}\), \(\varvec{{v}}\in V\) with a constant \(K>0\). Subsequently, we prove that B is completely continuous, that is, it maps weakly convergent sequences in V to strongly convergent ones in \(V^*\). Let \(\{ {\varvec{u}}^k \} \subset V\) be such that \({\varvec{u}}^k \rightharpoonup {\varvec{u}}\) in V, as \(k\rightarrow \infty \). Therefore, we have \({\varvec{u}}^k \rightarrow {\varvec{u}}\) in \(L^4(\Omega ; {\mathbb {R}}^d)\). From the estimate

$$\begin{aligned}&\Vert u_i^ku_j^k-u_i u_j\Vert _{L^2(\Omega )} = \Vert (u_i^k-u_i)(u_j^k-u_j) -u_j(u^k_i-u_i)-u_i(u^k_j-u_j)\Vert _{L^2(\Omega )} \\&\quad \le \Vert (u_i^k-u_i)(u_j^k-u_j)\Vert _{L^2(\Omega )} +\Vert u_j(u^k_i-u_i)\Vert _{L^2(\Omega )} +\Vert u_i(u^k_j-u_j)\Vert _{L^2(\Omega )}\rightarrow 0, \end{aligned}$$

it follows that

$$\begin{aligned} u^k_i u_j^k \rightarrow u_i u_j \ \ \text{ in } \ \ L^2(\Omega ) \ \ \text{ for } \text{ all } \ \ i, j=1,\ldots , d. \end{aligned}$$
(36)

Let \(\varvec{{v}}\in V\). We apply the Hölder inequality to obtain

$$\begin{aligned}&|\langle B{\varvec{u}}^k - B{\varvec{u}}, \varvec{{v}}\rangle | = |\int _\Omega \left( ({\varvec{u}}^k \otimes {\varvec{u}}^k) - ({\varvec{u}}\otimes {\varvec{u}}) \right) : {\mathbb {D}}(\varvec{{v}}) \, dx| \\&\quad \le | \sum _{i,j}^d \int _{\Omega } (u^k_i u^k_j - u_i u_j) \, {\mathbb {D}}_{ij}(\varvec{{v}}) \, dx| \le \sum _{i,j}^d \Vert u^k_i u^k_j - u_i u_j \Vert _{L^2(\Omega )} \, \Vert {\mathbb {D}}_{ij}(\varvec{{v}}) \Vert _{L^2(\Omega ;{\mathbb {M}}^d)} \\&\qquad \le c \, \Big ( \sum _{i,j}^d \Vert u^k_i u^k_j - u_i u_j \Vert _{L^2(\Omega )} \Big ) \, \Vert \varvec{{v}}\Vert _V \end{aligned}$$

with a constant \(c>0\) independent of k. Hence, by (36), \(\Vert B {\varvec{u}}^k - B {\varvec{u}}\Vert _{V^*} \rightarrow 0\), as \(k\rightarrow \infty \) which proves that B is the completely continuous operator.

Using, for instance, the characterization of pseudomonotonicity of [26, Proposition 3.66], we know that every bounded and completely continuous map is pseudomonotone. We deduce that \(B :V \rightarrow V^*\) is a pseudomonotone operator.

\(\bullet \) The sum of two pseudomonotone operators remains pseudomonotone, see, e.g., [26, Proposition 3.70], so the operator \(A+B :V \rightarrow V^*\) is pseudomonotone.

\(\bullet \) Consider now the multivalued operator \(T_1 :V \rightarrow 2^{V^*}\) defined by

$$\begin{aligned} T_1 {\varvec{u}}= A{\varvec{u}}+ B{\varvec{u}}+ M^* \partial J(M{\varvec{u}}, M{\varvec{u}}) \ \ \text{ for } \ \ {\varvec{u}}\in V. \end{aligned}$$

We will prove that \(T_1\) is a bounded, multivalued pseudomonotone operator.

First, we note that by hypotheses H(j), H(k), and [26, Theorem 3.47(v), (vi)], we obtain

$$\begin{aligned} \partial J({\varvec{w}}, {\varvec{u}}) \subseteq \int _{\Gamma _1} k({\varvec{w}}) \partial j({\varvec{u}}) \, d\Gamma \ \ \text{ for } \text{ all } \ \ {\varvec{w}}, {\varvec{u}}\in X, \end{aligned}$$

which entails

$$\begin{aligned} \Vert \partial J ({\varvec{w}}, {\varvec{u}}) \Vert _{X^*} \le c_3 + c_4 \Vert {\varvec{u}}\Vert _X \ \ \text{ for } \text{ all } \ \ {\varvec{w}}, {\varvec{u}}\in X, \end{aligned}$$
(37)

where \(c_3\), \(c_4 \ge 0\). In particular, this estimate shows that the operator \(\partial J:X \times X \rightarrow 2^{X^*}\) is bounded, recall that the subdifferential is taken with respect to the last variable of a function.

Second, we show that \(T_1\) is a bounded operator. Let \({{\mathcal {D}}}\) be a bounded subset of V. Since \(A+B\) and M are bounded operators, it is clear that \((A+B)\varvec{{v}}\) and \(M\varvec{{v}}\) for \(\varvec{{v}}\in {{\mathcal {D}}}\) are bounded in \(V^*\) and X, respectively. The boundedness of \(\partial J\) implies that \(\partial J (M\varvec{{v}}, M\varvec{{v}})\) for \(\varvec{{v}}\in {{\mathcal {D}}}\) remains in a bounded subset of \(X^*\), and subsequently \(M^* \partial J (M\varvec{{v}}, M\varvec{{v}})\) for \(\varvec{{v}}\in {{\mathcal {D}}}\) is bounded in \(V^*\). Thus \(T_1\) is a bounded operator as a sum of bounded operators.

Third, we note that \(T_1 \varvec{{v}}\) is a nonempty, closed and convex subset of \(V^*\) for all \(\varvec{{v}}\in V\). Indeed, we know, see [26, Proposition 2.23(iv)], that for all \({\varvec{w}}\in X\), the generalized gradient \(\partial J({\varvec{w}}, {\varvec{w}})\) is a nonempty, convex and weakly compact set in \(X^*\). Hence, this set is also weakly closed, and by convexity, it is closed in \(X^*\). Therefore, the set \(M^* \partial J (M\varvec{{v}}, M\varvec{{v}})\) is nonempty, closed and convex for all \(\varvec{{v}}\in V\), and the values of \(T_1\) are nonempty, closed and convex too.

Fourth, in order to prove that \(T_1\) is pseudomonotone, by Proposition [10, Proposition 1.3.66], it is enough to show that \(T_1\) is generalized pseudomonotone. To this end, let \(\varvec{{v}}_n \in V\), \(\varvec{{v}}_n \rightharpoonup \varvec{{v}}\) in V, \(\varvec{{v}}_n^* \in T_1 \varvec{{v}}_n\), \(\varvec{{v}}_n^* \rightharpoonup \varvec{{v}}^*\) in \(V^*\) and \(\limsup \langle \varvec{{v}}_n^*, \varvec{{v}}_n - \varvec{{v}}\rangle \le 0\). We will prove that \(\varvec{{v}}^* \in T_1 \varvec{{v}}\) and \(\langle \varvec{{v}}_n^*, \varvec{{v}}_n \rangle \rightarrow \langle \varvec{{v}}^*, \varvec{{v}}\rangle \). Since \(\varvec{{v}}_n^* \in T_1 \varvec{{v}}_n\), we have \(\varvec{{v}}_n^* = (A + B) \varvec{{v}}_n + M^* {\varvec{w}}_n\) with \({\varvec{w}}_n \in \partial J(M\varvec{{v}}_n, M\varvec{{v}}_n)\). It is obvious that \(\{ \varvec{{v}}_n \}\) is bounded in V and \(\{ M \varvec{{v}}_n \}\) is bounded in X. Using the boundedness of \(\partial J\), we obtain that \(\{ {\varvec{w}}_n \}\) remains in a bounded subset of \(X^*\). Hence, by passing to a subsequence, if necessary, we may assume that

$$\begin{aligned} {\varvec{w}}_n \rightharpoonup {\varvec{w}}\ \ \text{ in } \ \ X^* \ \text{ with } \ {\varvec{w}}\in X^*. \end{aligned}$$
(38)

On the other hand, by the compactness of operator M, see (9), we infer that

$$\begin{aligned} M \varvec{{v}}_n \rightarrow M \varvec{{v}}\ \ \text{ in } \ \ X. \end{aligned}$$
(39)

Exploiting (38), (39) and the closedness of the graph of \(\partial J\), see [26, Proposition 3.23(v)], from \({\varvec{w}}_n \in \partial J(M\varvec{{v}}_n, M\varvec{{v}}_n)\), we get \({\varvec{w}}\in \partial J(M\varvec{{v}}, M\varvec{{v}})\). Next, by passing to the \(\limsup \) in the following equality

$$\begin{aligned} \langle \varvec{{v}}_n^*, \varvec{{v}}_n -\varvec{{v}}\rangle = \langle (A+B) \varvec{{v}}_n, \varvec{{v}}_n - \varvec{{v}}\rangle + \langle {\varvec{w}}_n, M\varvec{{v}}_n - M \varvec{{v}}\rangle _{X^* \times X}, \end{aligned}$$

we get

$$\begin{aligned} \limsup \langle (A+B) \varvec{{v}}_n, \varvec{{v}}_n -\varvec{{v}}\rangle = \limsup \langle \varvec{{v}}_n^*, \varvec{{v}}_n - \varvec{{v}}\rangle \le 0. \end{aligned}$$

Since \(A+B\) is pseudomonotone, we have

$$\begin{aligned} (A+B) \varvec{{v}}_n \rightharpoonup (A+B) \varvec{{v}}\ \ \text{ in } \ V^* \ \ \text{ and } \ \ \langle (A+B) \varvec{{v}}_n, \varvec{{v}}_n -\varvec{{v}}\rangle \rightarrow 0. \end{aligned}$$

Finally, we are able to pass to the limit in the equality \(\varvec{{v}}_n^* = (A+B) \varvec{{v}}_n + M^* {\varvec{w}}_n\) to get

$$\begin{aligned} \varvec{{v}}^* = (A+B) \varvec{{v}}+ M^* {\varvec{w}}\in (A+B)\varvec{{v}}+ M^* \partial J(M\varvec{{v}}, M\varvec{{v}}) = T_1 \varvec{{v}}. \end{aligned}$$

Moreover, we obtain

$$\begin{aligned}&\langle \varvec{{v}}_n^*, \varvec{{v}}_n \rangle = \langle (A+B) \varvec{{v}}_n, \varvec{{v}}_n -\varvec{{v}}\rangle + \langle (A+B) \varvec{{v}}_n, \varvec{{v}}\rangle + \langle {\varvec{w}}_n, M \varvec{{v}}_n \rangle _{X^* \times X} \\&\qquad \qquad \quad \rightarrow \langle (A+B) \varvec{{v}}, \varvec{{v}}\rangle + \langle {\varvec{w}}, M \varvec{{v}}\rangle _{X^* \times X} = \langle (A+B) \varvec{{v}}+ M^* {\varvec{w}}, \varvec{{v}}\rangle = \langle \varvec{{v}}^*, \varvec{{v}}\rangle . \end{aligned}$$

This completes the proof that the operator \(T_1\) is pseudomonotone.

\(\bullet \) Consider now the multivalued operator \(T_2 :V \rightarrow 2^{V^*}\) defined by

$$\begin{aligned} T_2 {\varvec{u}}= \partial \varphi ({\varvec{u}}) \ \ \text{ for } \ \ {\varvec{u}}\in V \end{aligned}$$

with \(\varphi :V \rightarrow {\mathbb {R}}\) given by (29). We will prove that \(T_2\) is a bounded, multivalued pseudomonotone operator. Let \(I :L^2(\Omega ;{{\mathbb {M}}^d}) \rightarrow {\mathbb {R}}\) be the functional defined by

$$\begin{aligned} I({\varvec{w}}) = \int _{\Omega } g \, \Vert {\varvec{w}}\Vert _{{\mathbb {M}}^d}\, dx \ \ \text{ for } \ \ {\varvec{w}}\in L^2(\Omega ;{{\mathbb {M}}^d}). \end{aligned}$$

Note that, by H(fg), the integrand \(h:\Omega \times {\mathbb {M}}^d \rightarrow {\mathbb {R}}\), \(h({\varvec{x}}, {\varvec{w}}) = g({\varvec{x}}) \Vert {\varvec{w}}\Vert _{{\mathbb {M}}^d}\) has the following properties: \(h(\cdot , {\varvec{w}})\) is measurable for all \({\varvec{w}}\in {\mathbb {M}}^d\), \(h({\varvec{x}}, \cdot )\) is convex and continuous for a.e. \({\varvec{x}}\in \Omega \), and

$$\begin{aligned} |I({\varvec{w}}) | \le \Vert g \Vert _{L^2(\Omega )} \Vert {\varvec{w}}\Vert _{L^2(\Omega ;{{\mathbb {M}}^d})} \ \ \text{ for } \text{ all } \ \ {\varvec{w}}\in L^2(\Omega ;{{\mathbb {M}}^d}). \end{aligned}$$

We apply [3, Proposition 2.53] to deduce that the functional I is convex, lsc, everywhere finite, and

$$\begin{aligned} \partial I({\varvec{w}}) = \left\{ \, {\varvec{w}}^* \in L^2(\Omega ;{{\mathbb {M}}^d}) \mid {\varvec{w}}^*({\varvec{x}}) \in \partial h({\varvec{x}}, {\varvec{w}}({\varvec{x}})) \ \text{ a.e. } \ {\varvec{x}}\in \Omega \, \right\} \end{aligned}$$

for all \({\varvec{w}}\in L^2(\Omega ;{{\mathbb {M}}^d})\). The latter implies that \(\Vert \partial I({\varvec{w}}) \Vert _{L^2(\Omega ;{{\mathbb {M}}^d})} \le \Vert g \Vert _{L^2(\Omega )}\) for all \({\varvec{w}}\in L^2(\Omega ;{{\mathbb {M}}^d})\). Then \(\varphi ({\varvec{u}}) = I({\mathbb {D}}{\varvec{u}})\) for \({\varvec{u}}\in H^1(\Omega ;{\mathbb {R}}^d)\) is convex, lsc, and finite. Since \({\mathbb {D}}:H^1(\Omega ;{\mathbb {R}}^d) \rightarrow L^2(\Omega ;{{\mathbb {M}}^d})\) is linear and continuous, we use the chain rule to obtain

$$\begin{aligned} \Vert \partial \varphi ({\varvec{u}}) \Vert _{V^*} \le c_5 \, \Vert g \Vert _{L^2(\Omega )} \ \ \text{ for } \ \ {\varvec{u}}\in V \ \ \text{ with } \ \ c_5 > 0. \end{aligned}$$
(40)

From [10, Theorem 1.3.19] and (40), it follows that \(\partial \varphi :V \rightarrow 2^{V^*}\) is a maximal monotone and bounded operator, and \(D(\partial \varphi ) = V\). In consequence, by [10, Corollary 1.3.67], the operator \(T_2 :V \rightarrow 2^{V^*}\) is bounded and pseudomonotone.

\(\bullet \) Since \(T_1\), \(T_2 :V \rightarrow 2^{V^*}\) are bounded and pseudomonotone operators, from [10, Proposition 1.3.68], we deduce that \(T = T_1 + T_2\) is also bounded and pseudomonotone operator.

\(\bullet \) We will prove that \(T :V \rightarrow 2^{V^*}\) is coercive, i.e.,

$$\begin{aligned} \frac{\inf \{ \langle {\varvec{u}}^*, {\varvec{u}}\rangle \mid {\varvec{u}}^* \in T{\varvec{u}}\}}{\Vert {\varvec{u}}\Vert _V} \rightarrow +\infty , \ \ \text{ as } \ \ \Vert {\varvec{u}}\Vert _V \rightarrow \infty . \end{aligned}$$

Let \({\varvec{u}}\in V\). From hypothesis \(H(\mu )\)(i), we have

$$\begin{aligned} \langle A {\varvec{u}}, {\varvec{u}}\rangle = \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert )\, {\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} {\varvec{u}}\, dx \ge \mu _0 \int _{\Omega } \Vert {{\mathbb {D}}}{\varvec{u}}\Vert ^2 \, dx = \mu _0 \Vert {\varvec{u}}\Vert ^2_V. \end{aligned}$$
(41)

Next, we use the integration by parts formula, see [26, Theorem 2.23], and the conditions \(\mathop {\mathrm{div}}\nolimits {\varvec{u}}= 0\) in \(\Omega \), \({\varvec{u}}= 0\) on \(\Gamma _0\), and \(u_\nu = 0\) on \(\Gamma _1\) to get

$$\begin{aligned} \langle B {\varvec{u}}, {\varvec{u}}\rangle&= - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}({\varvec{u}}) \, dx = - \frac{1}{2} \sum _{i,j=1}^d \int _\Omega \left( u_i u_j \frac{\partial u_i}{\partial x_j} + u_i u_j \frac{\partial u_j}{\partial x_i} \right) \, dx \nonumber \\[1mm]&= - \frac{1}{2} \sum _{i,j=1}^d \int _\Omega \left( u_j \, \frac{\partial }{\partial x_j} \Big (\frac{u_i^2}{2}\Big ) + u_i \, \frac{\partial }{\partial x_i} \Big (\frac{u_j^2}{2}\Big )\right) \, dx\nonumber \\&= \frac{1}{2} \left( \int _\Omega (\mathop {\mathrm{div}}\nolimits {\varvec{u}}) \Vert {\varvec{u}}\Vert ^2 \, dx + \int _{\Gamma } u_\nu \Vert {\varvec{u}}\Vert ^2 \, d\Gamma \right) = 0. \end{aligned}$$
(42)

Further, we profit from H(j)(iii), H(k)(iii) and (32), to deduce that for every \({\varvec{\xi }}\in \partial J(M{\varvec{u}}, M{\varvec{u}})\), it holds

$$\begin{aligned}&\langle M^* {\varvec{\xi }}, {\varvec{u}}\rangle \\&\quad = \langle {\varvec{\xi }}, M{\varvec{u}}\rangle _{X^* \times X} \le J^0(M{\varvec{u}}, M{\varvec{u}}; M{\varvec{u}}) \le \int _{\Gamma _1} k({\varvec{u}}_\tau ) j^0({\varvec{u}}_\tau ; {\varvec{u}}_\tau )\, d\Gamma \\&\quad \le k_1 \int _{\Gamma _1} \left( c_0({\varvec{x}}) + c_1 \Vert {\varvec{u}}\Vert _{{\mathbb {R}}^d} \right) \Vert {\varvec{u}}\Vert _{{\mathbb {R}}^d} \, d\Gamma = k_1 \Vert c_0\Vert _{L^2(\Gamma _1)} \Vert {\varvec{u}}\Vert _X + k_1 c_1 \Vert {\varvec{u}}\Vert ^2_X. \end{aligned}$$

Recalling that the trace operator \(\gamma :V \rightarrow X\) is continuous, we arrive at

$$\begin{aligned} |\langle M^* \partial J(M{\varvec{u}}, M{\varvec{u}}), {\varvec{u}}\rangle | \le k_1 \Vert c_0\Vert _{L^2(\Gamma _1)} \Vert \gamma \Vert \Vert {\varvec{u}}\Vert _V + k_1 c_1 \Vert \gamma \Vert ^2 \Vert {\varvec{u}}\Vert ^2_V. \end{aligned}$$
(43)

Also, for the convex subdifferential term, we use (40) to obtain

$$\begin{aligned} |\langle \partial \varphi ({\varvec{u}}), {\varvec{u}}\rangle | \le \Vert \partial \varphi ({\varvec{u}}) \Vert _{V^*} \Vert {\varvec{u}}\Vert _V \le c_5 \, \Vert g\Vert _{L^2(\Omega )} \Vert {\varvec{u}}\Vert _V. \end{aligned}$$
(44)

Combining (41), (42), (43) and (44), we have

$$\begin{aligned} \langle T {\varvec{u}}, {\varvec{u}}\rangle \ge (\mu _0 - k_1 c_1 \Vert \gamma \Vert ^2 ) \Vert {\varvec{u}}\Vert ^2_V - k_1 \Vert c_0\Vert _{L^2(\Gamma _1)} \Vert \gamma \Vert \Vert {\varvec{u}}\Vert _V - c_5 \, \Vert g\Vert _{L^2(\Omega )} \Vert {\varvec{u}}\Vert _V \end{aligned}$$

for all \({\varvec{u}}\in V\). Taking into account the smallness condition (25), we deduce that \(T :V \rightarrow 2^{V^*}\) is a coercive operator.

In conclusion, the operator \(T :V \rightarrow 2^{V^*}\) is pseudomonotone and coercive, and so it is surjective. Therefore, for each \({{{\widetilde{{\varvec{f}}}}}} \in V^*\), Problem 5 has a solution. The proof is now complete. \(\square \)

The existence results for Problem 3 in Theorem 4 extends [24, Theorem 4.1] obtained for the Oseen model under more restrictive (Clarke regularity) hypothesis on \(j({\varvec{x}}, \cdot )\) and on the viscosity function \(\mu \). Further, we have also neglected the hypothesis on the relaxed monotonicity of \(\partial j({\varvec{x}}, \cdot )\) for a.e. \({\varvec{x}}\in \Gamma _1\) which has been widely used for hemivariational inequalites, see [26, 27] and the references therein.

We state a version of Theorem 4 which provides a sufficient condition for existence of solution without the smallness condition (25).

Theorem 6

Assume the hypotheses \(H(\mu )\), H(fg), H(k), H(j) (i), (ii), and H(j) (iii) is replaced by

$$\begin{aligned} \Vert \partial j({\varvec{x}},{\varvec{\xi }}) \Vert _{{\mathbb {R}}^d} \le c_0 ({\varvec{x}}) + c_1 \, \Vert {\varvec{\xi }}\Vert ^r_{{\mathbb {R}}^d} \ \ \text{ for } \text{ all } \ \ {\varvec{\xi }}\in {\mathbb {R}}^d, \ \text{ a.e. } \ {\varvec{x}}\in \Gamma _1 \end{aligned}$$
(45)

with \(c_0 \in L^2(\Gamma _1)\), \(c_0\), \(c_1 \ge 0\), and \(0 \le r < 1\). Then Problem 3 has a solution.

Proof

Under the growth condition (45), the coercivity condition of the operator T reads as follows

$$\begin{aligned} \langle T {\varvec{u}}, {\varvec{u}}\rangle \ge \mu _0 \Vert {\varvec{u}}\Vert ^2_V - k_1 c_1 \Vert \gamma \Vert ^2 \Vert {\varvec{u}}\Vert ^{r+1}_V - k_1 \Vert c_0\Vert _{L^2(\Gamma _1)} \Vert \gamma \Vert \Vert {\varvec{u}}\Vert _V - c_5 \, \Vert g\Vert _{L^2(\Omega )} \Vert {\varvec{u}}\Vert _V \end{aligned}$$

for all \({\varvec{u}}\in V\). The rest of the proof is similar to the one of Theorem 4. \(\square \)

We complete this section by showing a possible passage from Problem 3 to Problem 1 in the case of regular solutions and \(g=0\).

Proposition 7

Let \({\varvec{u}}\in H^2(\Omega ; {\mathbb {R}}^d)\) be the solution to Problem 3 with \(g = 0\). Then there is \(p\in L^2(\Omega )\) with \(\int _{\Omega } p \, dx = 0\) such that \(({\varvec{u}}, p)\) satisfies conditions of Problem 1.

Proof

Let \({\varvec{u}}\in V\) be a smooth solution to Problem 3 for \(g = 0\). We recover the pressure in Problem 3. Let \({\varvec{\psi }}\in C^\infty _0(\Omega ;{\mathbb {R}}^d)\) be such that \(\mathop {\mathrm{div}}\nolimits {\varvec{\psi }}= 0\) in \(\Omega \). Define \(\varvec{{v}}:= {\varvec{u}}\pm {\varvec{\psi }}\). Then \(\varvec{{v}}\in V\), so we choose it as a test function in Problem 3 to get

$$\begin{aligned} \int _{\Omega } \big ( \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} {\varvec{\psi }}- ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}} {\varvec{\psi }}\big ) \, dx = \int _{\Omega } {\varvec{f}}\cdot {\varvec{\psi }}\, dx . \end{aligned}$$

This means that the linear functional \(L :H^1_0(\Omega ;{\mathbb {R}}^d) \rightarrow {\mathbb {R}}\) defined by

$$\begin{aligned} \langle L, {\varvec{\psi }}\rangle = \int _{\Omega } \big ( \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} {\varvec{\psi }}\, dx - ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}} {\varvec{\psi }}- {\varvec{f}}\cdot {\varvec{\psi }}\big ) \, dx \end{aligned}$$

satisfies \(L \in H^{-1}(\Omega ;{\mathbb {R}}^d)\) and \(\langle L, {\varvec{\psi }}\rangle = 0\) for all \({\varvec{\psi }}\in C^\infty _0(\Omega ;{\mathbb {R}}^d)\) with \(\mathop {\mathrm{div}}\nolimits {\varvec{\psi }}= 0\) in \(\Omega \). We are in a position to apply [35, Lemma 2.2.2, p.75] to deduce that that there exists a unique \(p \in L^p(\Omega )\) such that \(\int _{\Omega } p \, dx =0\) and \(L = \nabla p\) in the sense of distributions. By the definition of the divergence operator, we have

$$\begin{aligned} \langle L, {\varvec{\psi }}\rangle = \int _{\Omega } \left( -\mathop {\mathrm{Div}}\nolimits \big ( \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}- ({\varvec{u}}\otimes {\varvec{u}}) \big ) + \nabla p - {\varvec{f}}\right) \cdot {\varvec{\psi }}\, dx = 0. \end{aligned}$$

Since \(\psi \) is arbitrary, by the variational lemma, see [37, Proposition 18.6], we deduce

$$\begin{aligned} - \mathop {\mathrm{Div}}\nolimits \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}+ \mathop {\mathrm{Div}}\nolimits ({\varvec{u}}\otimes {\varvec{u}}) + \nabla p = {\varvec{f}}\ \ \text{ in } \ \ \Omega , \end{aligned}$$
(46)

which proves that (1) with condition (2) holds (recall that, for convenience, we have assumed that \(\rho \equiv 1\)).

Subsequently, we interpret the boundary condition (5). We multiply the equation (46) by a suitable test function and compare the result with the inequality in Problem 3. Let \({\varvec{\psi }}\in C^\infty (\Omega ;{\mathbb {R}}^d)\) with \(\mathop {\mathrm{div}}\nolimits {\varvec{\psi }}= 0\) in \(\Omega \), \({\varvec{\psi }}= 0\) on \(\Gamma _0\) and \(\psi _\nu = 0\) on \(\Gamma _1\). We choose \(\varvec{{v}}:= {\varvec{u}}+ {\varvec{\psi }}\), we use the properties of \({\varvec{\psi }}\) and the fact \({\varvec{u}}\in V\) to infer that \(\varvec{{v}}\in V\). We multiply (46) by \(\varvec{{v}}- {\varvec{u}}\) and get

$$\begin{aligned}&\int _{\Omega } (-\mathop {\mathrm{Div}}\nolimits \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}) \cdot {\varvec{\psi }}\, dx + \int _\Omega \mathop {\mathrm{Div}}\nolimits ({\varvec{u}}\otimes {\varvec{u}}) \cdot {\varvec{\psi }}\, dx + \int _{\Omega } \nabla p \cdot {\varvec{\psi }}\, dx\nonumber \\&= \int _{\Omega } {\varvec{f}}\cdot {\varvec{\psi }}\, dx. \end{aligned}$$
(47)

Analogously as in the weak formulation (see Sect. 3), we have

$$\begin{aligned}&\int _{\Omega } (-\mathop {\mathrm{Div}}\nolimits \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}) \cdot {\varvec{\psi }}\, dx\nonumber \\&\quad = \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {\mathbb {D}} {\varvec{\psi }}\, dx - \int _{\Gamma _1} {\varvec{\tau }}_\tau ({\varvec{u}}) \cdot {\varvec{\psi }}_\tau \, d\Gamma \end{aligned}$$
(48)

and

$$\begin{aligned} \int _\Omega \nabla p \cdot {\varvec{\psi }}\, dx = 0. \end{aligned}$$
(49)

Inserting (48) and (49) into (47), we obtain

$$\begin{aligned}&\int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {\mathbb {D}} {\varvec{\psi }}\, dx - \int _{\Gamma _1} {\varvec{\tau }}_\tau ({\varvec{u}}) \cdot {\varvec{\psi }}_\tau \, d\Gamma +\int _\Omega \mathop {\mathrm{Div}}\nolimits ({\varvec{u}}\otimes {\varvec{u}}) \cdot {\varvec{\psi }}\, dx\nonumber \\&\quad = \int _{\Omega } {\varvec{f}}\cdot {\varvec{\psi }}\, dx. \end{aligned}$$
(50)

On the other hand, Problem 3 for \(\varvec{{v}}= {\varvec{u}}+ {\varvec{\psi }}\) implies

$$\begin{aligned} \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {\mathbb {D}} {\varvec{\psi }}\, dx - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}} {\varvec{\psi }}\, dx + \int _{\Gamma _1} k({\varvec{u}}_\tau ) j^0({\varvec{u}}_\tau ; {\varvec{\psi }}_\tau ) \, d\Gamma \ge \int _\Omega {\varvec{f}}\cdot {\varvec{\psi }}\, dx.\nonumber \\ \end{aligned}$$
(51)

Combining (50) with (51), we get

$$\begin{aligned} \int _{\Gamma _1} -{\varvec{\tau }}_\tau ({\varvec{u}}) \cdot {\varvec{\psi }}_\tau \, d\Gamma \le \int _{\Gamma _1} k({\varvec{u}}_\tau ) j^0({\varvec{u}}_\tau ; {\varvec{\psi }}_\tau ) \, d\Gamma . \end{aligned}$$

Note that \(\int _\Omega \mathop {\mathrm{Div}}\nolimits ({\varvec{u}}\otimes {\varvec{u}}) \cdot {\varvec{\psi }}\, dx= - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}{\varvec{\psi }}\, dx\). Since \({\varvec{\psi }}\) is arbitrary, choosing \({\varvec{\psi }}_\tau = {\varvec{\xi }}\) on \(\Gamma _1\) with \({\varvec{\xi }}\in {\mathbb {R}}^d\), we obtain

$$\begin{aligned} -\frac{{\varvec{\tau }}_\tau ({\varvec{u}})}{k({\varvec{u}}_\tau )} \cdot {\varvec{\xi }}\le j^0({\varvec{u}}_\tau ; {\varvec{\xi }}) \ \ \text{ on } \ \ \Gamma _1 \ \ \text{ for } \text{ all } \ \ {\varvec{\xi }}\in {\mathbb {R}}^d . \end{aligned}$$

We use the definition of the generalized gradient to have \(-{\varvec{\tau }}_\tau ({\varvec{u}}) \in k({\varvec{u}}_\tau ) \partial j({\varvec{u}}_\tau )\) on \(\Gamma _1\). Since \({\varvec{u}}\in V\), we know that \(u_\nu =0\) on \(\Gamma _1\). This proves the boundary condition (5) and completes the proof. \(\square \)

It is an interesting open problem whether the regularity of the solution in Proposition 7 is the minimal hypothesis. Also, it would be interesting to show a similar result to Proposition 7 for \(g \not = 0\).

5 Optimal Control Problem and Convergence Result

The purpose of this section is twofold. First, we analyse the class of optimal control problems for a system governed by the variational–hemivariational inequality in Problem 3. Second, we study the dependence of the solution on the yield limit g.

Let \({{\mathscr {U}}} = H = L^2(\Omega ; {\mathbb {R}}^d)\) be the space of controls representing external body forces. For every \({\varvec{f}}\in {{\mathscr {U}}}\), let \(S({\varvec{f}}) \subset V\) be the solution set to Problem 3 corresponding to \({\varvec{f}}\). From the previous section, we know that under the hypotheses of Theorem 4 or of Theorem 6, the set \(S({\varvec{f}})\) is nonempty. Consider the following optimal control problem. Given a nonempty subset \({\mathscr {U}}_{ad}\) of \({\mathscr {U}}\) denoting the set of admissible controls, and an objective functional \(F :{\mathscr {U}} \times V \rightarrow {\mathbb {R}}\), \(F=F({\varvec{f}}, {\varvec{u}})\), find \({\varvec{f}}^* \in {\mathscr {U}}_{ad}\) and a state \({\varvec{u}}^* \in S({\varvec{f}}^*)\) such that

$$\begin{aligned} F({\varvec{f}}^*, {\varvec{u}}^*) = \inf \{ \, F({\varvec{f}}, {\varvec{u}}) \mid {\varvec{f}}\in {{\mathscr {U}}}_{ad}, \, {\varvec{u}}\in S({\varvec{f}}) \, \}. \end{aligned}$$
(52)

A pair which solves (52) is called an optimal solution. We will show the existence of optimal solutions to (52) using a result on the continuous dependence of solutions on the right hand side.

Theorem 8

Let the hypotheses of Theorem 4 or of Theorem 6 hold, \({\varvec{f}}_n\), \({\varvec{f}}\in {{\mathscr {U}}}\) and \({\varvec{f}}_n \rightharpoonup {\varvec{f}}\) in \({{\mathscr {U}}}\). Then, for every sequence \(\{{\varvec{u}}_n \}\) of solutions to Problem 3 corresponding to \({\varvec{f}}_n\), one can find a subsequence \(\{{\varvec{u}}_{n_k} \} \subset \{ {\varvec{u}}_n \}\) such that \({\varvec{u}}_{n_k} \rightharpoonup {\varvec{u}}\) in V and \({\varvec{u}}\) is a solution to Problem 3 corresponding to \({\varvec{f}}\).

Proof

Let \({\varvec{f}}_n\), \({\varvec{f}}\in {{\mathscr {U}}}\) and \({\varvec{f}}_n \rightharpoonup {\varvec{f}}\) in \({{\mathscr {U}}}\). Suppose that the hypotheses of Theorem 4 hold. The proof under the hypotheses of Theorem 6 is analogous. We infer that there exists \({\varvec{u}}_n \in V\) such that (see (31))

$$\begin{aligned}&\langle (A + B){\varvec{u}}_n, \varvec{{v}}- {\varvec{u}}_n \rangle +\int _{\Gamma _1} k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; \varvec{{v}}_\tau -{\varvec{u}}_{n\tau }) \, d\Gamma \nonumber \\&\quad + \, \varphi (\varvec{{v}}) - \varphi ({\varvec{u}}_n) \ge \int _{\Omega } {\varvec{f}}_n \cdot (\varvec{{v}}-{\varvec{u}}_n)\, dx \end{aligned}$$
(53)

for all \(\varvec{{v}}\in V\), where A, B and \(\varphi \) are defined by (26), (27) and (29), respectively. First we show the a priori bound for the solutions. To this end, we take \(\varvec{{v}}={\varvec{0}}\) as a test function in (53) to obtain

$$\begin{aligned} \langle (A + B) {\varvec{u}}_n, {\varvec{u}}_n \rangle \le \int _{\Gamma _1} k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau };-{\varvec{u}}_{n\tau }) \, d\Gamma + \varphi ({\varvec{0}}) - \varphi ({\varvec{u}}_n) + \int _{\Omega } {\varvec{f}}_n \cdot {\varvec{u}}_n \, dx . \end{aligned}$$

From (41) and (42), we know that \(\langle (A + B) {\varvec{u}}_n, {\varvec{u}}_n \rangle \ge \mu _0\, \Vert {\varvec{u}}_n\Vert _V^2\). By H(j)(iii), H(k)(iii), the continuity of \(\gamma \), see (9), and [26, Proposition 3.23(iii)], we have

$$\begin{aligned}&\int _{\Gamma _1} k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau };-{\varvec{u}}_{n\tau }) \, d\Gamma \le k_1 \int _{\Gamma _1} \Vert \partial j ({\varvec{u}}_{n\tau }) \Vert \Vert {\varvec{u}}_{n\tau }\Vert \, d\Gamma \\&\quad \le k_1 \int _{\Gamma _1} (c_0({\varvec{x}}) + c_1 \Vert {\varvec{u}}_{n\tau }\Vert ) \Vert {\varvec{u}}_{n\tau }\Vert \, d\Gamma \le k_1 \Vert c_0\Vert _{L^2(\Gamma _1)} \Vert \gamma \Vert \Vert {\varvec{u}}_n\Vert _V + k_1 c_1\Vert \gamma \Vert ^2\Vert {\varvec{u}}_n\Vert ^2_V. \end{aligned}$$

We also get \(\varphi ({\varvec{0}}) = 0\) and \(\varphi ({\varvec{u}}_n) \ge 0\). By the properties stated above, it follows that

$$\begin{aligned} (\mu _0 - k_1 c_1\Vert \gamma \Vert ^2)\Vert {\varvec{u}}_n\Vert _V^2 \le (k_1 \Vert c_0\Vert _{L^2(\Gamma _1)} \Vert \gamma \Vert + \Vert {\varvec{f}}_n \Vert _H) \, \Vert {\varvec{u}}_n\Vert _V, \end{aligned}$$

and next by the smallness condition (25), we deduce

$$\begin{aligned} \Vert {\varvec{u}}_n\Vert _V \le M \ \ \text{ with } \ \ M := \frac{k_1 \Vert c_0\Vert _{L^2(\Gamma _1)}\Vert \gamma \Vert + \Vert {\varvec{f}}_n\Vert _H}{\mu _0-k_1c_1\Vert \gamma \Vert ^2}. \end{aligned}$$
(54)

From (54) we deduce that \(\{{\varvec{u}}_n\}\) is bounded in V uniformly with respect to n. Hence, by passing to a subsequence if necessary, we may assume that there exists \({\varvec{u}}\in V\) such that

$$\begin{aligned} {\varvec{u}}_n \rightharpoonup {\varvec{u}}\ \ \text{ in } \ \ V. \end{aligned}$$
(55)

It remains to show that \({\varvec{u}}\in S({\varvec{f}})\). From the compactness of the trace operator \(\gamma \), see (9), and [26, Theorem 2.39], we may suppose that at least for a next subsequence, we have \({\varvec{u}}_n \ \rightarrow {\varvec{u}}\) in \(L^2(\Gamma _1;{\mathbb {R}}^d)\) and

$$\begin{aligned} {\varvec{u}}_n({\varvec{x}}) \rightarrow {\varvec{u}}({\varvec{x}}) \ \ \text{ a.e. } \ \ {\varvec{x}}\in \Gamma _1, \ \ \Vert {\varvec{u}}_n ({\varvec{x}})\Vert _{{\mathbb {R}}^d} \le \rho ({\varvec{x}}) \ \ \text{ a.e. } \ {\varvec{x}}\in \Gamma _1 \end{aligned}$$
(56)

with \(\rho \in L^2(\Gamma _1)\).

Let \(\varvec{{v}}={\varvec{u}}\) be a test function in (53). We have

$$\begin{aligned}&\langle (A +B){\varvec{u}}_n, {\varvec{u}}_n - {\varvec{u}}\rangle \le \int _{\Gamma _1} k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \nonumber \\&\quad + \varphi ({\varvec{u}}) - \varphi ({\varvec{u}}_n) + \int _{\Omega } {\varvec{f}}_n \cdot ({\varvec{u}}_n-{\varvec{u}}) \, dx. \end{aligned}$$
(57)

We use [26, Proposition 3.23(ii)] and convergence (56) to get

$$\begin{aligned} \limsup j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{n\tau }-{\varvec{u}}_{\tau }) \le j^0({\varvec{u}}_{n\tau }; {\varvec{0}}) = 0 \ \ \text{ for } \text{ a.e. } \ \ {\varvec{x}}\in \Gamma _{1}. \end{aligned}$$

Combining the latter and the estimate

$$\begin{aligned} k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau } -{\varvec{u}}_{n\tau }) \le k_1 (c_0({\varvec{x}}) + c_1 \rho ({\varvec{x}})) (\Vert {\varvec{u}}({\varvec{x}}) \Vert + \rho ({\varvec{x}})) =: q({\varvec{x}}) \end{aligned}$$

with \(q \in L^1(\Gamma _{1})\), by Fatou’s lemma, we get

$$\begin{aligned}&\limsup \int _{\Gamma _1} k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \le \int _{\Gamma _1} \limsup k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \nonumber \\&\quad \le \int _{\Gamma _1} \limsup (k({\varvec{u}}_{n\tau })-k({\varvec{u}}_\tau )) q \, d\Gamma + \int _{\Gamma _1} k({\varvec{u}}_{\tau }) \limsup j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \le 0. \end{aligned}$$
(58)

We use the compactness of the embedding \(V \subset H\) and (55) to get \({\varvec{u}}_n \rightarrow {\varvec{u}}\) in H. Note that the functional \(\varphi \) is weakly lsc on V, being convex and lsc on V, which entails \(\limsup \, (\varphi ({\varvec{u}}) - \varphi ({\varvec{u}}_n)) \le 0\). The latter together with (58) allows to deduce from (57) that \(\limsup \langle (A+B) {\varvec{u}}_n, {\varvec{u}}_n - {\varvec{u}}\rangle \le 0\). Recalling that the operator \(A+B :V \rightarrow V^*\) is pseudomonotone, by a result in Sect. 2, we have

$$\begin{aligned} (A+B) {\varvec{u}}_n \rightharpoonup (A+B) {\varvec{u}}\ \ \text{ in } \ V^* \ \ \text{ and } \ \ \langle (A+B) {\varvec{u}}_n, {\varvec{u}}_n -{\varvec{u}}\rangle \rightarrow 0. \end{aligned}$$
(59)

Next, let \(\varvec{{v}}\in V\) be arbitrary. We will pass to the limit in the inequality (53). From (59), it follows

$$\begin{aligned} \lim \langle (A +B){\varvec{u}}_n, \varvec{{v}}- {\varvec{u}}_n \rangle = \langle (A +B){\varvec{u}}, \varvec{{v}}- {\varvec{u}}\rangle . \end{aligned}$$
(60)

Analogously to (58), we obtain

$$\begin{aligned} \limsup \int _{\Gamma _1} k({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; \varvec{{v}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \le \int _{\Gamma _1} k({\varvec{u}}_{\tau }) j^0({\varvec{u}}_{\tau }; \varvec{{v}}_{\tau }-{\varvec{u}}_{\tau }) \, d\Gamma , \end{aligned}$$
(61)

and

$$\begin{aligned} \limsup \, (\varphi (\varvec{{v}}) - \varphi ({\varvec{u}}_n)) \le \varphi (\varvec{{v}}) - \varphi ({\varvec{u}}). \end{aligned}$$
(62)

Exploiting (60), (61) and (62) from (53) in the limit as \(n \rightarrow \infty \), we obtain that \({\varvec{u}}\in V\) solves Problem 3. Hence \({\varvec{u}}\in S({\varvec{f}})\) which completes the proof. \(\square \)

Next, we consider a distributed parameter control problem which can be realized through electromagnetic forces, see [19] and the references therein. We need the following hypotheses.

\(\underline{H({{\mathscr {U}}}_{ad})}:\)    \({{\mathscr {U}}}_{ad}\) is a bounded and weakly closed subset of \({{\mathscr {U}}}\).

\(\underline{H(F)}:\)    F is lower semicontinuous on \({{\mathscr {U}}} \times V\) endowed with the weak topology.

Theorem 9

Under the hypotheses of Theorems 4 or 6, \(H({{\mathscr {U}}}_{ad})\) and H(F), the problem (52) has an optimal solution.

Proof

Let \(\{ ({\varvec{f}}_n, {\varvec{u}}_n)\}\) be a minimizing sequence for problem (52). This means that \({\varvec{f}}_n \in {{\mathscr {U}}}_{ad}\), \({\varvec{u}}_n \in S({\varvec{f}}_n)\) and

$$\begin{aligned} \lim F({\varvec{f}}_n, {\varvec{u}}_n) = \inf \{ \, F({\varvec{f}}, {\varvec{u}}) \mid {\varvec{f}}\in {{\mathscr {U}}}_{ad}, \, {\varvec{u}}\in S({\varvec{f}}) \, \} = : q. \end{aligned}$$

From \(H({{\mathscr {U}}}_{ad})\) it is clear that the sequence \(\{ {\varvec{f}}_n \}\) stays in a bounded subset of V. We use the reflexivity of V and we may assume, by passing to a subsequence if necessary, that \({\varvec{f}}_n \rightharpoonup {\varvec{f}}^*\) in \({{\mathscr {U}}}\). The weak closedness of the set \({{\mathscr {U}}}_{ad}\) entails \({\varvec{f}}^* \in {{\mathscr {U}}}_{ad}\). By Theorem 8 we infer, by passing to a subsequence, that \({\varvec{u}}\rightharpoonup {\varvec{u}}^*\) in V, where \({\varvec{u}}^* \in S({\varvec{f}}^*)\). By the weak lower semicontinuity of F, we have \(q \le F({\varvec{f}}^*, {\varvec{u}}^*) \le \liminf F({\varvec{f}}_n, {\varvec{u}}_n) = q\). This completes the proof. \(\square \)

Applying [10, Theorem 4.6.7] and the continuity of the map \(V \ni {\varvec{u}}\mapsto {{\mathbb {D}}}{\varvec{u}}\in L^2(\Omega ;{\mathbb {M}}^d)\), we deduce that if \(L:\Omega \times {\mathbb {R}}^d \times {{\mathbb {M}}}^d \times {\mathbb {R}}^d \rightarrow {\mathbb {R}}\) is a measurable function such that

  1. (a)

    \(L({\varvec{x}}, \cdot , \cdot , \cdot )\) is lower semicontinuous for a.e. \({\varvec{x}}\in \Omega \),

  2. (b)

    \(L({\varvec{x}}, {\varvec{u}}, \cdot , \cdot )\) is convex for all \({\varvec{u}}\in {\mathbb {R}}^d\), a.e. \({\varvec{x}}\in \Omega \),

  3. (c)

    \(L({\varvec{x}},{\varvec{u}}, {\varvec{\xi }},{\varvec{f}}) \ge \alpha ({\varvec{x}}) - c \, (\Vert {\varvec{u}}\Vert + \Vert {\varvec{\xi }}\Vert + \Vert {\varvec{f}}\Vert )\) for all \(({\varvec{u}},{\varvec{\xi }}, {\varvec{f}}) \in {\mathbb {R}}^d \times {{\mathbb {M}}}^d \times {\mathbb {R}}^d\),

      a.e. \({\varvec{x}}\in \Omega \) with \(c > 0\) and \(\alpha \in L^1(\Omega )\),

then the integral cost

$$\begin{aligned} F({\varvec{f}}, {\varvec{u}}) = \int _{\Omega } L(x, {\varvec{u}}(x), {{\mathbb {D}}} {\varvec{u}}(x), {\varvec{f}}(x)) \, dx \ \ \text{ for } \ \ ({\varvec{f}}, {\varvec{u}}) \in {{\mathscr {U}}} \times V \end{aligned}$$

satisfies H(F). A typical example of the cost is as follows

$$\begin{aligned} F({\varvec{f}}, {\varvec{u}}) = \int _{\Omega } \big ( w_1 \, \Vert {\varvec{u}}- {\varvec{u}}_d \Vert ^2 + w_2 \, \Vert \mathrm{rot} {\varvec{u}}\Vert ^2 + + w_3 \, \Vert {{\mathbb {D}}} {\varvec{u}}\Vert ^2 + w_4 \, \Vert {\varvec{f}}- {\varvec{f}}_0 \Vert ^2 \big )\, dx, \end{aligned}$$

where \(w_i \in L^\infty (\Omega )\), \(i=1\), \(\ldots \), 4 represent nonnegative weights. This objective functional is a compromise to minimize: the quadratic velocity tracking term, where \({\varvec{u}}_d\) is a given desired velocity profile, the vorticity of the optimal fluid flow responsible for turbulence, the rate dissipation energy (drug), and the control force to be as close as possible to a given element \({\varvec{f}}_0\). Here \(\mathrm{rot} {\varvec{u}}\) denotes the vorticity of the vector function

$$\begin{aligned} \mathrm{rot} {\varvec{u}}= {\left\{ \begin{array}{ll} \displaystyle \left( \frac{\partial u_3}{\partial x_2}- \frac{\partial u_1}{\partial x_3}, \frac{\partial u_1}{\partial x_3}- \frac{\partial u_3}{\partial x_1}, \frac{\partial u_2}{\partial x_1}- \frac{\partial u_1}{\partial x_2} \right) &{}\mathrm{if} \, d=3, \\ \displaystyle \frac{\partial u_2}{\partial x_1}- \frac{\partial u_1}{\partial x_2} &{}\mathrm{if} \, d =2. \end{array}\right. } \end{aligned}$$

Finally, we deal with the upper semicontinuity property of the solution set to Problem 3 with respect to the functions \((k, {\varvec{f}})\) and the yield limit g of the fluid which is assumed to be a nonnegative constant. Consider a sequence of inequalities in Problem 3, that is, find a velocity \({\varvec{u}}\in V\) such that

$$\begin{aligned} P(f_n, g_n, k_n)\qquad \left\{ \begin{array}{lll} \displaystyle \int _{\Omega } \mu (\Vert {{\mathbb {D}}}{\varvec{u}}\Vert ){\mathbb {D}}{\varvec{u}}: {{\mathbb {D}}} (\varvec{{v}}-{\varvec{u}}) \, dx - \int _\Omega ({\varvec{u}}\otimes {\varvec{u}}) : {\mathbb {D}}(\varvec{{v}}-{\varvec{u}}) \, dx\\ \displaystyle \quad + \, g_n \int _\Omega (\Vert {\mathbb {D}}\varvec{{v}}\Vert - \Vert {\mathbb {D}} {\varvec{u}}\Vert )\, dx +\int _{\Gamma _1} k_n({\varvec{u}}_\tau ) j^0({\varvec{u}}_\tau ; \varvec{{v}}_\tau -{\varvec{u}}_\tau ) \, d\Gamma \\ \displaystyle \qquad \ge \int _\Omega {\varvec{f}}_n \cdot (\varvec{{v}}-{\varvec{u}}) \, dx \ \ \text{ for } \text{ all } \ \ \varvec{{v}}\in V. \end{array} \right. \end{aligned}$$

Theorem 10

Let hypotheses \(H(\mu )\), H(j) and (25) hold, and k, \(k_n\) satisfy H(k) uniformly with respect to n. We assume

$$\begin{aligned}&g_n \ge 0, \ g_n \rightarrow g, \\&\quad {\varvec{f}}, {\varvec{f}}_n \in H, \ {\varvec{f}}_n \rightharpoonup {\varvec{f}}\ \ \text{ in } \ \ H, \\&\quad k_n({\varvec{x}}, {\varvec{\xi }}_n) \rightarrow k({\varvec{x}}, {\varvec{\xi }}) \ \ \text{ for } \text{ all } \ \ {\varvec{\xi }}_n \rightarrow {\varvec{\xi }}, \ \text{ a.e. } \ {\varvec{x}}\in \Gamma _1, \end{aligned}$$

as \(n \rightarrow \infty \). Let \(\{ {\varvec{u}}_n \} \subset V\) be a sequence of solutions to problem \(P(f_n, g_n, k_n)\). Then, there is a subsequence of \(\{ {\varvec{u}}_n \}\) converging weakly in V to an element \({\varvec{u}}\), where \({\varvec{u}}\in V\) is a solution to problem P(fgk).

Proof

Under the notation (26), (27) and (30), problem \(P(f_n, g_n, k_n)\) for the sequence \(\{{\varvec{u}}_n\}\) takes the form

$$\begin{aligned}&\langle (A +B){\varvec{u}}_n, \varvec{{v}}- {\varvec{u}}_n \rangle +\int _{\Gamma _1} k_n({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; \varvec{{v}}_\tau -{\varvec{u}}_{n\tau }) \, d\Gamma \nonumber \\&\quad + \, g_n (I(\varvec{{v}}) - I({\varvec{u}}_n)) \ge \int _\Omega {\varvec{f}}_n \cdot (\varvec{{v}}-{\varvec{u}}_n) \, dx \end{aligned}$$
(63)

for all \(\varvec{{v}}\in V\), where \(I :V \rightarrow {\mathbb {R}}\) is defined by

$$\begin{aligned} I(\varvec{{v}}) = \int _\Omega \Vert {\mathbb {D}} \varvec{{v}}\Vert \, dx, \ \ \varvec{{v}}\in V. \end{aligned}$$

For the sequence \(\{ {\varvec{u}}_n \} \subset V\), we obtain an estimate analogous to (54) with \(M > 0\) independent of n, since the sequence \(\{ {\varvec{f}}_n \}\) is uniformly norm bounded.

Then, we may assume that there exists \({\varvec{u}}\in V\) such that, at least for a subsequence, we have

$$\begin{aligned} {\varvec{u}}_n \rightharpoonup {\varvec{u}}\ \ \text{ in } \ \ V, \end{aligned}$$
(64)

and by the compactness of the trace operator, \({\varvec{u}}_n \ \rightarrow {\varvec{u}}\) in \(L^2(\Gamma _1;{\mathbb {R}}^d)\) and

$$\begin{aligned} {\varvec{u}}_n({\varvec{x}}) \rightarrow {\varvec{u}}({\varvec{x}}) \ \ \text{ a.e. } \ \ {\varvec{x}}\in \Gamma _1, \ \ \Vert {\varvec{u}}_n ({\varvec{x}})\Vert _{{\mathbb {R}}^d} \le \rho ({\varvec{x}}) \ \ \text{ a.e. } \ {\varvec{x}}\in \Gamma _1 \end{aligned}$$
(65)

with \(\rho \in L^2(\Gamma _1)\). It remains to show that \({\varvec{u}}\) is a solution to problem P(fgk). The passage to the limit, as \(n \rightarrow \infty \), in (63) can be done as in the proof of Theorem 8.

First, due the weak lower semicontinuity on V of the functional I, by (64), we have

$$\begin{aligned} I({\varvec{u}}) \le \liminf I({\varvec{u}}_n). \end{aligned}$$
(66)

Since \(g_n \ge 0\), \(I({\varvec{u}}_n) \ge 0\), by (66), it follows

$$\begin{aligned}&\limsup \, (g_n (I(\varvec{{v}}) - I({\varvec{u}}_n)) \le \limsup \, (g_n I(\varvec{{v}})) - \liminf (g_n I({\varvec{u}}_n)) \nonumber \\&\quad \le g I(\varvec{{v}}) - (\lim g_n)(\liminf I({\varvec{u}}_n)) \le g (I(\varvec{{v}}) - I({\varvec{u}})) \ \ \text{ for } \text{ all } \ \ \varvec{{v}}\in V, \end{aligned}$$
(67)

which, in particular, for \(\varvec{{v}}= {\varvec{u}}\) entails

$$\begin{aligned} \limsup \, (g_n (I({\varvec{u}}) - I({\varvec{u}}_n)) \le 0. \end{aligned}$$
(68)

Next, we proceed in two steps. First, we show that

$$\begin{aligned} (A+B) {\varvec{u}}_n \rightharpoonup (A+B) {\varvec{u}}\ \ \text{ in } \ V^* \ \ \text{ and } \ \ \langle (A+B) {\varvec{u}}_n, {\varvec{u}}_n -{\varvec{u}}\rangle \rightarrow 0. \end{aligned}$$
(69)

To this end, we choose \(\varvec{{v}}= {\varvec{u}}\) as a test function in problem (63) to get

$$\begin{aligned}&\langle (A +B){\varvec{u}}_n, {\varvec{u}}- {\varvec{u}}_n \rangle +\int _{\Gamma _1} k_n({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_\tau -{\varvec{u}}_{n\tau }) \, d\Gamma \nonumber \\&\quad + \, g_n (I({\varvec{u}}) - I({\varvec{u}}_n)) \ge \int _\Omega {\varvec{f}}_n \cdot (\varvec{{v}}-{\varvec{u}}_n) \, dx. \end{aligned}$$
(70)

By H(k)(iii), similarly as in (58), by (65) and Fatou’s lemma, we have

$$\begin{aligned}&\limsup \int _{\Gamma _1} k_n({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \le \int _{\Gamma _1} \limsup k_n({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \nonumber \\&\quad \le \int _{\Gamma _1} \limsup |k_n({\varvec{u}}_{n\tau })-k({\varvec{u}}_\tau )| q \, d\Gamma + \int _{\Gamma _1} k({\varvec{u}}_{\tau }) \limsup j^0({\varvec{u}}_{n\tau }; {\varvec{u}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \le 0, \nonumber \\ \end{aligned}$$
(71)

where \(q \in L^1(\Gamma _{1})\). We use (68), (71) and

$$\begin{aligned} \int _{\Omega } {\varvec{f}}_n \cdot {\varvec{u}}_n \, dx \rightarrow \int _{\Omega } {\varvec{f}}\cdot {\varvec{u}}\, dx, \end{aligned}$$
(72)

and then from (70) we deduce (69).

Subsequently, let \(\varvec{{v}}\in V\) be arbitrary. We shall pass to the limit in (63). To this end, we exploit (67), (69), (72) and the inequality

$$\begin{aligned} \limsup \int _{\Gamma _1} k_n({\varvec{u}}_{n\tau }) j^0({\varvec{u}}_{n\tau }; \varvec{{v}}_{\tau }-{\varvec{u}}_{n\tau }) \, d\Gamma \le \int _{\Gamma _1} k({\varvec{u}}_{\tau }) j^0({\varvec{u}}_{\tau }; \varvec{{v}}_{\tau }-{\varvec{u}}_{\tau }) \, d\Gamma . \end{aligned}$$
(73)

Finally, we conclude that \({\varvec{u}}\in V\) is a solution to problem P(fgk) which concludes the proof. \(\square \)

As a corollary from Theorem 10, we deduce that for\(g_n \rightarrow 0\) and \(\mu (r) = \mu _0\), \(k_n = k\), \({\varvec{f}}_n = {\varvec{f}}\), the Bingham fluid tends to behave as a Newtonian one.

6 Final Comments

The variational formulation in Problem 3 does not involve the fluid pressure due to a property of the test functions. For a non zero yield limit, it is an interesting and challenging problem to recover the pressure from this variational formulation. An attempt to retrieve the pressure for the non-Newtonian fluid with Dirichlet boundary conditions can be found in [36]. Related results on this topic involve [34] and a more recent paper [31]. The problems for the time-dependent Bingham type models via the variational–hemivariational inequalities have not been considered in literature and are worth studying in the future.

It would be also of interest to extend the results of this paper to the case of a more general law \({{\mathbb {S}}} = {{\mathbb {G}}}({{\mathbb {D}}})\) with a given function \({{\mathbb {G}}} :{{\mathbb {M}}}^d \rightarrow {{\mathbb {M}}}^d\) instead of \({{\mathbb {D}}} \mapsto \mu (\Vert {{\mathbb {D}}}\Vert ) {{\mathbb {D}}}\) in the Bingham constitutive relation (2). This will open a way to deal with various phenomena for fluids that the model (2) could not cope with, see, e.g., [6], and the references therein.