Past exam of the mathematics course of the University of Cambridge / 2015 / iii / Paper 74 (Fluid Dynamics of the Environment)

Use the inviscid, nonrotating Boussinesq approximation, with a stable background mass density $\bar{\rho }(z)$ and reference density ${\rho }_{\ast }$. Write $B=-g{\rho }^{\prime }/{\rho }_{\ast }$ for the perturbation buoyancy and $\pi =p^{\prime }/{\rho }_{\ast }$ for the kinematic pressure. The buoyancy frequency is $N_0^2=-(g/{\rho }_{\ast }){\bar{\rho }}_z>0$. Dropping products of perturbations in the Boussinesq equations gives the Linearized Boussinesq equationsThe last equation expresses incompressible flow; the second follows by advecting the background mass density gradient. Neglecting rotation and viscosity is part of this internal gravity wave model.

For a plane internal gravity wave proportional to $\exp [i({\mathbf{k}}_h\cdot {\mathbf{x}}_h+mz-\omega t)]$, put $\kappa =|{\mathbf{k}}_h|>0$ and $K^2={\kappa }^2+m^2$. Eliminating the horizontal velocity, pressure and buoyancy from the linear equations givesFor example, the horizontal momentum and continuity equations give $\pi =-\omega mw/{\kappa }^2$; substituting $B=-i{N}_0^{2}w/\omega$ into vertical momentum yields the displayed dispersion relation. Thus an IGW has $0<|\omega |\le N_0$ in this model.

On the positive-frequency branch, the phase velocity normal to a constant-phase plane and the group velocity areConsequently $({\mathbf{k}}_h,m)\cdot {\mathbf{c}}_g=0$, and phase and group velocity are perpendicular. Equivalently, the dispersion relation is homogeneous of degree zero in the wave vector, so differentiating with respect to its scale proves the same orthogonality. Energy travels with the group velocity, along the phase planes. At the degenerate limit $m=0$, $\omega =N_0$ and the group velocity vanishes; the orthogonality statement then has this limiting interpretation.

For a two-dimensional internal gravity wave with fixed horizontal wave number $k\ne 0$, the vertical velocity amplitude satisfiesThis follows from the Linearized Boussinesq equations just as in part (a), allowing the background buoyancy frequency to vary. The WKB approximation gives locallyIt requires slowly varying stratification compared with the local vertical wavelength, in particular $|m_z|\ll m^2$, and a small-amplitude wave in the inviscid, nonrotating regime. Propagation requires $N>|\omega |$; where $N<|\omega |$, the local vertical wave number is imaginary and the solution is evanescent. Turning levels and abrupt changes need a separate connection calculation. Finite viscosity, background flow, rotation or nonlinear effects require modifications of this model.

The ray tracing equations use the local dispersion relation as a Hamiltonian:The stationary, horizontally uniform background conserves the frequency and horizontal wave number. In two dimensions the group velocity gives $dx/dz=-m/k$. For the specified buoyancy frequency and $\omega =1$, $m^2=256k^2{z}^2$. Hence for $z>0$ the two geometric internal-wave rays through the specified point areEach curve has both directed propagation senses, giving four local directed rays. Their formal intersections with $z=0$ are $x=-8$ and $x=8$.

At $z=0$, $N=\omega$ and $m=0$: this is a double-zero internal-wave turning level, since $m^2$ touches zero without changing sign. The ray becomes vertical geometrically, but its group velocity tends to zero. Taking $k>0$ without loss of generality, the downward branch hasThe formal travel time diverges logarithmically as $z\to 0$. Meanwhile $|m_z|/m^2=1/(16|k|z^2)$ grows, so the WKB approximation fails before that limit; its predicted $|m{|}^{-1/2}$ amplitude cannot be extrapolated to infinity. The spatial ray reaches a finite limiting position, but finite-wavelength behaviour near the zero must be found beyond WKB. There is no evanescent half-space supplied by this particular $N(z)$, so ordinary simple-turning-point reflection cannot be assumed solely from the local ray construction.

Use the spatial internal-wave ray tracing model, with $dx/dz=\pm 16z$, and reflect the rays so that the frequency and tangential wave number are preserved. At a vertical wall the horizontal energy direction reverses; at a horizontal wall the vertical energy direction reverses. On the sloping wall the reflection depends on internal-wave slope criticality. The three-reflection orbit below encounters that wall at a height above $3/16$, so the wall reverses the horizontal energy direction there.

First take the packet to leave the bottom upward and to the left, with $0<x_0<6$. The four segments of the loop areThe left-wall reflection is at $z=\sqrt{x_0/8}$ and the top reflection at $x=8-x_0$. The slope intersection solves $8z_s^2+3z_s+x_0-10=0$, giving the internal-wave ray return mapThus the next formal bottom reflection is at $(R(x_0),0)$. On $0<x_0<6$, the slope heights are between approximately $0.544$ and $0.946$, consistent with the assumed reflection type.

An internal-wave attractor is a stable periodic ray, so solve $R(x_{\ast })=x_{\ast }$. This gives $z_s=2/3$, and thereforeThe successive vertices areThis is clockwise circulation in the $(x,z)$ plane. Differentiating the internal-wave ray return map showsNeighbouring bottom intersections converge to the periodic orbit, geometrically focusing the energy beam.

For completeness, the opposite upward launch direction has the reversed three-reflection itinerary: slope, top, left. In the interval where this itinerary is valid, putThe three-reflection reverse map is $R^{-1}$, with domain $R(0)<x_0<R(6)$. Its fixed point is the same, but $R_{\mathrm{rev}}^{\prime }(x_{\ast })=41/23>1$. Counterclockwise packets defocus from that orbit. Except for the exactly periodic reverse ray, they eventually leave this itinerary. In the full geometric reflection construction, omitted-side or additional-slope reflections redirect generic rays into the clockwise attracting branch. Tangency and exactly critical reflections are singular exceptions requiring finite-wavelength or dissipative treatment.

The bottom reflection is a formal spatial-ray idealization: part (b) shows that the WKB approximation itself fails near $z=0$, so the diagram does not imply arrival and reflection there in finite WKB travel time. In a real fluid, viscous attenuation of an internal-wave beam, finite wavelength, internal-wave breaking and ensuing turbulence and mixing limit the beam width and energy density. Unlimited focusing is a property of the ideal geometric model, not a prediction for the resolved physical wave.

The Boussinesq approximation requires $\delta =(\rho -{\rho }_0)/{\rho }_0=\varphi ({\rho }_p/{\rho }_0-1)\ll 1$. The reference mass density ${\rho }_0$ can then be used in inertia, while the small density excess is retained in buoyancy. Dilute particle volume fraction $\varphi \ll 1$ alone does not suffice if ${\rho }_p/{\rho }_0$ is exceptionally large. Define $G=g({\rho }_p-{\rho }_0)/{\rho }_0$, so the reduced gravity is $g^{\prime }=G\varphi$.

The shallow water equations require depth small compared with the horizontal evolution scale, weak vertical acceleration and approximately hydrostatic pressure. Here the ambient is deep and quiescent, so its leading contribution is a hydrostatic reference pressure; the excess pressure in the layer is ${\rho }_{0}G\varphi (h-z)$. The model also needs nearly uniform streamwise velocity across the section and the stated absence of secondary circulation. The advancing gravity current nose is a separate region where these assumptions can fail.

A nearly uniform particle volume fraction can be maintained by turbulent mixing with eddy diffusivity $K_z$, provided $h^2/K_z\ll h/w$ and the mixing time is short compared with horizontal advection and bulk evolution. Equivalently, $wh/K_z\ll 1$ keeps the bulk settling-induced gradient small. The particles should respond rapidly enough to follow the mixing motions. A thin boundary layer adjacent to the absorbing bed or walls need not be well mixed.

Uniform bulk concentration does not imply zero sedimentation. The relative downward particle velocity still supplies a particle deposition flux of approximately $w\varphi$ to an absorbing boundary. Mixing redistributes the remaining particles and maintains the bulk concentration profile while its overall level decreases; it need not suspend every particle indefinitely. Deposition requires negligible resuspension in the model.

The channel is prismatic: its width depends on height, not on $x$. A layer of depth $h$ has cross-sectional area $A={\int }_0^h\beta zdz=\beta h^2/2$. Volume conservation gives $A_t+(Au{)}_x=0$. With $g^{\prime }=G\varphi$, the integrated excess hydrostatic pressure force isThus the streamwise momentum balance is $(Au{)}_t+(A{u}^2+G\varphi \beta h^3/6{)}_x=0$. Under the specified vertical-settling approximation, the horizontal projection of the depositional boundary has width $\beta h$, so the particle balance is $(A\varphi {)}_t+(Au\varphi {)}_x=-\beta hw\varphi$, where $w=w_{0}f(\varphi )>0$ is the downward speed magnitude. The prismatic triangular-channel shallow water equations are thereforeThe pressure coefficient $h/3$ is the section-weighted mean of $h-z$. The factor two in deposition comes from top width divided by area. There is no streamwise widening term, because $\beta$ is constant along the channel.

Define the positive characteristic speed scale $c=\sqrt{G\varphi h/2}$. In variables $(u,\varphi ,c)$ the equations becomeThe coefficient matrix of this hyperbolic system has eigenvalues $u,u+c,u-c$. Hence the three characteristic families areAlong $dx/dt=u$, the concentration equation is the ordinary differential equation $d\varphi /dt=-2w\varphi /h$. Along $dx/dt=u\pm c$, the left eigenvectors $(1,\mp 4c/(3\varphi ),\pm 4)$ give the sedimenting triangular-channel characteristic compatibility equationsHere every derivative in a given equation follows that characteristic family, and $h=2c^2/(G\varphi )$. These three compatibility equations form the characteristic description; they are not three independent closed equations for all fields on any one curve. In particular, $u\pm 4c$ are not conserved Riemann invariants when the concentration varies: its differential and the deposition source must both be retained. The description assumes $h,\varphi >0$; the dry or zero-buoyancy limit is degenerate.

At the gravity current nose, the depth and velocity change over a distance comparable with the depth, with appreciable vertical acceleration and mixing. The hydrostatic approximation and the depth-uniform interior model do not resolve that structure. A suitable gravity-current front condition iswith a positive, order-one Froude number determined by the nose closure and geometry. It is not fixed by the interior equations alone. This definition uses $\sqrt{g^{\prime }h}$; if the front condition is expressed using $c=\sqrt{g^{\prime }h/2}$, its numerical coefficient changes by $\sqrt{2}$.

For a late-time gravity-current box model, take the depth and particle volume fraction to be uniform over $0<x<X(t)$, and neglect entrainment and resuspension. The initial volume is $\mathcal{V}=\beta h_0^{2}L/2$, soThe total particle volume is $\mathcal{V}\varphi$. Deposition over the projected boundary width gives $d(\mathcal{V}\varphi )/dt=-\beta hXw\varphi$. With the specified hindered settling factor, the prismatic triangular-channel gravity-current box model isInitializing this approximate model with $X(0)=L$ and $\varphi (0)={\varphi }_0$ matches the original volume and concentration; it does not assert that the immediate dam-break flow is a uniform box.

Eliminate time to obtainSince the left primitive is $2\mathrm{artanh}\sqrt{\varphi }$, the hindered-settling runout invariant in a prismatic triangular channel isThis relation, the volume constraint and $\dot{X}$ determine the late-time evolution. If desired, the time is given by the quadrature $t={\int }_L^{X(t)}[\mathrm{Fr}\sqrt{G\varphi (s)h_0\sqrt{L/s}}{]}^{-1}ds$. The physical branch has nonnegative right-hand side in the boxed relation; it ends when that expression reaches zero.

Thus the limiting runout length of a gravity current, measured from the closed end, isThe advance beyond the dam is $X_{\infty }-L$. In the dilute limit, $\mathrm{artanh}\sqrt{{\varphi }_0}\simeq \sqrt{{\varphi }_0}$. The ideal model approaches its finite runout length of a gravity current only as $t\to \infty$: near the limit $\varphi$ decays exponentially with rate $2w_0/h_{\infty }$ and the front speed tends to zero. The exponent $7/4$ uses the prismatic triangular geometry; a channel widening with $x$ has a different volume constraint and runout law.

A square circle packing of equal radii has one circle per cell of side $2r_1$. Its packing fraction is thereforeFor the alternating-radius square circle packing shown in the source figure, let $d$ be the nearest unlike-centre separation. A two-circle repeat cell has area $2d^2$ and solid area $\pi (r_1^2+r_2^2)$. Unlike circles require $d\ge r_1+r_2$; like circles on the two sublattices require $\sqrt{2}d\ge 2r_1$ and $\sqrt{2}d\ge 2r_2$. ThusAll three nonoverlap constraints matter when varying the radius ratio.

Take $r_1\ge r_2$ as in the illustrated arrangement and set $s=r_2/r_1$. If $0\le s\le \sqrt{2}-1$, the large circles determine the spacing and $\eta =(\pi /4)(1+s^2)$ increases with $s$. If $\sqrt{2}-1\le s\le 1$, unlike contacts determine it and $\eta =(\pi /2)(1+s^2)/(1+s{)}^2$, whose derivative is proportional to $s-1$ and is nonpositive. The maximum is therefore at the joining point:The large circles then touch their large neighbours, while each small circle exactly fills a square interstice. Interchanging the species gives the reciprocal radius ratio. This is the highest packing fraction within the specified alternating square arrangement; it is not a claim about every possible unequal-circle packing.

Let the particle and carrier-fluid mass densities be ${\rho }_p$ and ${\rho }_f$, with identical density contrast for both species. At small particle Reynolds number, balance Stokes drag law against submerged weight:ConsequentlyAt large particle Reynolds number, use approximately constant drag coefficient $C_D$ and quadratic drag. Balancing $\frac{1}{2}{C}_D{\rho }_f\pi r^2{W}_s^2$ against the same submerged weight gives $W_s^2=8({\rho }_p-{\rho }_f)gr/(3C_D{\rho }_f)$, soBoth particles must lie in the constant-$C_D$ regime for this simple square-root ratio. More generally it contains an additional factor $\sqrt{C_D(r_2)/C_D(r_1)}$. The speeds are positive downward magnitudes; the signs of the shock velocities below depend on the chosen spatial coordinate.

Use downward distance $y=H-z$, with $z$ the height above the bottom, and write $W=W_s(a)>0$. This convention will also reproduce the printed negative deposit-front velocity in part (d). The ideal simple-cubic sphere packing has packing fraction $\pi /6$; leaving ${\varphi }_{\max }$ symbolic keeps the jump calculation independent of its value. In suspension, the kinematic sedimentation equation isIntegrate this conservation law across a moving discontinuity to obtain the Rankine-Hugoniot conditionThe brackets denote values on the increasing-$y$ side minus those on the other side. A sedimentation shock is a concentration jump moving at this secant slope of the particle flux. Deposited material is stationary and has zero flux.

For the upper clearing front, the states are $0$ and ${\varphi }_{\max }/8$, giving $U_1=7W/8$ downward. For the lower deposition front, the states are ${\varphi }_{\max }/8$ and ${\varphi }_{\max }$, giving $U_2=-W/8$, upward. ThusIn terms of height above the bottom, the shock paths areThey meet when $H=Wt$. Therefore complete settling occurs atThe final height also follows directly from particle volume fraction conservation, ${\varphi }_{\max }{h}_c=({\varphi }_{\max }/8)H$. The stationary deposit after the meeting carries no particle flux.

The diagram plots time vertically against height, as requested. Its shock slopes have the opposite spatial sign to $U_1,U_2$, because those velocities were defined in the downward coordinate.

Continue with downward coordinate $y=H-z$ and $W=W_s(a)>0$. By Stokes drag law, the small-particle isolated settling velocity is $W/4$. Normalize concentrations by ${\varphi }_{\max }$: write $q_i={\varphi }_i/{\varphi }_{\max }$. In suspension the specified independent hindered settling laws give fluxes $j_1/{\varphi }_{\max }=W{q}_1(1-q_1)$ and $j_2/{\varphi }_{\max }=(W/4)q_2(1-q_2)$. The initial state is $(q_1,q_2)=(1/8,1/8)$.

Two upper sedimentation shocks clear the species separately. The faster front removes large particles, leaving $(0,1/8)$ behind it; the small concentration is unchanged across this front because the species fluxes are independent. The slower front separates that small-only suspension from clear fluid. Their downward velocities areThere is also a growing mixed deposit at the bottom. Let its normalized solid concentrations be $(p_1,p_2)$, with $p_1+p_2=1$, and let its upward speed magnitude be $s=-{\hat{U}}_3>0$. Applying the Rankine-Hugoniot condition to each species givesAdding these relations provesThe deposited flux is zero even though each species occupies only part of the packing: the settling law applies to the suspension, while the deposit is a separate stationary state. These jump relations determine the deposit composition from sedimentation jump conditions.

The large-particle clearing front has height $H-7Wt/8$ and meets the mixed-deposit front at height $35Wt/192$. ThusAt this merger all large particles have entered the deposit. Only the small particles remain suspended, at normalized concentration $1/8$. The new material deposited above the mixed layer is pure small-particle solid, at total concentration ${\varphi }_{\max }$. Its sedimentation shock therefore has downward speedFor $t>t_1$, the deposit top and the remaining clearing front have heightsTheir meeting, or the total solid-volume balance, givesThis is two-stage bidisperse batch sedimentation: the fast species finishes first, followed by deposition of the remaining slow species.

The final solid composition by volume, and also by mass because the particle densities are identical, isThe upper pure-small layer has thickness $9H/116$. As checks, $(29/40)h_1=H/8$ accounts for all large-particle solid volume, and $(11/40)h_1+(h_2-h_1)=H/8$ accounts for all small-particle solid volume, in units of ${\varphi }_{\max }$. If “percentage of particles” means number rather than volume, a large sphere occupies eight times the volume of a small sphere: the mixed bottom layer contains $29/117\simeq 24.8\%$ large particles by number and $88/117\simeq 75.2\%$ small. The upper layer remains entirely small under either convention.

At a fixed height, first determine the stored quantities per unit vertical height and per unit length along the line source. With $\Delta \rho =\rho -{\rho }_0$, the triangular-profile line plume has the useful profile integralThe geometric volume is $\mathcal{A}=2b$; the mass is $\mathcal{D}={\rho }_0{I}_0+\Delta \rho I_1=b(\rho +{\rho }_0)$. The upward momentum is $\mathcal{P}=W({\rho }_0{I}_1+\Delta \rho I_2)$, while the buoyancy force is $\mathcal{B}=-g\Delta \rho I_1$. ThusHere positive buoyancy corresponds to a rising light plume. If instead buoyancy is defined kinematically, divide $\mathcal{B}$ by ${\rho }_0$; the density-weighted convention here is the one matching the printed equations.

Multiply each local density by the vertical velocity to obtain the fluxes through a horizontal section, still per unit source length:In particular, the stored momentum equals the mass flux $Q$, and the buoyancy flux is $F=g({\rho }_{0}V-Q)$. The different triangular-profile factors in storage and flux must not be replaced by top-hat coefficients.

The Batchelor entrainment hypothesis models turbulent ingestion of ambient fluid through the plume edges. A self-similar turbulent plume has a spreading rate proportional to its axial velocity scale; the hypothesis represents the unresolved mixing by an inward edge speed proportional to that scale. For the prescribed centreline velocity, the Boussinesq approximation givesThe entrainment coefficient $\alpha$ is dimensionless and depends on the chosen velocity/profile convention and mixing regime. Here $e$ is the inward entrainment speed relative to the moving plume boundary, and the factor two counts its two sides.

A suitable density-corrected extension is $e=\alpha W\sqrt{\rho /{\rho }_0}$, representing the reduction in ambient ingestion for a much lighter non-Boussinesq plume. It reduces to the boxed closure when $|\rho -{\rho }_0|\ll {\rho }_0$. This is a phenomenological model, not an exact identity for arbitrary turbulent flow. Woods's non-Boussinesq plume model uses this density correction. The balances in part (c) can be written with either closure before taking their common Boussinesq limit.

Apply volume conservation, mass conservation and vertical momentum flux balance to a slice of the triangular-profile line plume. Ambient ingestion contributes $2e$ of volume and $2{\rho }_{0}e$ of mass per unit height and source length. The ambient is quiescent, so it supplies no leading vertical momentum. The integrated driving force is the stored buoyancy $\mathcal{B}$. Using the quantities from part (a), the non-Boussinesq triangular-profile line-plume balances areTogether with Batchelor entrainment, these are three equations for $W,\rho ,b$. They neglect viscous boundary stresses and streamwise pressure-force corrections within the integral-plume approximation.

Multiplying volume conservation by $g{\rho }_0$ and subtracting $g$ times mass conservation cancels the ambient sources exactly, givingThis is buoyancy conservation in the homogeneous ambient; it is a consequence of the first two balances, not an additional independent equation.

Now take the Boussinesq approximation in inertia and entrainment, retaining the small density deficit in buoyancy. Then $Q={\rho }_{0}bW$, $M=(2/3){\rho }_{0}b{W}^2$ and $F=(2/3)gbW({\rho }_0-\rho )$. ConsequentlyThe mass source is $2\alpha {\rho }_{0}W=3\alpha {\rho }_{0}M/Q$. Substitution into the mass conservation, momentum and buoyancy balances gives the unsteady Boussinesq triangular-profile line-plume balancesThese use density-weighted mass and momentum fluxes. Removing ${\rho }_0$ from some flux definitions while retaining it in the source would mix incompatible conventions.

Seek a nonzero-buoyancy, separable similarity solution of the unsteady Boussinesq triangular-profile line-plume balances, using powers of distance and time. Let $\tau =t+t_{\ast }>0$, allowing a time origin, and initially take distance $z$ from the virtual origin. Balancing the powers in the momentum and buoyancy equations and the entrainment balance gives the formOne way to recover the exponents is to let $Q\propto z^a{\tau }^p$ and $M\propto z^b{\tau }^q$. Matching the spatial powers in ${\partial }_t(Q^2/M)$, $Q_z$ and $M/Q$ gives $a=2,b=3$; matching the time powers gives $p=-1,q=-2$. The momentum equation then gives the powers of $F$. In the resulting solution, $Q^2/M$ is time-independent, so the storage term in the first equation vanishes exactly.

Direct substitution yields the coefficient equationsFor a rising buoyant solution $C>0$, the last equation gives $B=2A/3$. The first then gives $A=\alpha {\rho }_0$, and the second gives $C=B$. The separable decaying line-plume similarity is thereforeReconstructing the physical fields confirmsFor example, ${\partial }_t(FQ/M)=-2\alpha {\rho }_0{z}^2/{\tau }^3$ and $F_z=2\alpha {\rho }_0{z}^2/{\tau }^3$ cancel; the momentum balance similarly gives $-\alpha {\rho }_0{z}^2/{\tau }^2+2\alpha {\rho }_0{z}^2/{\tau }^2=FQ/M$.

The unshifted form uses $t_{\ast }=0$ and is defined for $t>0$, with a singular zero-time limit. To give a finite solution for all $t\ge 0$, choose $t_{\ast }>0$. The equations also permit replacing $z$ by $z+z_{\ast }>0$, a plume virtual origin; at a finite physical source distance this family has flux histories proportional to ${\tau }^{-1},{\tau }^{-2},{\tau }^{-3}$. The Boussinesq approximation requires $(z+z_{\ast })/(g{\tau }^2)\ll 1$ in the region modelled. This is a decaying similarity family; its time and virtual origins require source or initial data. The PDF supplies neither, so a unique forced-startup history cannot be selected from it. In particular, the singular unshifted field should not be presented as a regular solution at $t=0$.