Stability Analysis of Navier–Stokes–Voigt Fluids in Porous Media with Slippery Effect

TL;DR

This study examines how slip conditions and porous media affect the stability of fluid flow, showing that symmetric slip and strong damping can prevent instability.

Contribution

The paper introduces a stability analysis of NSV fluids in porous media with asymmetric slip, revealing how symmetry breaking and material properties influence flow stability.

Findings

Symmetric slip conditions lead to optimal flow stability, while asymmetric slip causes significant destabilization.

Increasing porous medium permeability or Voigt regularization can counteract slip-induced instability.

Manufacturing defects in slip conditions can be mitigated by optimizing the porous medium matrix.

Abstract

This paper investigates the linear stability of Navier–Stokes–Voigt (NSV) fluid flow in a channel filled with a homogeneous porous medium under general asymmetric slip boundary conditions. This study bridges the research gap between idealized theoretical models (uniform coating) and realistic engineering surfaces in superhydrophobic channels. In practice, manufacturing defects often lead to non-uniform slip distributions. By solving the generalized eigenvalue problem using the Chebyshev spectral collocation method, we quantify the sensitivity of the critical Reynolds number to symmetry breaking. The results reveal that symmetric slip achieves optimal stability, whereas symmetry breaking causes a significant destabilizing effect. Energy analysis clarifies the physical origin of this instability. Furthermore, we find that increasing the porous medium permeability parameter or the Voigt regularization parameter effectively counteracts the slip-induced instability. Specifically, flow stability can be restored even under highly asymmetric slip conditions if the porous damping or the viscoelastic regularization effect is sufficiently strong. This implies that inevitable manufacturing defects in engineering can be compensated for by optimizing the porous medium matrix.

Figures

Full text

1. Introduction

For decades, fluid system stability research has remained a highly compelling and prominent direction in fluid mechanics. Specifically, the stability of incompressible channel flow is a fundamental problem in fluid mechanics, which has been extensively studied from a theoretical perspective [1]. Beyond its theoretical significance, this problem has gained increasing attention in various engineering applications, such as microfluidic transport, filtration processes in porous media, and lubrication systems [2]. The roots of this research trace back to Reynolds’ pioneering experiment in 1883 [3]. This work marked the start of modern fluid stability studies and inspired Orr’s [4] and Sommerfeld’s [5] theoretical investigations. Working independently, they considered small-amplitude traveling-wave perturbations imposed on steady parallel flows. From this, they derived the classic Orr–Sommerfeld equation. Subsequent scholars conducted in-depth studies on various flow configurations. These include planar Poiseuille flow (PPF), planar Couette flow (PCF), Hagen–Poiseuille flow (HPF), and boundary layer flow [6]. These efforts have progressively refined the theoretical framework of fluid stability.

However, traditional fluid stability studies often rely on the no-slip boundary condition assumption. Research on small-scale fluid flows shows that the commonly assumed “smooth surface” fails to meet hydrodynamic smoothness criteria. This leads to significant deviations when the classical no-slip model describes micro- and nano-scale flows. As early as 1823 [7], Navier first proposed a mathematical model for slip flow:

$$[eqn]$$

In this equation, $[eqn]$ denotes the slip length, $[eqn]$ is the slip velocity component in the flow direction, $[eqn]$ indicates the wall position, and y represents the wall normal direction. This condition establishes that the tangential velocity near the surface is proportional to the surface shear rate. When $[eqn]$ , it corresponds to the no-slip condition (i.e., $[eqn]$ ). When $[eqn]$ , it represents an ideal slip surface. Subsequently, after a century of theoretical exploration, Vinogradova [8] provided a comprehensive review and further quantified the slip phenomena in uniform high-hydrophobicity systems through experimental comparison. Experimental observations showed a slip length of $[eqn]$ for surfaces with contact angles $[eqn]$ . With advances in materials science, research on nanoscale superhydrophobic surfaces has deepened. Tretheway and Meinhart [9] used micrometer-resolution particle image velocimetry ( $[eqn]$ -PIV) to measure significant apparent slip in hydrophobic microchannels coated with a single layer of OTS. The measured slip length reached $[eqn]$ . Choi’s [10] precise experiments demonstrated that slip length on hydrophobic surfaces increases approximately linearly with shear rate.

Additionally, under identical pressure differentials, flow rates in hydrophobic channels consistently exceed those in hydrophilic channels. This provides direct evidence of the positive correlation between slip length and drag reduction effects. These properties have proven valuable in practical applications, such as efficient transport in microfluidic chips and sweat dissipation in wearable electronics [11]. Regarding flow stability under slip boundary conditions, academic debate once existed. Chu [12] incorrectly concluded that critical Reynolds number decreases with increasing slip length. This error stemmed from improper boundary condition treatment. In contrast, Lauga and Cossu [13] and Spille [14] obtained opposite results. Later, Webber [15] used the accurate Chebyshev tau numerical method to verify the stabilizing effect of slip boundaries. This finding aligned with Lauga et al.’s conclusions.

Notably, most previous channel flow stability studies focused on Newtonian fluids. Non-Newtonian fluids—such as polymer solutions, blood, and emulsions—are now widely used in engineering. As a result, the Newtonian fluid assumption can no longer meet practical demands. Research on non-Newtonian slip flow has thus emerged as a hot topic. For example, Patlazhan and Vagner [16] investigated shear-thinning fluid flow on superhydrophobic surfaces. Their work provided a theoretical basis for designing high-efficiency drag-reducing microfluidic devices. Ferrás et al. [17] systematically analyzed viscoelastic fluid behavior under slip boundary conditions. Their findings have been applied in polymer processing, lubrication, and coating technologies.

The classical Navier–Stokes equations serve as the core theoretical framework for describing fluid motion, but their three-dimensional well-posedness remains an unsolved challenge. Potential singularities severely restrict theoretical analysis and numerical computation, prompting mathematical regularization as a critical research direction and leading to the development of the Navier–Stokes–Voigt (NSV) model. Notably, the NSV model is a simplified version of the second-grade fluid model [18,19], retaining key viscous and elastic stress terms related to strain rate and its time derivative while omitting higher-order nonlinear terms. This simplification ensures mathematical tractability while capturing core dissipative and regularization effects, making it suitable for stability analysis under complex boundary conditions. In his pioneering work, Oskolkov first proposed the NSV equations to describe the motion of viscoelastic incompressible fluids following the Kelvin–Voigt model. He also established a comprehensive framework covering several important viscoelastic fluid models, including Oldroyd fluids and Maxwell fluids. Unlike classical models such as Oldroyd-B (derived from polymer dynamics), the Navier–Stokes–Voigt (NSV) model introduces an elastic stress term proportional to the time derivative of the strain rate into its constitutive equations. This creates a highly tractable mathematical framework. While this simplification prevents the model from predicting typical viscoelastic phenomena (e.g., shear thinning and normal stress differences), it effectively suppresses singularities inherent in the original Navier–Stokes equations. Subsequently, Oskolkov and Shadiev [20] rigorously proved the existence of global classical solutions for initial-boundary value problems of Oldroyd and Kelvin–Voigt fluids over infinite time intervals. They systematically summarized the motion equations for these viscoelastic fluids, laying the theoretical foundation for NSV model applications. Although Oskolkov’s work initially confirmed the NSV model’s mathematical regularization properties, its connection to the original Navier–Stokes equations required explicit verification. In 2012, Berselli and Bisconti [21] conducted the first systematic analysis of the limit behavior of Euler–Voigt and Navier–Stokes–Voigt (NSV) models. Through rigorous proof, they showed that under appropriate initial conditions, as the Voigt parameter approaches zero, the model solutions strongly converge to those of the original Euler or Navier–Stokes equations. This conclusion provides core theoretical support for applying Voigt-type models in turbulence simulation and computational fluid dynamics. It also transformed these models from purely mathematical tools to practical research instruments for flow problems.

Building on Oskolkov’s foundational work and subsequent theoretical advancements, the NSV model has gradually been applied to stability analysis in complex scenarios—such as flow through porous media. Here, the influence of slip boundary conditions has become a key research focus. Straughan et al. [22] conducted pioneering studies on the effect of slip boundary conditions on Poiseuille flow instability in Brinkman-type porous media. Based on the Brinkman porous medium model, they derived flow control equations incorporating slip boundary conditions. They introduced two key parameters: a dimensionless slip coefficient (characterizing slip length) and a Darcy friction coefficient (characterizing porous medium resistance). To ensure the accuracy of critical parameter calculations, they validated their approach using two Chebyshev configuration methods, the finite element method, and the QZ algorithm for solving eigenvalue problems. Their study yielded several important findings: With fixed porous medium parameters, increasing the slip coefficient significantly raises the critical Reynolds number. This indicates that slip boundaries exert a strong stabilizing effect on flow, delaying laminar–turbulent transition. With fixed slip coefficient, increasing the porous medium resistance parameter (reducing porosity and increasing resistance) further increases the critical Reynolds number. This synergizes with the stabilizing effect of slip boundaries. As slip length increases, the critical wave number decreases, indicating sparser periodic structures in flow disturbances. As the porous medium coefficient increases, the critical wave number first decreases and then increases, with an inflection point in the 2–3 range. This research clarified the stabilizing effects of increased Darcy friction coefficient and slip length on flow stability. The Brinkman model is suitable for describing flow in large-pore porous media, while the Darcy model applies to small-pore scenarios. Together, they provide a foundation for boundary condition analysis in Navier–Stokes–Voigt fluid stability studies within porous media.

While the stability of Navier–Stokes–Voigt (NSV) fluids has been studied in various contexts, most research focuses on clear fluid channels. Recently, Shankar and Shivakumara [23] extended this analysis to porous media; nevertheless, their investigation was limited to idealized symmetric slip conditions ( $[eqn]$ ). In reality, this assumption of perfect symmetry often fails to capture engineering complexities for two key reasons. First, unavoidable manufacturing tolerances in superhydrophobic coatings—such as uneven spray thickness—inevitably result in non-uniform slip lengths. This fabrication-induced asymmetry is a phenomenon that symmetric models cannot accurately describe. Second, in advanced microfluidic design, intentionally using asymmetric coatings (e.g., $[eqn]$ ) allows for the “directional regulation” of flow stability. This approach offers significantly greater flexibility than symmetric configurations. Consequently, the impact of slip asymmetry on flow stability remains a critical but unresolved issue.

To bridge the gap between idealized theory and these practical engineering scenarios, this study performs a comprehensive linear stability analysis of NSV fluid flow under general asymmetric slip conditions. The novel contributions of this work are threefold. First, instead of assuming perfect symmetry, we quantify the sensitivity of the critical Reynolds number to slip asymmetry using the Chebyshev spectral collocation method. This establishes a theoretical basis for evaluating the allowable range of manufacturing defects in engineering design. Second, we investigate the coupling mechanism between the fluid’s viscoelastic regularization (Voigt term) and the Darcy drag of the porous medium in the presence of asymmetry. Finally, we propose a strategy to counteract the instability caused by manufacturing-induced asymmetry. This is achieved by optimizing the porous medium permeability parameter (M) or the Voigt parameter ( $[eqn]$ ).

The remainder of this paper is organized as follows: Section 2 establishes the mathematical framework for NSV fluid flow in porous media. Specifically, Section 2.1 and Section 2.2 present the governing equations and the steady base flow solutions under varying slip conditions; Section 2.3 derives the modified Orr–Sommerfeld perturbation equation; Section 2.4 discusses the analogy between asymmetric slip flow and Poiseuille–Couette flow; and Section 2.5 formulates the disturbance energy budget equation. Finally, Section 2.6 justifies the selection of physical parameters based on experimental references. Section 3 discusses the numerical results and physical mechanisms. It begins with the solution of the generalized eigenvalue problem via the Chebyshev spectral collocation method and code validation. Subsequently, the neutral stability curves and critical Reynolds numbers are analyzed with respect to key parameters, including slip asymmetry, porous medium permeability, and the Voigt parameter. Finally, based on the global energy budget analysis, the physical origins of slip-induced instability and the corresponding compensation mechanisms are thoroughly investigated. Section 4 summarizes the conclusions.

2. Mathematical Formulation

Consider a layer of incompressible Navier–Stoke–Voigt fluid confined between two infinitely long rigid parallel plates. The distance between the plates is $[eqn]$ , and the flow is driven by a constant pressure gradient, with the system embedded in a uniform porous medium. According to Squire’s theorem, for the basic flow of incompressible laminar flow, if a two-dimensional unstable disturbance exists, there must be a three-dimensional unstable disturbance with a lower critical Reynolds number. Consequently, our stability analysis only needs to consider the two-dimensional model. A Cartesian coordinate system is adopted, where the x-axis is parallel to the flow direction, the y-axis points upward, and the origin is located at the center of the channel, as illustrated in Figure 1.

Consider a layer of incompressible Navier–Stokes–Voigt fluid confined between two infinitely long rigid parallel plates. The distance between the plates is $[eqn]$ , and the flow is driven by a constant pressure gradient, with the system embedded in a uniform porous medium. It is well known that Squire’s theorem does not universally apply to all non-Newtonian fluids. However, because the viscoelastic Voigt regularization term is mathematically linear, recent rigorous derivations by Shankar and Shivakumara [23] have established that the classical Squire’s transformations remain perfectly valid for the NSV fluid model. According to this validated Squire’s theorem, if a three-dimensional unstable disturbance exists, there must be a two-dimensional unstable disturbance with a lower critical Reynolds number. Consequently, our stability analysis rigorously only needs to consider the two-dimensional model. A Cartesian coordinate system is adopted, where the x-axis is parallel to the flow direction, the y-axis points upward, and the origin is located at the center of the channel, as illustrated in Figure 1.

Herein, we generalize the model for the flow of viscoelastic fluids in channels: the fluid exhibits a specific slip length $[eqn]$ at the upper plate and a specific slip length $[eqn]$ at the lower plate. When $[eqn]$ , the flow corresponds to Poiseuille flow; when $[eqn]$ , it represents symmetric flow with slip conditions; and the configurations $[eqn]$ and $[eqn]$ are physically equivalent. To avoid redundancy and maintain mathematical generality, we only discuss the general asymmetric case where $[eqn]$ .

The mathematical formulation of the above model is given by the continuity equation and the Navier–Stokes equations:

$$[eqn]$$ $$[eqn]$$

$[eqn]$ is the stress tensor of the zeroth-order Navier–Stokes–Voigt fluid as given by Oskolkov [24] and Straughan [25].

$$[eqn]$$

Here, $[eqn]$ is the unit tensor and $[eqn]$ denotes the transpose of the indicated tensor. In addition, we re-scale the elastic parameter and we replace $[eqn]$ simply with $[eqn]$ . From the Dupuit–Forchheimer relation $[eqn]$ [26], we get

$$[eqn]$$

Herein, $[eqn]$ denotes the velocity vector, $[eqn]$ represents the fluid density, $[eqn]$ is the Darcy velocity, t stands for time, and p denotes the pressure. The permeability k characterizes the conductivity of the porous medium: a larger k corresponds to a smaller resistance to fluid flow, while a smaller k leads to a larger flow resistance. The porosity $[eqn]$ of the porous medium refers to the ratio of the pore volume to the total volume of the porous medium; a smaller $[eqn]$ results in narrower effective flow channels. $[eqn]$ is defined as the Kelvin–Voigt parameter. Under all considered cases, the half-width h of the channel is selected as the length scale, and the maximum fluid velocity $[eqn]$ is adopted as the velocity scale to normalize $[eqn]$ and k. The appropriate dimensionless treatment is given as follows:

$$[eqn]$$

Substituting Equation (5) back into Equation (4) and omitting the asterisks, the component forms of the nondimensionalized Navier–Stokes equations and continuity equation are obtained as:

$$[eqn]$$

Herein, $[eqn]$ denotes the dynamic viscosity coefficient of the fluid. $[eqn]$ represents the Reynolds number. It roughly describes the ratio of inertial forces to viscous forces in the flow. Its specific expression is given by $[eqn]$ . As a fundamental indicator, this parameter is crucial for determining the flow regime and evaluating flow stability. M is defined as the permeability parameter. It is specifically used to quantify the fluid’s ability to penetrate the porous medium. The magnitude of M is positively correlated with permeability: a larger M indicates stronger resistance for the fluid to penetrate the pores of the porous medium, while a smaller M means weaker permeation resistance for the fluid. Mathematically, it is expressed as $[eqn]$ , where h is the channel half-height, $[eqn]$ is the porosity of the porous medium, and k is the permeability of the porous medium. This parameter directly reflects the inhibitory or facilitatory effect of the porous medium’s structure on fluid flow. $[eqn]$ refers to the Kelvin–Voigt parameter. Its intrinsic physical meaning is to describe the ratio of elastic forces to viscous forces during the flow of the Navier–Stokes–Voigt (NSV) fluid. As a key dimensionless number characterizing the viscoelastic properties of the fluid, its mathematical expression is $[eqn]$ , where $[eqn]$ is the dimensional Kelvin–Voigt elastic parameter and h is the channel half-height. The magnitude of $[eqn]$ directly affects the fluid’s response to perturbations. Generally, a larger $[eqn]$ enhances the fluid’s ability to suppress unstable flow modes, thereby exerting a significant regulatory effect on flow stability.

2.1. Base Flow Analysis

To establish the stability equations under asymmetric slip boundaries, determining the base flow velocity profile for the NSV fluid is essential. The detailed derivation proceeds as follows. Since the Navier–Stokes–Voigt terms in Equation (6) do not contribute to the base flow solution, and the base flow is a fully developed, unidirectional and steady laminar flow, its velocity vector can be expressed as $[eqn]$ . Here, the subscript “b” denotes the basic state. Substituting this velocity vector into Equation (6) yields a second-order ordinary differential equation.

$$[eqn]$$

Let $[eqn]$ represent the pressure of the undisturbed fluid motion. To solve for the base flow, Equation (7) is solved simultaneously with the boundary conditions:

$$[eqn]$$ $$[eqn]$$

The resulting solution for the base flow is as follows:

$$[eqn]$$

where

$$[eqn]$$ $$[eqn]$$

2.2. Boundary Condition

2.2.1. Asymmetric Slip Conditions

When $[eqn]$

$$[eqn]$$

When $[eqn]$

$$[eqn]$$

When $[eqn]$ and $[eqn]$

$$[eqn]$$

When $[eqn]$

$$[eqn]$$

2.2.2. Symmetric Slip Conditions

If $[eqn]$ , then $[eqn]$ takes the form:

$$[eqn]$$

When $[eqn]$

$$[eqn]$$

However, if $[eqn]$ , $[eqn]$ simplifies to:

$$[eqn]$$

When $[eqn]$

$$[eqn]$$

Thus, Figure 2 presents a comparative illustration of the base flows corresponding to the three aforementioned slip scenarios of Newtonian fluids.

2.3. Linear Stability Analysis

Since Squire’s theorem is applicable to the slip flow modeled in this study, the following analysis is restricted to two-dimensional disturbances. To investigate the hydrodynamic stability of the problem, the instantaneous velocity and stress are expressed as the sum of the base state and the perturbation quantities. A linear stability analysis is performed on small perturbations to the basic flow:

$$[eqn]$$

where $[eqn]$ . Substituting Equation (19) into Equation (6) and linearizing with respect to $[eqn]$ , the linearized Navier–Stokes equations are obtained:

$$[eqn]$$

where $[eqn]$ . For simplicity, the hat symbols on the perturbation variables have been omitted, and primes denote differentiation with respect to y. A stream function $[eqn]$ is introduced such that:

$$[eqn]$$

Assuming a wave-like solution of the form:

$$[eqn]$$

For the streamwise wave number $[eqn]$ and phase velocity $[eqn]$ . Take the partial derivatives of the first two equations in Equation (20) with respect to x and y respectively, then subtract the resulting equations to eliminate the pressure-related terms. Finally, substitute Equations (21) and (22) into the pressure-eliminated equation, and Equation (23) can be obtained. When $[eqn]$ , the equation reduces to the well-known Orr–Sommerfeld equation.

$$[eqn]$$

The eigenvalue problem, Equation (23), is then solved subject to an equivalent set of boundary conditions:

$$[eqn]$$ $$[eqn]$$

Herein, $[eqn]$ denotes the base flow, and $[eqn]$ represents the second derivative with respect to y. Consequently, in the temporal stability analysis with a real wave number $[eqn]$ , the complex phase velocity $[eqn]$ serves as the eigenvalue of the problem. Correspondingly, the flow exhibits linear stability when $[eqn]$ and linear instability when $[eqn]$ .

The neutral curve constitutes a fundamental analytical tool in fluid stability studies, rooted in linear stability theory. This approach involves solving the eigenvalue problem of the Orr–Sommerfeld equation to determine the relationship between perturbation growth rates and flow parameters. In the ( $[eqn]$ ) plane, the neutral curve defined by the critical condition of zero growth rate ( $[eqn]$ ) serves as a demarcation boundary separating stable ( $[eqn]$ , decaying perturbations) from unstable ( $[eqn]$ , growing perturbations) regimes. Systematic analysis of the neutral curve’s morphology and position enables determination of: (i) critical flow instability thresholds, (ii) dominant instability modes and their physical mechanisms, and (iii) laminar–turbulent transition criteria. Furthermore, the curve’s evolution reveals how stability characteristics are influenced by flow parameters including wall slip conditions, porosity, and viscoelastic coefficients, providing crucial theoretical insights for flow control strategies.

Numerically, the generalized eigenvalue problem is solved using the Chebyshev spectral collocation method. The interval $[eqn]$ is discretized into N + 1 collocation points, and the disturbance velocity is expanded in a Chebyshev polynomial series. Substituting into the modified Orr–Sommerfeld equation leads to the discrete matrix system $[eqn]$ , where: $[eqn]$ ; $[eqn]$ . The eigenvalue problem is solved using the QZ algorithm in MATLAB R2024b.

2.4. Analogy to the Poiseuille–Couette Flow

The Couette flow model can be regarded as a variant of asymmetric slip. After conducting a stability analysis on the base flow Equation (15), we solve the Orr–Sommerfeld equation Equation (23) under the condition of no applied slip.

$$[eqn]$$ $$[eqn]$$

Pure Couette flow is characterized by a zero pressure drive, $[eqn]$ . Poiseuille–Couette flow is distinguished by a constant pressure, $[eqn]$ .

This leads to the base flow equation:

$$[eqn]$$

When $[eqn]$

$$[eqn]$$

where U is the velocity of the channel wall (in this particular model, the upstream channel wall is stationary). By comparing Equation (29) with the asymmetric base flow Equation (15),

$$[eqn]$$

The following relationship between the wall velocity U and the slip length $[eqn]$ is derived:

$$[eqn]$$

This approach reveals a certain connection between this specific problem and Poiseuille–Couette flow. Thus, the linear stability of Poiseuille–Couette flow is equivalent to analyzing the asymmetric slip flow Equation (15) with the application of no-slip boundary conditions to the Orr–Sommerfeld Equation (23).

2.5. Disturbance Energy Budget Equation

To elucidate the physical mechanisms governing stability enhancement, we employ the classical Reynolds–Orr energy analysis method [1,23]. We take the inner product of the linearized perturbation momentum equation (Equation (20)) with the perturbation velocity vector $[eqn]$ , where u and v denote the streamwise and wall-normal disturbance components, respectively. Then, we integrate over the interval $[eqn]$ in the y-direction. This yields the rate of change in the total disturbance kinetic energy, defined as $[eqn]$ . The derived energy budget equation for the NSV fluid in a porous medium is:

$$[eqn]$$

Here, the terms on the right-hand side correspond to distinct physical processes:

$$[eqn]$$

This term represents the rate at which the perturbation extracts energy from the base flow shear via Reynolds stress. Here, $[eqn]$ denotes the real part of a complex quantity, and the overbar $[eqn]$ denotes the complex conjugate. A positive value ( $[eqn]$ ) serves as the primary source of instability.

$$[eqn]$$

This term accounts for the kinetic energy loss caused by Newtonian viscosity. Mathematically, $[eqn]$ represents the total intensity of fluid deformation. It is strictly negative, indicating a stabilizing effect.

$$[eqn]$$

Unique to this study, this term quantifies the Darcy drag exerted by the porous matrix. It functions as an additional energy sink, uniformly suppressing disturbances across the channel.

$$[eqn]$$

This term captures the contribution of the Voigt viscoelastic term.

In linear stability analysis, it acts as a regularization mechanism that smooths out high-frequency oscillations. Under the framework of linear stability theory, the disturbance growth rate $[eqn]$ is determined by the net energy budget. Instability occurs only when the production $[eqn]$ exceeds the sum of all dissipation mechanisms.

2.6. Physical Parameter Selection

The key dimensionless parameters in this study are selected based on physical feasibility and experimental references:

Voigt parameter ( $[eqn]$ ): Defined as $[eqn]$ ( $[eqn]$ is the dimensional Kelvin–Voigt elastic parameter, h is the channel half-width). For typical viscoelastic fluids (e.g., PEO/PAM aqueous solutions) with relaxation time $[eqn]$ – $[eqn]$ s [27], and microchannel half-width $[eqn]$ – $[eqn]$ m, $[eqn]$ ranges from $[eqn]$ to $[eqn]$ .

Permeability parameter (M): Expressed as $[eqn]$ ( $[eqn]$ is porosity, k is permeability). Covering scenarios from pure fluids ( $[eqn]$ , $[eqn]$ ) to low-permeability porous media ( $[eqn]$ ), it corresponds to porosity $[eqn]$ – $[eqn]$ and permeability $[eqn]$ – $[eqn]$ , typical of soil, sandstone, and shale [28].

Slip length ( $[eqn]$ ): Dimensionless slip length is the ratio of physical slip length to h. Physical slip lengths of 1–10 μm (consistent with experimental measurements [9,10]) correspond to $[eqn]$ = 0–0.05 for $[eqn]$ m.

3. Discussion of Result

3.1. Numerical Convergence and Validation

Next, a convergence test via eigenvalue solution is performed for the Navier–Stokes–Voigt fluid under the conditions of $[eqn]$ , $[eqn]$ , $[eqn]$ , and $[eqn]$ . For the selected parameters ( $[eqn]$ , $[eqn]$ ), the convergence test for the least stable eigenmode ( $[eqn]$ ) is presented with an increasing number of collocation points. Table 1 first conducts the convergence test for the least stable eigenmode of Newtonian fluid under no-slip boundary conditions. Subsequently, the convergence test is also carried out for the case of $[eqn]$ , $[eqn]$ , and $[eqn]$ . This table confirms that when the number of collocation points $[eqn]$ , the numerical scheme achieves good convergence with at least six decimal places. These convergence results lay a benchmark for the subsequent comprehensive analysis. Furthermore, in comparison with the data of Newtonian fluid under no-slip boundary conditions ( $[eqn]$ , $[eqn]$ ) reported by Ceccacci et al. [29], the least stable eigenmode obtained via the numerical method in this study is $[eqn]$ , which is in perfect agreement with their results up to six decimal places. Compared with the least stable eigenmode of $[eqn]$ reported by Arjun et al. [30] for Navier–Stokes–Voigt fluid with $[eqn]$ , $[eqn]$ , $[eqn]$ ( $[eqn]$ =10,000, $[eqn]$ ), our result also shows complete consistency at the six-decimal-place level.

3.2. Compare with Poiseuille–Couette Flow

When $[eqn]$ and $[eqn]$ , the fluid behaves as a Newtonian fluid. Figure 3 illustrates the Poiseuille–Couette flow with one-sided slip under the conditions of $[eqn]$ and $[eqn]$ . Values of $[eqn]$ are taken as $[eqn]$ and the corresponding comparisons with the Couette flow satisfying Equation (30) are plotted. As $[eqn]$ decreases, the two curves tend to coincide. With increasing $[eqn]$ , although U still varies proportionally, a certain discrepancy in the critical Reynolds number emerges. Figure 4 compares the two flow models for a Navier–Stokes–Voigt fluid with $[eqn]$ and $[eqn]$ , under the condition that the relationship between the slip length $[eqn]$ and the Couette plate velocity U, as previously described, remains approximately valid. The discrepancy between the two models again increases with increasing $[eqn]$ . It is found that this approach is only suitable for small slip-length approximations. By performing a Taylor expansion of Equation (30), it is observed that both models exhibit identical linearized behavior to first-order approximation. However, when $[eqn]$ , terms of order $[eqn]$ and higher become non-negligible. For large $[eqn]$ , the discretization error of Chebyshev collocation points near the boundary is amplified, and the condition number of the generalized eigenvalue problem deteriorates. Therefore, the optimal range of agreement is found to be $[eqn]$ . This study is also relevant to stability analysis in the context of superhydrophobic surfaces: when the material slip length is small, the stability of such flows may be approximated by that of Couette flow.

3.3. Eigenspectrum

Figure 5 compares the eigenvalue spectra for Newtonian fluids ( $[eqn]$ ) at $[eqn]$ and $[eqn]$ under both no-slip and slip conditions. As shown in Figure 5a, when no-slip boundaries are examined, the classic Newtonian spectrum exhibits a characteristic Y-shaped structure. This structure consists of three distinct branches—A, P, and S—collectively referred to as the APS branches [31]. In Figure 5b, we introduce slip boundary conditions. The results demonstrate that introducing wall slip causes the A-branch to shift downward into the stable half-plane ( $[eqn]$ ), which aligns with the stabilizing effect reported in previous literature. The APS branches are three types of modes in the spectrum of Newtonian fluids, classified based on the distribution of perturbation energy, the trend of phase velocity, and stability behavior. The A branch corresponds to wall modes, which are perturbation modes dominated by the shear effect of the flow wall. The phase velocity c tends to 0, so the perturbation propagates extremely slowly, almost “attached” near the wall. They are usually damped modes $[eqn]$ , but under specific flow conditions, they may transform into growing modes, triggering flow instability near the wall. Changing the magnitude of some influencing factors can easily make them unstable or restore their stability. The P branch corresponds to central modes, which are perturbation modes dominated by the core region of the flow, with energy mainly concentrated in the central region where the flow velocity is relatively high. The phase velocity tends to the maximum velocity of the base flow, so the perturbation propagates at a speed close to the mainstream velocity. They may exhibit as growing modes $[eqn]$ or damped modes, and are key modes governing the global stability of the flow. The S branch corresponds to damped modes, which is the third family of eigenvalues derived by Schensted [32]. This family has an infinite range, with two asymptotic expressions for its eigenvalues. One is the formula proposed by Schensted: $[eqn]$ (where $[eqn]$ is unspecified), and the other is the formula put forward by Grosch and Salwen [33] in their extensive study of the spectrum of plane Poiseuille flow: $[eqn]$ . According to Grosch and Salwen, the asymptotic value of $[eqn]$ is $[eqn]$ of the average velocity of the mean flow in the channel.

Figure 6a presents a configuration with an upper-wall slip length of $[eqn]$ , a lower-wall slip length of $[eqn]$ , $[eqn]$ = 11,000, and $[eqn]$ . As the viscoelastic parameter $[eqn]$ increases from 0 to $[eqn]$ , the locally magnified view reveals that the most unstable mode crosses the $[eqn]$ -axis, indicating enhanced flow stability with increasing $[eqn]$ . However, it is observed that when $[eqn]$ increases to $[eqn]$ , the S-branch begins to shift leftward. A further increase in $[eqn]$ to $[eqn]$ results in a more pronounced leftward translation of the S-branch. Examination of the local detailed view shows that the most unstable mode now lies below zero. Moreover, even if slip and porosity effects are reduced, the flow remains stable under these conditions. This behavior is consistent with the findings reported by Shankar [23] et al.

Figure 6b explores the eigenvalue spectrum of NSV fluids with larger Voigt parameters (up to $[eqn]$ ). In some non-Newtonian models, increasing viscoelastic parameters triggers severe instability. However, the results for the NSV model are different. As $[eqn]$ increases, the traditional Y-shaped eigenvalue branches collapse toward the real axis. These branches remain strictly within the stable half-plane ( $[eqn]$ ). This confirms that the Voigt term acts as a consistent stabilizing factor.

Figure 6c investigates the effect of the permeability coefficient of porous media on the asymmetric slip flow of viscoelastic fluids. With $[eqn]$ , $[eqn]$ , $[eqn]$ = 11,000, and $[eqn]$ kept constant, as M increases from 0 to $[eqn]$ , it is noted that the unstable mode smoothly moves below $[eqn]$ , as shown in the local enlarged view. This confirms the stabilizing effect of increasing the porous medium parameter M on the flow.

The influence of slip length on the flow of Navier–Stokes–Voigt fluids is illustrated in Figure 6d. With $[eqn]$ and $[eqn]$ fixed, the lower plate has a slip length of $[eqn]$ , while the upper plate slip length $[eqn]$ is varied from $[eqn]$ to $[eqn]$ . The results clearly demonstrate that increasing the slip length plays a significant role in enhancing the stability of this viscoelastic fluid. In summary, the unstable modes in all the aforementioned cases belong to the A-branch.

Finally, a direct comparison between Figure 5b and Figure 6 robustly validates our core conclusions. In the purely Newtonian regime without a porous medium ( $[eqn]$ , as shown in Figure 5b), the specific asymmetric slip configuration ( $[eqn]$ ) yields a maximum eigenvalue $[eqn]$ , indicating an unstable flow. However, as demonstrated in Figure 6a,c, introducing and increasing either the viscoelastic parameter $[eqn]$ or the porous medium parameter M under the identical asymmetric slip condition systematically shifts the most unstable mode back into the stable half-plane ( $[eqn]$ ). This comparison firmly establishes the powerful stabilizing roles of both the Voigt regularization and the Darcy dissipation in further enhancing the overall stability of the fluid system.

3.4. Neutral Stability Curves

To verify the reliability of the numerical scheme, we compare the predicted critical wavenumber ( $[eqn]$ ) and critical Reynolds number ( $[eqn]$ ) for Poiseuille flow in a porous medium under symmetric slip conditions with results from the literature. For validation, parameters $[eqn]$ and $[eqn]$ are used. As shown in Figure 7, as the porous medium parameter M increases from 0 to 2, the computed values of $[eqn]$ and $[eqn]$ (e.g., $[eqn]$ 10,158 at $[eqn]$ ) agree well with those reported by Straughan and Harfash. This confirms the accuracy and reliability of the numerical implementation.

We next investigate the influence of various parameters on flow stability when the top and bottom wall slip lengths are asymmetric. Figure 8 shows the variation in the critical Reynolds number with $[eqn]$ for $[eqn]$ , $[eqn]$ , and $[eqn]$ . The results indicate that $[eqn]$ increases with $[eqn]$ .

Figure 9 presents the variation in $[eqn]$ and $[eqn]$ with M for $[eqn]$ , $[eqn]$ , and $[eqn]$ . When $[eqn]$ , the critical Reynolds number exceeds $[eqn]$ , which is two orders of magnitude higher than the reference case with $[eqn]$ , $[eqn]$ , and no-slip boundaries ( $[eqn]$ ). This indicates that the combined effect of wall slip and porous media significantly raises the instability threshold of NSV fluid flow. Numerical results show that increasing M (which inherently couples porosity and permeability) leads to a significant rise in the critical Reynolds number $[eqn]$ . This indicates that Darcy dissipation in porous media is a powerful mechanism for suppressing disturbance growth. However, it should be noted that both M and $[eqn]$ are dimensionless parameters. Their relative importance depends heavily on specific physical applications and the chosen parameter ranges. Therefore, this comparison only highlights the qualitative differences between the two mechanisms rather than ranking their importance.

Finally, Figure 10 presents the neutral stability curves for $[eqn]$ , $[eqn]$ , and unequal slip lengths ( $[eqn]$ ). The numerical results indicate that increasing the slip length on either boundary fundamentally stabilizes the flow. For example, when $[eqn]$ and $[eqn]$ , $[eqn]$ = 81,024.5. Additionally, we found that for a fixed total slip length, a more symmetric slip distribution results in higher flow stability. For instance, the symmetric case with $[eqn]$ yields a higher critical Reynolds number ( $[eqn]$ ) than the asymmetric case with $[eqn]$ and $[eqn]$ .

3.5. Energy Budget Analysis

To elucidate the physical mechanisms governing the competition between destabilizing and stabilizing forces, we performed a global energy budget analysis.

Figure 11a,b present the normalized contributions of the Reynolds stress production ( $[eqn]$ ), viscous dissipation ( $[eqn]$ ), porous damping ( $[eqn]$ ), and Voigt regularization work ( $[eqn]$ ) to the time rate of change in the disturbance kinetic energy ( $[eqn]$ ). As shown in Figure 11a, the flow exhibits a distinct energy imbalance at a lower porous parameter of $[eqn]$ . The Reynolds stress production (orange bar, normalized to 1.0) extracts energy from the asymmetric mean shear. Its magnitude exceeds the sum of all dissipative mechanisms. Specifically, the viscous dissipation (blue bar, $[eqn]$ ), combined with the minor porous damping and Voigt terms, is insufficient to counteract the energy production. This results in a net energy surplus ( $[eqn]$ ). This result physically confirms the occurrence of linear instability driven by the asymmetric boundary conditions. In contrast, increasing the porous parameter to $[eqn]$ significantly alters the energy balance, as illustrated in Figure 11b. The stronger interaction between the fluid and the porous matrix modifies the perturbation mode shape. This leads to a substantial enhancement in viscous dissipation (magnitude increases to $[eqn]$ ). Consequently, the total dissipation overcomes the energy production, resulting in a net energy decay ( $[eqn]$ ). This demonstrates that the porous medium acts as a critical stabilizing factor. When M is sufficiently large, it suppresses the instability induced by the slip. It is noteworthy that the Voigt term consistently makes a negative contribution in both stable and unstable regimes, despite its small magnitude (order of $[eqn]$ ). This confirms that the Voigt term functions primarily as an auxiliary energy sink (dissipative regularization) rather than an elastic energy source, regardless of the flow’s stability status.

To uncover the physical origin of the destabilizing effect observed in the neutral curves, we examine the spatial distribution of the energy production term, $[eqn]$ . This term represents the rate of energy transfer from the mean flow to the disturbance.

Figure 12 compares the production profiles between the symmetric case ( $[eqn]$ ) and the asymmetric case ( $[eqn]$ ) at $[eqn]$ = 11,000 and $[eqn]$ . In the symmetric configuration (blue solid line), the energy production exhibits a balanced, bi-modal distribution. Two equal peaks are located near the critical layers at both walls. This symmetry indicates that the destabilizing potential is evenly distributed across the channel. However, the introduction of slip asymmetry (red dashed line) drastically alters this balance. The profile becomes highly skewed. It is characterized by a dominant peak near the upper wall ( $[eqn]$ ), where the slip length is smaller ( $[eqn]$ ). Fluid tends to flow more easily over the wall with higher slip ( $[eqn]$ ). This causes the maximum velocity of the base flow to shift towards the lower wall. Consequently, to satisfy mass conservation and boundary conditions, the mean velocity gradient (shear rate, $[eqn]$ ) steepens significantly near the opposite wall (the upper wall). As evidenced in Figure 12, the peak production in the asymmetric case significantly exceeds the maximum value observed in the symmetric case. This localized surge acts as a powerful source of instability. It injects energy into the disturbance field at a rate that the global viscous and Voigt dissipation mechanisms cannot effectively counteract. This triggers the onset of instability at lower Reynolds numbers, as predicted by the neutral curves.

As shown in the Figure 13, the Voigt regularization work $[eqn]$ is not uniformly distributed across the channel; instead, it is highly concentrated near the walls, targeting the regions with steep velocity gradients. Consistent with the mathematical definition of the Voigt term, as the parameter $[eqn]$ increases, the amplitude of this local regularization work grows significantly, resulting in sharper and more pronounced dissipative peaks near the boundaries (e.g., for $[eqn]$ ). As $[eqn]$ increases, the Voigt term acts as a targeted energy sink, inherently magnifying the localized dissipation precisely at the near-wall regions where the asymmetric slip induces the maximum shear and Reynolds stress production. Ultimately, this enhanced, localized energy sink becomes strong enough to fully suppress the instability driven by the asymmetric shear, stabilizing the global flow. Crucially, as shown in the global energy budget analysis (Figure 11), the net integration of this profile yields a negative value. This confirms that despite local fluctuations, the global effect of the localized Voigt term is to dissipate kinetic energy. It effectively suppresses wall-bounded disturbances, thereby stabilizing the flow.

4. Conclusions

This study systematically investigates the linear stability of Navier–Stokes–Voigt (NSV) fluid flow in porous media under general asymmetric slip boundary conditions. By extending the classical analysis of symmetric slip boundaries to asymmetric ones, and supported by both numerical evidence and energy budget analysis, the core conclusions are summarized as follows:

First, key parameters demonstrate a consistent stabilizing effect on the flow. Increasing the Voigt parameter $[eqn]$ or the porous medium parameter M significantly raises the critical Reynolds number $[eqn]$ . For a high-permeability medium ( $[eqn]$ ), the stability threshold exceeds $[eqn]$ . This is two orders of magnitude higher than the critical value for Newtonian fluids under no-slip conditions ( $[eqn]$ ). Notably, the flow becomes unconditionally stable when $[eqn]$ . Furthermore, comparing the stability characteristics of the NSV model with recent non-Newtonian models provides new insights. For instance, Sun [34] demonstrated a specific phenomenon in certain non-Newtonian fluids, such as those governed by second-order nonlinear spatial viscosity terms. In these fluids, increasing the non-Newtonian parameter beyond a critical threshold causes entire eigenvalue branches to enter the unstable region. This triggers extreme flow complexities and turbulence. In stark contrast, our eigenspectrum (Figure 6b) and energy analyses (Figure 11) confirm that the NSV fluid behaves in a completely opposite manner. Applying an extremely large $[eqn]$ does not induce turbulence. Instead, it unconditionally maintains flow stability. Energy analysis reveals the underlying physical mechanism. The Voigt term consistently acts as a dissipative energy sink (a regularization term) rather than an elastic energy source. Working together with the Darcy drag from the porous medium, it effectively suppresses disturbances near the wall.

Second, an important finding of this study is the impact of wall slip on stability. Introducing or increasing wall slip clearly stabilizes the overall flow. However, the system is relatively sensitive to how the slip is distributed. It must be emphasized that independently increasing the slip length on either wall enhances fluid stability. For example, in Figure 10, the neutral curves shift to the right when $[eqn]$ is fixed and $[eqn]$ is increased. However, our research also offers an optimal solution for a practical engineering problem: how to maximize flow stability when applying a fixed amount of superhydrophobic coating to a microchannel. A fixed amount of coating material represents a constant total slip in our model (i.e., $[eqn]$ ). When evaluating the slip distribution under this premise, symmetric slip provides the optimal stability margin. For instance, as also shown in Figure 10, for a total slip length of $[eqn]$ , the critical Reynolds number $[eqn]$ for the symmetric distribution ( $[eqn]$ ) is significantly higher than that for the asymmetric counterpart ( $[eqn]$ ). Energy analysis clarifies the physical root of this phenomenon. Asymmetry distorts the base flow shear. This creates a localized surge in Reynolds stress energy production. Consequently, any deviation from symmetry prevents the system from fully maximizing the stabilizing potential of a given slip budget. Thus, instability is triggered earlier than in a symmetric configuration with the exact same total slip. Furthermore, as the degree of asymmetry increases, the trend of stability reduction gradually flattens. This indicates a saturation of this destabilizing effect.

Furthermore, as the degree of asymmetry increases, the trend of stability reduction gradually flattens, indicating a saturation of the destabilizing effect. Third, from an engineering perspective, these findings establish theoretical bounds for manufacturing tolerances in superhydrophobic channels. While uniform coating is ideal, our results indicate that performance degradation caused by unavoidable manufacturing asymmetry can be compensated. This is achieved by slightly increasing the porous medium parameter (M) or the Voigt regularization parameter ( $[eqn]$ ). This offers a robust design strategy for microfluidic chips and hemodialyzers, ensuring flow stability even when coating uniformity is imperfect. Fourth, we clarified the applicability boundaries of the simplified analogy between asymmetric slip flow and Poiseuille–Couette flow. This analogy is valid only for small slip lengths ( $[eqn]$ , first-order approximation). For larger slip lengths ( $[eqn]$ ), higher-order terms become non-negligible, making the analogy suitable only for qualitative verification rather than quantitative prediction. Additionally, the numerical scheme demonstrates high robustness, achieving high accuracy (error $[eqn]$ ) with $[eqn]$ collocation points. However, for the large slip regime ( $[eqn]$ ), we recommend using $[eqn]$ to address numerical sensitivity.

Finally, regarding contributions and future outlook, this study advances existing research by moving beyond the idealized symmetric assumption of Shankar and Shivakumara [23]. It also extends the analytical scope of Baranovskii [18,19] by incorporating stability criteria in porous media. Future research should extend the model to full second-grade fluid or turbulent flow systems. This will allow for exploring the nonlinear coupling between slip asymmetry and viscoelasticity, which should be validated through experimental measurements.