Past exam of the mathematics course of the University of Cambridge / 2006 / iii / Paper 75 (Solidification of Fluids)

Take $z$ upward from the unperturbed interface in the apparatus frame, with liquid at $z>0$ and pulling velocity $-V{\mathbf{e}}_z$. The primary ice contains no salt. We use the dilute, equal-density approximation and neglect fluid motion relative to the mould. The initial salt inventory per unit horizontal area, with a common density factor suppressed, is $N=h{C}_0$. Its conservation is essential: the far liquid has zero concentration, so an imposed nonzero far-field composition would describe a different experiment.

The steady advection-diffusion equation is $-V{C}^{\prime }=D{C}^{\prime \prime }$. Decay into the pure water and interfacial salt rejection giveSince ${\int }_0^{\infty }C(z)dz=D{C}_I/V=N$, the finite-inventory directional salt rejection solution isThe decay length is $D/V$. Local liquidus equilibrium givesThus the thermal apparatus selects the interface position at this isotherm. The stated speed inequality is precisely $C_I<C_E$, so the interface is warmer than the eutectic temperature, $T_I>T_E=-m{C}_E$. These are steady values after the initial finite saline layer has rearranged; the final exponential distribution is not confined to the original depth $h$.

For $Vh{C}_0/D>C_E$, a steady pure-ice front would require an interfacial composition beyond the eutectic composition. That branch is unavailable. During the transient, rejected brine concentrates until it reaches $C_E$; secondary salt-bearing crystals can then form and the residual brine freezes as an eutectic system. Salt is incorporated into the product rather than remaining indefinitely in a super-eutectic exponential layer. The detailed pattern can include brine pockets and a nonplanar freezing region; it depends on the thermal history and the morphological stability. The speed inequality alone does not decide that latter stability.

For the stability calculation use an interface displacement $\eta =\hat{\eta }{e}^{ikx+st}$ and a liquid perturbation $c=\hat{c}{e}^{ikx+st-qz}$. Here $k>0$ is the wavenumber, $s$ is the temporal rate and $\Re q>0$ gives spatial decay. The frozen thermal field is $T_I+Gz$. Define $\Gamma ={\gamma }_{sl}{T}_m/(\rho L)>0$, where ${\gamma }_{sl}$ is solid-liquid surface energy, $T_m$ is absolute pure-solvent melting temperature, $\rho$ is mass density and $L$ is specific latent heat. For a solid protrusion into liquid the curvature is $\mathcal{K}=-{\eta }_{xx}=k^2\eta$, and the Gibbs-Thomson effect gives local equilibrium $T=-mC-\Gamma \mathcal{K}$. We assume isotropic surface energy and instantaneous interface kinetics, so no kinetic undercooling term is present.

Linearizing the equilibrium condition at the displaced interface givesThe perturbation advection-diffusion equation givesThe local solute balance at a moving graph, to first order, isIts perturbation is $(V-Dq)\hat{c}+s{C}_I\hat{\eta }=0$, because the terms $D\eta C^{\prime \prime }(0)$ and $V\eta C^{\prime }(0)$ cancel. Consequently the full frozen-temperature solutal interface dispersion isNo quasistatic approximation to solute diffusion was made here. For a growing mode the unique admissible liquid root isWriting $b=(G+\Gamma k^2)/(m{C}_I)$ and eliminating $s$ gives $(q-V/D)(q-V/D+b)=k^2$. The root continuously connected to the growing interface branch is thereforeIn particular, unstable wavelengths exist exactly when $mV{C}_I/D>G$, and their band isThis is the constitutional supercooling criterion: immediately above the plane, the local liquidus rises at rate $-m{C}^{\prime }(0)=mV{C}_I/D$, so it overtakes the imposed thermal gradient when that rate exceeds $G$. The Gibbs-Thomson effect cuts off short waves. At long wavelength $s\sim D(V/D-b)k^2/b$ for positive $G$. The exactly uniform inventory-changing displacement is not an instability of a fixed-inventory experiment. On stable branches the full diffusion problem also admits relaxation transients; the root selection above is sufficient to identify every exponentially growing interface mode.

If one additionally neglects the perturbation time derivative in the liquid diffusion equation, the alternative reduced rate isIt has the same onset and cutoff, but need not have the same rate away from the slow-growth regime. The frozen-temperature approximation by itself does not justify this extra reduction.

After onset, protrusions reach more weakly saline liquid and advance faster, while grooves accumulate rejected brine. The planar interface develops cells and, as growth becomes sufficiently strong, crystal dendrites with concentrated liquid between them. Nonlinear competition changes the selected spacing; linear rates alone cannot fix the final geometry. The interdendritic liquid forms a mushy layer, concentrates toward the eutectic composition, and ultimately freezes into a salt-bearing eutectic product when it reaches $T_E$. This replaces the pure planar salt-rejection picture with trapping and freezing of brine between crystals.

Work in the frame of the steadily advancing mush, so the rigid primary-crystal skeleton moves with $-V{\mathbf{e}}_z$. Let $\varphi$ be solid fraction, $C$ the pore-liquid concentration, $\mathbf{q}=(u,w)$ the Darcy flux relative to that skeleton, and $K$ the common thermal conductivity. Set $\kappa =K/(\rho c_p)$, the thermal diffusivity. The actual pore velocity is $\mathbf{q}/(1-\varphi )-V{\mathbf{e}}_z$. Primary crystals are taken to exclude solute; the residual pore liquid freezes into a two-solid eutectic product at the bottom. These are the standard binary mushy layer assumptions implicit in using a single liquidus and an eutectic floor.

It is convenient to use hydraulic pressure $P=p+\rho gz$, absorbing the constant-density gravitational term. With equal phase properties, no solute diffusion, negligible fluid inertia and local liquidus equilibrium, a complete mush model isThe last equation follows from conservation of energy with mixture specific enthalpy $c_{p}T-L\varphi$. Sensible heat is transported by the total volume flux $\mathbf{q}-V{\mathbf{e}}_z$; the primary solid carries the latent-enthalpy flux $VL\varphi {\mathbf{e}}_z$. In particular the sign of the latent-heat term is fixed by the downward-moving skeleton.

At the liquid roof, imposeHere $P_0=p_0+\rho gh$. The gradient condition matches heat flux into the imposed liquid field, since there is no phase-fraction jump at this roof and the phase conductivities are equal. The pore velocity there equals the relative Darcy velocity. At the floor,There is no liquid penetration into fully solid material. The unspecified solid fraction at the floor is determined by the solute balance, not set to one: remaining eutectic liquid freezes as it crosses the floor. If $C_s$ denotes the bulk concentration of the final solid product, the floor solute balance and Stefan condition areMaintaining the floor at $T_E$ permits the solid-side bath to remove this final latent heat. Equivalently its required gradient is $T_z(0^{-})=T_z(0^{+})+\rho LV[1-\varphi (0)]/K$. No additional imposed solid-side gradient is specified. Laterally we take the quadratic-pressure strain solution, without sidewalls; finite sidewalls would require additional flow conditions and generally break this horizontal invariance.

For uniform permeability of a porous medium, incompressible flow gives ${\nabla }^{2}P=0$. The roof pressure and zero normal flow at the floor are satisfied byHorizontal outflow is balanced by downward inflow. The parameter $\alpha$ is an inverse time; $E^2$ in the printed pressure condition already has dimensions of pressure divided by length squared, so no extra density factor belongs in this definition.

With $T=T(z)$, $C=C(z)$ and $\varphi =\varphi (z)$, the $x$-dependence disappears from both scalar balances. Although $u$ is nonzero, its horizontal divergence cancels $w_z$ in $\nabla \cdot (\mathbf{q}C)$. The exact strain-driven steady mush equations becomeTogether with $T(0)=T_E$, $T(h)=T_0$, $T^{\prime }(h)=G$ and $\varphi (h)=0$, they determine the scalar fields and the unknown thickness. This reduction proves that the equations admit horizontally invariant scalar solutions despite the imposed straining flow. An equivalent useful form, obtained by differentiating the solute equation, isThis exact form also shows why one must not silently discard latent heat in an arbitrary strong-flow calculation.

Define $\Delta T=T_0-T_E>0$, $R=\alpha h/V$ and $S=L/(c_p\Delta T)$. For the stated regular ordering $R\gg S\gg 1$, strain transport dominates skeleton transport and the bounded phase-change enthalpy contribution. More precisely, at fixed finite thermal strain number $\alpha h^2/\kappa$, the dimensionless energy equation isso its skeleton and integrated latent-heat corrections are small in $1/R$ and $S/R$ respectively. The leading strain-dominated mush temperature profile obeys $\kappa T^{\prime \prime }+\alpha z{T}^{\prime }=0$. Put $\delta =\sqrt{2\kappa /\alpha }$ and $H=h/\delta$. Integrating first for $T^{\prime }$ and then for $T$ givesThe error function profile follows from integrating a Gaussian, with $\mathrm{erf}x=(2/\sqrt{\pi }){\int }_0^x{e}^{-s^2}ds$. Matching the roof gradient gives the requested transcendental thickness relationFor prescribed $G>0$ and $\Delta T>0$ the right side decreases continuously from infinity to zero as $H$ increases, so it fixes one positive leading-order thickness. The resulting $h$ must still satisfy the assumed ordering and the phase-fraction constraints.

To find the product composition, integrate the exact solute balance from the roof to an arbitrary depth. With $f=1-\varphi$,At $z=0$, the exact floor composition is thereforeThe minus sign has a direct physical meaning: the horizontal outflow removes more salt than arrives at the same depth from a uniform-$C_0$ liquid, leaving a depleted solid product. Using the leading error function profile andyieldsThe full solid fraction follows by substituting this same profile into the preceding integral formula for $f(z)$. It has $\varphi (h)=0$ and $1-\varphi (0)=C_s/C_E$. A physical branch requires $0\le C_s\le C_E$ and $0\le \varphi \le 1$. In this monotone profile, $C_s\ge 0$ ensures a nonnegative liquid fraction everywhere; it is not permissible to interpret a negative formal result as a material composition.

The error-function reduction is a regular leading approximation, not a uniform statement for every possible scaling of the thermal strain number. If $\alpha h^2/\kappa$ becomes singularly large along the ordering limit, a thin basal layer can amplify the discarded terms. The printed two inequalities alone do not specify that extra scaling. In such a regime the displayed exact scalar equations and boundary conditions, rather than an unqualified zero-latent-heat profile, provide the appropriate thickness and composition calculation.

Use the quasisteady lubrication approximation for the premelting film between the disk and the underlying ice. The ice is anchored to the lower plane; new ice formed from water supplied through this upper film raises the disk. Let $U$ be its upward separation rate, $d(r)$ the premelted thickness, $\mu$ the liquid dynamic viscosity, $\rho$ a common ice/water mass density, $L$ the specific latent heat of fusion, and $T_m$ the absolute bulk melting temperature. As usual in this leading model, neglect the ice–water density difference, pressure-dependent changes in $T_m$ due to a common phase pressure, inertia and edge corrections. The imposed thermal field is taken as maintained by an external heat supply/removal system, so no additional heat-conduction bottleneck is introduced.

Write $\Delta =T_m-T_0>0$, $K_T=\rho L/T_m$ and $\Delta (r)=\Delta (1-r^2/a^2)$. The Clapeyron pressure relation for equal-density phases makes the pressure difference across the solid–liquid interfaceThis thermomolecular pressure is the repulsive stress transmitted through the film. To calculate a definite rate one also needs a film-thickness constitutive law. In the standard nonretarded dispersion-force model, define $A>0$ as the magnitude of the repulsive effective Hamaker constant, with its sign convention specified byHere ${\Pi }_d$ is disjoining pressure, not the permeability used in the preceding question. The sign of an effective Hamaker coefficient depends on the three materials; the positive parameter here denotes an interaction that stabilizes a premelted film. The problem does not specify $A$ numerically, so it must remain a defined material parameter.

Take the reservoir pressure at the rim to be $p_{\infty }$. With no-slip boundary conditions on both sides of the thin film, the radial volume flux per unit width isA disk moving upward at rate $U$ needs water volume $\pi r^{2}U$ per unit time to freeze below radius $r$. Regularity on the axis and mass conservation therefore give $2\pi r{q}_r=-\pi r^{2}U$, or $q_r=-Ur/2$. ThusDefine $B=A/(6\pi K_T)$, so $d^3=B/[\Delta (1-r^2/a^2)]$. Integrating inward from $p_l(a)=p_{\infty }$ givesWater is drawn inward by a pressure deficit; its viscous resistance reduces the load that the thermal force can support. The net normal stress on the disk is $p_l-p_{\infty }+P_T$. Let $W$ be the downward apparent load after subtracting ordinary buoyancy. For an ideal thin disk with negligible buoyant volume, $W=Mg$, where $g$ is gravitational acceleration. If a finite displaced volume is known, use its submerged weight instead; the stated mass and radius alone do not determine that volume.

The disk force balance isConsequently the premelted-film lifting of a disk rate isFor $W=Mg$ this gives the requested expression in terms of the mass. Every additional symbol in this formula has been defined above. The growing-ice lifting branch exists only when the bracket is positive.

At the threshold the required water flux and viscous pressure deficit vanish. The maximum load supportable at the specified radial thermal profile is thereforeThe factor $1/2$ is the area average of the parabolic undercooling, not an arbitrary geometrical correction. At equality the lift rate is zero; strictly positive lift requires $W<W_{\max }$. There is no radius-only, temperature-independent maximum load in this ideal model: one must specify the imposed undercooling or a physically admissible range of it.

For a fixed positive load define ${\Delta }_c=2W{T}_m/(\pi a^2\rho L)$ and $U_{\infty }=A/(6\pi \mu a^2)$. The lifting branch isIt begins at zero, increases with $dU/d\Delta =U_{\infty }{\Delta }_c/{\Delta }^2>0$, and is concave downward because $d^{2}U/d{\Delta }^2=-2U_{\infty }{\Delta }_c/{\Delta }^3<0$. Below threshold no positive steady heave is possible: the thermomolecular force is too weak to support the load. The mathematical continuation to negative $U$ is not a prediction of upward ice-driven separation there.

Increasing undercooling increases the available thermal pressure, but it also thins the premelting film. In the inverse-cube disjoining pressure model, $d^3\propto {\Delta }^{-1}$, so the hydraulic conductance decreases in exact inverse proportion to the driving scale. The remaining load correction tends to zero and the rate approaches a plateau:For positive load this is approached from below, not attained at finite undercooling. With zero load the reduced model gives the same constant rate wherever its thin-film assumptions hold.

Two limits of the calculation are worth making explicit. First, the formal $d(r)$ diverges at the warm rim, so a narrow outer connection to bulk water replaces the thin-film solution there; this leading model assumes that its contribution to the integrals is negligible. Second, very large undercooling produces a molecularly thin film and eventually invalidates continuum lubrication theory and the simple dispersion-force law. Thus the plateau is the maximum within the stated ideal constitutive model, not an unrestricted extrapolation for real water.

If a different molecular interaction is intended, the result is not universal. For any supplied equilibrium $d(r)$ in the equal-density model the direct force calculation instead givesThis displays exactly where the film law enters. The closed formula and plateau above use the conventional repulsive inverse-cube law; premelting alone, without that material law, would not determine a unique rate curve.