2.1.1. Outer region

The outer region of the domain is where $X=O(1)$, $\tilde {y}=\mathrm {ord}(1)$, i.e. strictly of order unity as $\epsilon \rightarrow 0$. It is sufficiently far from the mixed boundary condition imposed at $\tilde {y}=0$ that the velocity field is one-dimensional, i.e. $\tilde {w} = \tilde {w}(\tilde {y})$, as justified henceforth. Equation (2.11) reduces, to leading order, to $\partial ^{2}\tilde {w}/\partial X^{2}=0$. Integrating, satisfying the symmetry boundary condition in $X$ and integrating again shows that, to leading order, $\tilde {w}=\tilde {w}(\tilde {y})$. Thus, the $\partial ^2\tilde {w}/\partial \tilde {y}^2$ and $-2$ terms in (2.11) balance, implying that $\tilde {w}=O(1)$. We thus, as an ansatz, expand the velocity as a regular power series in $\epsilon$ as per

(2.16)\begin{equation} \tilde{w}=\sum_{n=0}^{\infty}\epsilon^{n}\tilde{w}_{n}+ \mathrm{e.s.t.},\quad\epsilon \rightarrow 0,\end{equation}

where $\tilde {w}_{n}=O(1)$ for all $n\geqslant 0$ and e.s.t. denotes exponentially small terms that are smaller than any algebraic order of $\epsilon$. We do not attempt to solve for these terms, and only include them for error quantification. Substituting (2.16) into (2.11) yields an equation at each (algebraic) order in $\epsilon$ as per

(2.17)$$\begin{gather} O\left(\epsilon^{{-}2}\right):\quad \frac{\partial^{2} \tilde{w}_{0}}{\partial X^2}=0 \end{gather}$$

(2.18)$$\begin{gather}O\left(\epsilon^{{-}1}\right):\quad \frac{\partial^{2} \tilde{w}_{1}}{\partial X^2}=0 \end{gather}$$

(2.19)$$\begin{gather}O\left(\epsilon^{0}\right):\quad \frac{\partial^{2} \tilde{w}_{2}}{\partial X^2}+\frac{\partial^{2} \tilde{w}_{0}}{\partial \tilde{y}^{2}}={-}2 \end{gather}$$

(2.20)$$\begin{gather}O\left(\epsilon^{n}\right):\quad \frac{\partial^{2} \tilde{w}_{n+2}}{\partial X^2}+\frac{\partial^{2} \tilde{w}_{n}}{\partial \tilde{y}^{2}}=0,\quad n\geqslant 1. \end{gather}$$

The symmetry conditions at $X=0,1$ apply to all orders of $\tilde {w}$. Integrating the $O(\epsilon ^{-1})$ problem and applying them shows that $\tilde {w}_{1}$, like $\tilde {w}_{0}$, is purely a function of $\tilde {y}$. The $O(\epsilon ^{0})$ problem shows that $\partial ^{2}\tilde {w}_{2}/\partial X^{2}$ is purely a function of $\tilde {y}$ and, on account of the symmetry conditions, so is $\tilde {w}_{2}$. By induction in $n$, the remaining problems show that $\tilde {w}_{n}=\tilde {w}{}_{n}(\tilde {y})$ for $n\geqslant 0$. Thus, (2.19)–(2.20) collapse to

(2.21)$$\begin{gather} \frac{\mathrm{d}^{2}\tilde{w}_{0}}{\mathrm{d}\tilde{y}^{2}}={-}2 \end{gather}$$

(2.22)$$\begin{gather}\frac{\mathrm{d}^{2}\tilde{w}_{n}}{\mathrm{d}\tilde{y}^{2}}=0,\quad n\geqslant 1. \end{gather}$$

Integrating and applying the symmetry condition at the channel centreline as per (2.14) yields

(2.23)$$\begin{gather} \tilde{w}_{0}\left(\hat{y}\right)={-}\tilde{y}^{2}+2\tilde{y}+C_{0}, \end{gather}$$

(2.24)$$\begin{gather}\tilde{w}_{n}\left(\tilde{y}\right)=C_{n}, \end{gather}$$

where $C_{n}$ are constants. Thus the solution to the outer problem is

(2.25)\begin{equation} \tilde{w}\sim{-}\tilde{y}^{2}+2\tilde{y}+C\left(\epsilon\right)+\mathrm{e.s.t.}, \end{equation}

where

(2.26)\begin{equation} C\left(\epsilon\right)=\sum_{n=0}^{\infty}C_{n}\epsilon^{n}.\end{equation}

We turn to the inner problem to satisfy the mixed boundary condition at the composite interface, i.e. the meniscus and the solid–liquid interface.

2.1.2. Inner region

The inner region of the domain is where $\tilde {x} \sim \tilde {y} =O(\epsilon )$, i.e. they are on the scale of the pitch of the ridges, as $\epsilon \rightarrow 0$. This is sufficiently close to the mixed boundary condition imposed at $\tilde {y}=0$ that the velocity field is two-dimensional, i.e. $\tilde {w} = \tilde {w}(\tilde {x},\tilde {y})$. Our stretched coordinates are $X=\tilde {x}/\epsilon$, $Y=\tilde {y}/\epsilon$ such that the inner region corresponds to $X\sim Y=O(1)$ as $\epsilon \rightarrow 0$ as per figure 2(a). Since the $X$ and $Y$ scales are of the same order, the viscous stress (Laplacian) terms are of the same order in the momentum equation. The velocity scale for the inner problem follows from the outer solution as $\tilde {y}$ is decreased to $O(\epsilon )$. Substituting $\tilde {y}=\epsilon Y$ into (2.25) yields

(2.27)\begin{equation} \tilde{w}\sim{-}\epsilon^{2}Y^{2}+2\epsilon Y+C\left(\epsilon\right)+\mathrm{e.s.t.} \end{equation}

Hence, when $C(\epsilon )=O(1)$, i.e. $C_{0}\neq 0$ as required by (2.11), $\tilde {w}={O(1)}$ in the inner region.

Expressing the inner velocity field as $\tilde {w}=W(X,Y)$, the problem becomes

(2.28)\begin{equation} \frac{\partial^{2}W}{\partial X^{2}}+\frac{\partial^{2}W}{\partial Y^{2}}={-}2\epsilon^{2},\end{equation}

subject to

(2.29)$$\begin{gather} \frac{\partial W}{\partial Y}=0\quad \mathrm{at}\ Y=0\quad \mathrm{for}\ 0< X<\delta \end{gather}$$

(2.30)$$\begin{gather}W=0\quad \mathrm{at}\ Y=0\quad \mathrm{for}\ \delta< X<1 \end{gather}$$

(2.31)$$\begin{gather}W\sim{-}\epsilon^{2}Y^{2}+2\epsilon Y+C(\epsilon)+\mathrm{e.s.t.}\quad \mathrm{as}\ Y\rightarrow\infty\ \mathrm{for}\ 0< X<1 \end{gather}$$

(2.32)$$\begin{gather}\frac{\partial W}{\partial X}=0\quad \mathrm{at}\ X=0\quad \mathrm{and}\quad X=1\quad \mathrm{for}\ Y>0, \end{gather}$$

where (2.31) is the (Van Dyke) matching condition. Note that the condition is written as an infinite series here (since $C(\epsilon )$ is (2.26)) and the inner solution $W$ has not yet been expanded in $\epsilon$, but the matching can be done up to a finite number of terms in the usual way – see Appendix A. Expressing $W$ as

(2.33)\begin{equation} W={-}\epsilon^{2}Y^{2}+2\epsilon\hat{W},\end{equation}

$\hat {W}$ must satisfy

(2.34)\begin{equation} \nabla^{2}\hat{W}=0, \end{equation}

subject to

(2.35)$$\begin{gather} \frac{\partial\hat{W}}{\partial Y}=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ 0< X<\delta \end{gather}$$

(2.36)$$\begin{gather}\hat{W}=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ \delta< X<1 \end{gather}$$

(2.37)$$\begin{gather}\hat{W}\sim Y+\frac{C\left(\epsilon\right)}{2\epsilon}+\mathrm{e.s.t.}\quad\mathrm{as}\ Y\rightarrow\infty \end{gather}$$

(2.38)$$\begin{gather}\frac{\partial\hat{W}}{\partial X}=0\quad\mathrm{at}\ X=0\quad\mathrm{and}\quad X=1\quad\mathrm{for}\ 0< Y<\infty. \end{gather}$$

This $\hat {W}$ problem is the superposition of a one-dimensional linear-shear flow over a smooth surface and a perturbation to it which accommodates a mixed boundary condition and manifests itself with a constant mean velocity over the width of the domain. It has been solved using a conformal map by Philip (Reference Philip1972a,Reference Philipb). The solution up to the exponentially small error in the matching condition is

(2.39)\begin{equation} \hat{W}={\rm Im}\left\{ \frac{2}{\rm \pi}\cos^{{-}1}\left[\frac{\mathrm{cos}\left({\rm \pi} \varTheta_{{\parallel}}/2\right)}{\cos\left({\rm \pi}\delta/2\right)}\right]\right\}+\mathrm{e.s.t.}, \end{equation}

where $\varTheta _{\parallel }=X+\mathrm {i}Y$, and

(2.40)\begin{equation} \hat{W}\sim Y+\lambda+\mathrm{e.s.t.}\quad\mathrm{as}\ Y\rightarrow\infty, \end{equation}

where

(2.41)\begin{equation} \lambda=\frac{2}{\rm \pi}\ln[\mathrm{sec}\left(\delta{\rm \pi}/2\right)].\end{equation}

It follows that

(2.42)\begin{equation} C\left(\epsilon\right)=2\epsilon \lambda\end{equation}

and the inner and outer solutions are

(2.43)\begin{gather} W\sim\overbrace{2\,{\rm Im}\left\{ \frac{2}{\rm \pi}\cos^{{-}1}\left[\frac{\cos\left({\rm \pi} \varTheta_{{\parallel}}/2\right)}{\mathrm{cos}\left({\rm \pi}\delta/2\right)}\right]\right\}}^{W_1}\epsilon+\overbrace{({-}Y^{2})}^{W_2}\epsilon^{2}+\mathrm{e.s.t.} \end{gather}

(2.44)\begin{gather}\tilde{w}\sim \overbrace{-\tilde{y}^{2} +2\tilde{y}}^{\tilde{w}_0}+\overbrace{2\lambda}^{\tilde{w}_1}\epsilon+\mathrm{e.s.t.}, \end{gather}

respectively. Note that we did not have to expand $W$ (or $\hat {W}$) in powers of $\epsilon$ to arrive at the above inner solution and it has been checked that matching order by order (as in Appendix A) gives the same result.

2.1.3. Composite solution

The solutions for the inner and outer regions are in agreement in the overlap region, where the outer one keeps its form. Therefore, the inner solution, (2.43), which is accurate to all algebraic orders in $\epsilon$, is uniformly valid throughout the domain, i.e. it is the composite solution as per, in outer variables,

(2.45)\begin{equation} \tilde{w}_{comp}\sim{-}\tilde{y}^2 +2\epsilon\,{\rm Im}\left\{\frac{2}{\rm \pi}\cos^{{-}1}\left[\frac{\cos\left({\rm \pi} \left(\tilde{x}/\epsilon + \mathrm{i}\tilde{y}/\epsilon\right)/2\right)}{\cos\left({\rm \pi}\delta/2\right)}\right]\right\}+\mathrm{e.s.t}. \end{equation}

2.1.4. Slip length

The slip length quantifies the flow enhancement provided by texturing a microchannel with a SH surface(s). It is related to the one-dimensional (Navier-slip) problem, which does not resolve the local velocity field in the inner region, and is governed by

(2.46)\begin{equation} \frac{\mathrm{d}^2 \tilde{w}_{{1d}}}{\mathrm{d} \tilde{y}^2}={-}2, \end{equation}

where ${\tilde {w}_{1{d}}}$ is the dimensionless velocity averaged over a ridge period. It is subject to the symmetry condition at $\tilde {y} = 1$ and a Navier-slip boundary condition on the SH surface as per

(2.47)\begin{equation} \tilde{w}_{1{d},\tilde{y}=0}=\tilde{b}\left.\frac{\mathrm{d}\tilde{w}_{{1d}}}{\mathrm{d} \tilde{y}}\right|_{\tilde{y}=0}, \end{equation}

where

(2.48)\begin{equation} \tilde{b}=\frac{b}{h}, \end{equation}

is the dimensionless slip length. It follows that

(2.49)\begin{equation} \tilde{w}_{{1d}}={-}\tilde{y}^2+2\tilde{y}+2\tilde{b}, \end{equation}

and its mean value is

(2.50)\begin{equation} \bar{\tilde{w}}_{{1d}}=2\tilde{b} +\tfrac{2}{3}. \end{equation}

Following Lauga & Stone (Reference Lauga and Stone2003), we compute the slip length by equating the mean velocities from the composite and one-dimensional (Navier-slip) velocity profiles. As with the flow configurations resolved by Philip (Reference Philip1972a,Reference Philipb), the velocity field for our configuration is the sum of that for a one-dimensional, Poiseuille flow in a microchannel with smooth walls and a velocity-increment field. The latter is not subject to any forcing term and obeys the shear-free boundary condition at $\tilde {y}=1$. Therefore, there is no net momentum transferred across any plane normal to the composite interface in the velocity-increment problem and thus the mean velocity increment over the width of the domain is a constant, i.e. $2\epsilon \lambda +\mathrm {e.s.t.}$ Thus, the outer solution is valid throughout the domain insofar as computing the slip length and it follows from (2.25), (2.42) and (2.49) that

(2.51)\begin{equation} \tilde{b}\sim \epsilon \lambda+\mathrm{e.s.t}. \end{equation}

2.2.1. Formulation

The dimensional form of the thermal energy equation is

(2.53)\begin{equation} w\frac{\partial T}{\partial z}=\alpha\left(\frac{\partial^{2}T}{\partial x^{2}}+\frac{\partial^{2}T}{\partial y^{2}}\right),\end{equation}

where axial conduction is absent because the temperature field is fully developed and $\alpha =k/(\rho c_{p})$ is the thermal diffusivity of the liquid, where $k$ is its thermal conductivity, $\rho$ is its density and $c_{p}$ is its specific heat at constant pressure. The boundary conditions on (2.53) are

(2.54)$$\begin{gather} \frac{\partial{T}}{\partial{y}}=0\quad\mathrm{at}\ {y}=0\quad\mathrm{for}\ 0<{x}< a \end{gather}$$

(2.55)$$\begin{gather}\frac{\partial{T}}{\partial{y}}={-}\frac{q''_{sl}}{k}\quad\mathrm{at}\ {y}=0\quad\mathrm{for}\ a<{x}< d \end{gather}$$

(2.56)$$\begin{gather}\frac{\partial{T}}{\partial{y}}=0\quad\mathrm{at}\ {y}=H\quad\mathrm{for}\ 0<{x}< d \end{gather}$$

(2.57)$$\begin{gather}\frac{\partial{T}}{\partial{x}}=0\quad\mathrm{at}\ {x}=0\quad\mathrm{and}\quad{x}=a\quad\mathrm{for}\ 0<{y}< H, \end{gather}$$

where $q''_{sl}$ is the (constant) heat flux prescribed along the solid–liquid interface. Moreover, again, because of the thermally fully developed assumption, along with the constant heat flux boundary condition, we have

(2.58)\begin{equation} \frac{\partial T}{\partial z}=\frac{q_{sl}''(d-a)}{Q\rho c_{p}}, \end{equation}

where $Q$ is the volumetric flow rate of liquid through a half-period.

Henceforth, we proceed with the non-dimensional form of the problem, where temperature is non-dimensionalized (in a similar way as was done by Kirk et al. Reference Kirk, Hodes and Papageorgiou2017) by subtracting the (bulk) mean liquid temperature, $T_{m}$, and dividing by a characteristic temperature scale of $q_{sl}''H/k$ such that

(2.59)\begin{equation} \tilde{T}=\frac{k\left(T-T_{m}\right)}{q_{sl}''H}, \end{equation}

where

(2.60)\begin{equation} T_{m}=\frac{\displaystyle\int\nolimits_0^H \int_0^dwT\,\mathrm{d}\,x\,\mathrm{d}y}{\displaystyle\int\nolimits_0^H \int_0^dw\,\mathrm{d}\,x\,\mathrm{d}y}.\end{equation}

The thermal energy equation becomes

(2.61)\begin{equation} \frac{\epsilon \phi}{\tilde{Q}}\tilde{w}=\frac{\partial^{2}\tilde{T}}{\partial\tilde{x}^{2}}+\frac{\partial^{2}\tilde{T}}{\partial\tilde{y}^{2}}, \end{equation}

where $\phi =1-\delta$ is the solid fraction of the ridges and $\tilde {Q} = 2\mu Q/[(-\mathrm {d}p/\mathrm {d}z) H^4]$. The forcing term, $\epsilon \phi \tilde {w}/\tilde {Q}$, is known (to all algebraic orders) from the solution to the hydrodynamic problem.

Noting that

(2.62)\begin{equation} \tilde{Q}\sim \left(2\epsilon \lambda +2/3\right)\epsilon+\mathrm{e.s.t.},\quad\epsilon \rightarrow 0, \end{equation}

the dimensionless thermal energy equation (written in terms of $X$ such that the domain boundaries are independent of $\epsilon$) becomes

(2.63)\begin{equation} \frac{\phi}{2/3+2\epsilon \lambda}\tilde{w}+\mathrm{e.s.t.}=\frac{1}{\epsilon^{2}}\frac{\partial^{2}\tilde{T}}{\partial X^{2}}+\frac{\partial^{2}\tilde{T}}{\partial\tilde{y}^{2}},\end{equation}

subject to

(2.64)$$\begin{gather} \frac{\partial{\tilde{T}}}{\partial{\tilde{y}}}=0\quad\mathrm{at}\ \tilde{y}=0\quad\mathrm{for}\ 0<{X}<\delta \end{gather}$$

(2.65)$$\begin{gather}\frac{\partial{T}}{\partial{\tilde{y}}}={-}1\quad\mathrm{at}\ \tilde{y}=0\quad\mathrm{for}\ \delta<{X}< 1 \end{gather}$$

(2.66)$$\begin{gather}\frac{\partial\tilde{T}}{\partial\tilde{y}}=0\quad\mathrm{at}\ \tilde{y}=1\quad\mathrm{for}\ 0<{X}<1 \end{gather}$$

(2.67)$$\begin{gather}\frac{\partial\tilde{T}}{\partial{X}}=0\quad\mathrm{at}\ {X}=0\quad\mathrm{and}\quad X=1. \end{gather}$$

2.2.2. Outer region

The forcing term in the above problem is $O(1)$ and thus, by the same reasoning used to justify that $\tilde {w} =\tilde {w}(\tilde {y}) = O(1)$, we also have $\tilde {T}=\tilde {T}(\tilde {y}) = O(1)$ to leading order. We thus expand the temperature field as

(2.68)\begin{equation} \tilde{T}=\sum_{n=0}^{\infty}\epsilon^{n}\tilde{T}_{n}+\mathrm{e.s.t.},\quad\epsilon \rightarrow 0,\end{equation}

where $\tilde {T}_{n}=O(1)$ for all $n\geqslant 0$.

Substituting (2.68) and the solution for the outer velocity profile as per (2.44) into (2.63), we find that, at the various orders of $\epsilon$,

(2.69)$$\begin{gather} O(\epsilon^{{-}2}):\quad\frac{\partial^{2} \tilde{T}_{0}}{\partial X^2}=0 \end{gather}$$

(2.70)$$\begin{gather}O(\epsilon^{{-}1}):\quad\frac{\partial^{2} \tilde{T}_{1}}{\partial X^2}=0 \end{gather}$$

(2.71)$$\begin{gather}O(\epsilon^{0}):\quad\frac{\partial^{2} \tilde{T}_{2}}{\partial X^2}+ \frac{\partial^{2} \tilde{T}_{0}}{\partial \tilde{y}^{2}}=\frac{3\phi}{2}\left(-\tilde{y}^{2}+2\tilde{y}\right) \end{gather}$$

(2.72)$$\begin{gather}O(\epsilon^{n}):\quad\frac{\partial^{2} \tilde{T}_{n+2}}{\partial X^2}+ \frac{\partial^{2} \tilde{T}_{n}}{\partial \tilde{y}^{2}}=\frac{3\phi}{2}[ ({-}3\lambda)^n(-\tilde{y}^2+2\tilde{y})+({-}3\lambda)^{n-1}2\lambda],\quad n\geqslant 1. \end{gather}$$

Following the same approach as in the outer hydrodynamic problem shows that $\tilde {T}_{n}=\tilde {T}{}_{n}(\tilde {y})$ for all $n$. The solution to the outer problem, which satisfies all 3 symmetry conditions, follows as

(2.73)\begin{equation} \tilde{T}=\frac{\phi}{2/3+2\epsilon \lambda}\left[-\frac{\left(\tilde{y}-1\right)^{4}}{12}+\left(\frac{1}{2}+\epsilon \lambda\right)\left(\tilde{y}-1\right)^{2}+D\left(\epsilon\right)\right] +\mathrm{e.s.t.},\quad\epsilon \rightarrow 0, \end{equation}

where

(2.74)\begin{equation} D(\epsilon)=\sum_{n=0}^{\infty}\epsilon^{n}D_{n}. \end{equation}

2.2.3. Inner region

Denoting the inner temperature profile as $\tilde {T}=\theta (X,Y)$, the thermal energy equation becomes

(2.75)\begin{align} \frac{\partial^{2}\theta}{\partial X^{2}}+\frac{\partial^{2}\theta}{\partial Y^{2}}=\frac{\phi}{2/3+2\epsilon \lambda}\left(-\epsilon^{4}Y^2+2\epsilon^{3}{\rm Im}\left\{ \frac{2}{\rm \pi}\cos^{{-}1}\left[\frac{\cos\left({\rm \pi} Z/2\right)}{\cos\left({\rm \pi}\delta/2\right)}\right]\right\} \right)+\mathrm{e.s.t.}\end{align}

The forcing is $O(\epsilon ^{3})$; therefore, neglecting terms of this order, it reduces to Laplace's equation as per

(2.76)\begin{equation} \frac{\partial^{2}\theta}{\partial X^{2}}+\frac{\partial^{2}\theta}{\partial Y^{2}}= O\left(\epsilon^3\right), \end{equation}

subject to

(2.77)\begin{gather} \frac{\partial\theta}{\partial Y}=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ 0< X<\delta \end{gather}

(2.78)\begin{gather}\frac{\partial\theta}{\partial Y}={-}\epsilon\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ \delta< X<1 \end{gather}

(2.79)\begin{gather}\theta\sim{-}\epsilon\phi Y+\phi\left[\frac{1}{2}+\frac{1/12+D(\epsilon)}{2/3+2\epsilon \lambda}\right]+O \left(\epsilon^{3}\right)\quad\mathrm{as}\ Y\rightarrow\infty \end{gather}

(2.80)\begin{gather}\frac{\partial\theta}{\partial X}=0\quad\mathrm{at}\ X=0\quad\mathrm{and}\quad X=1\quad\mathrm{for}\ Y>0, \end{gather}

where (2.79), which follows from (2.73) with $\tilde {y} = \epsilon Y$, is the matching condition. As per the solution by Mikic (Reference Mikic1957) to this problem in the context of (thermal) spreading resistance

(2.81)\begin{align} \theta ={-}\epsilon \phi Y - \frac{2\epsilon}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} X\right)\mathrm{e}^{{-}n {\rm \pi}Y}}{n^2}+\phi\left[\frac{1}{2}+\frac{1/12+D(\epsilon)}{2/3+2\epsilon \lambda}\right]+O\left(\epsilon^{3}\right). \end{align}

We solve for $D(\epsilon )$ by enforcing that the dimensionless mean liquid temperature, as defined by (2.59) and (2.60), is zero, i.e. using outer variables

(2.82)\begin{equation} \int_{0}^{1}\int_0^{\epsilon}\tilde{w}\tilde{T}\,\mathrm{d}\,\tilde{x}\,\mathrm{d}\tilde{y}=0. \end{equation}

In Appendix B, we show that we can enforce this condition using the outer solutions only, up to an error of $O(\epsilon ^3)$, and the final temperature along the composite interface is found to be

(2.83)\begin{equation} \theta \left(X,0\right) = \phi \hat{D}\left(\epsilon \lambda\right) - \frac{2\epsilon}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} X\right)}{n^2}+O\left(\epsilon^{3}\right), \end{equation}

where $\hat {D}(\epsilon \lambda )$ (which is related to $D(\epsilon )$) is given by the explicit formula (B9).

2.2.4. Composite solution

The inner solution is not the composite one as in the hydrodynamic problem because (2.79) is consistent with (2.73) only up to $O(\epsilon ^2)$. Rather

(2.84)\begin{equation} \tilde{T}_{comp} = \tilde{T}_{out} + \tilde{T}_{in} - \tilde{T}_{overlap}; \end{equation}

consequently

(2.85)\begin{align} \tilde{T}_{comp} &\sim \frac{\phi}{\dfrac{2}{3}+2\epsilon \lambda}\left[-\frac{\left( \tilde{y}-1\right)^{4}}{12}+\left(\frac{1}{2}+\epsilon \lambda\right)\left(\tilde{y}-1\right)^{2}+ \breve{D}\left(\epsilon\right)\right]\nonumber\\ &\quad - \frac{2\epsilon}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} \dfrac{\tilde{x}}{\epsilon}\right)\exp\left({-n {\rm \pi}\dfrac{\tilde{y}}{\epsilon}}\right)} {n^2}+O\left(\epsilon^3\right). \end{align}

2.2.5. Nusselt numbers

We define the local Nusselt number along the solid–liquid interface as

(2.86)\begin{equation} {Nu}(X) = \frac{4h(X)H}{k}, \end{equation}

where $h(X)$ is the (local) heat transfer coefficient, i.e. $q''_{sl}/(T_{sl}-T_{m})$, which is finite along the solid–liquid interface and zero along the meniscus. Then

(2.87)\begin{equation} {Nu}(X) = \begin{cases} 0 & \mathrm{for}\ 0\leqslant X < \delta\\ \dfrac{4}{\theta\left(X,0\right)} & \mathrm{for}\ \delta < X \leqslant 1. \end{cases} \end{equation}

The minimum local Nusselt number (${Nu_{min}}$) occurs at the centre of a ridge, i.e. $X=1$, as per

(2.88)\begin{equation} {Nu_{min}} = \frac{4}{\theta\left(1,0\right)}. \end{equation}

This is an important engineering parameter because the magnitude of the maximum temperature difference between the ridge and bulk liquid follows from it. From (2.83), the minimum local Nusselt number (${Nu_{min}}$) may be expressed as

(2.89)\begin{equation} {Nu_{min}} = \frac{4}{ \phi \left[ \hat{D}\left(\epsilon \lambda\right)+\dfrac{2\epsilon}{\phi} \sum_{n=1}^{\infty}\dfrac{\sin\left(n{\rm \pi}\phi\right)}{\left(n{\rm \pi}\right)^2}\right]}+O\left(\epsilon^{3}\right). \end{equation}

More generally, the local Nusselt number anywhere along the ridge is

(2.90)\begin{align} {Nu}(X) = \frac{4}{ \phi \left( \hat{D}\left(\epsilon \lambda\right) - \dfrac{2\epsilon}{\phi} \sum_{n=1}^{\infty}\dfrac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} X\right)}{n^2{\rm \pi}^2}\right) }+O\left(\epsilon^{3}\right), \quad \delta < X \leqslant 1. \end{align}

This was derived in the limit of small $\epsilon$, but the subsequent behaviour in the small solid fraction limit, $\phi \to 0$, is also of interest. This limit is physically relevant, but care should be taken when making further approximations in $\epsilon$ that the correct asymptotic behaviour in the secondary limit $\phi \to 0$ is preserved. If it is not preserved, e.g. if one simply expands (2.90) directly in powers of $\epsilon$, the accuracy and validity of such an expression is significantly reduced. This is elucidated in Appendix C, where we instead expand in powers of $\mathcal {E}(\epsilon )= \epsilon / (\hat {D} + \epsilon \lambda )$, yielding

(2.91)\begin{align} {Nu}(X) &= \frac{4}{\phi \left[\hat{D}\left(\epsilon \lambda\right)+\epsilon \lambda\right]} \left\{1+ \mathcal{E}(\epsilon)\left[ \lambda+S(X)\right] + \mathcal{E}(\epsilon)^2\left[\lambda+S(X)\right]^2\right\} \nonumber\\ &\quad + O\left(\epsilon^3\right)\quad \mathrm{for}\ \delta < X \leqslant 1. \end{align}

The quantity $\mathcal {E}$ remains bounded even if $\phi \to 0$, and hence the expansion stays well ordered. The $O(\epsilon ^3)$ error term accounts for the $O(\mathcal {E}^3)$ terms since $\mathcal {E} = O(\epsilon )$. We note that, by evaluating Nu in the limit when the solid fraction is 1, i.e. for a smooth channel such that $\lambda \rightarrow 0$, we recover the well-known result that $Nu = 140/17$. (The $O(\epsilon ^3)$ error term can be shown to vanish in this limit.)

The average Nusselt number for the composite interface is defined in the manner of Maynes & Crockett (Reference Maynes and Crockett2014) and Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) as

(2.92)\begin{equation} \overline{Nu} = \frac{4\bar{h}H}{k}, \end{equation}

where $\bar {h} = \int _a^d h\,\mathrm {d}\,x/d$. Thus

(2.93)\begin{equation} \overline{Nu}=\int_{\delta}^{1}{Nu}(X)\,\mathrm{d}X. \end{equation}

Performing the integration on (2.91), and substituting the series $S(X)$, we find that

(2.94)\begin{align} \overline{Nu} &\sim \frac{4}{\hat{D}\left(\epsilon \lambda\right)+\epsilon \lambda}\left\{ 1+\left[\lambda-\sum_{n=1}^{\infty}\frac{2\sin^2\left(n{\rm \pi}\delta\right)}{\phi^2\left(n{\rm \pi}\right)^3} \right]\mathcal{E}(\epsilon)\right. \nonumber\\ &\quad \left.+ \left[\lambda^2-2\lambda\sum_{n=1}^{\infty}\frac{2\sin^2\left(n{\rm \pi}\delta\right)}{\phi^2\left(n{\rm \pi}\right)^3}+\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}L_{m,n}\right]\mathcal{E}(\epsilon)^2\right\} +O\left(\epsilon^3\right), \end{align}

where

(2.95)\begin{align} L_{mn} = \frac{4}{\phi^3}\frac{\sin\left(n{\rm \pi}\delta\right) \sin\left(m{\rm \pi}\delta\right)}{\left(n{\rm \pi}\right)^2\left(m{\rm \pi}\right)^2} \begin{cases} \dfrac{1}{2}\left[1-\delta - \dfrac{\sin\left(2n{\rm \pi}\delta\right)}{2n{\rm \pi}}\right] & \mathrm{for}\ m=n\\[10pt] -\dfrac{\sin\left[\left(m+n\right){\rm \pi}\delta\right]}{2{\rm \pi}\left(m+n\right)} - \dfrac{\sin\left[\left(m-n\right){\rm \pi}\delta\right]}{2{\rm \pi}\left(m-n\right)} & \mathrm{for}\ m \neq n. \end{cases} \end{align}

Alternatively, a simpler expression valid for any solid fraction with no sums to evaluate, but with an $O(\epsilon )$ error term, follows from the leading term of (2.94)

(2.96)\begin{equation} \overline{Nu} \sim \frac{4}{\hat{D}\left(\epsilon \lambda\right)+\epsilon \lambda} + O\left(\epsilon\right). \end{equation}

An average Nusselt number may also be defined in the manner of Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) as

(2.97)\begin{equation} \overline{Nu}' =\frac{4\bar{h}'H}{k}, \end{equation}

by defining the average heat transfer coefficient as $\bar {h}'= \phi q''_{sl}/(\overline {T_{sl}-T_{m}})$ such that it is based upon the mean heat flux over the composite interface ($\phi q''_{sl}$) and the mean driving force for heat transfer over the solid–liquid one ($\overline {T_{sl}-T_{m}}$). Then, it is given by

(2.98)\begin{equation} \overline{Nu}' = \frac{4\phi^2}{\displaystyle\int\nolimits_{\delta}^1 \theta \left(X,0\right) \mathrm{d}X}. \end{equation}

Performing the integration we find that

(2.99)\begin{align} \overline{Nu}'&= 140\left(1+3\epsilon \lambda\right)^2\left/\left\{ 17+ \left[84\lambda + \frac{70}{\phi^2{\rm \pi}^3}\sum_{n=1}^{\infty}\frac{\sin^2\left(n{\rm \pi}\delta\right)}{n^3}\right]\epsilon \right. \right. \nonumber\\ &\quad \left. +\left[ 105\lambda^2 +\frac{420\lambda}{\phi^2{\rm \pi}^3}\sum_{n=1}^{\infty}\frac{\sin^2\left(n{\rm \pi}\delta\right)}{n^3} \right]\epsilon^2+\left[ \frac{630\lambda^2}{\phi^2{\rm \pi}^3}\sum_{n=1}^{\infty}\frac{\sin^2\left(n{\rm \pi}\delta\right)}{n^3}\right]\epsilon^3\right\} +O\left(\epsilon^3\right). \end{align}

There is no need to expand this result further in $\epsilon$, and it also vanishes as expected when $\lambda \rightarrow \infty$ (i.e. $\phi \rightarrow 0$). Also, the last term in the denominator is retained as, depending on the relative values of $\lambda$ and $\epsilon$, it may contribute significantly to the result. Finally, we note that, regardless of which definition of the mean Nusselt number is used, the local one is given by (2.86). In Appendix D, we validate the asymptotic behaviour of this result against the exact solution of Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017).

2.2.6. Prior results

The prior studies by Maynes & Crockett (Reference Maynes and Crockett2014) and Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) provide the local Nusselt number as per (2.86) and, by implication, the minimum one as per (2.88), as well as the mean Nusselt Number $\overline {Nu}$ as per (2.92). The mean Nusselt number $\overline {Nu}'$ as per (2.97) also follows from the temperature profiles developed in these studies. Moreover, the studies by Enright et al. (Reference Enright, Eason, Dalton, Hodes, Salamon, Kolodner and Krupenkin2006) and Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) provide $\overline {Nu}'$ and machinery that may be adopted to find the local Nusselt numbers. We proceed to discuss and contrast all of them and our own results.

Maynes & Crockett (Reference Maynes and Crockett2014) used a Navier-slip velocity profile in the thermal energy equation. An analytical solution for the temperature profile was found by representing the discontinuous (Neumann) boundary condition along the composite interface as a Fourier series and using separation of variables. Their temperature profile (in our notation) along the composite interface is

(2.100)\begin{align} \theta_{M}\left(X,0\right) &= \frac{17}{35}\phi -\frac{2\epsilon}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\delta\right)\coth\left(n{\rm \pi}/\epsilon\right) \cos\left(n{\rm \pi} X\right)}{n^2} -\frac{18}{35}\frac{\phi \epsilon \lambda}{1+3\epsilon \lambda}\nonumber\\ &\quad +\frac{6}{35}\frac{\phi\epsilon^2 \lambda^2}{\left(1+3\epsilon \lambda\right)^2}, \end{align}

where the subscript ‘M’ denotes a result by Maynes & Crockett (Reference Maynes and Crockett2014). Closed-form expressions for ${Nu}_{M}$ and ${Nu}_{{min,M}}$ follow and their behaviour is discussed by Maynes & Crockett (Reference Maynes and Crockett2014). A closed-form expression for $\overline {Nu}'_{M}$ also follows from integrating (2.100). Numerical integration was used by Maynes & Crockett (Reference Maynes and Crockett2014) to compute $\overline {Nu}_{M}$. Although Maynes & Crockett (Reference Maynes and Crockett2014) did not quantify the error in their Nusselt number expressions due to using a Navier-slip velocity profile, we have shown that it is $O(\epsilon ^3)$. Relatedly, upon noting that $\coth (n{\rm \pi} /\epsilon ) \sim 1 + \mathrm {e.s.t.}$ as $\epsilon \rightarrow 0$, an expansion of (2.100) is consistent with (2.83) and all of the Nusselt numbers become the same as in our analysis. We note, however, that, since Maynes & Crockett (Reference Maynes and Crockett2014) did not do this, their results differ slightly from ours.

Insofar as $\overline {Nu}'$, Enright et al. (Reference Enright, Eason, Dalton, Hodes, Salamon, Kolodner and Krupenkin2006) were the first to compute it and, as here, solid–liquid interfaces were isoflux. They too utilized a Navier-slip velocity profile as per (2.50), or the corresponding expression when one side of the microchannel is textured, but did not capture the discontinuous thermal boundary condition along the composite interface. Consequently, the Nusselt number was solely dependent upon the slip length non-dimensionalized by the channel height for the texture of interest (parallel or transverse ridges, pillars, etc.). Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) refined this approach by further imposing an apparent temperature jump along the composite (c) interfaces on one or both sides of a microchannel as per

(2.101)\begin{equation} \bar{T}_{sl} - \bar{T}_{c} ={-}b_{t} \left.\frac{\partial \bar{T}}{\partial n}\right|_{c}, \end{equation}

where $b_{t}$ is the apparent thermal slip length (also referred to as the temperature jump length) and $n$ the direction normal to a composite interface and pointing into the liquid. Using the approach of Nield (Reference Nield2004, Reference Nield2008) to accommodate asymmetrical boundary conditions, closed-form expressions for the Nusselt number were developed as a function of arbitrary values of $b$ and $b_{t}$ imposed on each side of the microchannel and, additionally, the ratio of heat fluxes averaged over the composite interfaces ($\phi q''_{sl}$) bounding the domain (liquid). In the case of a symmetric channel, as in this study, the results by Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) reduce to those developed by Inman (Reference Inman1964) 50 years earlier in the context of molecular slip effects in gas flows, where the surface boundary conditions are mathematically equivalent. Kane & Hodes (Reference Kane and Hodes2019) extended the Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) analysis to a combined Poiseuille and Couette flow.

The (one-dimensional) temperature field provided by the foregoing studies in the case of parallel ridges is that averaged over the width of the domain. Insofar as the definition of dimensionless temperature adopted here as per (2.59), only the composite interface temperature is relevant; consequently, the thermal slip length is irrelevant. The resulting temperature profile, which we denote by $\tilde {T}_{E}$, is

(2.102)\begin{align} \tilde{T}_{E} = \frac{\phi}{2/3+2\epsilon \lambda}\left( -\frac{\tilde{y}^4}{12} + \frac{\tilde{y}^3}{3} + \epsilon \lambda \tilde{y}^2 -2\epsilon \lambda \tilde{y} - \frac{2}{3} \tilde{y} +\frac{2}{105}\frac{17+84\epsilon \lambda + 105 \epsilon^2 \lambda^2}{1+3\epsilon \lambda} \right). \end{align}

This expression is identical to that given by (2.73) with $D(\epsilon )$ given by (B7). Terms of $O(\epsilon ^3)$ were neglected in our derivation and thus, albeit not pointed out by Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014), it is only accurate to $O(\epsilon ^3)$.

In the foregoing studies (except in Enright et al. Reference Enright, Eason, Dalton, Hodes, Salamon, Kolodner and Krupenkin2006), the difference between the mean temperatures of the solid–liquid and composite interfaces follows by imposing an apparent thermal slip boundary condition as per (2.101). Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) utilized the spreading resistance solution by Mikic (Reference Mikic1957) to evaluate $b_{t}$ such that

(2.103)\begin{equation} b_{t}=\frac{2d}{{\rm \pi}^3\phi^2} \sum_{n=1}^{\infty} \frac{\sin^2\left(n{\rm \pi}\phi\right)}{n^3}. \end{equation}

The Nusselt number $\overline {Nu}'_{E}$, where the subscript ‘E’ indicates a result from Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014), then follows from the definition of the bulk temperature and it does not capture the error term as in the present analysis. Consequently, the relationship between the present result and that by Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) is

(2.104)\begin{equation} \overline{Nu}'=\overline{Nu}'_{E} + O(\epsilon^3). \end{equation}

Replacing $b_{t}$ in the Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) analysis with its maximum value, $b_{{t,max}}$, to define the apparent temperature jump along the composite interface in terms of the temperature of the centre of the ridge, (2.103) becomes

(2.105)\begin{equation} T_{sl}\left(x=d\right) - \bar{T}_{c} ={-}b_{{t,max}} \left.\frac{\partial \bar{T}}{\partial n}\right|_{c}. \end{equation}

It follows from the Mikic (Reference Mikic1957) solution that

(2.106)\begin{equation} b_{{t,max}}=\frac{2d}{{\rm \pi}^2\phi} \sum_{n=1}^{\infty} \frac{\sin\left(n{\rm \pi}\phi\right)}{n^2}. \end{equation}

Then, the Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) approach yields

(2.107)\begin{align} {Nu}_{{min,E}} &= 140\left(1+3\epsilon \lambda\right)^2\left/\left\{ 17\phi+ \left[84\lambda\phi + \frac{70}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\phi\right)}{n^2}\right]\epsilon \right. \right. \nonumber\\ &\quad \left. +\left[ 105\lambda^2\phi +\frac{420\lambda}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\phi\right)}{n^2} \right]\epsilon^2 +\left[ \frac{630\lambda^2}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\phi\right)}{n^2}\right]\epsilon^3\right\}. \end{align}

Upon rearrangement, this result may be shown to be identical to our own result as per (2.89), except that it does not quantify the error. More generally, the thermal slip length in the Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) formulation may be based upon any location along the ridge and thus the variation of the local Nusselt number along it determined. Then, the local Nusselt number $\overline {Nu}$ follows. We compare the value of $\overline {Nu}'_{E}$ (or, equivalently, $\overline {Nu}'$ from this study) with its exact value as per the Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) study in figure 6. Of course, the agreement only breaks down at sufficiently large values of $\epsilon$.

Only the aforementioned exact results by Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) are valid for arbitrary $\epsilon$. Also, by taking the limit of the dual-series equations as $\epsilon \rightarrow 0$, and the same limit in the temperature problem, an expression for $\overline {Nu}$ with an error term of $O(\epsilon ^2)$ was found by Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017). It breaks down for sufficiently small solid fraction. We summarize the Nusselt number developed in this study in table 1.

Table 1. Nusselt number results for parallel ridges.

2.2.7. Comparisons

A plot of ${Nu}_{min}$ vs $\phi$ for selected values of $\epsilon$ based on the present result and the aforementioned studies is shown in figure 3. As expected, it asymptotes to $\infty$ and 140/17 as $\phi$ approaches 0 and 1, respectively. When $\epsilon$ gets sufficiently large, ${Nu}_{min}$ no longer monotonically decreases with $\phi$. This is because there is essentially no flow over the middle of the ridge such that ${Nu}_{min} \rightarrow 0$; however, the aforementioned behaviour of it as $\phi \rightarrow 0,1$ remains valid. When the domain is square ($\epsilon = 1$), our solution and the Maynes & Crockett (Reference Maynes and Crockett2014) one overpredict ${Nu_{min}}$ by maximum amounts of 21.6 % and 21.7 %, respectively, at $\phi$ = 0.64. When $\epsilon = 2$, a local minimum in ${Nu}_{min}$ occurs at approximately $\phi = 1/2$, presumably because the flow velocity is rather low above the centre of the ridge relative to above the triple contact line (and meniscus). Using only the exact solution from Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017), in figure 4, we show the same results up to $\epsilon = 16$. The value of ${Nu}_{min}$ becomes very low as $\epsilon$ is sufficiently increased, except as $\phi \rightarrow 0$ or $\phi \rightarrow 1$, because most of the flow in the domain is between opposing menisci rather than between opposing ridges.

Figure 3. Value of ${Nu}_{min}$ vs $\phi$ for selected values of $\epsilon$ from present result (or, equivalently, that which follows from the approach of Enright et al. Reference Enright, Hodes, Salamon and Muzychka2014) per (2.89) shown by dashed curves, Maynes & Crockett (Reference Maynes and Crockett2014) expression for composite interface temperature per (2.100) shown by dot-dashed curves and exact solution by Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) shown by solid curves.

Figure 4. Value of ${{Nu_{min}}}$ vs $\phi$ for selected values of relatively large $\epsilon$ calculated using the (exact) method of Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017).

A semilog plot of $\overline {Nu}$ vs $\phi$ for selected values of $\epsilon$ based on our results and those from the aforementioned studies is shown in figure 5. When $\epsilon = 0.5$, a substantive value, the maximum discrepancies between $\overline {Nu}$ from Maynes & Crockett (Reference Maynes and Crockett2014), the present result with an error term of $O(\epsilon ^3)$, (2.94), and that with an error term of $O(\epsilon )$, (2.96), compared with the exact solution are only 2.6 %, 1.2 % and 5.2 %, respectively. When $\epsilon$ is increased to 2, they become 36.6 %, 38.7 % and 46.5 %, respectively, occurring near a solid fraction of 0.4. An extensive discussion of the physics governing the behaviour of the curves in figure 5, except that given by (2.94), is given by Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017). The behaviour of the latter is analogous to that of the asymptotic result in Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017), except that the error is $O(\epsilon ^3)$ rather than $O(\epsilon ^2)$ and it does not break down at small $\phi$. Finally, figure 6 shows the comparison between the Enright et al. result (equivalently, (2.99)) and the exact solution, which is rather good up to $\epsilon = 2$.

Figure 5. Semilog plot of $\overline {Nu}$ vs $\phi$ for selected values of $\epsilon$ from present result per (2.94) shown by purple curves, simplified form of present result per (2.96) shown by cyan curves, numerically evaluated integral from Maynes & Crockett (Reference Maynes and Crockett2014) shown by red curves, asymptotic expression by Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) shown by green curves and exact solution by Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) shown by black curves.

Figure 6. Value of $\overline {Nu}'$ vs $\phi$ for $\epsilon$ = 1/20, 0.1, 0.5 and 2 from Enright et al. or, equivalently, the present study, (2.99), shown by purple curves and exact solution by Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) shown by black curves.

2.2.8. Other relevant studies

Beyond the problem considered here, Kirk et al. (Reference Kirk, Hodes and Papageorgiou2017) accounted for meniscus curvature, assuming a small deflection from flat, by using a boundary perturbation. Also, the numerical results by Game et al. (Reference Game, Hodes, Kirk and Papageorgiou2018) consider arbitrary meniscus curvature, so long as it is a circular arc. Finally, Karamanis et al. (Reference Karamanis, Hodes, Kirk and Papageorgiou2018) solved the extended Graetz–Nusselt problem in the presence of viscous dissipation and uniform volumetric heat generation for a flat meniscus assuming isothermal, parallel ridges. Nusselt number results for transverse ridges are provided below. Those for (square) pillars using the Enright et al. (Reference Enright, Hodes, Salamon and Muzychka2014) approach are compared with numerical results in that study. Finally, very recently, Sharma et al. (Reference Sharma, Gaddam, Agrawal, Joshi and Dimov2020) provided numerical results for square, triangular and herringbone pillars in regular and staggered arrays.

2.2.9. Alternate form of inner solution and corresponding spreading resistances

Since the heat flux is discontinuous along the composite interface ($Y = 0$), the series in (2.81) converges non-uniformly in the domain ($Y > 0$) and requires many terms to accurately represent the temperature and heat flux profiles near the composite interface. Therefore, we ‘sum the series’ to obtain the exact solution in closed form. The dimensionless temperature field is a harmonic function, which can be written in terms of the imaginary part of a complex potential as per

(2.108)\begin{equation} \theta={\rm{Im}}\left[f\left(\varTheta_{{\parallel}}\right)\right], \end{equation}

where, recall, $\varTheta _{\parallel } = X + \mathrm {i}Y$. Judicious use of Euler's formula (as in Bazant (Reference Bazant2004), § V.G.) shows that

(2.109)\begin{align} f\left(\varTheta_{{\parallel}}\right)={-}\epsilon \left\{\phi \varTheta_{{\parallel}}+\frac{{Li}_{2}\left(M_+\right)-{Li}_{2}\left(M_{-}\right)}{{\rm \pi}^{2}}\right\} + \mathrm{i}\phi\left[\frac{1}{2}+\frac{1/12+D(\epsilon)}{2/3+2\epsilon \lambda}\right]+O\left(\epsilon^{3}\right), \end{align}

where the polylogarithm (special) function is

(2.110)\begin{equation} {Li}_{s}\left(M\right)=\sum_{n=1}^{\infty}\frac{M^{n}}{n^{s}},\quad\left|M\right|<1 ,\end{equation}

and

(2.111)\begin{equation} M_{{\pm}}=\exp\left({\mathrm{i}{\rm \pi}\left(\varTheta_{{\parallel}}\pm\delta\right)}\right). \end{equation}

Appealingly, in that no special functions are required, the dimensionless heat flux vector follows as

(2.112)\begin{equation} -\boldsymbol{\nabla} \theta ={-}\left(\frac{\partial \theta}{\partial X} + \mathrm{i}\frac{\partial \theta}{\partial Y}\right) = \frac{1}{\mathrm{i}}\overline{f'\left(\varTheta_{{\parallel}}\right)}, \end{equation}

where, using the relation that $\mathrm {d}{Li}_2( M_{\pm })/\mathrm {d}M_{\pm } = -\ln (1-M_{\pm })/M_{\pm }$,

(2.113)\begin{equation} f'\left(\varTheta_{{\parallel}}\right) ={-}\epsilon \phi +\frac{\mathrm{i}\epsilon}{\rm \pi}\left[\ln\left(1-M_+\right)-\ln\left(1-M_-\right)\right]. \end{equation}

The preceding result is also useful in the context of spreading resistance, $R_{sp}$, which we define conventionally. Hence, it is the increase in the temperature difference driving heat transfer, i.e. the mean heat source temperature minus the far-field temperature, beyond that of the corresponding one-dimensional problem, divided by the heat rate. Thus, the (dimensional) spreading resistance based upon the heat rate per unit depth of the domain, i.e. in units of ${\rm W}\ ({\rm m}\ {\rm K})^{-1}$, is

(2.114)\begin{equation} R'_{sp} = \frac{\displaystyle\int\nolimits_a^d T\left(x,0\right)\mathrm{d}\,x/\left(d-a\right) - \lim\nolimits_{y\rightarrow \infty}T}{q''_{sl}\left(d-a\right)} - \frac{T_{{1d}}\left(0\right) - \lim\nolimits_{y\rightarrow \infty}T_{{1d}}}{q''_{sl}\left(d-a\right)}, \end{equation}

where $T_{{1d}}(y)$ is, to within an additive constant, $-q''_{sl}(d-a)y/(kd)$. Non-dimensionalizing the spreading resistance by $1/k$, the well-known result by Mikic (Reference Mikic1957) is

(2.115)\begin{equation} \tilde{R}'_{sp} = \frac{2}{{\rm \pi}^3\phi^2} \sum_{n=1}^{\infty} \frac{\sin^2\left(n{\rm \pi}\phi\right)}{n^3}, \end{equation}

or, albeit not previously realized, in closed form

(2.116)\begin{equation} \tilde{R}'_{sp} = \frac{1}{{\rm \pi}^3\phi^2}{\rm Re}\left[{Li}_3\left(1\right) - {Li}_3\left(\mathrm{e}^{\mathrm{i}2{\rm \pi}\phi}\right) \right]. \end{equation}

This result is an exact result as the spreading resistance problem is one of pure diffusion in a semi-infinite domain. It is twice the value reported by Mikic (Reference Mikic1957) because the width of his domain corresponds to one ridge pitch (i.e. the full width of the heat source) and ours is half of that. (The spreading resistance based on the far-field heat flux, $R''_{sp}$, does not depend on the domain width and, non-dimensionalizing it by $(d-a)/k$, it follows that $\tilde {R}''_{sp}= \tilde {R}'_{sp}$.) More conservatively, the corresponding dimensionless spreading resistance based on the maximum source temperature, $\tilde {R}'_{{sp,max}}$ (or, equivalently, $\tilde {R}''_{{sp,max}}$), is

(2.117)\begin{equation} \tilde{R}'_{{sp,max}} = \frac{2}{{\rm \pi}^2\phi}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\phi\right)}{n^2}, \end{equation}

or, in closed form,

(2.118)\begin{equation} \tilde{R}'_{{sp,max}} =\frac{2}{{\rm \pi}^2\phi}{\rm Im}\left[{Li}_2\left(\mathrm{e}^{{\rm i}{\rm \pi}\phi}\right)\right]. \end{equation}

In the context of the asymptotics performed here, if a constant heat flux is applied along the top of a finite rather than infinite height domain, the error in the spreading resistance is exponentially small as $\epsilon \rightarrow 0$. For example, (2.115) would become

(2.119)\begin{equation} \tilde{R}'_{{sp,a}} \sim \frac{2}{{\rm \pi}^3\phi^2} \sum_{n=1}^{\infty} \frac{\sin^2\left(n{\rm \pi}\phi\right)}{n^3} + \mathrm{e.s.t.},\quad\epsilon \rightarrow 0, \end{equation}

where the subscript ‘sp,a’ signifies this is an asymptotic limit. Clearly, the foregoing results may be used to eliminate some of the infinite summations in the results that we provided. Finally, we note that Hodes, Kirk & Crowdy (Reference Hodes, Kirk and Crowdy2018) discusses analogies between the slip length for a linear-shear flow and spreading and contact resistances.

3.1.1. Outer region

The outer region is where $Z=O(1)$ and $\tilde {y}=\mathrm {ord}(1)$ as $\epsilon \rightarrow 0$. Substituting (3.10)–(3.12) into (3.2)–(3.4) yields, at the algebraic orders of $\epsilon$,

(3.13)\begin{gather} O\left(\epsilon^{{-}2}\right):\ \frac{\partial^{2} \tilde{v}_{0}}{\partial Z^2}=0 \end{gather}

(3.14)\begin{gather}\frac{\partial^{2} \tilde{w}_{0}}{\partial Z^2}=0 \end{gather}

(3.15)\begin{gather} O\left(\epsilon^{{-}1}\right):\ \frac{\partial \tilde{w}_{0}}{\partial Z}=0 \end{gather}

(3.16)\begin{gather}{Re}\left(\tilde{w}_{0}\frac{\partial \tilde{v}_{0}}{\partial Z}\right)=\frac{\partial^{2}\tilde{v}_{1}}{\partial Z^{2}} \end{gather}

(3.17)\begin{gather}{Re}\left(\tilde{w}_{0}\frac{\partial \tilde{w}_{0}}{\partial Z}\right)=2\frac{\partial \tilde{p}_{{p,}0}}{ \partial Z}+\frac{\partial^{2}\tilde{w}_{1}}{\partial Z^{2}} \end{gather}

(3.18)\begin{gather} O\left(\epsilon^{0}\right): \frac{\partial \tilde{v}_{0}}{\partial\tilde{y}}+\frac{\partial \tilde{w}_{1}}{\partial Z}=0 \end{gather}

(3.19)\begin{gather}{Re}\left(\tilde{v}_{0}\frac{\partial\tilde{v}_{0}}{\partial\tilde{y}}+\tilde{w}_{0} \frac{\partial \tilde{v}_{1}}{\partial Z}+\tilde{w}_{1}\frac{\partial\tilde{v}_{0}}{\partial Z}\right)= 2\frac{\partial\tilde{p}_{{p,}0}}{\partial\tilde{y}}+\frac{\partial^{2}\tilde{v}_{0}}{\partial\tilde{y}^{2}}+ \frac{\partial^{2}\tilde{v}_{2}}{\partial Z^{2}} \end{gather}

(3.20)\begin{gather}{Re}\left(\tilde{v}_{0}\frac{\partial\tilde{w}_{0}}{\partial\tilde{y}}+\tilde{w}_{0} \frac{\partial \tilde{w}_{1}}{\partial Z}+\tilde{w}_{1}\frac{\partial\tilde{w}_{0}}{\partial Z}\right)= 2+2\frac{\partial\tilde{p}_{{p,}1}}{\partial Z}+\frac{\partial^{2}\tilde{w}_{0}}{\partial \tilde{y}^{2}}+\frac{\partial^{2}\tilde{w}_{2}}{\partial Z^{2}} \end{gather}

(3.21)\begin{gather} O\left(\epsilon^{n}\right): \frac{\partial\tilde{v}_{n}}{\partial\tilde{y}}+\frac{\partial\tilde{w}_{n+1}}{\partial Z}=0\quad \mathrm{for}\ n\geqslant 0 \end{gather}

(3.22)\begin{gather}{Re}\left(\sum_{m=0}^n \tilde{v}_m \frac{\partial \tilde{v}_{n-m}}{\partial \tilde{y}} + \sum_{m=0}^{n+1} \tilde{w}_m \frac{\partial \tilde{v}_{n+1-m}}{\partial Z} \right)=2\frac{\partial\tilde{p}_{{p,}n}}{\partial\tilde{y}}+\frac{\partial^{2}\tilde{v}_{n}} {\partial\tilde{y}^2}+\frac{\partial^{2}\tilde{v}_{n+2}}{\partial Z^{2}} \quad\mathrm{for}\ n\geqslant 1 \end{gather}

(3.23)\begin{gather}{Re}\left(\sum_{m=0}^n \tilde{v}_m \frac{\partial \tilde{w}_{n-m}}{\partial \tilde{y}} + \sum_{m=0}^{n+1} \tilde{w}_m \frac{\partial \tilde{w}_{n+1-m}}{\partial Z}\right)= 2\frac{\partial\tilde{p}_{{p,}n+1}}{\partial Z}+\frac{\partial^{2}\tilde{w}_{n}}{\partial\tilde{y}^2}+ \frac{\partial^{2}\tilde{w}_{n+2}}{\partial Z^{2}}\quad\mathrm{for} \ n\geqslant 1. \end{gather}

The boundary conditions as per (3.6)–(3.9) apply at all orders.

Integrating the leading-order momentum equations, i.e. (3.13) and (3.14), twice and applying periodicity shows that $\tilde {v}_0 = \tilde {v}_0(\tilde {y})$ and $\tilde {w}_0 = \tilde {w}_0(\tilde {y})$. Then, the leading-order continuity equation, (3.18), becomes $\partial \tilde {w}_{1}/\partial Z = -\mathrm {d}\tilde {v}_0/\mathrm {d}\tilde {y}.$ Integrating it with respect to $Z$ and applying periodicity on $\tilde {w}_1$ shows that $\mathrm {d}\tilde {v}_{0}/\mathrm {d}\tilde {y}=0$ and thus, from the mirror symmetry condition at $\tilde {y}=1$, that $\tilde {v}_{0}=0$, and $\tilde {w}_{1} = \tilde {w}_{1}(\tilde {y})$. Relatedly, (3.16) and (3.17) collapse to $\partial ^2 \tilde {v}_1/\partial Z^2 = 0$ and $\partial \tilde {p}_{{p,}0}/\partial Z =0$, respectively. By implication, $\tilde {v}_{1}=\tilde {v}_{1}(\tilde {y})$ and $\tilde {p}_{{p,}0}=\tilde {p}_{{p,}0}(\tilde {y})$. In turn, (3.21) along with the boundary conditions implies that $\tilde {v}_{1}=0$ and $\tilde {w}_{2} = \tilde {w}_{2}(\tilde {y})$. Equations (3.19) and (3.20) then reduce to $\partial ^{2}\tilde {v}_{2}/\partial Z^{2}=-2\,\mathrm {d}\tilde {p}_{{p,}0}/\mathrm {d}\tilde {y}$ and $\partial \tilde {p}_{{p},1}/\partial Z = -\mathrm {d}^{2}\tilde {w}_{0}/\mathrm {d}\tilde {y}^{2}/2-1$, respectively. Integrating the former with respect to $Z$ and applying periodicity shows that $\mathrm {d}\tilde {p}_{{p,}0}/\mathrm {d}\tilde {y} = 0$ and thus $\tilde {p}_{{p,}0}$ is a constant ($E_0$), and doing so again that $\tilde {v}_{2}=\tilde {v}_{2}(\tilde {y})$. Integrating the latter with respect to $Z$ and applying periodicity shows that $\tilde {p}_{{p,}1}=\tilde {p}_{{p,}1}(\tilde {y})$ and $\mathrm {d}^{2}\tilde {w}_{0}/\mathrm {d}\tilde {y}^{2}=-2$. Therefore, per the mirror symmetry condition,

(3.24)\begin{equation} \tilde{w}_{0}={-}\tilde{y}^{2}+2\tilde{y}+G_{0}, \end{equation}

where $G_0$ is a constant. In turn, the $O(\epsilon ^2)$ continuity equation (not shown) along with the boundary conditions implies that $\tilde {v}_{2}=0$ and $\tilde {w}_{3} = \tilde {w}_{3}(\tilde {y})$. Then, (3.22) and (3.23) reduce to $\partial ^{2}\tilde {v}_{3}/\partial Z^{2}=-2\,\mathrm {d}\tilde {p}_{{p,}1}/\mathrm {d}\tilde {y}$ and $\partial \tilde {p}_{{p},2}/\partial Z=-\mathrm {d}^{2}\tilde {w}_{1}/\mathrm {d}\tilde {y}^{2}/2$, respectively. The former leads to $\tilde {p}_{{p,}1} = E_1$ and $\tilde {v}_3=\tilde {v}_3(\tilde {y})$. The latter leads to $\tilde {p}_{{p,}2}=\tilde {p}_{{p,}2}(\tilde {y})$ and $\widetilde {w_1} = G_1$. Generalizing the foregoing process for $n \geqslant 2$, the $O(\epsilon ^n)$ continuity equation leads to $\tilde {v}_n = 0$ and, in turn, the $O(\epsilon ^{n-1})$ transverse- and streamwise-momentum equations lead to $\tilde {p}_{{p},n-1} = E_{n-1}$ and $\tilde {w}_{n-1} = G_{n-1}$, respectively. This leads to a unidirectional velocity profile to all algebraic orders in the outer region as per

(3.25)\begin{gather} \tilde{v} \sim \mathrm{e.s.t.}, \end{gather}

(3.26)\begin{gather}\tilde{w}\sim{-}\tilde{y}^{2}+2\tilde{y}+G\left(\epsilon\right)+ \mathrm{e.s.t.}, \end{gather}

where

(3.27)\begin{equation} G\left(\epsilon\right)=\sum_{n=0}^{\infty}G_{n}\epsilon^{n}. \end{equation}

The corresponding periodic component of the pressure field is $p_{p} = E(\epsilon )$, where $E(\epsilon ) = \sum _{n=0}^{\infty }E_{n}\epsilon ^{n}$, but it need not be resolved to find the slip length and Nusselt number.

3.1.2. Inner region

As in the case of parallel ridges, the inner region is where $Y\sim Z=O(1)$ as $\epsilon \rightarrow 0$, but, unlike in the case of parallel ridges, the flow is bidirectional. Using notation such that $\tilde {v}=V(Y,Z)$, $\tilde {w}=W(Y,Z)$ and $\tilde {p}_{p}={P}_{p}(Y,Z)$ in this region, the problem becomes

(3.28)\begin{gather} \frac{\partial V}{\partial Y}+\frac{\partial W}{\partial Z}=0 \end{gather}

(3.29)\begin{gather}\frac{Re}{\epsilon}\left(V\frac{\partial V}{\partial Y}+W\frac{\partial V}{\partial Z}\right)=\frac{2}{\epsilon}\frac{\partial P_{p}}{\partial Y}+\frac{1}{\epsilon^{2}}\left(\frac{\partial^{2}V}{\partial Y^{2}}+\frac{\partial^{2}V}{\partial Z^{2}}\right) \end{gather}

(3.30)\begin{gather}\frac{Re}{\epsilon}\left(V\frac{\partial W}{\partial Y}+W\frac{\partial W}{\partial Z}\right)=2+\frac{2}{\epsilon}\frac{\partial P_{p}}{\partial Z}+\frac{1}{\epsilon^{2}}\left(\frac{\partial^{2}W}{\partial Y^{2}}+\frac{\partial^{2}W}{\partial Z^{2}}\right), \end{gather}

subject to

(3.31)\begin{gather} V=\frac{\partial W}{\partial Y}=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ 0<\vert{Z}\vert<\delta \end{gather}

(3.32)\begin{gather}V=W=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ \delta<\vert{Z}\vert<1 \end{gather}

(3.33)\begin{gather} V\sim \mathrm{e.s.t.},\quad W\sim{-}\epsilon^{2}Y^{2}+2\epsilon Y+G(\epsilon)+\mathrm{e.s.t.},\quad P_{p}\sim E\left(\epsilon \right)+\mathrm{e.s.t.}\nonumber\\ \quad \mathrm{as}\ Y\rightarrow\infty\quad\mathrm{for}\ 0<\vert{Z}\vert<1 \end{gather}

(3.34)\begin{gather}\chi(Y,-1)=\chi(Y,1)\quad\mathrm{where}\ \chi=V,W,{P}_{p}, \frac{\partial V}{\partial Z},\frac{\partial W}{\partial Z}, \end{gather}

where (3.33) is the Van Dyke matching condition. We expand our velocities and the periodic component of the pressure field as

(3.35)\begin{gather} V\sim \sum_{n=0}^{\infty}\epsilon^{n}V_{n}+\mathrm{e.s.t.} \end{gather}

(3.36)\begin{gather}W\sim \sum_{n=0}^{\infty}\epsilon^{n}W_{n}+\mathrm{e.s.t.} \end{gather}

(3.37)\begin{gather}P_{p}\sim \sum_{n=0}^{\infty}\epsilon^{n} P_{{p},n}+\mathrm{e.s.t.} \end{gather}

Substituting the expansions into the momentum equations shows that, at leading (algebraic) order, i.e. $O(\epsilon ^{-2})$, $\nabla ^{2}V_{0}=0$ and $\nabla ^{2}W_{0}=0$. There is no shear rate at this order as per the matching condition; consequently, $V_{0}=W_{0}=0$ (and $V$ and $W$ are $o(1)$) and thus $G_{0}=0$ (and $G(\epsilon ) = o(1)$). The momentum equations become

(3.38)\begin{gather} \nabla^{2}V ={-}2\epsilon\frac{\partial P_{p}}{\partial Y}+O\left(\epsilon^{3}{Re}\right) \end{gather}

(3.39)\begin{gather}\nabla^{2}W ={-}2\epsilon^{2}-2\epsilon\frac{\partial P_{p}}{\partial Z}+O\left(\epsilon^{3}{Re}\right), \end{gather}

which shows that neglecting the inertial terms in them introduces an $O(\epsilon ^3{Re})$ error.

Writing

(3.40)\begin{gather} V = 2\epsilon\hat{V} \end{gather}

(3.41)\begin{gather}W ={-}\epsilon^{2}Y^{2}+2\epsilon\hat{W}, \end{gather}

we find that $\hat {V}$ and $\hat {W}$ satisfy

(3.42)\begin{gather} \frac{\partial \hat{V}}{\partial Y} + \frac{\partial \hat{W}}{\partial Z}=0 \end{gather}

(3.43)\begin{gather}\nabla^{2}\hat{V} ={-}\frac{\partial P_{p}}{\partial Y}+O\left(\epsilon^{2}{Re}\right) \end{gather}

(3.44)\begin{gather}\nabla^{2}\hat{W} ={-}\frac{\partial P_{p}}{\partial Z}+O\left(\epsilon^{2}{Re}\right), \end{gather}

subject to

(3.45)\begin{gather} \hat{V}=\frac{\partial \hat{W}}{\partial Y}=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ 0<\vert{Z}\vert<\delta \end{gather}

(3.46)\begin{gather}\hat{V}=\hat{W}=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ \delta<\vert{Z}\vert<1 \end{gather}

(3.47)\begin{gather}\hat{V}\sim \mathrm{e.s.t},\quad\hat{W}\sim Y+\frac{G(\epsilon)}{2\epsilon}+\mathrm{e.s.t.},\quad {P}_{p}\sim E\left(\epsilon\right)+\mathrm{e.s.t.}\ \mathrm{as}\ Y\rightarrow\infty \end{gather}

(3.48)\begin{gather}\chi(Y,-1)=\chi(Y,1)\quad\mathrm{where}\ \chi=\hat{V},\hat{W},{P}_{p},\frac{\partial \hat{V}}{\partial Z}, \frac{\partial \hat{W}}{\partial Z}. \end{gather}

This problem has been solved by Philip (Reference Philip1972a) by using conformal maps to find the streamfunction. The resulting velocity profile is

(3.49)\begin{gather} \hat{V}={-}\frac{\partial}{\partial Z}\left(Y\,{\rm Im}\left\{ \frac{1}{\rm \pi}\cos^{{-}1}\left[\frac{\cos \left({\rm \pi}\varTheta_{{\perp}}/2\right)}{\cos\left({\rm \pi}\delta/2\right)}\right]\right\} \right)+O \left(\epsilon^{2}{Re}\right) \end{gather}

(3.50)\begin{gather}\hat{W}=\left(Y\frac{\partial}{\partial Y}+1\right){\rm Im}\left\{ \frac{1}{\rm \pi}\cos^{{-}1}\left[\frac{\cos \left({\rm \pi}\varTheta_{{\perp}}/2\right)}{\cos\left({\rm \pi} \delta /2 \right)}\right]\right\}+O\left(\epsilon^{2} {Re}\right), \end{gather}

where $\varTheta _{\perp }=Z+\mathrm {i}Y$. Moreover, it follows from the far-field behaviour documented by Philip (Reference Philip1972a) for his result that

(3.51)\begin{equation} \hat{V}=O\left(\epsilon^{2}{Re}\right),\quad\hat{W}=Y+\frac{\lambda}{2} +O\left(\epsilon^{2}{Re}\right)\quad\mathrm{as}\ Y\rightarrow\infty, \end{equation}

such that

(3.52)\begin{equation} G\left(\epsilon\right)=\epsilon \lambda +O\left(\epsilon^3{Re}\right). \end{equation}

The complete inner velocity profile becomes

(3.53)\begin{gather} V={-}2\epsilon\frac{\partial}{\partial Z}\left(Y\,{\rm Im}\left\{ \frac{1}{\rm \pi}\cos^{{-}1}\left[ \frac{\cos\left({\rm \pi}\varTheta_{{\perp}}/2\right)}{\cos\left({\rm \pi}\delta/2\right)}\right]\right\} \right)+O \left(\epsilon^{3}{Re}\right) \end{gather}

(3.54)\begin{gather}W={-}\epsilon^{2}Y^{2}+2\epsilon\left(Y\frac{\partial}{\partial Y}+1\right){\rm Im}\left\{ \frac{1}{\rm \pi}\cos^{{-}1}\left[\frac{\cos\left({\rm \pi}\varTheta_{{\perp}}/2\right)}{\cos \left({\rm \pi}\delta/2\right)}\right]\right\} +O\left(\epsilon^{3}{Re}\right). \end{gather}

The periodic component of the pressure field is not required to resolve the thermal problem and is thus not developed here.

3.1.3. Composite solution

As in the case of parallel ridges, the outer solution keeps its form in the overlap region. Therefore, the inner solution is the composite one as per

(3.55)\begin{gather} \tilde{v}_{comp}={-}2\epsilon\frac{\partial}{\partial \tilde{z}}\left[\tilde{y}\,{\rm Im}\left\{ \frac{1}{\rm \pi} \cos^{{-}1}\left[\frac{\cos\left({\rm \pi} \left(\tilde{z}/\epsilon + \mathrm{i}\tilde{y}/\epsilon\right)/2\right)}{\cos \left({\rm \pi}\delta/2\right)}\right]\right\} \right]+O\left(\epsilon^{3}{Re}\right) \end{gather}

(3.56)\begin{gather}\tilde{w}_{comp}={-}\tilde{y}^{2}+2\epsilon\left(\tilde{y}\frac{\partial}{\partial \tilde{y}}+1\right) {\rm Im}\left( \frac{1}{\rm \pi}\cos^{{-}1}\left\{\frac{\cos\left[{\rm \pi} \left(\tilde{z}/\epsilon + \mathrm{i} \tilde{y}/\epsilon\right)/2\right]}{\cos\left({\rm \pi}\delta/2\right)}\right\}\right) +O\left(\epsilon^{3}{Re}\right). \end{gather}

3.1.4. Slip length

The procedure for finding the slip length by equating $\tilde {w}_{1{d}}$ as per (2.49) and the outer velocity profile in the case of parallel ridges also holds for transverse ones; consequently

(3.57)\begin{equation} \tilde{b} = \epsilon \lambda/2 + O\left(\epsilon^3{Re}\right); \end{equation}

half of its value for parallel ridges. Teo & Khoo (Reference Teo and Khoo2009) obtained this result by taking the limit of the dual-series equations accommodating the mixed boundary condition along the composite interface as $\epsilon \rightarrow 0$, but did not provide an error term of $O(\epsilon ^3{Re})$. Davies et al. (Reference Davies, Maynes, Webb and Woolford2006) numerically resolved the full problem, including the effect of viscous shear by the gas phase on the liquid. Then, the dimensionless slip length ($\tilde {b}$) further depends on the non-dimensional depth of the cavity, Reynolds number, etc. ($\epsilon$ need not be small in the numerical study by Davies et al. (Reference Davies, Maynes, Webb and Woolford2006).).

3.2.1. Formulation

The temperature field is decomposed into

(3.58)\begin{equation} T=\gamma{z}+T_{p}\left({y},{z}\right), \end{equation}

where $\gamma z$ is the linear temperature gradient in the liquid when a uniform heat flux of $\phi q''_{sl}$ is applied over the composite interface and $T_{p}({y},{z})$ is its periodic component, which repeats itself over a period of $2d$. It follows from an energy balance that

(3.59)\begin{equation} \gamma=\frac{q''_{sl}\phi}{\rho Q' c_{p}}, \end{equation}

where $Q'$ is the volumetric flow rate of liquid per unit depth of the domain. The dimensional form of the thermal energy equation is

(3.60)\begin{equation} v\frac{\partial T_{p}}{\partial y}+w\left(\gamma+\frac{\partial T_{p}}{\partial z}\right)=\alpha\left( \frac{\partial^2 T_{p}}{\partial y^2}+\frac{\partial ^2T_{p}}{\partial z^2}\right). \end{equation}

Unlike in the case of parallel ridges, the axial conduction term is finite. Moreover, we must retain it because, as shown below, it is essential in the inner region. Also, rather than being constant in the domain, the bulk temperature varies along its streamwise direction per

(3.61)\begin{equation} T_{m}\left(z\right)=\frac{\displaystyle\int\nolimits_0^H wT\,\mathrm{d}y}{\displaystyle\int\nolimits_0^Hw\,\mathrm{d}y}. \end{equation}

Consequently, the non-dimensional temperature is defined relative to the (bulk) mean temperature at the domain inlet as per

(3.62)\begin{equation} \tilde{T}=\frac{k\left[T-T_{m}\left({z}=0\right)\right]}{q_{sl}''H}. \end{equation}

We, again, work in terms of $Z$ such that our domain boundaries are independent of $\epsilon$. The dimensionless form of the thermal energy equation governing the periodic component of the temperature field is

(3.63)\begin{equation} {Pe}\left(\tilde{v}\frac{\partial \tilde{T}_{p}}{\partial \tilde{y}}+ \frac{\tilde{w}}{\epsilon}\frac{\partial \tilde{T}_{p}}{\partial Z}\right) +\frac{\tilde{w}\phi}{\tilde{Q}'}= \frac{\partial^2 \tilde{T}_{p}}{\partial \tilde{y}^2}+\frac{1}{\epsilon^2} \frac{\partial^2 \tilde{T}_{p}}{\partial Z^2}, \end{equation}

where the Péclet number Pe equals RePr, where the Prandtl number ${Pr}$ equals $c_{p}\mu /k$, i.e.

(3.64)\begin{equation} {Pe} = \frac{-\rho \beta H^3 c_{p}}{2\mu k}, \end{equation}

and $\tilde {Q}' = 2\mu Q'/(-\beta H^3)$. It is subject to the discontinuous (Neumann) boundary condition along the composite interface, symmetry along the channel centreline and periodicity on the upstream and downstream faces of the domain as per, respectively,

(3.65)\begin{gather} \frac{\partial \tilde{T}_{p}}{\partial{\tilde{y}}}=0\quad\mathrm{at}\ {\tilde{y}}=0\quad\mathrm{for}\ \vert Z\vert<{\delta} \end{gather}

(3.66)\begin{gather}\frac{\partial \tilde{T}_{p}}{\partial{\tilde{y}}}={-}1\quad\mathrm{at}\ {\tilde{y}}=0\quad\mathrm{for}\ \vert \delta \vert <\vert Z\vert<1 \end{gather}

(3.67)\begin{gather}\frac{\partial \tilde{T}_{p}}{\partial{\tilde{y}}}=0\quad\mathrm{at}\ \tilde{y}=1\quad\mathrm{for}\ \vert Z\vert<1 \end{gather}

(3.68)\begin{gather}\chi(\tilde{y},-1)=\chi(\tilde{y},1),\quad\mathrm{where}\ \chi=\tilde{T}_{p}, \frac{\partial \tilde{T}_{p}}{\partial \tilde{z}}. \end{gather}

3.2.2. Outer region

Recalling that, in the outer region, $\tilde {v}$ is exponentially small as per (3.25), the thermal energy equation becomes

(3.69)\begin{equation} {Pe} \frac{\tilde{w}}{\epsilon}\frac{\partial \tilde{T}_{p}}{\partial Z}+\frac{\tilde{w}\phi}{\tilde{Q}'} = \frac{\partial^2 \tilde{T}_{p}}{\partial \tilde{y}^2}+\frac{1}{\epsilon^2} \frac{\partial^2 \tilde{T}_{p}}{\partial Z^2} +\mathrm{e.s.t.} \end{equation}

Defining

(3.70)\begin{equation} \hat{T}_{p} =\tilde{Q}' \tilde{T}_{p}, \end{equation}

and recalling that the outer (streamwise) velocity profile is $\tilde {w} = -\tilde {y}^2 + 2\tilde {y} + G(\epsilon ) + \mathrm {e.s.t.}$, it becomes

(3.71)\begin{equation} \frac{Pe}{\epsilon}\left[-\tilde{y}^2 + 2\tilde{y} +G\left(\epsilon\right) \right] \frac{\partial \hat{T}_{p}}{\partial Z} +\phi \left[-\tilde{y}^2 + 2\tilde{y} + G\left(\epsilon \right)\right] = \frac{\partial^2 \hat{T}_{p}}{\partial \tilde{y}^2}+\frac{1}{\epsilon^2} \frac{\partial^2 \hat{T}_{p}}{\partial Z^2} + \mathrm{e.s.t.} \end{equation}

The boundary conditions are given by (3.67) and, for $\epsilon \ll \tilde {y} <1$, by (3.68), where $\tilde {T}_{p}$ is replaced by $\hat {T}_{p}$.

We, as an ansatz, assume that the periodic component of the temperature field is $O(1)$ and thus asymptotically expand it as

(3.72)\begin{equation} \hat{T}_{p} \sim \sum_{n=0}^{\infty}\epsilon^{n} \hat{T}_{{p},n}+\mathrm{e.s.t.}, \end{equation}

where $\hat {T}_{{p,}n}=O(1)$ for $n\geqslant 0$. Then, at the various orders of $\epsilon$, the thermal energy equation is

(3.73)\begin{gather} O\left(\epsilon^{{-}2}\right):\ \frac{\partial^{2}\hat{T}_{{p},0}}{\partial Z^{2}}=0 \end{gather}

(3.74)\begin{gather}O\left(\epsilon^{{-}1}\right):\ {Pe}\,\tilde{w}_0\frac{\partial \hat{T}_{{p},0}}{\partial Z} =\frac{\partial^{2}\hat{T}_{{p},1}}{\partial Z^{2}} \end{gather}

(3.75)\begin{gather}O\left(\epsilon^n\right):\ {Pe}\sum_{m=0}^{n+1}\tilde{w}_{n+1-m}\frac{\partial \hat{T}_{{p},m}}{\partial Z} +\phi \tilde{w}_n = \frac{\partial^{2}\hat{T}_{{p},n}}{\partial \tilde{y}^{2}} +\frac{\partial^{2}\hat{T}_{{p},n+2}}{\partial Z^{2}}\quad\mathrm{for}\ n\geqslant 0, \end{gather}

where the $\tilde {w}_n$ are independent of $Z$ as per (3.26). The periodicity boundary condition applied, in succession, to (3.73) and (3.74) shows that $\hat {T}_{{p},n} = \hat {T}_{{p},n}(\tilde {y})$ for $n =0,1$. Continuing this process for the thermal energy equation for $O(\epsilon ^n)$, where $n \geqslant 0$, shows that

(3.76)\begin{equation} \frac{\partial^{2}\hat{T}_{{p,}n+2}}{\partial Z^{2}} =\frac{\partial^{2}\hat{T}_{{p,}n+2}}{\partial Z^{2}}\left({Pe},\phi,\tilde{w}_0,\ldots,\tilde{w}_{n+1}, \hat{T}_{{p,}0},\ldots,\hat{T}_{{p,}n+1} \right). \end{equation}

Moreover, once the thermal energy equation for $O(\epsilon ^n)$ has been reached, $\hat {T}_{{p},i}$ for $i \leqslant n+1$ has been shown to be only a function of $\tilde {y}$. Therefore, the periodicity boundary condition implies that $\hat {T}_{{p},n}$ is independent of $Z$ (and thus Pe) for all $n$. Making further use of (3.26), the thermal energy equation becomes

(3.77)\begin{gather} O\left(\epsilon^0\right): \frac{\mathrm{d}^{2}\hat{T}_{{p},0}}{\mathrm{d}\tilde{y}^{2}}=\phi \left(-\tilde{y}^2 + 2\tilde{y} + G_0\right) \end{gather}

(3.78)\begin{gather}O\left(\epsilon^n\right):\ \frac{\mathrm{d}^{2}\hat{T}_{{p},n}}{\mathrm{d}\tilde{y}^{2}}=\phi G_n,\quad n \geqslant 1. \end{gather}

Integrating these equations, applying the symmetry boundary condition at $\tilde {y} =1$ and integrating them again yields

(3.79)\begin{gather} O\left(\epsilon^0\right):\ \hat{T}_{{p},0} =\phi \left[ -\frac{\left(\tilde{y}-1\right)^4}{12} +\frac{\left(\tilde{y}-1\right)^2}{2} + G_0 \frac{\left(\tilde{y}-1\right)^2}{2} \right] + \alpha_0 \end{gather}

(3.80)\begin{gather}O\left(\epsilon^n\right):\ \hat{T}_{{p},n} = \phi G_n \frac{\left(\tilde{y}-1\right)^2}{2} + \alpha_n,\quad n \geqslant 1, \end{gather}

where the $\alpha _n$ are constants. Consequently

(3.81)\begin{equation} \hat{T}_{p} = \phi \left[ -\frac{\left(\tilde{y}-1\right)^4}{12} +\frac{\left(\tilde{y}-1\right)^2}{2} + G\left(\epsilon\right)\frac{\left(\tilde{y}-1\right)^2}{2} \right]+\sum_{n=0}^{\infty}\alpha_n \epsilon^n +\mathrm{e.s.t.}\end{equation}

Finally, it follows from the slip length given by (3.57) that $\tilde {Q}' = 2/3 + \epsilon \lambda + O(\epsilon ^3{Re})$ and thus

(3.82)\begin{equation} \frac{1}{\tilde{Q}'} = \frac{1}{2/3 +\epsilon \lambda} + O\left(\epsilon^3{Re}\right), \end{equation}

and, as per (3.52), $G(\epsilon )=\epsilon \lambda +O(\epsilon ^3{Re})$. Therefore, reverting from $\hat {T}_{p}$ back to $\tilde {T}_{p}$, (3.81) is expressed as

(3.83)\begin{equation} \tilde{T}_{p}=\frac{\phi}{2/3+\epsilon \lambda}\left[-\frac{\left(\tilde{y}-1\right)^{4}}{12}+\left(\frac{1}{2}+\frac{\epsilon \lambda}{2}\right)\left(\tilde{y}-1\right)^{2}+H\left(\epsilon\right)+\mathrm{e.s.t.}\right] +O\left(\epsilon^3{Re}\right), \end{equation}

where $H_n = \alpha _n/\phi$ and

(3.84)\begin{equation} H\left(\epsilon\right)=\sum_{n=0}^{\infty} H_n \epsilon^n.\end{equation}

This outer temperature profile is the same as that for parallel ridges as per (2.73), except that the (asymptotic limit of the) slip length for transverse ridges ($\epsilon \lambda /2$) replaces that for parallel ones ($\epsilon \lambda$) and the error is $O(\epsilon ^3{Re})$ rather than exponentially small.

3.2.3. Inner region

Denoting $\tilde {T}_{p}$ by ${\theta }_{p}(Y,Z)$ in the inner region, the thermal energy equation is

(3.85)\begin{equation} \frac{Pe}{\epsilon}\left(V\frac{\partial{\theta}_{p}}{\partial Y} +W\frac{\partial{\theta}_{p}}{\partial Z}\right)+ \frac{W\phi}{\tilde{Q}'} =\frac{1}{\epsilon^{2}}\left(\frac{\partial^{2}{\theta}_{p}}{\partial Y^{2}}+\frac{\partial^{2}{\theta}_{p}}{\partial Z^{2}}\right),\end{equation}

where $V$ and $W$ are given by (3.53) and (3.54), respectively, and are $O(\epsilon )$. The boundary conditions are

(3.86)\begin{gather} \frac{\partial\theta_{p}}{\partial Y}=0\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ 0< Z<\delta \end{gather}

(3.87)\begin{gather}\frac{\partial\theta_{p}}{\partial Y}={-}\epsilon\quad\mathrm{at}\ Y=0\quad\mathrm{for}\ \delta< Z<1 \end{gather}

(3.88)\begin{gather}\theta_{p}\sim{-}\epsilon\phi Y+\phi\left[\frac{1}{2}+\frac{1/12+H(\epsilon)}{2/3+\epsilon \lambda}\right] +O\left(\epsilon^{3},\epsilon^3{Re}\right)\quad\mathrm{as}\ Y\rightarrow\infty \end{gather}

(3.89)\begin{gather}\chi(\tilde{y},Z)=\chi(\tilde{y},Z+2),\quad\mathrm{where}\ \chi=\theta_{p}, \frac{\partial\theta_{p}}{\partial \tilde{z}}, \end{gather}

where (3.88), which follows from the manipulation of (3.83), is the matching condition. As an ansatz, we assume that the periodic component of the dimensionless, inner temperature field is $O(1)$ and thus asymptotically expand it as

(3.90)\begin{equation} {\theta}_{p} \sim \sum_{n=0}^{\infty}\epsilon^{n} {\theta}_{{p},n}+\mathrm{e.s.t.} \end{equation}

The thermal energy equation reduces to Laplace's equation at leading order, i.e. $O(\epsilon ^{-2})$; therefore,

(3.91)\begin{equation} \frac{\partial \theta_{{p},0}}{\partial Y} \sim \frac{\partial \theta_{{p},0}}{\partial Z}. \end{equation}

These quantities do not exceed their values near the ridge and are thus $O(\epsilon )$ as per the boundary condition along it as per (3.87) and $V,W=O(\epsilon )$; consequently, the thermal energy equation reduces to

(3.92)\begin{equation} \frac{\partial^{2}\theta_{p}}{\partial Y^{2}}+\frac{\partial^{2}\theta_{p}}{\partial Z^{2}}= O\left(\epsilon^3{Pe},\epsilon^3\right). \end{equation}

Accepting an accuracy of $O(\epsilon ^{2})$ for $\theta _{p}$, the periodic boundary conditions along the streamwise borders of the domain manifest themselves as symmetry conditions on a diffusion problem such that (3.89) is replaced by

(3.93)\begin{equation} \frac{\partial\theta_{p}}{\partial Z}=0\quad\mathrm{at}\ \left|Z\right|=1. \end{equation}

Relatedly, since the outer problem has no dependence on $Z$, the temperature field and thus local Nusselt number are symmetric about $Z = 0$ and we henceforth need only consider our domain to extend from 0 to $1$. The solution given by Mikic (Reference Mikic1957) is again employed such that the dimensionless, inner periodic temperature profile is given by

(3.94)\begin{align} \theta_{p}&={-}\epsilon \phi Y - \frac{2\epsilon}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} Z\right)\mathrm{e}^{{-}n {\rm \pi}Y}}{n^2}+\phi\left[\frac{1}{2}+\frac{1/12+H(\epsilon)}{2/3+\epsilon \lambda}\right]\nonumber\\ &\quad + O\left(\epsilon^3,\epsilon^3{Re},\epsilon^3{Pe}\right). \end{align}

Turning our attention to finding $H(\epsilon )$, the details are provided in Appendix E. The temperature along the composite interface becomes

(3.95)\begin{align} \theta \left(X,0\right)&= \phi \hat{H}\left(\epsilon \lambda\right) - \frac{2\epsilon}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} Z\right)}{n^2}\nonumber\\ &\quad + O\left( \frac{\epsilon^3}{Pe}, \frac{\epsilon^3}{Pr},\epsilon^3,\epsilon^3{Re},\epsilon^5{Pe}, \epsilon^7{Re}\,{Pe} \right), \end{align}

where $\hat {H}(\epsilon \lambda )=\hat {D}(\epsilon \lambda /2 )$.

3.2.4. Composite solution

The composite solution follows from (2.85) and includes the linear component of the temperature field. It is

(3.96) \begin{align} \tilde{T}_{comp} &= \frac{\phi}{2/3+\epsilon \lambda}\left[\frac{\tilde{z}}{Pe} -\frac{\left(\tilde{y}-1\right)^{4}}{12}+\left(\frac{1}{2}+\frac{\epsilon\lambda}{2}\right) \left(\tilde{y}-1\right)^{2}+\breve{H}\left(\epsilon\right)\right] \nonumber\\ &\quad- \frac{2\epsilon}{{\rm \pi}^2}\sum_{n=1}^{\infty}\frac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} \tilde{z}/\epsilon\right) \mathrm{e}^{{-}n {\rm \pi} \tilde{y}/\epsilon}}{n^2}\nonumber\\ &\quad +O\left( \frac{\epsilon^3}{Pe}, \frac{\epsilon^3}{Pr},\epsilon^3,\epsilon^4{Re},\epsilon^3{Pe}, \epsilon^7{Re}\,{Pe} \right). \end{align}

3.2.5. Nusselt numbers

The Nusselt number along the ridge follows from (2.87) and (3.95) as

(3.97)\begin{align} {Nu} &= \frac{4}{\phi \hat{D}\left(\epsilon \lambda/2\right) - \dfrac{2\epsilon}{{\rm \pi}^2} \sum_{n=1}^{\infty}\dfrac{\sin\left(n{\rm \pi}\delta\right)\cos\left(n{\rm \pi} Z\right)}{n^2}}\nonumber\\ &\quad +O\left(\frac{\epsilon^3}{Pe}, \frac{\epsilon^3}{Pr},\epsilon^3,\epsilon^4{Re},\epsilon^5{Pe}, \epsilon^7{Re}\,{Pe} \right), \end{align}

and, as in the case of parallel ridges, approaches $\infty$ as solid fraction approaches 0. Regularization of this expression and, subsequently, an expansion of it in terms of $\epsilon /[\hat {D}(\epsilon \lambda /2)+\epsilon \lambda ]$ follows in the same manner as in the case of parallel ridges. The result is that the mean Nusselt number $\overline {Nu}$ is given by (2.94) when $\hat {D} (\epsilon \lambda )$ is replaced by $\hat {D}(\epsilon \lambda /2)$ and the error term is replaced by that in (3.97). For the alternate definition of the mean Nusselt number, it follows, by replacing $\epsilon \lambda$ with $\epsilon \lambda /2$ in (2.99) and adjusting the error term, that

(3.98)\begin{align} \overline{Nu}'&= 140\left(1+\frac{3\epsilon \lambda}{2}\right)^2\left/\left\{ 17+ \left[42\lambda + \frac{70}{\phi^2{\rm \pi}^3}\sum_{n=1}^{\infty}\frac{\sin^2\left(n{\rm \pi}\delta\right)}{n^3}\right] \epsilon \right. \right. \nonumber\\ &\quad \left.+\left[ \frac{105\lambda^2}{4} +\frac{210\lambda}{\phi^2{\rm \pi}^3}\sum_{n=1}^{\infty}\frac{\sin^2\left(n{\rm \pi}\delta\right)}{n^3} \right]\epsilon^2 +\left[ \frac{630\lambda^2}{4\phi^2{\rm \pi}^3}\sum_{n=1}^{\infty}\frac{\sin^2\left(n{\rm \pi}\delta\right)}{n^3}\right]\epsilon^3\right\} \nonumber\\ &\quad + O\left( \frac{\epsilon^3}{Pe}, \frac{\epsilon^3}{Pr},\epsilon^3,\epsilon^4{Re},\epsilon^5{Pe}, \epsilon^7{Re}\,{Pe} \right). \end{align}

In summary, table 1 applies transverse ridges when $\hat {D} (\epsilon \lambda )$ is replaced by $\hat {D}(\epsilon \lambda /2)$ in the Nusselt number expressions and the $O(\epsilon ^3)$ error terms are replaced by $O( {\epsilon ^3}/{Pe}, {\epsilon ^3}/{Pr},\epsilon ^3,\epsilon ^4{Re},\epsilon ^5{Pe}, \epsilon ^7{Re}\,{Pe} )$.