Viscous fingering, a classic hydrodynamic instability, is governed by the the competition between destabilising viscosity ratios and stabilising surface tension or thermal diffusion. We show that the channel confinement can induce ‘diffusion’-like stabilising effects on viscous fingering even in the absence of interfacial tension and thermal diffusion, when a clear oil invades the mixture of the same oil and non-colloidal particles. The key lies in the generation of long-range dipolar disturbance flows by highly confined particles that form a monolayer inside a Hele-Shaw cell. We develop a coarse-grained model whose results correctly predict universal fingering dynamics that is independent of particle concentrations. This new mechanism offers insights into manipulating and harnessing collective motion in non-equilibrium systems.

The impact of several ‘flavours’ of free-stream turbulence (FST) on the structural response of a cantilever cylinder, subjected to a turbulent cross-flow is investigated. At high enough Reynolds numbers, the cylinder generates a spectrally rich turbulent wake that contributes significantly to the experienced loads. The presence of FST introduces additional complexity through two primary mechanisms: directly, by imposing a fluctuating velocity field on the cylinder’s surface, and indirectly, by altering the vortex shedding dynamics, modifying the experienced loads. We employ concurrent temporally resolved particle image velocimetry and distributed strain measurements using Rayleigh backscattering fibre optic sensors to instrument the surrounding velocity field and the structural strain respectively. By using various turbulence-generating grids, and manipulating their distance to the cylinder, we assess a broad FST parameter space allowing us to explore individually the influence of the transverse integral length scale ($\mathcal{L}_{13}/D$) and turbulence intensity of the FST on the developing load dynamics. The FST enhances the magnitude of the loads acting on the cylinder. This results from a decreased vortex formation length, increased coherence of regular vortex shedding, and energy associated with this flow structure in the near wake. The cylinder’s structural response is driven mainly by the vortex shedding dynamics, and its modification induced by the presence of FST, i.e. the indirect effect outweighs the direct effect. From the explored FST parameter space, turbulence intensity was seen to be the main driver of enhanced loading conditions, presenting a positive correlation with the fluctuating loads magnitude at the root.

We study theoretically and experimentally the propagation of two bubbles in a Hele-Shaw cell under a uniform background flow. We consider the regime where the bubbles are large enough to be flattened by the cell walls into a pancake-like shape, but small enough such that each bubble remains approximately circular when viewed from above. In a system of two bubbles of different radii, if the smaller bubble is in front, it will be overtaken by the larger bubble. Under certain circumstances, the bubbles may avoid collision by rolling over one another while passing. We find that, for a given ratio of the bubble radii, there exists a critical value of a dimensionless parameter (the Bretherton parameter) above which the two bubbles will never collide, regardless of their relative size and initial transverse offset, provided they are initially well separated in the direction of the background flow. Additionally, we determine the corrections to the bubble shape from circular for two bubbles aligned with the flow direction. We find that the front bubble flattens in the flow direction, while the rear bubble elongates. These shape changes are associated with changes in velocity, which allow the rear bubble to catch the bubble in front even when they are of the same size.

From anthropogenic litter carried by ocean currents to plant stems travelling through the atmosphere, geophysical flows are often seeded with elongated, fibre-like particles. In this study, we used a large-scale laboratory model of a tidal current – representative of a widespread class of geophysical flows – to investigate the tumbling motion of long, slender and floating fibres in the complex turbulence generated by flow interactions with a tidal inlet. Despite the non-stationary, non-homogeneous and anisotropic nature of this turbulence, we find that long fibres statistically rotate at the same frequency as eddies of similar size, a phenomenon called scale selection, which is known to occur in ideal turbulence. Furthermore, we report that the signal of the instantaneous transverse velocity difference between the fibre ends changes significantly from the signal produced by the flow in the fibre surroundings, although the two are statistically equivalent. These observations have twofold implications. On the one hand, they confirm the reliability of using the end-to-end velocity signal of rigid fibres to probe the two-point transverse statistics of the flow, even under realistic conditions: oceanographers could exploit this observation to measure transverse velocity differences through elongated floats in the field, where superdiffusion complicates collecting sufficient data to probe two-point turbulence statistics at a fixed separation effectively. On the other hand, by addressing the dynamics of inertial range particles floating in the coastal zone, these observations are crucial to improving our ability to predict the fate of meso- and macro-litter, a size class that is currently understudied.

Complex materials with internal microstructure such as suspensions and emulsions exhibit time-dependent rheology characterised by viscoelasticity and thixotropy. In many large-scale applications such as turbulent pipe flow, the elastic response occurs on a much shorter time scale than the thixotropy, hence these flows are purely thixotropic. The fundamental dynamics of thixotropic turbulence is poorly understood, particularly the interplay between microstructural state, rheology and turbulence structure. To address this gap, we conduct direct numerical simulations (DNS) of fully developed turbulent pipe flow of a model thixotropic (Moore) fluid as a function of the thixoviscous number $\Lambda$, which characterises the thixotropic kinetic rate relative to turbulence eddy turnover time, ranging from slow ($\Lambda \ll 1$) to fast ($\Lambda \gg 1$) kinetics. Analysis of DNS results in the Lagrangian frame shows that, as expected, in the limits of slow and fast kinetics, these time-dependent flows behave as time-independent purely viscous (generalised Newtonian) analogues. For intermediate kinetics ($\Lambda \sim 1$), the rheology is governed by a path integral of the thixotropic fading memory kernel over the distribution of Lagrangian shear history, the latter of which is modelled via a simple stochastic model for the radially non-stationary pipe flow. The DNS computations based on this effective viscosity closure exhibit excellent agreement with the fully thixotropic model for $\Lambda =1$, indicating that the purely viscous (generalised Newtonian) analogue persists for arbitrary values of $\Lambda \in (0,\infty ^+)$ and across nonlinear rheology models. These results significantly simplify our understanding of turbulent thixotropic flow, and provide insights into the structure of these complex time-dependent flows.

Convective boundary layers are governed by an interplay of vertical turbulent convection and shear-driven turbulence. Here, we investigate vertical velocity and buoyancy fields in convective boundary layers for varying atmospheric conditions by combining probability density function methods and direct numerical simulations. The evolution equations for the probability density functions of vertical velocity and buoyancy contain unclosed terms in the form of conditional averages. We estimate these terms from our direct numerical simulations data, and discuss their physical interpretation. Furthermore, using the method of characteristics, we investigate how these unclosed terms jointly determine the average evolution of a fluid element in a convective boundary layer, and how it relates to the evolution of the probability density functions of vertical velocity and buoyancy as a function of height. Thereby, our work establishes a connection between the turbulent dynamics of convective boundary layers and the resulting statistics.

Predicting the temperature distribution in laminar two-phase flows is essential in a wide range of engineering applications, like heat dissipation of electronic equipment and thermal design of biological reactors. Motivated by this, we extend the classical Graetz problem, studying the heat transfer between two flowing phases in a core-annular flow configuration. Using a rigorous two-scale asymptotic analysis, we derived two coupled one-dimensional advection–diffusion heat-transfer equations (one for each phase) embedding the effects of advection, diffusion (both axial and transverse) and viscous dissipation. Specifically, the heat-transfer mechanisms are described through effective velocity and effective diffusion coefficients, while the interaction between the phases is accounted for via ad hoc coupling and source terms, respectively. The dynamics of the problem is controlled by seven dimensionless groups: the Péclet and Brinkman numbers, the heat flux, the viscosity, thermal diffusivity and thermal conductivity ratios, and the volume fraction. Our analysis reveals the existence of two main regimes, depending on the disparity in thermal conductivity between the phases. When the conductivity ratio is of order one, the problem is strongly coupled; otherwise, the phases are thermally decoupled. Interestingly, we investigate the evolution of the heat-transfer coefficient in the thin-film limit, shedding light on the most common assumptions underlying extensively used models in the context of film flows. Finally, we derived closed-form scaling laws for the Nusselt number clarifying the impact of the phases topology on heat-transfer dynamics. Since our model has been derived by first principles, we hope that it will improve the understanding of two-phase forced convection.

Power minimisation in branched fluidic networks has gained significant attention in biology and engineering. The optimal network is defined by channel radii that minimise the sum of viscous dissipation and the volumetric energetic cost of the fluid. For limit cases including laminar flows, high-Reynolds-number turbulence or smooth-channel approximations, optimal solutions are known. However, current methods do not allow optimisation for a large intermediate part of the parameter space which is typically encountered in realistic fluidic networks that exhibit turbulent flow. Here, we present a unifying optimisation approach based on the Darcy friction factor, which has been determined for a wide range of flow regimes and fluid models and is applicable to the entire parameter space: (i) laminar and turbulent flows, including networks that exhibit both flow types, (ii) non-Newtonian fluids (powerlaw, Bingham and Herschel–Bulkley) and (iii) networks with arbitrary wall roughness, including non-uniform relative roughness. The optimal channel radii are presented analytically and graphically. All existing limit cases are recovered, and a concise framework is presented for systematic optimisation of fluidic networks. Finally, the parameter $x$ in the optimisation relationship $Q\propto R^{x}$, with $Q$ the flow rate and $R$ the channel radius, was approximated as a function of the Reynolds number, revealing in which case the entire network can be optimised based on one optimal channel radius, and in which case all radii must be optimised individually. Our approach can be extended to a wide range of fluidic networks for which the friction factor is known, such as different channel curvatures, bubbly flows or specific wall slip conditions.

The scaling of pressure and vorticity in aquatic swimming can provide insights into the mechanisms of propulsion. This is investigated through self-propelled, wall-resolved, large-eddy simulations of a lamprey (an anguilliform swimmer) and a mackerel (a carangiform swimmer) using the curvilinear immersed boundary method. It is observed that the pressure around the swimmers scales with theoretical fluid acceleration, which includes both local body and the convective acceleration, for anguilliform swimmers, whereas it scales with both acceleration and the angle of attack (AoA) for carangiform swimmers. This indicates that the main mechanism for propulsion in anguilliform swimmers is added mass (unsteady), whereas both lift-based (steady) and added mass (unsteady) are at play for carangiform swimmers. Furthermore, it is observed that the vorticity in the boundary layer of the swimmer initially follows the body rotation at low speeds but not at high speeds during the quasisteady swimming. This is explained by identifying the scaling of vorticity components: one due to body rotation and the other due to shear, which scale with Strouhal number ($St$) and Reynolds number ($\sqrt {Re}$), respectively. Here $St$ (body rotation) dominates at low speeds, but $\sqrt {Re}$ (shear) dominates at high speeds. Finally, it is observed that the pressure decreases as the swimming speed increases. This counterintuitive observation is explained by showing that both fluid acceleration and AoA decrease as swimming speed increases. This suggests that for efficient swimming, the pressure difference across the body should be minimised, but high enough to overcome the viscous drag.

The instability characteristics and laminar–turbulent transition of a series of laminar separation bubbles (LSBs) formed due to a single sinusoidal surface waviness are investigated in the absence of external disturbances or forcing. A scaling based on the geometrical parameters of the waviness and flow Reynolds number is found that enables the prediction of flow separation on the wall leeward side. The analysis of three-dimensional instabilities of two-dimensional base flows reveals a relation between the number of changes in the curvature sign of the recirculating streamlines and the number of unstable centrifugal modes that coexist for the same flow. When multiple curvature changes occur, in addition to the usual steady mode reported for two-dimensional recirculation bubbles, a new self-excited mode with a higher growth rate emerges, localised near the highest streamline curvature, close to the reattachment point. A detailed analysis of the mode growth and saturation using DNS reveals that the localised mode only disturbs the LSB locally, while the usual one leads to a global distortion of the bubble in the spanwise direction; this has a distinctive impact on the self-excited secondary instabilities. Then, the complete transition scenario is studied for two selected LSB cases. The first one only presents an unstable eigenmode, namely the usual centrifugal mode in recirculating flows. The second case presents three unstable eigenmodes: two centrifugal eigenmodes (the usual and the localised ones) and a two-dimensional eigenmode associated with the self-sustained Kelvin–Helmholtz waves. These results show how completely different transition scenarios can emerge from subtle changes in the LSB characteristics.

Plugging of a hydraulic fracture because of particle bridging in the fracture channel is ubiquitous in drilling operations and reservoir stimulations. The particles transported in the fluid and fracture can aggregate under certain conditions and finally form a plug. The plug reduces the permeability of the flow channel and blocks the fluid pressure from reaching the fracture front, leading to fracture arrest. In this paper, a numerical model is developed to describe the plugging process of a hydraulic fracture driven by a slurry of solid particles in a viscous fluid while accounting for the rock deformation, slurry flow in the fracture channel, fracture propagation, particle transport and bridging. Three dimensionless numbers are derived from the governing equations, which reveal two length scales that control the fracture propagation and particle transport behaviour, respectively. The difference in magnitude between the two length scales implies three limiting regimes for fracture propagation, i.e. static regime, fluid-driven regime and slurry-driven regime, which correspond to fracture arrest, fracture driven by clean fluid, and fracture driven by slurry, respectively. Numerical results show that the fracture will sequentially transition through the static regime, fluid-driven regime and slurry-driven regime as the fracture length increases. The transition between regimes is controlled by the ratio between the two length scales. Simulation results also reveal two plugging modes, with the plug located near the fracture tip region and at the fracture inlet. The transition between the two plugging modes is controlled by the ratio of the length scales and the injected particle concentration.

The dispersion behaviour of solutes in flow is crucial to the design of chemical separation systems and microfluidics devices. These systems often rely on coupled electroosmotic and pressure-driven flows to transport and separate chemical species, making the transient dispersive behaviour of solutes highly relevant. However, previous studies of Taylor dispersion in coupled electroosmotic and pressure-driven flows focused on the long-term dispersive behaviour and the associated analyses cannot capture the transient behaviour of solute. Further, the radial distribution of solute has not been analysed. In the current study, we analyse the Taylor dispersion for coupled electroosmotic and pressure-driven flows across all time regimes, assuming a low zeta potential (electric potential at the shear plane), the Debye–Hückel approximation and a finite electric double layer thickness. We first derive analytical expressions for the effective dispersion coefficient in the long-time regime. We also derive an unsteady, two-dimensional (radial and axial) solute concentration field applicable in the latter regime. We next apply Aris’ method of moments to characterise the unsteady propagation of the mean axial position and the unsteady growth of the variance of the solute zone in all time regimes. We benchmark our predictions with Brownian dynamics simulations across a wide and relevant dynamical regime, including various time scales. Lastly, we derive expressions for the optimal relative magnitudes of electroosmotic versus pressure-driven flow and the optimum Péclet number to minimise dispersion across all time scales. These findings offer valuable insights for the design of chemical separation systems, including the optimisation of capillary electrophoresis devices and electrokinetic microchannels and nanochannels.

This study presents an automatic differentiation (AD)-based optimisation framework for flow control in compressible turbulent channel flows. Using a differentiable solver, JAX-Fluids, we designed fully differentiable boundary conditions that allow for the precise calculation of gradients with respect to boundary control variables. This facilitates the efficient optimisation of flow control methods. The framework’s adaptability and effectiveness are demonstrated using two boundary conditions: opposition control and tunable permeable walls. Various optimisation targets are evaluated, including wall friction and turbulent kinetic energy (TKE), across different time horizons. In each optimisation, there were around $4\times 10^4$ control variables and $3\times 10^{9}$ state variables in a single episode. Results indicate that TKE targeted opposition control achieves a more stable and significant reduction in drag, with effective suppression of turbulence throughout the channel. In contrast, strategies that focus directly on minimising wall friction were found to be less effective, exhibiting instability and increased turbulence in the outer region. The tunable permeable walls also show potential to achieve stable drag reduction through a ‘flux-inducing’ mechanism. This study demonstrates the advantages of AD-based optimisation in complex flow control scenarios and provides physical insight into the choice of the quantity of interest for improved optimisation performance.

A new lattice Boltzmann model (LBM) is presented to describe chemically reacting multicomponent fluid flow in homogenised porous media. In this work, towards further generalising the multicomponent reactive lattice Boltzmann model, we propose a formulation which is capable of performing reactive multicomponent flow computation in porous media at the representative elementary volume (REV) scale. To that end, the submodel responsible for interspecies diffusion has been upgraded to include Knudsen diffusion, whereas the kinetic equations for the species, the momentum and the energy have been rewritten to accommodate the effects of volume fraction of porous media through careful choice of the equilibrium distribution functions. Verification of the mesoscale kinetic system of equations by a Chapman–Enskog analysis reveals that at the macroscopic scale, the homogenised Navier–Stokes equations for compressible multicomponent reactive flows are recovered. The dusty gas model (DGM) capability hence formulated is validated over a wide pressure range by comparison of experimental flow rates of component species counter diffusing through capillary tubes. Next, for developing a capability to compute heterogeneous reactions, source terms for maintaining energy and mass balance across the fluid phase species and the surface adsorbed phase species are proposed. The complete model is then used to perform detailed chemistry simulations in porous electrodes of a solid oxide fuel cell (SOFC), thereby predicting polarisation curves which are of practical interest.

The spatiotemporal dynamics of a turbulent boundary layer subjected to an unsteady pressure gradient are studied. A dynamic sequence of favourable to adverse pressure gradients (FAPGs) is imposed by deforming a section of the wind tunnel ceiling, transitioning the pressure gradient from zero to a strong FAPG within 0.07 s. At the end of the transient, the acceleration parameter is $K$ = $6 \times 10^{-6}$ in the favourable pressure gradient (FPG) region and $K$ = $-4.8 \times 10^{-6}$ in the adverse pressure gradient (APG) region. The resulting unsteady response of the boundary layer is compared with equivalent steady pressure gradient cases in terms of turbulent statistics and coherent structures. While the steady FAPG effects, as shown by Parthasarathy & Saxton-Fox (2023), caused upstream stabilisation in the FPG, a milder APG response downstream, and the formation of an internal layer, the unsteady case presented in this paper shows a reduced stabilisation in the FPG region, a stronger APG response and a weaker internal layer. This altered response is hypothesised to stem from the different spatiotemporal pressure gradient histories experienced by turbulent structures when the pressure gradient changes at a time scale comparable to their convection.

We solve a Bayesian inverse Navier–Stokes (N–S) problem that assimilates velocimetry data by jointly reconstructing a flow field and learning its unknown N–S parameters. We devise an algorithm that learns the most likely parameters of a Carreau shear-thinning viscosity model, and estimates their uncertainties, from velocimetry data of a shear-thinning fluid. We conduct a magnetic resonance velocimetry experiment to obtain velocimetry data of an axisymmetric laminar jet in an idealised medical device (US Food and Drug Administration’s benchmark nozzle) for a blood analogue fluid. The algorithm successfully reconstructs the flow field and learns the most likely Carreau parameters. Predictions from the learned model agree well with rheometry measurements. The algorithm accepts any differentiable algebraic viscosity model, and can be extended to more complicated non-Newtonian fluids (e.g. Oldroyd-B fluid if a viscoelastic model is incorporated).

This paper presents a peridynamics-based computational approach for modelling coupled fluid flow and heat transfer problems. A new thermo-hydrodynamic peridynamics model is formulated with the semi-Lagrangian scheme and non-local operators. To enhance accuracy and numerical stability, a multi-horizon scheme is developed to introduce distinct horizons for the flow field and thermal field. The multi-horizon scheme helps to capture the convective zone and complex thermal flow pattern while effectively mitigating possible oscillations in temperature. We validate the computational approach using benchmarks and numerical examples including heat conduction, natural convection in a closed cavity, and Rayleigh–Bénard convection cells. The results demonstrate that the proposed method can accurately capture typical thermal flow behaviours and complex convective patterns. This work offers a new foundation for future development of a unified peridynamics framework for robust, comprehensive multi-physics analysis of thermal fluid–solid interaction problems with complex evolving discontinuities in solids.

The influence of parametric forcing on a viscoelastic fluid layer, in both gravitationally stable and unstable configurations, is investigated via linear stability analysis. When such a layer is vertically oscillated beyond a threshold amplitude, large interface deflections are caused by Faraday instability. Viscosity and elasticity affect the damping rate of momentary disturbances with arbitrary wavelength, thereby altering the threshold and temporal response of this instability. In gravitationally stable configurations, calculations show that increased elasticity can either stabilize or destabilize the viscoelastic system. In weakly elastic liquids, higher elasticity increases damping, raising the threshold for Faraday instability, whereas the opposite is observed in strongly elastic liquids. While oscillatory instability occurs in Newtonian fluids for all gravity levels, we find that parametric forcing below a critical frequency will cause a monotonic instability for viscoelastic systems at microgravity. Importantly, in gravitationally unstable configurations, parametric forcing above this frequency stabilizes viscoelastic fluids, until the occurrence of a second critical frequency. This result contrasts with the case of Newtonian liquids, where under the same conditions, forcing stabilizes a system for all frequencies below a single critical frequency. Analytical expressions are obtained under the assumption of long wavelength disturbances predicting the damping rate of momentary disturbances as well as the range of parameters that lead to a monotonic response under parametric forcing.

Feigenbaum universality is shown to occur in subcritical shear flows. Our testing ground is the counter-rotation regime of the Taylor–Couette flow, where numerical calculations are performed within a small periodic domain. The accurate computation of up to the seventh period-doubling bifurcation, assisted by a purposely defined Poincaré section, has enabled us to reproduce the two Feigenbaum universal constants with unprecedented accuracy in a fluid flow problem. We have further devised a method to predict the bifurcation diagram up to the accumulation point of the cascade based on the detailed inspection of just the first few period-doubling bifurcations. Remarkably, the method is applicable beyond the accumulation point, with predictions remaining valid, in a statistical sense, for the chaotic dynamics that follows.

Miscible Rayleigh–Taylor (RT) turbulence exhibits a wide range of length scales in both the velocity and density fields, leading to complex deformations of isoscalar surfaces and enhanced mixing due to nonlinear interactions among different scales. Through high-resolution numerical simulations and a coarse-graining analysis, we demonstrate that the variance of the heavy fluid concentration, initially maximised by the unstable stratification, progressively cascades from larger to smaller scales, eventually dissipates at the smallest scale. The transfer of scalar variance, $\Pi ^Y$, primarily governed by the filtered strain rate tensor, is effectively captured by a nonlinear model that links $\Pi ^Y$ to the isoscalar surface stretching. On the other hand, the backscatter of scalar variance transfer, represented by the negative component of $\Pi ^Y$, is influenced by the filtered vorticity field. Furthermore, we examine the directional anisotropy of scalar transfer in RT turbulence, enhancing the accuracy of the nonlinear model by separating the horizontal mean of the mass fraction from its fluctuating part.