An analytical methodology to characterizing the effects of heat transport in internal laminar flows over ridged patterns, mimicking superhydrophobic surfaces, is indicated. The finite slip velocity on such surfaces and the thermal conductivity characteristics of the constituent material are both shown to modify the convective heat transport in the fluid. We use an effective medium approach to model the lowered thermal conductivity caused by the presence of air in the ridge interstices. The proposed analytical solutions for fully developed flow were verified through comparison with numerical simulations for a periodically ridged geometry in laminar flow. While the convective heat transport and the Nusselt Number (*Nu*) increase due to the modified fluid velocity profile on superhydrophobic surfaces, the decrease in the thermal conductivity of the substrate may play a larger role in determining the overall heat transfer in the channel.

## I. INTRODUCTION

Quantitative understanding and control of heat transfer in micro- and nano-fluidics are essential for a variety of physical, chemical, and biological applications,^{1} such as temperature gradient focusing for electrophoresis,^{2} mixing,^{3,4} and polymerase chain reaction for DNA amplification.^{5,6} A key to efficient thermal control of many such phenomena is the ability to adjust effective convective and conductive heat transfer coefficients. In this context, micro- and nano-patterned surfaces provide a promising way to modulate heat transfer in flowing fluid and the containing solid channel. Air pockets trapped in the valleys of a superhydrophobic (SH) patterned surface reduce interfacial frictional drag and induce a finite slip velocity at the interface^{7–11}—assuming that the viscosity of the air (*μ*_{air}) is sufficiently low compared to that of the fluid. While significant attention has been given to static phenomena, the understanding of dynamic fluid flow, from the laminar to the turbulent regimes, needs more exploration.^{11–13} A related aspect that also deserves study is the effect of the SH surface on the thermal transport across the surface/fluid interface. It may, for example, be assumed that the trapped air would affect heat transfer through both (a) modifying the flow dynamics, as well as (b) reducing the heat conduction.

SH characteristics generally require a large (typically greater than 150°) contact angle between the fluid and the surface and a low contact angle hysteresis.^{12,14,15} Consequently, when a wetting liquid does not penetrate the interstices of the roughness, e.g., under relatively low-pressure conditions so as not to exceed the Laplace pressure, the resultant state has been termed the Cassie state. For fluid in the regions of direct contact with the solid surface, the no-slip boundary condition is pertinent,^{16} while in regions with a fluid/air boundary, a shear-free condition applies. When averaged over the entire composite area, the net enhanced velocity has been considered through an apparent slip^{11} and quantified through the slip velocity (*v _{s}*).

With the overall goal of understanding heat transport phenomena in multiphase flow over structured surfaces, and to apply such understanding for practical applications requiring strict thermal control, we aim to show the sensitivity of the macroscopic heat transfer to temperature gradients in both the fluid (e.g., water will be the prototype) and solid. It will be indicated that tuning the air fraction of a ridged surface as well as the conductivity of the substrate could in concert modulate the *net* heat transfer rate.

We analyze the heat transport through a model SH surface by treating the air/solid composite through an effective medium approximation (EMA) based approach. The EMA is employed since the roughness length scales, as related to the air pocket dimensions, are typically small and the averages over the channel length, such as velocity (as related to the flow rate) and heat transfer coefficients, are measured at a much larger scale. Specifically, we will consider, in addition to the average slip velocity at the fluid/surface boundary, an effective thermal conductivity, *k*_{eff}, for the entire SH composite region. Consequently, we seek to decouple the hydrodynamic effects related to an altered fluid velocity profile and related convective heat transport, from the thermal effects of the trapped air pockets. Our approach enables the obtaining of an analytical solution for the temperature profile, heat transfer coefficients, and the Nusselt number (*Nu*).

Previous heat transport studies on SH surfaces in micro-channels focused on solving, through numerical or analytical methods, the coupled mass, momentum, and energy equations in the fluid, assuming constant temperature^{17} or heat flux^{18} at the substrate boundary. While it was indicated that the heat transfer could be reduced, an averaging over the SH surface was not employed, and relatively little attention was paid to specific solid characteristics. We aim to extend such studies through explicitly incorporating the thermal conductivity of the solid and the surface roughness needed to achieve SH characteristics. This will serve as the basis for analyzing more generalized cases of heat flux variation within the flowing fluid and enable insight into the pathways of heat flow and loss.

We focus on laminar flow through rectangular shaped micro-channels. We will, in the present work, extend our study to cases where the axial conduction through the micro-channel walls does not have a large effect on the fluid temperature profile and heat flux. Essentially, it was assumed that the heat transfer is primarily in the transverse (wall-normal) direction to the liquid flow. Such a consideration is parameterized through a low ratio of the transverse thermal resistance in the liquid to the lengthwise thermal resistance in the channel walls, through^{19} *β* = (*H*/*k*_{w})(*k*_{s}/*L*_{tot}) where *H* is the half-height of the channel, *k*_{w} and *k*_{s} are the thermal conductivities of the fluid (water) and the substrate, respectively, and *L*_{tot} is the total length of the channel. In commonly fabricated micro-channels with a small *H*/*L*_{tot} ratio, the *β* would generally be small (typical values considered here are in the range of *β* = 0.001 − 0.1, assuming an *L*_{tot} = 10 cm). The rationale for studying low *β* flows is that at larger *β* values, the heat transfer in the substrate and ridged region would be a significant fraction of the heat transported in the fluid. Consequently, the influence of the SH surface on fluid flow and heat transport would not be as well defined.

## II. METHODS

The investigated model system is comprised of an effectively two-dimensional fluid flow channel (of height: 2*H*, see Fig. 1) between large parallel plates. The outer surfaces of the plates are exposed to ambient air (with convective heat transfer coefficient: *h*_{ext}), at a constant temperature: *T*_{∞}. The inner surfaces at the top and bottom are ridged (see Fig. 1), so as to mimic SH-like character.^{20,21}

A periodic ridge pattern was used (see Fig. 2), with ridge height *d*, with air (of thermal conductivity *k*_{a} = 0.028 W/mK) filling the interstices between the ridges, and assuming a flat air meniscus.^{18,22} A linear air fraction was defined through *ϕ*_{a} = *l*_{a}/*L* (where *l*_{a} is the width of the air pocket and *L* is the combined width of the solid rib of length *l*_{s} and the air pocket, i.e., *L* = *l*_{s} + *l*_{a}). The substrate is of thickness *D* and thermal conductivity *k*_{s}. In the periodically rough region corresponding to the SH surface, an effective thermal conductivity, *k*_{eff}, is defined by considering the solid ridges and air pockets as thermal resistances to heat loss as follows:

Pressure driven laminar flow, of the Poiseuille type, of water (of thermal conductivity *k*_{w} = 0.65 W/mK with a temperature at the channel input of *T* > *T*_{∞}) was considered. Generally, the heat flux along the channel would be variable, since as the water flows over the SH surface, the temperature could drop yielding a decaying heat flux.^{23}

For laminar incompressible flow, in the steady state at constant density and viscosity, the solution of the Navier-Stokes equation: $ \u2202 p \u2202 x =\mu \u2202 2 u \u2202 y 2 $, for the velocity *u*, with slip boundary conditions corresponding to the Navier condition ($ v s =b \u2202 u \u2202 y $, with *b* as the slip length) at the bottom and top SH surfaces was derived to be

For the no-slip case, *b* = 0 and the standard Poiseuille flow solution can be obtained. Alternately, we indicate a normalized velocity (also see Subsection 1 of the Appendix)

where *y*^{∗} = *y*/*H*, *b*^{∗} = *b*/*H*, and $ u \u0304 $ is the average velocity in the channel,

In Eq. (3), *A* = 1 + 2*b*^{∗} can be considered equivalent to the square of an effective height as, per the slip length definition, the velocity would tend to zero at one slip length away from the wall, while the coefficient *B* = 1 + 3*b*^{∗} reflects the increase in the average velocity.

To obtain an average *b*, for square ridges oriented perpendicular to the fluid flow, we use a previously obtained relation by Philip,^{24} where

Such an expression was obtained for alternating bands of no-slip and shear-free boundaries in Stokes flow and would be applicable for low Reynolds number (*Re*) flow and *μ*_{air} ≪ *μ*_{w}. In the case of higher *Re*, we could employ a modified expression such as the one developed by Woolford *et al.*^{25} showing a small decrease in the *b* with increasing *Re* (yielding a value of ∼0.98 of that predicted by Eq. (5) at *Re* = 100 to 0.74 at *Re* = 1000 for the length scales considered here). Equation (5) shows that for a channel height 2*H* of 300 *μ*m and a horizontal roughness length scale *L* of 80 *μ*m, *b*^{∗} varies from 0.0 to 0.2 as *ϕ _{a}* increases from 0.0 to 0.94, resulting in values of

*v*

_{s}of up to $0.3 u \u0304 $.

The temperature distribution, *T*(*x*, *y*), in the fluid was obtained through solving the thermally fully developed form of the energy equation^{26}

where $ d T d x $ indicates the stream-wise temperature gradient and *α* is the thermal diffusivity. The viscous dissipation term has been neglected, since the Eckert number ($Ec= u \u0304 2 / C p \Delta T$, with *C _{p}* as the specific heat capacity of water and Δ

*T*as the characteristic temperature difference) was determined to be 10

^{−7}. The axial heat conduction in the fluid has also been neglected due to the Péclet number (

*Pe*)—a measure of the thermal energy convected by the fluid to that conducted within the fluid—being large,

^{26,27}such that the effect was small for

*Pe*(=

*RePr*) > 10, where

*Pr*is the Prandtl number (for water, 2 <

*Pr*< 7).

We used a normalized temperature, defined per,^{23,26}

where *T*_{∞} was the ambient temperature and *T*_{CL}, the centerline temperature (Fig. 2). Additionally, we use a non-dimensional lengthwise temperature gradient

It has been previously discussed^{23} that for a heated flow over a surface, *λ* is also a measure of the characteristic heat flux (*q*) decay constant, i.e., *q* ∼ *e*^{−λx}. Using *x*^{∗} = *x*/(*H* × *Pe*) with $Pe= u \u0304 2H/\alpha $, the energy equation was rewritten as (see Subsection 2 of the Appendix)

The differential equation was solved with boundary conditions (i) assuming a symmetric temperature profile of the liquid around the centerline, i.e., $ \u2202 \theta \u2202 y \u2217 | y \u2217 = 0 =0$, along with (ii) *θ*|_{y∗=0} = 1, following from the definition of *θ*. We assumed a power series solution for *θ* of the form (see Subsection 3 of the Appendix)

where *C*_{0} = 1, *C*_{1} = − 3*λA*/8*B*, and *C _{n}* = 3

*λ*(

*C*

_{n−2}−

*AC*

_{n−1})/(4

*B*(2

*n*)(2

*n*− 1)). The series solution for

*θ*, to the fourth order, is

The *λ* was obtained through equating the conductive heat flux at the fluid-wall interface, to the convective loss, from the liquid, through^{23}

where *h*_{w} encompasses the *net* conductance from the inner surface of the wall (*y*^{∗} = 1) to the ambient, see Fig. 2, and was calculated from the series sum of the thermal resistances between the inner wall surface and ambient from

Consequently, using *R*^{∗} = *R*_{w}*k*_{w}/*H*, we obtain the following:

The solution of Eq. (14) by evaluating the temperature and temperature gradient at the wall from Eq. (11) results in a value for *λ* of order 1/*R*^{∗} (see Subsection 4 of the Appendix).

We have shown previously through Eq. (3) that the change in average velocity and flow rate due to the slip was *B*(=1 + 3*b*/*H*). For *b* in the range of tens of micrometers, for a significant impact on the flow conditions, *H* would need to be of the order of hundreds of micrometers or less (in our study, we use an *H* = 150 *μ*m). From Eq. (13), for typical *d* = 50 *μ*m, *D* = 100 *μ*m, and *k _{s}* = 0.5 W/mK, we note that

*h*

_{ext}(typically 1 W/m

^{2}K - 250 W/m

^{2}K)

^{28}dominates

*R*

_{w}. Even for a large

*h*

_{ext}of 250 W/m

^{2}K and

*H*= 800

*μ*m,

*R*

^{∗}> 3.3, which we use as a benchmark value. Such

*R*

^{∗}values result in

*λ*< 1, with the implication that the temperature series expression in Eq. (11) can be truncated to the fourth order with minimal loss to accuracy. The form in Eq. (11) was verified against the numerical integration of Eq. (9) for a

*R*

^{∗}= 4 with the obtained

*θ*

_{w}and

*dθ*/

*dy*accurate to within 0.5%.

A heat transfer coefficient: *h _{c}* for the internal flows was then defined as

and can be obtained from Eq. (11) by calculating a mass-averaged mean temperature,

and through evaluating ∂*θ*/∂*y*^{∗} from Eq. (11). Consequently,

A corresponding Nusselt number: $Nu= h c D h k w $, where *D*_{h} = 4*H* is the hydraulic diameter, can be obtained from Eq. (17). The *Nu* resulting from Eq. (17), with *b*^{∗} = 0 (no-slip) and *λ* = 0 (no heat flux decrease in the horizontal direction), reduces to the typically quoted textbook value^{26} of *Nu* = 8.235 for parallel plates with a constant heat flux from fluid to ambient.

However, such a conventional definition of *Nu* in internal flows incorporates the effects of the altered velocity profile due to the patterned surface and only indirectly considers the modified thermal conductance of the channel wall (Eq. (14)). To further understand and quantify the net impact of the SH composite surface on the heat transfer from the liquid to the ambient, we define another heat transport coefficient: *h*_{sh} as follows:

The *h*_{sh} incorporates the change in the convective heat transport coefficient *h*_{c} as well as the influence of the air pockets on the conductive heat transport.

It has been assumed that the *k*_{eff} would be the only pertinent parameter for heat transport through the ridged region, as the convective heat transfer within the interstices can be neglected. This was justified on the basis that the equivalent resistance per unit area, =*d*/*k*_{a} of the air pocket, with, e.g., *d* = 40 *μ*m would be ∼0.0014 m^{2} K/W. This would correspond to an equivalent heat transfer coefficient of ∼665 W/m^{2} K, a relatively high value for air, and requires a large air velocity. Such conditions may only be important when the *Re* for the fluid flow in the channel corresponds to a turbulent regime.

We carried out numerical simulations on a representative case, using COMSOL Multiphysics (computational fluid dynamics coupled with heat transfer solver), to corroborate our analysis. We chose a representative case with a *Pe* = 60, within the limits of applicability. The fluid inflow temperature was 40 K over ambient which results in an average *Pr* = 3. The water flow was driven by a pressure gradient of 1000 Pa/m, resulting in an approximate velocity of 1.8 cm/s, and a *Re* of 20. The SH surface ridge parameters (Fig. 2) were in the range of 0 *μ*m > *l _{a}* > 80

*μ*m and

*d*= 40

*μ*m, chosen to generate sufficient slip yet within the practical limitations imposed by air pocket stability. The

*h*

_{ext}was chosen as 80 W/m

^{2}K, an approximate value for forced air convection, and a

*k*

_{s}= 0.5 W/mK corresponding to a thermally non-conductive substrate such as polydimethylsiloxane (PDMS) (

*β*= 0.001). Additionally, a thermally conductive substrate case with a

*k*= 45 W/mK, as for steel (

_{s}*β*= 0.1), was simulated.

## III. RESULTS

It was seen that the main variables defining the heat transfer were the air fraction *ϕ*_{a} and the corresponding normalized slip length *b*^{∗}, with the normalized thermal resistance *R*^{∗} and *λ* having a small effect. Fig. 3 indicates the variation of the *ϕ*_{a} and the resulting *k*_{eff} for different *b*^{∗} considered through Eqs. (1) and (5). It was evident that a large *ϕ*_{a} yields large slip for the given conditions. However, the heat transfer is impeded by the concomitantly lowered *k*_{eff}, reducing the conductivity of the composite region from 0.5 W/mK to as low as 0.06 W/mK at *ϕ _{a}* = 0.94.

The temperature and velocity profiles in the fluid obtained from Eq. (11) with representative *R*^{∗} and *b*^{∗} values are shown in Figure 4, showing the decrease in temperature from the maximum at the centerline. The temperature gradient in the fluid was higher for small *R*^{∗} values, as may be deduced from Eq. (14). A small (/large) thermal resistance between the heated fluid and the outside air would result in a large (/small) temperature change within the fluid. Increasing the *b*^{∗} at a constant *R*^{∗} makes the temperature profile in the fluid more uniform due to an increased *v*_{s}. Such enhanced uniformity mimics the change in the velocity profiles (Fig. 4 inset), which becomes more plug-like in flow over SH surfaces. Interestingly, while the overall temperature profile in the fluid increases for increasing *b*^{∗}, the mass-averaged temperature *θ*_{m} actually decreases, due to an increased flow close to the wall at higher slip.

The temperature profiles in the fluid and the channel walls are shown in Figure 5 for the case of *k*_{s} = 0.5 W/mK with the typical length scales used in the numerical analysis and simulations. An increased relative temperature in the fluid was again seen with increasing *b*^{∗}. However, a progressively larger temperature gradient in the SH surface composite region was also observed. Consequently, there was a marked decrease in temperature through the SH zone for increasing *b*^{∗}. A concomitant increase in *R*^{∗}, through the decrease in *k*_{eff} from the increased air fraction, is indicated in the inset of the figure. The temperature gradient in the SH composite region from the simulations was used to estimate an effective thermal conductivity for each surface; Table I compares the *k*_{eff} used in the analytical derivation from Eq. (1), with the estimate from the simulations. A close match to within 12% indicates the validity of the EMA based approach for understanding heat transfer through the composite region.

Slip length . | Air fraction . | Analytical . | Simulation . |
---|---|---|---|

b* . | ϕ_{a}
. | k_{eff} (W/mK)
. | |

0.00 | 0.000 | 0.500 | 0.500 |

0.05 | 0.625 | 0.205 | 0.180 |

0.10 | 0.801 | 0.123 | 0.110 |

0.15 | 0.891 | 0.080 | 0.074 |

Slip length . | Air fraction . | Analytical . | Simulation . |
---|---|---|---|

b* . | ϕ_{a}
. | k_{eff} (W/mK)
. | |

0.00 | 0.000 | 0.500 | 0.500 |

0.05 | 0.625 | 0.205 | 0.180 |

0.10 | 0.801 | 0.123 | 0.110 |

0.15 | 0.891 | 0.080 | 0.074 |

An increase in the *Nu* ($\u223c q T m \u2212 T w $) with increasing *b*^{∗} may be initially expected, as convection was enhanced due to the increased slip velocity, as illustrated in Fig. 6. Also, from Eq. (15), the *Nu* was seen to be a function of both the temperature gradient at the wall (in the numerator), as well as the temperature difference between fluid and wall (in the denominator). However, there seems to be a subtle interplay; it was observed that the temperature difference *θ*_{m} − *θ*_{w} decreases much more (of the order of 20% between *b*^{∗} = 0 and *b*^{∗} = 0.2) compared to the decrease in the temperature gradient ($ d \theta d y \u2217 | y \u2217 = 1 $, of around 10%) resulting in an overall increase in the *Nu* (Fig. 6) from 8.21 to 9.42. The increase in *Nu* arises primarily from the increase in slip velocity. Increasing air fraction also works to increase *R*^{∗} (Fig. 5, inset) leading to a decrease in *λ* and an additional increase in the *Nu*. However, the effect of *λ* was very small (Eq. (17)) for our representative case, and generally for the considered values (i.e., with *R*^{∗} > 3.3).

A decrease in the *h*_{sh} (as defined in Eq. (18)) with *b*^{∗} is shown in Fig. 7. The variation is in comparison to the no-slip value/baseline value (where *b*^{∗} = 0). While the *h*_{c}/*Nu* was increased following the previous discussion, the *k*_{eff} was involved in the overall decrease of *h*_{sh}. A decrease in the net heat transfer then seems to be the rule for thermally non-conductive substrates. However, the magnitude of the decrease was dependent on the specific values of the *k*_{s} and *d*; e.g., with *k*_{s} = 0.5 W/mK and *d* = 40 *μ*m, the increased convective heat transfer was offset by the *k*_{eff} resulting in a decreased *h*_{sh}, yet the heat transfer rate decrease could be limited by a smaller *d* or higher conductivity substrate *k _{s}*.

The results for the thermally conductive substrate (*k*_{s} = 45 W/mK) are shown in Figs. 8 and 9. There was a sharper decrease in the *k*_{eff} with *b*^{∗} in the simulation results than was predicted through Eq. (1), shown in Fig. 8. The variation in *k*_{eff} was presumably from the additional effect of conduction in the streamwise direction in the composite roughness zone, which is more significant with increasing |*k*_{s} − *k*_{a}|.^{28} In this case, the effective thermal conductivity of the interface will be between the Wiener bounds, of the sum of the resistances in parallel as in Eq. (1) (heat flow in the transverse direction) and that in series (heat flow in the streamwise direction).

The results for *h*_{sh} for the thermally conductive substrate are shown in Fig. 9. The *h*_{sh} variation in *b*^{∗} should be understood as due to the opposing effects of an increase in *h*_{c} and a decrease in *k*_{eff}. With an increase in the thermal conductivity of the substrate, the relative magnitude of *d*/*k*_{eff} is reduced, such that the increase in *h*_{c} should maintain *h*_{sh} relatively unchanged (see Fig. 9). However, as the simulations show, the increased reduction in *k*_{eff} due to the streamwise effects results in a decreased heat transfer rate.

A coefficient of performance (COP) was postulated as the ratio of the relative change in the heat transfer coefficient compared to a non-SH surface (no-slip) and the relative change in the pumping power as compared to a non-SH surface as follows:

where Δ*P* is the pressure drop across the channel, and *Q* is the volumetric flow rate. The COP variation for the *k*_{s} = 0.5 W/mK and *k*_{s} = 45 W/mK cases is shown in Fig. 10, both resulting in decreasing COP with increasing *b*^{∗}.

## IV. CONCLUSIONS

The EMA based approach yields physical insights into the nature of thermal transport (including temperature profiles and heat flux) through SH patterned surfaces in laminar flow. The analytical methodology is capable of being extended to the study of other surfaces (ridges parallel to flow, posts, etc.) through the use of modified slip length and effective thermal conductivity relations. We have indicated the subtle balancing of the effects of an increased velocity (which enhances the *h*_{c}) with the thermal insulation provided by air. For low *k*_{s} (of ∼1 W/mK, as typical to polymeric substrates), the effect of a SH surface seems to be the reduction of the overall heat transfer in micro-channels, regardless of the air fraction. For substrates with a larger *k*_{s}, the net heat transport also decreases with air fraction, though to a lesser degree. However, in this case, there was a deviation from the form of Eq. (1) (for the EMA *k*_{eff}) due to the contributions from the multidimensional effects (e.g., both transverse and streamwise). The study of such influences with respect to the roughness characteristics would be important for practical microfluidics applications and forms the basis for a fertile field of inquiry.

## Acknowledgments

D. Moreira appreciates the support from the Jacobs Fellowship program at UC, San Diego, Professor K. Seshadri for insight into convective transport, and Dr. S. Hobbs for PDE help. We are also thankful for financial support (through CMMI 1246800) from the National Science Foundation (NSF).

### APPENDIX: DERIVATION OF THE FLUID AND HEAT TRANSPORT EQUATIONS

#### 1. Non-dimensionalized velocity

In order to obtain the non-dimensionalized form of the velocity expression (Eq. (2)), we obtain the average velocity in the channel,

The normalized velocity, $ u \u2217 =u/ u \u0304 $, becomes

which can be written in normalized variables *y*^{∗} = *y*/*H*, *b*^{∗} = *b*/*H* as

Introducing two coefficients, *A* = 1 + 2*b*^{∗} and *B* = 1 + 3*b*^{∗}, we can write *u*^{∗} as

#### 2. Non-dimensionalized energy equation

Starting with the energy equation in dimensional form (Eq. (6)),

and substituting using the non-dimensional variables $ u \u2217 =u/ u \u0304 $, *y*^{∗} = *y*/*H*, and *x*^{∗} = *x*/(*H* × *Pe*), where $Pe=Re\xd7Pr=2H u \u0304 /\alpha $, *θ* = (*T* − *T*_{∞})/(*T*_{CL} − *T*_{∞}), we obtain the following:

Noting that for fully developed flow, *θ* = *f*(*y*^{∗}) only and *T _{CL}* =

*f*(

*x*

^{∗}) only, the partial derivative becomes

Substituting the definition of *λ* = − 1/(*T*_{CL} − *T*_{∞})(∂*T*_{CL}/∂*x*^{∗}), we obtain the final, non-dimensionalized form

#### 3. Series solution for the energy equation

The solution to the energy equation (Eq. (9)) can be represented by a power series in the form

and the second derivative with respect to y becomes

Substituting into Eq. (9) results in

Shifting the coefficients in order to write all terms as summations of *y ^{i}* and starting the summation at the same starting index, we obtain the following:

Requiring the coefficients of each power of y to be equal to zero results in

The power series solution becomes

Applying the boundary condition

results in all odd terms being zero and *a*_{1} = 0. From the definition of *θ*, *a*_{0} = 1, such that the final form is

Since only even terms are left, we change the summation index from *i* = 2*n* to get the final form

where *C*_{0} = 1, *C*_{1} = − 3*λA*/8*B*, and *C _{n}* = 3

*λ*(

*C*

_{n−2}−

*AC*

_{n−1})/(4

*B*(2

*n*)(2

*n*− 1)).

#### 4. Derivation of *λ*

To obtain *λ*, the boundary condition at the wall equating the heat flux conducted through the fluid with the heat flux to the ambient temperature, Eq. (12),

In non-dimensionalized form, with *θ _{w}* = (

*T*−

_{w}*T*

_{∞})/(

*T*

_{CL}−

*T*

_{∞}), we obtain the form

Using *R*^{∗} = *k*_{w}/(*h*_{w}*H*), a non-dimensionalized resistance from wall to ambient, we obtain the following:

From the temperature expression Eq. (11), we can obtain both the temperature gradient at the wall (LHS of Eq. (D3)) and *θ _{w}* (RHS of Eq. (D3)). These are

Substituting these values into the non-dimensionalized form (Eq. (D3)) results in

and solving the resultant quadratic equation for *λ*, we obtain the following:

where *E* = 4*B*(6*A* − 1 + 8*BR*^{∗}), *F* = 16*B*^{2}[(1 − 6*A* − 8*BR*^{∗})^{2} − 8(3*A*^{2} + 12*A*^{2}*R*^{∗})], and *G* = 3*A*^{2} + 12*A*^{2}*R*^{∗}.