Unsteady Pressure-Driven Electrokinetic Slip Flow and Heat Transfer of Power-Law Fluid through a Microannulus

To guarantee the transporting efficiency of microdevices associated with fluid transportation, mixing, or separation and to promote the heat transfer performance of heat exchangers in microelectronics, the hydrodynamic behaviors at unsteady and steady states, as well as the thermal characteristics at the steady state in a pressure-driven electrokinetic slip flow of power-law fluid in a microannulus are studied. To present a more reliable prediction, the slip phenomenon at walls and nonlinear rheology of liquid are incorporated. The modified Cauchy momentum equation applicable to all time scales and energy equations, are analytically solved in the limiting case of a Newtonian fluid and numerically solved for power-law fluids. The transient velocity profile, time evolution of flow rate, temperature profile, and heat transfer rate are computed at different flow behavior indices, electrokinetic width, slip lengths, and Brinkman numbers, thereby, the coupling effect of nonlinear rheology, slip hydrodynamics, and annular geometry on flow and thermal behaviors is explored. The unsteady flow takes a longer time to achieve the steady state for shear thinning fluids or greater slip lengths. The flow behavior index and slip length play a significant role in the flow rate and heat transfer performance. The relevant discussion can serve as a theoretical guide for the operation and thermal management of annular geometry-related flow actuation systems.


Introduction
In microflows, the surface effect predominates, and the contact between the charged surface of the channel wall and ionic liquid leads to the redistribution of nearby ions, inducing the electric double layer (EDL) [1]. The application of an external electric field tangential to the liquid in the microchannel results in the migration of excess counter ions in the EDL, finally leading to electroosmotic flow (EOF) under the viscous drag force, which is generally termed electrokinetic flow. With the development of the micro-electronicmechanical system (MEMS), and due to its desirable attributes, electroosmosis is used as a flow actuation mechanism in chemical and biomedical analysis, membrane separation, soil remediation, and the thermal management of microelectronic systems [2,3].
To meet the growing demand for electroosmosis actuation systems, a vast majority of the literature has been reported from different aspects, such as the influence of different geometries on EOF and theoretical surveys on two-layer electroosmosis systems [4][5][6][7]. Among them, the electrokinetic flow in an annular geometry remains an extensive topic of scientific and technological interest. From a scientific point of view, the cylindrical (rectangular) annulus model is of more universal significance, and can be treated as a parallel plate or circular cylinder (rectangular microchannel) in the limiting cases of the ratio of the inner radius to the outer radius [8,9], and can also serve as a basis for the development of two-layer fluid systems [9,10]. From a practical point of view, in studying the chemical remediation of contaminated soil, treating the pores inside the soils as microannulus rather To the best of the authors' knowledge, there is still a lack of comprehensive investigation on the unsteady slippery EOF of power-law fluids through an annulus. However, the corresponding analysis not only reveals the underlying mechanism that how the fluid is initiated but is also of practical use in the operation of bio-MEMS and remediation of contaminated soil. And the search for effective ways to drive various fluids into microchannels of different geometries persists. Bearing this in mind, this paper aims to explore the slip hydrodynamics at unsteady and steady states, as well as thermal characteristics in the pressure-driven electrokinetic flow of power-law fluids through a microannulus. The fundamental understanding of this phenomenon is crucial for guaranteeing the transporting efficiency of microdevices associated with fluid transportation, mixing or separation, and promoting the heat transfer performance of heat exchangers in microelectronics. For the sake of generality, the flow of the power-law fluid is driven by the pressure gradient and electric field. In addition, the slip phenomenon at the walls is incorporated in mathematical formulation to obtain a more realistic microflow model. The novelty of the present work is the investigation of the flow mechanism for all time scales under the nonlinear coupling of power-law rheology, slippery interfacial hydrodynamics, and annular geometrical effect, as well as the evaluation of the heat transfer performance of a slip flow system under the combined effect of viscous dissipation and Joule heating. The governing equations are analytically solved in the limiting case of a Newtonian fluid and are numerically solved for power-law fluids. The velocity and temperature distributions, flow rate, and heat transfer rate are evaluated at different pertinent parameters.

Modified Cauchy Momentum Equation
The continuity equation and Cauchy momentum equation for an incompressible laminar flow are given as [4].
where V is the velocity vector, p * is the pressure, τ is the stress for power-law fluid, F is the body force acting on liquid, and ρ is the liquid density. Here, the flow is driven under the combined effect of electroosmotic force arising from the applied electric field, pressure gradient, and slip hydrodynamics at the walls, as sketched in Figure 1. The microannulus is characterized by the inner radius r 1 * and outer radius r 2 * in which the slip lengths near the inner and outer walls are l * 1 and l * 2 , respectively, and the strength of the electric field is E 0 . To facilitate the analysis and capture the primary physics of slip flow, the following assumptions are applied: Equation (1) the thickness of EDLs at two walls is far less than r i * , and thus the EDLs will not overlap, the zeta potentials near the inner and outer walls ζ i * are constant; Equation (2) in microchannels, the gravity of fluid is negligible; Equation (3) since the length of the microchannel is far longer than the characteristic radius of microannulus, the radial velocity can be ignored and only the axial velocity v * is considered; and Equation (4) the slip phenomenon is represented by Navier's slip model [36]. Applying the assumptions above to Equations (1) and (2), one has the modified Cauchy momentum equation, the slip conditions at the channel walls, and the initial condition [36,37].  In Equation (3), η * = mn|∂v * /∂r * | n−1 represents the apparent viscosity of the power-law fluid where n denotes the flow behavior index of the unit [1], n < 1, n = 0, and n > 1 physically represent the shear thinning, Newtonian, and shear thickening fluids, respectively, and mn denotes the flow consistency index of the unit [Nm −2 s n ] which reduces to the dynamic viscosity of the Newtonian fluid, namely, m0 of the unit [Nm −2 s], when n = 1. In Equation (3), the local net charge density is subject to Boltzmann distribution based on the assumption of a local thermodynamic equilibrium ρe = −2ezn0 sinh[(ezψ * )/(kBTa)]. The electric potential distribution ψ * for a symmetric electrolyte due to the presence of the EDL is determined by the well-known Poisson-Boltzmann (P-B) equation [1,4].
where e denotes the elementary charge, z denotes the ion valence, n0 is the ion concentration in liquid, kB is the Boltzmann number, and Ta is the absolute temperature. In Equation (4), l * 1 and l * 2 are the property parameters of the solid wall and working liquid, representing the constant slip lengths near the inner and outer walls, respectively. Physically, l * i indicates the thickness beyond the solid-liquid interface where the velocity extrapolates to zero, as shown in Figure 1. According to the Debye-Hückel linear approximation [1], the hyperbolic term sinh[(ezψ * )/(kBTa)] can be linearized as [(ezψ * )/(kBTa)] by expanding the hyperbolic function up to first order, and it can work well when [(ezψ * )/(kBTa)] < 1, physically for small zeta potentials, as compared to the thermal energy of ions, kBTa/e < 25.7 mV [4]. Consequently, the electric potential distribution is governed by the following linearized P-B equation [15].
Introducing the dimensionless variables r = r * /R, ψ * = zeψ * /(kBTa), v = v * /Vhs, and t = t * m0/(ρR 2 ) to Equations (3), (4), (6) and (7) produces the dimensionless governing equations [36,37]. In Equation (3), η * = m n |∂v * /∂r * | n−1 represents the apparent viscosity of the powerlaw fluid where n denotes the flow behavior index of the unit [1], n < 1, n = 0, and n > 1 physically represent the shear thinning, Newtonian, and shear thickening fluids, respectively, and m n denotes the flow consistency index of the unit [Nm −2 s n ] which reduces to the dynamic viscosity of the Newtonian fluid, namely, m 0 of the unit [Nm −2 s], when n = 1. In Equation (3), the local net charge density is subject to Boltzmann distribution based on the assumption of a local thermodynamic equilibrium ρ e = −2ezn 0 sinh[(ezψ * )/(k B T a )]. The electric potential distribution ψ * for a symmetric electrolyte due to the presence of the EDL is determined by the well-known Poisson-Boltzmann (P-B) equation [1,4].
where e denotes the elementary charge, z denotes the ion valence, n 0 is the ion concentration in liquid, k B is the Boltzmann number, and T a is the absolute temperature. In Equation (4), l * 1 and l * 2 are the property parameters of the solid wall and working liquid, representing the constant slip lengths near the inner and outer walls, respectively. Physically, l * i indicates the thickness beyond the solid-liquid interface where the velocity extrapolates to zero, as shown in Figure 1. According to the Debye-Hückel linear approximation [1], the hyperbolic term sinh[(ezψ * )/(k B T a )] can be linearized as [(ezψ * )/(k B T a )] by expanding the hyperbolic function up to first order, and it can work well when [(ezψ * )/(k B T a )] < 1, physically for small zeta potentials, as compared to the thermal energy of ions, k B T a /e < 25.7 mV [4]. Consequently, the electric potential distribution is governed by the following linearized P-B equation [15].
Using Equations (8) and (9), one has the dimensionless unsteady flow rate

Energy Equation
The temperature distribution T * for the thermally fully developed steady flow is governed by the energy equation [26] ( where c p denotes the specific heat at constant pressure, k is the thermal conductivity, σ is the electric conductivity of the liquid, v s * is the velocity distribution of steady flow, q s denotes the constant wall heat flux, and T * w denotes the wall temperature. Further, since the constant heat flux boundary condition q s = const is adopted, one has ∂[(T * w − T * )/(T * w − T * m )]/∂r * = 0, thereby, ∂T * /∂r * = dT * w /dr * = dT * m /dr * ≡const in which T * m denotes the mean temperature. Imposing the overall energy balance condition over an elemental control volume yields [26] dT * m dr * = 2πr * 1 q s + 2πr * 2 q s + σE 0 2 π(r * 2 2 − r * 2 where v * sm indicates the mean temperature of steady flow. Introducing the dimensionless group T = k(T * − T * w )/(q s R), v s = v s * /V hs , F 1 = r 2 r 1 v s rdr, F 2 = r 2 r 1 |∂v s /∂r| n−1 (∂v s /∂r) 2 rdr, S = σE 2 0 R/q s and Br = m n V 2 hs /(q s R)(V hs /R) n−1 into Equations (13) and (14), and combining with Equation (15), one has the dimensionless energy equations where the first term in the left-hand side (LHS) of Equation (16) denotes the heat generation caused by heat diffusion, and the remaining terms in LHS denote the heat generations arising from axial conduction, Joule heating, and viscous dissipation, respectively. To note, S is the Joule heating parameter and Br is the Brinkman number, representing the viscous dissipation effect. To provide an in-depth insight into the heat transfer performance of electrokinetic slip flow, the first law analysis is conducted, i.e., the Nusselt number implying the heat transfer efficiency is deduced [26,38]. According to Nu = 2q s (r * 2 − r * 1 )/(T * w − T * m )/k and the mean temperature one has the Nusselt number for the slip flow in a microannulus

Solution Method and Validation
First, P-B Equations (10) and (11) are solved where I 0 and K 0 imply the 0-th order modified Bessel functions of the first kind and second kind, respectively, and the coefficients are given as A =

For Newtonian Fluid (n = 1) and without Consideration of Viscous Dissipation (Br = 0)
In the limiting case of a Newtonian fluid, Equations (8) and (9) where v N is the velocity of a Newtonian fluid. Since it is inhomogeneous, Equation (21) is homogenized as Using the method of variable separation, Equations (23) and (24) The coefficient C m and the basis function M 0 (λ m r) are presented in the Appendix A for conciseness and readability. Applying Duhamel's principle and combining Equation (25), the solution to Equations (21) and (22) As the slip flow evolves to the steady state, the temporal term in Equation (21) vanishes and the steady velocity is solved from Equations (21) and (22) as where the coefficients D and E are provided in the Appendix A. The first and second terms in the RHS of Equation (26) represent the steady velocity and the transient one, respectively. When t→∞, the second term in Equation (26) vanishes and thus, v N →v s N .
Therefore, in evaluating the unsteady velocity, to eliminate the oscillation caused by a truncated error in series solution (26), the first term in Equation (26) is replaced with v s N , given by Equation (27), and the final version of velocity for unsteady slip flow reads Similarly, as t→∞, v N automatically becomes v s N . With v s N obtained, for the Newtonian fluid, Equations (16) and (17) are analytically solved in the absence of viscous dissipation (Br = 0) where the coefficients F, C 0 , T 0 , and the intermediate function Q 0 (r) are presented in the Appendix A. The symbolic simplification of several coefficients is carried out with the help of Maple, and the computation and graphical presentation of Equations (28) and (29) are carried out using Matlab.

For Power-Law Fluids (n = 1) and with Consideration of Viscous Dissipation (Br = 0)
For power-law fluids and when considering the viscous dissipation effect, the momentum Equations (8) and (9) and the energy Equations (16) and (17) are solved using the finite difference method. Let t l = l∆t, r j = j∆r, v l s,j = v s (j∆r,l∆t), v l j = v(j∆r,l∆t), ψ j = ψ(j∆r), and T l j = T(j∆r,l∆t), l = 1,2, . . . ,L and j = 1,2, . . . ,J. The explicit finite difference scheme is used for time and the central difference scheme is adopted for space. According to Equation (9), the numerical velocity at the boundaries is v l . The bulk liquid velocity is computed by following the numerical algorithm v l+1 , j = 2,3, . . . ,J -1. When t→∞, the transient velocity approaches steady velocity v l s , i.e., v l − v l+1 < err with err being a specified criterion. Consequently, the flow rates can be computed from Equation (12) by the composite trapezoidal integration method. The temporal term ∂T/∂t has been introduced in the RHS of Equation (16), which vanishes when t→∞, and T becomes the fully developed temperature. The numerical velocity and temperature are given as vectors . The Crank-Nicolson scheme is used for time and the central difference scheme is adopted for space. Combining with the boundary conditions (17), the numerical algorithm is Here, Λ = tridiag(I,A, means a sparse matrix denoting that, except for x p,q = u, the remaining elements equal zero. The inhomogeneous term in Equation (13) is discretized as With the numerical temperature T 1 obtained, the mean temperature and Nusselt number can be computed from Equations (15) and (16) by using the composite trapezoidal integration method.
In Figure 2a, in the limiting case of a Newtonian fluid, the numerical velocity computed based on Equation (30) is compared with the analytical velocity obtained from Equation (28) and the existing data presented in the study [28] by Yavari et al. and [13] by Tsao. In Figure 2b, for power-law fluids, the numerical velocities computed from Equation (30) are compared with the existing data from [30] by Shamshiri et al. In Figure 2c, the average temperature T m is numerically computed for different flow behavior indices n and grid numbers in terms of space J. Irrespective of the value of n, when c = 0.016 and J ≥ 151, the curve of T m with J shows little change, meaning that J ≥ 151 is enough to obtain steady numerical solutions. To be economical, the grid number J is chosen as 201. Figure 2a,b shows when J = 201, the numerical velocities agree well with the analytical velocity or the existing data, meaning that the numerical algorithm is feasible.

Results
The unsteady hydrodynamics and heat transfer in the electrokinetic slip flow of power-law fluids in a microannulus are investigated by evaluating the unsteady velocity field, flow rate, fully developed temperature field, and Nusselt number at different parameters. The physical parameters take the given values in Nomenclature based on practical uses [4][5][6]39]. To note, ζi * = −0.025V with i = 1,2, and the reason for assuming a small zeta potential is to use the Debye-Hückel linearization and the fact that the small zeta potential is physically acceptable [1]. To obtain more realistic predictions, it is essential to Figure 2. (a) The comparison among the analytical velocity, numerical velocity, and existing data from the studies [13,28] without consideration of slip velocity at the channel walls, (b) the comparison between numerical velocity and existing data in [30] different n, and (c) the grid independence study of the mean temperature for shear thinning fluid and shear thickening fluid when Br = 0.005, S = 2, and ζ 1 = ζ 2 = −1.

Results
The unsteady hydrodynamics and heat transfer in the electrokinetic slip flow of powerlaw fluids in a microannulus are investigated by evaluating the unsteady velocity field, flow rate, fully developed temperature field, and Nusselt number at different parameters. The physical parameters take the given values in Nomenclature based on practical uses [4][5][6]39]. To note, ζ i * = −0.025V with i = 1,2, and the reason for assuming a small zeta potential is to use the Debye-Hückel linearization and the fact that the small zeta potential is physically acceptable [1]. To obtain more realistic predictions, it is essential to present the permissible ranges of governing parameters. Based on the well-established electroosmosis theory of power-law fluid, the flow behavior index n ranges from 0.6 to 1.4 [23] and the electrokinetic width W ranges from 10 to 100 [4,23]. The order of the Brinkman number Br can be of O(10 −2 ) since the orders of reference velocity, apparent viscosity, and reference radius are identical to that from the studies [7,39]. According to the study [36], the slip lengths l i range from 0 to 0.1 and the dimensionless pressure gradient ∇p is of the order O(1), which is chosen as 5. Figure 3 presents the unsteady velocity distribution over a cross-sectioned area of a microannulus at t = 0.006, t = 0.06, and t = 0.6 for the shear thinning, Newtonian, and shear thickening fluids. From Figure 3a,d,g, irrespective of the fluid type, at first the fluid at and near the channel walls is set in motion under the electroosmotic force and slippery effect, which then drags the bulk fluid in the central area forward under the shear stress and pressure gradient, forming the pressure driven electrokinetic slip flow. Figure 3a-c shows that for the shear thinning fluid, as time evolves from 0.006 to 0.6, the velocity grows obviously, and the magnitude of velocity eventually goes beyond 3. In contrast, Figure 3g-i shows that for the shear thickening fluid, from t = 0.006 to t = 0.06, the velocity distribution evolves, however, it changes little when time lapses from t = 0.06 to t = 0.6. The comparison among Figure 3a-c,d-f,g-i indicates that as time evolves, the shear thinning fluid attains the greatest velocity which is two times more than that of Newtonian and shear thickening fluids. Moreover, as time goes beyond 0.06, the change in velocity becomes smaller when the fluid changes from a shear thinning to a shear thickening fluid, meaning that the shear thickening fluid achieves the steady state earlier than the shear thinning and Newtonian fluids. In addition, since l i = 0.01 (i = 1, 2), the fluid at the channel walls is not static and slips over the solid walls.
As shown in Figure 4, the variation of the velocity profile with slip length ratio l r (l 2 /l 1 ) at different times is presented for the shear thinning, Newtonian, and shear thickening fluids. In terms of the influence of the slip length ratio, from Figure 4a,d,g,j,m, no matter what type of fluid is considered, since the slip flow is initiated close to the channel walls at first and l 1 is fixed, the flow near the channel walls is accelerated with the increase in l r . As time evolves, with the bulk liquid set in motion, the influence of the change in the slip length ratio extends from the outer channel wall to the bulk liquid. When l r = 1, the velocity shows a nearly symmetric profile, while the asymmetric velocity profile is observed when l r = 1. To note, when l r = 0, the fluid at the outer channel wall is static all the time, corresponding to the no-slip boundary condition. In addition, different from that under the no-slip condition, the fluid under the slip condition is not subject to resistance at the walls and slips over the solid walls, leading to a higher velocity for the bulk liquid. In terms of the influence of fluid type, namely, n, as shown in Figure 4a shows that for the shear thickening fluid, from t = 0.006 to t = 0.06, the velocity distribution evolves, however, it changes little when time lapses from t = 0.06 to t = 0.6. The comparison among Figure 3a-c,d-f,g-i indicates that as time evolves, the shear thinning fluid attains the greatest velocity which is two times more than that of Newtonian and shear thickening fluids. Moreover, as time goes beyond 0.06, the change in velocity becomes smaller when the fluid changes from a shear thinning to a shear thickening fluid, meaning that the shear thickening fluid achieves the steady state earlier than the shear thinning and Newtonian fluids. In addition, since li = 0.01 (i = 1, 2), the fluid at the channel walls is not static and slips over the solid walls. As shown in Figure 4, the variation of the velocity profile with slip length ratio lr (l2/l1) at different times is presented for the shear thinning, Newtonian, and shear thickening fluids. In terms of the influence of the slip length ratio, from Figure 4a  As time evolves, with the bulk liquid set in motion, the influence of the change in the slip length ratio extends from the outer channel wall to the bulk liquid. When lr = 1, the velocity shows a nearly symmetric profile, while the asymmetric velocity profile is observed when lr ≠ 1. To note, when lr = 0, the fluid at the outer channel wall is static all the time, corresponding to the no-slip boundary condition. In addition, different from that under the noslip condition, the fluid under the slip condition is not subject to resistance at the walls and slips over the solid walls, leading to a higher velocity for the bulk liquid. In terms of the influence of fluid type, namely, n, as shown in Figure 4a   In Figure 5, the variation of the velocity profile with electrokinetic width W at a different time is presented for the shear thinning, Newtonian, and shear thickening fluids. When t = 0.006, from Figure 5a,d,g,j,m, since li ≠ 0, the fluid at and near the channel walls attains velocity first, increasing with W and decreasing with n. When the time increases to 0.06, as shown in Figure 5b,e,h,k,n, the fluid near the channel walls drives the fluid in the central area to move forward. The velocity profile of the shear thinning fluid is still developing and the shear thickening fluid flow approaches the developed state. This is especially evident with the greater value of electrokinetic width W. From Figure 5c,f,i,l,o, the velocity profile becomes steady, which shows a plug-like pattern when the electrokinetic width changes from 10 to 100 or when the fluid changes from shear thickening to shear thinning. Compared to the Newtonian and shear thickening fluids, the slip flow of a shear thinning fluid is much more sensitive to the change in electrokinetic width W. In Figure 5, the variation of the velocity profile with electrokinetic width W at a different time is presented for the shear thinning, Newtonian, and shear thickening fluids. When t = 0.006, from Figure 5a,d,g,j,m, since l i = 0, the fluid at and near the channel walls attains velocity first, increasing with W and decreasing with n. When the time increases to 0.06, as shown in Figure 5b,e,h,k,n, the fluid near the channel walls drives the fluid in the central area to move forward. The velocity profile of the shear thinning fluid is still developing and the shear thickening fluid flow approaches the developed state. This is especially evident with the greater value of electrokinetic width W. From Figure 5c,f,i,l,o, the velocity profile becomes steady, which shows a plug-like pattern when the electrokinetic width changes from 10 to 100 or when the fluid changes from shear thickening to shear thinning. Compared to the Newtonian and shear thickening fluids, the slip flow of a shear thinning fluid is much more sensitive to the change in electrokinetic width W.  In Figure 5, the variation of the velocity profile with electrokinetic width W at a different time is presented for the shear thinning, Newtonian, and shear thickening fluids. When t = 0.006, from Figure 5a,d,g,j,m, since li ≠ 0, the fluid at and near the channel walls attains velocity first, increasing with W and decreasing with n. When the time increases to 0.06, as shown in Figure 5b,e,h,k,n, the fluid near the channel walls drives the fluid in the central area to move forward. The velocity profile of the shear thinning fluid is still developing and the shear thickening fluid flow approaches the developed state. This is especially evident with the greater value of electrokinetic width W. From Figure 5c,f,i,l,o, the velocity profile becomes steady, which shows a plug-like pattern when the electrokinetic width changes from 10 to 100 or when the fluid changes from shear thickening to shear thinning. Compared to the Newtonian and shear thickening fluids, the slip flow of a shear thinning fluid is much more sensitive to the change in electrokinetic width W.     Figure 7 shows that no matter what type of fluid is considered, with the increase in the slip length ratio lr, the left side of the temperature profile augments, and the right side of the temperature profile shows a slight decrement, in the meantime, the minimum value of temperature shifts from  Figure 6 illustrates the time evolution of the unsteady flow rate ratio for different flow behavior indices n at (a) l 1 = l 2 = 0.01 and (b) l 1 = l 2 = 0.05. The unsteady flow rate ratio is taken as the ratio of the unsteady flow rate to the respective steady flow rate to capture the temporal physical picture of the unsteady slip flow for different types of fluids. As shown in Figure 6a, when n ranges from 0.6 to 1.4, the fluid changes from shear thinning to shear thickening, and the slip flow achieves the steady state earlier and earlier. This is due to the fact that the higher apparent viscosity of the shear thickening fluid results in a stronger resistance to the flow and enables the unsteady motion to arrive at the steady state faster. Compared to Figure 6a,b shows that the increase in slip lengths l i intensifies the driving force near the channel walls; as a result, the unsteady flow needs a longer time to arrive at the steady state.  Figure 6 illustrates the time evolution of the unsteady flow rate ratio for different flow behavior indices n at (a) l1 = l2 = 0.01 and (b) l1 = l2 = 0.05. The unsteady flow rate ratio is taken as the ratio of the unsteady flow rate to the respective steady flow rate to capture the temporal physical picture of the unsteady slip flow for different types of fluids. As shown in Figure 6a, when n ranges from 0.6 to 1.4, the fluid changes from shear thinning to shear thickening, and the slip flow achieves the steady state earlier and earlier. This is due to the fact that the higher apparent viscosity of the shear thickening fluid results in a stronger resistance to the flow and enables the unsteady motion to arrive at the steady state faster. Compared to Figure 6a,b shows that the increase in slip lengths li intensifies the driving force near the channel walls; as a result, the unsteady flow needs a longer time to arrive at the steady state.   Figure 7 shows that no matter what type of fluid is considered, with the increase in the slip length ratio lr, the left side of the temperature profile augments, and the right side of the temperature profile shows a slight decrement, in the meantime, the minimum value of temperature shifts from  Figure 7 shows the variation of the temperature profile with slip length ratio l r for the shear thinning fluid, Newtonian fluid, and shear thickening fluid. Figure 7 shows that no matter what type of fluid is considered, with the increase in the slip length ratio l r , the left side of the temperature profile augments, and the right side of the temperature profile shows a slight decrement, in the meantime, the minimum value of temperature shifts from left to right. In addition, in the limiting case of the no-slip condition at one channel wall (l r = 0), the temperature difference between the bulk liquid and channel walls is the widest, decreasing when the slip length ratio is enhanced, meaning that the presence of slip velocity at solid walls promotes the heat transfer of the electrokinetic slip flow. Furthermore, the influence of n and that of l r will not interact with each other. left to right. In addition, in the limiting case of the no-slip condition at one channel wall (lr = 0), the temperature difference between the bulk liquid and channel walls is the widest, decreasing when the slip length ratio is enhanced, meaning that the presence of slip velocity at solid walls promotes the heat transfer of the electrokinetic slip flow. Furthermore, the influence of n and that of lr will not interact with each other.  Figure 8 shows the variation of the temperature profile with the slip length ratio lr for the shear thinning fluid, Newtonian fluid, and shear thickening fluid. It is obvious that the temperature profile decreases; namely, the temperature difference is widened with the increase in Brinkman number Br. The greater Brinkman number implies a stronger viscous dissipation effect which inevitably hinders the heat transfer of the electrokinetic flow; as a result, the widened temperature difference is observed. When the fluid changes from shear thinning to shear thickening, the influence of the Brinkman number Br on the temperature profile is reduced, since the corresponding velocity gradient of electrokinetic flow becomes smaller.   Figure 8 shows the variation of the temperature profile with the slip length ratio l r for the shear thinning fluid, Newtonian fluid, and shear thickening fluid. It is obvious that the temperature profile decreases; namely, the temperature difference is widened with the increase in Brinkman number Br. The greater Brinkman number implies a stronger viscous dissipation effect which inevitably hinders the heat transfer of the electrokinetic flow; as a result, the widened temperature difference is observed. When the fluid changes from shear thinning to shear thickening, the influence of the Brinkman number Br on the temperature profile is reduced, since the corresponding velocity gradient of electrokinetic flow becomes smaller. left to right. In addition, in the limiting case of the no-slip condition at one channel wall (lr = 0), the temperature difference between the bulk liquid and channel walls is the widest, decreasing when the slip length ratio is enhanced, meaning that the presence of slip velocity at solid walls promotes the heat transfer of the electrokinetic slip flow. Furthermore, the influence of n and that of lr will not interact with each other.  Figure 8 shows the variation of the temperature profile with the slip length ratio lr for the shear thinning fluid, Newtonian fluid, and shear thickening fluid. It is obvious that the temperature profile decreases; namely, the temperature difference is widened with the increase in Brinkman number Br. The greater Brinkman number implies a stronger viscous dissipation effect which inevitably hinders the heat transfer of the electrokinetic flow; as a result, the widened temperature difference is observed. When the fluid changes from shear thinning to shear thickening, the influence of the Brinkman number Br on the temperature profile is reduced, since the corresponding velocity gradient of electrokinetic flow becomes smaller.  The variation in the Nusselt number Nu with the flow behavior index n is presented in Figure 9 for different slip lengths l i . When considering the shear thinning fluid (n < 1), the Nusselt number Nu augments with the flow behavior index n, and the increasing rate is enhanced with slip lengths l i . When considering the shear thickening fluid (n > 1), the Nusselt number Nu shows little change with n, though it still increases with l i , meaning that the influence of n on the Nusselt number Nu becomes smaller; in other words, Nu is not as susceptible as that of the shear thinning fluid to the change in the flow behavior index n. The curve of Nu with n increases as a whole with the slip lengths l i . Although the flows of both shear thinning and shear thickening fluids are accelerated by the greater slip lengths, the velocity distribution becomes more uniform for shear thickening fluids; in contrast, a wider velocity difference between the channel walls and central region occurs for shear thinning fluids, which partly suppresses the heat transfer of the bulk liquid. Therefore, the increase in Nu with slip lengths is especially significant for shear thickening fluids. Consequently, in engineering practice, when driving the shear thickening fluids in the microannulus, the heat transfer performance can be promoted by adjusting the slip length. The variation in the Nusselt number Nu with the flow behavior index n is presented in Figure 9 for different slip lengths li. When considering the shear thinning fluid (n < 1), the Nusselt number Nu augments with the flow behavior index n, and the increasing rate is enhanced with slip lengths li. When considering the shear thickening fluid (n > 1), the Nusselt number Nu shows little change with n, though it still increases with li, meaning that the influence of n on the Nusselt number Nu becomes smaller; in other words, Nu is not as susceptible as that of the shear thinning fluid to the change in the flow behavior index n. The curve of Nu with n increases as a whole with the slip lengths li. Although the flows of both shear thinning and shear thickening fluids are accelerated by the greater slip lengths, the velocity distribution becomes more uniform for shear thickening fluids; in contrast, a wider velocity difference between the channel walls and central region occurs for shear thinning fluids, which partly suppresses the heat transfer of the bulk liquid. Therefore, the increase in Nu with slip lengths is especially significant for shear thickening fluids. Consequently, in engineering practice, when driving the shear thickening fluids in the microannulus, the heat transfer performance can be promoted by adjusting the slip length.

Conclusions
The temporal physical picture of the unsteady pressure-driven electrokinetic slip flow for power-law fluids in a microannulus is provided, based on which, the heat transfer of a steady flow is analyzed for different governing parameters including the flow behavior index (n), slip lengths (li), electrokinetic width (W), and Brinkman number (Br). In terms of the temporal hydrodynamical behavior of power-law fluids, the flow is accelerated with the increase in slip lengths and electrokinetic width and the flow of shear thinning fluids is greater and much more sensitive to the change of the above parameters than that of Newtonian and shear thickening fluids. When the fluid changes from shear thinning to shear thickening or the slip lengths are reduced, the flow reaches the steady state earlier. In terms of the thermal behavior of the steady flow, the temperature profile increases with the slip length ratio and flow behavior index, which decreases with the Brinkman number, meaning that the more uniform the velocity distribution, the narrower the temperature difference between the channel walls and bulk liquid, and the more intense the heat transfer performance. Therefore, in practical uses, the slip surface can be reliably engineered to achieve a higher flow rate and promote heat transfer. The analysis above can serve as a theoretical guidance for the design and optimum operation of relevant microfluidic devices.

Conclusions
The temporal physical picture of the unsteady pressure-driven electrokinetic slip flow for power-law fluids in a microannulus is provided, based on which, the heat transfer of a steady flow is analyzed for different governing parameters including the flow behavior index (n), slip lengths (l i ), electrokinetic width (W), and Brinkman number (Br). In terms of the temporal hydrodynamical behavior of power-law fluids, the flow is accelerated with the increase in slip lengths and electrokinetic width and the flow of shear thinning fluids is greater and much more sensitive to the change of the above parameters than that of Newtonian and shear thickening fluids. When the fluid changes from shear thinning to shear thickening or the slip lengths are reduced, the flow reaches the steady state earlier. In terms of the thermal behavior of the steady flow, the temperature profile increases with the slip length ratio and flow behavior index, which decreases with the Brinkman number, meaning that the more uniform the velocity distribution, the narrower the temperature difference between the channel walls and bulk liquid, and the more intense the heat transfer performance. Therefore, in practical uses, the slip surface can be reliably engineered to achieve a higher flow rate and promote heat transfer. The analysis above can serve as a theoretical guidance for the design and optimum operation of relevant microfluidic devices.   with J ν and Y ν (i = 0,1) denoting the ν-th order Bessel functions of the first kind and second kind. The coefficients in Equation (27) are given as D = ζ 2 + ∇pr 2 2 /4 + l 1 [ψ (r 1 ) + ∇pr 2 /2] − ζ 1 − ∇pr 2 1 /4 + l 2 [ψ (r 2 ) + ∇pr 2 /2] ln r 1 − l 2 /r 2 − ln r 2 − l 1 /r 1 (A13) E = l 1 [ψ (r 1 ) + ∇pr 1 /2 + D/r 1 ] − (ζ 1 + ∇pr 2 1 /4 + D ln r 1 ) The coefficients and intermediate functions P(r) and Q 0 (r) in Equation (29) are C 0 = P(r 2 ) − P(r 1 ) + Q 0 (r 1 )(ln r 1 − ln r 2 ) ln r 1 + ln r 2 (A16) F = r 1 + r 2 + S(r 2 2 − r 2 1 ) P 0 (r 2 ) − P 0 (r 1 ) (2E − D)r 2 + D 2 r 2 ln r (A20)