The used simulation model is based on the open source software package CFDEMcoupling^® (Goniva et al., 2012). This package combines the DEM code LIGGGHTS^® (Kloss et al., 2012) and the CFD code OpenFOAM^® (version 4) for the simulation of multiphase flows. For LIGGGHTS^® and CFDEMcoupling, the academic version published by the department of particulate flow modeling (PFM, JKU Linz, Austria) was used. The transient turbulent flow inside the shake flask was simulated with the volume-of-fluid (VOF) solver “cfdemSolverMultiphase” (Vångö et al., 2018). In short, the VOF approach uses a single momentum and mass balance equation for all phases (Hirt and Nichols, 1981; Brackbill et al., 1992). To track the phase interface, the distribution of the liquid and gas phase fractions is calculated by a transport equation of the phase void fraction (Gueyffier et al., 1999). A Reynolds-averaged Navier-Stokes (RANS) approach was used in combination with a renormalization group (RNG) k-epsilon model for the modeling of the turbulence (Li et al., 2013). The centrifugal forces acting inside a ShF on orbital shaker (Certomat S ll, Sartorius AG, Göttingen, Germany) with an orbital amplitude of 50 mm were considered by imposing a centrifugal acceleration on the fluid phase and the beads. For this purpose, the gravity field of fluid solver was extended in x- and y-direction to include the centrifugal acceleration (Schrader et al., 2019). In the unresolved CFD-DEM coupling, the CFD cells are too large to directly calculate the flow between the particles and the resulting forces. For this reason, the momentum exchange between the fluid and solid phases is realized by force models for drag, pressure gradient, surface tension, and viscous forces (Bérard et al., 2020). In detail, the Di Felice model was used for the drag force (Di Felice, 1994). The void fraction field was calculated using the divided model approach of the of the CFDEM framework and smoothed using a diffusion based equation (Pirker et al., 2011). On the DEM side, the Hertz-Mindlin model was used to calculate contact forces in the normal and tangential direction (Tsuji et al., 1992; Di Renzo and Di Maio, 2004; Di Renzo and Di Maio, 2005). The constant directional torque (cdt) model was applied as rolling friction model (Ai et al., 2011).

In the previous simulations (Schrader et al., 2019), the ShF geometry had to be slightly extended in both horizontal planes to ensure that the beads were always inside the CFD mesh. As a result, some beads, especially near the curved sidewall, were no longer inside the first boundary cells. Therefore, the open-source meshing tool cfMesh (Juretić, 2014) was used to generate a new mesh without artificially increasing the dimensions. This CFD mesh was able to describe the curved flask geometry more accurately while maintaining cell dimensions similar to the first mesh. The velocity boundary condition was set to noSlip at the wall of the shaking flask and pressureInletOutletVelocity at the opening of the shaking flask. The pressure boundary condition was set to zeroGradient at the wall and totalPressure at the opening. For a more accurate description of the wetting behavior of the liquid on the shake flask wall, the contact angles listed in Table 1 were used instead of the 90° contact angle used in the past.

Previously, the hydrodynamic interactions between an approaching and separating pair of beads near the contact were considered by applying a lubrication force model. To account for sliding forces and moments as well as torsional moments in addition to normal forces, the corresponding model equations of Kroupa et al. (2016) were implemented in the DEM code. Due to the divergence of the lubrication forces and moments at an infinitesimally small separation distance, a limit was set at a separation distance (slip length) of 10% of the bead radius. Additional simulation parameters are shown in Table 1.

The coefficients of friction for particle-particle and particle-wall contacts were assumed to be identical because the wall and beads are made of glass. Ceramic and glass contacts were also assumed to have the same friction coefficient as ceramic-ceramic contacts. In addition to the fluid properties for the cultivation medium, the values for water are given, because water was used for validation.

Based on the macroparticle-enhanced cultivations of Schrader et al. (2019) and Schrinner et al. (2021) several simulation series (see Table 2) were performed with different combinations of shaking frequency f , bead diameter d b , bead density ρ b and bead volume concentration c v , b . In the first series, different bead sizes (0.540, 0.658, 0.969, 1.183, 1.513, 1.746, and 1.932 mm) were tested at different shaking frequencies and constant bead concentration. Since the number of beads N b changes when the bead concentration is constant and the bead diameter varies, different bead sizes and a constant number of 4,200 beads were used in another study. Glass and ceramic bead mass concentrations were varied to study the effect of bead number at constant bead size. To improve the comparability of the bead concentrations, the bead volume concentrations are also given in Table 2.

The stress energy (SE) according to the analytical model for stirred media mills (Kwade, 2003) is defined as S E ∝ S E b = d b 3 ∙ v t 2 ∙ ρ b (1)

In this study, the maximum achievable SEs for different contact types were calculated analogously to the energy equations described by Beinert et al. (2015). A distinction is made between normal, shear, torsional and rolling bead-bead or bead-wall contacts. Depending on the type of contact, the corresponding kinetic energy is calculated using the translational or rotational relative velocity instead of the stirrer tip velocity v t . These kinetic energies correspond to the maximum possible SE in the cultivation process. The SF is defined as the number of bead-bead or bead-wall contacts per second. Instead of a constant time interval for the evaluation of the bead contacts, a fixed number of six shaker rotations was chosen. This avoids potential over- or under-evaluation of the temporally fluctuating SE when the shaking frequency is changed. Within the six shaker rotations, all new bead contacts that occurred in a time step were saved. A cumulative distribution of the SE was then calculated from the list of all contacts of for six rotations. Therefore, the distribution is not a time average. Instead, it represents the SE spectrum within the time period under consideration. The translational shear velocity of the beads was calculated from the translational shear energy E t s based on Eq. 2 (Beinert et al., 2015). E t s = 1 2 ∙ π 6 ∙ d b 3 ∙ v b , t s 2 ∙ ρ b (2)

The translational bead velocities, used in the validation (see chapter 3.5), were read out at 0.01 s intervals for a period of five shaker rotations and saved in a list. This list was then used to calculate a cumulative distribution of bead velocity based on the evaluated DEM data for five shaker rotations.

The mean volumetric power consumption P v in the cultivation medium, which also depends on the shaking frequency, was calculated from the turbulent energy dissipation rate ϵ and the medium density ρ m according to the following equation (Li et al., 2013): P v = ϵ ∙ ρ m (3)

Data were analyzed every 0.01 s to calculate the time averaged power consumption for the last five shaker rotations.

To track the complex motion of the liquid and the beads from below the ShF, an experimental setup was constructed (Supplementary Figure S1). An aluminum plate with four vertical threaded rods was screwed onto the base plate of the orbital shaker. Afterwards, a 250 mL ShF with four baffles, already filled with 50 mL of water and beads, was clamped between two transparent polymer plates at a height of about 40 cm. To increase the stability of the experimental setup, two 3D-printed brackets and guy ropes (for 160 min⁻¹ only) were installed. A high-speed camera (Promon 501, AOS Technologies AG, Switzerland) mounted vertically below the center of the ShF recorded images (896 × 920 pixels) at a frame rate of 200 s⁻¹. Optimal illumination of the ShF was crucial to minimize disturbing image artifacts caused by light refraction and reflection on the water surface. On the one hand, a square 30 cm LED panel was placed statically above the orbital shaker. On the other hand, a circular LED strip was installed at the level of the bottom of the ShF. A cardboard cylinder, covered on the inside with aluminum foil, was used for homogeneous illumination. Glass (1.183, 1.513, and 1.932 mm) and ceramic (0.918 mm) beads (Sigmund Lindner GmbH, Warmensteinach, Germany) were used in the experiments at shaking frequencies of 100, 120, and 160 min⁻¹. In contrast to the simulations, the minimum bead size was limited to approximately 1 mm, as smaller beads could not be reliably detected in the image analysis. The transparent glass beads were difficult to identify in the images during pre-processing. Therefore, a mixture of 26% black coated marker beads and 74% transparent beads of the same size was chosen. In the case of white ceramic beads, it was not possible to use beads of a different color, so the amount of beads was reduced from 40 mL L⁻¹ (glass) to 2.63 mL L⁻¹ to avoid overlapping beads in the images.

Prior to evaluation, the images were pre-processed using MATLAB^® to separate the beads from the background structures of the liquid-gas phase interface. In the case of glass beads, the imported images were converted to grayscale and smoothed with a flat-field correction. In the second step, a binary image was generated using an adaptive threshold. Then, the centers of circular objects within a defined radius range were detected in order to draw a new image with only the detected beads as white circles on a black background. In the case of white ceramic beads, an adapted image pre-processing was necessary, since no black marker beads could be used. In addition to image optimization and filtering steps, a marker-controlled watershed separation was applied to separate the partially overlapping beads. In addition, the centers of the beads were mapped to generate a new image with non-overlapping circles for PTV analysis. The actual calculation of bead velocities out of the pre-processed images was performed using the open-source MATLAB^® tool PTVlab (version 1.0.0.0) (Brevis et al., 2011; Patalano and Wernher, 2013). The number of images used to calculate the velocity distribution was the equivalent number of images taken during five rotations of the shaker. Therefore, the distribution is not a time average. Instead, it represents the velocity spectrum within the time period under consideration. In detail, the Gaussian mask option was selected for particle detection and the cross-correlation method for particle tracking. By default, PTVlab provides an interpolated velocity field whose expansion is determined by the spatial distribution of the beads. As a result, bead velocities that do not exist were assumed in zones without beads. For this reason, the PTVlab source code was extended to output the discrete velocity distribution without spatial interpolation.

Shaking frequency is the driving force for the beads movement. The bead-wall shear velocity for all bead diameters increases steadily with increasing shaking frequency (Figure 2). It is noticeable that the shear velocity initially increases gradually more strongly up to a shaking frequency of 140 min⁻¹ and then the slope decreases. Overall, larger beads reach a higher velocities than smaller ones. However, this is not true for the two largest bead diameters. In the high shaking frequency range, the shear velocities were slightly below the values for a bead size of 1.513 mm.

Forces affecting the bead motion pattern must be considered to explain the observed trends in ShF. For simplicity, the beads are assumed to move with velocity v b on a circular path with radius r r o t (Eq. 4). This radius is a function of the shaking frequency, because the centrifugal forces F a , caused by the centrifugal acceleration a c , push the beads outward (Eq. 5). The ShF radius limits the theoretical maximum rotation radius. However, gravity counteracts the outward movement of the bead as soon as a bead enters the curved area of the ShF bottom. This balance of forces, which is also influenced by frictional, fluid-particle interaction and particle contact forces, determines the maximum bead rotation radius. If the centrifugal influence is considered in Eq. 4, this results in a theoretical increase in bead velocity with the cube of the shaker frequency in Eq. 6. Besides, larger beads experience higher centrifugal forces, which can additionally affect the bead motion. Overall, this explains the initial large increase in velocity at the lower shaking frequencies in Figure 2. v b = r r o t F a c ⋅ 2 π ⋅ f (4) r r o t ∝ F a ∝ a c ⋅ m b ∝ r o r b ⋅ f 2 ⋅ m b (5) v b ∝ r o r b ⋅ f 3 ⋅ m b (6)

As the shaking frequency increases, the centrifugal effect gradually decreases because most of the beads are centrifuged to the outer zone of the ShF. As a result, it follows that the bead velocity increases linearly with the shaking frequency.

In Figure 3A, the bead diameter was varied from 0.540 to 1.932 mm at a constant bead volume concentration of 40 mL L⁻¹ to determine the dependence of bead-wall shear velocity on the bead diameter. For the shaking frequency of 100 min⁻¹, only a small linear increase in velocity is observed. Whereas at higher shaking frequencies, the shear velocity increases progressively more strongly up to a bead size of 1.183 mm. For 120 min⁻¹ and 140 min⁻¹, the shear velocity starts to increase only slightly from a bead diameter of 1.183 mm, while a decrease occurs at even higher shaking frequencies and larger bead diameters. The increase in velocity can be related to the influence of centrifugal forces, which depend on bead mass and shaking frequency (see Eq. 4). In contrast, other bead size dependent forces (e.g., drag, pressure gradient, lubrication, friction forces) increase with increasing diameter. As a result, additional forces dampen the increase in shear velocity with the bead diameter.

Furthermore, for a given bead volume concentration, the number of beads increases with decreasing diameter. As a result, more beads are distributed on the bottom of the flask and new layers of beads are formed on top of the bottom layer (Supplementary Figure S3 and Supplementary Video S3). Consequently, more bead-bead and bead-baffle interactions make a free bead motion more difficult. In addition, new shear planes are created between the bead layers, so that the shear velocity between the beads increases as the number of beads increases (Supplementary Figure S2A).

At higher shaking frequencies, the beads are forced intensively against the curved bottom of the ShF. On the one hand, the beads largely accumulate in layers with increasing rotational speed, resulting in an increased shear velocity between the beads (Supplementary Figure S2A). On the other hand, the normal forces and thus the friction forces on the beads increase. As a result, the beads begin to roll faster instead of sliding over the glass surface of the ShF. This reduces the translational shear velocity portion and conversely increases the translational rolling velocity portion. To illustrate the effect, the average velocity of the beads is shown in Figure 3B. The spatially averaged bead velocity was taken directly from the x- and y-velocity components of the DEM simulations. For example, for 200 min⁻¹ it can be seen that the spatially averaged bead velocity stagnates from a bead size of 1.2 mm. In contrast, the shear velocity decreases (Figure 3A). Consequently, if the spatially averaged bead velocity is constant, the rotational velocity of the beads must be increasing. The rotational acceleration of a sphere is inversely proportional to its rotational inertia and rolling friction, according to the equation of rotational motion. Due to the proportionality of the rotational inertia to the mass and size of a bead, larger beads change their rotational speed more slowly than smaller or lighter beads because of the applied moments. Furthermore, tangential forces on a bead cause moments that accelerate or decelerate the rotation depending on the direction of the tangential force. The tangential forces depend on the normal bead forces and the coefficient of friction between the beads. In addition, beads in layers hamper each other’s rotation. In sum, the intensity of translational shear is determined by the rolling behavior of a bead. The complex interplay of influencing factors illustrates that a direct prediction of shear stress by beads is extremely complex or even impossible without simulations.

When comparing the absolute velocities (Figure 3B), it is noticeable that the bead velocities are higher than the shear velocities. Therefore, it should be noted that rolling of the beads may be hindered by the biomass in the cultivation, so that in reality deviating shear velocities can occur.

In addition to this effect, in Figure 3A, at a shaking frequency of 120 min⁻¹ and a constant bead number the shear velocity decreases above 1.183 mm. In this case, the underlying cause cannot only be increased bead rolling, as the average bead velocity also decreases in Figure 3B. In all other simulations with constant bead mass, the number of beads decreased with increasing diameter, so that the reduced number of beads compensated for the increased space requirement per bead. Instead, with a constant number of beads, the ShF bottom fills up more at large bead diameters, making bead motion more difficult. As a result, the velocity of the beads in Figure 4 decreases above a volume concentration of approximately 73 mL L⁻¹.

In addition to shaking frequency and bead size, the bead volume concentration was varied at a constant shaking frequency of 120 min⁻¹ (Figure 4). For a constant glass bead size of 0.969 mm, the bead-wall shear velocity first increases slightly from 0.12 to almost 0.13 m s⁻¹ and then decreases with increasing bead volume concentration. However, in contrast to glass beads, an increase in the number of ceramic beads led to a significantly higher increase in shear velocity. Accordingly, at a bead volume concentration of around 40 mL L⁻¹, approximately 43% higher shear velocities were achieved between the beads and the wall of the ShF. Interestingly, the velocity of the ceramic beads decreased only slightly with increasing concentration.

The simulations illustrate (Supplementary Figures S4, S5 and Supplementary Videos S4, S5) that ceramic and glass beads formed layers of beads with increasing bead concentration. Compared to the ceramic beads, the glass beads spread over a larger area of the ShF bottom than the ceramic beads. The higher velocities and greater proportion of bead layers for the ceramic beads can be explained by their 52% higher density, which results in increased centrifugal forces. In addition, the coefficients of friction and rolling friction for ceramics are lower than for glass beads, which in turn can lead to a variation in the rolling and shearing velocity proportions.

In the study by Schrader et al. (2019), only the liquid level at the ShF was considered for the initial validation of the simulations. In the present study, experimental images of the beads at the bottom of the ShF were taken to verify the bead velocities calculated by the CFD-DEM simulations. Figure 8 shows the simulated and experimentally determined images for the largest glass bead diameter of almost 2 mm for three shaking frequencies between 100 and 160 min^-1 used during the cultivations (note that corresponding videos to the images can be found in the Supplementary Materials). Next to the almost transparent glass beads are the black marker glass beads used to obtain the bead velocities using PTV. At the lowest shaking frequency of 100 min⁻¹, the beads are distributed over a large area of the ShF bottom, but centrifugal forces pushed the beads further out as the speed increased. At a shaking frequency of 160 min⁻¹, the beads were in a tightly packed bed in the curved part of the ShF bottom. This results in greater interaction between the baffles and the beads, further increasing the stresses. Moreover, the formation of multiple layers of beads made it difficult to detect of individual marker beads. Light reflections from the surface created dark artifacts that made image analysis even more complex. By using multiple light sources and a cardboard cylinder with aluminum foil surrounding the ShF, the light artifacts were reduced to a minimum.

The bottom and side views of the simulations show the glass beads and the fluid in the interface, colored according to their velocity. As in the experiments, the beads in the simulation are pushed further towards the outer diameter of the ShF as the shaking frequency increases. This significantly increases the velocity of both the beads and the liquid (note that the colors change with shaking frequency). In addition, the side view shows that the liquid is increasingly forced against the side walls. While at 100 min⁻¹ there is almost exclusively a monolayer of beads on the bottom, at higher velocities, layers of beads are formed as already seen in the experiments. At a shaking frequency of 160 min⁻¹, some beads are deflected at the baffles and moved upwards so that the beads are mixed more intensively. In summery, the simulations show a good qualitative agreement with the experiments. However, the beads are more densely packed in the experiments than in the simulations, especially at higher shaking frequencies. Nevertheless, the CFD-DEM simulations are able to describe the motion behavior of the beads in the ShF with sufficient accuracy.

In addition the qualitative comparison in section 3.5.1, the velocities determined by PTV are compared with the results of simulations. Therefore, Figures 9A–C shows the cumulative bead velocities determined experimentally and numerically for three bead diameters and shaking frequencies. In general, analogous to the discussion in section 3.2.2, the bead velocity is higher than the shear velocity. This is because the overall bead velocity includes the translational rolling velocity in addition to the translational shear velocity. At the lowest shaking frequency of 100 min⁻¹, the experimental distributions for different bead diameters are similar, indicating that bead size has little influence on bead velocity. At higher shaking frequencies, the distributions differ depending on bead size, especially in the lower velocity range. At the same time, the maximum velocities achieved are almost identical. A higher bead mass with larger bead diameters results in stronger centrifugal forces, so that more beads are transported to the edge of the ShF. Furthermore, a constant mass of beads was used, which means that the number of beads decreases with bead diameter to the power of three as the bead size increases. In addition, the area required per bead increases as the projection area increases with bead diameter to the power of two. Assuming a single layer of beads, the required area of all beads is the product of the number of beads and the projection area. Consequently, the area requirement increases inversely with decreasing bead diameter for a constant total bead mass concentration. This, in turn, results in more beads moving at lower velocity in the central region of the ShF for smaller bead sizes, and finally in a broader velocity distribution. In reality, smaller beads form vertical layers, which partially reduces the effect discussed before.

Overall, there is good agreement between the experimental and numerical distributions for the two lowest shaking frequencies of 100 and 120 min⁻¹. However, the experimental distributions for 100 min⁻¹ flatten out in the upper velocity range. It is noticeable that at a shaking frequency of 120 min⁻¹, the two smallest bead sizes have a bimodal velocity distribution, while this effect is somewhat more pronounced for smaller beads. A likely cause is that the simulation data includes velocities from all vertical bead layers. In contrast, PTV captures only a two-dimensional velocity field without the frequency of bead velocities in the vertical direction. Accordingly, this method cannot capture the frequency of bead velocities in different bead layers. Smaller bead diameters tend to have more layers of beads, so the effect described above becomes more pronounced. Somewhat larger deviations between experiment and simulation are seen at the highest speed. The main reason for the deviations might be that the beads are now arranged very densely. This increases the likelihood of incorrect mappings of particle pairs between two images in the PTV analysis. Figure 8 shows that the beads at 160 min⁻¹ are slightly more compressed in the experiment than in the simulation, so this could be another reason for the observed differences.

In addition to the glass beads, white ceramic beads were also investigated. The experimental and simulated mean and maximum velocities are close. However, the experimental distribution shows slightly lower velocities in the lower velocity range. Overall, it is shown that the presented experimental method is suitable to determine the velocities of the beads in the ShF. Moreover, the comparison between experiment and simulation shows a high degree of agreement. Accordingly, the chosen simulation approach is a suitable tool for the analysis and prediction of the stresses caused by the beads.

However, biomass was not considered in the simulations and experiments. The cultivation broth becomes turbid to non-transparent after a certain cultivation time. Thus, reliable identification and tracking of the beads on the camera images becomes much more difficult due to poorer illumination conditions. Depending on the morphological state of the culture and the biomass concentration, the viscosity of the culture broth increases and the viscosity behavior often becomes non-Newtonian (Olsvik and Kristiansen, 1994). The use of transparent model fluids may be a feasible approach to achieve realistic viscosity behavior while maintaining optical accessibility (Bliatsiou et al., 2020). It is expected that the general motion behavior of the beads is not significantly affected by collisions with pellets due to the higher inertia of the beads. However, biomass entrapment between the beads may alter (rolling) friction. This, in turn, would affect the intensity of bead stress.