Novel Aspects of Cilia-Driven Flow of Viscoelastic Fluid through a Non-Darcy Medium under the Influence of an Induced Magnetic Field and Heat Transfer

Abstract

The spontaneous movement of natural motile cilia in the form of metachronal waves is responsible for fluid transport. These cilia, in particular, play important roles in locomotion, feeding, liquid pumping, and cell delivery. On the other hand, artificial cilia can be useful in lab-on-a-chip devices for manipulation processes. In this study, a novel model for the ciliated tapered channel in Sutterby fluid flow under the impact of an induced magnetic field and heat transport is proposed. The Darcy–Brinkman–Forchheimer law for porous media with a viscous dissipation function is considered. With the help of lubrication theory, the simplified non-linear form of the leading equation with cilia-oriented boundary conditions is achieved. The analytical results of differential equations are based on the topological perturbation approach. The numerical simulation is performed to elaborate on the physical interpretations of emerging parameters through computer software.

1. Introduction

Nature has created distinct mechanisms for producing fluid flow, most of which are related to animal propulsion, such as swimming or flying. Significant examples at large scales are the waving tails of fish and the flapping wings of birds. Flapping wings can be found in insects at smaller scales as well. A mechanism of fluid manipulation used in nature at tiny scales, usually sub-millimeter sizes, is that of flagella or cilia.

Cilia look like small hairs or flexible rods and are present in almost all mammalian groups. They are usually divided into two more significant kinds: “motile” cilia and “immotile” cilia. The two main functions of motile and immotile cilia are fluid flow generation and detection, respectively. The asymmetric movements of motile cilia in a concentrated fashion (back and forth) are instrumental in generating fluid flow. In this study, we focused on motile cilia, which perform many physiological functions inside the human body. For instance, these motile cilia are present in the fallopian tube, where they are very crucial for ovum transport; in the ventricles and spinal cord of the adult brain, where they help in the movement of polarized fluid, which is essential for neuronal migration and the circulation of cerebrospinal fluid; and also in the airway, where they are necessary for mucus clearance. In 1962, Sleigh [1] also elaborated on the structure of cilia and the diverse patterns of metachronal waves caused by the well-coordinated movement of cilia. In the literature, the sublayer and envelope models are the two effective approaches used to analyze fluid transport due to the systematic beating of cilia. The sublayer model is useful for the computation of the average or instantaneous velocity at any point outside or inside the cilia layer. Still, it is effective only for widely spaced cilia. Blake [2,3] was the first to take the initiative to consider the sublayer model for the infinite plan wall and also extend his work to two parallel walls. In contrast, the envelope model assumes an envelope of cilia tip profiles, in which the cilia are closely packed together and very effective at determining the velocity field outside the cilia layer. This model is very useful in a variety of biological systems. Shack and Lardner [4] used an envelope-model approach to investigate fluid transport in the ductus efferentes of the male reproductive tract due to the symplectic pattern of a metachronal wave. Akbar and Akhtar [5] recently computed mathematical results for the cilia-driven flow of the Phan–Thien–Tanner (PTT) fluid model inside a channel.

Heat transport is the most exciting phenomenon in cilia-mediated physiological systems. Heat-transfer effects may also substantially influence the viscoelastic fluid cilia motion across a porous material. Convective heat transfer may be induced by the temperature differential across the porous medium, which can impact the fluid flow and cilia motility. Moreover, heat transmission may modify the viscoelastic fluid’s rheological parameters, which can affect the cilia motility. Baetjer [6] investigated the impact of vapor pressure and temperature on the clearance rate of the cilia mucus. The rate of cilia beating, mucus viscosity, cilia density, temperature, and airway humidity are the key players in mucociliary transport [7]. Imran et al. [8] developed a new mathematical model for the heat transport of cilia-oriented nanofluids through a ciliated channel. The leading partial differential equations of the proposed mathematical model have been found after utilizing the concept of boundary-layer approximation. They also found the exact solution and represented it graphically. Further, their results are very effective in helping to understand the ductus efferentes of the human male reproductive tract. A recent article on Casson fluid and the constant temperature at the ciliated walls of the elliptic duct was presented by Fuzhang et al. [9]. They noticed that a Newtonian fluid moves faster than a non-Newtonian fluid in this elliptic duct with ciliated walls. Moreover, several recent studies about cilia-driven flow can be seen in [10,11,12].

The application of a magnetic field has been reported to cause ciliary motility. The interaction between the magnetic field and the microtubules inside the cilia causes magneto-mechanical coupling, which causes cilia motion. Both in vitro and in vivo tests have shown this phenomenon [13,14]. Microfluidic mixing and pharmaceutical fields are two uses of magnetic-field-induced cilia motility. In microfluidic mixing, magneto-mechanical coupling has been employed to create ciliary motion, resulting in fluid mixing in a microfluidic channel. This technology can potentially be used in quick and effective fluid mixing in various biological and chemical applications. Magnetic-field-induced cilia motility has been investigated as a strategy for targeted medication delivery in drug delivery. This application has the prospective to be used in treating a wide range of disorders, including cancer, where customized medication delivery to the afflicted region may considerably increase therapy effectiveness. Numerous scientists have investigated the potential uses of magnetic-field-induced cilia motility. Siddique et al. [15] gave the theoretical results of cilia-oriented flow in a tube. They chose a Newtonian fluid under the influence of MHD and later on extended their investigations, particularly for the non-Newtonian Casson fluid. Akbar et al. [16] discussed the creeping flow of copper–water nanofluids in a ciliated parallel plate channel. They included the effects of a high Reynolds number and explained their physical significance with the help of graphical results. Javid et al. [17] consider the ciliated micro-channel (asymmetric) and investigated the effects of thermal slip, porosity, magnetic field, and electro-osmosis on the nanofluids. They attained the exact solutions of mathematically modeled equations and produced three-dimensional graphs for the diverse emerging parameters of the considered problem for further discussion. Many researchers have utilized MHD [18,19,20,21] in various flows and geometries. On the other hand, biomimetic artificial cilia are critical players in lab-on-a-chip devices for fluid manipulation. They have several microfluidic applications, including particle manipulation, microsensors, bacterium collection, droplets, microfluidic pumping, and self-cleaning microfluidic mixing. Hanasoge et al. [22] elaborated on the motion of artificial magnetic cilia in a microchannel loop that is crucial for flow visualization. Toonder and Mayne [23] also utilized the microfluidic channel and discussed in detail the power of magnetic cilia in bacterium collection and fluid mixing.

Many transportation processes via porous media are important in mechanical industries, physiological processes, and engineering [24,25]. In the literature, different models such as Darcy, Darcy–Forchecimer, and Brinkman-extended Darcy are utilized by researchers for the study of porous media [26]. Akbar et al. [27] used the Darcy law for the creeping flow inside the slippery ciliated tube. Manzoor et al. [28] investigated the effective properties of the viscoelastic model (Jeffrey fluid) via a ciliated channel. They also included the effects of MHD, porous media, and convective heat transfer, and they employed the Adomian decomposition method (ADM) for mathematical solutions. Furthermore, Javid et al. [29] also considered a porous medium with Jeffrey fluids and found the numerical solution for the divergent ciliated channel with the aid of the BVP4C technique. Recently, Aich et al. [30] discussed the porosity effects of hybrid nanofluids through a cilia-oriented cylindrical tube. Moreover, advancement in the thermal and mechanical properties of human tissue that govern the biological process is due to the utilization of fundamental engineering principles in the study of many heat- and mass-transport applications in biology and medicine. Great interest has emerged in bioheat-transfer phenomena during the last twenty years, with a particular focus on therapeutic and diagnostic applications. Depending on advanced computational techniques, the establishment of mathematical models has increased the ability to study many types of the bioheat-transfer process [31,32,33,34,35]. Usman et al. [36] considered the Williamson fluid through a ciliated channel and utilized the concept of Darcy–Forchheimer for porous media with heat transfer.

The primary attention is to study the cilia-oriented, induced MHD flow of a Sutterby fluid via an asymmetric tapered channel, which is yet to be available in the literature. This current effort is devoted to filling this gap. In addition, the Darcy–Brinkman–Forchheimer law has been utilized for a porous medium with a viscous dissipation function. Temperature maintenance has a significant importance in physiological flows, especially in cilia-oriented flows. Moreover, the investigation of heat transport is also part of this study. The homotopy perturbation solutions of non-dimensional equations have been found with the support of “Mathematica” software. The pictorial results of different emerging parameters are also elaborated along with their physical interpretations.

2. Modeling of Sutterby Fluid in a Cilia-Oriented Asymmetric Tapered Channel

Assume that an incompressible Sutterby fluid is moving through a non-Darcy porous medium in the presence of a magnetic field. The geometry of the problem has a specific significance in physiological flow. Therefore, we consider an asymmetric tapered channel whose walls are cilia-oriented and moving towards the right as a metachronal wave with a constant speed c. Moreover, Cartesian coordinates ($\tilde{X}$, $\tilde{Y}$) are assumed, the direction of flow in a tapered channel is parallel to the $\tilde{X}$ axis, and the $\tilde{Y}$ axis is perpendicular to the flow direction. The cilia-oriented asymmetric tapered channel’s lower and upper walls with temperatures ${\tilde{T}}_0$ and ${\tilde{T}}_1$ are sustained, respectively, as presented in Figure 1.

The following is a well-defined Cauchy stress tensor for a Sutterby fluid:

$$\tilde{\mathit{\tau }}=-\tilde{p}\tilde{\mathit{I}}+\tilde{\mathit{S}},$$

(1)

in which $\tilde{\mathit{S}}$ is the extra stress tensor. For a Sutterby fluid its mathematical form is presented as [37]

$$\tilde{\mathit{S}}=\tilde{\mu }{\left[\frac{{\sinh }^{-1}\left(D\dot{\mathit{\Omega }}\right)}{D\dot{\mathit{\Omega }}}\right]}^m\tilde{\mathbf{\Delta }}\cong \tilde{\mu }\left[1-\frac{m}{6}{\left(D\dot{\mathit{\Omega }}\right)}^2\right]\tilde{\mathbf{\Delta }},$$

(2)

$$\tilde{\mathbf{\Delta }}=\left(\nabla \tilde{\mathit{V}}\right)+{\left(\nabla \tilde{\mathit{V}}\right)}^t,$$

(3)

$$\dot{\mathit{\Omega }}=\sqrt{\frac{1}{2}tr\left({\tilde{\mathbf{\Delta }}}^2\right)},$$

(4)

here $\tilde{p}$, $\tilde{\mathit{I}}$, $\tilde{\mu }$, $D$, and $m$ are the symbols that signify the pressure, identity matrix, viscosity, and material constants, respectively. Equations (3) and (4) define the deformation tensor and the second invariant strain tensor, respectively. In the case of two-dimensional flow, Equation (2) in component form becomes:

$${\tilde{S}}_{\tilde{X}\tilde{X}}=2\tilde{\mu }\stackrel{⌣}{\chi }\left(\frac{\partial \tilde{U}}{\partial \tilde{X}}\right),$$

(5)

$${\tilde{S}}_{\tilde{X}\tilde{Y}}=\tilde{\mu }\stackrel{⌣}{\chi }\left(\frac{\partial \tilde{V}}{\partial \tilde{X}}+\frac{\partial \tilde{U}}{\partial \tilde{Y}}\right)={\tilde{S}}_{\tilde{Y}\tilde{X}},$$

(6)

$${\tilde{S}}_{\tilde{Y}\tilde{Y}}=2\tilde{\mu }\stackrel{⌣}{\chi }\left(\frac{\partial \tilde{V}}{\partial \tilde{Y}}\right),$$

(7)

$$\mathrm{where} \stackrel{⌣}{\chi }=1-\frac{m{D}^2}{6}\left\{{\left(\frac{\partial \tilde{V}}{\partial \tilde{X}}+\frac{\partial \tilde{U}}{\partial \tilde{Y}}\right)}^2+2\left\{{\left(\frac{\partial \tilde{V}}{\partial \tilde{Y}}\right)}^2+{\left(\frac{\partial \tilde{U}}{\partial \tilde{X}}\right)}^2\right\}\right\}.$$

(8)

The rheological equations that govern the Sutterby flow in the existence of a body force are stated below.

$$\nabla \cdot \tilde{\mathit{V}}=0,$$

(9)

$${\tilde{\phi }}^{-1}\rho \left(\frac{\partial \tilde{\mathit{V}}}{\partial \tilde{t}}+{\tilde{\phi }}^{-1}\left(\tilde{\mathit{V}}\cdot \nabla \right)\tilde{\mathit{V}}\right)+\nabla \tilde{p}=\nabla \cdot \tilde{\mathit{S}}-\tilde{\mathit{V}}\left(\left(1-\frac{m}{6}{\left(D\dot{\mathit{\Omega }}\right)}^2\right)\frac{\tilde{\mu }}{\tilde{K}}+\frac{\rho \tilde{C}\left|\tilde{\mathit{V}}\right|}{{\tilde{K}}^{\frac{1}{2}}}\right)+\tilde{\mathit{J}}\times \tilde{\mathit{B}},$$

(10)

$$\rho \tilde{\xi }\frac{d\tilde{\mathit{T}}}{d\tilde{t}}=\tilde{\tau }\cdot \mathit{L}-\nabla \cdot \tilde{\mathit{\vartheta }}+\left(\frac{\tilde{\mu }}{\tilde{K}}\left(1-\frac{m}{6}{\left(D\dot{\mathit{\Omega }}\right)}^2\right)+\frac{\rho \tilde{C}\left|\tilde{\mathit{V}}\right|}{{\tilde{K}}^{\frac{1}{2}}}\right)\tilde{\mathit{V}}\cdot \tilde{\mathit{V}},$$

(11)

$$\begin{array}{l}0=\nabla \cdot \tilde{\mathit{B}}, \tilde{\mathit{J}}=\tilde{\sigma }\left(\tilde{\mathit{E}}+\tilde{\mathit{V}}\times \tilde{\mathit{B}}\right)\\ \stackrel{⌢}{\mu }\tilde{\mathit{J}}=\nabla \times \tilde{\mathit{B}},-\frac{\partial \tilde{\mathit{B}}}{\partial \tilde{t}}=\nabla \times \tilde{\mathit{E}}\end{array}\},$$

(12)

where $\tilde{\mathit{V}}=[\tilde{U}(\tilde{X},\tilde{Y},\tilde{t}),\tilde{V}(\tilde{X},\tilde{Y},\tilde{t}),0]$ is the velocity field, $\rho$ is the density of the fluid, $\tilde{\phi }$ is the porosity parameter, $\tilde{t}$ is the time, $\tilde{C}$ is the Forchheimer coefficient, $\tilde{K}$ is the permeability, $\tilde{\xi }$ is the specific heat, $\tilde{\mathit{\vartheta }}$ is the heat flux vector, which is defined as $\tilde{\mathit{\vartheta }}\left(=-\widehat{k}d\mathit{T}/d\tilde{t}\right)$, $\tilde{\sigma }$ is the electric conductivity, $\stackrel{⌢}{\mu }$ is the magnetic permeability, $\tilde{\mathit{J}}={\tilde{\mathit{J}}}_1+{\tilde{\mathit{J}}}_0$ is the total current density, which is the sum of the induced and applied current density, and similarly, $\tilde{\mathit{B}}={\tilde{\mathit{B}}}_0+{\tilde{\mathit{B}}}_1$ and $\tilde{\mathit{E}}={\tilde{\mathit{E}}}_0+{\tilde{\mathit{E}}}_1$ are the total magnetic and electric fields, respectively. Now, by considering the high Reynolds number approximation for the incorporation of an induced magnetic field, which is due to the induced currents ${\tilde{\mathit{J}}}_1$, then the advection of the magnetic field is also dominant and thus ${\tilde{\mathit{E}}}_1={\tilde{\mathit{B}}}_1={\tilde{\mathit{J}}}_1\ne 0$ and the induction equations can be written as follows:

$$\nabla \times \left(\tilde{\mathit{V}}\times \tilde{\mathit{B}}\right)+\tilde{\zeta } {\nabla }^2\tilde{\mathit{B}}=\frac{\partial \tilde{\mathit{B}}}{\partial \tilde{t}},$$

(13)

here $\tilde{\zeta }=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\stackrel{⌢}{\mu }\tilde{\sigma }$}\right.$ is the magnetic diffusivity. By utilizing the value of ${\tilde{\mathit{J}}}_1$ with the use of Ampere’s law rather than Ohm’s law, we obtain

$$\tilde{\mathit{J}}=\frac{1}{\stackrel{⌢}{\mu }}\left(\nabla \times {\tilde{\mathit{B}}}_1\right)+\tilde{\sigma }\left({\tilde{\mathit{E}}}_0+\tilde{\mathit{V}}\times {\tilde{\mathit{B}}}_0\right),$$

(14)

where the transverse constant applied magnetic field for two-dimensional flow yields

$$\begin{array}{l}{\tilde{\mathit{B}}}_0=\left[0,{\tilde{B}}_0,0\right], {\tilde{\mathit{E}}}_0=\left[0,0,\tilde{E}\right],\\ {\tilde{\mathit{B}}}_1=\left[{\tilde{B}}_1\left(\tilde{X},\tilde{Y},\tilde{t}\right),{\tilde{B}}_2\left(\tilde{X},\tilde{Y},\tilde{t}\right),0\right],\\ {\tilde{\mathit{E}}}_1=\left[0,0,{\tilde{E}}_3\left(\tilde{X},\tilde{Y},\tilde{t}\right)\right],\end{array}\}.$$

(15)

Thus, the expression for the Lorentz force takes the form [38,39]:

$$\tilde{\mathit{J}}\times \tilde{\mathit{B}}=\left[\begin{array}{c}-\left({\tilde{B}}_0+{\tilde{B}}_2\right)\times \left\{\tilde{\sigma }\left(\tilde{E}+\tilde{U}{\tilde{B}}_0\right)+\frac{1}{\stackrel{⌢}{\mu }}\left(\frac{\partial {\tilde{B}}_2}{\partial \tilde{X}}-\frac{\partial {\tilde{B}}_1}{\partial \tilde{Y}}\right)\right\},\\ +{\tilde{B}}_1\times \left\{\tilde{\sigma }\left(\tilde{E}+\tilde{U}{\tilde{B}}_0\right)+\frac{1}{\stackrel{⌢}{\mu }}\left(\frac{\partial {\tilde{B}}_2}{\partial \tilde{X}}-\frac{\partial {\tilde{B}}_1}{\partial \tilde{Y}}\right)\right\}\end{array}\right].$$

(16)

By ignoring the ${\tilde{\mathit{B}}}_1$ (i.e., the induced magnetic field), Equation (16) can be recast as

$$\tilde{\mathit{J}}\times \tilde{\mathit{B}}=[\tilde{\sigma }\left(\tilde{E}+\tilde{U}{\tilde{B}}_0\right),0].$$

(17)

The cilia-assisted flow has been analyzed by utilizing the sublayer model or the envelope model. Here, we consider the envelope model for the tapered channel; thus, the mathematical expression for the cilia tip envelope can be defined as [40]

$$\begin{array}{c}\tilde{Y}={\tilde{H}}_l={\tilde{F}}_1\left(\tilde{X},\tilde{t}\right)=-\left(a+\tilde{X}\tan \gamma \right)+a\widetilde{ϵ}\cos 2\left[-\frac{\phi }{4}+\frac{\pi }{\lambda }\left(\tilde{X}-c\tilde{t}\right)\right]\\ \tilde{Y}={\tilde{H}}_{up}={\tilde{F}}_2\left(\tilde{X},\tilde{t}\right)=a+\tilde{X}\tan \gamma +a\widetilde{ϵ}\cos 2\left[\frac{\phi }{4}+\frac{\pi }{\lambda }\left(\tilde{X}-c\tilde{t}\right)\right] \end{array}\},$$

(18)

where ${\tilde{H}}_l$ and ${\tilde{H}}_{up}$ represent the lower and upper walls, respectively, of the cilia-oriented asymmetric tapered channel, $\gamma$ is an inclined angle for unperturbed walls, $c$ is the speed, $\lambda$ is the wavelength of the metachronal wave, $\widetilde{ϵ}$ is the non-dimensional measure concerning the cilia length, and $\phi$ is the phase difference. As reported by Sleigh, the movement of the cilia tips in an elliptical path is such that the mathematical expression for the horizontal position of the cilia tips can be represented by

$$\tilde{X}=\tilde{\omega }\widetilde{ϵ}a\sin \left[\left(\tilde{X}-c\tilde{t}\right)\frac{2\pi }{\lambda }\right]+{\tilde{X}}_0={\tilde{F}}_0\left({\tilde{X}}_0,\tilde{X},\tilde{t}\right),$$

(19)

where $\tilde{\omega }$ and ${\tilde{X}}_0$ represent the measure of the eccentricity of the elliptic path and the reference position of the cilia, respectively. The horizontal and vertical velocities of the cilia positions at ${\tilde{X}}_0$ are introduced as follows

$$\tilde{U}=\frac{\partial {\tilde{F}}_0}{\partial \tilde{t}}{|}_{{\tilde{X}}_0}=\frac{\partial {\tilde{F}}_0}{\partial \tilde{X}}\frac{\partial \tilde{X}}{\partial \tilde{t}}+\frac{\partial {\tilde{F}}_0}{\partial \tilde{t}}=\frac{\partial {\tilde{F}}_0}{\partial \tilde{X}}\tilde{U}+\frac{\partial {\tilde{F}}_0}{\partial \tilde{t}},$$

(20)

$$\begin{array}{c}\tilde{V}=\frac{\partial {\tilde{F}}_1}{\partial \tilde{t}}{|}_{{\tilde{X}}_0}=\frac{\partial {\tilde{F}}_1}{\partial \tilde{X}}\frac{\partial \tilde{X}}{\partial \tilde{t}}+\frac{\partial {\tilde{F}}_1}{\partial \tilde{t}}=\frac{\partial {\tilde{F}}_1}{\partial \tilde{X}}\tilde{U}+\frac{\partial {\tilde{F}}_1}{\partial \tilde{t}},\\ \tilde{V}=\frac{\partial {\tilde{F}}_2}{\partial \tilde{t}}{|}_{{\tilde{X}}_0}=\frac{\partial {\tilde{F}}_2}{\partial \tilde{X}}\frac{\partial \tilde{X}}{\partial \tilde{t}}+\frac{\partial {\tilde{F}}_2}{\partial \tilde{t}}=\frac{\partial {\tilde{F}}_2}{\partial \tilde{X}}\tilde{U}+\frac{\partial {\tilde{F}}_2}{\partial \tilde{t}} \end{array}\}.$$

(21)

Substituting Equation (19) into (20) and Equation (18) into (21) provides the alternative form of the velocity components, as follows:

$$\tilde{U}=\frac{2ac\tilde{\omega }\widetilde{ϵ}\pi \cos \left[\frac{2\pi }{\lambda }\left(\tilde{X}-c\tilde{t}\right)\right]}{2a\tilde{\omega }\widetilde{ϵ}\pi \cos \left[\frac{2\pi }{\lambda }\left(\tilde{X}-c\tilde{t}\right)\right]-\lambda },$$

(22)

$$\tilde{V}=\pm \frac{2ac\widetilde{ϵ}\pi \left(\sin \left[\frac{\phi }{2}\pm \frac{2\pi }{\lambda }\left(\tilde{X}-c\tilde{t}\right)\right]+\tilde{\omega }\cos \left[\frac{2\pi }{\lambda }\left(\tilde{X}-c\tilde{t}\right)\right]\tan \gamma \right)}{2a\tilde{\omega }\widetilde{ϵ}\pi \cos \left[\frac{2\pi }{\lambda }\left(\tilde{X}-c\tilde{t}\right)\right]-\lambda }.$$

(23)

Now, we are introducing some dimensionless terms as follows:

$$\begin{array}{l}x=\frac{\tilde{X}}{\lambda },y=\frac{\tilde{Y}}{a},t=\frac{c\tilde{t}}{\lambda },p_m=\frac{a^2{\tilde{p}}_m}{c\lambda \mu },u=\frac{\tilde{U}}{c},S=\frac{a}{\tilde{\mu }c}\tilde{S},v=\frac{\tilde{V}}{\delta c},\delta =\frac{a}{\lambda },h_x=\frac{{\tilde{B}}_1}{{\tilde{B}}_0},\\ h_y=\frac{{\tilde{B}}_2}{{\tilde{B}}_0},E=\frac{\tilde{E}}{{\tilde{B}}_{0}c},\theta =\frac{\tilde{T}-{\tilde{T}}_s}{{\tilde{T}}_0-{\tilde{T}}_s},h_l=\frac{{\tilde{H}}_l}{a},h_{up}=\frac{{\tilde{H}}_{up}}{a},\end{array}\}.$$

(24)

Here, ${\tilde{p}}_m$ is the modified pressure, which is defined as ${\tilde{p}}_m=\tilde{p}+\left\{{\left({\tilde{B}}_1\right)}^2+{\left({\tilde{B}}_2\right)}^2\right\}/2\tilde{\mu }$. Now, with the help of Equation (24), the leading Equations (9)–(12) in dimensionless form can be inscribed as:

$$\begin{array}{l}{\tilde{\phi }}^{-1}Re\delta \left(\frac{\partial u}{\partial t}+{\tilde{\phi }}^{-1}u\frac{\partial u}{\partial x}+{\tilde{\phi }}^{-1}v\frac{\partial u}{\partial y}\right)\\ =-\frac{\partial p_m}{\partial x}+\delta \frac{\partial }{\partial x}\left(S_{xx}\right)+\frac{\partial }{\partial y}\left(S_{xy}\right)-\left(F_r\sqrt{u^2+{\delta }^2{v}^2}\right)u\\ \hspace{1em}-\frac{u}{Da}\left\{1-\eta \left({\left(\frac{\partial u}{\partial y}+{\delta }^2\frac{\partial v}{\partial x}\right)}^2+2\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]{\delta }^2\right)\right\}\\ \hspace{1em}-M^2\left(1+h_y\right)\left\{E+u+\frac{1}{R_m}\left(\delta \frac{\partial h_y}{\partial x}-\frac{\partial h_x}{\partial y}\right)\right\},\end{array}$$

(25)

$$\begin{array}{l}{\tilde{\phi }}^{-1}Re{\delta }^2\left(\frac{\partial v}{\partial t}+{\tilde{\phi }}^{-1}u\frac{\partial v}{\partial x}+{\tilde{\phi }}^{-1}v\frac{\partial v}{\partial y}\right)\\ =-\frac{\partial p_m}{\partial y}+{\delta }^2\frac{\partial }{\partial x}\left(S_{yx}\right)+\delta \frac{\partial }{\partial y}\left(S_{yy}\right)-{\delta }^2\left(F_r\sqrt{u^2+{\delta }^2{v}^2}\right)v\\ \hspace{1em}-\frac{{\delta }^{2}v}{D_a}\left\{1-\eta \left({\left(\frac{\partial u}{\partial y}+{\delta }^2\frac{\partial v}{\partial x}\right)}^2+2\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]{\delta }^2\right)\right\}\\ \hspace{1em}-M^2\delta h_x\left\{E+u+\frac{1}{R_m}\left(\delta \frac{\partial h_y}{\partial x}-\frac{\partial h_x}{\partial y}\right)\right\},\end{array}$$

(26)

$$\begin{array}{l}{\tilde{\phi }}^{-1}RePr\delta \left(\frac{\partial \theta }{\partial t}+{\tilde{\phi }}^{-1}u\frac{\partial \theta }{\partial x}+{\tilde{\phi }}^{-1}v\frac{\partial \theta }{\partial y}\right)\\ =Br\left[{\left(\frac{\partial u}{\partial x}\right)}^2{\delta }^2+\frac{1}{2}{\left(\frac{\partial v}{\partial x}{\delta }^2+\frac{\partial u}{\partial y}\right)}^2+{\delta }^2{\left(\frac{\partial v}{\partial y}\right)}^2\right]\\ \hspace{1em}\times \left[\left\{1-\eta \left({\left(\frac{\partial v}{\partial x}{\delta }^2+\frac{\partial u}{\partial y}\right)}^2+2\left[{\left(\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial u}{\partial x}\right)}^2\right]{\delta }^2\right)\right\}\right]\\ \hspace{1em}+\left(\frac{{\delta }^2{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}\right)+Br{F}_r\sqrt{u^2+{\delta }^2{v}^2}\left(u^2+{\delta }^2{v}^2\right)\\ \hspace{1em}+\frac{Br}{Da}\left\{1-\eta \left(2{\delta }^2\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]+{\left(\frac{\partial u}{\partial y}+{\delta }^2\frac{\partial v}{\partial x}\right)}^2\right)\right\}\left(u^2+{\delta }^2{v}^2\right),\\ \end{array}$$

(27)

$$\frac{\partial }{\partial y}\left(u\left(h_y+1\right)-\delta v{h}_x\right)+\frac{1}{R_m}\left(\frac{{\partial }^2{h}_x}{\partial y^2}+{\delta }^2\frac{{\partial }^2{h}_x}{\partial x^2}\right)=\delta \frac{\partial h_x}{\partial t},$$

(28)

$$\frac{\partial }{\partial x}\left(u\left(h_y+1\right)-\delta v{h}_x\right)-\frac{1}{R_m}\left(\frac{{\partial }^2{h}_x}{\partial y^2}+{\delta }^2\frac{{\partial }^2{h}_x}{\partial x^2}\right)=-\frac{\partial h_y}{\partial t},$$

(29)

where $h_x$ and $h_y$ denote the magnetic force functions, which can be written as: $h_x=\frac{\partial \varphi }{\partial y},$ $h_y=-\delta \frac{\partial \varphi }{\partial x}$, respectively, and some new dimensionless parameters that are introduced in Equations (25)–(27) are:

$$\begin{array}{l}Re=\frac{\rho ac}{\tilde{\mu }},F_r=\frac{\rho Fc{a}^2}{\sqrt{\tilde{K}}\tilde{\mu }},Da=\frac{\tilde{K}}{a^2},M^2=\frac{\tilde{\sigma }{a}^2{\tilde{B}}^2{}_0}{\tilde{\mu }}\\ \eta =\frac{m{D}^2{c}^2}{6a^2},Pr=\frac{\tilde{\mu }\tilde{\xi }}{\widehat{k}},Ec=\frac{c^2}{\tilde{\xi }({\tilde{T}}_0-{\tilde{T}}_s)},Br=PrEc\end{array}\}.$$

(30)

Here, $Re$ is the Reynolds number, $F_r$ is the Forchheimer number, $Da$ is the Darcy number, $M$ is the Hartman number, $\eta$ is the Sutterby fluid parameter, $Pr$ is the Prandtl number, $Ec$ is the Eckert number, and $Br$ is the Brinkmann number.

Finally, we employed the lubrication approach in Equations (25)–(29), which leads to:

$$\frac{\partial p_m}{\partial x}=\frac{\partial }{\partial y}\left(S_{xy}\right)-F_r{u}^2-\frac{u}{Da}\left\{1-\eta {\left(\frac{\partial u}{\partial y}\right)}^2\right\}-M^2\left\{E+u-\frac{1}{R_m}\left(\frac{\partial h_x}{\partial y}\right)\right\},$$

(31)

$$\frac{\partial p_m}{\partial y}=0,$$

(32)

$$u=-\frac{1}{R_m}\left(\frac{\partial h_x}{\partial y}\right),$$

(33)

$$\begin{array}{l}\frac{{\partial }^2\theta }{\partial y^2}+Br{\left(\frac{\partial u}{\partial y}\right)}^2\left\{1-\eta {\left(\frac{\partial u}{\partial y}\right)}^2\right\}+\frac{Br}{Da}\left\{1-\eta {\left(\frac{\partial u}{\partial y}\right)}^2\right\}{u}^2+Br{F}_r{u}^3=0,\\ \end{array}$$

(34)

where

$$\begin{array}{l}{S}_{xx}=0,\\ S_{xy}=S_{yx}=\frac{\partial u}{\partial y}\left\{1-\eta {\left(\frac{\partial u}{\partial y}\right)}^2\right\}\\ S_{yy}=0,\end{array}\}.$$

(35)

Further, the boundary condition in dimensionless form can be written as:

$$\begin{array}{c}u{|}_{y=h_l}=u{|}_{y=h_{up}}=\frac{2a\tilde{\omega }\widetilde{ϵ}\pi \delta \cos [2\pi \left(x-t\right)]}{2a\tilde{\omega }\widetilde{ϵ}\pi \delta \cos [2\pi \left(x-t\right)]-1}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}v{|}_{ y=h_l}=\frac{\partial h_l}{\partial t},v{|}_{ y=h_{up}}=\frac{\partial h_{up}}{\partial t}\end{array}\},$$

(36)

$$\theta {|}_{y=h_l}=1,\theta {|}_{y=h_{up}}=\Gamma \},$$

(37)

where $\Gamma$ denotes the temperature ratio.

3. Analytical Results

Now, we will determine the series solutions of Equation (31) by utilizing the well-known homotopy perturbation method [41]. For Equation (31), the homotopy can be constructed as:

$$\stackrel{⌣}{H}\left(v,\widehat{\beta }\right)=\left(1-\widehat{\beta }\right)\left({\widehat{L}}_r\left(v\right)-{\widehat{L}}_r\left(v_0\right)\right)+\widehat{\beta }\left[\begin{array}{l}{\widehat{L}}_r\left(v\right)-3\eta \frac{{\partial }^{2}u}{\partial y^2}{\left(\frac{\partial u}{\partial y}\right)}^2-F_r\left(u^2\right)\\ -\frac{u}{Da}\left\{1-\eta {\left(\frac{\partial u}{\partial y}\right)}^2\right\}-M^2\left(E+2u\right)-\frac{\partial p_m}{\partial x}\end{array}\right],$$

(38)

here $\widehat{\beta }\in \left[0,1\right]$ is the embedding parameter and ${\widehat{L}}_r$ is the linear operator, which is defined as:

$${\widehat{L}}_r=\frac{{\partial }^2}{\partial y^2}.$$

(39)

Furthermore, their initial guess and the perturbation expansion are defined as follows:

$$v_0=\left(y-h_l\right)\left(y-h_{up}\right)+N, \mathrm{where} N=\frac{2a\tilde{\omega }\widetilde{ϵ}\pi \delta \cos [2\pi \left(x-t\right)]}{2a\tilde{\omega }\widetilde{ϵ}\pi \delta \cos [2\pi \left(x-t\right)]-1},$$

(40)

$$v={\stackrel{⌢}{v}}_0+\widehat{\beta }{\stackrel{⌢}{v}}_1+{\widehat{\beta }}^2{\stackrel{⌢}{v}}_2+\dots$$

(41)

Substituting Equation (41) in Equation (31) gives us the system of differential equations with their boundary conditions as follows:

a.

The zero-order system and its solution:

$${\widehat{L}}_r\left({\stackrel{⌢}{v}}_0\right)-{\widehat{L}}_r\left(v_0\right)=0, \mathrm{also} {\stackrel{⌢}{v}}_0=N \mathrm{at} y=h_l,y=h_{up}.$$

(42)

$${\stackrel{⌢}{v}}_0=\left(y-h_1\right)\left(y-h_2\right)+N$$

(43)

b.

The first-order system and its solution:

$$\begin{gathered} L\left({\stackrel{⌢}{v}}_1\right)+2-3\gamma \frac{{\partial }^2{\stackrel{⌢}{v}}_0}{\partial y^2}{\left(\frac{\partial {\stackrel{⌢}{v}}_0}{\partial y}\right)}^2-\frac{{\stackrel{⌢}{v}}_0}{Da}\left\{1-\eta {\left(\frac{\partial {\stackrel{⌢}{v}}_0}{\partial y}\right)}^2\right\}-F_r{\stackrel{⌢}{v}}^2{}_0-2M^2{\stackrel{⌢}{v}}_0-M^{2}E-\frac{\partial p_m}{\partial x}=0, \\ \mathrm{also}, {\stackrel{⌢}{v}}_1=0 \mathrm{at} y=h_l,y=h_{up}. \end{gathered}$$

(44)

$$\begin{array}{l}{\stackrel{⌢}{v}}_1=\frac{1}{60Da}\left(h_l-y\right)\left(h_{up}-y\right)\\ \times \left[\begin{array}{l}-5\left(h_l^2-3h_l{h}_{up}+h_l^2-6N+\left(h_l+h_{up}\right)y-y^2\right)\\ +\\ \left(\begin{array}{l}{h}_l^4-5h_l^3{h}_{up}-5h_l{h}_{up}^3+h_{up}^4-10h_l^{2}N\\ -10h_{up}^{2}N+\left(h_l+h_{up}\right)\left(h_l^2+14h_l{h}_{up}\right)\\ +h_{up}^2+20Ny-\left(3\left(3h_l+h_{up}\right)\left(h_l+3h_{up}\right)+20N\right)\\ y^2+16\left(h_l+h_{up}\right)y^3-8y^4\end{array}\right)\eta +\\ Da\left(F_r\left(\begin{array}{l}{h}_l^4-5h_l^3{h}_{up}+10h_l^2{h}_{up}^2\\ -5h_l{h}_{up}^3+h_{up}^4-10h_l^{2}N+30h_l{h}_{up}-10h_{up}^{2}N+30N^2+\\ \left(h_l+h_{up}\right)\left(h_l^2-6h_l{h}_{up}+h_{up}^2-10N\right)y+\left(h_l^2+10h_l{h}_{up}+h_{up}^2+10N\right)y^2-\\ 4\left(h_l+h_{up}\right)y^3+2y^4+10\left(\begin{array}{l}3\left(-2+\frac{d{p}_m}{dx}\right)+M^2\left(\begin{array}{l}3E-h_l^2+3h_l{h}_{up}-h_{up}^2\\ +6N-\left(h_l+h_{up}\right)y^2\end{array}\right)+\\ 6\left(h_l^2+h_{up}^2-2\left(h_l+h_{up}\right)y+2y^2\right)\end{array}\right)\end{array}\right)\right)\end{array}\right]\}.\end{array}$$

(45)

c.

The second-order system and its solution:

$$\begin{gathered} L\left({\stackrel{⌢}{v}}_2\right)-3\eta \frac{\partial }{\partial y}\left[{\left(\frac{\partial {\stackrel{⌢}{v}}_0}{\partial y}\right)}^2\frac{\partial {\stackrel{⌢}{v}}_1}{\partial y}\right]-\left(\frac{{\tilde{v}}_1}{Da}\left\{1-\eta {\left(\frac{\partial {\stackrel{⌢}{v}}_0}{\partial y}\right)}^2\right\}-\frac{1}{Da}2\eta \frac{\partial {\stackrel{⌢}{v}}_0}{\partial y}\frac{\partial {\stackrel{⌢}{v}}_1}{\partial y}\right)-2F_r{\stackrel{⌢}{v}}_0{\stackrel{⌢}{v}}_1-2M^2{\stackrel{⌢}{v}}_1=0, \\ \mathrm{with}, {\stackrel{⌢}{v}}_2=0 \mathrm{at} y=h_l,y=h_{up}, \end{gathered}$$

(46)

$${\stackrel{⌢}{v}}_2={\stackrel{⌢}{v}}_{2,0}+{\stackrel{⌢}{v}}_{2,1}y+{\stackrel{⌢}{v}}_{2,2}{y}^2+{\stackrel{⌢}{v}}_{2,3}{y}^3+{\stackrel{⌢}{v}}_{2,4}{y}^4+{\stackrel{⌢}{v}}_{2,5}{y}^5+{\stackrel{⌢}{v}}_{2,6}{y}^6+{\stackrel{⌢}{v}}_{2,7}{y}^7+{\stackrel{⌢}{v}}_{2,8}{y}^8+{\stackrel{⌢}{v}}_{2,9}{y}^9+{\stackrel{⌢}{v}}_{2,10}{y}^{10},$$

(47)

and so on.

Finally, we achieved the desired solution of the differential Equations (42)–(47) up to the second order by taking $\widehat{\beta }\to 1,$ which leads to

$$u=v={\stackrel{⌢}{v}}_0+{\stackrel{⌢}{v}}_1+{\stackrel{⌢}{v}}_2+\dots$$

(48)

Now, Equation (48) is utilized in Equation (34) and we have

$$\theta ={\psi }_0+{\psi }_{1}y+{\psi }_2{y}^2+\dots +{\psi }_{40}{y}^{40},$$

(49)

where (${\stackrel{⌢}{v}}_{2,0}-{\stackrel{⌢}{v}}_{2,10}$) in Equation (47) and (${\psi }_0-{\psi }_{40}$) in Equation (49) are the list of constants, which can also be found by utilizing the “Mathematica” software.

In fixed frame, the mathematical expression for the instantaneous flow rate is given below

$$Q={\displaystyle {\int }_{h_l}^{h_{up}}u(x,y,t)} dy.$$

(50)

In addition, the non-dimensional form of the current density in terms of the magnetic force function can be inscribed as:

$$J_z+\frac{{\partial }^2\varphi }{\partial y^2}=R_m\left(E+u\right),$$

(51)

with the following boundary conditions

$$\varphi {|}_{y=h_l}=0, \varphi {|}_{y=h_{up}}=0.$$

(52)

The expression of the pressure-rise term can be evaluated numerically as:

$$\Delta p_m=\underset{0}{\overset{1}{{\displaystyle \int }}}\frac{d{p}_m}{dx}dx.$$

(53)

4. Pictorial Representation and Discussion

In this section, we discuss the effects of different physical parameters that are based on the graphical results of the axial velocity ($u$), temperature ($\theta$), magnetic force function ($\varphi$), induced magnetic field ($h_x$),magnitude of current density ($J_z$), and pressure rise ($\Delta p$). Moreover, a further discussion on trapping phenomena is also part of this section.

4.1. Velocity Profiles

The influences of $M$, $\eta$, $\gamma$, $\tilde{\omega }$, $Da$, $F_r$, and $\widetilde{ϵ}$ on the axial velocity are depicted in Figure 2a–g, respectively. Figure 2a illustrates that as we increase the influence of the Hartman number ($M$) by taking some non-negative values, the velocity of the cilia-driven flow slowly declines in the central area of the cilia-oriented asymmetric tapered channel and rises closer to the walls. It is noticeable because ($M$) is the ratio of electromagnetic forces to viscous forces, which creates a drag-like resistive force to slow down the velocity of the fluid. Consequently, the velocity profile shows a decreasing behavior. Moreover, in Figure 2b–d, the same situation occurred for diverse values of the Sutterby fluid parameter ($\eta$), inclined angle ($\gamma$), and eccentricity parameter ($\tilde{\omega }$). From Figure 2e, it is noticed that the velocity increases in the core region and slows down near the walls of a ciliated tapered channel for several values of the Darcy number ($Da$). Accordingly, we can say that large values of ($Da$) create less resistance, which means more permeability, which is in favor of the fluid flow. Likewise, in Figure 2f the velocity of the fluid is enhanced when the Forchheimer number ($F_r$) is elevated. In Figure 2g, when we elevate the cilia length ($\widetilde{ϵ}$)the velocity profile shows some hindrance in the middle of the cilia-oriented asymmetric tapered channel but close to the walls, its behavior is quite different.

4.2. Magnetic Force Function

Variations in the magnitude of the magnetic force function for all emerging parameters are demonstrated in Figure 3a–h. It can be noticed that in the central area of the cilia-oriented asymmetric tapered channel, the magnetic force function progressively decreases due to the non-negative values of the Hartman number ($M$), Sutterby fluid parameter ($\eta$), and eccentricity parameter ($\tilde{\omega }$), as shown in Figure 3a–c. On the other hand, quite different behavior is observed for the magnetic force function when the magnetic Reynolds number ($R_m$) is elevated, and it is defined as the ratio of magnetic advection to magnetic diffusion. Therefore, higher ($R_m$) increase the magnetic permeability, which is why the graph of magnetic force function goes up, as exposed in Figure 3d. Similarly, in Figure 3e–f we can easily observe that in the core region of the ciliated channel the magnetic force function increases with increasing Darcy number ($Da$) and Forchheimer number ($F_r$). Figure 3g–h illustrate diverse values of the inclined angle ($\gamma$) and cilia length ($\widetilde{ϵ}$). We found that large values of inclined angle ($\gamma$) also increase the amount of magnetic force function, as shown in Figure 3g, but the quite opposite behavior is noticed in the case of ($\widetilde{ϵ}$).

4.3. Magnetic Field Characteristics

The magnetic field characteristics for altered values of the Hartman number ($M$) and magnetic Reynolds number ($R_m$) are exposed in Figure 4a,b. It is observed that the induced magnetic field decreases in the region ($y<0$) but increases in the region ($y>0$) when ($M$) increases. However, in Figure 4b, the behavior of the induced magnetic field increases in the region ($y\prec 0$) and decreases in the region ($y\succ 0$) when ($R_m$) is elevated. Moreover, these types of dual behavior of the induced magnetic field are noticed for the Sutterby fluid ($\eta$), Darcy number ($Da$), Forchheimer number ($F_r$), eccentricity ($\tilde{\omega }$), and inclined angle ($\gamma$), as shown in Figure 4c–g. From Figure 4h, it can be seen that the power of an induced magnetic field decreases over the whole domain when the cilia length ($\widetilde{ϵ}$) increases.

4.4. Current Density

The variation in the current density profile against several emerging parameters are depicted in Figure 5a–i. These graphs have a parabolic shape and show almost the same trend as we already discussed for the velocity profile. Thus, we can see that the enormous value of the Hartman number ($M$), the Sutterby parameter ($\eta$), and the inclined angle ($\gamma$) are opposing the current density near the central area of the channel but are supportive near the channel boundaries, as displayed in Figure 5a–c. We can notice that the rising values of the eccentricity ($\tilde{\omega }$) parameter have no significant impacts on the profile of current density, as depicted in Figure 5d. From Figure 5e,f, it can be realized that enormous values of the Darcy number ($Da$) and Forchheimer number ($F_r$) both support the magnitude of the current density at the core part of the channel while behaving differently at the boundaries of the channel. Due to the elevated values of the electric field ($E$) and the magnetic Reynolds number ($R_m$) the profiles of the current density are raised dramatically, as demonstrated in Figure 5g,h. Figure 5i illustrates the influence of the cilia length ($\widetilde{ϵ}$) on the current density. It can be noticed that as we increase the value of ($\widetilde{ϵ}$)the magnitude of the current density shows hindrance near the maximal portion of the channel but the lower and upper portion of the channel behave completely differently from each other.

4.5. Pressure Rise

In this subsection, we analyze the pumping characteristics of the entire flow region. Figure 6 shows that the relationship between the pressure differential and volumetric flow rate is inversely proportional, which means that the pressure increases yield bigger values for lower volumetric rates and vice versa. The pressure differential is likewise greater for the absent time-averaged flow rate. All of the deviations, however, are linearly dependent. Figure 6a–g is illustrated to check the influences of diverse physical parameters such as the Hartman number ($M$), Sutterby fluid parameter ($\eta$), inclined angle ($\gamma$), Darcy number ($Da$), cilia length ($\widetilde{ϵ}$), Forchheimer number ($F_r$), and eccentricity ($\tilde{\omega }$) on the pressure-rise profile. From Figure 6a, the impact of the Hartman number ($M$) for different values shows that the pumping rate increases in the region ($-0.4<Q\le -1$) but declines in the region ($-0.4<Q\le -1$). On the other hand, the pumping rate rises in the region ($-1\le Q<0$) but is opposite in the region ($0<Q\le 1$) for diverse values of the Sutterby fluid parameter ($\eta$), as shown in Figure 6b. Figure 6c also shows the identical behavior of pressure rise for different values of the inclined angle ($\gamma$). In Figure 6d,e, we note the uniformly decreasing behavior of the pressure rise for numerous values of the Darcy number ($Da$) and cilia length ($\widetilde{ϵ}$) throughout the pumping region. Moreover, Figure 6f,g denote the increasing behavior of the pressure rise for higher values of the Forchheimer number ($F_r$) and eccentricity parameter ($\tilde{\omega }$), respectively.

4.6. Temperature Profiles

Figure 7a–e, are plotted to evaluate the effects of the Brinkman number ($Br$), Darcy number ($Da$), Sutterby fluid parameter ($\eta$), Forchheimer number, and temperature ratio on the temperature profile. It can be noticed in Figure 7a that as we increase the Brinkman number values, the temperature enriches throughout the entire region. As the Brinkman number is the product of the Prandtl number and the Eckert number, its higher values weaken the conduction process due to viscous dissipation, and the temperature profile gradually looks like a parabolic shape. The variation in the temperature profile is notable for diverse values of the Darcy number ($Da$) in Figure 7b. It has been found that less heat is generated for large values of ($Da$). This behavior is physically justified because when we progress from a porous medium toward a non-porous medium, there will be less friction between the fluid and the porous medium. Consequently, the temperature profile declined. Furthermore, in Figure 7c, a similar situation of temperature profile has been found for dissimilar values of the Sutterby fluid parameter ($\eta$). From Figure 7d,e, it can be observed that the temperature profile is enriched with a higher Forchheimer number ($F_r$) and temperature ratio ($\Gamma$), respectively.

4.7. Trapping Phenomena

In this subsection, we discuss the very fascinating phenomenon of trapped boluses of fluid with the help of streamlines, as depicted in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14. The streamlines have been plotted with the help of the stream function $\overline{\psi }$ and they are defined as:

$$\tilde{U}=\frac{\partial \overline{\psi }}{\partial \tilde{Y}},\tilde{V}=-\frac{\partial \overline{\psi }}{\partial \tilde{X}}.$$

(54)

This is the most fascinating mechanism in peristaltic phenomena. It is the creation of the fluid’s internally circulating bolus circumscribed by the streamlines. The confined bolus is pushed by the peristaltic waves. The formation of wall eddies as the width of the channel rises can be easily observed. When the peristaltic feature is weak, the eddies begin to disengage from the walls, relying only on kinematic principles. The obtained flow pattern is determined by the direction of wave propagation. As the parameters are changed, the pressure field varies, causing variations and the formation of multiple eddies. For instance, when the pressure field and the peristaltic motion converge in the same direction, the eddies grow and approach the centerline, giving birth to the trapping mechanism. In Figure 8, it can be observed that the boluses are initially expanded in size when the Hartman number ($M$) is elevated. In Figure 9 and Figure 10, conflicting behavior is observed for diverse values of the Sutterby fluid ($\eta$) and Forchheimer ($F_r$) numbers. In Figure 11, we observed that the number of trapped boluses ascended in size and number for diverse values of the Darcy number ($Da$). Figure 12 shows that the number of boluses diminished with the inclined angle ($\gamma$), but in Figure 13 there are no significant changes found in trapped boluses when the eccentricity ($\tilde{\omega }$) is prominent. Lastly, we can observe the number of boluses increases when the cilia length ($\widetilde{ϵ}$) is elevated, as depicted in Figure 14.

5. Conclusions

In this theoretical study, we consider the asymmetric tapered-channel geometry with ciliated walls. The study is based on the cilia-driven MHD flow of a Sutterby fluid with non-Darcy resistance as well as heat-transfer phenomena. After some suitable dimensionless parameters, the leading equations in dimensionless form have been found. Further, lubrication theory is introduced for mathematical simplifications. A very interesting procedure, the homotopy perturbation method (HPM), has been adopted for the solutions of non-dimensional governing equations with appropriate boundary conditions. Throughout the study, the computational software “Mathematica” has been utilized for the mathematical, numerical, and pictorial results. Some useful outcomes we achieved from this study are given below:

• An increasing behavior of the velocity, current density, and magnetic force function profiles is observed due to the elevated values of the Darcy and Forchheimer numbers, while the Sutterby parameter and Hartman number show the converse trends.
• The increasing values of the eccentricity parameter have almost negligible influences on the velocity, current density, and magnetic force function profiles.
• The velocity and current density profiles gradually slow down in the center for rising values of the inclined angle, while the opposite pattern in the case of the magnetic force function throughout the region has been observed.
• The increasing values of cilia length slow down the magnetic force function effects and the velocity and current density show some hindrance in the central regions but show different behavior near the walls.
• The magnitude of the magnetic force function grows more quickly when the magnetic Reynolds number is elevated.
• The current density magnitude is found to increase with higher values of the electric field and magnetic Reynolds number.
• Every graph of the induced magnetic field has dual behavior or two opposite trends for all physical parameters except the cilia length parameter.
• With increasing values of the Brinkman number, temperature ratio, and Forchheimer number, an increase in temperature is observed, while a decrease in the temperature magnitude is observed for both the Darcy and Sutterby parameters.
• From streamlined patterns, it can be seen that boluses diminish randomly for the Sutterby fluid parameter, inclined angle, and Hartman number but are enriched by the influence of the cilia length and Darcy number.
• No significant changes occur in streamlines for the eccentricity parameter, and boluses expand for elevated values of the Forchheimer number.
• In addition, we can also deduce the mathematical result for the Newtonian fluid by taking m = 0.