Past exam of the mathematics course of the University of Cambridge / 2013 / iii / Paper 69 (Fluid Dynamics of Energy Systems)

Use a constant porosity $\varphi$ and let $v$ denote pore velocity. Darcy law gives the Darcy velocity $q=\varphi v$. If the velocity convention has already absorbed porosity, set $\varphi =1$ in the following formulas. Assume positive permeability of a porous medium throughout the layer, so $k_0>0$ and $k_0+k_1>0$.Here $\mu$ is dynamic viscosity. The flux-weighted residence time in a porous layer is the pore volume divided by throughput. Equivalently, weighting each streamline's transit time $L/v(y)$ by its inlet flux $v(y)dy$ givesThis is also the advective travel scale after cross-layer mixing has homogenized the concentration. For resident-weighted transit time in a layered flow, sampling the initial fluid uniformly by resident volume, while suppressing cross-stream exchange, is a different experiment: its mean iswith the continuous $k_1=0$ limit $\mu \varphi L^2/(\Delta P{k}_0)$. Specifying the sampling convention avoids confusing the reciprocal of the mean speed with the mean reciprocal speed.

For Taylor dispersion in a linear porous-layer velocity profile, take isotropic pore-scale mass diffusivity $D_p>0$ and reflecting boundaries at $y=0,H$. Write $v-U=\gamma (y-H/2)$, where $\gamma =\Delta P{k}_1/(\mu \varphi LH)$. Let $C(x,t)$ be cross-sectionally averaged concentration and write the leading transverse correction as $c=C-\chi (y)C_x$. Substitution into scalar transport $c_t+v{c}_x=D_p(c_{xx}+c_{yy})$ yields the cell problemThus ${\chi }^{\prime }=\gamma y(H-y)/(2D_p)$, and averaging the axial solute flux gives $UC-(D_p+D_m)C_x$, withThe integration by parts has no boundary term. It also proves nonnegativity, even when $k_1<0$. The total effective mass diffusivity is $\mathcal{D}=D_p+D_m$. For anisotropic pore dispersion, the denominator uses the transverse coefficient and the added molecular term uses the axial coefficient.

There are two separate tests for significance. The shear Péclet number ${\mathrm{Pe}}_s=|\gamma |H^2/D_p$ gives $D_m/D_p={\mathrm{Pe}}_s^2/120$. A large value implies substantial enhancement of axial spreading. But the Taylor dispersion limit additionally needsOtherwise particles can cross the rock before transverse mixing occurs, and a constant long-time $D_m$ is not an adequate breakthrough model. The relative front width is of order $\sqrt{\mathcal{D}/(UL)}$ once the long-time model applies.

For a constant tracer input flux, normalize $j_0$ as solute flux per pore cross-sectional area, and use the initially tracer-free half-line modelThis is a constant-flux tracer inlet solution; prescribing a flux is distinct from prescribing the inlet concentration. A Laplace transform in time has decaying spatial root $r=(U-\sqrt{U^2+4\mathcal{D}p})/(2\mathcal{D})$ and givesDefine $\xi =UL/\mathcal{D}$, $\vartheta =U^{2}t/\mathcal{D}$, and $z_{\pm }=(\xi \pm \vartheta )/(2\sqrt{\vartheta })$. Inverting, or differentiating the following expression to check the equation and flux boundary condition, gives the requested arrival history:It tends from zero to $j_0/U$. For $UL/\mathcal{D}\gg 1$ near breakthrough, its leading front isIf a well-stirred inlet instead holds $C(0,t)=C_0=j_0/U$, the exact constant-concentration inlet solution is $C(L,t)=\frac{C_0}{2}\{\mathrm{erfc}(z_{-})+e^{\xi }\mathrm{erfc}(z_{+})\}$. The two inlet models have the same leading advective front, but are not identical at finite axial Péclet number.

The required boundary model is a fracture that drains after the storm, with $h(0,t)=0$, and a current of finite extent whose flux vanishes at its advancing nose. Introduce mass density $\rho$, gravity $g$, constant porosity $\varphi$ and dynamic viscosity $\mu$. Hydrostatic Darcy law gives the depth-integrated discharge $q=-k\rho gh{h}_x/\mu$. Fluid-volume conservation then gives the Boussinesq equation for an unconfined aquiferThe fluid volume per fracture length is $V=\varphi {\int }_0^{\infty }hdx$; omitting the constant $\varphi$ does not change its fractional decay rate.

The conserved first moment of a draining porous current follows directly:The finite outlet discharge does not contribute to $[xh{h}_x]$ at $x=0$. This establishes $D=\int xhdx$ constant. An unspecified nonzero fracture head would instead give $\dot{D}=\kappa h(0,t{)}^2/2$; the drained boundary is essential.

For the dipole similarity solution of a draining porous current, set $h=a{t}^{\alpha }f(\eta )$ and $\eta =x/(d{t}^{\beta })$. The conserved first moment requires $\alpha +2\beta =0$, while matching the time powers in the evolution equation gives $\alpha -1=2\alpha -2\beta$. HenceChoose the nose at $\eta =1$ and the normalization $\kappa a/d^2=1$. The profile equation becomesFor $f=b{\eta }^2+c{\eta }^{1/2}$, the two sides have coefficients $-b,-5c/8$ and $6b^2,15bc/4$, respectively. The nonzero positive profile therefore has $b=-1/6$; its zero at $\eta =1$ gives $c=1/6$. Its normalization is determined byThus a complete choice of the constants isSet $h=0$ beyond the nose. The square-root behavior at the outlet allows finite drainage: $f{f}^{\prime }\to 1/72$ as $\eta \downarrow 0$, while $f{f}^{\prime }\to 0$ at the nose.

Since ${\int }_0^{1}fd\eta =1/18$, the volume isThe drainage volume decay in a dipole porous current shows that the printed positive rate cannot describe drainage and its factor is also inconsistent with the conserved-moment scaling. This solution supplies a direct counterexample to that rate and the corrected one. If the initial volume $V_0$ is prescribed at the start of a similarity phase, replace $t$ by $t+t_0$ with $t_0=(\varphi ad/(18V_0){)}^4$. A general post-storm initial shape need not be exactly self-similar; this is the dipole similarity profile and its long-time scaling, not an asserted exact profile for every initial condition. The initial volume alone does not specify $D$.

Take $x$ increasing upslope and let $u$ be the superficial Darcy velocity, so the mobile discharge is $uh$. Its pore transport speed is $u/\varphi$. This interpretation is required by the factors of porosity in the printed equations; a literal pore-speed convention would instead replace $u$ there by $\varphi u$.

For capillary residual trapping, let $M(x,t)={\max }_{0\le \tau \le t}h(x,\tau )$ be the maximum invaded thickness, including the initial state. The stored carbon-dioxide volume per plan area isDuring advance, $M=h$ and new pores fill with mobile carbon dioxide; during recession, $M$ is fixed and the newly vacated pores retain saturation $s$. Conservation of mobile plus trapped fluid isConsequentlyAssume $a,L,u,\varphi >0$ and $0<s<1$. The two speeds differ because the advancing front fills a whole pore volume while the receding tail removes only its mobile fraction.

The method of characteristics keeps height constant on $x=x_0+v_{A}t$ on the leading face and $x=x_0+v_{R}t$ on the trailing face. Matching the two linear profiles gives the triangular current with capillary retention:Before extinction, the rear, front, crest position and crest height areEvery positive height has one trailing and one leading characteristic; their intersections trace the crest. The mobile current disappears at

The final trapped region is the maximum-invasion envelope of a retained current. For $0<x<L$, the initial height is already the maximum, $M_{\infty }=ax$. For $L<x<x_{\ast }$, the maximum occurs when the crest passes $x$, at $t_m=2(x-L)/(v_R+v_A)$. Substituting into the trailing profile givesIf $z$ is distance into the aquifer measured normally from its upper boundary, trapped carbon dioxide occupies $0<z<M_{\infty }(x)$ at residual pore saturation $s$. The occupied pore volume per unit transverse width isexactly the initial mobile volume. This is an independent mass-conservation check of both the envelope and extinction distance. In the limit $s=0$, the current translates without shrinking and nothing is trapped; the finite-extinction formula is not used in that limit.

Write $G=g^{\prime }$ for the positive reduced gravity of dense fluid. With $z$ upwards, $G_z>0$ places denser fluid above lighter fluid and permits convective turbulence. A displacement over the radius scale $d$ samples a buoyancy difference of order $d{G}_z$; balancing kinetic energy per mass with buoyant work over $d$ givesThe buoyancy-gradient mixing-length closure then gives an eddy diffusivity $K\sim ud=d^2\sqrt{G_z}$. We use the paper's unit-prefactor convention for this dimensional closure and its effective cross-sectional area $d^2$; order-one mixing constants and the cylinder's geometric factor $\pi$ are suppressed in its displayed evolution equation.

The downward reduced-gravity transport magnitude and its convergence are different quantities:The printed expression labelled a flux is actually the second quantity, a flux convergence per unit height. The dimensions distinguish them: $\mathcal{F}$ has units $L^4{T}^{-3}$, whereas $F_{\mathrm{conv}}$ has units $L^3{T}^{-3}$. A constant positive gradient has nonzero transport but zero convergence, providing a direct counterexample to identifying the derivative with the transport itself.

The upward total buoyancy flux is $QG-\mathcal{F}$. Conservation therefore yieldsThis flux convergence in turbulent buoyancy mixing obtains the displayed evolution equation with the correct transport interpretation.

For the arrested cubic buoyancy profile, an undisturbed region below $-H$ has $G=0$ and no turbulent buoyancy flux. The joining condition is $G(-H)=0$ andIn a steady state, integrating the conservation equation from that front gives $QG=d^4(G_z{)}^{3/2}$. For a nonzero mixed region,and henceThe associated gradient is $Q^2(z+H{)}^2/(9d^8)$, so it matches the required zero-gradient state continuously. We take the mixed interval to start immediately above $-H$; adding a further zero interval would merely relocate the arrest front.

The source buoyancy flux is $B_s=Q_s{g}_s^{\prime }$, with $g_s^{\prime }=g\Delta \rho /\rho$ in the Boussinesq approximation. The stated top-source convention means that this entire flux enters the downward turbulent transport there: $\mathcal{F}(0)=B_s$. Since steady transport also satisfies $\mathcal{F}(0)=QG(0)$,This is the extent in the paper's dimensional-closure normalization. It presupposes positive $Q$; with no upflow there is no finite steady arrest depth. For a mixture of the two original fluids, $G(0)\le g_s^{\prime }$ is necessary, so an extrapolation to $Q_s/Q>1$ cannot represent a physical concentration profile under this top-flux boundary model. A dilute-source regime avoids appreciable source-volume corrections.

Assuming the same scalar-mixing coefficient for a passive tracer, its eddy diffusivity in the mixed region isIt vanishes at the arrest front and increases linearly towards the source. A different turbulent scalar diffusivity ratio would multiply this result by its corresponding constant.

The geometric normalization can be restored without changing the reasoning. With actual area $S=\pi d^2$ and $K=C_m{d}^2\sqrt{G_z}$, put $R=S{C}_m{d}^2$. Then $G=Q^2(z+H{)}^3/(27R^2)$, $H=3R^{2/3}{B}_s^{1/3}/Q$, and $K=Q(z+H)/(3S)$. These formulas show which numerical prefactors are convention-dependent.

Use the well-mixed natural ventilation model for a distributed floor source. Let the building height be $H_b$, its volume be $\mathcal{V}$, its temperature excess be $\theta$, and the exterior mass density be $\rho$. Let ${\alpha }_T$ be the coefficient of thermal expansion, so reduced gravity is $g^{\prime }=g{\alpha }_T\theta$. Convert the total heat input $Q$ to the temperature-volume flux $\mathcal{H}=Q/(\rho c_p)$, where $c_p$ is specific heat capacity. If the quoted heat flux is per floor area, first multiply it by that area.

For effective opening area for pressure-driven ventilation, define the equal-opening hydraulic coefficient $K_A$ such that the signed throughflow $q$ obeys $q|q|=K_A\Pi$, with $\Pi$ the total driving pressure divided by mass density. In the effective-area convention of question 6, each opening has discharge $A\sqrt{\delta p/\rho }$, so $K_A=A^2/2$. The usual physical-area orifice law instead gives $K_A=C_d^2{A}^2$; that prefactor does not affect the stability conclusions.

Set $W=\Delta P/\rho >0$ and $G_T=g{\alpha }_T{H}_b$. Positive $q$ denotes entry at the floor and exit at the roof. hydrostatic pressure and wind combine to giveBecause the inflowing air is at exterior temperature in either direction, the heat budget isThe well-mixed ventilation temperature balance at steady state is thereforeAll roots must satisfy $\theta >0$ and their appropriate flow-direction inequality; squaring a pressure relation without imposing those conditions could add spurious branches.

The assisting case has a unique positive solution and it is stable. For opposing wind, define ${\theta }_w=W/G_T$. There is always one solution with $\theta >{\theta }_w$, a buoyancy-dominated upward flow. Below ${\theta }_w$, the heat-removal curve is $R(\theta )=\sqrt{K_A}\theta \sqrt{W-G_T\theta }$. It is zero at both endpoints, and its maximum occurs at $\theta =2{\theta }_w/3$:For $0<\mathcal{H}<R_{\max }$ there are two downward-flow roots as well as the upward-flow root. Equivalently, setting $y=G_T\theta /W$ and $\epsilon =G_T\mathcal{H}/(\sqrt{K_A}{W}^{3/2})$, the equilibrium curve isFor linear stability analysis of the opposing-wind ventilation bistability, linearize the heat budget: the eigenvalue is $-R^{\prime }(\theta )/\mathcal{V}$. The cooler downward root has $0<\theta <2{\theta }_w/3$ and $R^{\prime }>0$, so it is stable. The warmer downward root has $2{\theta }_w/3<\theta <{\theta }_w$ and $R^{\prime }<0$, so it is unstable and separates the two basins. The upward root has $R^{\prime }>0$ and is stable. At $\mathcal{H}=R_{\max }$ the two downward roots merge at a saddle-node bifurcation; above that value only the upward equilibrium remains.

For ventilation switching under wind and heating changes, increasing wind pressure or decreasing heat input both reduce $\epsilon$, so both can create the pair of wind-driven equilibria. They nevertheless have different physical and transient effects. On the hot upward branch, implicit differentiation of $\mathcal{H}=\sqrt{K_A}\theta \sqrt{G_T\theta -W}$ shows that increasing $W$ increases $\theta$, while decreasing $\mathcal{H}$ decreases $\theta$. Both reduce the upward ventilation rate, but only a wind change shifts ${\theta }_w$.

Under slow parameter variation, a state follows its current stable branch while that branch persists; the upward branch does not disappear when wind increases or heating decreases. A sufficiently large abrupt wind increase can raise ${\theta }_w$ beyond the current temperature and place that state below the new unstable threshold, causing flow reversal and cooling towards the wind-dominated equilibrium. A jump insufficient to cross its basin boundary returns to the upward equilibrium.

At fixed $W$, reducing a still-positive heat input cannot force an initially upward-flow state through ${\theta }_w$: at that boundary the temperature balance has $\dot{\theta }=\mathcal{H}/\mathcal{V}>0$. Thus heat reduction alone does not cause that reversal in this model. Conversely, decreasing wind or increasing heat from the stable downward branch eventually destroys it at the fold and forces a switch to the upward branch, giving history dependence. Wind strengthening and heat reduction have the same effect on the steady dimensionless control parameter, but not the same temperature change or switching dynamics. These stability statements use quasi-steady opening flow and a single well-mixed temperature; other stratification or airflow-inertia models require their own stability analysis.

Let $q_b$ be floor inflow and $q_t$ roof outflow. Use effective-area discharge $A\sqrt{\delta p/\rho }$, matching the factors in the supplied relations. With negligible source volume, both equal $q$. The lower region has exterior mass density while the warm upper layer has reduced gravity $g^{\prime }$ and depth $H-h$. hydrostatic pressure therefore supplies a total opening pressure head $g^{\prime }(H-h)$. Equal openings sharing equal discharge take half each, givingSteady lower-layer volume conservation sets this ventilation rate equal to the entrained turbulent plume flux at the interface, $q=\lambda B^{1/3}{h}^{5/3}$. Upper-layer buoyancy conservation gives $B=q{g}^{\prime }$. These are the required balances for displacement ventilation:Eliminating $g^{\prime }$ and $q$ gives the displacement-ventilation interface height equationThe left side is strictly increasing from zero to infinity, so there is one interface height. It depends on opening geometry and plume entrainment, not on $B$; the resulting $g^{\prime }=B^{2/3}/(\lambda h^{5/3})$ and throughflow do depend on $B$.

For finite source volume $Q$, conservation instead gives $q_t=q_b+Q$. The two pressure drops are no longer equal:Let $q_p(h;Q)$ be a physically entraining source plume: $q_p(0;Q)=Q$ and $q_p(h;Q)>Q$ for $h>0$. Its interface balance is $q_p=q_t=Q+q_b$. Therefore, as $q_b$ reaches zero, a positive-depth lower layer has no replenishment and cannot remain steady: this is no steady displacement layer without ambient supply. The limiting interface is $h=0$. At the threshold $q_t=Q$ and $g^{\prime }=B/Q$, and the entire hydrostatic head drives the roof discharge. HenceThis is source-volume blocking of displacement ventilation. The half-head factor of the negligible-source state must not be kept after the floor inflow vanishes. In the usual physical-area convention $q=A_{\mathrm{phys}}{C}_d\sqrt{2\delta p/\rho }$, the equivalent formula is $Q_c=(2C_d^2{A}_{\mathrm{phys}}^{2}BH{)}^{1/3}$.

The quoted pure-plume law cannot be used unmodified down to the origin for a source with nonzero volume: it would incorrectly give zero flux there. A compatible finite-source example uses a virtual origin,It recovers the supplied pure-plume limit when $Q\to 0$, and gives the same blocking threshold. The complete finite-$Q$ interface trajectory depends on that plume model, but the threshold needs only source-volume conservation and positive entrainment. For $Q>Q_c$, the assumed steady regime with a lower inflow no longer exists.