We found a multi-scale steady solution of the Boussinesq equations for Rayleigh–Bénard convection in a three-dimensional periodic domain between horizontal plates with a constant temperature difference. This was realised using a homotopy from the wall-to-wall optimal transport solution provided by Motoki et al. (J. Fluid Mech., vol. 851, 2018, R4). A connected steady solution, which is a consequence of bifurcation from a thermal conduction state at Rayleigh number $Ra\sim 10^{3}$, is tracked up to $Ra\sim 10^{7}$ using a Newton–Krylov iteration. The three-dimensional exact coherent thermal convection exhibits a scaling of $Nu\sim Ra^{0.31}$ (where $Nu$ is the Nusselt number) as well as multi-scale thermal plume and vortex structures, which are quite similar to those in turbulent Rayleigh–Bénard convection. The mean temperature profiles and the root-mean-square of the temperature and velocity fluctuations are in good agreement with those of the turbulent states. Furthermore, the energy spectrum follows Kolmogorov's $-5/3$ scaling law with a consistent prefactor, and the energy transfer to small scales with a nearly constant flux in the wavenumber space is in accordance with the turbulent energy transfer.

We analyse the model problem of rigid, spherical pendula impacting on rubber membranes at different Reynolds numbers to understand the contact dynamics in deformable bodies in a viscous fluid. We have investigated the problem both by laboratory and numerical experiments and a new contact model has been developed to perform the simulations. We have found that the collision dynamics depends on many parameters, the most important ones being the impact Stokes number ($St$) and the ratio of the membrane thickness ($e$) to sphere diameter ($2R$). For $e/(2R)=0.1$ the pendulum rebounds for $St \geqslant 18$ and this threshold increases as the membrane thickness decreases. Also, the membrane surface affects the dynamics of the impact by delaying the rebound for smaller membranes.

Celebrating the life and work of Professor George K. Batchelor FRS Founding Editor of the Journal of Fluid Mechanics.

Graphical Abstract © Mark Hallworth. Reproduced with Permission.

We present a numerical study of the dynamics of an elastic fibre in a shear flow at low Reynolds number, and seek to understand several aspects of the fibre's motion using the equations for slender-body theory coupled to the elastica. The numerical simulations are performed in the bead-spring framework including hydrodynamic interactions in two theoretical schemes: the generalized Rotne–Prager–Yamakawa model and a multipole expansion corrected for lubrication forces. In general, the two schemes yield similar results, including for the dominant scaling features of the shape that we identify. In particular, we focus on the evolution of an initially straight fibre oriented in the flow direction and show that the time scales of fibre bending, curling and rotation, which depend on its length and stiffness, determine the overall motion and evolution of the shapes. We document several characteristic time scales and curvatures representative of the shape that vary as power laws of the bending stiffness and fibre length. The numerical results are further supported by an interpretation using an elastica model.

Flowing granular materials exhibit several features that distinguish them from molecular fluids. A prominent feature is dilatancy, or volume deformation caused by shear deformation. Its significance in sustained flow has not been much appreciated, as its effect was thought to be confined to thin shear layers. However, it was demonstrated recently by Krishnaraj & Nott (Nat. Commun., vol. 7, 2016, 10630) that dilatancy drives a large-scale secondary flow in a cylindrical Couette device. They hypothesized that the combination of shear and gravity, when their directions are non-collinear, is necessary for the occurrence of the secondary flow. In this paper we investigate the phenomenon by considering a more complex primary flow generated in a split-bottom Couette device, wherein a part of the base of the container filled with a granular material is translated or rotated. It is known from previous studies that the height to which the material is filled determines the shape and extent of the shearing region in the primary flow. We show that the fill height also determines the form of the secondary flow, and argue that the two are intricately coupled and evolve together. Though the secondary flow is more complex than in a cylindrical Couette device, the mechanism is indeed what Krishnaraj & Nott hypothesized: the combined effect of dilatancy driven expansion, and gravity-driven flow down regions of low density. Unlike fluid instabilities that are typically driven by inertia, the secondary flow occurs at arbitrarily low shear rate and appears to be an integral part of the kinematic response.

Capsules, composed of a liquid core protected by a thin deformable membrane, offer high-potential applications in many fields of industry such as bioengineering. One of their limitations comes from the absence of models of capsule damage and/or rupture when they are subjected to an external flow. To assess when rupture is initiated, we develop a fluid–structure interaction (FSI) numerical model of a capsule in Stokes flow that accounts for potential damage of the capsule membrane. We consider the framework of continuum damage mechanics and model the membrane with an isotropic brittle damage model, in which the membrane damage state depends on the history of loading. The FSI problem is solved by coupling the finite element method, to solve for the membrane deformation, with the boundary integral method, to solve for the inner and outer fluid flows. The model is applied to an initially spherical capsule subjected to a simple shear flow. Damage initiates at a critical value of the capillary number, ratio of the fluid viscous forces to the membrane elastic forces and rupture occurs at a higher capillary number, when it reaches a threshold value. The material parameters introduced in the damage model do not influence the mode of damage but only the values of the critical and threshold capillary numbers. When the capillary number is larger than the critical value, damage develops in the two symmetric central regions containing the vorticity axis. It is indeed in these regions that the internal tensions are the highest on the membrane.

Turbulent mixing exerts a significant influence on many physical processes in the ocean. In a stably stratified Boussinesq fluid, this irreversible mixing describes the conversion of available potential energy (APE) to background potential energy (BPE). In some settings the APE framework is difficult to apply and approximate measures are used to estimate irreversible mixing. For example, numerical simulations of stratified turbulence often use triply periodic domains to increase computational efficiency. In this set-up, however, BPE is not uniquely defined and the method of Winters et al. (J. Fluid Mech., vol. 289, 1995, pp. 115–128) cannot be directly applied to calculate the APE. We propose a new technique to calculate APE in periodic domains with a mean stratification. By defining a control volume bounded by surfaces of constant buoyancy, we can construct an appropriate background buoyancy profile $b_\ast (z,t)$ and accurately quantify diapycnal mixing in such systems. This technique also permits the accurate calculation of a finite-amplitude local APE density in periodic domains. The evolution of APE is analysed in various turbulent stratified flow simulations. We show that the mean dissipation rate of buoyancy variance $\chi$ provides a good approximation to the mean diapycnal mixing rate, even in flows with significant variations in local stratification. When quantifying measures of mixing efficiency in transient flows, we find significant variation depending on whether laminar diffusion of a mean flow is included in the kinetic energy dissipation rate. We discuss how best to interpret these results in the context of quantifying diapycnal diffusivity in real oceanographic flows.

Coating of filtration membranes with a polymer brush holds great promise for efficiently preventing the deposition of fouling particles. The polymer chains are compressed by incoming particles, carried with the permeation flow towards the membrane, and consequently exert a repulsive force that acts to keep the particles away from the membrane surface. Here, we theoretically investigate the effect of a polymer brush coating on the permeation-induced hydrodynamic force, $F_{h}$, pulling a particle towards the membrane, and its balance with the steric repulsion exerted by the compressing brush, resisting the particle's approach. Lubrication theory yields an ordinary differential equation for the pressure, from which $F_{h}$ is calculated numerically. Further, an asymptotic analysis is performed for the limiting cases of a dilute or dense brush, providing analytic expressions that demonstrate how brush properties affect $F_{h}$. Finally, the equilibrium position of a particle is evaluated by considering a balance between the opposing forces. The results provide an upper boundary for the brush properties, beyond which the brush is barely compressed under conditions typical of membrane filtration processes. Further increasing the brush density or thickness only decreases the total system permeance, resulting in increased energy consumption. The results shed light on the mechanisms by which a polymer brush affects the forces acting on a foulant particle, providing quantitative measures for assessing the potential efficacy of brush coatings.

A concentrated, vertical monolayer of identical spherical squirmers, which may be bottom heavy, and which are subjected to a linear shear flow, is modelled computationally by two different methods: Stokesian dynamics, and a lubrication-theory-based method. Inertia is negligible. The aim is to compute the effective shear viscosity and, where possible, the normal stress differences as functions of the areal fraction of spheres $\phi$, the squirming parameter $\beta$ (proportional to the ratio of a squirmer's active stresslet to its swimming speed), the ratio $Sq$ of swimming speed to a typical speed of the shear flow, the bottom-heaviness parameter $G_{bh}$, the angle $\alpha$ that the shear flow makes with the horizontal and two parameters that define the repulsive force that is required computationally to prevent the squirmers from overlapping when their distance apart is less than a critical value. The Stokesian dynamics method allows the rheological quantities to be computed for values of $\phi$ up to $0.75$; the lubrication-theory method can be used for $\phi > 0.5$. For non-bottom-heavy squirmers, which are unaffected by gravity, the effective shear viscosity is found to increase more rapidly with $\phi$ than for inert spheres, whether the squirmers are pullers ($\beta > 0$) or pushers ($\beta < 0$); it also varies with $\beta$, although not by very much. However, for bottom-heavy squirmers the behaviour for pullers and pushers as $G_{bh}$ and $\alpha$ are varied is very different, since the viscosity can fall even below that of the suspending fluid for pushers at high $G_{bh}$. The normal stress differences, which are small for inert spheres, can become very large for bottom-heavy squirmers, increasing with $\beta$, and varying dramatically as the orientation $\alpha$ of the flow is varied from 0 to ${\rm \pi} /2$. A major finding is that, despite very different assumptions, the two methods of computation give overlapping results for viscosity as a function of $\phi$ in the range $0.5 < \phi < 0.75$. This suggests that lubrication theory, based on near-field interactions alone, contains most of the relevant physics, and that taking account of interactions with more distant particles than the nearest is not essential to describe the dominant physics.

We consider the solidification of idealised two-component mixtures comprising a solvent or suspending fluid and dissolved solute molecules or suspended colloidal particles, each considered as hard spheres. We review some fundamental thermodynamic ideas regarding relative motion between species and phase equilibria in such mixtures to show how the related solid–liquid phase diagrams depend on the size of the spheres. Using similarity solutions, we first describe freezing of the solvent to form a pure solid (here referred to as ‘ice’), with the solute rejected from the solid forming a boundary layer or dense particle layer ahead of the freezing front. We extend ideas of constitutional supercooling to the case of colloidal suspensions and show that, for a given temperature difference driving solidification, constitutional supercooling occurs only for an intermediate range of particle sizes. Constitutional supercooling promotes the formation of a mushy layer in which segregated ice separates regions of concentrated solute or particles on the microscale. We formulate a continuum model of the mushy layer that relies on a key observation that the regelative motion of concentrated clusters of particles is independent of the size and geometry of the cluster. Our modelling begins with a description of relative motion as a Fickian diffusive process. However, at high particle concentrations, we show that it is more convenient and more computationally tractable to use an equivalent formulation in terms of Darcy flow of the solvent. Within a mushy layer these diffusive fluxes correspond directly to the regelative flux of particle clusters at a rate determined by the local temperature and temperature gradient.

Motivated by the desire to understand complex transient behaviour in fluid flows, we study the dynamics of an air bubble driven by the steady motion of a suspending viscous fluid within a Hele-Shaw channel with a centred depth perturbation. Using both experiments and numerical simulations of a depth-averaged model, we investigate the evolution of an initially centred bubble of prescribed volume as a function of flow rate and initial shape. The experiments exhibit a rich variety of organised transient dynamics, involving bubble breakup as well as aggregation and coalescence of interacting neighbouring bubbles. The long-term outcome is either a single bubble or multiple separating bubbles, positioned along the channel in order of increasing velocity. Up to moderate flow rates, the life and fate of the bubble are reproducible and can be categorised by a small number of characteristic behaviours that occur in simply connected regions of the parameter plane. Increasing the flow rate leads to less reproducible time evolutions with increasing sensitivity to initial conditions and perturbations in the channel. Time-dependent numerical simulations that allow for breakup and coalescence are found to reproduce most of the dynamical behaviour observed experimentally, including enhanced sensitivity at high flow rate. An unusual feature of this system is that the set of steady and periodic solutions can change during temporal evolution because both the number of bubbles and their size distribution evolve due to breakup and coalescence events. Calculation of stable and unstable solutions in the single- and two-bubble cases reveals that the transient dynamics is orchestrated by weakly unstable solutions of the system that can appear and disappear as the number of bubbles changes.

Direct numerical simulations have been performed for turbulent thermal convection between horizontal no-slip, permeable walls with a distance $H$ and a constant temperature difference $\Delta T$ at the Rayleigh number $Ra=3\times 10^3\text {--}10^{10}$. On the no-slip wall surfaces $z=0$, $H$, the wall-normal (vertical) transpiration velocity is assumed to be proportional to the local pressure fluctuation, i.e. $w=-\beta p'/\rho$, $+\beta p'/\rho$ (Jiménez et al., J. Fluid Mech., vol. 442, 2001, pp. 89–117). Here $\rho$ is mass density, and the property of the permeable wall is given by the permeability parameter $\beta U$ normalised with the buoyancy-induced terminal velocity $U=(g\alpha \Delta TH)^{1/2}$, where $g$ and $\alpha$ are acceleration due to gravity and volumetric thermal expansivity, respectively. The critical transition of heat transfer in convective turbulence has been found between the two $Ra$ regimes for fixed $\beta U=3$ and fixed Prandtl number $Pr=1$. In the subcritical regime at lower $Ra$ the Nusselt number $Nu$ scales with $Ra$ as $Nu\sim Ra^{1/3}$, as commonly observed in turbulent Rayleigh–Bénard convection. In the supercritical regime at higher $Ra$, on the other hand, the ultimate scaling $Nu\sim Ra^{1/2}$ is achieved, meaning that the wall-to-wall heat flux scales with $U\Delta T$ independent of the thermal diffusivity, although the heat transfer on the wall is dominated by thermal conduction. In the supercritical permeable case, large-scale motion is induced by buoyancy even in the vicinity of the wall, leading to significant transpiration velocity of the order of $U$. The ultimate heat transfer is attributed to this large-scale significant fluid motion rather than to transition to turbulence in boundary-layer flow. In such ‘wall-bounded’ convective turbulence, a thermal conduction layer still exists on the wall, but there is no near-wall layer of large change in the vertical velocity, suggesting that the effect of the viscosity is negligible even in the near-wall region. The balance between the dominant advection and buoyancy terms in the vertical Boussinesq equation gives us the velocity scale of $O(U)$ in the whole region, so that the total energy budget equation implies the Taylor dissipation law $\epsilon \sim U^3/H$ and the ultimate scaling $Nu\sim Ra^{1/2}$.

In this work, model closures of the multiphase Reynolds-averaged Navier–Stokes (RANS) equations are developed for homogeneous, fully developed gas–particle flows. To date, the majority of RANS closures are based on extensions of single-phase turbulence models, which fail to capture complex two-phase flow dynamics across dilute and dense regimes, especially when two-way coupling between the phases is important. In the present study, particles settle under gravity in an unbounded viscous fluid. At sufficient mass loadings, interphase momentum exchange between the phases results in the spontaneous generation of particle clusters that sustain velocity fluctuations in the fluid. Data generated from Eulerian–Lagrangian simulations are used in a sparse regression method for model closure that ensures form invariance. Particular attention is paid to modelling the unclosed terms unique to the multiphase RANS equations (drag production, drag exchange, pressure strain and viscous dissipation). A minimal set of tensors is presented that serve as the basis for modelling. It is found that sparse regression identifies compact, algebraic models that are accurate across flow conditions and robust to sparse training data.

Electrohydrodynamics of drops is a classic fluid mechanical problem where deformations and microscale flows are generated by application of an external electric field. In weak fields, electric stresses acting on the drop surface drive quadrupolar flows inside and outside and cause the drop to adopt a steady axisymmetric shape. This phenomenon is best explained by the leaky-dielectric model under the premise that a net surface charge is present at the interface while the bulk fluids are electroneutral. In the case of dielectric drops, increasing the electric field beyond a critical value can cause the drop to start rotating spontaneously and assume a steady tilted shape. This symmetry-breaking phenomenon, called Quincke rotation, arises due to the action of the interfacial electric torque countering the viscous torque on the drop, giving rise to steady rotation in sufficiently strong fields. Here, we present a small-deformation theory for the electrohydrodynamics of dielectric drops for the complete Melcher–Taylor leaky-dielectric model in three dimensions. Our theory is valid in the limits of strong capillary forces and highly viscous drops and is able to capture the transition to Quincke rotation. A coupled set of nonlinear ordinary differential equations for the induced dipole moments and shape functions are derived whose solution matches well with experimental results in the appropriate small-deformation regime. Retention of both the straining and rotational components of the flow in the governing equation for charge transport enables us to perform a linear stability analysis and derive a criterion for the applied electric field strength that must be overcome for the onset of Quincke rotation of a viscous drop.

For a small sessile or pendant droplet it is generally assumed that gravity does not play any role once the Bond number is small. This is even assumed for evaporating binary sessile or pendant droplets, in which convective flows can be driven due to selective evaporation of one component and the resulting concentration and thus surface tension differences at the air–liquid interface. However, recent studies have shown that in such droplets gravity indeed can play a role and that natural convection can be the dominant driving mechanism for the flow inside evaporating binary droplets (Edwards et al., Phys. Rev. Lett., vol. 121, 2018, 184501; Li et al., Phys. Rev. Lett., vol. 122, 2019, 114501). In this study, we derive and validate a quasi-stationary model for the flow inside evaporating binary sessile and pendant droplets, which successfully allows one to predict the prevalence and the intriguing interaction of Rayleigh and/or Marangoni convection on the basis of a phase diagram for the flow field expressed in terms of the Rayleigh and Marangoni numbers.

Circular milling, a stunning manifestation of collective motion, is found across the natural world, from fish shoals to army ants. It has been observed recently that the plant-animal worm Symsagittifera roscoffensis exhibits circular milling behaviour, both in shallow pools at the beach and in Petri dishes in the laboratory. Here we investigate this phenomenon experimentally and theoretically, from a fluid dynamical viewpoint, focusing on the effect that an established circular mill has on the surrounding fluid. Unlike systems such as confined bacterial suspensions and collections of molecular motors and filaments that exhibit spontaneous circulatory behaviour, and which are modelled as force dipoles, the front–back symmetry of individual worms precludes a stresslet contribution. Instead, singularities such as source dipoles and Stokes quadrupoles are expected to dominate. We analyse a series of models to understand the contributions of these singularities to the azimuthal flow fields generated by a mill, in light of the particular boundary conditions that hold for flow in a Petri dish. A model that treats a circular mill as a rigid rotating disc that generates a Stokes flow is shown to capture basic experimental results well, and gives insights into the emergence and stability of multiple mill systems.

We revisit the dynamics of a thick liquid film flowing down a vertical fibre. Instead of deriving a long-wave model, we directly solve the Navier–Stokes equations using a domain mapping technique and the exact steady travelling wave solutions are explored using a dynamical systems theory approach. Three distinct flow regimes, labelled ‘a’, ‘b’ and ‘c’, observed in previous experiments (Kliakhandler et al., J. Fluid Mech., vol. 429, 2001, pp. 381–390) are investigated. Flow regime ‘a’ refers to a steady flow state in which large droplets are separated by a long thin film. Flow regime ‘b’ is a necklace-like flow. In flow regime ‘c’, a cyclic process of droplet coalescence and breakup was observed. By matching the mean flow rates of the travelling wave solutions and experimental data, our travelling wave solutions show an excellent agreement with flow regimes ‘a’ and ‘b’. The time-periodic flow regime ‘c’ is compared with direct numerical simulation of the Navier–Stokes equations. A snapshot of the simulation shows a remarkable similarity to an experimental image and the discrepancy of mean wave speed and maximal wave height between our numerical simulation and experimental data is negligible.

Rotating Rayleigh–Bénard convection is investigated numerically with the use of an asymptotic model that captures the rapidly rotating, small Ekman number limit, $Ek \rightarrow 0$. The Prandtl number ($Pr$) and the asymptotically scaled Rayleigh number ($\widetilde {Ra} = Ra Ek^{4/3}$, where $Ra$ is the typical Rayleigh number) are varied systematically. For sufficiently vigorous convection, an inverse kinetic energy cascade leads to the formation of a pair of large-scale vortices of opposite polarity, in agreement with previous studies of rapidly rotating convection. With respect to the kinetic energy, we find a transition from convection dominated states to a state dominated by large-scale vortices at an asymptotically reduced (small-scale) Reynolds number of $\widetilde {Re} \approx 6$ ($\widetilde {Re} = Re Ek^{1/3}$, where $Re$ is the Reynolds number associated with vertical flows) for all investigated values of $Pr$. The ratio of the depth-averaged kinetic energy to the kinetic energy of the convection reaches a maximum at $\widetilde {Re} \approx 24$, then decreases as $\widetilde {Ra}$ is increased. This decrease in the relative kinetic energy of the large-scale vortices is associated with a decrease in the convective correlations with increasing Rayleigh number. The scaling behaviour of the convective flow speeds is studied; although a linear scaling of the form $\widetilde {Re} \sim \widetilde {Ra}/Pr$ is observed over a limited range in Rayleigh number and Prandtl number, a clear departure from this scaling is observed at the highest accessible values of $\widetilde {Ra}$. Calculation of the forces present in the governing equations shows that the ratio of the viscous force to the buoyancy force is an increasing function of $\widetilde {Ra}$, that approaches unity over the investigated range of parameters.

The phenomena of self-sustained shock-wave oscillations over conical bodies with a blunt axisymmetric base subject to uniform high-speed flow are investigated in a hypersonic wind tunnel at Mach number $M = 6$. The flow and shock-wave dynamics is dictated by two non-dimensional geometric parameters presented by the three length scales of the body, two of which are associated with the conical forebody and one with the base. Time-resolved schlieren imagery from these experiments reveals the presence of two disparate states of shock-wave oscillations in the flow, and allows for the mapping of unsteadiness boundaries in the two-parameter space. Physical mechanisms are proposed to explain the oscillations and the transitions of the shock-wave system from steady to oscillatory states. In comparison with the canonical single-parameter problem of shock-wave oscillations over spiked-blunt bodies reported in literature, the two-parameter nature of the present problem introduces distinct elements to the flow dynamics.

The dynamics of intense vorticity is investigated by means of synchronisation experiments in direct numerical simulations of isotropic turbulence. By imposing similar dynamics above the dissipative range, the same structures of intense vorticity appear in two independent turbulent flows, showing that intense vorticity synchronises to large-scale dynamics. Remarkably, this synchronisation takes place despite the presence of chaos, and affects mostly the intense vorticity, but not so much the weak vorticity background, which remains comparatively asynchronous. These results pinpoint the role of large-scale dynamics in the formation of intense vorticity structures, the so-called ‘worms’, and rule out the possibility that they emerge primarily due to interactions within the dissipative range, and then grow or coalesce into elongated structures. The stretching of the vorticity vector by the large-scale rate-of-strain tensor is identified as the mechanism responsible for the synchronisation of intense vorticity, supporting the extended view of vortex stretching as a fundamental inter-scale mechanism in turbulence.