Direct numerical simulation (DNS) is performed to explore turbulent Rayleigh–Bénard convection in spherical shells. Our simulations cover six distinct values of radius ratio, $\eta = r_i/r_o = 0.2$, 0.3, 0.4, 0.5, 0.6 and 0.8, under the assumption of a centrally condensed mass with gravity profile $g \sim 1/r^{2}$; where $r_i$, $r_o$ and $r$ denote the inner shell radius, the outer shell radius and the local radial coordinate, respectively. The Prandtl number is kept constant at unity while the Rayleigh number ($Ra$) is varied from $3 \times 10^{3}$ to $5 \times 10^8$. Our primary aim is to analyze how the radius ratio influences the global transport properties and flow physics. To gain insights into the scaling behaviour of the Nusselt number ($Nu$) and the Reynolds number ($Re$) with respect to $Ra$ and $\eta$, we apply the Grossmann–Lohse (GL) theory (Grossmann & Lohse, J. Fluid Mech., vol. 407, 2000, pp. 27–56) to the system. It is observed that the scaling exponents for $Nu$ and $Re$ in relation to $Ra$ are more significant for smaller $\eta$ values, suggesting that the simulations with smaller $\eta$ reach the classical $Nu\sim Ra^{1/3}$ regime at a relatively lower $Ra$. This observation could also imply the systems with smaller $\eta$ might transition to the ultimate regime earlier at a smaller $Ra$. Based on our extensive DNS data, we establish that the thickness of the inner thermal boundary, $\lambda _{\vartheta }^{i}$, follows a scaling relationship of $\lambda _{\vartheta }^{i} \sim \eta ^{1/2}$. This relationship, in turn, leads to a scaling law for $Nu$ in the form of $Nu \sim f(\eta ) Ra^{\gamma }$, where the function $f(\eta )$ is defined as $f(\eta ) = {\eta ^{1/2}}/{(1+\eta ^{4/3})}$, and the exponent $\gamma$ depends on both $Ra$ and $\eta$. Additionally, we characterize and explain the asymmetry in the velocity field by introducing the separate Reynolds numbers for the inner and outer shells. The asymmetry of the kinetic and thermal energy dissipation rates in the inner and outer boundary layers (BLs) is also quantified.

By comparing the budget of a data-driven quasi-linear approximation (DQLA) (Holford, Lee & Hwang, J. Fluid Mech., vol. 980, 2024, A12) and direct numerical simulation (DNS) (Lee & Moser, J. Fluid Mech., vol. 860, 2019, pp. 886–938), the energetics of linear models for wall-bounded turbulence are assessed. The DQLA is implemented with the linearised Navier–Stokes equations with a stochastic forcing term and an eddy viscosity diffusion model. The self-consistent nature of the DQLA allows for a global comparison across all wavenumbers to assess the role of the various terms in the linear model in replicating the features present in DNS. Starting from the steady-state second-order statistics of a Fourier mode, a spectral budget equation is derived, connecting Lyapunov-like equations to the transport budget equations obtained from DNS. It is found that the DQLA and DNS are in good qualitative agreement for the streamwise-elongated structures present in DNS, comparing well for production, viscous transport and wall-normal turbulent transport. However, the DQLA does not have an energy-conservative nonlinear term. This results in no dissipation under molecular viscosity, with energy instead being dissipated locally through the eddy viscosity model, which models the energy removal by the nonlinear term at integral length scales. Comparison of the pressure–strain statistics also highlights the absence of the streak instability, with production and forcing mainly being retained in the streamwise and wall-normal components or shifted to the spanwise component. It is demonstrated that the eddy viscosity diffusion term locally enforces a self-similar budget, making the model for the nonlinear term self-consistent with a logarithmic mean profile. Implications and recommendations to improve the current eddy viscosity enhanced linear models are also discussed concerning the comparison with DNS, as well as considerations with regard to pressure statistics to mimic the role of the streak instability through colour of turbulence models.

This DNS study considers transition to turbulence in plane Couette flow (pCf) with a rough stationary wall and a smooth moving wall. The roughness elements are square ribs of height $k=0.2h$ (where $h$ is the half-channel height). Two different pitch separations, $\lambda =2k$ and $10k$, are considered, i.e. d-type and k-type roughness, respectively. The transition in both rough pCf cases takes place through a stage of alternate laminar–turbulent bands aligned in an oblique fashion. However, roughness causes a shift in the transitional Reynolds number ($Re$) range. In the k-type roughness, stable bands are observed in the range $Re \in [300, 325]$, which is a downwards shift from the transitional $Re$ range for the smooth pCf ($Re \in [325,400]$). The d-type roughness, on the other hand, surprisingly shifts the transitional $Re$ range upwards to $Re \in [350,425]$. This peculiar behaviour is associated with the ability of the ribs to shed and regenerate vorticity. Large-scale components extracted using a filtering process relate to the transition bands and flow parallel to the oblique laminar–turbulent boundaries. Counter-rotating vortices are present in the turbulent regions of the flow field, which exist in tandem with the high- and low-velocity streaks. Another interesting observation is the secondary Reynolds shear stresses, $-\overline {v^{\prime }u^{\prime }}$ and $-\overline {w^{\prime }v^{\prime }}$, which are non-zero in the transitional regime, in contrast to the turbulent regime where they are negligible.

Extensive three-dimensional boundary-integral simulations are presented for the steady-state, low-Reynolds-number motion of a non-wetting deformable drop in another liquid on an inclined solid wall. The drop remains separated from the wall by a lubricating film. The boundary-integral formulation is based on the half-space Green function. The focus is on the challenging case of small tilt angles $\theta$ combined with high drop-to-medium viscosity ratios $\lambda$, when the drop travels with strong hydrodynamical interaction very close to the wall. Simulations to steady state have required ultrahigh drop surface resolutions (to $3\times 10^5$ boundary elements) achieved through multipole acceleration and combined with novel regularization to fully eliminate near-singular behaviour of the double-layer integrals due to small clearance. Non-dimensional drop speed $U$ is presented for $\theta \geq 7.5^\circ$, $\lambda \leq 300$ and in a broad range of Bond numbers $B$, covering from nearly spherical to strongly pancaked drops. The results are consistent with published experiments on liquid–liquid systems. At small $\theta$ and $\lambda \gg 1$, $U$ is a strong, decreasing function of $B$; the asymptotic regime $U\to 0$ at $B\to 0$ is not observed in the simulated range. For small $B$, the dimpled thin-film geometry is insensitive to $\lambda =1\unicode{x2013}300$. For pancaked drops, the lubrication film is much thicker for $\lambda =1$ than for $\lambda \gg 1$ drops. Approximate thin-film uniformity in the drop motion direction is confirmed for pancaked, but not for $B\ll 1$, drops. Kinematics of drop motion shows that neither perfect tank treading, nor perfect rolling can be approached for liquid–liquid systems in the purely hydrodynamical formulation. The methodology is applicable to other problems and can allow for direct inclusion of short-range colloidal forces in three-dimensional boundary-integral simulations.

We present herein the derivation of a lubrication-mediated (LM) quasi-potential model for droplet rebounds off deep liquid baths, assuming the presence of a persistent dynamic air layer which acts as a lubricating pressure transfer. We then present numerical simulations of the LM model for axisymmetric rebounds of solid spheres and compare quantitatively to current results in the literature, including experimental data in the low-speed impact regime. In this regime the LM model has the advantage of being far more computationally tractable than direct numerical simulation (DNS) and is also able to provide detailed behaviour within the micro-metric thin lubrication region. The LM system has an interesting mathematical structure, with the lubrication layer providing a free-boundary elliptic problem mediating the drop and bath free-boundary evolutionary equations.

The effect of microbial activity on buoyancy-driven flow within a porous layer is analysed. The input fluid provides an energy source for the growth of biofilms on the porous rock. At each location within the porous layer, the porosity and permeability begin to decrease once the input fluid has invaded. This leads to an evolving rock heterogeneity that depends on the passing time of the input fluid. Hence, the evolution of the flow is partly controlled by its own history. We present an axisymmetric gravity current model, accounting for this effect. In general, a reduction in permeability leads to the flow having a lesser extent in the radial direction and greater thickness (extent in the cross-flow direction), whilst a reduction in porosity has negligible effect on the thickness but leads to a much greater radial extent. The flow is fastest near the free surface where the permeability is greatest. In the case where the porosity and permeability reduce as power-law functions of fluid residence time, the evolution of the flow and the rock properties are self-similar. Consumption of the input fluid by the microbes is also incorporated in the model and it generally leads to flows with lesser radial extent but little change in the thickness. The three impacts of microbial growth (volume loss owing to consumption and the reduction in permeability and porosity) each influence the flow in substantially different ways and the interplay is analysed. A motivation of the study, the underground storage of hydrogen, is briefly discussed.

Simulations of elastic turbulence, the chaotic flow of highly elastic and inertialess polymer solutions, are plagued by numerical difficulties: the chaotically advected polymer conformation tensor develops extremely large gradients and can lose its positive-definiteness, which triggers numerical instabilities. While efforts to tackle these issues have produced a plethora of specialized techniques – tensor decompositions, artificial diffusion, and shock-capturing advection schemes – we still lack an unambiguous route to accurate and efficient simulations. In this work, we show that even when a simulation is numerically stable, maintaining positive-definiteness and displaying the expected chaotic fluctuations, it can still suffer from errors significant enough to distort the large-scale dynamics and flow structures. We focus on two-dimensional simulations of the Oldroyd-B and FENE-P equations, driven by a large-scale cellular body forcing. We first compare two positivity-preserving decompositions of the conformation tensor: symmetric square root (SSR) and Cholesky with a logarithmic transformation (Cholesky-log). While both simulations yield chaotic flows, only the latter preserves the pattern of the forcing, i.e. its fluctuating vortical cells remain ordered in a lattice. In contrast, the SSR simulation exhibits distorted vortical cells that shrink, expand and reorient constantly. To identify the accurate simulation, we appeal to a hitherto overlooked mathematical bound on the determinant of the conformation tensor, which unequivocally rejects the SSR simulation. Importantly, the accuracy of the Cholesky-log simulation is shown to arise from the logarithmic transformation. We also consider local artificial diffusion, a potential low-cost alternative to high-order advection schemes. Unfortunately, the artificially enhanced diffusive smearing of polymer stress in regions of intense stretching substantially modifies the global dynamics. We then show how the spurious large-scale motions, identified here, contaminate predictions of scalar mixing. Finally, we discuss the effects of spatial resolution, which controls the steepness of gradients in a non-diffusive simulation.

We propose a theoretical method to decompose the solution of a Stokes flow past a body immersed in a confined fluid into two simpler problems, related separately to the two geometrical elements of these systems: (i) the body immersed in the unbounded fluid (represented by its Faxén operators); and (ii) the domain of the confinement (represented by its Stokesian multipoles). Specifically, by using a reflection method, and assuming linear and reciprocal boundary conditions (Procopio & Giona, Phys. Fluids, vol. 36, issue 3, 2024, 032016), we provide the expression for the velocity field, the forces, torques and higher-order moments acting on the body in terms of: (i) the volume moments of the body in the unbounded ambient flow; (ii) the multipoles in the domain of the confinement; (iii) the collection of all the volumetric moments on the body immersed in all the regular parts of the multipoles considered as ambient flows. A detailed convergence analysis of the reflection method is developed. In light of the practical applications, we estimate the truncation error committed by considering only the lower-order moments (thus, truncating the matrices) and the errors associated with the approximated expressions available in the literature for force and torques. We apply the theoretical results to the archetypal hydrodynamic system of a sphere with Navier-slip boundary conditions near a plane wall with no-slip boundary conditions, to determine forces and torques on a translating and rotating sphere as a function of the slip length and of the distance of the sphere from the plane. The hydromechanics of a spheroid is also addressed.

Nonlinear electrokinetic phenomena, where electrically driven fluid flows depend nonlinearly on the applied voltage, are commonly encountered in aqueous suspensions of colloidal particles. A prime example is the induced-charge electro-osmosis, driven by an electric field acting on diffuse charge induced near a polarizable surface. Nonlinear electrohydrodynamic flows also occur in non-polar fluids, driven by the electric field acting on space charge induced by conductivity gradients. Here, we analyse the flows about a charge-neutral spherical solid particle in an applied uniform electric field that arise from conductivity dependence on local field intensity. The flow pattern varies with particle conductivity: while the flow about a conducting particle has a quadrupolar pattern similar to induced-charge electro-osmosis, albeit with opposite direction, the flow about an insulating particle has a more complex structure. We find that this flow induces a force on a particle near an electrode that varies non-trivially with particle conductivity: while it is repulsive for perfectly insulating particles and particles more conductive than the suspending medium, there exists a range of particle conductivities where the force is attractive. The force decays as the inverse square of the distance to the electrode and thus can dominate the dielectrophoretic attraction due to the image dipole, which falls off with the fourth power with the distance. This electrohydrodynamic lift opens new possibilities for colloidal manipulation and driven assembly by electric fields.

First predicted by Richtmyer in 1960 and experimentally confirmed by Meshkov in 1969, the Richtmyer–Meshkov instability (RMI) is crucial in fields such as physics, astrophysics, inertial confinement fusion and high-energy-density physics. These disciplines often deal with strong shocks moving through condensed materials or high-pressure plasmas that exhibit non-ideal equations of state (EoS), thus requiring theoretical models with realistic fluid EoS for accurate RMI simulations. Approximate formulae for asymptotic growth rates, like those proposed by Richtmyer, are helpful but rely on heuristic prescriptions for compressible materials. These prescriptions can sometimes approximate the RMI growth rate well, but their accuracy remains uncertain without exact solutions, as the fully compressible RMI growth rate is influenced by both vorticity deposited during shock refraction and multiple sonic wave refractions. This study advances previous work by presenting an analytic, fully compressible theory of RMI for reflected shocks with arbitrary EoS. It compares theoretical predictions with heuristic prescriptions using ideal gas, van der Waals gas and three-term constitutive equations for simple metals, the latter being analysed with detailed and simplified ideal-gas-like EoS. We additionally offer an alternative explicit approximate formula for the asymptotic growth rate. The comprehensive model also incorporates the effects of constant-amplitude acoustic waves at the interface, associated with the D'yakov–Kontorovich instability in shocks.

In a vertical channel driven by an imposed horizontal temperature gradient, numerical simulations (Gao et al., Phys. Rev. E, vol. 88, 2013, 023010; Phys. Rev. E, vol. 91, 2015, 013006; Phys. Rev. E, vol. 97, 2018, 053107) have previously shown steady, time-periodic and chaotic dynamics. We explore the observed dynamics by constructing invariant solutions of the three-dimensional Oberbeck–Boussinesq equations, characterizing the stability of these equilibria and periodic orbits, and following the bifurcation structure of the solution branches under parametric continuation in Rayleigh number. We find that in a narrow vertically periodic domain of aspect ratio 10, the flow is dominated by the competition between three and four co-rotating rolls. We demonstrate that branches of three- and four-roll equilibria are connected and can be understood in terms of their discrete symmetries. Specifically, the $D_4$ symmetry of the four-roll branch dictates the existence of qualitatively different intermediate branches that themselves connect to the three-roll branch in a transcritical bifurcation due to $D_3$ symmetry. The physical appearance, disappearance, merging and splitting of rolls along the connecting branch provide a physical and phenomenological illustration of the equivariant theory of $D_3$–$D_4$ mode interaction. We observe other manifestations of the competition between three and four rolls, in which the symmetry in time or in the transverse direction is broken, leading to limit cycles or wavy rolls, respectively. Our work highlights the interest of combining numerical simulations, bifurcation theory and group theory, in order to understand the transitions between and origin of flow patterns.

Wind tunnel measurements of the incident turbulent velocity fields and axial forces on a horizontal axis turbine and porous disc analogues are reported. The models were tested in both a simulated atmospheric boundary layer (ABL) and in grid turbulence, allowing for a range of turbulence length scale to rotor diameter ratios to be considered. A theoretical framework to account for the combined effect of distortion and potential flow blocking in the induction zone is presented. In the case of very large length-scale turbulence to diameter ratios, where distortion effects are minimal, a quasi-steady approach is adopted for the effect of blocking. For the small length-scale ratio limit, the method is developed from the classical analyses for rapid distortion of turbulence and blockage from flow through a porous sheet of resistance. For general length-scale ratios, an efficient prediction method based on interpolation between the two length-scale ratio extremes is established. For very large length-scale ratios, a quasi-steady theory without distortion is appropriate for a rotor or disc in a simulated ABL. The small length-scale theory is applicable for tests conducted in grid turbulence. The results of the study can inform the prediction and interpretation of typical measurements of turbulence within the induction zone and the fluctuating loads on a rotor, at both prototype and full scale. This is of particular importance to fatigue load assessments.

The causal relevance of local flow conditions in open-channel turbulence is analysed using ensembles of interventional experiments in which the effect of perturbing the flow within a small cell is monitored at some future time. When this is done using the relative amplification of the perturbation energy, causality depends on the flow conditions within the cell before it is perturbed, and can be used as a probe of the flow dynamics. The key scaling parameter is the ambient shear, which is also the dominant diagnostic variable for wall-attached perturbations. Away from the wall, the relevant variables are the streamwise and wall-normal velocities. Causally significant cells are associated with sweeps that carry the perturbation towards the stronger shear near the wall, whereas irrelevant ones are associated with ejections that carry it towards the weaker shear in the outer layers. Causally significant and irrelevant cells are themselves organised into structures that share many characteristics with classical sweeps and ejections, such as forming spanwise pairs whose dimensions and geometry are similar to those of classical quadrants. At the wall, this is consistent with causally significant configurations in which a high-speed streak overtakes a low-speed one, and causally irrelevant ones in which the two streaks pull apart from each other. It is argued that this is probably associated with streak meandering.

During oscillatory wetting, a phase retardation emerges between contact angle variation and contact line velocity, presenting as a hysteresis loop in their correlation – an effect we term dynamic hysteresis. This phenomenon is found to be tunable by modifying the surface with different molecular layers. A comparative analysis of dynamic hysteresis, static hysteresis and contact line friction coefficients across diverse substrates reveals that dynamic hysteresis is not a result of dissipative effects but is instead proportionally linked to the static hysteresis of the surface. In the quest for appropriate conditions to model oscillatory contact line motion, we identify the generalized Hocking's linear law and modified generalized Navier boundary condition as alternative options for predicting realistic dynamic hysteresis.

Modular floating solar farms exhibit periodic open surface coverages due to the strip configuration of floating modules that support the photovoltaic (PV) panels on top. The associated modulations in the surface boundary layer and its turbulence characteristics are investigated in the present study under fully developed open channel flows. Different coverage percentages of 100 % (i.e. full cover), 60 %, 30 % and 0 % (i.e. open surface) were tested and measurements were obtained using particle image velocimetry. The results showed that the turbulence statistics are similar when the coverage decreases from 100 % to 60 %. However, with 30 %, both the turbulence intensities and Reynolds stresses increase substantially, reaching up to 50 % higher compared with the 100 % coverage, and the boundary layer thickness increases by more than 25 %. The local skin friction beneath the openings increases by 50 %. Analysis of spanwise vortices and premultiplied spectra indicates that the periodic coverage elongates the hairpin vortex packets and reduces their inclination angle, imposing limitations on sustainable coherent structures. At 30 %, flow detachment and smaller-scale vortices become dominant, reducing the mean velocities and increasing the turbulence intensities. Decreasing coverage percentage with flow detachment also shifts the energy transfer to higher wavenumbers, increasing energy dissipation and decreasing the bulk flow velocity. The kinetic energy and Reynolds stress carried by very large-scale motions decreases from 40 %–50 % with the 100 % and 60 % coverage to around 30 %–40 % with the 30 % coverage. Further research studies involving spanwise heterogeneity, higher Reynolds number and varying submergence of PV modules are needed for environmental considerations.

We study the effect of surface texture on an overlying turbulent flow for the case of textures made of an alternating slip/no-slip pattern, a common model for superhydrophobic surfaces, but also a particularly simple form of texture. For texture sizes $L^+ \gtrsim 25$, we have previously reported that, even though the texture effectively imposes homogeneous slip boundary conditions on the overlying, background turbulence, this is not its sole effect. The effective conditions only produce an origin offset on the background turbulence, which remains otherwise smooth-wall-like. For actual textures, however, as their size increases from $L^+ \gtrsim 25$ the flow progressively departs from this smooth-wall-like regime, resulting in additional shear Reynolds stress and increased drag, in a non-homogeneous fashion that could not be reproduced by the effective boundary conditions. This paper focuses on the underlying physical mechanism of this phenomenon. We argue that it is caused by the nonlinear interaction of the texture-coherent flow, directly induced by the surface topology, and the background turbulence, as it acts directly on the latter and alters it. This does not occur at the boundary where effective conditions are imposed, but within the overlying flow itself, where the interaction acts as a forcing on the governing equations of the background turbulence, and takes the form of cross-advective terms between the latter and the texture-coherent flow. We show this by conducting simulations where we remove the texture and introduce additional, forcing terms in the Navier–Stokes equations, in addition to the equivalent homogeneous slip boundary conditions. The forcing terms capture the effect of the nonlinear interaction on the background turbulence without the need to resolve the surface texture. We show that, when the forcing terms are derived accounting for the amplitude modulation of the texture-coherent flow by the background turbulence, they quantitatively capture the changes in the flow up to texture sizes $L^+ \approx 70{-}100$. This includes not just the roughness function but also the changes in the flow statistics and structure.

This paper explores active wake-flow control on a notchback Ahmed body using genetically inspired optimization. Hotwire and particle image velocimetry measurements record velocity data and flow structures in the wake. Pulsed jets at four actuation slots (two at the roof trailing edge, two at the side trailing edges) dynamically control the wake to minimize aerodynamic drag. The study achieves up to 9.2 % (without consideration of energy consumption) drag reduction, primarily by manipulating vortices from the roof rear end. The paper elucidates the underlying flow mechanism and evaluates various actuation strategies, highlighting how optimal control leads to reattachment of wake separation at the rear slant, diminishing the slant bubble and promoting downstream reattachment for enhanced drag reduction.

The resolvent analysis reveals the worst-case disturbances and the most amplified response in a fluid flow that can develop around a stationary base state. The recent work by Padovan et al. (J. Fluid Mech., vol. 900, 2020, A14) extended the classical resolvent analysis to the harmonic resolvent analysis framework by incorporating the time-varying nature of the base flow. The harmonic resolvent analysis can capture the triadic interactions between perturbations at two different frequencies through a base flow at a particular frequency. The singular values of the harmonic resolvent operator act as a gain between the spatiotemporal forcing and the response provided by the singular vectors. In the current study, we formulate the harmonic resolvent analysis framework for compressible flows based on the linearized Navier–Stokes equation (i.e. operator-based formulation). We validate our approach by applying the technique to the low-Mach-number flow past an airfoil. We further illustrate the application of this method to compressible cavity flows at Mach numbers of 0.6 and 0.8 with a length-to-depth ratio of $2$. For the cavity flow at a Mach number of 0.6, the harmonic resolvent analysis reveals that the nonlinear cross-frequency interactions dominate the amplification of perturbations at frequencies that are harmonics of the leading Rossiter mode in the nonlinear flow. The findings demonstrate a physically consistent representation of an energy transfer from slow-evolving modes toward fast-evolving modes in the flow through cross-frequency interactions. For the cavity flow at a Mach number of 0.8, the analysis also sheds light on the nature of cross-frequency interaction in a cavity flow with two coexisting resonances.

This work reports an experimental study of the turbulent entrainment into the planar wake of a circular cylinder, exposed to various turbulent backgrounds, from the near- to the far-field. The background turbulence features independently varying turbulence intensity and integral length scale, thereby rendering different turbulent/turbulent interfaces (TTIs) between the background and the primary flow (wake). Combined, simultaneous particle image velocimetry and planar laser induced fluorescence measurements were conducted to quantify the entrainment characteristics across these various TTIs at an inlet Reynolds number of 3800. The primary focus was on understanding how turbulent entrainment evolves spatially in conjunction with the rapid development of the large-scale coherent vortices in the planar wake, and how such evolution is affected by the background turbulence. It is found that TTIs can establish two layers when the background turbulence is sufficiently intense, which distinguishes TTIs from the turbulent/non-turbulent interface (TNTI). The two layers are underpinned by different physical mechanisms but have the same thickness and appear to scale with the local Kolmogorov length scale after the wake spreading transition position (Chen & Buxton, J. Fluid Mech., vol. 969, 2023, A4). It is also found that the probability density functions of the entrainment velocity for both TTIs and a TNTI display power law tails, which are associated with extremely large entrainment velocities occurring more frequently than for a Gaussian process. These intermittent, extreme entrainment velocities make a remarkable contribution to the mean entrainment velocity, particularly in the near wake, which leads to a much higher mean entrainment velocity than farther downstream, for both a TNTI and the TTIs. Conditionally averaged analysis reveals that these extreme events of the entrainment velocity are directly associated with intense enstrophy structures close to the interface.

We consider the steady flow of a viscoelastic film over an inclined plane featuring periodic trenches normal to the main flow direction. The trenches have a square cross-section and side length 5–8 times the capillary length. Owing to the orientation of the substrate, the film fails to coat the topographical feature entirely, forming a second gas–liquid interface inside the trench and two three-phase contact lines at the points where the free surface meets the wall of the trench. The volume of entrapped air depends on material and flow parameters and geometric conditions. We develop a computational model and carry out detailed numerical simulations based on the finite element method to investigate this flow. We solve the two-dimensional mass and momentum conservation equations using the exponential Phan-Thien & Tanner constitutive model to account for the rheology of the viscoelastic material. Due to the strong nonlinearity, multiple steady solutions possibly connected by turning points forming hysteresis loops, transcritical bifurcations and isolated solution branches are revealed by pseudo-arc-length continuation. We perform a thorough parametric analysis to investigate the combined effects of elasticity, inertia, capillarity and viscosity on the characteristics of each steady flow configuration. The results of our simulations indicate that fluid inertia and elasticity may or may not promote wetting, while shear thinning and hydrophilicity always promote the wetting of the substrate. Interestingly, there are conditions under which the transition to the inertial regime is not smooth, but a hysteresis loop arises, signifying an abrupt change in the film hydrodynamics. Additionally, we investigate the effect of the geometrical characteristics of the substrate, and our results indicate that there is a unique combination of the geometry and viscoelastic properties that either maximizes or minimizes the wetting lengths.