Variational Integrators and Fluid‐Structure‐Interaction at Low Reynolds‐Number

For the special case of incompressible and highly viscous fluids, the interaction with a rigid body can be collected in a damping matrix, relating the velocities and angular velocities of the body with the fluid force and torque. This damping matrix (a.k.a. viscous resistance matrix) depends exclusively on the geometry and needs to be computed only once. We consider a rigid body moving in an unbound fluid. The generalized dissipative forces from the fluid onto the body enter the time discretization via the discrete D'Alembert Principle. As generally large rotations may occur, we chose quaternions for a singularity‐free description of the body orientation. The corresponding holonomic constraint of a unit quaternion is enforced on the position and momentum level by the RATTLE algorithm. The problem of Stokes drag on a sedimenting particle serves as an example.

For the special case of incompressible and highly viscous fluids, the interaction with a rigid body can be collected in a damping matrix, relating the velocities and angular velocities of the body with the fluid force and torque. This damping matrix (a.k.a. viscous resistance matrix) depends exclusively on the geometry and needs to be computed only once. We consider a rigid body moving in an unbound fluid. The generalized dissipative forces from the fluid onto the body enter the time discretization via the discrete D'Alembert Principle. As generally large rotations may occur, we chose quaternions for a singularity-free description of the body orientation. The corresponding holonomic constraint of a unit quaternion is enforced on the position and momentum level by the RATTLE algorithm. The problem of Stokes drag on a sedimenting particle serves as an example.

Problem formulation
Variational Integrators [1], originally developed for mechanical systems, have been successfully extended to further domains and coupled problems. The dynamics of inviscid incompressible fluids and their interaction with structures fit well into the variational framework [2]. For the other extreme of highly viscous fluids (Re < 1) there is a variational formulation [3], too Indeed, this principle is very similar to the minimum of potential energy in statics, of the velocities v as displacements, the viscosity µ as elasticity, and of the pressure p as Lagrange multiplier to enforce incompressibility. However, eq. (1) is not compatible to Hamilton's Principle, which is already apparent by its different units. Although the resulting fluid equations are linear, analytical solutions are only available for simple geometries, such as for two planar eccentric circles [4]. For general problems we need to resort to numerical methods. A popular choice are finite elements [3], leading to a high dimensional discretization of the fluid domain Ω. An approximation in the sense of Trefftz [5], which fulfills the field equations exactly and minimizes the error at the boundary Γ, leads astray, since trial functions fulfilling ∇ · σ = 0 do not generate forces or torques around any closed contour. We decide for another approach, a boundary integral formulation. It is promising, as it reduces the geometric dimension by one; hence, in three dimensions only a discretization of the boundary surface is needed.
Having multi-body systems in mind, we consider here the model problem of a sedimenting particle in an unbound fluid, which turns out easy, as the fluid-structure interaction reduces to the classical linear damping matrix with constant entries. For brevity, we use a non-dimensional formulation with the reference length L, reference speed U , and consequently reference time T = L/U , reference acceleration G = U 2 L , reference pressure P = µU L and reference mass M = µ L 2 U .

Fluid flow and interface forces
We define the stress σ = −pI + ∇v + ∇ T v to write Stokes equations, governing the fluid dynamics [6], compactly assuming b = 0. The linearity of the field equations leads to a linear force-velocity relation for the immersed body [7] f where the damping matrix D is constant in the body-fixed frame (X-Y -Z). This matrix depends only on the geometry of the body and will be calculated next. We use Green's functions for Stokes flow, a.k.a. Stokeslet, with J(r) = 1 8π I r + rr This integral equation can be solved for the unknown surface tractions t(r p ) by a quadrature rule, made simplest by assuming constant values on the finite boundary elements. However, an apparent difficulty is that the integrands at y = r are singular. As a workaround, we follow Power and Miranda [8] to rewrite them as proper integrals and set up a linear equation system for the surface tractions at the meshed particle boundary [9]. Once the discretized surface tractions t d (r p ) are found, fluid force and torque result from numerical integration of expressions like F X = Γp t d X dΓ and M X = Γp t d Z Y − t d Y Z dΓ on the particle boundary. We solve for the six degrees of freedom of the particle separately, meaning we set one velocity component to unity while the others are zero, from which we assemble the damping matrix D. Note that the asymptotic behavior at the outer boundary v(r) → 0 and p(r) → 0 for |r| → ∞ is already implied by the Stokeslet solution.
3 Dynamics of an immersed particle in an unbound fluid The motion of the particle is conveniently described by its center of mass (q x , q y , q z ) and an unit quaternion (q 0 , q 1 , q 2 , q 3 ) for its orientation. The formulation of the Lagrangian L = 1 2q T Mq − gq z is straightforward. The virtual work of the nonconservative forces δW nc =q T R T DRδq reflects the fact, that the time derivativesq = [q x ,q y ,q z ,q 0 ,q 1 ,q 2 ,q 3 ] T expressed in the space-fixed frame (x-y-z) have to be converted to the velocities and angular velocities expressed in the body-fixed frame by the relation v = Rq. The holonomic constraint φ(q) = q 2 0 + q 2 1 + q 2 2 + q 2 3 − 1 = 0 of the unit quaternion is included by the RATTLE algorithm [10], written in the customary notation of Variational Integrators [1] Indeed, RATTLE corresponds to a Variational Integrator, enforcing the holonomic constraint on position and momentum level by the Lagrange multipliersλ andλ, respectively. Observe that, the mass matrix M (non-dimensional) consists of a diagonal block M 3 = diag(1, 1, 1) for the translational part of the kinetic energy and an augmented mass matrix for the rotational part M 4 [11], which is invertible and thus allows the momentum definition p = Mq to be inverted, namelyq = M −1 p. Fig.1 shows the results for the sedimentation of an ellipsoid (a = 2, b = 1, c = 1.5) subjected to gravity (g = −ge z ). Its body-fixed frame is initially aligned with the space-fixed frame and starts with v x = 10, v y = v z = 0, ω x = ω y = 0 and ω z = 10. The damping matrix was found to be D = diag(26.13, 29.92, 27.68, 129. 2, 194.6, 161.9).
Future work is related to wall effects, multi-particle systems and non-vanishing velocity at infinity. Although the fluidstructure interaction will be more complex then, it still fits into the discrete D'Alembert Principle.