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

Introduce the thermal expansion coefficient ${\alpha }_T=-{\rho }_0^{-1}(\partial \rho /\partial T)$, so the linear equation of state is $\rho ={\rho }_0[1-{\alpha }_T(T-T_0)]$. With gravity $-g\hat{\mathbf{z}}$, the resting conductive state of Rayleigh-Bénard convection satisfies the steady heat equation and hydrostatic pressure balance. Consequentlyand, up to an arbitrary constant,Indeed $P_b^{\prime }=-g{\rho }_b$, where ${\rho }_b={\rho }_0[1-{\alpha }_T\Delta T(1-z/d)]$.

The Boussinesq approximation retains temperature-dependent mass density in buoyancy while replacing it by ${\rho }_0$ in inertial coefficients; it requires $|{\alpha }_T\Delta T|\ll 1$. Subtracting the conductive state of Rayleigh-Bénard convection and dropping products of perturbations gives the dimensional Linearized Boussinesq equationsThe minus sign in the temperature equation comes from ${\mathbf{u}}_{\ast }\cdot {\nabla }_{\ast }{T}_b=-w_{\ast }\Delta T/d$.

Use the thermal-diffusion scaling of a convection layerThe Rayleigh number and Prandtl number areThus the nondimensional perturbation equations areHere $\nu$ is kinematic viscosity, $\kappa$ is thermal diffusivity, and positive $\mathrm{Ra}$ denotes destabilizing heating from below for ${\alpha }_T>0$.

A direct elimination of pressure gives the same result as taking a double curl of the momentum equation. Incompressibility gives ${\nabla }^{2}p=\mathrm{Ra}\mathrm{Pr}{\partial }_z\theta$. Apply ${\nabla }^2$ to the vertical momentum equation and subtract this pressure derivative to obtainHere ${\nabla }_h^2={\partial }_x^2+{\partial }_y^2$ is the horizontal Laplacian. At $z=0,1$, impermeability gives $w=0$ and fixed temperature gives $\theta =0$. The stress-free boundary condition gives ${\partial }_{z}u={\partial }_{z}v=0$, since horizontal derivatives of the boundary value $w=0$ vanish. Differentiating incompressibility once in $z$ therefore givesFor a horizontal Fourier mode $X=e^{i(k_{x}x+k_{y}y)}$, ${\nabla }_h^{2}X=-{\lambda }^{2}X$ with real $\lambda =\sqrt{k_x^2+k_y^2}$. More generally the nonnegative spectrum of $-{\nabla }_h^2$ follows by integration by parts for periodic or square-integrable horizontal disturbances. Substituting the separated normal mode gives, with $D=d/dz$,The boundary conditions admit a complete Fourier sine series in $z$. For a nonzero horizontal wavenumber, take $W=W_n\sin (n\pi z)$, $\Theta ={\Theta }_n\sin (n\pi z)$, and set $a_n^2={\lambda }^2+n^2{\pi }^2$. The two amplitude equations areThus the stress-free convection growth-rate polynomial isIts discriminant is $(\mathrm{Pr}-1{)}^2{a}_n^4+4\mathrm{Ra}\mathrm{Pr}{\lambda }^2/a_n^2$, so both growth rates are real for $\mathrm{Ra}>0$. For $\mathrm{Ra}<0$, its two coefficients after the leading term are positive: their sum gives ${\sigma }_1+{\sigma }_2=-(1+\mathrm{Pr})a_n^2<0$, and their product is $\mathrm{Pr}(a_n^4-\mathrm{Ra}{\lambda }^2/a_n^2)>0$. Real roots are both negative; nonreal roots have real part $-(1+\mathrm{Pr})a_n^2/2<0$. Hence stable heating from above damps every mode. Horizontally uniform thermal modes and uncoupled horizontal viscous modes are also purely damped; impermeability and incompressibility exclude a horizontally uniform nonzero $w$.

Because the unstable-side roots are real, the onset is stationary, an exchange of stabilities in stress-free convection. Setting $\sigma =0$ gives the free-slip convection neutral curveWith $x={\lambda }^2$, its derivative has the sign of $2x-n^2{\pi }^2$. Its minimum is at $x=n^2{\pi }^2/2$, with value $27n^4{\pi }^4/4$. Minimizing also over $n\ge 1$ gives the critical values

Use the Boussinesq approximation with a constant inertial reference mass density ${\rho }_0$, keeping the background mass density gradient in buoyancy. Write $V(z)=\overline{U}(z)-c$ and $P=\hat{p}/{\rho }_0$. For $k\ne 0$, the linearized normal mode equations areThe reference mass density is chosen close to the background values, whose relative variation must be small. Hydrostatic background pressure satisfies ${\overline{p}}^{\prime }=-g\overline{\rho }$.

From incompressibility, $\hat{u}=i{\hat{w}}^{\prime }/k$. The horizontal momentum equation and the density equation then giveSubstitute these into the vertical momentum equation and multiply by $ik$:For $V\ne 0$, rearrangement gives the Taylor–Goldstein equationThe buoyancy frequency $N$ is real for stable stratification, where ${\overline{\rho }}^{\prime }<0$. Points where $\overline{U}=c$ are critical levels of an internal gravity wave; the derivation there is understood through a limiting or piecewise formulation. The printed density gradient is ${\overline{\rho }}^{\prime }$, not the pressure gradient accidentally substituted in the local TeX.

Let $\hat{\eta }$ be the displacement of a material interface and use $[f{]}_{-}^{+}=f_{+}-f_{-}$. On each side, the kinematic boundary condition isThere is one displaced material interface, so its displacement is continuous. Thus the first jump condition isIn particular, vertical velocity need not be continuous when the background velocity jumps.

With no surface tension, pressure continuity applies on the displaced interface, rather than at its undisplaced height. Linearizing gives $[\hat{p}+\hat{\eta }{\overline{p}}^{\prime }{]}_{-}^{+}=0$. Since ${\overline{p}}^{\prime }=-g\overline{\rho }$, this becomes $[\hat{p}-g\overline{\rho }\hat{\eta }{]}_{-}^{+}=0$. The horizontal momentum calculation in part (a) givesSubstitute this and $\hat{\eta }=\hat{w}/[ik(\overline{U}-c)]$ to obtain the second jump condition:These jump conditions for stratified inviscid shear flow use one-sided values of ${\overline{U}}^{\prime }$ and apply to jumps of mass density, vorticity, or velocity. If $\overline{U}$ is continuous, the first condition reduces to continuity of $\hat{w}$; a jump of vorticity generally still prevents continuity of ${\hat{w}}^{\prime }$. This derivation avoids multiplying singular distributional derivatives of a discontinuous velocity in the Taylor–Goldstein equation.

Take $k>0$, without loss for this symmetric real base state, and write $s=1+J$, $e=e^{-2\alpha }$, and $y={\tilde{c}}^2$. The given quartic becomes a real quadratic $y^2+Ay+B=0$. Its constant term factors asTherefore $B<0$ exactly whenA negative product means the two quadratic roots are real and of opposite sign: the discriminant is $A^2-4B>0$. The negative root gives a pair $\tilde{c}=\pm i\sqrt{-y_{-}}$. One sign has positive temporal growth rate $k{c}_i$ and is unstable. Using the definitions of the hyperbolic functions gives the required unstable band:This proves the unstable band of a three-layer stratified shear flow directly, without needing to solve all four quartic roots.

For the wave interpretation, consider either interface in isolation with decaying normal modes $\hat{w}=w_0{e}^{-k|z-z_i|}$. Put $S=\Delta U/h$, $U_i=\overline{U}(z_i)$, and $v=c-U_i$, the intrinsic phase velocity. At either interface $[\overline{\rho }]=-\Delta \rho /2$, while $[{\overline{U}}^{\prime }]=-S$ at the upper interface and $+S$ at the lower one. The jump conditions for stratified inviscid shear flow give the gravity-vorticity interface wave relationFor stable density jumps, $J\ge 0$, both intrinsic branches are real. In units $C=2c/\Delta U$, the isolated wave speeds areThe upper wave travelling against its positive background current has $C_u=1-[1+\sqrt{1+8\alpha J}]/(4\alpha )$. The lower wave travelling against its negative background current has the opposite speed, $C_l=-C_u$. Their counterpropagating wave resonance in a three-layer shear flow occurs at $C_u=C_l=0$, which givesAt large wavenumber, coupling across the separation $h$ is exponentially weak, of order $e^{-kh}=e^{-2\alpha }$. The unstable band becomeswith full width $4\alpha e^{-2\alpha }+O(\alpha e^{-6\alpha })$. It is a counterpropagating wave instability: the two waves propagate oppositely relative to their local currents but have almost equal laboratory phase velocities, allowing weak coupling to lock their phases and extract mean-flow energy. The unstable band of a three-layer stratified shear flow thus becomes narrowly concentrated around this resonance. For fixed $J$, arbitrarily large $k$ lies outside that band; large-$k$ resonance requires $J$ to grow with $\alpha$.

Under the convention $e^{i(kx-\omega t)}$, a mode's amplitude grows as $e^{{\omega }_{i}t}$. The linear stability analysis of real-wavenumber disturbances is therefore summarized byAt ${\omega }_{i,\max }=0$ there is marginal or neutral hydrodynamic stability: no normal mode grows exponentially, but some may fail to decay. If “linearly stable” includes neutral modes, its condition is ${\omega }_{i,\max }\le 0$. These modal statements do not rule out transient growth from non-normal modes or algebraic growth.

The fixed observer is the ray $V=0$. The physically selected zero-group velocity saddle defines the absolute wavenumber $k_0$ and absolute frequency ${\omega }_0$:For a smooth frequency branch with $D_{\omega }\ne 0$, the saddle condition is equivalently $D_k(k_0,{\omega }_0)=0$. The absolute growth rate isHere $k_0$ and ${\omega }_0$ can be complex. A physical saddle point is selected by the localized impulse response and deformation of its Fourier transform contour; an arbitrary algebraic stationary point need not determine absolute hydrodynamic instability.

For a temporally unstable base flow, the relevant saddle growth rate along a ray distinguishes two cases in the chosen observation frame:A convective hydrodynamic instability amplifies a travelling wave packet but lets a localized disturbance decay at any fixed point. An absolute hydrodynamic instability grows at a fixed point. The separating case $\sigma (0)=0$ is marginal in exponential rate and can have an algebraic prefactor. In a frame moving with $V$, replace $\sigma (0)$ by $\sigma (V)=\mathrm{Im}({\omega }_{\star }-V{k}_{\star })$: this distinction depends on the observation frame.

In the linear complex Ginzburg-Landau equation, take $U,c_d,\mu$ real. $U$ gives advection, carrying the envelope downstream at speed $U$. $c_d$ gives wave dispersion: the imaginary part of the second-derivative coefficient changes wave phase without directly dissipating a real-wavenumber Fourier mode. $\mu$ gives local linear production or damping: a spatially uniform perturbation grows for $\mu >0$ and decays for $\mu <0$.

The real part of the diffusion coefficient is the scaled value one. It smooths the envelope and damps a real-wavenumber Fourier mode at rate $k^2$. Thus the linear complex Ginzburg-Landau equation models competing advection, smoothing through a diffusion equation, wave dispersion and local growth, rather than nonlinear saturation.

For a normal mode, ${\partial }_t\mapsto -i\omega$, ${\partial }_x\mapsto ik$, and ${\partial }_x^2\mapsto -k^2$. Substitution into the linear complex Ginzburg-Landau equation givesTherefore the dispersion relation isFor real $k$, its temporal growth rate is ${\omega }_i=\mu -k^2$, with unique maximum $\mu$ at $k=0$. Its real angular frequency is $Uk+c_d{k}^2$, displaying the separate advection and wave dispersion terms.

The dispersion relation gives the zero-group velocity saddleThus the absolute growth rate is $\sigma (0)=\mu -U^2/[4(1+c_d^2)]$. More generally the saddle growth rate along a ray isso a frame travelling at $V=U$ sees the maximal rate $\mu$.

For $U>0$, the three open regions areThe line $\mu =0$ is the temporal marginal boundary; the parabola $\mu =U^2/[4(1+c_d^2)]$ is the convective-to-absolute boundary. The Green function of the linear complex Ginzburg-Landau equation has the prefactor $t^{-1/2}$, so zero exponential rate on the latter boundary still allows algebraic decay of the impulse response. This is the stability diagram of the linear complex Ginzburg-Landau equation.

Impermeability requires $\hat{w}(-1)=\hat{w}(1)=0$. Let $c=\omega /\alpha$ with real $\alpha \ne 0$, and set $a=|\alpha |$. If $c$ is nonreal, or real outside the channel, the Rayleigh equation for inviscid shear flow gives ${\hat{w}}^{\prime \prime }-a^2\hat{w}=0$ everywhere. The two wall conditions then force $\hat{w}=0$. In fact there is no nonzero globally smooth eigenfunction for any $c$: the complete family is the inviscid Couette continuous spectrum.

For each $\xi \in (-1,1)$, permit a localized vorticity sheet. Let $G_a$ be the Dirichlet Green function satisfying $(D^2-a^2)G_a(z,\xi )=\delta (z-\xi )$. An explicit generalized eigenfunction iswhere $z_{<}=\min (z,\xi )$ and $z_{>}=\max (z,\xi )$. It vanishes at both walls, is continuous at $z=\xi$, and has derivative jump $[{\partial }_z{G}_a{]}_{{\xi }^{-}}^{{\xi }^{+}}=1$. Consequentlyusing the Dirac delta multiplication identity. These are vorticity-sheet eigenfunctions of inviscid Couette flow, understood as generalized eigenfunctions, not as a discrete smooth Sturm-Liouville eigenfunction expansion.

To see completeness, define the vorticity variable $q=(D^2-a^2)w$. Its evolution is $({\partial }_t+i\alpha z)q=0$, so for any admissible initial vorticity $q_0$,The homogeneous Dirichlet problem for $D^2-a^2$ has only the zero solution, so this inversion reconstructs every initial vertical-velocity field in its usual function space. The generalized frequencies fill the interval with endpoints $\pm \alpha$; the endpoint values are understood as the closure of the continuous spectrum.

Every generalized frequency in the inviscid Couette continuous spectrum is real: ${\omega }_{\xi }=\alpha \xi$ for real $\alpha ,\xi$. Thus every time factor has unit modulus and there is no exponentially growing normal mode:The vorticity equation also gives $q(z,t)=q_0(z)e^{-i\alpha zt}$, preserving its $L^2$ norm. For fixed nonzero $\alpha$, the bounded Dirichlet Green function inversion therefore gives a uniform bound on the reconstructed velocity in terms of the initial vorticity norm. Superpositions can still rearrange their velocity energy and display transient growth from non-normal modes; neutral hydrodynamic stability does not require every velocity component to decrease monotonically. Treating only smooth discrete modes would miss this continuous neutral family.

Use the kinetic-energy inner product $⟨v,u⟩=\int v^{\ast }\cdot udV$, with perturbation energy $\Vert u{\Vert }^2/2$. On divergence-free fields the primal generator, before the pressure projection, isThe base state is also divergence-free. The periodic boundary terms in $x,z$ cancel, and decay at $|y|\to \infty$ removes the remaining boundary terms. Integration by parts changes $-\mathbf{U}\cdot \nabla$ to $+\mathbf{U}\cdot \nabla$; the Laplacian is self-adjoint under these conditions. Writing $J_{ij}={\partial }_j{U}_i$, the shear term $-Ju$ has adjoint operator $-J^{T}v$. A pressure gradient pairs to zero with a divergence-free field. Hencewhere the last term enforces the divergence-free projection.

To conserve $⟨v(t),u_p(t)⟩$ over the finite interval, the physical-time adjoint obeys $-{\partial }_{t}v=A^{†}(t)v$. Put $\tau =-t$ and $u_d(\tau )=v(-\tau )$; then ${\partial }_{\tau }{u}_d=A^{†}(-\tau )u_d$. This reverses the interval from $t=T$ towards $t=0$, or $\tau =-T$ towards $\tau =0$.

For $\mathbf{\Omega }=\nabla \times \mathbf{U}$, the vorticity cross-product identity and the curl of a cross producthold when both fields are divergence-free. Their difference equals the transport and shear terms in the adjoint operator. Therefore the adjoint linearized Navier-Stokes evolution isThe Leray-Helmholtz projection determines the adjoint pressure. The time-dependent base state is evaluated at $-\tau$; no extra time derivative of that base state appears in the spatial adjoint operator.

The matrix $N$ is a skew-symmetric matrix, so $x^{T}Nx=0$ for real $x$. Therefore the energy-neutral rotational nonlinearity contributes nothing directly to the energy balance:The nonlinear term rotates the state without changing its Euclidean norm. It can still change the orientation, and thereby affect later energy growth through the linear term. Thus energy cancellation does not mean the nonlinear and linear trajectories have identical energy histories.

For $x_1\ne 0$, the energy balance isThe concave quadratic changes sign at the roots of ${\rho }^2-4R\rho +4=0$. For $R>1$, the neutral orientations and the growth interval areOutside this orientation interval for transient energy growth, the energy decays; at either endpoint it is instantaneously stationary. The direction $x_1=0$ also has $\dot{E}=-x_2^2/(4R)<0$.

In the linear problem, the orientation evolves according to $\dot{\rho }=\rho [3/(4R)-\rho ]$. Since both ${\rho }^{\pm }>3/(4R)$, a trajectory traversing the positive growth sector has decreasing orientation: growth begins at ${\rho }^{+}$ and ends at ${\rho }^{-}$, as in the question.

For the full nonlinear model, however,The added rotation can reverse the direction of crossing. Thus the two roots always delimit energy growth, but the printed “begins” and “ends” labels are conditional on a decreasing-orientation passage, and are not universal for arbitrary nonlinear initial conditions.

In the linear problem, ${\dot{x}}_2=-x_2/(4R)$, so $x_2(t)=x_{20}{e}^{-t/(4R)}$. Solve ${\dot{x}}_1+x_1/R=x_2(t)$ with an integrating factor to obtain the complete solution:Equivalently the matrix exponential isBoth eigenvalues of $L$ are negative, so each solution eventually decays, despite the possible transient growth generated by $b$.

For the linear evolution from part (iii), maximize the energy ratio at each time. The optimal energy amplification of a linear system isThus the optimizing initial state is a right singular vector, not generally an eigenvector of $L$.

There is a convenient exact way to find the maximizing time. Let $x_0$ be a unit optimizing initial vector and $y=M(t)x_0$. The singular value decomposition gives $M^{T}y=G_{\mathrm{opt}}{x}_0$. For $t>0$, the top singular value is simple. At an interior time maximum,Since $LM=ML$, this also equals $G_{\mathrm{opt}}{x}_0^{T}L{x}_0$. Thus both initial and final orientations are neutral energy-growth directions. The positive entries of $M^{T}M$ let the top right singular vector be chosen with positive components. Its orientation decreases through the growth sector, so it must start at ${\rho }^{+}$ and end at ${\rho }^{-}$.

Put $q=a/d=e^{-3t/(4R)}$ and $K=4R/3$. The solution givesSetting this equal to ${\rho }^{-}$ yieldsTherefore the exact optimal time for $R>1$ isHere and below $\log$ is the natural logarithm. This is the unique positive stationary time: $q$ decreases monotonically from one to zero, and the neutral-direction condition determines one $q_{\ast }\in (0,1)$. Also $G_{\mathrm{opt}}(0)=1$, $G_{\mathrm{opt}}(\infty )=0$, and $R>1$ admits initial energy growth, so this stationary point is the global maximum. For $R\le 1$, the instantaneous energy-growth criterion for a linear system gives no amplification, and the maximum is one at time zero.

At the optimum, $d=q_{\ast }^{1/3}$ and the initial direction is $(1,{\rho }^{+})$. Since the final direction has ratio ${\rho }^{-}$,Using ${\rho }^{\pm }=2R(1\pm s)$, $R^2(1-s^2)=1$, and $q_{\ast }=(5-3s)/(5+3s)$ simplifies this to the exact maximum gain:As $R\to \infty$, $s=1-1/(2R^2)+O(R^{-4})$, so $q_{\ast }\to 1/4$ and $(1+s)/(1-s)\sim 4R^2$. Consequently the requested large-Reynolds-number transient growth scaling isThese are the optimal time and gain of a Reynolds-scaled triangular model. The leading gain can also be read from the off-diagonal matrix exponential entry: with $\tau =t/R$, $b=(4R/3)(e^{-\tau /4}-e^{-\tau })$ is maximized at $\tau =4\log 4/3$ and has leading maximum $R{4}^{-1/3}$.

Both eigenvalues of the linear generator, $-1/R$ and $-1/(4R)$, are negative. Nevertheless the non-normal matrix transfers the slowly decaying second component into the first before viscous decay wins. Modal decay coexists with energy amplification of order $R^2$ over a time of order $R$.

This illustrates transient growth from non-normal modes: large Reynolds number permits much greater disturbance amplification even though the linear stability analysis finds no unstable eigenvalue. The energy-neutral rotational nonlinearity can redirect the amplified state and affect its subsequent evolution. The linear gain alone establishes neither sustained nonlinear instability nor a nonlinear transition threshold.