2.1 Governing equations in the moving frame

We consider a two-dimensional (2-D) horizontal channel of height $h$ , as illustrated in figure 1. A Cartesian coordinate system $(x^{\ast },y^{\ast })$ is defined at the bottom wall, with $x^{\ast }$ pointing horizontally along the channel and $y^{\ast }$ pointing upward. A travelling sinusoidal temperature wave is imposed at the bottom wall. Since the temperature variation is small compared with the mean temperature, we apply the Boussinesq approximation, and the flow and temperature fields obey

(2.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}}+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\ast }\boldsymbol{u}^{\ast }=-{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}\unicode[STIX]{x1D735}^{\ast }p^{\ast }+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{\ast 2}\boldsymbol{u}^{\ast }+\unicode[STIX]{x1D6FD}g\left(T^{\ast }-T_{f}\right)\boldsymbol{j}, & \displaystyle\end{eqnarray}$$

(2.1b ) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}^{\ast }\boldsymbol{\cdot }\boldsymbol{u}^{\ast }=0, & \displaystyle\end{eqnarray}$$

(2.1c ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}T^{\ast }}{\unicode[STIX]{x2202}t^{\ast }}}+\boldsymbol{u}^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\ast }T^{\ast }=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FB}^{\ast 2}T^{\ast }, & \displaystyle\end{eqnarray}$$

(2.1d ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}^{\ast }(t^{\ast },x^{\ast },y^{\ast }=0)=\boldsymbol{u}^{\ast }(t^{\ast },x^{\ast },y^{\ast }=h)=0, & \displaystyle\end{eqnarray}$$

(2.1e ) $$\begin{eqnarray}\displaystyle & \displaystyle T^{\ast }(t^{\ast },x^{\ast },y^{\ast }=0)=T_{f}+A\sin [k(x^{\ast }-ct^{\ast })],T^{\ast }(t^{\ast },x^{\ast },y^{\ast }=h)=T_{f}, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D735}^{\ast }=\boldsymbol{i}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x^{\ast }+\boldsymbol{j}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y^{\ast }$

$\boldsymbol{u}^{\ast }=u^{\ast }\boldsymbol{i}+v^{\ast }\boldsymbol{j}$

$\unicode[STIX]{x1D70C}$

$\unicode[STIX]{x1D708}$

$\unicode[STIX]{x1D6FD}$

$g$

$T^{\ast }$

$\unicode[STIX]{x1D6FC}$

$T_{f}$

$A$

$k$

$c$

$x^{\prime }=x^{\ast }-ct^{\ast }$

$y^{\prime }=y^{\ast }$

$t^{\prime }=t^{\ast }$

$u^{\prime }=u^{\ast }$

$v^{\prime }=v^{\ast }$

$p^{\prime }=p^{\ast }$

$T^{\prime }=T^{\ast }$

(2.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle (\boldsymbol{u}^{\prime }-c\,\boldsymbol{i})\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }\boldsymbol{u}^{\prime }=-{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}}}\unicode[STIX]{x1D735}^{\prime }p^{\prime }+\unicode[STIX]{x1D708}\unicode[STIX]{x1D6FB}^{\prime 2}\boldsymbol{u}^{\prime }+\unicode[STIX]{x1D6FD}g\left(T^{\prime }-T_{f}\right)\boldsymbol{j},\\[2.0pt] \unicode[STIX]{x1D735}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime }=0,\\ (\boldsymbol{u}^{\prime }-c\,\boldsymbol{i})\boldsymbol{\cdot }\unicode[STIX]{x1D735}^{\prime }T^{\prime }=\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FB}^{\prime 2}T^{\prime },\\ \boldsymbol{u}^{\prime }(x^{\prime },y^{\prime }=0)=\boldsymbol{u}^{\prime }(x^{\prime },y^{\prime }=h)=0,\\ T^{\prime }(x^{\prime },y^{\prime }=0)=T_{f}+A\sin (kx^{\prime }),T^{\prime }(x^{\prime },y^{\prime }=h)=T_{f},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D735}^{\prime }=\boldsymbol{i}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x^{\prime }+\boldsymbol{j}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y^{\prime }$ and $\boldsymbol{u}^{\prime }=u^{\prime }\boldsymbol{i}+v^{\prime }\boldsymbol{j}$ . By defining $\boldsymbol{u}^{\prime }=\boldsymbol{u}^{\ast }$ , we are solving the absolute velocity in the moving reference frame. Consequently, the boundary conditions for $\boldsymbol{u}^{\prime }$ remain homogeneous, but the momentum and heat equations gain new convective terms. The solution is periodic in $x^{\prime }$ because of the periodic boundary conditions. There are sufficient differential equations and boundary conditions to solve the problem. The induced axial flow rate along the channel can then be found as

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle Q^{\prime }=\int _{0}^{h}u^{\prime }\,\text{d}y^{\prime }. & \displaystyle\end{eqnarray}$$

The main purpose of this paper is to find $Q^{\prime }$ .

2.2 Non-dimensionalization

The channel height $h$ provides a natural length scale. A dominant balance between the driving buoyancy term and the viscous term in the vertical momentum equation yields a velocity scale:

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle V=\unicode[STIX]{x1D6FD}gAh^{2}/\unicode[STIX]{x1D708}. & \displaystyle\end{eqnarray}$$

Based on these scales, a set of dimensionless variables is defined:

(2.5a-g ) $$\begin{eqnarray}\displaystyle x={\displaystyle \frac{x^{\prime }}{h}},\quad y={\displaystyle \frac{y^{\prime }}{h}},\quad u={\displaystyle \frac{u^{\prime }}{V}},\quad v={\displaystyle \frac{v^{\prime }}{V}},\quad p={\displaystyle \frac{hp^{\prime }}{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D708}V}},\quad T={\displaystyle \frac{T^{\prime }-T_{f}}{A}},\quad Q={\displaystyle \frac{Q^{\prime }}{Vh}}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The governing equations and boundary conditions are non-dimensionlized:

(2.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle (Re\,\boldsymbol{u}-Rc\,\boldsymbol{i})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+T\boldsymbol{j}, & \displaystyle\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

(2.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle (Re\,\boldsymbol{u}-Rc\,\boldsymbol{i})\boldsymbol{\cdot }\unicode[STIX]{x1D735}T=Pr^{-1}\unicode[STIX]{x1D6FB}^{2}T, & \displaystyle\end{eqnarray}$$

(2.6d ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}(x,y=0)=\boldsymbol{u}(x,y=1)=0, & \displaystyle\end{eqnarray}$$

(2.6e ) $$\begin{eqnarray}\displaystyle & \displaystyle T(x,y=0)=\sin (Kx),T(x,y=1)=0, & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D735}=\boldsymbol{i}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}x+\boldsymbol{j}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}y$

$\boldsymbol{u}=u\boldsymbol{i}+v\boldsymbol{j}$

$Re=Vh/\unicode[STIX]{x1D708}$

$V$

$Rc=ch/\unicode[STIX]{x1D708}$

$c$

$Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D6FC}$

$K=kh$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle Q=\int _{0}^{1}u\,\text{d}y. & \displaystyle\end{eqnarray}$$

In the next section, we seek a solution in the limit $Re\rightarrow 0$ because $Re$ is small in microchannels. A full numerical solution will be presented in § 4.

3.1 Leading order

At the leading order, the governing equations and boundary conditions are

(3.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle -Rc{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{\mathbf{0}}}{\unicode[STIX]{x2202}x}}=-\unicode[STIX]{x1D735}p_{0}+\unicode[STIX]{x1D735}^{2}\boldsymbol{u}_{\boldsymbol{ 0}}+T_{0}\,\boldsymbol{j},\\[8.0pt] \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\mathbf{0}}=0,\\[8.0pt] \displaystyle -Rc{\displaystyle \frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}x}}=Pr^{-1}\unicode[STIX]{x1D6FB}^{2}T_{0},\\[8.0pt] \displaystyle \boldsymbol{u}_{\mathbf{0}}(x,0)=\boldsymbol{u}_{\mathbf{0}}(x,1)=0,\quad T_{0}(x,0)=\sin (Kx),\quad T_{0}(x,1)=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The pressure term is then eliminated from the momentum equation by cross-differentiation. The horizontal velocity $u_{0}$ is eliminated from the combined momentum equation using the continuity equation, and we are left with

(3.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle Rc\unicode[STIX]{x1D6FB}^{2}{\displaystyle \frac{\unicode[STIX]{x2202}v_{0}}{\unicode[STIX]{x2202}x}}+\unicode[STIX]{x1D6FB}^{4}v_{0}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}T_{0}}{\unicode[STIX]{x2202}x^{2}}}=0,\\[8.0pt] \displaystyle Rc{\displaystyle \frac{\unicode[STIX]{x2202}T_{0}}{\unicode[STIX]{x2202}x}}+Pr^{-1}\unicode[STIX]{x1D6FB}^{2}T_{0}=0,\\[8.0pt] \displaystyle v_{0}(x,0)=v_{0}(x,1)={\displaystyle \frac{\unicode[STIX]{x2202}v_{0}}{\unicode[STIX]{x2202}y}}(x,0)={\displaystyle \frac{\unicode[STIX]{x2202}v_{0}}{\unicode[STIX]{x2202}y}}(x,1)=0,\\[8.0pt] T_{0}(x,0)=\sin (Kx),\quad T_{0}(x,1)=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Since the thermal boundary condition is periodic in $x$ , we take

(3.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle v_{0}=\text{Im}[F(y)\text{e}^{\text{i}Kx}], & \displaystyle\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle T_{0}=\text{Im}[G(y)\text{e}^{\text{i}Kx}], & \displaystyle\end{eqnarray}$$

$\text{Im}$

(3.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\text{d}^{4}F}{\text{d}y^{4}}}+(\text{i}RcK-2K^{2}){\displaystyle \frac{\text{d}^{2}F}{\text{d}y^{2}}}-K^{2}(\text{i}RcK-K^{2})F=K^{2}G, & \displaystyle\end{eqnarray}$$

(3.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\text{d}^{2}G}{\text{d}y^{2}}}+(\text{i}RcKPr-K^{2})G=0, & \displaystyle\end{eqnarray}$$

(3.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle F(0)=F(1)={\displaystyle \frac{\text{d}F}{\text{d}y}}(0)={\displaystyle \frac{\text{d}F}{\text{d}y}}(1)=0, & \displaystyle\end{eqnarray}$$

(3.6d ) $$\begin{eqnarray}\displaystyle & \displaystyle G(0)=1,G(1)=0. & \displaystyle\end{eqnarray}$$

(3.7) $$\begin{eqnarray}\displaystyle & \displaystyle G(y)={\displaystyle \frac{\sin [\sqrt{\text{i}RcKPr-K^{2}}(1-y)]}{\sin [\sqrt{\text{i}RcKPr-K^{2}}]}}. & \displaystyle\end{eqnarray}$$

Equation (3.6a ) then yields an analytic solution for $F(y)$ , which is listed in appendix A because of its length. The horizontal velocity $u_{0}$ is found from $v_{0}$ as

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle u_{0}=\text{Im}\left[{\displaystyle \frac{\text{i}}{K}}{\displaystyle \frac{\text{d}F}{\text{d}y}}\text{e}^{\text{i}Kx}\right]. & \displaystyle\end{eqnarray}$$

Thus, the leading-order volume flow rate is found to be zero:

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{0}=\int _{0}^{1}u_{0}\,\text{d}y=0. & \displaystyle\end{eqnarray}$$

3.2 First order

The first-order governing equations and boundary conditions are

(3.10a ) $$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D735}p_{1}+\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}_{\mathbf{1}}+T_{1}\,\boldsymbol{j}+Rc{\displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}_{\mathbf{1}}}{\unicode[STIX]{x2202}x}}=\boldsymbol{u}_{\mathbf{0}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{\mathbf{0}}, & \displaystyle\end{eqnarray}$$

(3.10b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{\mathbf{1}}=0, & \displaystyle\end{eqnarray}$$

(3.10c ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{1}{Pr}}\unicode[STIX]{x1D6FB}^{2}T_{1}+Rc{\displaystyle \frac{\unicode[STIX]{x2202}T_{1}}{\unicode[STIX]{x2202}x}}=\boldsymbol{u}_{\mathbf{0}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}T_{0}, & \displaystyle\end{eqnarray}$$

(3.10d ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{u}_{\mathbf{1}}(x,0)=\boldsymbol{u}_{\mathbf{1}}(x,1)=T_{1}(x,0)=T_{1}(x,1)=0. & \displaystyle\end{eqnarray}$$

$Q$

$u_{1}$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle -{\displaystyle \frac{\unicode[STIX]{x2202}p_{1}}{\unicode[STIX]{x2202}x}}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}u_{1}}{\unicode[STIX]{x2202}x^{2}}}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}u_{1}}{\unicode[STIX]{x2202}y^{2}}}+Rc{\displaystyle \frac{\unicode[STIX]{x2202}u_{1}}{\unicode[STIX]{x2202}x}}=u_{0}{\displaystyle \frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}x}}+v_{0}{\displaystyle \frac{\unicode[STIX]{x2202}u_{0}}{\unicode[STIX]{x2202}y}}. & \displaystyle\end{eqnarray}$$

Averaging (3.11) in the $x$ direction over the thermal wavelength $2\unicode[STIX]{x03C0}/K$ yields

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\text{d}^{2}\overline{u}_{1}}{\text{d}y^{2}}}={\displaystyle \frac{\text{d}\overline{u_{0}v_{0}}}{\text{d}y}}, & \displaystyle\end{eqnarray}$$

where the overline means horizontal average over the thermal wavelength. All other terms vanish owing to periodicity in $x$ . The equation states that the averaged horizontal viscous shear stress is balanced by the Reynolds stress. Solving (3.12) together with the boundary conditions $\overline{u}_{1}(0)=\overline{u}_{1}(1)=0$ gives

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{u}_{1}(y)=\int _{0}^{y}\overline{u_{0}v_{0}}\,\text{d}\unicode[STIX]{x1D709}-y\int _{0}^{1}\overline{u_{0}v_{0}}\,\text{d}\unicode[STIX]{x1D709}=\int _{0}^{y}\left(\overline{u_{0}v_{0}}-\langle \overline{u_{0}v_{0}}\rangle \right)\,\text{d}\unicode[STIX]{x1D709}. & \displaystyle\end{eqnarray}$$

The angle brackets denote vertical average over the entire interval $[0,1]$ , so this term is a constant. The first-order volume flow rate $Q_{1}$ is found as

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{1}=\int _{0}^{1}\overline{u}_{1}(y)\,\text{d}y=\int _{0}^{1}\int _{0}^{y}\left(\overline{u_{0}v_{0}}-\langle \overline{u_{0}v_{0}}\rangle \right)\,\text{d}\unicode[STIX]{x1D709}\,\text{d}y. & \displaystyle\end{eqnarray}$$

If $F(y)$ is decomposed into real and imaginary parts as $F_{R}(y)+\text{i}F_{I}(y)$ , $v_{0}$ and $u_{0}$ can be found by substituting it into (3.5a ) and (3.8), respectively. Then, the Reynolds stress is found as

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{u_{0}v_{0}}={\displaystyle \frac{1}{2K}}\left(F_{I}{\displaystyle \frac{\text{d}F_{R}}{\text{d}y}}-F_{R}{\displaystyle \frac{\text{d}F_{I}}{\text{d}y}}\right). & \displaystyle\end{eqnarray}$$

The integrals in (3.14) are evaluated numerically to an accuracy of seven significant digits (appendix B). Thus, the volume flow rate $Q$ is

(3.16) $$\begin{eqnarray}\displaystyle Q(Re,Rc,Pr,K)=ReQ_{1}(Rc,Pr,K). & & \displaystyle\end{eqnarray}$$

This concludes the asymptotic solution of the volume flow rate. In the next section, we present the numerical method used to solve (2.6), and two typical solutions.

5.1 Reynolds number $Re$

The asymptotic solution in (3.16) gives $Q=ReQ_{1}(Rc,Pr,K)$ . Thus $Q/Re$ is plotted against $Re$ in figure 5 for $Pr=1$ , $K=2\unicode[STIX]{x03C0}$ and various $Rc$ . We expect that the numerical solutions would agree with the asymptotic predictions at $Re\ll 1$ , and the two would diverge starting at $Re\sim 1$ . Surprisingly, the two maintain agreement all the way to $Re=1000$ for all $Rc$ cases studied. We also try other $Pr$ and $K$ values and observe the same agreement. This is unexpected, especially considering that the velocity fields already start to separate at $Re=1000$ in figure 3. We compute the $Rc=10$ case for $Re>10^{3}$ , and find that the asymptotic prediction begins to deviate from the numerical solution at $Re=10^{4}$ . We double and redouble the mesh resolutions and conclude that the $Q$ values for these three high $Re$ cases are accurate to at least six significant digits. Thus, the deviation is physical and not caused by numerical errors. Similar deviations are observed for other $Pr$ and $K$ values.

Figure 5.

Ratio of volume flow rate $Q$ to $Re$ versus $Re$ for $K=2\unicode[STIX]{x03C0}$ , $Pr=1$ and various $Rc$ values. The $Rc=10$ case is pushed to extra high $Re$ to show the divergence of the numerical and asymptotic solutions.

The horizontally averaged velocity $\bar{u}(y)$ is studied for $Re=1$ and $1000$ . Since $Q$ is calculated by a vertical integration of $\bar{u}(y)$ , the numerical and asymptotic solutions of $\bar{u}(y)$ will help to explain why the two solutions of $Q$ agree well at $Re=10^{4}$ , whereas those of the velocity field start to deviate at $Re=10^{3}$ . The velocity $\bar{u}(y)$ is found by averaging $u(x,y)$ over the wavelength $2\unicode[STIX]{x03C0}/K$ in the horizontal direction. Figure 6 shows $\bar{u}(y)$ for the two cases studied in figures 2 ( $Re=1$ ) and 3 ( $Re=1000$ ). The leading-order asymptotic solution $u_{0}(x,y)$ in (3.8) consists of simple $\sin (x)$ and $\cos (x)$ functions and averages to zero. The first-order solution $u_{1}(x,y)$ contains $\sin ^{2}(x)$ and $\cos ^{2}(x)$ terms and averages to non-zero. The good agreement between $Re\,\bar{u}_{1}(y)$ and the numerical solution for $Re=1000$ shows that the differences between the velocity fields in figure 3 mostly cancel themselves when averaged. Since $Q$ is the vertical integration of $\bar{u}(y)$ , and can be considered as an average of an average, the asymptotic solution of $Q$ is accurate at much higher $Re$ than that of $\boldsymbol{u}_{0}$ .

Figure 6.

The horizontally averaged velocity profile $\bar{u}(y)$ . The cases in figures 2 and 3 are plotted in (a) and (b), respectively. Shown in each panel are the numerical solution (open circles), the leading-order asymptotic solution $\bar{u}_{0}$ (dashed line) and the first-order solution $Re\,\bar{u}_{1}$ (solid line).

Since $Re$ is not high for flow in microchannels, we study $Q$ only for $Re\leqslant 10^{3}$ , and rely solely on the asymptotic solution $Q_{1}$ to explore the effects of other parameters ( $Rc$ , $K$ and $Pr$ ).

5.2 Reynolds number $Rc$ based on the wave speed

Figure 5 shows that $Q$ varies non-monotonically with $Rc$ ; $|Q|$ increases then decreases with $Rc$ , and vanishes in the small and large $Rc$ limits. The limiting behaviour can be understood as follows. When $Rc=0$ , the flow consists of convection cells that span across the channel and therefore $Q=0$ ; when $Rc\rightarrow \infty$ , both flow velocity and $Q$ vanish owing to the fast alternation of heating and cooling at the bottom wall. The asymptotic solutions of $Q_{1}$ in the limits $Rc\rightarrow 0$ and $\infty$ are derived next.

5.2.1 Asymptotic solution of $Q_{1}$ in the limit $Rc\rightarrow 0$

Instead of expanding the complex solution of $F(y)$ in appendix A, we find it more revealing to start from the governing equation (3.6a ), and expand $F$ in an asymptotic series in $Rc$ :

(5.1) $$\begin{eqnarray}\displaystyle F(y)=F_{0}(y)+RcF_{1}(y)+\cdots \,. & & \displaystyle\end{eqnarray}$$

The solution of $G(y)$ in (3.7) is expanded in the limit $Rc\rightarrow 0$ as

(5.2) $$\begin{eqnarray}\displaystyle G(y) & = & \displaystyle \text{csch}(K)\sinh [K(1-y)]\nonumber\\ \displaystyle & & \displaystyle +\,{\textstyle \frac{1}{4}}\text{i}PrRc\,\text{csch}^{2}(K)\{(y-2)\sinh (Ky)-y\sinh [K(y-2)]\}+O(Rc^{2}).\end{eqnarray}$$

Note that the leading-order term of $G(y)$ is non-zero which is responsible for generating convection cells when $Rc=0$ , but it cannot generate axial flow in the channel. When (5.1) and (5.2) are substituted into (3.6a,c ), we obtain

(5.3a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\displaystyle \frac{\text{d}^{4}F_{0}}{\text{d}y^{4}}}-2K^{2}{\displaystyle \frac{\text{d}^{2}F_{0}}{\text{d}y^{2}}}+K^{4}F_{0}=K^{2}\text{csch}(K)\sinh [K(1-y)],\end{eqnarray}$$

(5.3b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\displaystyle \frac{\text{d}^{4}F_{1}}{\text{d}y^{4}}}-2K^{2}{\displaystyle \frac{\text{d}^{2}F_{1}}{\text{d}y^{2}}}+K^{4}F_{1}=-\text{i}K{\displaystyle \frac{\text{d}^{2}F_{0}}{\text{d}y^{2}}}+\text{i}K^{3}F_{0}\nonumber\\ \displaystyle & & \displaystyle \quad \displaystyle +\,\frac{1}{4}\text{i}K^{2}Pr\,\text{csch}^{2}(K)\{(y-2)\sinh (Ky)-y\sinh [K(y-2)]\},\end{eqnarray}$$

(5.3c ) $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle F_{i}(0)=F_{i}(1)={\displaystyle \frac{\text{d}F_{i}}{\text{d}y}}(0)={\displaystyle \frac{\text{d}F_{i}}{\text{d}y}}(1)=0\quad \text{for}~i=0,1.\end{eqnarray}$$

Since the source term in (5.3a ) is real, $F_{0}$ has a real solution. Thus, $F_{1}$ has a pure imaginary solution, according to (5.3b ). These solutions are readily obtained with any computer algebra system, and omitted here due to their length. Substitution of $F$ in (5.1) into (3.15) gives

(5.4) $$\begin{eqnarray}\displaystyle & \displaystyle \overline{u_{0}v_{0}}={\displaystyle \frac{Rc}{2\text{i}K}}\left({\displaystyle \frac{\text{d}F_{0}}{\text{d}y}}F_{1}-F_{0}{\displaystyle \frac{\text{d}F_{1}}{\text{d}y}}\right)+O(Rc^{2}). & \displaystyle\end{eqnarray}$$

The first-order volume flow rate is then calculated from (3.14):

(5.5) $$\begin{eqnarray}\displaystyle \displaystyle Q_{1} & = & \displaystyle -{\displaystyle \frac{(1+Pr)Rc}{645120K^{7}\left(\text{e}^{2K}-1\right)^{2}\left[(\text{e}^{2K}-1)^{2}-4K^{2}\text{e}^{2K}\right]}}\left\{315\left[\text{e}^{8K}(K-5)+K+5\right]\right.\nonumber\\ \displaystyle & & \displaystyle +\,2\text{e}^{6K}\left(12K^{7}+70K^{6}-420K^{5}+2520K^{4}-3360K^{3}+3150K^{2}+2520K+1575\right)\nonumber\\ \displaystyle & & \displaystyle +\,2\text{e}^{4K}\left(8K^{9}+256K^{7}+1680K^{5}-5880K^{3}-5355K\right)+2\text{e}^{2K}\left(12K^{7}-70K^{6}\right.\nonumber\\ \displaystyle & & \displaystyle \left.\!\left.-\,420K^{5}-2520K^{4}-3360K^{3}-3150K^{2}+2520K-1575\right)\!\right\}.\end{eqnarray}$$

This is plotted in figure 7 and agrees with the full asymptotic solution when $Rc$ is small.

Figure 7.

First-order asymptotic flow rate $Q_{1}$ versus $Rc$ for $K=2\unicode[STIX]{x03C0}$ and various $Pr$ . The solid lines are the full solution of $Q_{1}$ calculated by the method described in appendix B. The dashed lines are the small and large $Rc$ asymptotic solutions (5.5) and (5.14), respectively.

5.2.2 Asymptotic solution of $Q_{1}$ in the limit $Rc\rightarrow \infty$

Equation (3.6a ) is singular in this limit. The solution is expected to have a boundary layer at the bottom wall because as the wave speed increases, the heating and cooling by the wave at the bottom wall will cancel away from the wall. A dominant balance between the first and second terms of (3.6a ) suggests that the boundary-layer thickness scales as $1/\sqrt{Rc}$ . Thus, we define a small parameter $\unicode[STIX]{x1D716}=1/\sqrt{Rc}$ , and (3.6a ) becomes

(5.6) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D716}^{2}{\displaystyle \frac{\text{d}^{4}F}{\text{d}y^{4}}}+K(\text{i}-2K\unicode[STIX]{x1D716}^{2}){\displaystyle \frac{\text{d}^{2}F}{\text{d}y^{2}}}-K^{3}(\text{i}-K\unicode[STIX]{x1D716}^{2})F=\unicode[STIX]{x1D716}^{2}K^{2}G(y). & & \displaystyle\end{eqnarray}$$

Equation (5.6) is singular in the limit $\unicode[STIX]{x1D716}\rightarrow 0$ , and is solved by the method of matched asymptotic expansions.

In the outer region, $G(y)$ expands to $\text{e}^{-(1-\text{i})\sqrt{KPr/2}\,y/\unicode[STIX]{x1D716}}$ , which is exponentially small in the limit $\unicode[STIX]{x1D716}\rightarrow 0$ and can be safely dropped. Substituting $F=F_{0}(y)+\unicode[STIX]{x1D716}F_{1}(y)+\unicode[STIX]{x1D716}^{2}F_{2}(y)+\unicode[STIX]{x1D716}^{3}F_{3}(y)+\unicode[STIX]{x1D716}^{4}F_{4}(y)+\cdots \,$ into (5.6), we find

(5.7a ) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\displaystyle \frac{\text{d}^{2}F_{0}}{\text{d}y^{2}}}-K^{2}F_{0}=0,\\[10.0pt] \displaystyle {\displaystyle \frac{\text{d}^{2}F_{1}}{\text{d}y^{2}}}-K^{2}F_{1}=0,\\[7.0pt] \displaystyle {\displaystyle \frac{\text{d}^{2}F_{2}}{\text{d}y^{2}}}-K^{2}F_{2}={\displaystyle \frac{\text{i}}{K}}\frac{\text{d}^{4}F_{0}}{\text{d}y^{4}}-2\text{i}K\frac{\text{d}^{2}F_{0}}{\text{d}y^{2}}+\text{i}K^{3}F_{0},\\[10.0pt] \displaystyle {\displaystyle \frac{\text{d}^{2}F_{3}}{\text{d}y^{2}}}-K^{2}F_{3}={\displaystyle \frac{\text{i}}{K}}\frac{\text{d}^{4}F_{1}}{\text{d}y^{4}}-2\text{i}K\frac{\text{d}^{2}F_{1}}{\text{d}y^{2}}+\text{i}K^{3}F_{1},\\[10.0pt] \displaystyle {\displaystyle \frac{\text{d}^{2}F_{4}}{\text{d}y^{2}}}-K^{2}F_{4}={\displaystyle \frac{\text{i}}{K}}\frac{\text{d}^{4}F_{2}}{\text{d}y^{4}}-2\text{i}K\frac{\text{d}^{2}F_{2}}{\text{d}y^{2}}+\text{i}K^{3}F_{2}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Thus, the fourth-order ordinary differential equation (5.6) has been reduced by two orders, and only two boundary conditions are needed. We apply the no-penetration condition at the top wall:

(5.7b ) $$\begin{eqnarray}\displaystyle \displaystyle F_{i}(1)=0\quad \text{for}~i=0,1,2,3,\ldots . & & \displaystyle\end{eqnarray}$$

The other boundary condition will come from matching with the inner solution. We find $F_{0}=F_{1}=F_{2}=F_{3}=0$ because their corresponding inner solutions are zero (see below). Thus,

(5.8) $$\begin{eqnarray}\displaystyle & \displaystyle F_{4}=C\text{e}^{-ky}\left(\text{e}^{2ky}-\text{e}^{2k}\right), & \displaystyle\end{eqnarray}$$

where $C$ is an integration constant to be determined by matching.

In the inner region, we define $\unicode[STIX]{x1D702}=y/\unicode[STIX]{x1D716}$ and $f(\unicode[STIX]{x1D702})=F(y)$ , and (5.6) becomes

(5.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\text{d}^{4}f}{\text{d}\unicode[STIX]{x1D702}^{4}}}+K\left(\text{i}-2K\unicode[STIX]{x1D716}^{2}\right){\displaystyle \frac{\text{d}^{2}f}{\text{d}\unicode[STIX]{x1D702}^{2}}}-\unicode[STIX]{x1D716}^{2}K^{3}\left(\text{i}-K\unicode[STIX]{x1D716}^{2}\right)f=K^{2}\unicode[STIX]{x1D716}^{4}\text{e}^{-(1-\text{i})\sqrt{KPr/2}\,\unicode[STIX]{x1D702}}. & \displaystyle\end{eqnarray}$$

The inner dependent variable is expanded as $f=f_{0}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}f_{1}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}^{2}f_{2}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}^{3}f_{3}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}^{4}f_{4}(\unicode[STIX]{x1D702})+\cdots \,$ . Since the source term ${\sim}\unicode[STIX]{x1D716}^{4}$ , $f_{0}=f_{1}=f_{2}=f_{3}=0$ and $f_{4}$ obeys

(5.10) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\text{d}^{4}f_{4}}{\text{d}\unicode[STIX]{x1D702}^{4}}}+\text{i}K{\displaystyle \frac{\text{d}^{2}f_{4}}{\text{d}\unicode[STIX]{x1D702}^{2}}}=K^{2}\text{e}^{-(1-\text{i})\sqrt{KPr/2}\,\unicode[STIX]{x1D702}}, & & \displaystyle\end{eqnarray}$$

which is subject to the boundary conditions $f_{4}(0)=\text{d}f_{4}/\text{d}y(0)=0$ , and the requirement that $f_{4}(\unicode[STIX]{x1D702}\rightarrow \infty )$ remains finite. These conditions are sufficient to determine all the integration constants, and we have

(5.11) $$\begin{eqnarray}\displaystyle & \displaystyle f_{4}(\unicode[STIX]{x1D702})={\displaystyle \frac{1}{(\sqrt{Pr}+1)Pr}}\left[{\displaystyle \frac{\sqrt{Pr}\,\text{e}^{-(1-\text{i})\sqrt{K/2}\,\unicode[STIX]{x1D702}}-\text{e}^{-(1-\text{i})\sqrt{KPr/2}\,\unicode[STIX]{x1D702}}}{\sqrt{Pr}-1}}-1\right]. & \displaystyle\end{eqnarray}$$

From the asymptotic matching principle, $\lim _{y\rightarrow 0}F=\lim _{\unicode[STIX]{x1D702}\rightarrow \infty }f\equiv F_{m}$ , we get

(5.12a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle C={\displaystyle \frac{1}{(\sqrt{Pr}+1)Pr}}\cdot {\displaystyle \frac{1}{(\text{e}^{2K}-1)}},\quad F_{m}=-{\displaystyle \frac{\unicode[STIX]{x1D716}^{4}}{(\sqrt{Pr}+1)Pr}}. & \displaystyle\end{eqnarray}$$

The composite solution is

(5.13) $$\begin{eqnarray}\displaystyle & \displaystyle F_{c}(y)=f(y/\unicode[STIX]{x1D716})+F(y)-F_{m}, & \displaystyle\end{eqnarray}$$

which is separated into real and imaginary parts and substituted into (3.15) to give the Reynolds stress. The leading-order volume flow rate is found from (3.14) as

(5.14) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{1}=-{\displaystyle \frac{1}{8KPr^{2}(1+Pr)}}Rc^{-4}+O(Rc^{-9/2}). & \displaystyle\end{eqnarray}$$

This flow rate comes from the inner region; the outer and matching regions make no contribution to $Q_{1}$ at the leading order (appendix C.1). Equation (5.14) is plotted in figure 7. Since the error term in (5.14) is $O(Rc^{-1/2})$ compared with the leading-order term, the asymptotic solution agrees with the full analytic solution only at extremely large $Rc$ in figure 7. Note that $Q_{1}$ varies differently with $Pr$ in the small and large $Rc$ limits; $|Q_{1}|$ increases with $Pr$ when $Rc$ is small, but decreases with $Pr$ when $Rc$ is large. This is why the curves in figure 7 cross and switch order as $Rc$ increases. We will further discuss this in the next section.

While the temperature $G(y)$ is confined within the boundary layer (exponentially small in the outer region), the vertical velocity $F(y)$ extends to the outer region. This is confirmed by a numerical simulation using a relatively large $Rc$ value of $500$ (appendix C.3). Because of the non-zero outer solution, a top wall boundary layer is needed to satisfy the no-slip condition. This boundary layer has negligible contribution to the net flow rate, and is presented in appendix C.2.

Figure 8.

First-order volume flow rate $Q_{1}$ versus wavenumber $K$ for $Rc=1$ and various $Pr$ . The solid curves are the full solution of $Q_{1}$ calculated by the method described in appendix B. The dashed lines are the small and large $K$ asymptotic solutions (5.15) and (5.16), respectively. The open circles are the small $Rc$ asymptotic solution (5.5).

5.3 Wavenumber K

In figure 8, $Q_{1}$ is plotted as a function of $K$ for $Rc=1$ and various $Pr$ . It shows that $|Q_{1}|$ first increases with $K$ until $K\approx 3$ before decreasing. The behaviour of $Q_{1}$ as $K\rightarrow 0$ or $\infty$ makes physical sense. As $K\rightarrow 0$ , the thermal wavelength goes to infinity. In this limit, the thermal boundary condition at the bottom wall reduces to zero as shown by (2.6e ). Thus, there is no driving force and the axial volume flow rate vanishes. In the limit $K\rightarrow \infty$ , the thermal wavelength goes to zero, and the close proximity of temperature difference cancels itself. Thus, the driving force and therefore the axial volume flow rate again reduce to zero. The asymptotic analysis in the small and large limits of $K$ are similar to that of $Rc$ , so we only discuss the main points, leaving the mathematical details to appendix D.

In the limit $K\rightarrow 0$ , the expansion of $F(y)$ in (3.6a ) is regular, and the axial volume flow rate is found to be

(5.15) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{1}=-{\displaystyle \frac{229(1+Pr)RcK^{4}}{435\,891\,456\,000}}+O(K^{5}). & \displaystyle\end{eqnarray}$$

This is plotted in figure 8 and agrees well with the full solution of $Q_{1}$ for $K\ll 1$ .

In the limit $K\rightarrow \infty$ , the expansion of $F(y)$ is singular. There is a boundary layer of thickness $O(K^{-1})$ near $y=0$ . The temperature variation and buoyancy driven flow are confined within this boundary layer. Outside the boundary layer, the temperature is uniform and the fluid is stationary. The axial volume flow rate is

(5.16) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{1}=-{\displaystyle \frac{(1+Pr)Rc}{2048K^{6}}}+O(K^{-7}). & \displaystyle\end{eqnarray}$$

This is plotted in figure 8 and agrees well with the full solution of $Q_{1}$ for $K\gg 1$ . Since the boundary-layer thickness is $O(K^{-1})$ and the axial volume flow within the boundary layer is $Q_{1}\sim K^{-6}$ , it would be difficult for a numerical method to resolve this boundary layer, and to accurately calculate $Q$ .

5.4 Prandtl number $Pr$

Figure 9 shows the dependence of $Q_{1}$ on $Pr$ . Since $Pr$ is a physical property of the fluid, it cannot be varied widely, and therefore the limit of $Pr\rightarrow 0$ or $\infty$ has little practical value. The limits also do not yield simpler equations. For these reasons, we do not seek asymptotic solutions and only briefly discuss the limiting behaviour.

The limit $Pr\rightarrow 0$ differs from the small $K$ or $Rc$ limit in that $Q$ approaches a non-zero constant, i.e. an axial flow exists for infinite thermal diffusivity (figure 9). The repetitive appearance of the relation $Q_{1}\propto (1+Pr)$ in previous sections suggests its wide applicability. Indeed, $Q_{1}=Q_{1}(Pr=0)~(1+Pr)$ fits the numerical solution very well, where the proportionality constant $Q_{1}(Pr=0)$ is the axial flow rate at $Pr=0$ and is calculated numerically. We note that for the curves with $Rc=0.1$ and $1$ in figure 9, $Q_{1}(Pr=0)$ can be calculated by the analytic solution in (5.5); for the curve with $Rc=10$ , equation (5.5) is not accurate. The large $Pr$ limit implies infinite viscosity or zero thermal diffusivity, and $Q_{1}$ vanishes. The scaling $Q_{1}\sim Pr^{-3}$ taken from the large $Rc$ asymptotic solution (5.14) captures the dependence quite well, although the proportionality constant is not straightforward to derive.

In figure 9, the curves with different $Rc$ cross as $Pr$ increases. This implies different dependence of $Q_{1}$ on $Rc$ as $Pr$ varies from infinitesimal to infinity. Figure 9 appears to suggest that $Q_{1}$ is proportional to $Rc$ as $Pr\rightarrow 0$ , and to $Rc^{-3}$ as $Pr\rightarrow \infty$ . A detailed analysis would be needed to validate this observation.

Figure 9 also shows that $|Q_{1}|$ reaches a maximum before decreasing as $Pr$ increases for fixed $K$ and $Rc$ . Similar behaviour of $Q_{1}$ is observed in figures 7 and 8 as $Rc$ and $K$ vary. We will take a close look at this in § 6.

Figure 9.

First-order asymptotic volume flow rate $Q_{1}$ versus Prandtl number for $K=2\unicode[STIX]{x03C0}$ and various $Rc$ . The solid curves are the full solution of $Q_{1}$ calculated by the method described in appendix B. The dashed lines are the scaling $1+Pr$ as $Pr\rightarrow 0$ and $Pr^{-3}$ as $Pr\rightarrow \infty$ . All curves shown exhibit a maximum in $|Q_{1}|$ .

5.5 Discussion and comparison with previous work

Equations (5.15) and (5.16) describe how the axial flow rate $Q_{1}$ behaves in the limit the thermal wavenumber $K\rightarrow 0$ and $\infty$ , respectively. In each equation, the scaling with $K$ shows that the dependence is quite sensitive. It is interesting that $Q_{1}\propto (1+Pr)Rc$ in both small and large $K$ limits. For example, the large $K$ or ‘short wave limit’ corresponds to a thermal wave travelling beneath an infinitely deep channel. The axial flow rate is linearly proportional to how fast the thermal wave moves.

The asymptotic solutions (5.15) and (5.16) agree with the asymptotic solution (5.5) obtained in the limit $Rc\rightarrow 0$ when it is further expanded in the limits $K\rightarrow 0$ and $\infty$ . This is surprising because in deriving (5.5), $Rc$ is assumed to be small in the beginning, while in (5.15) and (5.16) no assumption is made about the magnitude of $Rc$ . In fact, a plot of (5.5) in figure 8 shows that it captures the variation of $Q$ for the entire range of $K$ when $Pr$ is small (0.1 and 1). For the larger $Pr=10$ value, it fails at intermediate $K$ values, where the large and small $K$ asymptotic solutions also do not hold.

While we could derive (5.15) and (5.16) from (5.5), we could not reverse the direction and derive the small and large $Rc$ asymptotic solutions from the small $K$ solution (5.15) by expanding it in the small and large $Rc$ limits, respectively. Mathematically this states the non-trivial relation among the different limits of parameters; while the limit $Rc\rightarrow 0$ can be swapped either with the limit $K\rightarrow 0$ or $K\rightarrow \infty$ , the limit $Rc\rightarrow \infty$ can neither be swapped with the limit $K\rightarrow 0$ nor with $K\rightarrow \infty$ .

The dependence of $Q$ on the thermal wave velocity has been well studied in the long wave limit (Stern Reference Stern1959; Davey Reference Davey1967). After translating their dimensionless groups to ours, the net volume flow rate given by Stern (Reference Stern1959) is

(5.17) $$\begin{eqnarray}\displaystyle Q_{1}=\left\{\begin{array}{@{}ll@{}}\displaystyle -5\times 10^{-7}RcK^{4} & KRc\,\text{small},\\[10.0pt] \displaystyle {\displaystyle \frac{1}{48\sqrt{2}}}Rc^{-5/2}K^{-1/2} & KRc\,\text{large}.\end{array}\right. & & \displaystyle\end{eqnarray}$$

The second expression says that $Q$ reverses its direction to be the same as the thermal wave, provided that $KRc$ is large enough (the threshold was estimated to be $10^{2}$ ). Instead of assuming a volumetric temperature field, Davey (Reference Davey1967) used a more realistic approach and obtained $T(x,y)$ as a solution of the heat diffusion equation. Consequently,  $Q_{1}$ was found to be

(5.18) $$\begin{eqnarray}\displaystyle Q_{1}=\left\{\begin{array}{@{}ll@{}}\displaystyle -{\displaystyle \frac{1+Pr}{12!}}RcK^{4} & KPrRc\,\text{small},\\[10.0pt] \displaystyle -{\displaystyle \frac{1}{4Pr^{2}(1+Pr)}}Rc^{-4}K^{-1} & KPrRc\,\text{large}.\end{array}\right. & & \displaystyle\end{eqnarray}$$

The dependence of $Q_{1}$ on $Rc$ and $K$ in the small $Rc$ limit differs from that of Stern (Reference Stern1959) only in the numerical coefficient. The dependence on $Pr$ is interesting as it shows that $Pr=0$ is a regular limit, which explains why the (volumetric temperature field) assumption made by Stern (Reference Stern1959) works. The large $Rc$ limit of $Q_{1}$ flows opposite to the thermal wave direction, and is proportional to $Rc^{-4}$ instead of $Rc^{-5/2}$ . These differences highlight how the initial simplification could qualitatively affect the result.

We do not make any initial simplification on $Rc$ , $K$ and $Pr$ . Our result has the same $(1+Pr)Rc$ dependence as the one obtained by Davey (Reference Davey1967). In the limit $K\rightarrow 0$ , equation (5.5) reduces to

(5.19) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{1}=-{\displaystyle \frac{229(1+Pr)RcK^{4}}{435\,891\,456\,000}}. & \displaystyle\end{eqnarray}$$

Equation (5.19) results from taking two limits, first in $Rc\rightarrow 0$ and followed by $K\rightarrow 0$ . It is similar to the one by Davey (Reference Davey1967), with about $1/4$ the numerical coefficient value. This means that in the small $Rc$ and small $K$ limits, applying thermal waves at both walls generate four times greater axial flow rate, compared with applying the thermal wave at only one wall.

If we take $K\rightarrow \infty$ , equation (5.5) becomes

(5.20) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{1}=-{\displaystyle \frac{(1+Pr)Rc}{2048K^{6}}}, & \displaystyle\end{eqnarray}$$

which differs from the large $KRc$ and large $KPrRc$ limits obtained by Stern (Reference Stern1959) and Davey (Reference Davey1967), respectively. This is expected as $K\rightarrow \infty$ is a ‘short wave’ limit, and is incompatible with the long wave assumption made by these authors. It may also be argued that the double limits of $Rc\rightarrow 0$ followed by $K\rightarrow \infty$ do not define the limit of their product $KRc$ . However, this is not a valid argument because (5.16), which is the same as (5.20), is derived in the limit $K\rightarrow \infty$ without any assumption on the magnitude of $Rc$ .

The large $Rc$ limit (5.14) is exactly half that of Davey (Reference Davey1967) reproduced above. This is because we apply the thermal wave only at the bottom wall, while in Davey (Reference Davey1967) the thermal wave is applied at both bottom and top walls. It is interesting that in the large $Rc$ limit, applying the thermal wave at both walls only double the axial flow rate, apparently due to all thermal effects being confined within the boundary layer.