Steady-state non-uniform Poiseuille shear flows with Navier boundary condition

Full Text

Introduction

The classical Poiseuille exact solution \(\boldsymbol{V}(z) = \bigl(V_x(z), 0, 0\bigr)\) describes a layered unidirectional gradient flow of a viscous incompressible fluid with a parabolic profile \(V_x = C_1 z^2 + C_2 z + C_3\), governed by three integration constants [1, 2, 4–7]. This exact solution serves as a basis for solving boundary value problems with diverse boundary conditions across various fluid types [6–13]. These boundary conditions are typically formulated either kinematically, in terms of velocities (the no-slip condition), or dynamically, in terms of stresses (the slip condition) [6–13].

Historically, the no-slip condition has been employed in the vast majority of hydrodynamic flow studies [5–7]. In recent years, however, the slip of a fluid at channel boundaries has regained significant attention in computational and theoretical analyses [8–13]. Poiseuille flow, subject to either no-slip or Navier slip conditions, has been examined not only for Newtonian fluids but also for fluids with shear-dependent viscosity and for complex media with internal microstructure, such as micropolar fluids [14–19]. Moreover, the exact Poiseuille solution plays a fundamental role in hydrodynamic stability theory [5, 8, 20–22]. The stability of this background flow is investigated for various classes of perturbations (or secondary flows) in both classical and non-classical, including rheological, media [5, 8, 20–22].

Therefore, finding exact solutions of the Navier–Stokes equations that describe gradient non-one-dimensional flows of incompressible fluids remains an important and current challenge in hydrodynamics. To develop new ansatzes for these equations, we observe that the exact Poiseuille solution satisfies an overdetermined system of partial differential equations [5–7, 23]. Specifically, it fulfills an additional constraint, which is the continuity equation [5–7, 14–18, 23].

A natural progression is to consider exact solutions of the form \(\boldsymbol{V}(z) = \bigl(V_x(z), V_y(z), 0\bigr)\). While this solution formally corresponds to a two-dimensional velocity field, a simple rotation transformation can reduce it to an equivalent unidirectional flow [5–7, 19, 23]. To capture nonlinear effects, one must construct exact Poiseuille-type solutions with a velocity vector of the form \(\boldsymbol{V}(x,y,z) = \bigl(V_x(x, y, z), V_y(x, y, z), 0\bigr)\). For such shear flows, the governing system reduces to an overdetermined set of equations [6, 7, 19, 24]. The first nontrivial exact solution to this system was derived within the framework now known as the Lin–Sidorov–Aristov class [24–26].

Nontrivial exact solutions for Couette-type flows under various boundary conditions have been constructed and analyzed for an infinite horizontal layer [6–10, 23, 24]. Similar analyses have been extended to two-layer fluid systems [27, 28]. Investigations into non-uniform Couette-type flows, however, have primarily been confined to the no-slip condition [6, 7, 19, 24, 27, 28]. For instance, the study in [23] explored the velocity field structure for gradient motions and flows down an inclined plane (Nusselt-type flows), yet exclusively for the no-slip case. It was demonstrated that isobaric flows exhibit a single stagnation point within the fluid layer, whereas gradient flows can develop up to four distinct flow zones, corresponding to as many as three stagnation points [29].

The primary motivation for the present work is to investigate non-uniform shear flows of the Poiseuille type within the Lin–Sidorov–Aristov class, with Navier slip conditions imposed on one boundary of an infinite horizontal fluid layer.

1. Problem statement

The steady isothermal flow of a viscous incompressible fluid is described by the Navier–Stokes equations:
\[\begin{equation} \tag{1}
(\boldsymbol{V}\cdot\boldsymbol{\nabla} ) \boldsymbol{V}=- \boldsymbol{\nabla} P +\nu\Delta \boldsymbol{V},
\end{equation}\]
\[\begin{equation}
\tag{2}
\boldsymbol{\nabla}\cdot\boldsymbol{V}=0.
\end{equation}\]
The following notation is used in (1), (2): $\boldsymbol{V}(x,y,z)= (V_x,V_y,V_z)$ is the velocity vector; $P$ is the pressure deviation from the hydrostatic level, normalized by the mean density; $\nu$ is the kinematic viscosity of the fluid; $\boldsymbol{\nabla}=\bigl\{\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z}\bigr\}$ is the nabla operator; $\Delta=\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$ is the Laplace operator.

Projecting the invariant form of equations (1), (2 onto the axes of a Cartesian coordinate system yields the following system of nonlinear partial differential equations:
\[\begin{equation}
\tag{3}
\begin{array}{l}
V_x\dfrac{\partial V_x}{\partial x} + V_y\dfrac{\partial V_x}{\partial y} + V_z\dfrac{\partial V_x}{\partial z} =
-\dfrac{\partial P}{\partial x}+\nu\Bigl(\dfrac{\partial^2 V_x}{\partial x^2} + \dfrac{\partial^2 V_x}{\partial y^2} + \dfrac{\partial^2 V_x}{\partial z^2}\Bigr),
\\
V_x\dfrac{\partial V_y}{\partial x} + V_y\dfrac{\partial V_y}{\partial y} + V_z\dfrac{\partial V_y}{\partial z} = 
-\dfrac{\partial P}{\partial y}+\nu\Bigl(\dfrac{\partial^2 V_y}{\partial x^2} + \dfrac{\partial^2 V_y}{\partial y^2} + \dfrac{\partial^2 V_y}{\partial z^2}\Bigr),
\\
V_x\dfrac{\partial V_z}{\partial x} + V_y\dfrac{\partial V_z}{\partial y} + V_z\dfrac{\partial V_z}{\partial z} = 
-\dfrac{\partial P}{\partial z}+\nu\Bigl(\dfrac{\partial^2 V_z}{\partial x^2} + \dfrac{\partial^2 V_z}{\partial y^2} + \dfrac{\partial^2 V_z}{\partial z^2}\Bigr),
\\
\dfrac{\partial V_x}{\partial x} + \dfrac{\partial V_y}{\partial y} + \dfrac{\partial V_z}{\partial z} = 0.
\end{array}
\end{equation}\]
We seek an exact solution to system (3) in the form [6, 7, 19, 23, 24, 27, 28]:
\[\begin{equation}
\tag{4}
\begin{array}{l}
V_x (x,y,z )=U(z)+xu_1(z)+yu_2(z),
\\
V_y(x,y,z)=V(z)+xv_1(z)+yv_2(z),
\\
V_z(x,y,z)=w(z).
\end{array}
\end{equation}\]
The pressure is assumed to be a linear function of the spatial coordinates [6, 7, 19, 23, 24, 27, 28]:
\[\begin{equation}
\tag{5}
P(x,y,z)=P_0(z)+P_1(z)x+P_2(z)y.
\end{equation}\]
Substituting the ansatz (4), (5) into (3) and projecting the results onto the $Ox$, $Oy$, and $Oz$ axes gives:
\[\begin{multline*}
(U+u_1x+u_2y)u_1 + (V+v_1x+v_2y )u_2 + w\Bigl(\frac{\partial U}{\partial z}+\frac{\partial u_1}{\partial z}x+\frac{\partial u_2}{\partial z}y\Bigr)={}
\\
{}=-P_1+\nu\Bigl(\frac{\partial^2 U}{\partial z^2}+\frac{\partial^2 u_1}{\partial z^2}x+\frac{\partial^2 u_2}{\partial z^2}y\Bigr),
\end{multline*}\]
\[\begin{multline}
\tag{6}
(U+u_1x+u_2y)v_1 + (V+v_1x+v_2y)v_2 + w\Bigl(\frac{\partial V}{\partial z}+\frac{\partial v_1}{\partial z}x+\frac{\partial v_2}{\partial z}y\Bigr)={}
\\
{} =-P_2+\nu\Bigl(\frac{\partial^2 V}{\partial z^2}+\frac{\partial^2 v_1}{\partial z^2}x+\frac{\partial^2 v_2}{\partial z^2}y\Bigr),
\end{multline}\]
\[
w\frac{\partial w}{\partial z}=-\Bigl(\frac{\partial P_0}{\partial z}+\frac{\partial P_1}{\partial z}x+\frac{\partial P_2}{\partial z}y\Bigr)
+\nu\frac{\partial^2 w}{\partial z^2},
\quad 
u_1+v_2+\frac{\partial w}{\partial z}=0.
\]
Both sides of each equation in (6) are linear in the coordinates $x$ and $y$. Equating the coefficients of like terms yields a system of ten ordinary nonlinear differential equations for the ten unknown functions:
\[\begin{equation}
\tag{7}
\begin{array}{c}
Uu_1+Vu_2+w\dfrac{\partial U}{\partial z}=-P_1+\nu\dfrac{\partial^2 U}{\partial z^2},
\\
u_1^2+v_1u_1+w\dfrac{\partial u_1}{\partial z}=\nu\dfrac{\partial^2 u_1}{\partial z^2},
\quad
u_1u_2+v_2u_2+w\dfrac{\partial u_2}{\partial z}=\nu\dfrac{\partial^2 u_2}{\partial z^2},
\\
Uv_1+Vv_2+w\dfrac{\partial V}{\partial z}=-P_2+\nu\dfrac{\partial^2 V}{\partial z^2},
\\
u_1v_1+v_1v_2+w\dfrac{\partial v_1}{\partial z}=\nu\dfrac{\partial^2 v_1}{\partial z^2},
\quad
u_2v_1+v_2^2+w\dfrac{\partial v_2}{\partial z}=\nu\dfrac{\partial^2 v_2}{\partial z^2},
\\
w\dfrac{\partial w}{\partial z}=-\dfrac{\partial P_0}{\partial z}+\nu\dfrac{\partial^2 w}{\partial z^2},
\quad
\dfrac{\partial P_1}{\partial z}=0,
\quad
\dfrac{\partial P_2}{\partial z}=0,
\\
u_1+v_2+\dfrac{\partial w}{\partial z}=0.
\end{array} 
\end{equation}\]
All derivatives in (7) are taken with respect to $z$. Denoting $d/dz$ by a prime, the system becomes
\[\begin{equation}
\tag{8}
\begin{array}{c}
\nu u_1^{\prime\prime}-wu_1^{\prime}-u_1^2-v_1u_1=0,
\quad
\nu u_2^{\prime\prime}-wu_2^{\prime}-u_1u_2-v_2u_2=0,
\\
\nu v_1^{\prime\prime}-wv_1^{\prime}-u_1v_1-v_1v_2=0,
\quad
\nu v_2^{\prime\prime}-wv_2^{\prime}-u_2v_1-v_2^2=0,
\\
w^{\prime}+u_1+v_2=0,
\\
\nu U^{\prime\prime}-wU^{\prime}-Uu_1-Vu_2=P_1,
\quad
\nu V^{\prime\prime}-wV^{\prime}-Uv_1-Vv_2=P_2,
\\
P_0^{\prime}=\nu w^{\prime\prime}-ww^{\prime},
\quad
P_1^{\prime}=0,
\quad
P_2^{\prime}=0.
\end{array}
\end{equation}\]
We now non-dimensionalize the system. The horizontal coordinates $x$ and $y$ are scaled by a length $l$, and the vertical coordinate $z$ by the layer thickness $h$:
\[
X= {x}/{l},
\quad
Y={y}/{l},
\quad
Z={z}/{h}.
\]
Two characteristic lengths are introduced: $h$ for the thin-layer (shear-flow) approximation, and $l$ to model flows in laterally bounded domains where a permeability condition holds on the sidewalls.

Correspondingly, two characteristic velocities are defined: $[V]$ for the longitudinal motion and $[w]$ for the transverse motion:
\[\begin{gather*}
\overline{V}_x= {V_x}/{[V]},
\quad
\overline{V}_y= {V_y}/{[V]},
\quad
\overline{V}_z= {V_z}/{[w]}.
\end{gather*}\]
Applying this scaling to the representation (4) gives
\[
\overline{V}_x=\frac{V_x}{[V]}=\frac{U+Xl\cdot u_1+Yl\cdot u_2}{[V]}=\frac{U}{[V]}+X\frac{u_1l}{[V]}+Y\frac{u_2l}{[V]}.
\]
Hence,
\[
\overline{U}= {U}/{[V]},
\quad
\overline{u}_1= {u_1l}/{[V]},
\quad
\overline{u}_2= {u_2l}/{[V]}.
\]
Similarly,
\[
\overline{V}= {V}/{[V]},
\quad
\overline{v}_1= {v_1l}/{[V]},
\quad
\overline{v}_2= {v_2l}/{[V]}.
\]
A characteristic pressure scale $[P_0]$ is introduced as
\[
\overline{P}_0= {P_0}/{[P_0]},
\quad
\overline{P}_1= {P_1l}/{[P_0]},
\quad
\overline{P}_2= {P_2l}/{[P_0]}.
\] 
We proceed to non-dimensionalize system (8). Starting from the continuity equation $w^{\prime}+u_1+v_2=0$,
\[\begin{gather*}
\frac{dw}{dz}+u_1+v_2=0,
\quad
\frac{d\overline{w}[w]}{d (Z\cdot h )}+\overline{u}_1\frac{[V]}{l}+\overline{v}_2\frac{[V]}{l}=0,
\\
\overline{w}^{\prime}\frac{[w]}{h}+\left(\overline{u}_1+\overline{v}_2\right)\frac{[V]}{l}=0,
\quad
\overline{w}^{\prime}+ (\overline{u}_1+\overline{v}_2 )\frac{[V]h}{[w]l}=0,
\\
\overline{w}^{\prime}+ (\overline{u}_1+\overline{v}_2 )\Bigl(\frac{h}{l}\Bigr)^2\cdot\frac{[V]l}{\nu}\cdot\frac{\nu}{[w]h}=0,
\quad 
\overline{w}^{\prime}+ (\overline{u}_1+\overline{v}_2 )\delta^2\cdot\frac{\textsf{Re}_V}{\textsf{Re}_w}=0.
\end{gather*}\]
Here the prime denotes $d/dZ$, $\textsf{Re}_V= {[V]l}/{\nu}$, $\textsf{Re}_w= {[w]h}/{\nu}$, $\delta= {h}/{l}$ are the Reynolds numbers (the subscript refers to the reference velocity) and the geometric anisotropy parameter.

Applying the same scaling to the remaining equations of (8) yields the dimensionless system
\[\begin{equation}
\tag{9}
\begin{array}{c}
\overline{u}_1^{\prime\prime}-\mathsf{Re}_w\overline{w}\,\overline{u}_1^{\prime}-\delta^2\mathsf{Re}_V (\overline{u}_1^2+\overline{v}_1\overline{u}_1 )=0,
\\
\overline{u}_2^{\prime\prime}- \mathsf{Re}_w\overline{w}\,\overline{u}_2^{\prime}-\delta^2 \mathsf{Re}_V  (\overline{u}_1\overline{u}_2+\overline{v}_2\overline{u}_2 )=0,
\\
\overline{v}_1^{\prime\prime}- \mathsf{Re}_w\overline{w}\,\overline{v}_1^{\prime}-\delta^2 \mathsf{Re}_V  (\overline{u}_1\overline{v}_1+\overline{v}_1\overline{v}_2 )=0,
\\
\overline{v}_2^{\prime\prime}-\mathsf{Re}_w\overline{w}\,\overline{v}_2^{\prime}-\delta^2\mathsf{Re}_V (\overline{u}_2\overline{v}_1+\overline{v}_2^2 )=0,
\\
\overline{w}^{\prime}+ (\overline{u}_1+\overline{v}_2 )\delta^2\dfrac{\mathsf{Re}_V}{\mathsf{Re}_w}=0,
\\
\overline{U}^{\prime\prime}-\mathsf{Re}_w\overline{w}\,\overline{U}^{\prime}-\delta^2 \mathsf{Re}_V  (\overline{U}\,\overline{u}_1+\overline{V}\,\overline{u}_2 )=\delta^2\mathsf{Fr}_V\overline{P}_1,
\\
\overline{V}^{\prime\prime}-\mathsf{Re}_w\overline{w}\,\overline{V}^{\prime}-\delta^2 \mathsf{Re}_V  (\overline{U}\,\overline{v}_1+\overline{V}\,\overline{v}_2 )=\delta^2 \mathsf{Fr}_V\overline{P}_2,
\\
\overline{P}_0^{\prime}=\dfrac{1}{\mathsf{Fr}_w} (\overline{w}^{\prime\prime}-\mathsf{Re}_w\overline{w}\,\overline{w}^{\prime} ),
\quad
\overline{P}_1^{\prime}=0,
\quad
\overline{P}_2^{\prime}=0.
\end{array}
\end{equation}\]
Here $\mathsf{Fr}_V= {[P_0]l}/({\nu [V]})$ and $\mathsf{Fr}_w= {[P_0]l}/({\nu [w]})$ are the Froude numbers based on the velocities $V$ and $w$, respectively.

The equations $\overline{P}_1^{\prime}=0$ and $\overline{P}_2^{\prime}=0$ show that the pressure gradients are constant. The background pressure $\overline{P}_0$ is obtained from the vertical velocity $\overline{w}$. In the subsequent analysis we restrict ourselves to shear flows by setting $\overline{w}=0$.

2. Shear flows

For shear flows, the system (9) simplifies considerably:
\[\begin{equation*}
\begin{array}{c}
\overline{u}_1^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{u}_1^2 + \overline{v}_1\overline{u}_1 ) = 0,
\quad   
\overline{u}_2^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{u}_1\overline{u}_2 + \overline{v}_2\overline{u}_2 ) = 0,
\\
\overline{v}_1^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{u}_1\overline{v}_1 + \overline{v}_1\overline{v}_2 ) = 0, 
\quad
\overline{v}_2^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{u}_2\overline{v}_1 + \overline{v}_2^2 ) = 0,
\\
\overline{u}_1 + \overline{v}_2 = 0,
\\
\overline{U}^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{U}\,\overline{u}_1 + \overline{V}\,\overline{u}_2 ) = \delta^2 \mathsf{Fr}_V \overline{P}_1,
\\
\overline{V}^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{U}\,\overline{v}_1 + \overline{V}\,\overline{v}_2 ) = \delta^2 \mathsf{Fr}_V \overline{P}_2 .
\end{array}
\end{equation*}\]

With $\overline{w}=0$, the pressure field in representation (5) reduces to a linear function with constant coefficients:
\[\begin{equation}
\tag{10}
\overline{P} = \overline{P}_0 + \overline{P}_1 X + \overline{P}_2 Y, \qquad \overline{P}_i = \text{const},
\end{equation}\]
where the constants may also be determined from the boundary conditions.

Note further that, owing to the continuity condition $\overline{u}_1 + \overline{v}_2 = 0$, the nonlinear terms in the equations for $\overline{u}_2$ and $\overline{v}_1$ vanish, leaving
\[
\overline{u}_2^{\prime\prime}=0, \qquad \overline{v}_1^{\prime\prime}=0.
\]
These equations integrate to linear functions of the dimensionless vertical coordinate $Z$.

Thus, the governing system for shear flows becomes
\[\begin{equation}
\tag{11}
\begin{array}{c}
\overline{u}_1^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{u}_1^2 + \overline{v}_1\overline{u}_1 ) = 0,
\quad 
\overline{u}_2^{\prime\prime} = 0,\quad \overline{v}_1^{\prime\prime} = 0,
\\
\overline{v}_2^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{u}_2\overline{v}_1 + \overline{v}_2^2 ) = 0,
\quad \overline{u}_1 + \overline{v}_2 = 0,
\\
\overline{U}^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{U}\,\overline{u}_1 + \overline{V}\,\overline{u}_2 ) = \delta^2 \mathsf{Fr}_V \overline{P}_1,
\\
\overline{V}^{\prime\prime} - \delta^2 \mathsf{Re}_V  (\overline{U}\,\overline{v}_1 + \overline{V}\,\overline{v}_2 )  = \delta^2 \mathsf{Fr}_V \overline{P}_2 .
\end{array}
\end{equation}\]

3. Generating solution

Consider a generating solution of the form [6, 7, 19, 24, 27, 29]:
\[\begin{equation} 
\tag{12}
\overline{V}_x = \overline{U} + Y\overline{u}_2, \qquad \overline{V}_y = \overline{V}.
\end{equation}\]
This solution is a special case of representation (4) with $\overline{u}_1 = \overline{v}_1 = \overline{v}_2 = 0$. The choice of class (12) is motivated by two factors. First, the incompressibility condition is identically satisfied, thereby eliminating the issue of overdeterminacy in the governing system.

Second, by rotating the coordinates in the horizontal $XY$-plane, one can transform class (12) into the general form (4). This is achieved by the simultaneous replacements [6, 7, 19, 29]:
\[
X \; \longrightarrow \; X\cos\psi + Y\sin\psi \equiv \widetilde{X},
\qquad
Y \; \longrightarrow \; -X\sin\psi + Y\cos\psi \equiv \widetilde{Y},
\]
and
\[
\overline{V}_x \; \longrightarrow \; \overline{V}_x\cos\psi + \overline{V}_y\sin\psi \equiv \widetilde{V}_x,
\qquad
\overline{V}_y \; \longrightarrow \; -\overline{V}_x\sin\psi + \overline{V}_y\cos\psi \equiv \widetilde{V}_y .
\]

Indeed, using the coordinate transformation and solving for the old coordinates gives
\[
\begin{pmatrix} X \\ Y \end{pmatrix}
= \begin{pmatrix} \cos\psi & -\sin\psi \\ \sin\psi & \cos\psi \end{pmatrix}
  \begin{pmatrix} \widetilde{X} \\ \widetilde{Y} \end{pmatrix},
\]  
so that
\[
X = \widetilde{X}\cos\psi - \widetilde{Y}\sin\psi,\quad
Y = \widetilde{X}\sin\psi + \widetilde{Y}\cos\psi .
\]

Substituting these relations into the velocity transformation yields
\[\begin{multline*} 
\widetilde{V}_x =  (\overline{U} + \overline{u}_2 Y )\cos\psi + \overline{V}\sin\psi  = \\ =
 \bigl[\overline{U} + \overline{u}_2 (\widetilde{X}\sin\psi + \widetilde{Y}\cos\psi )\bigr]\cos\psi + \overline{V}\sin\psi = 
\\
=  (\overline{U}\cos\psi + \overline{V}\sin\psi )
   +  (\overline{u}_2\sin\psi\cos\psi )\widetilde{X}
   +  (\overline{u}_2\cos^{2}\psi )\widetilde{Y} \equiv 
\widetilde{U} + \widetilde{u}_1\widetilde{X} + \widetilde{u}_2\widetilde{Y},
\end{multline*}\]
and analogously,
\[\begin{multline*} 
\widetilde{V}_y = - (\overline{U} + \overline{u}_2 Y )\sin\psi + \overline{V}\cos\psi = \\ =
-\bigl[\overline{U} + \overline{u}_2 (\widetilde{X}\sin\psi + \widetilde{Y}\cos\psi )\bigr]\sin\psi + \overline{V}\cos\psi = 
\\
= (-\overline{U}\sin\psi + \overline{V}\cos\psi)
   + (-\overline{u}_2\sin^{2}\psi )\widetilde{X}
   + (-\overline{u}_2\sin\psi\cos\psi )\widetilde{Y}  \equiv 
   \widetilde{V} + \widetilde{v}_1\widetilde{X} + \widetilde{v}_2\widetilde{Y}.
\end{multline*}\]

Returning to system (11) and specializing it to class (12) gives
\[\begin{equation} \tag{13}
\overline{u}_2^{\prime\prime}=0,\qquad
\overline{V}^{\prime\prime} - \delta^{2}\mathsf{Re}_V\,\overline{V}\,\overline{u}_2 = \delta^{2}\mathsf{Fr}_V\overline{P}_1,\qquad
\overline{V}^{\prime\prime} = \delta^{2}\mathsf{Fr}_V\overline{P}_2 .
\end{equation}\]

The general solution of (13), together with (10), can be written as
\[\begin{equation} \tag{14}
\begin{array}{c}
\overline{P}_0 = C_1,\quad 
\overline{P}_1 = C_2,\quad 
\overline{P}_2 = C_3,
\\
\overline{u}_2 = C_4 Z + C_5,\quad
\overline{V} = \dfrac{\delta^{2}}{2}\,\mathsf{Fr}_V C_3 Z^{2} + C_6 Z + C_7,
\\
\overline{U} = \dfrac{\mathsf{Fr}_V\delta^{2}}{2}
\Bigl[C_2 + \dfrac{1}{60}\,\mathsf{Re}_V\delta^{2}C_3
                (3C_4 Z + 5C_5 )Z^{2}\Bigr]Z^{2} +{}
\hspace{2cm} ~
\\
+ \dfrac{1}{6}\,\mathsf{Re}_V\delta^{2}
                \Bigl[\dfrac{1}{2}Z  (C_4 Z + 2C_5 )C_6
                +  (C_4 Z + 3C_5 )C_7\Bigr]Z^{2}
                + C_8 Z + C_9 .
\end{array}
\end{equation}\]

To obtain a particular solution corresponding to specific flow conditions, one must impose a set of boundary conditions that determine the integration constants $C_1$, $C_2$, $\dots$, $C_9$ in (14).

4. Boundary conditions

Consider a flow in a horizontal layer of thickness $h$. We assume that the Navier slip condition is imposed at the lower boundary $z = 0$ [30]:
\[
\Bigl[
    \alpha \frac{\partial V_x}{\partial z} - V_x
\Bigr]_{z=0}=0, \qquad
\Bigl[
    \alpha \frac{\partial V_y}{\partial z} - V_y
\Bigr]_{z=0}=0.
\]

The parameter $\alpha$ is the slip length. Rewriting these conditions in terms of the normalized variables gives
\[\begin{equation}
\tag{15}
\begin{array}{ll}
\Bigl[
    \alpha \dfrac{\partial  (\overline{V}_x\cdot [V] )}{\partial  (Z\cdot h )} - \overline{V}_x\cdot [V]
\Bigr]_{z=0}=0, &
\Bigl[
    \alpha \dfrac{\partial (\overline{V}_y\cdot [V] )}{\partial (Z\cdot h)} - \overline{V}_y\cdot [V]
\Bigr]_{z=0}=0;
\\
\Bigl[
    \overline{\alpha} \dfrac{\partial \overline{V}_x}{\partial Z} - \overline{V}_x
\Bigr]_{Z=0}=0,  &
\Bigl[
    \overline{\alpha} \dfrac{\partial \overline{V}_y}{\partial Z} - \overline{V}_y
\Bigr]_{z=0}=0.
\end{array}
\end{equation}\]
In (15) we have introduced the dimensionless slip length $\overline{\alpha}={\alpha}/{h}$.

At the upper boundary of the layer $z = h$ (i.e., $Z=1$) we prescribe the pressure distribution
\[
P|_{z=h}=S_0+S_1x+S_2y.
\]
In the dimensionless variables this becomes
\[\begin{gather} 
\tag{16}
\overline{P}|_{Z=1}=\overline{S}_0+\overline{S}_1 X+\overline{S}_2 Y,
\end{gather}\]
where $\overline{S}_0=S_0/[P_0]$, $\overline{S}_1=S_1l/[P_0]$, and $\overline{S}_2=S_2l/[P_0]$.

We also specify the velocity distribution at the upper boundary. For the chosen solution structure ($u_1=v_1=v_2=0$) it reduces to
\[\begin{gather*}
V_x|_{z=h}=[V]\cos\varphi+\Omega y, \quad V_y|_{z=h}=[V]\sin\varphi.
\end{gather*}\]
Thus, the characteristic horizontal velocity $[V]$ is identified with the uniform component of the velocity measured at one boundary, while $\Omega$ quantifies the spatial gradient of $V_x$ (i.e., the flow inhomogeneity).

In dimensionless form these conditions are
\[\begin{gather} 
\tag{17}
\overline{V}_x|_{Z=1}=\cos\varphi+\frac{\Omega l}{[V]} Y, \quad 
\overline{V}_y|_{Z=1}=\sin\varphi.
\end{gather}\]
The dimensionless group ${\Omega l}/{[V]}$ can be transformed as follows:
\[\begin{gather*}
\frac{\Omega l}{[V]}=\frac{\nu}{[V]l}\cdot\frac{2\Omega l^2}{\nu}\cdot\frac{1}{2}=\frac{\mathsf{Ta}}{2\mathsf{Re}_V},
\end{gather*}\]
where $\mathsf{Ta} = 2\Omega l^2/\nu$ is a modified Taylor number¹ that reflects the ratio of centrifugal to viscous forces. This modification stems from the fact that the spatial acceleration acts as a local angular velocity of the fluid [6, 9, 32]. Using this notation, condition (17) can be recast as
\[\begin{gather} 
\tag{18}
\overline{V}_x|_{Z=1}=\cos\varphi+\frac{\mathsf{Ta}}{2\mathsf{Re}_V} Y, \quad 
\overline{V}_y|_{Z=1}=\sin\varphi.
\end{gather}\]

The solution to the boundary‑value problem (14)–(16), (18) reads
\[\begin{gather*}
\overline{P}_0=\overline{S}_0,
\quad
\overline{P}_1=\overline{S}_1,
\quad
\overline{P}_2=\overline{S}_2,
\\
\overline{u}_2=\frac{\mathsf{Ta} ( Z+\overline{\alpha } )}{2\mathsf{Re}_V ( 1+\overline{\alpha } )},
\quad
\overline{V}=\frac{\delta^{2}}{2}\mathsf{Fr}_V\overline{S}_2 Z^{2}
            +\frac{ ( Z+\overline{\alpha } )
                   ( 2\sin\varphi -\delta^{2}\mathsf{Fr}_V\overline{S}_2 )}
                   {2 ( 1+\overline{\alpha }  )},
\end{gather*}\]

\[\begin{multline}
\tag{19}
\overline{U}= \frac{\delta^{4}\mathsf{Fr}_V\overline{S}_2\mathsf{Ta}}
                   {80 ( 1+\overline{\alpha }  )}\,Z^{5}
            +\bigl[ \delta^{2}\mathsf{Fr}_V\overline{S}_2
                    ( \overline{\alpha}^{2}+\overline{\alpha}-1  )
                   +2\sin\varphi \bigr]
                   \frac{\delta^{2}\mathsf{Ta}\,Z^{4}}
                        {48 ( 1+\overline{\alpha }  )^{2}} -
\\
            -\frac{\overline{\alpha}\delta^{2}\mathsf{Ta}
                    ( \delta^{2}\mathsf{Fr}_V\overline{S}_2-2\sin\varphi )}
                   {12 ( 1+\overline{\alpha } )^{2}}\,Z^{3} +
\\
            +\bigl[ 4 ( 1+\overline{\alpha } )^{2}\mathsf{Fr}_V\overline{S}_1
                   -\overline{\alpha}^{2}\mathsf{Ta}
                    ( \delta^{2}\mathsf{Fr}_V\overline{S}_2-2\sin\varphi ) \bigr]
                   \frac{\delta^{2}Z^{2}}
                        {8 ( 1+\overline{\alpha } )^{2}}+
\\
            +\Bigl[
                   \frac{\delta^{2}\mathsf{Ta} (
                         \delta^{2}\mathsf{Fr}_V\overline{S}_2
                         ( 2+12\overline{\alpha}+25\overline{\alpha}^{2} )
                         -10 ( 1+4\overline{\alpha}+6\overline{\alpha}^{2} )
                         \sin\varphi )}
                         {240 ( 1+\overline{\alpha }  )^{3}} +
\\
            + \frac{2\cos\varphi-\delta^{2}\mathsf{Fr}_V\overline{S}_1}
                        {2 ( 1+\overline{\alpha } )}
            \Bigr] ( Z+\overline{\alpha } ).
\end{multline}\]

Solution (19) describes three distinct regimes at the lower wall: complete adhesion ($\overline{\alpha}=0$), finite slip ($\overline{\alpha}$ a finite positive number), and perfect slip ($\overline{\alpha}\to\infty$). The case $\overline{\alpha}=0$ has been analyzed in detail in [24]; we therefore concentrate on the remaining two cases.

5. Perfect Slip

Consider the limiting case of perfect slip, which corresponds to the limit $\overline{\alpha} \to \infty$ in solution (19). The limit is readily evaluated because it involves only rational functions of $\overline{\alpha}$:
\[\begin{gather*}
    \overline{u_2}=\frac{\mathsf{Ta}}{2\mathsf{Re}_V},
\\
    \overline{V}=\frac{\delta^{2}}{2}\mathsf{Fr}_V\overline{S}_2 Z^{2}
                +\frac{ ( 2\sin\varphi -\delta^{2}\mathsf{Fr}_V\overline{S}_2 )}{2}
                =\frac{\delta^{2}}{2}\mathsf{Fr}_V\overline{S}_2 ( Z^{2}-1 )+\sin\varphi,
\\
    \overline{U}= \frac{\mathsf{Fr}_V\overline{S}_2 \delta^{4}\mathsf{Ta}\,Z^{4}}{48}
                +\Bigl[ \mathsf{Fr}_V\overline{S}_1
                      -\frac{\mathsf{Ta} ( \delta^{2}\mathsf{Fr}_V\overline{S}_2 -2\sin\varphi )}{4}
                      \Bigr]\frac{\delta^{2}Z^{2}}{2}+~~~~~
\\
 ~\hspace{4cm} +\frac{\delta^{2}\mathsf{Ta} ( 5\delta^{2}\mathsf{Fr}_V\overline{S}_2-12\sin\varphi )}{48}
                +\frac{2\cos\varphi -\delta^{2}\mathsf{Fr}_V\overline{S}_1}{2}.
\end{gather*}\]
The pressure components remain unchanged: $\overline{P}_0=\overline{S}_0$, $\overline{P}_1=\overline{S}_1$, $\overline{P}_2=\overline{S}_2$.

We begin the analysis with the velocity $\overline{V}$. Rewriting it as
\[\begin{gather*}
\overline{V}=\frac{\delta^{2}}{2}\mathsf{Fr}_V\overline{S}_2 Z^{2}
            +\Bigl( \sin\varphi-\frac{\delta^{2}}{2}\mathsf{Fr}_V\overline{S}_2 \Bigr)
\end{gather*}\]
makes its behaviour more transparent. In the flow domain this function is strictly monotonic (we discard the trivial case $\mathsf{Fr}_V\overline{S}_2=0$, where the leading coefficient vanishes). Consequently, it can possess at most one zero inside the layer. Such a zero exists if
\[\begin{gather*}
\sin\varphi\Bigl( \sin\varphi-\frac{\delta^{2}}{2}\mathsf{Fr}_V\overline{S}_2 \Bigr)<0,
\end{gather*}\]
which simply means that $\overline{V}$ takes opposite signs at the two boundaries.

Next we analyse the velocity $\overline{U}$. It is a linear combination of three terms: two strictly monotonic functions (a quartic and a quadratic) and a constant. The coefficients of this combination offer considerable freedom. Hence, referring to [31], we conclude that $\overline{U}$ can have at most two zeros inside the layer (Fig. 1)—a result consistent with the fact that a biquadratic equation admits at most two positive roots.

We now ask whether a zero of $\overline{V}$ can simultaneously be a zero of $\overline{U}$. Let $Z_*$ denote a point where $\overline{V}(Z_*)=0$; then
\[
\overline{V} (Z_* )=
\frac{\delta^{2}}{2}\mathsf{Fr}_V\overline{S}_2 (Z_*^{2}-1 )+\sin\varphi=0,
\quad
Z_*^{2}=1-\frac{2\sin\varphi}{\delta^{2}\mathsf{Fr}_V\overline{S}_2}.
\]
For $Z_*$ to lie inside the layer we require
\[\begin{gather*}
0<Z_*^{2}<1,
\quad\text{or equivalently}\quad
0<1-\frac{2\sin\varphi}{\delta^{2}\mathsf{Fr}_V\overline{S}_2}<1,
\end{gather*}\]
which in particular implies
\[
0<\frac{\sin\varphi}{\mathsf{Fr}_V\overline{S}_2}<\frac{\delta^{2}}2 .
\]

Substituting $Z_*$ into the expression for $\overline{U}$ and simplifying yields
\[
\overline{U}(Z_*)=
\frac{12\mathsf{Fr}_V\overline{S}_2\cos\varphi
     -12\mathsf{Fr}_V\overline{S}_1\sin\varphi
     +2\delta^{2}\mathsf{Fr}_V\overline{S}_2\mathsf{Ta}\sin\varphi
     -5\mathsf{Ta}\sin^{2}\varphi}
     {12\mathsf{Fr}_V\overline{S}_2} .
\]
Notice that when $\sin\varphi=0$ (i.e., when the flow at the upper boundary is unidirectional), $\overline{U}(Z_*)=\pm1$; therefore, in that case $\overline{U}$ and $\overline{V}$ share no common zero inside the layer.

 

Figure 1. Profile of the velocity $\overline{U}$ in the perfect slip case

 

6. Slippage at the Lower Boundary

We now analyse the general case of slip at the lower boundary, where $\overline{\alpha}$ is a finite positive number. For this purpose we examine solution (19) in more detail.

We begin with the spatial acceleration $\overline{u}_2$:
\[
\overline{u}_2 =\frac{\mathsf{Ta} (Z+\overline{\alpha} )}{2\mathsf{Re}_V (1+\overline{\alpha} )}.
\]
As in the perfect slip limit, $\overline{u}_2$ has no zeros inside the layer; however, its magnitude now varies with the vertical coordinate $Z$. This component of the velocity field keeps the same sign throughout the layer, the sign being determined by the ratio of the dimensionless numbers $\mathsf{Ta}$ and $\mathsf{Re}_V$.

Next we consider the background velocity $\overline{V}$:
\[
\overline{V}=\frac{\delta^2}{2}\mathsf{Fr}_V\overline{S}_2 Z^2+
\frac{ (Z+\overline{\alpha} ) (2\sin\varphi-\delta^2\mathsf{Fr}_V\overline{S}_2 )}{2 (1+\overline{\alpha})}.
\]
Clearly, this is a quadratic polynomial; therefore, when the leading coefficient is non‑zero ($\mathsf{Fr}_V\overline{S}_2\neq0$), it can possess at most two zeros.

If $\mathsf{Fr}_V\overline{S}_2=0$, there are no zeros at all. In that case $\overline{V}$ is strictly monotonic and retains the sign of the difference $(2\sin\varphi-\delta^2\mathsf{Fr}_V\overline{S}_2 )$. 

Now assume $\mathsf{Fr}_V\overline{S}_2\neq0$. The velocity $\overline{V}$ vanishes only when the equation
\[
Z^2=k\left(Z+\overline{\alpha}\right),
\]
with
\[
k=-\frac{ (2\sin\varphi-\delta^2\mathsf{Fr}_V\overline{S}_2 )}
       {\delta^2\mathsf{Fr}_V\overline{S}_2 (1+\overline{\alpha} )},
\]
admits a root in the interval $(0, 1)$. The definition of $k$ is meaningful only for $k>0$ (i.e., when $2\sin\varphi\big/ [\delta^2\mathsf{Fr}_V\overline{S}_2 (1+\overline{\alpha})]<1$).

Both sides of the equation are strictly increasing functions. Consequently, they can intersect at most once. Hence, $\overline{V}$ can have at most one zero inside the fluid layer.

For the velocity $\overline{U}$ it is not feasible to obtain similarly simple analytical bounds on the number of zeros, because the polynomial that describes it is of high order and its coefficients are not homogeneous. Invoking the fundamental theorem of algebra together with [31], we can only state that the number of zeros of $\overline{U}$ cannot exceed five.

In practice, the attainable number of zeros is somewhat lower. The six coefficients of the polynomial $\overline{U}$ are determined by seven independent parameters. Some of these coefficients are strongly correlated (for instance, the coefficients of the linear and constant terms), while some of the parameters are restricted to specific ranges (e.g., $\overline{\alpha}>0$, ${\delta>0}$) or contribute only weakly to the overall shape of the profile (e.g., when $\delta\ll 1$ or ${|\sin\varphi|\leqslant 1}$). A qualitative example of a velocity profile possessing three internal zeros is displayed in Fig. 2}.

 

Figure 2. Profile of the velocity $\overline{U}$ in the slip case

 

Conclusion

This paper presents an exact solution of an overdetermined Navier–Stokes system describing steady non‑uniform shear flows of a vertically swirling fluid. The non‑uniform Poiseuille flow of “two‑and‑a‑half” dimensions is analysed for a boundary condition that incorporates Navier slip and, as a limiting case, perfect slip at the substrate (the lower boundary of an infinite horizontal layer). To integrate the overdetermined system, a generating solution belonging to the Lin–Sidorov–Aristov class is employed. A polynomial exact solution of the boundary‑value problem is derived. An analysis of the polynomial roots reveals the existence of three stagnation points, which stratifies the velocity field into four zones with counter‑directed flows.

Competing interests. The authors declare no competing interests.
Authors' contributions. All authors contributed equally to the conceptualization, analysis, and writing of the manuscript. All authors read and approved the final version.
Funding. This work was supported by the Russian Science Foundation (project no. 25-29-00339, https://rscf.ru/en/project/25-29-00339/).

---

¹The traditional Taylor number is $\textsf{Ta} = (2\Omega l^2/\nu)^2$; the present choice is necessary to distinguish between positive and negative spatial velocity gradients [6, 9, 32]. 

About the authors

Natalya V. Burmasheva

Ural Federal University named after the first President of Russia B. N. Yeltsin; Institute of Engineering Science, Ural Branch of RAS

Author for correspondence.
Email: nat_burm@mail.ru
ORCID iD: 0000-0003-4711-1894
Scopus Author ID: 57193346922
ResearcherId: E-3908-2016
https://www.mathnet.ru/eng/person52636

Cand. Tech. Sci.; Associate Professor; Dept. of Information Technology and Automation; Senior Researcher; Sect. of Nonlinear Vortex Hydrodynamics

Russian Federation, 620002, Ekaterinburg, Mira st., 19; 620049, Ekaterinburg, Komsomolskaya st., 34

Evgeniy Yu. Prosviryakov

Ural Federal University named after the first President of Russia B. N. Yeltsin; Institute of Engineering Science, Ural Branch of RAS

Email: evgen_pros@mail.ru
ORCID iD: 0000-0002-2349-7801
SPIN-code: 3880-5690
Scopus Author ID: 57189461740
ResearcherId: E-6254-2016
https://www.mathnet.ru/eng/person41426

Dr. Phys. & Math. Sci.; Professor; Dept. of Information Technology and Automation; Head of Sector; Sect. of Nonlinear Vortex Hydrodynamics

Russian Federation, 620002, Ekaterinburg, Mira st., 19; 620049, Ekaterinburg, Komsomolskaya st., 34

Mikhail Yu. Alies

Udmurt Federal Research Center of the Ural Branch of the Russian Academy of Sciences

Email: aliesmy@mail.ru
ORCID iD: 0000-0001-8853-5365
https://www.mathnet.ru/eng/person37555

Dr. Phys. & Math. Sci., Professor; Director

Russian Federation, 426067, Izhevsk, T. Baramzina str., 34

References

Supplementary files