Influence of Stress Jump Condition at the Interface Region of a Two-Layer Nanofluid Flow in a Microchannel with EDL Effects

Abstract

The influence of stress jump conditions on a steady, fully developed two-layer magnetohydrodynamic electro-osmotic nanofluid in the microchannel, is investigated numerically. A nanofluid is partially filled into the microchannel, while a porous medium, saturated with nanofluid, is immersed into the other half of the microchannel. The Brinkmann-extended Darcy equation is used to effectively explain the nanofluid flow in the porous region. In both regions, electric double layers are examined, whereas at the interface, Ochoa-Tapia and Whitaker’s stress jump condition is considered. The non-dimensional velocity, temperature, and volume fraction of the nanoparticle profiles are examined, by varying physical parameters. Additionally, the Darcy number, as well as the coefficient in the stress jump condition, are investigated for their profound effect on skin friction and Nusselt number. It is concluded that, taking into account the change in shear stress at the interface has a significant impact on fluid flow problems.

1. Introduction

Two-layer flow in a microchannel is essential in practical applications like crude-oil extraction, thermal insulation, solidification of castings, and several other geophysical applications. Another example is the design of micro-electromechanical systems (MEMS). In addition, fluid flow properties depict unusual behaviors in a microchannel compared to a macro-scale channel. Consequently, it is of significant importance to scientifically study the two-layer microchannel flow, particularly taking into account the possible effect of EDL. Due to this reason, many research studies have been conducted on flows through a microchannel, considering the electric double layer effects for Newtonian fluids [1,2,3,4], and non-Newtonian fluids [5,6,7,8,9,10]. However, most of the works mentioned above are connected to single-layer flow. The flow attributes of immiscible liquid are noticeable in the biochemical and biological investigation processes [11]. A laminar fluid interface is rendered when two or more immiscible liquids stream in microfluidic devices. In most cases, the influences of the fluid interface are noteworthy and cannot be neglected in the investigation of biological sample separation. Some research studies that have investigated this correlation include the work Gao et al. [12], who obtained theoretical and experimental results to investigate the two-fluid electro-osmotic flow in microchannels, but the Maxwell stress balance condition at the interface was not taken in account. Later, Gao et al. [13] modified the interface condition, by including the shear stress balances, that result in a jump at the interface resulting from the specific surface charge density. Some more exciting work [14,15,16] includes the investigation of two-layer microchannel flow along with the electro-osmotic effect, and using the shear stress balance interface condition. Recently, Niazi and Xu [17] used nanofluids to assess the electro-osmotic effect in two-layer microchannel flow. They used Buongiorno’s model [18] to construct a mathematical model, and obtained the analytic solutions for their problem. Mainly, they concluded that the flow behavior was altered dramatically in the presence of Brownian diffusion, thermophoresis diffusion, and viscosity. H. Tahir et al. [19] used the optimal homotopy approach, to analyze the performance of a hybridized two-phase ferromagnetic nanofluid of ferrite nanoparticles, and their effects on heat transmission in the flow of the hybrid nanofluid. Based on their investigation, it can be concluded that the thermophysical characteristics and Curie temperature with two or more ferrites suspended in two or more base fluids, can be enhanced. In-depth analysis by Hammad et al. [20], covers the numerous uses of nanofluids, as well as the implications of variables such as nanoparticle type and size, which may open up new prospects for commercial applications.

Porous media are also critical for exploring the applications described above. For example, thickening alloys do not have a eutectic composition, resulting in the separation of the frozen and liquid portions of the casting. In this instance, the partially frozen areas can be thought of as a porous medium with varying permeability. While porous media have been used for a wide range of commercial and geological purposes, there are opportunities to investigate alternative uses, particularly for energy systems, such as compact heat exchangers, heat pipes, electronic cooling, and solar collectors, by exploiting porous media. For certain applications, it is not necessary to entirely fill the system with the porous medium; partial filling is adequate. In comparison to a system that is totally filled with porous media, partial filling reduces the pressure drop. In addition, partial filling prevents contact between the porous material and the surface, reducing heat loss from the porous material to the surface. Such a criterion is necessary in a system where the primary objective is to improve the thermal coupling between the porous medium and fluid flow, and reduce the system’s high thermal coupling with the surrounding environment. For instance, the objective of Mohamad’s [21] solar air heater, was to increase the rate of heat transfer from the porous medium, which is heated by solar radiation, to air, while minimizing heat loss to the ambient environment. In addition, partial filling helps to decrease the pressure drop. A partial filling of a channel with porous media drives the flow to exit from the core area to the outer region, depending on the permeability of the medium. This decreases the thickness of the boundary layer and therefore increases the rate of heat transfer. The porous medium also alters the effective thermal conductivity and heat capacity of the flow, and the solid matrix increases the rate of radiative heat transfer in a gas-based system. Hence, increases in heat transfer occur through three mechanisms: flow redistribution, thermal conductivity adjustment, and medium radiative property modification. Beavers and Joseph [22] pioneered this type of study, by modeling flow in a porous material using Darcy’s law. The effects of the interfacial layer on fluid mechanics and heat transmission are discussed in further detail in [23,24]. These articles investigate non-Darcian effects in flow in porous media, via the Brinkman–Forchheimer-extended Darcy equation. In [24], the authors presented a precise solution for the flow field at the interface. The fluid layer is located among a semi-infinite porous object and an impermeable outer border, in their proposed model. Nield [25] demonstrated that velocity shear is continuous across the porous part of the contact. This is not always the case for solid sections, as the averaged velocity shears do not always coincide. Then, Kuznetsov [26], and Ochoa-Tapia and Whitaker [27,28], developed the strategy for comparing the Brinkman-extended Darcy law to Stokes’ equations, that need a discontinuity in the stress but retain continuity in the fluid flow. They determined that solving at the interface utilizing the Ochoa-Tapia and Whitaker conditions, resolves the over-determination problem demonstrated in Nield [25].

We intend to investigate fluid flow in a microchannel half filled with porous media, in light of the practical implications of two-layer fluid flow in a microchannel. In the region with a porous layer, the Brinkmann-extended Darcy’s law is used to mathematically predict fluid flow, whereas Buongiorno’s model is used in the other zone. For this topic, we used the stress jump boundary condition at the interface, which had been overlooked in prior research, as well as the impacts of the electric double layer (EDL) and magnetic field. Utilizing the interface stress jump condition, it is possible to correct for the overestimation of the physical parameters involved in the problem. The Darcy number and stress jump condition variations are critical in analyzing heat and mass transport in this two-layer fluid flow problem.

2. Problem Formulation

We analyze the flow of an electro-osmotic fluid within a microchannel divided into two distinct regions (I and II). The elongated rectangular microchannel is horizontally positioned, with a width W, that is adequately greater than its height H ($W/H>4$; see Dauenhauer and Majdalani [29]). The length of the microchannel, L, is believed to be sufficient to prevent the apertures at the end from having an effect. $H_1+H_2=H$ are the height of the lower and upper layers, respectively. The interface among immiscible fluids is planar, based on the aforementioned assumptions. The parallel flow proposition can also be used to reduce the dimensions of the problem to two $\left(2D\right)$. In Figure 1, the Cartesian coordinate system $(x,y,z)$ is used, with x along the streamwise direction, y is parallel to the surfaces and normal to x, and z is perpendicular to the plates, parallel to each other. The lower and top walls have zeta potentials, temperature and nanoparticle volumes are represented as ${\overline{\zeta }}_1,T_w,C_w$ and ${\overline{\zeta }}_2,T_w,C_w$. Region I receives nanofluid containing Al₂O₃ nanoparticles, whereas Region II has porous media saturated with TiO₂.

Table 1. Numeric values for the physical attributes of the fluid and nanoparticles.

Physical Characteristic	${\mathbf{H}}_2\mathbf{O}$	${\mathbf{Al}}_2{\mathbf{O}}_3$	${\mathbf{TiO}}_2$
$c_p(\mathrm{J}·{\mathrm{kg}}^{-1}·{\mathrm{K}}^{-1})$	4179.0	765.0	686.2
$\rho (\mathrm{kg}·{\mathrm{m}}^{-3})$	997.1	3970.0	4250.0
$k(\mathrm{W}·{\mathrm{m}}^{-1}·{\mathrm{K}}^{-1})$	0.6130	40.0	8.9538
$\alpha \times {10}^{-7}({\mathrm{m}}^2·{\mathrm{s}}^{-1})$	1.47	131.70	30.70
$\beta \times {10}^{-5}\left({\mathrm{K}}^{-1}\right)$	21.00	0.85	0.90

lists the physical parameters of the fluid and nanoparticles. The Buongiorno model is used to simulate nanofluid flow in Region I. The Brinkmann-extended Darcy law is employed to illustrate the flow of nanofluids in a porous layer region. The steady-state laminar flow is considered to be one-dimensional, owing to the significant presence of an electric field. due to the presence of an electric double layer (EDL) and the applied pressure.

The governing equations are modeled after the Navier-Stokes equations, with the driving force deriving from the electric and magnetic field, along with a pressure gradient. The mathematical models representing the physical phenomena in both the regions are as follows:

• Region I:

$$\begin{array}{cc}\hfill {\nabla }^2{\overline{\psi }}_1=-& \frac{{\overline{\rho }}_{e1}}{{\epsilon }_0{\epsilon }_{R_1}},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \nabla ·{\mathbf{V}}_{\mathbf{1}}=& 0,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\left({\rho }_1\right)}_f({\mathbf{V}}_{\mathbf{1}}·\nabla ){\mathbf{V}}_{\mathbf{1}}=& -\nabla \overline{P}+{\mu }_1{\nabla }^2{\mathbf{V}}_{\mathbf{1}}+{\mathbf{F}}_{\mathbf{1}},\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill {\left({\rho }_1{c}_{p1}\right)}_f({\mathbf{V}}_{\mathbf{1}}·\nabla ){\overline{T}}_1=& k_{nf1}{\nabla }^2{\overline{T}}_1+{\left({\rho }_1{c}_{p1}\right)}_s[D_{B_1}\nabla {\overline{T}}_1·\nabla {\overline{C}}_1\hfill \\ \hfill & +\frac{D_{T_1}}{T_0}\nabla {\overline{T}}_1·\nabla {\overline{T}}_1]+{\mu }_1\overline{{\Phi }_1},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill ({\mathbf{V}}_{\mathbf{1}}·\nabla ){\overline{C}}_1=& D_{B_1}{\nabla }^2{\overline{C}}_1+\frac{D_{T_1}}{T_0}{\nabla }^2{\overline{T}}_1.\hfill \end{array}$$

(5)

Region II:

$$\begin{array}{cc}\hfill {\nabla }^2{\overline{\psi }}_2=-& \frac{{\overline{\rho }}_{e2}}{{\epsilon }_0{\epsilon }_{R_2}},\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \nabla ·{\mathbf{V}}_{\mathbf{2}}=& 0,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{{\left({\rho }_2\right)}_f}{{ϵ}^2}({\mathbf{V}}_{\mathbf{2}}·\nabla ){\mathbf{V}}_{\mathbf{2}}=& -\nabla \overline{P}+\frac{{\mu }_2}{ϵ}{\nabla }^2{\mathbf{V}}_{\mathbf{2}}-\frac{{\mu }_2}{\kappa }{\mathbf{V}}_{\mathbf{2}}+{\mathbf{F}}_{\mathbf{2}},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\left({\rho }_2{c}_{p2}\right)}_f{\mathbf{V}}_{\mathbf{2}}·\nabla {\overline{T}}_2=& k_{nf2}(\nabla ·\nabla {\overline{T}}_2)+ϵ{\left({\rho }_2{c}_{p2}\right)}_s[D_{B_2}\nabla {\overline{T}}_2·\nabla {\overline{C}}_2\hfill \\ \hfill & +\frac{D_{T_2}}{T_0}\nabla {\overline{T}}_2·\nabla {\overline{T}}_2]+\overline{{\Phi }_2},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{1}{ϵ}({\mathbf{V}}_{\mathbf{2}}·\nabla ){\overline{C}}_2=& D_{B_2}{\nabla }^2{\overline{C}}_2+\frac{D_{T_2}}{T_0}{\nabla }^2{\overline{T}}_2.\hfill \end{array}$$

(10)

Here, ${\overline{\psi }}_1$ and ${\overline{\psi }}_2$ represent the dimensional electrostatic potentials in the two regions, and ${\mathrm{\Phi }}_1$ and ${\mathrm{\Phi }}_2$ are the viscous dissipation factor in two regions. The general forms of ${\mathrm{\Phi }}_1$ and ${\mathrm{\Phi }}_2$ are as follows:

$$\begin{array}{cc}\hfill {\overline{\mathrm{\Phi }}}_1=& 2\left({\left(\frac{\partial {\overline{u}}_1}{\partial \overline{x}}\right)}^2+{\left(\frac{\partial {\overline{v}}_1}{\partial \overline{y}}\right)}^2+{\left(\frac{\partial {\overline{w}}_1}{\partial \overline{z}}\right)}^2\right)+{\left(\frac{\partial {\overline{v}}_1}{\partial \overline{x}}+\frac{\partial {\overline{u}}_1}{\partial \overline{y}}\right)}^2+{\left(\frac{\partial {\overline{w}}_1}{\partial \overline{y}}+\frac{\partial {\overline{v}}_1}{\partial \overline{z}}\right)}^2\hfill \\ \hfill & +{\left(\frac{\partial {\overline{u}}_1}{\partial \overline{z}}+\frac{\partial {\overline{w}}_1}{\partial \overline{x}}\right)}^2-\frac{2}{3}{\left(\frac{\partial {\overline{u}}_1}{\partial \overline{x}}+\frac{\partial {\overline{v}}_1}{\partial \overline{y}}+\frac{\partial {\overline{w}}_1}{\partial \overline{z}}\right)}^2,\hfill \end{array}$$

(11)

$${\overline{\mathrm{\Phi }}}_2=\frac{{\mu }_2}{\kappa }\mathbf{V}·\mathbf{V}+{\mu }_{eff}\mathbf{V}·{\nabla }^2\mathbf{V}$$

(12)

2.1. Problem Statement and Assumptions

• The direction of the flow is assumed to be along the x-axis.
• The flow velocity in the $\overline{z}$-direction is negligible, since the length of microchannel L is much larger than its height H. Hence ${\overline{w}}_i\approx 0$,
• The velocity component in the $\overline{y}$-direction is considered to be zero, i.e., ${\overline{v}}_i=0$,
• The flow is assumed to uni-directional along the $\overline{x}$-axis but its properties changes with respect to the $\overline{z}$-axis, hence ${\mathbf{V}}_{\mathbf{i}}=({\overline{u}}_i\left(\overline{z}\right),0,0)$,
• The body force, ${\mathbf{F}}_{\mathbf{i}}={\mathrm{\ae }}_{\mathbf{ei}}\mathbf{E}+J_i\times \mathbf{B}$, represents the sum of electro-osmosis and the electromagnetic forces, where $\mathbf{E}=(E_x,E_y,0)$ is the electric field, $\mathbf{B}=(0,0,B_0)$ is the applied magnetic field, and $J_i={\sigma }_i(\mathbf{E}+{\mathbf{V}}_i\times B)$ is the current density of the ion.
• The inertial effects in the porous region of the microchannel (Region II) are negligible.
• Region I of the channel is filled with nanofluid, while the channel’s Region II is filled with the porous medium saturated with nanofluid, having uniform permeability only.
• Proceeding from the analysis presented in [26], the stress jump condition is utilized at the interface. Simultaneously, the electric potential, temperature, nanoparticle concentration, and flux at the interface are presumed to be continuous. Finally, the no-slip condition is applied to the velocity boundaries, while the temperature and nanoparticle concentration are assumed to have a constant distribution on the boundaries.

In light of the above assumptions, Equations (1)–(10) now take the form,

• Region I: ($-H_1\le \overline{z}\le 0$)

$$\begin{array}{cc}\hfill \frac{{\partial }^2{\overline{\psi }}_1}{\partial {\overline{z}}^2}=& -\frac{{\overline{\rho }}_{e1}\left(\overline{z}\right)}{{\epsilon }_0{\epsilon }_{R_1}}\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\mu }_1\frac{{\partial }^2{\overline{u}}_1}{\partial {\overline{z}}^2}=& \frac{\partial \overline{P}}{\partial \overline{x}}-E_x{\overline{\rho }}_{e1}-{\sigma }_1{B}_0(E_y-B_0{\overline{u}}_1)\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\overline{u}}_1\frac{\partial {\overline{T}}_1}{\partial \overline{x}}+{\overline{w}}_1\frac{\partial {\overline{T}}_1}{\partial \overline{z}}=& {\alpha }_1\frac{{\partial }^2{\overline{T}}_1}{\partial {\overline{z}}^2}+{\tau }_1[D_{B_1}\frac{\partial {\overline{T}}_1}{\partial \overline{z}}\frac{\partial {\overline{C}}_1}{\partial \overline{z}}\hfill \\ \hfill & +\frac{D_{T_1}}{T_0}{\left(\frac{\partial {\overline{T}}_1}{\partial \overline{z}}\right)}^2]+\frac{{\mu }_1}{{\left({\rho }_1{c}_{p1}\right)}_f}{\left(\frac{\partial {\overline{u}}_1}{\partial \overline{z}}\right)}^2\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\overline{u}}_1\frac{\partial {\overline{C}}_1}{\partial \overline{x}}+{\overline{w}}_1\frac{\partial {\overline{C}}_1}{\partial \overline{z}}=& D_{B_1}\frac{{\partial }^2{\overline{C}}_1}{\partial {\overline{z}}^2}+\frac{D_{T_1}}{T_0}\frac{{\partial }^2{\overline{T}}_1}{\partial {\overline{z}}^2}.\hfill \end{array}$$

(16)

Region II: ($0\le \overline{z}\le H_2$)

$$\begin{array}{cc}\hfill \frac{{\partial }^2{\overline{\psi }}_2}{\partial {\overline{z}}^2}=& -\frac{{\overline{\rho }}_{e2}\left(\overline{z}\right)}{{\epsilon }_0{\epsilon }_{R_2}},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\mu }_{eff}\frac{{\partial }^2{\overline{u}}_2}{\partial {\overline{z}}^2}=& \frac{\partial \overline{P}}{\partial \overline{x}}-E_x{\overline{\rho }}_{e2}-{\sigma }_2{B}_0(E_y-B_0{\overline{u}}_2)-\frac{{\mu }_2}{\kappa }{\overline{u}}_2,\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill {\mu }_{eff}\frac{\partial {\overline{T}}_2}{\partial \overline{x}}+{\overline{w}}_2\frac{\partial {\overline{T}}_2}{\partial \overline{z}}=& {\alpha }_2\frac{{\partial }^2{\overline{T}}_2}{\partial {\overline{z}}^2}+{\tau }_2[D_{B_2}\frac{\partial {\overline{T}}_2}{\partial \overline{z}}\frac{\partial {\overline{C}}_2}{\partial \overline{z}}\hfill \\ \hfill & +\frac{D_{T_2}}{T_0}{\left(\frac{\partial {\overline{T}}_2}{\partial \overline{z}}\right)}^2]+\frac{1}{{\left({\rho }_2{c}_{p2}\right)}_f}\frac{{\mu }_2}{\kappa }{\overline{u}}_2^2,\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\overline{u}}_2\frac{\partial {\overline{C}}_2}{\partial \overline{x}}+{\overline{w}}_2\frac{\partial {\overline{C}}_2}{\partial \overline{z}}=& D_{B_2}\frac{{\partial }^2{\overline{C}}_2}{\partial {\overline{z}}^2}+\frac{D_{T_2}}{T_0}\frac{{\partial }^2{\overline{T}}_2}{\partial {\overline{z}}^2}.\hfill \end{array}$$

(20)

Here, ${\alpha }_i=k_{n{f}_i}/{\left({\rho }_i{c}_{p_i}\right)}_f$ is the thermal diffusivity, with $i=1,2$ representing Region I and Region II, ${\tau }_1={\left({\rho }_1{c}_{p_1}\right)}_s/{\left({\rho }_1{c}_{p_1}\right)}_f$, ${\tau }_2=ϵ{\left({\rho }_2{c}_{p_2}\right)}_s/{\left({\rho }_2{c}_{p_2}\right)}_f$ is the heat capacity ratio between the two regions of the microchannel, and ${\mu }_{eff}={\mu }_2/ϵ$, where $ϵ$ is the porosity. The boundary conditions for the above stated governing equations in the two regions are as follows:

• when $z=-H_1$:

$${\overline{\psi }}_1={\overline{\zeta }}_1,{\overline{u}}_1=0,{\overline{T}}_1=T_w,{\overline{C}}_1=C_w,$$

(21)

when $z=0$:

$$\begin{array}{ccc}\hfill {\overline{\psi }}_1={\overline{\psi }}_2,{\overline{u}}_1& =& {\overline{u}}_2,{\overline{T}}_1={\overline{T}}_2,{\overline{C}}_1={\overline{C}}_2,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\epsilon }_1\frac{\partial {\psi }_1}{\partial \overline{z}}={\epsilon }_2\frac{\partial {\psi }_2}{\partial \overline{z}},{\mu }_{eff}\frac{\partial {\overline{u}}_2}{\partial \overline{z}}-{\mu }_1\frac{\partial {\overline{u}}_1}{\partial \overline{z}}& =& \frac{\beta {\mu }_1}{\sqrt{\kappa }}\overline{u_2},k_{n{f}_1}\frac{\partial {\overline{T}}_1}{\partial \overline{z}}=k_{n{f}_2}\frac{\partial {\overline{T}}_2}{\partial \overline{z}},\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \frac{D_{T_1}}{T_0}\frac{\partial {\overline{T}}_1}{\partial \overline{z}}+D_{B_1}\frac{\partial {\overline{C}}_1}{\partial \overline{z}}& =& \frac{D_{T_2}}{T_0}\frac{\partial {\overline{T}}_2}{\partial \overline{z}}+D_{B_2}\frac{\partial {\overline{C}}_2}{\partial \overline{z}},\hfill \end{array}$$

(24)

when $z=H_2$:

$${\overline{\psi }}_2={\overline{\zeta }}_2,{\overline{u}}_2=0,{\overline{T}}_2=T_w,{\overline{C}}_2=C_w.$$

(25)

where $\beta$ is the adjustable stress jump coefficient and $\kappa$ is the permeability of the porous medium. The Poisson–Boltzmann equation [29], simplifies the relationship between the electrostatic potential ${\overline{\psi }}_i$, near the surface and the cumulative number of electrical charges for each unit of volume ${\rho }_{ei}$, at any point in the fluid.

$$\frac{{\partial }^2{\overline{\psi }}_i}{\partial {\overline{z}}^2}=\frac{2n_0\widehat{z}e}{{\epsilon }_0{\epsilon }_{R_i}}sinh\left(\frac{\widehat{z}e\overline{{\psi }_i}\left(\overline{z}\right)}{k_B\widehat{T}}\right),$$

(26)

When the electrical potential is sufficiently small in comparison to the thermal energy of the ions, the Debye–Huckel linear approximation holds true, i.e., $|k_B\widehat{T}|\ll |\widehat{z}e{\overline{\psi }}_i\left(\overline{z}\right)|$, Equation (26) is reduced to

$$\frac{{\partial }^2{\overline{\psi }}_i}{\partial {\overline{z}}^2}-\left(\frac{2n_0{\widehat{z}}^2{e}^2}{{\epsilon }_0{\epsilon }_{R_i}{k}_B\widehat{T}}\right){\overline{\psi }}_i=0,i=1,2.$$

(27)

2.2. Problem Non-Dimensionalization

To transform the governing equations to dimensionless forms, we introduce the following similarity transformations:

$$\begin{array}{ccc}& & \eta =\frac{\overline{z}}{H},{\mathrm{\Psi }}_i=\frac{\widehat{z}e\overline{{\psi }_i}}{k_B\widehat{T}},u_i=\frac{{\overline{u}}_i}{U_{ai}},{\theta }_i=\frac{{\overline{T}}_i-T_0}{T_w-T_0},{\varphi }_i=\frac{{\overline{C}}_i-C_0}{C_w-C_0},\hfill \end{array}$$

(28)

Substituting the non-dimensional variables defined in Equation (28), the fluid flow region is changed to $[-h_1,h_2]$, with $h_1=H_1/H,h_2=H_2/H$, and the governing equations now take the form:

• Region I $(-h_1\le \eta \le 0)$:

$$\begin{array}{ccc}& & {\mathrm{\Psi }}_1^{″}-k_1^2{\mathrm{\Psi }}_1=0,\hfill \end{array}$$

(29)

$$\begin{array}{ccc}& & u_1^{″}-H{a}_1^2{u}_1+S_{e_1}H{a}_1+{\mathrm{\Gamma }}_1+k_1^2{\mathrm{\Psi }}_1=0,\hfill \end{array}$$

(30)

$$\begin{array}{ccc}& & {\theta }_1^{″}+N_{B1}{\theta }_1^{\prime }{\varphi }_1^{\prime }+N_{T1}{\left({\theta }_1^{\prime }\right)}^2+B{r}_1{\left(u_1^{\prime }\right)}^2=0,\hfill \end{array}$$

(31)

$$\begin{array}{ccc}& & {\varphi }_1^{″}+\frac{N_{T1}}{N_{B1}}{\theta }_1^{″}=0.\hfill \end{array}$$

(32)

Region II $(0\le \eta \le h_2)$:

$$\begin{array}{ccc}& & {\mathrm{\Psi }}_2^{″}-k_2^2{\mathrm{\Psi }}_2=0,\hfill \end{array}$$

(33)

$$\begin{array}{ccc}& & u_2^{″}+\frac{1}{{\gamma }^2}(-H{a}_2^2{u}_2+S_{e_2}H{a}_2+{\mathrm{\Gamma }}_2+k_2^2{\mathrm{\Psi }}_2-\frac{1}{Da}{u}_2)=0,\hfill \end{array}$$

(34)

$$\begin{array}{ccc}& & {\theta }_2^{″}+N_{B2}{\theta }_2^{\prime }{\varphi }_2^{\prime }+N_{T2}{\left({\theta }_2^{\prime }\right)}^2+\frac{B{r}_2}{Da}{u}_2^2=0,\hfill \end{array}$$

(35)

$$\begin{array}{ccc}& & {\varphi }_2^{″}+\frac{N_{T2}}{N_{B2}}{\theta }_2^{″}=0.\hfill \end{array}$$

(36)

The corresponding boundary conditions are reduced to,

• when $\eta =-h_1$:

$${\mathrm{\Psi }}_1={\zeta }_1,u_1=0,{\theta }_1=1,{\varphi }_1=1,$$

(37)

when $\eta =0$:

$$\begin{array}{ccc}\hfill {\mathrm{\Psi }}_1={\mathrm{\Psi }}_2,{\mathrm{\Psi }}_1^{\prime }={\lambda }_{\epsilon }{\mathrm{\Psi }}_2^{\prime }& ,\hfill & \hfill u_1=\frac{{\lambda }_{\epsilon }}{{\lambda }_{\mu }}{u}_2,\\ \hfill {\gamma }^2{\lambda }_{\epsilon }{u}_2^{\prime }-u_1^{\prime }=\frac{\beta }{\sqrt{Da}}\frac{{\lambda }_{\epsilon }}{{\lambda }_{\mu }}{u}_2,{\theta }_1={\theta }_2& ,\hfill & \hfill {\theta }_1^{\prime }={\lambda }_{nf}{\theta }_2^{\prime },\\ \hfill N_{B1}({\varphi }_1^{\prime }-{\lambda }_{D_B}{\varphi }_2^{\prime })+N_{T1}({\theta }_1^{\prime }-{\lambda }_{D_T}{\theta }_2^{\prime })=0& ,\hfill & \hfill {\varphi }_1={\varphi }_2\end{array}$$

(38)

when $\eta =h_2$:

$${\mathrm{\Psi }}_2={\zeta }_2,u_2=0,{\theta }_2=1,{\varphi }_2=1.$$

(39)

where ${\zeta }_i=\widehat{z}{e}_0{\overline{\zeta }}_i/(k_B\widehat{T})$ is the zeta potential, $k_i$ the electro-osmotic parameter.

$$\begin{array}{ccc}& & k_i=\widehat{z}eH\sqrt{\frac{2n_0}{{\epsilon }_0{\epsilon }_{Ri}{k}_B\widehat{T}}},H{a}_i=B_{0}H\sqrt{\frac{{\sigma }_i}{{\mu }_i}},S_{e_i}=\frac{E_{y}H}{U_{ai}}\sqrt{\frac{{\sigma }_i}{{\mu }_i}},\hfill \\ & & {\mathrm{\Gamma }}_i=-\frac{H^2}{{\mu }_i{U}_{ai}}\frac{d\overline{P}}{d\overline{x}},U_{ai}=\frac{{\epsilon }_0{\epsilon }_{Ri}{E}_x{k}_B\widehat{T}}{{\mu }_i\widehat{z}e},Da=\frac{\kappa }{H^2},\gamma =\sqrt{\frac{{\mu }_{eff}}{{\mu }_2}}\hfill \\ & & N_{Bi}=\frac{{\tau }_i{D}_{Bi}(C_w-C_0)}{{\alpha }_i},N_{Ti}=\frac{{\tau }_i{D}_{Ti}(T_w-T_0)}{{\alpha }_i{T}_0},\hfill \\ & & B{r}_i=\frac{{\mu }_i{U}_{ai}^2}{k_{nfi}(T_w-T_0)}\hfill \end{array}$$

(40)

To measure the difference of physical properties, the following ratios are defined:

$$\begin{array}{ccc}& & {\lambda }_{\epsilon }=\frac{{\epsilon }_2}{{\epsilon }_1},{\lambda }_{nf}=\frac{k_{nf2}}{k_{nf1}},{\lambda }_{D_B}=\frac{D_{B2}}{D_{B1}},{\lambda }_{\mu }=\frac{{\mu }_2}{{\mu }_1}\hfill \\ & & {\lambda }_{\sigma }=\frac{{\sigma }_2}{{\sigma }_1},{\lambda }_{D_t}=\frac{D_{T2}}{D_{T1}},{\lambda }_{\alpha }=\frac{{\alpha }_2}{{\alpha }_1},{\lambda }_{\tau }=\frac{{\tau }_2}{{\tau }_1}\hfill \end{array}$$

(41)

where the physical parameters of the two regions are related as follows:

$$\begin{array}{ccc}& & U_{a2}=\frac{{\lambda }_{\epsilon }}{{\lambda }_{\mu }}{U}_{a1},S_{e_2}=\frac{\sqrt{{\lambda }_{\mu }{\lambda }_{\sigma }}}{{\lambda }_{\epsilon }}{S}_{e_1},B{r}_2=\frac{{\lambda }_{\epsilon }^2}{{\lambda }_{\mu }{\lambda }_{nf}}B{r}_1,{\mathrm{\Gamma }}_2=\frac{1}{{\lambda }_{\epsilon }}{\mathrm{\Gamma }}_1,\hfill \\ & & k_2=\frac{1}{\sqrt{{\lambda }_{\epsilon }}}{k}_1,N_{B2}=\frac{{\lambda }_{\tau }{\lambda }_{D_B}}{{\lambda }_{\alpha }}{N}_{B1},N_{T2}=\frac{{\lambda }_{\tau }{\lambda }_{D_T}}{{\lambda }_{\alpha }}{N}_{T1},\hfill \\ & & H{a}_2=\sqrt{\frac{{\lambda }_{\sigma }}{{\lambda }_{\mu }}}H{a}_1\hfill \end{array}$$

(42)

The required ratios are calculated using the values from Table 1. Such as ${\lambda }_{nf}=1$, ${\lambda }_{\mu }=1$, ${\lambda }_{\alpha }=1$, and ${\lambda }_{\tau }=0.96$, and other ratios are chosen as ${\lambda }_{\sigma }=1.2$, ${\lambda }_{\epsilon }=1.2$, ${\lambda }_{D_B}=1.2$, ${\lambda }_{D_T}=1.2$, and ${\lambda }_{\rho }=1$.

2.3. Skin Friction Coefficient and Nusselt Number

For the heat and mass transfer analyses we calculated the skin friction coefficient and Nusselt number as follows:

$$C_{fi}=\frac{{\tau }_{wi}}{\frac{1}{2}{\rho }_i{U}_{ai}^2},N{u}_i=\frac{H_i{q}_{wi}}{k_{nfi}(T_w-T_0)}$$

(43)

where $i=1,2$ denotes regions I and II, ${\tau }_{wi}$ denotes the shear stress, and $q_{wi}$ denotes the heat flux, which can be calculated using

$${\tau }_{wi}={\mu }_i{\left.\frac{\partial \overline{u_i}}{\partial \overline{z}}\right|}_{\overline{z}={(-1)}^i{H}_i},q_{wi}=-k_{nfi}{\left.\frac{\partial \overline{T_i}}{\partial \overline{z}}\right|}_{\overline{z}={(-1)}^i{H}_i}.$$

(44)

Substituting Equations (28) and (44) into Equation (43), we get

$$C_{fi}=\frac{2}{R{e}_i}{u}_i^{\prime }({(-1)}^i{h}_i),N{u}_i=-\frac{h_i}{H}{\theta }_i^{\prime }({(-1)}^i{h}_i)$$

(45)

where $R{e}_i=\frac{H{\rho }_i{U}_{ai}}{{\mu }_i}$ is the Reynolds number. The relationship between two regions’ Reynolds numbers is defined by

$$R{e}_2=\frac{{\lambda }_{\rho }{\lambda }_{\epsilon }}{{\lambda }_{\mu }^2}R{e}_1$$

(46)

where ${\lambda }_{\rho }={\rho }_2/{\rho }_1$.

3. Problem Solution

Because Equations (29), (30), (33) and (34), are all second-order ODEs, they can be solved for their exact solutions. The Matlab built-in function $dsolve$ is used to solve Equations (29), (30), (33), and (34), with respect to the boundary conditions Equations (37), (38) and (39) for ${\mathrm{\Psi }}_1$, $u_1$, ${\mathrm{\Psi }}_2$, and $u_2$. Because Equations (31), (32), (35), and (36) are nonlinear ODEs and exact solutions are difficult to obtain; thus, numerical simulations based on FDM are performed, as described in [30]. The nonlinearity of the equations is dealt with by Picard’s method. The iterative procedure stops once the following criterion is met.

$$\frac{{\sum }_{j=1}^N|F^{i+1}\left(\eta \right)-F^i\left(\eta \right)|}{{\sum }_{j=1}^N\left|F^{i+1}\left(\eta \right)\right|}\le {10}^{-8},i\ge 1,$$

(47)

where F represents the field variables $\psi$, u, $\theta$, and $\varphi$, and N represents the number of grid points. Here, ${10}^{-8}$ is the predefined tolerance error. To confirm our findings, we replicated those of Niazi and Xu [17], by setting $\beta =0$, $Da=1$, and $\gamma =1$. The comparisons for the velocity and temperature profiles are shown in Figure 2, which validates the results of the current problem. In this analysis, values of the parameters are selected based on the properties of the nanofluid given in Table 1. These values can vary depending upon the values of other parameters, to keep the system stable. For some values of parameters, such as $\beta$ and $S{e}_1$, we have referred to the papers by Kuznetsov [26] and Niazi et al. [17].

Figure 3 illustrates the velocity profiles calculated for various different values of $\beta$ and $Da$. The interface is located at $\eta =0$. According to the analysis in [26], we chose values for $\beta$ that range between $-0.8$ and $+0.8$, and a Darcy number of the order of ${10}^{-1}$ or less. Figure 3a demonstrates that a change in stress can fundamentally alter the velocity profiles. When the stress at the interface increases, the slope of the tangent to the velocity distribution at $\eta =0$ changes dramatically. By gradually increasing the coefficient $\beta$, which accounts for the stress jump, the velocity is reduced noticeably. When $\beta$ is negative, this apparent impact is particularly strong. Additionally, when the stress jump coefficient is varied, there is no change in the velocity profile near the upper wall. The Darcy number’s effect on the velocity profile is depicted in Figure 3b. The curve computed for $Da={10}^{-2}$ contains three segments. One portion is contained within the momentum boundary layer adjacent to the boundary at $\eta =-1$, while the other portion is contained within the momentum boundary layer adjacent to the interface at $\eta =0$. As per the classical Darcy law, the fluid velocity increases inbetween two boundary layers, but stays unchanged in the porous layer. Additionally, as it enters the porous layer, the velocity decreases more rapidly in this third section. Similarly, the curves corresponding to $Da={10}^{-3}$ and $Da={10}^{-4}$ are almost identical; the difference is simply not visible due to the low velocity in the porous layer. There is no point on the curve $Da={10}^{-1}$ where the velocity is constant. This is because, as the Darcy number increases, the width of the momentum boundary layers decreases.

The influences of the physical ratios on the flow characteristics are displayed in Figure 4. It is observed in Figure 4a that, with the increase in the ratio of electric conductivity (${\lambda }_{\epsilon }$), the average velocity increases in Region I. At the same time, it decreases in Region II with a porous medium. The fluid in Region I conducts electricity better than in Region II. Further, in Figure 4b, an increase in the viscosity ratio (${\lambda }_{\mu }$) decreases the velocity throughout the channel. In the case of a larger viscosity ratio, the velocity is smaller in Region I and Region II. The reason for these phenomena, is that when the viscosity ratio ${\lambda }_{\mu }>1$, the fluid viscosity in Region I is greater than that in Region II, resulting in a larger value of the average velocity in Region I.

The significant influence of stress jump condition coefficient and Darcy number on the temperature profile was examined, and is shown in Figure 5. It can be observed that the temperature decreases throughout the channel for larger values of $\beta$, as displayed in Figure 5a. This figure demonstrates the significant effect of the stress jump on the non-dimensional temperature profile. A peak in the $\theta \left(\eta \right)$ is observed for a considerably smaller value of $\beta =-0.8$, and this curve flattens for a very large value of $\beta$. This illustrates that the increase in the stress jump coefficient reduces the temperature throughout the channel. While an adverse behavior in $\theta \left(\eta \right)$ is seen in the case of the Darcy number, as shown in Figure 5b. For lower values of the Darcy number, the temperature profile decreases significantly in the two regions. It is also observed that the temperature profile is indeed more significant in Region II than in Region I, which in turn depicts that the heat transfer rate is higher in the porous layer. Figure 6 shows the variation in the temperature profile for distinct values of the Brinkman number and viscosity ratio. An increase in Brinkman number $B{r}_1$, tends to increase the temperature profile, as in Figure 6a. A higher value of $B{r}_1$, slows the conduction of heat produced and hence the temperature rise is more considerable. The viscosity ratio shows an opposite trend on $\theta \left(\eta \right)$ as compared to the Brinkman number, as given in Figure 6b. As the value of the viscosity ratio, ${\lambda }_{\mu }$, increases, the temperature profile decreases throughout the channel. Physically, as the value of $\lambda \mu$ increases, so does the amount of molecular conduction in the second region. As a result, the temperature of Region II decreases, which results in a decrease in the temperature throughout the channel, as illustrated in Figure 6b. However, when the viscosity ratio, ${\lambda }_{\mu }$, is minor (${\lambda }_{\mu }\le 1$), the position of the maximal value for $\theta \left(\eta \right)$ shifts towards Region II. For larger values of ${\lambda }_{\mu }$, the position shifts towards Region I. This shift occurs because the fluid interface must satisfy the boundary condition for continuous thermal flux.

Figure 7 illustrates the evolution of the nanoparticle volume fraction $\varphi \left(\eta \right)$ as the stress jump coefficient $\left(\beta \right)$ and Darcy number $\left(Da\right)$ increase. As shown in Figure 7a, the volume fraction of nanoparticles decreases rapidly as $\beta$ decreases, particularly for negative values of $\beta$. Nevertheless, an adverse behavior is seen in the $\varphi \left(\eta \right)$ profile for the case of the Darcy number, as given in Figure 7b. The increase in Darcy number causes a significant decrease in the nanoparticle profile. In contrast, smaller values of the Darcy number have a minimal impact on the $\varphi \left(\eta \right)$ and give almost a flat curve for $Da={10}^{-3}$ and $Da={10}^{-4}$. The impacts of the physical ratios and Brinkman number on the nanoparticle volume fraction, are shown in Figure 8. It is observed that the influence of $B{r}_1$ and ${\lambda }_{\mu }$ on $\varphi \left(\eta \right)$, show opposite trends. Figure 7a shows that the $\varphi \left(\eta \right)$ decreases as the Brinkman number increases, due to the increased fluid viscosity; while increasing the value of ${\lambda }_{\mu }$ accelerates the movement of nanoparticles toward the upper wall, this results in a decrease in the nanoparticle volume fraction, as illustrated in Figure 8b.

The variations in skin friction coefficient, $C_{fi}$, with ${\kappa }_1$, $\beta$, and $Da$ are presented in Figure 9 and Figure 10, respectively. Since it has hitherto been observed that the $\beta$ and $Da$ produce a significant effect on the velocity profile, given in Figure 3, this pattern also holds true for the fluctuation of $C_{fi}$, but the orientation on the upper wall is the inverse of the direction on the lower wall. When $\beta <0$, the local skin friction ($C_{f1}$) rises, for increasing values of ${\kappa }_1$, and this increase becomes increasingly evident for larger values of ${\kappa }_1$. On the other hand, for $\beta \ge 0$, there is a slight increase in skin friction ($C_{f1}$) for larger values of ${\kappa }_1$ but for smaller values of ${\kappa }_1$ this increase is negligible. At the upper wall, variation in the $\beta$ parameter has no effect on the skin friction coefficient ($C_{f2}$) when the electro-osmotic parameter is changed. Figure 10 illustrates the increase in the skin friction coefficient on the lower wall and the decrease on the upper wall, as ${\kappa }_1$ increases. This is because the coordinates are set in such a way that at the interface between the two layers, the signs on the top and bottom walls are reversed. An increase in the Darcy number increases the $C_{f1}$ significantly, but for smaller values of the Darcy number, this increase is not prominent. While at the upper wall, the decrease in $C_{f2}$ is evident for the increasing values of Darcy number. The influence of Nusselt number ($N{u}_i$) with $\beta$ and $Da$, for several values of ${\kappa }_1$, are shown in Figure 11 and Figure 12, respectively. As can be seen from these figures, increasing the value of ${\kappa }_1$ decreases the Nusselt number on the top wall, while increasing it on the bottom wall. Physically, the increase in ${\kappa }_1$ reduces the EDL effect, which leads to the enhancement of the fluid motion. Thus, it causes more heat conduction than heat convection, causing the decrease in the Nusselt number at the upper wall. On the other hand, for larger values of $\beta$ and $Da$, the effect on $N{u}_i$ is more evident, both at the upper and lower wall. As a result, calculations that do not take the increase in stress into account, may suffer a significant loss of accuracy.

4. Conclusions

In the EDL effect, physical analysis was performed on two-layer nanofluids in a microchannel partially filled with porous medium. The mathematical model for the two-layer nanofluid flow is developed using Buongiorno’s model. At the interface, Ochoa-Tapia and Whitaker’s proposed jump boundary condition is used, to map the Brinkman-extended Darcy equation to the Stokes’ equation. The finite difference method is employed to solve and investigate the nonlinear system of differential equations. Exact solutions are obtained for the electrostatic potential and velocity. Different from the work of Niazi and Xu [17], we have considered a partially filled microchannel with a porous medium and used the stress jump condition at the interface. Two momentum boundary layers are formed in the porous region, to ensure uniform permeability of the porous medium, as illustrated graphically. One of these boundary layers forms immediately adjacent to the impermeable boundary, while the other forms immediately adjacent to the interface. The fluid velocity is constant between these momentum boundary layers. For this reason, the Darcy number and stress jump coefficient show significant variations with velocity, temperature, and nanoparticle concentration. In addition, increasing the physical ratio, i.e., viscosity, decreases the velocity and temperature profiles. Simultaneously, an adverse trend is observed in the volume fraction profile of the nanoparticles. Additionally, the effect of the Darcy number and stress jump coefficient on the skin friction, and the Nusselt number on the upper and lower microchannel walls, are visible. As a result, it is concluded that the stress jump boundary condition is critical for solving fluid flow problems in a wide variety of practical applications.