Numerical Analysis and Entropy Generation of Electrokinetic Flow of Power-Law Fluid in a Microtriangular Prism

Abstract

This research aims to study the characteristics of thermal transport and analyse the entropy generation of electroosmotic flow of power-law fluids in a microtriangular prism in the presence of pressure gradient. Considering a fully developed flow subject to constant wall heat flux, the nonlinear electric potential, momentum, and linear heat transfer equations are solved numerically by developing an iterative finite difference method with a nonuniform grid. The thermal efficiency of the model is explored under the light of the second law of thermodynamics. Effect/impact of governing physical parameters on velocity, temperature, Nusselt number, and entropy distributions is studied, and the results are demonstrated graphically; we found that the Nusselt number decreases with the increase of power-law index, and average entropy generation increases with power-law index. We believe that the obtained result in the present study shall be useful for design of energy efficient microsystems which utilize the dual electrokinetic and centrifugal pumping effects.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Discovery of microfluidic devices requires designing a suitable pumping system which is of considerable importance. Among various different techniques developed such as magnetohydrodynamics, piezoelectrics, and electrohydrodynamics [1–3], electroosmosis [4] has received special attention for lab-on-a-chip microfluidic devices because of the simple design and easier fabrication of electroosmotic pumping systems. The aforementioned techniques are generally conducted in micron-sized ducts with arbitrary cross section; therefore, the name microfluidics. There are several advantages to the use of microfluidic devices, such as low energy and material consumption, higher accuracy, easier control, and automation. These are the main reasons that made microfluidics one of the most attractive research fields in recent years.

Since the generation of fluid flow is crucially important at microscale, a classical pressure gradient-driven mechanism which contains moving components is extremely difficult to be designed and manufactured at microscales and is mostly prone to mechanical failure. In this respect, generation of the electroosmotic flow (EOF) system is the only unique alternative, indeed, establishing fluid flow called electroosmosis. Moreover, since electroosmotic pumping involves only application of the electric field on a duct with no moving components, electroosmotical design, and fabrication-actuated microfluidic devices, this process is much easier than their pressure-based counterparts. In 1809, Reuss [5] was first to report the study of EOF. He showed that by the application of suitable electric voltage, water can flow through a plug of clay. This work was followed by the theoretical study of Helmholtz in 1879 [6] on the electric double layer (EDL). In the early 1900s, Smoluchowski [7] investigated electrokinetically driven flows under the conditions where the EDL thickness is much smaller than the channel height.

Regarding to the research on hydrodynamics of electroosmotic flow, Burgreen and Nakatche [8] and Rice and Whitehead [9] considered the electrokinetic flow for slit and cylindrical capillaries, respectively, where they assumed low values of the zeta potential and used Debye–Hückel linearization. Later, Levine et al. [10] used analytical approximation solution of Poisson–Boltzmann equation with high zeta functions. When we have the situation that two cylindrical walls carry high zeta potentials, Kang et al. [11] analytically analysed the electroosmotic flow through an annulus. Yang [12] and Wang et al. [13] obtained analytical solutions for fully developed and pressure effects on electroosmotic flow in rectangular and semicircular microchannels. Recently, Azari et al. [14] considered Graetz problem analytically for the electroosmotic flow in microchannels and pressure-driven with distributed wall heat flux; their results indicate that the average Nusselt number is a decreasing function of pressure-driven velocity and the electric double layer thickness, regardless of the wall heat flux distribution.

Biofluids encountered in many lab-on-chip devices show non-Newtonian behavior; their rheological properties significantly differ from the Newtonian law of viscosity, and usually, their viscosities are dependent on the rate of shear. Hence, analysing such fluid under the action of the electroosmotic force is necessary for design the lab-on-chip devices. Literature reviews indicate that there is a growing interest in modelling of electroosmotic flow of non-Newtonian fluid. Under the sole influence of electrokinetic forces, in a rectangular microchannel, Das and Chakraborty [15] investigated the transport characteristics of a non-Newtonian fluid flow. They used Debye–Hückel linearization and obtained an analytical solution for velocity and temperature field. Their results reveal that significant reductions in species concentration levels characterized by more significant viscous effects may be achieved by higher hematocrit fraction on account of stronger dispersions in the velocity profiles. This work was followed by Zhao and coworkers [16], first in [17], and they found the analytical solution for power-law fluids in a slit channel and small zeta function; later in [18], this group obtained the expression for the general Smoluchowski velocity by solving the same problem without restriction on zeta function for electroosmosis of power-law fluids in a similar fashion to the classic Smoluchowski velocity for Newtonian fluids. The group later extended their study to the EOF of non-Newtonian power-law fluids in a cylindrical microchannel and found the exact solution for power-law index $n=1\text{\hspace{0.17em}and\hspace{0.17em}}0.5$ and approximate solutions obtained for arbitrary values of fluid behavior index. The steric effect induced alterations in streaming potential, and energy transfer efficiency of power-law fluid is investigated analytically by Bandopadhyay and Chakraborty [19], and they found that giant augmentations in the energy transfer efficiency may be caused by confluence of the steric interactions with the non-Newtonian transport characteristics for shear thickening fluids under appropriate conditions. In rectangular microchannels, Vakili et al. [20] analysed electrokinetically driven fluidic transport of power-law fluids, and their result indicates that Poiseuille number is an increasing function of the zeta potential, the flow behavior index, channel aspect ratio, and the dimensionless Debye–Hückel parameter. Regarding to triangular geometry, Chaves et al. [21] used generalized integral transform technique and employed the laminar forced convection (hydrodynamically fully developed and thermally developing laminar flow which was generated by pressure gradient) of power-law non-Newtonian fluids inside ducts with arbitrary-shaped cross sections. We have also more recent study of Mukherjee et al. [22]; they used COMSOL commercial package programme and employed the laminar forced convection in power-law and Bingham plastic fluids in ducts of semicircular and other cross-section problems where the flow is generated by the gradient of pressure. They show that their results reduce to the power-law fluid prediction for zero and infinite shear viscosities. Effects of homogeneous-heterogeneous reactions in MHD flow of Casson fluid flow over a stretched surface is considered by Khan et al. [23].

Combination of electroosmotic and pressure-driven force may be involved in many practical applications. However, there are significant differences between the hydrodynamic characteristics of combined pressure-driven and electroosmotic flow (PDEOF) from those of both conventional pressure gradient-driven flow (PDF) and pure electroosmotic flow (EOF). Literature is very rich also in this topic. But, in this study, we concern in power-law fluid; in this regard, we have the study of Vakili et al. [24] in which they consider pressure effects on electroosmotic flow of power-law fluids in rectangular microchannels, and they later studied thermal transport characteristic as well in [25]. Thermally fully developed case for electroosmotic flow of power-law nanofluid in a rectangular microchannel is investigated by Deng in [26]; their result indicates an increase in the Nusselt number with the flow behavior index and with electrokinetic width when considering the surface heating effect, which decreases with the Joule heating parameter.

Recently, there is a growing interest of scientist on the physical importance for the entropy generation analysis. This concept explored in the pioneering work of Bejan [27] where the author considered the thermodynamic second law features of heat transfer by forced convection on four different flow configurations: pipe flow, flow in a rectangular duct, boundary layer over flat plate, and cross-flow; then, he analysed the irreversibility due to heat transfer through finite temperature gradient and irreversibility due to the viscous effect. In a series of research studies, Khan and his co-authors investigated the entropy generation for a different flow and heat transfer problem. Khan et al. in [28] studied second order velocity slip flow of viscous fluid by a variable thickened stretchable surface of disk problem and obtained the numerical results for entropy generation. We also note following studies on entropy generation analysis for the flow and heat transfer problem in [29–31].

Microchannel has close to a rectangular shape cross section in most lab-on-chip systems [32, 33]; but due to the limitation of space in some cases, a microchannel of triangular cross-section must be used. Thus, the flow of fluids in a triangular microchannel has received special attention [34]. Literature reviews reveal that no study has yet explored the electroosmotic flow thermal features of power-law fluids in a microtriangular duct. This provides enough motivation for the current study. The numerical method is developed to explore the thermal characteristics and entropy generation of mixed electroosmotic and pressure-driven flow of power-law fluids in a microtriangular duct by considering the viscous dissipation effects. The triangular contour is chosen because it is the cross section representing the largest deviation from the circular in the family of axially symmetric tube contours.

The organization of this paper is as follows. Formulation of novel problem is given in Section 2. Governing equations for our novel flow and heat transfer problem are derived in Section 3. We studied the numerical method in Section 4, Section 5 devoted to the introduction of the numerical method, and entropy generation analysis is given Section 6. Grid independence of the numerical method is studied in Section 7, and numerical results are discussed in Section 8. New finding documented in the final conclusion section where we state new result that the Nusselt number decreases with increase of power-law index and average entropy generation increases with power-law index and discussion of the results.

2. Formulation of Problem

We consider the heat transfer associated with a mixed electroosmotic and pressure gradient-driven flow of a power-law fluid in a left triangular duct with the dimension given in Figure 1. The flow considered to be steady laminar, both hydrodynamically and thermally, developed. We assume that fluid has constant thermophysical properties, and it contains an ideal solution of fully dissociated salt. The EDL is assumed to have a constant zeta potential on the wall (Stern layer).

[figure omitted; refer to PDF]

The triangular duct wall is a subject to zeta potential which uniformly distributed over the walls. We also assume that EDLs formed on the triangular wall do not overlap. Furthermore, when we calculate the potential field and the charge density, we assume that the temperature variations over the triangular cross section are negligible as compared to the absolute temperature. We calculate the charge density based on the velocity-weighted average (bulk) temperature. Therefore, Boltzmann equation can be used to describe the spatial distribution of the electric charge density.

3. Governing Equations

The following Poisson’s equation will be used to describe the electrical potential distribution within the microtriangular:$$\begin{array}{}\text{(1)}& {\nabla }^2\phi =-\frac{{\rho }_e}{\epsilon },\end{array}$$where $\epsilon$ represents the fluid permittivity, and ${\rho }_e$ is the net electric charge. The potential $\varphi$ is due to the combination of externally imposed field and EDL potential; therefore, equation (1) can be written as$$\begin{array}{}\text{(2)}& \phi =\varphi +\psi .\end{array}$$

The electric charge density for an ideal solution of fully dissociated salt is given by [24]$$\begin{array}{}\text{(3)}& {\rho }_e=-2n_{0}eZ\sinh \frac{eZ\psi }{k_{B}T},\end{array}$$where $n_0,e,Z,T,\text{and\hspace{0.17em}}{k}_B$ is the ion density, valence of ions, proton charge, absolute temperature, and Boltzmann constant, respectively. Since our flow problem hydrodynamically fully developed, EDL potential is reduced to $\psi =\psi x,y$, and the external potential gradient here is in the axial direction only which means $\varphi =\varphi z$, where $z$ stands for the axial coordinate. For a constant voltage gradient in the $z$-direction, equation (1) becomes$$\begin{array}{}\text{(4)}& \frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}=\frac{2n_{0}eZ}{\epsilon }\sinh \frac{eZ\psi }{k_{B}T}.\end{array}$$

The effect of temperature on the potential distributions is shown to be negligible [24]; hence, potential distribution can be calculated on the basis of average temperature over the triangular region. Now, we write the equation (4) in dimensionless form as$$\begin{array}{}\text{(5)}& \frac{{\partial }^2{\psi }^{\ast }}{\partial x^{\ast 2}}+\frac{{\partial }^2{\psi }^{\ast }}{\partial y^{\ast 2}}=K^2\sinh {\psi }^{\ast },\end{array}$$where ${\psi }^{\ast }=\epsilon Z\psi /k_B{T}_{av},x^{\ast }=x/H,y^{\ast }=y/H$, and $K$ represents the dimensionless Debye–Hückel parameter given by $K=H/{\lambda }_D$ with ${\lambda }_D={\epsilon k_B{T}_{av}/2n_0{e}^2{Z}^2}^{1/2}$ (we omit the star expression). It is apparent from the definition that a large value of $K$ denotes a relatively thin electric double layer and vice versa. In this study, we consider equation (1) with the following boundary conditions:$$\begin{array}{c}\text{(6)}& {\psi }^{\ast }1,y={\zeta }^{\ast },{\psi }^{\ast }x,0={\zeta }^{\ast },\\ {\psi }^{\ast }x,x={\zeta }^{\ast }\text{,}\\ \text{or}\\ {\psi }^{\ast }y,y={\zeta }^{\ast },\hspace{1em}0\le x,y\le 1.\end{array}$$

The governing equation of flow field for the power-law fluid is the continuity, and Cauchy momentum equations are given as$$\begin{array}{}\text{(7)}& \nabla \cdot \text{u}=0,\text{(8)}& \rho \frac{\text{Du}}{\text{Dt}}=-\nabla p+\nabla \cdot \mathrm{\tau }+F,\end{array}$$where $\text{u}$ is the velocity vector, $p$ is the pressure, $\rho$ is the density, $\tau$ is the total stress tensor (Cauchy extrastress tensor), and $F$ is the body force. The total stress tensor can be related to the strain rate tensor (generalized Newtonian fluid) as$$\begin{array}{}\text{(9)}& \tau =2\mu \dot{\gamma }\dot{\gamma },\end{array}$$where $\dot{\gamma }=\nabla \text{u}+\nabla {\text{u}}^T/2$, and $\mu \dot{\gamma }$ is the effective viscosity with $\dot{\gamma }$ as the magnitude of the strain tensor, which is defined as$$\begin{array}{}\text{(10)}& \dot{\gamma }={\frac{1}{2}\dot{\gamma }\cdot \dot{\gamma }}^{1/2}.\end{array}$$

The effective viscosity for power-law fluid is modelled as$$\begin{array}{}\text{(11)}& \mu \dot{\gamma }=\eta {2\dot{\gamma }}^{n-1},\end{array}$$where $\eta$ is the flow consistency index, and $n$ is the flow behavior index. We note that power-law fluid is subdivided into three different types of fluids based on the flow behavior index that $n<1$ corresponds to shear thinning $n=1$ Newtonian and $n>1$ shear thickening fluids, respectively. Since we also assumed that the flow is fully developed, in this case, velocity vector becomes $\text{u}=\mathrm{0,0},wx,y$, then continuity equation is automatically satisfied, and momentum equation (8) is reduced to$$\begin{array}{}\text{(12)}& 0=-\frac{\text{d}p}{\text{d}z}+\frac{\partial }{\partial x}{\tau }_{xz}+\frac{\partial }{\partial y}{\tau }_{yz}+F_z.\end{array}$$

The z-component of the electric body force, generated due to the interaction of the electric charge density in EDL and the electric field, equals ${\rho }_e{E}_z$ with $E_z=-\text{d}/\text{d}z$ representing the electric field in the axial direction. By substituting ${\rho }_e$ from equation (3), the $z$-component of the electric body force becomes$$\begin{array}{}\text{(13)}& F_z=-2n_{0}Z{E}_z\sinh \frac{eZ\psi }{k_B{T}_{\text{av}}}.\end{array}$$Magnitude of rate of strain and effective viscosity is obtained as$$\begin{array}{}\text{(14)}& \dot{\gamma }=\frac{1}{2}{{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{1/2},\text{(15)}& \mu \dot{\gamma }=\eta {{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n-1/2}.\end{array}$$

Now, substituting (14) and (15) in (9), the shear stresses reduce to$$\begin{array}{}\text{(16)}& {\tau }_{xz}=\eta {{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n-1/2}\frac{\partial w}{\partial x},\text{(17)}& {\tau }_{yz}=\eta {{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n-1/2}\frac{\partial w}{\partial y}.\end{array}$$

Substituting (13), (16), and (17) into the momentum equation (12), we obtain$$\begin{array}{}\text{(18)}& \frac{\partial }{\partial x}\eta {{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n-1/2}\frac{\partial w}{\partial x}+\frac{\partial }{\partial y}\eta {{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n-1/2}\frac{\partial w}{\partial y}-\frac{dp}{dz}-2n_{0}Z{E}_z\sinh \frac{eZ\psi }{k_B{T}_{av}}=0.\end{array}$$

After performing the differentiation, we obtain$$\begin{array}{}\text{(19)}& \eta n{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2\frac{{\partial }^{2}w}{\partial x^2}+\eta n{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2\frac{{\partial }^{2}w}{\partial z^2}+2n-1\eta \frac{\partial w}{\partial x}\frac{\partial w}{\partial y}\frac{{\partial }^{2}w}{\partial x\partial y}-{{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{3-n/2}\frac{dp}{dz}+2n_{0}Z{E}_z\sinh \frac{eZ\psi }{k_B{T}_{av}}=0.\end{array}$$

The Helmholtz–Smoluchowski electroosmotic velocity is considered here as the reference velocity, $U_{\text{HS}}$, and this for power-law fluids at small zeta potentials becomes [24]$$\begin{array}{}\text{(20)}& U_{\text{HS}}=n{\lambda }_D^{n-1/n}{-\frac{\epsilon \zeta E_z}{\eta }}^{1/n}.\end{array}$$

After defining dimensionless velocity as $w^{\ast }=w/U_{\text{HS}}$, the dimensionless form of momentum equation is given by (again omitting star expression)$$\begin{array}{c}\text{(21)}& n{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2\frac{{\partial }^{2}w}{\partial x^2}+n{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2\frac{{\partial }^{2}w}{\partial z^2}+2n-1\frac{\partial w}{\partial x}\frac{\partial w}{\partial y}\frac{{\partial }^{2}w}{\partial x\partial y}=-{{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{3-n/2}{\frac{n+1}{n}}^n\Gamma +\frac{K^{n+1}}{n^n\zeta }\sinh \psi ,\end{array}$$where $\mathrm{\Gamma }=U_{\text{PD}}^n/U_{\text{HD}}^n$, and $U_{\text{HD}}^n={n+1/n}^n-1/\eta \text{d}p/\text{d}z{H}^{n+1}$. This is the subject to the following boundary conditions:$$\begin{array}{c}\text{(22)}& w1,y=0,wx,0=0,\\ wx,y=0\text{\hspace{0.17em}on\hspace{0.17em}}x=y.\end{array}$$

Fundamental physical law of energy conservation will be used to obtain temperature distribution. Considering the effect of viscous dissipation and Joule heating, the energy equation is written as$$\begin{array}{}\text{(23)}& \rho c_p\frac{DT}{Dt}=\nabla \cdot k\cdot \nabla T+s+\mathrm{\tau }:\nabla w,\end{array}$$where $k$ is the thermal conductivity, $c_p$ is the specific heat, $s$ is the Joule heating term, and $\tau :\nabla w$ represents the volumetric heat generation caused by viscous dissipation. The term $s=E_z^2/\delta T$ represents the Joule heating, where $\delta T$ is the liquid electrical resistivity given by$$\begin{array}{}\text{(24)}& \delta T=\frac{{\delta }_0}{\cosh eZ\psi /k_B{T}_{av}},\end{array}$$where ${\delta }_0$ is the neutral liquid electrical resistivity. The term representing viscous dissipation can be expressed as$$\begin{array}{}\text{(25)}& \mathrm{\tau }:\nabla w={\tau }_{xz}\frac{\partial w}{\partial x}+{\tau }_{yz}\frac{\partial w}{\partial y}.\end{array}$$

Using (16) and (17) energy equations for steady fully developed flow of power-law fluid becomes$$\begin{array}{}\text{(26)}& \rho c_{p}wx,y\frac{\partial T}{\partial z}=k\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}T+\eta \Phi +\frac{E_z^2}{{\delta }_0}\cosh \frac{eZ\psi }{k_B{T}_{\text{av}}},\end{array}$$where$$\begin{array}{}\text{(27)}& \Phi ={{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n-1/2}{\frac{\partial w}{\partial x}}^2+{{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n-1/2}{\frac{\partial w}{\partial y}}^2.\end{array}$$

Since the flow is assumed to be thermally fully developed, this gives$$\begin{array}{}\text{(28)}& \frac{\partial }{\partial z}\frac{T_w-T}{T_w-T_m}=0,\end{array}$$where $T_m$ and $T_w$ are the bulk and wall temperatures, respectively. The convection rate equation could be written as$$\begin{array}{}\text{(29)}& q_w=\frac{k}{H}{T}_w-T_m,\end{array}$$representing a constant wall heat flux $q_w$ boundary condition (invariant of axial coordinates), and then from (28), we have$$\begin{array}{}\text{(30)}& \frac{\partial T}{\partial z}=\frac{\text{d}{T}_m}{\text{d}z}=\frac{\text{d}{T}_w}{\text{d}z}=\text{constant}.\end{array}$$

Substituting the above equation into (26), we find$$\begin{array}{}\text{(31)}& \rho c_{p}wx,y\frac{\text{d}{T}_m}{\text{d}z}=k\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}T+\eta \Phi +\frac{E_z^2}{{\delta }_0}\cosh \frac{eZ\psi }{k_B{T}_{\text{av}}},\end{array}$$where$$\begin{array}{}\text{(32)}& \Phi ={{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n+1/2}.\end{array}$$

The following relation represents the energy balance accounting for heat fluxes on walls of triangle, Joule heating on an elemental control volume, and viscous dissipation:$$\begin{array}{}\text{(33)}& \rho c_{p}wx,y\frac{\text{d}{T}_m}{\text{d}z}=\frac{4q_w}{d_h}+\frac{1}{H^2}\eta \Phi +\frac{E_z^2}{{\delta }_0}\cosh \frac{eZ\psi }{k_B{T}_{\text{av}}},\end{array}$$where the angle bracket denotes an average taken over the cross section, and $d_h=4A/P$ is the hydraulic diameter of the left angular duct.

For the fully developed flow, the dimensionless temperature is defined as$$\begin{array}{}\text{(34)}& \theta x,y=\frac{T-T_{w}z}{q_{w}H/k}.\end{array}$$

Substituting (33) in (31) and using dimensionless temperature function and dimensionless parameters which are defined before, we obtain$$\begin{array}{}\text{(35)}& \frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}=\frac{w}{w_m}\frac{4}{{\widehat{d}}_h}+S{I}_1+\frac{S{S}_v}{K^2}{I}_2-\frac{S{S}_v}{K^2}{{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n+1/2}-S\cosh \psi ,\end{array}$$where the dimensionless parameters are$$\begin{array}{c}\text{(36)}& S=\frac{E_z^{2}H}{q_w{\delta }_0},S_v=\frac{\eta {\delta }_0{U}_{HS}^{n+1}}{E_z^2{\lambda }_D^2{H}^{n-1}},\\ w_m={\displaystyle {\int }_0^1{\displaystyle {\int }_0^{x}wx,y\text{d}y\text{d}x}},\\ I_1={\displaystyle {\int }_0^1{\displaystyle {\int }_0^x\cosh \psi \text{d}y\text{d}x}},\\ I_2={\displaystyle {\int }_0^1{\displaystyle {\int }_0^x{{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n+1/2}\text{d}y\text{d}x}}.\end{array}$$

The dimensionless boundary conditions for the energy equation are$$\begin{array}{c}\text{(37)}& \theta x,0=0,\theta 0,y=0,\\ \theta x,y,z=0,\text{\hspace{0.17em}on\hspace{0.17em}}x-y=0.\end{array}$$

Based on the hydraulic diameter of the trigonometric geometry, the Nusselt number Nu can be represented by$$\begin{array}{}\text{(38)}& Nu=-\frac{{\widehat{d}}_h}{{\theta }_m},\end{array}$$where$$\begin{array}{}\text{(39)}& {\theta }_m=\frac{{\displaystyle {\int }_0^1{\displaystyle {\int }_0^{x}wx,y\theta x,y\text{d}y\text{d}x}}}{{\displaystyle {\int }_0^1{\displaystyle {\int }_0^{x}wx,y\text{d}y\text{d}x}}}=\frac{w\theta }{w}.\end{array}$$

4. Numerical Method

In this work, we studied the finite difference procedure with a nonuniform grid as described in [35] for our problem in (5), (21), and (35). We used following stretching functions to generate the nonuniform grid:$$\begin{array}{c}\text{(40)}& x_i=\frac{1}{2}-\frac{{\sin }^{-1}-\alpha \cos \pi i/N}{2{\sin }^{-1}\alpha },\hspace{1em}i=1,\dots ,N,y_j=\frac{1}{2}-\frac{{\sin }^{-1}-\alpha \cos \pi j/N}{2{\sin }^{-1}\alpha },\hspace{1em}j=1,\dots ,N.\end{array}$$

We now use finite difference approximation with a nonuniform grid for equations (5)–(21) and (35) as$$\begin{array}{c}\text{(41)}& \frac{1}{{\Delta x_i}^2}+\frac{1}{\Delta x_i\Delta x_{i-1}}+\frac{1}{{\Delta y_j}^2}+\frac{1}{\Delta y_j\Delta y_{j-1}}{\psi }_{i,j}=\frac{1}{{\Delta x_i}^2}{\psi }_{i+1,j}+\frac{1}{\Delta x_i\Delta x_{i-1}}{\psi }_{i-1,j}+\frac{1}{{\Delta y_j}^2}{\psi }_{i,j+1}+\frac{1}{\Delta y_j\Delta y_{j-1}}{\psi }_{i,j-1}+K^2\sinh {\psi }_{i,j},\text{(42)}& {{\frac{w_{i+1,j}-w_{i,j}}{\Delta x_i}}^2+{\frac{w_{i,j+1}-w_{i,j}}{\Delta j_j}}^2}^{n-3/2}n\frac{w_{i+1,j}-w_{i,j}}{\Delta x_i}+{\frac{w_{i,j+1}-w_{i,j}}{\Delta j_j}}^2+\frac{1}{{\Delta x_i}^2}{w}_{i+1,j}-\frac{1}{{\Delta x_i}^2}+\frac{1}{\Delta x_i\Delta x_{i-1}}{w}_{i,j}+\frac{1}{\Delta x_i\Delta x_{i-1}}{w}_{i-1,j}{\frac{w_{i+1,j}-w_{i,j}}{\Delta x_i}}^2+n{\frac{w_{i,j+1}-w_{i,j}}{\Delta j_j}}^2\\ \frac{1}{{\Delta j_j}^2}{w}_{i,j+1}-\frac{1}{{\Delta y_j}^2}+\frac{1}{\Delta y_j\Delta y_{j-1}}{w}_{i,j}+\frac{1}{\Delta y_j\Delta y_{j-1}}{w}_{i,j-1}\\ 2n-1\frac{w_{i+1,j}-w_{i,j}}{\Delta x_i}\frac{w_{i,j+1}-w_{i,j}}{\Delta j_j}\frac{w_{i+1,j+1}-w_{i+1,j-1}-w_{i-1,j+1}+w_{i-1,j-1}}{x_{i+1}-x_{i-1}{y}_{i+1}-y_{i-1}}\\ ={\frac{n+1}{n}}^n\Gamma +\frac{K^{n+1}}{n^n\zeta }\sinh {\psi }_{i,j},\\ \text{(43)}& -\frac{1}{{\Delta x_i}^2}+\frac{1}{\Delta x_i\Delta x_{i-1}}+\frac{1}{{\Delta y_j}^2}+\frac{1}{\Delta y_j\Delta y_{j-1}}{\theta }_{i,j}+\frac{1}{{\Delta x_i}^2}{\theta }_{i+1,j}+\frac{1}{\Delta x_i\Delta x_{i-1}}{\theta }_{i-1,j}+\frac{1}{{\Delta y_j}^2}{\theta }_{i,j+1}+\frac{1}{\Delta y_j\Delta y_{j-1}}{\theta }_{i,j-1}=\frac{w_{i,j}}{w_m}\frac{4}{{\widehat{d}}_h}+S{I}_1+\frac{S{S}_v}{K^2}{I}_2-\frac{S{S}_v}{K^2}{{\frac{w_{i+1,j}-w_{i,j}}{\Delta x_i}}^2+{\frac{w_{i,j+1}-w_{i,j}}{\Delta y_j}}^2}^{n+1/2}S\cosh {\psi }_{i,j}.\end{array}$$

Implementing boundary conditions (6), (22), and (36) easily and applying the predictor-corrector scheme as usual for nonlinear equations (41) and (42), we obtain following system of equations for each step:$$\begin{array}{}\text{(44)}& \frac{1}{{\Delta x_i}^2}+\frac{1}{\Delta x_i\Delta x_{i-1}}+\frac{1}{{\Delta y_j}^2}+\frac{1}{\Delta y_j\Delta y_{j-1}}{\psi }_{i,j}^k=\frac{1}{{\Delta x_i}^2}{\psi }_{i+1,j}^k+\frac{1}{\Delta x_i\Delta x_{i-1}}{\psi }_{i-1,j}^k+\frac{1}{{\Delta y_j}^2}{\psi }_{i,j+1}^k+\frac{1}{\Delta y_j\Delta y_{j-1}}{\psi }_{i,j-1}^k+K^2\sinh {\psi }_{i,j+1}^k,\hspace{1em}i,j=1,\dots ,N,\text{(45)}& A_{ij}^{k-1}\frac{1}{{\Delta x_i}^2}{w}_{i+1,j}^k-\frac{1}{{\Delta x_i}^2}+\frac{1}{\Delta x_i\Delta x_{i-1}}{w}_{i,j}^k+\frac{1}{\Delta x_i\Delta x_{i-1}}{w}_{i-1,j}^k+\frac{1}{{\Delta y_j}^2}{w}_{i,j+1}^k-\frac{1}{{\Delta y_j}^2}+\frac{1}{\Delta y_j\Delta y_{j-1}}{w}_{i,j}^k+\frac{1}{\Delta y_j\Delta y_{j-1}}{w}_{i,j-1}^k+C_{i,j}^{k-1}{D}_{i,j}^{k-1}\frac{1}{{\Delta x_i}^2}{w}_{i+1,j}^k-\frac{1}{{\Delta x_i}^2}+\frac{1}{\Delta x_i\Delta x_{i-1}}{w}_{i,j}^k\frac{1}{\Delta x_i\Delta x_{i-1}}{w}_{i-1,j}^k+C_{i,j}^{k-1}{E}_{i,j}^{k-1}\frac{1}{{\Delta y_j}^2}{w}_{i,j+1}^k-\frac{1}{{\Delta y_j}^2}+\frac{1}{\Delta y_j\Delta y_{j-1}}{w}_{i,j}^k+\frac{1}{\Delta y_j\Delta y_{j-1}}{w}_{i,j-1}^k+F_{i,j}^{k-1}\frac{w_{i+1,j}^k-w_{i,j}^k}{\Delta x_i}+G_{i,j}^{k-1}\frac{w_{i,j+1}^k-w_{i,j}^k}{\Delta y_i}={\frac{n+1}{n}}^n\Gamma +\frac{K^{n+1}}{n^n\zeta }\sinh {\psi }_{i,j}.\end{array}$$

First guest is taken for (45) from Newtonian case $n=1$ in which case exact analytical solution is possible; for (44), we can use Debye–Hückel linearization. We then apply the successive overrelaxation method to solve equations (44) and (45). In each iteration, the value of ${\psi }_{i,j}^{\text{prev}}$ is used in one previous iteration and updated for each iteration for (44). We continue until we achieve the required overall relative error of ${10}^{-6}$. Similar logic is used to solve equation (45). One’s solutions are obtained for potential and velocity field. Energy equation in (43) was solved by the SOR method again, and results are given in the following section.

5. Entropy Generation Analysis

In Sections 2 and 3, we discussed thermodynamical behaviors of the system; it is well known that entropy generation plays an important role in this field. We now study analysis of entropy generation for our current flow problem along the same line as in Bejan [27]. The generation of total entropy is the sum of viscous dissipation, Joule heating, and thermal gradients, and this expression can be mathematically formulated as$$\begin{array}{}\text{(46)}& S_G^{\ast }=S_T^{\ast }+S_F^{\ast }+S_J^{\ast },\end{array}$$where, in the above right-hand side terms, $S_T^{\ast }$ is the entropy generation rate because of heat transfer, $S_F^{\ast }$ is the fluid friction, and $S_J^{\ast }$ represents the Joule heating irreversibilities, and these terms can be expressed by$$\begin{array}{c}\text{(47)}& S_T^{\ast }=\frac{k}{T_0^2}{\nabla T}^2,S_F^{\ast }=\frac{1}{T_0}\tau :\nabla w,\\ S_J^{\ast }=\frac{Q_j}{T_0}.\end{array}$$

The dimensionless form of the total entropy generation can be derived as $S_G=S_T+S_F+S_J$, where$$\begin{array}{c}\text{(48)}& S_T={\displaystyle {\int }_0^1{\displaystyle {\int }_0^x{S}_T^{c}x,y\text{d}y\text{d}x}},S_F={\displaystyle {\int }_0^1{\displaystyle {\int }_0^x{S}_F^{c}x,y\text{d}y\text{d}x}},\\ S_J={\displaystyle {\int }_0^1{\displaystyle {\int }_0^x{S}_J^{c}x,y\text{d}y\text{d}x}},\end{array}$$whereas$$\begin{array}{c}\text{(49)}& S_T^c={\frac{\partial \theta }{\partial x}}^2+{\frac{\partial \theta }{\partial y}}^2,S_F^c=\frac{S{S}_v}{R{K}^2}{{\frac{\partial w}{\partial x}}^2+{\frac{\partial w}{\partial y}}^2}^{n+1/2},\\ S_J^c=\frac{S}{R}\cosh \psi ,\end{array}$$where $R=q_{w}H/T_{0}k$.

6. Method Validation

Excellent convergence was achieved for all the results. The following algorithm was used:

Step 1 Specify values of $\zeta ,n,\text{and}K$

Step 2 Solve equation (44) by using a predictor-corrector along with the SOR method

Step 3 Substitute the solution of equation (44) into (45), solve equation (45) for $n=1$, and then, use this solution as an initial velocity solution for $n\ne 1$ where we used again a predictor-corrector along with the SOR method

Step 4 Substituting the solution of equation (45) into (43), solve equation (43) with the SOR method

Step 5 Calculate the entropy generation with a variety of parameters

In order to the reliability of the numerical result, we need to validate numerical result. We first check the grid dependency of the numerical result for velocity result. As usual, we first obtained approximate solution for the grid system $101\times 101$, and compared with the gridding $151\times 151$ and $201\times 201$, we have seen that differences between obtained results for mean velocity is less than $1\times {10}^{-5}$; hence, we use system $151\times 151$ for all our computation in the numerical results.

Besides the analysis of grid dependency, we need to compare our numerical result with the benchmark data. The results are compared with two limiting cases: one of them is the Debye–Hückel linearization which can be applied to the equation (5) which is $\sinh {\psi }^{\ast }\approx {\psi }^{\ast }$; hence, equation (5) becomes linear. We then use the eigen function expansion as in the work of Wang [36] to find the approximate analytical solution of (5), and we first change the dependent variable as $\psi ={\zeta }^{\ast }-Q$, and then$$\begin{array}{c}\text{(50)}& \frac{{\partial }^{2}Q}{\partial x^2}+\frac{{\partial }^{2}Q}{\partial y^2}\approx K^2{\zeta }^{\ast }-Q,\hspace{1em}\text{with\hspace{0.17em}}Q1,y=0,Qx,0=0,Qx,x=0.\end{array}$$

We now express the unknown functions and unity in terms of Helmholtz eigen functions as$$\begin{array}{c}\text{(51)}& Qx,y={\displaystyle \sum _j{A}_j{\phi }_j},1={\displaystyle \sum _j{a}_j{\phi }_j}.\end{array}$$

Substituting (51) into (50), we obtain$$\begin{array}{}\text{(52)}& A_j=\frac{K^2{\zeta }^{\ast }{a}_j}{K^2+{\lambda }_j},\hspace{1em}j=1,\dots ,N.\end{array}$$

Eigen values and eigen functions of Helmholtz eigenvalue problem for our geometry are given by$$\begin{array}{c}\text{(53)}& {\lambda }_{mn}={\pi }^2{m}^2+n^2,{\phi }_{mn}=\sin n\pi x\sin m\pi y-\sin n\pi y\sin m\pi x.\end{array}$$

We now compare the exact solution (51) with the result of the finite difference solution (44) for ${\zeta }^{\ast }=1\text{\hspace{0.17em}and\hspace{0.17em}}K=5$ which is shown in Figure 2(a); we see that differences between the exact and numerical solution are negligible $\text{less\hspace{0.17em}than\hspace{0.17em}}{10}^{-4}$ for $151\times 151$ number of grid points. In order to check the grid-independent result, we increase the number of grids to $201\times 201$, and we see that differences between the two numerical solution are less than ${10}^{-5}$. Therefore, we conclude that the result of $151\times 151$ number of grid points is enough for the grid-independent result.

[figures omitted; refer to PDF]

The other one is the simplification of the momentum equation for the Newtonian case. In the case of pressure gradient-driven flow and left triangular geometry, the exact solution is given by Pipkin and Rivlin (1963) in [37] after transformation as$$\begin{array}{}\text{(54)}& wx,y=-\frac{3}{2}yx-yx-1.\end{array}$$

We compared this exact analytical solution with the numerical solution equation (45) with $n=1$ in Figure 2(b) which shows the satisfactory result.

7. Results and Discussion

In this section, numerical results are presented to show the effect of parameters on the velocity and temperature field where we assumed that the EOF (electroosmotic flow) field is hydrodynamically and thermally fully developed. Furthermore, we discuss the effect of triangular geometry and flow behavior concerning to the non-Newtonian rheology on the velocity and temperature field in a detailed analysis. The electroosmotic flow is driven by the combined mechanism of the pressure gradient and axially imposed electric field. The dynamic behavior of the power-law fluid is studied for combined mechanism EOF pressure-driven flow with no-slip conditions at the boundary of triangle. We note here that power-law fluid depending on power-law index exhibits shear thinning or shear thickening behaviors. In this study, the characteristic of fluid flow and heat transfer are studied numerically by using the finite difference method with a nonuniform grid; when the triangle length $L=5\times {10}^{-3}\text{m}$, height $H=10.0\mu \text{m}$, electron charge $e=1.6\times {10}^{-19}\text{C}$, field strength of electricity $E_0=1.0\times 105\text{V}/\text{m}$, Boltzmann constant $k_B=1.3805\times {10}^{-23}\text{J}/\text{K}$, $T_{ab}=300\text{K}$ which is absolute temperature, we considered valence of ion $z$ as 1, fluid permittivity $\epsilon =6.95039\times {10}^{-10}\text{C}/\text{Vm}$, heat flux of wall $q_w=1.5\times {10}^3\text{W}/{\text{m}}^2$, electrical conductivity $\sigma =0.001\text{S}/\text{m}\text{\hspace{0.17em}to\hspace{0.17em}}1.0\text{S}/\text{m}$, specific heat $c_p=3760\text{J}/\text{kg\hspace{0.17em}K}$, and thermal conductivity $k=613\times {10}^{-3}\text{W}/\text{mK}$ [38].

There are several methods to solve Poisson–Boltzmann equation which was outlined by Lu et al. [39]. In this work, we chose the finite difference method with a nonuniform grid. Next, we investigate the impact of power-law index on the centerline velocity field. Here, we have investigated three different extreme cases, which correspond to the three different values of the scale ratio which are $\mathrm{\Gamma }=-1$ corresponding to pressure opposing flow, $\mathrm{\Gamma }=0$ corresponding to purely electrokinetic flow, and $\mathrm{\Gamma }=1$ corresponding to pressure assisting flow. For $n=0.6,$ it can be seen in Figure 3 that velocity increases as $\mathrm{\Gamma }$ increases. Figures 4 and 5 show the dimensionless velocity distribution at different values of $\Gamma$ and $n$ for ${\zeta }^{\ast }=1\text{\hspace{0.17em}and\hspace{0.17em}}K=10$; we see that maximum occurs at the center of the triangle. These are true relatively small values of $\Gamma$ in comparison to Debye–Hückel parameter, where $n=0.6\text{\hspace{0.17em}and\hspace{0.17em}}1.2$ corresponding to shear thinning and shear thickening behavior of power-law fluid, respectively. Note that the center of contour line velocity distributions for $n=1.2$ is much smaller than velocity distributions for $n=0.6$. It is clear from Figures 6 and 7 that the Joule heating effect cannot be ignored, and it depends on the value of $S$. Also, there is decrease or increase of the fluid temperature at a uniform rate. We note here that fluid temperature is negative for the case of heat generation $S>0$. However, it is positive for the heat absorption process $S<0$. We also bring about another important observation, which is the parabolic shape of the temperature profile; during heat generation, it is concave and convex during heat absorption in all of the cases examined here. Figures 6 and 7 are explained in detail. Figure 6 explores the effect of zeta potential on the temperature field where power-law index is considered as $n=0.6\text{\hspace{0.17em}and\hspace{0.17em}}1.2,$ respectively. In these figures, we take the other dimensionless parameters, which are used in calculations to be $K=10,S=-4,\text{\hspace{0.17em}and\hspace{0.17em}SV}=5.$ We observe that an increase in the value of dimensionless zeta potential would cause the increase in the temperature difference profiles inside the triangle. Increasing the value of $\zeta$, for which the amount of Joule heating in the proximity of the wall of the triangle increases, results in higher wall temperature which leads to a higher temperature difference. It is important to note that we used $81\times 81$ grid points for more visibility to sketch these figures. For the present problem, variation of Nusselt number with power-law index is given in Figure 8; we can clearly observe that the Nusselt number decreases with increase in the power-law index, in the case of shear thinning fluid, and becomes smaller for a shear thickening fluid. We note here that this decrease is not linear, and further increase in the value of $\zeta$ reduces the value of the Nusselt number as seen in the same figure. Again in the case of heat absorption characterized by negative values of the Joule heating parameter $S=-4$, the Nusselt number increases as the power-law index increases, which is not reported before in any other geometry. In the initial phase, the velocity and the temperature field in the flow of a fluid are determined properly with momentum and energy equations along with the microduct system; following this aspect, we tried to obtain the thermodynamic irreversibilities because of heat conduction, fluid friction, and Joule heating, and we derived explicit expressions concerning entropy generation analysis. As we stated earlier, temperature variation in the microtriangular duct system is essentially dominated by the Joule heating effect because of the external applied electric field. Therefore, the effect of the power-law effect on characteristics of the entropy production is consisted by the conductive heat transfer and viscous dissipation effects. Increasing the value of power-law index and $\zeta$ increases the average entropy, which is expected, because we obtained large gradients of the flow velocity and temperature in the microtriangular duct. Comparing entropy generation with shear thickening fluid, we can see that entropy generation is more pronounced for shear thickening fluid than that of the shear thinning fluid, which is seen clearly in Figure 9 for higher values of wall zeta function. We obtained larger average entropy generation as seen in Figure 9. We also examined the average entropy generation for negative $S$, and we found that average entropy generation increases with power-law index. Finally, we tried to find the effect of the Debye–Hückel parameter on temperature, Nusselt number, and average entropy generation while keeping the other parameters fixed. Figure 10 shows the effect of centerline temperature distribution, and we can see that the value of temperature distribution increases with increasing the Debye–Hückel parameter. Figure 11 represents the effect of Debye–Hückel parameter on the Nusselt number, and we see that the Nusselt number increases with the increase in the value of the Debye–Hückel parameter. Finally, Figure 12 shows the effect of the Debye–Hückel parameter on the average entropy generation, and we can see that the average entropy generation decreases with the increase in the value of the Debye–Hückel parameter as expected.

[figure omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figure omitted; refer to PDF]

[figures omitted; refer to PDF]

[figures omitted; refer to PDF]

[figure omitted; refer to PDF][figure omitted; refer to PDF][figure omitted; refer to PDF]

8. Conclusions

In this study, flow of power-law fluid, convective heat transfer, and entropy generation characteristics are considered for the microtriangular prism for both hydrodynamically and thermally developed EOF. We developed numerical algorithm based on the finite difference method with a nonuniform grid. The effects of key parameters such as power-law index in particular, viscous dissipation, Joule heating parameter on the flow, heat transfer, and entropy production are examined for the triangular geometry.

The following result can be drawn from the current research:

(a) A new finite difference algorithm with a nonuniform grid is developed for our flow problem

(b) Suitable adjustment of the Joule heating may control the fluid temperature distribution

(c) Nusselt number decreases with the increase of power-law index, and average entropy generation increases with power-law index

(d) We could use algorithm developed here for the simulation of electrokinetically modulated transport of blood at the microcirculatory system; one example concerning this system is where the human body is subjected to modern diagnostic tests such as MRI

The results of the numerical algorithm undertaken in this research are useful to verify theoretical/experimental models that may be involved with more complex three-dimensional electrokinetic flows. Also, we believe that concepts of the present research may be used to develop more realistic ionic tissue models to investigate the electrokinetic blood flow.

References

[1] A. V. Lemoff, A. P. Lee, "An AC magnetohydrodynamic micropump," Sensors and Actuators B: Chemical, vol. 63 no. 3, pp. 178-185, DOI: 10.1016/s0925-4005(00)00355-5, 2000.

[2] H. T. G. van Lintel, F. C. M. van de Pol, S. Bouwstra, "A piezoelectric micropump based on micromachining of silicon," Sensors and Actuators, vol. 15 no. 2, pp. 153-167, DOI: 10.1016/0250-6874(88)87005-7, 1988.

[3] A. Richter, A. Plettner, K. A. Hofmann, H. Sandmaier, "A micromachined electrohydrodynamic (EHD) pump," Sensors and Actuators A: Physical, vol. 29 no. 2, pp. 159-168, DOI: 10.1016/0924-4247(91)87118-m, 1991.

[4] S. Arulanandam, D. Li, "Liquid transport in rectangular microchannels by electroosmotic pumping," Colloids and Surfaces A: Physicochemical and Engineering Aspects, vol. 161 no. 1, pp. 89-102, DOI: 10.1016/s0927-7757(99)00328-3, 2000.

[5] F. F. Reuss, "Sur un nouveleffet de I’electricitegalvanique," Memoires de la SocieteImperiale des Naturalistes de Moscou, vol. 2, pp. 327-337, 1809.

[6] H. Helmholtz, "Studien uber electrische Grenzschichten," Annalen der Physik und Chemie, vol. 243 no. 7, pp. 337-382, DOI: 10.1002/andp.18792430702, 1879.

[7] M. Smoluchowski, ElektrischeEndosmose und Stromungsstr, 1914.

[8] D. Burgreen, F. R. Nakache, "Electrokinetic flow in ultrafine capillary Slits1," The Journal of Physical Chemistry, vol. 68 no. 5, pp. 1084-1091, DOI: 10.1021/j100787a019, 1964.

[9] C. L. Rice, R. Whitehead, "Electrokinetic flow in a narrow cylindrical capillary," The Journal of Physical Chemistry, vol. 69 no. 11, pp. 4017-4024, DOI: 10.1021/j100895a062, 1965.

[10] S. Levine, J. R. Marriott, G. Neale, N. Epstein, "Theory of electrokinetic flow in fine cylindrical capillaries at high zeta-potentials," Journal of Colloid and Interface Science, vol. 52 no. 1, pp. 136-149, DOI: 10.1016/0021-9797(75)90310-0, 1975.

[11] Y. Kang, C. Yang, X. Huang, "Dynamic aspects of electroosmotic flow in a cylindrical microcapillary," International Journal of Engineering Science, vol. 40 no. 20, pp. 2203-2221, DOI: 10.1016/s0020-7225(02)00143-x, 2002.

[12] D. Y. Yang, "Analytical solution of mixed electroosmotic and pressure-driven flow in rectangular microchannels," Key Engineering Materials, vol. 483, pp. 679-683, DOI: 10.4028/www.scientific.net/kem.483.679, 2011.

[13] C. Y. Wang, Y. H. Liu, C. C. Chang, "Analytical solution of electro-osmotic flow in a semicircular microchannel," Physics of Fluids, vol. 20 no. 6,DOI: 10.1063/1.2939399, 2008.

[14] M. Azari, A. Sadeghi, S. Chakraborty, "Graetz problem for combined pressure-driven and electroosmotic flow in microchannels with distributed wall heat flux," International Journal of Heat and Mass Transfer, vol. 128, pp. 150-160, DOI: 10.1016/j.ijheatmasstransfer.2018.08.106, 2019.

[15] S. Das, S. Chakraborty, "Analytical solutions for velocity, temperature and concentration distribution in electroosmotic microchannel flows of a non-Newtonian bio-fluid," Analytica Chimica Acta, vol. 559 no. 1, pp. 15-24, DOI: 10.1016/j.aca.2005.11.046, 2006.

[16] C. Zhao, E. Zholkovskij, J. H. Masliyah, C. Yang, "Analysis of electroosmotic flow of power-law fluids in a slit microchannel," Journal of Colloid and Interface Science, vol. 326 no. 2, pp. 503-510, DOI: 10.1016/j.jcis.2008.06.028, 2008.

[17] C. Zhao, C. Yang, "Nonlinear Smoluchowski velocity for electroosmosis of Power-law fluids over a surface with arbitrary zeta potentials," Electrophoresis, vol. 31 no. 5, pp. 973-979, DOI: 10.1002/elps.200900564, 2010.

[18] C. Zhao, C. Yang, "Electroosmotic flows of non-Newtonian power-law fluids in a cylindrical microchannel," Electrophoresis, vol. 34 no. 5, pp. 662-667, DOI: 10.1002/elps.201200507, 2013.

[19] A. Bandopadhyay, S. Chakraborty, "Steric-effect induced alterations in streaming potential and energy transfer efficiency of non-Newtonian fluids in narrow confinements," Langmuir, vol. 27 no. 19, pp. 12243-12252, DOI: 10.1021/la202273e, 2011.

[20] M. A. Vakili, A. Sadeghi, M. H. Saidi, A. A. Mozafari, "Electrokinetically driven fluidic transport of power-law fluids in rectangular microchannels," Colloids and Surfaces A: Physicochemical and Engineering Aspects, vol. 414, pp. 440-456, DOI: 10.1016/j.colsurfa.2012.07.030, 2012.

[21] C. L. Chaves, J. N. N. Quaresma, E. N. Macedo, L. M. Pereira, J. A. Lima, "Forced convection heat transfer to power-law non-Newtonian fluids inside triangular ducts," Heat Transfer Engineering, vol. 25 no. 7, pp. 23-33, DOI: 10.1080/01457630490495814, 2004.

[22] S. Mukherjee, A. K. Gupta, R. P. Chhabra, "Laminar forced convection in power-law and Bingham plastic fluids in ducts of semi-circular and other cross-sections," International Journal of Heat and Mass Transfer, vol. 104, pp. 112-141, DOI: 10.1016/j.ijheatmasstransfer.2016.08.007, 2017.

[23] M. I. Khan, M. Waqas, T. Hayat, A. Alsaedi, "A comparative study of Casson fluid with homogeneous-heterogeneous reactions," Journal of Colloid and Interface Science, vol. 498, pp. 85-90, DOI: 10.1016/j.jcis.2017.03.024, 2017.

[24] M. A. Vakili, A. Sadeghi, M. H. Saidi, "Pressure effects on electroosmotic flow of power-law fluids in rectangular microchannels," Theoretical and Computational Fluid Dynamics, vol. 28 no. 4, pp. 409-426, DOI: 10.1007/s00162-014-0325-6, 2014.

[25] M. A. Vakili, M. H. Saidi, A. Sadeghi, "Thermal transport characteristics pertinent to electrokinetic flow of power-law fluids in rectangular microchannels," International Journal of Thermal Sciences, vol. 79, pp. 76-89, DOI: 10.1016/j.ijthermalsci.2013.11.016, 2014.

[26] S. Deng, Q. An, M. Li, "The effect of streaming potential and viscous dissipation in the heat transfer characteristics of power-law nanofluid flow in a rectangular microchannel," Micromachines (Basel), vol. 11 no. 4,DOI: 10.3390/mi11040421, 2020.

[27] S. Deng, "Thermally fully developed electroosmotic flow of power-law nanofluid in a rectangular microchannel," Micromachines, vol. 10 no. 6,DOI: 10.3390/mi10060363, 2019.

[28] A. Bejan, "A study of entropy generation in fundamental convective heat transfer," Journal of Heat Transfer, vol. 101 no. 4, pp. 718-725, DOI: 10.1115/1.3451063, 1979.

[29] M. I. Khan, F. Alzahrani, A. Hobiny, Z. Ali, "Fully developed second order velocity slip Darcy-Forchheimer flow by a variable thicked surface of disk with entropy generation," Int. Commun. Heat Mass Transfer, vol. 117,DOI: 10.1016/j.icheatmasstransfer.2020.104778, 2020.

[30] J. Wang, R. Muhammad, M. I. Khan, W. A. Khan, S. Z. Abbas, "Entropy optimized MHD nanomaterial flow subject to variable thicked surface," Computer Methods and Programs in Biomedicine, vol. 189,DOI: 10.1016/j.cmpb.2019.105311, 2020.

[31] R. Muhammad, M. I. Khan, N. B. Khan, M. Jameel, "Magnetohydrodynamics (MHD) radiated nanomaterial viscous material flow by a curved surface with second order slip and entropy generation," Computer Methods and Programs in Biomedicine, vol. 189,DOI: 10.1016/j.cmpb.2019.105294, 2020.

[32] R. Muhammad, M. I. Khan, M. Jameel, N. B. Khan, "Fully developed Darcy-Forchheimer mixed convective flow over a curved surface with activation energy and entropy generation," Computer Methods and Programs in Biomedicine, vol. 188,DOI: 10.1016/j.cmpb.2019.105298, 2020.

[33] X. Y. Chen, K. C. Toh, J. C. Chai, C. Yang, "Developing pressure-driven liquid flow in microchannels under the electrokinetic effect," International Journal of Engineering Science, vol. 42 no. 5-6, pp. 609-622, DOI: 10.1016/j.ijengsci.2003.07.008, 2004.

[34] H. A. Stone, A. D. Stroock, A. Ajdari, "Engineering flows in small devices," Annual Review of Fluid Mechanics, vol. 36 no. 1, pp. 381-411, DOI: 10.1146/annurev.fluid.36.050802.122124, 2004.

[35] D. Maynes, B. W. Webb, "The effect of viscous dissipation in thermally fully-developed electro-osmotic heat transfer in microchannels," International Journal of Heat and Mass Transfer, vol. 47 no. 5, pp. 987-999, DOI: 10.1016/j.ijheatmasstransfer.2003.08.016, 2004.

[36] R. K. Shukla, X. Zhong, "Derivation of high-order compact finite difference schemes for non-uniform grid using polynomial interpolation," Journal of Computational Physics, vol. 204 no. 2, pp. 404-429, DOI: 10.1016/j.jcp.2004.10.014, 2005.

[37] C.-Y. Wang, C.-F. Kung, C.-C. Chang, "Approach to analytic solutions for electroosmotic flow in micro-ducts by eigenfunctions of the Helmholtz equation," Microfluid Nanofluid, vol. 20,DOI: 10.1007/s10404-016-1764-8, 2016.

[38] A. C. Pipkin, R. S. Rivlin, "Normal stresses in flow through tubes of non-circular cross-section," Zeitschrift für angewandte Mathematik und Physik, vol. 14 no. 6, pp. 738-742, DOI: 10.1007/bf01601165, 1963.

[39] A. K. Nayak, A. Haque, B. Weigand, W. Steve, "Thermokinetic transport of dilatant/pseudoplastic fluids in a hydrophobicpatterned micro-slit," Physics of Fluids, vol. 32,DOI: 10.1063/5.0011716, 2020.

[40] B. Z. Lu, Y. C. Zhou, M. J. Holst, J. A. Mc Cammon, "Recent Progress in Numerical Methods for the PoissonBoltzmann Equation in Biophysical Applications," Communication in Computational Physics, vol. 3 no. 5, pp. 973-1009, 2008.