Numerical Investigation of Thermally Developing and Fully Developed Electro-Osmotic Flow in Channels with Rounded Corners

: Electro-osmotic ﬂow, that is, the motion of a polar ﬂuid in microducts induced by an external electric ﬁeld, is one micro-effect which allows ﬂuid circulation without the use of mechanical pumping. This is of interest in the thermal management of electronic devices, as microchannels with cross sections of almost arbitrary shape can easily be integrated on the chips. It is therefore important to assess how the geometry of the channel inﬂuences the heat transfer performance. In this paper, the thermal entry region and the fully developed electro-osmotic ﬂow in a microchannel of rectangular cross section with smoothed corners is investigated for uniform wall temperature. For the fully developed region, correlations for the Poiseuille and Nusselt numbers considering the aspect ratio and nondimensional smoothing radius are given, which can be used for practical design purposes. For thermally developing ﬂow, it is highlighted how smoothing the corners increases the value of the local Nusselt number, with increases up to 18% over sharp corners, but that it also shortens the thermal entry length. It is also found that Joule heating in the ﬂuid may cause a reversal of the heat ﬂux, and that the thermal entry length has a linear dependence on the Reynolds number and the hydraulic diameter and on the logarithm of the nondimensional Joule heating.


Introduction
Fundamental research on microchannels has given rise to a vast body of studies ever since the seminal work by Tuckerman and Pease almost four decades ago [1]. Knowledge acquired over the years allowed microchannels to establish themselves as the building blocks of so-called micro-flow devices (MFDs), as can be realized by comparing older reviews on the subject [2,3], to more recent ones [4][5][6], a fact also recognized by the Tuckerman and Pease themselves in a retrospective article published in 2012 [7].
Currently, there is still an amount of research addressing fundamental subjects [8][9][10][11][12][13][14][15][16], but MFDs have now found applications in several fields [17], e.g., micro heat exchangers (MHXs) are employed in air conditioning systems [18] and heat pumping equipment [19]. The use of microjets and two-phase flow as cooling means has been investigated and applied too [20][21][22]. Microchannel heat sinks, a subset of MHXs, have great applicative potential where removal of high heat fluxes is in demand. This makes them particularly interesting for the thermal management of electronic devices, which are steadily growing in compactness and, as a consequence, in power density. Pressure drop across such devices, however, may be significant, especially when liquids are employed as coolants, with the reduction in hydraulic diameter of the ducts quickly leading to viscous heating of the fluid circulated [23] which reverses the direction of the heat flow.
to obtain them are reported in Appendix A. The equations are written for a Cartesian coordinate system, which has x as the axial coordinate (flow direction), while y and z are the coordinates of the cross section (orthogonal to flow direction), as illustrated in Figure 1.

Electric Potential Field
The nondimensional formulation of the equation governing the electric potential distribution is shown below, where k D is the Debye-Hückel parameter (m −1 ): The value of the potential at the Stern plane ζ is used as a boundary condition over the channel wall: As T-boundary conditions are applied to the temperature field, the nondimensional electric potential Ψ, the Debye-Hückel parameter k D , and the nondimensional Stern potential Ψ 0 , Equations (2) and (3) are defined with respect to the imposed wall temperature T w , which is chosen as the characteristic temperature. Accordingly, variations of T w and, subsequently, of the inlet temperature of the fluid, T i , when the temperature drop T w − T i is fixed, determine the range of values taken by the Debye-Hückel parameter in this work.

Velocity Field
A steady, hydrodynamically developed laminar flow of a Newtonian, incompressible fluid with uniform thermophysical properties through a microchannel of uniform cross section with rigid, non-porous wall is considered.
Starting from the governing equation and introducing suitable normalising quantities, as reported in Appendix A.2, the component along x axis of Equation (A5) can be written in nondimensional form: The no-slip boundary condition is imposed through the solid wall, U| ∂Ω = 0. The only velocity component in the case of fully-developed flow is along the X direction; therefore, equations for the other two velocity components are not shown. Once the velocity field U(Y, Z) is computed, the corresponding Poiseuille number can be estimated, details are given in Appendix A.2.

Temperature Field
The nondimensional temperature distribution Θ(X, Y, Z) is obtained from the solution of the nondimensional energy equation below, where Pr is the Prandtl number of the fluid, while U b and Re are the nondimensional bulk velocity and the Reynolds number, The temperature T w is imposed (T-boundary conditions) through the channel wall, and a uniform temperature profile is chosen over the inlet section X 0 = 0 of the channel: The choice of the boundary conditions imposed to the temperature field through the inlet section will be discussed further in Section 3. Once the temperature field Θ(X, Y, Z) is computed, the average Nusselt number along channel axis can be estimated, as explained in Appendix A.3.

Numerical Modeling
The numerical solution of the governing equations-Equations (1), (4) and (5)-was obtained using the open source software GNU Octave. Additional libraries, provided by Octave Forge, were used: the msh package, which relies on the open source software gmsh, allows to generate triangular and tetrahedral unstructured meshes for Finite Element Method (FEM) or Finite Volume Method (FVM) solvers, and the bim package, which implements both the Finite Volume Scharfetter-Gummel (FVSG) method and the Finite Element Method (FEM), allows to solve Diffusion Advection Reaction (DAV) partial differential equations. The functions bim2a_laplacian, bim2a_reaction, and bim2a_rhs belonging to bim library allow to assemble the finite element discretised operators for solving partial differential equations of the form on a two-dimensional (2D) unstructured, triangular mesh. The 2D computational domain, consisting of the channel cross section geometry, is discretized through gmsh software, using a triangular mesh with increased accuracy imposed close to the cross section perimeter, where higher gradients for both the electric potential and the velocity fields are expected. Equation (1) is numerically solved first over the 2D mesh of the channel cross section to compute the electric potential Ψ(Y, Z). In order to reduce Equation (1) to the same form of Equation (8), linearization of the reaction term, (k D D h ) 2 sinh Ψ, is carried out using a first-order Taylor polynomial and the numerical solution of Equation (1) is iterated until convergence of the electric potential field is obtained. Then, the numerical solution of Equation (4), which requires the knowledge of Ψ, is computed and the developed velocity field U(Y, Z) is solved over the channel cross section.
As the energy equation, Equation (5), has the form of a parabolic partial differential equation, the three-dimensional physical problem can be reduced to a 2D mathematical problem. Implementing the θ−weighted scheme for marching over channel axis, Equation (5) is discretized as where n denotes the integration step along the channel axis. At each step n, Equation (9), which has the same form as Equation (8), is numerically solved over the 2D mesh of the channel cross section through bim library functions.

Problem Setup
The electro-osmotic flow over a channel of length L with constant cross section is investigated. The geometrical configuration, defined by a rectangular cross section with rounded corners as shown in Figure 1, is characterized by the aspect ratio β = b/a ∈ [0, 1] and the nondimensional smoothing radius γ = 2 r/b ∈ [0, 1]. The corresponding hydraulic diameter can be calculated as a function of the length of the longer side, a, and of the nondimensional parameters β and γ: Taking advantage of the problem symmetry, only one-quarter of the whole channel section is investigated, in order to reduce the computational cost of simulations, and symmetry conditions are imposed through the lines Y = 0 and Z = 0, ∇Ψ ·n = 0, ∇U ·n = 0, ∇Θ ·n = 0 (11) Dirichlet boundary conditions are imposed along the heated perimeter Γ, which correspond to imposed Stern potential, no-slip condition, and prescribed wall temperature respectively. In addition, an inlet temperature profile must be imposed in order to solve the energy equation, Equation (5). As pointed out in [48,49], when a heat source is present, an adiabatic preparation of the fluid (at least over the hydrodynamically developing region) should be provided in order to ensure physical consistency of the boundary conditions at X = 0. Yet, under certain flow conditions the increase in the bulk temperature ∆T 0 due to Joule heating over the adiabatic section can be neglected if compared to the prescribed temperature difference between wall and inlet. This makes the uniform inlet temperature profile a still consistent approximation. In fact, writing the macroscopic energy balance over the adiabatic section, and assuming that L 0 equals the hydrodynamically developing length, L hy = 0.05 Re D h , the increment ∆T 0 reduces to If water at ambient temperature T amb = 20 • C is considered, i.e., Pr ≈ 7, ∆Θ 0 < 7.2% is ensured for Q ≤ 10, which corresponds to the maximum value of the heat source investigated in this paper. The Graetz problem for an electro-osmotic flow, which consists of the numerical solution of Equations (1), (4) and (5) with the prescribed boundary conditions, is defined by the following parameters: • aspect ratio β of the channel cross section; • nondimensional smoothing radius γ; • nondimensional Stern potential Ψ 0 ; • nondimensional parameter k D D h , defined by Equation (2), which affects the electric potential field; • nondimensional electric fieldẼ = E D h /ζ, which causes fluid motion; • nondimensional heat source Q due to Joule heating; and • reference Peclet number, where u 0 is the reference velocity, calculated through Equation (A6).
Note that the reference Peclet number is not fully representative of the flow characteristics. The effective value of the Peclet number, is verified a posteriori to vary in the range 2 × 10 1 < Pe < 4 × 10 2 ; this makes the assumption of negligible axial conduction a consistent approximation.

Grid Independence Analysis
A grid independence analysis (GDA) was also conducted before running several parametric computations. Excluding simulations involving model verification and comparison with literature results, the values of the Debye-Hückel parameter and the nondimensional Stern potential were set to k D D h ≤ 9.77 and Ψ 0 ≤ 3.89 respectively for most of the parametric computations. As mentioned at the end of Section 2.1, these are determined by the choice of the wall and fluid inlet temperatures; the remaining quantities which compose the Debye-Hückel parameter and the Stern potential have the values of deionized, ultra-filtered water, and are often adopted in the literature [28,[30][31][32][33]44,45,50]. The highest values were considered for the GDA, while the channel cross section aspect ratio and the smoothing radius were set to β = 3/4 and γ = 1/2.
Equations (3), (4) and (A15) were numerically solved for increasing mesh accuracy, which is defined by the imposed grid spacing parameter l c . Both the Poiseuille and Nusselt numbers of the fully developed flow were computed through Equations (A7) and (A14). Results are listed in Table 1, where the corresponding numbers of mesh nodes and triangular elements are also reported, and are plotted in Figure 2. As a discrepancy of less than 0.35 was observed for Po reducing l c from 5 × 10 −3 to 2.5 × 10 −3 , a value of l c = 5 × 10 −3 was chosen for all the numerical computations. Figure 2 also shows the triangular mesh, generated through the gmsh software, for l c = 4 × 10 −2 : it can be noticed that mesh refinement occurs along the heated perimeter, where higher gradients of the electric potential Ψ(Y, Z), velocity U(Y, Z) and developed temperature Θ ∞ (Y, Z) fields are expected.
When the Graetz problem is solved, the discretization step ∆X over the channel axis direction is imposed in order to ensure at least 200 iterations before thermally developed condition is reached, where Nu(X) and Nu ∞ are computed via Equations (A12) and (A14).

Model Validation
Comparison with numerical results in the literature was first conducted in order to check the numerical solution of both the Poisson-Boltzmann equation and the momentum equation, Equations (1) and (4), which were both numerically solved in [45] by means of COMSOL Multiphysics ® . Considering the same test case as in [45], k D D h = 9.85 and Ψ 0 = 7.92, several geometrical configurations were simulated and the corresponding Poiseuille number computed. As high gradients ∇Ψ ·n, due to the high Stern potential, were expected close to the channel wall, the grid spacing parameter was set to l c = 2.5 × 10 −3 , leading to meshes with more than 10 5 nodes. A good agreement with the numerical results of [45] can be observed in Table 2, with differences below 1% on the computed Poiseuille number always obtained. The electric potential field was also qualitatively verified with the analytical solution, which is available in the limiting case of a cylindrical duct for Stern potential approaching zero. Solving the axisymmetric Poisson-Boltzmann equation under Debye-Hückel linearization, i.e., lim Ψ→0 (sinh Ψ) → Ψ, gives [51], where R ∈ [0, 1/2] is the nondimensional radial coordinate; I 0 and K 0 are the first and second kind of modified Bessel functions; and C 1 and C 2 are constants defined by the boundary conditions imposed. As a fully developed flow subject to H-boundary condition was investigated in [45], comparison was not possible. However, the developed temperature field under T-boundary condition was quantitatively verified with the analytical solution available for a cylindrical duct. Solving the axisymmetric energy equation for a fully developed flow under T-boundary conditions gives the analytical temperature profile, where R is the nondimensional radial coordinate. Both normalized, numerical electric potential profiles for different values of the imposed Stern potential Ψ 0 and normalized, analytical electric potential, valid in the limit of Ψ 0 approaching zero, are plotted in Figure 3 over the cylindrical cross section: note that the numerical profile approaches the analytical one when low values of Ψ 0 are imposed. The numerical velocity profiles for different values of Ψ 0 are shown in Figure 3: as mentioned, the higher the value of the Stern potential, the higher the velocity gradient close to the channel wall, resulting in a flow pattern similar to plug flow. The developed temperature profile, which is also shown in Figure 3 and does not depend on the velocity profile, was computed both via the numerical solution of Equation (A15) and as the asymptotic behavior for X → ∞ from the solution of Equation (5). Thus, the numerical temperature profile was quantitatively verified with the analytical solution, Equation (19), computing the mean squared error (MSE) at mesh nodes, (20) showing perfect agreement. The numerical model was also validated with experimental results from the literature. First, the test case of [52] is replicated. Thus, the bulk flow velocity of an EOF is numerically computed and compared with experimental data of [52], where the Authors used nanoparticle image velocimetry technique to measure the mean EOF velocity over 4.9-24.7 µm thick channels. The physical properties of a dilute aqueous solution of sodium tetraborate (borax) were chosen according to the method in [52]. The ζ potential required for the boundary condition on the potential field is derived from the experimental values of the mobility factor µ eo via the Helmholtz-Smoluchowski formula: In Figure 4a, numerical results are compared with the experimental data in [52], showing a good agreement for all the concentrations investigated. A small discrepancy is observed at higher values of the electric field, since a linear relationship between bulk velocity and electric field is predicted by the numerical model. However, a linear relationship was also reported according to [26].    In order to validate the temperature field solution, the experimental setup in [26] was also simulated. In [26], the Authors manufactured and tested a micro-heat sink with the cooling fluid serving both as a heat exchanger and an integrated pump. Twenty parallel 0.3 × 0.1 × 30 mm 3 microchannels were etched on a silicon substrate and covered by a low thermal conductivity PDMS plate. A heater was attached to the bottom. The borax concentration was set to 0.4 mol/m 3 , different values of the voltage driving the EOF were imposed and heating powers up to 4 W were investigated in [26]. Equations (1), (4) and (5) were therefore numerically solved and the resulting bulk temperature at both the middle section (x = 1.5 cm) and the outlet section (x = 3 cm) compared with the experimental data of [26], for different values of the imposed heating power. The required wall potential was again estimated through Equation (21). Adiabatic conditions were imposed at the channel wall facing the PDMS cover, whilst H1 thermal boundary conditions were applied to the heated perimeter. The required heat flux per unit length was derived from the data supplied in [26]. The electric field was set to 1.33 × 10 4 V/m. As shown by Figure 4b, numerical temperatures are in good agreement with experimental measurements over the whole heating power range, with a perfect matching observed when the heater is switched off (temperature rise along the channel due solely to Joule heating). The resulting temperature profile through the outlet section is also provided by Figure 4b.
The electric potential field Ψ, the developed velocity field U and the developed temperature field Θ ∞ resulting from a single computation are shown in Figure 5. As pointed out for the cylindrical duct case, significant gradients of both the electric potential and the velocity can be observed close to the cross section wall. The magnitude of both U and Θ linearly depends on the imposed values ofẼ and Q. The local Nusselt number was also numerically estimated as over the heated perimeter of the cross section, Γ. As demonstrated by Figure 6, the highest temperature gradient occurs at the midpoint of the longer side of the rectangular section, whereas smoother slopes are observed close to the rounded corners. Averaging the local Nusselt number along the heated perimeter leads to the mean Nusselt number Nu = 6.662, which can be computed through Equation (A14). The results from 110 computations are shown in Figure 7, where the computed Poiseuille and Nusselt numbers are plotted as a function of the channel section geometry, defined by β and γ. The highest values for both Po and Nu were found at β = 1/10 and γ = 1, whilst the minimum Po, corresponding to highest efficiency of the micropump, was found for β = 1 and γ = 0. Following the approach in [45], the trends of Po and Nu can be approximated, for a given aspect ratio, via a polynomial correlation as a function of the smoothing radius. However, a single, 3rd-order polynomial correlation taking into account the influence of both β and γ is provided here for Po and Nu: Po = a 00 + 3 ∑ n=1 a n0 β n + 3 ∑ n=1 a 0n γ n + a 11 a 11 β γ + a 12 β γ 2 + a 21 β 2 γ (23) The corresponding fitting coefficients are listed in Table 3. The Root Mean Square Error, also listed in Table 3, demonstrates that an accurate estimate of both Po and Nu is ensured.

Graetz Problem Solution
Finally, the Graetz problem, defined by Equations (1), (4) and (5), was solved, in order to determine the thermal entry length under different operating conditions. The aspect ratio was fixed to β = 3/4 and four different, nondimensional radii of curvature, namely, γ = 0, 1/4, 1/2, 1 were investigated. The Debye-Hückel parameter and the nondimensional Stern potential were set to k D D h = 8.45 and Ψ 0 = 2.92 respectively. The nondimensional electric field and the reference Peclet number wereẼ = 3 × 10 2 and Pe = 1.15. The effect of changing the heat source due to Joule heating, Q, was investigated: several values ranging between 10 −3 and 10 were simulated and the influence on the thermal entry length assessed. Note that a variation of the heat source Q at fixed k D D h , Ψ 0 ,Ẽ and Pe, which are the other model input parameters, can be physically obtained varying the reference temperature difference T w − T i and the electric conductivity σ (so that Pe is unaffected).
The resulting temperature profile over the channel sections corresponding to Y = 0, Z = 0 and X ∈ [0, 2] is shown in Figure 8 for a computation at very high heat generation (Q = 8.04) and low electric field (Ẽ = 2.25 × 10 −2 ). Even if the fully developed temperature profile is reached shortly after the channel entrance because of the high value of Q, important considerations can be made: • close to the entry section, the contribution of Joule heating is still negligible and the heat transfer process is mainly driven by conduction and convection, which is the most profitable operating condition for a microchannel heat sink; • looking at the mean Nusselt number along the channel axis in Figure 8, a vertical asymptote at X 0.535 is observed, corresponding to the change in sign of the fluid bulk temperature; • for X > 0.535, the heat generation due to Joule heating drives the heat transfer process and the microchannel heat sink is no longer able to dissipate heat from the walls; and • the thermal entry length L th is reached at X = 1.37.
Assessing the existence of a change in the direction of heat transfer is obviously important in practical applications, but knowledge of the extension of the thermal entry length is also useful. In fact, if L th occupies a significant portion of the channel, correlations for the Nusselt number that assume fully developed flow cannot be recommended. Therefore, the rest of the section is devoted to the investigation of L th . Results from parametric computations are shown in Figure 9, where the thermal entry length is plotted as a function of the nondimensional source term Q and the smoothing radius γ. It can be pointed out that L th decreases for increasing Joule heating following a logarithmic behavior, when the other model input parameters (namely, k D D h , Ψ 0 ,Ẽ, Pe) are fixed: The mean Nusselt number along channel axis is also plotted in Figure 9 for different γ and fixed Q = 1, showing that • smoothing the corners shortens the thermal entry length, and, as a consequence, the efficacious channel length (i.e., useful to remove heat from the walls) at high Joule heating, that is for high values of Q; • heat transfer performance decreases when sharp corners are considered, since lower temperature gradients are observed close to them, meaning that lower Nusselt numbers over the entry region are reached. As an example, in Figure 9 for x/D h = 0.9 Nu goes from 3.95 for γ = 1 (maximum smoothing) to 3.38 for γ = 0 (sharp corners); • small changes in Nu(x) correspond to significant variations of L th at different values of the smoothing radius γ.
In order to extend the study to a wider range of practical applications, the analysis is expanded to different hydraulic diameters of the channel cross section. Thus, the nondimensional input parameters of the model are scaled with the hydraulic diameter according to their definitions: Starting from the set of input parameters investigated and fixing the cross section geometry to β = 3/4 and γ = 1/2, computations are run at increasing hydraulic diameters. Thus, beginning at D 0 = 3 µm, the hydraulic diameter is progressively increased up to 15 µm, which corresponds to an increase of the nondimensional Debye-Hückel parameter k D D h from 8.45 to 42.25. Three different values of the nondimensional heat source, Q = 10 −1 , 10 −2 , and 10 −3 , are chosen. Both the thermal entry length and the Reynolds number, the latter resulting from velocity field solution and depending on the imposed hydraulic diameter as well, are computed. Results are plotted in Figure 10. The nondimensional thermal entry length linearly depends on the hydraulic diameter, with values up to L th = 40 D h observed when D h /D 0 approaches 5. This trend can be easily explained looking at Equation (5): the advection term is proportional to Re and Pr, both acting as a scaling parameter for the temperature field evolution along the channel axis, as long as Q is kept constant. On the other hand, this means that L th /D h weakly depends on the shape of the normalized velocity profile U/U b , which tends to plug flow for increasing D h but does not modify the linear correlation between L th /D h and Re. As discussed earlier for Figure 9, the lower the magnitude of heat source term Q, the longer the nondimensional thermal entry region L th /D h . The dependence of the Reynolds number on D h /D 0 is also reported in Figure 10, and is clearly linear. As the magnitude of the dimensional bulk velocity u b mainly depends on the imposed electric field and the Stern potential but weakly on the hydraulic diameter, as long as k D D h 1, the product of the Prandtl number, which depends on fluid properties, times the Reynolds number Re = ρ u b D h /µ is simply proportional to the hydraulic diameter. As a consequence, it can be stated that the nondimensional length of the thermal entry region L th /D h linearly depends on the normalized hydraulic diameter D h /D 0 too. The influence of wall temperature, which was chosen as the characteristic fluid temperature, was also investigated. Therefore, T w was varied between 20 and 100 • C and three different values were again considered for the nondimensional heating source, Q = 10 −1 , 10 −2 , and 10 −3 . The channel section geometry, the wall potential and the nondimensional electric field were set to β = 3/4, γ = 1/2, Ψ 0 = 2.92, andẼ = 3 ×   Figure 11a shows the dependence of the nondimensional thermal entry length on the imposed wall temperature. In accordance with Re Pr, L th /D h increases slightly nonlinearly with the imposed temperature T w : this is readily explained if one considers that the kinematic viscosity in the expressions of the Prandtl and Reynolds numbers cancel out.

Conclusions
Although not the main focus of the paper, the fully developed EOF under uniform wall temperature was first investigated for channels with sharp and smoothed corners and rectangular cross sections, after validation with two different experimental works and verification with analytical and numerical cases. The numerical findings highlighted that smoothing increases both the Poiseuille and Nusselt numbers, i.e., friction and heat transfer, respectively. Two correlations which capture simultaneously the influence of the nondimensional aspect ratio β and radius of curvature γ on Po and Nu were obtained. These are recommended for design purposes, when the flow can be considered as fully developed for most of the channel length. It was also highlighted that neither the magnitude ofẼ (nondimensional electric field) nor Q (nondimensional Joule heating) influence the form of the velocity and temperature profiles. Similar results have not been previously reported in the literature, to the best of the Authors' knowledge.
The Graetz problem under T boundary condition was then investigated, using input parameters representative of practical applications, such as micropumps and micro heat sinks. This subject is even less explored, owing to the experimental difficulties in obtaining temperature distribution in the fluid and at the walls of the channels. Numerical studies also appear to be limited to fully-developed cases and very few are concerned with channels with smoothed corners. The investigation highlighted that high values of Q cause a reversal in the direction of the heat flux, which renders the channel unfit for its purpose. An estimate of the usable channel length under such circumstances was also provided via the solution of the temperature field over the thermal entry region under different operating conditions.
It was then demonstrated that smoothing the corners shortens L th , which may significantly reduce the usable length of the channel as well, but the Nusselt number in the entry region increases with the radius of curvature. It can therefore be concluded that, when Q is sufficiently low that no heat flux reversal occurs in the channel, smoothing the corners increases the Nusselt number (for the cases investigated an improvement up to around 18% was achieved).
Finally, the thermal entry length was analyzed for different values of the model input parameters. The magnitude of the heat source term and of the variation in hydraulic diameter were investigated, with a wide range of operating conditions covered, and it was found that the nondimensional thermal entry length decreases linearly with the logarithm of Q, while scaling linearly to D h . The dependence of L th on the Reynolds and Prandtl numbers was also studied. If the Prandtl number is constant, the thermal entry length is proportional to the Reynolds number; when both vary because of changes in the reference temperature, L th slightly deviates from proportionality to Re Pr because of the dependence of thermal diffusivity on temperature, whereas the kinematic viscosity cancels out.
The results presented are crucial for the design of microchannel heat sinks working with an electro-osmotic flow; moreover, knowledge of the Poiseuille number and Nusselt number over the channel length under different operating conditions is required in order to determine, e.g., the best microchannel configuration via a combined first-and secondlaw-based optimization under typical PEC constraints [40,43,53,54]. In the next future the Authors plan to investigate different geometries, defined by the aspect ratio β and the smoothing radius γ, and the optimal configuration determined for different PECs and different operating conditions (hydraulic diameter D h , electric double layer thickness k D −1 , electric potential E x , source heat power Q, and Stern potential ζ). Moreover, the effect of axial conduction through the fluid can eventually be accounted for in the numerical model, in order to properly investigate electro-osmotic flows characterized by lower Reynolds numbers and, thus, lower Peclet numbers. Note that the developed temperature field Θ ∞ (Y, Z) can be either estimated as Θ(X, Y, Z)| X→∞ or via the solution over the 2D channel cross section of the following partial differential equation, which Equation (5) reduces to for X → ∞.