*Article* **Transient Two-Layer Electroosmotic Flow and Heat Transfer of Power-Law Nanofluids in a Microchannel**

**Shuyan Deng \* and Tan Xiao**

Institute of Architecture and Civil Engineering, Guangdong University of Petrochemical Technology, Maoming 525011, China; xiaotan@gdupt.edu.cn

**\*** Correspondence: sydeng4-c@my.cityu.edu.hk; Tel.: +86-173-7689-1017

**Abstract:** To achieve the optimum use and efficient thermal management of two-layer electroosmosis pumping systems in microdevices, this paper studies the transient hydrodynamical features in twolayer electroosmotic flow of power-law nanofluids in a slit microchannel and the corresponding heat transfer characteristics in the presence of viscous dissipation. The governing equations are established based on the Cauchy momentum equation, continuity equation, energy equation, and power-law nanofluid model, which are analytically solved in the limiting case of two-layer Newtonian fluid flow by means of Laplace transform and numerically solved for two-layer power-law nanofluid fluid flow. The transient mechanism of adopting conducting power-law nanofluid as a pumping force and that of pumping nonconducting power-law nanofluid are both discussed by presenting the two-layer velocity, flow rates, temperature, and Nusselt number at different power-law rheology, nanoparticle volume fraction, electrokinetic width and Brinkman number. The results demonstrate that shear thinning conducting nanofluid represents a promising tool to drive nonconducting samples, especially samples with shear thickening features. The increase in nanoparticle volume fraction promotes heat transfer performance, and the shear thickening feature of conducting nanofluid tends to suppress the effects of viscous dissipation and electrokinetic width on heat transfer.

**Keywords:** transient two-layer flow; electroosmotic flow; power-law nanofluid; heat transfer; Laplace transform; nanoparticle volume fraction

#### **1. Introduction**

It is well known that in microchannels the contact between the electrolyte solution and the solid surface of the channel wall leads to the rearrangement of charged ions, inducing an electric double layer (EDL) near the channel wall. In the presence of EDL, a layer of conducting fluid under a tangentially-applied electric field moves forward, forming electroosmotic flow (EOF); this phenomenon is called electroosmosis. Due to such favorable attributes as its ease of integration, plug-like profile, and the independence of its non-mechanical parts, the electroosmosis pumping mechanism has become a common transport phenomena in microfluidic devices [1]. In order to meet the growing demand for electroosmosis-based applications, a large number of works have theoretically studied EOF from different point of view. The transport characteristics of EOF in containers with different geometries, including slit microchannels [2], microtubes [3,4], rectangular microchannels [5], elliptic microchannels [6], and T-shaped microchannels [7] have been investigated. Several working liquids in microdevices, such as biomedical lab-on-a-chip, show nonlinear rheological behaviors, which is where the use of non-Newtonian fluid modeling becomes relevant. The nonlinear relationship between the shear stress and shear rate has been carefully treated using the power-law model [8,9], Casson model [10], Maxwell model [11], Carreau model [12], etc. The power-law model was first proposed by Das and Chakraborty [8] to describe the rheological behavior of blood, which has received great attention due to its wide coverage and the simple rheological relation [2,13]. The

**Citation:** Deng, S.; Xiao, T. Transient Two-Layer Electroosmotic Flow and Heat Transfer of Power-Law Nanofluids in a Microchannel. *Micromachines* **2022**, *13*, 405. https://doi.org/10.3390/ mi13030405

Academic Editors: Jin-yuan Qian and Kwang-Yong Kim

Received: 19 January 2022 Accepted: 26 February 2022 Published: 1 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

various aspects of power-law fluid for EOF have been discussed; the power-law model incorporates the shear thinning rheological behavior encountered in DNA solutions, and shear thickening rheological behavior encountered in cornstarch solution [14]. Recently, external environmental effects on EOF, such as a rotating frame or peristalsis, have been considered. In microchannel flow, a rotating environment induces Coriolis force, which causes a secondary flow; this is applied in biofluid transportation, drug delivery, and DNA analysis [15]. The theoretical analysis of rotating EOF was first studied by Chang and Wang [16], and has subsequently been extended in literature [17,18]. In order to obtain a comprehensive understanding of the intricate mechanism of rotational flow for biofluid flow, Kaushik et al., engaged in a transient analysis of the rotational electrohydrodynamics of power-law fluids under the effect of EDL [19]. Moreover, in the application of EOF to biomedical and biochemical analysis, peristalsis is introduced to assist in the EOF of biofluids, and thus electroosmosis-modulated peristaltic flow has recently become a frequently-studied research topic [20].

Owing to the application of external electric fields in electroosmosis-driven flow systems, the applied electric voltage leads to an inherent byproduct of the Ohmic resistance of electrolytes, namely, the Joule heating effect. Joule heating-induced heat transfer exerts an influence on transport performance by altering the electric properties of working liquids, especially for certain thermally-liable samples, which has been widely discussed in numerous works [21–25]. To optimize the hydrodynamic transport process and minimize the Joule heating effect, combined electric and magnetic fields are applied to working liquids in order to improve the actuation mechanism in microfluidics, which has the advantage of lower voltage operation, convenient manufacture, and the independence of moving parts [26,27]. On the other, to promote heat transfer and reduce entropy generation in the heat exchange equipment of microfluidics, nanofluid is created by adding nanosized metal particles, which possesses boosted thermal conductivity compared to conventional pure fluids, to the sample [28–30]. Nanofluid flow has been extensively applied in different fields, as it has none of the usual drawbacks such as sedimentation, blockage, and pressure drop [31–33]. AI2O3–water Nanofluid is used for cooling microprocessors or other microelectronic components due to its enhanced thermal conductivity [34], which exhibits shear thinning rheological behavior in certain ranges of nanoparticle volume fraction [35]. Moreover, in order to provide a better understanding of blood flow and other non-Newtonian biological flows in biomicrofluidic chips (such as pseudoplastic aqueous nanoliquid flow driven by electroosmosis and peristalsis and Cu/CuO–blood microvascular nanoliquid flow under thermal, microrotation, and electromagnetic field effects) were studied in [36,37], respectively. Carboxymethyl cellulose (CMC) water with γ-AI2O3, TiO<sup>2</sup> and CuO particles has been experimentally investigated for the achievement of efficient thermal management in microelectronics [38]. A comprehensive literature survey of topics in nanofluid flow indicates that the viscosity of several nanofluids mentioned above shows a nonlinear dependence on the shear rate and volume fraction of nanoparticles; thus, the power-law nanofluid model is proposed to precisely describe the rheological behavior of such nanofluids [39–42]. Regarding the power-law nanofluid, the heat transfer characteristics in magneto-hydrodynamic flow [43], convection flow [44], and EOF [45–47] have been extensively investigated.

The great advent of technologies in microelectrical mechanical systems (MEMS) enables the delivery, mixing, or separation of multi-liquids at microscales. However, certain liquids, for instance oil, blood, and ethanol, have low electrical conductivity (<10−<sup>6</sup> S m−<sup>1</sup> ) and are defined as nonconducting fluids [48], which fail to be driven by electroosmosis. Furthermore, for certain liquids the applied electric voltage leads to undesirable problems such as the generation of gases, fluctuation of PH value, or electrochemical decomposition. In this context, Brask et al. developed a two-layer flow system where the conducting fluid driven by electroosmotic force is adopted as driving mechanism to drag the nonconducting fluid; this has gained great attention in recent decades [49]. The steady hydrodynamic behaviors of two-Newtonian fluid EOF in a rectangular microchannel [50], two-powerlaw fluid combined electroosmotics with pressure driven flow in a microtube [51], and Newtonian–Casson fluid EOF in a microtube [48] have all been theoretically studied in this context. In terms of transient hydrodynamical behaviors, Gao et al. characterized the transient two-layer EOF of Newtonian fluids in a rectangular microchannel by presenting analytical velocities and flow rates at different viscosities and different electroosmotic properties [52]. Su et al., presented semi-analytical velocities for two-layer combined electroosmotic and pressure driven flow of Newtonian fluids in a slit microchannel at different electric and hydrodynamic parameters [53]. Time periodic transport characteristics of two-Newtonian liquid combining electroosmotic and pressure driven flows in a microtube have been studied numerically [54]. In order to improve transport efficiency and reduce the Joule heating effect in a two-layer pumping system, a magnetic field can be applied in addition to the pressure gradient to form a magneto-hydrodynamic EOF; the corresponding entropy generation analysis has been conducted as well [55,56]. Two-layer EOF assisted by peristalsis force was proposed by Ranjit et al., who analyzed entropy generation and heat transfer in such cases [57]. External factors such as a rotating environment [58] or varying wall shapes together with zeta potential [59] have been considered in two-layer electroosmotic systems, and the resulting influence on hydrodynamic behavior has been discussed. Furthermore, because of its desirable thermal conductivity properties, nanofluids can be applied in two-layer mixed convection flows, which are characterized by the power-law nanofluid model; the outcomes can help with the promotion of heat transfer performance [60]. Entropy generation and heat transfer in immiscible EOF of two conducting power-law nanofluid flows through a microtube have been analyzed by computing their temperature, Nusselt number, and entropy generation at different nanoparticle volume fractions and different rheological and electroosmotic properties [61]; the rheological effect of the peripheral fluid plays a dominant role in thermal performance as compared to that of inner fluid.

The application of chemical mixing/separation in thermofluidic micropumps has become increasingly frequent; therefore, the corresponding microscale cooling and heat exchangers need to be carefully designed as working liquids combined with nanoparticles show nonlinear rheological behavior in nature. To the authors' best knowledge, the underlying mechanism of the transient transport process as it develops from an unsteady to a steady state in two-power-law nanofluid EOF in a slit microchannel remains to be discovered. Therefore, this paper studies transient two-layer flow with one layer of conducting powerlaw nanofluid and one layer of nonconducting power-law nanofluid in a slit microchannel, with consideration of Joule heating and viscous dissipation. The governing equations are established based on the Cauchy momentum equation, continuity equation, energy equation, and power-law nanofluid rheological relation, which are analytically solved for two-Newtonian fluid flow and numerically solved for two-power-law nanofluid flow. For the hydrodynamic aspect, the mechanisms involved in using power-law nanofluid as a pumping force as well as those of pumping power-law nanofluid as the system develops from unsteady to a steady state are carefully discussed by presenting the time evolution of velocity and flow rate at different parameters. For the thermal aspect, in order to guarantee efficient thermal performance, the heat transfer characteristics arising from the interplay of the nanoparticles, power-law rheological behavior, viscous dissipation, and electrokinetic effects in a two-layer system are analyzed theoretically. The results are relevant for assisting in determining the operating parameters for optimal performance of microdevices characterized by multi-fluid delivery, mixing, or cooling.

### **2. Problem Formulation**

#### *2.1. Electric Potential Distribution*

A two-layer immiscible EOF of power-law nanofluids in a slit microchannel is considered where the power-law nanofluid in layer I is conducting and the power-law nanofluid in layer II is nonconducting (Figure 1). The AI2O<sup>3</sup> nanoparticles are suspended and uniformly distributed in a two-layer power-law base fluid system where carboxymethyl

cellulose–water (CMC-water) can be represented by shear thinning fluid [38] and, without loss of generality, the base fluids fall into the categories of a shear thinning fluid and a shear thickening one. According to previously published work [50], the zeta potential difference near the two-liquid interface (*y* \* = 0) is negligible. The heights of layer I and layer II are represented by *H*. In the two-layer system, the EDL forms near the channel wall within the region of layer I, which creates the electric potential distribution *ψ* \* . When a uniform electric field *E* ∗ *z* is tangentially exerted across layer I, the conducting power-law nanofluid moves forward under the electroosmotic force due to the existence of EDL near the lower channel wall, which drags the nonconducting power-law nanofluid in layer II via the interfacial viscous stress, eventually forming a two-layer EOF. potential distribution *ψ*\* . When a uniform electric field \* *Ez* is tangentially exerted across layer I, the conducting power-law nanofluid moves forward under the electroosmotic force due to the existence of EDL near the lower channel wall, which drags the nonconducting power-law nanofluid in layer II via the interfacial viscous stress, eventually forming a two-layer EOF.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 4 of 27

A two-layer immiscible EOF of power-law nanofluids in a slit microchannel is

considered where the power-law nanofluid in layer I is conducting and the power-law

nanofluid in layer II is nonconducting (Figure 1). The AI2O3 nanoparticles are suspended

and uniformly distributed in a two-layer power-law base fluid system where

carboxymethyl cellulose–water (CMC-water) can be represented by shear thinning fluid

[38] and, without loss of generality, the base fluids fall into the categories of a shear

thinning fluid and a shear thickening one. According to previously published work [50],

the zeta potential difference near the two-liquid interface (*y*\* = 0) is negligible. The

heights of layer I and layer II are represented by *H*. In the two-layer system, the EDL

forms near the channel wall within the region of layer I, which creates the electric

**2. Problem Formulation**

*2.1. Electric Potential Distribution*

**Figure 1.** Schematic of two-layer EOF of power-law nanofluids in a slit microchannel. **Figure 1.** Schematic of two-layer EOF of power-law nanofluids in a slit microchannel.

It is assumed that the zeta potential *ζ*\* is small and the EDL thickness is far less than It is assumed that the zeta potential *ζ* \* is small and the EDL thickness is far less than the height of microchannel; thus, the EDLs near the channel walls will not overlap. Eventually, based on electrostatic theory, the electric potential distribution *ψ* \* is governed by the well-known Poisson–Boltzmann (P-B) equations [50]

$$\frac{d^2\phi^\*}{dy^{\*2}} = -\frac{\rho\_t}{\varepsilon} \tag{1}$$

ρ

ε

0 0

According to the established model in the scientific literature [32,62], when

nanoparticles with an order of nm and with a volume fraction of *ϕ* ≤ 10% are distributed

in the channel with μm-sized height, it is physically reasonable to assume that the EDLs

around the nanoparticles are rather small compared to the EDLs near the channel walls,

and can thus be neglected, as well as that there is no electrophoretic force on the

nanoparticles. Correspondingly, the nanoparticles have no influence on the local electric

charge density *ρe*, which can be governed by the P-B equations, and their influence on

liquid property can be incorporated into the effective viscosity and effective thermal

conductivity of the nanofluid [28–30] such that their importance can be attached to the

0

*z e*

*B*

*k T*

\*

0

ψ

$$\rho\_{\varepsilon} = -2n\_0 z\_0 \varepsilon\_0 \text{sech}\left(\frac{z\_0 e \psi^\*}{k\_B T\_0}\right) \tag{2}$$

2 \* \*2 *d <sup>e</sup> dy* ψ = − (1) 2 sinh *<sup>e</sup> nze* ρ = − (2) According to the established model in the scientific literature [32,62], when nanoparticles with an order of nm and with a volume fraction of *φ* ≤ 10% are distributed in the channel with µm-sized height, it is physically reasonable to assume that the EDLs around the nanoparticles are rather small compared to the EDLs near the channel walls, and can thus be neglected, as well as that there is no electrophoretic force on the nanoparticles. Correspondingly, the nanoparticles have no influence on the local electric charge density *ρe* , which can be governed by the P-B equations, and their influence on liquid property can be incorporated into the effective viscosity and effective thermal conductivity of the nanofluid [28–30] such that their importance can be attached to the effect of power-law

feature and the effect of nanoparticles on the hydrodynamic and thermal characteristics of transient two-layer flow.

The P-B equations are subject to the following boundary conditions:

$$\frac{d\psi^\*}{dy^\*}\Big|\_{y^\*=0} = 0, \left.\psi^\*\right|\_{y^\*=H} = \zeta^\* \tag{3}$$

where *ε* denotes the dielectric constant, *n*<sup>0</sup> denotes the ion density, *z*<sup>0</sup> denotes the valence, *e* is the electron charge, and *k<sup>B</sup>* and *T*<sup>0</sup> represent the Boltzmann constant and the absolute temperature, respectively.

Introducing the nondimensional variables *ψ = ez*0*ψ* \*/(*kBT*0), *ζ* = *ez*0*ζ* \*/(*kBT*0), and *y* = *y* \*/*H*, *K* = *κH* with *κ* <sup>2</sup> = 2*e* <sup>2</sup>z<sup>0</sup> <sup>2</sup>*n*0/(*εkBT*0), and applying Debye–Hückel approximation (sinh*ψ* ≈ *ψ* [50]) to Equations (1)–(3) when |*ζ* \*<sup>|</sup> <sup>≤</sup> 0.025 V, the nondimensional versions of the P-B equations can be rewritten as

$$\frac{d^2\psi}{dy^2} = \mathcal{K}^2 \psi$$

$$\left. \frac{d\psi}{dy} \right|\_{y=0} = 0, \left. \psi \right|\_{y=1} = \zeta \tag{5}$$

Solving Equations (4) and (5) yields the electric potential distribution, as below:

$$\psi = \frac{\zeta \cosh(Ky)}{\cosh(K)}\tag{6}$$

With Equation (6), the electroosmotic force driving the conducting nanofluid can be obtained.

#### *2.2. Two-Layer Velocity Distribution and Flow Rates*

Focusing on the hydrodynamical aspects of transient two-layer EOF of power-law nanofluids in a slit microchannel, the governing equations for velocity distribution are represented by the Cauchy momentum equation and the continuity equation, as below:

$$\nabla \cdot \stackrel{\rightarrow}{v}^\* = \mathbf{0} \tag{7}$$

$$
\rho \left[ \frac{\partial \overrightarrow{\boldsymbol{v}}^{\*}}{\partial t^{\*}} + \left( \overrightarrow{\boldsymbol{v}}^{\*} \cdot \nabla \right) \overrightarrow{\boldsymbol{v}}^{\*} \right] = \nabla \cdot \overrightarrow{\boldsymbol{\tau}} + \overrightarrow{\boldsymbol{f}} - \nabla p \tag{8}
$$

where <sup>→</sup> *v* is the velocity vector, *ρ* is the liquid density, *t* \* is the time, τ is the shear stress, → *f* denotes the body force vector, and ∇*p* denotes the pressure gradient. The channel is openended and no pressure gradient is induced. The following assumptions are made for the purpose of analysis: (1) the properties of the nanofluids are independent of the external electric field, ion concentration, and temperature [48,50]; (2) the two-layer flow is immiscible, laminar, and incompressible, and the two-liquid interface remains distinguishable [48,50]; and (3) the gravity force and buoyancy force of the nanofluids are neglected [32]. As a result, there is only velocity component along *z* \* direction *v \* i* (*y* \* ,*t* \* ), with *i* = 1,2 and the body force equal to the electroosmotic force, *f<sup>z</sup>* = *Ezρ<sup>e</sup>* . In a system with two power-law nanofluids flowing through a slit microchannel, the shear stress of a power-law nanofluid is *τ<sup>i</sup>* = *ηeffi*·*∂v \* <sup>i</sup>*/*∂y* \* , where *ηeffi* implies the effective dynamic viscosity of the power-law nanofluid, which nonlinearly depends on the shear rate *∂v \* <sup>i</sup>*/*∂y* \* and the nanoparticle volume fraction *φ*, namely, *ηe f f i* = *m*<sup>∗</sup> 0 /(1 + *φ*) 5/2 · *∂v* ∗ *i* /*∂y* ∗ *ni*−1 [13,28,43,51], where *m*<sup>0</sup> \* is the consistency viscosity coefficient, *n<sup>i</sup>* is the flow behavior index, the subscript *i* = 1 represents the conducting nanofluid in layer I, and *i* = 2 represents the nonconducting nanofluid in layer II. Note that *n<sup>i</sup>* < 1 corresponds to shear thinning base fluid, *n<sup>i</sup>* = 1 corresponds to Newtonian base fluid and *n<sup>i</sup>* > 1 corresponds to shear thickening base fluid. Accordingly, the modified Cauchy momentum equation for the conducting power-law nanofluid in layer I under the electroosmotic force is expressed as

$$\frac{\partial v\_1^\*}{\partial t^\*} = \frac{m\_0^\*}{\left(1+\phi\right)^{5/2}\rho\_1} \frac{\partial}{\partial y^\*} \left( \left| \frac{\partial v\_1^\*}{\partial y^\*} \right|^{n\_1-1} \frac{\partial v\_1^\*}{\partial y^\*} \right) + \frac{1}{\rho\_1} \rho\_\varepsilon E\_z^\* \text{ for } 0 \le y^\* \le H \tag{9}$$

The modified Cauchy momentum equation of nonconducting power-law nanofluid in layer II is expressed as

$$\frac{\partial v\_2^\*}{\partial t^\*} = \frac{m\_0^\*}{(1+\phi)^{5/2}\rho\_2} \frac{\partial}{\partial y^\*} \left( \left| \frac{\partial v\_2^\*}{\partial y^\*} \right|^{n\_2-1} \frac{\partial v\_2^\*}{\partial y^\*} \right) \text{ for } -H \le y^\* \le 0 \tag{10}$$

The velocities and shear stresses of the two-layer liquid satisfy the matching conditions, the velocities at the channel walls satisfy the no-slip condition, and the two-layer flow is initially set as motionless [55], thusly:

$$\left.v\_1^\*\right|\_{y^\*=0} = \left.v\_2^\*\right|\_{y^\*=0'} \eta\_{eff1} \cdot \frac{\partial v\_1^\*}{\partial y^\*}\Big|\_{y^\*=0} = \eta\_{eff2} \cdot \frac{\partial v\_2^\*}{\partial y^\*}\Big|\_{y^\*=0} \left.v\_1^\*\right|\_{y^\*=H} = 0 \,, \left.v\_1^\*\right|\_{y^\*=-H} = 0 \,, \left.v\_i^\*\right|\_{i^\*=0} = 0 \tag{11}$$

With the introduction of the nondimensional variables *t* = *t* \**m*0/(*ρ*1*H*<sup>2</sup> ), *v<sup>i</sup>* = *v* \* *<sup>i</sup>*/*U*, and *y* = *y* \*/*H*, by replacing Equations (2) and (6) with Equations (9)–(11), the nondimensional versions of the governing equations can be obtained as follow

$$\frac{\partial v\_1}{\partial t} = \frac{m\_1}{m\_0} \frac{\partial}{\partial y} \left( \left| \frac{\partial v\_1}{\partial y} \right|^{n\_1 - 1} \frac{\partial v\_1}{\partial y} \right) - \frac{GE\_2 \zeta}{\cosh(K)} \cosh(Ky) \text{ for } 0 \le y \le 1 \tag{12}$$

$$\frac{\partial v\_2}{\partial t} = \rho\_r \frac{m\_2}{m\_0} \frac{\partial}{\partial y} \left( \left| \frac{\partial v\_2}{\partial y} \right|^{n\_2 - 1} \frac{\partial v\_2}{\partial y} \right) \text{ for } -1 \le y \le 0 \tag{13}$$

$$\left|v\_{1}\right|\_{y=1} = v\_{2}\big|\_{y=-1} = 0, \left.v\_{1}\right|\_{y=0} = v\_{2}\big|\_{y=0'} m\_{1} \left|\frac{\partial v\_{1}}{\partial y}\right|^{n\_{1}-1} \cdot \left.\frac{\partial v\_{1}}{\partial y}\right|\_{y=0} = m\_{2} \left|\frac{\partial v\_{2}}{\partial y}\right|^{n\_{2}-1} \cdot \left.\frac{\partial v\_{2}}{\partial y}\right|\_{y=0}, v\_{1}|\_{t=0} = 0\tag{14}$$

where *m<sup>i</sup>* = (*U*/*H*) *ni*−1 *m*0/(1 − *φ*) 5/2 , *ρ<sup>r</sup>* = *ρ*2/*ρ*1, *G* = 2*zen*0*ζ* \*/(*ρ*1*U*<sup>2</sup> ), and *E<sup>z</sup>* = *E* \* *<sup>z</sup>HRe*/*ζ* \* , with *Re* = *ζ* \**UH*/*m*0, *U* the reference velocity.

With the nondimensional transient velocities *v*<sup>1</sup> and *v*<sup>2</sup> solved as in Equations (12)–(14), the transient flow rates are defined as follows:

$$Q\_1 = \int\_0^1 v\_1 dy \,\,\, Q\_2 = \int\_{-1}^0 v\_2 dy \,\,\tag{15}$$

As time elapses, the transient velocities for layer I and layer II, namely, *v*<sup>1</sup> and *v*2, reach steady status, and are then expressed as *vs*1(*y*) = lim *t*→∞ *v*1(*y*, *t*) and *vs*2(*y*) = lim *t*→∞ *v*2(*y*, *t*). To compare the flow rate of the conducting nanofluid in layer I (flow rate I) and the nonconducting nanofluid in layer II (flow rate II) with different parameters, the steady flow rate ratio is defined as

$$Q\_r = \frac{\int\_0^1 v\_{s1} dy}{\int\_{-1}^0 v\_{s2} dy} \tag{16}$$

#### *2.3. Two-Layer Temperature Distribution and Heat Transfer Rate*

With the steady velocity distribution obtained from Equations (12)–(14), the temperature distribution for the thermally fully developed two-layer flow can be determined from the following energy equation:

$$(\rho c\_p)\_{eff} \left(\frac{\partial T}{\partial t^\*} + \overrightarrow{\boldsymbol{v}}\_s^\* \cdot \nabla T\right) = k\_{eff} \nabla^2 T + \lambda \boldsymbol{E}\_z^{\*2} + \eta\_{eff} \boldsymbol{\Phi}^\* \tag{17}$$

where *T* denotes the temperature distribution, <sup>→</sup> *v <sup>s</sup>* means the steady velocity vector, *c<sup>p</sup>* means the specific heat at constant pressure, *k* is the thermal conductivity, *λ* is the electric conductivity of base fluid, Φ\* measures the viscous dissipation effect, and the subscript *eff* means the nanofluid.

The assumption that the two-layer flow is fully thermally developed leads to the vanishing of the unsteady part of Equation (17), *∂T*/*∂t* \* , hence producing the following energy equations for the conducting nanofluid and nonconducting nanofluid, respectively, along with their corresponding boundary conditions [26,55]:

$$\left(\rho c\_p\right)\_{eff} v\_{s1}^\* \frac{\partial T\_1}{\partial z^\*} = k\_{eff} \frac{d^2 T\_1}{dy^{\*2}} + \lambda E\_z^2 + \eta\_{eff1} \left(\frac{dv\_{s1}^\*}{dy^\*}\right)^2 \text{ for } 0 \le y^\* \le H \tag{18}$$

$$\left(\rho c\_p\right)\_{eff} 2v\_{s2}^\* \frac{\partial T\_2}{\partial z^\*} = k\_{eff} 2 \frac{d^2 T\_2}{dy^{\*2}} + \eta\_{eff2} \left(\frac{dv\_{s2}^\*}{dy^\*}\right)^2 \text{ for} - H \le y^\* \le 0 \tag{19}$$

$$\left.T\_1\right|\_{y^\*=0} = T\_2\big|\_{y^\*=0} \\ k\_{eff1} \left.\frac{dT\_1}{dy^\*}\right|\_{y^\*=0} = k\_{eff2} \left.\frac{dT\_2}{dy^\*}\right|\_{y^\*=0} \\ \left.T\_1\right|\_{y^\*=H} = T\_w \\ \left.T\_2\right|\_{y^\*=-H} = T\_w \\ \tag{20}$$

where *T<sup>w</sup>* means the temperature at the channel wall, subscript *i =* 1 implies the conducting nanofluid, and *i =* 2 implies the nonconducting nanofluid. Regarding the thermal properties of power-law nanofluids, the model of Choi and Yu has been applied [30,63] as it is capable of predicting the thermal conductivity of nanoliquids suspended with various kind of nonspherical nanoparticles, namely, (*ρcp*)*effi*= *φ*(*ρcp*)*<sup>p</sup>* + (1 − *φ*)(*ρcp*)*<sup>b</sup>* and *ke f f <sup>i</sup>* = *kp*+2*kbi*+2(*kp*−*kbi* )(1+*ω*) 3 *φ kp*+2*kbi*−2(*kp*−*kbi* )(1+*ω*) 3 *φ kbi* , where *ω* represents the ratio of nanolayer thickness to the original particle radius and the subscripts *p* and *b* mean nanoparticles and base fluid, respectively. The left-hand side of Equation (18) measures the heat generation due to axial conduction, while the right-hand sides of Equation (18) measure the heat generation caused by heat diffusion, heat generation from Joule heating, and heat generation caused by viscous dissipation. Imposing the constant heat flux boundary condition *qw*≡*const* for the fully thermally developed flow above, namely, *∂*[*(T<sup>w</sup>* − *T<sup>i</sup>* )/(*T<sup>w</sup>* − *Tm*)]/*∂y* \* = 0, leads to *∂T*1/*∂y* \* = *∂T*2/*∂y* \* <sup>=</sup> *dTw*/*dy*\* <sup>=</sup> *dTm*/*dy*\*≡*const*, in which *<sup>T</sup><sup>m</sup>* implies the mean temperature [26,55]. Furthermore, the overall energy balance condition over an elemental control volume results in

$$\frac{dT\_m}{dz^\*} = \frac{2q\_w + \lambda H E\_z^2 + \frac{m\_0}{(1-\rho)^{5/2}} \left( \int\_0^H \left| \frac{\partial v\_{s1}^\*}{\partial y^\*} \right|^{n\_1 - 1} \frac{\partial v\_{s1}^\*}{\partial y^\*} dy^\* + \int\_{-H}^0 \left| \frac{\partial v\_{s2}^\*}{\partial y^\*} \right|^{n\_2 - 1} \frac{\partial v\_{s2}^\*}{\partial y^\*} dy^\* \right)}{H \left( \rho c\_p \right)\_{eff1} v\_{m\text{s} 1}^\* + H \left( \rho c\_p \right)\_{eff2} v\_{m\text{s} 2}^\*} \tag{21}$$

where *v* ∗ *ms*<sup>1</sup> = R *<sup>H</sup>* 0 *v* ∗ *s*1 *dy*∗/*H* and *v* ∗ *ms*<sup>2</sup> = R 0 −*H v* ∗ *s*2 *dy*∗/*H* are the dimensional average velocities in layer I and layer II, respectively. Introducing the nondimensional temperature variable *θ<sup>i</sup>* = *kf*1(*T<sup>i</sup>* − *Tw*)/(*qwH*) and placing Equation (21) into Equations (18)–(20) yields

$$\frac{d^2\theta\_1}{dy^2} + \frac{k\_{f1}}{k\_{eff1}}(\mathcal{S} + Br\Phi\_1) = \frac{k\_{f1}}{k\_{eff1}} \frac{(2 + \mathcal{S} + Br\Gamma\_1 + m\_I Br\Gamma\_2)}{v\_{\text{ms1}} + (\rho c\_p)\_r v\_{\text{ms2}}} v\_{\text{s1}} \text{ for } 0 \le y \le 1 \tag{22}$$

$$\frac{d^2\theta\_2}{dy^2} + \frac{k\_{f1}}{k\_{eff2}} m\_r Br \Phi\_2 = \frac{k\_{f1}}{k\_{eff2}} \frac{(2 + S + Br\Gamma\_1 + m\_r Br\Gamma\_2)}{v\_{\text{ms1}}/(\rho c\_p)\_r + v\_{\text{ms2}}} v\_{\text{s2}} \text{ for} - 1 \le y \le 0 \tag{23}$$

$$\left. \theta\_1 \right|\_{y=1} = \theta\_2 \vert\_{y=-1} = 0 \; \text{\textquotedblleft}\_1 \vert\_{y=0} = \theta\_2 \vert\_{y=0'} \; k\_{\ell f1} \frac{d\theta\_1}{dy} \bigg|\_{y=0} = k\_{\ell f2} \frac{d\theta\_2}{dy} \bigg|\_{y=0} \tag{24}$$

where *vms*<sup>1</sup> = R 1 0 *vs*1*dy*, *vms*<sup>2</sup> = R 0 −1 *vs*2*dy*, Φ*<sup>i</sup>* = *dvsi*/*dy <sup>n</sup>i*−1(*dvsi*/*dy*) 2 , Γ<sup>1</sup> = R 1 <sup>0</sup> Φ1*dy*, Γ<sup>2</sup> = R 0 <sup>−</sup><sup>1</sup> <sup>Φ</sup>2*dy*, (*ρcp*)*r*= (*ρcp*)*eff*2/(*ρcp*)*eff*1, *<sup>m</sup><sup>r</sup>* <sup>=</sup> *<sup>m</sup>*2/*m*1, *Br* <sup>=</sup> *<sup>m</sup>*1*U*2/(*qwH*) is the Brinkman number, and *S* = *λHE<sup>z</sup>* <sup>2</sup>/*q<sup>w</sup>* is the Joule heating parameter.

With the two-layer temperature distributions solved from Equations (22)–(24), the heat transfer performance can be examined by evaluating heat transfer rate as represented by the Nusselt number, *Nu = hD<sup>h</sup> /keff.* With the convective heat transfer coefficient at the channel surface *h = qw*/(*T<sup>w</sup>* − *Tm*) and the characteristic height *D<sup>h</sup>* = 2*H* [55,57], the further rearrangement produces

$$Nu = -\frac{2k\_{f1}}{k\_{eff1}\theta\_m} \tag{25}$$

where *<sup>θ</sup><sup>m</sup>* = (R <sup>1</sup> 0 *vs*1*θ*1*dy* + R 0 −1 *vs*2*θ*2*dy*)/( R 1 0 *vs*1*dy* + R 0 −1 *vs*2*dy*) is the mean temperature.

#### *2.4. Entropy Generation Analysis*

Based on the second law of thermodynamics, a certain amount of energy is inevitably destroyed during the heat transfer process, that is, thermal irreversibility inherently accompanies the thermal behaviors and reduces system efficiency. This irreversibility is represented by the entropy generation rate, and the thermal performance of system is thereby assessed. The local entropy generation over a given cross-section of the microchannel can be given for two-layer flow as

$$S\_{l}^{\*}1(y^{\*}) = \frac{k\_{eff}}{T\_{1}^{2}} \left(\frac{dT\_{1}}{dy^{\*}}\right)^{2} + \frac{\sigma E\_{z}^{2}}{|T\_{1}|} + \frac{m\_{eff}}{|T\_{1}|} \left|\frac{dv\_{s1}^{\*}}{dy^{\*}}\right|^{n\_{i}-1} \left(\frac{dv\_{s1}^{\*}}{dy^{\*}}\right)^{2} \text{ for } 0 \le y \le 1 \tag{26}$$

$$S\_{l}^{\*}(y^{\*}) = \frac{k\_{eff}}{T\_{2}^{2}} \left(\frac{dT\_{2}}{dy^{\*}}\right)^{2} + \frac{m\_{eff2}}{|T\_{2}|} \left|\frac{dv\_{s2}^{\*}}{dy^{\*}}\right|^{n\_{2}-1} \left(\frac{dv\_{s2}^{\*}}{dy^{\*}}\right)^{2} \text{ for } -1 \le y \le 0 \tag{27}$$

in which the right-hand side of Equation (26) implies entropy generation caused by heat transfer, Joule heating and viscous dissipation, respectively. With the introduction of the nondimensional groups *Sli* = *H2S \* li*/*kf1* and Θ = *kf1Tw*/(*qwH*), the respective nondimensional local entropy generation is obtained as follows

$$S\_{l1}(y) = \frac{k\_{eff1}}{k\_{f1}(\theta\_1 + \Theta)^2} \left(\frac{d\theta\_1}{dy}\right)^2 + \frac{S}{|\theta\_1 + \Theta|} + \frac{Br}{|\theta\_1 + \Theta|} \left|\frac{dv\_{s1}}{dy}\right|^{n\_1 - 1} \left(\frac{dv\_{s1}}{dy}\right)^2 \text{ for } 0 \le y \le 1\tag{28}$$

$$S\_{f2}(y) = \frac{k\_{eff2}}{k\_{f1}(\theta\_2 + \Theta)^2} \left(\frac{d\theta\_2}{dy}\right)^2 + \frac{Br}{|\theta\_2 + \Theta|} \left|\frac{dv\_{s2}}{dy}\right|^{\theta\_2 - 1} \left(\frac{dv\_{s2}}{dy}\right)^2 \text{ for} - 1 \le y \le 0 \tag{29}$$

The total entropy generation can be computed by integrating Equations (28) and (29) over the relevant cross-section area of microchannel:

$$S\_l = \int\_0^1 S\_{l1} dy + \int\_{-1}^0 S\_{l2} dy \tag{30}$$

#### *2.5. Solutions of Modelling and Validation*

2.5.1. In the Case of Newtonian Fluids

The analytical velocities and analytical temperatures of a two-layer transient EOF of Newtonian fluids are first solved, and can then be employed to validate the numerical algorithm proposed for power-law nanofluid flow. Specifically, when *n<sup>i</sup>* = 1 and *φ* = 0, Equations (12)–(14) reduce to

$$\frac{\partial v\_1^N}{\partial t} = \frac{\partial^2 v\_1^N}{\partial y^2} - \frac{GE\_z \zeta}{\cosh(K)} \cosh(ky) \text{ for } 0 \le y \le 1 \tag{31}$$

$$\frac{\partial v\_2^N}{\partial t} = \frac{\partial^2 v\_2^N}{\partial y^2} \text{ for } -1 \le y \le 0 \tag{32}$$

$$v\_1^N \Big|\_{y=0} = v\_2^N \Big|\_{y=0} \cdot \frac{\partial v\_1^N}{\partial y} \Big|\_{y=0} = \frac{\partial v\_2^N}{\partial y} \Big|\_{y=0} \,, \ v\_1^N \Big|\_{y=1} = 0 \,, \ v\_1^N \Big|\_{y=-1} = 0 \,, \ v\_1^N \Big|\_{t=0} = 0 \tag{33}$$

where the superscript *N* denotes the special case of Newtonian fluid. Using the method of Laplace transformation, the analytical expressions of the transient velocities are obtained as follows:

$$\begin{split} v\_{1}^{N} &= \frac{GE\_{z}l}{K^{2}} \left[ \frac{1-\cosh(K)}{2\cosh(K)} (y-1) + \frac{\cosh(Ky)}{\cosh(K)} - 1 \right] \\ &+ 8GE\_{z}l \sum\_{P=1}^{\infty} \frac{(-1)^{P+l} - (2^{p-1})^{2} \pi^{2/4}}{(2P-1)\pi [4K^{2} + (2P-1)^{2}\pi^{2}]} \cos\left[\frac{(2P-1)\pi y}{2}\right] \\ &+ 6GE\_{z}l \sum\_{P=1}^{\infty} \frac{[1/\cosh(K) - (-1)^{P+l}]^{+} e^{-(P\pi)^{2}t}}{P\pi [K^{2} + (P\pi)^{2}]} \sin(P\pi y) \text{ for } 0 \le y \le 1 \end{split} \tag{34}$$
 
$$\begin{split} v\_{2}^{N} &= \frac{GE\_{z}l}{K^{2}} \left[ \frac{1-\cosh(K)}{2\cosh(K)} (y-1) + \frac{1}{\cosh(K)} - 1 \right] \\ &+ 8GE\_{z}l \sum\_{P=1}^{\infty} \frac{(-1)^{P+l} - (2^{p-1})^{2}\pi^{2/4} 0}{(2P-1)\pi [4K^{2} + (2P-1)^{2}\pi^{2}]} \cos\left[\frac{(2P-1)\pi y}{2}\right] \\ &+ 6GE\_{z}l \sum\_{P=1}^{\infty} \frac{[1/\cosh(K) - (-1)^{P+l}]e^{-(P\pi)^{2}t}}{Pr\{K^{2} + (P\pi^{2})\}} \sin(P\pi y) \text{ for } -1 \le y \le 0 \end{split}$$

for which the solving procedure is presented in the Appendix A in the interest of conciseness and readability.

As time tends to infinity, the two-fluid velocities reach a steady status as follows:

$$v\_{s1}^N(y) = \lim\_{t \to \infty} v\_1^N(y, t) = \frac{GE\_2\zeta}{K^2} \left[ \frac{1 - \cosh(K)}{2\cosh(K)}(y - 1) + \frac{\cosh(Ky)}{\cosh(K)} - 1 \right] \text{ for } 0 \le y \le 1 \tag{36}$$

$$v\_{s2}^{N}(y) = \lim\_{t \to \infty} v\_2^{N}(y, t) = \frac{GE\_2\zeta}{K^2} \left[ \frac{1 - \cosh(K)}{2\cosh(K)}(y - 1) + \frac{1}{\cosh(K)} - 1 \right] \text{ for } -1 \le y \le 0 \tag{37}$$

which are exactly the solutions to Equations (31)-(33) when the temporal term *∂v*i/*∂t* vanishes.

Replacing Equations (36) and (37) with Equations (22)–(24), the analytical temperature distributions are solved by integrating Equations (22) and (23) twice and combining them with Equation (24):

$$\theta\_1^N = A\_1 y^3 + A\_2 \cosh(\mathcal{K}y) + A\_3 y^2 + A\_4 \sinh(\mathcal{K}y) + A\_5 \left[ \frac{\cosh(\mathcal{K}y)^2}{\mathcal{K}^2} - y^2 \right] + D\_1 y + D\_2 \tag{38}$$

$$
\theta\_2^N = B\_1 y^3 + B\_2 y^2 + D\_3 y + D\_4 \tag{39}
$$

with the coefficients *A*, *B*, and *D* presented in the Appendix A for conciseness.

#### 2.5.2. In the Case of Power-Law Nanofluids

Due to the nonlinearity of a two-layer EOF of power-law nanofluids, Equations (12)– (14) and (22)–(24) are numerically solved based on the explicit finite difference method. At first, the following discretization is introduced: *t<sup>l</sup>* = *l*∆*t*, *y<sup>j</sup>* = *j*∆*y*, *v l <sup>i</sup>*,*<sup>j</sup>* = *vi*(*j*∆*y*, *l*∆*t*), and *θ l <sup>i</sup>*,*<sup>j</sup>* = *θi*(*j*∆*y*, *l*∆*t*), *l* = 1, 2, . . . , *L* and *j* = 1, 2, . . . , *J*.

At first, the numerical velocities at the channel walls are easily determined by discretizing the no slip conditions in Equation (14):

$$v\_{1,l}^l = v\_{2,1}^l = 0 \tag{40}$$

The numerical velocities at the two-liquid interface are computed using the bisection method from the discretized version of the two-liquid interface matching conditions in Equation (14)

$$v\_{1,1}^l = v\_{2,I}^l \, m\_1 \left(\frac{-3v\_{1,1}^l + 4v\_{1,2}^l - v\_{1,3}^l}{2\Delta y}\right)^{n\_1} = m\_2 \left(\frac{v\_{2,I-2}^l - 4v\_{2,I-1}^l + 3v\_{2,I}^l}{2\Delta y}\right)^{n\_2} \tag{41}$$

Note that for the proposed numerical algorithm, when *n<sup>i</sup>* ≥ 1, let the initial two-layer velocities be *v* 0 <sup>1</sup> = *v* 0 <sup>2</sup> = 0, and when *n<sup>i</sup>* < 1, to avoid the singularity caused by the zero denominator, let the initial two-layer velocities be *v* 0 <sup>1</sup> = *v* 0 <sup>2</sup> = *nonzero*.

Then, the velocities of the bulk liquid can be computed by means of the following numerical algorithm:

$$v\_{1,j}^{l+1} = v\_{1,j}^l + [\Lambda\_{1,j}^l + GE\_{\overline{z}}\zeta \cosh(Ky\_j)/\cosh(K)] \cdot \Delta t \tag{42}$$

$$v\_{2,j}^{l+1} = v\_{2,j}^l + \Lambda\_{2,j}^l \cdot \Delta t \tag{43}$$

where Λ *l <sup>i</sup>*,*<sup>j</sup>* = *mi m*<sup>0</sup> (*g l i*,*j* ) *<sup>n</sup>i*−<sup>1</sup> *<sup>v</sup> l <sup>i</sup>*,*j*+1−2*v l <sup>i</sup>*,*j*+*v l i*,*j*−1 ∆*y* <sup>2</sup> + (*n<sup>i</sup>* − 1)(*g l i*,*j* ) *<sup>n</sup>i*−<sup>2</sup> *<sup>g</sup> l <sup>i</sup>*,*j*+1−*g l i*,*j*−1 2∆*y v l <sup>i</sup>*,*j*+1−*v l i*,*j*−1 2∆*y* , *g l <sup>i</sup>*,*<sup>j</sup>* = *v l <sup>i</sup>*,*j*+1−*v l i*,*j*−1 2∆*y* , *i* = 1, 2, and *j* = 2, 3, . . . , *J* − 1. When *t* is great enough, the tran-

 sient velocities reach steady status, i.e., they are examined by k*v <sup>l</sup>* <sup>−</sup> *<sup>v</sup> <sup>l</sup>*+1<sup>k</sup> <sup>&</sup>lt; *Er*, with *Er* being a specified criterion, and the steady velocities *vs*<sup>1</sup> and *vs*<sup>2</sup> are solved numerically. Consequently, the flow rates can be computed from Equations (15) and (16) by means of the Simpson composite integration method [61].

The temporal term *∂θi*/*∂t* is introduced to Equations (22) and (23) to allow the same numerical algorithm to be used to solve both velocity distribution and temperature distribution. As time approaches to infinity, the fully thermally developed temperature distribution can be obtained. More specifically, at first the numerical temperatures at the boundaries and initial condition can be easily determined from

$$
\theta\_{1,l}^{l} = \theta\_{2,1}^{l} = 0,\\
\theta\_1^0 = \theta\_2^0 = 0 \tag{44}
$$

$$
\theta\_{1,1}^{l} = \theta\_{2,I'}^{l} \, k\_{\ell f1} \frac{-3\theta\_{1,1}^{l} + 4\theta\_{1,2}^{l} - \theta\_{1,3}^{l}}{2\Delta y} = k\_{\ell f2} \frac{\theta\_{2,I-2}^{l} - 4\theta\_{2,I-1}^{l} + 3\theta\_{2,I}^{l}}{2\Delta y} \tag{45}
$$

Then, the temperature distributions of the bulk liquid are computed based on the following numerical algorithm:

$$
\theta\_{1,j}^{l+1} = \theta\_{1,j}^{l+1} + (\Pi\_{1,j}^l + k\_{f1} \mathbb{S}/k\_{\varepsilon ff1}) \cdot \Delta t \tag{46}
$$

$$
\theta\_{2,j}^{l+1} = \theta\_{2,j}^l + (\Pi\_{2,j}^l) \cdot \Delta t \tag{47}
$$

2+*S*+*Br*Γ *l* *l*

where Π*<sup>l</sup>*

*<sup>i</sup>*,*<sup>j</sup>* = *<sup>i</sup>*,*j*+1−2*θ <sup>i</sup>*,*j*+*θ i*,*j*−1 ∆*y* <sup>2</sup> + *a<sup>i</sup> k <sup>f</sup>* <sup>1</sup> *ke f f <sup>i</sup>* Φ *l i*,*j* − *k <sup>f</sup>* <sup>1</sup> *ke f f <sup>i</sup>* <sup>1</sup>+*mrBr*Γ 2 *biv l ms*1+*civ l ms*2 *v l si* , Φ *l <sup>i</sup>*,*<sup>j</sup>* = (*g l i*,*j* ) *ni*−1 *v l si*,*j*+1−*v l si*,*j*−1 2∆*y* 2 , *a*<sup>1</sup> = *b*<sup>1</sup> = *c*<sup>2</sup> = 1, *a*<sup>2</sup> = *m<sup>r</sup>* , *b*<sup>2</sup> = 1/(*ρcp*) *r* , *c*<sup>1</sup> = (*ρcp*) *r* , *i* = 1, 2, *j* = 2, 3, . . . , *J –* 1, and the integrals Γ *l* <sup>1</sup> = R 1 <sup>0</sup> Φ *l* 1 *dy*, Γ *l* <sup>2</sup> = R 0 <sup>−</sup><sup>1</sup> <sup>Φ</sup> *l* 2 *dy*, *v l ms*<sup>1</sup> = R 1 0 *v l s*1 *dy*, *v l ms*<sup>2</sup> = R 0 −1 *v l s*2 *dy* are computed using the Simpson composite integration method. As time grows, if k*θ <sup>l</sup>* <sup>−</sup> *<sup>θ</sup> <sup>l</sup>*+1<sup>k</sup> <sup>&</sup>lt; *Er*, the thermally fully developed numerical temperature

*l*

*l*

 *θ l* distribution is obtained, based on which the Nusselt number in Equation (25) and total entropy generation in Equation (30) can be computed and the heat transfer analysis is conducted.

The physical parameters regarding the hydrodynamic properties of two power-law nanofluids are: *m*<sup>0</sup> \* = 9 <sup>×</sup> <sup>10</sup>−<sup>4</sup> Nm−<sup>2</sup> *s n* , *<sup>H</sup>* = 1 <sup>×</sup> <sup>10</sup>−<sup>4</sup> m, and *<sup>U</sup>* = 1 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>m</sup>·<sup>s</sup> −1 , which is of the same order as the Helmholtz–Smoluchowski velocity [52,53,61]. The physical parameters regarding the electric properties of two power-law nanofluids are : *<sup>ε</sup>*<sup>2</sup> = 7.08 <sup>×</sup> <sup>10</sup>−<sup>10</sup> F/m, *<sup>e</sup>* = 1.6 <sup>×</sup> <sup>10</sup>−<sup>19</sup> C, *<sup>z</sup>*<sup>0</sup> = 1, *<sup>E</sup><sup>z</sup>* \* = 1 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>V</sup>·m−<sup>1</sup> , *<sup>k</sup><sup>B</sup>* = 1.38 <sup>×</sup> <sup>10</sup>−<sup>23</sup> <sup>J</sup>·K−<sup>1</sup> , and *ζ* \* <sup>=</sup> <sup>−</sup>0.025 V [61]. To facilitate discussion, let the thermophysical properties of the fluids remain constant; the thermal conductivity of two base fluids are *kb*<sup>1</sup> = *kb*<sup>2</sup> = 0.618 Wm−1K−<sup>1</sup> , the thermal conductivity of AI2O<sup>3</sup> is *k<sup>p</sup>* = 40 Wm−1K−<sup>1</sup> , *T*<sup>0</sup> = 293 K, andthe Joule heating parameter *S* = 1 [61]. Furthermore, it is essential to provide the ranges of important nondimensional governing parameters based on practical physical uses. From the well-established electroosmotic theory of power-law modeling, *n<sup>i</sup>* = 0.6~1.4 [23] and the width of EDL is far less than the channel height; thus, *K* = 10~100 [23,52]. Based on the given order of reference (velocity, viscosity, and channel height), the Brinkman number ranges from *Br* = 0~0.1 [31,55], and the nanoparticle volume fraction *φ* = 0~0.09 is suitable for the chosen model of effective viscosity and effective thermal conductivity [44]. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 12 of 27 thermal conductivity of two base fluids are *kb*<sup>1</sup> = *kb*<sup>2</sup> = 0.618 Wm−1K−1, the thermal conductivity of AI2O3 is *kp* = 40 Wm−1K−1, *T*<sup>0</sup> = 293 K, andthe Joule heating parameter *S* = 1 [61]. Furthermore, it is essential to provide the ranges of important nondimensional governing parameters based on practical physical uses. From the well-established electroosmotic theory of power-law modeling, *ni* = 0.6~1.4 [23] and the width of EDL is far less than the channel height; thus, *K* = 10~100 [23,52]. Based on the given order of reference (velocity, viscosity, and channel height), the Brinkman number ranges from *Br* 

Concerning the numerical schemes, let *Er* = 1 <sup>×</sup> <sup>10</sup>−<sup>8</sup> in order to assure that the velocity and temperature reach steady status; then, a test grid dependence is conducted, with a grid system of 1 <sup>×</sup> <sup>10</sup><sup>3</sup> chosen. The good agreement between the analytical solutions and numerical solutions depicted in Figure 2 shows that the numerical algorithm proposed here is feasible for computing the two-layer velocity distribution and two-layer temperature distribution and for carrying out the heat transfer analysis. = 0~0.1 [31,55], and the nanoparticle volume fraction *ϕ* = 0~0.09 is suitable for the chosen model of effective viscosity and effective thermal conductivity [44]. Concerning the numerical schemes, let *Er* = 1 × 10−<sup>8</sup> in order to assure that the velocity and temperature reach steady status; then, a test grid dependence is conducted, with a grid system of 1 × 103 chosen. The good agreement between the analytical solutions and numerical solutions depicted in Figure 2 shows that the numerical algorithm proposed here is feasible for computing the two-layer velocity distribution and two-layer temperature distribution and for carrying out the heat transfer analysis.

**Figure 2.** Comparison of analytical solutions and numerical solutions for (**a**) velocity and (**b**) temperature. **Figure 2.** Comparison of analytical solutions and numerical solutions for (**a**) velocity and (**b**) temperature.

#### **3. Results and Discussion 3. Results and Discussion**

The transient hydrodynamic behavior of two-layer power-law nanofluid EOF is discussed by evaluating the transient two-layer velocities at different times and the flow rates for different governing nondimensional parameters. Then, with the steady velocities computed, the heat transfer analysis and entropy generation analysis are conducted by presenting the two-layer temperature profiles, Nusselt number and entropy generation rate at different governing nondimensional parameters. Specifically, *(ρCp*)*<sup>r</sup>* = 1 and *ρr* = 1, such that their importance can be attached to the effects of power-law rheology (represented by *ni*), the electroosmotic property (represented by electrokinetic width *K*), the nanoparticle volume fraction, *ϕ*, and the viscous dissipation effect (represented by the Brinkman number, *Br*) on hydrodynamic and thermal behaviors. The transient hydrodynamic behavior of two-layer power-law nanofluid EOF is discussed by evaluating the transient two-layer velocities at different times and the flow rates for different governing nondimensional parameters. Then, with the steady velocities computed, the heat transfer analysis and entropy generation analysis are conducted by presenting the two-layer temperature profiles, Nusselt number and entropy generation rate at different governing nondimensional parameters. Specifically, *(ρCp*)*<sup>r</sup>* = 1 and *ρ<sup>r</sup>* = 1, such that their importance can be attached to the effects of power-law rheology (represented by *ni* ), the electroosmotic property (represented by electrokinetic width *K*), the nanoparticle volume fraction, *φ*, and the viscous dissipation effect (represented by the Brinkman number, *Br*) on hydrodynamic and thermal behaviors.

#### The time evolutions of two-layer velocities with different types of conducting *3.1. Flow Characteristics in Two-Layer*

nanofluid, i.e., flow behavior index *n*<sup>1</sup> at different electrokinetic width *K*, are presented in Figure 3. It can be seen that at first the conducting nanofluid near the channel wall is set in motion due to the electroosmotic force, and as time elapses the bulk conducting nanofluid attains velocity. As opposed to the single-layer EOF, the flow of bulk The time evolutions of two-layer velocities with different types of conducting nanofluid, i.e., flow behavior index *n*<sup>1</sup> at different electrokinetic width *K*, are presented in Figure 3. It can be seen that at first the conducting nanofluid near the channel wall is set in motion

> conducting fluid fails to form a plug-like profile, as the delivery of momentum via the bulk conducting nanofluid is dissipated through interfacial viscous stress, causing a

*3.1. Flow Characteristics in Two-Layer*

due to the electroosmotic force, and as time elapses the bulk conducting nanofluid attains velocity. As opposed to the single-layer EOF, the flow of bulk conducting fluid fails to form a plug-like profile, as the delivery of momentum via the bulk conducting nanofluid is dissipated through interfacial viscous stress, causing a deviated parabolic profile in layer I. Accordingly, the nonconducting nanofluid is driven by interfacial shear stress; thus, in layer II, the closer to the two-liquid interface the flow, the higher the velocity. No matter what value *K* or *n*<sup>1</sup> takes, the conducting nanofluid near the wall at first moves forward and drags the nonconducting fluid through the interfacial hydrodynamical shear stress, and as *t* increases to 5, the transient two-layer velocity reaches its steady state. The comparison among Figure 3a–c describes the way in which, when the conducting nanofluid is shear thinning, the transient two-layer velocity is augmented with the increase of *K*, while as shown in Figure 3d–f, when the conducting nanofluid is shear thickening, the transient two-layer velocity shows abatement with the increase of *K*. The profiles of the steady velocities in Figure 3b,c show that the remarkable increase in the velocity of the conducting nanofluid with *K* results in a subtle turning of the velocity near the two-liquid interface. Moreover, when pumped by shear thinning nanofluid, the two-layer fluid is more sensitive to changes in the electrokinetic width *K* than when pumped by shear-thickening nanofluid. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 13 of 27 deviated parabolic profile in layer I. Accordingly, the nonconducting nanofluid is driven by interfacial shear stress; thus, in layer II, the closer to the two-liquid interface the flow, the higher the velocity. No matter what value *K* or *n*<sup>1</sup> takes, the conducting nanofluid near the wall at first moves forward and drags the nonconducting fluid through the interfacial hydrodynamical shear stress, and as *t* increases to 5, the transient two-layer velocity reaches its steady state. The comparison among Figure 3a–c describes the way in which, when the conducting nanofluid is shear thinning, the transient two-layer velocity is augmented with the increase of *K*, while as shown in Figure 3d–f, when the conducting nanofluid is shear thickening, the transient two-layer velocity shows abatement with the increase of *K*. The profiles of the steady velocities in Figure 3b,c show that the remarkable increase in the velocity of the conducting nanofluid with *K* results in a subtle turning of the velocity near the two-liquid interface. Moreover, when pumped by shear thinning nanofluid, the two-layer fluid is more sensitive to changes in the electrokinetic width *K* than when pumped by shear-thickening nanofluid.

**Figure 3.** Transient two-layer velocities at different electrokinetic widths *K* and different flow behavior indexes of conducting fluid *n*<sup>1</sup> when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *ϕ* = 0.03. (**a**) *K* = 10, *n*<sup>1</sup> = 0.8; (**b**) *K* = 50, *n*<sup>1</sup> = 0.8; (**c**) *K* = 100, *n*<sup>1</sup> = 0.8; (**d**) *K* = 10, *n*<sup>1</sup> = 1.2; (**e**) *K* = 50, *n*<sup>1</sup> = 1.2; (**f**) *K* = 100, *n*<sup>1</sup> = 1.2. **Figure 3.** Transient two-layer velocities at different electrokinetic widths *K* and different flow behavior indexes of conducting fluid *n*<sup>1</sup> when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *φ* = 0.03. (**a**) *K* = 10, *n*<sup>1</sup> = 0.8; (**b**) *K* = 50, *n*<sup>1</sup> = 0.8; (**c**) *K* = 100, *n*<sup>1</sup> = 0.8; (**d**) *K* = 10, *n*<sup>1</sup> = 1.2; (**e**) *K* = 50, *n*<sup>1</sup> = 1.2; (**f**) *K* = 100, *n*<sup>1</sup> = 1.2.

When the conducting nanofluid and nonconducting nanofluid change from shear thinning to shear thickening, the time evolutions of the two-layer velocities are as

When the conducting nanofluid and nonconducting nanofluid change from shear thinning to shear thickening, the time evolutions of the two-layer velocities are as delineated in Figure 4. No matter what type of pumped nonconducting nanofluid is considered, the decrease in the flow behavior index *n*1, namely, the shear thinning feature of conducting fluid, accelerates the flow and consequently enhances the two-layer velocity. From the comparison between Figures 4a–e, the change in nonconducting nanofluid type from shear thinning to shear thickening exerts a slight influence on the two-layer flow near the twoliquid interface, and shows little influence on that far away from the two-liquid interface. The influence of *n*<sup>1</sup> on two-layer flow far outweighs that of *n*<sup>2</sup> on two-layer flow. considered, the decrease in the flow behavior index *n*1, namely, the shear thinning feature of conducting fluid, accelerates the flow and consequently enhances the two-layer velocity. From the comparison between Figure 4a–c and Figure 4d–e, the change in nonconducting nanofluid type from shear thinning to shear thickening exerts a slight influence on the two-layer flow near the two-liquid interface, and shows little influence on that far away from the two-liquid interface. The influence of *n*<sup>1</sup> on two-layer flow far outweighs that of *n*<sup>2</sup> on two-layer flow.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 14 of 27

**Figure 4.** Transient two-layer velocities at different flow behavior indexes *n*<sup>1</sup> and different flow behavior indexes *n*<sup>2</sup> when *Br* = 0.02, *S* = 1, *ζ* = −1, *K* = 30 and *ϕ* = 0.03. (**a**) *n*<sup>1</sup> = 0.8, *n*<sup>2</sup> = 0.6; (**b**) *n*<sup>1</sup> = 1, *n*<sup>2</sup> = 0.6; (**c**) *n*<sup>1</sup> = 1.2, *n*<sup>2</sup> = 0.6; (**d**) *n*<sup>1</sup> = 0.8, *n*<sup>2</sup> = 1.4; (**e**) *n*<sup>1</sup> = 1, *n*<sup>2</sup> = 1.4; (**f**) *n*<sup>1</sup> = 1.2, *n*<sup>2</sup> = 1.4. **Figure 4.** Transient two-layer velocities at different flow behavior indexes *n*<sup>1</sup> and different flow behavior indexes *n*<sup>2</sup> when *Br* = 0.02, *S* = 1, *ζ* = −1, *K* = 30 and *φ* = 0.03. (**a**) *n*<sup>1</sup> = 0.8, *n*<sup>2</sup> = 0.6; (**b**) *n*<sup>1</sup> = 1, *n*<sup>2</sup> = 0.6; (**c**) *n*<sup>1</sup> = 1.2, *n*<sup>2</sup> = 0.6; (**d**) *n*<sup>1</sup> = 0.8, *n*<sup>2</sup> = 1.4; (**e**) *n*<sup>1</sup> = 1, *n*<sup>2</sup> = 1.4; (**f**) *n*<sup>1</sup> = 1.2, *n*<sup>2</sup> = 1.4.

Because the flow behavior index, *n*1, plays a crucial role in two-layer velocity, the time evolutions of flow rates at different flow behavior indexes of conducting nanofluid *n*<sup>1</sup> are presented in Figure 5 when (a) *ϕ* = 0 and (b) *ϕ* = 0.03. It is obvious that the unsteady flow rate of conducting nanofluid takes a longer time to reach steady status than that of nonconducting nanofluid; as the conducting nanofluid near the wall is set in motion first, the bulk conducting nanofluid moves forward, and eventually the nonconducting nanofluid is dragged by the bulk conducting fluid. Irrespective of the value of the nanoparticle volume fraction *ϕ*, the shear thinning feature of the conducting fluid enhances the flow rates of both the conducting and nonconducting nanofluids. A Because the flow behavior index, *n*1, plays a crucial role in two-layer velocity, the time evolutions of flow rates at different flow behavior indexes of conducting nanofluid *n*<sup>1</sup> are presented in Figure 5 when (a) *φ* = 0 and (b) *φ* = 0.03. It is obvious that the unsteady flow rate of conducting nanofluid takes a longer time to reach steady status than that of nonconducting nanofluid; as the conducting nanofluid near the wall is set in motion first, the bulk conducting nanofluid moves forward, and eventually the nonconducting nanofluid is dragged by the bulk conducting fluid. Irrespective of the value of the nanoparticle volume fraction *φ*, the shear thinning feature of the conducting fluid enhances the flow rates of both the conducting and nonconducting nanofluids. A comparison between Figure 5a,b

comparison between Figure 5a,b shows that the addition of the nanoparticle to the bulk

shows that the addition of the nanoparticle to the bulk nanofluid reduces the two-layer flow rate by improving the viscosity of the two fluids, as per Equations (9) and (10). nanofluid reduces the two-layer flow rate by improving the viscosity of the two fluids, as per Equations (9) and (10). nanofluid reduces the two-layer flow rate by improving the viscosity of the two fluids, as per Equations (9) and (10).

**Figure 5.** Time evolutions of flow rates of two power-law nanofluids at different flow behavior indexes *n*<sup>1</sup> and different volume fractions of nanoparticle *ϕ* when *Br* = 0.02, *S* = 1, *ζ* = −1, *K* = 10 and *n*<sup>2</sup> = 1.2. (**a**) *ϕ* = 0; (**b**) *ϕ* = 0.03 **Figure 5.** Time evolutions of flow rates of two power-law nanofluids at different flow behavior indexes *n*<sup>1</sup> and different volume fractions of nanoparticle *φ* when *Br* = 0.02, *S* = 1, *ζ* = −1, *K* = 10 and *n*<sup>2</sup> = 1.2. (**a**) *φ* = 0; (**b**) *φ* = 0.03. **Figure 5.** Time evolutions of flow rates of two power-law nanofluids at different flow behavior indexes *n*<sup>1</sup> and different volume fractions of nanoparticle *ϕ* when *Br* = 0.02, *S* = 1, *ζ* = −1, *K* = 10 and *n*<sup>2</sup> = 1.2. (**a**) *ϕ* = 0; (**b**) *ϕ* = 0.03

In Figure 6, the time evolutions of flow rates at different electrokinetic widths *K* are presented when (a) *n*<sup>1</sup> = 0.8 and (b) *n*<sup>2</sup> = 1.2. From Figure 6a, when the pumping conducting nanofluid is shear thinning, the flow rates of the two power-law nanofluids show augmentation with the electrokinetic width *K*, being consistent with Figure 3; in contrast, when it is shear thickening, the flow rates of both fluids show abatement with *K*, meaning that the change in fluid type of the pumping nanofluid from shear thinning to shear thickening can reverse the effect of the electrokinetic width. In addition, the flow rates show a noticeable reduction when *n*<sup>1</sup> increases, irrespective of what value *K* takes. In Figure 6, the time evolutions of flow rates at different electrokinetic widths *K* are presented when (a) *n*<sup>1</sup> = 0.8 and (b) *n*<sup>2</sup> = 1.2. From Figure 6a, when the pumping conducting nanofluid is shear thinning, the flow rates of the two power-law nanofluids show augmentation with the electrokinetic width *K*, being consistent with Figure 3; in contrast, when it is shear thickening, the flow rates of both fluids show abatement with *K*, meaning that the change in fluid type of the pumping nanofluid from shear thinning to shear thickening can reverse the effect of the electrokinetic width. In addition, the flow rates show a noticeable reduction when *n*<sup>1</sup> increases, irrespective of what value *K* takes. In Figure 6, the time evolutions of flow rates at different electrokinetic widths *K* are presented when (a) *n*<sup>1</sup> = 0.8 and (b) *n*<sup>2</sup> = 1.2. From Figure 6a, when the pumping conducting nanofluid is shear thinning, the flow rates of the two power-law nanofluids show augmentation with the electrokinetic width *K*, being consistent with Figure 3; in contrast, when it is shear thickening, the flow rates of both fluids show abatement with *K*, meaning that the change in fluid type of the pumping nanofluid from shear thinning to shear thickening can reverse the effect of the electrokinetic width. In addition, the flow rates show a noticeable reduction when *n*<sup>1</sup> increases, irrespective of what value *K* takes.

**Figure 6.** Time evolutions in flow rates of two power-law nanofluids at different flow behavior indexes *n*<sup>1</sup> and different electrokinetic widths *K* when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2 and *ϕ* = 0.03. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1.2. **Figure 6.** Time evolutions in flow rates of two power-law nanofluids at different flow behavior indexes *n*<sup>1</sup> and different electrokinetic widths *K* when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2 and *ϕ* = 0.03. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1.2. **Figure 6.** Time evolutions in flow rates of two power-law nanofluids at different flow behavior indexes *n*<sup>1</sup> and different electrokinetic widths *K* when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2 and *φ* = 0.03. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1.2.

In Figure 7, the variation in the ratio of flow rate I to flow rate II is plotted for a flow behavior index of conducting nanofluid *n*<sup>1</sup> at different nanoparticle volume fractions *ϕ*. When the pumping conducting nanofluid is shear thinning, that is, *n*<sup>1</sup> < 1, flow rate I is more than three times higher than flow rate II and two-layer flow is evidently affected by the change in the nanoparticle volume fraction, *ϕ*. In contrast, when the pumping In Figure 7, the variation in the ratio of flow rate I to flow rate II is plotted for a flow behavior index of conducting nanofluid *n*<sup>1</sup> at different nanoparticle volume fractions *ϕ*. When the pumping conducting nanofluid is shear thinning, that is, *n*<sup>1</sup> < 1, flow rate I is more than three times higher than flow rate II and two-layer flow is evidently affected by the change in the nanoparticle volume fraction, *ϕ*. In contrast, when the pumping conducting nanofluid is shear thickening (i.e., *n*<sup>1</sup> > 1), the two-layer flow is insensitive to In Figure 7, the variation in the ratio of flow rate I to flow rate II is plotted for a flow behavior index of conducting nanofluid *n*<sup>1</sup> at different nanoparticle volume fractions *φ*. When the pumping conducting nanofluid is shear thinning, that is, *n*<sup>1</sup> < 1, flow rate I is more than three times higher than flow rate II and two-layer flow is evidently affected by the change in the nanoparticle volume fraction, *φ*. In contrast, when the pumping

conducting nanofluid is shear thickening (i.e., *n*<sup>1</sup> > 1), the two-layer flow is insensitive to

flow.

conducting nanofluid is shear thickening (i.e., *n*<sup>1</sup> > 1), the two-layer flow is insensitive to the change in flow behavior index *n*<sup>1</sup> and nanoparticle volume fraction *φ*. Compared to the effect of the nanoparticle volume fraction, the effect of *n*<sup>1</sup> dominates in the two-layer flow. the change in flow behavior index *n*<sup>1</sup> and nanoparticle volume fraction *ϕ*. Compared to the effect of the nanoparticle volume fraction, the effect of *n*<sup>1</sup> dominates in the two-layer flow. the change in flow behavior index *n*<sup>1</sup> and nanoparticle volume fraction *ϕ*. Compared to the effect of the nanoparticle volume fraction, the effect of *n*<sup>1</sup> dominates in the two-layer

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 16 of 27

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 16 of 27

**Figure 7.** Variation of the ratio of flow rates with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *ϕ* when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *K* = 30. **Figure 7.** Variation of the ratio of flow rates with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *φ* when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *K* = 30. **Figure 7.** Variation of the ratio of flow rates with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *ϕ* when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *K* = 30.

In Figure 8, the variation in the ratio of flow rate I to flow rate II is presented for the flow behavior index of nonconducting nanofluid *n*<sup>2</sup> at different electrokinetic widths *K*. A higher flow behavior index of nonconducting nanofluid triggers stronger viscosity, and thus smaller flow rate of nonconducting nanofluid flow, finally leading to an increase in the flow rate ratio. When the type of conducting nanofluid is fixed, although a thinner EDL length improves velocities of both fluids (as shown in Figure 3a–c), the increase in flow rate I is much greater than that of flow rate II, and thus the increasing trend of flow rate ratio with *K* is observed, which is especially obvious for a shear thickening nonconducting nanofluid. Therefore, when dragging a shear thickening fluid, increasing electrokinetic width exerts more influence on flow rate I than on flow rate II. Furthermore, in comparison with Figure 7, the effect of the flow behavior index, *n*2, on the flow rate ratio is almost linear, and the increasing rate is enhanced for greater values of the electrokinetic width, *K*. In Figure 8, the variation in the ratio of flow rate I to flow rate II is presented for the flow behavior index of nonconducting nanofluid *n*<sup>2</sup> at different electrokinetic widths *K*. A higher flow behavior index of nonconducting nanofluid triggers stronger viscosity, and thus smaller flow rate of nonconducting nanofluid flow, finally leading to an increase in the flow rate ratio. When the type of conducting nanofluid is fixed, although a thinner EDL length improves velocities of both fluids (as shown in Figure 3a–c), the increase in flow rate I is much greater than that of flow rate II, and thus the increasing trend of flow rate ratio with *K* is observed, which is especially obvious for a shear thickening nonconducting nanofluid. Therefore, when dragging a shear thickening fluid, increasing electrokinetic width exerts more influence on flow rate I than on flow rate II. Furthermore, in comparison with Figure 7, the effect of the flow behavior index, *n*2, on the flow rate ratio is almost linear, and the increasing rate is enhanced for greater values of the electrokinetic width, *K*. In Figure 8, the variation in the ratio of flow rate I to flow rate II is presented for the flow behavior index of nonconducting nanofluid *n*<sup>2</sup> at different electrokinetic widths *K*. A higher flow behavior index of nonconducting nanofluid triggers stronger viscosity, and thus smaller flow rate of nonconducting nanofluid flow, finally leading to an increase in the flow rate ratio. When the type of conducting nanofluid is fixed, although a thinner EDL length improves velocities of both fluids (as shown in Figure 3a–c), the increase in flow rate I is much greater than that of flow rate II, and thus the increasing trend of flow rate ratio with *K* is observed, which is especially obvious for a shear thickening nonconducting nanofluid. Therefore, when dragging a shear thickening fluid, increasing electrokinetic width exerts more influence on flow rate I than on flow rate II. Furthermore, in comparison with Figure 7, the effect of the flow behavior index, *n*2, on the flow rate ratio is almost linear, and the increasing rate is enhanced for greater values of the electrokinetic width, *K*.

**Figure 8.** Variation of the ratio of volumetric flow rates with flow behavior index *n*<sup>2</sup> at different nanoparticle volume fractions when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>1</sup> = 0.8, and *φ* = 0.03.

#### *3.2. Heat Transfer Characteristics in Two-Layer Flow 3.2. Heat Transfer Characteristics in Two-Layer Flow* In Figure 9, the temperature profiles at different flow behavior indexes *n*1 are In Figure 9, the temperature profiles at different flow behavior indexes *n*1 are described when (a) *ϕ* = 0, (b) *ϕ* = 0.03, and (c) *ϕ* = 0.06. It is noted that the temperature

*3.2. Heat Transfer Characteristics in Two-Layer Flow*

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 17 of 27

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 17 of 27

nanoparticle volume fractions when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>1</sup> = 0.8, and *ϕ* = 0.03.

nanoparticle volume fractions when *Br* = 0.02, *S* = 1, *ζ* = −1, *n*<sup>1</sup> = 0.8, and *ϕ* = 0.03.

In Figure 9, the temperature profiles at different flow behavior indexes *n*<sup>1</sup> are described when (a) *φ* = 0, (b) *φ* = 0.03, and (c) *φ* = 0.06. It is noted that the temperature profiles of the two-layer flow are asymmetric around the two-liquid interface, in which the minimum value occurs within layer I. The temperature difference between the channel wall and bulk fluid is reduced, with the conducting nanofluid changing from a shear thinning fluid to a shear thickening one; therefore, it can be seen that the smaller the velocity gradient of two-layer flow is, the better the heat transfer performance becomes. On the other hand, the increase in the nanoparticle volume fraction, *φ*, leads to an overall increase in the temperature profiles for all types of conducting nanofluid, meaning that the introduction of nanoparticles truly improves the heat transfer in two-layer flow. described when (a) *ϕ* = 0, (b) *ϕ* = 0.03, and (c) *ϕ* = 0.06. It is noted that the temperature profiles of the two-layer flow are asymmetric around the two-liquid interface, in which the minimum value occurs within layer I. The temperature difference between the channel wall and bulk fluid is reduced, with the conducting nanofluid changing from a shear thinning fluid to a shear thickening one; therefore, it can be seen that the smaller the velocity gradient of two-layer flow is, the better the heat transfer performance becomes. On the other hand, the increase in the nanoparticle volume fraction, *ϕ*, leads to an overall increase in the temperature profiles for all types of conducting nanofluid, meaning that the introduction of nanoparticles truly improves the heat transfer in two-layer flow. profiles of the two-layer flow are asymmetric around the two-liquid interface, in which the minimum value occurs within layer I. The temperature difference between the channel wall and bulk fluid is reduced, with the conducting nanofluid changing from a shear thinning fluid to a shear thickening one; therefore, it can be seen that the smaller the velocity gradient of two-layer flow is, the better the heat transfer performance becomes. On the other hand, the increase in the nanoparticle volume fraction, *ϕ*, leads to an overall increase in the temperature profiles for all types of conducting nanofluid, meaning that the introduction of nanoparticles truly improves the heat transfer in two-layer flow.

**Figure 8.** Variation of the ratio of volumetric flow rates with flow behavior index *n*<sup>2</sup> at different

**Figure 8.** Variation of the ratio of volumetric flow rates with flow behavior index *n*<sup>2</sup> at different

**Figure 9.** Variation of temperature at different flow behavior indexes *n*<sup>1</sup> and different nanoparticle volume fractions *ϕ* when *K* = 30, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02. (**a**) *ϕ* = 0; (**b**) *ϕ* = 0.03; (**c**) *ϕ* = 0.06. **Figure 9.** Variation of temperature at different flow behavior indexes *n*<sup>1</sup> and different nanoparticle volume fractions *φ* when *K* = 30, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02. (**a**) *φ* = 0; (**b**) *φ* = 0.03; (**c**) *φ* = 0.06. volume fractions *ϕ* when *K* = 30, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02. (**a**) *ϕ* = 0; (**b**) *ϕ* = 0.03; (**c**) *ϕ* = 0.06.

As shown in Figure 10, the temperature profiles at different electrokinetic widths *K* are plotted when (a) *n*<sup>1</sup> = 0.8, (b) *n*<sup>1</sup> = 1, and (c) *n*<sup>1</sup> = 1.2 in order to demonstrate the influence of *n*<sup>1</sup> on the effect of *K*. The increase in *K*, namely a thinner EDL thickness, enlarges the temperature difference between the channel wall and bulk liquid and shifts the minimum value of the temperature to the right. This can be explained by the fact that the increase in *K* enhances the two-layer velocity and causes a drastic change in velocity near the channel wall, as described in Figure 3. The comparison among Figure 10a–c indicates that the shear thickening feature of the conducting nanofluid tends to suppress the effect of the electrokinetic width, *K*, on temperature profile. As shown in Figure 10, the temperature profiles at different electrokinetic widths *K* are plotted when (a) *n*<sup>1</sup> = 0.8, (b) *n*<sup>1</sup> = 1, and (c) *n*<sup>1</sup> = 1.2 in order to demonstrate the influence of *n*<sup>1</sup> on the effect of *K*. The increase in *K*, namely a thinner EDL thickness, enlarges the temperature difference between the channel wall and bulk liquid and shifts the minimum value of the temperature to the right. This can be explained by the fact that the increase in *K* enhances the two-layer velocity and causes a drastic change in velocity near the channel wall, as described in Figure 3. The comparison among Figure 10a–c indicates that the shear thickening feature of the conducting nanofluid tends to suppress the effect of the electrokinetic width, *K*, on temperature profile. As shown in Figure 10, the temperature profiles at different electrokinetic widths *K* are plotted when (a) *n*<sup>1</sup> = 0.8, (b) *n*<sup>1</sup> = 1, and (c) *n*<sup>1</sup> = 1.2 in order to demonstrate the influence of *n*<sup>1</sup> on the effect of *K*. The increase in *K*, namely a thinner EDL thickness, enlarges the temperature difference between the channel wall and bulk liquid and shifts the minimum value of the temperature to the right. This can be explained by the fact that the increase in *K* enhances the two-layer velocity and causes a drastic change in velocity near the channel wall, as described in Figure 3. The comparison among Figure 10a–c indicates that the shear thickening feature of the conducting nanofluid tends to suppress the effect of the electrokinetic width, *K*, on temperature profile.

**Figure 10.** Variation of temperature at different flow behavior indexes *n*<sup>1</sup> and different electrokinetic widths *K* when *ϕ* = 0.03, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1; (**c**) **Figure 10.** Variation of temperature at different flow behavior indexes *n*<sup>1</sup> and different electrokinetic widths *K* when *ϕ* = 0.03, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1; (**c**) *n*<sup>1</sup> = 1.2. **Figure 10.** Variation of temperature at different flow behavior indexes *n*<sup>1</sup> and different electrokinetic widths *K* when *φ* = 0.03, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1; (**c**) *n*<sup>1</sup> = 1.2.

As shown in Figure 11, the temperature profiles with different Brinkman numbers *Br* are plotted when (a) *n*<sup>1</sup> = 0.8, (b) *n*<sup>1</sup> = 1, and (c) *n*<sup>1</sup> = 1.2 in order to show the combined effect of flow behavior index *n*<sup>1</sup> and Brinkman number *Br*. As the Brinkman number *Br* increases, viscous dissipation in the two-layer flow is improved, which further enlarges

*n*<sup>1</sup> = 1.2.

the temperature difference between the channel wall and bulk liquid, meaning that the consideration of viscous dissipation evidently hinders the heat transfer in two-layer flow. Moreover, as obvious changes in temperature are observed in Figure 11a–c, the Brinkman number *Br* has an important effect on the temperature distribution for all types of conducting nanofluid; in other word, no matter what type of fluid the two-layer flow is driven by, the effect of *Br* on temperature cannot be neglected. the temperature difference between the channel wall and bulk liquid, meaning that the consideration of viscous dissipation evidently hinders the heat transfer in two-layer flow. Moreover, as obvious changes in temperature are observed in Figure 11a–c, the Brinkman number *Br* has an important effect on the temperature distribution for all types of conducting nanofluid; in other word, no matter what type of fluid the two-layer flow is driven by, the effect of *Br* on temperature cannot be neglected. *Br* are plotted when (a) *n*<sup>1</sup> = 0.8, (b) *n*<sup>1</sup> = 1, and (c) *n*<sup>1</sup> = 1.2 in order to show the combined effect of flow behavior index *n*<sup>1</sup> and Brinkman number *Br*. As the Brinkman number *Br*  increases, viscous dissipation in the two-layer flow is improved, which further enlarges the temperature difference between the channel wall and bulk liquid, meaning that the consideration of viscous dissipation evidently hinders the heat transfer in two-layer flow. Moreover, as obvious changes in temperature are observed in Figure 11a–c, the Brinkman number *Br* has an important effect on the temperature distribution for all

As shown in Figure 11, the temperature profiles with different Brinkman numbers

As shown in Figure 11, the temperature profiles with different Brinkman numbers *Br* are plotted when (a) *n*<sup>1</sup> = 0.8, (b) *n*<sup>1</sup> = 1, and (c) *n*<sup>1</sup> = 1.2 in order to show the combined effect of flow behavior index *n*<sup>1</sup> and Brinkman number *Br*. As the Brinkman number *Br*  increases, viscous dissipation in the two-layer flow is improved, which further enlarges

(**a**) (**b**) (**c**)

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 18 of 27

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 18 of 27

**Figure 11.** Variation of temperature at different flow behavior indicese *n*<sup>1</sup> and different Brinkman numbers *Br* when *ϕ* = 0.03, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1 and, *K* = 30. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1; (**c**) *n*<sup>1</sup> = 1.2. **Figure 11.** Variation of temperature at different flow behavior indicese *n*<sup>1</sup> and different Brinkman numbers *Br* when *φ* = 0.03, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1 and, *K* = 30. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1; (**c**) *n*<sup>1</sup> = 1.2. **Figure 11.** Variation of temperature at different flow behavior indicese *n*<sup>1</sup> and different Brinkman numbers *Br* when *ϕ* = 0.03, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1 and, *K* = 30. (**a**) *n*<sup>1</sup> = 0.8; (**b**) *n*<sup>1</sup> = 1; (**c**) *n*<sup>1</sup> = 1.2.

The temperature profiles at different flow behavior indices of nonconducting nanofluid *n*<sup>2</sup> are described in Figure 12. When the pumped nonconducting nanofluid changes from a shear thinning fluid to a shear thickening one, the temperature in the vicinity of the two-liquid interface slightly increases, while far away from the two-liquid interface it shows little change. This is because the increase in the flow behavior index, *n*2, causes a slight change in the two-layer velocity near the two-liquid interface. Therefore, compared to the influence of the fluid type of the pumping conducting nanofluid, that of the fluid type of the pumped nonconducting nanofluid (*n*2) on The temperature profiles at different flow behavior indices of nonconducting nanofluid *n*<sup>2</sup> are described in Figure 12. When the pumped nonconducting nanofluid changes from a shear thinning fluid to a shear thickening one, the temperature in the vicinity of the two-liquid interface slightly increases, while far away from the two-liquid interface it shows little change. This is because the increase in the flow behavior index, *n*2, causes a slight change in the two-layer velocity near the two-liquid interface. Therefore, compared to the influence of the fluid type of the pumping conducting nanofluid, that of the fluid type of the pumped nonconducting nanofluid (*n*2) on temperature profile is weak. The temperature profiles at different flow behavior indices of nonconducting nanofluid *n*<sup>2</sup> are described in Figure 12. When the pumped nonconducting nanofluid changes from a shear thinning fluid to a shear thickening one, the temperature in the vicinity of the two-liquid interface slightly increases, while far away from the two-liquid interface it shows little change. This is because the increase in the flow behavior index, *n*2, causes a slight change in the two-layer velocity near the two-liquid interface. Therefore, compared to the influence of the fluid type of the pumping conducting nanofluid, that of the fluid type of the pumped nonconducting nanofluid (*n*2) on temperature profile is weak.

temperature profile is weak.

**Figure 12.** Variation of temperature at different flow behavior indexes *n*<sup>2</sup> when *ϕ* = 0.03, *n*<sup>1</sup> = 0.8, *ζ*  **Figure 12.** Variation of temperature at different flow behavior indexes *n*<sup>2</sup> when *ϕ* = 0.03, *n*<sup>1</sup> = 0.8, *ζ*  = −1, *S* = 1, and *Br* = 0.02. **Figure 12.** Variation of temperature at different flow behavior indexes *n*<sup>2</sup> when *φ* = 0.03, *n*<sup>1</sup> = 0.8, *ζ* = −1, *S* = 1, and *Br* = 0.02.

= −1, *S* = 1, and *Br* = 0.02. The variation in the Nusselt number *Nu* with flow behavior index *n*<sup>2</sup> at different electrokinetic widths *K* is presented in Figure 13. The Nusselt number, *Nu*, shows a slightly ascending trend with flow behavior index *n*<sup>2</sup> no matter the value of the electrokinetic width. On the other hand, the increase in electrokinetic width, *K*, leads to the decrease in Nusselt number, *Nu.* This is because the rapid change in velocity profile caused by a lower EDL thickness triggers a widened temperature difference by suppressing the heat transfer performance of the two-layer flow. It should be noted that the effect of electroosmotic property of the fluid and the effect of the type of pumped nonconducting nanofluid do not interact with each other. caused by a lower EDL thickness triggers a widened temperature difference by suppressing the heat transfer performance of the two-layer flow. It should be noted that the effect of electroosmotic property of the fluid and the effect of the type of pumped nonconducting nanofluid do not interact with each other.

The variation in the Nusselt number *Nu* with flow behavior index *n*<sup>2</sup> at different electrokinetic widths *K* is presented in Figure 13. The Nusselt number, *Nu*, shows a slightly ascending trend with flow behavior index *n*<sup>2</sup> no matter the value of the electrokinetic width. On the other hand, the increase in electrokinetic width, *K*, leads to the decrease in Nusselt number, *Nu.* This is because the rapid change in velocity profile

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 19 of 27

**Figure 13.** Variation of Nusselt number with flow behavior index *n*<sup>2</sup> at different electrokinetic widths *K* when *ϕ* = 0.03, *n*<sup>1</sup> = 0.8, *ζ* = −1, *S* = 1, and *Br* = 0.02. **Figure 13.** Variation of Nusselt number with flow behavior index *n*<sup>2</sup> at different electrokinetic widths *K* when *φ* = 0.03, *n*<sup>1</sup> = 0.8, *ζ* = −1, *S* = 1, and *Br* = 0.02.

In Figure 14, the variation of the Nusselt number *Nu* with (a) electrokinetic width *K* and (b) Brinkman number *Br* is presented when choosing different types of conducting nanofluid. As shown in Figure 14a, a greater value of the electrokinetic width *K* reduces the heat transfer in two-layer flow, causing the dramatic change in velocity and widening temperature difference shown in Figures 3 and 10. Figure 14b shows that the Nusselt number *Nu* decreases with the Brinkman number *Br*, as the stronger viscous dissipation effect represented by greater values of *Br* enlarges the temperature difference between the channel wall and the bulk liquid, impeding heat transfer in two-layer flow. As opposed to the interaction between the effects of *K* and *n*2 predicted in Figure 13, the descending trend of the Nusselt number *Nu* with *K* or *Br* becomes smaller when the conducting nanofluid changes from a shear thinning fluid to a shear thickening one, meaning that when two-layer flow is driven by a shear thinning fluid the heat transfer performance is much more susceptible to changes in electrokinetic width, *K*, or the Brinkman number, *Br*. In Figure 14, the variation of the Nusselt number *Nu* with (a) electrokinetic width *K* and (b) Brinkman number *Br* is presented when choosing different types of conducting nanofluid. As shown in Figure 14a, a greater value of the electrokinetic width *K* reduces the heat transfer in two-layer flow, causing the dramatic change in velocity and widening temperature difference shown in Figures 3 and 10. Figure 14b shows that the Nusselt number *Nu* decreases with the Brinkman number *Br*, as the stronger viscous dissipation effect represented by greater values of *Br* enlarges the temperature difference between the channel wall and the bulk liquid, impeding heat transfer in two-layer flow. As opposed to the interaction between the effects of *K* and *n*<sup>2</sup> predicted in Figure 13, the descending trend of the Nusselt number *Nu* with *K* or *Br* becomes smaller when the conducting nanofluid changes from a shear thinning fluid to a shear thickening one, meaning that when twolayer flow is driven by a shear thinning fluid the heat transfer performance is much more susceptible to changes in electrokinetic width, *K*, or the Brinkman number, *Br*. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 20 of 27

**Figure 14.** Variation of the Nusselt number *Nu* with (a) electrokinetic width *K* and (b) Brinkman number *Br* at different flow behavior indexes *n*<sup>1</sup> when *ϕ* = 0.03, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *S* = 1. (**a**) *Br* = 0.02; (**b**) *K* = 30. **Figure 14.** Variation of the Nusselt number *Nu* with (a) electrokinetic width *K* and (b) Brinkman number *Br* at different flow behavior indexes *n*<sup>1</sup> when *φ* = 0.03, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *S* = 1. (**a**) *Br* = 0.02; (**b**) *K* = 30.

nanofluid changes from a shear thinning fluid to a shear thickening one, the Nusselt number *Nu* is augmented, hence the heat transfer performance of two-layer flow is enhanced. Furthermore, the ascending trend of *Nu* with *n*<sup>1</sup> becomes more slight. The increment in the nanoparticle volume fraction *ϕ* tends to improve the variation of *Nu* with *n*<sup>1</sup> as a whole, thus promoting the heat transfer of two-layer flow for all types of conducting nanofluid. This means that the addition of nanoparticles has little influence on the effect of the conducting nanofluid type on heat transfer (represented by *Nu*). According to Equation (25), growth in mean temperature with *ϕ* leads to the increase in *Nu*, implying that the temperature difference between the bulk fluid and channel walls is narrowed. Therefore, in practical engineering terms, when a wide temperature difference between the channel wall and the centerline occurs and might cause undesirable results (such as electrochemical decomposition or fluctuation of PH value in working liquids), the introduction of nanoparticles can intensify heat transfer

**Figure 15.** Variation of Nusselt number *Nu* with flow behavior index *n*<sup>1</sup> at different volume

fractions of nanoparticles *ϕ* when *K* = 30, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02.

Variation of the Nusselt number *Nu* with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *ϕ* is shown in Figure 15, revealing the interaction between

performance and help to avoid the problems mentioned above.

Variation of the Nusselt number *Nu* with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *φ* is shown in Figure 15, revealing the interaction between the effect of the conducting nanofluid and that of the nanoparticles. As the conducting nanofluid changes from a shear thinning fluid to a shear thickening one, the Nusselt number *Nu* is augmented, hence the heat transfer performance of two-layer flow is enhanced. Furthermore, the ascending trend of *Nu* with *n*<sup>1</sup> becomes more slight. The increment in the nanoparticle volume fraction *φ* tends to improve the variation of *Nu* with *n*<sup>1</sup> as a whole, thus promoting the heat transfer of two-layer flow for all types of conducting nanofluid. This means that the addition of nanoparticles has little influence on the effect of the conducting nanofluid type on heat transfer (represented by *Nu*). According to Equation (25), growth in mean temperature with *φ* leads to the increase in *Nu*, implying that the temperature difference between the bulk fluid and channel walls is narrowed. Therefore, in practical engineering terms, when a wide temperature difference between the channel wall and the centerline occurs and might cause undesirable results (such as electrochemical decomposition or fluctuation of PH value in working liquids), the introduction of nanoparticles can intensify heat transfer performance and help to avoid the problems mentioned above. nanoparticle volume fractions *ϕ* is shown in Figure 15, revealing the interaction between the effect of the conducting nanofluid and that of the nanoparticles. As the conducting nanofluid changes from a shear thinning fluid to a shear thickening one, the Nusselt number *Nu* is augmented, hence the heat transfer performance of two-layer flow is enhanced. Furthermore, the ascending trend of *Nu* with *n*<sup>1</sup> becomes more slight. The increment in the nanoparticle volume fraction *ϕ* tends to improve the variation of *Nu* with *n*<sup>1</sup> as a whole, thus promoting the heat transfer of two-layer flow for all types of conducting nanofluid. This means that the addition of nanoparticles has little influence on the effect of the conducting nanofluid type on heat transfer (represented by *Nu*). According to Equation (25), growth in mean temperature with *ϕ* leads to the increase in *Nu*, implying that the temperature difference between the bulk fluid and channel walls is narrowed. Therefore, in practical engineering terms, when a wide temperature difference between the channel wall and the centerline occurs and might cause undesirable results (such as electrochemical decomposition or fluctuation of PH value in working liquids), the introduction of nanoparticles can intensify heat transfer performance and help to avoid the problems mentioned above.

**Figure 14.** Variation of the Nusselt number *Nu* with (a) electrokinetic width *K* and (b) Brinkman number *Br* at different flow behavior indexes *n*<sup>1</sup> when *ϕ* = 0.03, *ζ* = −1, *n*<sup>2</sup> = 1.2, and *S* = 1. (**a**) *Br* =

Variation of the Nusselt number *Nu* with flow behavior index *n*<sup>1</sup> at different

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 20 of 27

**(a) (b)**

0.02; (**b**) *K* = 30.

**Figure 15.** Variation of Nusselt number *Nu* with flow behavior index *n*<sup>1</sup> at different volume fractions of nanoparticles *<sup>ϕ</sup>* when *K* <sup>=</sup> 30, *n*<sup>2</sup> = 1.2, *ζ* <sup>=</sup> <sup>−</sup>1, *S* <sup>=</sup> 1, and *Br* <sup>=</sup> 0.02. **Figure 15.** Variation of Nusselt number *Nu* with flow behavior index *<sup>n</sup>*<sup>1</sup> at different volume fractions of nanoparticles *φ* when *K* = 30, *n*<sup>2</sup> = 1.2, *ζ* = −1, *S* = 1, and *Br* = 0.02.

The variation of total entropy generation *S<sup>t</sup>* with (a) electrokinetic width *K* and (b) Brinkman number *Br* is presented in Figure 16 for different types of conducting nanofluids. An almost linear increase of total entropy generation *S<sup>t</sup>* with electrokinetic width *K* and Brinkman number *Br* can be observed in Figure 16. From Figure 16a, it can be inferred that a thinner EDL length (represented by a greater value of *K*) results in a wider velocity gap between the channel wall and the bulk fluid; therefore, the heat transfer is retarded, and entropy generation improves accordingly. Figure 16b implies that the viscous dissipation effect grows stronger with the Brinkman number, leading to an increase in entropy generation. In addition, the increasing trend with *K* and *Br* becomes smaller when nanofluid I changes from shear thinning to shear thickening.

**Figure 16.** Variation of total entropy generation with (**a**) electrokinetic width *K* and (**b**) Brinkman number *Br* at different flow behavior indexes *n*<sup>1</sup> when *ϕ* = 0.03, *n*<sup>2</sup> = 1.2, and *S* = 1. (**a**) *Br* = 0.02; (**b**) *K*  = 30. **Figure 16.** Variation of total entropy generation with (**a**) electrokinetic width *K* and (**b**) Brinkman number *Br* at different flow behavior indexes *n*<sup>1</sup> when *φ* = 0.03, *n*<sup>2</sup> = 1.2, and *S* = 1. (**a**) *Br* = 0.02; (**b**) *K* = 30.

The variation of total entropy generation *St* with (a) electrokinetic width *K* and (b) Brinkman number *Br* is presented in Figure 16 for different types of conducting nanofluids. An almost linear increase of total entropy generation *St* with electrokinetic width *K* and Brinkman number *Br* can be observed in Figure 16. From Figure 16a, it can be inferred that a thinner EDL length (represented by a greater value of *K*) results in a wider velocity gap between the channel wall and the bulk fluid; therefore, the heat transfer is retarded, and entropy generation improves accordingly. Figure 16b implies that the viscous dissipation effect grows stronger with the Brinkman number, leading to an increase in entropy generation. In addition, the increasing trend with *K* and *Br* becomes smaller when nanofluid I changes from shear thinning to shear thickening.

The variation of total entropy generation with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *ϕ* is presented in Figure 17. The entropy generation decreases when the conducting nanofluid changes from a shear thinning type to a shear thickening one. The shear thinning feature of conducting fluid accelerates the two-layer flow and enhances velocity distribution, while the increased velocity gradient and temperature gradient result in greater entropy generation as per Equations (28) and (29). Furthermore, entropy generation in two-layer flow driven by shear thinning nanofluid is more sensitive to changes in nanoparticle volume fraction compared to that driven by shear thickening nanofluid. In practical terms, when the two-layer flow is pumped by The variation of total entropy generation with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *φ* is presented in Figure 17. The entropy generation decreases when the conducting nanofluid changes from a shear thinning type to a shear thickening one. The shear thinning feature of conducting fluid accelerates the two-layer flow and enhances velocity distribution, while the increased velocity gradient and temperature gradient result in greater entropy generation as per Equations (28) and (29). Furthermore, entropy generation in two-layer flow driven by shear thinning nanofluid is more sensitive to changes in nanoparticle volume fraction compared to that driven by shear thickening nanofluid. In practical terms, when the two-layer flow is pumped by shear thinning fluid, it is more likely to utilize nanoparticles to adjust the thermal performance of the system. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 22 of 27

shear thinning fluid, it is more likely to utilize nanoparticles to adjust the thermal

**Figure 17.** Variation of total entropy generation with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *ϕ* when *K* = 30, *n*<sup>2</sup> = 1.2, *S* = 1, and *Br* = 0.02. **Figure 17.** Variation of total entropy generation with flow behavior index *n*<sup>1</sup> at different nanoparticle volume fractions *φ* when *K* = 30, *n*<sup>2</sup> = 1.2, *S* = 1, and *Br* = 0.02.

#### **4. Conclusions 4. Conclusions**

*K*.

weak.

generation.

	- When driven by shear thinning nanofluid, the two-layer flow accelerates for thinner EDL thicknesses and decelerates when driven by shear thickening • When driven by shear thinning nanofluid, the two-layer flow accelerates for thinner EDL thicknesses and decelerates when driven by shear thickening nanofluid.

nanofluid. The change in fluid type of pumped nonconducting nanofluid

nanofluid, the fluid type of the pumping conducting nanofluid plays a dominant role in two-layer flow and alters the effect of the electrokinetic width,

 In practical terms, the selection of a conducting nanofluid is crucial, as is the use of electrokinetic width to adjust two-layer flow for different types of

 As opposed to the variation of the flow rate ratio with *n*2, the variation with *n*<sup>1</sup> is nonlinear, and the flow rate of two-layer flow driven by shear thinning nanofluid is more sensitive to changes in the nanoparticle volume fraction. (2) With steady two-layer velocity obtained, the thermally developed heat transfer characteristics were discussed by presenting the temperature distribution, Nusselt

 The fluid type of the pumping conducting nanofluid, Brinkman number, nanoparticle volume fraction, and electrokinetic width all play important roles in the temperature profile, Nusselt number, and total entropy generation; in contrast, the influence of the type of pumped nonconducting nanofluid is

 In terms of the interactive influence of the governing parameters, shear thickening feature of the conducting nanofluid tends to suppress the effects of the Brinkman number and electrokinetic width on heat transfer and entropy

 No matter what type of conducting nanofluid is considered, increasing the nanoparticle volume fraction within a specified range truly enhances the heat

number, and total entropy generation at different parameters.

transfer performance of two-layer flow.

conducting nanofluid.

The change in fluid type of pumped nonconducting nanofluid exerts only a slight influence on velocity near the two-liquid interface. It is concluded that compared to the fluid type of pumped nonconducting nanofluid, the fluid type of the pumping conducting nanofluid plays a dominant role in two-layer flow and alters the effect of the electrokinetic width, *K*.

	- The fluid type of the pumping conducting nanofluid, Brinkman number, nanoparticle volume fraction, and electrokinetic width all play important roles in the temperature profile, Nusselt number, and total entropy generation; in contrast, the influence of the type of pumped nonconducting nanofluid is weak.
	- In terms of the interactive influence of the governing parameters, shear thickening feature of the conducting nanofluid tends to suppress the effects of the Brinkman number and electrokinetic width on heat transfer and entropy generation.
	- No matter what type of conducting nanofluid is considered, increasing the nanoparticle volume fraction within a specified range truly enhances the heat transfer performance of two-layer flow.
	- Entropy generation in two-layer flow driven by shear thinning nanofluid is more sensitive to changes in electrokinetic width, Brinkman number, and nanoparticle volume fraction.

**Author Contributions:** Conceptualization, S.D.; Data curation, S.D.; Formal analysis, S.D.; Funding acquisition, S.D. and T.X.; Investigation, S.D.; Methodology, S.D.; Project administration, S.D.; Resources, S.D.; Software, S.D.; Supervision, S.D.; Validation, S.D.; Visualization, S.D.; Writing—original draft, S.D.; Writing—review & editing, S.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Guangdong Basic and Applied Basic Research Foundation (Grant numbers 2021A1515012371 and 2020A1515011241), the National Natural Science Foundation of China (Grant number 11902082), and the Scientific Research Foundation of Universities in Guangdong Province for Young Talents (Grant number 2018KQNCX165). The authors are grateful to the Reviewers and Editors for their valuable comments and suggestions.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

To obtain the analytical solutions for two-layer EOF in the case of Newtonian fluid, Laplace transform to two-layer velocity distributions *v<sup>i</sup>* is introduced:

$$V\_i(s) = \mathcal{L}[v\_i(t)] = \int\_0^\infty v\_i^N(y\_\prime t) e^{-st} dt \tag{A1}$$

Transforming Equations (31)–(33), and letting *ρ<sup>r</sup>* = 1, the ordinary differential equations (ODEs) with respect to *s* can be obtained as follows:

$$\frac{d^2V\_1}{dy^2} - sV\_1 = \frac{1}{s}GE\_2\zeta \frac{\cosh(Ky)}{\cosh(K)} \text{ for } 0 \le y \le 1\tag{A2}$$

$$\frac{d^2V\_2}{d\overline{y}^2} - sV\_2 = 0 \text{ for } -1 \le y \le 0 \tag{A3}$$

$$\left.V\_1\right|\_{y=0} = \left.V\_2\right|\_{y=0^\nu} \left.\frac{dV\_1}{dy}\right|\_{y=0} = \left.\frac{dV\_2}{dy}\right|\_{y=0} \left.V\_2\right|\_{y=-1} = 0 \,\,\,\, V\_1|\_{y=1} = 0\tag{A4}$$

The solution to the above ODEs (A2)–(A3) takes the following forms:

$$V\_1 = \mathbb{C}\cosh(ky) + \mathbb{E}\cosh(\sqrt{s}y) + \mathbb{P}\sinh(\sqrt{s}y) \text{ for } 0 \le y \le 1\tag{A5}$$

$$V\_2 = M \cosh(\sqrt{s}y) + F \sinh(\sqrt{s}y) \text{ for } -1 \le y \le 0 \tag{A6}$$

Combining the boundary Equation (A4),

$$\begin{aligned} \mathcal{C} &= \frac{GE\_z \zeta}{s(K^2 - s)\cosh(K)}, E = \frac{-GE\_z \zeta}{2s(K^2 - s)\cosh\sqrt{s}} \left( 1 + \frac{\cosh\sqrt{s}}{\cosh(K)} \right) \\\\ F &= \frac{GE\_z \zeta}{2s(K^2 - s)\sinh\sqrt{s}} \left( \frac{\cosh\sqrt{s}}{\cosh(K)} - 1 \right), M = \frac{GE\_z \zeta}{s(K^2 - s)\cosh(K)} + E \end{aligned} \tag{A7}$$

Applying the Laplace inverse transform to Equations (A5) and (A6) yields

$$v\_i^N = \mathcal{L}^{-1}[V\_i] = \frac{1}{2\pi I} \int\_{\sigma - I\infty}^{\sigma + I\infty} V\_i e^{st} ds\tag{A8}$$

in which *I* is the imaginary unit, *σ* is a real number satisfying *σ* > 0, and the integral path is presented in Figure A1. Based on the residue theorem, Equation (A8) equals

$$v\_1^N = \sum\_k \text{Res}[V\_1(s)e^{st}, s\_k], \text{ for } 0 \le y \le 1\tag{A9}$$

$$\boldsymbol{v}\_{2}^{N} = \sum\_{k} \text{Res}[\boldsymbol{V}\_{2}(\boldsymbol{s})\boldsymbol{e}^{\boldsymbol{st}}, \boldsymbol{s}\_{k}]\_{\prime} \text{ for } -1 \le \boldsymbol{y} \le \boldsymbol{0} \tag{A10}$$

where Res[*V<sup>i</sup>* (s)*e st, s<sup>k</sup>* ] implies the residue of *V<sup>i</sup>* (*s*)*e st* at the isolated singularities, *s<sup>k</sup>* . The isolated singularities *s<sup>k</sup>* are enclosed in the simple closed curve *C<sup>R</sup>* + *L*, as shown in Figure A1. According to Equations (A5)–(A7), *s<sup>k</sup>* = 0, *K* 2 , –(2*P*–1)2π <sup>2</sup>/4, and –(*Pπ*) 2 , with *k* = 1, 2, 3, 4 and *P* = 1, 2, . . . . Consequently, the residues of *V<sup>i</sup>* (*s*)*e st* are evaluated as follows:

$$\operatorname{Res}[V\_1(s)e^{\xi t}, 0] = \lim\_{s\_1 \to 0} \frac{d[s^2 \cdot V\_1(s) \cdot e^{\xi t}]}{ds} = \frac{GE\_2\zeta}{K^2} \left[ \frac{1 - \cosh(K)}{2\cosh(K)}(y - 1) + \frac{\cosh(Ky)}{\cosh(K)} - 1 \right] \tag{A11}$$

$$\operatorname{Res}[V\_2(s)e^{st},0] = \lim\_{s\_1 \to 0} \frac{d[s^2 \cdot V\_2(s) \cdot e^{st}]}{ds} = \frac{GE\_2\zeta}{K^2} \left[ \frac{1 - \cosh(K)}{2\cosh(K)}(y-1) + \frac{1}{\cosh(K)} - 1 \right] \tag{A12}$$

$$\text{Res}[V\_l(s)e^{st}, K^2] = \lim\_{s\_2 \to K^2} (s - K^2) \cdot V\_l(s) \cdot e^{st} = 0 \tag{A13}$$

$$\text{Res}\left[V\_{l}(s)e^{st}, -\frac{\left[(2P-1)\pi\right]^{2}}{4}\right] = \lim\_{s\_{3P}\to -\frac{\left[(2P-1)\pi\right]^{2}}{4}} \left[s + \frac{\left[(2P-1)\pi\right]^{2}}{4}\right] \cdot V\_{l}(s) \cdot e^{st} \tag{A14}$$
 
$$\text{sCF:} \, \_{l}(-1)^{P+1} s^{-(2P-1)^{2}\pi^{2}t/4} \qquad \left[\text{(2P-1)\pi\nu}\right]$$

$$\begin{aligned} \gamma &= \frac{8GE\_5\zeta(-1)^{p+1}\epsilon^{-(2P-1)^2\pi^2t/4}}{(2P-1)\pi[4K^2+(2P-1)^2\pi^2]}\cos\left[\frac{(2P-1)\pi y}{2}\right] \\\\ \dots &\dots & \dots & \dots & \dots \end{aligned}$$

$$\begin{aligned} \text{Res}\left[V\_i(s)e^{st}, -\left(P\pi\right)^2\right] &= \lim\_{s\_{4p}\to -\left(P\pi\right)^2} \left[s + \left(P\pi\right)^2\right] \cdot V\_i(s) \cdot e^{st} \\\\ &= \frac{GE\_5[1/\cosh(K) - (-1)^{P+1}]e^{-\left(P\pi\right)^2 t}}{P\pi \left[K^2 + \left(P\pi^2\right)\right]} \sin(P\pi y) \end{aligned} \tag{A15}$$

with *i* = 1, 2. According to Equations (A9) and (A10), summarizing Equations (A11)–(A15) produces the analytical velocities for two-layer Newtonian fluid EOF, which are exactly expressed by Equations (34) and (35). with *i* = 1, 2. According to Equations (A9) and (A10), summarizing Equations (A11)– (A15) produces the analytical velocities for two-layer Newtonian fluid EOF, which are exactly expressed by Equations (34) and (35).

<sup>2</sup> 1 ()

π

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 24 of 27

*k*

*k*

in which *I* is the imaginary unit, *σ* is a real number satisfying *σ* > 0, and the integral path

where Res[*Vi*(s)*est, sk*] implies the residue of *Vi*(*s*)*est* at the isolated singularities, *sk*. The isolated singularities *sk* are enclosed in the simple closed curve *CR* + *L*, as shown in Figure A1. According to Equations (A5)–(A7), *sk* = 0, *K*2, –(2*P*–1)2π2/4, and –(*Pπ*)2, with *k* = 1, 2, 3,

*k*

*k*

*GEz <sup>K</sup> Ky <sup>y</sup> KK K*

<sup>−</sup> <sup>=</sup> − + <sup>−</sup>

*<sup>y</sup> KK K*

 π

π

π

*v V se s* ∑ , for 0 1 ≤ ≤ *<sup>y</sup>* (A9)

*v V se s* ∑ , for −≤ ≤ 1 0 *<sup>y</sup>* (A10)

1 cosh( ) cosh( ) ( 1) <sup>1</sup> 2cosh( ) cosh( )

1 cosh( ) <sup>1</sup> ( 1) <sup>1</sup> 2cosh( ) cosh( )

= − ⋅ ⋅= (A13)

<sup>−</sup> <sup>=</sup> − + <sup>−</sup> (A12)

(A11)

(A14)

(A15)

is presented in Figure A1. Based on the residue theorem, Equation (A8) equals 1 1 = Res[ ( ) , ] *<sup>N</sup> st*

2 2 = Res[ ( ) , ] *<sup>N</sup> st*

4 and *P* = 1, 2, …. Consequently, the residues of *Vi*(*s*)*est* are evaluated as follows:

*GEz K*

 π

ζ

ζ

2

*st st*

π

*P Pt GE e <sup>z</sup> P y*

 π

π

2 2

1

1

*s ds V s e V se* <sup>→</sup> *ds*

[ () ] Res[ ( ) , 0] lim

*s ds V s e V se* <sup>→</sup> *ds*

[ () ] Res[ ( ) , 0] lim

<sup>1</sup> <sup>0</sup>

<sup>2</sup> <sup>0</sup>

*st*

*st*

2

2

π

ζ

π

π

π

ζ 1

⋅ ⋅ <sup>=</sup> <sup>2</sup>

2

⋅ ⋅ <sup>=</sup> <sup>2</sup>

3

*P KP*

[(2 1) ] [(2 1) ] Res ( ) , lim ( ) 4 4 *<sup>P</sup>*

Res ( ) , ( ) lim ( ) ( ) *P*

2 2

*i i s P V se P s P Vs e* π

> π

+ − − − <sup>=</sup> <sup>+</sup>

→− − = + ⋅⋅

*i i <sup>P</sup> <sup>s</sup> P P V se <sup>s</sup> Vs e*

<sup>−</sup> →− − − − = + ⋅⋅

*st*

*st*

<sup>2</sup> <sup>2</sup> 2 2 Res[ ( ) , ] lim ( ) ( ) 0 *st st i i s K V se K s K V s e* <sup>→</sup>

> [(2 1) ] 4

+− − <sup>−</sup> <sup>−</sup> <sup>=</sup> − +−

2 2 ( )

[1/ cosh( ) ( 1) ] sin( ) [ ( )] *P Pt GE K e <sup>z</sup> P y PK P*

*st st*

2 2 1 (2 1) /4

2 2 2 8 ( 1) (2 1) cos (2 1) [4 (2 1) ] 2

<sup>2</sup> <sup>4</sup>

**Figure A1.** Schematic of the simple closed curve *CR* + *L* in the *s*-plane. **Figure A1.** Schematic of the simple closed curve *C<sup>R</sup>* + *L* in the *s*-plane.

In addition, the coefficients in Equations (38) and (39) are presented as follows:

*A*<sup>1</sup> = *k <sup>f</sup>* <sup>1</sup>*GE<sup>z</sup> ζ*[1−cosh(*K*)] 12*K*<sup>2</sup> cosh(*K*)*ke f f* <sup>1</sup> · (2+*S*+*Br*Γ1+*mrBr*Γ2) *vms*1+(*ρcp*) *r vms*<sup>2</sup> , *A*<sup>2</sup> = *k <sup>f</sup>* <sup>1</sup>*GE<sup>z</sup> ζ K*<sup>4</sup> cosh(*K*)*ke f f* <sup>1</sup> · (2+*S*+*Br*Γ1+*mrBr*Γ2) *vms*1+(*ρcp*) *r vms*<sup>2</sup> , *A*<sup>3</sup> = − *k <sup>f</sup>* <sup>1</sup> 2*ke f f* <sup>1</sup> *GEz ζ*[1+cosh(*K*)] 2*K*<sup>2</sup> cosh(*K*) · (2+*S*+*Br*Γ1+*mrBr*Γ2) *vms*1+(*ρcp*) *r vms*<sup>2</sup> + *S* + *Br*(*GEz ζ*) 2 [1−cosh(*K*)]<sup>2</sup> 4*K*<sup>4</sup> cosh (*K*) 2 *A*<sup>4</sup> = − *k <sup>f</sup>* <sup>1</sup>*Br*(*GE<sup>z</sup> ζ*) 2 [1−cosh(*K*)] *K*<sup>5</sup> cosh (*K*) 2 *ke f f* <sup>1</sup> , *A*<sup>5</sup> = − *k <sup>f</sup>* <sup>1</sup>*Br*(*GE<sup>z</sup> ζ*) 2 4*K*<sup>2</sup> cosh (*K*) 2 *ke f f* <sup>1</sup> , *B*<sup>1</sup> = *k <sup>f</sup>* <sup>1</sup>*GE<sup>z</sup> ζ*[1−cosh(*K*)] 12*K*<sup>2</sup> cosh(*K*)*ke f f* <sup>1</sup> · (2+*S*+*Br*Γ1+*mrBr*Γ2) *vms*1/(*ρcp*) *r* +*vms*<sup>2</sup> , *B*<sup>2</sup> = *k <sup>f</sup>* <sup>1</sup> 2*ke f f* <sup>1</sup> h *GE<sup>z</sup> ζ*[1−cosh(*K*)] 2*K*<sup>2</sup> cosh(*K*) · (2+*S*+*Br*Γ1+*mrBr*Γ2) *vms*1+*vms*2(*ρcp*) *r* − *mrBr*Φ<sup>2</sup> i , *D*<sup>1</sup> = − *ke f f* <sup>2</sup> (*d*1−*d*2) *ke f f* <sup>1</sup>+*ke f f* <sup>2</sup> , *D*<sup>2</sup> = − *ke f f* <sup>2</sup> *ke f f* <sup>1</sup>+*ke f f* <sup>2</sup> *ke f f* <sup>1</sup> *ke f f* <sup>2</sup> *d*<sup>1</sup> + *d*<sup>2</sup> , *D*<sup>3</sup> = *ke f f* <sup>1</sup> *ke f f* <sup>2</sup> (*KA*<sup>4</sup> + *D*1), *D*<sup>4</sup> = (*A*<sup>2</sup> + *A*5/*K* <sup>2</sup> + *D*2), *d*<sup>1</sup> = *A*<sup>1</sup> + *A*<sup>2</sup> cosh(*K*) + *A*<sup>3</sup> + *A*4sinh(*K*) + *A*<sup>5</sup> cosh (*K*) 2 *<sup>K</sup>*<sup>2</sup> − 1 *d*<sup>2</sup> = −*B*<sup>1</sup> + *B*<sup>2</sup> − *ke f f* <sup>1</sup> *ke f f* <sup>2</sup> *KA*<sup>4</sup> + *A*<sup>2</sup> + *A*5 *K*2 .

## **References**

