Analysis of electroosmotic flow and Joule heating effect in a hydrophobic channel

Accepted Manuscript Analysis of electroosmotic flow and Joule heating effect in a hydrophobic channel A.K. Nayak, A. Haque, B. Weigand PII: DOI: Reference:

S0009-2509(17)30623-1 https://doi.org/10.1016/j.ces.2017.10.014 CES 13846

To appear in:

Chemical Engineering Science

Received Date: Accepted Date:

26 June 2017 11 October 2017

Please cite this article as: A.K. Nayak, A. Haque, B. Weigand, Analysis of electroosmotic flow and Joule heating effect in a hydrophobic channel, Chemical Engineering Science (2017), doi: https://doi.org/10.1016/j.ces. 2017.10.014

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Analysis of electroosmotic flow and Joule heating effect in a hydrophobic channel A. K. Nayaka,∗, A. Haquea , B. Weigandb,∗∗ a

Department of Mathematics, Indian Institute of Technology Rorkee, Roorkee, Uttarakhand, 247667, India b Institut f¨ ur Thermodynamik der Luft- und Raumfahrt, Pfaffenwaldring 31 70569 Stuttgart

Abstract In this paper, a mathematical model has been developed to analyze the fluid flow and heat transfer effect in a hydrophobic microchannel filled with a power law fluid. The effects of Joule heating, thermal radiation and velocity slip boundary conditions are analyzed by considering different slip parameters, EDL thickness, pressure gradient and flow behavior indices. The analytical expression for fluid flow and heat transfer have been derived in terms of the flow governing parameters based on the Debye-Huckel linearizing principle. Depending on the experimental existing flow behavior index, the analytical solutions are obtained in closed form where as numerical results are presented for general parametric values. The impact of slip velocity parameters in terms of a flow enhancement factor (Ef b ), is studied to obtain the average flow velocity variation in a hydrophobic microchannel compared to a plane microchannel. The pressure assisted flow for pseudo-plastic (shear thinning) fluids achieve maximum velocity as compared to dilatant (shear thickening) fluids. The study finds that increase in thermal radiation minimizes heat transfer rate close to the hydrophobic wall, plays a vital role for the therapeutic treatment of hyperthermia (to understand the effect of heat transfer due to electric potential). It is also observed that Joule heating parameters enhance the heat transfer rate for classical Newtonian/non-Newtonian fluids for decrease in power indices and pressure gradient. Keywords: Electroosmotic flow, Power-law fluid, Hydrophobic wall, Joule heating.

Dr. A. K. Nayak B. Weigand Email addresses: [email protected] (A. K. Nayak), [email protected] (B. Weigand) ∗

∗∗

Nomenclature cp ,

specific heat at constant pressure, Jkg −1 K −1

e,

electron charge, C

Ex ,

external electric field, V m−1

H,

height of the channel, m

kB ,

Boltzmann contant, V m−1

kp

flow consistency index, kgm−1 sn−2

kT h ,

thermal conductivity, W m−1 K −1

L,

length of the Channel, m

n,

flow behavior index

n0 ,

ion density, m−3

N r,

thermal radiation parameter

Re,

Reynolds number

S,

Joule heating parameter

T,

temperature, K

Tab ,

absolute temperature, K

u,

axil velocity, ms−1

V

velocity vector

(x, y),

cartesian coordinate

zi ,

valence of species ‘i’

Greek letters β,

slip parameter

,

electrical permittivity of the solution , CV −1 m−1

η,

dimensionless vertical coordinate

Γ,

dimensionless pressure gradient

κ,

inverse Debye length (1/λD ), m−1

ψ,

induced potential, V

ω,

electroosmotic parameter, (λD /H)

Φ

total potential, V

ρe ,

net charge density, Cm−3

ρ,

fluid density, kgm−3

θm

dimensionless mean temperature

σ,

electrical conductivity, Sm−1

τ

shear stress, N m−2

1. Introduction A distinguishing feature of the research area in non-Newtonian electro-osmotic flow is the ample applicability in many industrial and biochemical applications such as electronic cooling, chemical reaction and drug delivery in small domains [Babaie et al. 2011]. Mostly charged surfaces induce electric charge when they are in contact with the electrolyte solvent [Hunter 1981]. The counterions are attracted by the charged surface to form a thin layer of mobile charges adjacent to it, termed as electric double layer (EDL) [Bera and Bhattacharyya 2013], that consists of two layers known as stern layer and diffuse layer. The EDL thickness is characterized by the Debye length [Chen 2012], is the distance from the surface to the region over which the net charge magnitude is zero along the bulk region. The ions can move and a flow can be generated when an external electric force will act inside the domain due to the effect of viscosity, which will pull the surrounding liquid with them [Ren and Li 2001]. Reuss [Reuss 1809] initiated the study of electroosmotic flow experimentally on clay diaphragms and observed that clay particles dispersed in water migrate under the influence of electric field. In 19th century, Wiedemann repeated the experimental studies and formally introduced the mathematical theory behind it [Wiedemann 1852]. The electro-kinetic phenomena has been developed for more than a century but the application of electrokinetic transport in microfluidics is invented about thirty years ago when micro fabricated systems were developed at Stanford University (gas chromatography) and at IBM (Nozzle diffuser). Recently it is utilized in various fields of microelectomechanical systems (MEMS) used for thermal, magnetic, fluidic and optical devices [Zheng et al. 2003]. The flow can also be possible in capillary using micropumping techniques and the flow rate of pressure driven flow is Qp ∼

W h3 ∆p, µL

where L, W and h

are length, width and height of the channel respectively (fabricated by Sandia National Lab [Paul et al. 1998]). Based on the above studies it can be predicted that the flow can be possible in small channels using electroosmotic or pressure driven technique or by the combined effect of both. Furthermore, the most important field of MEMS is BioMEMS, which attracts much attention in the beginning of the new millennium. In case of analyzing BioMEMS, most fluids are considered as non-Newtonian where Newtons law of viscosity is insufficient. BioMEMS technology offers solutions for higher throughput and specific analysis for cell culture technique [Vladisavljevic et al. 2013], DNA cloning [Li and Harrison 1997], hematology, detection of bacteria [Cady et al. 2005]-[Gascoyne et al. 2004] and mixing operations [Huang et al. 2014]. The tasks are implemented in the study of biofluid whose

behaviors are considered as non-Newtonian fluid such as blood [kang and Lee 2013], saliva [Stokes and Davies 2007] and synovia fluid [Fung 1993]. A great deal of research has been made on electroosmotic flow both numerically and experimentally to visualize the flow behavior using Newtonian fluids, both electrolytes and salt free solution. However, very few studies are reported in the literature for electroosmotic flow with non-Newtonian fluid which are of great importance in BioMEMS and Lab on Chip (LOC) operation. In this context, Das and Chakraborty [Das and Chakraborty 2006] obtained an analytical solution for steady electroosmotic flow in a capillary micro-channels filled with non-Newtonian inelastic power law fluids to investigate the effects of flow velocity, temperature and concentration on a Debye length variation with channel height. Meanwhile, Chakraborty [Chakraborty 2007] developed a semi analytical mathematical model for transport of power law fluids in a rectangular micro-channel using electroosmotic effects . Also, Zhao and Yang [Zhao and Yang 2011] obtained an analytical solution for the velocity distribution of a power law fluid in a microchannel using higher range of EDL thickness. The velocity and shear stress distribution of a power law fluid in a slit microchannel are studied both numerically and analytically by Zhao et al. [Zhao et al. 2008] using a hyperbolic cosine function for the electric potential distribution. They reported that shear stress is independent of fluid behavior index. The transient electroosmotic flow of viscoelastic fluids in a rectangular microchannel is put forward by Zhao and Yang [Zhao and Yang 2009]. They showed the flow behavior of viscoelastic fluids in BioMEMS and Lab-on-a-chip for the design of microfludic devices. Electroosmotic flow of non-Newtonian power law fluids in a cylindrical microchannel is studied by Zhao and Yang [Zhao and Yang 2013], who derived an analytical solution for a specified flow behavior index in terms of generalized hypergeometric functions. They represented some interesting phenomena and counterintuitive effects for nonlinear coupling among the electrostatics channel geometry and non-Newtonian hydrodynamics which is very useful for BioMEMS design. On the other hand a lot of research concerning combined effect of pressure driven and electroosmosis is conducted to characterize the dispersion control in micromixing with the variation of the wall zeta potential [Lee et al. 2004] - [Hadigol et al. 2011]. Mixing performance of a Joule heated high viscous fluid is carried by Gopalakrishnan et al. [Gopalakrishnan et al. 2010] to study the time dependency in the flow for the optimization of the electrode configuration and the heat supply pattern to obtain a better thermal and chemical homogenization of the solution. A theoretical study on Joule heating for electroos-

motic flow in a micro channel is performed by Chang and Homsy [Chang and Homsy 2005] to analyze the linear stability using time-modulated electric field. It is observed that the Joule heating is a key parameter to stabilize the flow field for a wide range of other flow depending parameters. The flow field and heat transfer characteristics is studied for viscous fluid between two stretching disks in presence of Joule heating and viscous dissipation by Khan et al. [Khan et al. 2015]. Homotopy analysis has been used in their study to get the analytical solution for a steady MHD flow. The combined effects of EOF and pressure driven flow reported by Zade et al. [Zade et al. 2007] found the superposition of flow velocity over the independent flow by EOF or mechanical driving force. To characterize the heat transfer effects in electrokinetics, Chen [Chen 2009] studied the Joule heating effect of combined EOF and pressure driven flow, since heat transfer plays an important role for biological flows under laser irradiation [He et al. 2006]. The work is extended by Sadeghi and Saidi [Sadeghi and Saidi 2010] for a viscous heating effect. They found that viscous heating increases when the Debye Huckel parameter increases with a fixed Brinkman number. The radiation effect on biological flow such as blood flow is an important aspect of research in the field of electrokinetics due to its rapid applications in biomedical engineering [Abe and Hiraoka 1985] and several medical treatment methods [Hamann and Bickhardt 1983] - [Ismail and Jaafar 2011], such as thermal therapeutic procedures. The technique of infrared radiation is frequently used for thermal treatment of various parts of the human body [Horng et al. 2007]. This technique is mostly used in heat therapy by using infrared radiation and it is possible to heat the blood capillaries of the affected area of the human body directly. When the temperature of the surroundings is below 200 C, the human body loses heat by conduction and radiation both, but if the temperature exceeds 200 C, heat transfer takes place from the surface of the skin by the process of evaporation through sweating [Shit et al. 2016]. Prakash and Makinde [Prakash and Makinde 2011] analyzed the effect of radiative heat transfer on blood flow in a stenosed artery and found that increasing thermal radiation absorption reduces the blood flow resistance in presence of a magnetic field due to the external source and induced a pressure gradient due to stenosis. The electromagnetohydrodynamic flow of blood in a capillary having thermal radiation effect is extensively studied by Sinha and Shit [Sinha and Shit 2015] in order to find the effect of Joule heating on heat transfer coefficient. They observed that heat transfer rate is decreased with increasing Joule heating parameter. Also Shit et al. [Shit et al. 2016] studied the effect of thermal radiation and velocity slip condition in a hydrophobic chan-

nel filled with a power law fluid using electrokinetics effects. But the radiation and Joule heating effect in combined EOF/pressure driven flow are not considered. Many microfluidic devices exploit electroosmotic flows to transport the chemical species to a target position through microchannels. Recently hydrophobic materials are widely used for the fabrication of microchannels to improve a significant slip velocity at the wall.The velocity-slip effects are especially depend on the hydrophobic polymers. When the velocity-slip occurs at the walls of the microchannel, it affects the volumetric flow rate under a given strength of the electric field or pressure gradient which extend the chemical reaction and the dispersion effects of solutes inside the microchannel [Park and Kim 2009]. Hence it is utmost important to take velocity slip into account for a reliable design and operation of microfluidic devices made of hydrophobic materials. The pressure-driven transient flow of an incompressible Newtonian fluid through a circular microtube with a Navier slip boundary condition is studied by Wu et al [Wu et al. 2008]. They have found that the influences of boundary slip on the flow behavior are different for different types of pressure driven flows. A transient electroosmotic flow in micro channel is conducted experimentally by Yan et al [Yan et al. 2007]. using Nacl solution. They have implemented a phase locking based micro PIV technique to study the time dependent flow behavior using the zeta potential which is obtained by using the Smoluchowski equation with the measured steady EOF velocity in the open end of the channel without considering the Joule heating effect. They have also used the slip velocity approach for the theoretical study of EOF to validate the experimental results at the closed-end micro channel. However, there is a lack of studies on the combined effect of electroosmotic and pressure driven flow with Navier velocity slip condition. To the best of authors knowledge no research has been conducted to perform transport of electrically charged non-Newtonian fluid in microchannel with heat transfer and radiation effect due to combined EOF pressure driven flow. Thus the objective of the present work is to obtain both theoretical and numerical solution of EOF pressure driven flow using a thermally fully developed flow condition to obtain heat transfer and radiation effect of a non-Newtonian power law fluid through a microchannel with the effect of velocity slip subject to constant heat flux along the boundary. The microchannel surface is considered as hydrophobic in nature and the boundary condition is taken as a hydrodynamic Navier-slip boundary condition. The effect of thermal radiation on heat transfer along with flow variation is analyzed numerically and validated with the theoretical observations. Since, it has relevance in several biological devices used for cancer treatment, cryosurgery, laser surgery and heat therapeutic. In the first phase,

we have developed a mathematical model for combined pressure driven EOF with Joule heating effects in microchannel and the complete analytical solution is presented. In the second phase the analytical solution is compared with the developed simulated results with the variation of power law indices, slip parameters, Debye length and dimensionless pressure gradient. The work emphasizes to investigate the effects of hydrophobicity due to slip boundary conditions in presence of surface zeta potential to probe the effect of Debye length in fluid flow and heat transfer. 2. Physical configuration Consider a very large and narrow channel of rectangular cross section whose height is micro-constrained and considered to be 2H and length L shown in Fig. 1. The walls of the channel are assumed to be negatively charged and hydrophobic in nature and the flow is in x− direction. The channel is filled with a non-Newtonian power law fluid with flow consistency index kp and flow behavior index ‘n’. The electrodes are placed at the inlet and outlet of the channel and the imposed external electric field Ex is along the x− direction (see Fig. 1).

+ +

+

+ + + +

+ + +

L

+ + + +

+

b

Hydrophobic wall

2H

Y

y = +H

++ + + ++++ + +++++++ +++ + + + + + + + + X

Ex

++ + + + -1 + + + + + + + + + + + ++ + + + + + +++ κ

y = -H

L

Figure 1: Schematic diagram of problem.

3. Mathematical Formulation The governing equations for combined pressure driven EOF flow with Joule heating effect are composed by the equations for the conservation of mass, momentum, equation for potential field in the electric double layer (EDL) and energy equation. The detailed equations along with the boundary conditions are expressed as follows:

3.1. Potential Equation Electrokinetic phenomenons are based on the electrostatic charges placed at the solidliquid interface and the electro kinetic effect is found due to the formation of the EDLs. This gives rise to a finite volumetric charge density (ρe ) inside the channel nearby the charged walls. For a symmetric electrolyte solution i.e. binary fluid (z+ = −z− = z) the

volumetric charge density is given by [Probstein 2003], ρe = eΣzi ni ,

i = −1, +1.

(1)

Here ni represents the number ionic concentration and zi is the valence of ith ion. The channel is considered to be a microchannel and the EDL thickness is assumed to be small compared to the channel height and the core neutrality can be preserved. Thus, the number ionic concentration ni follows the Boltzmann distribution, which can be written as,   zi eΦ∗ , ni = n0 exp − kB Tab

(2)

where n0 represents the bulk ionic concentration, e is the electron charge, Φ∗ is the potential distribution, kB is the Boltzmann constant and Tab is the absolute temperature. According to the theory of electrostatics, the electric potential distribution due to the presence of EDL in a microchannel is described by the Poisson equation such as, ∇ · (∇Φ∗ ) = −

ρe , 

(3)

where  is the permittivity of the medium. 3.2. Flow governing equation for velocity The basic difference of non-Newtonian fluid to a Newtonian fluid is that the viscous stress is not a linear function of the rate of strain and the variations of apparent viscosity with the rate of strain is described by a number of empirical expression. According to Hadigol et al. [Hadigol et al. 2011], the relation q between the rate of strain tensor and 1 ˙ where γ˙ is the rate of strain (γ˙ : γ), its magnitude for the power-law fluid is γ˙ = 2 rate tensor and γ˙ its magnitude. The apparent viscosity can be expressed as a function of γ˙ and is denoted by µ(γ). ˙ For the power law fluid µ is given by, µ(γ) ˙ = kp (2γ) ˙ n−1 , where kp is a constant relative to the properties of the fluid and n is the flow behavior index which classify different types of fluids i.e. shear thinning (or pseudo-plastic n < 1), shear thickening (or dilatant n > 1) and Newtonian (n = 1). For shear thinning (pseudoplastic) fluid the viscosity decreases with the increase of rate of deformation and for a shear thickening (dilatant) fluid the viscosity increases with the increase of rate of deformation.

The general relationship between the stress tensor and the strain tensor rate can be written as, !
τ = 2µ(γ) ˙ γ˙ = µ(γ) ˙ ∇V + (∇V )T ,

(4)

where ∇V is the velocity gradient tensor and (∇V )T represent the transpose of the velocity gradient tensor.

The flow governed by the combined effect of pressure and electroosmotic forces along the hydrophobic microchannel and the momentum equation along with electroosmotic force term for non-Newtonian fluid can be written as [Babaie et al. 2011], D (ρV ) = −∇P + ∇ · τ + ρe E, Dt

(5)

where ρ is the density, V is the velocity vector, P is the pressure and E is the electric field strength. The continuity equation is, ∂ρ + ∇ · (ρV ) = 0. ∂t

(6)

3.3. Energy Equation In the presence of thermal radiation, viscous dissipation and Joule heating effects, the governing equation for thermal energy transport can be written as [Das and Chakraborty 2006, Babaie et al. 2012], D (ρcp T ∗ ) = ∇ · (kT h ∇T ∗ ) + ∇V : τ − ∇ · qr + q, ˙ Dt

(7)

where cp is the specific heat at constant pressure, kT h is the thermal conductivity, ∇V : τ

is the volumetric heat generation rate due to viscous dissipation, qr is the radiative heat flux and q˙ is the heat generation per unit volume due to Joule heating. The volumetric Joule heating induced by conduction current may be accurately modeled using Ohm’s law [Burgreen and Nakache 1964] for low zeta potential (which is considered here) the current

density is uniform across the channel [Levine et al. 1975] for which the energy disipation will be uniformly distributed across the microchannel. According to Ohm’s law the charge density ie = σE, where σ is the electrical conductivity of the fluid and the volumetric heat generation due to Joule heating can be written as q˙ = σ(E · E).

1

One dimensional analytical solution Two dimensional numerical solution 0.8

1 0.8

ψ

0.6

ψ

0.6 0.4

0.4 0.2 2 1

0.2

0

0 1

0

0

0.2

0.4

0.6

0.8

1

0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

0.5

0

-0.5

-1 -1-2

ξ

η

η

(b)

(a)

Figure 2: (a) Comparison of one dimensional and two dimensional potential distribution for half channel and (b) two dimensional surface potential distribution when κH = 10, where ξ, η, ψ are scaled parameters.

4. Simplified mathematical model The electric field satisfying the Maxwell’s equations is governed by the electric potential term Φ∗ which can be written as Φ∗ = φ∗ + ψ ∗ , where φ∗ is the potential due to the external electric field and ψ ∗ is the induced potential due to the presence of EDL. The applied external electric field is assumed to be weak i.e.

∆φ∗ L

<<

|ζ| , κ−1

where ∆φ∗ is the

potential difference due the external electric field applied in the axial direction. Therefore, the change of the potential difference does not depend on the external electric field, that changes is found from the wall zeta potential. Thus the potential φ∗ can be written in the form φ∗ = φ0 − xEx (i.e. φ∗ = φ∗ (x)). The Poisson equation for the induced potential can

be rewritten in a modified form as,

ρe d2 ψ ∗ d2 ψ ∗ + = − dx2 dy 2 

(8)

The result for the induced potential is obtained numerically using Boltzmann distribution with Debye-Huckel approximation (which is appropriate for small co-ions and counter-ions concentration differences [Conlisk 2005] ) i.e. the electric potential is assumed to be small compared to the thermal energy of the ions (|zeψ ∗ | < |kB Tab |) and is presented in Fig. 2. From the figure it is observed that the induced potential can be treated as a function of y only i.e. ψ ∗ = ψ ∗ (y) [Probstein 2003].

The flow is assumed to occur in a straight two dimensional microchannel. The following assumptions are considered: 1. Fluid is non-Newtonian and incompressible. 2. Fluid is laminar and steady state. 3. Physical and thermo-physical properties are constant. 4. Flow is hydrodynamically and thermally fully developed. 5. Wall zeta potential is constant. 6. Electrolyte is symmetric. 7. Ions are point charges. 8. Permittivity of the fluid is constant and not affected by overall electric field strength. 9. The Debye-Huckel linearization principle holds. Under the above assumptions, the governing equations reduced to, d2 ψ ∗ ρe = − dy 2  " # du∗ n−1 du∗ d dP + ρe Ex = 0 kp − dy dy dy dx n−1  ∗ 2 2 ∗ ∗ du∗ ∂ T ∂qr du ∂T = kT h + k + σEx2 . − ρcp u∗ p 2 ∂x ∂y dy dy ∂y

(9) (10) (11)

In the fully developed region the constant heat flux condition considered along the wall yields

∂T ∗ ∂x

=

∂Tw ∂x

=

∂Tm ∂x

= constant, where Tw is the initial surface temperature and Tm is

the bulk temperature. The centerline symmetry is considered and the boundary conditions are used as dψ ∗ = 0 at y = 0 and ψ ∗ = ζ at y = H dy du∗ du∗ = 0 at y = 0 and u∗ + b = 0 at y = H dy dy ∂T ∗ ∂T ∗ ∗ = 0 at y = 0 and T = Tw + kj = 0 at y = H, ∂y ∂y

(12) (13) (14)

where b denotes the the slip length or slip coefficient and kj is the temperature jump factor at the wall of the microchannel.

5. Analytical Estimation The flow governing equations represent a fully developed flow without any heterogeneity along the hydrophobic surface. An attempt has been made to find the analytical solutions based on different approximations. 5.1. Induced potential The Poisson Boltzmann Eq. (9) is approximated by the Debye-Huckel theory and is represented in the form, d2 ψ ∗ = κ2 ψ ∗ , 2 dy where κ =

(15)

p 2n0 z 2 e2 /kB Tav is the Debye-Huckel parameter and κ−1 is called the Debye

length (λD ) which indicate the EDL thickness of the channel.

The solution of Eq. (15) subject to the boundary condition (12) is given by ψ ∗ (y) = ζ

cosh(κy) . cosh(κH)

(16)

In dimensionless form Eq. (16) can be written as, ψ(η) = where ψ(η) =

ψ ∗ (y) , ζ

η =

y H

and ω =

λD H

cosh(η/ω) , cosh(1/ω)

(17)

is the electroosmotic parameter (ratio of Debye

length to half channel height). 5.2. Axial velocity (u) We can rewrite the momentum Eq. (10) in terms of exclusively of the velocity field in the transverse coordinates which leads to "  n−1 ∗ # du∗ dP du d kp − − + ρe Ex = 0. dy dy dy dx

(18)

The solution of Eq. (18) subject to the boundary condition (13) yields, 1   1 κEx ζ n sinh(κy) 1 dP n du∗ =− − + y dy kp cosh(κH) κEx ζ dx

(19)

The dimensionless form of Eq. (19) is given by,  n1  1 sinh(η/ω) du =− + 2Γωη , dη nω cosh(1/ω)

(20)

where Γ =

H 2 dP 2ζEx dx

velocity and uHS

is the dimensionless pressure gradient, u = uuHS is the dimensionless   n1 1−n ζE x , which is known as the generalized Helmholtz− kp = nκ n ∗

Smoluchowski velocity distribution. The corresponding dimensionless boundary condition is u+β where β =

b H

du = 0 at η = 1, dη

(21)

is the slip parameter.

Solution of Eq. (20) with the prescribed boundary condition (21) yields,  n1     n1 Z 1  sinh(χ/ω) 1 u(η) = + 2Γωχ dχ . β tanh(1/ω) + 2Γω + nω cosh(1/ω) η

(22)

Equation (22) suggests that the integration can be possible only for specific values of the power law index n (or flow behavior index) such as n = 1,

1 1 , , 2 3

etc.

For n = 1, one obtains 

   cosh(η/ω) 2 + Γ 1 − η + A, u(η) = 1 − cosh(1/ω)

(23)

  β A= tanh(1/ω) + 2Γω . ω

(24)

where,

For n = 21 , one obtains   2 ω sinh(2η/ω) − 2η 4Γω 2 u(η) = B− η cosh(η/ω) − − ω cosh(1/ω) 4 cosh2 (1/ω)   4 2 2 3 ω sinh(η/ω) − Γ ω η , 3

(25)

where,    2 4 2 2 ω sinh(2/ω) − 2 2 + 4Γω 1 − ω tanh(1/ω) + Γ ω + β tanh(1/ω) + 2Γω . (26) B= 3 4 cosh2 (1/ω) and for n = 31 , one obtains   3ωΓ ω cosh(3η/ω) − 9ω cosh(η/ω) 3 − C− 2ωη sinh(2η/ω) − u(η) = ω 12 cosh3 (1/ω) 4 cosh2 (1/ω)   12ω 3 Γ2 2 2 ω cosh(2η/ω) − 2η − η 2 cosh(η/ω) − 2ωη sinh(η/ω) + cosh(1/ω)   2 3 3 4 (27) 2ω cosh(η/ω) − 2ω Γ η

where,   3ωΓ ω cosh(3/ω) − 9ω cosh(1/ω) 2 + C= 2ω sinh(2/ω) − ω cosh(2/ω) − 2 + 12 cosh3 (1/ω) 4 cosh2 (1/ω)    3 3 2 2 3 3 12ω Γ 1 − 2ω tanh(1/ω) + 2ω + 2ω Γ + β tanh(1/ω) + 2Γω . (28) The important physical quantity of interest is the friction factor, which is defined as fr = R 2τw 1 ,where u = udAc is the bulk mean fluid velocity over the channel cross section. av 2 ρu Ac Ac av

The skin-friction Cf can be defined as [Chen 2012],   n+1 4 tanh(1/ω) + 2Γω  n , Cf = fr Re = R1 2 nω 0 u(η)dη

(29)

where, Re = ρu2−n av Dh /kp is represents the Reynolds number and Dh = 4H is the hydrodynamic diameter of the microchannel. The skin-friction for n = 1 is 

24 tanh(1/ω) + 2Γω Cf =





ω 3 + A + 2Γ − ω tanh(1/ω)

,

(30)

for n = 21 , 



4 tanh(1/ω) + 2Γω   Cf =  ω 2 cosh(2/ω)−2 B − 8 cosh2 (1/ω) − 4Γω 3 tanh(1/ω) − 2ω − 31 Γω+

 1/2 ω2 2 , 1 − 64Γω cosh(1/ω) 8 cosh2 (1/ω)

(31)

and for n = 31 , √ 3 Cf = 

C−

ω 2 sinh(3/ω)−27ω 2 sinh(1/ω) 36 cosh3 (1/ω) 2

− 12Γ ω

4





32 tanh(1/ω) + 2Γω   3ωΓ − 36 cosh ω 2 cosh(2/ω) − ω 3 sinh(2/ω) − 32 3 (1/ω)





2 tanh(1/ω) − 4ω + 6ω tanh(1/ω) − Γ3 ω 3 5 4

1/3

. (32)

5.3. Thermal transport equation The radiative heat flux qr is approximated by the Rosseland approximation [Shit et al. 2016] and can be written as qr = −

4σ1 ∂T ∗ 4 , 3k1 ∂y

(33)

where σ1 is the Stefan Boltzmann constant and k1 is the mean absorption coefficient. Assuming that the temperature difference within the channel is small such that T ∗ 4 can be expressed in terms of linear combination of the temperature difference. Therefore T ∗ 4 can be expressed by a Taylor’s series form around the wall temperature Tw as, T ∗ 4 = Tw4 + 4Tw3 (T ∗ − Tw ) + 6Tw2 (T ∗ − Tw )2 ...

(34)

Since the temperature differences within the channel are small, so neglecting the higher powers of (T ∗ − Tw ), the expression for T ∗ 4 is reduced in the form, T ∗ 4 = Tw4 + 4Tw3 (T ∗ − Tw ) .

(35)

With the aid of Eq. (33) and Eq. (35) with Eq. (11) we can rewrite the energy equation as,  n−1  ∗ 2 ∂ 2T ∗ du∗ du + σEx2 , = kT h (1 + N r) 2 + kp − ρcp u ∂x ∂y dy dy ∗ ∂T

where N r =

3 16σ1 Tw 3k1 kT h

∗

(36)

is the thermal radiation parameter.

Here the volumetric heat source is obtained due to the presence of Joule heating and the viscous dissipation components. It is observed that for a wall zeta potential of −25mV

with a channel height larger than 0.1µm, the contribution of viscous dissipation is much less than the Joule heating component [Chen 2009] and in that scenario we can neglect the viscous dissipation term. Also we can neglect the Steric effects for small wall zeta potential [Jimenez et al. 2017, Ranjit and Shit 2017]. The overall performance of the energy balance term on the flow field leads to, −(1 + N r)q + σEx2 H ∂T ∗ = , ∂x ρcp uHS H

(37)

is the wall heat flux. where q = −kT h ∂T ∂y ∗

Combining Eq. (36) and Eq. (37) results the following dimensionless equation; (1 + N r)

∂ 2T + S = (S − N r − 1)u, ∂η 2

(38)

where S =

σEx2 H q

is the dimensionless volumetric heat generation or heat absorption i.e.

either (S > 0) or (S < 0) due to the effect of Joule heating. The dimensionless form of temperature T can be written as T =

T ∗ −Tw . qH/kT h

The boundary conditions in dimensionless form are, ∂T ∂T = 0 at η = 0 and T = St at η = 1, ∂η ∂η where St =

kj H

(39)

is the dimensionless temperature jump factor.

Solution of Eq. (38) subject to the boundary conditions (39) yields,   Z 1Z % S 2 (1 + N r)T (η) = u(χ)dχd% + 1 − η − (S − N r − 1) 2 0 η   Z 1 u(η)dη − S . St (S − N r − 1)

(40)

0

The temperature distribution for the flow behavior index n = 1, is given by

where,

     1 S 2 1 − η + (S − N r − 1) 1 + Γ + A η2 − (1 + N r)T (η) = 2 2  cosh(η/ω) Γ 4 η −ω + D, 12 cosh(1/ω) 

(41)



  2Γ D = St (S − N r − 1) 1 − ω tanh(1/ω) + +A −S + 3 (N r + 1 − S) (6A + 5Γ − 12ω 2 + 6). 12 For n = 21 , the temperature distribution     2 B 2 3ω 3 sinh(2η/ω) − 4η 3 S 2 1 − η + (S − N r − 1) η − (1 + N r)T (η) = 2 ω 2 48cosh2 (1/ω)     4Γω 4 1 2 2 5 − η cosh(η/ω) − 3ω sinh(η/ω) − Γ ω η + Eη − F + G, cosh(1/ω) 15

(42)

(43)

where,

  ω 2 1 − 64Γω cosh(1/ω) E= 4 cosh2 (1/ω)     1 2 2 2 B 3ω 3 sinh(2/ω) − 4 4 − 4Γω 1 − 3ω tanh(1/ω) − Γ ω + E − F = ω 2 15 48 cosh2 (1/ω)      2 2 ω cosh(2/ω) − 2 3 G = St (S − N r − 1) − 4Γω tanh(1/ω) − 2ω − B− ω 8 cosh2 (1/ω)    1 2 2 (44) Γ ω +E −S . 3

and for n = 13 , one obtains     3 C 2 ω 3 cosh(3η/ω) − 81ω 3 cosh(η/ω) S 2 η − (1 + N r)T (η) = 1 − η + (S − N r − 1) 2 ω 2 108cosh3 (1/ω)    3 12ω 5 Γ2 3ω 4 1 4 ω 3ωΓ η sinh(2η/ω) − cosh(2η/ω) − η − η 2 cosh(η/ω) − − 4 6 cosh(1/ω) 4 cosh2 (1/ω) 2    1 3 3 6 2 (45) 6ωη sinh(η/ω) + 12ω cosh(η/ω) − ω Γ η − I + J 15 where,   3 3 C ω 3 cosh(3/ω) − 81ω 3 cosh(1/ω) 3ωΓ ω I= − − sinh(2/ω) − 2 ω 2 108cosh3 (1/ω) 4 cosh (1/ω) 2    1 12ω 5 Γ2 3ω 4 2 cosh(2/ω) − − cosh(1/ω) − 6ω sinh(1/ω) + 12ω cosh(1/ω) − 4 6 cosh(1/ω)  1 3 3 ω Γ 15   3 ω 2 sinh(3/ω) − 27ω 2 sinh(1/ω) − J = St (S − N r − 1) C− ω 36 cosh3 (1/ω)    3ωΓ 2 4 2 3 − 12ω tanh(1/ω) − 4ω + ω cosh(2/ω) − ω sinh(2/ω) − 3 4 cosh2 (1/ω)    2 3 3 2 6ω tanh(1/ω) − Γ ω −S . 5 

(46)

The heat transfer rates are expressed as Nusselt number which is given by Nu =

q(4H) 4 hDh = =− , kT h kT h (Tw − Tm ) θm

(47)

where h is the heat transfer coefficient and θm is the dimensionless mean temperature which is defined by, θm =

R1 0

u(η)T (η)dη . R1 u(η)dη 0

(48)

The analytical expression for dimensionless mean temperature for a Newtonian fluid (n = 1) is θm =

(A + 1)A1 +

4ΓS 15

+ (S − N r − 1)ΓA2 + A3 − (S − N r − 1)A4 − Dω tanh(1/ω) (49) A5

where, A1 A2 A3 A4

A5

  S A + 1 9Γ 2 = + (S − N r − 1) + − ω tanh(1/ω) + D 3 6 60 A + 1 13Γ − − 2ω 3 + 2ω 4 tanh(1/ω) = 15 210  2DΓ 2 = − Sω 1 − ω tanh(1/ω) 3    Γω (1 + A + Γ)ω 2 tanh(1/ω) − 2ω + 2ω tanh(1/ω) − tanh(1/ω) − 4ω = 2 12    ω 2 3 4 + 12ω tanh(1ω) − 24ω + 24ω tanh(1/ω) − 1 + ω sinh(2/ω) 4 cosh(1/ω)   2Γ − ω tanh(1/ω) . = (1 + N r) 1 + A + 3

The expressions of θm for n =

1 2

and n =

1 3

are quite large and complicated, so these are

not presented here. 6. Results and Discussions In the present study a binary fluid is driven by the combined pressure driven EOF inside a hydrophobic microchannel. The considered channel is fabricated from a silicon material where the wall exhibit a negative zeta potential. The characteristics of fluid flow and heat transfer are analyzed both analytically and numerically with the variation of dimensionless pressure gradient (0 ≤ Γ ≤ 10), slip parameter (0 ≤ β ≤ 0.05), power law index (0 <

n < 2), Joule heating parameter (−10 ≤ S ≤ 10), electrical conductivity (0 < σ ≤ 1)

and Debye length. The integrations involved in Eqs. (22), (29), (40) and (47) are carried

out using a composite trapezoidal rule which gives an approximate solution. The step size of the space variable is considered to be 0.0001 for the numerical solution. The working fluid is considered as a power-law fluid containing binary symmetric electrolytes such as KCl or NaCl. The employed parameters are: the wall zeta potential ζ = −25mV , channel

length L = 0.005m, half channel height H = 40µm, electron charge e = 1.6 × 10−19 C, electric field strength Ex = 2 × 104 V m−1 , Boltzmann constant kB = 1.3805 × 10−23 J/K,

absolute temperature Tab = 300K, valence of ion z = 1, permittivity of the fluid  = 5.3 × 10−10 C/V m, wall heat flux q = 1500W/m2 , electrical conductivity σ = 10−3 S/m to

1.0S/m, specific heat cp = 3760J/kgK, thermal conductivity kT h = 0.613W/mK, Stefan Boltzmann constant σ1 = 5.67 × 10−8 W/m2 K 4 and absorption coefficient k1 = 10−1 m−1 .

1

1

ω = 0.1 0.8

ω = 0.2

0.8

Γ = 0, 0.3, 0.5

ω=1

0.6

Analytical solution Numerical solution

η

η

0.6

Analytical solution

0.4

0.4

Numerical solution 0.2

0

0.2

0

0.5

1

u

1.5

0

2

(a)

0

0.2

0.4

u

0.6

0.8

1

(b)

Figure 3: Velocity comparison of analytical solution with numerical solution for (a) different values of pressure gradient when ω = 0.3, β = 0 and n =

1 2

(b) different values of ω when Γ = 0, β = 0 and n = 21 .

6.1. Velocity Distribution The velocity distribution along the channel is computed numerically as well as analytically for a non hydrophobic microchannel (i.e. β = 0) and is presented in Fig. 3(a). It is clearly emphasized that, for a plane microchannel (i.e. β = 0) Eq. (22) reduces to the same expression for axial velocity as derived by Chen [Chen 2012]. In this figure the dimensionless velocity profiles are compared for different values of the dimensionless pressure gradient with the flow behavior index n = 12 . From this figure it is observed that composite Trapezoidal rule provides a good comparative solution with the analytical results. For a purely electrosmotic flow the fluid flows vary slowly as compared to the pressure driven flow, but the flow velocity is increased if the flow is driven due to the combined effect of electroosmosis and pressure. It is clearly visualized that the flow is increased with the increment of dimensionless pressure gradient as expected, and is presented in Fig. 3(a). But the flow speed is further increased for the flow over hydrophobic walls. These results are demonstrated later. A comparison of velocity distribution is also made for a variation of EDL thickness and is presented in Fig. 3(b). The results are obtained for a purely one dimensional EOF (i.e. Γ = 0) in a microchannel and it is observed that the velocity profile achieve a top-hat profile when the values of ω is decreased. This is due to thin EDL, which drags all the fluid layers due to the effect of viscosity. For the absence of pressure gradient

electrically driven force is applied to move the fluid due the movement of EDL. The velocity distribution of combined EOF with a pressure driven flow for power law fluids are obtained both analytically and numerically from Eq. (22). The numerical solutions are in good agreement with the analytical solution for pseudoplastic fluids (n < 1). The axial velocity distribution for different values of slip parameters (β) and flow behavior index (n) are shown in Figs. 4 (a), (b) & (c) for a fixed value of Γ = 0.5 and ω = 0.3. Figs. 4 (a), (b) and (c) represent the velocity profiles for β = 0.001, β = 0.01 and β = 0.05 respectively. From this figures it is found that the velocity is increasing with decreasing values of the power law index ‘n’ for each values of Γ. For a dilatant fluid (n > 1), the shear thickening nature enhances the maximum velocity at the center line and produces a parabolic profile. Conversely, for a pseudo plastic fluid (n < 1), the shear thinning behavior reduces the optimum velocity but increases the velocity near the channel wall, leading to a plug flow nature. It is clearly observed from the figures that the flow profiles show opposite behavior along the center line of the channel in case of shear thinning and shear thickening fluid. From the figures it is noticeable that the velocity slip bears the potential at the channel wall to alter the axial velocity distribution significantly. Also it is observed that the velocity is increased with the increase of the slip factor for a fixed dimensionless pressure gradient. Hence it is concluded that the hydrodynamic slip factor provides a significant contribution to the velocity distribution along the wall. For a particular value of the electric field strength or pressure gradient the velocity slip at the walls of the electrokinetic microchannel has a significant effect on the volumetric flow RH rate. The average axial velocity is determined by using the relation; uav = H1 0 u(y)dy

and the volumetric flow rate is defined as Qf = Huav . The variation of volumetric flow rate Qf due to pressure gradient Γ for various values of slip parameters are shown in Figs.

4 (d), (e) & (f) which represent the volumetric flow rate for different power law indices n = 1,

1 2

and 31 . From Fig. 4(d) it reveals that the volumetric flow rate is a linear function

of pressure gradient in case of a Newtonian fluid (n = 1). But this variation tends to a nonlinear function of Γ for a pseudo plastic fluid (n < 1) and the flow rate is increased with a decrease in the power law index, because the shear thinning effect reduces the maximum velocity for which the average velocity is gradually increased. In these figures the velocity slip parameter has a significant contribution for enhancing the flow rate distribution. When the slip parameter is increased the wall slip impact is increased and the flow velocity is increased and the increment is 1.2 times for Newtonian fluid (n = 1) and 1.5 times for pseudo plastic fluid (n = 13 ) when the slip parameter is increased from β = 0 to β = 0.05.

1

1

n = 1/3

n = 1/3

0.8

1

0.8

n = 1/2 n=1

n = 1.2

n = 1.2

n = 1/2

n=1 0.6

0.4

0.4

0

Analytical solution 0.2

Numerical solution

0

0.5

n = 1.2

0.4

Analytical solution 0.2

n=1

η

η

0.6

η

0.6

n = 1/3

0.8

n = 1/2

1

1.5

2

0

2.5

0.2

Numerical solution

0

0.5

1

1.5

2

0

2.5

Analytical solution Numerical solution 0

0.5

1

1.5

2

u

u

u

(a)

(b)

(c)

2.5

3

3.5

7 1.6

3.5

β=0 β = 0.001 β = 0.01 β = 0.05

1.5 1.4

β=0 β = 0.001 β = 0.01 β = 0.05

3

1.3

β=0 β = 0.001 β = 0.01 β = 0.05

6

5

2.5

Qf

Qf

Qf

1.2 4

1.1 2 3

1 0.9

1.5 2

0.8 1

0.7

1 0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

Γ

Γ

Γ

(d)

(e)

(f)

0.8

1

Figure 4: Velocity distribution for different values of the power law index when (a) β = 0.001 (b) β = 0.01 & (c) β = 0.05. Variation of volumetric flow rate with dimensionless pressure gradient for different values of slip parameter when (d) n = 1 (e) n =

1 2

& (f) n = 13 .

When the power law fluid is passing through a microchannel the frictional force is developed on to the surface walls which resists the further movement of the fluid into the channel and this friction is called skin friction. In the present study the variation of friction coefficient with the flow behavior index for different values of slip parameters are shown in Figs. 5 (a), (b) & (c), when ω = 0.2. The friction coefficient variations are also shown for various values of the dimensionless pressure gradient. In general for a fixed values of β the friction coefficient is increased with an increase of the power law index (n) for any values of Γ. Since shear thinning fluids are inherently more amenable to the flow field, such type of fluid exhibit less friction as compared to the shear thickening fluid. In case of pure EOF (i.e. Γ = 0) in a plane microchannel (i.e. β = 0) produces the maximum friction coefficient which is clearly shown in Fig. 5(a) but the friction coefficient is decreasing with the increase of either Γ or β which are assisted factors for the fluid flow. It is obvious to conclude that an increase in β enhances the slip velocity along the wall. Figure 5((d), (e) & (f)) implies the results in terms of a flow enhancement factor Ef b where Ef b measures the average flow velocity in the microchannel when β > 0 (i.e. hydrophobicity is present at the walls) divided by the average flow velocity under the same conditions when β = 0 (i.e. without hydrophobicity). Therefore, Ef b measures the increment in average flow velocity by using a hydrophobic structure in the channel walls. The results provide to study the influence of the slip parameter on the flow field. The flow enhancement factor is plotted as a function of the flow behavior index for various values ω and is presented in Figs. 5 (d), (e) & (f) in case of a purely EOF (i.e. Γ = 0). Since shear thinning fluids have more amenable properties due to which the flow enhancement factor gradually decreases with the increment of flow behavior index (n). As the EDL thickness increases i.e. for weak electrolytes the flow in the bulk region is less pronounced and Ef b is reduced due to the combined effects of movement of diffuse layer and diffusion property of the fluid. But when the thickness of the EDL is very thin, the opposite ions are rigidly bounded by the counter ions, fluid movement is more dominated by EDL compared to the diffusion property produce due to the movement of the bulk fluid. The impact of the slip parameter is clearly visualized from Figs. 5 (d), (e) & (f). One sees that the average flow velocity is approximately increased two times and six times for ω = 0.3 when the slip parameter increases from 0.001 to 0.05. It is obvious that the average flow velocity is increased by increasing the pressure gradient but the Navier-velocity slip is more sensitive in case of pure EOF to get a maximum average flow velocity as expected.

300 200

β=0 β=0.001 β=0.01 β=0.05

150

Cf

100

Cf

200

β=0 β=0.001 β=0.01 β=0.05

150

Cf

250

150

β=0 β=0.001 β=0.01 β=0.05

100

100

50 50

50

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

0.2

1

1.2

(a)

(b)

(c)

ω = 0.05 ω = 0.1 ω = 0.2 ω = 0.3

1.8

4

1.04

1.4

3

1.02

1.2

2

0.9

1.2

1.5

1.8

1.4

1.6

1.8

ω = 0.05 ω = 0.1 ω = 0.2 ω = 0.3

5

Efb

1.6

Efb

1.06

0.6

0.8

n

ω = 0.05 ω = 0.1 ω = 0.2 ω = 0.3

0.3

0.6

n

1.08

1

0.4

n

Efb

0.2

1

0.3

0.6

0.9

1.2

n

n

(d)

(e)

1.5

1.8

1

0.3

0.6

0.9

1.2

1.5

1.8

n

(f)

Figure 5: Variation of friction coefficient with power law index (n) for different values of β when (a) Γ = 0 (b) Γ = 1 & (c) Γ = 5. Variation of flow enhancement factor with flow behavior index n for different average velocity at β=0.01 average velocity at β=0.05 average velocity at β=0.001 (e) & (f) values of ω (d) average velocity at β=0 average velocity at β=0 average velocity at β=0 when Γ = 0.

1

1 n = 1/3 n = 1/2

0.8

n=1

n = 1/3

0.8 n=1

n = 1/2

0.6

0.6

η

η

Γ = 0.3

0.4

Analytical solution Numerical solution

β = 0.001

0.4

β = 0.05 Analytical solution Numerical solution

0.2

0.2 Γ = 0.5

0 -0.3

-0.2

-0.1

0

0

0.1

-0.3

T

-0.2

-0.1

0

0.1

T

(a)

(b)

Figure 6: Validation of analytical and numerical results of temperature distribution (a) for different values of Γ at β = 0.001 and (b) for different values of β at Γ = 0 with the flow behavior index n when ω = 0.2, S = 2.0, St = 0.3, N r = 0.

6.2. Temperature Distribution The isothermal variation across the channel is computed along the upper half of the channel by assuming symmetry at the centreline. The analytical and numerical solutions of the temperature distribution are compared and represented in Figs. 6(a) and 6(b) for three different flow behavior indices i.e. n = 13 ,

1 2

and 1. It is clearly observed that for different

values of Γ and β the numerical solutions agree excellently with the analytical solutions. Figure 6(a) reveals that the dimensionless temperature profiles show a large variation when dimensionless pressure gradient is increased and the temperature lines move towards the positive direction along the channel wall. Since the temperature difference between the bulk region and surface wall is increased with the increasing value of dimensionless pressure gradient. The temperature lines achieve maximum values at the centerline and decrease towards the channel wall for Γ = 0.3. But there is an interesting observation that the maximum temperature for Γ = 0.5 and n =

1 3

appears at the channel wall, since the

increasing value of Γ reduced the maximum temperature along the centerline of the channel. A similar behavior also occurred when β is increasing i.e. temperature is increasing and vice versa. This is presented in Fig. 6(b). The radiation effect on heat transfer are the electromagnetic radiation effects generated

1

1

1

0.8

0.8

0.8

Nr = 0 Nr = 0.5 Nr = 1 Nr = 1.5 Nr = 2.0

0.4

0.6

0.4

0.2

0

0.6

η

Nr = 0 Nr = 0.5 Nr = 1 Nr = 1.5 Nr = 2

η

η

0.6

0.2

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0 -0.4

Nr = 0 Nr = 0.5 Nr = 1 Nr = 1.5 Nr = 2

0.4

0.2

-0.3

-0.2

-0.1

0

0.1

0.2

0 -0.4

-0.3

-0.2

-0.1

0

T

T

T

(a)

(b)

(c)

0.1

0.2

Figure 7: Temperature distribution for different values of thermal radiation parameter when ω = 0.1, Γ = 1.0, β = 0.001, S = 2.0, St = 0.3 and (a) n = 0.5 (b) n = 1.0 (c) n = 1.5.

due to the thermophoretic transport of charged particles which is displayed in Figs.7 (a), (b) and (c) for different power law indices i.e. n = 0.5, 1 & 1.5 respectively. From these figures it is observed that, for a fixed power law index, the temperature increases with the increment of the radiation parameter (N r) which occurs along the centerline of the channel but a reversed trend is followed along the channel wall. This reversed behavior of the temperature happens due to the presence of the temperature jump factor which controls the wall temperature. Also, the temperature increment depends on the mean absorption coefficient (Eq. 33) which produces the maximum variation of temperature along the centerline of the channel. This effect is most efficiently used for thermal radiation effect on blood flow in medical treatment of heat therapy, that produces heat on the affected areas through blood capillaries [Prakash and Makinde 2011]. This maximum variation of temperature is obtained in case of a shear thinning fluid which is presented in Fig.7(a). The Joule heating parameter is the most important effective parameter for heat transfer and its influence on thermophoretic transport is presented in Figs. 8 (a), (b), (c), (d), (e) & (f). Figs. 8 (a), (b) & (c) reveals the temperature distribution for pseudoplastic fluid (n = 0.5), Newtonian fluid (n = 1) and dilatant fluid (n = 1.5) in the presence of radiation effect of pressure assisted flow for different values of the Joule heating parameter (S). These values of temperature in the plots are computed in the absence of a temperature jump factor, where S > 0 indicates heat generation and S < 0 heat absorption. It is clearly visualized from Fig. 8 (a), (b) & (c) that the increment values of S, provides a higher value of the temperature difference for all three types of fluid (i.e. shear thinning, Newtonian &

1

1

S = -10 S = -5 S = -1 S=0 S=1 S=5 S = 10

0.8

1

S = -10 S = -5 S = -1 S=0 S=1 S=5 S = 10

0.8

0.8

η

0.6

η

0.6

η

0.6

0.4

0.4

0.4

0.2

0.2

0.2

-0.5

0

0.5

1

1.5

2

2.5

0 -0.5

3

0

0.5

1

1.5

2

0

2.5

0.5

1

T

T

(a)

(b)

(c)

1

S = -5 S = -1 S=0 S=1 S=5 S = 10

S = -5 S = -1 S=0 S=1 S=5 S = 10

0.8

0.4

0.4

0.2

0.2

0.2

0

0.2

0.4

0.6

0.8

1

S = -5 S = -1 S=0 S=1 S=5 S = 10

0.8

0.4

-0.8 -0.6 -0.4 -0.2

2

0.6

0.6

η

η

0.6

1.5

1

1

0.8

0

0

T

η

0

S = -10 S = -5 S = -1 S=0 S=1 S=5 S = 10

0

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

0

-0.6

-0.4

-0.2

0

0.2

T

T

T

(d)

(e)

(f)

0.4

0.6

0.8

Figure 8: Variation of temperature profiles with different values of Joule heating parameter (S) when ω = 0.1, Γ = 1.0, β = 0.001, N r = 2.0 & (a) n = 0.5 (b) n = 1.0 (c) n = 1.5 for St = 0.0 and (d) n = 0.5 (e) n = 1.0 (f) n = 1.5 for St = 0.3.

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

St = 0 St = 0.1 St = 0.3 St = 0.5

0.2

0

η

1

η

1

η

1

-0.6

-0.4

-0.2

0.4

St = 0 St = 0.1 St = 0.3 St = 0.5

0.2

0

0.2

0.4

0.6

0

-0.6

-0.4

-0.2

0

St = 0 St = 0.1 St = 0.3 St = 0.5

0.2

0.2

0.4

0.6

0

-0.4

-0.2

0

T

T

T

(a)

(b)

(c)

0.2

0.4

0.6

Figure 9: Temperature profiles for different values of St when ω = 0.1, Γ = 1.0, β = 0.001, S = 2.0 and N r = 2.0, (a) n = 0.5 (b) n = 1.0 (c) n = 1.5.

shear thickening). Since the temperature difference between the channel surface and the bulk region of the fluid decreases, i.e. the relative magnitude of Joule heating and constant surface heat-flux is increased. The influence of Joule heating is more prominent for higher values of S of all three types of fluid and the maximum temperature is observed along the centerline. The temperature profiles change its sign from negative to positive when S increases negatively. This signifies that for certain values of S the temperature variation will be zero, which will generate a singularity in the Nusselt number represented in Eq. (47). It can also be observed that

∂T ∗ ∂x

= 0 for S − N r = 1 (Eq. 37). In this case the

temperature profiles can be uniquely determined from Eq. (40) which is propotional to S (1 − η 2 ) − S 2

St and will represent a parabolic profile for all three types of fluid. In Figs.8

(d), (e) & (f) the result includes all the same parameters used in Figs. 8 (a), (b) & (c) with an extra parameter temperature jump factor (St = 0.3) to find the relative effect of St for different flow behavior indices, n = 0.5, 1 and 1.5 at particular values of S. From these figures it can be observed that the Joule heating parameter enhances the temperature upto a certain height of the channel and then decreases in all three type of fluids (i.e. shear thinning, Newtonian and shear thickening). This type of behavior is controlled by the temperature jump factor which decreases the temperature along the channel wall. One may also observe that the temperature profiles change its sign from negative to positive as the Joule heating parameter is increased negatively. This type of behavior is responsible for the heat absorption coefficient. Figures 9 (a), (b) and (c) represent the temperature distribution for different power

law indices n = 0.5, 1 and 1.5 at selected values of St . It is observed that the temperature distribution gradually decreases throughout the channel as the temperature jump factor (St ) increases for all three types of fluid. From the figures, it is visualized that the temperature variation is dependent on the flow behavior index (n). In case of St = 0 the temperature distribution is positive throughout the entire section of the channel but this variation is slowly shipped towards the negative temperature with increasing St and the maximum temperature is obtained along the centerline of the channel. From all the figures, it is concluded that St plays a vital role to control the temperature distribution along the channel wall. The Nusselt number generally depends on the mean temperature (see Eq. (47)) on the flow field and its variation with the dimensionless pressure gradient (Γ). For a constant surface heat flux boundary conditions, it is represented in Fig. 10 (a). In general the magnitude of the Nusselt number decreases with the increment of Γ for all three types of fluid. This behavior depends on the temperature difference in the flow field which increases as the dimensionless pressure gradient increases and it is easily understandable from Fig. 6(a). The rate of heat transfer is maximum for the shear thinning fluid and minimum for the shear thickening fluid when Γ < 3 but it is conversed for Γ > 3. Therefore, at Γ = 3 the heat transfer rate does not depend on the flow behavior index. The absolute Nusselt number variation with the flow behavior index (n) are presented in Fig. 10 (b) for different slip parameters (β) and it is observed that the rate of heat transfer gradually decreases with the increment of power law index. Also, a significant change occurs according to the slip parameter, as the temperature variation depends on the slip parameter which is presented in Fig.6(b) i.e. the temperature difference between the bulk region and the surface increases as slip parameter increases. For a fully developed flow, the Nusselt number is mainly dependent on the dimensionless Joule heating parameters (S) as well as on the thermal radiation parameter (N r) for the uniform heat flux boundary condition. The variation of the absolute Nusselt number with the flow behavior indices for different values of the Joule heating parameter is presented in Fig. 10 (c), it is found that the absolute Nusselt number strongly depends on S and indicating more significant contribution on heat transfer. The effect of Joule heating on heat transfer rate is recognized to be maximum and minimum at n = 0.3 and n = 1.8 respectively for S < 2.5. The next immediate value of S giving maximum heat transfer rate is at n = 0.2. A weak variation of heat transfer is observed in case of a shear thickening fluid up to S = 1, since the temperature difference is decreased between the surface and

60 50

β=0 β = 0.001 β = 0.01

50

n = 0.8 n=1 n = 1.2

45

|Nu|

|Nu|

40

40

30

35

20

30

10 0

2

4

6

8

10

0.2

0.4

0.6

0.8

1

Γ

1.2

1.4

1.6

1.8

n

(a)

(b)

120

S = -3 S = -1 S=0 S=1 S=2 S = 2.5 S=3

|Nu|

60

100

n = 0.8 n=1 n = 1.2

80

|Nu|

80

60 40

40 20

0.2

0.4

0.6

0.8

1

n

(c)

1.2

1.4

1.6

1.8

20

2

4

6

8

10

Nr

(d)

Figure 10: Variation of Nusselt number with (a) dimensionless pressure gradient for different flow behavior index (n) (b) flow behavior index (n) for different values of slip parameter (c) flow behavior index (n) for different values of Joule heating parameter & (d) thermal radiation parameter (N r) for different flow behavior index.

the fluid which is clearly recognized from Fig. 8. It is also observed that the temperature profile becomes negative when S increases positively and the heat transfer is decreasing for all three types of fluid, however the magnitude of the Nusselt number is a function of power law index as well as Joule heating parameter. Figure 10 (d) exhibits the magnitude of Nusselt number variations with the change of thermal radiation parameter for different values of the power law index. The magnitude of heat transfer rate is gradually decreasing with increasing the radiation parameter N r and this is clearly visualized from Fig. 10 (d). In case of a pseudo-plastic fluid the decay rate of the heat transfer is more prominent as compared to other fluids, because the temperature is decreasing with an increasing flow behavior index (see Fig. 7). The maximum variation of the Nusselt number is found for 0 < N r ≤ 6 and after that the variation is very weak or almost constant for all three types

of fluid.

7. Conclusions In the present study, fluid flow and heat transfer characteristics are analyzed for a combined pressure driven EOF of a power law fluid in a hydrophobic microchannel with Joule heating effects. A constant heat flux boundary condition with thermal radiation effects is considered to obtain a closed form analytical solution in terms of flow governing parameters, where as the numerical results are determined for a wide range of the flow governing parameters. The concluded remarks are as follows: 1. The fully developed nature of EOF is reported when the length scale ratio (ω) is very small i.e. the velocity achieves a top-hat profile, since the coefficient of friction is high. 2. The velocities are increasing with the increase of pressure gradient and slip parameter but behaves oppositely for power law index. The velocity slip bears the potential along the channel wall to extend the flow velocity, and the volumetric flow rate is increased with increasing Γ and β (which is shown in Fig. 4 (d), (e) & (f)). 3. The shear thinning fluid, which is intrinsically more amenable, reduces the friction coefficient in a flow field and the friction coefficient is rapidly decreased with the increment of dimensionless pressure gradient as well as the slip parameter (in Fig. 5 (a), (b) & (c)). 4. Increasing the slip parameter enhances the flow enhancement factor which reflects the ratio of average flow velocity through a hydrophobic channel with the average flow

velocity through a plane channel. The flow enhancement decays with the increment of the power law index. 5. For all the discussed cases in the result section the highest heat transfer is always found along the centerline of the channel and decreases towards the channel wall. However, thermal radiation effect produces a significant change in the temperature field along the boundary layer and controlling the temperature in the flow field by increasing/decreasing the radiation parameter. Also a crucial effect on temperature distribution along the channel wall is produced by the temperature jump factor. 6. The temperature distribution is proportional to

S (1 2

− η 2 ) − S St for the specific

values of Joule heating and radiation parameters which does not depend on the flow behavior index, pressure gradient and slip parameters. The temperature profile shows maximum variation in case of shear thinning fluid with the increase of the Joule heating parameter. 7. It is found that the magnitude of the Nusselt number varies inversely with dimensionless pressure gradient, slip parameter, power law index and radiation parameter, but increases with the increment of Joule heating parameter for a fixed power law index. Acknowledgement The authors are thankful to the reviewers for their valuable comments and suggestion, which enabled to present an improved version of the paper. Author Dr. A. K. Nayak is thankful to SERB (under DST financed by Project no. SR/S4/MS: 765/12) India for providing financial support during the preparation of this manuscript.

References Abe, M., Hiraoka, H., 1985. Localized hyperthermia and radiation in cancer therapy, Int. J. Radiat Biol. Relat. Stud. Phys. Chen. Med. 47, 347-359. Babaie, A., Sadeghi, A., Saidi, M. H., 2011. Combined electroosmotically and pressure driven flow of power-law fluids in a slit microchannel, J. Non-Newtonian Fluid Mech. 166, 792-798. Babaie, A., Saidi, M. H., Sadeghi, A., 2012. Heat transfer characteristics of mixed electroosmotic and pressure driven flow of power-law fluids in a slit microchannel, Int. J. Therm. Sci. 53, 71-79. Bera, S., Bhattacharyya, S., 2013. On mixed electroosmotic-pressure driven flow and mass transport in microchannels, Int. J. Eng. Sci. 62, 165-176. Burgreen, D., Nakache, F. R., 1964. Electrokinetic Flow in Ultrafine Capillary Slits, J. Phys. Chem. 68, 1084-1091. Cady, N. C., Stelick, S., Kunnavakkam, M.V., Batt, C. A., 2005. Real-time PCR detection of Listeria monocytogenes using an integrated microfluidics platform, Sensors Actuators B 107, 332-341. Chakraborty, S., 2007. Electroosmotically driven capillary transport of typical non-Newtonian biofluids in rectangular microchannels, Anal. Chim. Acta 605, 175-184. Chang, M. H., Homsy, G. M., 2005. Effects of Joule heating on the stability of time-modulated electroosmotic flow, Phys. Fluids, 17(7), 074107. Chen, C. H., 2012. Fully-developed thermal transport in combined electroosmotic and pressure driven flow of power-law fluids in microchannels, Int. J. Heat Mass Transfer 55, 2173-2183. Chen, C. H., 2009. Thermal transport charecteristics of mixed pressure and electroosmotically driven flow in micro- and nanochannels with Joule heating, ASME J. Heat Transfer 131, 022401. Conlisk, A. T., 2005. The DebyeHuckel Approximation: Its Use in Describing Electro-Osmotic Flow in Micro- and Nano-Channels, Electrophoresis 26, 1896-1912. Das, S., Chakraborty, S., 2006. Analytical solutions for velocity, temperature and concentration distribution in electroosmotic microchannel flows of a non-Newtonian bio-fluid, Anal. Chim. Acta 559, 15-24. Fung, Y. C., 1993. Biomechanics: Mechanical Properties of Living Tissues, second ed., Springer. Gascoyne, P., Satayavivad, J., Ruchirawat, M., 2004. Microfluidic approaches to malaria detection, Acta Trop. 89, 357-369. Gopalakrishnan, S., Thess, A., Weidmann, G., Lange, U., 2010. Chaotic mixing in a Joule-heated glass melt, Phys. Fluids, 22, 013101.

Hadigol, M., Nosrati, R., Nourbakhsh, A., Raisee, M., 2011. Numerical study of electroosmotic micromixing of non-Newtonian fluids, J. Non-Newtonian Fluid Mech. 166, 965-971. Hadigol, M., Nosrati, R., Raisee, M., 2011. Numerical analysis of mixed eletroosmotic/pressure driven flow of power-law fluids in microchannels and micropums, J. Colloids Surf. A: Physicochem. Eng. Aspects 374, 142-153. Hamann, H. J., Bickhardt, K., 1983. Radiation effects in blood used for the in vitro feeding of tsetse flies, Int. J. Appl. Radiat. Isot. 34, 851-854. He, Y., Shirazaki, M., Liu, H., Himeno, R., Sun, Z., 2006. A numerical coupling model to analyze the blood flow, temperature, and oxygen transport in human breast tumor under laser irradiation, Comput. Biol. Med. 36, 1336-1350. Horng, T. L., Lin, W. L., Liauh, C. T., Shih, T. C., 2007. Effects of pulsatile blood flow inlarge vessels on thermal dose distribution during thermal therapy, Med. Phys. 34, 1312-1320. Huang, K. R., Hong, Z. H., Chang, J. S., 2014. Microfluidic mixing on application of traveling wave electroosmosis, Eur. J. Mech. B Fluids 48, 153-164. Hunter, R. J., 1981. Zeta Potential in Colloid Science, Academic Press, New York. Ismail, A. H., Jaafar, Md. S., 2011. Interaction of low-intensity nuclear radiation dose with the human blood: Using the new technique of CR-39NTDs for an vitro study, Appl. Radiat. Isot. 69, 559-566. Jimenez, E., Escand´on, J., M´endez, F., Bautista, O., 2017. Combined viscoelectric and stric effects on the electroosmotic flow in a microchannel under induced high zeta potentials, Colloids and Surfaces A, 513, 221-233. Kang, Y. J., Lee, S. J., 2013. Blood viscoelasticity measurement using steady and transient flow controls of blood in a microfluidic analogue of Wheastonebridge channel, Biomicrofluidics 7, 054122. Khan, N., Sajid, M., Mahmood, T., 2015. Heat transfer analysis for magnetohydrodynamics axisymmetric flow between stretching disks in the presence of viscous dissipation and Joule heating, AIP Advances, 5, 057115. Lee, G. B., Fu, L. M., Lin, C. H., Lee, C. Y., Yang, R. J., 2004. Dispersion control in microfluidic chips by localized zeta potential variation using the field effect, Electrophoresis 25, 1879-1887. Levine, S., Marriott, J. R., Neale, G., Epstein, N., 1975. Theory of ElectrokineticFlow in Fine Cylindrical Capillaries at High Zeta-Potentials, J. Colloid Interface Sci. 52, 136-149. Li, P. C. H., Harrison, D. J., 1997. Transport, manipulation, and reaction of biological cells on-chip using electrokinetic effects, Anal. Chem. 69, 1564-1568.

Park, H. M., Kim, T. W., 2009. Extension of the HelmholtzSmoluchowski velocity to the hydrophobic microchannels with velocity slip, Lab Chip 9, 291-296. Paul, P. H., Garguilo, M. G., Rakestraw, D. J., 1998. Imaging of pressure and electrokinetically driven flows through open capillaries, Anal. Chem. 70, 2459-2467. Prakash, J., Makinde, O. D., 2011. Radiative heat transfer to blood flow through a stenotic artery in the presence of magnetic field, Lat. Am. Appl. Res. 41, 273-277. Probstein, R. F., 2003. Physicochemical Hydrodynamics: An Introduction, second ed., WileyInterscience. Ranjit, N. K., Shit, G. C., 2017. Joule heating effects on electromagnetohydrodynamic flow through a peristaltically induced micro-channel with different zeta potential and wall slip, Physica A, 482, 458-476. Ren, L., Li, D., 2001. Electrosmotic flow in heterogeneous microchannels, J. Colloid Interface Sci. 243, 255-261. Reuss, F. F., 1809. Charge-induced flow, proc. Imp. Soc. Nat. Mosc. 3, 327-344. Sadeghi, A., Saidi, M. H., 2010. Viscous dissipation effects on thermal transport characteristics of combined pressure and electroosmotically driven flow in microchannels, Int. J. Heat Mass Transfer 53, 3782-3791. Shit, G. C., Mondal, A., Sinha, A., Kundu, P. K., 2016. Electro-osmotic flow of power-law fluid and heat transfer in a micro-channel with effects of Joule heating and thermal radiation, Physica A 462, 1040-1057. Shit, G. C., Mondal, A., Sinha, A., Kundu, P. K., 2016. Electro-osmotically driven MHD flow and heat transfer in micro-channel, Physica A 449, 437-454. Sinha, A., Shit, G. C., 2015. Electromagnetohydrodynamic flow of blood and heat transfer in a capillary with thermal radiation, J. Magn. Magn. Mater. 378, 143-151. Stokes, J. R., Davies, G. A., 2007. Viscoelasticity of human whole saliva collected after acid and mechanical stimulation, Biorheology 44, 141-160. Vladisavljevic, G. T., Khalid, N., Neves, M.A., Kuroiwa, T., Nakajima, M., Uemura, K., Ichikawa, S., Kobayashi, I., 2013. Industrial lab-on-a-chip: Design, applications and scale-up for drug discovery and delivery, Adv. Drug Deliv. Rev. 65, 1626-1663. Wiedemann, G., 1852. Pogg. Ann. 87, 321. Wu, Y. H., Wiwatanapataphee, B., Hu, M., 2008. Pressure-driven transient flows of Newtonian fluids through microtubes with slip boundary, Physica A 387, 5979-5990.

Yan, D., Yang, C., Nguyen, N. T., Huang, X., 2007. Diagnosis of transient electrokinetic flow in microfluidic channels, Phys. Fluids, 19(1), 017114. Zade, A. Q., Manzari, M. T., Hannani, S. K., 2007. An analytical solutiopn for thermally fully devoloped combined pressure-electroosmotically driven flow in microchannels, Int. J. Heat Mass Transfer 50, 1087-1096. Zhao, C., Yang, C., 2013. Electroosmotic flow of non-Newtonian power law fluids in a cylindrical microchannel, Electrophoresis 34, 662-667. Zhao, C., Yang, C., 2011. An exact solution for electroosmosis of non-Newtonian fluids in microchannels, J. Non-Newtonian Fluid Mech. 166, 1076-1079. Zhao, C., Yang, C., 2009. Exact solution for electroosmotic flow of viscoelastic fluids in rectangular micro-channels, Appl. Math. Comput. 211, 502-509. Zhao, C., Zholkovskij, E., Masliyah, J., Yang, C., 2008. Analysis of electroosmotic flow of power-law fluids in a slit microchannel, J. Colloid Interface Sci. 326, 503-510. Zheng, Z., Hansford, D. J., Conlisk, A. T., 2003. Effect of multivalent ions on electroosmotic flow in micro and nanochannels, Electrophoresis 24, 3006-3017.

*Highlights

Highlights

Flow and heat transfer effect of power law fluid is studied in hydrophobic microchannel. Joule heating and thermal radiation effect is evaluated using the velocity slip boundary conditions. Both analytical and numerical results are validated by Debye-Huckel linearizing principle. Flow enhancement factors are distinguished in shear thinning and shear thickening fluid. Increase in thermal radiation minimizes heat transfer rate close to the hydrophobic wall.