Numerical Simulation of the Thermally Developed Pulsatile Flow of a Hybrid Nanoﬂuid in a Constricted Channel

: Heat transfer analysis of the pulsatile ﬂow of a hybrid nanoﬂuid through a constricted channel under the impact of a magnetic ﬁeld and thermal radiation is presented. Hybrid nanoﬂuids form a new class of nanoﬂuids, distinguished by the thermal properties and functional utilities for improving the heat transfer rate. The behaviors of a water-based copper nanoﬂuid and water-based copper plus a single-wall carbon nanotube, i.e., ( Cu – SWCNT /water), hybrid nanoﬂuid over each of velocity, wall shear stress, and temperature proﬁles, are visualized graphically. The time-dependent governing equations of the incompressible ﬂuid ﬂow are transformed to the vorticity-stream function formulation and solved numerically using the ﬁnite difference method. The laminar ﬂow simulations are carried out in 2D for simplicity as the ﬂow proﬁles are assumed to vary only in the 2D plane represented by the 2D Cartesian geometry. The streamlines and vorticity contours are also shown to demonstrate the ﬂow behviour along the channel. For comparison of the ﬂow characteristics and heat transfer rate, the impacts of variations in Hartmann number, Strouhal number, Prandtl number, and the thermal radiation parameter are analyzed. The effects of the emerging parameters on the skin friction coefﬁcient and Nusselt number are also examined. The hybrid nanoﬂuid is demonstrated to have better thermal characteristics than the traditional one.


Introduction
Blood flow (pulsatile flow) in the arteries exhibits a periodically echoing time scale that affects the flow-induced mass transfer. This makes the modeling of the physiological processes challenging because the time scales of the characteristics are larger than the pulse period. To address the challenge, researchers look at performing time-averaging over the period to reduce the flow variations within the time period so as to capture a characteristic flow profile. In addition, for discussing the blood flow characteristics in arteries, researchers used constricted channels instead of straight channels because of their geometrical advantages. The phenomena of blood flow through arteries have been associated as an essential factor in hemostasis. For example, many events such as atherogenesis, platelet adhesion, thrombosis, and red blood cell lysis have been connected to various hemodynamic factors.
A nanofluid (NF) is formed by an engineered suspension of nanoparticles (NPs) of a highly conducting substance (such as copper, aluminum, carbon nanotubes, etc.) in a base fluid such as water. The nanofluids are used extensively for their higher thermal conductivity than the base fluids alone. NPs can be categorized into various types by scale, morphology, and physical and chemical properties. As explained by Saidur et al. [1], NFs are in use as heat transfer fluids, biological fluids, pharmaceutical fluids, environmental fluids, etc. Researchers have been analyzing the characteristics of NFs in various applications both experimentally and numerically for even further enhancement in thermal conductivity. Choi [2] was the first one to study the enhancement of thermal conductivity in NFs. Buongiorno [3] simplified the mathematical formulation of the NFs. He described the basic mechanics that lead to the enhancement of the thermal characteristics. This analysis was extended by Kuznetsov and Nield [4] for the boundary layer model of free convective NF flows. Wong and Leon [5] researched NFs as suspensions of NPs in fluids showing substantial improvement in their properties at moderate NP concentrations. Akbar [6] examined the peristaltic flow of incompressible viscous fluid containing metallic NPs in an irregular duct. The combined impact of Lorentz and Kelvin forces in two dimensional Fe 3 O 4 /water NF flow (NFF) was studied by Sheikoleslami et al. [7]. Hayat et al. [8] studied the structures of chemical species on NFF impacted by the rotation of a disc. They reported that the amount of surface shear stress rises with the solid volume fraction of NPs. Paul and Mandal [9] investigated the stretched flow of radiative MHD micropolar NF. Alizadeh et al. [10] studied thermally developed MHD micropolar NFF. They reported that the Nusselt number rises with radiation strength and volume fraction. Dib et al. [11], reported that the presence of NPs enhances heat transfer for their problem of setting the flow between two plates. Cacua et al. [12], demonstrated that NPs improve the thermal characteristics. In practice, NFs might have severe stability issues due to the agglomeration and sedimentation of NPs with time, which limits the applicability of the NFs. To cover the stability issue, Cacua. et al. [12] added surfactants for producing electrostatic or steric repulsion between NPs. Choudhary et al. [13] presented the stability analysis of a (Al 2 O 3 )/water NF. They investigated the stability with the help of zeta potential and visual inspection methods. Lim et al. [14,15] used electroosmotic flow with two or more fluids in various microfluidic applications. Shadloo et al. [16] estimated the convection heat transfer coefficient NFF via circular pipes making use of the Support Vector Machine method. For the forced convection flow of TiO 2 /water NF, Rashidi et al. [17] compared single-phase and two-phase flows in a horizontal tube with a steady wall heat flux as the boundary condition.
A novel class of working NFs, immersed with NPs of two or more solid materials, has come into use during recent years. Such fluids are termed hybrid nanofluids (HNFs). Usually, a material alone does not have all the desirable characteristics needed for a specific purpose; it may have either good thermal or rheological properties. It may be appropriate, however, to have a trade-off between several properties in several practical applications. HNF might serve the desired purpose by combining the properties of its constituent materials. Due to the synergistic impact, HNFs are found to exhibit better thermal conductivity compared with individual NFs [18]. The examples of HNF in use include those formed by the mixing of NPs of two or more of copper (Cu), silver (Ag), copper oxide (CuO), aluminium oxide (Al 2 O 3 ), Titanium oxide (TiO 2 ), single wall carbon nanotubes (SWCNTs), and multi wall carbon nanotube (MWCNTs). Although HNFs have been used in several applications, a lot of research work is still needed for their preparation, characterization, mixing-ratios, and stability [19]. Ahmad and Nadeem [20] discussed the impact of Hall and ion slip with heating for micropolar HNF involving NPs of SWCNT and MWCNT. They also studied the effects of radiation, Darcy-Forchheimer, viscous dissipation, and variable viscosity. Khan et al. [21] numerically studied the impact variable viscosity in inclined MHD Williamson NFF over a nonlinearly stretching sheet. They reported that the velocity profile is declined with the rise of inclination angle, Hartmann number, and variable viscosity. The 3D flow of Cu-AI 2 O 3 /water HNF on an expanding surface was studied by Devi and Anjali [22] using the RK-Fehlberg integration technique. The numerical results demonstrated higher heat transfer rate for the HNF than the copper-based NF. Suleman et al. [23] analyzed the effects of Newtonian heating and thermal radiation in chemically reacting silver/water NFF. Nadeem et al. [24] studied the heat transfer feature with heat generation in SWCNT-MWCNT/water HNF flow. Some examples of relevant studies can be found in [25][26][27][28]. Eshgarf et al. [29] explained that making use of NPs in the conventional fluids to enhance the performance of thermal devices saves energy. Vasua et al. [30] simulated the 2D rheological laminar blood flow via a tapered artery representing mild stenosis. They discussed the effect of various NPs in the blood for pharmacological contexts. Priyadharshini and Ponalagusamy [31] carried out a similar study along with the influence of magnetic field and periodic body acceleration.
In view of the above-mentioned studies, it is notable that no adequate investigations are available in the literature regarding the pulsating MHD HNF flow through a constricted channel subject to thermal radiation. Therefore, to address this research concern, the present study is carried out to analyze the heat transfer of HNF in a constricted channel. The HNF is considered to have NPs of Cu and SWCNT in water. The objective is to examine the cumulative impact of the applied magnetic field and thermal radiation on the wall shear stress (WSS), velocity, and temperature profiles. The effects of various emerging parameters, namely, the Hartman number (M), Strouhal number (St), Prandtl number (Pr), and radiation parameter (Rd) on the HNF are discussed. The effects of Rd, Pr, and the Reynolds number (Re) on the Nusselt number (Nu) and skin friction coefficient (C f ) profiles are also examined. The streamlines and contour graphs are also shown to demonstrate the streamline pattern and flow behviour along the channel. Various comparisons are presented in the present work to observe that which of the NF and HNF possess better thermophysical characteristics. The present study has furnished new insights for the pulsatile flow of a hybrid NF in a constricted channel. The results would provide insight to the relevant researchers and engineers to design concerning products based on the findings. The simulation of the pulsatile flow with nanoparticles is motivated by drug delivery (pharmacology) applications.
The remaining parts of the article are organized in the following manner. The mathematical modeling and its transformation into a solvable form, specifically using the vorticity-stream function formulation, is presented in Section 2. Section 3 discusses the findings and related discussions. Section 4 contains the concluding remarks.

Mathematical Model
We consider the pulsatile flow of an incompressible electrically conductive HNF through a rectangular channel. The HNF, having NPs of Cu and SWCNT in water, is assumed to be a uniform and stable fluid. A uniform-magnetic field B is applied perpendicular to the channel walls. We consider a Cartesian coordinate system ( x, y) such that the flow direction and the direction of B is along the x-axis and y-axis, respectively. The channel walls have a pair of symmetrical constrictions. In the transformed system (x, y), to be discussed later on, the symmetrical constrictions are placed in [−x 0 , x 0 ], as depicted in Figure 1. We assume that the magnetic-Reynolds number for the flow is very low so that the induced electric and magnetic fields are negligible. The lower and upper walls are defined, respectively, as: Here, h 1 and h 2 denote the respective constriction heights. The unsteady incompressible Navier-Stokes equations, in ( x, y) system, are given by: The continuity equation: The momentum equation: The energy equation: ∂ y is the radiative heat flux. Expanding T 4 about T ∞ and ignoring higher-order terms, we get where u and v represent the velocity components along xand y-axes, respectively, hn f in the subscript represents the nanofluid, p the pressure, ρ the density, ν the kinematic viscosity, T the temperature, k the thermal conductivity, C p the specific heat, J ≡ J x , J y , J z the current density, B ≡ (0, B 0 , 0) the magnetic field with the uniform strength B 0 across the flow, σ the electric conductivity, µ the dynamic viscosity, and µ m the magnetic permeability of the medium. If E ≡ E x , E y , E z represents the electric field such that the electric current flows along the normal to the plane of the flow, then E ≡ (0, 0, E z ). Next, from Ohm's law For the steady flow, Maxwell's equation ∇ × E = 0 implies that E z is constant. For the current study, we suppose that E z = 0. Then, Equation (7) gives, J z = σ f uB 0 . Therefore, The following quantities are added for obtaining the non-dimensional form of the governing equations: Here, L is the maximum width of the channel, T is the period of the flow pulsation, U is the characteristic flow velocity, and k * is the mean absorption coefficient.
The effective thermophysical properties of HNF [10,22,32] are presented as follows: The HNF's effective dynamic viscosity µ hn f is given as Here, φ 1 and φ 2 show the volume fractions of NPs Cu and SWCNT, respectively. The HNF's effective density ρ hn f is given as Here, s 1 and s 2 in subscripts correspond to the property of Cu and SWCNT, respectively. Moreover, hn f , n f and f in the subscripts correspond to the property of the hybrid-nanofluid, nanofluid, and base fluid water, respectively. The HNF's heat capacitance ρC p hn f is given as The HNF's effective thermal conductivity k hn f is given as where The thermophysical properties of the base fluid water and the NPs under consideration are shown in Table 1. The properties in Table 1 are given at 24.6 • C. The viscosity of water at 24.6 • C is approximately 0.0000009009 m 2 /s. Thus, according to Equation (10), Using the quantities from Equations (9) and (12) in Equations (2)- (4) and (6) gives ∂u ∂x

The Boundary Conditions
For the steady case of the flow problem under consideration, the boundary condition, obtained by solving the dimensionless form of Equation (8) and performing some manipulations using Equation (10), are given as follows, When M = 0, the u-velocity at the inlet takes the form: The conditions of the outlet boundary are taken as for fully developed flow for all the variables. The flow is presumably sinusoidal for the pulsatile flow: Further, u = 0 and v = 0 (i.e., no-slip conditions) are assumed on the walls. In the non-dimensional form, the temperature condition at the lower wall, θ = 1 and at the upper wall, θ = 0.

The Vorticity-Stream Function Formulation
The transformation from the primitive variables, u and v, to the vorticity-stream functions, ψ and ω, is defined as: The relation for computing the WSS is given by Here,n and ω wall denote the normal vector and vorticity at the wall, respectively. As the vorticity and WSS are orthogonal to each other, the WSS can also be determined by the vorticity expression [33].
We differentiate the momentum Equations (16) and (17) Additionally, the stream function ψ equation (Poisson equation) is given as

Transformation of Coordinates
By incorporating the following coordinate transformation, the constricted portion of the channel is treated as a straight one: In the (ξ, η) system, η = 0 and η = 1 represent the lower and upper walls of the channel, respectively. Moreover, Equations (18) and (24), (25) assume the form as where P = P(ξ, η) = The velocity components u and v in terms of (ξ, η) take the forms In the new coordinate system, the initial conditions for ψ and ω are given by ψ = 0, for all t ≤ 0 ω = 0, for all t ≤ 0 Furthermore, the boundary conditions for the walls in (ξ, η) system for ψ and ω are given as where denotes the pulsating amplitude. The steady and pulsatile flow conditions are given by = 0 and = 1, respectively. For the temperature, the boundary conditions after the transformation are given as θ = 1 at η = 0 and θ = 0 at η = 1. The other concerning nondimensional physical quantities include the skin friction coefficient and Nusselt number, defined by where τ w and q w are defined as After using dimensionless variables from (9) and the coordinate transformation from (26), we get

Results and Discussion
The problem Equations (27)-(29) subject to boundary conditions, including (32)- (34), are solved numerically using FDM, as in [34][35][36][37][38][39]. The domain {(ξ, η)|ξ ∈ [−x 1 , x 1 ] and η ∈ [0, 1]} is discretized by choosing the grid points ξ i , η j with the step-sizes ∆ξ and ∆η in ξ-direction and η-direction, respectively. For time integration of the solution, a fixed step size ∆t is used. The solution at time level l is used to obtain the solution is at time level l + 1, for l = 0, 1, 2, · · · . Firstly, Equation (28) is solved for ψ = ψ(ξ, η). The linear system obtained through the finite difference discretizations is solved using the Tri-Diagonal Matrix Algorithm (TDMA). Secondly, Equations (27) and (29) are solved using the ADI (Alternating Direction Implicit) method for ω = ω(ξ, η) and θ = θ(ξ, η), respectively. The linear systems obtained through the finite difference discretizations are solved using the TDMA at each of the two half-steps of the ADI method. A detailed description of the algorithm can be found in [38]. The values of the step sizes are taken as follows to ensure the consistency of the numerical scheme: ∆t = 0.00005 for the pulsatile flow with ∆ξ = 0.05, ∆η = 0.02, and Re = 700. The computations for the present study are performed in a sequential fashion. The results can be found by parallel computing for time-efficient solutions [40].
In the current study, the computational domain, considered as −10 ≤ ξ ≤ 10 and 0 ≤ η ≤ 1, is discretized by a grid of 400 × 50 elements. The length of the constriction is considered as 2x 0 = 4. The constrictions on the lower and upper walls are assumed to have heights h 1 = h 2 = 0.35; hence, at the throat of the constriction, the channel width remains 30% of the maximum channel width.
The pulsatile motion is modeled by adding in the inflow boundary condition the sinusoidal time-dependent function sin(2πt). For a comparison of the behaviors of Cu-SWCNT/water HNF with Cu/water NF, the effects of the physical parameters in-cluding M, St, Pr, and Rd are studied on the non-dimensional streamwise velocity (u) and temperature (θ) profiles. We performed simulations for a long enough time that, in most cases, the results are displayed at selected time instants and axial locations, especially x = 0 (throat of the constriction) and x = 2 (in the lee of the constriction) where the fluid enters in the low-pressure zone from the high-pressure zone. For a pulse cycle, 0 < t < 0.25 is the acceleration phase and 0.25 < t < 0.75 is the deceleration phase.
For validation, the present results for the pulsatile flow of the base fluid (i.e., with φ 1 = 0, φ 2 = 0) are compared with those obtained by Bandyopadhyay and Layek [35] in the case of Newtonian fluid without the heat effect. Figure 2 shows a good agreement of the present results, specifically the wall shear stresses, with [35] for M = 5 and 10 at t = 0.25.  Figure 3. It is seen that during a complete period, the wall shear stress varies considerably, especially in the vicinity and downstream of the constriction. However, no major changes are observed upstream of the constriction on the WSS. Over the duration 0 < t < 0.25, the peak value of WSS decreases and reaches its highest value at t = 0.25 at peak flow time. The flow separation area rises during this period of the time cycle. After that, during 0.25 < t < 0.75, the flow starts to decelerate, and the WSS also decreases, but the area of the flow separation tends to increase. The WSS is oscillating, and its sign changes at zero net flow rate, i.e., at t = 0.75. Figure 4 presents a comparison among the pure water, Cu/water NF, and Cu-SWCNT /water HNF by fixing the parameters as φ = 0.02, M = 5, St = 0.02, Pr = 1, and Rd = 0.2. The velocity profile for pure water is observed to be the highest, followed by the NF and then the HNF. For the HNF, the particle concentration in the fluid is even higher than that of the base fluid and the NF, which causes enhancement of the density and dynamic viscosity. The flow velocity experiences reduction accordingly. The peak value of u tends to rise, and the curves appear to be parabolic for 0 ≤ t ≤ 0.5. At t = 0.75, u is dropped substantially, and backflow occurs near the walls. Figure 5    A comparison among the pure water, Cu/water NF, and Cu-SWCNT/water HNF for the temperature profiles computed with φ = 0.02, M = 5, St = 0.02, Pr = 1, and Rd = 0.2. In Figure 6, the profiles are shown at x = 0 for a complete pulse cycle. However, in Figure 7, the profiles are shown at four different x-locations at t = 0.25. It is noticed that the HNF hits higher temperatures compared to the NF. A rapid rise in temperature is due to the hybrid nature of the nanofluid since it increases the thermal conductivity. Figures 8 and 9 illustrate the behavior of u and θ profiles, respectively, for the variation of M to make a comparison between the NF and the HNF. The other parameters are fixed as φ = 0.02, St = 0.02, Pr = 1, and Rd = 0.2. It is observed that the velocity is greater for the NF than the HNF. The velocity increases as the magnetic field strength is increased. That is, the viscous boundary layer thickness is decreased for higher values of M. This is due to the production of a resistive Lorentz force, whose strength increases with increasing values of M. The HNF temperature profile is slightly higher than that of the NF. The thermal boundary layer thickness decreases, resulting in a decline of the temperature profile with enhancement in the magnetic field strength. The appreciation in M maximizes the velocity gradient and minimizes the temperature gradient at the boundary. It is observed that the velocity is greater for the NF than the HNF. The velocity decreases with increment in St, i.e., the viscous boundary layer thickness increases with St. The HNF temperature profile is slightly higher than that of the NF. The thermal boundary layer thickness decreases, resulting in a decrease of the temperature profile with increasing values of St. The velocity exhibits a parabolic profile at x = 0. The profiles are not parabolic as some backflow in the vicinity of the walls is observed at x = 2. The backflow reduces with an increase in St. Figure 12 illustrates the behavior of the temperature profiles for the variation of Pr to compare the NF and HNF. The other parameters are fixed as: φ = 0.02, M = 5, St = 0.02, and Rd = 0.2. The HNF temperature profile is slightly higher than that of the NF. The thermal boundary layer thickness decreases as the Pr increases, resulting in a decline in the temperature profile at both locations. The rising values of Pr cause a reduction in thermal diffusivity. Consequently, the temperature of the fluid is lowered.  Figure 13 illustrates the behavior of the temperature profiles for the variation of Rd to compare the NF and HNF. The other parameters are fixed as φ = 0.02, M = 5, St = 0.02, and Pr = 1. Appreciation in Rd enhances the temperature of the fluids. It is noticed that the temperature profile of the HNF is slightly higher than that of the NF. The thermal boundary layer thickness increases, resulting in the rise of temperature as Rd is increased. Physically it is due to the fact that the mean absorption coefficient decreases, which boosts Rd, thus the radiative heat transfer rate is increased. Figures 14 and 15 illustrate the behavior of the velocity and temperature profiles, respectively, for the variation of solid volume fraction φ to compare the NF and HNF. The other parameters are fixed as M = 5, St = 0.02 Pr = 1, and Rd = 0.2. The velocity is greater for the NF than that of the HNF. We further observe that the velocity decreases for enhancing the solid volume fraction φ, i.e., the viscous boundary layer thickness increases with increasing values of φ. The temperature distribution and momentum boundary layer thickness show a rising trend with the solid volume fraction. It is also clear from the figure that the HNF exhibits a higher temperature profile than the NF. The velocity exhibits a parabolic profile at x = 0. The profiles are not parabolic as some backflow in the vicinity of the walls is observed at x = 2.    The results are shown in Figure 16, and the streamlines distribution plots are shown Figure  17 for t = 0.25. As expected, the magnitudes of the peak values of the profiles are lower for the lower Re. The profile wave and the eddies downstream of the constriction are weaker for the lower Reynold number.
The effects of Rd, Pr, and Re on the Nusselt number and skin friction profiles are shown in Figures 18 and 19, respectively. It is observed that the Nusselt number increases for higher values of Pr and Re, but decreases for increasing Rd. Additionally, the skin friction coefficient increases for an increase in the values of Re. However, it remains unchanged for Pr and Re.          Figure 20 for comparison among pure water, Cu/water NF and Cu-SWCNT/water HNF. The formation of vertical eddies in the vicinity of the walls can be observed. Over time, the eddies grow for the HNF and slowly occupy a major part of the channel downstream of the constriction. The inclusion of multiple types of nanoparticles adds more energy. This results in higher increments in the temperature as well as the thickness of the thermal boundary layer. A rapid rise in the temperature is due to the HNF since it increases the thermal conductivity.

Concluding Remarks
In this research, the comparison between the behavior of a traditional NF and an emerging HNF in the presence of viscous incompressible MHD pulsatile fluid flow over a rectangular channel is numerically investigated. The comparison between the behavior of the pure water, Cu/water NF, and Cu-SWCNT/water HNF over the velocity, temperature, and WSS distributions are visualized graphically. The impacts of each of M, St, Pr, and Rd on the flow profiles are studied. The summary of the impacts is as follows: • The effects of the pure water, Cu/water, and Cu-SWCNT/water on the WSS at both walls are computed. The flow separation area rises during the duration 0 < t < 0.25. After that, the flow starts to decelerate, during 0.25 < t < 0.75, and the WSS also decreases, but the area of the flow separation tends to increase. • Density and dynamic viscosity are the factors of fluid velocity reduction, which is enhanced by hybridity as more immense particles are introduced, so the velocity reduces when comparing the pure water, the NF, and the HNF. • A rapid rise in temperature is observed due to the HNF since it increases the thermal conductivity when comparing the pure water, the NF, and the HNF. Hence, it is interpreted that the HNF hits higher temperatures compared to the NF.

•
The increasing values of M reduce the overall density, which raised the temperature of the HNF within the boundary layer. Finally, the appreciation in M maximizes the velocity gradient and minimizes the temperature gradient at the boundary.

•
The HNF temperature profile is slightly higher than that of the NF in the case of varying St and Pr. The thermal boundary layer thickness decreases, resulting in a decrease of the temperature profile with increasing values of St and Pr. The higher value of Pr number causes a reduction in the thermal diffusivity.

•
The thermal boundary layer thickness increases, resulting in the rise of temperature as Rd is increased. Temperature profile and thermal boundary layer thickness are enhanced by increasing the values of Rd. • Velocity decreases for enhancing the solid volume fraction φ, i.e., the viscous boundary layer thickness increases with increasing values of φ. Additionally, the HNF exhibits a higher temperature profile than the NF.

•
The HNF is observed to be a better thermal conductor as compared with the traditional NF.

•
The Nusselt number increases for higher values of Pr and Re, but decreases for increasing Rd.

•
The skin friction coefficient increases for an increase in the values of Re; however, it remains unchanged for Pr and Re.