Unsteady MHD Boundary Layer Similarity Solutions

1. Introduction

The diffusion-time similarity transformation proposed by Sun [1] reduces the two-dimensional unsteady laminar boundary layer equations from three independent variables to two. Writing the stream function as $\psi =U\left(x\right)\delta \left(x\right)f(\eta ,\tau )$ with $\eta =y/\delta \left(x\right)$ and $\tau =\nu t/{\delta }^2\left(x\right)$, the streamwise coordinate is eliminated entirely provided the coefficients of the reduced equation are constant, and a single two-dimensional problem then describes the flow at every station. This is a partial reduction rather than a reduction to an ordinary differential equation: the symmetry algebra of the unsteady boundary layer equations, classified by Ma and Hui [2], is too small to permit the latter for a general outer stream. The construction has since been applied to conducting fluids by Fu, Ni and Zhang [3], who obtained closed-form approximations for stagnation and converging flows under a spanwise field.

The reduction raises a question that does not arise for the classical steady similarity variables, and that appears not to have been addressed. Because $\tau$ depends on x through $\delta \left(x\right)$, surfaces of constant $\tau$ are not in general transverse to the direction along which the parabolic boundary layer equations propagate information. One should therefore ask whether the reduced equation may be integrated as an initial-value problem in $\tau$ at all. We show below that it may not: the equation changes type at a critical value of $\tau$ that depends on the outer-flow exponent, and beyond that value forward marching is ill posed in the sense of Hadamard.

The threshold turns out to have a transparent physical meaning. For a power-law outer stream $U=C{x}^m$ the parabolicity parameter reduces to $\mathrm{\Lambda }=(1-m)ut/x$, so that for a flat plate the reduced problem is well posed precisely in the region $t<x/u$ that has not yet been reached by the starting front propagating from the leading edge. This is the same region in which Stewartson’s small-time solution [4] holds, recovered here from the type of the transformed equation rather than from a matched asymptotic expansion. For accelerated outer flows, $m\ge 1$, the obstruction is absent and the march is unconditionally well posed; the stagnation-point case $m=1$ is distinguished twice over, since it is also the only case in which a spatially uniform applied magnetic field is compatible with similarity.

We work throughout with an electrically conducting fluid at low magnetic Reynolds number, ${\mathit{R}}_m={\mu }_0\sigma UL\ll 1$, the regime relevant to liquid metals in ducts, moulds and fusion blankets [5,6]. The induced field is then negligible and the electromagnetic problem reduces to the quasi-static determination of the electric potential. The magnetic case is of interest here not only for its own sake but because the Lorentz force provides an independent control on the velocity profile: it alters $f_{\eta }$, hence $\mathrm{\Lambda }$, hence the location of the threshold, and it supplies a second verification route through the strong-field asymptotics developed in Sec. Section 4.3.

The paper is organised as follows. Section 2 sets out the governing equations, the quasi-static closure of the Lorentz force, and the similarity reduction, together with the conditions under which the coefficients are constant. Section 3 contains the central result: the parabolicity condition, its physical interpretation, and the consequences for the marching direction. Section 4 collects the limiting solutions used for verification, including the two-term strong-field expansion. Section 5 describes the numerical scheme and its verification, and Sec. Section 6 reports the numerical test of the threshold and its physical consequences.

2. Formulation and Similarity Reduction

2.1. Governing Equations and the Lorentz Closure

Consider the two-dimensional unsteady flow of an incompressible, viscous, electrically conducting fluid past a plane wall, with a magnetic field $\mathit{B}=B_0\left(x\right){\widehat{\mathit{e}}}_y$ applied normal to the wall (Figure 1). At low ${\mathit{R}}_m$ the current follows from Ohm’s law for a moving conductor, $\mathit{J}=\sigma (-\nabla \varphi +\mathit{u}\times \mathit{B})$, with $\nabla ·\mathit{J}=0$.

The force involves two successive cross products, and it is worth writing both out. The motional term is

$$\mathit{u}\times \mathit{B}=u{B}_0\left({\widehat{\mathit{e}}}_x\times {\widehat{\mathit{e}}}_y\right)+v{B}_0\left({\widehat{\mathit{e}}}_y\times {\widehat{\mathit{e}}}_y\right)=u{B}_0{\widehat{\mathit{e}}}_z,$$

(1)

so the wall-normal velocity does not contribute and the induced current is purely spanwise, $J_z=\sigma \left(E_z+u{B}_0\right)$. This spanwise current then crosses the wall-normal field a second time, and since ${\widehat{\mathit{e}}}_z\times {\widehat{\mathit{e}}}_y=-{\widehat{\mathit{e}}}_x$, the force is directed along the wall:

$$\mathit{J}\times \mathit{B}=J_z{B}_0\left({\widehat{\mathit{e}}}_z\times {\widehat{\mathit{e}}}_y\right)=-\sigma B_0\left(E_z+u{B}_0\right){\widehat{\mathit{e}}}_x.$$

(2)

A spanwise current in a wall-normal field therefore produces a streamwise force; this is the mechanism responsible for the Hartmann layer [16]. Note that the result depends on $B_0^2$, so reversing the field direction, $\mathit{B}=-B_0{\widehat{\mathit{e}}}_y$, leaves (2) unchanged once $E_z$ is reversed with it.

The value of $E_z$ in (2) is fixed by the electrical boundary conditions. Applying the outer-flow momentum balance at the edge of the layer,

$${\displaystyle \frac{\partial U}{\partial t}}+U{\displaystyle \frac{dU}{dx}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}+{\displaystyle \frac{1}{\rho }}{(\mathit{J}\times \mathit{B})}_x{|}_{u=U},$$

(3)

shows that the short-circuited choice $E_z=0$, which gives ${(\mathit{J}\times \mathit{B})}_x=-\sigma B_0^{2}u$, is incompatible with a steady outer stream: with $-{\rho }^{-1}\partial p/\partial x=UdU/dx$ it requires $\sigma B_0^{2}U=0$. It is admissible only if the free stream is allowed to decay as $U\propto exp(-\sigma B_0^{2}t/\rho )$, which destroys the separability the reduction requires. We therefore take $E_z=-U{B}_0$, so that the current vanishes in the outer stream and

$${(\mathit{J}\times \mathit{B})}_x=-\sigma B_0^2\left(u-U\right),$$

(4)

in agreement with the closure used by Chiam [7] and Parand et al. [8]. Under the boundary layer approximation the governing equations are then

$${\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}=0,$$

(5)

$${\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+v{\displaystyle \frac{\partial u}{\partial y}}=U{\displaystyle \frac{dU}{dx}}+\nu {\displaystyle \frac{{\partial }^{2}u}{\partial y^2}}-{\displaystyle \frac{\sigma B_0^2}{\rho }}\left(u-U\right),$$

(6)

with $u=v=0$ at $y=0$, $u\to U\left(x\right)$ as $y\to \infty$, and an impulsive start $u(x,y,0)=U\left(x\right)$ for $y>0$.

Remark 1

(Why a wall-normal field). The alternative orientation is a spanwise field $\mathit{B}=B_0{\widehat{\mathit{e}}}_z$, for which $\mathit{u}\times \mathit{B}=-u{B}_0{\widehat{\mathit{e}}}_y+v{B}_0{\widehat{\mathit{e}}}_x$ lies in the plane of the flow. The streamwise force is then ${(\mathit{J}\times \mathit{B})}_x=J_y{B}_0=\sigma B_0(E_y-u{B}_0)$, which for $E_y=U{B}_0$ reduces to $-\sigma B_0^2(u-U)$: thesameexpression as (4). The two orientations therefore give identical momentum equations, and differ only in the electrostatic problem that must be solved to obtain them.

That difference is the reason for the present choice. With $\mathit{B}=B_0{\widehat{\mathit{e}}}_y$ the motional term (1) is purely spanwise, so in two dimensions $\nabla ·\mathit{J}=\partial J_z/\partial z=0$ holds identically: no potential is induced, and $E_z$ is fixed by the external circuit alone. With $\mathit{B}=B_0{\widehat{\mathit{e}}}_z$, by contrast, $\nabla ·(\mathit{u}\times \mathit{B})\ne 0$ in general, ϕ must be obtained from a Poisson problem, and for an insulating wall the potential contributes a further half of the total force, doubling the effective coefficient [3]. The wall-normal orientation is thus the one for which the two-dimensional low-${\mathit{R}}_m$ problem closes without an auxiliary electrostatic equation; all results below carry over to the spanwise case under $M\to 2M$.

2.2. The Diffusion-Time Transformation

Introducing the stream function $u={\psi }_y$, $v=-{\psi }_x$, Eqs. (5) and (6) collapse to

$${\psi }_{ty}+{\psi }_y{\psi }_{xy}-{\psi }_x{\psi }_{yy}=U{\displaystyle \frac{dU}{dx}}+\nu {\psi }_{yyy}-{\displaystyle \frac{\sigma B_0^2}{\rho }}\left({\psi }_y-U\right),$$

(7)

and we apply the transformation of Ref. [1],

$$\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 }^2\left(x\right)}}.$$

(8)

The derivatives of the new variables are ${\eta }_x=-\eta {\delta }_x/\delta$, ${\eta }_y=1/\delta$, ${\tau }_x=-2\tau {\delta }_x/\delta$ and ${\tau }_t=\nu /{\delta }^2$, whence

$${\psi }_y=U{f}_{\eta },{\psi }_{yy}={\displaystyle \frac{U}{\delta }}{f}_{\eta \eta },{\psi }_{yyy}={\displaystyle \frac{U}{{\delta }^2}}{f}_{\eta \eta \eta },$$

(9)

$${\psi }_{ty}={\displaystyle \frac{U\nu }{{\delta }^2}}{f}_{\eta \tau },{\psi }_x={\left(U\delta \right)}_{x}f-U\eta {\delta }_x{f}_{\eta }-2U\tau {\delta }_x{f}_{\tau },$$

(10)

$${\psi }_{xy}=U_x{f}_{\eta }-{\displaystyle \frac{\eta U{\delta }_x}{\delta }}{f}_{\eta \eta }-{\displaystyle \frac{2U\tau {\delta }_x}{\delta }}{f}_{\eta \tau }.$$

(11)

The terms proportional to $\eta$ in ${\psi }_y{\psi }_{xy}-{\psi }_x{\psi }_{yy}$ cancel identically, and dividing (7) by $U\nu /{\delta }^2$ gives

$$\begin{array}{cc}\hfill f_{\eta \eta \eta }+af{f}_{\eta \eta }& +b\left(1-f_{\eta }^2\right)+M\left(1-f_{\eta }\right)\hfill \\ \hfill & =f_{\eta \tau }+c\tau \left(f_{\tau }{f}_{\eta \eta }-f_{\eta }{f}_{\eta \tau }\right),\hfill \end{array}$$

(12)

with

$$\begin{array}{cccc}\hfill a& ={\displaystyle \frac{\delta {\left(U\delta \right)}_x}{\nu }},\hfill & \hfill b& ={\displaystyle \frac{{\delta }^2{U}_x}{\nu }},\hfill \\ \hfill c& ={\displaystyle \frac{2U\delta {\delta }_x}{\nu }},\hfill & \hfill M& ={\displaystyle \frac{\sigma B_0^2{\delta }^2}{\rho \nu }}.\hfill \end{array}$$

(13)

Setting $M=0$ recovers the equation of Ref. [1], and the form (12) agrees with that of Ref. [3] under the correspondence $(\alpha ,\beta ,\gamma )\leftrightarrow (a,b,c)$. The boundary conditions become

$$f(0,\tau )=0,f_{\eta }(0,\tau )=0,f_{\eta }\to 1(\eta \to \infty ).$$

(14)

Remark 2.

Had the short-circuited closure been used, the magnetic term in (12) would read $-M{f}_{\eta }$. Evaluating the equation as $\eta \to \infty$, where $f_{\eta }\to 1$ and all derivatives vanish, then leaves the residual $-M$, so no solution satisfying (14) exists for $M\ne 0$. With (4) the same limit gives $b(1-1)+M(1-1)=0$ identically. This compatibility check is worth performing on any similarity formulation of a damped boundary layer.

2.3. Complete Determination of the Admissible Flows

Equation (12) is a genuine two-variable reduction if and only if a, b, c and M are independent of x. Rather than postulate a form for U and $\delta$ and verify constancy afterwards, we take constancy as the defining condition and integrate the resulting system, which determines the admissible flows completely.

Writing $S\equiv {\delta }^2$, the definitions (13) of b and c are two first-order ordinary differential equations for the pair $\left(U,S\right)$,

$$S{U}_x=b\nu ,U{S}_x=c\nu ,$$

(15)

while a is not independent, since

$$a={\displaystyle \frac{{\delta }^2{U}_x}{\nu }}+{\displaystyle \frac{U\delta {\delta }_x}{\nu }}=b+{\displaystyle \frac{c}{2}}.$$

(16)

The condition on M separates from (15): from (13), $M=\sigma B_0^2{\delta }^2/\left(\rho \nu \right)$ is constant if and only if

$$\begin{array}{c}{B}_0\left(x\right)\delta \left(x\right)=\mathrm{const},\\ \mathrm{that}\mathrm{is},M={\mathit{Ha}}_{\delta }^2,{\mathit{Ha}}_{\delta }\equiv B_0\delta \sqrt{{\displaystyle \frac{\sigma }{\rho \nu }}}.\end{array}$$

(17)

The requirement is therefore simply that the Hartmann number formed with the local layer thickness be the same at every station; the applied field must scale inversely with the layer it damps. This statement is independent of the solution of (15) and holds for every case below.

Lemma 1

(Admissible outer flows). The system (15) with b, c constant admits exactly three families, distinguished by the value of $K\equiv b+c=2a-b$.

(A) $K\ne 0$: power-law flows.Dividing the two equations of (15) gives $dlnS/dlnU=c/b$, whence $S=A{U}^{c/b}$; substituting back and integrating yields $U^{K/b}\propto x-x_0$, that is

$$U\left(x\right)=C{\left(x-x_0\right)}^m,{\delta }^2\left(x\right)={\displaystyle \frac{K\nu \left(x-x_0\right)}{U\left(x\right)}},$$

(18)

with $m=b/K$ and

$$b=mK,a={\displaystyle \frac{(1+m)K}{2}},c=(1-m)K.$$

(19)

Constancy of M then requires, by (17), $B_0\propto {(x-x_0)}^{(m-1)/2}$.

(B) $K=0$, $b\ne 0$: exponential flows.Here $c=-b$, $a=b/2$, and $S=A/U$, so that (15) gives $A{U}_x/U=b\nu$ and

$$U\left(x\right)=U_0{e}^{kx},{\delta }^2\left(x\right)={\displaystyle \frac{A}{U_0}}{e}^{-kx},b={\displaystyle \frac{Ak}{\nu }},$$

(20)

with $B_0\propto e^{kx/2}$ for constant M.

(C) $b=c=0$: uniform flows.Then U and δ are both constant, $a=0$, and $B_0$ is uniform. This is the degenerate case of Sec. Section 4.2.

Three points deserve emphasis. First, family (A) contains a virtual origin $x_0$ that the usual statement of the similarity conditions suppresses; it is a genuine free parameter, and for the flat plate ($b=0$, U const) it appears as ${\delta }^2=c\nu (x-x_0)/U$, the familiar leading-edge offset.

Second, family (B) is not a power law and is not reached by the ansatz ${\delta }^2=K\nu x/U$, which degenerates when $K=0$. It is nevertheless an exact similarity family: an exponentially accelerating outer stream over an exponentially thinning layer, with an exponentially strengthening field. Its coefficients satisfy $c=-b$ identically, so by the criterion of Sec. Section 3 an accelerating member ($k>0$, hence $b>0$, $c<0$) is unconditionally well posed, while a decelerating one ($k<0$) is conditional. For this family the parabolicity parameter takes the particularly simple form $\mathrm{\Lambda }=-kut$.

Third, the residual freedom in K within family (A) is a gauge choice for $\delta$, not additional physics. We fix it by the Falkner–Skan convention $K=2/(1+m)$, so that

$$a=1,b={\displaystyle \frac{2m}{1+m}},c={\displaystyle \frac{2(1-m)}{1+m}}=2(1-b),$$

(21)

$\delta =\sqrt{2\nu x/\left[\right(1+m\left)U\right]}$ with $x_0=0$, and b is the Hartree pressure-gradient parameter. Table 1 summarises the classification; the sign of c, which will be shown in Sec. Section 3 to govern the type of the reduced equation, is determined entirely by whether the outer flow accelerates faster or slower than linearly.

3. Change of Type and Well-Posedness

3.1. The Parabolicity Condition

The terms of (12) containing $\tau$-derivatives may be grouped as

$$\left(1-c\tau f_{\eta }\right)f_{\eta \tau }+c\tau f_{\eta \eta }{f}_{\tau },$$

(22)

so that the coefficient multiplying the highest mixed derivative is not unity but $1-\mathrm{\Lambda }$, where

$$\mathrm{\Lambda }(\eta ,\tau )\equiv c\tau f_{\eta }(\eta ,\tau ).$$

(23)

This single observation determines the character of the reduced problem.

Lemma 2

(Parabolicity condition). Let f solve (12) and consider a short-wave perturbation $f\mapsto f+ϵe^{\mathrm{i}k\eta +\lambda \tau }$ with $k\to \infty$ at frozen coefficients. Retaining the highest derivative on each side, the principal balance is $-\mathrm{i}{k}^3=(1-\mathrm{\Lambda })\mathrm{i}k\lambda$, whence

$$\lambda =-{\displaystyle \frac{k^2}{1-\mathrm{\Lambda }}}.$$

(24)

Forward evolution in τ is therefore parabolic and well posed where $\mathrm{\Lambda }<1$, and backward-parabolic—hence ill posed as an initial-value problem in increasing τ—where $\mathrm{\Lambda }>1$.

Because $f_{\eta }\le 1$, with equality attained only in the free stream, the critical condition is met first, and simultaneously across the whole outer part of the layer, at

$${\tau }_{\mathrm{crit}}={\displaystyle \frac{1}{c}}(c>0).$$

(25)

For $c\le 0$, that is for $m\ge 1$, one has $1-\mathrm{\Lambda }\ge 1$ everywhere and the forward march is unconditionally well posed. Equation (24) also shows that the effective diffusivity in $\tau$ is ${(1-\mathrm{\Lambda })}^{-1}$, so that the approach to the threshold is accompanied by a progressive stiffening of the problem even before the type changes.

3.2. Physical Interpretation: The Starting Front

The parameter $\mathrm{\Lambda }$ is not an artefact of the algebra. Using (18) and (19), $\tau =\nu t/{\delta }^2=Ut/\left(Kx\right)$ and $c=(1-m)K$, so that

$$\mathrm{\Lambda }=(1-m){\displaystyle \frac{ut}{x}}.$$

(26)

For a flat plate, $m=0$, the condition $\mathrm{\Lambda }<1$ is exactly $t<x/u$: the reduction may be marched forwards in $\tau$ precisely in the region that has not yet been reached by the starting front propagating downstream from the leading edge. This is the classical criterion delimiting the validity of Stewartson’s small-time solution for the impulsively started plate [4], obtained here from the type of the transformed equation rather than from a matched asymptotic expansion. The last column of Table 1 lists ${\tau }_{\mathrm{crit}}$ for each case; decelerated flows reach the threshold sooner than the flat plate, and accelerated flows with $m>1$ never reach it. For the exponential family (B) of Lemma 1 the same parameter is $\mathrm{\Lambda }=-kut$, so that an accelerating exponential stream is well posed for all $\tau$.

The magnetic field enters this criterion indirectly but not negligibly. Since the Lorentz force makes the profile fuller, raising $f_{\eta }$ at fixed $\eta$, it raises $\mathrm{\Lambda }$ throughout the layer and therefore brings the interior of the layer to the threshold earlier, even though ${\tau }_{\mathrm{crit}}$ itself, fixed by the free stream, is unchanged.

3.3. The Correct Marching Direction

The loss of well-posedness for $\mathrm{\Lambda }>1$ is not a defect of the reduction but a statement about the direction in which information travels. At fixed t, decreasing $\tau$ corresponds to increasing x. The region $\tau >1/c$ is therefore the already-established region near the leading edge, and there the boundary layer equations are parabolic in x with the marching direction downstream, that is, towards decreasing $\tau$. The properly posed formulation for $m<1$ and $\tau >1/c$ is thus a terminal-value problem, integrated backwards in $\tau$ from the steady Falkner–Skan state at $\tau \to \infty$, and (24) with $\mathrm{\Lambda }>1$ gives decay rather than growth in that direction. Attempting instead to march forwards through $\tau =1/c$ produces the grid-dependent divergence documented in Sec. Section 6.

Two consequences follow. First, any statement about the transient of a decelerated or zero-pressure-gradient layer obtained by forward marching in $\tau$ beyond $\tau =1/c$ is not supported by the equation. Second, the stagnation-point case $m=1$ is the natural setting for unsteady computations within this framework, being both unconditionally well posed and compatible with a uniform applied field.

3.4. Relation to the Group-Theoretic Classification

It is worth locating this obstruction within the Lie-group classification of similarity reductions of the unsteady boundary layer equations given by Ma and Hui [2]. Their analysis shows that the symmetry algebra of the unsteady equations admits only a small number of families of invariant solutions, distinguished by the dependence of the free stream on x and t; it is too small to reduce the general problem to an ordinary differential equation. The transformation (8) is a partial reduction that sidesteps this by eliminating x at the cost of constraining U and $B_0$. The obstruction identified above is the further, and to our knowledge previously unremarked, price of that construction: because $\tau$ depends on x, the level surfaces of $\tau$ cross the characteristic surfaces of the original parabolic system, and the reduced problem inherits a change of type at the crossing. A characterisation of this obstruction directly in terms of the characteristic surfaces of the unreduced system would be a natural extension of Ref. [2].

4. Limiting Solutions

Three limits of (12) admit closed-form or classical results. The first two are known and are recorded here because they supply the initial condition and the verification cases for the numerical scheme; the third is derived below and appears to be new.

4.1. Steady Limit

Setting $f_{\tau }=f_{\eta \tau }=0$ removes the c-group entirely and leaves the MHD Falkner–Skan equation

$$f^{\prime \prime \prime }+af{f}^{\prime \prime }+b\left(1-f^{\prime 2}\right)+M\left(1-f^{\prime }\right)=0,$$

(27)

with $f\left(0\right)=f^{\prime }\left(0\right)=0$, $f^{\prime }(\infty )=1$, studied by Chiam [7], Abbasbandy and Hayat [9], and Parand et al. [8], whose parameter $M^2$ is the present M. For $M=0$ it is the classical Falkner–Skan equation [10,11,12].

4.2. Small-Time Limit

As $\tau \to 0$ the physical layer thickness $2\sqrt{\nu t}$ is much smaller than $\delta$, so $\eta =O\left(\sqrt{\tau }\right)$, $f_{\eta \eta \eta }$ and $f_{\eta \tau }$ are $O\left({\tau }^{-1}\right)$, while $af{f}_{\eta \eta }$, $b(1-f_{\eta }^2)$ and $M(1-f_{\eta })$ are all $O\left(1\right)$. The leading balance is

$$f_{\eta \tau }=f_{\eta \eta \eta }+M\left(1-f_{\eta }\right)+O\left(\tau \right),$$

(28)

which is exact for all $\tau$ in the degenerate case $a=b=c=0$. With $u=f_{\eta }$ and $w=1-u$ this is $w_{\tau }=w_{\eta \eta }-Mw$ subject to $w(0,\tau )=1$, $w(\infty ,\tau )=0$, $w(\eta ,0)=0$, whose solution is the MHD Rayleigh–Stokes profile of Chang and Yen [15],

$$\begin{array}{cc}\hfill w=\frac{1}{2}[& e^{-\eta \sqrt{M}}\mathrm{erfc}\left(\frac{\eta }{2\sqrt{\tau }}-\sqrt{M\tau }\right)\hfill \\ \hfill & +e^{+\eta \sqrt{M}}\mathrm{erfc}\left(\frac{\eta }{2\sqrt{\tau }}+\sqrt{M\tau }\right)],\hfill \end{array}$$

(29)

with wall shear

$$f_{\eta \eta }(0,\tau )=\sqrt{M}\mathrm{erf}\left(\sqrt{M\tau }\right)+{\displaystyle \frac{e^{-M\tau }}{\sqrt{\pi \tau }}}.$$

(30)

For $M=0$ these reduce to $u=\mathrm{erf}(\eta /2\sqrt{\tau })$ and ${\left(\pi \tau \right)}^{-1/2}$; as $\tau \to \infty$ they tend to the exponential Hartmann profile of thickness $M^{-1/2}$ and to $\sqrt{M}$. Equation (29) supplies the starting profile at $\tau ={\tau }_0\ll 1$ with error $O\left({\tau }_0\right)$, and together with (30) provides an exact benchmark for the unsteady solver.

4.3. Strong-Field Expansion

For $M\gg 1$ the velocity defect is confined to a Hartmann layer of thickness $M^{-1/2}$. Introducing $\zeta =\sqrt{M}\eta$ and $f=\epsilon F\left(\zeta \right)$ with $\epsilon =M^{-1/2}$, Eq. (27) becomes

$$F^{\prime \prime \prime }+{\epsilon }^2\left[aF{F}^{\prime \prime }+b\left(1-F^{\prime 2}\right)\right]+\left(1-F^{\prime }\right)=0.$$

(31)

At leading order $F_0^{\prime \prime \prime }+(1-F_0^{\prime })=0$ with $F_0\left(0\right)=F_0^{\prime }\left(0\right)=0$, $F_0^{\prime }(\infty )=1$, giving

$$F_0^{\prime }=1-e^{-\zeta },F_0=\zeta -1+e^{-\zeta }.$$

(32)

Writing $F=F_0+{\epsilon }^2{F}_1$ and $W_1=F_1^{\prime }$,

$$W_1^{\prime \prime }-W_1=-a\zeta e^{-\zeta }+(a-2b)e^{-\zeta }-(a-b)e^{-2\zeta },$$

(33)

whose bounded solution with $W_1\left(0\right)=0$ is

$$W_1=\left({\displaystyle \frac{a}{4}}{\zeta }^2+\left(b-{\displaystyle \frac{a}{4}}\right)\zeta \right)e^{-\zeta }+{\displaystyle \frac{a-b}{3}}\left(e^{-\zeta }-e^{-2\zeta }\right).$$

(34)

Since $F_1^{\prime \prime }\left(0\right)=W_1^{\prime }\left(0\right)=\frac{2}{3}b+\frac{1}{12}a$, the wall shear is

$$f_{\eta \eta }\left(0\right)\simeq M^{1/2}+\left({\displaystyle \frac{a}{12}}+{\displaystyle \frac{2b}{3}}\right)M^{-1/2}+O\left(M^{-3/2}\right),$$

(35)

and the integral thicknesses satisfy ${\delta }^{*}/\delta \simeq M^{-1/2}$, $\theta /\delta \simeq \frac{1}{2}{M}^{-1/2}$, so that the shape factor tends to the value of a purely exponential profile,

$$H={\displaystyle \frac{{\delta }^{*}}{\theta }}\to 2(M\to \infty ).$$

(36)

Equation (35) confirms the $\sqrt{M}$ friction scaling reported in Ref. [3] and fixes its prefactor: the leading coefficient is exactly unity in the present normalisation, whereas the closed-form estimate of Ref. [3] gives $1/\sqrt{3}$. The discrepancy originates in the near-wall Taylor truncation of $f{f}^{\prime \prime }$ and $f^{\prime 2}$ used there, which preserves the exponent but not the coefficient. The accuracy of (35) is assessed against independent spectral data in Sec. Section 6.

5. Numerical Method and Verification

The steady problem (27) is solved by fourth-order collocation with residual-based mesh adaptation, tolerance ${10}^{-10}$, and truncation length ${\eta }_{\infty }=max(6,min(26,14/\sqrt{M}))$; selected cases were reproduced by shooting with an eighth-order Runge–Kutta integrator and a secant iteration on $f^{\prime }\left({\eta }_{\infty }\right)-1$, the two agreeing to twelve significant figures for $b\ge 0$. For $M=0$ the scheme returns $f^{\prime \prime }\left(0\right)=0.4696000$ for the flat plate and $1.2325877$ for plane stagnation flow, and the separation eigenvalue $b_{\mathrm{sep}}=-0.198838$, reproducing the classical values [11,13] to better than $6\times {10}^{-7}$.

For the unsteady problem we take $u=f_{\eta }$ as the primary unknown on a uniform grid ${\eta }_j=jh$, $h={\eta }_{\infty }/N$, recovering f by cumulative trapezoidal quadrature. Since the solution spans several decades of $\tau$, we march in $\xi =ln\tau$, so that $\tau {\partial }_{\tau }={\partial }_{\xi }$ and (12) reads

$$\begin{array}{cc}\hfill \left(e^{-\xi }-cu\right)u_{\xi }+c{f}_{\xi }{u}_{\eta }=& u_{\eta \eta }+af{u}_{\eta }\hfill \\ \hfill & +b\left(1-u^2\right)+M\left(1-u\right).\hfill \end{array}$$

(37)

Spatial derivatives use second-order central differences and ${\partial }_{\xi }$ the BDF2 formula, with implicit Euler on the first step. The quadratic term is treated by Newton linearisation and the non-local couplings through f and $f_{\xi }$ by Picard iteration, so each iterate requires only a tridiagonal solve; iteration stops when successive iterates differ by less than ${10}^{-11}$, typically after three to six sweeps. The scheme is fully implicit and imposes no stability restriction on $\Delta \xi$ within the well-posed region $\mathrm{\Lambda }<1$—a point of some importance for the test in Sec. Section 6, since it means that any divergence observed there cannot be attributed to a step-size restriction. The initial condition at ${\tau }_0={10}^{-3}$ is (29).

Three verification tests were performed. First, with $a=b=c=0$ the computed solution is compared with (29); Figure 2(a) shows the relative error in the wall shear for $M=0$–25, below ${10}^{-5}$ for $\tau \gtrsim 0.1$ and largest at the earliest times, where the layer thickness $2\sqrt{{\tau }_0}\approx 0.06$ is resolved by only a few cells. Second, Table 2 and Figure 2(b) report a refinement study with spatial and temporal resolution varied independently; the observed spatial order is $2.00$, $2.00$, $1.99$, $1.95$ over four refinements. Third, for $m=1$ the value of $f_{\eta \eta }(0,\tau )$ obtained by marching to $\tau =60$ agrees with the independently computed solution of (27) to $8\times {10}^{-7}$ at $M=0$, $7\times {10}^{-6}$ at $M=4$ and $1\times {10}^{-4}$ at $M=25$, the degradation reflecting the thinning of the Hartmann layer at fixed h.

6. Results

6.1. Test of the Strong-Field Expansion

Table 3 compares the computed wall shear with the independent spectral results of Ref. [8] and with the expansion (35), for a strongly decelerated ($b=-3$) and an accelerated ($b=4/3$) outer flow. Agreement between the present computation and the published values is to seven or eight significant figures across four decades in M, which validates the solver; the interest here is the last two columns. The two-term expansion is accurate to $3\times {10}^{-3}$ already at $M=25$ and to better than ${10}^{-8}$ at $M\ge 2500$, and the convergence is consistent with the predicted $O\left(M^{-3/2}\right)$ remainder: between $M=2500$ and $M={10}^4$ at $b=4/3$ the relative error falls by a factor of 16, against the factor of 8 expected from the remainder alone and 32 from a naive $M^{-2}$ estimate.

Figure 3 shows the accompanying prediction for the integral quantities. The wall shear follows $M^{1/2}$ and the displacement thickness $M^{-1/2}$, and the shape factor collapses onto $H=2$ from above, reaching $H<2.006$ at $M=100$ for every b examined. The last result gives a compact diagnostic: since H is bounded below by 2 and takes values near $2.6$ for Blasius flow, the departure of H from 2 measures how far the magnetic field is from having taken over the momentum balance entirely.

6.2. Numerical Test of the Parabolicity Threshold

Table 4 and Figure 4 test Lemma 2 directly. For the flat plate, $c=2$ and (25) predicts ${\tau }_{\mathrm{crit}}=0.5$. On four successive meshes spanning a factor of eight in h, the computed wall shear at $\tau =0.20$, $0.35$ and $0.45$ agrees to six significant figures, with differences decreasing as $h^2$; the solution below the threshold is thus fully grid-converged.

Beyond $\tau =0.5$ the behaviour changes qualitatively. The computed values diverge, and—this is the essential point—they diverge differently on each mesh, non-monotonically in N: at $\tau =1$ the four grids return $1.15\times {10}^{58}$, $1.37\times {10}^{32}$, $5.00\times {10}^1$ and $1.22\times {10}^{83}$. Refinement does not restore convergence, and the coarsest grid is not the worst. This is precisely what (24) predicts for $\mathrm{\Lambda }>1$: the amplification rate is proportional to $k^2$ and hence unbounded as $h\to 0$, so the outcome is governed by whatever short-wave content the particular mesh admits. Because the scheme is fully implicit, the divergence cannot be attributed to a time-step restriction.

Figure 4(b) shows the complementary case. For $m=1$, 2 and 3, all of which have $c\le 0$, the march proceeds smoothly over five decades in $\tau$ and converges to the respective steady values. The short-time behaviour is common to all cases, as it must be, since $c\tau \to 0$ as $\tau \to 0$ irrespective of the sign of c; the histories separate only once $\left|c\right|\tau =O\left(1\right)$.

6.3. The Unconditionally Well-Posed Case

For $m=1$ the reduction is well posed for all $\tau$, and the complete transient may be computed. Figure 5 shows the evolution from the Rayleigh–Stokes state to the steady MHD stagnation-point profile: at $\tau =0.005$ the computed profile is indistinguishable from (29), and by $\tau \approx 1$ the steady state has essentially been attained.

Figure 6(a) shows that at small $\tau$ all wall-shear histories collapse onto the field-free Rayleigh singularity ${\left(\pi \tau \right)}^{-1/2}$, since the magnetic term is $O\left(M\right)$ while the diffusive terms are $O\left({\tau }^{-1}\right)$; the field is felt only once $M\tau =O\left(1\right)$. Thereafter each curve relaxes to its own plateau. Defining the relaxation time ${\tau }_{1\%}$ as the value of $\tau$ beyond which $f_{\eta \eta }(0,\tau )$ remains within $1\%$ of its steady value, we obtain

The product $M{\tau }_{1\%}$ tends to a constant near $1.8$, so that ${\tau }_{1\%}\sim 1.8/M$ at strong field: the magnetic damping time $\rho /\left(\sigma B_0^2\right)$ replaces the viscous diffusion time ${\delta }^2/\nu$ as the controlling timescale, and the field accelerates the approach to the steady state as well as altering the state itself.

Figure 6(c) shows the shape factor. Every case starts from the value $H={(\sqrt{2}-1)}^{-1}=2.4142$ of the error-function profile, independently of M, and relaxes monotonically to its steady value—$2.216$ at $M=0$, $2.079$ at $M=4$, $2.018$ at $M=25$—approaching the limit (36). The monotonicity indicates that the profile becomes fuller throughout the transient, with no intermediate overshoot.

7. Conclusions

We have extended the diffusion-time similarity reduction of Sun [1] to the unsteady boundary layer of a conducting fluid at low magnetic Reynolds number and examined the type of the reduced equation. The main conclusions are the following.

The reduced problem changes type. Grouping the time-derivative terms shows that the coefficient of the highest mixed derivative is $1-\mathrm{\Lambda }$ with $\mathrm{\Lambda }=c\tau f_{\eta }$, and a frozen-coefficient analysis gives the growth rate $\lambda =-k^2/(1-\mathrm{\Lambda })$. Forward marching in the diffusion time is therefore well posed only where $\mathrm{\Lambda }<1$, and for $c>0$ the threshold is reached at ${\tau }_{\mathrm{crit}}=1/c$.

The threshold is the starting front. For a power-law outer stream $\mathrm{\Lambda }=(1-m)ut/x$, so the condition $\mathrm{\Lambda }<1$ is, for a flat plate, exactly $t<x/u$. The restriction on the reduction thus coincides with the classical criterion delimiting Stewartson’s small-time solution, recovered here from the type of the equation rather than from a matched expansion. Since decreasing $\tau$ corresponds to increasing x, the correct formulation beyond the threshold is a terminal-value problem integrated backwards in $\tau$, which is downstream marching.

Numerical experiment confirms the threshold sharply. For $c=2$ the wall shear is grid-converged to six figures up to $\tau =0.45$ on four meshes and diverges non-monotonically in the mesh beyond $\tau =0.5$, over eighty orders of magnitude, despite a fully implicit scheme. For $m\ge 1$ ($c\le 0$) the march is unconditionally stable, and the stagnation-point case $m=1$ is doubly distinguished, being also the only case compatible with a uniform applied field.

A two-term strong-field expansion, $f_{\eta \eta }\left(0\right)\simeq M^{1/2}+(a/12+2b/3)M^{-1/2}$, reproduces independent spectral benchmarks to eight significant figures and fixes the prefactor of the $\sqrt{M}$ friction law at unity, in place of the $1/\sqrt{3}$ obtained in Ref. [3] from a near-wall truncation. In the same limit ${\delta }^{*}/\delta \to M^{-1/2}$ and $H\to 2$, and the relaxation time of the transient scales as $M^{-1}$.

Finally, a formulation point of general applicability: the quasi-static Lorentz force must be closed as $-\sigma B_0^2(u-U)$ if the outer stream is steady, so that it enters the reduced equation as $+M(1-f_{\eta })$. The alternative closure leaves a residual $-M$ in the free stream and admits no solution satisfying $f_{\eta }\to 1$ (Remark 2).

Two extensions suggest themselves. The first is to solve the terminal-value problem for $m<1$ and so close the flat-plate transient within the similarity framework, which would permit a direct quantitative comparison with Stewartson’s matched asymptotic solution across $\mathrm{\Lambda }=1$. The second is to characterise the change of type directly in terms of the characteristic surfaces of the unreduced system, within the group-theoretic setting of Ref. [2].