Heat and mass transfer characteristics in flow of bi-viscosity fluid through a curved channel with contracting and expanding walls: A finite difference approach

Abstract

This article investigates the heat and mass transfer in flow of bi-viscosity fluid through a porous saturated curved channel with sinusoidally deformed walls. The magnetic field and Joule heating effects are also taken into account. The equations describing the flow and heat/mass transfer are developed using curvilinear coordinates. A reduction of these equations is made based on lubrication approximation. The reduced linear ordinary differential equations are integrated numerically using an implicit finite difference scheme. It is observed that, the bi-viscosity fluid parameter, permeability parameter, and Hartmann number have analogous effects on the longitudinal velocity. Moreover, temperature of the fluid, heat coefficient, and mass concentration increase by increasing bi-viscosity fluid parameter, Brinkmann number, and Hartmann number. Further, mass concentration increases by increasing the rate of chemical reaction and bi-viscosity parameter. The size of circulating roll in lower half of the channel boosts up with larger variation of bi-viscosity parameter and permeability parameter. The flow patterns in the channel illustrating the effects of bi-viscosity parameter, permeability parameter, and Hartmann number are also displayed.

Introduction

Flow due to sinusoidally deformed vessel walls is commonly known as peristaltic flow. Peristaltic flows are caused by rhythmic contractions and relaxations of smooth muscles of the vessel. Peristalsis is a major mechanism found in the human body for physiological fluid transport. Examples of physiological fluids transported by peristaltic activity are urine, blood, chyme, spermatic fluid, etc. Modern machines which operate on the principle of peristalsis are heart-lung machine, diabetic pumps, roller and finger pumps, pharmacological delivery systems, etc.

The fluid dynamics of peristalsis are extensively investigated in the literature in past few decades. The first theoretical and experimental investigation on peristaltic flows was presented by Shapiro et al. ¹ based on the absence of curvature and inertial features. This prestigious investigation has been further extended for nonlinear fluid models,^2–9 curvature consequences,^10–14 magnetohydrodynamic prospective,^15,16 electro-osmotic effects, ¹⁷ and slip effects.^18–20 The important phenomena of heat and mass transfer have also been investigated by researchers for peristaltic transport as it involves motivated significance in area of laser therapy, hemodialysis, chemotherapy, and cryosurgery.^21–27 More recently peristaltic flow of nanofluids is an active area of research because of novel applications in micropumps, drug delivery systems, and pharmacological engineering. ²⁸

Due to multidisciplinary applications of non-Newtonian materials, researchers have concentrated on explanation of complex physical properties in recent days. It is well justified fact that all physical properties of non-Newtonian liquids cannot be completely predicted by intruding a single mathematical relation. Following such distinguish rheological features, bi-viscosity fluid characterizes to non-Newtonian fluids because of complex rheological properties. The relation for shear stress and deformation rate is extensively different from traditional viscous fluid model. The unique feature associated with this fluid model behaves like solid when yield stress is larger as compared to shear stress and deforms continuously when amount of stress becomes dominant with respect to yield stress. The motivated application encountered by bi-viscosity liquid reports physical applications in food processing, metallurgy and drilling processes. Biviscosity fluid materials include priceless applications in many biological fluids, blood, synovial liquids, plasma, lubricants, many pharmaceutical products, proteins, sewage sludge, china clay, and many more. Number of researchers investigated rheological impact of this model in various flow scenarios. For instance Dash et al. ²⁹ described the flow characteristics of bi-viscosity fluid through a homogeneous porous medium bounded by a circular tube. The mixed convection and variable thermal features for dissipative flow of bi-viscosity nanofluid have been intended by Hayat et al. ³⁰ The effects of chemical reaction on flow of biviscosity fluid are thoroughly studied by Shehzad et al. ³¹ and Reddy et al. ³² Further, studies regarding flows of bi-viscosity fluids can be found in references.^33–35

The porous medium is a stuff that encompasses spaces between solid areas through which fluid can pass. In recent years, considerable interest has been noticed to examine the flow through saturated porous space due to its fundamental applications in petroleum engineering, soil sickness, geophysical processes, and medical applications. Some fundamental applications for saturated porous media flow include filtration and extraction of oil from wells, water draining applications for irrigation, oil flow in porous rocks, energy extraction associated with geo-thermal systems, the fluid flow in ion-exchange beds, and oil-spills cleaning. The foremost contribution on flows through porous media was made by Brinkman, ³⁶ who calculated the viscous force exerted by a flowing fluid through a porous media. Peristaltic flows through porous medium are also an active area of research. Several authors contributed to this area. Mention may be made to the work of El Shehawey, ³⁷ Mekheimer, ³⁸ Kothandapani and Srinivas, ³⁹ and Hayat et al. ⁴⁰ in various scenarios. Javed et al. ⁴¹ studied the effect of mixed convection through micropolar nanofluid under the impact of magnetohydrodynamics. In a separate study Javed et al. ⁴² examined the effect of heat transfer in water–Cu nanofluid saturated porous medium through trapezoidal cavities under the influence of magnetic field. Javed et al. ⁴³ investigated effect of natural convection in a ferro-fluid filled square cavity in the presence of uniform magnetic field. Mehmood et al. ⁴⁴ investigated the combined effect of MHD-Mixed convection flow in a lid driven trapezoidal cavity in presence of heated bottom wall. In a separate study Mehmood ⁴⁵ simulated numerically the effect of mixed thermomagnetic convection through two lid-driven entrapped trapezoidal cavities enclosing ferrofluid saturated porous medium. Mehmood and Javed ⁴⁶ considered the impact of MHD on energy transmission through ferrofluid saturating porous medium contained in a lid-driven trapezoidal container with triangular notched heater configuration.

The phenomenon of heat transfer in flowing materials involves a variety of applications in many chemical and mechanical industries, thermal engineering, electronic devices, nuclear thermal-hydraulic processes, etc. Owing to such interesting applications, many authors analyzed the heat transfer characteristics in flow of different fluids. For instance, Belhocine ⁴⁷ examined the heat transfer analysis in flow developed flow configured by a circular tube. The work of Belhocine and Omar ⁴⁸ deals with study of convective heat transfer due to circular tube where analytical solution were computed for formulated flow problem. In another useful contribution, Belhocine and Abdullah ⁴⁹ deals the flow and heat transfer of viscous fluid by using constant wall temperature assumptions, usually called the Graetz problem. The thermal boundary layer analysis for modified Levèque problem has been performed by Belhocine and Abdullah. ⁵⁰ Joule heating appears due to the resistance offered to electric current passing through some conductive material. There are number of systems in which Joule heating effect has vital role such as dielectrophoretic trapping, electric fuses, PCR reactors, hot plate, microvalves for fluid control, electric heaters and stoves, thermistors and soldering irons, etc. The effect of Joule heating in peristalsis has also been investigated. For detail see references.^51–55

From literature cited above, some investigations regarding the peristaltic flow of bi-viscosity model encountered by a porous saturated straight vessel have been performed but no continuation has been utilized to explore the peristaltic phenomenon in a curved porous-saturated channel. Such motivated applications insist us to report the heat and mass transportation in magnetized bi-viscosity fluid configured by a curved channel in presence of porous medium. The objective here is to report the consequences of formulated dimensionless parameters namely Brinkmann number, bi-viscosity fluid constant, radius of curvature, Hartmann number, and porous medium permeability on combined heat/mass transportation. It is remarked that although variety of studies are available in the literature regarding the peristaltic flow of non-Newtonian fluids in curved channel but analysis regarding peristaltic transport of magnetized bi-viscosity fluid in presence of Joule heating, porous medium and heat and mass transfer phenomenon has not be investigated yet. The aim of current contribution is to fulfill this research gap. The governing equations for the concerned flow model are constituted under the assumption of dominant viscous effects as compared to the inertial features. The solution for formulated flow problem has been suggested by employing implicit finite difference scheme with desirable accuracy. Later on, a comprehensive graphical analysis for each flow parameters is reported with relevant physical consequences.

Modeling of flow problem

We assume a curved porous channel having width $2w$ coiled a circle with center $O$ along with radius $R_0.$ The curved channel is filled with incompressible bi-viscosity liquid. The flow has been induced sinusoidal type channel walls deformation where impact of magnetic force is employed radially. Let $T_0$ and $C_0$ are being upper wall temperature and concentration respectively. The upper wall maintained constant temperature and concentration represented by $T_1$ and $C_1,$ respectively. The flow pattern has been followed by using curvilinear coordinate system $(R,\chi ,Z)$ in which $R$ is specified in radial direction, $\chi$ is taken the direction fluid flow while $Z$ is taken perpendicular plane associated with $R$ and $\chi$ The physical problem is described by sketching Figure 1. Equations representing both walls shape are illustrated as follows^12,22,23,25:

$$H_1(\chi ,t)=w+a\sin \left(\left(\frac{2\pi }{{\lambda }^{*}}\right)(\chi -\mathrm{ct})\right),\mathrm{Upper}\mathrm{wall},$$

(1)

$$H_2(\chi ,t)=-w-b\sin \left(\left(\frac{2\pi }{{\lambda }^{*}}\right)(\chi -\mathrm{ct})\right),\mathrm{Lower}\mathrm{wall},$$

(2)

where $t,c,{\lambda }^{*}$ and $(a,b)$ are being time, wave speed, wavelength and amplitudes, respectively. For bi-viscosity fluid, the governing flow equations regarding heat/mass transfer phenomenon are given by^12,23,25:

$$\nabla ·\mathbf{U}=0,\left(ContinuityEquation\right)$$

(3)

$$\rho \frac{d\mathbf{U}}{dt}=\nabla ·\tau -\frac{\mu }{k^{\ast }}\mathbf{U}+\mathbf{J}\times \mathbf{B},\left(MomentumEquation\right)$$

(4)

$$\rho c_p\frac{dT}{dt}=k{\nabla }^{2}T+\mu \mathrm{\Phi }+\frac{\sigma B^{\ast 2}{\tilde{R}}^2}{{\left(R+\tilde{R}\right)}^2}{U}^2,\left(EnergyEquation\right)$$

(5)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{dC}{\mathrm{dt}}=D{\nabla }^{2}C+\frac{D{K}_T}{T_m}{\nabla }^{2}T-k_{1}C,}\\ \multicolumn{1}{c}{\left(\mathrm{Mass}\mathrm{Concentration}\mathrm{Equation}\right)}\end{array}$$

(6)

where $\mathbf{\tau }$ notify the Cauchy stress tensor, $\mu$ determine the viscosity, $\beta$ symbolized the bi-viscosity fluid parameter, T is for temperature, $k$ is the thermal conductivity, $c_p$ noted specific heat at constant pressure, $\mathbf{U}$ is the velocity, $k^{*}$ is the permeability of porous medium, $C$ is the mass concentration, $\rho$ expresses the fluid density, $D$ mass diffusivity coefficient, $K_{\mathit{T}}$ is thermal diffusivity, $T_m$ depicts the mean temperature, $\mathrm{\Phi }$ represents the relation for dissipation function, $\sigma$ is the Stefan-Boltzman constant while $k_1$ is the rate of chemical reaction. The viscosity of a fluid is a measure of its resistance to deformation at a given rate while porosity of the medium is a measure of the void spaces in a material.

Figure 1.

Geometry of the problem.

For bi-viscosity fluid model ²⁹

$$\overline{\tau }=\{\begin{array}{c}({\mu }_b+\frac{P_y}{\sqrt{2\mathrm{\Pi }}}){\overline{A}}_1,\mathrm{when}\mathrm{\Pi }>{\mathrm{\Pi }}_C\hfill \\ \multicolumn{1}{c}{({\mu }_b+\frac{P_y}{\sqrt{2\mathrm{\Pi }}}){\overline{A}}_1,\mathrm{when}\mathrm{\Pi }<{\mathrm{\Pi }}_C}\end{array}$$

(7)

where $P_y$ is the yield stress of the fluid given by

$$P_y=\frac{{\mu }_b\sqrt{2\mathrm{\Pi }}}{\beta },$$

(8)

in which ${\mu }_b$ is the plastic dynamic viscosity, $\mathrm{\Pi }$ is second invariant of ${\overline{A}}_1$ , ${\mathrm{\Pi }}_C$ is the critical value based on non-Newtonian model. For bi-viscosity fluid flow, where $\mathrm{\Pi }>{\mathrm{\Pi }}_C$ , it is possible to write

$$\mu ={\mu }_b+\frac{P_y}{\sqrt{2\mathrm{\Pi }}}.$$

(9)

Using equations (8) and (9), we get stress tensor for bi-viscosity fluid flow

$$\overline{\tau }={\mu }_b(1+\frac{1}{\beta }){\overline{A}}_1,\mathrm{when}\mathrm{\Pi }>{\mathrm{\Pi }}_C.$$

(10)

In view of equation (10), equation (4) becomes ²⁹

$$\rho \frac{d\mathbf{U}}{\mathrm{dt}}={\mu }_b(1+\frac{1}{\beta })\nabla ·{\overline{A}}_1-\frac{\mu }{k^{*}}\mathbf{U}+\mathbf{J}\times \mathbf{B}.$$

(11)

The form of the radial magnetic field B is given by ²²

$$\mathbf{B}=\left(\frac{B^{*}\tilde{R}}{R+\tilde{R}}\right){\mathit{e}}_{\mathit{R}},$$

(12)

where $B^{*}$ is the characteristic magnetic induction in the limit $\tilde{R}\to \infty$ and ${\mathit{e}}_{\mathit{R}}$ is the unit vector along radial direction. It is remarked that expression for magnetic force defined in equation (12) is solenoidal and satisfies Maxwell’s equations.

Using equation (12), the term $\mathbf{J}\times \mathbf{B}$ in equation (4) is given by ¹³

$$\mathbf{J}\times \mathbf{B}=-\left(\frac{\sigma {B^{*}}^2{U}_2{\tilde{R}}^2}{{(R+\tilde{R})}^2}\right){\mathit{e}}_{\mathit{X}},$$

(13)

where ${\mathit{e}}_{\mathit{X}}$ azimuthal direction unit vector.

Assuming $\mathbf{U}=[U_1(\chi ,R,t),U_2(\chi ,R,t),0],T=T(\chi ,R,t),C=C(\chi ,R,t),$ equations (3) to (6) yield^12,22,23,25

$$\frac{\partial }{\partial R}\left\{(R+\tilde{R})U_1\right\}+\tilde{R}\frac{\partial U_2}{\partial \chi }=0,$$

(14)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial U_1}{\partial t}+(U_1.\nabla )U_1-\frac{{U_2}^2}{R+\tilde{R}}=-\frac{1}{\rho }\frac{\partial P}{\partial R}+\frac{{\mu }_b}{\rho }(1+\frac{1}{\beta })}\\ \multicolumn{1}{c}{[{\nabla }^2{U}_1-\frac{U_1}{{(R+\tilde{R})}^2}-\frac{2\tilde{R}}{{(R+\tilde{R})}^2}\frac{\partial U_2}{\partial \chi }]-\frac{\nu }{k^{*}}{U}_1,}\end{array}$$

(15)

$$\begin{array}{c}\frac{\partial U_2}{\partial t}+(U_1.\nabla )U_2+\frac{U_1{U}_2}{R+\tilde{R}}=-\frac{\tilde{R}}{\rho (R+\tilde{R})}\frac{\partial P}{\partial \chi }+\frac{{\mu }_b}{\rho }(1+\frac{1}{\beta })\hfill \\ \multicolumn{1}{c}{[{\nabla }^2{U}_2-\frac{U_2}{{(R+\tilde{R})}^2}+\frac{2\tilde{R}}{{(R+\tilde{R})}^2}\frac{\partial U_1}{\partial \chi }]-\frac{\nu }{k^{*}}{U}_2}\\ -\frac{\sigma {B^{*}}^2{\tilde{R}}^2}{\rho {(R+\tilde{R})}^2}{U}_2,\hfill \end{array}$$

(16)

$$\rho c_p[\frac{\partial T}{\partial t}+(U_1.\nabla )T]=k{\nabla }^{2}T+\mu \mathrm{\Phi }+\frac{\sigma {B^{*}}^2{\tilde{R}}^2{U}_2^2}{{(R+\tilde{R})}^2},$$

(17)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial C}{\partial t}+(U_1.\nabla )C=D{\nabla }^{2}C+\frac{D{K}_T}{T_m}}\\ \multicolumn{1}{c}{(\frac{{\partial }^{2}T}{\partial R^2}+\frac{\tilde{R}}{R+\tilde{R}}\frac{\partial T}{\partial R}+{\left(\frac{\tilde{R}}{R+\tilde{R}}\right)}^2\frac{{\partial }^{2}T}{\partial {\chi }^2})-k_{1}C,}\end{array}$$

(18)

where

$$U_1.\nabla =U_1\frac{\partial }{\partial R}+\frac{\tilde{R}{U}_2}{R+\tilde{R}}\frac{\partial }{\partial \chi },$$

(19)

$${\nabla }^2=\frac{1}{R+\tilde{R}}\frac{\partial }{\partial R}\left\{(R+\tilde{R})\frac{\partial }{\partial R}\right\}+{\left(\frac{\tilde{R}}{R+\tilde{R}}\right)}^2\frac{{\partial }^2}{\partial {\chi }^2},$$

(20)

$$\begin{array}{l}\mathrm{\Phi }=[2{\left(\frac{\partial U_1}{\partial R}\right)}^2+\left(\frac{\tilde{R}}{R+\tilde{R}}\frac{\partial U_1}{\partial \chi }-\frac{U_2}{R+\tilde{R}}\right)\\ \left(\frac{\partial U_2}{\partial R}+\frac{\tilde{R}}{R+\tilde{R}}\frac{\partial U_1}{\partial \chi }-\frac{U_2}{R+\tilde{R}}\right)+\frac{\partial U_2}{\partial R}\\ [\frac{\partial U_2}{\partial R}+\frac{\tilde{R}}{R+\tilde{R}}\frac{\partial U_1}{\partial \chi }-\frac{U_2}{R+\tilde{R}}]+2\left(\frac{\tilde{R}}{R+\tilde{R}}\frac{\partial U_2}{\partial \chi }+\frac{U_1}{R+\tilde{R}}\right)].\end{array}$$

(21)

The boundary conditions associated with equations (14) to (18) are ²³

$$U_2=0,U_1=\frac{\partial H_1}{\partial t},T=T_o,C=C_o\mathrm{at}R=H_1(\chi ,t),$$

(22)

$$U_2=0,U_1=\frac{\partial H_2}{\partial t},T=T_1,C=C_1\mathrm{at}R=H_2(\chi ,t).$$

(23)

In order to shift from fixed frame $(R,\chi )$ to wave frame $(r,x)$ , we employ the transformations^22,23

$$\begin{array}{c}\multicolumn{1}{c}{x=\chi -ct,r=R,p=P,u_1=U_1,u_2}\\ \multicolumn{1}{c}{=U_2-c,T=T.}\end{array}$$

(24)

The governing equations are transformed to wave frame by using above transformations. These equations after defining the dimensionless variables^22,23

$\overline{x}=\frac{2\pi }{{\lambda }^{*}}x,\eta =\frac{r}{a_1},$ ${\overline{u}}_1=\frac{u_1}{c_1},$ ${\overline{u}}_2=\frac{u_2}{c_1},\overline{p}=\frac{2\pi {a_1}^2}{{\lambda }^{*}{\mu }_{b}c}p,$ $\theta =\frac{T-T_1}{T_0-T_1},$ $\varphi =\frac{C-C_1}{C_0-C_1},$ $\delta =\frac{2\pi a_1}{{\lambda }^{*}},$ $\mathrm{Re}=\frac{\rho c{a}_1}{{\mu }_b},$ and invoking lubrication approximations ( $\delta \approx 0,\mathrm{Re}\approx 0$ ) reduce to:

$$\frac{\partial p}{\partial \eta }=0,$$

(25)

$$\begin{array}{c}-\frac{\partial p}{\partial x}+\frac{1}{\gamma }\{-\frac{\partial }{\partial \eta }\left((\eta +\gamma )\frac{{\partial }^2\psi }{\partial {\eta }^2}\right)-\frac{1}{\eta +\gamma }(1-\frac{\partial \psi }{\partial \eta })\}\hfill \\ \multicolumn{1}{c}{-\frac{\eta +\gamma }{K^{*}\gamma }(1-\frac{\partial \psi }{\partial \eta })}\\ -\frac{\gamma H{a}^2}{\eta +\gamma }(1-\frac{\partial \psi }{\partial \eta })=0,\hfill \end{array}$$

(26)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2\theta }{\partial {\eta }^2}+\frac{1}{(\eta +\gamma )}\frac{\partial \theta }{\partial \eta }}\\ \multicolumn{1}{c}{+\mathrm{Br}{(-\frac{1}{\eta +\gamma }(1-\frac{\partial \psi }{\partial \eta })-\frac{{\partial }^2\psi }{\partial {\eta }^2})}^2}\\ \multicolumn{1}{c}{+\frac{\mathrm{BrH}{a}^2{\gamma }^2}{{(\eta +\gamma )}^2}{(1-\frac{\partial \psi }{\partial \eta })}^2=0,}\end{array}$$

(27)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2\varphi }{\partial {\eta }^2}+\frac{1}{(\gamma +\eta )}\frac{\partial \varphi }{\partial \eta }-R_c\varphi }\\ \multicolumn{1}{c}{=-\mathrm{SrSc}(\frac{{\partial }^2\theta }{\partial {\eta }^2}+\frac{1}{(\gamma +\eta )}\frac{\partial \theta }{\partial \eta })+R_c.}\end{array}$$

(28)

In above equations Re, δ, $\gamma$ , $K^{*}$ , and $R_c$ represent the Reynolds number, the wave number, the dimensionless radius of curvature, dimensionless permeability parameter and dimensionless rate of chemical reaction respectively. The stream function $\psi$ and velocity components $u_1$ and $u_2$ are related through the expressions ²²

$$u_1=\delta \frac{\gamma }{\eta +\gamma }\frac{\partial \psi }{\partial x},u_2=-\frac{\partial \psi }{\partial \eta }.$$

(29)

After eliminating pressure between equations (25) and (26), we get

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\partial }{\partial \eta }[\frac{1}{\gamma }(\frac{\partial }{\partial \eta }\{(\eta +\gamma )\frac{{\partial }^2\psi }{\partial {\eta }^2}\}+\frac{1}{\eta +\gamma }(1-\frac{\partial \psi }{\partial \eta }))}\\ \multicolumn{1}{c}{+\frac{\eta +\gamma }{K^{*}}\gamma (1-\frac{\partial \psi }{\partial \eta })+\frac{\gamma {\mathrm{Ha}}^2}{\eta +\gamma }(1-\frac{\partial \psi }{\partial \eta })]=0.}\end{array}$$

(30)

The transmuted boundary conditions (22) and (23) are

$$\begin{array}{c}\multicolumn{1}{c}{\psi =-\frac{q}{2},\frac{\partial \psi }{\partial \eta }=1,\theta =0,\varphi =0,\mathrm{at}\eta =h_1}\\ \multicolumn{1}{c}{=1+\lambda \sin x,}\end{array}$$

(31)

$$\begin{array}{c}\multicolumn{1}{c}{\psi =\frac{q}{2},\frac{\partial \psi }{\partial \eta }=1,\theta =1,\varphi =1,\mathrm{at}\eta }\\ \multicolumn{1}{c}{=h_2=-1-\lambda \prime \sin x,}\end{array}$$

(32)

where $\lambda =a/w$ and $\lambda \prime =b/w$ stands for amplitude ratios. We need to compute the solution of equations (27), (28), and (30) by using boundary conditions in equations (31) and (32). Following expressions are initiated for pressure per wavelength $\mathrm{\Delta }p,$ local Nusselt number at both the wall $z_i(i=1,2)$ and Sherwood number $S{h}_i(i=1,2)$ ^15,22,25

$$\mathrm{\Delta }p=\underset{0}{\overset{2\pi }{\int }}\frac{\mathrm{dp}}{\mathrm{dx}}\mathrm{dx},$$

(33)

$$z_i=\frac{\partial h_i}{\partial x}{\frac{\partial \theta }{\partial \eta }|}_{\eta =h_i},i=1,2.$$

(34)

$$\mathrm{Sh}=\frac{\partial h_i}{\partial x}{\frac{\partial \varphi }{\partial \eta }|}_{\eta =h_i},i=1,2.$$

(35)

Finite difference scheme

The set of equations (27), (28), and (30) subject to boundary conditions given in equations (31) and (32) are quite complex and highly nonlinear for which analytical solution is not possible. On this end, we follow finite difference approach to simulate the numerical solution of such highly complex equations. This method successfully tacked nonlinear differential equations with excellent accuracy. This method is also sufficient to compute the solution of problem associated with curved boundaries and complex boundary conditions. The numerical computations are performed by using FORTRAN software. This method is based on following four steps:

• (i) First step is to construct an iterative technique to allow the new nonlinear boundary value problem (BVP) to be transformed into a linear one at the (m + 1)th iterative stage. For this specific problem, the subsequent iterative technique is suggested:

$$\begin{array}{c}\frac{{\partial }^4{\psi }^{(m+1)}}{\partial {\eta }^4}+\frac{2}{(\eta +\gamma )}\frac{{\partial }^3{\psi }^{(m+1)}}{\partial {\eta }^3}\hfill \\ \multicolumn{1}{c}{-\{\frac{1}{{(\eta +\gamma )}^2}+\frac{{\gamma }^{2}H{a}^2}{{(\eta +\gamma )}^2}+\frac{1}{K^{*}}\}\frac{{\partial }^2{\psi }^{(m+1)}}{\partial {\eta }^2}}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{+\{\frac{1}{{(\eta +\gamma )}^3}-\frac{{\gamma }^{2}H{a}^2}{{(\eta +\gamma )}^3}-\frac{1}{K^{*}(\eta +\gamma )}\}}\\ \frac{\partial {\psi }^{(m+1)}}{\partial \eta }-\frac{1}{{(\eta +\gamma )}^3}+\frac{{\gamma }^{2}H{a}^2}{{(\eta +\gamma )}^3}+\frac{1}{K^{*}(\eta +\gamma )}=0.\hfill \end{array}$$

(36)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2{\theta }^{(m+1)}}{\partial {\eta }^2}+\frac{1}{\gamma +\eta }\frac{\partial {\theta }^{(m+1)}}{\partial \eta }}\\ \multicolumn{1}{c}{=-\mathrm{Br}{(-\frac{1}{\gamma +\eta }(1-\frac{\partial {\psi }^{(m)}}{\partial \eta })-\frac{{\partial }^2{\psi }^{(m)}}{\partial {\eta }^2})}^2}\\ \multicolumn{1}{c}{-\frac{\mathrm{BrH}{a}^2{\gamma }^2}{{(\gamma +\eta )}^2}{(1-\frac{\partial {\psi }^{(m)}}{\partial \eta })}^2}\end{array}$$

(37)

$$\begin{array}{c}\multicolumn{1}{c}{\frac{{\partial }^2{\varphi }^{(m+1)}}{\partial {\eta }^2}+\frac{1}{\eta +\gamma }\frac{\partial {\varphi }^{(m+1)}}{\partial \eta }-R_c{\varphi }^{(m+1)}}\\ \multicolumn{1}{c}{=-\mathrm{SrSc}(\frac{{\partial }^2{\theta }^{(m)}}{\partial {\eta }^2}+\frac{1}{(\gamma +\eta )}\frac{\partial {\theta }^{(m)}}{\partial \eta })+R_c,}\end{array}$$

(38)

$$\begin{array}{c}\multicolumn{1}{c}{{\psi }^{m+1}=-\frac{q}{2},\frac{\partial {\psi }^{m+1}}{\partial \eta }=1,{\theta }^{(m+1)}=0,{\varphi }^{(m+1)}=0,}\\ \multicolumn{1}{c}{\mathrm{at}\eta =h_1,}\end{array}$$

(39)

$$\begin{array}{c}\multicolumn{1}{c}{{\psi }^{(m+1)}=\frac{q}{2},\frac{\partial {\psi }^{(m+1)}}{\partial \eta }=1,{\theta }^{(m+1)}=1,{\varphi }^{(m+1)}=1,}\\ \multicolumn{1}{c}{\mathrm{at}\eta =h_2.}\end{array}$$

(40)

Here, the index (m) shows the iterative step. It is now clear that above BVP is linear in ${\psi }^{(m+1)}$ .

• (ii) In this step, we insert finite difference approximations of ${\psi }^{(m+1)}$ , ${\theta }^{(m+1)}$ , ${\varphi }^{(m+1)}$ and their derivatives into equations (36) to (38). In this way, we get a system of linear algebraic equations at each iterative step.

• (iii) In third step, the system of algebraic equations obtained in previous step are solved at each cross-section to get numerical results of ${\psi }^{(m+1)}$ , ${\theta }^{(m+1)}$ , and ${\varphi }^{(m+1)}$ . Evidently, suitable initial guesses are required for ${\psi }^{\left(m\right)}$ , ${\theta }^{\left(m\right)}$ , and ${\varphi }^{\left(m\right)}$ at each cross-section to start the iterative procedure. The iterative procedure at cross-section is carried out until a convergent solution is reached.

• (iv) In this step, we employed method of successive under-relaxation to get rapid convergent solution. In this method the values of ${\tilde{\psi }}^{(m+1)},$ ${\tilde{\theta }}^{(m+1)}$ , and ${\tilde{\varphi }}^{(m+1)},$ at (m+1)th iterative step are refined as

$${\psi }^{(m+1)}={\psi }^{(m)}+\tau ({\tilde{\psi }}^{(m+1)}-{\psi }^{(m)}),$$

(41)

$${\theta }^{(m+1)}={\theta }^{(m)}+\tau ({\tilde{\theta }}^{(m+1)}-{\theta }^{(m)}),$$

(42)

$${\varphi }^{(m+1)}={\varphi }^{(m)}+\tau ({\tilde{\varphi }}^{(m+1)}-{\varphi }^{(m)}),$$

(43)

where $\tau$ is under relaxation parameter.

Usually, $\tau$ is chosen small for rapid convergence. In present computation, the iterative procedure is terminated after achieving the values of $\psi ,\theta$ , and $\varphi$ convergent to ${10}^{-8}.$ Some authors have already been used for instance, Wang et al., ¹⁵ Hayat et al., ¹⁹ Ali et al., ²³ and Abbas et al. ²⁷ The complete iterative process in represented in Figure 2.

Figure 2.

The flow chart of the numerical method.

Verification of solution

The numerical solution computed in previous section is compared for velocity profile with work of Ali et al. ²² as a limiting case $(K^{*}\to \infty ,\mathrm{Ha}=0,\gamma =2)$ . On this end, Figure 3 is prepared for which superimposed line show the solution reported by Ali et al. ²² while solid line represents the current numerical solution. An excellent accuracy between solutions is noted.

Figure 3.

Comparison of solution for velocity: superimposed line denotes the work of Ali et al. ¹¹ while solid line represents numerical solution.

Results and discussion

In this section, it is intended to analyze the important features of peristaltic motion such as velocity, pressure rise, temperature distribution, mass concentration and trapping for various values of the curvature parameter $(\gamma )$ , bi-viscosity parameter $(\beta ),$ Brinkman number $(\mathrm{Br})$ , rate of chemical action ( $R_c$ ), permeability parameter $(K^{*})$ , and Hartmann number $(Ha).$ The heat transfer coefficient at upper wall is also discussed. To this end, Figures 4–28 are plotted.

Figure 4.

$u_2\left(\eta \right)$ for $\beta$ with $\gamma$  = 2.5, $K^{*}$  = 0.2, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 5.

$u_2\left(\eta \right)$ for $K^{*}$ with $\gamma$  = 2.5, $\beta$  = 0.2, Ha = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 6.

$u_2\left(\eta \right)$ for $\mathrm{Ha}$ with $\gamma$  = 2.5, $K^{*}$  = 0.2, $\beta$  = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 7.

$u_2\left(\eta \right)$ for $\gamma$ with $\beta$  = 0.2, $K^{*}$  = 0.2, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 8.

$\mathrm{\Delta }p$ for $\beta$ with $\gamma$  = 2.5, $K^{*}$  = 0.2, Ha =0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 9.

$\mathrm{\Delta }p$ for $K^{*}$ with $\gamma$  = 2.5, Ha = 0.5, $\beta$  = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 10.

$\mathrm{\Delta }p$ for $\mathrm{Ha}$ with $\gamma$  = 2.5, $K^{*}$  = 0.2, $\beta$  = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 11.

$\mathrm{\Delta }p$ for $\gamma$ with $K^{*}$  = 0.2, $\beta$  = 0.2, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 12.

$\theta \left(\eta \right)$ for $\beta$ with Br = 0.2, $\gamma$  = 2, $K^{*}$  = 1, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 13.

$\theta \left(\eta \right)$ for Br with $\beta$ =0.2, $\gamma$  = 2, $K^{*}$  = 1.5, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 14.

$\theta \left(\eta \right)$ for $K^{*}$ with Br = 0.2, $\gamma$  = 2, Ha = 0.5, $\beta$  = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 15.

$\theta \left(\eta \right)$ for Ha with Br = 0.2, $\gamma$  = 2, $K^{*}$  = 1.5, $\beta$  = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 16.

z at upper wall for $\beta$ with Br = 0.2, $\gamma$  = 2, $K^{*}$  = 1, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 17.

z at upper wall for Br with $\beta$  = 0.2, $\gamma$  = 2, $K^{*}$  = 1.5, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 18.

z at upper wall $K^{*}$ with Br = 0.2, $\gamma$  = 2, Ha = 0.5, $\beta$  = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 19.

z at upper wall for $\mathrm{Ha}$ with Br = 0.2, $\gamma$  = 2, $K^{*}$  = 1.5, $\beta$  = 0.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 20.

ϕ $\left(\eta \right)$ for $\beta$ with Br = 2, $\gamma$  = 2, $R_c=0.2,$ Sr = 1.5, Sc = 1.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 21.

ϕ $\left(\eta \right)$ for Br with $\beta$  = 0.4, $\gamma$  = 2, $R_c=0.2,$ Sr = 1.5, Sc = 1.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 22.

ϕ $\left(\eta \right)$ for $R_c$ with $\beta$  = 0.4, Br = 2, $\gamma$  = 2, Sr = 1.5, Sc = 1.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 23.

ϕ $\left(\eta \right)$ for Ha with $\beta$  = 0.4, Br = 2, $\gamma$  = 2, $R_c$  = 0.2, Sr = 1.5, Sc = 1.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 24.

ϕ $\left(\eta \right)$ for Sc with $\beta$  = 0.4, Br = 2, $\gamma$  = 2, $R_c=0.2,$ Sr = 1.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 25.

ϕ $\left(\eta \right)$ for Sr with $\beta$  = 0.4, Br = 2, $\gamma$  = 2, $R_c=0.2,$ Sc = 1.2, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 26.

Streamlines in wave frame for (a) $\beta =0.1,$ (b) $\beta =0.5,$ (c) $\beta =1$ and (d) $\beta =1.5$ . The other parameters chosen are $\gamma$  = 2, $K^{*}$  = 1, Ha = 0.5, $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 27.

Streamlines in wave frame for (a) $K^{*}=0.1,$ (b) $K^{*}=0.15,$ (c) $K^{*}=0.2,$ and (d) $K^{*}=0.3.$ The other parameters chosen are $\gamma$  = 2, Ha = 0.5, $\beta =0.5,$ $\lambda =0.4,$ and $\Theta$  = 1.5.

Figure 28.

Streamlines in wave frame for (a) $\mathrm{Ha}=1,$ (b) $\mathrm{Ha}=5,$ (c) $\mathrm{Ha}=10,$ and (d) $\mathrm{Ha}=15.$ The other parameters chosen are $\gamma$  = 2, $K^{*}$  = 0.2, $\beta =0.5,$ $\lambda =0.4,$ and $\theta$  = 1.5.

The axial velocity distribution for some specific values of bi-viscosity parameter $(\beta ),$ permeability parameter $(K^{*})$ , Hartmann number $(Ha)$ , and curvature parameter $(\gamma )$ is shown in Figures 4–7, respectively. Figure 4 depicts that velocity profile increases by increasing bi-viscosity parameter. The truncation in velocity component $u_2$ for permeability parameter $\left(K^{*}\right)$ is verified by illustrating Figure 5. The weaker porous medium consequences are deliberated due to lower $\left(K^{*}\right)$ values in contrast to higher values which denote the case when resistance corresponds to porous matrix is more dominant. Moreover, the amplitude of velocity due to variation of $\left(K^{*}\right)$ slow down the velocity amplitude the velocity is shift from lower channel wall to central line $(\eta =0).$ Figure 6 exposes the consequences of applied magnetic field on the flow velocity. Here it is quite obvious that flow velocity exhibits boundary layer character for larger values of Hartmann number. In fact for large values of Hartmann number, the disturbance in flow velocity is confined in thin layer near both upper and lower walls. The fluid outside the boundary layer moves with a velocity which varies linearly with radial coordinate $\eta$ . In contrast, the fluid outside the boundary layers moves with constant velocity in a straight channel. The consequences of radius of curvature $(\gamma )$ on velocity have been predicted in Figure 7. The velocity distribution is observed to be asymmetric about $\eta =0$ associated with lower values of $\gamma .$ It is further noticed that for larger values of $\gamma ,$ velocity recovers its symmetry about $\eta =0.$ This is expected because for large values of $\gamma$ channel becomes straight. The graphical visualization of pressure rise per wavelength against diverse flow parameters has been reported in Figures 8–11. Figure 8 shows the effects of bi-viscosity parameter $\beta$ on $\mathrm{\Delta }p$ . It is observed that $\mathrm{\Delta }p$ is decreasing by increasing $\beta$ in the peristaltic pumping region. The profile of pressure rise per wavelength for different values of $K^{*}$ (permeability parameter) is shown in Figure 9. A declining trend in $\mathrm{\Delta }p$ due to enhanced porous medium permeability has been seen in the region of peristaltic pumping $(\Theta >0,\mathrm{\Delta }p>0).$ Thus, the normal fluid flow inside the channel impedes due to porous medium and peristalsis play a valuable role against the pressure rise to maintain the same flux as in the case of clear medium inside the channel. This fact reduced the pumping efficiency subsequently. Figures 10 and 11 are plotted to show the variation of $\mathrm{\Delta }p$ against dimensionless mean flow rate $\Theta$ for various values of Ha and $\gamma$ , respectively. The maximum variation in Ha and $\gamma$ , an improved pressure rise per wavelength has been resulted in the pumping region. However, for case of free pumping ( $\mathrm{\Delta }p$  = 0) and co-pumping case $(\Theta >0,\mathrm{\Delta }p<0),$ this situation is quite different. A lower trend in $\mathrm{\Delta }p$ is observed due to Ha and $\gamma$ . Figure 11 revealed that change in $\mathrm{\Delta }p$ is maximum in case of straight channel as compared to curved configuration up to a certain values of $\Theta$ . However, opposite observations are referred above this critical numerical value. Similar reverse styles encountered in region of co-pumping and free pumping. The alteration in radial temperature distribution inside the channel against flow parameters like bi-viscosity parameter $(\beta ),$ permeability parameter $(K^{*})$ Hartmann number (Ha) and Brinkmann number $\left(\mathrm{Br}\right)$ are executed in Figures 12–15. Figures 12, 13, and 15, respectively, report that $\theta (\eta )$ increases by increasing $\beta ,\mathrm{Br}$ and Ha inside the channel. In contrast, Figure 14 shows a decrease in temperature with increasing permeability parameter. It is due to the fact that heat transportation appeared from the walls to the fluid is reduced for clear medium and as a result, a decrease in temperature is noted inside the channel. Figures 16–19 indicate change in heat transfer coefficient $z$ at the upper channel wall for $K^{*},\beta ,\mathrm{Br},$ and Ha. The distribution of $z$ oscillates periodically which is associated with the oscillatory nature of wall surfaces. The results reveal that a relatively improved oscillatory profile is referred to the maximum values of $\beta ,\mathrm{Br}$ and Ha. However, a declining periodically oscillating phenomenon is evaluated with enlarge values of $K^{*}.$

The mass concentration variation inside the channel for several values of bi-viscosity parameter, Brinkmann number, rate of chemical reaction, Hartmann number, Schmidt number and Soret number are shown in Figures 20–25. It is observed that mass concentration inside the channel increases with increasing each of these parameters. The streamlines of flow inside the channel for different values of bi-viscosity parameter $(\beta ),$ permeability parameter $(K^{*})$ , Hartmann number (Ha) and curvature parameter $(\gamma )$ are shown in Figures 26–29. The aim of these curves is to highlight the trapping phenomenon where fluid volume known as bolus is trapped within closed streamlines. Figure 26 shows two nearly symmetric circulating rolls in lower and upper halves of the channel. However, lower roll increases in size while the upper one shrinks with increasing $\beta .$ Figure 27 illustrates the effects of $K^{*}$ on trapping phenomenon. For smaller values of $K^{*}$ i.e. for strong permeability effects the circulating roll in lower half contains two small eddies. The size of lower half is greater than the upper one. With increasing $K^{*},$ the two eddies inside the lower roll merge into a single cell. Moreover, the size of upper roll reduces with increasing $K^{*}$ . The formation of two circulating eddies in the lower roll for small values of $K^{*}$ is never reported in the earlier available literature on the peristaltic flows through porous-saturated curved channel. Figure 28 demonstrates the effects of applied magnetic field on trapping phenomenon. It is observed that the role of applied magnetic field is to increase the size of circulating roll. The streamlines highlighting the effects of dimensionless radius of curvature on trapping has been conveyed through Figure 29. The graphical analysis revealed that both circulating roll become in size with increasing $\gamma .$

Figure 29.

Streamlines in wave frame for (a) $\gamma =2,$ (b) $\gamma =3.5$ and (c) $\gamma \to \infty .$ The other parameters chosen are Ha = 0.5, $K^{*}$  = 0.2, $\beta =0.5,$ $\lambda =0.4,$ and $\theta$  = 1.5.

Concluding remarks

The peristaltic flow of bi-viscosity fluid in a curved porous channel with extended and contracting walls has been studied under the assumptions of long wavelength. The channel with sinusoidally deformed walls is assumed to be porous. Rheological non-Newtonian material characteristics are evaluated by using bi-viscosity fluid model. The heat and mass transfer features are examined in presence of Joule heating effects. The numerical simulations based on finite difference scheme are performed with desirable accuracy. The novel characteristics of pressure, velocity, temperature, concentration and streamlines are graphically underlined. This analysis predicts following interesting observations:

• The flow velocity is observed up to maximum level with variation of $\beta .$

• Pressure rise per wavelength decreases with increasing $\beta$ and $K^{*}.$

• A declining temperature variation has been observed due to impact of $K^{*}.$ However, it increases with increasing $\beta .$

• The mass concentration increases with increasing $R_c$ and $\beta .$

• The size of circulating roll in lower half of the channel boosts up with larger variation of $\beta$ and $K^{*}.$

• The symmetry in velocity and streamlines pattern is observed when $\gamma \to \infty$ .

• The results for Newtonian fluid are achieved for larger values of $\beta .$

Supplemental Material

Appendix – Supplemental material for Heat and mass transfer characteristics in flow of bi-viscosity fluid through a curved channel with contracting and expanding walls: A finite difference approach

Supplemental material, Appendix for Heat and mass transfer characteristics in flow of bi-viscosity fluid through a curved channel with contracting and expanding walls: A finite difference approach by Raheel Ahmed, Nasir Ali, Sami Ullah Khan, Ali Chamkha and Iskander Tlili in Advances in Mechanical Engineering