Analytical Solution of the Navier–Stokes Equations for Describing Inhomogeneous Couette Flow with a Quadratic Velocity Profile in a Layer with Permeable Boundaries

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.

1. Biomedical microfluidics and so-called labs-on-a-chip. Efficient mixing of reagents or targeted cell delivery require, firstly, fluid injection through porous membranes (V_w ≠ 0). Secondly, heterogeneous shear flows generated by wall microtopography (modelled by boundary coefficients A and B) are needed [3].
2. Energy-efficient membrane technologies (desalination, gas separation). In this area, performance is determined by the interaction between the longitudinal flow and transverse filtration through a permeable wall, as well as by the viscous properties of the medium being separated [4].
3. Tribology of porous and textured bearings. Pressure distribution, friction, and wear directly depend on the lubricant flow in a microscopic gap with complex topography (A, B) and potential filtration (V_w) through a porous liner [5].

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:

• dynamic viscosity (μ);
• permeability (α = ρV_w/μ);
• linear (A) и quadratic (B) corrections to the boundary velocity profile.

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.

1. Output and rigorous verification of the analytical solution, including checking the limiting transitions to known special cases.
2. Investigation of the effect of the sign and magnitude of the boundary coefficients A and B on the spatial distribution of velocity and vorticity.
3. Quantitative analysis of the effect of the dynamic viscosity of real fluids (water, engine oil, petroleum) on the velocity profile and shear layer thickness.
4. Discussion of the practical significance of the results, model limitations, and promising directions for the model generalization.

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 (V_w) is constant over the entire surface of each plate and is the same for both plates. If V_w > 0, fluid enters the channel through the bottom plate and exits through the top plate; if V_w < 0, the flow direction through the walls is reversed.

The longitudinal velocity component V_x 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), u₁(z) and u₂(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), u₁(z) and u₂(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 V_x is independent of x, V_y = 0 and V_z = V_w = 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 α = ρV_w / μ 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 u₁(z), the conditions u₁(0) = 0 and u₁(h) = A yield:

(11)

Similarly, for u₂(z) from the conditions u₂(0) = 0 and u₂(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 u₂(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 V_w → 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:

• verification and analysis of the obtained solution;
• investigation of the effect of physical parameters on the flow structure.

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:

• visual representation of velocity fields;
• quantitative comparison of flow characteristics for various fluid types and boundary conditions.

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

The inputs are:

• physical properties of the fluid (density ρ and dynamic viscosity μ);
• geometric parameters of the problem (distance between the plates h);
• kinematic characteristics (W, A, B, V_w);
• z-coordinate grid.

The output is the basic functions U(z), u₁(z), u₂(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 V_x(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 300F¹ 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 = μ∂V_x/∂z and the vorticity ω_y = −∂V_x/∂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:

• verify the correctness of the analytical solution;
• conduct an in-depth study of the effect of the physical properties of the working medium and complex boundary conditions on the flow hydrodynamics;
• obtain extensive visual and quantitative data for further analysis.

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, V_w = 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, V_w = 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 V_x 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 V_x for B > 0 and B < 0 have different curvatures.

Figure 5 shows the profiles V_x(y) at different levels, i.e., at z = 0, h, h/4, h/2 and 3h/4. At z = 0, the velocity V_x is zero. As z increases, the variation V_x 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, V_w = 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 V_x(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 V_x(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 ∂V_x/∂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 (∂²V_x/∂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 μ (α = ρV_w / μ). 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.