12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

ACTIVE TURBULENCE AND ITS MODELING

Benjamin Deußen
Martin Oberlack
Chair of Fluid Dynamics
Chair of Fluid Dynamics
Technical Univerity of Darmstadt
Technical Univerity of Darmstadt
Otto-Berndt-Str. 2, D-64287 Darmstadt, Germany Otto-Berndt-Str. 2, D-64287 Darmstadt, Germany
deussen@fdy.tu-darmstadt.deoberlack@fdy.tu-darmstadt.deYongqi Wang
Chair of Fluid Dynamics
Technical Univerity of Darmstadt
Otto-Berndt-Str. 2, D-64287 Darmstadt, Germany
wang@fdy.tu-darmstadt.de

ABSTRACT

the suspension (Baskaran & Marchetti, 2009), the shape of
the particles (Wensink et al., 2012; Be’er et al., 2020) and theenclosure of the entire system (Wu et al., 2017; Dhar et al.,2020). For high particle densities the flow pattern shows similarities to hydrodynamic turbulence, which leads to the occasionally used term ”active turbulence” (Bratanov et al., 2015;
Alert et al., 2022). While both phenomena - active and hydrodynamic turbulence - show similarities, the differences are evident. Hydrodynamic turbulence only occurs in systems with avery large Reynolds number, whereas active suspensions withmicroscopic particles exist in a regime of very low Reynolds
numbers Re ≪ 1. Hence, hydrodynamic turbulent flows aregoverned by inertia and active suspensions by viscosity.

Two dimensional active suspensions are examined. The
suspension consists of a Newtonian fluid and multiple identical active Janus particles. Janus particles are a model for biological or artificial particles capable of driving themselves forward. The length of such a particle lies in the range of a few
micrometers, hence, the Reynolds number Re ≪ 1 is small.
Due to the chaotic behaviour of the particles in suspensions,some authors have used the term active turbulence to describe
this behaviour.
The similarities and differences between active and hydrodynamic turbulence are examined. We will use statistical
moment equations and direct numerical simulations (DNS) to
obtain results. The results are quantified with the help of statistical measures, such as the probability density function (PDF)
and the mean values of the physical quantities.

GOVERNING EQUATIONS
A general, incompressible, Newtonian fluid is describedby the Navier-Stokes equation

INTRODUCTIONF F

Active suspension consists of a fluid, which is usually but
not necessarily Newtonian, and a large number N of active particles. The latter term describes objects, which are capable ofdriving themselves forward (Bechinger et al., 2016). In general, there is no size restriction, i.e. flocks of birds or swarming fish can be described as active particles as well as bacteria or artificial microscopic objects, however, our focus lies
on microscopic active particles. The simplest realisation ofan artificial active particle is the Janus particle (Sharan et al.,2020), which is driven forward by the chemical potential created by two different surface coatings. More complex structures are described as microrobots (Zhang et al., 2019; Yigit
et al., 2020). Possible fields of application are environmental monitoring and cleaning (Soler & Sánchez, 2014; JuradoSánchez et al., 2015; Orozco et al., 2016) and self-healing
electronics (Li et al., 2015).
Colonies of motile strains of Escherichia coli or Bacillus
subtilis are often used as model organisms to examine the behaviour of active particles (Wu & Libchaber, 2000; Drescher
et al., 2011; Sokolov & Aranson, 2012; Dombrowski et al.,2004). Important parameters for the characteristics of the collective behaviour of an active suspension are the density of∂ ui u j∂ uFi∂ pF
1 ∂ 2 uFi
+
=−
+
,
∂t
∂xj
∂ xi
Re ∂ x2j(1)

where uFi is the fluid velocity and pF the fluid pressure. Due tothe low Reynolds number the convective term ∂ uFi uFj /∂ x j issmall compared to the viscous term, hence, it can be neglected,
leading to the unsteady Stokes equation
∂ uFi
∂ pF
1 ∂ 2 uFi
=−
.
+
∂t
∂ xi
Re ∂ x2j(2)

Applying the rotation operator to (2) delivers the equation for
the vorticity ω F∂ ωiF
1 ∂ 2 ωiF=.
∂t
Re ∂ x2j(3)


The vortex stretching term ω Fj ∂ uFi /∂ x j does not occur in(3) as it is a result of the convective term in (1). Subsequently,

1

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
must be fulfilled. Furthermore, a mechanical stress is applied,which permits a slip velocity in the tangential direction. The
tangential vector ti is defined as

eP

n



ti = εi jk ε jlm ePl nm nk ,twhere ePi is the orientation vector of the particle, see Fig. 1.
Due to its definition the tangential vector points always in the
backwards direction ePi ti ≤ 0 of the particle. Hence, the active
stress is

Γ pa

b
a

fiac = Ac ti

Γacthe energy cascade of three dimensional hydrodynamic turbulence will not occur in an low-Reynolds active suspension regardless of the spatial dimension. Furthermore, as the injectionof energy into the system is on the same scale as the dissipation, which dominates all scales of an active suspension, noenergy cascade is expected to occur (Alert et al., 2022).
The motion of the particles is described by the NewtonEuler equations for rigid particles, where we will restrict ourselves to two dimensions

m

=

Z

τiFj n j ds,

dt
ΓP
Z


P
dω
ΘP
=
ε3i j xi − χiP τ Fjk nk ds.
dt
ΓP

To get a first impression of the behaviour of an active suspension we use volume-averaged equations. The underlying
idea of Eulerian spatial averaging theory, see e.g. Drew (1983),
is to express all physical variables in terms of their averages
Z
DE
1ψ P (xx,t) =ψ P (rr ,t) γ P (rr ,t) dV,
V V

(4)

∂ uFj
∂ uFi
+
∂xj
∂ xi

γ (xx,t) =!
(6)∀xx ∈ Γac .

0 x∈/ V P,

α S + α F = 1,

(12)

(13)

S

∀xx ∈ Γ pa .

where α S = VV is the volume fraction of the particle phase and(7)F

α F = VV is the volume fraction of the fluid phase. Using therespective volume fractions, we introduce the average over thevolume occupied by a single phase ψ P = ψ P /α P , which
we call true physical average.
Using the volume average (11) to derive balance equations for the mass, linear and angular momentum for both

The surface of the particle is impermeable, thus, at the active
part of the surface Γac the conditionuFi ni = uPi ni(
1 x ∈ VP

where V P is the volume occupied by the phase P.
Due to the average each phase is present at each spacetime point (xx,t), hence it is necessary to solve the equations
for both phases simultaneously. Furthermore, we have the saturation conditionis the source of momentum for the the particles.
The equations (4) and (5) are valid for active and passive
particles. The capability of self-driven motion is introduced
into our model by the means of an active boundary condition.
The particle surface is split into two parts by the shorter half
axis b (Fig. 1), i.e. our model is a so called Janus particle.
At the passive part Γ pa of the boundary a no-slip condition isapplied. Hence, the local particle and fluid velocity are equaluFi = uPi = vPi + εi3 j ω P x j − χ Pj

(11)

where P indicates an arbitrary phase, ψ P an arbitrary physical variable of the given phase and V the respective averaging
volume. Note that V needs to be sufficiently larger than theparticle volume but also much smaller than the entire computational domain. The function γ P is the phase indicator function
which is defined as

(5)

P

1
Re

(10)

Eulerian spatial averaging theory

Here, vPi is the translational velocity of the particle P, ω P is therotational velocity, χiP is the position of the centre of mass, mP
is the mass and ΘP the moment of inertia and ni is the outward
pointing normal vector at the particle surface ΓP . The fluid
stress tensor
τiFj = −pF δi j +∀xx ∈ Γac ,where Ac is the constant non-dimensional active stress magnitude. Subsequently, the fluid at the active boundary is accelerated in the backwards direction. Due to momentum conservation the resulting reaction force of the fluid will accelerate theparticle in the direction of ePi by the means of the integral termin (4).

Figure 1. A two dimensional Janus particle. The surface ofthe particle is split by the shorter axis b into a passive Γ pa andan active part Γ pa .

P
P dvi

(9)

(8)

2

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
(b)

(a)
x2 1

(c)

1
oben rechts

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0.5

0.6

0.7

0
αS

1

0

0.5

1

0
v1S

0

0.5

1

uF
1

Figure 2. Results obtained for a pressure and active stress driven channel flow. (a) Profiles of the the solid phase volume fraction
α S . (b) Averaged particle translational x1 -velocity vS1 . (c) Averaged fluid x1 -velocity uF1 . On the vertical axis the remaining spatial), Ac = 500 () and Ac = 1000 ( ). Figure
coordinate vertical to the channel is denoted. Results are plotted for Ac = 0 (reproduce from Deußen et al. (2022), Fig. 8a-c.
(d)

(e)

x2 an1 equation for the averaged orientation of1the parphases andticle phase delivers in two dimensions(f)

Results

1

The model is used to obtain results for a planar
0.8
Poiseuille channel flow.
0.8
0.8To obtain results MATLAB (Version
9.2.0.556344, R2017a) is used. The equations are solved us0.6
0.6
0.6 equations solver pdepe.m, employing the partial differential∂ α S ∂ α S vSi∂ α F ∂ α F uFi
ing
finitedifferences(Skeel & Berzins, 1990) for the spa+
= 0,
+
=0
(14)
∂t 0.4 ∂ xi
∂t
∂ xi
0.4
tial discretisation and 0.4
numerical differentiation formulas for
F uF uF
the time discretisation, implemented by ode15s.m (Shampine
F
F pF
∂
α
∂ α F u0.2
∂
α
i j
i
0.2
0.2
+
=−
& Reichelt, 1997; Shampine
et al., 1999). The respective
∂t
∂xj
∂ xi
MATLAB-scripts are publicly available at TUdatalib (Deußen00
0 were, together with further setups,
F
The results∂ uFj
0
1 −1∂ S  ∂ u−0.5S−20 et al., 2021).
20
0.2 0.4 0.6 0.8
1Si + ⟨F−2i ⟩w=
1 0,
ω3Sin Deußen et al. (2022).
e
α
+
(15)
+ F
first
presented
Re ∂ x j∂xj∂ xi ⋅10
The results presented in Fig. 2 show a quasi steady state,
where the averaged vertical velocities are negligibly small∂ α S vSi vSj 2 ∂ α S K∂ α S vSi
compared to the velocities parallel to the channel wall and are
+
+
D
∂t
∂xj
3 ∂ xinot displayed here. The flow is driven by a pressure gradient∂ pF /∂ x1 = −106 and the active force of the particles. Three∂ vSj∂  S S  ∂ vSi
numerical experiments are carried out with varying averaged − ⟨Fi ⟩ = 0,
α µ
+
(16)
−
∂xj
∂ x(a)∂ xi(c)jc = 1000.
(b) active forces Ac = 0, Ac = 500 and A
In case of the passive
suspension
Ac = 0 the parabolic ve" x2 1
1
1
S S
∂ α S ω S vSi
2 ∂ α SW ilocity profile of an Newtonian fluid is retrieved. Due to the low
S ∂α ω
Θ++
Reynolds number the∂t∂ xi
3 ∂ xi0.8 particles follow the stream lines of the
0.8
0.8
!#
fluid. Negligible forces in the vertical direction are observed,∂∂ ωS
hence, the volume 0.6
fraction α S does not change in space, see
0.6− 0.6 α S µ S− ⟨T ⟩ = 0,
(17)
∂ xi
∂ xi
Fig. 2a.
Despite the expectedchaotic behaviour clear structures
0.4
0.4
0.4
S
S
S
∂ α ei v j∂ α S eSi
∂ α S Ac aReM
emerge in results for the averaged quantities of the active sus+−
∂t 0.2∂xj∂ xi 96πλ 20.2 in Fig. 2a material accumulates at thepensions. As visible
0.2



channelwalls
foractivesuspensions. The reason for this acS∂ eS ∂ e j∂
0 For once, active particles have an in − α S εi3k ω S e0S = 0.
cumulation are twofold.
− 0 α S µ S  i +
−5
0
k 0
0.5
1
1.5
∂ x j 0.5
∂ x0.6∂ xi 0.7 αS
w1S et al.,jherenttendencyuF
1 to accumulate at solid walls (Bechinger⋅10−3
2016; Berke et al., 2008). This is represented by the term K(18)
in the equation for the averaged linear momentum (16) of theparticles, see Deußen et al. (2022) for a more detailed discussion and derivation. The second reason for the accumulationis the averaged Saffman-force, acting on the particles
perpenHere, uFi and pF are the mean fluid velocity and pressure, vSi ,
S = uF − vS and rotation
S
S
diculartotherelative
velocityw
ω and ei are the mean translational and rotational particle veiii
locity and the particle orientation. Additional terms arise dueSFωi = 0.5ωi − ωi between the two phases (Deußen et al.,to the interactions between the phases; ⟨Fi ⟩ represents inter2022)
action forces and ⟨T ⟩ interaction torques. The two quantities
K and W represent interactions between the particles and ber
DE1atween particles and thewall based on the averaged active forceLSSqω |− 2 εi jk wSj ωk ,|ωF=
1.615α(19)iM
2aAc R πRe222π sin(σ )A =ε sin (σ ) + cos (σ )dσ . An in depth
c

εV P

2

derivation of the given equation can be found in Deußen et al.
(2022).

The averaged relative velocity wSi is caused by the active force,3

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022

Figure 3. Snapshots of a DNS with 211 active particles and a single passive particle. The left figure shows the situation at t = 20,
the right figure shows the same domain at t = 100. The particles form clusters cluster alongside the wall and to a lesser extent in theproximity of the passive disk-shaped particle.

fdthus, this effect is not visible in a passive suspension. In fact
there are several studies (Leighton & Acrivos, 1985; Derksen & Larsen, 2011; Agudo et al., 2017), which show that the
Saffman force is negligible for the case of passive spheres.
The profiles for the velocity of the two phases in case of
the active suspension show an area of reduced velocity gradients at the channel walls. This is caused by the accumulationof material. In the limiting case, where α S reaches the maximal packing fraction, the suspension will behave like a solid,
leading to a vanishing velocity gradient. The overall velocityof the active suspension is reduced, compared to the passive
suspension. The main reason is the additional stress between
the two phases. This effect may be somewhat resembling the
behaviour of a turbulent flow, where increased stresses provide
a flatter velocity profile. Obviously, however, the cause for thesimilar behaviour is different. In the case of hydrodynamic turbulence, one has to deal with increased internal stresses, whichbecome visible as Reynolds stresses in the averaged NavierStokes equations. In contrast, the active stresses are a material
parameter of the active suspension.

4

2

0

0

0.5

1

1.5

2
d min

Figure 4. The PDF f d for the minimal distance between two
particles is marked by (). The lower limit for the PDF
is the particle dimension plus the minimal distance threshold,marked by ( ). In case that all particles would be equallydistributed in the domain, the distance between each particle
would be d = 1.24 ().

DIRECT NUMERICAL SIMULATIONS
The model system consisting of the unsteady Stokes equation (2) and the Newton-Euler equations (4) and (5) is implemented into the Bounded Support Spectral Solver (BoSSS)
framework 1 (Kummer et al., 2021). This framework provides
high-order discontinuous Galerkin (DG) methods to solve thefluid phase of the suspension. Within this framework it is possible to represent the particle surface with sub-cell accuracy byemploying a level-set method, implemented by the extend discontinuous Galerkin method (XDG) (Kummer, 2017; Kummer
et al., 2018).
We present results for a closed system of 211 active particles and a single passive particle. The length of the active
particles is aac = 1, the aspect ratio is ε = 5/2 and the active
stress magnitude Ac = 10. The passive particle is a disk with adiameter a pa = 2. The simulation is carried out in a 18×18 domain with the mesh size h = 1/12. Fig. 3 shows the simulation
progress between t = 20 and t = 100. From an approximatly

uniform distribution at the first time the particle form clusters
at the second time. These clusters mainly form at the domain
walls and are, thus, similar to the results of the volume averaged simulation. However, they also appear within the domain
and in proximity to the passive particle.
We analyse the behaviour of the particles using probability density functions (PDF), obtained by using the MATLAB
function fitdist. As a measure of the clustering phenomenon, a
PDF for the minimal distance f d between two particle centresof mass is introduced, see Fig. 4. This PDF returns the probability of finding a second particle at a certain distance froma probing particle. It clearly shows that the highest probability is in the closest proximity of the probing particle. Hence,the visual impression of cluster formation obtained from thesnapshots in Fig. 3 can be confirmed by the PDF. In case ofa hypothetical uniform distribution, the distance between theparticles would overall be larger, as indicated by the dashedline in Fig. 4. Hence, clear structures are formed in contrast to

1 Openly available at https://github.com/FDYdarmstadt/BoSSS

4

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
f x 0.1

DISCUSSION
Structures are made visible with both models. Active particles have a tendency to attach to walls, as is visible in thecase of both, the volume-averaged model and the DNS. This
attraction is due solely to the active stress at the surface, as it
does not occur in the case of passive suspension. Furthermore,loose clusters are also formed in the domain, as can be seenfrom the snapshots of the DNS. In this sense, active suspensions thus behave in a predictable and less chaotic manner.
The observation of the velocity PDF shows a nonGaussian behaviour. This shows similarities with the nonGaussianity in the case of hydrodynamic turbulence and possibly indicates the phenomenon of intermittency. Further analysis in this regard is needed. Intermittency, i.e. the sudden
occurrence of a chaotic state in an otherwise ordered system,
thus indicates a behaviour similar to hydrodynamic turbulence.
It should be noted, however, that in an active system there isno classical energy cascade, as can be easily shown by considering the vorticity equation. Considering the energy cascade
as defining property of hydrodynamic turbulence, the term active turbulence might be misleading from the standpoint of
fluid mechanics. Nevertheless, applying techniques known
from hydrodynamic turbulence research to active suspension
might lead to interesting new results as shown for example byDeußen et al. (2021).

8 · 10−2
6 · 10−2
4 · 10−2
2 · 10−2
0

−5

0

x1

5

Figure 5. The PDF f x for the x1 -coordinate of the positionof the particles within the domain. The results are preliminary.
With a longer run-time of the simulation and thus additionaldata, it is to be expected that f x is constant within the domain and only two peaks form on the walls, which are alreadyrecognisable in the available smaller data set.

fv
10

REFERENCES
Agudo, JR, Illigmann, C, Luzi, G, Laukart, A, Delgado, A& Wierschem, A 2017 Shear-induced incipient motion of asingle sphere on uniform substrates at low particle reynolds
numbers. J. Fluid Mech. 825, 284–314.
Alert, R., Casademunt, J. & Joanny, J-F 2022 Active turbulence. Annual Review of Condensed Matter Physics 13,
143–170.
Baskaran, A. & Marchetti, M. C. 2009 Statistical mechanics and hydrodynamics of bacterial suspensions. Proc. Natl.
Acad. Sci. USA 106 (37), 15567–15572.
Bechinger, C., Di Leonardo, R., Löwen, H., Reichhardt, C.,
Volpe, G. & Volpe, G. 2016 Active particles in complex
and crowded environments. Rev. Mod. Phys. 88, 045006.
Be’er, A., Ilkanaiv, B., Gross, R., Kearns, D. B., Heidenreich,S., Bär, M. & Ariel, G. 2020 A phase diagram for bacterial
swarming. Commun. Phys. 3 (1), 66.
Berke, A. P., Turner, L., Berg, H. C. & Lauga, E. 2008 Hydrodynamic attraction of swimming microorganisms by surfaces. Phys. Rev. Lett. 101 (3), 038102.
Bratanov, V., Jenko, F. & Frey, E. 2015 New class of turbulence in active fluids. Proc. Natl. Acad. Sci. USA 112 (49),15048–15053.
Derksen, JJ & Larsen, RA 2011 Drag and lift forces on random assemblies of wall-attached spheres in low-reynoldsnumber shear flow. J. Fluid Mech. 673, 548–573.
Deußen, B., Oberlack, M. & Wang, Y. 2021 Probability theoryof active suspensions. Phys. Fluids 33 (6), 061902.
Deußen, Benjamin, Wang, Yongqi & Oberlack, Martin 2022
A deterministic two-phase model for an active suspension
with non-spherical active particles using the eulerian spatial
averaging theory. Physics of Fluids 34 (2), 023302.
Deußen, B., Wang, Y. & Oberlack, M. 2021 A deterministic two-phase model for an active suspension with non-spherical active particles using eulerian spatial averaging theory. TUdatalib, Dataset,
https://doi.org/10.48328/tudatalib-781.
Dhar, A., Burada, P. S. & Sekhar, G. P. Raja 2020 Hydrody-

5

0

−0.2

0

0.2

x1

Figure 6. The PDF f v for the particle translation x1 -velocity(). While the PDF might look Gaussian, the comparison
with a normal distribution () reveals that especially forlarger possible velocities it is non-Gaussian.

the idea of a chaotic behaviour.
The increase of α S in the volume averaged simulation at
the channel walls can also be observed in the DNS-simulation.
The PDF f x (Fig. 5) gives the probability of finding a particle centre of mass at a certain x1 -coordinate of the domain. Itshould be noted, that the data shown in Fig. 5 are preliminary, asteady state has not jet been reached. Nevertheless, it alreadyindicate what is to expected from the final data. Within thedomain the distribution will be constant, whereas material accumulated at the domain wall, leading to the two peaks visible
in f x .
While the PDFs for the minimal distance f d and positionx
f indicate a structured behaviour, i.e. the formation of clusters, the analysis of the velocity PDF reveals a more chaotic
behaviour. In Fig. 6 the PDF for the v1 -velocity of the particles is plotted. Non-Gaussian ”heavy tails” are clearly visible, resembling the behaviour of a turbulent flow, see for example Wilczek & Friedrich (2009). This non-Gaussian behaviour could indicate that effects such as intermittency are
occurring; a phenomenon that is often studied in turbulence
research based on vorticity (Li & Meneveau, 2005).

5

12th International Symposium on Turbulence and Shear Flow Phenomena (TSFP12)
Osaka, Japan (Online), July 19-22, 2022
Shampine, Lawrence F & Reichelt, Mark W 1997 The matlab
ode suite. SIAM J. Sci. Comput. 18 (1), 1–22.
Shampine, Lawrence F, Reichelt, Mark W & Kierzenka,
Jacek A 1999 Solving index-1 daes in matlab and simulink.
SIAM Rev. 41 (3), 538–552.
Sharan, Priyanka, Postek, Witold, Gemming, Thomas,
Garstecki, Piotr & Simmchen, Juliane 2020 Study of active
janus particles in the presence of an engineered oil–waterinterface. Langmuir .
Skeel, Robert D & Berzins, Martin 1990 A method for the spatial discretization of parabolic equations in one space variable. SIAM J. Sci. Comput. 11 (1), 1–32.
Sokolov, Andrey & Aranson, Igor S. 2012 Physical propertiesof collective motion in suspensions of bacteria. Phys. Rev.
Lett. 109, 248109.
Soler, Lluı́s & Sánchez, Samuel 2014 Catalytic nanomotors for environmental monitoring and water remediation.
Nanoscale 6 (13), 7175–7182.
Wensink, Henricus H., Dunkel, Jörn, Heidenreich, Sebastian,Drescher, Knut, Goldstein, Raymond E., Löwen, Hartmut& Yeomans, Julia M. 2012 Meso-scale turbulence in living
fluids. Proc. Natl. Acad. Sci. USA 109 (36), 14308–14313.
Wilczek, Michael & Friedrich, Rudolf 2009 Dynamical origins
for non-gaussian vorticity distributions in turbulent flows.
Physical Review E 80 (1), 016316.
Wu, K.-T., Hishamunda, J.B., Chen, D.T.N., DeCamp, S.J.,
Chang, Y.-W., Fernandez-Nieves, A., Fraden, S. & Dogic,
Z. 2017 Transition from turbulent to coherent flows in confined three-dimensional active fluids. Science 355.
Wu, Xiao-Lun & Libchaber, Albert 2000 Particle diffusion
in a quasi-two-dimensional bacterial bath. Phys. Rev. Lett.
84 (13), 3017.
Yigit, Berk, Alapan, Yunus & Sitti, Metin 2020 Cohesive selforganization of mobile microrobotic swarms. Soft Matter
16 (8), 1996–2004.
Zhang, Zhuoran, Wang, Xian, Liu, Jun, Dai, Changsheng &Sun, Yu 2019 Robotic micromanipulation: Fundamentals
and applications. Annu. rev. control robot. auton. syst. 2,
181–203.

namics of active particles confined in a periodically taperedchannel. Phys. Fluids 32 (10), 102005.
Dombrowski, C., Cisneros, L., Chatkaew, S., Goldstein, R. E.
& Kessler, J. O. 2004 Self-concentration and large-scale coherence in bacterial dynamics. Phys. Rev. Lett. 93, 098103.
Drescher, K., Dunkel, J., Cisneros, L. H., Ganguly, S. & Goldstein, R. E. 2011 Fluid dynamics and noise in bacterial cell–
cell and cell–surface scattering. Proc. Natl. Acad. Sci. USA
108 (27), 10940–10945.
Drew, D. A. 1983 Mathematical modeling of two-phase flow.
Annu. Rev. Fluid Mech. 15 (1), 261–291.
Jurado-Sánchez, B., Sattayasamitsathit, S., Gao, W., Santos,
L., Fedorak, Y., Singh, V. V., Orozco, J., Galarnyk, M. &Wang, J. 2015 Self-propelled activated carbon janus micromotors for efficient water purification. Small 11 (4), 499–
506.
Kummer, F. 2017 Extended discontinuous galerkin methods
for two-phase flows: the spatial discretization. International
Journal for Numerical Methods in Engineering 109 (2),259–289.
Kummer, F., Müller, B. & Utz, T. 2018 Time integration forextended discontinuous galerkin methods with moving domains. Int. J. Numer. Methods Eng. 113 (5), 767–788.
Kummer, F., Weber, J. & Smuda, M. 2021 Bosss: A package for multigrid extended discontinuous galerkin methods.
Comput. Math. with Appl. 81, 237–257.
Leighton, David & Acrivos, Andreas 1985 The lift on a small
sphere touching a plane in the presence of a simple shear
flow. Z. Angew. Math. Phys. 36 (1), 174–178.
Li, Jinxing, Shklyaev, Oleg E, Li, Tianlong, Liu, Wenjuan,Shum, Henry, Rozen, Isaac, Balazs, Anna C & Wang,
Joseph 2015 Self-propelled nanomotors autonomously seek
and repair cracks. Nano Letters 15 (10), 7077–7085.
Li, Yi & Meneveau, Charles 2005 Origin of non-gaussianstatistics in hydrodynamic turbulence. Physical review letters 95 (16), 164502.
Orozco, Jahir, Mercante, Luiza A, Pol, Roberto & Merkoçi,
Arben 2016 Graphene-based janus micromotors for the dynamic removal of pollutants. J. Mater. Chem. A 4 (9), 3371–
3378.

6