Past exam of the mathematics course of the University of Cambridge / 2018 / iii / Paper 331 (Hydrodynamic Stability)

For irrotational flow, ${\mathbf{u}}_j=\nabla {\varphi }_j$. An impermeable stationary horizontal wall has zero normal velocity, so the kinematic boundary condition isNo tangential no-slip condition is imposed on these inviscid fluids.

The interface is material: fluid particles on $F=z-\eta (x,t)=0$ remain there. Its material derivative therefore vanishes, giving $w_j={\eta }_t+u_j{\eta }_x=D_j\eta /Dt$ on the interface. Since $w_j={\partial }_z{\varphi }_j$, this is the kinematic boundary condition for a free-surface graph on each side. The material derivatives use the respective tangential fluid velocities; they enforce the same moving normal boundary.

The Unsteady Bernoulli equation gives $p_j+{\rho }_j[{\partial }_t{\varphi }_j+|\nabla {\varphi }_j{|}^2/2+gz]=C_j(t)$. Absorb the spatially constant functions $C_j$ into the gauges of the velocity potentials. Evaluating at $z=\eta$ and subtracting the equations yields the stated dynamic boundary condition for an inviscid interface. Surface tension supplies the pressure jump, rather than pressure continuity. Thus the first condition enforces material motion, and the second enforces the Bernoulli pressure balance.

Use the linearized boundary condition: evaluate at $z=0$ and discard products of perturbations. The kinematic boundary condition and Young–Laplace equation becomeIncompressible flow and irrotational flow give the Laplace equation for each potential. For a normal mode, the wall conditions selectWith $\eta =B{e}^{i(kx-\omega t)}$, the kinematic conditions give $A_1=i\omega B/[k\sinh (k{L}_1)]$ and $A_2=-i\omega B/[k\sinh (k{L}_2)]$. Hence ${\hat{\varphi }}_1(0)=i\omega B\coth (k{L}_1)/k$ and ${\hat{\varphi }}_2(0)=-i\omega B\coth (k{L}_2)/k$. Substitution into the dynamic condition gives the finite-depth Rayleigh-Taylor dispersion relation:Its denominator is positive. Assuming $g>0$ and $\gamma >0$, negative ${\omega }^2$ therefore produces a growing branch $\omega =i\sigma$ with $\sigma >0$, as in Rayleigh-Taylor instability.

Put $A=g({\rho }_1-{\rho }_2)>0$. In the deep-layer limit, the finite-depth Rayleigh-Taylor dispersion relation gives ${\omega }^2=(\gamma k^3-Ak)/({\rho }_1+{\rho }_2)$. Its sign changes atThe Rayleigh-Taylor instability band is $0<k<k_c$; $k=k_c$ is neutral and $k>k_c$ is oscillatory.

For the growing normal mode in two deep layers,It vanishes at the endpoints of the unstable band and has its unique interior maximum where $A-3\gamma k^2=0$. Using $k_c=\sqrt{A/\gamma }$ givesThis is the maximum exponential growth rate, not the real oscillation frequency.

The stationary-point condition for the squared growth rate is $A-3\gamma k_m^2=0$. Its second derivative is $-6\gamma k_m/({\rho }_1+{\rho }_2)<0$, so

For $L_2\to \infty$, the denominator is ${\rho }_1\coth (k{L}_1)+{\rho }_2>0$. The numerator is unchanged, so the cutoff is exactly independent of depth:This remains exact before taking the thin-layer limit in the finite-depth Rayleigh-Taylor dispersion relation.

Write $A=g({\rho }_1-{\rho }_2)$ and $\epsilon =k_c{L}_1\ll 1$. Throughout the unstable band, $k{L}_1\le \epsilon$ and $\coth (k{L}_1)=1/(k{L}_1)+O(k{L}_1)$. Since ${\rho }_2<{\rho }_1$,The leading polynomial is maximized at $k^2=A/(2\gamma )$, where it equals $L_1{A}^2/(4{\rho }_1\gamma )$. Thus the thin-layer Rayleigh-Taylor growth asymptotics giveThis is a leading asymptotic expression; the exact maximum for a small but nonzero layer depth has corrections.

Maximize the leading thin-layer squared growth rate, $(L_1/{\rho }_1)(A{k}^2-\gamma k^4)$. Its positive stationary point obeys $2Ak-4\gamma k^3=0$ and is the unique interior maximum. ThereforeThe maximizing point stays within the thin-layer regime because $k_m{L}_1\sim k_c{L}_1/\sqrt{2}\ll 1$.

Use superscripts $T$ and $D$ for the thin-upper-layer and two-deep-layer cases. The cutoff is identical, whileThus confinement suppresses the Rayleigh-Taylor instability growth strongly, but moves the most amplified disturbance to a somewhat larger wavenumber. The cutoff depends on the gravity-capillary balance, whereas growth depends also on the inertia of the moving layers.

First suppose $c_i\ne 0$ and put $F=\hat{w}/(U-c)$. The wall condition gives $F(\pm L)=0$. Substitution into the Rayleigh equation for inviscid shear flow givesMultiply by $\overline{F}$ and use integration by parts. The resulting weighted identity for Howard's semicircle theorem isIts imaginary part gives $\int UQ=c_{r}I$, and its real part then gives $\int U^{2}Q=(c_r^2+c_i^2)I$. Set $a=U_{\min }$ and $b=U_{\max }$. Since $(U-a)(b-U)\ge 0$,Completing the square proves Howard's semicircle theorem:In particular every growing normal mode lies in the upper semicircle for $k>0$. If $c$ is real and outside $[a,b]$, $F$ is again regular, but the same identity has a strictly positive integrand unless the mode vanishes. Such a neutral mode is impossible. Consequently neutral regular modes also satisfy the bound; no division by $c_i$ is needed for this last conclusion. The nonreal proof avoids critical-level singularities altogether.

Normalize $U$ and $c$ by $U_{\max }$ and write $h={\zeta }_L\in (0,1)$, $\alpha =kL>0$. In each open region $U^{\prime \prime }=0$, so the Rayleigh equation for inviscid shear flow reduces to ${\hat{w}}^{\prime \prime }-{\alpha }^2\hat{w}=0$ for the nonreal modes of interest. Incorporating the wall conditions from the start, the bounded piecewise-linear shear layer hasThese are four independent regional amplitudes before matching. The same expressions extend to noncritical neutral modes and their limits.

Impermeability gives $\hat{w}(\pm 1)=0$. At $\zeta =\pm h$, the base velocity is continuous but its derivative jumps. The normal velocity and pressure must remain continuous. From incompressibility and the streamwise linearized momentum equation,Therefore, with derivatives now taken in $\zeta$, the vorticity-jump matching for an inviscid shear flow isAt the upper jump $[U^{\prime }]=-1/h$ and at the lower jump $[U^{\prime }]=1/h$, where brackets mean the limit from larger $\zeta$ minus that from smaller $\zeta$. Equivalently $(U-c)[{\hat{w}}^{\prime }]=[U^{\prime }]\hat{w}$ away from $U=c$. Setting $[{\hat{w}}^{\prime }]=0$ would discard the vorticity jumps and give the wrong eigenproblem.

Let $q=\alpha h$, $d=\alpha (1-h)$, and use the four regional amplitudes above. Continuity of $\hat{w}$ at the two jumps givesPressure continuity, using vorticity-jump matching for an inviscid shear flow, gives the other two equations:These are the required homogeneous four equations. As a useful reduction, put $X=\tanh q$, $Y=\tanh d$, $H=q(1+XY)$ and $J=q(X+Y)$. Eliminating $A,D$ leavesThe vanishing determinant gives $c^2=(H-Y)(J-XY)/(HJ)$, which expands to the printed dispersion relation.

For fixed $h\in (0,1)$, $X=\alpha h+O({\alpha }^3)$ and $Y=\alpha (1-h)+O({\alpha }^3)$. The factored dispersion relation therefore giveswhere the latter follows from $X,Y\to 1$ and $c^2\sim (1-1/(2\alpha h){)}^2$. A real negative $c^2$ yields a growing branch $c=i\sqrt{-c^2}$ for positive $k$. Thus $h<1/2$ gives instability at sufficiently small nonzero wavenumber.

For the converse, monotonicity toward the limiting value one would suffice, as permitted in the question. In fact a direct sign argument avoids this extra assumption. With $q=\alpha h>0$, $X=\tanh q<q$, so $J-XY=qX+(q-X)Y>0$. If $h\ge 1/2$, then $q\ge d=\alpha (1-h)$ and $Y=\tanh d<d\le q$. Hence $H-Y=q(1+XY)-Y>0$, so $c^2>0$ for every $\alpha >0$. We obtain the exact instability threshold for a bounded piecewise-linear shear layer:At ${\zeta }_L=1/2$ the long-wave limiting value is zero, but every finite positive wavenumber has $c^2>0$.

Choose a length $\ell$, velocity $U_s$, time $\ell /U_s$, and pressure $\rho U_s^2$, giving Reynolds number $Re=U_s\ell /\nu$. Write perturbation velocity as $(u,v,w)$ and let $\mathcal{L}={\partial }_t+U{\partial }_x-R{e}^{-1}{\nabla }^2$. Linearizing the incompressible Navier-Stokes equations givesTaking the divergence gives ${\nabla }^{2}p=-2U^{\prime }{v}_x$. Applying ${\nabla }^2$ to the wall-normal momentum equation and usingthen cancels the pressure derivatives and yields $\mathcal{L}{\nabla }^{2}v-U^{\prime \prime }{v}_x=0$. Define wall-normal vorticity $\eta =u_z-w_x$. Applying ${\partial }_z$ to the streamwise equation minus ${\partial }_x$ to the spanwise equation gives $\mathcal{L}\eta =-U^{\prime }{v}_z$.

For the assumed normal modes, set $D=d/dy$ and ${\kappa }^2={\alpha }^2+{\beta }^2$. The two equations become the Orr-Sommerfeld equation and Squire equation:At rigid no-slip walls the conditions are $\hat{v}=D\hat{v}=\hat{\eta }=0$ for nonzero horizontal wavenumber.

An Orr-Sommerfeld mode has $\hat{v}\ne 0$ satisfying the Orr-Sommerfeld equation, with accompanying vorticity satisfying the forced Squire equation. A Squire mode has $\hat{v}=0$ and $\hat{\eta }\ne 0$, so its Squire equation is homogeneous.

For a Squire mode multiply that homogeneous equation by $\overline{\eta }$ and integrate over the channel. Under the no-slip condition $\eta =0$ at the walls, integration by parts givesSince $U,\alpha$ are real, taking the real part provesThe numerator is strictly positive for a nontrivial finite-channel mode, and $Re>0$. Thus Squire modes are exponentially damped in the convention $e^{-i\omega t}$. This does not rule out transient growth from non-normal modes in the coupled system.

With $U=U_0$ and $\beta =0$, the Orr-Sommerfeld equation factors into constant-coefficient operators:For distinct nonzero characteristic roots, the general solution isEither square-root choice for $\gamma$ gives the same solution space. The no-slip boundary condition and impermeability impose $\hat{v}=D\hat{v}=0$ at each wall.

The displayed exponential basis needs the usual repeated-root qualification. If ${\gamma }^2={\alpha }^2\ne 0$, replace it by $(a_1+a_{2}y)e^{\alpha y}+(a_3+a_{4}y)e^{-\alpha y}$. If just one root pair is zero, that pair contributes $1,y$; if both pairs are zero, the basis is $1,y,y^2,y^3$. These are the complete limiting cases, not four independent copies of a repeated exponential.

Take $k^{\ast }>0$ and positive shear scale $S^{\ast }$, with kinematic viscosity ${\nu }^{\ast }$. ChooseThe velocity scale is $S^{\ast }/k^{\ast }$. A reversed shear can be treated by reversing the corresponding coordinate orientation and mode convention. The dimensional Orr-Sommerfeld equation, with $U^{\ast \prime \prime }=0$, isSubstitution of the scales and cancellation gives

In two dimensions at nondimensional streamwise wavenumber one, incompressibility gives $i\hat{u}+D\hat{v}=0$, so $\hat{u}=iD\hat{v}$. The spanwise vorticity is the $z$ component of the curl:ConsequentlyThe symbol $\hat{\omega }$ here denotes a vorticity amplitude, not the temporal frequency used elsewhere. In dimensional variables, ${\hat{\zeta }}^{\ast }=-i(D_{\ast }^2-(k^{\ast }{)}^2){\hat{v}}^{\ast }/k^{\ast }$, with the same sign convention.

Put $W=(D^2-1)\hat{v}$ and ${\delta }^3=R{e}^{-1}$. The constant-shear Orr-Sommerfeld equation givesLet $a=e^{i\pi /6}$ and $Z=a(y-c-i{\delta }^3)/\delta$. Since $d/dy=(a/\delta )d/dZ$ and $a^3=i$, the equation becomesThis is the Airy ordinary differential equation, yielding the Airy reduction of the Orr-Sommerfeld equation in constant shear. It is distinct from the evolution PDE represented by the existing Airy equation article.

Write $W(s)=a_1\mathrm{Ai}(Z(s))+a_2\mathrm{Bi}(Z(s))$, using the independent Airy functions. To solve $(D^2-1)\hat{v}=W$, the variation of parameters formula, equivalently a direct convolution, gives the velocity reconstruction from the Orr-Sommerfeld vorticity:The reference point $y_0$ is arbitrary. Twice differentiating the integral produces $W(y)$ because $\sinh 0=0$ and its derivative at zero is one, verifying the formula. The four constants are independent: applying $D^2-1$ to an identically zero combination first forces $a_1=a_2=0$, then the complementary solutions force $C_{+}=C_{-}=0$. The boundary conditions would be $\hat{v}(y_j)=D\hat{v}(y_j)=0$, $j=1,2$. Their compatibility selects the eigenvalues; they need not be imposed here.

The matrix exponential is the absolutely convergent power seriesIt converges in any finite-dimensional operator norm, since $\Vert L^m\Vert \le \Vert L{\Vert }^m$. Termwise differentiation gives $B^{\prime }=LB$, $B(0)=I$, so the solution is $q(t)=B(t)q_0$. Also $e^{Lt}{e}^{-Lt}=I$, proving invertibility without requiring $L$ itself to be invertible.

The question's expansion in eigenvectors tacitly needs an eigenbasis: an invertible matrix alone need not be a diagonalizable matrix. The exponential solution and the singular-value arguments remain valid without that assumption.

The relevant norm is the operator norm induced by the complex Euclidean norm, also called the matrix 2-norm:The maximum exists because the unit sphere in finite dimensions is compact. It obeys $\Vert Bq{\Vert }_2\le \Vert B{\Vert }_2\Vert q{\Vert }_2$ and is the norm relevant to the specified perturbation energy. Other matrix norms do not in general give this gain.

A normal matrix admits unitary diagonalization of a normal matrix: $L=V\mathrm{diag}({\lambda }_j)V^{†}$ with orthonormal eigenvectors. ThusFor $t\ge 0$, the largest weight belongs to an eigenvalue with maximal real part. The weighted average is bounded by that weight and attains it on the corresponding eigenspace. Therefore the optimal energy amplification of a linear system isIf the largest real part is repeated, any nonzero combination within that maximal eigenspace is optimal. For negative times the ordering reverses; the displayed formula concerns forward evolution.

The singular value decomposition of $B$ is $B=U\Sigma V^{†}$, where $U,V$ are unitary and $\Sigma =\mathrm{diag}({\sigma }_1,\dots ,{\sigma }_n)$, ${\sigma }_1\ge \cdots \ge {\sigma }_n>0$ since $B$ is invertible. Its right singular vector $v_j$ and left singular vector $u_j$ are the corresponding unit columns of $V$ and $U$, withEquivalently $B^{†}B{v}_j={\sigma }_j^2{v}_j$ and $B{B}^{†}{u}_j={\sigma }_j^2{u}_j$. The positive numbers ${\sigma }_j$ are the singular values. Degenerate singular subspaces allow any orthonormal choice of paired vectors.

For $q(t)=B{q}_0$, the energy ratio is the Rayleigh quotient of the Hermitian positive-definite matrix $B^{†}B$:Expanding $q_0$ in the orthonormal right singular vectors makes this a weighted average of ${\sigma }_j^2$. HenceAny nonzero vector in the largest right-singular subspace is optimal if ${\sigma }_1$ is repeated. Its amplified state is proportional to the paired left singular vector. For a non-normal matrix $L$, this optimal initial direction need not be an eigenvector of $L$, explaining transient growth from non-normal modes.

A normal mode $x\propto e^{-i\omega t}$ has eigenvalue $\lambda =-i\omega$ of $A$. Setting $d_{\omega }=\omega -{\omega }_1$, the characteristic equation isThusOne branch then has positive imaginary frequency and hence positive exponential growth rate. If $4ac<b^2$, both frequencies are real and the distinct-eigenvalue system has bounded oscillatory motion. At equality there is no positive exponential growth rate, but a nontrivial Jordan block can cause algebraic growth, so exponential neutrality must not be confused with Lyapunov stability. The parameter $d$ listed in the PDF does not appear in the printed matrix and has no effect on this calculation.

Put $q=10/Re>0$. The frequency formula becomesFor large positive Reynolds number, $\sqrt{q-q^2/4}=\sqrt{q}[1-q/8+O(q^2)]$. ThereforeMore precisely the real correction is $-5/Re$, and the next imaginary correction is $O(R{e}^{-3/2})$. The growing branch has a rate tending to zero as $Re\to \infty$, despite being unstable at every sufficiently large finite $Re$.

For $q=10/Re$, the eigenvalues areFor physical $Re>0$, this gives the following classification. If $Re>5/2$, then $0<q<4$ and one eigenvalue has positive real part: exponential instability. If $0<Re<5/2$, then $q>4$, both eigenvalues are distinct and imaginary, and solutions are bounded oscillatory combinations. At $Re=5/2$ there is a repeated imaginary eigenvalue ${\lambda }_0=-i{\omega }_1+2i$. The matrix $N=A-{\lambda }_{0}I=\left(\begin{array}{cc}-2i& 1\\ 4& 2i\end{array}\right)$ is nonzero and satisfies $N^2=0$, so $e^{At}=e^{{\lambda }_{0}t}(I+tN)$ has algebraically growing trajectories.

Thus the exponential instability threshold is $Re>5/2$, consistently with $4ac>b^2$, which here is $4q>q^2$. At the threshold the matrix is defective, and the eigenvector-basis assumption from part (a) cannot be used; the matrix-exponential and singular-value formulas still apply.

A normal matrix satisfies $A{A}^{†}=A^{†}A$. Direct multiplication, with all parameters real, gives the normality criterion for a two-mode shear model:For $a=1$, $b=c=q=10/Re$, both $1-q^2=0$ and $q(1+q)=0$ are required. Their only common solution is $q=-1$. HenceThe original PDF really has these signs and $b=c=10/Re$. Thus parts (v) and (vi) have no physical positive-Reynolds-number case as printed. We give their algebraic continuation at $Re=-10$ below, without silently changing the matrix. The tempting value $Re=10$ makes the off-diagonal magnitudes equal but fails the off-diagonal commutator condition.

At the formal value $Re=-10$, $q=-1$ and $A=-i{\omega }_{1}I+\left(\begin{array}{cc}0& 1\\ -1& -i\end{array}\right)$ is a skew-Hermitian matrix. Put $r_{\pm }=(-1\pm \sqrt{5})/2$. Its two eigenvalues and normalized eigenvectors areSubstitution uses $r_{\pm }^2+r_{\pm }-1=0$ and verifies $A{v}_{\pm }={\lambda }_{\pm }{v}_{\pm }$. Their Hermitian inner product issince $r_{+}{r}_{-}=-1$. This provides the requested orthogonality in the formal continuation. There is no corresponding physical $R{e}_N>0$ for the printed system.