Thermally Fully Developed Electroosmotic Flow of Power-Law Nanofluid in a Rectangular Microchannel

The hydrodynamic and thermal behavior of the electroosmotic flow of power-law nanofluid is studied. A modified Cauchy momentum equation governing the hydrodynamic behavior of power-law nanofluid flow in a rectangular microchannel is firstly developed. To explore the thermal behavior of power-law nanofluid flow, the energy equation is developed, which is coupled to the velocity field. A numerical algorithm based on the Crank–Nicolson method and compact difference schemes is proposed, whereby the velocity, temperature, and Nusselt number are computed for different parameters. A larger nanoparticle volume fraction significantly reduces the velocity and enhances the temperature regardless of the base fluid rheology. The Nusselt number increases with the flow behavior index and with electrokinetic width when considering the surface heating effect, which decreases with the Joule heating parameter. The heat transfer rate of electroosmotic flow is enhanced for shear thickening nanofluids or at a greater nanoparticle volume fraction.


Introduction
The development of the micro-electric mechanical system (MEMS) has received great attention because of its innovative application in chemical, medical, and biological-related industries. For instance, Lab-on-a-chip, which conducts experiments by processing bioliquids or chemical solutions on a microchip, has become increasingly popular [1]. Electroosmotic flow (EOF) being a major electrokinetic effect refers to the motion of an ionized liquid inside microchannels with respect to the stationary charged channel walls under an external electric field applied tangentially along the microchannel [2]. An electroosmosis-based micropump has been successfully developed [3] and became one of the most important components in lab-on-a-chip since it has advantages, such as the production of pulse-free and plug-like flow, and the dependence on non-mechanical parts [4] over the conventional pressure-driven flow (PDF). Therefore, a fundamental understanding of the mechanism of EOF is vital for the precise control and optimal design of microfluidic devices.
The hydrodynamic behavior of EOF of different fluids is investigated in different microchannels, such as slit [5], rectangular [6,7], elliptic [8], and circular microchannels [9]. The core of the attention on EOF has shifted to the heat transfer characteristics of EOF [9][10][11], two-layer EOF [12], rotating EOF [13,14], and pressure effects on EOF [15]. In biological and chemical industries, biofluids, such as blood and DNA, solutions manipulated in microfluidic devices show nonlinear rheological behavior, such as the viscosity dependent on shear rate, which cannot be modeled by the linear constitutional relation. Therefore, the theory of non-Newtonian fluids offers researchers in the fields of mathematics and physics challenges in developing solutions to the nonlinear governing equations. The EOF of Jeffrey fluid [16] and Maxwell fluid [17] under an alternating current (AC) electric field are investigated in terms of the influence of nonlinear rheological behavior on flow performance. The corresponding fluid flow. Therefore, in microscale devices with characteristics of compact structure, to achieve the efficient thermal management of the non-Newtonian fluid, the power-law nanofluid was developed, where nanoparticles are added, and the base fluid rheology is described by the power-law model. The EOF of power-law nanofluid in a parallel plate microchannel was studied where one-dimensional momentum equation and energy equation were analytically solved [36]. A mixed convection flow and the corresponding heat transfer of shear thinning power-law nanofluid was respectively investigated by Ellahi et al. [37] and by Si et al. [38]. It turns out that the solid volume fractions exert special effects on both the velocity field and temperature field of fluid flow.
An up-to-date literature review indicates that for non-Newtonian fluids, the studies on the thermal behavior of EOF lack detailed insights into the substance of microstructure, namely the role that nanoparticles play, and the influence of channel geometry, owing to the complexity of governing equations for non-Newtonian nanofluid EOF in complex microchannels. Moreover, the thermal behavior of pure EOF of a power-law fluid in complex microchannels has not been studied yet. Due to the practical and fundamental significance of EOF of power-law nanofluid in a complex microchannel, such as a rectangular channel, and motivated by the works reviewed above, it is necessary to provide a fundamental understanding of the Joule heating induced heat transfer characteristics. This paper develops a model for thermally fully developed EOF of power-law nanofluid in a rectangular microchannel. The viscosity of power-law nanofluid depends both on the flow behavior index and the volume fraction of nanoparticles. By numerically solving the modified momentum equation and the energy equation, the velocity distribution, temperature distribution, and Nusselt number under the influence of the Joule heating effect are evaluated for different nanoparticle volume fraction, electrokinetic width, flow behavior index. The results are useful for the optimal design and active control of microfluidic systems.

Mathematical Modeling
As sketched in Figure 1, a rectangular microchannel of width 2b and height 2a is considered. The channel is filled with power-law nanofluid with a dielectric constant ε and is subject to constant heat flux q s . The channel wall is uniformly charged with zeta potential ξ, and the electric double layers (EDLs) will not overlap. As an electric field is applied tangentially along the channel (z-direction), the liquid inside is set in motion due to the existence of free ions in EDLs. Owing to the symmetry, the following analysis is limited to a quarter of the channel, Ω, and the coordinate system used is shown in Figure 1. developed, where nanoparticles are added, and the base fluid rheology is described by the powerlaw model. The EOF of power-law nanofluid in a parallel plate microchannel was studied where onedimensional momentum equation and energy equation were analytically solved [36]. A mixed convection flow and the corresponding heat transfer of shear thinning power-law nanofluid was respectively investigated by Ellahi et al. [37] and by Si et al. [38]. It turns out that the solid volume fractions exert special effects on both the velocity field and temperature field of fluid flow. An up-to-date literature review indicates that for non-Newtonian fluids, the studies on the thermal behavior of EOF lack detailed insights into the substance of microstructure, namely the role that nanoparticles play, and the influence of channel geometry, owing to the complexity of governing equations for non-Newtonian nanofluid EOF in complex microchannels. Moreover, the thermal behavior of pure EOF of a power-law fluid in complex microchannels has not been studied yet. Due to the practical and fundamental significance of EOF of power-law nanofluid in a complex microchannel, such as a rectangular channel, and motivated by the works reviewed above, it is necessary to provide a fundamental understanding of the Joule heating induced heat transfer characteristics. This paper develops a model for thermally fully developed EOF of power-law nanofluid in a rectangular microchannel. The viscosity of power-law nanofluid depends both on the flow behavior index and the volume fraction of nanoparticles. By numerically solving the modified momentum equation and the energy equation, the velocity distribution, temperature distribution, and Nusselt number under the influence of the Joule heating effect are evaluated for different nanoparticle volume fraction, electrokinetic width, flow behavior index. The results are useful for the optimal design and active control of microfluidic systems.

Mathematical Modeling
As sketched in Figure 1, a rectangular microchannel of width 2b and height 2a is considered. The channel is filled with power-law nanofluid with a dielectric constant ε and is subject to constant heat flux qs. The channel wall is uniformly charged with zeta potential ξ, and the electric double layers (EDLs) will not overlap. As an electric field is applied tangentially along the channel (z-direction), the liquid inside is set in motion due to the existence of free ions in EDLs. Owing to the symmetry, the following analysis is limited to a quarter of the channel,  , and the coordinate system used is shown in Figure 1.

Electric Potential Distribution
Based on the electrostatics theory, the relationship between the electric potential distribution ѱ(x, y) and the local net charge density ρe is expressed by the Poisson equation. (1)

Electric Potential Distribution
Based on the electrostatics theory, the relationship between the electric potential distribution ψ(x, y) and the local net charge density ρ e is expressed by the Poisson equation.
The local net charge density ρ e is subject to Boltzmann distribution.
where χ denotes the valence of ions, e denotes the elementary charge, n 0 is the ion density of the bulk liquid, k b is the Boltzmann constant, and T 0 implies the absolute temperature. Substituting Equation (2) to Equation (1) yields the two-dimensional Poisson-Boltzmann (P-B) equation for the electric potential distribution of EDL in the rectangular microchannel, The corresponding boundary conditions are given as To facilitate the discussion, the following scales are introduced x = x/a, y = y/a, ξ = χeξ/(k b T 0 ), α = b/a, K = κa, with K defined as the dimensionless electrokinetic width and 1/κ as the thickness of EDL, where κ = [2χe 2 n 0 /(εk b T 0 )] 1/2 . Accordingly, the dimensionless P-B equation and the boundary conditions are obtained ∂ψ/∂x x=0 = 0 , ∂ψ/∂y y=0 = 0 , ψ x=α = ξ , ψ y=1 = ξ.

Velocity Distribution
The vector formation of Cauchy momentum equation reads γ = (2e kl e kl ) 1/2 . µ eff denotes the effective viscosity coefficient. For a fully developed, incompressible and laminar EOF of power-law nanofluid, it is assumed that there is only an axial velocity w = w(x, y) [7,20], accordingly, the effective viscosity of the power-law nanofluid is µ eff = µ f / (1 − φ) 5/2 . The viscosity of the base fluid is µ f = η(|∂w/∂x| n − 1 , |∂w/∂y| n − 1 ) which is a function of strain rate, flow consistency index η of dimension (Nm −2 s n ) [7,20], and volume fraction of nanoparticle φ [39]. Finally, one can express the constitutive equation of power-law nanofluid as where τ is the shear stress, and n is the flow behavior index. n = 1 represents a Newtonian nanofluid, n < 1 and n > 1 represents the shear thinning nanofluids (pseudoplastic nanofluids) and the shear thickening nanofluids (dilatant nanofluids), respectively. Consequently, the modified Cauchy momentum equation for power-law nanofluid EOF is The corresponding boundary conditions are Introducing the following scales w = w/W with W = −εk b ET 0 /(µ 0 χe) where µ 0 is the viscosity coefficient of Newtonian fluid of dimension (Nm −2 s), one can obtain the dimensionless modified Cauchy momentum equation and the boundary conditions: ∂w/∂x x=0 = 0 , ∂w/∂y y=0 = 0 , w x=α = 0 , w y=1 = 0.
When considering Newtonian nanofluid, i.e., n = 1, and using the known Debye-Hückel principle, an analytical solution can be obtained by the method of variable separation and the method of constant variation [40], which is expressed as where β I = (2I − 1)π/(2α) with I = 1,2,3, . . . , In the case of non-Newtonian power-law nanofluid flow (n 1), a numerical algorithm is proposed to obtain the velocity field, which is presented in the next section.

Temperature Distribution
Assuming the thermophysical properties are constant, the heat generation induced by the viscous dissipation is much smaller than that induced by Joule heating [23,41] and thus, is neglected in the present analysis. The energy equation for the flow field in a rectangular microchannel is [32,42]. T denotes temperature field, ω implies the ratio of the nanolayer thickness to the original particle radius. σ, k, and (ρc p ) denote electrical conductivity, thermal conductivity, and heat capacity of power-law nanofluid at the reference pressure, respectively. The subscripts s, f, and eff denote the solid nanoparticles, base fluid, and nanofluid, respectively. The second term at the right-hand side of Equation (13) denotes the volumetric heat generation owing to the Joule heating effect. The corresponding boundary conditions are where q s = h(T w − T m ) implies the constant wall heat flux and h is the convective heat transfer coefficient. Subscripts w and m denote the local wall value and mean value, respectively. As assumed in the preceding section, for thermally fully developed EOF, one has ∂T/∂z = dT m /dz = dT w /dz = 0 and thus ∂ 2 T/∂z 2 = 0. It is noted that the temperature variation with respect to z-direction cannot be neglected and thus, the temperature on the wall is a function of z. Plus, an overall energy balance at an elemental control volume on a length of duct dz gives: Introducing T = k f (T − T w )/(q s a) and the Joule heating parameter S = σE 2 a/q s , namely, the ratio of Joule heating to the heat flux from the channel wall, the dimensionless energy equation and the corresponding boundary conditions are obtained as Micromachines 2019, 10, 363 In the case of Newtonian fluid and using the known Debye-Hückel linearization principle, an analytical solution of temperature distribution is obtained by the method of variable separation and the method of constant variation, where The detailed procedure is presented in Appendix A. For non-Newtonian power-law nanofluid, the energy equation can be solved using the numerical algorithm developed later.

Entropy Generation Analysis
With the temperature distribution obtained, Nusselt number indicating heat transfer rate from the channel wall to the liquid inside the channel is deduced as: where the mean temperature is and the convective heat transfer coefficient The entropy generation rate per unit volume is presented, which evaluates the amount of useful energy destroyed during the process, namely, the thermodynamic irreversibility of EOF [35,43,44].
Since the heat generated by viscous dissipation is negligible [23,41], the irreversibility of the system equals the summation of local volumetric entropy generation rate due to heat transfer and Joule heating. The corresponding nondimensional form is where Θ = T w k f /q w a [44]. The global total entropy generation rate can be obtained by integrating Equation (22) over the entire domain Ω. Then an important dimensionless number, Bejan number, representing the relative dominance of entropy generation due to heat transfer and Joule heating [45], is presented as It is noted that Be > 0.5 denotes that heat conduction irreversibility dominates and Be < 0.5 denotes the Joule heating irreversibility dominates.

Numerical Algorithm
Due to the nonlinearity of the governing equations above, the following numerical algorithm is developed to solve the velocity and temperature. Let t l = l∆t, x i = i∆x, y j = j∆y, where ∆t is the temporal step, ∆x and ∆y are the spatial steps with l = 1, 2, . . . , L, and i = 1, 2, . . . , I×J represents the matrix obtained by discretizing a variable f over the points above.
A numerical algorithm is proposed using the time splitting method in terms of time and the compact difference method in space. The complete P-B equation is numerically carried out based on the compact difference schemes, and the corresponding discretized form can be found in [20] and not presented here for conciseness.
To obtain the discretized velocity fieldw I×J by using the iterative method, one can introduce ∂/∂t on the left-hand side of Equations (11) and (16) and set a nonzero initial value for them, More specifically, each time step is split into two steps. In the first step, from t l to t * , one can focus on the third term on the right-hand side of Equations (19) and (20), and rewrite them as Accordingly, the velocity field at t = t * is given asw * Micromachines 2019, 10, 363 The subscript m implies the mean value of w l i,j . Setting d l = c 1 c l + c 2 and the discretized temperature at t = t l as T l , the Runge-Kutta method is applied to Equation (23): where In the second step, from t * to t l+1 , one can focus on the second-order derivatives in Equations (19) and (20), which are generally expressed as To note, as solving w l i,j , , the corresponding discretized forms have been given in Appendix B. When solving T l i,j , C 1 = 1 and C 2 = 1. According to the Crank-Nicolson technique, the discretized form of Equation (25) is Using the compact difference schemes and following the procedure presented in Appendix B, f i,j l+1 can be obtained as where and h x = ∆t ∆x 2 , h y = ∆t ∆y 2 , a 1 = 1 12 − . As a result, replacing f with w and T, the discretized solutions w l i,j and T l i,j can be computed. Finally, a specified criterion ∆L is given to identify that if the velocity is fully developed, i.e., w l − w l+1 < ∆L. Then, the fully developed velocity w i,j and temperature T i,j are obtained based on the iterative method above.

Results
A parametric study for velocity, temperature, and Nusselt number is conducted in the case of different flow behavior index n, electrokinetic width K, volume fraction of nanoparticles φ, and Joule heating parameter S. The nanoparticle is regarded as aluminum oxide [38]. The other physical parameters of the base fluid come from [9] and [20]. The typical values are presented in Table 1 below.

Parameters (Notation) Value (unit)
The fluid permittivity ε 7.08 × 10 − For validating the numerical algorithm proposed above, a test of grid dependence is conducted. As a result, the volumetric domain Ω is discretized to a grid system of 121 × 121. The numeric algorithm shows that the numerical solutions can achieve the accuracy of O(∆t 2 + ∆x 4 + ∆y 4 ), which is demonstrated in numerical computation by obtaining the numerical error of the order 10 −8 between the numerical velocity and analytical velocity of Newtonian nanofluid EOF. Furthermore, the iterative criterion is given as ∆L < 10 −10 to make sure the EOF is fully developed. For the EOF of Newtonian nanofluid, in Figure 2a, the numerical velocity profile at y = 0 is compared with the analytical velocity obtained from Equation (13)  Valence of ions χ 1 Electrokinetic width K 5-50 Flow behavior index n 0.7-1.4 Nanoparticle volume fraction ϕ 0.0-0.06 Joule heating Parameter S −5-5 Ratio of nanolayer thickness to original particle radius ω 1.1 For validating the numerical algorithm proposed above, a test of grid dependence is conducted. As a result, the volumetric domain Ω is discretized to a grid system of 121 × 121. The numeric algorithm shows that the numerical solutions can achieve the accuracy of which is demonstrated in numerical computation by obtaining the numerical error of the order 10 −8 between the numerical velocity and analytical velocity of Newtonian nanofluid EOF. Furthermore, the iterative criterion is given as ΔL < 10 −10 to make sure the EOF is fully developed. For the EOF of Newtonian nanofluid, in Figure 2a Figure 3b,d,f show the velocity distributions of shear thinning nanofluid, Newtonian nanofluid, and shear thickening nanofluid at φ = 0.06. It is noted from these two groups of figures that the velocity shows homogeneous abatement across the channel as the volume fraction of nanoparticles increases no matter what fluid rheology is considered. This is because the nanoparticles in fluid enhance the effective viscosity of power-law nanofluid in response to the sear rate and cause greater dispersion of the velocity distribution. As observed in [7], although considering the EOF of power-law nanofluid, the shear thinning feature of base fluid leads to greater velocity compared to the shear thickening feature of the base fluid.   Without consideration of nanoparticles in fluid (φ = 0) and when K = 10, the temperature distributions of shear thinning fluid, Newtonian fluid, and shear thickening fluid at S = 0 are displayed in Figure 4a,c,e. For the same prescribed conditions, the temperature distributions at S = 3 are displayed in Figure 4b,d,f. The comparison between these two groups of figures shows that the variation of temperature along the channel becomes greater with the increase of the Joule heating parameter S irrespective of the fluid type. Furthermore, the dimensionless temperature increases with the flow behavior index n, indicating that the difference of heat transfer characters caused by different shear rates of fluid type is remarkable, especially for shear thinning base fluid.
Micromachines 2019, 10, x 12 of 22 0 y  is illustrated in Figure 6 for different volume fractions of the nanoparticles. From Figure 6 and the comparison between Figures 4 and 5, the dimensionless temperature distribution shows an increment when the volume fraction of nanoparticles ϕ increases, or the flow behavior index n increases. It means that the combined effects of larger thermal conductivity of nanoparticles and the change in shear stress leads to the increase in thermal diffusion effect and in further the temperature. To observe the influence of volume fraction of nanoparticles on the heat transfer characteristics, the temperature distributions are evaluated in Figure 5 when the nanoparticle volume fraction φ is increased to 0.06, and other parameters keep unchanged. The corresponding temperature profiles at y = 0 is illustrated in Figure 6 for different volume fractions of the nanoparticles. From Figure 6 and the comparison between Figures 4 and 5, the dimensionless temperature distribution shows an increment when the volume fraction of nanoparticles φ increases, or the flow behavior index n increases. It means that the combined effects of larger thermal conductivity of nanoparticles and the change in shear stress leads to the increase in thermal diffusion effect and in further the temperature. To observe the influence of the Joule heating parameter on the temperature, in Figure 7, the temperature profiles at 0 y = for different Joule heating parameters S are illustrated in the case of shear thinning nanofluid, Newtonian nanofluid, and shear thickening nanofluid when ϕ = 0.06 and K = 10. As seen in Figure 5, no matter what nanofluid is considered, as the Joule heating parameter increases the variation of temperature increases, namely, the temperature difference between the centerline and the channel wall enlarges. It implies that although the energy from Joule heating effect is uniformly enhanced, the heat transferred through convection near the channel wall is much less than near the centerline. And this effect is more notable when considering shear thinning base fluid. To observe the influence of the Joule heating parameter on the temperature, in Figure 7, the temperature profiles at y = 0 for different Joule heating parameters S are illustrated in the case of shear thinning nanofluid, Newtonian nanofluid, and shear thickening nanofluid when φ = 0.06 and K = 10. As seen in Figure 5, no matter what nanofluid is considered, as the Joule heating parameter increases the variation of temperature increases, namely, the temperature difference between the centerline and the channel wall enlarges. It implies that although the energy from Joule heating effect is uniformly enhanced, the heat transferred through convection near the channel wall is much less than near the centerline. And this effect is more notable when considering shear thinning base fluid.
The variations of the Nusselt number with the electrokinetic width are shown in Figure 8 for different fluid types and Joule heating parameters when φ = 0.06. The cases S < 0 and S > 0 correspond to the surface cooling effect and surface heating effect, respectively. Importantly, in the case of S < 0, the Nusselt number decreases as electrokinetic width K increases, which increases when S ≥ 0 and other parameters remain unchanged. Surface heating (cooling) has a uniform influence on the temperature distribution, which leads to the increasing (decreasing) tendency of Nusselt number with electrokinetic width. Furthermore, the decrease in flow behavior index results in the pronounced abatement of Nusselt number, meaning that the heat transfer rate of shear thinning nanofluid is much smaller than the heat transfer rate of Newtonian fluid and shear thickening fluid. This is consistent with the prediction in [22].
increases the variation of temperature increases, namely, the temperature difference between the centerline and the channel wall enlarges. It implies that although the energy from Joule heating effect is uniformly enhanced, the heat transferred through convection near the channel wall is much less than near the centerline. And this effect is more notable when considering shear thinning base fluid. The variations of the Nusselt number with the electrokinetic width are shown in Figure 8 for different fluid types and Joule heating parameters when ϕ = 0.06. The cases S < 0 and S > 0 correspond to the surface cooling effect and surface heating effect, respectively. Importantly, in the case of S < 0, the Nusselt number decreases as electrokinetic width K increases, which increases when S ≥ 0 and other parameters remain unchanged. Surface heating (cooling) has a uniform influence on the temperature distribution, which leads to the increasing (decreasing) tendency of Nusselt number with electrokinetic width. Furthermore, the decrease in flow behavior index results in the pronounced abatement of Nusselt number, meaning that the heat transfer rate of shear thinning nanofluid is much smaller than the heat transfer rate of Newtonian fluid and shear thickening fluid. This is consistent with the prediction in [22].   Figure 9 depicts the variations of the Nusselt number with the Joule heating parameter for different fluid types and different nanoparticle volume fraction when K = 10. It shows that the Nusselt number is inversely proportional to the Joule heating parameter, and the decreasing rate becomes smaller as nanofluid changes from shear thickening to shear thinning. Since it is considered that Joule heating has a spatially uniform effect, the greater Joule heating parameter corresponds to relatively uniform heat distribution. Consequently, the convective heat transfer rate is reduced as presented in Figure 7, and one can observe the decreasing trend of the Nusselt number, namely, heat transfer rate. Being consistent with the prediction in Figure 7, the smaller convective heat transfer rate of shear thinning nanofluid leads to a smaller Nusselt number. Furthermore, the comparison between Figure 9a  The variation of the Nusselt number Nu and the variation of Bejan number Be with flow behavior index n for different volume fractions of nanoparticles are depicted in Figure 10 when K = 20 and S = 3. The combined effects of nonlinear rheological behavior and nanoparticle volume fraction on Nusselt number are investigated. To note, the dependence of Nusselt number on the flow behavior index shows nonlinear. For shear thinning base fluid, the Nusselt number is more amenable to the change in flow behavior index n. The increase in viscosity for shear thinning fluid resulting from the decrease in flow behavior index reduces the Nusselt number, and thus, the heat transfer rate of shear thinning nanofluid is obviously less than that of shear thickening nanofluid. Furthermore, the variation of the Nusselt number with the volume fraction of nanoparticles is almost linear. It shows that the heat transfer performance of nanofluid is clearly enhanced compared to the conventional fluid regardless of the base fluid rheology. These results are also witnessed in Figure 10b. The increase in nanoparticle volume fraction and flow behavior index leads to the enhancement of heat transfer, which further reduces the "loss" of useful energy, i.e., the total irreversibility of the system. Moreover, the Be > 0.5 implies that the primary irreversibility is contributed to the heat conduction when S = 3. The variation of the Nusselt number Nu and the variation of Bejan number Be with flow behavior index n for different volume fractions of nanoparticles are depicted in Figure 10 when K = 20 and S = 3. The combined effects of nonlinear rheological behavior and nanoparticle volume fraction on Nusselt number are investigated. To note, the dependence of Nusselt number on the flow behavior index shows nonlinear. For shear thinning base fluid, the Nusselt number is more amenable to the change in flow behavior index n. The increase in viscosity for shear thinning fluid resulting from the decrease in flow behavior index reduces the Nusselt number, and thus, the heat transfer rate of shear thinning nanofluid is obviously less than that of shear thickening nanofluid. Furthermore, the variation of the Nusselt number with the volume fraction of nanoparticles is almost linear. It shows that the heat transfer performance of nanofluid is clearly enhanced compared to the conventional fluid regardless of the base fluid rheology. These results are also witnessed in Figure 10b. The increase in nanoparticle volume fraction and flow behavior index leads to the enhancement of heat transfer, which further reduces the "loss" of useful energy, i.e., the total irreversibility of the system. Moreover, the Be > 0.5 implies that the primary irreversibility is contributed to the heat conduction when S = 3. thinning nanofluid is obviously less than that of shear thickening nanofluid. Furthermore, the variation of the Nusselt number with the volume fraction of nanoparticles is almost linear. It shows that the heat transfer performance of nanofluid is clearly enhanced compared to the conventional fluid regardless of the base fluid rheology. These results are also witnessed in Figure 10b. The increase in nanoparticle volume fraction and flow behavior index leads to the enhancement of heat transfer, which further reduces the "loss" of useful energy, i.e., the total irreversibility of the system. Moreover, the Be > 0.5 implies that the primary irreversibility is contributed to the heat conduction when S = 3.

Conclusions
In the present paper, a numerical investigation for hydrodynamically and thermally fully developed EOF of power-law nanofluid in a rectangular microchannel has been conducted. The interactive influences of electrokinetic width K, flow behavior index n, nanoparticle volume fraction

Conclusions
In the present paper, a numerical investigation for hydrodynamically and thermally fully developed EOF of power-law nanofluid in a rectangular microchannel has been conducted. The interactive influences of electrokinetic width K, flow behavior index n, nanoparticle volume fraction φ, and Joule heating parameter S on the velocity field, temperature field, Nusselt number, and Bejan number have been studied. The following conclusions can be drawn:

•
The velocity is reduced due to the increase in viscosity of nanofluid resulting from the addition of nanoparticles.

•
The volume fraction of nanoparticles φ truly enhances the dimensionless temperature and Nusselt number by changing the thermal conductivity and shear stress of velocity, no matter what type of base fluid is considered. Accordingly, the Bejan number decreases with increasing nanoparticle volume fraction and flow behavior index.

•
Nusselt number is an increasing function of the electrokinetic width K for S ≥ 0 and a decreasing function for S < 0.

•
The Nusselt number shows abatement, and the temperature exhibits greater variation along the channel as Joule heating parameter S increases. The heat transfer rate of shear thickening nanofluid is much larger than that of shear thinning fluid.
To sum up, nanoparticle volume fraction and Joule heating effect show significant influence on the thermal behavior of EOF of power-law nanofluid. It should be noted that the dependence of the Nusselt number on the flow behavior index is nonlinear. The Nusselt number becomes a weak function of electrokinetic width as the electrokinetic width gets larger. The heat transfer rate of shear thinning nanofluid is more amenable to the change in nanofluid rheology. Therefore, the heat transfer of power-law nanofluid has to be carefully treated, and one should make a distinction when manipulating different types of nanofluid in practical problems. This understanding of the mechanism of thermally fully developed EOF can be referred to in more efficient control of fluids and optimization of microfluidic devices.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A
The velocity distribution of the Newtonian nanofluid is solved following the procedure in [40] and is presented as Equation (13). As the nanoparticle volume fraction vanishes, the velocity of Equation (13) becomes Equation (14) in [21].
The temperature distribution of Newtonian nanofluid EOF is solved as follows: in the first step, based on the method of variable separation [40], the temperature can be expressed as with β I = (2I − 1)π/(2α), I = 1,2,3, . . . . In the second step, Y I (y) needs to be determined. Substituting Equation (A1) into Equation (17) yields χ I (y) can be expressed from Equation (A1a) using the orthogonality of Fourier series, namely, where k 1 = In the third step, with the expression of χ I (y), the ODE of Equation (A1b) is solved adopting the method of constant variation, the solution form is given as Y IG (y) = H 1 e C I y + H 2 e −C I y , Consequently, the analytical solution of temperature for Newtonian nanofluid flow is carried out.

Appendix B
In the case of f = w, the compact difference schemes for computing C 1 and C 2 are given as 1 6 F y,i,j+1 (w) + 2 3 F y,i,j (w) + 1 6 F y,i,j−1 (w) = (w i,j+1 − w i, j−1 )/2∆y, in which F x,i,j (w) and F y,i,j (w) are the difference approximations of ∂w/∂x and ∂w/∂y, respectively.