The interaction of a swimmer with unsteady vortices in complex flows remains a topic of interest and open discussion. The present study, employing the immersed boundary method with a flexible fin model, explores swimming behaviours behind a circular cylinder with vortex-induced vibration (VIV). Five distinct swimming modes are identified on the $U_r$–$G_0$ plane, where $U_r$ denotes the reduced velocity, and $G_0$ represents the fin’s initial position. These modes include drifting upstream I/II (DU-I/II), Karman gait I/II (KG-I/II), and large oscillation (LO), with the DU-II, KG-II and LO modes being newly reported. The fin can either move around or cross through the vortex cores in the KG-I and KG-II modes, respectively, for energy saving and maintaining a stable position. When the upstream cylinder vibrates with its maximum amplitude, a double-row vortex shedding forms in the wake, allowing the DU-II mode to occur with the fin to achieve high-speed locomotion. This is attributed to a significant reduction in the streamwise velocity caused by vortex-induced velocity. Furthermore, a symmetry breaking is observed in the fin’s wake in the DU-II mode, potentially also contributing to high-speed locomotion. Overall, compared to the case without an upstream cylinder, we demonstrate that a self-propelled fin gains hydrodynamic advantages with various swimming modes in different VIV wakes. Interestingly, increased power transferred from flows by the oscillating cylinder leads to a more favourable environment for the downstream fin’s propulsion, indicating that a fin in VIV wakes obtains more advantages compared to the vortex street generated by a stationary cylinder.

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.

Turbulent mixing driven by the reshocked Richtmyer–Meshkov (RM) instability plays a critical role in numerous natural phenomena and engineering applications. As the most fundamental physical quantity characterizing the mixing process, the mixing width transitions from linear to power-law growth following the initial shock. However, there is a notable absence of quantitative models for predicting the pronounced compression of initial interface perturbations or mixing regions at the moment of shock impact. This gap has restricted the development of integrated algebraic models to only the pre- and post-shock evolution stages. To address this limitation, the present study develops a predictive model for the compression of the mixing width induced by shocks. Based on the general principle of growth rate decomposition proposed by Li et al. (Phy. Rev. E, vol. 103, issue 5, 2021, 053109), two distinct types of shock-induced compression processes are identified, differentiated by the dominant mechanism governing their evolution: light–heavy and heavy–light shock-induced compression. For light–heavy interactions, both stretching (compression) and penetration mechanisms are influential, whereas heavy–light interactions are governed predominantly by the stretching (compression) mechanism. To characterize these mechanisms, the average velocity difference between the extremities of the mixing zone is quantified, and a physical model of RM mixing is utilized. A quantitative theoretical model is subsequently formulated through the independent algebraic modelling of these two mechanisms. The proposed model demonstrates excellent agreement with numerical simulations of reshocked RM mixing, offering valuable insights for the development of integrated algebraic models for mixing width evolution.

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.

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.

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.

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.

We perform direct numerical simulations of centrifugal convection with an oscillating rotational velocity of small amplitude to study the effects of oscillatory boundary motion. The oscillation period is the main control parameter, with its range corresponding to a Womersley number in the range $1\lt Wo\lt 300$. Oscillating boundaries generate a circumferential shear flow, which significantly inhibits heat transfer, with maximum suppression $87\,\%$ observed in the present parameter space. Through analysis of the background flow, we find that as the oscillation period increases, the increasing penetration depth of the oscillation and weakening local shear strength result in non-monotonic changes in heat transfer. Under high-frequency oscillation, the characteristic length scale of the viscous layer induced by the oscillation is smaller than the convective length scale, and shear manifests primarily as a continuous suppression of the boundary layer. In contrast, under low-frequency oscillation, the shear flow covers the entire region but with weak strength. The suppression effect of such shear flow exhibits periodicity, leading to alternating phases of convection inhibition and convection generation. The present findings explore the physical mechanisms behind the suppression of convective heat transfer by oscillation, and offer a new strategy for controlling convection systems, with potential implications for both fundamental research and industrial applications.

A generalised multiparameter model for linear modal stability and sensitivity analysis is developed. The stability and sensitivity equations are derived from a generalised vector-form governing equation comprised of multiple dimensionless parameters that represent different physical forces affecting the system’s stability. By introducing adjoint variables and constructing the Lagrangian identity, a differential relationship between the eigenvalue of the perturbation mode and dimensionless parameters is determined and defined as the global sensitivity gradient. It provides the constraint that must be satisfied for changes in different dimensionless parameters along the isoeigenvalue curve, which aids in the fast computation of the neutral curve. Moreover, the global sensitivity gradient can directly and intuitively evaluate the competitive relationship among the influences of various parameters on system instability. Based on the global sensitivity gradient, an optimal stability control strategy for transitioning from an unstable state to a stable state is discussed. Additionally, the relative sensitivity function is also introduced to investigate the influence of relative parameter variations on instability. To demonstrate the effectiveness of this method, three applications are presented: two-dimensional flow around a circular cylinder with a single dimensionless parameter Re; three-dimensional axisymmetric magnetohydrodynamic (MHD) flow around a sphere with two parameters Re and $N$; and two-dimensional MHD mixed convection with three parameters Re, ${\textit{Gr}}$ and $\textit{Ha}$.

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 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.

This paper investigates the weakly nonlinear isotropic bidirectional Benney–Luke (BL) equation, which is used to describe oceanic surface and internal waves in shallow water, with a particular focus on soliton dynamics. Using the Whitham modulation theory, we derive the modulation equations associated with the BL equation that describe the evolution of soliton amplitude and slope. By analysing rarefaction waves and shock waves within these modulation equations, we derive the Riemann invariants and modified Rankine–Hugoniot conditions. These expressions help characterise the Mach expansion and Mach reflection phenomena of bent and reverse bent solitons. We also derive analytical formulae for the critical angle and the Mach stem amplitude, showing that as the soliton speed is in the vicinity of unity, the results from the BL equation align closely with those of the Kadomtsev–Petviashvili (KP) equation. Corresponding numerical results are obtained and show excellent agreement with theoretical predictions. Furthermore, as a far-field approximation for the forced BL equation – which models wave and flow interactions with local topography – the modulation equations yield a slowly varying similarity solution. This solution indicates that the precursor wavefronts created by topography moving at subcritical or critical speeds take the shape of a circular arc, in contrast to the parabolic wavefronts observed in the forced KP equation.

The recently proposed near-wall turbulence predictive model quantifies the degree of the superposition and the amplitude modulation exerted by large-scale coherent structures on small scales in the linear and nonlinear terms of the formula, respectively, and achieves the prediction of streamwise velocity in the inner region. However, the multiscale effect and the time shift confirmed in the amplitude modulation have not yet been simultaneously taken into account in the model, which could limit the prediction accuracy especially at high Reynolds numbers. In this study, the role of the nonlinear term in the model is clarified based on high-quality flow data obtained in atmospheric surface layers: it redistributes the energy of the universal signal in the time domain and determines the accuracy of the predictive odd moments. An analysis of the multiscale effect and the time shifts in the nonlinear term is subsequently conducted, followed by a demonstration of the refinement in the quality of the universal signal after separately incorporating them into the model. The amplitude modulation is revealed when the two factors are simultaneously considered, and profiles of the scales that dominate the modulation and time shifts with height is provided. Thus, the nonlinear term of the existing model is modified, proposing an polished scheme that can quantify the nonlinear modulation terms more accurately.

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.

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.