Custom 3D printed ridged geometries were fabricated on a Stratasys Objet30 Scholar (Stratasys, Ltd.) and a Formlabs 2 SLA (Formlabs, Inc.) printer to fit the Anton Paar rheometer shaft for disposable geometries.

Patterning with 3D printing allows to create well-defined geometries with good control to produce fine features with a dimensional accuracy of ∼100 μm. The benefit of using the 3D printed ridged geometries is that they can be easily incorporated into a flow imaging technique such as MRI [44, 45], which allows to fully quantify the flow profiles and extract the wall stresses present within the system.

The geometries with ridges perpendicular to the rotational axis were designed by making cutouts of various spacings (0.5–3 mm) and depths (0.5–3 mm) from the standard CC17 geometry with constant R _i = 8.5 mm and R _o = 9 mm. They were printed on the Stratasys Objet30 Scholar, using the Vero Black photopolymer (Stratasys, Ltd.) with water-soluble supports. The geometries were printed vertically with a layer thickness of 16 μm.

The geometries with ridges at an angle with respect to the rotational axis were designed with the angle ranging from 10° to 60° with 10° intervals. Geometries with ridge angles between 70° and 90° could not be printed due to limitations of the printing technique. The geometries were designed such that the shortest distance (spacing) between the ridges was 2 mm, with a ridge depth of 1 mm and ridge thickness of 0.3 mm. This was achieved by changing the width and the number of cutouts. The width and the number of cutouts were calculated using a Python script with the length and the diameter of the ridged geometry, ridge angle and shortest distance between the ridges as input. These geometries with angled ridges were fabricated on a Formlabs 2 SLA printer, using the glass reinforced Rigid resin (Formlabs, Inc.) with a layer height of 50 μm and xy-resolution of 140 μm. The support needed for the printing was only attached to the shaft and the top chamber to prevent the introduction of artefacts coming from small residuals of the support material.

The printed geometries were washed twice for 10–15 min in 95% isopropanol. After drying, the geometries were post-cured under a 366 nm lamp for 1.5–2 h. The geometries were rotated 180° halfway through to cure evenly. The support material was removed manually using flush side-clippers after curing. The shaft of all geometries was printed with a diameter of 8 mm and milled down on a lathe to 6.95 mm to fit the Anton Paar shaft for disposable geometries. We determined that this lathe post-processing ensured the best concentricity of the geometry with the Anton Paar rheometer shaft and cup.

Rheological experiments were performed on Anton Paar MCR 300 and MCR 501 rheometers, using the custom geometries in standard cups and cup holders (Anton Paar CC17). The experiments were performed at 20.5°C, using either a Peltier heat exchange element and waterbath at 20°C as a heat sink (for MCR 501) or a high flowrate waterbath at 20.5°C (for MCR 300). All the printed geometries were tested on a Newtonian fluid. We used castor oil as a Newtonian fluid, due to its high viscosity (∼1 Pa s at 20°C compared to 50 mPa s for most oils) and its stability against oxidation and high temperatures for an extended period of time. Therefore, its viscosity as a function of time at a given temperature can be considered constant. Unlike other viscous fluids such as glycerol, castor oil is non-hygroscopic and does not evaporate at room temperature, does not degrade over time, and does not swell nor degrade the 3D printed geometries, which sets it apart from most organic solvents.

Before each measurement, the geometry was placed in a container with the test fluid to eliminate air bubbles between the ridges. The geometries were then visually inspected for the presence of air bubbles, which, if present, were removed by rubbing the geometry against the cup while it was submerged in the fluid.

We performed steady shear measurements with average shear rates ranging from γ ̇ = 10^–3 to 10² s⁻¹, resulting in a range of Re from 2.5 × 10⁻⁵ to 2.5 for castor oil (only the laminar regime was probed). In these measurements, the average shear rate γ ̇ is defined as the surface velocity of the moving inner cylinder divided by the gap size of the CC17 geometry. Each measurement consisted of 26 datapoints, distributed logarithmically, and each datapoint started with a stepwise increase of γ ̇ , followed by an equilibration period of 2 min, after which the data was averaged for 5 s. To test the viability of using 3D printed geometries in rheological measurements, a solid CC geometry was designed with an outer diameter as close as possible to that of a commercially available smooth standard stainless steel CC geometry. These geometries were checked by comparing the flow curves of the test fluid. Both the 3D printed and the standard geometries yielded identical flow curves indicating that the accuracy of the 3D printer was high enough to produce complex rheological geometries and no deviations were expected in the results due to the intrinsic roughness of the 3D printed geometries.

Finite Element calculations were done using COMSOL Multiphysics^® [46] to complement the experiments. Performing numerical calculations in addition to experiments was necessary given the limitations inherent in the printing technique, as it was not feasible to fabricate all the different types of custom geometries and to experimentally vary all the geometric parameters. To simplify the calculations, the ridged geometries were approximated as infinitely large parallel plate geometries. A depiction of the computational domain is shown in Supplementary Figure S1 in Supplementary Material. We also performed mesh independence tests, see e.g. results in Supplementary Table S1. We performed simulations for varying S, D, and θ, while keeping the gap size R _gap,0 (i.e., the distance between the outer wall and the tip of the ridges) and ridge width W constant and equal to the experimentally used values. We then impose a constant velocity difference Δv between the inner and outer wall, and solve the Stokes equation (hence, ignoring inertial effects) numerically using quadratic quadrilateral elements for the velocity components and linear elements for the pressure. Since the flat plate limit does not account for the finite curvature of the actual CC geometry, we performed additional simulations for curved geometries with horizontal and vertical ridges to quantify the effects of curvature (see Supplementary Figures S2, S3). These simulations show that for the horizontal ridges, the effects of curvature are negligible for the used gap size, while curvature effects for the vertical ridges become significant at ridge depths D > 1 mm, and at larger ridge spacings, as depicted in Supplementary Figure S3. We expect that curvature effect for the angled ridges can be interpolated between these two limiting cases. In the current study, all the angled ridged geometries have a ridge depth D = 1 mm, which is sufficiently small to neglect curvature effects.

When a Newtonian fluid is sheared between two surfaces at a certain average shear rate γ ̇ , a shear stress σ is generated and it can be defined as σ = η γ ̇ , with η being the viscosity of the fluid. The fluid flow behaviour can be modified using different boundary conditions, imposed by employing different geometries in a rheometer. Concentric cylinder (CC) is one of the most commonly used geometries, where the fluid is sheared in the gap between the inner and outer cylinder of radius R _i and R _o , respectively. The inner cylinder is made to rotate at a rotational velocity Ω _i in revolutions per second (rps). The average shear rate γ ̇ , in s⁻¹, is defined as the surface velocity of the moving inner cylinder divided by the gap size R _gap. It should be noted that the shear rate is not constant within the gap.

Adhering to Feynman’s approach [47–49] for deriving the final expression for torque in a concentric cylinder configuration, the torque M acting across the cylindrical surface is given by Eq. 1, where ω _i is the angular velocity in rad s⁻¹. M = 4 π η ω i L R o R i 2 R o 2 − R i 2 where ω i = 2 π Ω i (1) The shear force in the azimuthal direction F _shear induced by the fluid can be expressed as the ratio of M and the radius of the inner cylinder R _i . Thus, the shear stress in the azimuthal direction σ can be expressed as F _shear divided by the surface area of the inner cylinder A _i (Eq. 2). σ = F s h e a r A i = M 2 π R i 2 L where F s h e a r = M R i , A i = 2 π R i L (2)

When roughness is introduced into the CC geometry in the form of ridges (as shown in Figure 1B), both the gap size and shear rate within the gap are not constant. Additionally, the flow profiles penetrate in the space between the ridges. We now make the step to quantify this flow profile penetration not via the flow field, but via the effect it has on the effective stress on the wall. We denote the extent of this penetration by a parameter called the penetration depth δ, similar to the protrusion height as mentioned in previous studies [19, 29, 34, 35, 50, 51]. Typically, the penetration depth is determined by inducing fluid flow over a flat surface with ridges on it. In their theoretical investigation, Bechert and coworkers [19] simulated velocity distributions for various ridge configurations and found that the apparent origin of the velocity profile in between the ridges lies in the gap between the inner flat surface and the tip of the ridge. Lee and Lee [51] experimentally investigated the flow structures inside the semi-circular ridges. They referred to the distance between this origin and the ridge tip as the protrusion height. Another equivalent way to determine penetration depth is to create ridges on a concentric cylinder geometry and introduce continuous flow by rotating the ridged geometry with respect to the fluid. We use the second method to generate different flow patterns. Due to the introduction of the ridged geometries, the effective gap size R _gap changes as a function of the ridge spacing S and the ridge depth D, as denoted in Eq. 3. R g a p S , D = R o − R i − δ (3) where R _i is assumed to be identical to the radius of the inner cylinder without the cutouts. Due to the penetration of the flow in between the ridges, the effective gap size increases. The ridged geometry can thus be considered similar to the standard solid CC geometry with an effective inner radius R _i,eff = R _i − δ (a similar equivalence was reported by Luchini et al. [34] for a flat corrugated wall).

For every geometry with a particular ridge spacing S and ridge depth D, the penetration depth δ can be measured from standard rheological experiments. Rewriting the torque-rotation rate conversion equation (Eq. 1), we get M = 4 π η ω i L R o R i , e ff 2 R o 2 − R i , e ff 2 where R i , e ff = R i − δ (4) In this equation, R _i , R _o , L are known from the geometry and Ω _i , M were measured from the rheological experiments. η was determined by fitting the flow curves in a non-ridged (solid) CC geometry with σ = η γ ̇ . In this way, the penetration depth for every rotation rate can be obtained by solving Eq. 4. For Newtonian fluids, the penetration depth does not depend on the rotation rate Ω _i , therefore, we averaged the penetration depths corresponding to ten higher rotation rates to give a single penetration depth value for a geometry with a specific ridge spacing and ridge depth.

From the simulations in the flat plate limit, the penetration depth was calculated from the average shear stress on the outer wall using σ = η γ ̇ = η Δ v / ( R g a p , 0 + δ ) .

We measured the flow curves from the rotational tests on multiple horizontal ridged geometries. On a log-log scale, the torque M increased linearly with the rotation rate Ω for over four orders of magnitude in rotation rate (Figure 3A), indicating that the presence of ridges did not change the Newtonian shear rate dependence of the flow. To investigate the variations due to the different roughness parameters, we looked more closely on a linear scale (Figure 3B) and found that the measured torque decreases with increasing ridge depths D for a given ridge spacing S. This is due to the influence of the flow profile penetration on the measured torque, i.e., the penetration depth δ increases with ridge depth, leading to an increase in the effective gap size (per Eq. 3) and a decrease in the measured torque. It should be noted that in all these experiments, the rotation rates were low enough to avoid Taylor vortices.

Figure 3C shows the penetration depth, obtained from the measured torques using Eq. 4, as a function of D for various ridge spacings S. We see that δ first increases with D, but then levels off to a plateau that depends on the value of S. Hence, the flow can penetrate between the ridges only up to a certain maximum depth, δ _∞ , which increases approximately linearly with the ridge spacing S (Figure 3D).

To complement the experiments and to gain more insight into the flow profile penetration, Finite Element simulations were performed. The geometries were recreated as described in Section 2.3, and the curvature of the geometry was neglected. For this geometry, the only non-zero velocity component is in the azimuthal (ϕ) direction, while all gradients in the ϕ-direction are zero. The pressure was also constant in this case. The simulated penetration depths as a function of different ridge depths and ridge spacings are plotted in Figure 3C as dashed lines. The experimental and simulation results match quite well, especially at larger D. Similar to the experiments, the δ values plateau after a certain ridge depth; and the plateau value δ _∞ exhibits a linear increase with S, see Figure 3D.

For a solid cylinder case (S = 0), δ _∞ must be zero as the fluid cannot penetrate into a solid cylinder, and therefore, we would normally expect δ _∞ to scale linearly with the ridge spacing and the linear plot (dashed line in Figure 3D) to extrapolate to the origin. However, this is not observed and δ _∞ seems to deviate from the linear dependence on S at small ridge spacings below 0.5 mm. This is because in this small S limit, the ridge spacing becomes comparable to the ridge width, which also starts to play a role in determining the penetration depth. Strikingly, the normalised penetration depth δ/δ _∞ for different ratios of D/S collapses onto a master curve (Figure 3E). The collapse to the master curve is observed to be satisfactory for the numerical data, whereas it is poor in case of the experimental data, particularly in the small D limit. The underlying reasons behind this observation are discussed in the subsequent paragraph. This master curve can be described by a simple exponential function: δ δ ∞ = 1 − e − k D / S with k ≈ 5.68 (5) The initial part of this master curve, corresponding to small D or large S, can be linearized to give δ/δ _∞ ≈ kD/S. It is expected that δ → D in this regime, so that we expect k = 1/k ₂, with k ₂ being the slope of the curve of δ _∞ versus S. We find 1/k ₂ ≈ 5.88 from Figure 3D, which is indeed close to the value found for k. This can also be understood as follows: a static boundary layer is created due to the presence of the small ridges. The thickness of this static boundary layer is D − δ and it goes to zero as the spacing between the ridges becomes extremely large. The prefactor k takes into account the other geometric effects, most likely the effect from the other length scale, i.e., the ridge width W. At a constant ridge spacing S, the thicker the ridges, the lower will be the flow penetration, hence, δ will be smaller. We expect that as the ridge width increases, δ will be small and hence, δ _∞ will also drop. If the normalised penetration depth δ/δ _∞ decreases, then k increases. But if both δ and δ _∞ decrease by the same factor, then δ/δ _∞ will not change, hence, k will not change. Therefore, k can be written as a function of normalized W as: k = f (W/S).

The deviations of the experimental results from that of simulations are more pronounced at smaller D, and this is a consistent deviation for all ridge spacings. In order to understand these deviations, we explored different possible sources of error. First, we tried to reduce the end effects from the top and bottom conical parts of the geometry. Their contribution to the torque was calibrated on a solid cylindrical geometry and the torques for the ridged geometries were corrected for these end effects. However, deviations were still observed, suggesting the presence of additional sources of error. For instance, small inaccuracies during the fabrication of the geometries have a bigger influence at smaller cutouts. Likewise, the measurement of the actual dimensions of the 3D printed geometries such as the inner and the outer radii, as well as the length of the geometries is challenging due to the small size of the geometries. Additionally, boundary conditions like wall slip are uncontrolled in the experiments, especially for large S and small D, which may further contribute to small errors. Another possibility is that small air bubbles may be trapped between the ridges, especially for small ridge spacings, and can lead to occasional deviations in the results. Finally, we want to highlight the additional physics that can come into play at ridge spacing below 0.5 mm (small S limit), where the ridge spacing becomes comparable to the ridge width. The ridge width is another important length scale that is relevant in case of ridges with small spacings. However, the analysis of the effect of different ridge width is beyond the scope of the present work. Due to these reasons, we took the plateau value of delta (δ _∞ ) in experimental results as the average of the last three δ values. By linking the experiments and the simulations, it is possible to determine how changing the D and S of the ridges changes the penetration depth δ.

The second limiting case for patterned walls is a cylindrical geometry with vertical ridges. Here, we show only the numerical results in the limit where the curvature can be neglected (i.e., for large cylinder radii), see Section 2.3. In this case, the fluid velocity has components both in the ϕ and the r-directions, and gradients in the z-direction are zero.

The penetration depth for vertical ridges with different ridge depths and ridge spacings were calculated as described in Section 2.4 and plotted in Figure 4A. The flow profiles penetrate more deeply in between the ridges with increasing ridge depth and spacing, therefore, δ increases. This behaviour is similar to what was observed in the case of horizontal ridges, however, the penetration depth values in vertical ridges are smaller than for the horizontal ridge case because the main flow direction is perpendicular to the ridges while for the horizontal ridges, the flow is along the ridges. In contrast to the horizontal ridge case, we see two (quasi-)plateau values of δ as a function of D: the first plateau value is at around D = 0.5 mm, where the ridge depth is close to the ridge width and most likely, the plateau value is attributed to the influence of the ridge width. The second plateau value is at large ridge depths similar to what was observed in the case of horizontal ridges. We consider the second plateau value to be δ _∞ .

δ _∞ increases linearly with the ridge spacing S (Figure 4B), but the slope is about 2.5 times smaller in comparison to the horizontal ridge case. Additionally, the linear plot of δ _∞ against S does not precisely extrapolate to the origin and there seems to be non-linear dependence at small ridge spacings below 0.5 mm.

The normalised penetration depth δ/δ _∞ for different ratios of D/S collapses onto a master curve (Figure 4C). However, a single exponential fit does not work so well compared to what was seen in the horizontal ridge case, suggesting a two-step function fit. This is most likely due to the presence of two plateaus in the δ vs. D variation (Figure 4A). For S ≫ D, i.e., for small values of D/S, we expect δ → D; the plot of δ/δ _∞ against D/S is linear in this region and the slope of this plot is inversely related to the slope of the variation of δ _∞ with S.

Finally, it is worth noting that the relationship between δ _∞ and S could potentially depend on the ridge width W and the gap size R _o − R _i . These are a few additional length scales to explore, which may give different regimes (for example: two plateau values in δ vs. D plot). The current plots shown are for S ≫ W. This is a substantially general rescaling that to our knowledge has not been reported in the literature, and it is helpful to get a better understanding of the angled ridge case, as discussed below.

Having studied the two limiting cases, we now investigate the influence of geometries with ridges at finite angles on the flow behaviour of a Newtonian fluid. All the angled geometries have the same ridge spacing and ridge depth, and the ridge angles vary from 0° to 60° with 10° intervals. For different angles, the torque required to sustain a certain Ω shows a linear dependence on the rotation rate even at higher rotation rates, indicating that the angled ridges do not affect the Newtonian shear rate dependence of the flow (Figure 5A). However, the slope of the torque versus rotation rate changes as a function of ridge angle (Figure 5B), indicating that the extent of flow profile penetration into the gap between the ridges changes with the angle. To examine the effect of different angles, we calculated the penetration depths from the M vs. Ω relationship as previously done for the limiting cases.

Interestingly, δ shows a monotonic dependence on the ridge angle (Figure 5C) and increases as θ increases, with the largest δ for the horizontal ridges (θ = 90°). There are a few outliers at smaller ridge angles. One of the plausible reasons is δ is derived from the assumption that there is a shear stress only in the azimuthal direction. In reality, there is a normal force (and hence, a shear stress in the axial direction) and it is not clear whether the normal force is completely decoupled from the azimuthal shear stress.

We compare our experimental findings with numerical results. In the case of angled ridges, the velocity has non-zero components in the ϕ, r, and z-directions, and also the velocity gradients in all three directions are non-zero, so, a full 3D calculation is necessary. As we restrict ourselves to a flat geometry (corresponding to large cylinder radii), we can make use of symmetries to limit the computational domain. In particular, we model a domain with five ridges (see Supplementary Figure S1A) and impose a zero-outflow condition at the bottom and top of the domain, and periodic boundary conditions along the azimuthal directions. To reduce the influence of end effects, we exclude the top and bottom regions when calculating the average shear stress on the wall.

As shown in Figure 5C, the simulated δ values for angled ridges in the flat plate regime exhibit a similar trend as the experimental values, with δ monotonically increasing as the angle increases from 0° (vertical) to 90° (horizontal). However, the experimental values are systematically higher than the simulated ones, with an offset of ∼ 0.1 mm. It is not clear what determines this offset; potential factors include partial wall slip (which would reduce the measured torque) or curvature effects. Although we have previously seen that the effect of curvature is negligible in horizontal and vertical ridges, it is highly unlikely that the curvature will have an effect in the case of angled ridges (refer to Section 2.3 for details). Likewise, the presence of even a small amount of wall slip suppresses the shear rate because in such a slip situation, there would be a large effective gap size which lowers the torque. A lower torque would then result in higher penetration depth, as observed in the case of experiments. It is important to note that the flows in two directions (azimuthal and axial) are independent of each other; the superimposition of an axial flow along the symmetry axis of the rotating inner cylinder would not change the expression for the torque M. Such axial flow components cannot be the source of the discrepancy either.