To simulate the flow past different Drusinae species, we use realistic body shapes and resolve the problem numerically in time and space using OpenFOAM. To that end, we employed micro-computer tomography (µ-CT) data of different species to construct the body shapes. Based on the larvae’s natural setting, a rectangular computational domain in the form of an open channel was created. The reconstructed bodies were then placed in the channel and exposed to a turbulent flow, based on flow velocities measured in situ at the larvae locations. The workflow is illustrated in Fig. 1 where the mean flow is in positive x direction (from left to right). We use a cartesian coordinate system which is aligned with the channel walls. The origin of the coordinate system is placed on the bottom wall exactly beneath the upstream-facing tip of the larva’s head (Fig. 1). While the channel height is fixed to H = 25 mm, the elevation of the body varies within 2.20 mm ≤ h ≤ 3.30 mm, depending on the respective larva considered.

The turbulent flow is fully described by the incompressible Navier-Stokes equations and suitable boundary conditions as indicated in Fig. 1. On the surface of the body, as well as on the bottom wall, no-slip boundary conditions are imposed. The top boundary is assumed to be stress-free, emulating a free surface. For computational economy, the channel’s side boundaries are considered periodic. The fluid enters the domain in a turbulent state through the inlet with a bulk Reynolds number

Reb=UbHν (1)

where U_b is the bulk (or mean) velocity and ν is the kinematic viscosity (ν = 1.31 × 10⁻⁶ m²/s for water at a temperature of 10.5 °C). At the outlet, we assume zero streamwise gradients of all flow velocities. The numerical solution yields the instantaneous velocity u→(x→,t) and pressure p(x→,t) fields. After a transient time, these quantities are averaged to obtain temporal mean values.

For the flow measurements, nine final instar larvae from Drusinae were selected in their characteristic habitats, representing the three clades of the subfamily (two shredders, four grazers, three filter feeders). The measurements were taken directly in front of the larvae using a tripod mounted Schiltknecht MiniWater 20 Micro velocimeter probe (resolution: 0.01 m/s, probe diameter H_p = 11 mm) with a sampling rate of 1 Hz. The measured velocities ranged from 0 to 0.93 m/s. The total water depths were between 10 mm and 30 mm. Typically, larvae were sitting on top of flattened sediment pebbles, aligned with the flow direction and head upstream. Further information and an analysis of the measured data can be found in Waringer et al. (2021). The present simulations were carried out for the two velocities out of the above range, U_b = 0.55 m/s and U_b = 0.40 m/s.

We employed micro-computer tomography (µ-CT) to create the external surface of five specimens according to the three different feeding modes prevailing in Drusinae: Drusus alpinus Meyer-Dür, 1875 (representing the shredder clade), D. bosnicus Klapálek, 1899 and D. monticola McLachlan, 1876 (representing the grazer clade) and D. discolor (Rambur, 1842) and C. nebulicola McLachlan, 1867 (representing the filtering carnivore clade). For µ-CT scans, the merged volume was exported as *.TXM file into Amira 2019.1 (FEI SAS, Mérignac, France, part of Thermo Fisher Scientific^TM). Image segmentation was achieved in Amira 6.5.0 (Visage Imaging, Inc., San Diego, CA, USA). Then we used the Amira Surface Generate tool to create the three-dimensional surface renderings. Further information regarding the technical aspects of tomographic scans and a discussion of the three-dimensional tomography reconstruction of the external surface can be found in Zittra et al. (2022).

The surface renderings were further processed using Blender to prepare 3D geometries for the numerical simulations. In Blender, the original finely triangularised geometry was smoothed and transformed into a symmetric quadrilateral surface mesh. One example of the re-topology is shown in Fig. 2 (head only).

After re-topology of the torso, the legs were moved to a position which resembles the clinging position of the larvae observed in the field. An example is shown in Fig. 3. During the whole processing of the µ-CT data, we tried to preserve the original scale of the animal.

Based on U_b and H from above, the simulation is carried out for bulk Reynolds numbers Re_b = 7634 and Re_b = 10496. Since the flow is expected to be turbulent at these Reynolds numbers, we use Large Eddy Simulation (LES). The turbulent inlet flow is generated by the divergence-free synthetic eddy method (DFSEM) (Poletto et al. 2013). DFSEM needs an auxiliary Reynolds-Averaged-Navier-Stokes simulation (RANS) of the channel flow without the body beforehand. The procedure is illustrated symbolically in Fig. 4 and the set-up follows that of Vieira et al. (2023) for the flow past a wall-mounted cuboid.

The simulations are carried out using OpenFOAM. This open-source software employs cell-centred finite volumes to discretise the governing equations. The mesh was generated with the tool snappyHexMesh, which requires a cubic hexahedral background mesh (without the body). This background mesh was created with cells of size ∆x⁺ = ∆y⁺ = ∆z⁺ = ∆⁺_max = 11.32, where the superscript ‘+’ indicates non-dimensional quantities (in wall units). They are based on the friction velocity uτ=τw¯/ρ , where τw¯ is the mean wall stress of the turbulent flow (here acquired from the auxiliary RANS) and ρ the fluid density. Hence, the non-dimensional velocity, length and time are u+→=u→/uτ, x+→=x→/(ν/uτ), and t⁺ = t / (H / u_τ), respectively. The pressure is scaled differently (see below). In a final step, the mesh was refined locally near the channel floor and the larva’s surface reaching ∆y⁺_min = 0.36. The maximum number of grid points was 6.49 × 10⁶ for D. monticola and the minimum was 5.03 × 10⁶ for C. nebulicola. A non-dimensional time of t⁺ = 0.44 was adequate for all field variables to reach a statistically steady state. After that, mean quantities were calculated by averaging over t⁺ ∈ [0.44,1.75].

To analyse the flow past the larva, it is useful to define a local velocity

UL=11.4HHL∫0HL∫-0.7H0.7Hu¯x0,y,zdydz, (3)

to which the larva is exposed. The velocity U_L is defined as the mean over the area [0,H_L] × [−0.7H,0.7H] at the position x = x₀ upstream of the larva. Here, we use the height H_L = 10 mm (0.4H) which is comparable to the height H_P = 11 mm of the velocity probe used to measure the flow velocity in front of the larvae, and x₀ = −H_L = 0.4H. Note that the mean elevation of the larvae from the ground is h = 2.86 mm which is only a fraction of H_L (Table 2).

Fig. 5 shows the mean velocity profiles u¯x0=-HL,y,z=0 in the symmetry plane z = 0 upstream of the larva for Re_b = 10496. The horizontal dashed line marks the average height h = 2.86 mm (0.286H) of the larva models. From Fig. 5 and Table 2, it is seen that the larvae are mainly found in the buffer region y⁺ ∈ [5,30] of the turbulent flow, such that both viscous shear stress and turbulent shear stress are significant to the results.

Fig. 6 shows contours of the streamwise component of the mean velocity u¯/UL near the larvae for Re_b = 10496 and z = 0. Generally, the streamwise velocity distribution does not vary much amongst the larvae. In the dark blue regions, the streamwise velocity is negative (u¯<0), indicating the existence of mean recirculation regions. A recirculation region upstream of the larva is a common feature amongst all species and is commonly associated with a horseshoe vortex. Additionally, we recognise a mean recirculation region between the head and the downstream part of the body, which is especially pronounced for D. bosnicus. Additional mean recirculation regions are depicted between the posterior area of the body and the bottom wall next to the downstream flow region. In this case, they are more pronounced for D. bosnicus and D. discolor than for D. alpinus.

For a 3D visualisation, we show in Fig. 7 the pressure distribution on the surface of the larvae for Re_b = 10496. A common property of all pressure distributions is the high pressure on the forehead of the larvae. In addition, hypersurfaces are shown in cyan on which u¯=0, marking back-flow regions. The extended cyan regions in Fig. 7 suggest a mean recirculating flow exists between the posterior ventral part of the body and the bottom wall for D. bosnicus, D. discolor and C. monticola (Fig. 7B, C, D). The reason seems to be related to the body shapes and postures. The back-flow region downstream of the bodies resembles a wake vortex, similar as for other bluff bodies. The wake vortex is not detected for D. nebulicola (Fig. 7E), probably because it is significantly smaller than the other species (Table 2) and, thus, is facing a smaller mean velocity over its body height. The body tapering off very gently also tends to suppress a wake vortex. Other smaller recirculation regions arise in body indentations and behind sharp edges, for example, the characteristic pronounced pronotum of D. bosnicus.

An important question is whether different species perceive different mean forces due to the flow. In the symmetric arrangement considered, only drag Fx¯ and lift forces (Fy¯) arise. They are usually expressed in terms of the drag and lift coefficients

Cd¯=Fx¯(ρ/2)UL2A, (4)

Cl¯=Fy¯(ρ/2)UL2A, (5)

which are just the scaled forces. As the scale, we employ the pressure rise (ρ / 2)U²_L in the forward stagnation point of a body in a homogeneous flow with velocity U_L multiplied by the streamwise projected area of the body (Table 3). Since the forces arise due to the pressure, as well as the viscous stress distributions on the body, we decompose the drag and lift coefficients into a pressure and a viscous part,

Cd¯=Cd.p¯+Cd,ν¯, (6)

Cl¯=Cl,p¯+Cl,ν¯, (7)

which are distinguished by the subscripts p (pressure part) and ν (viscous part). Note that the viscous normal stress is taken care of by the pressure part.

Fig. 8 shows the mean drag and lift coefficients for all larvae and Reynolds numbers simulated. The pressure and viscous contributions are shown in black and blue, respectively. For Re_b = 7634 the maximum overall drag is found in C. nebulicola with a difference of 41% compared to D. alpinus (which has the lowest drag value) for the same Reynolds number (Fig. 8A). In the case of Re_b = 10496, the difference between species decreases to ≈ 28% for the drag coefficient, with the lowest value present in D. alpinus and the highest in D. monticola. The lift coefficients are about 3.5 times smaller than the drag coefficients (Fig. 8B). The highest values are found consistently for D. bosnicus for the two Reynolds numbers tested. Additionally, the minimum lift values are found consistently for C. nebulicola.

Fig. 9 shows the Reynolds stress u'u'¯/uL2, describing the squared magnitude of the fluctuation velocity in the x direction. The Reynolds stress level for D. discolor (and for C. nebulicola, not shown) is found to be higher than for the other species. In channel flow, the production of turbulent kinetic energy has a maximum near y⁺ = 11 (Laadhari 2002; Monkewitz 2022). Accordingly, the streamwise velocity fluctuations, measured by u'u'¯, also reach their maximum in this region (Lee and Moser 2015). Additionally, the characteristic shape of their bodies, the concave face of D. discolor and the high number of edges of C. nebulicola might induce flow separation and recirculation, as well as shedding of vortices, which can increase the stresses.

The smallest overall Reynolds stresses u'u'¯/uL2 are found for Drusus bosnicus. A reason could be Drusus bosnicus has the smallest length-to-height ratio amongst the species. Furthermore, the tapered shape of the body may create a streamlined silhouette, which may be responsible for stress reduction. Regarding D. alpinus, the highest stresses arise past the dorsal line.

Fig. 10 shows the distribution of the streamwise component of the mean skin-friction coefficient Cf¯=2τw,x/ρUb2 for Re_b = 10496, where τw,x=τw→·ex→ is the x component of the vectorial wall shear stress. We noticed that the mean viscous streamwise surface stress is most pronounced near the edges of the bodies, i.e. at the edge of the head and the pronotum. These are regions from which vortices can be shed into the flow.