The Kelvin–Voigt–Brinkman–Forchheimer Equations with Non-Homogeneous Boundary Conditions

Abstract

We investigate the well-posedness of an initial boundary value problem for the Kelvin–Voigt–Brinkman–Forchheimer equations with memory and variable viscosity under a non-homogeneous Dirichlet boundary condition. A theorem about the global-in-time existence and uniqueness of a strong solution of this problem is proved under some smallness requirements on the size of the model data. For obtaining this result, we used a new technique, which is based on the operator treatment of the initial boundary value problem with the consequent application of an abstract theorem about the local unique solvability of an operator equation containing an isomorphism between Banach spaces with two kind perturbations: bounded linear and differentiable nonlinear having a zero Fréchet derivative at a zero element. Our work extends the existing frameworks of mathematical analysis and understanding of the dynamics of non-Newtonian fluids in porous media.

1. Introduction

Let $Q:=\mathsf{\Omega }\times S$ be a space-time cylinder, where $\mathsf{\Omega }$ is a bounded domain in space ${\mathbb{R}}^d$$(d=2,3)$ and S denotes a bounded time interval $(0,T)$ with fixed $T>0$. In this paper, we are concerned with the Kelvin–Voigt–Brinkman–Forchheimer equations, which describe unsteady flows of an incompressible viscoelastic fluid through a porous media:

(1)$\left\{\begin{array}{cc}{\partial }_t(\overrightarrow{u}-\alpha \mathsf{\Delta }\overrightarrow{u})+\delta \nabla ·(\overrightarrow{u}\otimes \overrightarrow{u})-\nabla ·\left(\eta \mathbb{D}\left(\overrightarrow{u}\right)\right)\hfill & \hfill \\ {\displaystyle -\underset{0}{\overset{t}{\int }}\xi (s,t)\mathsf{\Delta }\overrightarrow{u}(\overrightarrow{x},s)\mathrm{d}s+\lambda \overrightarrow{u}+\nabla \pi =\overrightarrow{f}}\hfill & \mathrm{i}\mathrm{n} Q,\hfill \\ \mathrm{div}\left(\overrightarrow{u}\right)=0\hfill & \mathrm{i}\mathrm{n} Q,\hfill \end{array}\right.$

where

• $\overrightarrow{u}$ is the velocity field, $\overrightarrow{u}:\overline{Q}\to {\mathbb{R}}^d$;

• $\pi$ is the pressure, $\pi :\overline{Q}\to \mathbb{R}$;

• $\overrightarrow{f}$ is the external force field; $\overrightarrow{f}:\overline{Q}\to {\mathbb{R}}^d$;

• $\delta$ is a dimensionless parameter included in the convective term; $\delta \in \{0,1\}$;

• $\alpha$ is a parameter characterizing the elasticity of the media, $\alpha \ge 0$ ;

• $\eta$ is the viscosity function, $\eta :\overline{Q}\to [0,\infty )$;

• $\xi$ is a function describing hereditary effects (memory), $\xi :\overline{S}\times \overline{S}\to \mathbb{R}$;

• $\lambda$ is a parameter characterizing the permeability of the media, $\lambda \ge 0$;

• the symbol ⊗ denotes the tensor product of vectors, that is, $\overrightarrow{x}\otimes \overrightarrow{y}:={\left(x_i{y}_j\right)}_{i,j=1}^d$ for any vectors $\overrightarrow{x},\overrightarrow{y}\in {\mathbb{R}}^d$;

• the symbol ∇ denotes for the gradient with respect to the space variables $x_1,\cdots ,x_d$, that is, $\nabla \pi :={\left({\displaystyle \frac{\partial \pi }{\partial x_1}},\cdots ,{\displaystyle \frac{\partial \pi }{\partial x_d}}\right)}^{\top }$;

• the differential operators div, $\mathsf{\Delta }$ and $\nabla ·$ are defined as follows:

  $\mathrm{div}\left(\overrightarrow{v}\right)=\sum _{i=1}^d{\displaystyle \frac{\partial v_i}{\partial x_i}},\mathsf{\Delta }\overrightarrow{v}:=\sum _{i=1}^d{\displaystyle \frac{{\partial }^2\overrightarrow{v}}{\partial x_i^2}},\nabla ·\mathbb{M}:={\left(\sum _{i=1}^d{\displaystyle \frac{\partial {\mathrm{M}}_{i1}}{\partial x_i}},\cdots ,\sum _{i=1}^d{\displaystyle \frac{\partial {\mathrm{M}}_{id}}{\partial x_i}}\right)}^{\top },$

  for a vector-valued function $\overrightarrow{v}:{\mathbb{R}}^d\to {\mathbb{R}}^d$ and a matrix-valued function $\mathbb{M}:{\mathbb{R}}^d\to {\mathbb{R}}^{d\times d}$.

As can be seen from Table 1, particular cases of system (1) appear in studying many important models for dynamics of incompressible fluids, including various models for viscoelastic media with memory. Moreover, the assumption that the viscosity is variable allows us to consider the important case of a mixture of two (or more) immiscible homogeneous viscous fluids.

In the present article, we will consider the so-called flow-through problem for the Kelvin–Voigt–Brinkman–Forchheimer equations, assuming that, for the flow velocity in system (1), the non-homogeneous Dirichlet boundary condition on the set $\partial \mathsf{\Omega }\times S$ is prescribed:

(2)$\overrightarrow{u}={\overrightarrow{u}}_b\mathrm{o}\mathrm{n} \partial \mathsf{\Omega }\times S.$

Moreover, we supplement this system with the initial condition

(3)$\overrightarrow{u}{|}_{t=0}={\overrightarrow{u}}_0\mathrm{i}\mathrm{n} \mathsf{\Omega }.$

In (2) and (3), ${\overrightarrow{u}}_b$ and ${\overrightarrow{u}}_0$ are given vector functions defined on $\partial \mathsf{\Omega }\times S$ and $\mathsf{\Omega }$, respectively.

Figure 1 shows an example of the flow configuration for a flow-through problem with

$\partial \mathsf{\Omega }={\mathrm{\Gamma }}_{\mathrm{in}}\cup {\mathrm{\Gamma }}_{\mathrm{out}}\cup {\mathrm{\Gamma }}_{\mathrm{lat}}^1\cup {\mathrm{\Gamma }}_{\mathrm{lat}}^2$

and

${\overrightarrow{u}}_b:=\left\{\begin{array}{cc}{\overrightarrow{u}}_{\mathrm{in}}\hfill & \mathrm{o}\mathrm{n} {\mathrm{\Gamma }}_{\mathrm{in}},\hfill \\ {\overrightarrow{u}}_{\mathrm{out}}\hfill & \mathrm{o}\mathrm{n} {\mathrm{\Gamma }}_{\mathrm{out}},\hfill \\ \overrightarrow{0}\hfill & \mathrm{o}\mathrm{n} {\mathrm{\Gamma }}_{\mathrm{lat}}^1\cup {\mathrm{\Gamma }}_{\mathrm{lat}}^2.\hfill \end{array}\right.$

Although various particular cases of the flow model (1) have been studied extensively by many researchers (see the books and the papers that are mentioned in Table 1 and the references in them), in the general case, the unique solvability of system (1)–(3) is a still open problem. One of the reasons is that non-zero Dirichlet boundary conditions produce serious difficulties in deriving a priori estimates for solutions and proving the well-posedness of corresponding boundary value problems for nonlinear governing equations, primarily in the case of three-dimensional flow problems through a domain with boundary having more than one connected component [36,37,38]. However, the proof of the unique solvability property is very important, since this is the first step in approbation of flow-through models, which are not merely of academic interest but has important consequences for engineering applications, for example, in modeling the time-dependent flows of fluids with complex rheology in tubes and pipeline networks.

It should be mentioned that the analysis of the well-posedness of model (1) with $\alpha >0$ and $\lambda =0$ under the no-slip boundary condition $\overrightarrow{u}{|}_{\partial \mathsf{\Omega }\times S}=\overrightarrow{0}$ was performed by Oskolkov [30]. More specifically, supposing that

$\begin{array}{c}\partial \mathsf{\Omega }\in C^{2,\mu },{\overrightarrow{u}}_0\in {\mathit{C}}^{2,\mu }(\overline{\mathsf{\Omega }})\cap {\mathit{H}}_0^1(\mathsf{\Omega }),\overrightarrow{f}\in {\mathit{L}}^{\infty }\left(S,{\mathit{C}}^{\mu }(\overline{\mathsf{\Omega }})\right)\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h} \mu \in (0,1),\\ {\partial }_t\overrightarrow{f}\in {\mathit{L}}^2\left(Q\right),\xi (t,s)\equiv k(t-s),\end{array}$

where k is a $C^1$-smooth function on $\overline{S}$, he has proved that the corresponding initial boundary value problem (IBVP) admits a unique solution $(\overrightarrow{u},\pi )$ satisfying the following inclusions:

$\overrightarrow{u}\in {\mathit{W}}^{1,\infty }(S,{\mathit{C}}^{2,\mu }(\overline{\mathsf{\Omega }})\cap {\mathit{H}}_0^1(\mathsf{\Omega })),{\displaystyle \frac{\partial \pi }{\partial x_i}}\in L^{\infty }(S,C^{\mu }(\overline{\mathsf{\Omega }})),\forall i=1,\cdots ,d.$

Di Plinio et al. [39] showed that the Kelvin–Voigt–Brinkman–Forchheimer system, where instantaneous viscosity is completely replaced by a memory term, is dissipative (in the sense of dynamical systems) and even admits exponential and global attractors of finite fractal dimension. Such properties of asymptotic well-posedness in the absence of instantaneous viscosity are rarely observed in the field of dynamical systems arising from fluid models. On the other hand, Yushkov [40] has discovered the blow-up effect for solutions of IBVP (1)–(3) with $\alpha =1$, $\delta =1$, $\eta \equiv 2$, $\xi (s,t)\equiv exp(-(t-s\left)\right)$, $\lambda =0$, $\overrightarrow{f}\equiv \overrightarrow{0}$ and ${\overrightarrow{u}}_b\equiv \overrightarrow{0}$ in the presence of the cubic source $-\overrightarrow{u}{\left|\overrightarrow{u}\right|}^2$. He has also found upper and lower bounds for the time of blowing up a solution.

This article is a continuation of the work [31], in which the well-posedness of an IBVP for flow model (1) with $\eta \equiv \mathrm{const}$ and $\lambda =0$ has been established for the case of small data. Our main aim is to investigate the existence and uniqueness of a regular-in-time strong solution to problem (1)–(3) under the assumption that the boundary velocity field ${\overrightarrow{u}}_b$ belongs to suitable fractional Sobolev space, for any $t\in \overline{S}$. We will construct solutions in a function class where the uniqueness and regularity properties are assured even for the 3D case.

The outline of our work is as follows. In Section 2, for the convenience of the reader, we give necessary preliminaries, including some valuable results from linear functional analysis (see Theorems 1 and 2) and one abstract theorem about the local well-posedness of a nonlinear operator equation involving an isomorphism between Banach spaces with Fréchet differentiable perturbations (see Theorem 3). Section 3 introduces the notion of a regular-in-time strong solution to the Kelvin–Voigt–Brinkman–Forchheimer system (see Definition 3). In Section 4, we state the main results of the article—the theorem about the unique solvability of problem (1)–(3) in the strong formulation, under appropriate smallness conditions for the model data $\xi$, $\overrightarrow{f}$, ${\overrightarrow{u}}_b$, ${\overrightarrow{u}}_0$ and the viscosity gradient $\nabla \eta$ (see Theorem 4). In Section 5, we formulate and prove some auxiliary propositions about the existence of a suitable divergence-free lifting operator and properties of two linear operators associated with the Kelvin–Voigt–Brinkman–Forchheimer equations (see Propositions 1–3). Next, using these propositions and Theorem 3, in Section 6, we prove the main results of this paper. Finally, Section 7 provides our concluding remarks.

2. Preliminaries

2.1. Isomorphisms and Some Related Results

We recall the notion of an isomorphism between Banach spaces.

For Banach spaces $\mathit{E}$ and $\mathit{F}$, by $\mathcal{L}(\mathit{E},\mathit{F})$ denote the collection of all continuous linear mappings from $\mathit{E}$ into $\mathit{F}$. As is known, $\mathcal{L}(\mathit{E},\mathit{F})$ is a Banach space with the operator norm ${\parallel ·\parallel }_{\mathcal{L}(\mathit{E},\mathit{F})}$ defined as follows:

${\parallel \mathit{L}\parallel }_{\mathcal{L}(\mathit{E},\mathit{F})}≔inf\left\{c\ge {0:\parallel \mathit{L}\mathit{v}\parallel }_{\mathit{F}}{\le c\parallel \mathit{v}\parallel }_{\mathit{E}}\mathrm{f}\mathrm{o}\mathrm{r} \mathrm{a}\mathrm{l}\mathrm{l} \mathit{v}\in \mathit{E}\right\},$

for any operator $\mathit{L}\in \mathcal{L}(\mathit{E},\mathit{F})$.

Definition 1.

Banach spaces $\mathit{E}$ and $\mathit{F}$ are said to be isomorphic if there exists a linear one-to-one (bijective) mapping $\mathit{A}:\mathit{E}\to \mathit{F}$ such that

(4) $\mathit{A}\in \mathcal{L}(\mathit{E},\mathit{F}),{\mathit{A}}^{-1}\in \mathcal{L}(\mathit{F},\mathit{E}).$

Definition 2.

A one-to-one operator $\mathit{A}$ satisfying both inclusions in (4) is called an isomorphism.

By ${\mathcal{L}}_{\mathrm{Isom}}(\mathit{E},\mathit{F})$ we denote the set of all isomorphisms from $\mathit{E}$ onto $\mathit{F}$.

If $\mathit{E}=\mathit{F}$, then, for brevity, we will write $\mathcal{L}\left(\mathit{E}\right)$ and ${\mathcal{L}}_{\mathrm{Isom}}\left(\mathit{E}\right)$ instead of $\mathcal{L}(\mathit{E},\mathit{E})$ and ${\mathcal{L}}_{\mathrm{Isom}}(\mathit{E},\mathit{E})$, respectively.

An important property of the isomorphisms set ${\mathcal{L}}_{\mathrm{Isom}}(\mathit{E},\mathit{F})$ is that this set is open in the space $\mathcal{L}(\mathit{E},\mathit{F})$. Namely, the following theorem holds.

Theorem 1

(see [41], Section 23). Let $\mathit{E}$ and $\mathit{F}$ be Banach spaces and $\mathit{A}\in {\mathcal{L}}_{\mathrm{Isom}}(\mathit{E},\mathit{F})$. Then, for an arbitrary operator $\mathit{B}\in \mathcal{L}(\mathit{E},\mathit{F})$ such that

${\parallel \mathit{B}\parallel }_{\mathcal{L}(\mathit{E},\mathit{F})}<{\displaystyle \frac{1}{\parallel {\mathit{A}}^{-1}{\parallel }_{\mathcal{L}(\mathit{F},\mathit{E})}}},$

the sum of the operators $\mathit{A}$ and $\mathit{B}$ is an isomorphism, that is,

$(\mathit{A}+\mathit{B})\in {\mathcal{L}}_{\mathrm{Isom}}(\mathit{E},\mathit{F}).$

Let us also give the formulation of the bounded inverse theorem (also called Banach isomorphism theorem).

Theorem 2

(see [42], Section 8.2). Let $\mathit{E}$ and $\mathit{F}$ be Banach spaces. Suppose that $\mathit{M}\in \mathcal{L}(\mathit{E},\mathit{F})$ and $\mathit{M}$ is a one-to-one mapping. Then, the inverse operator ${\mathit{M}}^{-1}$ is bounded, and hence the inclusion $\mathit{M}\in {\mathcal{L}}_{\mathrm{Isom}}(\mathit{E},\mathit{F})$ holds.

2.2. Local Solvability of Equations with Fréchet Differentiable Operators

Using the ideas from the work [31], we here present a result on the local unique solvability of a class of abstract nonlinear equations in Banach spaces.

Theorem 3

(Abstract theorem on local solvability). Suppose $\mathit{X}$ and $\mathit{Y}$ are isomorphic Banach spaces over the field $\mathbb{R}$ and ${\mathit{T}}_i:\mathit{X}\to \mathit{Y}$, where $i\in \{0,1,2\}$, are given operators. Moreover, let the five following conditions be satisfied:

• ($\mathcal{H}.$1)

  The inclusions ${\mathit{T}}_0\in {\mathcal{L}}_{\mathrm{Isom}}(\mathit{X},\mathit{Y})$ and ${\mathit{T}}_1\in \mathcal{L}(\mathit{X},\mathit{Y})$ hold.

• ($\mathcal{H}.$2)

  The operator norms $\parallel {\mathit{T}}_0^{-1}{\parallel }_{\mathcal{L}(\mathit{Y},\mathit{X})}$ and $\parallel {\mathit{T}}_1{\parallel }_{\mathcal{L}(\mathit{X},\mathit{Y})}$ obey the following relation

  (5)$\parallel {\mathit{T}}_0^{-1}{\parallel }_{\mathcal{L}(\mathit{Y},\mathit{X})}{\parallel {\mathit{T}}_1\parallel }_{\mathcal{L}(\mathit{X},\mathit{Y})}<1.$

• ($\mathcal{H}.$3)

  The mapping ${\mathit{T}}_2:\mathit{X}\to \mathit{Y}$ is continuously Fréchet differentiable.

• ($\mathcal{H}.$4)

  The equality ${\mathit{T}}_2\left(\mathbf{0}\right)=\mathbf{0}$ holds.

• ($\mathcal{H}.$5)

  The Fréchet derivative $\mathit{D}{\mathit{T}}_2\left(\mathbf{0}\right)$ is equal to the zero operator.

Then, there exist $ϵ>0$ and (open) neighborhood $\mathcal{U}$ of the zero element $\mathbf{0}$ in the Banach space $\mathit{X}$ such that, for any given $\mathit{h}\in \mathit{Y}$, the equation

(6) $({\mathit{T}}_0+{\mathit{T}}_1+{\mathit{T}}_2)\left(\mathit{u}\right)=\mathit{h}$

has a unique solution $\mathit{u}={\mathit{u}}_{\mathit{h}}$ in the neighborhood $\mathcal{U}$ provided that $\mathit{h}\in {\mathfrak{B}}_{ϵ}\left(\mathbf{0}\right)$, where

${\mathfrak{B}}_{ϵ}\left(\mathbf{0}\right):={\{\mathit{y}\in \mathit{Y}:\parallel \mathit{y}\parallel }_{\mathit{Y}}<ϵ\}.$

Proof. 

First, we introduce the auxiliary operator $\mathit{G}:\mathit{X}\to \mathit{X}$ by the formula

(7)$\mathit{G}:={\mathit{T}}_0^{-1}\circ {\mathit{T}}_1+{\mathit{T}}_0^{-1}\circ {\mathit{T}}_2.$

Clearly, this operator is continuously Fréchet differentiable.

Next, by applying ${\mathit{T}}_0^{-1}$ to the right-hand and left-hand sides of Equation (6) and taking into account (7), we obtain

(8)$({\mathbf{I}}_{\mathit{X}}+\mathit{G})\left(\mathit{u}\right)={\mathit{T}}_0^{-1}\mathit{h}.$

Here, ${\mathbf{I}}_{\mathit{X}}$ stands for the identity operator acting in the Banach space $\mathit{X}$.

Since Equations (6) and (8) are equivalent, it is sufficient to prove the theorem for Equation (8).

Due to the operators ${\mathit{T}}_0$ and ${\mathit{T}}_1$ are linear and condition ($\mathcal{H}.$4) is valid, we have

$\mathit{G}\left(\mathbf{0}\right)=({\mathit{T}}_0^{-1}\circ {\mathit{T}}_1)\left(\mathbf{0}\right)+({\mathit{T}}_0^{-1}\circ {\mathit{T}}_2)\left(\mathbf{0}\right)=\mathbf{0},$

and hence

(9)$({\mathbf{I}}_{\mathit{X}}+\mathit{G})\left(\mathbf{0}\right)=\mathbf{0}.$

Furthermore, taking into account relation (5) and condition ($\mathcal{H}.$5), we obtain the following estimate of the operator norm ${\parallel \mathit{D}\mathit{G}\left(\mathbf{0}\right)\parallel }_{\mathcal{L}\left(\mathit{X}\right)}$:

$\begin{array}{cc}\hfill {\parallel \mathit{D}\mathit{G}\left(\mathbf{0}\right)\parallel }_{\mathcal{L}\left(\mathit{X}\right)}& =\parallel {\mathit{T}}_0^{-1}\circ {\mathit{T}}_1+{\mathit{T}}_0^{-1}\circ \mathit{D}{\mathit{T}}_2{\left(\mathbf{0}\right)\parallel }_{\mathcal{L}\left(\mathit{X}\right)}\hfill \\ \hfill & =\parallel {\mathit{T}}_0^{-1}\circ {\mathit{T}}_1{\parallel }_{\mathcal{L}\left(\mathit{X}\right)}\hfill \\ \hfill & \le \parallel {\mathit{T}}_0^{-1}{\parallel }_{\mathcal{L}(\mathit{Y},\mathit{X})}{\parallel {\mathit{T}}_1\parallel }_{\mathcal{L}(\mathit{X},\mathit{Y})}\hfill \\ \hfill & <1,\hfill \end{array}$

whence one can conclude that the operator ${\mathbf{I}}_{\mathit{X}}+\mathit{D}\mathit{G}\left(\mathbf{0}\right)$ is an isomorphism, that is,

(10)$\left({\mathbf{I}}_{\mathit{X}}+\mathit{D}\mathit{G}\left(\mathbf{0}\right)\right)\in {\mathcal{L}}_{\mathrm{Isom}}\left(\mathit{X}\right).$

In view of relations (9) and (10), we can apply the inverse function theorem (see, for example, the monograph [42], Theorem 10.4) to Equation (8). This yields that there exist open neighborhoods ${\mathcal{U}}_1$ and ${\mathcal{U}}_2$ of the element $\mathbf{0}$ in the Banach space $\mathit{X}$ such that the restriction of the mapping ${\mathbf{I}}_{\mathit{X}}+\mathit{G}$ to the subset ${\mathcal{U}}_1$ satisfies the following property:

(11)$({\mathbf{I}}_{\mathit{X}}+\mathit{G}){|}_{{\mathcal{U}}_1}:{\mathcal{U}}_1\to {\mathcal{U}}_2 \mathrm{i}\mathrm{s} \mathrm{a} \mathrm{o}\mathrm{n}\mathrm{e}-\mathrm{t}\mathrm{o}-\mathrm{o}\mathrm{n}\mathrm{e} \mathrm{m}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{i}\mathrm{n}\mathrm{g}.$

Let us fix a positive number $ϵ$ such that, for any $\mathit{h}\in {\mathfrak{B}}_{ϵ}\left(\mathbf{0}\right)$, the inclusion ${\mathit{T}}_0^{-1}\mathit{h}\in {\mathcal{U}}_2$ holds.

Now, note that, to prove this theorem, it suffices to set $\mathcal{U}={\mathcal{U}}_1$. Indeed, let $\mathit{h}\in {\mathfrak{B}}_{ϵ}\left(\mathbf{0}\right)$. Then, taking into account property (11), we can define the element ${\mathit{u}}_{\mathit{h}}$ as follows:

${\mathit{u}}_{\mathit{h}}:=[{({\mathbf{I}}_{\mathit{X}}+\mathit{G})}^{-1}\circ {\mathit{T}}_0^{-1}]\mathit{h}.$

It can be directly checked that this element is a unique solution of Equation (8) in the neighborhood $\mathcal{U}$ of $\mathbf{0}\in \mathit{X}$. Thus, the proof of Theorem 3 is complete. □

2.3. Notation for Scalar Product and Euclidean Norm in ${\mathbb{R}}^d$ and ${\mathbb{R}}^{d\times d}$

Let $\overrightarrow{a},\overrightarrow{b}\in {\mathbb{R}}^d$ and $\mathbb{A},\mathbb{B}\in {\mathbb{R}}^{d\times d}$. We will use the following notation for the scalar product and the Euclidean norm in ${\mathbb{R}}^d$ and ${\mathbb{R}}^{d\times d}$, respectively:

$\begin{array}{cc}{\displaystyle \overrightarrow{a}·\overrightarrow{b}:=\sum _{i=1}^d{a}_i{b}_i,}\hfill & {\displaystyle |\overrightarrow{a}|:={\left(\sum _{i=1}^d{a}_i^2\right)}^{\frac{1}{2}},}\hfill \\ {\displaystyle \mathbb{A}:\mathbb{B}:=\sum _{i,j=1}^d{\mathrm{A}}_{ij}{\mathrm{B}}_{ij},}\hfill & {\displaystyle \left|\mathbb{A}\right|:={\left(\sum _{i,j=1}^d{\mathrm{A}}_{ij}^2\right)}^{\frac{1}{2}}.}\hfill \end{array}$

2.4. Spaces of Time-Independent Functions

We will use the following functional spaces:

• $C_c^{\infty }\left({\mathbb{R}}^d\right):=\left\{v\in C^{\infty }\left({\mathbb{R}}^d\right):\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\left(v\right) \mathrm{i}\mathrm{s} \mathrm{c}\mathrm{o}\mathrm{m}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{t}\right\}$;

• $C_c^{\infty }\left(\overline{\mathsf{\Omega }}\right):=\left\{{u|}_{\overline{\mathsf{\Omega }}}:u\in C_c^{\infty }\left({\mathbb{R}}^d\right)\right\}$;

• the Lebesgue space $L^p\left(\mathsf{\Omega }\right)$, $p\ge 1$ and the Sobolev space $H^n\left(\mathsf{\Omega }\right):=W^{n,2}\left(\mathsf{\Omega }\right)$, $n\in \mathbb{N}$.

The corresponding classes of functions with values in ${\mathbb{R}}^d$ are designated by the same symbols, but the first letter is highlighted in bold. For example,

$\begin{array}{cc}\hfill {\mathit{L}}^p\left(\mathsf{\Omega }\right)& :=\underset{d\mathrm{times}}{\underbrace{L^p\left(\mathsf{\Omega }\right)\times \cdots \times L^p\left(\mathsf{\Omega }\right)}} ,p\ge 1,\hfill \\ \hfill {\mathit{H}}^n\left(\mathsf{\Omega }\right)& :=\underset{d\mathrm{times}}{\underbrace{H^n\left(\mathsf{\Omega }\right)\times \cdots \times H^n\left(\mathsf{\Omega }\right)}} ,n\in \mathbb{N},\hfill \\ \hfill {\mathit{H}}^{n-{\displaystyle \frac{1}{2}}}\left(\partial \mathsf{\Omega }\right)& :=\underset{d\mathrm{times}}{\underbrace{H^{n-{\displaystyle \frac{1}{2}}}\left(\partial \mathsf{\Omega }\right)\times \cdots \times H^{n-{\displaystyle \frac{1}{2}}}\left(\partial \mathsf{\Omega }\right)}} ,n\in \mathbb{N}.\hfill \end{array}$

For handling boundary traces of functions belonging to $H^n\left(\mathsf{\Omega }\right)$, the fractional Sobolev space $H^{n-\frac{1}{2}}\left(\partial \mathsf{\Omega }\right)$ is used; see the books [9,43] for details.

By ${\gamma }_{\partial \mathsf{\Omega }}$ we denote the trace operator (see, for example, [9], Chapter III). This mapping is a surjective continuous linear operator from ${\mathit{H}}^n\left(\mathsf{\Omega }\right)$ into ${\mathit{H}}^{n-\frac{1}{2}}\left(\partial \mathsf{\Omega }\right)$ such that ${\gamma }_{\partial \mathsf{\Omega }}\overrightarrow{v}=\overrightarrow{v}{|}_{\partial \mathsf{\Omega }},$ for any $\overrightarrow{v}\in {\mathit{C}}_c^{\infty }\left(\overline{\mathsf{\Omega }}\right)$.

Furthermore, we define the subspace of ${\mathit{H}}^n\left(\mathsf{\Omega }\right)$ consisting of solenoidal (divergence-free) vector-valued functions:

${\mathit{H}}_{\mathrm{div}}^n\left(\mathsf{\Omega }\right)≔\left\{\overrightarrow{v}\in {\mathit{H}}^n\left(\mathsf{\Omega }\right):\mathrm{div}\left(\overrightarrow{v}\right)=0in\mathsf{\Omega }\right\},n\in \mathbb{N},$

and the subspace of ${\mathit{H}}^{n-\frac{1}{2}}\left(\partial \mathsf{\Omega }\right)$ consisting of boundary traces that satisfy the zero flux condition on the surface $\partial \mathsf{\Omega }$:

${\dot{\mathit{H}}}^{n-\frac{1}{2}}\left(\partial \mathsf{\Omega }\right)≔\left\{\overrightarrow{\omega }\in {\mathit{H}}^{n-\frac{1}{2}}\left(\partial \mathsf{\Omega }\right):\underset{\partial \mathsf{\Omega }}{\int }\overrightarrow{\omega }·\overrightarrow{\nu }\mathrm{d}\sigma =0\right\},n\in \mathbb{N}.$

In the last formula, $\overrightarrow{\nu }$ is the outward-pointing unit normal on $\partial \mathsf{\Omega }$ and the term $\mathrm{d}\sigma$ indicates an element of surface area on $\partial \mathsf{\Omega }$.

Now, we introduce the three spaces for the setting of flow problems in regions with impermeable solid boundaries:

$\begin{array}{cc}\hfill {\mathfrak{D}}_{\mathrm{div}}\left(\mathsf{\Omega }\right)& :=\left\{\overrightarrow{v}\in {\mathit{C}}^{\infty }\left(\mathsf{\Omega }\right)\cap {\mathit{H}}_{\mathrm{div}}^1\left(\mathsf{\Omega }\right):\mathrm{supp}\left(\overrightarrow{v}\right)\subset \mathsf{\Omega }\right\},\hfill \\ \hfill {\mathit{H}}_{0,\mathrm{div}}^1\left(\mathsf{\Omega }\right)& :=\mathrm{t}\mathrm{h}\mathrm{e} \mathrm{c}\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e} \mathrm{o}\mathrm{f} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{s}\mathrm{e}\mathrm{t} {\mathfrak{D}}_{\mathrm{div}}\left(\mathsf{\Omega }\right)\mathrm{i}\mathrm{n} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{S}\mathrm{o}\mathrm{b}\mathrm{o}\mathrm{l}\mathrm{e}\mathrm{v} \mathrm{s}\mathrm{p}\mathrm{a}\mathrm{c}\mathrm{e} {\mathit{H}}^1\left(\mathsf{\Omega }\right),\hfill \\ \hfill {\mathit{H}}_{0,\mathrm{div}}^2\left(\mathsf{\Omega }\right)& :={\mathit{H}}^2\left(\mathsf{\Omega }\right)\cap {\mathit{H}}_{0,\mathrm{div}}^1\left(\mathsf{\Omega }\right).\hfill \end{array}$

It is clear that a vector function belonging to ${\mathit{H}}_{0,\mathrm{div}}^1\left(\mathsf{\Omega }\right)$ or ${\mathit{H}}_{0,\mathrm{div}}^2\left(\mathsf{\Omega }\right)$ vanishes on the surface $\partial \mathsf{\Omega }$. More specifically, we have ${\gamma }_{\partial \mathsf{\Omega }}\overrightarrow{v}=\overrightarrow{0}$ for any $\overrightarrow{v}\in {\mathit{H}}_{0,\mathrm{div}}^l\left(\mathsf{\Omega }\right)$, where $l=1$ or 2.

2.5. Spaces of Time-Dependent Functions

Let $\mathit{E}$ be a Banach space and $S≔(0,T)$ with $0<T<\infty$.

We will consider spaces of functions defined on $\overline{S}$ with values in some functional space, for example, in a Lebesgue space or a Sobolev space.

By $\mathit{C}(\overline{S},\mathit{E})$ denote the space of all continuous functions from $\overline{S}$ into $\mathit{E}$ with the norm

${\parallel \mathit{v}\parallel }_{\mathit{C}(\overline{S},\mathit{E})}≔\underset{t\in \overline{S}}{max}{\parallel \mathit{v}\left(t\right)\parallel }_{\mathit{E}}.$

We will also use the space of all continuously differentiable functions

${\mathit{C}}^1(\overline{S},\mathit{E})≔\left\{\mathit{v}:\overline{S}\to \mathit{E}:\mathit{v}\in \mathit{C}(\overline{S},\mathit{E})\mathrm{a}\mathrm{n}\mathrm{d} {\mathit{v}}^{\prime }\in \mathit{C}(\overline{S},\mathit{E})\right\}$

with the norm

${\parallel \mathit{v}\parallel }_{{\mathit{C}}^1(\overline{S},\mathit{E})}≔{\parallel \mathit{v}\parallel }_{\mathit{C}(\overline{S},\mathit{E})}+{\parallel {\mathit{v}}^{\prime }\parallel }_{\mathit{C}(\overline{S},\mathit{E})}.$

Here and in the following discussion, the prime symbol ′ denotes the derivative with respect to time t.

2.6. Helmholtz–Weyl Decomposition and Leray Projection

By $\mathit{H}\left(\mathsf{\Omega }\right)$ we denote the closure of the set ${\mathfrak{D}}_{\mathrm{div}}\left(\mathsf{\Omega }\right)$ in the Lebesgue space ${\mathit{L}}^2\left(\mathsf{\Omega }\right)$.

Let

$\mathit{G}\left(\mathsf{\Omega }\right)≔\nabla H^1\left(\mathsf{\Omega }\right)=\left\{\nabla h:h\in H^1\left(\mathsf{\Omega }\right)\right\}.$

We will use the Helmholtz–Weyl decomposition of the space ${\mathit{L}}^2\left(\mathsf{\Omega }\right)$ into the solenoidal and gradient parts:

(12)${\mathit{L}}^2\left(\mathsf{\Omega }\right)=\mathit{H}\left(\mathsf{\Omega }\right)\oplus \mathit{G}\left(\mathsf{\Omega }\right);$

for details, see [9], Chapter IV. In relation (12), the symbol ⊕ denotes the orthogonal sum.

By ${\mathbf{P}}_{\mathit{H}\left(\mathsf{\Omega }\right)}$ we denote the orthogonal projection from the space ${\mathit{L}}^2\left(\mathsf{\Omega }\right)$ into its subspace $\mathit{H}\left(\mathsf{\Omega }\right)$. This projection is called the Leray projection.

2.7. Quotient of Sobolev Space $H^1\left(\mathsf{\Omega }\right)$ by Constants

In order to describe the pressure $\pi$, which is included in (1) via the gradient term $\nabla \pi$, it is convenient to introduce the equivalence relation “∼” on the Sobolev space $H^1\left(\mathsf{\Omega }\right)$ by

$f\sim g⟺\nabla (f-g)=\overrightarrow{0}\mathrm{a}.\mathrm{e}. \mathrm{i}\mathrm{n} \mathsf{\Omega },$

for any functions $f,g\in H^1\left(\mathsf{\Omega }\right)$.

Let

$R\left(\mathsf{\Omega }\right)≔\left\{\phi \in H^1\left(\mathsf{\Omega }\right):\nabla \phi =\overrightarrow{0}\mathrm{a}.\mathrm{e}. \mathrm{i}\mathrm{n} \mathsf{\Omega }\right\}.$

It is easily shown that

$\psi \in R\left(\mathsf{\Omega }\right)⟺\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} \mathrm{e}\mathrm{x}\mathrm{i}\mathrm{s}\mathrm{t}\mathrm{s} \mathrm{a} \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{t} c \mathrm{s}\mathrm{u}\mathrm{c}\mathrm{h} \mathrm{t}\mathrm{h}\mathrm{a}\mathrm{t} \psi \left(\overrightarrow{x}\right)=c\mathrm{a}.\mathrm{e}. \mathrm{i}\mathrm{n} \mathsf{\Omega }.$

By $H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)$ we denote the quotient of $H^1\left(\mathsf{\Omega }\right)$ by $R\left(\mathsf{\Omega }\right)$.

For an arbitrary function $\mathit{v}\in H^1\left(\mathsf{\Omega }\right)$, we introduce the equivalence class $⟨\mathit{v}⟩$ as follows:

$⟨\mathit{v}⟩≔\left\{w\in H^1\left(\mathsf{\Omega }\right):w\sim \mathit{v}\right\}\in H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right).$

Let $⟨f⟩\in H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)$. We define the gradient $\nabla ⟨f⟩$ of the equivalence class $⟨f⟩$ by

$\nabla ⟨f⟩≔\nabla f.$

Finally, let us introduce the norm in the space $H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)$ by the following formula

$\parallel ⟨f⟩{\parallel }_{H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)} ≔\parallel \nabla ⟨f⟩{\parallel }_{L^2\left(\mathsf{\Omega }\right)}.$

Taking into account Proposition 1.2 from the book [7] (see Chapter I, § 1), one can establish that the norm ${\parallel ·\parallel }_{H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)}$ is well defined.

3. Strong Formulation of Problem

Let us assume that the model data $\overrightarrow{f}$, ${\overrightarrow{u}}_b$, ${\overrightarrow{u}}_0$, $\xi$, $\eta$ satisfy the six conditions:

(13)$\begin{array}{ccc}\overrightarrow{f}\in \mathit{C}\left(\overline{S},{\mathit{L}}^2\left(\mathsf{\Omega }\right)\right),& {\overrightarrow{u}}_b\in {\mathit{C}}^1\left(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right)\right),& {\overrightarrow{u}}_0\in {\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right),\\ \xi \in C(\overline{S}\times \overline{S}),& \eta \in C\left(\overline{Q}\right),& \nabla \eta \in \mathit{C}\left(\overline{Q}\right).\end{array}$

Before presenting our main results, we introduce the notion of a strong solution to the IBVP under consideration.

Definition 3

(Strong solution). A pair $\left(\overrightarrow{u},⟨\pi ⟩\right)$ is called a strong solution of IBVP (1)–(3) if the following three conditions hold:

• $\overrightarrow{u}\in {\mathit{C}}^1\left(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)\right);$

• $⟨\pi ⟩\in C\left(\overline{S},H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)\right);$

• the vector function $\overrightarrow{u}$ and the equivalence class $⟨\pi ⟩$ satisfy the system

  $\left\{\begin{array}{cc}{\overrightarrow{u}}^{\prime }-\alpha \mathsf{\Delta }{\overrightarrow{u}}^{\prime }+\delta \nabla ·(\overrightarrow{u}\otimes \overrightarrow{u})-{(\nabla \eta )}^{\top }\mathbb{D}\left(\overrightarrow{u}\right)-\frac{\eta }{2}\mathsf{\Delta }\overrightarrow{u}\hfill & \\ {\displaystyle -\underset{0}{\overset{t}{\int }}\xi (s,t)\mathsf{\Delta }\overrightarrow{u}(\overrightarrow{x},s)\mathrm{d}s+\lambda \overrightarrow{u}+\nabla ⟨\pi ⟩=\overrightarrow{f},}\hfill & \hfill \\ {\gamma }_{\partial \mathsf{\Omega }}\overrightarrow{u}={\overrightarrow{u}}_b,\hfill & \hfill \\ \overrightarrow{u}\left(0\right)={\overrightarrow{u}}_0.\hfill \end{array}\right.$

The first equality in the last system is derived from (1) by taking into account the easy-to-verify relation

$\nabla ·\left(\eta \mathbb{D}\left(\overrightarrow{u}\right)\right)={(\nabla \eta )}^{\top }\mathbb{D}\left(\overrightarrow{u}\right)+{\displaystyle \frac{\eta }{2}}\mathsf{\Delta }\overrightarrow{u}.$

Note that the above definition differs in a fashion from the usual definition of strong solutions to IBVPs for Navier–Stokes–Voigt-type equations (cf. [11,44]). In particular, we have required $\overrightarrow{u}\in {\mathit{C}}^1\left(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)\right)$ instead of $\overrightarrow{u}\in {\mathbf{W}}^{1,2}\left(S,{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)\right)$, that is, the velocity field $\overrightarrow{u}$ is assumed to be $C^1$-smooth with respect to time t. We are able to construct such regular solutions due to the conditions given in (13). At the same time, the main results of this paper also hold for “standard” strong solutions subject to appropriate changes in the problem statement and in (13).

4. Main Results

Let ${\eta }_{\mathrm{aver}}$ be a constant, which represents a typical viscosity value in the flow model (1).

Theorem 4

(Existence and uniqueness of strong solutions). Suppose that

• Ω is a bounded domain in space ${\mathbb{R}}^d$, where $d=2$ or 3, and $\partial \mathsf{\Omega }\in C^2;$

• ${\eta }_{\mathrm{aver}}>0$, $\alpha >0$ and $\lambda \ge 0;$

• all inclusions in (13) are valid;

• the vector functions ${\overrightarrow{u}}_0$ and ${\overrightarrow{u}}_b$ satisfy the compatibility condition

(14)${\gamma }_{\partial \mathsf{\Omega }}{\overrightarrow{u}}_0={\overrightarrow{u}}_b\left(0\right).$

Then the following two statements hold:

• (a). There exist two positive constants ${\epsilon }_1(Q,\alpha ,{\eta }_{\mathrm{aver}},\lambda )$ and ${\epsilon }_2(Q,\alpha ,{\eta }_{\mathrm{aver}},\eta ,\xi ,\lambda )$ such that if the norms of η, ξ, $\overrightarrow{f}$, ${\overrightarrow{u}}_b$ and ${\overrightarrow{u}}_0$ satisfy the following estimates:

  (15)$\begin{array}{cc}\hfill {\parallel \nabla \eta \parallel }_{\mathit{C}(\overline{S},{\mathit{L}}^4\left(\mathsf{\Omega }\right))}+\parallel \eta -{\eta }_{\mathrm{aver}}{\parallel }_{C\left(\overline{Q}\right)}+{\parallel \xi \parallel }_{C(\overline{S}\times \overline{S})}& \le {\epsilon }_1(Q,\alpha ,{\eta }_{\mathrm{aver}},\lambda ),\hfill \end{array}$

  (16)$\begin{array}{cc}\hfill \parallel \overrightarrow{f}{\parallel }_{\mathit{C}(\overline{S},{\mathit{L}}^2\left(\mathsf{\Omega }\right))}+\parallel {\overrightarrow{u}}_b{\parallel }_{{\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}(\partial \mathsf{\Omega }))}+{\parallel {\overrightarrow{u}}_0\parallel }_{{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)}& \le {\epsilon }_2(Q,\alpha ,{\eta }_{\mathrm{aver}},\eta ,\xi ,\lambda ),\hfill \end{array}$

  then IBVP (1)–(3) has a unique strong solution in an open neighborhood $\mathcal{U}(Q,\alpha ,{\eta }_{\mathrm{aver}},\eta ,\xi ,\lambda )$ of the pair $\left(\overrightarrow{0},⟨0⟩\right)$ in the Cartesian product ${\mathit{C}}^1\left(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)\right)\times C\left(\overline{S},H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)\right)$.

• (b). If $\delta =0$ and estimate (15) is valid, then IBVP (1)–(3) is uniquely solvable in the strong formulation without the smallness requirement (16) on the norms of vector functions $\overrightarrow{f}$, ${\overrightarrow{u}}_b$, and ${\overrightarrow{u}}_0$.

The proof of Theorem 4 is given in Section 6. This proof is based on the operator treatment of IBVP (1)–(3) and appropriately applying the abstract result formulated in Theorem 3.

Remark 1.

Conditions (15) and (16) imply, in particular, the smallness of viscosity variations, external body forces and both boundary and initial velocities, which corresponds to the case of laminar (slow) flows.

5. Auxiliary Propositions

5.1. Boundary Trace and Divergence-Free Lifting

For $C^1$-smooth functions defined on the time interval $\overline{S}$ with values in the Hilbert space ${\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)$, one can introduce the boundary trace operator ${\mathbf{Tr}}_{\partial \mathsf{\Omega }}$ by

$\begin{array}{cc}\hfill & {\mathbf{Tr}}_{\partial \mathsf{\Omega }}:{\mathit{C}}^1\left(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)\right)\to {\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right)),\hfill \\ \hfill & \left[{\mathbf{Tr}}_{\partial \mathsf{\Omega }}\overrightarrow{u}\right]\left(t\right):={\gamma }_{\partial \mathsf{\Omega }}\left[\overrightarrow{u}\left(t\right)\right],\mathrm{f}\mathrm{o}\mathrm{r} \mathrm{a}\mathrm{n}\mathrm{y}\overrightarrow{u}\in {\mathit{C}}^1\left(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)\right) \mathrm{a}\mathrm{n}\mathrm{d} t\in \overline{S}.\hfill \end{array}$

Below, we will construct a continuous right inverse for the operator ${\mathbf{Tr}}_{\partial \mathsf{\Omega }}$.

Proposition 1

(Existence of a divergence-free lifting). Let $\partial \mathsf{\Omega }\in C^2$. Then, there exists a lifting operator ${\mathbf{Lif}}_{\mathsf{\Omega }}$ such that

(17)$\begin{array}{c}{\mathbf{Lif}}_{\mathsf{\Omega }}\in \mathcal{L}({\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right)),{\mathit{C}}^1(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right))),\end{array}$

(18)$\begin{array}{c}{\mathbf{Tr}}_{\partial \mathsf{\Omega }}\circ {\mathbf{Lif}}_{\mathsf{\Omega }}=\mathbf{I},\end{array}$

where $\mathbf{I}$ denotes the identity operator in the space ${\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right))$.

Proof. 

Let us fix a vector function $\overrightarrow{\phi }$ belonging to the space ${\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right)$ and consider a non-homogeneous Dirichlet boundary value problem for the Stokes-type system:

(19)$\left\{\begin{array}{cc}{\mathbf{P}}_{\mathit{H}\left(\mathsf{\Omega }\right)}\mathsf{\Delta }\overrightarrow{\mathit{v}}=\overrightarrow{0}\hfill & \mathrm{i}\mathrm{n} \mathsf{\Omega },\hfill \\ \mathrm{div}\left(\overrightarrow{\mathit{v}}\right)=0\hfill & \mathrm{i}\mathrm{n} \mathsf{\Omega },\hfill \\ \overrightarrow{\mathit{v}}=\overrightarrow{\phi }\hfill & \mathrm{o}\mathrm{n} \partial \mathsf{\Omega }.\hfill \end{array}\right.$

Here, ${\mathbf{P}}_{\mathit{H}\left(\mathsf{\Omega }\right)}:{\mathit{L}}^2\left(\mathsf{\Omega }\right)\to \mathit{H}\left(\mathsf{\Omega }\right)$ is the Leray projection (see Section 2.6).

Taking into account the known results about the well-posedness of the stationary Stokes equations with non-homogeneous Dirichlet boundary conditions (see, for example, the monograph [7], Chapter I), we conclude that, in the space ${\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)$, there exists a unique vector function $\overrightarrow{\mathit{v}}$ satisfying system (19).

Moreover, there exists a positive constant $c_0\left(\mathsf{\Omega }\right)$ such that

(20)$\parallel \overrightarrow{v}{\parallel }_{{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)}\le c_0\left(\mathsf{\Omega }\right){\parallel \overrightarrow{\phi }\parallel }_{{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right)}.$

For problem (19), we introduce the data-to-solution mapping $\mathit{A}$ as follows:

$\mathit{A}:{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right)\to {\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right),\mathit{A}\overrightarrow{\phi }:=\overrightarrow{v}.$

This operator is well defined and, due to estimate (20), we have the inclusion

$\mathit{A}\in \mathcal{L}({\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right),{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)).$

Now, we can define the operator ${\mathbf{Lif}}_{\mathsf{\Omega }}$ by

$\begin{array}{c}{\mathbf{Lif}}_{\mathsf{\Omega }}:{\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right))\to {\mathit{C}}^1(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)),\\ [{\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\psi }]\left(t\right):=\mathit{A}[\overrightarrow{\psi }\left(t\right)], \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{a}\mathrm{n}\mathrm{y} \overrightarrow{\psi }\in {\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right)) \mathrm{a}\mathrm{n}\mathrm{d} t\in \overline{S}.\end{array}$

It is easily shown that relations (17) and (18) hold for the operator ${\mathbf{Lif}}_{\mathsf{\Omega }}$. Thus, Proposition 1 is proved. □

We will refer to ${\mathbf{Lif}}_{\mathsf{\Omega }}$ as a divergence-free lifting operator.

Clearly, the operators ${\mathbf{Lif}}_{\mathsf{\Omega }}$ and ${\gamma }_{\partial \mathsf{\Omega }}$ satisfy the following equality:

${\gamma }_{\partial \mathsf{\Omega }}[{\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\psi }\left(t\right)]=\overrightarrow{\psi }\left(t\right), \mathrm{f}\mathrm{o}\mathrm{r} \mathrm{a}\mathrm{n}\mathrm{y} \overrightarrow{\psi }\in {\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}(\partial \mathsf{\Omega })) \mathrm{a}\mathrm{n}\mathrm{d} t\in \overline{S}.$

5.2. Two Linear Operators Associated with Kelvin–Voigt–Brinkman–Forchheimer System

Let us introduce the two main functional spaces:

$\begin{array}{cc}\hfill & \mathit{X}\left(Q\right):={\mathit{C}}^1\left(\overline{S},{\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right)\right)\times C\left(\overline{S},H^1\left(\mathsf{\Omega }\right)/R\left(\mathsf{\Omega }\right)\right),\hfill \\ \hfill & \mathit{Y}\left(Q\right):=\{(\overrightarrow{g},\overrightarrow{\varphi },\overrightarrow{a})\in \mathit{C}\left(\overline{S},{\mathit{L}}^2\left(\mathsf{\Omega }\right)\right)\times {\mathit{C}}^1(\overline{S},{\dot{\mathit{H}}}^{\frac{3}{2}}\left(\partial \mathsf{\Omega }\right))\times {\mathit{H}}_{\mathrm{div}}^2\left(\mathsf{\Omega }\right):\overrightarrow{\varphi }\left(0\right)={\gamma }_{\partial \mathsf{\Omega }}\overrightarrow{a}\}.\hfill \end{array}$

In order to analyze properties of the linear part of IBVP (1)–(3), we define two operators ${\mathit{T}}_0:\mathit{X}\left(Q\right)\to \mathit{Y}\left(Q\right)$ and ${\mathit{T}}_1:\mathit{X}\left(Q\right)\to \mathit{Y}\left(Q\right)$ as follows:

(21)$\begin{array}{cc}\hfill & {\mathit{T}}_0\left(\overrightarrow{u},⟨\pi ⟩\right):=\left({\overrightarrow{u}}^{\prime }-\alpha \mathsf{\Delta }{\overrightarrow{u}}^{\prime }-\frac{{\eta }_{\mathrm{aver}}}{2}\mathsf{\Delta }\overrightarrow{u}+\lambda \overrightarrow{u}+\nabla ⟨\pi ⟩,{\mathbf{Tr}}_{\partial \mathsf{\Omega }}\overrightarrow{u},\overrightarrow{u}\left(0\right)\right),\hfill \end{array}$

(22)$\begin{array}{cc}\hfill & {\mathit{T}}_1\left(\overrightarrow{u},⟨\pi ⟩\right):=\left(-{(\nabla \eta )}^{\top }\mathbb{D}\left(\overrightarrow{u}\right)-{\displaystyle \frac{\eta -{\eta }_{\mathrm{aver}}}{2}}\mathsf{\Delta }\overrightarrow{u}-\underset{0}{\overset{t}{\int }}\xi (s,t)\mathsf{\Delta }\overrightarrow{u}(\overrightarrow{x},s)\mathrm{d}s,\overrightarrow{0},\overrightarrow{0}\right).\hfill \end{array}$

Proposition 2.

Suppose that Ω is a bounded domain in space ${\mathbb{R}}^d$, $d=2,3$, with boundary $\partial \mathsf{\Omega }$ of class $C^2$ and, moreover, the inequalities ${\eta }_{\mathrm{aver}}>0$, $\alpha >0$ and $\lambda \ge 0$ hold. Then, the operator ${\mathit{T}}_0$ is an isomorphism, that is,

${\mathit{T}}_0\in {\mathcal{L}}_{\mathrm{Isom}}\left(\mathit{X}\left(Q\right),\mathit{Y}\left(Q\right)\right).$

Proof. 

We divide the proof of this proposition in three steps.

Step 1: Continuity property.

It is easy to show that the following estimate holds:

$\parallel {\mathit{T}}_0\left(\overrightarrow{u},⟨\pi ⟩\right){\parallel }_{\mathit{Y}\left(Q\right)}\le c(Q,\alpha ,{\eta }_{\mathrm{aver}},\lambda )\parallel \left(\overrightarrow{u},⟨\pi ⟩\right){\parallel }_{\mathit{X}\left(Q\right)},\forall \left(\overrightarrow{u},⟨\pi ⟩\right)\in \mathit{X}\left(Q\right),$

with some positive constant $c(Q,\alpha ,{\eta }_{\mathrm{aver}},\lambda )$. Therefore, the linear operator ${\mathit{T}}_0$ is continuous, that is, the inclusion

${\mathit{T}}_0\in \mathcal{L}\left(\mathit{X}\left(Q\right),\mathit{Y}\left(Q\right)\right)$

is valid.

Step 2: Injectivity property.

Now, we prove that ${\mathit{T}}_0$ is an injective mapping. Suppose that $\left({\overrightarrow{u}}_1,⟨{\pi }_1⟩\right)$ and $\left({\overrightarrow{u}}_2,⟨{\pi }_2⟩\right)$ are pairs such that

$\left({\overrightarrow{u}}_1,⟨{\pi }_1⟩\right)\in \mathit{X}\left(Q\right),\left({\overrightarrow{u}}_2,⟨{\pi }_2⟩\right)\in \mathit{X}\left(Q\right)$

and

${\mathit{T}}_0\left({\overrightarrow{u}}_1,⟨{\pi }_1⟩\right)={\mathit{T}}_0\left({\overrightarrow{u}}_2,⟨{\pi }_2⟩\right).$

Let us show that

(23)$\left({\overrightarrow{u}}_1,⟨{\pi }_1⟩\right)=\left({\overrightarrow{u}}_2,⟨{\pi }_2⟩\right).$

Clearly, if

$\left({\overrightarrow{u}}_{★},⟨{\pi }_{★}⟩\right):=\left({\overrightarrow{u}}_1-{\overrightarrow{u}}_2,⟨{\pi }_1⟩-⟨{\pi }_2⟩\right),$

then we have

${\mathit{T}}_0\left({\overrightarrow{u}}_{★},⟨{\pi }_{★}⟩\right)=(\overrightarrow{0},\overrightarrow{0},\overrightarrow{0}).$

The last relation yields that

(24)${\overrightarrow{u}}_{★}^{\prime }-\alpha \mathsf{\Delta }{\overrightarrow{u}}_{★}^{\prime }-{\displaystyle \frac{{\eta }_{\mathrm{aver}}}{2}}\mathsf{\Delta }{\overrightarrow{u}}_{★}+\lambda {\overrightarrow{u}}_{*}+\nabla ⟨{\pi }_{★}⟩=\overrightarrow{0}inQ.$

Multiplying both sides of (24) by the vector function ${\overrightarrow{u}}_{★}$ and integrating the domain $\mathsf{\Omega }$ lead us to

(25)$\begin{array}{cc}\hfill \underset{\mathsf{\Omega }}{\int }{\overrightarrow{u}}_{★}^{\prime }·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}& -\alpha \underset{\mathsf{\Omega }}{\int }\mathsf{\Delta }{\overrightarrow{u}}_{★}^{\prime }·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}-{\displaystyle \frac{{\eta }_{\mathrm{aver}}}{2}}\underset{\mathsf{\Omega }}{\int }\mathsf{\Delta }{\overrightarrow{u}}_{★}·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & +\lambda \underset{\mathsf{\Omega }}{\int }{\overrightarrow{u}}_{★}·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}+\underset{\mathsf{\Omega }}{\int }\nabla ⟨{\pi }_{★}⟩·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}=0,\forall t\in \overline{S}.\hfill \end{array}$

Taking into account the equalities $\mathrm{div}\left({\overrightarrow{u}}_{★}\right)=0$ in Q and ${\gamma }_{\partial \mathsf{\Omega }}{\overrightarrow{u}}_{★}=\overrightarrow{0}$ on $\partial \mathsf{\Omega }\times S$, by integration by parts, we obtain

(26)$\begin{array}{cc}\hfill \underset{\mathsf{\Omega }}{\int }\mathsf{\Delta }{\overrightarrow{u}}_{★}^{\prime }·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}& =\sum _{i,j=1}^d\underset{\mathsf{\Omega }}{\int }{\displaystyle \frac{{\partial }^2{u}_{★i}^{\prime }}{\partial x_j^2}}{u}_{★i}\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & =\sum _{i,j=1}^d\underset{\partial \mathsf{\Omega }}{\int }\underset{=0}{\underbrace{{\displaystyle \frac{\partial u_{★i}^{\prime }}{\partial x_j}}{u}_{★i}{n}_j}}\mathrm{d}\sigma -\sum _{i,j=1}^d\underset{\mathsf{\Omega }}{\int }{\displaystyle \frac{\partial u_{★i}^{\prime }}{\partial x_j}}{\displaystyle \frac{\partial u_{★i}}{\partial x_j}}\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & =-\underset{\mathsf{\Omega }}{\int }\nabla {\overrightarrow{u}}_{★}^{\prime }:\nabla {\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x},\hfill \end{array}$

(27)$\begin{array}{cc}\hfill \underset{\mathsf{\Omega }}{\int }\mathsf{\Delta }{\overrightarrow{u}}_{★}·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}& =\sum _{i,j=1}^d\underset{\mathsf{\Omega }}{\int }{\displaystyle \frac{{\partial }^2{u}_{★i}}{\partial x_j^2}}{u}_{★i}\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & =\sum _{i,j=1}^d\underset{\partial \mathsf{\Omega }}{\int }\underset{=0}{\underbrace{{\displaystyle \frac{\partial u_{★i}}{\partial x_j}}{u}_{★i}{n}_j}}\mathrm{d}\sigma -\sum _{i,j=1}^d\underset{\mathsf{\Omega }}{\int }{\left({\displaystyle \frac{\partial u_{★i}}{\partial x_j}}\right)}^2\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & =-\underset{\mathsf{\Omega }}{\int }{|\nabla {\overrightarrow{u}}_{★}|}^2\mathrm{d}\overrightarrow{x},\hfill \end{array}$

(28)$\begin{array}{cc}\hfill \underset{\mathsf{\Omega }}{\int }\nabla ⟨{\pi }_{★}⟩·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}& =\sum _{i=1}^d\underset{\mathsf{\Omega }}{\int }{\displaystyle \frac{\partial ({\pi }_1-{\pi }_2)}{\partial x_i}}{u}_{★i}\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & =\sum _{i=1}^d\underset{\partial \mathsf{\Omega }}{\int }\underset{=0}{\underbrace{({\pi }_1-{\pi }_2)u_{★i}{n}_i}}\mathrm{d}\sigma -\sum _{i=1}^d\underset{\mathsf{\Omega }}{\int }({\pi }_1-{\pi }_2){\displaystyle \frac{\partial u_{★i}}{\partial x_i}}\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & =\underset{\mathsf{\Omega }}{\int }({\pi }_1-{\pi }_2)\underset{=0}{\underbrace{\mathrm{div}\left({\overrightarrow{u}}_{★}\right)}}\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & =0.\hfill \end{array}$

Next, we substitute (26)–(28) into relation (25). This yields that

(29)$\begin{array}{cc}\hfill \underset{\mathsf{\Omega }}{\int }{\overrightarrow{u}}_{★}^{\prime }·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}+\alpha \underset{\mathsf{\Omega }}{\int }\nabla {\overrightarrow{u}}_{★}^{\prime }:\nabla {\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}& +{\displaystyle \frac{{\eta }_{\mathrm{aver}}}{2}}\underset{\mathsf{\Omega }}{\int }{|\nabla {\overrightarrow{u}}_{★}|}^2\mathrm{d}\overrightarrow{x}\hfill \\ \hfill & +\lambda \underset{\mathsf{\Omega }}{\int }{\left|{\overrightarrow{u}}_{★}\right|}^2\mathrm{d}\overrightarrow{x}=0,\forall t\in \overline{S}.\hfill \end{array}$

Using the two easy-to-verify relations:

$\underset{\mathsf{\Omega }}{\int }{\overrightarrow{u}}_{★}^{\prime }·{\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}={\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}\underset{\mathsf{\Omega }}{\int }{\left|{\overrightarrow{u}}_{★}\right|}^2\mathrm{d}\overrightarrow{x},\forall t\in \overline{S},$

and

$\underset{\mathsf{\Omega }}{\int }\nabla {\overrightarrow{u}}_{★}^{\prime }:\nabla {\overrightarrow{u}}_{★}\mathrm{d}\overrightarrow{x}={\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}\underset{\mathsf{\Omega }}{\int }{|\nabla {\overrightarrow{u}}_{★}|}^2\mathrm{d}\overrightarrow{x},\forall t\in \overline{S},$

we rewrite equality (29) as follows:

(30)${\displaystyle \frac{d}{dt}}\underset{\mathsf{\Omega }}{\int }(|{\overrightarrow{u}}_{★}{|}^2+\alpha {|\nabla {\overrightarrow{u}}_{★}|}^2)\mathrm{d}\overrightarrow{x}+{\eta }_{\mathrm{aver}}\underset{\mathsf{\Omega }}{\int }|\nabla {\overrightarrow{u}}_{★}{|}^2\mathrm{d}\overrightarrow{x}+\lambda \underset{\mathsf{\Omega }}{\int }{|{\overrightarrow{u}}_{★}|}^2\mathrm{d}\overrightarrow{x}=0,\forall t\in \overline{S}.$

Furthermore, we integrate both sides of (30) with respect to the variable t from 0 to $\tau$, where $\tau$ is an arbitrary point belonging to the time interval S. Taking into account the equality ${\overrightarrow{u}}_{★}\left(0\right)=\overrightarrow{0}$, we obtain

$\begin{array}{cc}\hfill \underset{\mathsf{\Omega }}{\int }\left(|{\overrightarrow{u}}_{★}{\left(\tau \right)|}^2+\alpha {|\nabla {\overrightarrow{u}}_{★}\left(\tau \right)|}^2\right)\mathrm{d}\overrightarrow{x}& +{\eta }_{\mathrm{aver}}\underset{0}{\overset{\tau }{\int }}\underset{\mathsf{\Omega }}{\int }{|\nabla {\overrightarrow{u}}_{★}\left(t\right)|}^2\mathrm{d}\overrightarrow{x}\mathrm{d}t\hfill \\ \hfill & +\lambda \underset{0}{\overset{\tau }{\int }}\underset{\mathsf{\Omega }}{\int }{\left|{\overrightarrow{u}}_{★}\right|}^2\mathrm{d}\overrightarrow{x}\mathrm{d}t=0,\forall \tau \in \overline{S}.\hfill \end{array}$

This implies that

(31)${\overrightarrow{u}}_{★}\left(\tau \right)=\overrightarrow{0},\forall \tau \in \overline{S},$

and hence ${\overrightarrow{u}}_1={\overrightarrow{u}}_2$.

Moreover, from relations (24) and (31), it follows that the equality $\nabla ⟨{\pi }_{★}⟩=\overrightarrow{0}$ holds. Due to this equality, we conclude that $⟨{\pi }_{★}⟩=⟨0⟩$, whence $⟨{\pi }_1⟩=⟨{\pi }_2⟩$. Thus, we have established the required equality (23), which means that ${\mathit{T}}_0$ is an injective operator.

Step 3: Surjectivity property.

Due to the properties of the operator ${\mathit{T}}_0$ that have been established in Steps 1 and 2 and Theorem 2, for completing of the proof of Proposition 2, it is sufficient to show that this operator is surjective.

Let us consider an arbitrary triple $(\overrightarrow{g},\overrightarrow{\varphi },\overrightarrow{a})$ belonging to the space $\mathit{Y}\left(Q\right)$. We will prove that the operator equation

(32)${\mathit{T}}_0(\overrightarrow{u},⟨\pi ⟩)=(\overrightarrow{g},\overrightarrow{\varphi },\overrightarrow{a})$

has a solution in the space $\mathit{X}\left(Q\right)$.

Clearly, Equation (32) is equivalent to the following IBVP:

(33)$\left\{\begin{array}{cc}{\displaystyle {\overrightarrow{u}}^{\prime }-\alpha \mathsf{\Delta }{\overrightarrow{u}}^{\prime }-\frac{{\eta }_{\mathrm{aver}}}{2}\mathsf{\Delta }\overrightarrow{u}+\lambda \overrightarrow{u}+\nabla ⟨\pi ⟩=\overrightarrow{g}}\hfill & \mathrm{i}\mathrm{n} Q,\hfill \\ \mathrm{div}\left(\overrightarrow{u}\right)=0\hfill & \mathrm{i}\mathrm{n} Q,\hfill \\ \overrightarrow{u}=\overrightarrow{\varphi }\hfill & \mathrm{o}\mathrm{n} \partial \mathsf{\Omega }\times S,\hfill \\ \overrightarrow{u}{|}_{t=0}=\overrightarrow{a}\hfill & \mathrm{i}\mathrm{n} \mathsf{\Omega }\hfill \end{array}\right.$

provided that it is considered in the strong formulation (in the sense of Definition 3).

Further, we introduce a new unknown vector function $\overrightarrow{v}$, which is connected with $\overrightarrow{u}$ by the relation $\overrightarrow{u}=\overrightarrow{v}+{\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\varphi }$. Using this notation, one can rewrite IBVP (33) as follows:

(34)$\left\{\begin{array}{cc}{\displaystyle {\overrightarrow{v}}^{\prime }-\alpha \mathsf{\Delta }{\overrightarrow{v}}^{\prime }-\frac{{\eta }_{\mathrm{aver}}}{2}\mathsf{\Delta }\overrightarrow{v}+\lambda \overrightarrow{v}+\nabla ⟨\pi ⟩={\overrightarrow{g}}_{*}}\hfill & \mathrm{i}\mathrm{n} Q,\hfill \\ \mathrm{div}\left(\overrightarrow{v}\right)=0\hfill & \mathrm{i}\mathrm{n} Q,\hfill \\ \overrightarrow{v}=\overrightarrow{0}\hfill & \mathrm{o}\mathrm{n} \partial \mathsf{\Omega }\times S,\hfill \\ \overrightarrow{v}{|}_{t=0}={\overrightarrow{a}}_{*}\hfill & \mathrm{i}\mathrm{n} \mathsf{\Omega },\hfill \end{array}\right.$

where

$\begin{array}{cc}\hfill {\overrightarrow{g}}_{*}& :=\overrightarrow{g}-{\mathbf{Lif}}_{\mathsf{\Omega }}{\overrightarrow{\varphi }}^{\prime }+\alpha \mathsf{\Delta }({\mathbf{Lif}}_{\mathsf{\Omega }}{\overrightarrow{\varphi }}^{\prime })+{\displaystyle \frac{{\eta }_{\mathrm{aver}}}{2}}\mathsf{\Delta }({\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\varphi })-\lambda {\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\varphi }\in \mathit{C}\left(\overline{S},{\mathit{L}}^2\left(\mathsf{\Omega }\right)\right),\hfill \\ \hfill {\overrightarrow{a}}_{*}& :=\overrightarrow{a}-[{\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\varphi }]\left(0\right)\in {\mathit{H}}_{0,\mathrm{div}}^2\left(\mathsf{\Omega }\right).\hfill \end{array}$

Applying a modified Faedo–Galerkin scheme (with the basis of eigenfunctions of the Stokes operator) as in the proof of the existence result obtained in [13], one can establish that IBVP (34) has a unique strong solution $\left({\overrightarrow{v}}_{\#},⟨{\pi }_{\#}⟩\right)$. Clearly, the pair $\left({\overrightarrow{v}}_{\#}+{\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\varphi },⟨{\pi }_{\#}⟩\right)$ is a strong solution of IBVP (33). Since this IBVP and Equation (32) are equivalent, we see that $\left({\overrightarrow{v}}_{\#}+{\mathbf{Lif}}_{\mathsf{\Omega }}\overrightarrow{\varphi },⟨{\pi }_{\#}⟩\right)$ is a solution of (32). In view of the arbitrariness of the triple $(\overrightarrow{g},\overrightarrow{\varphi },\overrightarrow{a})$, this yields that the operator ${\mathit{T}}_0$ is surjective. Thus, the proof of Proposition 2 is complete. □

An important consequence of Proposition 2 is formulated in the following statement.

Proposition 3.

Let all conditions of Proposition 2 be valid and, moreover,

$\xi \in C(\overline{S}\times \overline{S}),\eta \in C\left(\overline{Q}\right),\nabla \eta \in \mathit{C}\left(\overline{Q}\right).$

Then, the sum of the operators ${\mathit{T}}_0$ and ${\mathit{T}}_1$ is an isomorphism between the spaces $\mathit{X}\left(Q\right)$ and $\mathit{Y}\left(Q\right)$, that is,

$({\mathit{T}}_0+{\mathit{T}}_1)\in {\mathcal{L}}_{\mathrm{Isom}}\left(\mathit{X}\left(Q\right),\mathit{Y}\left(Q\right)\right),$

provided that the value of the sum of the norms

${\parallel \nabla \eta \parallel }_{\mathit{C}(\overline{S},{\mathit{L}}^4\left(\mathsf{\Omega }\right))}+\parallel \eta -{\eta }_{\mathrm{aver}}{\parallel }_{C\left(\overline{Q}\right)}+{\parallel \xi \parallel }_{C(\overline{S}\times \overline{S})}$

is sufficiently small.

The proof of the above proposition is based on applying Theorem 1.

6. Proof of Main Theorem

In this section, we prove Theorem 4, in which our main results are stated.

First, let us give the operator treatment IBVP (1)–(3). To this end, we need to introduce one extra operator, which is associated with the convective term $\nabla ·(\overrightarrow{u}\otimes \overrightarrow{u})$ and, unlike the operators ${\mathit{T}}_0$ and ${\mathit{T}}_1$, is nonlinear:

${\mathit{T}}_2:\mathit{X}\left(Q\right)\to \mathit{Y}\left(Q\right),{\mathit{T}}_2(\overrightarrow{u},⟨\pi ⟩):=(\nabla ·(\overrightarrow{u}\otimes \overrightarrow{u}),\overrightarrow{0},\overrightarrow{0}).$

It is easily seen that IBVP (1)–(3) in the strong formulation (in the sense of Definition 3) is equivalent to the operator equation

(35)$({\mathit{T}}_0+{\mathit{T}}_1+\delta {\mathit{T}}_2)(\overrightarrow{u},⟨\pi ⟩)=(\overrightarrow{f},{\overrightarrow{u}}_b,{\overrightarrow{u}}_0)$

with the parameter $\delta \in \{0,1\}$. Recall that the operators ${\mathit{T}}_0$ and ${\mathit{T}}_1$ are defined in formulas (21) and (), respectively.

Now, we separately consider the two cases: $\delta =1$ and $\delta =0$.

Case 1: $\delta =1$ (the problem with the convective term).

In view of $\partial \mathsf{\Omega }\in C^2$, ${\eta }_{\mathrm{aver}}>0$, $\alpha >0$ and $\lambda \ge 0$, by Proposition 2, we conclude that the operator ${\mathit{T}}_0$ is an isomorphism. Furthermore, it can easily be checked that there exists a positive constant ${\epsilon }_1(Q,\alpha ,{\eta }_{\mathrm{aver}},\lambda )$ such that, if condition (15) is valid, then

$\parallel {\mathit{T}}_0^{-1}{\parallel }_{\mathcal{L}\left(\mathit{Y}\right(Q),\mathit{X}(Q\left)\right)}{\parallel {\mathit{T}}_1\parallel }_{\mathcal{L}\left(\mathit{X}\right(Q),\mathit{Y}(Q\left)\right)}<1.$

Clearly, ${\mathit{T}}_2\left(\overrightarrow{0},⟨0⟩\right)=(\overrightarrow{0},\overrightarrow{0},\overrightarrow{0})$. We also note that ${\mathit{T}}_2$ is a continuously differentiable operator and its Fréchet derivative $\mathit{D}{\mathit{T}}_2$ can be calculated by the following formula

(36)$\left[\mathit{D}{\mathit{T}}_2\left(\overrightarrow{u},⟨\pi ⟩\right)\right]\left(\overrightarrow{w},⟨\sigma ⟩\right)=\left(\nabla ·(\overrightarrow{w}\otimes \overrightarrow{u})+\nabla ·(\overrightarrow{u}\otimes \overrightarrow{w}),\overrightarrow{0},\overrightarrow{0}\right),$

for any pairs $\left(\overrightarrow{u},⟨\pi ⟩\right)$ and $\left(\overrightarrow{w},⟨\sigma ⟩\right)$ belonging to the space $\mathit{X}\left(Q\right)$.

From equality (36) it follows that the Fréchet derivative of ${\mathit{T}}_2$ at the zero pair $\left(\overrightarrow{0},⟨0⟩\right)$ is the zero operator.

Taking into account the above observations, we can apply Theorem 3 to Equation (35) with $\delta =1$ and conclude that IBVP (1)–(3) has a unique strong solution $\left(\overrightarrow{u},⟨\pi ⟩\right)$ in an open neighborhood $\mathcal{U}(Q,\alpha ,{\eta }_{\mathrm{aver}},\eta ,\xi ,\lambda )$ of the pair of the zero functions $\left(\overrightarrow{0},⟨0⟩\right)$ in the space $\mathit{X}\left(Q\right)$ provided that conditions (14)–() hold.

Case 2: $\delta =0$ (the problem without the convective term).

In this case, IBVP (1)–(3) is linear and its unique solvability follows directly from Proposition 3, without any smallness assumptions on the model data $\overrightarrow{f}$, ${\overrightarrow{u}}_b$, and ${\overrightarrow{u}}_0$, but provided that compatibility condition (14) and estimate (15) are valid. Thus, Theorem 4 is completely proven.

Remark 2.

As can be seen from the above proof, the approach based on the use of Theorem 3 is quite universal and can be applied in the analysis of the local well-posedness of various other models for fluid motion [45,46,47,48,49,50,51,52,53,54] , heat and mass transfer models [55,56,57,58,59,60,61,62,63,64] and magnetohydrodynamics problems [65,66,67,68].

7. Final Comments

In this paper, we have established sufficient conditions for the unique solvability of the IBVP for the Kelvin–Voigt–Brinkman–Forchheimer equations describing unsteady flows of a variable viscosity fluid with memory through a porous media. More precisely, the existence and uniqueness of a regular-in-time strong solution satisfying non-zero Dirichlet boundary conditions have been proved for small model data, which corresponds to the case of laminar (slow) flows. Our approach is based on the application of a theorem about the local unique solvability of operator equations involving an isomorphism between Banach spaces with linear and nonlinear Fréchet differentiable perturbations. The present work can be considered as the first step in a comprehensive study of the inhomogeneous Kelvin–Voigt–Brinkman–Forchheimer system using methods of functional analysis. In this regard, open challenging problems related to proving the well-posedness and stability of IBVP (1)–(3) without any smallness assumptions for boundary, forcing and initial data, should be mentioned. Solving these problems is an important direction of future research and can ensure a deeper understanding of dynamics of non-Newtonian fluids with complex rheological properties.

Author Contributions

Conceptualization, E.S.B.; methodology, E.S.B.; investigation, E.S.B. and M.A.A.; writing—original draft preparation, E.S.B., M.A.A. and S.V.E.; visualization, E.S.B.; writing—review and editing, M.A.A., S.V.E. and A.V.Y. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figure and Table

Figure 1. Three-dimensional flow of a viscous fluid through a cylindrical domain [Forumla omitted. See PDF.] with an unmoved spherical obstacle contained inside this domain.

Table 1

Some particular cases of system (1).

Values of Parameters	Flow Model
$\delta =1$,  $\eta \equiv 0$,  $\alpha =0$,  $\xi \equiv 0$,  $\lambda =0$	The incompressible Euler equations for describing unsteady flows of an inviscid fluid [1,2,3,4,5]
$\delta =1$,  $\eta \equiv \mathrm{const}>0$,  $\alpha =0$,  $\xi \equiv 0$,  $\lambda =0$	The standard Navier–Stokes equations for describing Newtonian fluids motion [6,7,8,9]
$\delta =0$,  $\eta \equiv \mathrm{const}>0$,  $\alpha =0$,  $\xi \equiv 0$,  $\lambda =0$	The non-stationary Stokes equations for describing the creeping flows of an incompressible viscous fluid [6,7,8,9,10]
$\delta =1$,  $\eta \equiv \mathrm{const}>0$,  $\alpha >0$,  $\xi \equiv 0$,  $\lambda =0$	The incompressible Navier–Stokes–Voigt equations for describing the flows of Kelvin–Voigt-type viscoelastic fluids [11,12,13,14,15,16]
$\delta =1$,  $\eta \equiv \mathrm{const}>0$,  $\alpha >0$,  $\xi \equiv 0$,  $\lambda >0$	The Kelvin–Voigt–Brinkman–Forchheimer equations (without an integral memory term) for describing the dynamics of incompressible viscoelastic fluids in a porous medium [17,18,19,20,21,22,23]
$\delta =1$,  $\eta \equiv \mathrm{const}>0$,  $\alpha =0$,  $\xi (s,t)\equiv aexp\left(b\right(s-t))$,  $\lambda =0$	The model for describing flows of Jeffreys–Oldroyd-type viscoelastic fluids [24,25,26,27,28]
$\delta =1$,  $\eta \equiv \mathrm{const}>0$,  $\alpha >0$,  $\overline{)\xi }\cancel{\equiv }0$,  $\lambda =0$	The Oskolkov integro-differential system for describing the motion of non-Newtonian fluids with memory [29,30,31]
$\delta =1$,  $\eta \cancel{\equiv }\mathrm{const}$,  $\alpha =0$,  $\xi \equiv 0$,  $\lambda =0$	The generalized Navier–Stokes system for describing flows of non-homogeneous incompressible viscous fluids [32,33,34]
$\delta =0$,  $\eta \cancel{\equiv }\mathrm{const}$,  $\alpha =0$,  $\xi \equiv 0$,  $\lambda =0$	The generalized Stokes system for describing creeping flows of non-homogeneous incompressible viscous fluids [35]

References

1. Landau, L.D.; Lifshitz, E.M. Fluid Mechanics; 2nd ed. Pergamon Press: Oxford, UK, 1987; [DOI: https://dx.doi.org/10.1016/C2013-0-03799-1]

2. Chorin, A.J.; Marsden, J.E. A Mathematical Introduction to Fluid Mechanics; 3rd ed. Springer: New York, NY, USA, 1993; [DOI: https://dx.doi.org/10.1007/978-1-4612-0883-9]

3. Marchioro, C.; Pulvirenti, M. Mathematical Theory of Incompressible Nonviscous Fluids; Springer: New York, NY, USA, 1994; [DOI: https://dx.doi.org/10.1007/978-1-4612-4284-0]

4. Glass, O. Exact boundary controllability of 3-D Euler equation. ESAIM Control Optim. Calc. Var.; 2000; 5, pp. 1-44. [DOI: https://dx.doi.org/10.1051/cocv:2000100]

5. Ershkov, S.V.; Shamin, R.V. A Riccati-type solution of 3D Euler equations for incompressible flow. J. King Saud Univ.; 2020; 32, pp. 125-130. [DOI: https://dx.doi.org/10.1016/j.jksus.2018.03.010]

6. Ladyzhenskaya, O.A. The Mathematical Theory of Viscous Incompressible Flow; 2nd ed. Gordon and Breach: New York, NY, USA, 1969.

7. Temam, R. Navier–Stokes Equations: Theory and Numerical Analysis; North-Holland Publishing Co.: Amsterdam, The Netherlands, 1977.

8. Tartar, L. An Introduction to Navier–Stokes Equation and Oceanography; Springer: Berlin/Heidelberg, Germany, 2006; [DOI: https://dx.doi.org/10.1007/3-540-36545-1]

9. Boyer, F.; Fabrie, P. Mathematical Tools for the Study of the Incompressible Navier–Stokes Equations and Related Models; Springer: New York, NY, USA, 2013; [DOI: https://dx.doi.org/10.1007/978-1-4614-5975-0]

10. Varnhorn, W. The Stokes Equations; Akademie Verlag: Berlin, Germany, 1994.

11. Damazio, P.D.; Manholi, P.; Silvestre, A.L. L^q-theory of the Kelvin–Voigt equations in bounded domains. J. Differ. Equ.; 2016; 260, pp. 8242-8260. [DOI: https://dx.doi.org/10.1016/j.jde.2016.02.020]

12. Baranovskii, E.S. Mixed initial–boundary value problem for equations of motion of Kelvin–Voigt fluids. Comput. Math. Math. Phys.; 2016; 56, pp. 1363-1371. [DOI: https://dx.doi.org/10.1134/S0965542516070058]

13. Baranovskii, E.S. Strong solutions of the incompressible Navier–Stokes–Voigt model. Mathematics; 2020; 8, 181. [DOI: https://dx.doi.org/10.3390/math8020181]

14. Khompysh, K.; Kenzhebai, K. An inverse problem for Kelvin–Voigt equations perturbed by isotropic diffusion and damping. Math. Meth. Appl. Sci.; 2022; 45, pp. 3817-3842. [DOI: https://dx.doi.org/10.1002/mma.8018]

15. Straughan, B. Stability problems with generalized Navier–Stokes–Voigt theories. Ann. Univ. Ferrara; 2024; 70, pp. 1747-1766. [DOI: https://dx.doi.org/10.1007/s11565-024-00540-6]

16. Straughan, B. Double diffusion in a Navier–Stokes–Voigt fluid with a Christov heat law. Ann. Univ. Ferrara; 2025; 71, 23. [DOI: https://dx.doi.org/10.1007/s11565-024-00577-7]

17. Anh, C.T.; Trang, P.T. On the 3D Kelvin–Voigt–Brinkman–Forchheimer equations in some unbounded domains. Nonlinear Anal.; 2013; 89, pp. 36-54. [DOI: https://dx.doi.org/10.1016/j.na.2013.04.014]

18. Mohan, M.T. Global and exponential attractors for the 3D Kelvin–Voigt–Brinkman–Forchheimer equations. Discret. Contin. Dyn. Syst. Ser. B; 2020; 25, 9. [DOI: https://dx.doi.org/10.3934/dcdsb.2020067]

19. Gentile, M.; Straughan, B. Stabilization estimates for the Brinkman–Forchheimer–Kelvin–Voigt equation backward in time. Acta Mech.; 2023; 234, pp. 4001-4009. [DOI: https://dx.doi.org/10.1007/s00707-023-03592-5]

20. Badday, A.J.; Harfash, A.J. Thermosolutal convection in a Brinkman–Darcy–Kelvin–Voigt fluid with a bidisperse porous medium. Phys. Fluids; 2024; 36, 014119. [DOI: https://dx.doi.org/10.1063/5.0186934]

21. Caucao, S.; Yotov, I. Velocity-vorticity-pressure mixed formulation for the Kelvin–Voigt–Brinkman–Forchheimer equations. arXiv; 2024; [DOI: https://dx.doi.org/10.48550/arXiv.2406.16703] arXiv: 2406.16703

22. Shivaraj Kumar, D.L.; Basavaraj, M.S.; Aruna, A.S. Linear stability analysis of the viscoelastic Navier–Stokes–Voigt fluid model through Brinkman porous media: Modal and non-modal approaches. Int. J. Non-Linear Mech.; 2024; 167, 104885. [DOI: https://dx.doi.org/10.1016/j.ijnonlinmec.2024.104885]

23. Straughan, B. Kelvin–Voigt fluid models in double-diffusive porous convection. Transp. Porous Media; 2025; 152, 11. [DOI: https://dx.doi.org/10.1007/s11242-024-02147-z]

24. Oskolkov, A.P.; Shadiev, R. Towards a theory of global solvability on [0, ∞) of initial-boundary value problems for the equations of motion of Oldroyd and Kelvin–Voight fluids. J. Math. Sci.; 1994; 68, pp. 240-253. [DOI: https://dx.doi.org/10.1007/BF01249338]

25. Doubova, A.; Fernandez-Cara, E. On the control of viscoelastic Jeffreys fluids. Syst. Control Lett.; 2012; 61, pp. 573-579. [DOI: https://dx.doi.org/10.1016/j.sysconle.2012.02.003]

26. Zafar, A.A.; Fetecau, C.; Mirza, I.A. On the flow of Oldroyd-B fluids with fractional derivatives over a plate that applies shear stress to the fluid. Math. Rep.; 2016; 18, pp. 85-108.

27. Baranovskii, E.S.; Artemov, M.A. Global existence results for Oldroyd fluids with wall slip. Acta Appl. Math.; 2017; 147, pp. 197-210. [DOI: https://dx.doi.org/10.1007/s10440-016-0076-z]

28. Bir, B.; Goswami, D. On a three step two-grid finite element method for the Oldroyd model of order one. Z. Angew. Math. Mech.; 2021; 101, e202000373. [DOI: https://dx.doi.org/10.1002/zamm.202000373]

29. Oskolkov, A.P. Unsteady flows of viscoelastic fluids. Proc. Steklov Inst. Math.; 1984; 159, pp. 105-134.

30. Oskolkov, A.P. Theory of nonstationary flows of Kelvin–Voigt fluids. J. Math. Sci.; 1985; 28, pp. 751-758. [DOI: https://dx.doi.org/10.1007/BF02112340]

31. Baranovskii, E.S.; Artemov, M.A. Generalized Navier–Stokes equations with non-homogeneous boundary conditions. Fractal Fract.; 2022; 6, 373. [DOI: https://dx.doi.org/10.3390/fractalfract6070373]

32. Deteix, J.; Yakoubi, D. Improving the pressure accuracy in a projection scheme for incompressible fluids with variable viscosity. Appl. Math. Lett.; 2018; 79, pp. 111-117. [DOI: https://dx.doi.org/10.1016/j.aml.2017.12.004]

33. Plasman, L.; Deteix, J.; Yakoubi, D. A projection scheme for Navier–Stokes with variable viscosity and natural boundary condition. Int. J. Numer. Methods Fluids; 2020; 92, pp. 1845-1865. [DOI: https://dx.doi.org/10.1002/fld.4851]

34. Tartaglione, A. Regularity of weak solutions to the inhomogeneous stationary Navier–Stokes equations. Symmetry; 2021; 13, 1336. [DOI: https://dx.doi.org/10.3390/sym13081336]

35. Russo, R.; Tartaglione, A. The plane exterior boundary-value problem for nonhomogeneous fluids. J. Math. Fluid Mech.; 2020; 22, 14. [DOI: https://dx.doi.org/10.1007/s00021-019-0473-y]

36. Korobkov, M.V.; Pileckas, K.; Pukhnachov, V.V.; Russo, R. The flux problem for the Navier–Stokes equations. Russ. Math. Surv.; 2014; 69, pp. 1065-1122. [DOI: https://dx.doi.org/10.1070/RM2014v069n06ABEH004928]

37. Pukhnachev, V.V. Three dimensional flux problem for the Navier–Stokes equations. Vestn. Yuzhno-Ural. Univ. Ser. Mat. Model. Program.; 2015; 8, pp. 95-104. [DOI: https://dx.doi.org/10.14529/mmp150208]

38. Korobkov, M.; Pileckas, K.; Russo, R. The Steady Navier–Stokes System: Basics of the Theory and the Leray Problem; Birkhäuser: Cham, Switzerland, 2024; [DOI: https://dx.doi.org/10.1007/978-3-031-50898-1]

39. Di Plinio, F.; Giorgini, A.; Pata, V.; Temam, R. Navier–Stokes–Voigt equations with memory in 3D lacking instantaneous kinematic viscosity. J. Nonlinear Sci.; 2018; 28, pp. 653-686. [DOI: https://dx.doi.org/10.1007/s00332-017-9422-1]

40. Yushkov, E.V. On the blow-up of a solution of a non-local system of equations of hydrodynamic type. Izv. Math.; 2012; 76, pp. 190-213. [DOI: https://dx.doi.org/10.1070/IM2012v076n01ABEH002580]

41. Kolmogorov, A.N. , Fomin, S.V. Introductory Real Analysis; Dover Publications: New York, NY, USA, 1975.

42. Renardy, M.; Rogers, R. An Introduction to Partial Differential Equations; 2nd ed. Springer: New York, NY, USA, 2004.

43. Adams, R.A.; Fournier, J.J.F. Sobolev Spaces, Vol. 40 of Pure and Applied Mathematics; Elsevier: Amsterdam, The Netherlands, 2003.

44. Antontsev, S.N.; Khompysh, K. Inverse problems for a Boussinesq system for incompressible viscoelastic fluids. Math. Methods Appl. Sci.; 2023; 46, pp. 11130-11156. [DOI: https://dx.doi.org/10.1002/mma.9172]

45. Rajagopal, K.R. On some unresolved issues in non-linear fluid dynamics. Russ. Math. Surv.; 2003; 58, pp. 319-330. [DOI: https://dx.doi.org/10.1070/RM2003v058n02ABEH000612]

46. Fetecau, C.; Fetecau, C. The first problem of Stokes for an Oldroyd-B fluid. Int. J. Non-Linear Mech.; 2003; 38, pp. 1539-1544. [DOI: https://dx.doi.org/10.1016/S0020-7462(02)00117-8]

47. Pukhnachev, V.V. Mathematical model of an incompressible viscoelastic Maxwell medium. J. Appl. Mech. Tech. Phys.; 2010; 51, pp. 546-554. [DOI: https://dx.doi.org/10.1007/s10808-010-0071-5]

48. Baranovskii, E.S. On steady motion of viscoelastic fluid of Oldroyd type. Sb. Math.; 2014; 205, pp. 763-776. [DOI: https://dx.doi.org/10.1070/SM2014v205n06ABEH004397]

49. Artemov, M.A.; Baranovskii, E.S. Mixed boundary-value problems for motion equations of a viscoelastic medium. Electron. J. Differ. Equ.; 2015; 2015, 252.

50. Ershkov, S.V.; Christianto, V.; Shamin, R.V.; Giniyatullin, A.R. About analytical ansatz to the solving procedure for Kelvin–Kirchhoff equations. Eur. J. Mech. B Fluids; 2020; 79, pp. 87-91. [DOI: https://dx.doi.org/10.1016/j.euromechflu.2019.08.011]

51. Ye, Z. Global regularity of the high-dimensional Oldroyd-B model in the corotational case. J. Math. Anal. Appl.; 2020; 486, 123867. [DOI: https://dx.doi.org/10.1016/j.jmaa.2020.123867]

52. Fetecau, C.; Bridges, C. Analytical solutions for some unsteady flows of fluids with linear dependence of viscosity on the pressure. Inverse Probl. Sci. Eng.; 2021; 29, pp. 378-395. [DOI: https://dx.doi.org/10.1080/17415977.2020.1791109]

53. Mallea-Zepeda, E.; Ortega-Torres, E.; Villamizar-Roa, É.J. An optimal control problem for the Navier–Stokes-α system. J. Dyn. Control Syst.; 2023; 29, pp. 129-156. [DOI: https://dx.doi.org/10.1007/s10883-021-09566-x]

54. Sirota, D.D.; Gushchin, K.A.; Khan, S.A.; Kostikov, S.L.; Butov, K.A. Neural operators for hydrodynamic modeling of underground gas storage facilities. Russ. Tech. J.; 2024; 12, pp. 102-112. [DOI: https://dx.doi.org/10.32362/2500-316X-2024-12-6-102-112]

55. Mallea-Zepeda, E.; Lenes, E.; Valero, E. Boundary control problem for heat convection equations with slip boundary condition. Math. Probl. Eng.; 2018; 2018, 7959761. [DOI: https://dx.doi.org/10.1155/2018/7959761]

56. Brizitskii, R.V.; Saritskaya, Z.Y. Control problem for generalized Boussinesq model. Journal of Physics: Conference Series; IOP Publishing: Bristol, UK, 2019; Volume 1268, 012011. [DOI: https://dx.doi.org/10.1088/1742-6596/1268/1/012011]

57. Boldrini, J.L.; Mallea-Zepeda, E.; Rojas-Medar, M.A. Optimal boundary control for the stationary Boussinesq equations with variable density. Commun. Contemp. Math.; 2020; 22, 1950031. [DOI: https://dx.doi.org/10.1142/S0219199719500317]

58. Brizitskii, R.V.; Saritskaia, Z.Y. Analysis of inhomogeneous boundary value problems for generalized Boussinesq model of mass transfer. J. Dyn. Control Syst.; 2023; 29, pp. 1809-1828. [DOI: https://dx.doi.org/10.1007/s10883-023-09663-z]

59. Chebotarev, A.Y.; Pak, N.M.; Kovtanyuk, A.E. Analysis and numerical simulation of the initial-boundary value problem for quasilinear equations of complex heat transfer. J. Appl. Ind. Math.; 2023; 17, pp. 698-709. [DOI: https://dx.doi.org/10.1134/S1990478923040026]

60. Chebotarev, A.Y. Optimal control of quasi-stationary equations of complex heat transfer with reflection and refraction conditions. Comput. Math. Math. Phys.; 2023; 63, pp. 2050-2059. [DOI: https://dx.doi.org/10.1134/S0965542523110064]

61. Chebotarev, A.Y. Inhomogeneous problem for quasi-stationary equations of complex heat transfer with reflection and refraction conditions. Comput. Math. Math. Phys.; 2023; 63, pp. 441-449. [DOI: https://dx.doi.org/10.1134/S0965542523030053]

62. Baranovskii, E.S. Exact solutions for non-isothermal flows of second grade fluid between parallel plates. Nanomaterials; 2023; 13, 1409. [DOI: https://dx.doi.org/10.3390/nano13081409]

63. Brizitskii, R.V. Generalised Boussinesq model with variable coefficients. Sib. Electron. Math. Rep.; 2024; 21, pp. 213-227. [DOI: https://dx.doi.org/10.33048/semi.2024.21.015]

64. Baranovskii, E.S.; Brizitskii, R.V.; Saritskaia, Z.Y. Optimal control problems for the reaction–diffusion–convection equation with variable coefficients. Nonlinear Anal. Real World Appl.; 2024; 75, 103979. [DOI: https://dx.doi.org/10.1016/j.nonrwa.2023.103979]

65. Davidson, P.A. Introduction to Magnetohydrodynamics; 2nd ed. Cambridge University Press: Cambridge, UK, 2016; [DOI: https://dx.doi.org/10.1017/9781316672853]

66. Galtier, S. Introduction to Modern Magnetohydrodynamics; Cambridge University Press: Cambridge, UK, 2016; [DOI: https://dx.doi.org/10.1017/CBO9781316665961]

67. Xiros, N. Mathematical Magnetohydrodynamics; Springer: Cham, Switzerland, 2018; [DOI: https://dx.doi.org/10.1007/978-3-031-79601-2]

68. Brizitskii, R.V. Boundary value and control problems for the stationary magnetic hydrodynamic equations of heat conducting fluid with variable coefficients. J. Dyn. Control Syst.; 2024; 30, 39. [DOI: https://dx.doi.org/10.1007/s10883-024-09706-z]