1. Introduction
The development of the micro-electric mechanical system (MEMS) has received great attention because of its innovative application in chemical, medical, and biological-related industries. For instance, Lab-on-a-chip, which conducts experiments by processing bioliquids or chemical solutions on a microchip, has become increasingly popular [
1]. Electroosmotic flow (EOF) being a major electrokinetic effect refers to the motion of an ionized liquid inside microchannels with respect to the stationary charged channel walls under an external electric field applied tangentially along the microchannel [
2]. An electroosmosis-based micropump has been successfully developed [
3] and became one of the most important components in lab-on-a-chip since it has advantages, such as the production of pulse-free and plug-like flow, and the dependence on non-mechanical parts [
4] over the conventional pressure-driven flow (PDF). Therefore, a fundamental understanding of the mechanism of EOF is vital for the precise control and optimal design of microfluidic devices.
The hydrodynamic behavior of EOF of different fluids is investigated in different microchannels, such as slit [
5], rectangular [
6,
7], elliptic [
8], and circular microchannels [
9]. The core of the attention on EOF has shifted to the heat transfer characteristics of EOF [
9,
10,
11], two-layer EOF [
12], rotating EOF [
13,
14], and pressure effects on EOF [
15]. In biological and chemical industries, biofluids, such as blood and DNA, solutions manipulated in microfluidic devices show nonlinear rheological behavior, such as the viscosity dependent on shear rate, which cannot be modeled by the linear constitutional relation. Therefore, the theory of non-Newtonian fluids offers researchers in the fields of mathematics and physics challenges in developing solutions to the nonlinear governing equations. The EOF of Jeffrey fluid [
16] and Maxwell fluid [
17] under an alternating current (AC) electric field are investigated in terms of the influence of nonlinear rheological behavior on flow performance. The corresponding analytical solution for the velocity field has been developed by adopting the method of variable separation. Das and Charkraborty [
5] first proposed a non-Newtonian power-law model for the EOF of biofluids and acquired analytical solutions of velocity, temperature, and concentration distribution. The power-law model has received much more attention because of the concise expression and breadth coverage for working liquids in a MEMS. Zhao et al. highlighted the dynamic properties of EOF of power-law fluid under an AC electric field along a rectangular microchannel [
18]. My previous works conducted parametric studies on the EOF of power-law fluids in terms of the flow behavior index and electrokinetic width in a cylindrical circular microcapillary [
19] and in a rectangular microchannel [
20]. It showed that the nonlinear rheological behavior of power-law fluid has a significant influence on the flow pattern and volumetric flow rate of EOF. Particularly, the shear thinning fluid needs more developed temporal length and is more sensitive to the change in parameters than that of shear thickening fluid.
The thermal characteristics have to be carefully considered due to the involvement with fluid flow and the applications in the cooling systems of a MEMS. The Joule heating effect that is an inherent phenomenon in EOF resulting from the ohmic resistance of the electrolyte dominates the heat generation of fluid flow. It results in the change of temperature-dependent electrical/transport properties, such as viscosity and electric conductivity of the working liquids, which, in turn, influences the hydrodynamic behavior of EOF. Concerning the Joule heating-induced heat transfer of EOF, the models in a microcapillary [
9] and in a rectangular microchannel [
21] without consideration of viscous dissipation have been developed. A model in a slit microchannel with consideration of viscous dissipation and thermal radiation is developed by Shit et al. [
10], which has been extended to power-law fluid in a circular microchannel [
22]. The entropy generation rate for the EOF of a power-law fluid in a microchannel has been obtained, and it turns out that entropy generation is dominated by the Joule heating effect [
23]. In addition, the microchannel flow can be assisted by the mixed effects of electroosmosis and pressure gradient. The hydrodynamic and thermal behavior of combined EOF and PDF of Newtonian fluid in a microchannel was studied by solving the velocity distribution, temperature distribution, and Nusselt number [
24]. The model above has been extended to the case of power-law fluid in a slit microchannel [
25] and in a rectangular microchannel [
26], where the temperature distribution and Nusselt number have been numerically solved by taking into account the Joule heating effect and viscous dissipation. Lately, the effects of magnetic field on the heat generation of EOF has been given importance [
27].
Compared with the conventional fluid, the nanofluid containing nanoparticles with a diameter of 1 to 100 nm shows larger overall thermal conductivity. This is because the thermal conductivity of solid particles, such as Cu, is 700 times higher than the thermal conductivity of water, and the thermal conductivity of metal oxide, such as AI
2O
3, is much greater than that of a pure fluid, and hence, the thermal performance of nanofluid is dramatically enhanced in the presence of these nanoparticles. Furthermore, the nanofluid is stable, which has no additional problems like pressure drop, sedimentation, or blockage in microchannels. Choi firstly proposed the model of nanofluid [
28]. In the fields of energy, power, electronics, etc., the trend for using suspensions of nanoparticles in a fluid medium is growing [
29], which has led to several researches on the nanofluid flow [
30,
31]. The comparison among several existing models for thermal conductivity was carried out to assist in the accurate prediction of thermal behavior of nanofluid in [
32]. A more comprehensive survey on the turbulent heat transfer of nanofluid flow was conducted [
33]. Since the characteristics of nanofluid mentioned above meet the increasingly higher requirements in heat transfer of heat exchange equipment of MEMS, the nanofluid flow in microscale has been paid more and more attention recently. The electrokinetic effects of nanofluids induced by the streaming potential on PDF in microtubes were theoretically studied [
34]. The transient thermal characteristics of the combined EOF and PDF of nanofluids in a microchannel under the effect of a magnetic field were investigated [
35]. Furthermore, as heat transfer performance of fluid flow is increasingly encountered in micro electrics, the base fluid is not only limited to Newtonian fluid, but also the emphasis has been laid on the non-Newtonian fluid flow. Therefore, in microscale devices with characteristics of compact structure, to achieve the efficient thermal management of the non-Newtonian fluid, the power-law nanofluid was developed, where nanoparticles are added, and the base fluid rheology is described by the power-law model. The EOF of power-law nanofluid in a parallel plate microchannel was studied where one-dimensional momentum equation and energy equation were analytically solved [
36]. A mixed convection flow and the corresponding heat transfer of shear thinning power-law nanofluid was respectively investigated by Ellahi et al. [
37] and by Si et al. [
38]. It turns out that the solid volume fractions exert special effects on both the velocity field and temperature field of fluid flow.
An up-to-date literature review indicates that for non-Newtonian fluids, the studies on the thermal behavior of EOF lack detailed insights into the substance of microstructure, namely the role that nanoparticles play, and the influence of channel geometry, owing to the complexity of governing equations for non-Newtonian nanofluid EOF in complex microchannels. Moreover, the thermal behavior of pure EOF of a power-law fluid in complex microchannels has not been studied yet. Due to the practical and fundamental significance of EOF of power-law nanofluid in a complex microchannel, such as a rectangular channel, and motivated by the works reviewed above, it is necessary to provide a fundamental understanding of the Joule heating induced heat transfer characteristics. This paper develops a model for thermally fully developed EOF of power-law nanofluid in a rectangular microchannel. The viscosity of power-law nanofluid depends both on the flow behavior index and the volume fraction of nanoparticles. By numerically solving the modified momentum equation and the energy equation, the velocity distribution, temperature distribution, and Nusselt number under the influence of the Joule heating effect are evaluated for different nanoparticle volume fraction, electrokinetic width, flow behavior index. The results are useful for the optimal design and active control of microfluidic systems.
2. Mathematical Modeling
As sketched in
Figure 1, a rectangular microchannel of width 2
b and height 2
a is considered. The channel is filled with power-law nanofluid with a dielectric constant
ε and is subject to constant heat flux
qs. The channel wall is uniformly charged with zeta potential
ξ, and the electric double layers (EDLs) will not overlap. As an electric field is applied tangentially along the channel (
z-direction), the liquid inside is set in motion due to the existence of free ions in EDLs. Owing to the symmetry, the following analysis is limited to a quarter of the channel,
, and the coordinate system used is shown in
Figure 1.
2.1. Electric Potential Distribution
Based on the electrostatics theory, the relationship between the electric potential distribution
ψ(
x,
y) and the local net charge density
ρe is expressed by the Poisson equation.
The local net charge density
ρe is subject to Boltzmann distribution.
where
χ denotes the valence of ions,
e denotes the elementary charge,
n0 is the ion density of the bulk liquid,
kb is the Boltzmann constant, and
T0 implies the absolute temperature. Substituting Equation (2) to Equation (1) yields the two-dimensional Poisson–Boltzmann (P-B) equation for the electric potential distribution of EDL in the rectangular microchannel,
The corresponding boundary conditions are given as
To facilitate the discussion, the following scales are introduced
,
,
,
α = b/a,
, with
K defined as the dimensionless electrokinetic width and 1/
as the thickness of EDL, where
= [2
χe2n0/(
εkbT0)]
1/2. Accordingly, the dimensionless P-B equation and the boundary conditions are obtained
For small zeta potential (
ξ ≤ 0.025 V), the known Debye–Hückel linearization principle [
5]
is used to linearize the P-B equation. In this case, an analytical solution to Equations (5) and (6) can be readily solved by using variable separation method as
where
,
,
M,
N = 1,2,3,…
2.2. Velocity Distribution
The vector formation of Cauchy momentum equation reads
where
is the velocity field of EOF,
ρ represents the density of mass,
t denotes the variable of time,
is the pressure gradient,
is the electric field strength,
denotes the stress tensor which can be expressed by a second invariant of the strain rate tensor
eij = [(
∂ui/∂xj) + (
∂uj/∂xi)]/2 as
.
μeff denotes the effective viscosity coefficient. For a fully developed, incompressible and laminar EOF of power-law nanofluid, it is assumed that there is only an axial velocity
w =
w(
x, y) [
7,
20], accordingly, the effective viscosity of the power-law nanofluid is
μeff = μf/(
1 − ϕ)
5/2. The viscosity of the base fluid is
μf = η(
|∂w/∂x|n − 1, |∂w/∂y|n − 1) which is a function of strain rate, flow consistency index
η of dimension (Nm
−2s
n) [
7,
20], and volume fraction of nanoparticle
ϕ [
39]. Finally, one can express the constitutive equation of power-law nanofluid as
where τ is the shear stress, and
n is the flow behavior index.
n = 1 represents a Newtonian nanofluid,
n < 1 and
n > 1 represents the shear thinning nanofluids (pseudoplastic nanofluids) and the shear thickening nanofluids (dilatant nanofluids), respectively. Consequently, the modified Cauchy momentum equation for power-law nanofluid EOF is
The corresponding boundary conditions are
Introducing the following scales
with
where
μ0 is the viscosity coefficient of Newtonian fluid of dimension (Nm
−2s), one can obtain the dimensionless modified Cauchy momentum equation and the boundary conditions:
When considering Newtonian nanofluid, i.e.,
n = 1, and using the known Debye–Hückel principle, an analytical solution can be obtained by the method of variable separation and the method of constant variation [
40], which is expressed as
where
with I = 1,2,3,…,
,
,
,
,
, , , ,
.
In the case of non-Newtonian power-law nanofluid flow (n ≠ 1), a numerical algorithm is proposed to obtain the velocity field, which is presented in the next section.
2.3. Temperature Distribution
Assuming the thermophysical properties are constant, the heat generation induced by the viscous dissipation is much smaller than that induced by Joule heating [
23,
41] and thus, is neglected in the present analysis. The energy equation for the flow field in a rectangular microchannel is
with
and
[
32,
42].
T denotes temperature field,
ω implies the ratio of the nanolayer thickness to the original particle radius.
σ, k, and (
ρcp) denote electrical conductivity, thermal conductivity, and heat capacity of power-law nanofluid at the reference pressure, respectively. The subscripts
s,
f, and
eff denote the solid nanoparticles, base fluid, and nanofluid, respectively. The second term at the right-hand side of Equation (13) denotes the volumetric heat generation owing to the Joule heating effect. The corresponding boundary conditions are
where
qs = h(
Tw − Tm) implies the constant wall heat flux and
h is the convective heat transfer coefficient. Subscripts
w and
m denote the local wall value and mean value, respectively. As assumed in the preceding section, for thermally fully developed EOF, one has
and thus
. It is noted that the temperature variation with respect to
z-direction cannot be neglected and thus, the temperature on the wall is a function of
z. Plus, an overall energy balance at an elemental control volume on a length of duct
dz gives:
with the mean velocity
.
Introducing
and the Joule heating parameter
, namely, the ratio of Joule heating to the heat flux from the channel wall, the dimensionless energy equation and the corresponding boundary conditions are obtained as
In the case of Newtonian fluid and using the known Debye–Hückel linearization principle, an analytical solution of temperature distribution is obtained by the method of variable separation and the method of constant variation,
where
,
, , , and
.
The detailed procedure is presented in
Appendix A.
For non-Newtonian power-law nanofluid, the energy equation can be solved using the numerical algorithm developed later.
2.4. Entropy Generation Analysis
With the temperature distribution obtained, Nusselt number indicating heat transfer rate from the channel wall to the liquid inside the channel is deduced as:
where the mean temperature is
and the convective heat transfer coefficient is
h = qs/(
Tw − Tm).
The entropy generation rate per unit volume is presented, which evaluates the amount of useful energy destroyed during the process, namely, the thermodynamic irreversibility of EOF [
35,
43,
44].
Since the heat generated by viscous dissipation is negligible [
23,
41], the irreversibility of the system equals the summation of local volumetric entropy generation rate due to heat transfer and Joule heating. The corresponding nondimensional form is
where
[
44]. The global total entropy generation rate can be obtained by integrating Equation (22) over the entire domain
. Then an important dimensionless number, Bejan number, representing the relative dominance of entropy generation due to heat transfer and Joule heating [
45], is presented as
It is noted that Be > 0.5 denotes that heat conduction irreversibility dominates and Be < 0.5 denotes the Joule heating irreversibility dominates.
3. Numerical Algorithm
Due to the nonlinearity of the governing equations above, the following numerical algorithm is developed to solve the velocity and temperature. Let
,
,
, where
is the temporal step,
and
are the spatial steps with
l = 1, 2,…,
L, and
i = 1, 2, …,
I,
j = 1, 2,…,
J.
represents the matrix obtained by discretizing a variable
f over the points above. A numerical algorithm is proposed using the time splitting method in terms of time and the compact difference method in space. The complete P-B equation is numerically carried out based on the compact difference schemes, and the corresponding discretized form can be found in [
20] and not presented here for conciseness.
To obtain the discretized velocity field
with
and the discretized temperature field
by using the iterative method, one can introduce
on the left-hand side of Equations (11) and (16) and set a nonzero initial value for them,
More specifically, each time step is split into two steps. In the first step, from
to
, one can focus on the third term on the right-hand side of Equations (19) and (20), and rewrite them as
Accordingly, the velocity field at
is given as
. To obtain temperature field
at
, Equation (22) is discretized as
where
,
, and
. The subscript
m implies the mean value of
. Setting
and the discretized temperature at
as
, the Runge–Kutta method is applied to Equation (23):
where
α1 =
β1 = 1,
α2 = 3/4,
β2 = 1/4,
α3 = 1/3,
β2 = 2/3.
In the second step, from
to
, one can focus on the second-order derivatives in Equations (19) and (20), which are generally expressed as
To note, as solving
,
and
, the corresponding discretized forms have been given in
Appendix B. When solving
,
C1 = 1 and
C2 = 1. According to the Crank–Nicolson technique, the discretized form of Equation (25) is
Using the compact difference schemes and following the procedure presented in
Appendix B,
fi,jl+1 can be obtained as
where
,
,
,
,
, and
,
,
,
,
,
,
,
.
As a result, replacing f with and , the discretized solutions and can be computed. Finally, a specified criterion is given to identify that if the velocity is fully developed, i.e., . Then, the fully developed velocity and temperature are obtained based on the iterative method above.
4. Results
A parametric study for velocity, temperature, and Nusselt number is conducted in the case of different flow behavior index
n, electrokinetic width
K, volume fraction of nanoparticles
ϕ, and Joule heating parameter
S. The nanoparticle is regarded as aluminum oxide [
38]. The other physical parameters of the base fluid come from [
9] and [
20]. The typical values are presented in
Table 1 below.
For validating the numerical algorithm proposed above, a test of grid dependence is conducted. As a result, the volumetric domain
Ω is discretized to a grid system of 121 × 121. The numeric algorithm shows that the numerical solutions can achieve the accuracy of
, which is demonstrated in numerical computation by obtaining the numerical error of the order 10
−8 between the numerical velocity and analytical velocity of Newtonian nanofluid EOF. Furthermore, the iterative criterion is given as
ΔL < 10
−10 to make sure the EOF is fully developed. For the EOF of Newtonian nanofluid, in
Figure 2a, the numerical velocity profile at
is compared with the analytical velocity obtained from Equation (13) and the numerical temperature profile at
is compared with the analytical temperature obtained from Equation (19) when
K = 10,
ϕ = 0.06, and
S = 3. To clearly show the comparison, 25 grid points of the numerical solution are plotted in
Figure 2. It shows that the numerical solutions show good agreement with the analytical solutions and thus, the numerical algorithm is valid, which can be used to numerically solve velocity, temperature, Nusselt number, and Bejan number of power-law nanofluid EOF.
Figure 3 illustrates velocity distributions of shear thickening nanofluid, Newtonian nanofluid, and shear thinning nanofluid at different volume fractions when
K = 10.
Figure 3a,c,e show the velocity distributions of shear thinning fluid, Newtonian fluid, and shear thickening fluid without consideration of nanoparticles in fluid, namely at
ϕ = 0.
Figure 3b,d,f show the velocity distributions of shear thinning nanofluid, Newtonian nanofluid, and shear thickening nanofluid at
ϕ = 0.06. It is noted from these two groups of figures that the velocity shows homogeneous abatement across the channel as the volume fraction of nanoparticles increases no matter what fluid rheology is considered. This is because the nanoparticles in fluid enhance the effective viscosity of power-law nanofluid in response to the sear rate and cause greater dispersion of the velocity distribution. As observed in [
7], although considering the EOF of power-law nanofluid, the shear thinning feature of base fluid leads to greater velocity compared to the shear thickening feature of the base fluid.
Without consideration of nanoparticles in fluid (
ϕ = 0) and when
K = 10, the temperature distributions of shear thinning fluid, Newtonian fluid, and shear thickening fluid at
S = 0 are displayed in
Figure 4a,c,e. For the same prescribed conditions, the temperature distributions at
S = 3 are displayed in
Figure 4b,d,f. The comparison between these two groups of figures shows that the variation of temperature along the channel becomes greater with the increase of the Joule heating parameter
S irrespective of the fluid type. Furthermore, the dimensionless temperature increases with the flow behavior index
n, indicating that the difference of heat transfer characters caused by different shear rates of fluid type is remarkable, especially for shear thinning base fluid.
To observe the influence of volume fraction of nanoparticles on the heat transfer characteristics, the temperature distributions are evaluated in
Figure 5 when the nanoparticle volume fraction
ϕ is increased to 0.06, and other parameters keep unchanged. The corresponding temperature profiles at
is illustrated in
Figure 6 for different volume fractions of the nanoparticles. From
Figure 6 and the comparison between
Figure 4 and
Figure 5, the dimensionless temperature distribution shows an increment when the volume fraction of nanoparticles
ϕ increases, or the flow behavior index
n increases. It means that the combined effects of larger thermal conductivity of nanoparticles and the change in shear stress leads to the increase in thermal diffusion effect and in further the temperature.
To observe the influence of the Joule heating parameter on the temperature, in
Figure 7, the temperature profiles at
for different Joule heating parameters
S are illustrated in the case of shear thinning nanofluid, Newtonian nanofluid, and shear thickening nanofluid when
ϕ = 0.06 and
K = 10. As seen in
Figure 5, no matter what nanofluid is considered, as the Joule heating parameter increases the variation of temperature increases, namely, the temperature difference between the centerline and the channel wall enlarges. It implies that although the energy from Joule heating effect is uniformly enhanced, the heat transferred through convection near the channel wall is much less than near the centerline. And this effect is more notable when considering shear thinning base fluid.
The variations of the Nusselt number with the electrokinetic width are shown in
Figure 8 for different fluid types and Joule heating parameters when
ϕ = 0.06. The cases
S < 0 and
S > 0 correspond to the surface cooling effect and surface heating effect, respectively. Importantly, in the case of
S < 0, the Nusselt number decreases as electrokinetic width
K increases, which increases when S ≥ 0 and other parameters remain unchanged. Surface heating (cooling) has a uniform influence on the temperature distribution, which leads to the increasing (decreasing) tendency of Nusselt number with electrokinetic width. Furthermore, the decrease in flow behavior index results in the pronounced abatement of Nusselt number, meaning that the heat transfer rate of shear thinning nanofluid is much smaller than the heat transfer rate of Newtonian fluid and shear thickening fluid. This is consistent with the prediction in [
22].
Figure 9 depicts the variations of the Nusselt number with the Joule heating parameter for different fluid types and different nanoparticle volume fraction when
K = 10. It shows that the Nusselt number is inversely proportional to the Joule heating parameter, and the decreasing rate becomes smaller as nanofluid changes from shear thickening to shear thinning. Since it is considered that Joule heating has a spatially uniform effect, the greater Joule heating parameter corresponds to relatively uniform heat distribution. Consequently, the convective heat transfer rate is reduced as presented in
Figure 7, and one can observe the decreasing trend of the Nusselt number, namely, heat transfer rate. Being consistent with the prediction in
Figure 7, the smaller convective heat transfer rate of shear thinning nanofluid leads to a smaller Nusselt number. Furthermore, the comparison between
Figure 9a,b show that no matter what base fluid and what Joule heating parameter is taken, the increase in the volume fraction of nanoparticles
ϕ significantly improves the heat transfer rate.
The variation of the Nusselt number
Nu and the variation of Bejan number
Be with flow behavior index
n for different volume fractions of nanoparticles are depicted in
Figure 10 when
K = 20 and
S = 3. The combined effects of nonlinear rheological behavior and nanoparticle volume fraction on Nusselt number are investigated. To note, the dependence of Nusselt number on the flow behavior index shows nonlinear. For shear thinning base fluid, the Nusselt number is more amenable to the change in flow behavior index
n. The increase in viscosity for shear thinning fluid resulting from the decrease in flow behavior index reduces the Nusselt number, and thus, the heat transfer rate of shear thinning nanofluid is obviously less than that of shear thickening nanofluid. Furthermore, the variation of the Nusselt number with the volume fraction of nanoparticles is almost linear. It shows that the heat transfer performance of nanofluid is clearly enhanced compared to the conventional fluid regardless of the base fluid rheology. These results are also witnessed in
Figure 10b. The increase in nanoparticle volume fraction and flow behavior index leads to the enhancement of heat transfer, which further reduces the “loss” of useful energy, i.e., the total irreversibility of the system. Moreover, the
Be > 0.5 implies that the primary irreversibility is contributed to the heat conduction when
S = 3.