An adaptable estimation technique is presented to reconstruct time-evolving three dimensional (3-D) velocity fields from planar particle image velocimetry measurements. The methodology builds on the multi-time-delay estimation technique of Hosseini et al. (2015) by implementing the finite-impulse-response spectral proper orthogonal decomposition (FIR-SPOD) of Sieber et al. (2016). The candidate flow is the highly modulated turbulent near wake of a cantilevered square cylinder with a height-to-width ratio $h/d=4$, protruding a thin laminar boundary layer ($\delta /d=0.21$ with $\delta$ being the boundary layer thickness) at the Reynolds number $Re=10600$, based on d. The novelty of the estimation technique is in using the modal space obtained by FIR-SPOD to better isolate the spatio-temporal scales for correlating velocity and pressure modes. Using FIR-SPOD, irregular coherent contributions at frequencies centred at $f_{ac1}=(1\pm 0.05)f_s$ and $f_{ac2}=(1\pm 0.1)f_s$ (with $f_s$ the fundamental shedding frequency) could be separated, which was not possible using proper orthogonal decomposition. With the FIR-SPOD bases, the quality of the estimation improved significantly using only linear terms, and the correct phase relationships between pressure and velocity modes are retained, as is required for synchronizing coherent motions along the height of the obstacle. It is shown that a low-dimensional reconstruction of the flow field successfully captures the cycle-to-cycle variations of the dominant 3-D vortex shedding process, which give rise to vortex dislocation events. Thus, the present methodology shows promise in 3-D reconstruction of challenging turbulent flows, which exhibit non-periodic behaviour or contain multi-scale phenomena.

The demand for separating and analysing rare target cells is increasing dramatically for vital applications such as cancer treatment and cell-based therapies. However, there remains a grand challenge for high-throughput and label-free segregation of lesion cells with similar sizes. Cancer cells with different invasiveness usually manifest distinct deformability. In this work, we employ a hydrogel microparticle system with similar sizes but varied stiffness to mimic cancer cells and examine in situ their deformation and focusing under microfluidic flow. We first demonstrate the similar focusing behaviour of hydrogel microparticles and cancer cells in confined flow that is dominated by deformability-induced lateral migration. The deformation, orientation and focusing position of hydrogel microparticles in microfluidic flow under different Reynolds numbers are then systematically observed and measured using a high-speed camera. Linear correlations of the Taylor deformation and tilt angle of hydrogel microparticles with the capillary number are revealed, consistent with theoretical predictions. Detailed analysis of the dependence of particle focusing on the flow rate and particle stiffness enables us to identify a linear scaling between the equilibrium focusing position and the major axis of the deformed microparticles, which is uniquely determined by the capillary number. Our findings provide insights into the focusing and dynamics of soft beads, such as cells and hydrogel microparticles, under confined flow, and pave the way for applications including the separation and identification of circulating tumour cells, drug delivery and controlled drug release.

We report pattern formation in an otherwise non-uniform and unsteady flow arising in high-speed liquid entrainment conditions on the outer wall of a wide rotating drum. We show that the coating flow in this rotary dragout undergoes axial modulations to form an array of roughly vertical thin liquid sheets which slowly drift from the middle of the drum towards its sidewalls. Thus, the number of sheets fluctuates in time such that the most probable rib spacing varies ever so slightly with the speed, and a little less weakly with the viscosity. We propose that these axial patterns are generated due to a primary instability driven by an adverse pressure gradient in the meniscus region of the rotary drag-out flow, similar to the directional Saffman–Taylor instability, as is wellknown for ribbing in film-splitting flows. Rib spacing based on this mechanistic model turns out to be proportional to the capillary length, wherein the scaling factor can be determined based on existing models for film entrainment at both low and large capillary numbers. In addition, we performed direct numerical simulations, which reproduce the experimental phenomenology and the associated wavelength. We further include two numerical cases wherein either the liquid density or the liquid surface tension is quadrupled while keeping all other parameters identical with experiments. The rib spacings of these cases are in agreement with the predictions of our model.

Low-inertia pulsatile flows in highly distensible viscoelastic vessels exist in many biological and engineering systems. However, many existing works focus on inertial pulsatile flows in vessels with small deformations. As such, here we study the dynamics of a viscoelastic tube at large deformation conveying low-Reynolds-number oscillatory flow using a fully coupled fluid–structure interaction computational model. We focus on a detailed study of the effect of wall (solid) viscosity and oscillation frequency on tube deformation, flow rate, phase shift and hysteresis, as well as the underlying flow physics. We find that the general behaviour is dominated by an elastic flow surge during inflation and a squeezing effect during deflation. When increasing the oscillation frequency, the maximum inlet flow rate increases and tube distention decreases, whereas increasing solid viscosity causes both to decrease. As the oscillation frequency approaches either $0$ (quasi-steady inflation cycle) or $\infty$ (steady flow), the behaviours of tubes with different solid viscosities converge. Our results suggest that deformation and flow rate are most affected in the intermediate range of solid viscosity and oscillation frequency. Phase shifts of deformation and flow rate with respect to the imposed pressure are analysed. We predict that the phase shifts vary throughout the oscillation; while the deformation always lags the imposed pressure, the flow rate may either lead or lag depending on the parameter values. As such, the flow rate shows hysteresis behaviour that traces either a clockwise or counterclockwise curve, or a mix of both, in the pressure–flow rate space. This directional change in hysteresis is fully characterised here in the appropriate parameter space. Furthermore, the hysteresis direction is shown to be predicted by the signs of the flow rate phase shifts at the crest and trough of the oscillation. A distinct change in the tube dynamics is also observed at high solid viscosity which leads to global or ‘whole-tube’ motion that is absent in purely elastic tubes.

We investigate the fluid–solid interaction of suspensions of Kolmogorov-size spherical particles moving in homogeneous isotropic turbulence at a microscale Reynolds number of $Re_\lambda \approx 140$. Two volume fractions are considered, $10^{-5}$ and $10^{-3}$, and the solid-to-fluid density ratio is set to $5$ and $100$. We present a comparison between interface-resolved (PR-DNS) and one-way-coupled point-particle (PP-DNS) direct numerical simulations. We find that the modulated energy spectrum shows the classical $-5/3$ Kolmogorov scaling in the inertial range of scales and a $-4$ scaling at smaller scales, with the latter resulting from a balance between the energy injected by the particles and the viscous dissipation, in an otherwise smooth flow. An analysis of the small-scale flow topology shows that the particles mainly favour events with axial strain and vortex compression. The dynamics of the particles and their collective motion studied for PR-DNS are used to assess the validity of the PP-DNS. We find that the PP-DNS predicts fairly well both the Lagrangian and Eulerian statistics of the particle motion for the low-density case, while some discrepancies are observed for the high-density case. Also, the PP-DNS is found to underpredict the level of clustering of the suspension compared with the PR-DNS, with a larger difference for the high-density case.

We study the response of a flexible prism with a square cross-section placed in cross-flow through a series of experiments conducted at increasing flow velocities. We show that as the reduced velocity (a dimensionless flow velocity that also depends on the natural frequency of the structure) is increased, the prism undergoes vortex-induced vibration (VIV) in its first mode, which then transitions to VIV in the second mode and then third mode. In these ranges, the shedding frequency is synchronised with the oscillation frequency, and the oscillations are mainly in the transverse (cross-flow – CF) direction. As we keep increasing the reduced velocity, we observe a linear increase in the amplitude of the torsional oscillations of the prism, resembling a torsional galloping. This increase in the torsional oscillations then causes an increase in the amplitudes of the CF and inline (IL) oscillations while the third structural mode is still excited in the CF direction. A transition to oscillations in the fourth structural mode is observed at higher reduced velocities, which reduces the CF and IL amplitudes, while the torsional oscillations reach a plateau. After this plateau is reached in the torsional oscillations, galloping is observed in the CF oscillations of the response, which results in large-amplitude oscillations in both the CF and IL directions. The CF galloping response at these higher reduced velocities is accompanied by a torsional VIV response and the shedding frequency is synchronised with the frequency of the torsional oscillations.

In this paper, we develop an analytical model to investigate the generation of instability waves triggered by the upstream acoustic forcing near the nozzle lip of a supersonic jet. This represents an important stage, i.e. the jet receptivity, of the screech feedback loop. The upstream acoustic forcing, resulting from the shock-instability interaction (SII), reaches the nozzle lip and excites new shear-layer instability waves. To obtain the newly excited instability wave, we first determine the scattered sound field due to the upstream forcing using the Wiener–Hopf technique, with the kernel function factored using asymptotic expansions and overlapping approximations. Subsequently, the unsteady Kutta condition is imposed at the nozzle lip, enabling the derivation of the dispersion relation for the newly excited instability wave. A linear transfer function between the upstream forcing and the newly excited instability wave is obtained. We calculate the amplitude and phase delay in this receptivity process and examine their variations against the frequency. The analytically obtained phase delay enables us to evaluate the phase condition for jet screech and predict the screech frequency accordingly. The results show improved agreement with the experimental data compared with classical models. It is hoped that this model may help in developing a full screech model.

The bi-stable dynamics of a one-degree-of-freedom disk pendulum swept by a flow and allowed to rotate in the cross-flow direction is investigated experimentally. For increasing flow velocity, a subcritical bifurcation is observed from a Pendulum state, characterised by an increasing time-averaged pendulum angle with large amplitude fluctuations, to a rotating state with a non-zero mean rotation velocity at a critical free stream velocity $U_{P2W}$. The rotating state, referred to as Windmill state, presents a strong hysteresis: once initiated, it is sustained down to velocities $U_{W2P}\lt U_{P2W}$ before bifurcating towards the Pendulum state. A thorough experimental characterisation of the dynamical features of each state is reported, with a particular focus on the influence of the static yaw angle of the disk $\beta _0$ and the free stream velocity. In the Pendulum state, the system behaves differently depending on whether $\beta _0$ lies below or above the stall angle of the disk, with more regular dynamics below. We demonstrate that the bifurcation between the Pendulum state and the Windmill state is triggered by aerodynamic fluctuations, while the bifurcation between the Windmill state and the Pendulum state is deterministic. A stochastic model faithfully reproduces the dynamical features of both states, as well as the characteristics of the bifurcations.

When an oblate droplet translates through a viscous fluid under linear shear, it experiences a lateral lift force whose direction and magnitude are influenced by the Reynolds number, the droplet’s viscosity and its aspect ratio. Using a recently developed sharp interface method, we perform three-dimensional direct numerical simulations to explore the evolution of lift forces on oblate droplets across a broad range of these parameters. Our findings reveal that in the low-but-finite Reynolds number regime, the Saffman mechanism consistently governs the lift force. The lift increases with the droplet’s viscosity, aligning with the analytical solution derived by Legendre & Magnaudet (Phys. Fluids, vol. 9, 1997, p. 3572), and also rises with the droplet’s aspect ratio. We propose a semi-analytical correlation to predict this lift force. In the moderate- to high-Reynolds-number regime, distinct behaviours emerge: the $L\hbox{-}$ and $S\hbox{-}$mechanisms, arising from the vorticity contained in the upstream shear flow and the vorticity produced at the droplet surface, dominate for weakly and highly viscous droplets, respectively. Both mechanisms generate counter-rotating streamwise vortices of opposite signs, leading to observed lift reversals with increasing droplet viscosity. Detailed force decomposition based on vorticity moments indicates that in the $L\hbox{-}$mechanism-dominated regime for weakly to moderately viscous droplets, the streamwise vorticity-induced lift approximates the total lift. Conversely, in the $S\hbox{-}$mechanism-dominated regime, for moderately to highly viscous droplets, the streamwise vorticity-induced lift constitutes only a portion of the total lift, with the asymmetric advection of azimuthal vorticity at the droplet interface contributing additional positive lift to counterbalance the $S\hbox{-}$mechanism’s effects. These insights bridge the understanding between inviscid bubbles and rigid particles, enhancing our comprehension of the lift force experienced by droplets in different flow regimes.

Many problems in elastocapillary fluid mechanics involve the study of elastic structures interacting with thin fluid films in various configurations. In this work, we study the canonical problem of the steady-state configuration of a finite-length pinned and flexible elastic plate lying on the free surface of a thin film of viscous fluid. The film lies on a moving horizontal substrate that drives the flow. The competing effects of elasticity, viscosity, surface tension and fluid pressure are included in a mathematical model consisting of a third-order Landau–Levich equation for the height of the fluid film and a fifth-order Landau–Levich-like beam equation for the height of the plate coupled together by appropriate matching conditions at the downstream end of the plate. The properties of the model are explored numerically and asymptotically in appropriate limits. In particular, we demonstrate the occurrence of boundary-layer effects near the ends of the plate, which are expected to be a generic phenomenon for singularly perturbed elastocapillary problems.

A combination of physics-based and data-driven post-processing techniques is proposed to extract acoustic-related shear-layer perturbation responses directly from spatio-temporally resolved schlieren video. The physics-based component uses momentum potential theory to extract the irrotational (acoustic and thermal) component from density gradients embedded in schlieren pixel intensities. For the unheated shear layer, the method filters acoustic structures and tones not evident in the raw data. The filtered data are then subjected to an efficient data-driven dynamic mode decomposition reduced-order model, which provides the forced acoustic perturbation response for broad parameter ranges. A shear layer comprising Mach 2.461 and 0.175 streams, corresponding to a convective Mach number 0.88 and containing shocks, is adopted for illustration. The overall perturbation response is first obtained using an impulse forcing in the wall-normal direction of the splitter plate, extending into subsonic and supersonic streams. Subsequently, impulse and harmonic forcings are independently applied in a pixel-by-pixel manner for a precise receptivity study. The acoustic response shows a convective wavepacket and acoustic burst from the splitter plate. The interaction with the shock and associated wave dispersion emits a second, slower, acoustic wave. Harmonic forcing indicates higher frequency-dependent sensitivity in the supersonic stream, with the most sensitive location near the outer boundary-layer region, which elicits an order of magnitude larger acoustic response compared with disturbances in the subsonic stream. Some receptive forcing regions do not generate significant acoustic waves, which may guide excitation with low noise impact.

The real-fluid effect induced by large density variation at supercritical pressure (SCP) modulates the turbulent dynamics and heat transfer, and poses challenges to existing turbulence models that are based on ideal-gas conditions. This study conducts direct numerical simulations of fully developed channel flows at SCP, with the upper and lower channel walls being isothermally heated and cooled, respectively. Emphasis is placed on examining the effects of various levels of density variations on near-wall turbulence as well as turbulent heat transfer by changing wall temperatures. The results show that the density fluctuation significantly impacts both first-order and second-order turbulence statistics near the heated wall owing to the close vicinity of pseudo-boiling point. Such real-fluid impact increases substantially with increasing density ratio, and tends to weaken the turbulent kinetic energy by damping turbulence production, while simultaneously inducing an additional turbulent mass flux that partially offsets this reduction. Detailed quadrant analysis reveals that the ‘ejection’ events dominate diverse effects of density fluctuation on Reynolds shear stresses, with density fluctuation contributing positively on the cooled wall side, and negatively on the heated wall side. Regarding the turbulent heat transfer, density fluctuation enhances the enthalpy–pressure–gradient correlation, tending to weaken the turbulent heat flux, which is slightly compensated by additional terms induced by density fluctuations. The overall negative contribution of density fluctuation to turbulent heat flux stems primarily from ‘hot ejection’ motions. Instantaneous flow characteristics provide additional support for these findings. Additionally, the mechanisms by which density fluctuations affect Reynolds shear stress and turbulent heat flux could also be extended to the skin friction coefficient and Nusselt number, respectively.

Many particles, whether passive or active, possess elongated shapes. When these particles settle or swim in shear flows, they often form regions of accumulation and depletion. Additionally, the density contrast between the particles and the fluid can further alter the flow by increasing the local suspension density, resulting in a two-way buoyancy–flow coupling mechanism. This study investigates the buoyancy–flow coupled dispersion of active spheroids, examining the effects of elongation, orientation-dependent settling and gyrotaxis in a vertical pipe subjected to either downwards or upwards discharge. While the concentration and velocity profiles of passive settling spheroids and spherical gyrotactic swimmers can be analysed similarly to a recent study, notable differences in dispersion characteristics emerge due to different streamline-crossing mechanisms. For suspensions of elongated swimmers, the interplay between orientation-dependent settling, gyrotaxis-induced accumulation and shear-induced trapping results in distinct concentration and velocity distributions compared to those of neutrally buoyant particles and extremely dilute suspensions with negligible coupling effect. These differences further impact drift velocity, dispersivity, and the time elapsed to steady dispersion under varying flow rates. Interestingly, low-shear trapping of non-settling elongated swimmers around the centreline, commonly observed in planar Poiseuille flow, is absent in the vertical pipe due to the change of confinement from reflectional to rotational symmetry. However, elongated settling swimmers show a non-trivial concentration response to strong downwelling discharge. This phenomenon, linked to the centreline accumulation of passive settling spheroids, bears similarities to low-shear trapping observed in planar Poiseuille flow.

We study the dynamic deflation of a hydraulic fracture subject to fluid withdrawal through a narrow conduit located at the centre of the fracture. Recent work revealed a self-similar dipole-flow regime, when the influence of material toughness is negligibly small. The focus of the current work is on the influence of material toughness, which leads to an additional self-similar regime of fracture deflation with fixed frontal locations in the toughness-dominated regime. The two limiting regimes can be distinguished by a dimensionless material toughness $\Pi _k$, defined based on a comparison with the influence of the viscous thin film flow within the fracture: $\Pi _k \to 0$ indicates the dipole-flow regime, while $\Pi _k \to \infty$ indicates the fixed-length regime. For intermediate $\Pi _k$, the fracture’s front continues to propagate during an initial period of deflation before it remains pinned at a fixed location thereafter. A regime diagram is then derived, with key scaling behaviours for the frontal dynamics, pressure and volume evolution summarised in a table for the self-similar stage. A comparison is also attempted between theoretical predictions and available experimental observations of viscous backflows from transparent solid gelatins.

We study the dynamics of a thin liquid sheet that flows upwards along the sides of a vertically aligned, impacting plate. Upon impact of the vertical solid plate onto a liquid pool, the liquid film is ejected and subsequently continues to flow over the solid surface while the plate enters the water. With increasing impact velocity, the liquid film is observed to rise up faster and higher. We focus on the time evolution of the liquid film height and the thickness of its upper rim and discuss their dynamics in detail. Similar to findings in previous literature on sheet fragmentation during drop impact, we find the rim thickness to be governed by the local instantaneous capillary number based on gravity and the deceleration of the liquid sheet, showing that the retraction of the rim is primarily due to capillarity. In contrast, for the liquid film height, we demonstrate that the viscous dissipation in the thin boundary layer is the dominant factor for the vertical deceleration of the liquid sheet, by modelling the time evolution of the film height and showing that the influences of capillarity, gravity and deceleration due to the air phase are all negligible compared with the viscous term. Finally, we introduce characteristic viscous time and length scales based on the initial rim thickness and show that the maximum height of the film and the corresponding time can be determined from these viscous scales.

In this study, we conducted interface-capturing high-resolution simulations of a bubbly upflow in a vertical channel to investigate the bubble distribution and its interaction with surrounding turbulence, focusing on the effects of the density ratio. A bulk Reynolds number $Re_b=2300$ was used for all simulations. The influence of density ratio on vortex structures and turbulence statistics differed between the near-wall and core regions of the channel. Adding 5.43 $\%$ gas caused an increase in wall friction. By applying a generalised FIK identity to analyse wall friction, it was determined that the drag rise in the bubbly channel was mostly due to the near-wall region. Visualisation of the bubble and vortex structures showed that small bubbles near the wall induced larger magnitude of Reynolds shear stress and increased wall friction. Bubble behaviour near the wall region was similar for density ratios above 30, leading to wall friction saturation. In the core region, large deformable bubbles created wake vortices due to slip velocity between liquid and gas phases. Wake vortices help large bubbles absorb smaller bubbles and maintain their sizes. As the density ratio increased, the slip velocity increased owing to greater difference in the gravitational acceleration between liquid and gas phases, resulting in corresponding increase in wake intensity and velocity fluctuations. However, quadrant analysis showed that Q1 and Q3 events increased together with Q2 and Q4 events in the core region, cancelling out any net effect of wake vortices on Reynolds shear stress or wall friction.

We use direct numerical simulations to investigate the energy pathways between the velocity and the magnetic fields in a rotating plane layer dynamo driven by Rayleigh–Bénard convection. The kinetic and magnetic energies are divided into mean and turbulent components to study the production, transport and dissipation in large- and small-scale dynamos. This energy balance-based characterisation reveals distinct mechanisms for large- and small-scale magnetic field generation in dynamos, depending on the nature of the velocity field and the conditions imposed at the boundaries. The efficiency of a dynamo in converting the kinetic energy to magnetic energy, apart from the energy redistribution inside the domain, is found to depend on the kinematic and magnetic boundary conditions. In a small-scale dynamo with a turbulent velocity field, the turbulent kinetic energy converts to turbulent magnetic energy via small-scale magnetic field stretching. This term also represents the amplification of the turbulent magnetic energy due to work done by stretching the small-scale magnetic field lines owing to fluctuating velocity gradients. The stretching of the large-scale magnetic field plays a significant role in this energy conversion in a large-scale turbulent dynamo, leading to a broad range of energetic scales in the magnetic field compared with a small-scale dynamo. This large-scale magnetic field stretching becomes the dominant mechanism of magnetic energy generation in a weakly nonlinear dynamo. We also find that, in the weakly nonlinear dynamo, an upscale energy transfer from the small-scale magnetic field to the large-scale magnetic field occurs owing to the presence of a gradient of the mean magnetic field.

Recent experiments and simulations have sparked growing interest in the study of Rayleigh–Bénard convection in very slender cells. One pivotal inquiry arising from this interest is the elucidation of the flow structure within these very slender cells. Here we employ tomographic particle image velocimetry, for the first time, to capture experimentally the full-field three-dimensional and three-component velocity field in a very slender cylindrical cell with aspect ratio $\Gamma =1/10$. The experiments cover a Rayleigh number range $5.0 \times 10^8 \leqslant Ra \leqslant 5.0 \times 10^9$ and Prandtl number 5.7. Our experiments reveal that the flow structure in the $\Gamma =1/10$ cell is neither in the multiple-roll form nor in the simple helical form; instead, the ascending and descending flows can intersect and cross each other, resulting in the crossing events. These crossing events separate the flow into segments; within each segment, the ascending and descending flows ascend or descend side by side vertically or in the twisting manner, and the twisting is not unidirectional, while the segments near the boundary can also be in the form of a donut like structure. By applying the mode decomposition analyses to the measured three-dimensional velocity fields, we identified the crossing events as well as the twisting events for each instantaneous flow field. Statistical analysis of the modes reveals that as $Ra$ increases, the average length of the segments becomes smaller, and the average number of segments increases from 2.5 to 3.9 in the $Ra$ range of our experiments.

The acoustic field radiated by a system of contra-rotating propellers in wetted conditions (with no cavitation) is reconstructed by exploiting the Ffowcs Williams–Hawkings acoustic analogy and a database of instantaneous realizations of the flow. They were generated by high-fidelity computations using a large eddy simulation approach on a cylindrical grid of 4.6 billion points. Results are also compared against the cases of the front and rear propellers working alone. The analysis shows that the importance of the quadrupole component of sound, originating from wake turbulence and instability of the tip vortices, is reinforced, relative to the linear component radiated from the surface of the propeller blades. The sound from the contra-rotating propellers decays at a slower rate for increasing radial distances, compared with the cases of the isolated front and rear propellers, again due to the quadrupole component. The quadrupole sound is often neglected in the analysis of the acoustic signature of marine propellers, by considering the only linear component. In contrast, the results of this study point out that the quadrupole component becomes the leading one in the case of contra-rotating propulsion systems, due to the increased complexity of their wake. This is especially the result of the mutual inductance phenomena between the tip vortices shed by the front and rear propellers of the contra-rotating system.

The hydrodynamic behaviours of finite-size microorganisms in turbulent channel flows are investigated using a direct-forcing fictitious domain method. The classical ‘squirmer’ model, characterized by self-propulsion through tangential surface waves at its boundaries, is employed to mimic the swimming microorganisms. We adopt various simulation parameters, including a friction Reynolds number Reτ = 180, two squirmer volume fractions 𝜑0 = 12.7 % and 2.54 % and a blocking ratio (squirmer radius/half-channel width) κ = 0.125. Results show that pushers (propelled from the rear) induce a more pronounced decrease in the velocity profile than neutral squirmers and pullers (propelled from the front). This hindrance and the induced particle inner stress τpI positively correlate with the quantity of squirmers accumulated in the near-wall region. Notably, the increase in τpI primarily occurs at the expense of diminishing the fluid Reynolds stress τfR. Compared with passive spheres, a low volume fraction (𝜑0 = 2.54 %) of pullers results in a slightly enhanced velocity profile across the channel. In the near-wall region, the swimming direction of the squirmers shows no significant tendency with respect to the flow direction. In the bulk-flow region, pushers and neutral squirmers tend to align their axes more along the flow direction, whereas pullers exhibit a slight preference for alignment with the normal direction.