Swimming at Low Reynolds Number

We address the swimming problem at low Reynolds number. This regime, which is typically used for micro-swimmers, is described by Stokes equations. We couple a PDE solver of Stokes equations, derived from the Feel++ finite elements library, to a quaternion-based rigid-body solver. We validate our numerical results both on a 2D exact solution and on an exact solution for a rotating rigid body respectively. Finally, we apply them to simulate the motion of a one-hinged swimmer, which obeys to the scallop theorem.


Introduction
One of the first works in the field of microswimming was proposed by E. Purcell, who presented in 1977 his observations regarding micro-swimmers' swimming strategies. In [33] he noticed that, at the microscale, Reynolds number is small, which means that viscous dissipative forces overcome inertial ones. This observation has two consequences: firstly, Navier-Stokes fluid equations, at the swimmer's scale, can be approximated very well by Stokes equations. These have the characteristics of instantaneously diffusing any change in momentum, and they also guarantee the adherence of the fluid to the boundary of the swimmer at every instant; secondly, certain types of deformations, that usually at human scale produce a net displacement, do not give any motion at the micro-scale. The well-known Scallop theorem [13,33] specifies some of them.
Since Stokes equations are time-invariant, the important factors, in terms of swimming, are just the body configurations and not the speed of deformation.
A typical value of Reynolds number Re = U L ν at the microscale is 10 −5 . In this case the length-scale of microscopic swimmers is taken to be L ≈ 10 −6 m, their mean velocity U ≈ 30 µm/s and water kinematic viscosity ν ≈ 10 −6 m 2 /s. Mathematical analysis of microswimming has mainly developed in the field of control theory, due to its applications to micro-robotics. In fact, it is a matter of great importance to understand how the propelling mechanisms of a robot influence its motion and its controllability. As microrobots might be used as vessels for targeted drug administration or as microsurgery devices, the necessity of control theory becomes clear.
In [35, chapter 16] some applications are mentioned: they can range from medical imaging of regions that cannot

The model
In this section, the equations that are considered along the paper are introduced. Section 1.1 is devoted to the description of the fluid model. Section 1.2 describes how the rigid body motion can be recovered given an imposed swimmer body deformation.

The fluid equations
In our problem we consider a micro-swimmer immersed in an incompressible fluid which is modeled, in the Eulerian reference frame, by Navier-Stokes equations ∂u ∂t + (u · ∇)u = ν∆u − ∇p ρ + f, where u(t, x) : R × R d → R d , d = 2, 3, is the fluid velocity field, p(t, x) : R × R d → R is the fluid pressure field, ν is the fluid kinematic viscosity, ρ is the fluid density and f (t, x) : R × R d → R d is the distribution of external forces per unit mass. The Navier-Stokes system, in the case of micro-swimmers, can be simplified to the Stokes equations by neglecting the inertial contributions of the flow. This simplification is justified by the small value of the Reynolds number associated to the flow regime in which we are interested. This can be seen in the adimensionalized version of Navier-Stokes equations (here without external forces, for simplicity) where u = u U , p = pL µU , t = L U and U is the swimmer's characteristic swimming speed, L is its characteristic length and ν is the fluid kinematic viscosity. In the momentum conservation equations, Reynolds number multiplies the inertial terms, which can therefore be neglected, since in the case of micro-swimmers Re ≈ 10 −5 .
As a result, the set of equations we will be using to model the fluid behaviour is the Stokes system Denoting by σ(t, x) = −p(t, x)I + µ(∇u(t, x) + ∇u(t, x) T ) the Newtonian fluid stress tensor and considering a time-dependent family of open sets Ω t ⊆ R d ∀t, we can therefore write the fluid problem as follows: find u(t, x), p(t, x) such that where ∂Ω t D and ∂Ω t N are the portions of ∂Ω t where Dirichlet and, respectively, Neumann boundary conditions are imposed, and n(t, x) is the outward unit normal to ∂Ω t . In particular, one requires ∂Ω t = ∂Ω t D ∪ ∂Ω t N and ∂Ω t D ∩ ∂Ω t N = ∅. Conditions onū andḡ are imposed so that the problem has a unique solution in [10]. In the case of our interest, Ω t will be R d deprived of the region occupied by the swimmer, and the Dirichlet boundary conditions correspond to the velocity field on the swimmer boundary. Such vector field can be decomposed in two parts: one represents the rigid motion and the other comes from the deformation of the swimmer's bodyū(t, , where x 0 is the center of mass of the swimmer. When this latter is known in advance, it is possible to recover the linear velocity V and the angular velocity ω by solving a linear system as presented in [27]. This method exploits the linearity of Stokes equations with respect to pressure and velocity, producing a system of equations which is linear in (V, ω) and that depends on the stresses produced by a number of "fundamentally-driven" flows.

The rigid body dynamics
Following Newton's second law of dynamics and Euler's equations for angular momentum, we can write the equations describing the dynamics of the swimmer. However, in our case, the microscopic nature of the immersed object and the Stokes approximation of fluid equations lead us to neglect the inertial terms. As a result, the aforementioned equations read respectively (see [27]) The equations above are written in the swimmer's reference frame. As proposed and developed in [27], the linearity of the Stokes equation in pressure and velocity allows a decomposition of the two fields as a sum of solutions to fundamentally driven Stokes equations. Each of these problems is characterized by the boundary conditions on the swimmer's surface. If we denote the canonical basis of the swimmer's reference frame by and by v(t, x) the deformation velocity we already mentioned, the fundamental boundary conditions are of the formū(t, x) = e i for i ∈ {1, 2, 3},ū(t, x) = e i ∧ x for i ∈ {1, 2, 3} andū(t, x) = v(t, x). They provide seven different Stokes problems and allow an expansion of the velocity and pressure fields as where problems 1 to 3 are characterized byū(t, x) = e i , problems 4 to 6 are characterized byū(t, x) = e i ∧ x and problem 7 byū(t, x) = v(t, x). This decomposition carries onto the fluid stress tensors, since they are linear with respect to pressure and velocity: . Substituting this decomposition into (1.2) and appropriately rearranging the terms brings us to the solution of a linear system in (V, ω) of 6 × 6 matrix M and right-hand side N ∈ R 6 as follows Once the velocities (in the swimmer's reference frame) are recovered by solving the previous system, as done in [27], we can determine the position and orientation of the swimmer in the laboratory reference frame by solving a set of ODEs of the form where q is the unit quaternion representing the orientation of the body in the laboratory fixed frame, R(q) is the rotation matrix associated to the unit quaternion q, * denotes the quaternion product and X is the position of the swimmer's centroid in the fixed reference frame. Appropriate initial conditions for the ODEs have to be chosen. The derivation of the quaternion equation and definition of the quaternion product are described in Section 2.
Choosing unit quaternions to represent rotations in R d is dictated by the fact that they cover SO (3), sending the pair {q, −q} to a unique R(q) ∈ SO (3), and avoid the classical complications that arise when Euler angles are used [16,24].

The self-propulsion constraints
In the previous section we addressed the computation of the linear and angular velocity of the swimmer in its reference frame. In order to model the displacement of a swimmer, the deformation velocity v(t, x) has to satisfy two constraints ensuring that the motion is entirely due to internal forces and is not relying on external actions. In particular, these constraints express the fact that in the swimmer's reference frame, linear and angular momenta remain unchanged because no external forces are applied [12]. In the linear momentum case, such invariance writes Following [12], (5) and its angular momentum analogue can be rewritten as

Numerical method
In Subsection 2.1 we address the discretization of the Stokes system and in Subsection 2.2 we will address the derivation and discretization of the ODE system for rigid body motion.

The fluid model
This part presents the PDE discretization that is used in the rest of the paper to numerically recover the swimmer's motion using finite element method. The indices denoting the time dependence of the fluid domain are abandoned and the focus is placed on the stationary problem to be solved at each time instant. The solution of the Stokes equations was performed via mixed finite elements, where the approximation spaces were chosen in order to guarantee the inf-sup stability of the discrete problem. The variational formulation is the following: (Ω) such that the boundary conditions in (1) They are characterized by different approximation orders on u and p in order to guarantee the inf-sup stability of the discrete problem. Namely, a continuous, piecewise polynomial approximation of order k is performed on each velocity component and a continuous, piecewise polynomial approximation of order k − 1 is performed on the pressure.
In the following, we assume that the computational triangulated domain coincides with the continuous one.
. Dirichlet boundary conditions are applied strongly at the algebraic level.

The rigid body dynamics
As shown in Subsection 1.2, from the finite element solution of the fluid problem it is possible to recover the rigid body dynamics of the swimmer. Unit quaternions are a good option to represent the orientation of a rigid object, since they avoid parametrization issues that can affect other representations, like the gimbal lock for Euler angles [16,24].
A quaternion q is defined by a scalar q 0 and a 3-dimensional vector q v = (q v1 , q v2 , q v3 ), q = (q 0 , q v ). The conjugate of q will be denoted byq = (q 0 , −q v ) and the quaternion product is defined as . Each quaternion of unit norm represents a reflection in 3D space, and together with its conjugate concurs to describe a rotation in space. The rotation matrix linked to a unit quaternion q = (q 0 , q v ) is given by In what follows we present how to get the swimmer frame dynamics in the laboratory frame. Coordinates X in the laboratory frame are expressed in the swimmer frame by [0 in our case 2[0, v] = 2 dq dt * q = [0, R(q)ω]. Then, multiplying both sides of the equality by q gives the following dynamics for the quaternion dq dt

Implementation
The implementation part combined the FEEL++ finite element library [32] to solve the Stokes problems with an ODE solver for the solution of the rigid body motion. Matlab and Paraview were used for the post-processing.
Since the methodology giving the rigid body problem was posed in the swimmer's reference frame and the swimmer's motion in the computational domain produced strong deformations on the underlying triangular mesh, at every time-step the domain was displaced and remeshed. This procedure was case specific, since both the shape and the deformation velocity on the swimmer's boundary were externally imposed.
Thanks to FEEL++, a Stokes solver with a syntax closely resembling the finite element variational formulation is achieved. Moreover, implementation of the finite elements was strongly simplified: Taylor-Hood spaces of order P 2 /P 1 are chosen, but also spaces of higher order (ex. P 3 /P 2 ) could be supported easily.
The solution of the ODEs (4) describing rigid-body motion was embedded in the code through an external ODE solver in C++. A Runge-Kutta scheme was performed. At the end of each time-step, before the new rotation matrix is computed, the quaternion needs to be normalized to have unit norm: the propagation of numerical errors does not preserve the norm of the quaternion. Remark that, if this renormalization is not performed, the rotation matrix won't be orthogonal and the body would suffer a deformation.

Validation of fluid solver
To validate the fluid solver presented in this paper, we compared the analytical solution of Stokes problem presented in [7] to the numerical results given by our simulations. The analytical solution to Stokes problem on with forcing function Dirichlet boundary conditions derived from the exact solution presented below and zero-mean pressure is given by p(t, (x, y)) = (x − 0.5)(y − 0.5).

Figure 1 presents this latter Stokes solution.
We solve the Stokes equations using the pair (X 2 h , M 1 h ) for the mixed finite elements. Figure 2 shows a third order convergence of the L 2 norm of the velocity and a second order convergence for both the L 2 norm of the pressure and the H 1 norm of the velocity are obtained. These results are those we expected from the finite element method theory, see for instance [21,Thm 4

.3, pag 181].
Remark that, even though the solver is made for the 2D case, Feel++ allows to easily obtain a three dimensional Stokes solution when providing an appropriate 3D mesh.

Validation of rigid body solver
In order to validate the rigid body solver, we considered the results presented in [28]. The solution to complemented with adequate initial conditions, is searched. Let us remark that the exactness of the translational part of the rigid body motion is easy to check. As a result, we will consider only equations (9b) and (9d).
We validate the rotational part of the rigid body motion by extracting a set of sample points from the graphics of [28] and compared qualitatively such results with ours. The case we will consider is the following: the time-varying torque in the body frame will be given by The body reference frame at t = 0 coincides with the laboratory reference frame, which in terms of quaternions is expressed as q(0) = (1, 0, 0, 0). The components of the angular velocity ω at t = 0 along the three principal axes are (0, 0, 0.329) rad/s, while the inertia matrix is given by Since in [28] the results about the frame orientation are given in terms of Euler angles (φ x , φ y , φ z ), we use the following relationship to obtain them from a quaternion q = (q 0 , q 1 , q 2 , q 3 ):    Figure 3 shows that our results are in agreement with the numerical results in [28]. We remark that, even though we will be focusing on a two-dimensional simulation for the swimmer, the previous validation confirms that such solver can be used also on a 3D experiment. The resulting angular velocity is the one required in (4) to obtain the body orientation.

Validation of force and torque computations
Given (u, p) solution of Stokes equations, this part addresses the numerical computation of fluid forces and torques. Since the swimmer motion depends on them, it is fundamental to validate the techniques we use to that are subsequently approximated via 3D volumic quadratures. Feel++ library was used to perform these computations. The numerical results will be compared to approximate and exact solutions. We considered two test cases: the first one is the 3D Stokes example of a translating and rotating sphere, for which the approximate solutions are provided in [23, p.67-68], while the second one is a 3D example taken from [14]. In the moving sphere case, the experimental setting is composed of a sphere of radius r = 0.01 placed in the middle of a cubic box of size 1, µ = 1 and U = ω = 2. The approximated forces and torques, given by the formulas F = −6πµrU and T = −8πµr 3 ω, give a qualitative target for the numerically derived forces and torques.
The results regarding the 3D Stokes case are collected in Tables 1 to 4. Tables 1 and 2 show the numerical results that were obtained when using X 2 h /M 1 h finite elements on a piecewise linear approximation of the geometry: Table 1 shows that the surface quadratures give better approximation of forces while Table 2 show that torques are better approximated when using volumic quadratures. Tables 3 and 4 show the numerical results that were obtained when X 3 h /M 2 h finite elements on a piecewise quadratic approximation of the geometry were used: in this case volumic quadrature performs better than surface quadratures in both cases.  Table 1. Forces comparison for the translating and rotating sphere with piecewise linear approximation of the geometry and X 2 h /M 1 h finite elements. In the last column, the numerical result of the Stokes formula F = −6πµrU is reported.   Table 3. Forces comparison for the translating and rotating sphere with piecewise quadratic approximation of the geometry and X 3 h /M 2 h finite elements. In the last column, the numerical result of the Stokes formula F = −6πµrU is reported.  Table 4. Torques comparison for the translating and rotating sphere with piecewise quadratic approximation of the geometry and X 3 h /M 2 h finite elements. In the last column, the numerical result of the Stokes formula T = −8πµr 3 ω is reported.
The second benchmark that was considered regarded the exact computation of forces from an exact and fully three-dimensional solution to Navier-Stokes equations (where we neglected the time variation) [14] u = −a   e ax sin(ay + dz) + e az cos(ax + dy) e ay sin(az + dx) + e ax cos(ay + dz) e az sin(ax + dy) + e ay cos(az + dx) [e 2ax + e 2ay + e 2az + 2 sin(ax + dz) cos(az + dx)e ay+az + 2 sin(ay + dz) cos(ax + dy)e az+ax + 2 sin(az + dx) cos(ay + dz)e ax+ay ] (12) in the cube [−1, 1] 3 , with a = π 4 and d = π 2 . The exact forces were calculated using the symbolic mathematics package Sympy, and they were compared to forces issued from surface and volume quadratures. The results in Figure 4 show that, after fixing a high quadrature order for both surface and volume quadratures (13 in our case), forces are better approximated when using volume quadratures on (10) than surface quadratures. These two benchmarks encouraged us to use the volume formulation to compute forces and torques for the Scallop simulation in the next section.

Scallop theorem
In this section we will address the simulation of the so called "Scallop theorem". This theorem was first presented by Purcell in [33], and later generalized in [13] for other types of swimmers, highlighting the underlying geometrical properties of the problem. Swimming is characterized by periodical strokes, and the scallop theorem specifies when they do not provide a net advancement. Scallop theorem, as in [33], was specialized to Stokes flows, and stated that a swimmer which performs time-reversible shape changes will not be able to advance over a period. A particular case, from which the theorem takes his name, is the case of one-hinged swimmers. Having just one degree of freedom, swimmers like the scallop are forced to perform time-reversible shape changes and are not able to obtain a net advancement over a period of time. Figure 5 shows an example of these shape changes.
The aim of this part is to verify that our numerical solver is able to reproduce the scallop theorem. Problem (1) is solved and coupled with (4) to obtain the swimmer's displacement and change in orientation. Here Ω 0 is the complement in [0, 0.1] × [0, 0.1] of the scallop and the deformation velocity on the swimmer's boundary point x is given by v(t, x) = dA dt (t)x, where A(t) is the rotation matrix imposing the desired opening and closing of the valves.
We imposed on the boundary of the swimmer a displacement field which caused its hinges to open and close in a periodical fashion. We also allow a slight deformation of the body, to see whether any changes can be experienced with respect to the invariance predicted in [13,Section 3.5].
The computational domain coincides with Ω 0 , hence it consists in a pierced square whose hole has the shape of the swimmer. The mesh is conformal to the swimmer's geometry and the approximation is linear: interior and boundary edges are straight. Figure 6 shows that the scallop theorem is not exactly satisfied due to the domain discretization, but the errors over a period, and the errors cumulated between two consecutive periods, are negligible. We chose a period T = 2 s and a swimmer with total valve opening equal to 14 mm. The smallest triangulation edge size in is of 0.025 mm.  The figure on the right shows the swimmer's vertical displacement: it is periodic of period T = 2, and the net motion value on a period is also imputable to numerical errors.
We also located the scallop close to a no slip boundary of the computational domain. In that case the scallop theorem is not satsfied since the wall breaks the symmetry of the whole system. Figure 7 shows the displacement over a complete stroke period as a function of the proximity to the wall. In that case we noticed that the displacement at the end of the period was sensibly larger than what Figure 6 showed, and it exemplifies both the attractive effect that walls have on swimmers (Figure 7 right) and the enhanced displacement parallel to no slip boundaries (Figure 7 left). Let us quote the recent study [6] that overcomes the scallop theorem by opening and closing the shell with different speeds.
The pseudocode for the simulation of the swimmer motion is presented in Algorithm 1.

Conclusion and perspectives
In this note we present a solver able to provide the dynamics of a derforming swimmer moving into a Stokes flow. Here all the framework could fit a 3D model but the numerical results are given in 2D. All the numerical computations of this note relied on the finite element library Feel++. We applied a method proposed in [27] to extract rigid body motion from computation of fluid stresses. We validated our code and model using different benchmarks. The goal of this project was to lay the foundations of a numerical method which allowed the simulation of swimmers with more complex geometries. Displace the swimmer according to X, q, V, ω and Remesh for 1 ≤ i ≤ N Stokes do Solve Stokes problem i from Subsection 1.2 Compute surface forces and torques issued from Stokes problem i end for Find V, ω by using (2) and (3) Find X, q by solving (4) time ← time + ∆t end while The passage to several swimmers or to 3D are beyond the scope of this paper. The technique could be applied to several swimmers, but the basis onto which the Stokes solution has to be expanded is larger and the computational cost might increase noticeably. Another issue is that, in the Stokes regime, forces can become infinite when two swimmers are close, and collision avoidance strategies would need to be implemented.