Abstract

In this work, we study the linear temporal stability of a unidirectional Poiseuille blood flow under the influence of the magnetic parameter and the Reynolds number. Blood is modeled as a spatially inhomogeneous fluid, with its viscosity dependent on the red blood cell concentration. We focus on small blood vessels, such as the terminal branches of arteries, arterioles, and small veins. For the basic flow, assumed to be laminar and two-dimensional, the basic velocity profiles were obtained numerically, including the magnetic effect, using the BVP4C method in Matlab. Stability analysis was performed using classical linear temporal analysis in terms of normal modes, leading to a fourth-order eigenvalue problem solved numerically using the Chebyshev spectral collocation method. The results indicate that the flow is unconditionally unstable. However, increasing the magnetic parameter or the Reynolds number, coupled with decreasing red blood cell concentration, reduces the degree of flow instability. The stabilization of blood flow in microvessels is essential for protecting the vascular walls.

Keywords

Linear Temporal Stability, Blood Flow, Microvessels, BVP4C Technique, Spectral Collocation Method

Share and Cite:

1. Introduction

Newtonian fluid theory describes the mechanical behavior of many real fluids with good accuracy. However, numerous fluids cannot be adequately described by this theory and are generally referred to as non-Newtonian fluids. Blood is one such non-Newtonian fluid and is the focus of the present study. It is a complex and heterogeneous fluid composed of red blood cells suspended in a liquid medium known as plasma. Microcirculation refers to the flow of blood through the smallest blood vessels, including capillaries, arterioles, and venules. It plays a crucial role in maintaining the proper functioning of tissues and organs, and any disruption of microcirculation can lead to serious health complications.

Haynes [1] investigated the inhomogeneity of blood resulting from the non-uniform distribution of red blood cells across the cross-sections of vessels with diameters smaller than 300 μm. In other words, during microcirculation, red blood cells are not uniformly distributed over the vessel cross-section; instead, they tend to accumulate near the vessel axis, thereby forming a cell-rich core and a cell-free layer adjacent to the vessel wall. Moyers-Gonzalez et al. [2] modeled blood as a suspension of red blood cells (RBCs) in a liquid phase (plasma). During blood microcirculation, the volume fraction of cells (hematocrit) varies across the vessel cross-section, a phenomenon known as the Fahraeus-Lindqvist effect [3]. The inhomogeneous distribution of red blood cells has significant implications for the rheological properties of blood, particularly its viscosity, which depends on the spatial distribution of the cells. Fournier [4] showed that the Fahraeus-Lindqvist effect modifies blood viscosity as a function of vessel diameter. This variation in viscosity can induce instabilities in blood flow, making certain concentration profiles more stable than others [5]. Therefore, linear stability analysis is essential for understanding how these instabilities develop and for identifying the flow patterns that are most likely to persist under physiological conditions.

In this study, we investigate the linear stability of unidirectional flow in a channel filled with a fluid whose viscosity is a positive function of the hematocrit. When the hematocrit is uniform across the cross-section, the fluid exhibits the constitutive behavior of a Newtonian fluid. Conversely, when the hematocrit varies across the vessel cross-section, the flow becomes inhomogeneous. Relevant studies on stratified flows include those dealing with the so-called central annular flow, i.e., the parallel flow of two or more fluids with different viscosities. Blood flow in a microvessel, in connection with the Fahraeus-Lindqvist effect, is often modeled as a two-layer flow. Hickox [6] investigated the stability of this axisymmetric flow configuration, taking into account the effects of gravity and capillary forces acting at the interface between the two liquids. He showed that the steady Poiseuille flow of two immiscible fluids of different viscosities is unstable when the less viscous fluid is located at the center. Furthermore, Preziosi et al. [7] investigated more broadly the linear stability of the same basic flow by varying the viscosities, the volume ratios of the two fluids, and the Reynolds numbers. They showed that the flow is generally unstable, except when the less viscous fluid occupies the annular region, is sufficiently thin, and the Reynolds number lies within a limited range that depends on the fluid parameters. Anand and Rajagopal [5] demonstrated that inhomogeneous fluids, whose properties vary slightly from their mean values, can exhibit flow responses that differ by more than an order of magnitude from those associated with homogeneous fluids. Mohan Anand et al. [8] performed a linear stability analysis of a steady, fully developed flow of a shear-thinning fluid in a long cylindrical pipe. Using the shooting method, they found that all the computed eigenvalues are negative within the investigated Reynolds number range, indicating that the flow remains stable throughout this range. Miguel Moyers-Gonzalez et al. [2] applied the constitutive model originally proposed by Fang and Owens [9] and later developed by Owens [10] to steady, axisymmetric flow in rigid-walled tubes to describe non-homogeneous flows of healthy human blood. Their model accurately captures stress-induced cell migration in narrow tubes and predicts the Fahraeus-Lindqvist effect, according to which the apparent viscosity of healthy blood decreases with increasing tube diameter in sufficiently small vessels. Their numerical results show that this phenomenon is caused by the formation of a cell-depleted plasma layer near the vessel walls, which leads to a reduction in the tube hematocrit. Michela Ascolese et al. [3] showed that the marginal layer, although leading to a substantial reduction in flow resistance and an increase in discharge, does not decrease the energy dissipation rate. This result was obtained by considering six rheological models relating blood viscosity to hematocrit (the volume fraction occupied by erythrocytes). Lorenzo Fusi [11] studied the two-dimensional flow of a non-homogeneous incompressible fluid in a channel with a low aspect ratio. He extended the model to the case of a channel with variable thickness and compared his simulations with those obtained by Massoudi et al. [12]. Lorenzo Fusi et al. [13] studied the flow of a pressure-driven, inhomogeneous, incompressible thin film whose viscosity depends on its density. They showed that it is possible to determine an analytical solution to the problem when the boundary data are small perturbations of the homogeneous case. They then used this analytical solution to validate their numerical scheme. In [14], the authors illustrated the applications of mathematics to various physiological and artificial processes involving blood circulation, including hemorheology, microcirculation, coagulation, renal filtration, and dialysis, while providing a historical overview of each topic. They employed mathematical models to simulate processes occurring naturally in blood and to predict the effects of dysfunctions, such as coagulation disorders and renal failure, as well as the effects of therapies, to improve treatments. Khan A. et al. [15] presented a numerical study of the oscillatory motion of an Oldroyd-B fluid in a uniform magnetic field flowing through a small circular tube. First, they derived the orientation stress tensor by considering the Brownian force. Next, this tensor was incorporated into the Oldroyd-B model by considering Hookean dumbbells. The Oldroyd-B model was then reformulated by coupling it with the momentum equation and the total stress tensor. Finally, numerical simulations were performed to analyze the orientation stress tensor in the tube, showing that its effect is significant even when the Brownian force is sufficiently weak. According to Papalexandris [16] and Varsakelis et al. [17], the constitutive modeling of non-Brownian particle suspensions is considerably more complex. Furthermore, particle deformability must also be taken into account when modeling blood. Indeed, Goldsmith et al. [18] and Moyers-Gonzalez et al. [19] showed that, for Stokes flow in tube suspensions, particle deformability strongly influences the flow behavior. They described the steady Poiseuille flow of blood in a small tube using a two-layer fluid model consisting of an outer plasma layer and an inner core in which blood is treated as a suspension of rouleaux of different sizes represented by deformable dumbbells (Moyers-Gonzalez et al. [2]). Consequently, the total stress is viscoelastic and consists of a small Newtonian contribution from the plasma and an elastic contribution from the red blood cells. Such an approach was employed by Dimakopoulos et al. [20] to simulate blood flow in a stenosed vessel. Gwennou Coupier et al. [21] experimentally and numerically investigated the non-inertial transverse migration of a vesicle in a confined Poiseuille flow. They observed that the combined effects of the vessel walls and the curvature of the velocity profile induce migration toward the centerline of the channel.

In [22], L. Fusi and A. Farina studied the linear stability of unidirectional Poiseuille blood flow, modeling the fluid as spatially inhomogeneous, with viscosity depending on the concentration of red blood cells (RBCs). Their work focused on small vessels—such as terminal arterial branches, arterioles, or venules—where the inhomogeneity arises from the non-uniform distribution of red blood cells across the vessel’s cross-section. They found that distributions in which red blood cells are more concentrated near the vessel center are more “stable” than those in which red blood cells accumulate toward the vessel walls.

The objective and novelty of this work, compared to [22], lie in investigating the effect of a magnetic field and the Reynolds number on blood microcirculation through a linear temporal stability analysis, while—as in [22]—considering a hematocrit profile that depends only on the transverse coordinate, i.e., stratified flows.

The effect of the magnetic field on flow stability was investigated by analyzing the linear temporal stability of a viscous, incompressible, and electrically conductive fluid forming a dynamic laminar boundary layer over an impermeable horizontal flat magnetic plate, as presented by Lihonou et al. [23]. John et al. [24] demonstrated that increasing the Reynolds number, through either a larger cylinder radius or a greater sweep angle, in the presence of wall suction, can stabilize the boundary layer at the leading edges of swept wings. This counterintuitive mechanism challenges conventional expectations regarding flow transition, showing that increasing certain parameters may, under specific conditions, have a stabilizing effect on the flow. Lebbal [25] showed that, for pulsatile flow, increasing the Reynolds number can temporally stabilize the eTWF modes.

To achieve our objective, we consider two classes of hematocrit profiles: one that increases monotonically across the vessel cross-section and another that decreases monotonically. We focus on flows at low Reynolds numbers because our primary interest lies in microcirculation, namely arteries, terminal branches, arterioles, and venules (i.e., vessels with diameters ranging from 0.1 to 0.6 mm), where red blood cells are not uniformly distributed across the cross-section. To keep our analysis as general as possible, we employ three empirical laws relating blood viscosity to hematocrit. For each law, we determine the base velocity field corresponding to a prescribed transverse hematocrit profile. In particular, we examine two idealized cases representing the main hematocrit distributions: profiles that increase toward the vessel walls and profiles that decrease toward the vessel walls. Next, following the classical modal analysis of infinitesimal perturbations, we superimpose a small perturbation on the base flow and investigate its linear stability. The normal-mode formulation leads to a fourth-order eigenvalue problem, which is solved numerically using a Chebyshev polynomial method.

The remainder of this paper is organized as follows. In Section 2, we formulate the model for pressure-driven flow between two parallel plates. Section 3 presents the linear stability analysis, while Section 4 describes the numerical method used to solve the eigenvalue problem. The results and their discussion are presented in Section 5. Finally, concluding remarks are given in Section 6.

2. The Basic Flow

We consider a mechanically incompressible flow in a horizontal channel formed by two stationary parallel plates of length $L^o$ , located at $y^o=-H^o$ and $y^o=+H^o$ , giving a channel height of $2H^o$ . The flow is driven by a prescribed pressure gradient $G^o$ . A transverse magnetic field $B^o$ , directed toward the positive $y^o$ -axis, is applied across the channel. We exclude the cylindrical geometry, although it is physically more relevant, because it leads to a more involved numerical treatment. Likewise, the drag-driven flow case (i.e., flow induced by the motion of a plate in the absence of a pressure gradient) is not considered, as it is less relevant to blood flow applications. The present study is therefore restricted to a two-dimensional configuration.

We chose to model the microvessels as a two-dimensional parallel-plate channel because this geometry eliminates variations along the third dimension, radically simplifying the Navier-Stokes equations to yield exact solutions for the velocity profile. It facilitates the assessment of wall shear stresses. However, this choice represents a limitation when extrapolating the results to actual cylindrical vessels.

The superscript ${(.)}^o$ denotes dimensional quantities or variables. We denote by

$v^o=u^o{e}_x+v^o{e}_y$ (1)

the velocity field and we introduce the viscosity ${\mu }^o\left(\varphi \right)={\mu }_p^o\mu \left(\varphi \right)$ , where ${\mu }_p^o$ is a reference viscosity and $\mu \left(\varphi \right)$ is dimensionless function. Assuming that the Cauchy stress tensor is given by ${\mathbb{T}}^o=-p^o\mathbb{I}+2{\mu }^o\left(\varphi \right){\mathbb{D}}^o$ , the mathematical formulation of the problem reads

$\frac{\partial \varphi }{\partial t^o}+v^o\cdot {\nabla }^o\varphi =0,$ (2)

${\nabla }^o\cdot v^o=0,$ (3)

${\rho }^o\left(\frac{\partial v^o}{\partial t^o}+\left(v^o\cdot {\nabla }^o\right)v^o\right)=-{\nabla }^o{p}^o+{\mu }_p^o{\nabla }^o\cdot \left(2\mu \left(\varphi \right){\mathbb{D}}^o\right)+J^o\wedge B,$ (4)

${\nabla }^o\wedge B={\mu }_e\left(J^o+{\epsilon }_{e}E\right),$ (5)

${\nabla }^o\wedge E=-\frac{\partial B}{\partial t^o},$ (6)

${\nabla }^o\cdot B=0,$ (7)

${\nabla }^o\cdot E=0,$ (8)

${\nabla }^o\cdot J^o=0,$ (9)

Equation (2) is a simple advection equation. Equation (3) and Equation (4) are continuity and Newton’s second law respectively. Equations (5) - (9) are Ampere’s law, Faraday’s law, Maxwell’s law and Gauss law equations respectively, with

$J^o=\lambda \left(E+V^o\wedge B\right),$ (10)

where ${\rho }^o$ is the uniform fluid density, $E$ the electric field, $J^o$ the current density vector, ${\mu }_e$ the magnetic permeability, ${\epsilon }_e$ the absolute permittivity of the fluid, $t^o$ the time, $p^o$ the pressure, and $\lambda$ is the Stefan-Boltzmann constant.

We also considered the following

$B=\left(0,B_0,0\right),$ (11)

$E=\left(E_x,E_y,E_z\right),$ (12)

$J^o=\left(J_x^o,0,J_z^o\right),$ (13)

where $B_0$ is a constant. We assumed that no applied polarization voltage exists (i.e., $E=0$ ). Then Equation (10) and Equation (13) give

$J^o=\lambda B_0\left(-w^o,0,u^o\right)$ (14)

We rescale the problem with $x=\frac{x^o}{L^o}$ , $y=\frac{y^o}{L^o}$ ,  $H=\frac{H^o}{L^o}$ , $U=\frac{u^o}{U^o}$ , $V=\frac{v^o}{U^o}$ , $p=\frac{p^o}{{\rho }^o{U}^{o^2}}$ ,  $G=\frac{G^o{L}^o}{{\rho }^o{U}^{o^2}}$ , where $U^o$ is the characteristic velocity, still to be selected. The systems (2) - (4) becomes

$\{\begin{array}{l}\frac{\partial \varphi }{\partial t}+v\cdot \nabla \varphi =0,\hfill \\ \nabla \cdot v=0,\hfill \\ \mathrm{Re}\left(\frac{\partial v}{\partial t}+\left(v\cdot \nabla \right)v\right)=-\mathrm{Re}\nabla p+\nabla \cdot \left(2\mu \left(\varphi \right)\mathbb{D}\right)+J\wedge B.\hfill \end{array}$ (15)

In Cartesian coordinates, we have:

$\{\begin{array}{l}\frac{\partial \varphi }{\partial t}+U\frac{\partial \varphi }{\partial x}+V\frac{\partial \varphi }{\partial y}=0\hfill \\ \frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0\hfill \\ \mathrm{Re}\left(\frac{\partial U}{\partial t}+U\frac{\partial U}{\partial x}+V\frac{\partial U}{\partial y}\right)=-\mathrm{Re}\frac{\partial p}{\partial x}+\frac{\partial }{\partial x}\left(2\mu \left(\varphi \right)\frac{\partial U}{\partial x}\right)+\frac{\partial }{\partial y}\left(\mu \left(\varphi \right)\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)\right)-MU\hfill \\ \mathrm{Re}\left(\frac{\partial V}{\partial t}+U\frac{\partial V}{\partial x}+V\frac{\partial V}{\partial y}\right)=-\mathrm{Re}\frac{\partial p}{\partial y}+\frac{\partial }{\partial y}\left(2\mu \left(\varphi \right)\frac{\partial V}{\partial y}\right)+\frac{\partial }{\partial x}\left(\mu \left(\varphi \right)\left(\frac{\partial U}{\partial y}+\frac{\partial V}{\partial x}\right)\right)\hfill \end{array}$ (16)

where $\mathrm{Re}=L^o{U}^o/{\mu }_p^o$ is the Reynolds number and $M=\lambda B_0^2{L}^{o^2}/{\mu }_p^o$ , the magnetic parameter. We then look for a basic flow of the type $U=U\left(y\right)$ , $V=0$ , $p=p_o\left(x\right)$ , $\varphi ={\varphi }_o\left(y\right)$ , satisfying the boundary conditions

$U\left(H\right)=0,\frac{\partial U}{\partial y}\left(0\right)=0,p_o\left(0\right)=p_{in},p_o\left(1\right)=p_{in}-G,$ (17)

where $p_{in}$ is the dimensionless inlet pressure.

we set:

$\mu \left(\varphi \right)={\mu }_o\left(\varphi \right).$

So the system becomes:

$\{\begin{array}{l}0=-\mathrm{Re}\frac{\partial p_o}{\partial x}+\frac{\partial }{\partial y}\left({\mu }_o\left(\varphi \right)\left(\frac{\partial U}{\partial y}\right)\right)-MU\hfill \\ 0=\frac{\partial }{\partial x}\left({\mu }_o\left(\varphi \right)\left(\frac{\partial U}{\partial y}\right)\right)\hfill \end{array}$ (18)

The base flow rate is evaluated in the lower half of the channel, the upper half being obtained by symmetry.

Considering the first equation of system (18), we obtain

$\frac{\partial }{\partial y}\left({\mu }_o\left(\varphi \right)\frac{\partial U\left(y\right)}{\partial y}\right)-MU=\mathrm{Re}\frac{\partial P_o}{\partial x}$

considered the following

$p_o\left(x\right)=-Gx+p_{in},$

this leads $\frac{\partial P_o}{\partial x}=-G$

We therefore obtain:

$\frac{\partial }{\partial y}\left({\mu }_o\left(\varphi \right)\frac{\partial U\left(y\right)}{\partial y}\right)-MU=-\mathrm{Re}G.$ (19)

We considered:

$\mathrm{Re}G{\displaystyle {\int }_0^H\frac{\zeta }{{\mu }_o\left(\zeta \right)}\text{d}\zeta }=1,$ (20)

Thus, the Equation (19) can be written in the form:

${\mu }_o{U}^{″}+{\dot{\mu }}_o{\varphi }^{\prime }{U}^{\prime }-MU+\frac{1}{{\displaystyle {\int }_0^H\frac{y}{{\mu }_o\left(y\right)}\text{d}y}}=0,$(21)

where

$\mathrm{Re}G=\frac{1}{{\displaystyle {\int }_0^H\frac{y}{{\mu }_o\left(y\right)}\text{d}y}}$ (22)

${\frac{\text{d}{\mu }_o\left(\varphi \right)}{\text{d}\varphi }|}_{\varphi ={\varphi }_o\left(y\right)}={\dot{\mu }}_o\left(y\right),$

$\frac{\text{d}\varphi }{\text{d}y}={\varphi }^{\prime },$

The boundary conditions are:

$U\left(\pm H\right)=0,U^{\prime }\left(0\right)=0.$ (23)

Furthermore, we take:

$\varphi \left(y\right)=\frac{{\varphi }_M{\varphi }_m}{\left(1-y^2\right){\varphi }_m+y^2{\varphi }_M}$ (24)

In (24) ${\varphi }_m=0.1,{\varphi }_M=0.7$ are the minimum and maximum values attained by the hematocrit .

Thus, the hematocrit functions become:

$\varphi \left(y\right)=\frac{0.07}{\left(0.1+0.6y^2\right)}$

${\varphi }^{\prime }\left(y\right)=\frac{-0.084y}{{\left(0.1+0.6y^2\right)}^2}$

Such a choice of ${\varphi }_o\left(y\right)$ does not correspond to any real physiological situation; however, it provides a regular symmetric function bounded between two reasonable hematocrit values and attaining its maximum at ($y=0$ ). We recall that the hematocrit ${\varphi }_o$ varies locally, thereby affecting the viscous stress, while the fluid density remains unchanged.

In the literature, numerous empirical formulas have been proposed to relate viscosity to hematocrit (see, for instance, Fournier [4] and Hund et al. [26]). All these formulas share the common feature that ${\mu }_o$ is an increasing function of $\varphi$ . In the present work, we consider some of the most commonly used expressions reported in the literature, namely those proposed by Hatschek [27], Cokelet [28], and Nubar [29]:

$\mu \left(\varphi \right)=\frac{1}{1-{\varphi }^{1/3}}$ , (Hatschek [27])(25)

either:

$\dot{\mu }=\frac{1}{3{\left({\varphi }^{1/3}-{\varphi }^{2/3}\right)}^2},$ (26)

$\mu \left(\varphi \right)=\frac{1}{{\left(1-\varphi \right)}^{2.5}}$ , (Cokelet [28])(27)

either:

$\dot{\mu }=\frac{2.5{\varphi }^{1.5}}{{\left(1-\varphi \right)}^5},$ (28)

$\mu \left(\varphi \right)=\frac{0.75}{0.75-\varphi }$ , (Nubar [29])(29)

either:

$\dot{\mu }=\frac{0.75}{{\left(0.75-\varphi \right)}^2}.$ (30)

3. Linear Stability Analysis

To carry out a linear stability analysis of the basic flow (19), we introduce the following perturbed solution:

$\left[u-U\left(y\right),v,p-p_o\left(x\right),\varphi -{\varphi }_o\left(y\right)\right]=\left[\widehat{u}\left(y\right),\widehat{v}\left(y\right),\widehat{p}\left(y\right),\widehat{\varphi }\left(y\right)\right]{\text{e}}^{i\left(\alpha x-\omega t\right)}$(31)

e.i.

$\begin{array}{l}u\left(x,y,t\right)=U\left(y\right)+\widehat{u}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)},\\ v\left(x,y,t\right)=\widehat{v}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)},\\ p\left(x,y,t\right)=p_o\left(x\right)+\widehat{p}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)},\\ \varphi \left(x,y,t\right)={\varphi }_o\left(y\right)+\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)},\end{array}$

where $\left|\widehat{\cdot }\right|\ll 1$ , $\alpha \in ℝ$ is the wave number and $\omega \in ℂ$ is the frequency. Introducing

$c=\frac{\omega }{\alpha },\text{with}c\in ℂ,$

the perturbation phase can be rewritten as

$i\left(\alpha x-\omega t\right)=i\alpha \left(x-ct\right).$

Setting

${\frac{\text{d}\mu \left(\varphi \right)}{\text{d}\varphi }|}_{\varphi ={\varphi }_o\left(y\right)}={\dot{\mu }}_o\left(y\right),$(32)

By expanding the viscosity function in the vicinity of ${\varphi }_o$ up to the linear term, one successively obtains:

$\mu \left(\varphi \left(x,y,t\right)\right)=\mu \left({\varphi }_o\right)+{\frac{\left(\varphi -{\varphi }_o\right)}{1!}\frac{\text{d}\mu }{\text{d}\varphi }|}_{\varphi ={\varphi }_o}$

$\mu \left({\varphi }_o\right)={\mu }_o\left(y\right);{\frac{\text{d}\mu }{\text{d}\varphi }|}_{\varphi ={\varphi }_o}={\dot{\mu }}_o$ and $\varphi -{\varphi }_o=\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}$

Then, the viscosity function is expanded up to the linear term as follows:

$\mu \left(\varphi \left(x,y,t\right)\right)={\mu }_o\left(y\right)+{\dot{\mu }}_o\left(y\right)\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)},$(33)

By substituting Eqs. (31) into (15)₁, we have:

$\begin{array}{l}\frac{\partial }{\partial t}\left[\widehat{\varphi }{\text{e}}^{i\left(\alpha x-\omega t\right)}\right]+\left[U\left(y\right)+\widehat{u}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right]\frac{\partial }{\partial x}\left[{\varphi }_o\left(y\right)+\widehat{\varphi }{\text{e}}^{i\left(\alpha x-\omega t\right)}\right]\\ +\left(\widehat{v}{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)\frac{\partial }{\partial y}\left[{\varphi }_o\left(y\right)+\widehat{\varphi }{\text{e}}^{i\left(\alpha x-\omega t\right)}\right]=0.\end{array}$

Since $\frac{\partial \varphi }{\partial t}=\frac{\partial }{\partial t}\left[{\varphi }_o\left(y\right)+\widehat{\varphi }{\text{e}}^{i\left(\alpha x-\omega t\right)}\right]=\frac{\partial {\varphi }_o\left(y\right)}{\partial t}+\frac{\partial }{\partial t}\left(\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right)$ and $\frac{\partial {\varphi }_o\left(y\right)}{\partial t}=0$ , then

$\frac{\partial \varphi }{\partial t}=-i\omega \widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}$

Similarly, $\frac{\partial \varphi }{\partial x}=\frac{\partial }{\partial x}\left[{\varphi }_o\left(y\right)+\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right]=\frac{\partial \varphi \left(y\right)}{\partial x}+\frac{\partial }{\partial x}\left[\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right]$ with $\frac{\partial {\varphi }_o\left(y\right)}{\partial x}=0$ . Then

$\frac{\partial \varphi }{\partial x}=i\alpha \widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}$

At the same time, we have: $\frac{\partial \varphi }{\partial y}=\frac{\text{d}{\varphi }_o\left(y\right)}{\text{d}y}+\frac{\partial }{\partial y}\left[\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right]=\frac{\text{d}{\varphi }_o\left(y\right)}{\text{d}y}+\frac{\text{d}\widehat{\varphi }\left(y\right)}{\text{d}y}{\text{e}}^{i\left(\alpha x-\omega t\right)}$, e.i.:

$\frac{\partial \varphi }{\partial y}={{\varphi }^{\prime }}_o\left(y\right)+\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}$

Equation (15)₁ successively becomes:

$\begin{array}{l}i\omega \widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+\left[U\left(y\right)+\widehat{u}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right]\left[i\alpha \widehat{\varphi }\left(y\right){\text{e}}^{i\left(ax-\omega t\right)}\right]\\ +\widehat{v}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\left[{{\varphi }^{\prime }}_o\left(y\right)+{{\widehat{\varphi }}^{\prime }}_o\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right]=0\\ i\omega \widehat{\varphi }\left(y\right)+i\alpha \widehat{\varphi }\left(y\right)U\left(y\right)+i\alpha \widehat{u}\left(y\right)\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+\widehat{v}\left(y\right){{\varphi }^{\prime }}_o\left(y\right)\\ +\widehat{v}\left(y\right){{\widehat{\varphi }}^{\prime }}_o\left(y\right){\text{e}}^{i\left(ax-\omega t\right)}=0\\ i\left(\alpha U-\omega \right)\widehat{\varphi }\left(y\right)+\widehat{v}\left(y\right){{\varphi }^{\prime }}_o\left(y\right)+i\alpha \widehat{u}\left(y\right){\widehat{\varphi }}^{\prime }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+\widehat{v}\left(y\right){\widehat{\varphi }}^{\prime }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}=0\end{array}$

Neglecting the nonlinear terms, we obtain

$i\left(\alpha U-\omega \right)\widehat{\varphi }+\widehat{v}{{\varphi }^{\prime }}_o=0$ (34)

By substituting Eqs. (31) into (15)₂, we have:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=i\alpha \widehat{u}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+\frac{\text{d}\widehat{v}\left(y\right)}{\text{d}y}{\text{e}}^{i\left(\alpha x-\omega t\right)}=0$

with $\frac{\partial U\left(y\right)}{\partial x}=0$

Equation (15)₂ successively becomes:

$\begin{array}{l}i\alpha \widehat{u}\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+{\widehat{v}}^{\prime }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)}=0\hfill \end{array}$

${\widehat{v}}^{\prime }+i\alpha \widehat{u}=0,$ (35)

where

${(\cdot )}^{\prime }=\frac{\text{d}(\cdot )}{\text{d}y}.$

Substituting (31) and (33) into (15)₃ we obtain

$\mathrm{Re}\left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)=-\mathrm{Re}\frac{\partial p}{\partial x}+\frac{\partial }{\partial x}\left(2\mu \left(\varphi \right)\frac{\partial u}{\partial x}\right)+\frac{\partial }{\partial y}\left[\mu \left(\varphi \right)\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]-Mu$

and

$\mathrm{Re}\left(\frac{\partial v}{\partial t}+u\frac{\partial \nu }{\partial x}+v\frac{\partial v}{\partial y}\right)=-\mathrm{Re}\frac{\partial p}{\partial y}+\frac{\partial }{\partial y}\left[2\mu \left(\varphi \right)\frac{\partial v}{\partial y}\right]+\frac{\partial }{\partial x}\left[\mu \left(\varphi \right)\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]$

We can calculate:

$\frac{\partial \mu \left(\varphi \right)}{\partial x}=i\alpha {\dot{\mu }}_o\widehat{\varphi }\left(y\right){\text{e}}^{i\left(\alpha x-\omega t\right)},$

$\frac{\partial }{\partial x}\left[2\mu \left(\varphi \right)\frac{\partial u}{\partial x}\right]=2\frac{\partial \mu \left(\varphi \right)}{\partial x}\frac{\partial u}{\partial x}+2\mu \left(\varphi \right)\frac{{\partial }^{2}u}{\partial x^2}$

$\frac{\partial }{\partial x}\left[2\mu \left(\varphi \right)\frac{\partial u}{\partial x}\right]=-2{\alpha }^2{\dot{\mu }}_o\widehat{\varphi }\widehat{u}{\text{e}}^{2i\left(\alpha x-\omega t\right)}-2{\alpha }^2{\mu }_o\widehat{u}{\text{e}}^{i\left(\alpha x-\omega t\right)}-2{\alpha }^2{\dot{\mu }}_o\widehat{u}\widehat{\varphi }{\text{e}}^{2i\left(\alpha x-\omega t\right)},$

$\begin{array}{l}\frac{\partial }{\partial y}\left[\mu \left(\varphi \right)\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]\\ =\frac{\partial }{\partial y}\left[{\mu }_o{U}^{\prime }+\dot{\mu }{U}^{\prime }\widehat{\varphi }{\text{e}}^{i\left(\alpha x-\omega t\right)}+{\mu }_o\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+{\dot{\mu }}_o\widehat{\varphi }\left({\widehat{\mu }}^{\prime }+i\alpha \widehat{\nu }\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right]\end{array}$

$\frac{\partial }{\partial y}\left[2\mu \left(\varphi \right)\frac{\partial v}{\partial y}\right]=2\frac{\partial \mu \left(\varphi \right)}{\partial y}\frac{\partial v}{\partial y}+2\mu \left(\varphi \right)\frac{{\partial }^{2}v}{\partial y^2}$

$\begin{array}{c}\frac{\partial }{\partial y}\left[2\mu \left(\varphi \right)\frac{\partial v}{\partial y}\right]=2U^{\prime }{\widehat{v}}^{\prime }{\text{e}}^{i\left(\alpha x-\omega t\right)}+2{{\dot{\mu }}^{\prime }}_o{\widehat{v}}^{\prime }\widehat{\varphi }{\text{e}}^{2i\left(\alpha x-\omega t\right)}+2{\dot{\mu }}_o{\widehat{v}}^{\prime }{\widehat{\varphi }}^{\prime }{\text{e}}^{2i\left(\alpha x-\omega t\right)}\\ +2{\mu }_o{\widehat{v}}^{″}{\text{e}}^{i\left(\alpha x-\omega t\right)}+2{\dot{\mu }}_o\widehat{\varphi }{\widehat{v}}^{″}{\text{e}}^{i\left(\alpha x-\omega t\right)}\end{array}$

$\begin{array}{c}\frac{\partial }{\partial y}\left[2\mu \left(\varphi \right)\frac{\partial v}{\partial y}\right]=\frac{\text{d}}{\text{d}y}\left(2{\mu }_o{\widehat{v}}^{\prime }\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+2{\dot{\mu }}_o{\widehat{v}}^{\prime }\widehat{\varphi }{\text{e}}^{2i\left(\alpha x-\omega t\right)}\\ +2{\dot{\mu }}_o{\widehat{v}}^{\prime }{\widehat{\varphi }}^{\prime }{\text{e}}^{2i\left(\alpha x-\omega t\right)}+2{\dot{\mu }}_o\widehat{\varphi }{\widehat{v}}^{″}{\text{e}}^{i\left(\alpha x-\omega t\right)}\end{array}$

$\begin{array}{c}\frac{\partial }{\partial x}\left[\mu \left(\varphi \right)\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right]=\frac{\partial }{\partial x}\left({\mu }_o{U}^{\prime }\right)+[{\dot{\mu }}_o\widehat{\varphi }{U}^{\prime }{\text{e}}^{i\left(\alpha x-\omega t\right)}+{\mu }_o\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\\ +{\dot{\mu }}_o\widehat{\varphi }\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right){\text{e}}^{i\left(\alpha x-\omega t\right)}]\end{array}$

Thus, along (x) and (y), the equation (15)₃ becomes:

$\begin{array}{l}\mathrm{Re}\left[-i\omega \widehat{u}{\text{e}}^{i\left(\alpha x-\omega t\right)}+\left({\mu }_o+\widehat{u}{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)\left(i\alpha \widehat{u}{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)+\left(\widehat{v}{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)\left({{\mu }^{\prime }}_o+{\widehat{u}}^{\prime }{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)\right]\\ =-\mathrm{Re}\left(\frac{\partial p_o\left(x\right)}{\partial x}+i\widehat{p}{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)-2{\alpha }^2{\dot{\mu }}_o\widehat{u}\widehat{\varphi }{\text{e}}^{2i\left(ax-\omega t\right)}-2a^2{\mu }_o\widehat{u}{\text{e}}^{i\left(\alpha x-\omega t\right)}-2{\alpha }^2{\dot{\mu }}_o\widehat{u}\widehat{\varphi }{\text{e}}^{i\left(\alpha x-\omega t\right)}\\ +\frac{\partial }{\partial y}\left[{\mu }_o{U}^{\prime }+{\dot{\mu }}_o{U}^{\prime }\widehat{\varphi }{\text{e}}^{i\left(\alpha x-\omega t\right)}+{\mu }_o\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+{\dot{\mu }}_o\widehat{\varphi }\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right){\text{e}}^{2i\left(\alpha x-\omega t\right)}\right]\\ -MU-M\widehat{u}{\text{e}}^{i\left(\alpha x-\omega t\right)}\end{array}$

and

$\begin{array}{l}\mathrm{Re}\left[-i\omega \widehat{v}{\text{e}}^{i\left(\alpha x-\omega t\right)}+\left(U+\widehat{u}{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)\left(i\alpha \widehat{v}{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)\left({\widehat{v}}^{\prime }{\text{e}}^{i\left(\alpha x-\omega t\right)}\right)\right]\\ =-\mathrm{Re}{\widehat{p}}^{\prime }{\text{e}}^{i\left(\alpha x-\omega t\right)}+\frac{\text{d}}{\text{d}y}\left(2{\mu }_o{\widehat{v}}^{\prime }\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+2{{\dot{\mu }}^{\prime }}_o{\widehat{v}}^{\prime }\widehat{\varphi }{\text{e}}^{2i\left(\alpha x-\omega t\right)}\\ +2{\dot{\mu }}_o{\widehat{v}}^{\prime }{\widehat{\varphi }}^{\prime }{\text{e}}^{2i\left(\alpha x-\omega t\right)}+2{\dot{\mu }}_o\widehat{\varphi }{\widehat{v}}^{″}{\text{e}}^{i\left(\alpha x-\omega t\right)}+\frac{\partial }{\partial x}\left({\mu }_o{U}^{\prime }\right)\\ +i\alpha \left[{\dot{\mu }}_o\widehat{\varphi }{U}^{\prime }{\text{e}}^{i\left(\alpha x-\omega t\right)}+{\mu }_o\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right){\text{e}}^{i\left(\alpha x-\omega t\right)}+{\dot{\mu }}_o\widehat{\varphi }\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right){\text{e}}^{i\left(\alpha x-\omega t\right)}\right].\end{array}$

Based on the base flow, we have:

$-\mathrm{Re}\frac{\partial p_o\left(x\right)}{\partial x}+\frac{\partial }{\partial y}\left({\mu }_o{U}^{\prime }\right)-MU=0$

and $\frac{\partial }{\partial x}\left({\mu }_o{U}^{\prime }\right)=0$

Furthermore, by neglecting the non-linear terms and simplifying the exponential function in these preceding equations of motion, we obtain:

$\mathrm{Re}\left(-i\omega \widehat{u}+i\alpha U\widehat{u}+U^{\prime }\widehat{v}\right)=-i\alpha \mathrm{Re}\widehat{p}-2{\alpha }^2{\mu }_o\widehat{u}+\frac{\text{d}}{\text{d}y}\left[{\mu }_o\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right)+{\dot{\mu }}_o\widehat{\varphi }{U}^{\prime }\right]-M\widehat{u}$ (36)

and

$\mathrm{Re}\left(-i\omega \widehat{v}+i\alpha U\widehat{v}\right)=-\mathrm{Re}{\widehat{p}}^{\prime }+i\alpha \left({\mu }_o\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right)+{\dot{\mu }}_o\widehat{\varphi }{U}^{\prime }\right)+\frac{\text{d}}{\text{d}y}\left[2{\mu }_o{\widehat{v}}^{\prime }\right]$(37)

Exploiting (35) and introducing

$Q\left(y\right)={\mu }_o\left(y\right)-1$ ,

$P\left(y\right)={\dot{\mu }}_o\left(y\right){{\varphi }^{\prime }}_o\left(y\right)U^{\prime }\left(y\right)$ ,

Equation (36) and Equation (37) can be rewritten as

$\begin{array}{l}\mathrm{Re}\left(-i\omega \widehat{u}+i\alpha U\widehat{u}+U^{\prime }\widehat{v}\right)\\ =-i\alpha \mathrm{Re}\widehat{p}-2{\alpha }^{2}Q\widehat{u}+\left(\frac{{\text{d}}^2}{\text{d}{y}^2}-{\alpha }^2\right)\widehat{u}+\frac{\text{d}}{\text{d}y}\left[Q\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right)+P\frac{\widehat{\varphi }}{{{\varphi }^{\prime }}_o}\right]-M\widehat{u}\end{array}$(38)

$\mathrm{Re}\left(-i\omega \widehat{v}+i\alpha U\widehat{v}\right)=-\mathrm{Re}{\widehat{p}}^{\prime }+\left(\frac{{\text{d}}^2}{\text{d}{y}^2}-{\alpha }^2\right)\widehat{v}+i\alpha \left[Q\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right)+P\frac{\widehat{\varphi }}{{{\varphi }^{\prime }}_o}\right]+\frac{\text{d}}{\text{d}y}\left[2Q{\widehat{v}}^{\prime }\right]$ (39)

Eliminating the pressure between Equation (38) and Equation (39) and using Equation (35), we obtain

$\begin{array}{l}\mathrm{Re}\left[i\left(\alpha U-\omega \right)\left({\widehat{u}}^{\prime }-i\alpha \widehat{v}\right)+U^{″}\widehat{v}\right]\\ =\left(\frac{{\text{d}}^2}{\text{d}{y}^2}-{\alpha }^2\right)\left({\widehat{u}}^{\prime }-i\alpha \widehat{v}\right)+\left(\frac{{\text{d}}^2}{\text{d}{y}^2}+{\alpha }^2\right)\left[Q\left({\widehat{u}}^{\prime }+i\alpha \widehat{v}\right)+P\frac{\widehat{\varphi }}{{{\varphi }^{\prime }}_o}\right]\\ -2i\alpha \frac{\text{d}}{\text{d}y}\left[Q\left({\widehat{v}}^{\prime }-i\alpha \widehat{u}\right)\right]-M{\widehat{u}}^{\prime }\end{array}$(40)

Introducing the new variable $f\left(y\right)$

$\widehat{\varphi }={{\varphi }^{\prime }}_{o}f$ ,

implying that $\widehat{\varphi }$ is symmetric with respect to $y=0$ , Equation (34) and Equation (35) can be rewritten as

$\widehat{v}=i\alpha \left(c-U\right)f,\widehat{u}=-\frac{\text{d}}{\text{d}y}\left[\left(c-U\right)f\right]$.

We also have:

${\widehat{v}}^{\prime }=i\alpha \frac{\text{d}}{\text{d}y}\left[\left(c-U\right)f\right],{\widehat{u}}^{\prime }=-\frac{{\text{d}}^2}{\text{d}{y}^2}\left[\left(c-U\right)f\right]$.

On substituting the above into Equation (40) we obtain the fourth order eigenvalue problem

$\begin{array}{l}i\alpha \mathrm{Re}\left\{\left(U-c\right)\left({\text{D}}^2-{\alpha }^2\right)-U^{″}\right\}\left[\left(U-c\right)f\right]\\ ={\left({\text{D}}^2-{\alpha }^2\right)}^2\left[\left(U-c\right)f\right]+\left({\text{D}}^2+{\alpha }^2\right)\left\{Q\left({\text{D}}^2+{\alpha }^2\right)\left[\left(U-c\right)f\right]+Pf\right\}\\ -4{\alpha }^2\text{D}\left\{Q\text{D}\left[\left(U-c\right)f\right]\right\}-M{\text{D}}^2\left[\left(U-c\right)f\right]\end{array}$ (41)

where

${\text{D}}^k=\frac{{\text{d}}^k}{\text{d}{y}^k}$

and whose boundary conditions are

${f|}_{\pm H}=0,{f^{\prime }|}_{\pm H}=0.$ (42)

System (41) and (42) provides the eigenvalues $c\in ℂ$ , and the relative eigenfunctions $f$ , that allows one to establish if the basic flow, corresponding to the prescribed ${\varphi }_o$ and to the selected ${\mu }_o\left(\varphi \right)$ , is linearly stable or not. In particular, when $\mathrm{Im}\left(c\right)>0$ the system is unstable.

Remark 1. We observe that when $Q=P=0$ , i.e. when $\mu =1$ , setting $g=\left(U-c\right)f$ , the Equation (41) become:

$i\alpha \mathrm{Re}\left\{\left(U-c\right)\left({\text{D}}^2-{\alpha }^2\right)-U^{″}\right\}g={\left({\text{D}}^2-{\alpha }^2\right)}^{2}g-M{\text{D}}^{2}g$ .

Furthermore, when the magnetic parameter (M) is set to zero, the classical Orr-Sommerfeld equation is recovered

$i\alpha \mathrm{Re}\left\{\left(U-c\right)\left({\text{D}}^2-{\alpha }^2\right)-U^{″}\right\}g={\left({\text{D}}^2-{\alpha }^2\right)}^{2}g$ .

The aim here is not to solve the classical Orr-Sommerfeld equation, but rather to solve the complex Equation (41); we will first transform this into an eigenvalue equation, enabling a direct analysis of blood flow stability, simplifying the calculation of perturbation waves, and allowing the use of appropriate numerical methods.

To expand Equation (41) by letting the operator $\text{D}$ act on its contents, let us set:

${\Gamma }_L=i\alpha \mathrm{Re}\left\{\left(U-c\right)\left({\text{D}}^2-{\alpha }^2\right)-U^{″}\right\}\left[\left(U-c\right)f\right]$ , and

$\begin{array}{c}{\Gamma }_R={\left({\text{D}}^2-{\alpha }^2\right)}^2\left[\left(U-c\right)f\right]+\left({\text{D}}^2+{\alpha }^2\right)\left\{Q\left({\text{D}}^2+{\alpha }^2\right)\left[\left(U-c\right)f\right]+Pf\right\}\\ -4{\alpha }^2\text{D}\left\{Q\text{D}\left[\left(U-c\right)f\right]\right\}-M{\text{D}}^2\left[\left(U-c\right)f\right]\end{array}$

Thus, we have:

$\begin{array}{c}{\Gamma }_L=\left[i\alpha \mathrm{Re}U{U}^{″}f+i\alpha \mathrm{Re}{U}^2{\text{D}}^{2}f+2i\alpha \mathrm{Re}U{U}^{\prime }\text{D}f-i{\alpha }^3\mathrm{Re}{U}^{2}f-i\alpha \mathrm{Re}U{U}^{″}f\right]\\ +c[-i\alpha \mathrm{Re}U{\text{D}}^{2}f+2i{\alpha }^3\mathrm{Re}Uf-i\alpha \mathrm{Re}{U}^{″}f-i\alpha \mathrm{Re}U{\text{D}}^{2}f-2i\mathrm{Re}{U}^{\prime }\text{D}f\\ +i\alpha \mathrm{Re}{U}^{″}f]+c^2[i\alpha \mathrm{Re}{\text{D}}^{2}f\end{array}$

and

$\begin{array}{c}{\Gamma }_R=-4{\alpha }^2{Q}^{\prime }{U}^{″}f-4{\alpha }^2{Q}^{\prime }U{\text{D}}^{2}f-8{\alpha }^2{Q}^{\prime }{U}^{\prime }\text{D}f+4c{\alpha }^2{Q}^{\prime }\text{D}f-4{\alpha }^{2}Q{U}^{\left(3\right)}f\\ -4{\alpha }^{2}Q{U}^{″}\text{D}f-4{\alpha }^{2}Q{U}^{\prime }\text{D}-4{\alpha }^{2}QU{\text{D}}^{3}f-8{\alpha }^{2}Q{U}^{″}\text{D}f-8{\alpha }^{2}Q{U}^{\prime }{\text{D}}^{2}f\\ +4a^{2}Qc{\text{D}}^{2}f+{\alpha }^{4}Uf-{\alpha }^{4}cf+U{\text{D}}^{4}f+4U^{\prime }{\text{D}}^{3}f+4U^{\left(3\right)}\text{D}f+U^{\left(4\right)}f\\ -c{\text{D}}^{4}f-2{\alpha }^2{U}^{″}f-2{\alpha }^{2}U{\text{D}}^{2}f-4{\alpha }^2{U}^{\prime }\text{D}f+2{\alpha }^{2}c{\text{D}}^{2}f+{\alpha }^{2}Q{U}^{″}f\\ +{\alpha }^{2}QUf+2{\alpha }^{2}Q{U}^{\prime }\text{D}f-c{\alpha }^{2}Q{\text{D}}^{2}f+{\alpha }^{4}QUf-c{\alpha }^{4}Qf+{\alpha }^{2}Pf+P^{″}f\\ +P{\text{D}}^{2}f+2P^{\prime }\text{D}f-c{\alpha }^2{Q}^{″}f-c-2c{\alpha }^2{Q}^{\prime }Df-c{Q}^{″}{\text{D}}^{2}f-cQ{\text{D}}^{4}f\\ -c{Q}^{\prime }{\text{D}}^{3}f+Q^{″}{U}^{″}f+Q{U}^{\left(4\right)}f+Q{U}^{″}{\text{D}}^{2}f+2Q^{\prime }{U}^{\left(3\right)}f+2QU+2Q^{\prime }{U}^{″}\text{D}f\\ +Q^{″}U{\text{D}}^{2}f+Q{U}^{″}{\text{D}}^{2}f+QU{\text{D}}^{4}f+2Q^{\prime }{U}^{\prime }{\text{D}}^{2}f+2Q{U}^{\prime }{\text{D}}^{3}f+2Q^{\prime }U{\text{D}}^{3}f\\ +2Q^{″}{U}^{\prime }+2Q{U}^{\prime }{\text{D}}^{3}f+4Q^{\prime }{U}^{″}\text{D}f+4Q{U}^{″}+4Q^{\prime }{U}^{\prime }{\text{D}}^{2}f+\alpha Q^{″}Uf\\ +a^{2}Q{U}^{″}f+{\alpha }^{2}QU{\text{D}}^{2}f+2{\alpha }^2{Q}^{\prime }{U}^{\prime }+2{\alpha }^{2}Q{U}^{\prime }\text{D}f+2{\alpha }^2{Q}^{\prime }U\text{D}f-M{U}^{″}f\\ -2M{U}^{\prime }\text{D}f-MU{\text{D}}^{2}f+cM{\text{D}}^{2}f\end{array}$

By taking

$\begin{array}{l}{\Gamma }_L-{\Gamma }_R=0\hfill \end{array}$

and then factoring the $c$ , Equation (41) can be rewritten as

${\Gamma }_{0}f+c{\Gamma }_{1}f+c^2{\Gamma }_{2}f=0,$ (43)

where ${\Gamma }_j$ are the differential operators

${\Gamma }_j={\displaystyle \sum _{k=0}^4{\Gamma }_{jk}\left(y\right){\text{D}}^k}.$ (44)

The coefficients ${\Gamma }_{jk}\left(y\right)$ are reported in Appendix.

So,

${\Gamma }_0={\Gamma }_{00}+{\Gamma }_{01}\text{D}+{\Gamma }_{02}{\text{D}}^2+{\Gamma }_{03}{\text{D}}^3+{\Gamma }_{04}{\text{D}}^4,$ (45)

${\Gamma }_1={\Gamma }_{10}+{\Gamma }_{11}\text{D}+{\Gamma }_{12}{\text{D}}^2+{\Gamma }_{13}{\text{D}}^3+{\Gamma }_{14}{\text{D}}^4,$ (46)

${\Gamma }_2={\Gamma }_{20}+{\Gamma }_{21}\text{D}+{\Gamma }_{22}{\text{D}}^2+{\Gamma }_{23}{\text{D}}^3+{\Gamma }_{24}{\text{D}}^4.$ (47)

4. Numerical Solution Method

4.1. Basic Flow Solution Method

For the numerical solution of the basic flow, we employed the bvp4c solver implemented in MATLAB, following the approach of Shampine et al. [30]. This method has also been used by T. Jamir and H. Konwar (2022) [31] and by Lihonou et al. (2025) [32] [33]. Therefore, the highly coupled nonlinear ordinary differential Equation (21), together with the boundary conditions (23), is solved by introducing the following variables:

$A1=\left(Ip\ast gp\right)/I$ ;

$A2=M/I$ ;

$A3=RG/I$ ;

$U=Z_1$ ; $U^{\prime }=Z_2$ ;