Introduction. Precise control of flow structure in thin layers is critically important for modern technologies — from microfluidic devices to lubrication systems. However, in such cases, some classical hydrodynamic models fail to account for a complex set of technological parameters: wide range of working fluid viscosities, wall permeability, as well as their spatial inhomogeneity (roughness, distributed sources) [1]. An example of such an “incomplete” model is the Couette flow between impermeable plates [2]. The presented study is relevant for the advancement of the following applied areas.

The literature analysis made it possible to systematize approaches to describing the aforementioned factors and revealed a significant gap. Both for model [6] and applied problems [7], the effects of boundary permeability determining nonlinear (exponential) velocity profiles is investigated in detail. The effect of complex (including polynomial) boundary conditions on impermeable walls is well described, which is related to problems of stability, heat and mass transfer, and roughness modeling [8]. The effect of rheological properties, including the wide range of viscosities of real fluids, is deeply studied [9]. Modern numerical [10] and analytical [11] methods are actively used to solve conjugate problems. However, the aforementioned factors have been traditionally studied outside of a systems approach. Analytical solutions either describe permeability under homogeneous boundary conditions [12] or account for boundary nonlinearity only for impermeable walls [13]. Thus, the literature lacks a comprehensive analytical solution that explicitly combines the key dimensionless parameters:

It is this gap that prevents a direct parametric analysis of the synergistic influence of the aforementioned factors on the flow structure. Consequently, it is impossible to purposefully design devices for the applied areas listed above.

The objective of this work is to create and analyze an exact solution for generalized stationary Couette flow of a Newtonian fluid between permeable plates with a quadratic velocity profile at the boundary. The scientific novelty lies in obtaining a closed-form analytical solution that, for the first time, explicitly and comprehensively accounts for the synergy of the parameters α, A, B and μ. The primary advantage of this approach over numerical modeling is the ability to obtain a solution instantly and perform a direct analysis of physical dependences, unmediated by grid approximations. To achieve this objective, four tasks are addressed.

Materials and Methods. Let us consider a steady (stationary), laminar, isothermal flow of a Newtonian viscous and incompressible fluid in a channel formed by two infinite, parallel, flat plates (Fig. 1).

The distance between the plates is constant and equal to h. The bottom plate is stationary and located in the plane z = 0. The top plate is located in the plane z = h and moves with constant velocity W in the positive x direction.

A key feature of this problem is the permeability of both plates, which enables fluid filtration in the direction normal to their surfaces [6]. It is assumed that this normal velocity (Vw) is constant over the entire surface of each plate and is the same for both plates. If Vw > 0, fluid enters the channel through the bottom plate and exits through the top plate; if Vw < 0, the flow direction through the walls is reversed.

The longitudinal velocity component Vx is specified as a quadratic expansion in the transverse coordinate y. This generalization of the classical Couette profile allows for more complex boundary conditions on the top plate [5]:

(1)

Here, functions U(z), u1(z) and u2(z) represent the primary (axial), linear, and quadratic modes of the velocity profile, respectively, depending only on the coordinate z.

Boundary conditions are formulated on the basis of the no-slip condition. On the stationary bottom plate (z = 0), a complete absence of motion is specified:

(2)

On the top plate (z = h), a generalized velocity profile is specified:

(3)

where W — rate of translational motion of the plate; A — coefficient determining the linear velocity gradient along y; B — coefficient determining the profile curvature.

Unknown functions U(z), u1(z) and u2(z) can be found from the Navier–Stokes equation for the x‑component of velocity. Under the assumption of stationarity, it takes the form [2]:

(4)

The flow under consideration is a generalization of pure Couette flow, which is characterized by the absence of a pressure gradient along the channel (∂p/∂x = 0). Taking this condition into account, as well as the facts that Vx is independent of x, Vy = 0 and Vz = Vw = const [3]. In this case, equation (4) is simplified significantly:

(5)

Substituting (1) into (5) and computing the corresponding derivatives yields an equation that must hold for all values of y. Equating coefficients of like powers of y gives a system of three independent ordinary differential equations (ODE) [11]:

(6)

(7)

(8)

Expressions (6) and (7) are second-order linear homogeneous ordinary differential equations (ODE) with constant coefficients. To solve them, the dimensionless permeability parameter α = ρVw / μ is introduced, which physically represents the ratio of inertial forces to viscous forces in the normal direction. Dividing (6) and (7) by μ gives the canonical form:

(9)

The general solutions to these equations are:

(10)

Using boundary conditions (2) and (3) allows us to determine the integration constants. For function u1(z), the conditions u1(0) = 0 and u1(h) = A yield:

(11)

Similarly, for u2(z) from the conditions u2(0) = 0 and u2(h) = B, we find:

(12)

Final expressions for the linear and quadratic modes of the velocity profile:

(13)

(14)

Equation (8) for the fundamental mode U(z) is inhomogeneous. We substitute into it the expression found for u2(z) (14) and divide by μ:

(15)

For equation (15), the method of undetermined coefficients [13] is applied. Its general solution is the sum of the general solution to the homogeneous equation and a particular solution to the non-homogeneous equation. After determining the constants of integration from the boundary conditions (2) and (3), the final expression for the primary mode takes the form [14]:

(16)

Thus, substituting (13), (14), and (16) into the original ansatz (1) yields the full field of the longitudinal velocity component. The analytical solution obtained satisfies both the equation of motion (5) and all boundary conditions, i.e., (2) and (3).

To confirm the physical validity of the solution, let us consider the limiting transition Vw → 0, which is equivalent to α → 0. In this limit, using the Taylor series expansion of the exponential, we obtain for the quadratic mode:

(17)

Similarly for the linear mode:

(18)

For the fundamental mode:

(19)

Numerical modeling. At this stage, the following is performed:

To this end, the authors of the presented article conducted comprehensive numerical modeling in the MATLAB environment. The primary objectives of the software implementation were:

The central element of the code is solve_couette_full function, which implements formulas (13), (14), and (16).

The inputs are:

The output is the basic functions U(z), u1(z), u2(z) and their derivatives with respect to z for the subsequent calculation of such important hydrodynamic quantities as shear stress and vorticity.

The modeling was conducted in two stages. In the first stage, a basic calculation was performed for a standard set of parameters. This allowed for verification of the boundary conditions. The maximum deviation of the velocity from the specified values at the boundaries z = 0 and z = h did not exceed 10–12, which confirmed the high accuracy of the analytical approach and the correctness of the software implementation.

In the second stage, the code was extended to analyze two key aspects. First, the effect of the signs and magnitudes of the coefficients A and B in the boundary velocity profile on the top plate was investigated in detail. Four combinations were considered: A > 0, B > 0; A > 0, B < 0; A < 0, B > 0; A < 0, B < 0. This allowed for the simulation of several physical scenarios (from flow acceleration toward the channel center to its deceleration near the walls), as well as the transition from convex to concave velocity profiles. For each combination, two-dimensional maps of the velocity distribution Vx(y, z) and velocity profiles along the normal to the plates at fixed values of the transverse coordinate y were constructed. This provided a clear visual representation of how changing the sign of the quadratic term B altered drastically the shape of the velocity isolines, and how the sign of the linear coefficient A determined the asymmetry of the profile relative to the plane y = 0.

Secondly, the flow for various fluids was calculated. The model incorporated the physical properties of water SAE 300F1 engine oil, and petroleum at temperatures of 20°C, 40°C, or 60°C [14]. The use of reference data for dynamic viscosity and density made it possible to move from an abstract mathematical model to practically significant engineering calculations. For each fluid, the Reynolds number Re=Wh/ν was automatically determined, which made it possible to assess the flow regime. As expected, for water at room temperature, the Reynolds number reached Re≈1000. At the same time, the flow remained laminar, which was consistent with the known stability of classical Couette flow. For oil and petroleum, the flow remained deeply laminar (Re < 20).

Particular attention was paid to the calculation of derived quantities. For each value of the coordinate y, the shear stress τxz = μ∂Vx/∂z and the vorticity ωy = −∂Vx/∂z were computed analytically. This made it possible to avoid errors of numerical differentiation [15] and ensured high accuracy in constructing the distributions of the specified quantities both along the gap between the plates and as a function of the transverse coordinate y.

Thus, the program developed in MATLAB is a tool for the comprehensive analysis of generalized Couette flow. It can be used to:

Functions u₁(z), u₂(z), U(z) were parametrically investigated using the numerical modeling method in MATLAB based on the analytical solution derived from (13), (14), and (16). The calculations were performed with the following fixed parameters (Figs. 2–5): h = 1.0 m, W = 0.3 m/s, Vw = 0.001 m/s, μ = 0.01 Pa·s, for four combinations of signs: A = ±0.03 s–1 and B = ±0.005 m–1·s–1.

For real fluids, calculations were performed with h = 0.02 m, W = 0.05 m/s, A = 0.1 s–1, B = 0.02 m–1·s–1, Vw = 0.0005 m/s. The following were used: water at 20°C (μ = 1.002·10–3 Pa·s) and 40°C (μ = 0.653·10–3 Pa·s); petroleum at 20°C (μ = 9.15·10–3 Pa·s) and 40°C (μ = 4.72·10–3 Pa·s);); SAE 30 oil at 20°C (μ = 0.290 Pa·s) and 60°C (μ = 0.050 Pa·s).

Research Results. Figure 2 shows the vorticity fields ωy. For A > 0, the maximum positive values ωy are located in the region of positive coordinates y. For A < 0, maxima ωy are shifted to the region y < 0. The amplitude of the values ωy varies with changes in coefficient B.

Figure 3 shows the velocity magnitude field |V| and the vector field for two combinations of the signs of parameters A and B. For A > 0, the vectors are longer in the region y > 0, indicating a shift of the velocity maximum towards the positive transverse coordinate region. For A < 0, the opposite pattern is observed: the vectors are longer in the region y < 0, which indicates a shift of the maximum towards the negative region. This is a manifestation of flow asymmetry caused by the linear term A in the boundary condition. The shape of isolines |V| differs for the cases B > 0 and B < 0, which is related to the effect of the quadratic term B.

Figure 4 presents the isolines of the longitudinal velocity Vx and the streamlines. The streamlines deviate from the direction of the channel axis. The direction of the deviation is determined by the sign A. The isolines Vx for B > 0 and B < 0 have different curvatures.

Figure 5 shows the profiles Vx(y) at different levels, i.e., at z = 0, h, h/4, h/2 and 3h/4. At z = 0, the velocity Vx is zero. As z increases, the variation Vx with y increases and reaches the maximum at z = h, where the profile is determined by parameters A and B. At z = h/2 , the profile is close to linear.

The calculations were performed for real fluids with h = 0.02 m, W = 0.05 m/s, A = 0.1 s–1, B = 0.02 m–1·s–1, Vw = 0.0005 m/s (Figs. 6–8). The following were used: water at 20°C (μ = 1.002·10⁻³ Pa·s) and 40°C (μ = 0.653·10–3 Pa·s); petroleum at 20°C (μ = 9.15·10–3 Pa·s) and 40°C (μ = 4.72·10–3 Pa·s); SAE 30 oil at 20°C (μ = 0.290 Pa·s) and 60°C (μ = 0.050 Pa·s).

Figure 6 shows the profiles Vx(z) at y = 0. All profiles are plotted within the range z ∈ [0, 20] mm, which corresponds to the layer thickness.

Figure 7 displays the isolines |V| For water, the streamlines show minimal deviation from the horizontal, while for petroleum at 20°C, they are significantly curved near the upper wall.

Figure 8 shows two-dimensional vorticity fields ωy for various fluids.

Discussion. According to the data in Figure 6, for water at 20°C, the profile Vx(z) exhibits a weak deviation from the linear law, which is associated with its low viscosity. When the temperature increases to 40°C, this deviation becomes even less pronounced, and the overall slope of the profile slightly decreases. For SAE 30 oil (20°C), a moderate non-linearity of the profile is observed, which weakens significantly upon heating to 60°C. In this case, the profile approaches the linear one, which corresponds to a decrease in viscosity with increasing temperature. For petroleum (20°C), the profile is characterized by a pronounced non-linearity. When the temperature increases to 40°C, the curve becomes closer to the linear one, which is also consistent with a decrease in viscosity.

As shown in Figure 8, for water at 20°C and 40°C, the vorticity is negative and varies in the range from –35 to –5 s–1, which corresponds to weak vorticity throughout the channel. For SAE 30 oil at 20°C, a localized region of elevated vorticity magnitude (|ωy| > 2,3 с–1), is observed, concentrated near the upper wall. Upon heating to 60°C, the region of elevated vorticity expands, and its magnitude decreases, which is consistent with the reduction in viscosity. For petroleum at 20°C, vorticity reaches values up to –4.5 s–1, and at 40°C up to –1.5 s–1, which also corresponds to a decrease in viscosity.

The research results allow us to assert that the combined action of permeability, nonlinear boundary conditions and viscosity, determines qualitatively new flow regimes that cannot be reduced to a simple superposition of known effects [5].

The decisive influence of the sign of coefficient A on flow asymmetry (Figs. 2–5) is explained by its physical meaning as the velocity gradient ∂Vx/∂y at the boundary. This gradient, specified at the upper wall, generates a transverse flow component throughout the entire volume due to viscous diffusion, leading to a shift in streamlines and vorticity maxima. Coefficient B, governing the curvature of the profile (∂²Vx/∂y²), determines the distribution of this transverse flow across the channel width, which is reflected in the change in the shape of the velocity isolines (Figs. 3, 4). The weakening of the effect of the boundary coefficients with distance from the upper wall (Fig. 5) is quantitatively described by the exponential dependences (13) and (14), which is consistent with the concept of disturbance damping from the boundary in a viscous fluid [16].

The dimensionless permeability parameter α depends on the dynamic viscosity μ (α = ρVw / μ). This explains the pronounced difference between the flow of water and petroleum (Figs. 6, 8).

For high-viscosity fluids (low α values), the solution far from the boundaries tends to a linear profile. However, the inhomogeneity in equation (8) for U(z), caused by the term with B, leads to the concentration of shear in a thin near-wall layer, whose thickness is inversely proportional to |α|. Thus, viscosity is not merely a multiplicative factor but a parameter governing the spatial localization of shear deformation. This has significant practical implications, for example, in calculating friction in bearings with porous lubrication [17].

The obtained solution is verified, as it correctly reproduces known limiting cases. For α→0, it transforms into the solution for flow with a quadratic condition on an impermeable wall [18]. For A, B → 0, it reduces to the classical exponential profile for permeable walls [19]. The simultaneous transition α, A, B → 0 yields the linear profile of classical Couette flow [2].

The explicit form of the solution allows for direct assessment of the effect of each parameter on the flow field, which is valuable for engineering design. For example, in a microfluidic mixer, the microrelief of the wall can be varied (via A and B) to generate a specific vortex structure that enhances mixing [20]. In the problem of porous bearing lubrication, the model enables an analytical connection between oil viscosity, seepage rate, and surface roughness (via A, B), and the distribution of shear stress and energy dissipation [21].

The basic limitation of the model is its assumptions of steady, laminar, and Newtonian flow behavior. Extending the model to account for turbulence is beyond the scope of this research. This represents a separate complex problem, requiring a transition to the Reynolds-averaged Navier-Stokes equations (RANS) or Large Eddy Simulation (LES) models. [20]. Another promising direction is generalization for non-Newtonian fluids [21] and unsteady regimes [22]. This would allow for the modeling of a wider class of applied problems, such as pulsating flows in biomedical devices.

Conclusion. For the first time, an exact analytical solution to the stationary Navier-Stokes equations has been obtained, describing generalized Couette flow of a Newtonian fluid between permeable plates with a quadratic velocity profile at the boundary.

The parametric analysis has revealed that the linear coefficient A of the boundary condition determines the direction of asymmetry in the velocity and vorticity fields, while the quadratic coefficient B determines the degree of their spatial nonlinearity.

Dynamic viscosity is shown to be a key parameter controlling the thickness of the shear layer. For high-viscosity media, the main velocity gradient is localized in a thin near-wall region, whereas for low-viscosity fluids, the velocity profile is close to linear across the entire channel height.

The results of this research provide an analytical foundation for solving applied problems in microfluidics, membrane technologies, and tribology, where flow control under conditions of boundary permeability and complex boundary conditions is required.

Further research is related to refining and expanding the model. In perspective, it could be extended to account for the non-Newtonian properties of the fluid, unsteady regimes, and flow stability.

1. From the name of the American organization Society of Automotive Engineers — Association of Automotive Engineers.