Past exam of the mathematics course of the University of Cambridge / 2014 / iii / Paper 73 (Slow Viscous Flow)

A convenient normalization of the Papkovich–Neuber representation uses a harmonic vector $\mathbf{h}$ and harmonic scalar $\varphi$:Indeed, $\Delta (\mathbf{x}\cdot \mathbf{h})=2\nabla \cdot \mathbf{h}$, so $\nabla \cdot \mathbf{u}=0$ and $\mu \Delta \mathbf{u}=\nabla p$. This is a rescaling of the usual harmonic-potential representation of Stokes flow.

Measure $\mathbf{x}$ from the sphere center, set $r=|\mathbf{x}|$ and $\mathbf{n}=\mathbf{x}/r$. A translational vector harmonic $\mathbf{U}/r$ provides the decaying $1/r$ force field, while the scalar dipole $(\mathbf{U}\cdot \mathbf{x})/r^3$ adjusts the surface velocity without changing that leading far field. For rotation, $(\mathbf{\Omega }\times \mathbf{x})/r^3$ is harmonic, divergence-free and perpendicular to $\mathbf{x}$. Thus tryThe translational velocity is $\frac{A}{2r}(I+\mathbf{nn})\mathbf{U}-\frac{C}{2r^3}(I-3\mathbf{nn})\mathbf{U}$. Matching its independent tangential and radial components at $r=a$ gives $A=3a/2$, $C=-a^3/2$; matching rotation gives $B=a^3$. ConsequentlyThis superposes the translating sphere in Stokes flow and the rotating sphere in Stokes flow. It is exactly $\mathbf{U}+\mathbf{\Omega }\times \mathbf{x}$ at the sphere and tends to zero at infinity; the additive ambient pressure is set to zero.

Write $\mathbf{n}=\mathbf{R}/R$. At fixed geometry, Linearity of Stokes flow makes the translational velocities linear in the couples. With zero forces and couples perpendicular to $\mathbf{R}$, isotropy and reflection symmetry allow only terms proportional to ${\mathbf{G}}_j\times \mathbf{n}$ in each ${\mathbf{U}}_i$. There is no longitudinal term: reflect in the plane containing $\mathbf{R}$ and the common torque axis, taking account of the axial-vector transformation of torque, and use linearity under torque reversal. HenceThe direction of $\mathbf{R}$ can still change; constant separation does not mean stationary centers.

When the couples are equal, a rotation through $\pi$ about their common axis exchanges the identical spheres and reverses the in-plane translational velocities while leaving the axial angular velocities unchanged. Uniqueness of Stokes flow then givesThe midpoint is stationary and the pair can orbit it. The body spin need not equal the angular rate of that orbit. These symmetry statements hold at arbitrary noncontact separation, not just in a distant rotlet approximation.

The minimum dissipation theorem fixes the boundary velocities and the far-field velocity: among admissible incompressible velocity fields, the Stokes flow minimizesFor a competitor $\mathbf{u}+\mathbf{v}$ with zero boundary data for $\mathbf{v}$, integration by parts eliminates the cross term and leaves $\mathcal{D}[\mathbf{u}+\mathbf{v}]-\mathcal{D}[\mathbf{u}]=2\mu \int \mathbf{e}(\mathbf{v}):\mathbf{e}(\mathbf{v})dV\ge 0$.

Extend the actual two-sphere velocity rigidly through sphere 2. It is an admissible field in the one-sphere exterior, with the same translation and rotation of sphere 1. It is continuous across the filled boundary, incompressible, and adds zero strain dissipation inside. Applying the theorem to the exact isolated-sphere flow therefore givesWith all applied couples and the second force zero, the boundary-work identity is ${\mathcal{D}}_{\mathrm{pair}}={\mathbf{F}}_1\cdot {\mathbf{U}}_1$. Thus $6\pi \mu a|{\mathbf{U}}_1{|}^2\le {\mathbf{F}}_1\cdot {\mathbf{U}}_1\le |{\mathbf{F}}_1||{\mathbf{U}}_1|$, which provesThis fixed-force comparison of minimum viscous dissipation uses a comparison at fixed actual velocity first; directly comparing different-force or different-velocity solutions would not justify the result.

If ${\mathbf{G}}_1\ne 0$, the power is instead ${\mathbf{F}}_1\cdot {\mathbf{U}}_1+{\mathbf{G}}_1\cdot {\mathbf{\Omega }}_1$, so the preceding estimate no longer bounds the force contribution alone. The two-sphere hydrodynamic mobility matrix generally has a nonzero self translation-rotation coupling: the freely moving second sphere reflects the first sphere's torque field. In the planar geometry this coupling produces a velocity perpendicular to $\mathbf{R}$; by choosing the sign and magnitude of the couple when the force has a component in that direction, its contribution to ${\mathbf{F}}_1\cdot {\mathbf{U}}_1$ can exceed the isolated force-only value. There is no universal force-only inequality with an additional applied couple. Special geometries can eliminate the coupling, but do not restore a general theorem.

Let $\zeta =6\pi \mu a$, $ϵ=a/R$ and define the dimensionless leading cross-mobility tensorFor a given applied force, the outgoing Stokeslet is fixed by that force, regardless of the sphere's entrainment velocity. At the other center the Faxén translation law gives its force-only velocity plus the incident Stokeslet velocity. The finite-size dipole and the Laplacian Faxén correction first contribute $O({ϵ}^3|\mathbf{V}|)$; an emitted stresslet reflection is still smaller at this order. ThereforeThere is no second-order translational term for prescribed equal forces.

The vorticity of translating-sphere Stokes flow is $\mathbf{\omega }=3a\mathbf{V}\times \mathbf{r}/(2r^3)$; its potential dipole has zero curl. Apply the torque-free Faxén rotation law to the other sphere's field, using $\mathbf{r}=-\mathbf{R}$ at sphere 1 and $\mathbf{r}=\mathbf{R}$ at sphere 2. This yieldsthrough second order in the dimensionless spin $a\mathbf{\Omega }/|\mathbf{V}|$. Neglected spin terms are higher order. This is the method of reflections for Stokes flow, with the ambient field at each application excluding that sphere's own disturbance.

Put ${\mathbf{W}}_2={\mathbf{F}}_2/\zeta$. To the needed accuracy the Faxén translation law gives ${\mathbf{U}}_2={\mathbf{W}}_2+\mathsf{M}\mathbf{V}+O({ϵ}^3|\mathbf{V}|)$. Enforcing its prescribed velocity therefore givesHere the force on sphere 1 remains prescribed; no extra second-order inversion term is needed in ${\mathbf{F}}_2$. Substitute this force into the incident field at sphere 1:The negative second-order term reflects the reduced force required on sphere 2 because sphere 1 already entrains it.

Sphere 2 is now force-free and torque-free. Its translation and spin follow the incident rigid-motion components, so it emits neither a Stokeslet nor a rotlet. Its first disturbance is a stresslet from the incident rate of strain, of magnitude $\mu a^3(a|\mathbf{V}|/R^2)$. Its return velocity at sphere 1 is consequently $O(a^4|\mathbf{V}|/R^4)$:More precisely, that leading correction is $-(15/4){ϵ}^4(\mathbf{V}\cdot \mathbf{n})\mathbf{n}$. It vanishes for purely transverse forcing, for which the first translational correction is of order ${ϵ}^6$.

The first-return spin requires an additional cancellation. The leading incident strain at sphere 2 isIt is axisymmetric about the line of centers, so the reflected stresslet has zero vorticity at sphere 1 on that axis. A naive vorticity estimate of order $(|\mathbf{V}|/a){ϵ}^5$ therefore has zero coefficient. The finite-size correction to the incident strain and the next force-free scattering multipoles are two powers of $ϵ$ smaller; these generically giveThis suppressed spin from a force-free distant sphere is also consistent with the reciprocal torque-to-self-translation coupling. Purely longitudinal forcing gives zero spin exactly by axial symmetry. Distinguishing the stresslet's on-axis curl from its generic off-axis scale is essential here.

Use $y$ for height above the horizontal wall. In the lubrication approximation, vertical momentum is hydrostatic and the interfacial stress balance with variable surface tension givesThe horizontal equation is $\mu u_{yy}=B$, with $u(0)=0$ and $\mu u_y(h)={\gamma }_x=-A{C}_x$. ThusIntegrating across the film and adding surface diffusion to surfactant advection gives the dimensional thin-film mass flux and insoluble surfactant flux:These include Marangoni stress, hydrostatic leveling, capillarity and surface diffusion with their signs fixed by the surface traction. In particular the derivative acts on the product of surface tension and curvature, not just on curvature.

Choose $h_{\ast }=(A{C}_0/(\rho g){)}^{1/2}$ and $u_{\ast }=A{C}_0{h}_{\ast }/(\mu L)$. The dimensionless definitions areTogether with the stated concentration and horizontal scales, these giveIn steady flow $Q,J$ are constants by liquid and surfactant conservation.

With capillarity and diffusion neglected, solve these two linear equations for the gradients, in the region $H>0,\Gamma >0$:For positive $Q,J$, the phase plane nullclines are $H=3Q\Gamma /(2J)$ for ${\Gamma }_X=0$ and $H=2Q\Gamma /J$ for $H_X=0$. Above both lines, trajectories go left and upward; between them they go left and downward; below both they go right and downward. There is no positive-quadrant equilibrium. The axes are singular boundaries of this positive-flux reduction, not regular equilibria.

The phase portrait shows these trajectories for one choice of positive fluxes; the two nullcline slopes rescale with $Q/J$.

With zero surfactant flux and positive concentration, the nondiffusive surface velocity is zero. The reduced equations become ${\Gamma }_X=6Q/H^2$, $H_X=-12Q/H^3$. Integrating, with the prescribed initial height and concentration, givesAt the far end, the imposed concentration drop fixes $H(1{)}^2=H_0^2+4\delta$, henceThis steady surfactant film with zero surface flux rises to the right while its concentration decreases. The liquid flux is negative, despite the rightward Marangoni traction: the adverse hydrostatic gradient is strong enough to immobilize the surface and drive the interior leftwards.

In dimensionless height $Y=y/h_{\ast }$, the velocity profile isIt vanishes at both boundaries and is negative in the interior when $\delta >0$. Its integral is $Q=-H^3{H}_X/12$. For $\delta =0$, the nondegenerate film is flat and stationary. If its initial thickness is zero, the formal $X^{1/4}$ edge has an unbounded slope, so the lubrication approximation applies away from that edge rather than at the exact dry point.

For zero liquid flux the reduced gradients are ${\Gamma }_X=-4J/(H\Gamma )$ and $H_X=6J/(H^2\Gamma )$. Their ratio givesIntegrate $dX/dH=H^2\Gamma /(6J)$ to obtain the requested implicit height profile:The endpoint relation is $H_1^2=H_0^2+3\delta$, and the same expression at $X=1$ fixes $J$. For $H_0=0$ it simplifies tosoThis zero-liquid-flux surfactant film grows to the right from a formal dry edge, with $H\sim (18JX{)}^{1/3}$ there. Its dimensionless velocity isThe lower two thirds flow left, the upper third flows right, and the integrated liquid flux is zero. Surface advection carries the positive surfactant flux, with $U_s=H^2{H}_X/6=J/\Gamma$.

The flux curve obeys $dJ/d\delta =(\sqrt{3}/4)\sqrt{\delta }(1-\delta )$. It increases throughout the interval and reaches its formal maximum $J_{\max }=1/(5\sqrt{3})$ at $\delta =1$. Larger depletion strengthens the Marangoni driving, so increasing transport is plausible. But finite flux at a zero-concentration endpoint is a singular prediction: $U_s=J/\Gamma$ and the height gradient diverge there. Neglected diffusion, capillarity or an endpoint region must regularize the physical limit, and can alter its maximum. The dry initial edge is likewise outside a uniform small-slope approximation.

The profiles and velocity arrows illustrate the opposing interior-flow directions and the formal endpoint maximum; neither sketch treats the singular edges as resolved lubrication regions.

Take $z$ positive upward, and write $b=\tilde{\rho }g$. For a straight steady viscous buoyant conduit, the inner and outer axial equations areRegularity gives $w_i^{\prime }(0)=0$, the wall has $w_o(R)=0$, and velocity and tangential traction are continuous at $r=a_0$: $w_i=w_o$, $\lambda \mu w_i^{\prime }=\mu w_o^{\prime }$. With no interfacial tension specified, normal traction is continuous; for this straight steady flow that sets $P_i=P_o$ at the interface. Thus the modified gradients are common in this case.

When their value is zero, integration givesConsequentlyThe large center-to-interface speed ratio explains why the inner parabolic component dominates at small viscosity ratio.

For a slowly varying radius, set $B=b-P_{i,z}$. The local inner solution has $w_i=w_a+B(a^2-r^2)/(4\lambda \mu )$, soThe shear-driven outer interface speed has scale $a^{2}B\ln (R/a)/\mu$. The pressure-driven return flow needed for zero total flux adds a velocity of order $q/R^2$ and a shear-driven contribution of order $a^{2}B/\mu$. Relative to the inner parabolic velocity these terms are $O(\lambda \ln (R/a))$ and $O((a/R{)}^2)$, which are small under the stated assumptions. Axial viscous derivatives are also smaller because the plume is slender. ThereforeFor clarity, the outer pressure-gradient estimate must include both return flux and transmitted interface shear. Their flux scales are $R^4{P}_{o,z}/\mu$, $a^2{R}^{2}B/\mu$ and $q$. Zero total flux hence impliesThis justifies the negligible outer pressure gradient. We subsequently choose its negligible pressure reference as zero, as in the reduced model.

Near the moving interface the outer radial source flow is $u_r=a{a}_t/r$: the interface's axial advection is smaller because its axial velocity is small compared with the inner core speed. Its radial strain is ${\partial }_r{u}_r=-a{a}_t/r^2$. Continuity of normal traction at $r=a$, neglecting the much smaller inner normal viscous stress, givesWith $A=a^2/a_0^2$, this is $P_i=\mu A_t/A$ after dropping $P_o$. Conservation of plume volume is $\pi (a^2{)}_t+q_z=0$.

Let the axial and velocity scales beSince $\mu U_0/(b{\ell }^2)=1$, substitution gives the conduit equationThe buoyant flux and the viscous normal-stress pressure are both retained at this scaling.

Linearizing about unit area gives $a_T+2a_Z-a_{TZZ}=0$. ThusWave crests propagate upward. The group velocity is $2(1-k^2)/(1+k^2{)}^2$: it is downward for $|k|>1$, so the phase direction is not the direction of every wave packet.

For a travelling wave $f(\zeta )$, integrating once and using its uniform far field givesDivide by $f^2$, multiply by $f^{\prime }/f$ and integrate again. A convenient potential, defined up to an additive constant, isAt the positive crest $f=\alpha f_0$, the derivative is zero. Subtracting the far-field potential yieldsThis is the solitary-wave amplitude-speed relation for the conduit equation. A nontrivial elevation wave has $\alpha >1$ and $c>2f_0$; as $\alpha \downarrow 1$, the speed tends to the long-wave speed $2f_0$. Positive area is required throughout the derivation.

Let the actual dimensional thickness be $d=ϵah$, and let $\bar{z}$ increase downward. The leading curvature of the outer free surface isThe Young–Laplace equation gives $p_{\bar{z}}=-\gamma (d_{\bar{z}}/a^2+d_{\bar{z}\bar{z}\bar{z}})$. In the lubrication approximation, the downward velocity at distance $y$ from the wall is $(\rho g-p_{\bar{z}})(dy-y^2/2)/\mu$, from no slip and zero free-surface shear. Integrating this profile givesCircumferential corrections are higher order in $d/a$, so volume conservation is $d_{\bar{t}}+q_{\bar{z}}=0$ at leading order. DefineThen the thin liquid film on a vertical cylinder satisfiesThe circumferential-curvature term $h_z$ is retained along with axial curvature.

The perturbation calculation givesLong waves with $0<|k|<1$ grow, while $|k|>1$ decay; the fastest growth is $1/12$ at $|k|=1/\sqrt{2}$. Disturbances drift downward at speed $G$. The dimensional unstable wavelengths exceed $2\pi a$, with most unstable wavelength $2\pi \sqrt{2}a$. This is the thin-film version of the Rayleigh–Plateau instability driven by circumferential curvature and opposed by short-wave axial curvature.

For a travelling wave, one integration fixed by $H\to 1$ givesor $H^{\prime \prime \prime }+H^{\prime }=3c(H-1)/H^3-G(1-H^{-3})$. In the large positive-speed core, substitute the expansion using the supplied assumptions. The left side at order $c^{2/3}$ is $H_2^{\prime \prime \prime }+H_2^{\prime }$, while the gravity term is order one and $3c/H^2$ is only order $c^{-1/3}$. Thus $H_2^{\prime \prime \prime }+H_2^{\prime }=0$. Its general solution is a constant plus a sine and cosine. The vanishing leading height and slope at both endpoints selectThe core length is $2\pi$ in the axial coordinate; its amplitude is fixed by transition matching.

At either edge set $X=(3c{)}^{1/3}(x-x_e)$, with $x_e=0$ at the trailing edge and $x_e=2\pi$ at the leading edge. Keeping the height $H$ of order one, the exact wave equation becomesThus the transition equation is $H^3{H}_{XXX}=H-1$. Its linearization about $H=1$ has exponents $1$ and $-1/2\pm i\sqrt{3}/2$. Behind the pulse, approaching the uniform film as $X\to -\infty$ permits only the positive real exponent. The positive-growing, nontrivial branch has one amplitude that is removed by translation, giving the relevant unique $H_{-}$. In front, decay as $X\to +\infty$ permits the two real oscillatory amplitudes; translation removes one parameter and leaves a one-parameter family $H_{+}$. These are the unstable manifold and stable manifold dimensions for the transition system. The literal boundary condition also admits the constant solution and an opposite departing branch, so uniqueness here is specifically for the branch matching a raised positive pulse.

Match the supplied trailing quadratic asymptotic to the core near zero. It gives $K_2=\kappa 3^{2/3}$, because $\frac{1}{2}{c}^{2/3}{K}_2{x}^2=\frac{1}{2}\kappa (3c{)}^{2/3}{x}^2$. The leading maximum occurs at $x=\pi$, soThe front tail has the form $H-1\sim C{e}^{-X/2}\cos (\sqrt{3}X/2+\theta )$. Its capillary waves occur ahead of the descending pulse, not behind it; the trailing tail is monotone to leading order. Converting its frequency to physical distance gives

This large solitary pulse on a cylindrical film sketch shows the leading core and the front-tail eigenmodes; the edge transitions are schematic and the ringing is enlarged for visibility. The thin-film asymptotics require both $c\gg 1$ and $ϵc^{2/3}\ll 1$, rather than taking speed to infinity at fixed film slenderness.

At the next core order, $H_0^{\prime \prime \prime }+H_0^{\prime }=-G$. Integrating and imposing zero $H_0^{\prime }$ at both endpoints givesThe constants in the two quadratic transition asymptotics match the core endpoint heights: $A_0+B_0={\beta }_{-}$ and $A_0+B_0-2\pi G={\beta }_{+}$. Hence $2\pi G={\beta }_{-}-{\beta }_{+}$, and the dimensional uniform thickness selected by matching isWith the supplied constants this fixes $G=3.75/(2\pi )$; $B_0$ is not separately fixed by these height conditions. The thickness selection and capillary wavelength follow from matched asymptotic expansion, not from the linear instability calculation.

The raised-wave core expansion separates its leading height scale from the order-one gravitational correction.

The two short edge regions supply the curvature and height matching needed by the core approximation.

Order-one gravitational forcing is retained in the next core equation and selects the uniform film thickness.