Nusselt Numbers for Poiseuille Flow Over Isoflux Parallel Ridges for Arbitrary Meniscus Curvature

An mm-scale sessile droplet on a structured surface characterized by small periodic length scales compared to the capillary length is, on the macroscopic scale, in the shape of a spherical cap [1]. It may be stable in the wetted (Wenzel) state, where the liquid penetrates into the structures, or in the unwetted (Cassie) state, where solid–liquid contact is limited to the tips of the structures [2]. Heat transfer between a textured surface and a droplet is relevant to, e.g., enabling Wenzel to Cassie state transition [3], prevention of ice formation [4], and long-term inhibition of transition from dropwise to filmwise condensation [5].

A liquid flowing along a diabatic, structured surface may also be stable in the Cassie state. We assume this here and the necessary criteria are provided by Lam et al. [6]. Solid–liquid interfaces are subjected to the no-slip boundary condition, but lubrication is achieved because of the low shear stress at liquid–gas (and/or vapor) interfaces (menisci). The insulating gas causes an increased convective thermal resistance to the heat exchange between the liquid and the underlying substrate. Thus, the magnitude of the temperature difference between a solid–liquid interface and the bulk liquid is increased on account of the reduction in surface area for efficient heat transfer. The foregoing effects are referred to as apparent hydrodynamic slip, a (generally) desired effect, and apparent thermal slip (or temperature jump), a (generally) undesired effect, as discussed by, e.g., Rothstein [7] and Enright et al. [8], respectively. There are applications where insulating a flowing liquid from the ambient and thus apparent thermal slip are desirable effects [9].

Another potential application of apparent slip is the enhancement of (single phase) direct liquid cooling of microelectronics, where a liquid coolant is pumped through microchannels in the die of an operating semiconductor device [10]. Indeed, their small form factors, e.g., about 1.4 cm × 1.4 cm footprint × several hundred microns thick in the case of a microprocessor in a desktop computer [11], imply high flow and thus caloric resistances. This can be reduced by apparent hydrodynamic slip. Although this is accompanied by an increase in the aforementioned convective component of thermal resistance, it has been predicted that in the case of a liquid metal coolant, there is a reduction in total thermal resistance [6]. The same applies to indirect liquid cooling, where the liquid is pumped through the microchannels in a cold plate attached to the package of a microelectronics device. Not surprisingly, there has been an expanding body of literature providing Nusselt numbers (Nu) for convection in the presence of apparent slip, albeit restricted to the parallel-plate microchannel flow configuration thus far.

The Reynolds and Péclet numbers characterizing the flow and heat transfer, respectively, near a composite interface (solid–liquid interface and adjacent meniscus) may be low enough for transport to be diffusive. Simultaneously, the microchannel height to structure (ridge parallel or transverse to flow, square pillar, etc.) pitch may be high enough for advection to be relevant elsewhere in the microchannel as discussed by Enright et al. [8]. Then, assuming a flat, shear free and adiabatic meniscus, the apparent hydrodynamic and thermal slip lengths nondimensionalized by the structure pitch follow from existing solutions to the Stokes and Laplace equations, respectively, in semi-infinite domains and are a function of the solid fraction, and type, of the structures [8]. These parameters, in turn, manifest themselves in Robin boundary conditions imposed on the mean streamwise momentum and mean thermal energy equations, i.e., those where the dependent variable is the velocity and temperature averaged over the structure pitch, respectively [8]. In the case of hydrodynamically and thermally fully developed Poiseuille flow over isoflux structures of any type subjected to the foregoing constraints, Enright et al. [8] provide an expression for the Nusselt number. It is a function of the dimensionless apparent hydrodynamic and thermal slip lengths on each side of the channel and the ratio of heat fluxes imposed on the composite interfaces on opposite sides of the channel. Kane and Hodes [12] generalized this result to a combined Poiseuille and Couette flow.

Additional physical effects have been incorporated into more rigorous expressions for the apparent thermal slip length than those developed by Enright et al. [8]. Ng and Wang [13] relax the assumption of an adiabatic meniscus to capture conduction through the gas phase in the case of isothermal, ridge-type structures. The dimensionless apparent thermal slip length is then a function of the gas-to-liquid-phase thermal conductivity ratio and the ridge depth to ridge pitch ratio in addition to solid fraction. Hodes et al. [14] captured the effects of evaporation and condensation along the meniscus on account of its nonuniform temperature such that the dimensionless apparent thermal slip length is a function of a dimensionless heat transfer coefficient along the meniscus and solid fraction. Finally, Lam et al. [15] accounted for small meniscus curvature such that the thermal slip length is a function of both the solid fraction and the protrusion angle of the meniscus, i.e., that between the horizontal and the tangent to the meniscus emanating from the triple contact line, with positive and negative values of it corresponding to protrusion of the liquid phase into the channel and a cavity, respectively. More rigorous expressions for apparent hydrodynamic slip lengths than those available in Enright et al. [8] exist as well. For example, Sbragaglia and Prosperetti [16] captured the effects of small meniscus curvature on the apparent hydrodynamic slip length for ridges oriented parallel to the flow.

Mathematically, the Robin boundary conditions used by Enright et al. [8] are of the same forms as the simplest ones of those capturing the effects of molecular slip on flow and heat transfer when gases flow in the slip regime [17], i.e., when the Knudsen number (Kn = mean free path of gas to hydraulic diameter) is between about 0.001 and 0.1. In the case of Poiseuille flow in this regime, Colin [18] provides a critical review of relevant Nusselt number results. Subject to the foregoing constraints imposed by Enright et al. [8], many of them, e.g., those by Sparrow and Lin [19] for hydrodynamically and thermally fully developed flow between isothermal parallel plates, are directly applicable to flows exhibiting apparent rather than molecular slip. Conversely, the solution by Hooman and Ejlali [20] uses a perturbation method to capture the effects of (small) changes to viscosity and thermal conductivity due to their temperature dependence on Nusselt number. It would not be as accurate when utilized to find a Nusselt number for a flow exhibiting apparent slip while using the apparent slip length expressions documented by Enright et al. [8]. Indeed, the latter do not account for variations of these properties insofar as computing the velocity and temperature fields near the composite interface, an assumption that could be relaxed.

Caution must be exercised when selecting apparent slip length expressions for use in Nusselt number expressions that are developed by employing effective Robin boundary conditions, for example by Enright et al. [8]. This is because implicit in the analysis of Enright et al. is the assumption that the multidimensionality of the velocity and temperature fields is confined to an inner region close to the structuring, which requires the channel height to be large in comparison to the structure pitch. Analytical formulae for the apparent slip length of an unbounded shear-driven flow over the structures, such as Philip [21], can then be used in their Nusselt number expressions, representing the behavior in this inner region. But since Enright et al. did not account for curvature, only slip length formulae without curvature can be used with confidence. Using slip length formulae that account for curvature in an unbounded flow, such as those of Sbragaglia and Prosperetti [16] and Crowdy [22], would only account for changes to the velocity field, but not the change in cross-sectional area, which was shown to be an important contribution to the overall flow enhancement by Sbragaglia and Prosperetti [16] and Kirk et al. [1]. In particular, this contribution becomes more significant as the height of the channel decreases. The use of apparent thermal slip lengths that account for curvature, such as Lam et al. [15], may also introduce unforeseen complications. An alternative and unambiguous approach, which also relaxes the large channel height to structure pitch ratio constraint, is to solve for the velocity and temperature fields in the entirety of a domain and compute the Nusselt number directly. Subsequently, one may calculate the corresponding apparent thermal slip length from a Nusselt number expression based on Robin boundary conditions. However, since the Nusselt number is the engineering parameter of interest, this is superfluous, especially when the curvature of the meniscus is significant and the notions of apparent slip lengths become physically unrealistic for the bounded flows of interest here.

The first study on convection in the presence of apparent slip where the convective transfer equations were directly solved was a two-dimensional, numerical one for laminar, periodically hydrodynamically and thermally fully developed flow between parallel plates symmetrically textured with isothermal ridges oriented transverse to the flow by Maynes et al. [23]. The liquid and gas phases were water and air, respectively, and velocity, temperature, shear stress, and heat flux were matched along a flat meniscus. Inertial terms were retained in the streamwise and transverse momentum equations and the axial conduction term was retained in the thermal energy equation. No-slip and isothermal boundary conditions were imposed along the periphery of the wall of the cavity in which the air recirculates. The computed Nusselt numbers are dependent upon the solid fraction of the ridges, distance from the leading edge of a ridge, aspect ratio of the cavity, channel height to ridge pitch ratio and, due to inertial effects, Reynolds number of the water. When averaged over the composite interface, they do not exceed the value for a smooth channel (Nu = 7.54). However, local Nusselt numbers over the ridge always exceed this value.

A subsequent investigation by Maynes et al. [24] analytically considered the foregoing problem, but for isoflux rather than isothermal transverse ridges and the more general case of hydrodynamically developed and thermally developing flow, i.e., a Graetz-type problem. The meniscus is flat and, unlike in the foregoing study, shear free and adiabatic. A (unidirectional) Navier-slip velocity profile based on the apparent hydrodynamic slip length for creeping flow over transverse ridges in a circular tube in the limit of large tube radius to ridge pitch ratio developed by Lauga and Stone [25] multiplied by a term given by Woolford et al. [26] to account for inertial effects at sufficiently high Reynolds numbers was utilized in the thermal energy equation (Teo and Khoo [27] show that the Lauga and Stone [25] expression is valid for the parallel-plate geometry in the large channel height to ridge pitch ratio limit). The axial conduction term was neglected in the thermal energy equation, which was first solved for a constant heat flux boundary condition along the entire composite interface. Then, by embedding the resulting temperature field in Duhamel's integral, it was determined for arbitrary solid fraction (periodic heating). When averaged over the composite interface, Nusselt numbers in the periodically hydrodynamically and thermally fully developed flow limit do not exceed the value for a smooth channel (Nu = 8.235). Despite neglecting axial conduction, the Nusselt number is dependent on the Péclet number because it affects the thickness of the thermal boundary layers that develop along the ridges. (The Prandtl number is that of water and thus the Reynolds number is not an independent variable.) Additionally, Cowley et al. [28] considered this problem in the periodically hydrodynamically and thermally fully developed flow limit using the same Navier slip velocity, but numerically solved the thermal energy equation while retaining the axial conduction term. The mean temperature of the liquid was assumed to vary linearly above the ridge and to be constant above the meniscus. Alternatively, a rigorous (numerical) solution of the thermal energy equation may be obtained by separating the temperature into linear and periodic components using the formulation by Patankar et al. [29] as has been done in the studies on pillar-type structures discussed below. Cowley et al. [28] also performed a small number of simulations for isothermal rather than isoflux ridges in the context of the Graetz–Nusselt problem. No other study has used a brute-force numerical approach to resolve the temperature field in a domain longer than one structure pitch.

Maynes and Crockett [30] analytically determined the Nusselt number as a function of solid fraction and channel height to ridge pitch ratio for laminar, hydrodynamically and thermally fully developed flow in a parallel plate channel textured with isoflux, parallel ridges separated by flat, adiabatic menisci. A Navier-slip velocity profile was used, but the thermal energy equation, which is simpler than in the case of transverse ridges because there is no axial conduction, was directly solved. Kirk et al. [1] semi-analytically considered the same problem, but relaxed the assumption of a Navier-slip velocity profile and also captured the effects of small meniscus curvature by perturbing the domain boundary. The validity of the Navier-slip velocity profile assumption was shown to break down at sufficiently small channel height to ridge pitch ratios, even for a flat meniscus. Karamanis et al. [31] semi-analytically solved the Graetz–Nusselt problem for the case of one plate textured with isothermal, parallel ridges and the other one smooth and adiabatic or textured with either adiabatic or isothermal ridges either aligned or staggered relative to the ridges on the opposite plate. The full velocity and temperature profiles were resolved, but the meniscus was assumed to be flat, shear free, and adiabatic and axial conduction was neglected. A subsequent study by Karamanis et al. [32] considered the effects of axial conduction, i.e., the extended Graetz–Nusselt problem. Either both plates were textured with aligned, isothermal, parallel ridges or one plate was textured and the other one was smooth and adiabatic. Viscous dissipation and (uniform) volumetric heating were accounted for in this study.

Enright et al. [8] and, subsequently, Cowley et al. [33] numerically computed Nusselt numbers in the periodically hydrodynamically and thermally fully developed flow limit for a parallel-plate channel symmetrically textured with isoflux, square pillars. The meniscus was assumed to be flat, shear free, and adiabatic. However, the three-dimensional velocity and temperature fields were resolved by separating the pressure and temperature fields into linear and periodic components [29]. The Cowley et al. [33] results quantitatively delineate the ranges of height to structure pitch ratio, Reynolds number, and Péclet number for which the analytical results of Enright et al. for Poiseuille and Nusselt number are accurate. Heat transfer between ridges and a liquid flowing through a parallel plate channel also causes thermocapillary stress, and the potential to use this as a pumping mechanism and its effect on the apparent hydrodynamic slip length have been studied [34–36].

Except for the study by Kirk et al. [1], the effect of meniscus curvature on Nusselt number has not been directly considered. However, it is important in thermal management applications, where the full range of meniscus protrusion angles, i.e., between $−90 deg$ and $90 deg$⁠, is relevant. For example, the Lam et al. [6] study mentioned previously showed that Galinstan, a liquid metal at room temperature used as a nontoxic alternative to mercury, is better suited than water to enhanced microchannel cooling enabled by apparent slip. The thermal conductivity of Galinstan is 28 times that of water, whereas its volumetric heat capacity $(ρcp)$ is 45% that of water and its viscosity is 2.4 times that of water [37]. Consequently, the degradation in Nusselt number and thus convective resistance on account of using a textured rather than smooth surface may be more than compensated for by decreased flow and thus caloric resistance. Since the advancing contact angles of Galinstan on, e.g., Teflon^® and silicon nitride are $161.2 deg$ and $147.0 deg$⁠, respectively [38] (equivalently, protrusion angles of $−71.2 deg$ and $−57.0 deg,$ respectively), the effects of meniscus curvature on Nusselt number must be considered. It is also beneficial that the surface tension of Galinstan is 7.4 times that of water [38]. We also note that the advancing contact angle of water on Teflon is $110 deg$ (equivalently, a protrusion angle of $−20 deg$⁠). Thus, assuming the meniscus is pinned at the advancing contact angle at the channel inlet and is flat at its outlet, the Young–Laplace equation implies that the pressure difference across microchannel can be 20 times higher when Galinstan rather than water is the coolant in the Cassie state.

More generally, by texturing surfaces with re-entrant structures as previously done to stabilize droplets of low surface tension (wetting) liquids in the Cassie state [39,40], the protrusion angle of a Galinstan or water meniscus may reach $−90 deg$⁠. This is illustrated in Fig. 1, where the cross sections of the structures are tee-shaped. The top (straight) line represents the meniscus when there is no pressure difference across it. The middle (circular) arc represents it when this pressure difference equals that required to reach the advancing contact angle $(θa)$⁠. Then, the triple contact lines depin and advance down the vertical sides of the tops of the tees. Finally, the lower (circular) arc represents the meniscus when the contact angle relative to the horizontal surfaces on the bottom sides of the top of a tee is $90 deg$ (equivalently, a protrusion angle of $−90 deg$⁠). This configuration must be achieved before the triple contact lines depin from these surfaces and advance along them. Then, the surface tension force supporting the meniscus has only an (upward) vertical component; therefore, the pressure difference across it is the maximum the liquid can support. Again assuming a flat meniscus at the outlet of a microchannel, this pressure difference is also the maximum one available to drive liquid through it in the Cassie state. Then, the caloric resistance of enhanced microchannel cooling enabled by pumping liquid through a microchannel in the Cassie state is minimized. Regarding the top, rectangular portion of a tee-shaped structure, we note that its height $(rh)$ and the width it protrudes laterally outward from its stem $(rw)$ should be very small compared to the width of its stem $(sw)$⁠. Then, the robustness of the Cassie state of a re-entrant structure is realized, but insofar as lubrication and convection heat transfer, the structures are essentially parallel ridges. This is important because wetting of the sides of structures has been shown by Salamon et al. [41] to have a deleterious effect on apparent hydrodynamic slip and lateral temperature gradients in the top of a structure increase apparent thermal slip. Finally, it is plausible that the menisci may protrude into the channel rather than the cavities to further increase the pressure difference across the channel. Practically, this might be achieved by maintaining the liquid toward the outlet of the channel below atmospheric pressure, and adjusting the geometry of the outlet to prevent gas from escaping.

Having justified the need for Nusselt number expressions for arbitrary meniscus curvature, the rest of this paper is organized as follows: In Sec. 2, we formulate the Nusselt number problem for laminar, hydrodynamically and thermally fully developed Poiseuille flow through a parallel-plate channel symmetrically textured with isoflux parallel ridges. In Sec. 3 (and Appendix  A), we split the domain into two subdomains, both of which are then mapped onto square subdomains, and discuss the Chebyshev collocation (spectral) method used to solve the problem. In Sec. 4, we present and discuss our results. In particular, we calculate the average Nusselt number for a variety of parameter values, and compare our results with the asymptotic results of Kirk et al. [1]. In Sec. 5, we conclude our work with a summary.

The problem studied in this paper is analogous to that studied by Kirk et al. [1]. For completeness, the mathematical model is also included here (Reprinted with permission from Cambridge University Press).

We consider laminar, hydrodynamically and thermally fully developed Poiseuille flow through a parallel plate channel, with uniformly spaced parallel ridges aligned with the flow direction, as in Fig. 2. The top and bottom surfaces of the channel are structured identically; therefore, the flow and temperature fields are symmetric about the center plane of the channel. Surface tension maintains the liquid in the Cassie state, i.e., the liquid contacts the solid substrate only on the tops of the ridges. The grooves between ridges are filled with vapor (and perhaps also noncondensible gas). The pressure difference between the liquid and gas causes the liquid–gas interfaces (menisci) to protrude into (or out of) the cavities.

where μ is the dynamic viscosity of the liquid.

where n is the unit (outward) normal to the interface. The domain, equations, and boundary conditions for the w (and T) problem are summarized in Fig. 3.

where $Q*$ is the total volume flow rate through the cross section A.

and the adiabatic condition on the meniscus (7) is unchanged. Just as for the w problem, see Fig. 3 for a summary of the problem for T. The w problem can be solved in isolation, and then the solution w appears as a forcing in the Poisson problem for T. Note that this is possible for a meniscus of any shape or curvature.

Note that this is an average taken over the entire pitch (not just the ridge), but since there is no heat flux through the meniscus, the integrand is zero along it. The definition used by Kirk et al. [1] included a factor due to the increase in length of the composite interface; however, it was a higher order effect and did not appear in their small curvature results. Therefore, comparisons using definition (17) are valid.

In this paper, we do not provide results for the velocity problem, as there are other studies that provide these for small meniscus curvature [1,16], or for general meniscus curvature [43,44], from which one can infer the caloric resistance of such channels.

We consider the domain as being comprised of two distinct subdomains, henceforth referred to as subdomains 1 and 2 as in Fig. 4. These subdomains are separated by a vertical line originating from the point of intersection between the solid ridge and the meniscus. We map subdomain 1 (resp. 2) to [−1, 1] × [−1, 1] via transformation A (resp. B) as detailed in Appendix  A. The purpose of this is so that we are able to then solve the problem with a Chebyshev collocation approach, as detailed in Sec. 3.3.

Henceforth, we use the original derivative notation, e.g., $∂/∂x$⁠, with the understanding that this is shorthand for the transformed derivatives in terms of ξ_i and η_i as detailed in Appendix  A. In practice, for every x or y derivative that appears in the following, our algorithm applies the discretized form of the equations given in Appendix  A.

where $Hm(ξ1)=hm(x)$⁠. Note that Eq. (22) along with Eqs. (24)–(31) captures equations (2)–(5) for w, while Eq. (23) along with Eqs. (32)–(39) captures equations (7) and (13)–(16) for T.

In order to more accurately solve the system (22)–(39) with a Chebyshev collocation/domain decomposition method, it is necessary to remove the (first derivative) singularities at the triple contact points from the numerical problems. First, we carry out local analyses to establish the functional forms of the singularities which, once subtracted from the velocity and temperature (at the correct strength), give nonsingular solutions.

where a_n and b_n are constants.

where $λ=π/[2(π−θ)]$ corresponds to n when k = 0.

We use methods similar to those found in Peyret [45], and described in brief below, to remove these singularities from our problem since their form is known analytically. Ultimately, we are left with the discretization and computation of a singularity-free problem, with the singularity strengths computed as part of the solution. Thus, we are guaranteed increased accuracy of the Chebyshev collocation method along with the usual numerical convergence properties of nonsingular elliptic problems. It is important to note that using this treatment enables us to handle cases with small solid fractions and/or large (negative) meniscus protrusion angles where the singularities are strongest.

To generate a numerical solution to the system (22)–(39), we apply a Chebyshev collocation method to discretize the transformed equations and boundary conditions. We briefly describe the essential features of our approach here and refer the reader to standard texts [45–47] for more details.

In the following, we give a descriptive definition of each of the elements that appear in the preceding matrix equations. The constants α₁ and α₂ are the unknown strengths of the singularities to be solved as part of the system. The vectors $ui$ and $vi$ contain the values of the singular functions (42) and (43) at each of the Chebyshev grid points representing the ith subdomain, enumerated as described earlier. The unknown vectors $w̃1=w1−α1u1, w̃2=w2−α1u2, T̃1=T1−α2v1, T̃2=T2−α2v2$ represent the values of the velocity and temperature once the singularity is removed in each case. Here, $wi$ and $Ti$ are the unknown vectors containing the values of the velocity and temperature, respectively, at each of the grid points in the ith subdomain.

The 1 × (N + 1)² row vectors $ri$ impose a first-order derivative condition at the location of the singularity. Since the singularity causes first derivatives to become infinite, imposing any (nontrivial) finite first derivative condition should ensure the correct value of α_i is selected to remove the singularity. Numerical experimentation appeared to confirm the arbitrary nature of the choice of first derivative condition. In both problems, we choose to impose continuity of the x derivative of the perturbation functions at the singularity (i.e., we impose $∂xw̃1(1,0)=∂xw̃2(−1,0)$ and $∂xT̃1(1,0)=∂xT̃2(−1,0)$⁠). As before, we convert this x derivative into ξ and η derivatives before discretization. Making use of a derivative condition in this way differentiates this method from that found in Ref. [45], which is an iterative method.

The (N + 1)² × 1 vector E contains −1 in every position corresponding to an interior point of the Chebyshev grids (according to the aforementioned enumeration), and 0 elsewhere. The (N + 1)² × 1 vector G contains the corresponding value of $Pw1$ (which is known once the velocity problem is solved) in every position corresponding to an interior point of the Chebyshev grids, and 0 elsewhere. The (N + 1)² × 1 vector H contains −1 in every position corresponding to the last row of the Chebyshev grids, the corresponding value of $Pw2$ in every position corresponding to an interior point, and 0 elsewhere. The −1s correspond to the nondimensional pressure gradient and heat flux in the former and latter cases, respectively.

The following are all matrices of size (N + 1)² × (N + 1)² and consist mainly of zeros: A is the resultant matrix of applying the collocation method to the velocity (or temperature) problems in subdomain 1; B (resp. C) is the matrix that comprises the right-hand side of Eq. (25) (resp. (28)) using the values of w₂ (resp. w₁); D (resp. F) is the resultant matrix of applying the collocation method to the velocity (resp. temperature) problems in subdomain 2. It is worth noting that points on the corners of the grid could have one of two boundary conditions allocated to them by these matrices. In this case, it mostly does not matter which of the two are imposed, although derivative conditions are avoided at the singularity to prevent infinite values.

It is important to note that in Eqs. (45) and (46), the apparent matrix multiplication of $A,B,C,D$⁠, and F with $ui$ and $vi$ is a shorthand. In particular, these quantities are the result of applying the continuous differential operators (corresponding to each matrix, both in the interior and boundary of the relevant subdomain) to the singular functions analytically before discretization. Note also that, since these singular functions are harmonic by design, all nonzero elements result from the boundary condition operators. Since the singular functions are handled analytically, the accuracy of the remaining numerical solution, and hence the full solution, is greatly improved.

In our implementation, we build these matrices and vectors in MATLAB and solve the system directly using the backslash command. We select the number of Chebyshev nodes, and check that this is sufficient by running the simulation with 25% more nodes in each direction, and checking that the averaged Nusselt number, defined in Eq. (17), does not change by more than 0.2%. In the case of $δ→1$⁠, we also found that the convergence rate was greatly improved by carrying out a further domain decomposition by dividing the domain at a constant value of y, usually chosen to be a small positive number, e.g., 0.1. In this case, there are four subdomains instead of two; however, the details of implementing these extra domains are very similar to those already given. We can also vary the number of Chebyshev nodes in the x and y directions independently. This is advantageous in cases involving subdomains with very large aspect ratios, in particular when $δ→0$ or 1.

Nusselt numbers were computed for the parameter ranges 0.5 < h < 20, 0.01 < ϕ < 0.99, and $−87.1 deg<θ<87.1 deg$⁠. The limits to the range of ϕ and h were primarily because extreme aspect ratios can cause the matrix problem to become ill conditioned. The ranges could have been enhanced by using more Chebyshev nodes (or adjusting the number for each axis); however, these ranges were sufficient to fully demonstrate the flow behaviors of interest. No singular behavior was noticed in $Nu¯$ as θ tends to $±90 deg$⁠. The limits to the range of θ were primarily due to the fact that transformation A (given in Appendix  A) becomes infinite as $θ→±90 deg$⁠.

where $Nu¯(0)$ and $Nu¯(1)$ are given in Kirk et al. [1]. Note that we use the protrusion angle θ as a small parameter, rather than the curvature, and that values of $Nu¯approx$ appear as straight, dashed lines in Figs. 6, 7, and 8.

In the h = 0.5 graph, we have no values for the approximate range $π/4<θ<π/2$⁠. This is because at these values of θ and this value of h, the meniscus has protruded so far into the channel, that it has touched the symmetry line y = h. This can be physically interpreted as the two menisci protruding into the channel from opposite surfaces (y = 0 and $y=2h$⁠) coming into contact with each other. At this point, we expect the flow regime to change completely, and therefore we do not attempt to carry out calculations for this range of θ.

The case h = 1, ϕ = 0.01 in particular illustrates the dramatic dependence of $Nu¯$ on θ, with $Nu¯$ increasing by more than a factor of 7, from $Nu¯≈0.26$ at θ = π/2 (protruding into the channel) to $Nu¯≈1.94$ at θ = −π/2.

Considering any fixed value of θ, one can note that our results are consistent with the expected observation that $Nu¯→0$ as $h→0$⁠. This is partly due to the choice of using the same definition of hydraulic diameter as Maynes and Crockett [30] in the definition of $Nu¯$⁠, namely 4h. This was chosen for consistency; however, there is a case to be made for including the meniscus geometry in the definition of the hydraulic diameter, which would result in a hydraulic diameter that tends to a nonzero constant as h tends to zero, and hence $Nu¯$ would not approach zero as quickly. However, $Nu¯$ would certainly tend to zero as $h→0$ regardless of the definition of the hydraulic diameter. This is due to the much diminished advection in the liquid gap between the two solid boundaries, resulting from the very small velocity values. This creates a large temperature drop between the ridge surface and the liquid region above the meniscus in order to drive the diffusion of heat into this region, where the heat can be advected away. This temperature drop becomes infinite as $h→∞$⁠, and therefore $Nu¯→0$⁠.

With regard to the accuracy of the asymptotic approximation of Kirk et al. [1], we see that it performs remarkably well. It is noteworthy that as h decreases, the accuracy of this asymptotic approximation improves for positive meniscus protrusion angles (protrusion into the channel). The trend for negative protrusion angles is less clear, and the accuracy seems to depend more on the solid fraction in this case.

In Fig. 7, we again plot $Nu¯$ for $−87.1 deg<θ<87.1 deg$⁠, but instead using values of ϕ close to 1 (i.e., $ϕ=0.9,0.95,0.99$⁠). Hence, this is the large solid fraction case, where the shear-free meniscus comprises only a small part of the boundary. In general, there is far less variation in $Nu¯$ as θ varies in this case, as a percentage change or in absolute terms, due to the smaller influence of the meniscus on the flow and temperature profiles. Note that as $ϕ→1$ or $h→∞$⁠, we recover the well-known result from the case of parallel plates, $Nu¯→8.23$⁠. This result is perhaps more intuitive if we consider that increasing h is equivalent to reducing the period to channel height ratio. When the period approaches zero for a fixed height, the velocity and temperature fields become one-dimensional, and the well-known result follows.

From Fig. 6, we can see that the observation made in Kirk et al. [1] that $Nu¯$ decreases as θ increases (for $h≥0.5$⁠) generally holds true at larger angles. This is largely due to the additional liquid area into which the heat can diffuse. However, in the h = 0.5 case in Fig. 6, we see some interesting behavior in $Nu¯$ as θ approaches the critical value at which the meniscus touches the centerline. We see a local minimum, followed by a small range of θ in which $Nu¯$ increases with θ. We propose the following explanation. The velocity values above the meniscus are always larger than above the ridge, normally causing the (nondimensional) velocity-weighted mean temperature to have values closer those found above the meniscus. Once the area of the liquid above the meniscus becomes noticeably smaller, which occurs when θ is sufficiently large in the h = 0.5 case, the velocity weighted mean temperature becomes more similar to values seen over the ridge. The liquid temperature values above the ridge are closer to those at the surface of the ridge, and hence $Nu¯$ may begin to increase.

By increasing the solid fraction ϕ, we introduce two effects: (i) we increase the rate of heat flow into the channel since we are increasing the area of the ridge (ii) we decrease the lubrication. Effect (i) increases $Nu¯$⁠, but not by as much as one might assume. Since there is a larger rate of heat flow into the channel, a larger temperature drop is required to sustain it, partially counteracting the increase in $Nu¯$ due to the increased range of integration. Effect (ii) decreases $Nu¯$ since less heat can be advected away in the vicinity of the ridge. In Figs. 6 and 7, we can see the general trend that an increase in ϕ consistently leads to an increase in $Nu¯$⁠. It follows that effect (i) is more significant than effect (ii). However, the important role played by lubrication is indicated by the fact that a tenfold increase in ϕ (between $ϕ=0.01$ and ϕ = 0.1) leads to only a comparatively modest increase in $Nu¯$⁠.

Note that $Nu¯$ decreasing with ϕ is consistent with our previous observation: that textured microchannels have greater convective thermal resistances than parallel plate channels. It is worthwhile reiterating that this increase in convective thermal resistance can be more than offset by the decrease in caloric thermal resistance associated with the enhanced volumetric flow rate. As predicted by Lam et al. [6], this is particularly true in the case of liquid metals (due to large thermal conductivities). Therefore, for cooling applications, one might consider the results presented here in combination with those given in purely hydrodynamic studies [43,44] to acquire a predicted change in total thermal resistance. In applications for which it is desirable to retain the temperature of the liquid (for example, see Rosengarten et al. [9]), there is no compromise: both the flow rate enhancement and the decrease in convective heat transfer are beneficial.

In Table 1, we provide the range of θ for which the approximation given in Kirk et al. [1] is accurate to either 1% or 10%, compared to the values in the current work. This quantifies the remarkable accuracy of the asymptotic expansion. The values in the table were calculated using the same data as Figs. 6 and 7. In particular, $Nu¯$ was computed for 100 equally spaced values of θ ranging from $−87.1 deg$ to $87.1 deg$ (inclusive). Linear interpolation between these values of θ then gives the values found in the table (rounded to one decimal place). The elements given as 87.1 (resp. −87.1) indicate that the specified accuracy is achieved for all θ in the range $0≤θ≤87.1 deg$ (resp. $−87.1 deg≤θ≤0$⁠). The entries of the table given as N/A indicate that the parameter combination is impossible. This happens in the case of large positive θ, small ϕ, and small h as the meniscus touches the symmetry line.

In Fig. 8, we provide plots of $Nu¯$ as calculated in the current work alongside approximate values obtained using the $h→∞$ approximation¹ in Kirk et al. [1]. In Table 2, we provide the range of θ for which this approximation is accurate to 10%, compared to the values in the current work (calculated in the same way as in the previous table). This asymptotic solution is explicit, and therefore is much easier to implement; however, using this approximation instead of the solution for general h comes at the cost of some lost accuracy. Nonetheless, this table shows the remarkable accuracy of the $h→∞$ asymptotic expansion. As before, the elements given as 87.1 (resp. −87.1) indicate that the specified accuracy is achieved for all θ in the range $0 deg≤θ≤87.1 deg$ (resp. $−87.1 deg≤θ≤0 deg$⁠). The entries of this table given as N/A indicate that there is insufficient agreement at θ = 0 from which a range of accuracy can be meaningfully extended. This only happens at h = 5 and ϕ = 0.01. In general, the larger h and ϕ, the more accurate this approximation.

In Fig. 9, we give contour plots for the average Nusselt number as a function of ϕ and h at $θ=±87.1 deg$⁠. These values were chosen as they are the extremes of the values of θ for which $Nu¯$ was calculated. Note that, while close, these are not exactly at $±90 deg$⁠, because one of the coordinate transformations (transformation A in Appendix  A) becomes singular as θ approaches $90 deg$⁠. We do not encounter any singular behavior in $Nu¯$ as θ approaches $90 deg$⁠; hence, we view our results as representative of the behavior of $Nu¯$ at the extrema of the range of possible values of θ.

We have developed a numerical method to accurately compute the Nusselt numbers for Poiseuille flow over isoflux parallel ridges and arbitrary meniscus curvature. A primary objective of this paper was to use this method to assess the extent to which the asymptotic approximation of Kirk et al. [1] is accurate for protrusion angles which are not small.

We present our computed values of the averaged Nusselt number as a function of the meniscus protrusion angle θ, for various values of channel height h and solid fraction ϕ. We demonstrate the importance of meniscus curvature on heat transfer, particularly when the solid fraction is small; e.g., for $ϕ=0.01, Nu¯$ more than doubles as the protrusion angle decreases from $π/2$ to $−π/2$⁠, for the entire range of h studied.

These computations are given alongside the asymptotic approximations of Kirk et al. [1] for comparison. We show that this approximation agrees with the numerical work remarkably well over the range $−π/2<θ<π/2$⁠. In particular, 90% accuracy is typical for positive protrusion angles (into the channel) of less than $50 deg$⁠, whereas 90% accuracy is seen over the full range of negative protrusion angles (into the groove), for $h≥1$⁠.

We also compare our numerically derived results with the additional $h→∞$ approximation given in Kirk et al. [1]. This approximation is given explicitly, and is much easier to calculate. We find that for $h≳5$ (or h > 5 if ϕ is very small), this approximation typically has good accuracy for $θ⪅45 deg$⁠. Remarkably, for these parameter values, the asymptotic approximation remains good for the entire range of negative protrusion angles (into the groove).

• National Science Foundation, Division of Chemical, Bioengineering, Environmental, and Transport Systems (Grant No. 1402783).

• Engineering and Physical Sciences Research Council (Grant Nos. EP/K041134/1 and EP/L020564/1).

• Imperial College President's Ph.D. Scholarship.

A =

domain in which calculations are carried out

c_p =

specific heat of the liquid

h_m(x) =

dimensionless distance from the center-line of the channel (in the y direction) to the meniscus

h* =

dimensional channel height

k =

thermal conductivity of the liquid

n =

unit (outward) normal vector to the meniscus

$Nu¯$ =

average Nusselt number

P =

dimensionless constant defined as $ϕ/Q$

p* =

dimensional liquid pressure

$pg*$ =

dimensional gas of vapor pressure inside the grooves

Q =

dimensionless volume flow rate through A

Q* =

dimensional volume flow rate through A

$qsl″$ =

constant heat flux per unit area imposed on the solid–liquid interface

r =

radial coordinate in the polar coordinate system in which $(δ,0)$ is the origin

r_h =

height of the top of the tee-shaped structure comprising the ridge

r_w =

width of the top of the tee-shaped structure comprising the ridge

$R*$ =

dimensional radius of curvature of the meniscus

S =

the shear free liquid boundary, i.e., the meniscus

s_w =

width of the stem of the tee-shaped structure comprising the ridge

T =

dimensionless liquid temperature

T_i =

dimensionless temperature values given in terms of the transformed coordinates ξ_i and η_i in the ith transformed subdomain

$Ti$ =

vector containing the values of T_i at the Chebyshev points

T* =

dimensional liquid temperature

$Tm*$ =

dimensional mixed-mean liquid temperature

$Tsl*$ =

dimensional liquid temperature on the solid–liquid interface

u =

singular part of the velocity field

u_i =

singular part of the velocity field given in terms of the transformed coordinates ξ_i and η_i in the ith transformed subdomain

$ui$ =

vector containing the values of u_i at the Chebyshev points

v =

singular part of the temperature field

v_i =

singular part of the temperature field given in terms of the transformed coordinates ξ_i and η_i in the ith transformed subdomain

$vi$ =

vector containing the values of v_i at the Chebyshev points

w =

dimensionless streamwise velocity

w_i =

dimensionless streamwise velocity values given in terms of the transformed coordinates ξ_i and η_i in the ith transformed subdomain

$wi$ =

vector containing the values of w_i at the Chebyshev points

$w*$ =

dimensional streamwise velocity

x =

dimensionless spanwise coordinate

x* =

dimensional spanwise coordinate

y =

dimensionless transverse (wall-normal) coordinate

y_m(x) =

y coordinate of the meniscus

y* =

dimensional transverse (wall-normal) coordinate

z =

dimensionless streamwise coordinate

z* =

dimensional streamwise coordinate

2a* =

dimensional width of each groove

2d* =

dimensional period of the ridges

Greek Symbols

α =

thermal diffusivity of the liquid

α₁ =

strength of the velocity singularity

α₂ =

strength of the temperature singularity

γ =

constant surface tension of the meniscus

δ =

half of a dimensionless groove width referred to as liquid (or cavity) fraction

η_i =

vertical coordinate used in the ith transformed subdomain

θ =

protrusion angle of the meniscus into the channel

θ_a =

advancing contact angle

λ =

the power of the singularity in velocity field in the sense that $w∼rλ$ as $r→0$

μ =

dynamic stress viscosity of the liquid

ξ_i =

horizontal coordinate used in the ith transformed subdomain

ρ =

density of the liquid

ϕ =

solid fraction, defined as $1−δ$

$φ$ =

angular coordinate in the polar coordinate system in which $(δ,0)$ is the origin

where $Hm(ξ1)=hm(x)$⁠.