Thermal and Flow Analysis of Fully Developed Electroosmotic Flow in Parallel-Plate Micro- and Nanochannels with Surface Charge-Dependent Slip

This study analytically investigates the coupled effects of surface charge and boundary slip on the fully developed electroosmotic flow and thermal transfer in parallel plate micro and nanochannels under the high zeta potential. The electric potential, velocity, temperature, flow rate, and Nusselt number are obtained analytically. The main results are that the velocity of bulk flow is significantly reduced in the presence of the surface charge-dependent slip. Moreover, the maximum velocity at ζ = −125 mV is approximately twice as large as that at ζ = −25 mV. The velocity and dimensionless temperature increase as the zeta potential increases. The dimensionless temperature of the surface charge-dependent slip flow is larger than that of the surface charge-independent slip flow. For the surface charge-dependent slip flow, the maximum temperature at ζ = −125 mV is approximately four times larger than that at ζ = −25 mV. The Nusselt number decreases with Joule heating and increases with a positive heat transfer coefficient. The Nusselt number decreases as the electric field and the magnitude of the zeta potential increase. In the surface charge-dependent slip flows, the Nusselt number is smaller than that in the surface charge-independent slip flows.


Introduction
Over the last decades, the micro-and nanoscale fluid systems based on physics, medicine, chemistry and biology have been widely applied in various fields [1,2], such as the separation and detection of chemical and biological samples, the design of heat and mass transfer system, micro-and nanoflow control in lab-on-chip [2]. In these applications, the thermal properties of fluid flows are significant factors affecting the performance of micro-nanoscale fluid systems, including the temperature of fluid could affect the properties of biological and medical samples [3], as well as the formation and stability of nanobubbles at the solid-liquid interface [4]. Traditionally, compared with macroscopic flows, the reduced characteristic scale of microfluidic channels results in surface forces that are much larger than the volume force. Therefore, for certain aspects like surface force, capillary effect, slip effect, and rapid heat transfer effect [5], the micro-and nanofluidic mechanics exhibit new characteristics distinguished from the macro-scale fluid flow. At present, several mechanisms for fluid-driving in micro-and nanoflows have been put forward in the existing literature, such as the micropumps by pressure gradient [6], external electric field, external magnetic field [7][8][9][10], surface tension, and acoustic wave [11]. In most cases, a solid surface acquires a negative electric charge when brought into contact with an electrolyte solution. In return, the charged solid surface redistributes nearby ions in the surrounding solution forming a double-layer structure of ions in the vicinity of the surface of the solid, which is referred to as the electric double layer (EDL) [2,12]. The EDL is a region close to the charged surface, where the counter ions exceed the ions, thus neutralizing the surface charge. When an axial electric field is applied, the electric power generated by the interaction between the net charge density in EDL and the applied electric field can drive the flow. This kind of flow is known as electroosmotic flow (EOF) [2,12].
Due to its relative advantages, such as simple design requirements and no need for moving elements, the electroosmosis force has been extensively discussed as a more effective driving mechanism than a simple pressure driving mechanism. In this respect, much of the literature on the microchannel electroosmosis flow by Newtonian fluid and non-Newtonian fluid has been well studied [13][14][15][16][17][18][19][20][21][22][23][24][25]. Yang et al. [17] analyzed the EOF in an annulus with high zeta potentials. Kadet and Koryuzlov [22] have presented the results of experimental studies of EOF of two immiscible liquids in a thin slit channel. Through experimental, simulation and theoretical analyses in polymer-coated slits, Monteferrante and Melchionna [23] pointed out that electroosmotic suppression arises from the frictional forces acting on the ionic currents.
However, when the voltage is high, Joule heating is inevitable in the process of EOF. How to dramatically improve the performance of conventional heat transfer has become a major challenge for scientists and engineers. Therefore, heat transfer has become the focus of EOF research in recent years [26][27][28][29][30][31][32][33][34][35][36][37][38]. Under the condition of a constant heat flux boundary, Maynes and Webb [31] studied the thermal transport of EOF in parallel plates and circular microchannels and obtained the analytical expression of temperature distribution. In addition, they also studied the viscous dissipation effect [32]. Liechty et al. [33] investigated the fully-developed convection heat transfer of EOF in a circular microtube with arbitrary wall zeta potentials under imposed wall temperatures and imposed wall heat fluxes. Tang et al. [34] have obtained a numerically transient temperature field that takes into account Joule heating in a microchannel. Chakraborty [35] has gained analytical solutions of the temperature and the Nusselt number of the fully developed flow driven by both electroosmotic force and imposed pressure gradient. Sadeghi and Saidi [36] have investigated the influence of viscous dissipation on thermal transport characteristics of the fully developed combined electroosmotic and pressure-driven flow in a parallel-plate microchannel with uniform wall heat flux. The research results show that viscous dissipation and Joule heating dominate the heat transfer problem in EOF.
However, most studies on EOF problems are limited to no-slip boundary conditions. Actually, the experimental results show that when the fluid flows over the solid surface, there may be relative motion between the solid surface and the liquid near the solid surface. The velocity of the solid surface is called slip velocity u slip , which is related to slip length b. The slip length is usually in the range of several nanometers to tens of micrometers [39]. Previous studies have shown that EDL and slip can significantly affect micro-and nanoscale fluid behaviors [3,[40][41][42][43][44][45][46][47][48][49]. In the case of electroosmotic flow, Rojas et al. [19] conducted a dispersion analysis in the neutral non-reacting solute in the circular microchannel with boundary slip and concluded that the existence of slip can greatly magnify the dispersion increase caused by the induced pressure gradient caused by the non-uniform wall potential [20]. Kundu and Saha [24] discussed the thermal transport behavior of electrostatic-driven flow in microchannels based on the no-slip, first-and second-order slip at the boundary of velocity distribution and the no-jump, first-and second-order heat slip of thermal response under the condition of maintaining uniform wall heat flux. In a microchannel between two parallel plates with imposed heat flux, the liquid flows with the slip boundary condition have been numerically investigated by Ngoma and Erchiqui [38]. Park and Kim [46] studied the Electrokinetic slip flow through hydrophobic microchannels and obtained an analytical solution of volumetric flow rate and solute retention time.
Although most of the studies in the literature were made to investigate the effects of EOF and convective heat transfer on liquid flow, the impacts of the surface charge about boundary slip on heat transfer in electroosmosis flows were rarely concerned about [40,50].
Under the influence of the surface charge, EDL with an opposite net charge is formed in the solution near the surface of the wall. Then, this charge distribution near the wall results in an attractive electric force between the charged wall and the nearby fluid. Eventually, this interaction between the wall and the fluid further influences the slip boundary conditions [41,50]. Utilizing the molecular dynamics simulation method, Joly et al. [50] studied and used a mathematical model to describe the effect of surface charge on the boundary slip. Jing and Bhushan [41,51] used atomic force microscopy to study the relationship between the boundary slip and the surface charge and found the same relationship as that obtained by Joly et al. [50]. So far, although most researchers have found the dependence of boundary slip on a surface charge, the study of fluid types is mainly limited to Newtonian fluid [50][51][52][53][54][55][56]. Bhushan et al. [40,41] studied the effect of surface charge-dependent boundary slip on the electroviscous effect between microchannels. Pan et al. [57,58] further investigated the effect of surface charge-dependent slip on the electroviscous effect and heat transfer in a microchannel and obtained analytical expressions for velocity, temperature distributions, and Nusslt number.
The EOFs in slit channels have a variety of practical applications, such as micro/nanopumping, mixing, and oil extraction [22]. However, a large number of previous studies rarely discussed this problem theoretically. The purpose of this paper is to investigate the effects of the surface charge-dependent slip on the combined electroosmotic and pressuredriven flow in a parallel-plate micro and nanochannel with high zeta potential and constant wall heat flux. The paper is organized as follows: Section 2 presents the mathematical modeling. Results and comparisons are discussed in Section 3. The conclusions are drawn in Section 4.

Electric Potential Distribution
In this study, the steady combined electroosmotic and pressure-driven flow of a symmetric binary electrolyte between two parallel plates with surface charge-dependent slip and high zeta potential is studied with density ρ, dynamic viscosity µ, and electric conductivity σ of the fluid. A Cartesian orthonormal coordinate system (x, y, z) is used with the origin fixed at the middle of the channel, as shown in Figure 1. The channel width W in the z-direction and length L in the x-direction is much larger than the channel height 2H in the y-direction, i.e., 2H << L and W << L. EOF in the x-direction is generated when an electric field is applied along the x-direction. On the boundary, the effect of surface charge density σ s on the boundary slip b is considered [41,53]: where b 0 is the original slip length without considering surface charge effect, α~1 is a numerical factor, e is the elementary charge, d is the equilibrium distance of Lennard-Jones potential, and l B = e 2 /(4πεk B T av ), ε is the permittivity of the electrolyte solution, k B is the Boltzmann constant and T av is the average absolute temperature over the channel crosssection [53,54]. According to the EDL theory, the electric potential ψ and the net charge density ρ e satisfy Poisson equation [2,12] The ideal solution of fully dissociated symmetric salt, the electric charge density is given by [2,12]: where n 0 is the ion density, z v is the valence number of ions in solution. Using Equations (2) and (3), we can obtain the Poisson-Boltzmann equation The ideal solution of fully dissociated symmetric salt, the electric charge density is given by [2,12]: where n 0 is the ion density, z v is the valence number of ions in solution. Using Equations (2) and (3), we can obtain the Poisson-Boltzmann equation The following dimensionless parameters are introduced: where K is called the nondimensional electrokinetic width, κ = (2n 0 z v 2 e 2 /(εk B T av )) 1/2 is Debye-Hückel parameter and 1/κ is the characteristic thickness of the EDL. The dimensionless form of the Poisson-Boltzmann equation can be expressed as: Equation (6) is subject to the following boundary conditions: here, Z is the dimensionless zeta potential, i.e., Z = z v eζ/(k B T av ). From Equation (6) and boundary conditions (7), one obtains the dimensionless potential distribution:

Velocity Distribution
Since the velocity of EOF and the hydraulic diameter of micro/nano channels are very small, the Reynolds number is very low, and hence fully-developed laminar flow can be assumed. In the fully-developed parallel flow ( = 0, = 0, = 0), the Navier-Stokes equation becomes linear. Here, u, v, and w are the velocities in the x, y, and z-directions, respectively. The momentum equation in the x-direction under the combined action of pressure and electric field E0 is as follows: Equation (9) is subject to the following boundary conditions: The following dimensionless parameters are introduced: where K is called the nondimensional electrokinetic width, κ = (2n 0 z v 2 e 2 /(εk B T av )) 1/2 is Debye-Hückel parameter and 1/κ is the characteristic thickness of the EDL. The dimensionless form of the Poisson-Boltzmann equation can be expressed as: Equation (6) is subject to the following boundary conditions: here, Z is the dimensionless zeta potential, i.e., Z = z v eζ/(k B T av ). From Equation (6) and boundary conditions (7), one obtains the dimensionless potential distribution:

Velocity Distribution
Since the velocity of EOF and the hydraulic diameter of micro/nano channels are very small, the Reynolds number is very low, and hence fully-developed laminar flow can be assumed. In the fully-developed parallel flow ( ∂u ∂x = 0, v = 0, w = 0), the Navier-Stokes equation becomes linear. Here, u, v, and w are the velocities in the x, y, and z-directions, respectively. The momentum equation in the x-direction under the combined action of pressure and electric field E 0 is as follows: Equation (9) is subject to the following boundary conditions: where the slip length b is given by Equation (1).
According to the charge conservation law, the surface charge should be equal to the net charge in the fluid. Therefore To write the governing equations and the boundary conditions in dimensionless form, let: where U eo is the steady Helmholtz-Smoluchowski velocity. First, substituting Equation (12) in the governing equations and in their boundary conditions and then solving the corresponding dimensionless boundary value problem by integrating Equation (9) twice, we obtain in which u r = u max /U eo is the ratio of the velocity of pressure-driven flow and that of EOF between two parallel plates, u max = −H 2 (2µ) −1 dp/dx. Note that u r = 0 corresponds to a pure EOF.
According to the definition of the dimensional volumetric flow rate, we have Introduce a dimensionless volumetric flow rate Substituting Equation (13) into Equation (14), the dimensionless volumetric flow rate can be expressed as The integral in Equation (16) can be calculated numerically using Gaussian quadrature.

Temperature Distribution
Energy conservation, including the effects of Joule heating and viscous dissipation, can be written as where s and µ(du/dy) 2 denote the rate of volumetric heat generation due to Joule heating and viscous dissipation, respectively. s = i e 2 σ ρ , where i e is the current density and σ ρ is the liquid electrical resistivity. The current density is given by [33,36] where σ ρ −1 = 2z v 2 e 2 Dn 0 /(k B T av ), D is the diffusivity of ions in the electrolyte. The boundary conditions for the energy equation are as follows: The dimensionless temperature θ is introduced This depends only on y for the fully developed flow. Differentiating both sides of Equation (20) with respect to x gives in which T b is the bulk temperature. Considering an energy balance on a control volume with thickness dx along the centerline of the channel, we can write Equation (17) as: From Equation (22), we have where u m , β, γ are given by So the energy equation in dimensionless form can be written as where u m = u m /U eo is the dimensionless axial mean velocity, Br = µU 2 eo /qH is the Brinkman number, and S = E 0 2 H/σ ρ q is the dimensionless volumetric heat generation due to Joule heating. The thermal boundary conditions in the dimensionless form can be written as Integrating Equation (25) twice and using the boundary conditions (26), we obtain here In the above equation, for simplicity,ψ = ψ(ŷ). To obtain the Nusselt number, the dimensionless bulk temperature θ b must be calculated first, which is given by and Nusselt number is defined as

Results and Discussion
In this section, the discrete numerical solutions of flow rate Q, dimensionless temperature θ and Nusselt number are obtained by Gaussian integration in two parallel plates with surface charge-dependent slip, and constant heat flux, together with the analytical and semi-analytic solutions for velocity, flow rate, and temperature. In addition, a series of results are displayed in the form of graphs to intuitively express the coupling effect of surface charge and slip on the combined pressure and electroosmotic-driven flow in a parallel-plate micro-and nanochannel. In order to obtain veritable and effective results, the following typical parametric values are used [2,12,53,54]: d = 0.4 × 10 −9 m, α = 1, µ = 1.01 × 10 −3 Pa·s, ρ = 10 3 kgm −3 , e = 1.6 × 10 −19 C, z = 1, T av = 298 K, ε = 7 × 10 −10 Fm −1 , k B = 1.38 × 10 −23 JK −1 , D = 1.612 × 10 −9 m 2 s −1 , without special instructions dp/dx = 0 Pa/m.
In Figure 2, the present results in the case of charge-independent slip are compared with those obtained by Mondal et al. [37] and Ngoma et al. [38]. Ngoma and Erchiqui [38] investigated the fluid flow with the slip boundary condition in a microchannel between two parallel plates; that is, the case of low zeta potential and the slip length b = b 0 . In addition, Mondal and Shit [37] took into account the non-uniform walls varying sinusoidally, which can be simplified into the model in this paper in the limit case. The present results are found to be in agreement with those in Refs. [37,38]. Moreover, there is a slight error near the wall, which is mainly caused by the Debye-Hückel approximation used in Refs. [37,38]. Fm −1 , kB = 1.38 × 10 −23 JK −1 , D = 1.612 × 10 −9 m 2 s −1 , without special instructions dp/dx = 0 Pa/m. In Figure 2, the present results in the case of charge-independent slip are compared with those obtained by Mondal et al. [37] and Ngoma et al. [38]. Ngoma and Erchiqui [38] investigated the fluid flow with the slip boundary condition in a microchannel between two parallel plates; that is, the case of low zeta potential and the slip length b = b0. In addition, Mondal and Shit [37] took into account the non-uniform walls varying sinusoidally, which can be simplified into the model in this paper in the limit case. The present results are found to be in agreement with those in Refs. [37,38]. Moreover, there is a slight error near the wall, which is mainly caused by the Debye-Hückel approximation used in Refs. [37,38]. In Figure 3, the non-slip velocity profile is plotted as a function of the distance for different values of the zeta potential and slip length. It can be seen from Figure 3a that for large zeta potentials, the velocity distribution of no-slip flow varies with zeta potential in the thin layer near the solid walls. As shown in Figure 3, it is found that the velocities of both the no-slip and slip flows increase with the zeta potential. For the surface charge-dependent slip flow, the maximum velocity at ζ = −125 mV is approximately twice as large as that at ζ = −25 mV. The study also demonstrated that the velocity of the surface charge-independent slip flow is larger than that of the surface charge-dependent flow. The surface charge-dependent expresses the case of b = b0. In addition, the difference between the velocities of the surface charge-dependent and charge-independent slip flows increases with zeta potential, owing to the slip length decreasing with zeta potential, as shown in Equation (1). In Figure 3, the non-slip velocity profile is plotted as a function of the distance for different values of the zeta potential and slip length. It can be seen from Figure 3a that for large zeta potentials, the velocity distribution of no-slip flow varies with zeta potential in the thin layer near the solid walls. As shown in Figure 3, it is found that the velocities of both the no-slip and slip flows increase with the zeta potential. For the surface chargedependent slip flow, the maximum velocity at ζ = −125 mV is approximately twice as large as that at ζ = −25 mV. The study also demonstrated that the velocity of the surface charge-independent slip flow is larger than that of the surface charge-dependent flow. The surface charge-dependent expresses the case of b = b 0 . In addition, the difference between the velocities of the surface charge-dependent and charge-independent slip flows increases with zeta potential, owing to the slip length decreasing with zeta potential, as shown in Equation (1). In Figure 4, the effect of the surface charge-dependent slip on the flow rate is investigated. Here, the QID represents the flow rate without considering the effect of the surface charges on the slip length. The surface charge effect is increased when the magnitude of zeta potential ζ increases and the slip length b0 increases. Obviously, the slip length decreases with the zeta potential. A reasonable explanation is the increase of the attractive electrostatic force between the charged solid surface and the liquid near the solid-liquid interface. In addition, Figure 4 shows that the surface charge effect on the slip length in the microchannel is smaller than in the nanochannel. In Figure 4, the effect of the surface charge-dependent slip on the flow rate is investigated. Here, the Q ID represents the flow rate without considering the effect of the surface charges on the slip length. The surface charge effect is increased when the magnitude of zeta potential ζ increases and the slip length b 0 increases. Obviously, the slip length decreases with the zeta potential. A reasonable explanation is the increase of the attractive electrostatic force between the charged solid surface and the liquid near the solid-liquid interface. In addition, Figure 4 shows that the surface charge effect on the slip length in the microchannel is smaller than in the nanochannel.
In Figure 5, the temperature profile is plotted against the distance for different values of the zeta potential and the slip length. As seen in Figure 5, the dimensionless temperature in the nanochannel increases when the magnitude of the zeta potential increases. The reason is that as the magnitude of the zeta potential increases, the velocity and hence the effective Joule heating increases, and thus the heat transfer can be enhanced. Furthermore, it can be seen from Figure 5 that the dimensionless temperature in the surface charge-dependent slip flow is larger than that in the surface charge-independent slip flow. The viscous dissipation behaves like an energy source, increasing the temperature of the fluid, especially near the wall. The reason is that the highest shear rate occurs at this region while it is zero at the centerline. In addition, since the walls are cooling, heat is transferred from the interior to the exterior of the fluid, and fluid convection carries away a fraction of this heat, so the amplitude of the fluid velocity is large, further leading to a decrease in the amplitude of the temperature. For the surface charge-dependent slip flow, the maximum temperature at ζ = −125 mV is approximately four times larger than that at ζ = −25 mV. In Figure 5, the temperature profile is plotted against the distance for different values of the zeta potential and the slip length. As seen in Figure 5, the dimensionless temperature in the nanochannel increases when the magnitude of the zeta potential increases. The reason is that as the magnitude of the zeta potential increases, the velocity and hence the effective Joule heating increases, and thus the heat transfer can be enhanced. Furthermore, it can be seen from Figure 5 that the dimensionless temperature in the surface charge-dependent slip flow is larger than that in the surface charge-independent slip flow. The viscous dissipation behaves like an energy source, increasing the temperature of the fluid, especially near the wall. The reason is that the highest shear rate occurs at this region while it is zero at the centerline. In addition, since the walls are cooling, heat is transferred from the interior to the exterior of the fluid, and fluid convection carries away a fraction of this heat, so the amplitude of the fluid velocity is large, further leading to a decrease in the amplitude of the temperature. For the surface charge-dependent slip flow, the maximum temperature at ζ = −125 mV is approximately four times larger than that at ζ = −25 mV.
Nondimensionalizing this parameter by the total Joule heating per unit length in the case of low zeta potential, 2H (E 0 2 /σ ρ ), we obtain The total Joule heating per unit length is defined as Nondimensionalizing this parameter by the total Joule heating per unit length in the case of low zeta potential, 2H (E0 2 /σρ), we obtain The dimensionless total Joule heating per unit length of the channel Ee is plotted as a function of zeta potential ζ in Figure 6 for different values of H and n0. It is found that for low zeta potential (ζ < 25 mV), Ee = 1 at all values of H. The reason can be illustrated that for low zeta potential, cosh ψ~1. For a given zeta potential, the effective Joule heating Ee decreases as the half-height H of the micro and nanochannel increases. On the contrary, for a given H, the zeta potential strengths the effective Joule heating Ee, especially in nanochannels. By comparing Figure 6a,b, it can be found that ion density n0 can inhibit the increase of Joule heat. The dimensionless total Joule heating per unit length of the channel E e is plotted as a function of zeta potential ζ in Figure 6 for different values of H and n 0 . It is found that for low zeta potential (ζ < 25 mV), E e = 1 at all values of H. The reason can be illustrated that for low zeta potential, cosh ψ~1. For a given zeta potential, the effective Joule heating E e decreases as the half-height H of the micro and nanochannel increases. On the contrary, for a given H, the zeta potential strengths the effective Joule heating E e , especially in nanochannels. By comparing Figure 6a,b, it can be found that ion density n 0 can inhibit the increase of Joule heat.
For a given H, as the zeta potential ζ increases, the effective Joule heating E e increases, especially in nanochannels. In addition, as the ion density n 0 increases, the effective Joule heat E e decreases. The reason is that an increase in the channel height H, or in the ion density n 0 , implies an increase in the dimensionless electrokinetic width K; that is, the EDL becomes thinner. For a given H, as the zeta potential ζ increases, the effective Joule heating Ee increases, especially in nanochannels. In addition, as the ion density n0 increases, the effective Joule heat Ee decreases. The reason is that an increase in the channel height H, or in the ion density n0, implies an increase in the dimensionless electrokinetic width K; that is, the EDL becomes thinner. Figure 7 depicts the variations of the Nusselt number with zeta potential ζ for different values of the surface heat flux q and the electric field E0. Both positive and negative values of the wall heat flux are considered, where positive and negative values of q correspond to the case where the fluid is cooled and heated, respectively. In general, increasing the Br decreases the Nusselt number, so the Br decreases with the heat flow flux q. It can be seen that the Brinkman number has a greater influence on the Nusselt number at greater values of heat flux q (see Figure 7a). From Figure 7, we can see that the Nusselt number increases with the positive q because the convective heat transfer coefficient h increases with q. For a given constant wall heat flux q, the Nusselt number is large when E0 and ζ are small. In contrast to the wall-heating case, in the presence of viscous heating, which occurs at large velocity gradients, q has a greater effect on the Nusselt number. As a result of viscous heating, the bulk temperature is not affected significantly, but the wall temperature is greatly increased, and the Nusselt number will be much higher. Due to the large value of viscous heating, the wall temperature is much higher than the bulk temperature. This behavior of the Nusselt number is accompanied by the occurrence of a singularity in Nusselt number values of the wall heating cases; as a result of changing the sign of dimensionless bulk temperature from negative to positive values (see Figure 7b). Increasing the value of E0 leads to smaller values of the Nu, indicating that wall temperatures increase rather than bulk temperatures when Joule heating is applied. The reason is that although the distribution of energy generated by Joule heating is uniform throughout the channel cross-section, the energy transferred by convection decreases near the wall, and it equals zero at the wall (see Figure 7c). Moreover, Nu is smaller in the surface charge-dependent slip flow than in the surface charge-independent slip flow. The reason is that the convective heat transfer coefficient h is reduced by the decrease of the velocity on the boundary due to the surface charge effect on the slip.   Figure 7a). From Figure 7, we can see that the Nusselt number increases with the positive q because the convective heat transfer coefficient h increases with q. For a given constant wall heat flux q, the Nusselt number is large when E 0 and ζ are small. In contrast to the wall-heating case, in the presence of viscous heating, which occurs at large velocity gradients, q has a greater effect on the Nusselt number. As a result of viscous heating, the bulk temperature is not affected significantly, but the wall temperature is greatly increased, and the Nusselt number will be much higher. Due to the large value of viscous heating, the wall temperature is much higher than the bulk temperature. This behavior of the Nusselt number is accompanied by the occurrence of a singularity in Nusselt number values of the wall heating cases; as a result of changing the sign of dimensionless bulk temperature from negative to positive values (see Figure 7b). Increasing the value of E 0 leads to smaller values of the Nu, indicating that wall temperatures increase rather than bulk temperatures when Joule heating is applied. The reason is that although the distribution of energy generated by Joule heating is uniform throughout the channel crosssection, the energy transferred by convection decreases near the wall, and it equals zero at the wall (see Figure 7c). Moreover, Nu is smaller in the surface charge-dependent slip flow than in the surface charge-independent slip flow. The reason is that the convective heat transfer coefficient h is reduced by the decrease of the velocity on the boundary due to the surface charge effect on the slip.

Conclusions
In the present study, the fully developed combined pressure and electroosmotic driven flows in parallel plate micro and nanochannels with surface charge-dependent slip and high zeta potentials have been studied. The classical boundary conditions of surface charge-dependent slip and uniform wall heat flux are taken into account in the analysis. The effects of Joule heating and viscous dissipation are also considered. The electric potential in EDL, velocity, temperature, and Nusselt number have also been obtained analytically. The problem studied here is found to be governed by five parameters: the ratio of the velocity of pressure-driven flow and the velocity of EOF, the electrokinetic width (Debye-Hückel parameter), the zeta potential, Joule heating term, and Brinkman number. The main results of this study can be summarized as follows: • For a given large zeta potential, the EDL velocity variation is limited to the thin layer that clings to both solid walls.

•
Compared with the velocity of the surface charge-independent slip flow, that of the surface charge-dependent slip flow showed a significant reduction.

•
The increase in zeta potential causes an increase in the dimensionless temperature in the nanochannel.

•
The dimensionless temperature of the surface charge-dependent slip flow is larger than that of the surface charge-independent slip flow.

•
For a given zeta potential, the effective Joule heating Ee decreases as the channel height increases; for a given channel height, as the zeta potential ζ increases, the effective Joule heating Ee is increased, especially in the nanochannel.

•
Ion density n0 can inhibit the increase of Joule heat.

•
The Nu decreases with Joule heating S and increases with positive heat transfer coefficient q.

•
The Nu is large when the applied electric field E0 and zeta potential ζ are small.

•
The Nu in the surface charge-dependent slip flow is smaller than that in the surface charge-independent slip flow.

Conclusions
In the present study, the fully developed combined pressure and electroosmotic driven flows in parallel plate micro and nanochannels with surface charge-dependent slip and high zeta potentials have been studied. The classical boundary conditions of surface chargedependent slip and uniform wall heat flux are taken into account in the analysis. The effects of Joule heating and viscous dissipation are also considered. The electric potential in EDL, velocity, temperature, and Nusselt number have also been obtained analytically. The problem studied here is found to be governed by five parameters: the ratio of the velocity of pressure-driven flow and the velocity of EOF, the electrokinetic width (Debye-Hückel parameter), the zeta potential, Joule heating term, and Brinkman number. The main results of this study can be summarized as follows: • For a given large zeta potential, the EDL velocity variation is limited to the thin layer that clings to both solid walls.

•
Compared with the velocity of the surface charge-independent slip flow, that of the surface charge-dependent slip flow showed a significant reduction.

•
The increase in zeta potential causes an increase in the dimensionless temperature in the nanochannel.

•
The dimensionless temperature of the surface charge-dependent slip flow is larger than that of the surface charge-independent slip flow. • For a given zeta potential, the effective Joule heating E e decreases as the channel height increases; for a given channel height, as the zeta potential ζ increases, the effective Joule heating E e is increased, especially in the nanochannel. • Ion density n 0 can inhibit the increase of Joule heat.

•
The Nu decreases with Joule heating S and increases with positive heat transfer coefficient q.

•
The Nu is large when the applied electric field E 0 and zeta potential ζ are small.

•
The Nu in the surface charge-dependent slip flow is smaller than that in the surface charge-independent slip flow.