Temperature-Coupled Similarity Reductions for Unsteady Two-Dimensional Incompressible Navier–Stokes Equations

1. Introduction

Laminar forced convection is described, in almost every textbook and every engineering correlation, by the pair consisting of the Blasius momentum solution [2] and the Pohlhausen energy solution [9]. Both are steady. When the thermal problem is genuinely time dependent—a plate switched on, a wall temperature stepped, a start-up transient—the standard procedure retains $\partial T/\partial t$ in the energy equation while inserting the steady velocity field.

That procedure is an approximation, and often a good one: when the thermal layer is much thinner than the momentum layer, as at large Prandtl number, the velocity field evolves slowly on the scale that matters for T and the quasi-steady treatment is entirely adequate for engineering purposes. It is nevertheless not self-consistent—a fluid particle cannot convect a time-dependent temperature along a time-independent trajectory—and its accuracy has never been assessed against an exact coupled solution, because none was available. Our object is to supply one: not to displace the quasi-steady approximation, but to give it a benchmark. Exact solutions of the Navier–Stokes equations are rare [16,17,18], and unsteady two-dimensional ones rarer still; the classical analyses of Stokes [13], Rayleigh [10], Stewartson [12], Hall [5], Dennis [3] and Riley [11] produced asymptotic and numerical results, and the group classification of Ma and Hui [8] enumerated the similarity reductions of the boundary-layer system, but the thermal problem inherited the gap of necessity: an exact unsteady temperature field presupposes an exact unsteady velocity field to convect it.

Two developments make the coupled problem tractable. The first is the diffusion-time-scale transformation of Sun [15], which promotes time itself to a similarity variable by measuring it against the transverse diffusion time ${\delta }^2/\nu$ identified by Stuart [14]. The second is the observation, developed in a companion analysis, that the same transformation applied to the full equations—with no thin-layer approximation—closes only for particular members of the family, and that the stretching operator it generates annihilates a distinguished variable which reduces the problem to an ordinary differential equation.

The present paper carries that programme through for the temperature field. We work throughout with the complete equations: streamwise diffusion $\nu u_{xx}$, streamwise conduction $\chi T_{xx}$ and the full dissipation function are all retained, so that every statement is a statement about the Navier–Stokes system rather than about a truncation of it. Section 2 sets out the equations and §3 transforms the energy equation. Section 4 establishes the closure condition for the coupled system, which is narrower than for the momentum equation alone, and identifies dissipation as the obstruction. Section 5 obtains the exact solutions, §6 the heat-transfer characteristics, and §7 solves the stagnation-point member numerically. Section 8 concludes.

2. Formulation

A Newtonian fluid of constant density $\rho$, kinematic viscosity $\nu$, thermal diffusivity $\chi =k/\rho c_p$ and Prandtl number $\mathit{Pr}=\nu /\chi$ occupies a two-dimensional region bounded by a wall. The governing equations, complete, are

$$\begin{array}{c}{u}_x+v_y=0,\end{array}$$

(1)

$$\begin{array}{c}{u}_t+u{u}_x+v{u}_y=-{\displaystyle \frac{p_x}{\rho }}+\nu \left(u_{xx}+u_{yy}\right),\end{array}$$

(2)

$$\begin{array}{c}{v}_t+u{v}_x+v{v}_y=-{\displaystyle \frac{p_y}{\rho }}+\nu \left(v_{xx}+v_{yy}\right),\end{array}$$

(3)

$$\begin{array}{c}{T}_t+u{T}_x+v{T}_y=\chi \left(T_{xx}+T_{yy}\right)+{\displaystyle \frac{\mathsf{\Phi }}{\rho c_p}}.\end{array}$$

(4)

The dissipation function is $\mathsf{\Phi }=2\mu e_{ij}{e}_{ij}$ with $e_{ij}=\frac{1}{2}(u_{i,j}+u_{j,i})$. In two dimensions, using $u_x=-v_y$,

$${\displaystyle \frac{\mathsf{\Phi }}{\rho c_p}}={\displaystyle \frac{\nu }{c_p}}\left[4u_x^2+{\left(u_y+v_x\right)}^2\right],$$

(5)

of which the thin-layer approximation retains only $u_y^2$. We keep both terms; the first will turn out to control the closure of the coupled problem.

Continuity is satisfied by $u={\psi }_y$, $v=-{\psi }_x$. The wall is impermeable and at rest, and we consider the isothermal and adiabatic conditions

$$T(x,0,t)=T_w\left(x\right),\mathrm{or}{\left.T_y\right|}_{y=0}=0,$$

(6)

with $T\to T_{\infty }$ far from the wall.

3. Transformation of the Energy Equation

3.1. Variables and the Stretching Operator

Following Sun [15] we set

$$\psi =U\left(x\right)\delta \left(x\right)f(\eta ,\tau ),\eta ={\displaystyle \frac{y}{\delta \left(x\right)}},\tau ={\displaystyle \frac{\nu t}{\delta {\left(x\right)}^2}}.$$

(7)

Since $\eta \propto {\delta }^{-1}$ and $\tau \propto {\delta }^{-2}$, all streamwise chain rules collapse onto the stretching operator

$$\mathcal{D}\equiv \eta {\displaystyle \frac{\partial }{\partial \eta }}+2\tau {\displaystyle \frac{\partial }{\partial \tau }},$$

(8)

and for any $F(\eta ,\tau )$ and any $G\left(x\right)\propto x^q$,

$$\begin{array}{cc}\hfill & {\displaystyle \frac{\partial F}{\partial x}}=-{\displaystyle \frac{{\delta }^{\prime }}{\delta }}\mathcal{D}F,{\displaystyle \frac{\partial }{\partial x}}\left[GF\right]={\displaystyle \frac{G}{x}}\left(q-n\mathcal{D}\right)F,\hfill \\ \hfill & {\displaystyle \frac{{\partial }^2}{\partial x^2}}\left[GF\right]={\displaystyle \frac{G}{x^2}}\left[(q-1)-n\mathcal{D}\right]\left(q-n\mathcal{D}\right)F,\hfill \end{array}$$

(9)

with $n=x{\delta }^{\prime }/\delta$. The velocity components are

$$u=U{f}_{\eta },v={\displaystyle \frac{\nu }{\delta }}V,V={\displaystyle \frac{\gamma }{2}}\mathcal{D}f-\alpha f,$$

(10)

where, as in the momentum problem,

$$\alpha ={\displaystyle \frac{\delta }{\nu }}{\displaystyle \frac{\mathrm{d}\left(U\delta \right)}{\mathrm{d}x}},\beta ={\displaystyle \frac{{\delta }^2}{\nu }}{\displaystyle \frac{\mathrm{d}U}{\mathrm{d}x}},\gamma ={\displaystyle \frac{U}{\nu }}{\displaystyle \frac{\mathrm{d}{\delta }^2}{\mathrm{d}x}}=2(\alpha -\beta ),$$

(11)

and constancy of these forces the power-law family

$$U=C{x}^m,m={\displaystyle \frac{\beta }{2\alpha -\beta }},{\delta }^2={\displaystyle \frac{\nu (2\alpha -\beta )x}{U}}.$$

(12)

Three derived quantities recur below:

$$n\equiv {\displaystyle \frac{x{\delta }^{\prime }}{\delta }}={\displaystyle \frac{1-m}{2}},\sigma \equiv {\displaystyle \frac{1}{n}}=1+{\displaystyle \frac{2\alpha }{\gamma }},\mathit{Re}\equiv {\displaystyle \frac{U\left(x\right)x}{\nu }},$$

(13)

$\mathit{Re}$ being the local Reynolds number, together with $\epsilon \equiv {\delta }^{\prime 2}=n^2(2\alpha -\beta )/\mathit{Re}$, which measures the size of the streamwise-diffusion and streamwise-conduction terms.

3.2. The Wall Temperature Scale

Write

$$T=T_{\infty }+\Delta T\left(x\right)\theta (\eta ,\tau ),\Delta T=D{x}^k,\mathit{Ec}={\displaystyle \frac{U^2}{c_p\Delta T}}.$$

(14)

The Eckert number is constant if and only if $\Delta T\propto U^2$, that is

$$k=2m.$$

(15)

This is a genuine extension of the isothermal case: holding $T_w-T_{\infty }$ fixed forces $\mathit{Ec}\propto x^{2m}$ and restricts the reduction to the flat plate, whereas (15) admits the whole family. For $m=0$ it reduces to a uniform wall temperature.

3.3. The Cancellation

With T of the form (14) the two convective terms are, using (9),

$$u{T}_x={\displaystyle \frac{\Delta TU}{x}}{f}_{\eta }\left[k\theta -n\mathcal{D}\theta \right],v{T}_y={\displaystyle \frac{\nu \Delta T}{{\delta }^2}}V{\theta }_{\eta }.$$

(16)

Both contain $\eta f_{\eta }{\theta }_{\eta }$: the first through $-n\mathcal{D}\theta$, the second through the $\mathcal{D}f$ inside V. Since $n(2\alpha -\beta )=\gamma /2$ the two contributions are $\mp (\gamma /2)\eta f_{\eta }{\theta }_{\eta }$ and cancel identically. This is the same cancellation that makes the momentum equation reduce, and it is the step most easily lost: differentiating f with respect to x through $\eta$ only, and not through $\tau$, leaves a spurious $\frac{1}{2}\gamma \eta f_{\eta }{\theta }_{\eta }$ and destroys every reduction below.

3.4. The Transformed Equation

Multiplying (4) by $\mathit{Pr}{\delta }^2/\nu \Delta T$ and collecting gives the exact transformed energy equation

$$\begin{array}{cc}\hfill & \mathit{Pr}\left[{\theta }_{\tau }-\alpha f{\theta }_{\eta }+2\beta f_{\eta }\theta +\gamma \tau \left(f_{\tau }{\theta }_{\eta }-f_{\eta }{\theta }_{\tau }\right)\right]\hfill \\ \hfill & ={\theta }_{\eta \eta }+\epsilon \left(\mathcal{D}+4-\sigma \right)\left(\mathcal{D}+4-2\sigma \right)\theta +\mathit{Pr}\mathit{Ec}\mathfrak{D}.\hfill \end{array}$$

(17)

with the transformed dissipation

$$\mathfrak{D}={\left[f_{\eta \eta }-{\displaystyle \frac{n}{\mathit{Re}}}\left(1+\mathcal{D}\right)V\right]}^2+{\displaystyle \frac{4(2\alpha -\beta )}{\mathit{Re}}}{\left[\left(m-n\mathcal{D}\right)f_{\eta }\right]}^2.$$

(18)

The first group on the right of (17) is transverse conduction, the second streamwise conduction $\chi T_{xx}$, the third dissipation. The term $2\mathit{Pr}\beta f_{\eta }\theta$ is contributed by the varying wall temperature and vanishes for the flat plate. Setting $\epsilon \to 0$, $\beta =0$ and retaining only $f_{\eta \eta }^2$ in (18) recovers the thin-layer form of the transformed energy equation.

Two independent checks were made. At $\mathit{Pr}=1$, $\mathit{Ec}=0$, $\beta =0$ the streamwise velocity is itself a passive scalar, and setting $\theta =f_{\eta }$ in (17) returns the momentum equation exactly; with the spurious term retained it does not. Separately, (17) was verified against (4) in physical $(x,y,t)$ variables for $m=1$, using none of the algebra above (Appendix A).

4. Closure of the Coupled System

The parameters $\alpha ,\beta ,\gamma ,\sigma ,n,m$ are constants by construction. The reduction to two independent variables therefore requires that the remaining coefficients in (17)–(18) be free of x. There are three, and since $\mathit{Re}\propto x^{m+1}$ they are

$$\epsilon ={\displaystyle \frac{n^2(2\alpha -\beta )}{\mathit{Re}}},{\displaystyle \frac{n}{\mathit{Re}}},{\displaystyle \frac{2\alpha -\beta }{\mathit{Re}}},$$

(19)

all proportional to $x^{-(m+1)}$. The first two carry a factor n and therefore vanish identically when $n=0$, that is $m=1$; the third does not.

Theorem 1.

Under the transformation (7) with $\Delta T\propto x^{2m}$, the coupled system (1)–(4) reduces to two independent variables if and only if

1.

$m=-1$, when viscous dissipation is retained; or

2.

$m=1$ or $m=-1$, when dissipation is absent $(\mathit{Ec}=0)$.

Proof. 

Immediate from (19): $\epsilon$ and $n/\mathit{Re}$ are free of x if $n=0$ or $\mathit{Re}$ is constant, that is $m=1$ or $m=-1$; the coefficient $(2\alpha -\beta )/\mathit{Re}$ multiplying the streamwise-strain contribution to (18) is free of x only if $\mathit{Re}$ is constant, that is $m=-1$. □

The theorem is sharper than the corresponding statement for the momentum equation alone, which admits both $m=\pm 1$. The obstruction is the term $4u_x^2$ in (5), discarded in every thin-layer treatment: the streamwise strain heats the fluid at a rate that does not scale with the transverse shear, and only a geometry in which $\mathit{Re}$ is constant—a diverging or converging wedge, $U\delta =\mathrm{const}$—keeps the two in proportion.

The case singled out by the theorem deserves comment, since we do not solve it here. For $m=-1$ one has $\alpha =0$, $U\delta =\mathrm{const}$ and $\mathit{Re}$ constant: the flow occupies a wedge whose walls are surfaces of constant $\eta$, and the wall temperature varies as $\Delta T\propto x^{-2}$. The steady limit is Jeffery–Hamel flow [6,7] with heat transfer, and is an ordinary differential system. The unsteady problem is harder, and in a way that is structural rather than technical: with $\epsilon =1/\mathit{Re}$ finite, the streamwise second derivatives make the reduced system elliptic in $(\eta ,\tau )$, with principal symbol ${({\xi }_{\eta }^2+\epsilon s^2)}^2$ where $s=\eta {\xi }_{\eta }+2\tau {\xi }_{\tau }$. The similarity time is then not an evolution variable: data must be posed at both ends of the $\tau$ interval, the line $\tau =0$ is characteristic and admits fewer conditions than an elliptic boundary would, and marching is inconsistent. Whether the temperature field inherits the upstream influence of the velocity field, and how it behaves near the Jeffery–Hamel bifurcation where the steady solution loses uniqueness [1,4], are open questions that we regard as the natural continuation of this work.

Remark 1.

The case $m=0$, the flat plate, closes only in the thin-layer sense: the streamwise conduction and streamwise strain both survive with coefficients $\propto x^{-1}$. It is nevertheless the case for which exact solutions are available, because those solutions annihilate the offending groups individually (§5).

5. Exact Solutions

5.1. The Invariant Variable

Because $\mathcal{D}$ assigns weight 1 to $\eta$ and 2 to $\tau$, the ratio

$$\zeta \equiv {\displaystyle \frac{\eta }{\sqrt{\tau }}}={\displaystyle \frac{y}{\sqrt{\nu t}}}\mathrm{satisfies}\mathcal{D}\zeta =0,$$

(20)

independently of $\delta$ and hence of x. Temperature fields of the form $\theta =\mathsf{\Theta }\left(\zeta \right)$ are therefore group invariant, and three simplifications follow at once. The Jacobian term $\gamma \tau (f_{\tau }{\theta }_{\eta }-f_{\eta }{\theta }_{\tau })$ combines with $-\alpha f{\theta }_{\eta }$ into $-(\gamma /2)f_{\eta }\mathcal{D}\theta$ whenever the velocity field satisfies $V=0$, and then vanishes. The streamwise conduction group reduces to a multiple of $\theta$,

$$\epsilon \left(\mathcal{D}+4-\sigma \right)\left(\mathcal{D}+4-2\sigma \right)\mathsf{\Theta }=\epsilon \left(4-\sigma \right)\left(4-2\sigma \right)\mathsf{\Theta },$$

(21)

which vanishes if $\sigma =2$ or $\sigma =4$, that is $m=0$ or $m=\frac{1}{2}$. And the term $2\mathit{Pr}\beta f_{\eta }\theta$ carries an explicit $\tau$ unless $\beta =0$.

The three conditions are met simultaneously by the branch $\beta =0$, $\alpha =\gamma /2=1$, that is the flat plate, on which the companion momentum analysis gives the exact velocity field

$$f=\sqrt{\tau }G\left(\zeta \right),f_{\eta }=erfc\left({\displaystyle \frac{\zeta }{2}}\right),V=0,\mathcal{D}f=f,$$

(22)

the Rayleigh–Stokes layer, which satisfies the no-slip and outer conditions exactly at every finite $\tau$.

5.2. Reduction of the Energy Equation

Theorem 2.

On the velocity field (22), the transformed energy equation (17) is satisfied by $\theta =\mathsf{\Theta }\left(\zeta \right)$ for arbitrary $\mathit{Pr}$ and $\mathit{Ec}$ if and only if

$${\mathsf{\Theta }}^{\prime \prime }\left(\zeta \right)+{\displaystyle \frac{\mathit{Pr}\zeta }{2}}{\mathsf{\Theta }}^{\prime }\left(\zeta \right)+{\displaystyle \frac{\mathit{Pr}\mathit{Ec}}{\pi }}{e}^{-{\zeta }^2/2}=0.$$

(23)

Proof. 

From $\mathcal{D}f=f$ we have $2\tau f_{\tau }=f-\eta f_{\eta }$, so $-\alpha f{\theta }_{\eta }+\gamma \tau (f_{\tau }{\theta }_{\eta }-f_{\eta }{\theta }_{\tau })=-f_{\eta }\mathcal{D}\theta =0$ on invariant fields, and $\beta =0$ removes the remaining convective term. With $\sigma =2$, (21) vanishes. Since $f_{\eta \eta }=-{\left(\pi \tau \right)}^{-1/2}{e}^{-{\zeta }^2/4}$ and $V=0$, the dissipation (18) is $f_{\eta \eta }^2=e^{-{\zeta }^2/2}/\pi \tau$ up to the $O\left({\mathit{Re}}^{-1}\right)$ strain term, which is $4(2\alpha -\beta ){\mathit{Re}}^{-1}{\left[(m-n\mathcal{D})f_{\eta }\right]}^2=0$ because $m=0$ and $\mathcal{D}{f}_{\eta }=0$. Finally ${\theta }_{\tau }=-(\zeta /2\tau ){\mathsf{\Theta }}^{\prime }$ and ${\theta }_{\eta \eta }={\mathsf{\Theta }}^{\prime \prime }/\tau$; multiplying by $\tau$ gives (23). □

Equation (23) is the unsteady counterpart of the Pohlhausen equation. The correspondence is instructive: in the steady problem the coefficient of ${\mathsf{\Theta }}^{\prime }$ is $\frac{1}{2}\mathit{Pr}{f}_B\left({\eta }_B\right)$, the Blasius function, known only numerically, whereas here it degenerates to the linear function $\frac{1}{2}\mathit{Pr}\zeta$. That degeneracy is what makes the solutions below elementary while their steady analogues are not.

5.3. Solutions

Equation (23) is linear with integrating factor $e^{\mathit{Pr}{\zeta }^2/4}$, so

$${\mathsf{\Theta }}^{\prime }\left(\zeta \right)=e^{-\mathit{Pr}{\zeta }^2/4}\left[{\mathsf{\Theta }}^{\prime }\left(0\right)-{\displaystyle \frac{\mathit{Pr}\mathit{Ec}}{\pi }}{\int }_0^{\zeta }{e}^{(\mathit{Pr}-2)s^2/4}\mathrm{d}s\right].$$

(24)

Imposing $\mathsf{\Theta }\left(0\right)=1$, $\mathsf{\Theta }(\infty )=0$ and interchanging the order of integration gives the wall gradient

$$\begin{array}{cc}\hfill {\mathsf{\Theta }}^{\prime }\left(0\right)& =-\sqrt{{\displaystyle \frac{\mathit{Pr}}{\pi }}}+\mathit{Pr}\mathit{Ec}Q\left(\mathit{Pr}\right),\hfill \\ \hfill Q\left(\mathit{Pr}\right)& ={\displaystyle \frac{1}{\pi }}{\int }_0^{\infty }{e}^{-u^2/2}erfcx\left({\displaystyle \frac{\sqrt{\mathit{Pr}}u}{2}}\right)\mathrm{d}u,\hfill \end{array}$$

(25)

where $erfcx\left(x\right)=e^{x^2}erfc\left(x\right)$; the scaled form is used because the unscaled kernel overflows for $\mathit{Pr}>2$.

For $\mathit{Ec}=0$ the solution is a complementary error function of the group invariant,

$$\theta =erfc\left({\displaystyle \frac{\sqrt{\mathit{Pr}}}{2}}\zeta \right)=erfc\left({\displaystyle \frac{\sqrt{\mathit{Pr}}}{2}}{\displaystyle \frac{y}{\sqrt{\nu t}}}\right)$$

(26)

which satisfies $\theta (0,\tau )=1$ and $\theta (\infty ,\tau )=0$ exactly, at every finite τ. This is worth emphasising. Previous unsteady thermal solutions built on the Kummer-function branch of the velocity field satisfy their boundary conditions only as $\tau \to \infty$, with a defect of $O(1/2\tau )$; here there is no defect, because the underlying velocity field (22) has none.

At $\mathit{Pr}=2$ the integral in (24) is elementary and the full solution with dissipation closes:

$$\mathsf{\Theta }=\left(1-{\displaystyle \frac{2\mathit{Ec}}{\pi }}\right)erfc\left({\displaystyle \frac{\zeta }{\sqrt{2}}}\right)+{\displaystyle \frac{2\mathit{Ec}}{\pi }}{e}^{-{\zeta }^2/2},{\mathsf{\Theta }}^{\prime }\left(0\right)={\displaystyle \frac{\sqrt{2}(2\mathit{Ec}-\pi )}{{\pi }^{3/2}}}.$$

(27)

Profiles for general $\mathit{Pr}$ and $\mathit{Ec}$ are shown in Figure 1.

6. Heat Transfer Characteristics

The wall heat flux follows from ${\theta }_y{|}_0={\mathsf{\Theta }}^{\prime }\left(0\right)/\sqrt{\nu t}$, using $\delta \sqrt{\tau }=\sqrt{\nu t}$, so that with $x/\delta =\sqrt{{\mathit{Re}}_x/2}$,

$$\begin{array}{cc}\hfill {\mathit{Nu}}_x& ={\displaystyle \frac{q_{w}x}{k\Delta T}}={\displaystyle \frac{1}{\sqrt{2\pi }}}\left(1-\frac{1}{2}\mathit{Ec}r\right)\sqrt{{\displaystyle \frac{\mathit{Pr}{\mathit{Re}}_x}{\tau }}}\hfill \\ \hfill & =0.39894\left(1-\frac{1}{2}\mathit{Ec}r\right)\sqrt{{\displaystyle \frac{\mathit{Pr}{\mathit{Re}}_x}{\tau }}},\hfill \end{array}$$

(28)

or equivalently, in purely temporal form,

$$h={\displaystyle \frac{q_w}{\Delta T}}=\sqrt{{\displaystyle \frac{\mathit{Pr}}{\pi }}}\left(1-\frac{1}{2}\mathit{Ec}r\right){\displaystyle \frac{k}{\sqrt{\nu t}}},$$

(29)

which is independent of x: a signature of the Rayleigh-type structure of the unsteady layer, and a qualitative departure from the steady $h\propto x^{-1/2}$ law.

For an adiabatic wall, ${\mathsf{\Theta }}^{\prime }\left(0\right)=0$ in (24) with $\mathsf{\Theta }(\infty )=0$ gives the recovery factor

$$r\equiv {\displaystyle \frac{T_{aw}-T_{\infty }}{U^2/2c_p}}=2\sqrt{\pi \mathit{Pr}}Q\left(\mathit{Pr}\right),{\mathit{Ec}}^{*}={\displaystyle \frac{2}{r}},$$

(30)

${\mathit{Ec}}^{*}$ being the Eckert number at which (28) vanishes and the heat flux reverses. The identity ${\mathit{Ec}}^{*}=2/r$ states that reversal occurs when $T_w=T_{aw}$, as it must. At $\mathit{Pr}=1$ we find $Q=1/2\sqrt{\pi }$ exactly, hence

$$r\left(1\right)=1,{\mathit{Ec}}^{*}\left(1\right)=2,$$

(31)

recovering the classical unit recovery factor. Values are given in Table 1.

6.1. Limits of the Recovery Factor

Both limits of Q are analytic. As $\mathit{Pr}\to 0$ the kernel $erfcx(\sqrt{\mathit{Pr}}u/2)\to 1$ and

$$Q⟶{\displaystyle \frac{1}{\pi }}{\int }_0^{\infty }{e}^{-u^2/2}\mathrm{d}u={\displaystyle \frac{1}{\sqrt{2\pi }}},\mathrm{hence}r⟶\sqrt{2\mathit{Pr}},$$

(32)

so the classical $\sqrt{\mathit{Pr}}$ scaling is recovered with prefactor $\sqrt{2}$ in place of unity. As $\mathit{Pr}\to \infty$ the kernel behaves as $erfcx\left(z\right)\sim 1/z\sqrt{\pi }$, the integral is cut off at $u\sim 2/\sqrt{\mathit{Pr}}$, and the logarithmic divergence of $\int e^{-u^2/2}\mathrm{d}u/u$ gives

$$r\simeq {\displaystyle \frac{2}{\pi }}ln\mathit{Pr}+0.448.$$

(33)

Numerically the fitted slope over ${10}^3\le \mathit{Pr}\le {10}^6$ is $0.63605$ against $2/\pi =0.63662$, and (33) reproduces the computed values to $0.01$ over that range. The logarithmic growth is a genuinely unsteady effect: dissipation is confined to a Stokes sublayer of thickness $\sqrt{\nu t}$ while the thermal layer shrinks as ${\mathit{Pr}}^{-1/2}$, so their overlap saturates instead of following the steady $\sqrt{\mathit{Pr}}$ law. Both limits are shown in Figure 1(c).

Defining $99\%$ thicknesses through $erfc(\zeta /2)=0.01$ and $\theta =0.01$ gives, exactly and at all $\mathit{Pr}$,

$${\displaystyle \frac{{\delta }_T}{{\delta }_u}}={\mathit{Pr}}^{-1/2}.$$

(34)

This is worth contrasting with the steady Prandtl–Blasius layer, in which the corresponding ratio crosses over from ${\mathit{Pr}}^{-1/2}$ at small $\mathit{Pr}$ to ${\mathit{Pr}}^{-1/3}$ at large $\mathit{Pr}$. No crossover occurs here: both layers are diffusive rather than convective in structure, and the single exponent $-\frac{1}{2}$ registers that difference. The exponent itself is not new—the same ${\mathit{Pr}}^{-1/2}$ appears in the thermal analogues of the Rayleigh and Stokes problems [10,13], for the same reason, both layers there being purely diffusive. What is new is that it holds here for a genuinely two-dimensional flow with a nontrivial transverse velocity, where one would expect convection to reinstate the crossover; the reason it does not is that the convective operator is annihilated on $\mathcal{D}$-invariant fields, so the temperature never feels the velocity field except through the dissipation term.

7. The Stagnation-Point Member: A Numerical Example

The case $m=1$ is, by the theorem of §4, an exact reduction of the full system whenever dissipation is absent. Here $\alpha =\beta =1$, $\gamma =0$, $n=0$, $\delta ={(\nu /C)}^{1/2}$ is constant and $\tau =Ct$, the velocity field is $u=Cx{f}_{\eta }$, $v=-C\delta f$ with $u_{xx}=v_{xx}=0$ identically, and the wall temperature varies as $\Delta T\propto x^2$. Equation (17) collapses to

$$\mathit{Pr}\left[{\theta }_{\tau }-f{\theta }_{\eta }+2f_{\eta }\theta \right]={\theta }_{\eta \eta },$$

(35)

coupled to the momentum equation $f_{\eta \eta \eta }+f{f}_{\eta \eta }+1-f_{\eta }^2=f_{\eta \tau }$. Both are parabolic in $\tau$, because $\epsilon \equiv 0$ removes the elliptic contribution, so the pair may be marched.

We integrate from an impulsive start, $f_{\eta }=1$ and $\theta =0$ for $\eta >0$ at $\tau =0$, with $\theta (0,\tau )=1$. Physically this is the stagnation-point analogue of the impulsively started plate [5,12]: at $\tau =0$ the fluid is in the inviscid straining motion $u=Cx$, $v=-Cy$ everywhere and at the ambient temperature, and at $\tau =0^{+}$ the wall is brought to rest and to $T_w$ simultaneously. Both fields therefore begin with a discontinuity at the wall, and both initially diffuse one-dimensionally, which is the origin of the short-time laws below. It differs from the classical problem in that the outer flow here strains as well as convects, so the layer reaches a steady thickness instead of growing indefinitely. The integration uses second-order central differences in $\eta$ and implicit stepping in $\tau$, with Newton iteration for the momentum equation and a direct tridiagonal solve for (35). Results for air are shown in Figure 2.

Two limits are recovered. As $\tau \to 0$,

$$-{\theta }_{\eta }(0,\tau )\sim \sqrt{{\displaystyle \frac{\mathit{Pr}}{\pi \tau }}},$$

(36)

the one-dimensional conduction law, the computed value at $\tau ={10}^{-3}$ being $15.19$ against $15.03$; the convective and wall-temperature terms in (35) have not yet acted. As $\tau \to \infty$ the solution relaxes onto its steady state, with

$$-{\theta }_{\eta }(0,\infty )=0.85650,f_{\eta \eta }(0,\infty )=1.232588,$$

(37)

the latter reproducing the Hiemenz value to six figures and providing an independent check on the coupled solver. The relaxation is complete to three figures by $\tau \simeq 1$, that is after one diffusion time ${\delta }^2/\nu$; the thermal field lags the velocity field slightly, the ratio in Figure 2(b) settling later than either quantity separately.

8. The Fully Closing Case $\mathit{m}=-\mathbf{1}$

Theorem 1 singles out $m=-1$ as the only member for which the coupled system closes with dissipation retained, and it would be unsatisfactory to leave that case untouched. Here $\alpha =0$, $\beta =-1$, $\gamma =2$, $n=\sigma =1$, $U\delta =\mathrm{const}$ and $\mathit{Re}$ is a genuine constant; the walls are surfaces of constant $\eta$, so the flow occupies a wedge of half-angle ${\theta }_w={\eta }_w/{\mathit{Re}}^{1/2}$, and (15) requires a wall temperature $\Delta T\propto x^{-2}$. We take the steady problem, which is an ordinary differential system, and record what it implies for the unsteady one.

8.1. The Steady Thermal Problem

With ${\partial }_{\tau }=0$ the operator $\mathcal{D}$ becomes $\eta \mathrm{d}/\mathrm{d}\eta$, and (17) reduces to

$$\left(1+{\displaystyle \frac{{\eta }^2}{\mathit{Re}}}\right){\theta }^{″}+{\displaystyle \frac{6\eta }{\mathit{Re}}}{\theta }^{\prime }+\left({\displaystyle \frac{6}{\mathit{Re}}}+2\mathit{Pr}{f}^{\prime }\right)\theta +\mathit{Pr}\mathit{Ec}\mathfrak{D}=0,$$

(38)

with

$$\mathfrak{D}={\left[f^{\prime \prime }-{\displaystyle \frac{1}{\mathit{Re}}}\left(2\eta f^{\prime }+{\eta }^2{f}^{\prime \prime }\right)\right]}^2+{\displaystyle \frac{4}{\mathit{Re}}}{\left(f^{\prime }+\eta f^{\prime \prime }\right)}^2,$$

(39)

the second group being the streamwise-strain contribution that Theorem 1 identifies as the obstruction at $m=1$. The velocity field is the Jeffery–Hamel solution [6,7], obtained from $\mathsf{\Omega }=(1+{\eta }^2/\mathit{Re})f^{″}+(2\eta /\mathit{Re})f^{\prime }$ together with $(1+{\eta }^2/\mathit{Re}){\mathsf{\Omega }}^{″}+(6\eta /\mathit{Re}){\mathsf{\Omega }}^{\prime }+(6/\mathit{Re})\mathsf{\Omega }+2f^{\prime }\mathsf{\Omega }=0$. Equation (38) is linear in $\theta$, so it may be solved directly once f is known.

8.2. An Exact Correspondence at $\mathit{Pr}=1$

Comparison of (38) at $\mathit{Pr}=1$, $\mathit{Ec}=0$ with the vorticity equation above shows that the two are identical. The temperature and the vorticity therefore satisfy the same second-order operator, and $\mathsf{\Omega }$ is a known solution of the homogeneous thermal problem. Reduction of order then supplies the general solution by quadrature,

$$\theta =\mathsf{\Omega }\left(\eta \right)\left[A+B{\int }^{\eta }{\displaystyle \frac{\mathrm{d}s}{{\left(\mathit{Re}+s^2\right)}^3\mathsf{\Omega }{\left(s\right)}^2}}\right](\mathit{Pr}=1,\mathit{Ec}=0),$$

(40)

the weight following from $exp(-\int p)$ with $p=6\eta /(\mathit{Re}+{\eta }^2)$. Numerically, imposing $\mathsf{\Omega }$’s own boundary values reproduces $\mathsf{\Omega }$ to $2\times {10}^{-7}$, which is a stringent check on both solvers. Equation (40) is the $m=-1$ analogue of the unsteady Reynolds analogy: at unit Prandtl number the thermal field is slaved to the vorticity field, not to the velocity field, because the wall temperature and the vorticity decay downstream at the same rate. The quadrature is singular where $\mathsf{\Omega }$ vanishes, which for a diverging wedge occurs once inside the channel; the combination remains regular there, but the two branches must be matched across the zero.

8.3. Loss of Uniqueness at Critical Prandtl Numbers

The zeroth-order coefficient in (38) is $6/\mathit{Re}+2\mathit{Pr}{f}^{\prime }$, which is positive whenever the flow is unseparated. Equation (38) is therefore of Helmholtz type rather than of maximum-principle type, and the homogeneous problem $\theta \left(0\right)=\theta \left({\eta }_w\right)=0$ possesses nontrivial solutions at a discrete set of Prandtl numbers. Solving the generalised eigenvalue problem

$$\left(1+{\displaystyle \frac{{\eta }^2}{\mathit{Re}}}\right){\mathsf{\Theta }}^{″}+{\displaystyle \frac{6\eta }{\mathit{Re}}}{\mathsf{\Theta }}^{\prime }+{\displaystyle \frac{6}{\mathit{Re}}}\mathsf{\Theta }=-2{\mathit{Pr}}_c{f}^{\prime }\mathsf{\Theta },\mathsf{\Theta }\left(0\right)=\mathsf{\Theta }\left({\eta }_w\right)=0,$$

(41)

at $\mathit{Re}=10$, ${\eta }_w=2$ gives, grid-converged on 400, 800 and 1600 intervals,

$${\mathit{Pr}}_c=1.8225,9.7846,23.790,43.819,\dots$$

(42)

Below the first of these the boundary-value problem is uniquely solvable; at $\mathit{Pr}={\mathit{Pr}}_c$ it is solvable only if the dissipation term satisfies a compatibility condition, and the wall gradient diverges as $\mathit{Pr}$ approaches ${\mathit{Pr}}_c$ from below, taking the values $0.220$, $-0.524$, $-1.219$ and $-5.150$ at $\mathit{Pr}=0.1$, $0.71$, 1 and $1.5$.

The first critical value is strikingly insensitive to the Reynolds number, moving only from $1.750$ to $1.855$ as $\mathit{Re}$ varies from $2.5$ to 40 (Figure 3c); the higher modes fall slowly with $\mathit{Re}$. This mirrors, in the thermal problem, the non-uniqueness that Jeffery–Hamel flow itself exhibits above a critical Reynolds number [1,4], and it arises from the same feature of the geometry: the streamwise decay $\Delta T\propto x^{-2}$ required by (15) acts as a distributed source in the frame of the moving fluid, and at discrete $\mathit{Pr}$ that source resonates with the diffusive operator. Profiles below the first critical value, with and without dissipation, are shown in Figure 3(a,b).

8.4. The Unsteady Problem

For the unsteady problem at $m=-1$ the coefficient $\epsilon =1/\mathit{Re}$ is finite, so the streamwise second derivatives survive and the reduced system is elliptic in $(\eta ,\tau )$, with principal symbol ${({\xi }_{\eta }^2+\epsilon s^2)}^2$, $s=\eta {\xi }_{\eta }+2\tau {\xi }_{\tau }$, degenerating only on $\tau =0$. The similarity time is then not an evolution variable: data must be posed at both ends of the $\tau$ interval, the line $\tau =0$ is characteristic and admits fewer conditions than an elliptic boundary, and marching in $\tau$ is inconsistent. The temperature field inherits this structure through the same $\epsilon$ groups. Whether the critical Prandtl numbers (42) survive as genuine bifurcations of the unsteady problem, or are smoothed by the ellipticity, is the natural next question; we have not attempted it here.

9. Discussion and Conclusions

Applying the diffusion-time-scale transformation to the complete Navier–Stokes and energy equations gives (17), in which the energy equation carries the same operator $\mathcal{D}$ and the same parameters $\alpha ,\beta ,\gamma$ as the momentum equation, the terms in $\eta f_{\eta }{\theta }_{\eta }$ cancelling identically. A wall temperature $\Delta T\propto x^{2m}$ keeps the Eckert number constant and contributes one further term $2\mathit{Pr}\beta f_{\eta }\theta$.

The principal structural result is the closure theorem of §4. With dissipation retained the coupled system reduces to two variables only for $m=-1$; the obstruction is the streamwise strain $4u_x^2$, discarded in every thin-layer treatment, whose coefficient carries no factor of n and therefore survives at $m=1$. Dropping dissipation restores the second case. This is narrower than the closure condition for the momentum equation alone, and it identifies precisely which neglected term is responsible.

The principal explicit result is the exact solution (26)–(27). Because $\mathcal{D}$ annihilates $\zeta =y/\sqrt{\nu t}$, invariant temperature fields reduce the energy equation to the elementary linear equation (23), whose solution is a complementary error function without dissipation and closes in elementary functions at $\mathit{Pr}=2$. Unlike previous unsteady thermal solutions, it satisfies its boundary conditions exactly at every finite time rather than asymptotically, because the velocity field on which it is built does so. It yields a heat transfer coefficient independent of x, a recovery factor equal to unity at $\mathit{Pr}=1$, and a thickness ratio ${\mathit{Pr}}^{-1/2}$ exact at all Prandtl numbers, with none of the ${\mathit{Pr}}^{-1/3}$ crossover of the steady layer.

We would not present these results as a reason to abandon the quasi-steady treatment, which remains accurate wherever the velocity field evolves slowly on the thermal time scale, and which is far simpler to apply. Their value is as a benchmark: an exact, fully coupled solution of the unsteady problem against which the quasi-steady approximation, and unsteady numerical schemes, can be calibrated. Equation (28) in particular provides a closed-form target with no adjustable constants.

Three limitations should be stated. The solution (26) describes a layer that grows as $\sqrt{\nu t}$ and does not saturate at the Blasius thickness: it is an exact member of a class of unsteady layers, not a description of the complete Rayleigh-to-Blasius transition. The unsteady form of the case $m=-1$ is an elliptic boundary-value problem in $(\eta ,\tau )$ rather than an evolution problem, and only its steady limit is solved here. And the results of §6 assume constant properties throughout.

The mechanism identified in §5 is not special to this flow. Whenever a solution is invariant under a one-parameter group, the convective operator restricted to invariant scalars loses its dependence on the velocity field, and any passive scalar reduces to an ordinary differential equation with the flow entering only through the source term. That observation should apply equally to species transport, to the three-dimensional problem, and to any other reduction in which a distinguished invariant of the stretching operator can be identified.