The model

The FSI model is a coupling of partial differential equations describing fluid and solid mechanics.

We will use the subscript \(f\) to denote fluid-related quantities and the subscript \(s\) to denote the solid related ones. We will indicate with the superscript \( ^ * \) the variables and domains expressed in Lagrangian coordinates.

We use Navier-Stokes equations to model the fluid in its domain \(\Omega_f\), and they read as \[ \begin{aligned} &\rho_f \frac{\partial u_f}{\partial t} + \rho_f(u_f \cdot \nabla) u_f = \mu_f \Delta u_f - \nabla p_f + f_f \\ &\nabla \cdot u_f =0 \end{aligned} \] where \(u_f=u_f(t,x)\) is the velocity field, \(p_f=p_f(t,x)\) is the pressure field, \(\rho_f\) is the constant fluid density, \(\mu_f\) is the dynamic viscosity and \(f_f=f_f(t,x)\) is the resultant of the external volume forces. We point out that the Navier-Stokes equations are written in Eulerian coordinates.

This set of equations can be derived from Cauchy equations of continuum mechanics and mass balance equation in the hypothesis of incompressible Newtonian fluid, for which the fluid strain tensor \(\sigma_f = -p_f\mathbb{I} +2\mu_f D_f\) is linear with respect to the strain rate \(\nabla u\), with \[ D_f=\frac{1}{2} (\nabla u_f + \nabla u_f ^ T). \] This set of equations can be therefore rewritten as \[ \begin{aligned} &\rho_f \frac{\partial u_f}{\partial t} + \rho_f(u_f \cdot \nabla) u_f -\nabla \cdot \sigma_f= f_f \\ &\nabla \cdot u_f =0. \end{aligned} \] The solid swimmer is again modeled using Cauchy equations of continuum mechanics, with an appropriate constitutive model for the stress tensor. However, since solid mechanics is usually written in Lagrangian coordinates, we respect this convention and have, on \(\Omega_s ^ * \), the Cauchy equations which read as follows \[ \rho_s ^ * \frac{\partial ^ 2 \eta_s}{\partial t ^ 2} - \nabla \cdot (F_s\Sigma_s) = f_s \] where \(\rho_s ^ * \) is the solid density, \(\eta_s=\eta_s(t,x ^ * )\) is the displacement, defined as the difference \(\eta(t,x ^ * )=x(t,x ^ * )-x ^ * \) between the solid configuration at time \(t\) and the initial configuration, \(f_s=f_s(t,x ^ * )\) is the resultant of external volume forces and \(F_s=F_s(t,x ^ * )=\mathbb{I} + \nabla \eta_s(t,x ^ * )\) is the strain tensor.

The tensor \(\Sigma_s=\Sigma_s(t,x ^ * )\), called second Piola-Kirchhoff tensor, depends on the constitutive solid model: based on the swimmer geometrical characteristics and its swimming dynamics, which consist in moving the tail with ample motions, we chose a Saint-Venant-Kirchhoff hyper-elasticity model, for which \[ \Sigma_s(t,x ^ * ) = \lambda_s tr(E_s)\mathbb{I} + 2\mu_s E_s \] where \(E_s=E_s(t,x ^ * )=\frac{1}{2}(F_s ^ T F_s - \mathbb{I})\) is the deformation tensor and \[ \lambda_s = \frac{E\nu}{(1+\nu)(1-2\nu)} \quad \mu_s= \frac{E}{2(1+\nu)} \] are the Lamé coefficients, calculated using the Young modulus \(E\) and Poisson coefficients \(\nu\) of the solid. For hyper-elastic materials, \(\Sigma_s\) can be described in terms of an energy function \(W(E_s)\) as \(\Sigma_s(E_s) = \frac{\partial W}{\partial E}\), and in the Saint-Venant-Kirchhoff model this energy is given by \[ W(E_s) = \frac{\lambda_s}{2} tr(E_s) ^ 2 + \mu_s tr(E ^ 2) \] The fluid-solid interaction model is encoded in the continuity conditions that one imposes on the common boundaries. The first issue to address when coupling Navier-Stokes and Cauchy equations for solids is the different coordinate systems in which they are formulated. For this reason we define the ALE maps \(\mathcal{A}_f ^ t: \Omega_f ^ * \rightarrow \Omega_f ^ t\) which form a time dependent family of homeomorphisms linking the fluid volume at the initial time \(\Omega_f ^ * \) with the fluid volume at the current time \(\Omega_f ^ t\) and the family of maps \(\varphi_s ^ t : \Omega_s ^ * \rightarrow \Omega_s ^ t\) which do the same job on the solid domain. The two families map each point \(x ^ * \) in the reference domain to its current location \(x=x(t,x ^ * )\).

At this point, we have the two models, solid and fluid, expressed in two different coordinate systems: ALE formulation allows to maintain this difference, but we will rewrite the Eulerian time derivative in the Navier-Stokes equations. We will perform this later on, after having presented the ALE coupling conditions: \[ \begin{aligned} &\frac{\partial \eta_s}{\partial t} - u \circ \mathcal{A}_f ^ t = 0 \\ &F_s\Sigma_sn_s ^ * + F_{\mathcal{A}_f ^ t } ^ {-T} \hat{\sigma_f} n_f ^ * \det({\nabla_x \mathcal{A}_f ^ t } )= 0\\ & \varphi_s ^ t -\mathcal{A}_f ^ t = 0 \end{aligned} \] where \(F_{\mathcal{A}_f ^ t}=\nabla_x \mathcal{A}_f ^ t\).

The first equation imposes the velocity continuity at the interface \(\Gamma_{fsi} ^ * \), expressed in the Lagrangian reference frame; the second equation prescribes the continuity of the stresses and the third one the geometric continuity at the boundary \(\Gamma_{fsi} ^ * \).

The second equation, expressing the stress tensors in Lagrangian coordinates, needs to be developed: \(F_s\Sigma_sn_s ^ * \) is the stress tensor on the solid structure side, and has the same expression as in the Cauchy equation; \(J_{\mathcal{A}_f ^ t}F_{\mathcal{A}_f ^ t} ^ {-T}\hat{\sigma_f}n_f ^ * \) is the stress tensor in the fluid side: in Eulerian coordinates it would be just \(\sigma_fn_f\). To express it in Lagrangian reference frame, we need to do the following change of coordinates: \[ \begin{aligned} &\int_{\Omega_f ^ t} \sigma_f(t,x)n_f \mathrm{d}x= \int_{\Omega_f ^ * } F_{\mathcal{A}_f ^ t } ^ {-T} F_{\mathcal{A}_f ^ t } ^ {T} \sigma_f(t,\mathcal{A}_f ^ t(x ^ * )) F_{\mathcal{A}_f ^ t } F_{\mathcal{A}_f ^ t } ^ {-1} n_f \det(F_{\mathcal{A}_f ^ t })\mathrm{d}x ^ * =\\ &\int_{\Omega_f ^ * } F_{\mathcal{A}_f ^ t } ^ {-T} \hat{\sigma}_f(t,x ^ * )n_f ^ * \det(F_{\mathcal{A}_f ^ t })\mathrm{d}x ^ * \end{aligned} \] The Cauchy tensor, a 2-covariant tensor field, undergoes a change of coordinates which is the pull-back along the ALE map \(\mathcal{A}_f ^ t\).

As a last step, we need to rewrite Navier-Stokes equations in the Arbitrary-Lagrangian-Eulerian coordinate system: the only term we will change is the Eulerian temporal derivative. Given a function \(g:\Omega_f ^ t \times [0,T\) \rightarrow \mathbb{R}] defined in the Eulerian reference frame, we call \(\hat{g}=g(\mathcal{A}_f ^ T(x ^ * ),t)\) its equivalent in the ALE frame. We define the ALE temporal derivative as

\[ \frac{ \partial g }{ \partial t} |_{x ^ * } (x,t) = \frac{ \partial \hat{g} }{ \partial t} \Big( (\mathcal{A} ^ t) ^ {-1} (x),t \Big) \]

and the domain deformation velocity (ALE velocity) as

\[ w_f(x,t) = \frac{ \partial x}{ \partial t} |_{x ^ * } (x,t) = \frac{ \partial x}{ \partial t} \Big( (\mathcal{A} ^ t ) ^ {-1} (x),t \Big) \]

If we denote the Eulerian velocity as \(\frac{\partial}{\partial t}|_x\), we have that

\[ \frac{ \partial u}{ \partial t}|{x ^ * } = \frac{ \partial u}{ \partial t} |{x} + w_f \cdot \nabla_x u_f \]

and the Navier-Stokes equations can be rewritten as follows

\[ \begin{aligned} &\rho_f \frac{ \partial u_f}{ \partial t} |_{x ^ * } + \rho_f((u_f-w_f) \cdot \nabla_x) u_f -\nabla_x \cdot \sigma_f= f_f \\ &\nabla_x \cdot u_f =0. \end{aligned} \]

The reason why we have performed this substitution is to split the contributions due to the domain deformation and the "fluid particle" velocity.

When elasticity equations are complemented by pure Neumann boundary conditions, they are invariant with respect to infinitesimal rigid body motions, and a customary approach to solve them numerically is to look for a solution which belongs to the orthogonal space to such kernel.

Rigid body motion must be applied independently from the FSI model, and one can choose between two options: the Newton-Euler set of ordinary differential equations, or Lagrangian multipliers on the PDE model.

Newton-Euler differential equations describe the dynamics of a reference frame that is attached to the rigid body: its origin and the orientation of its axes with respect to the laboratory reference frame are the unknowns of the problem. The Newton equations, expressing the position of the center of mass of a rigid body, which coincides with the origin of the reference frame that we have mentioned above, can be written as follows:

\[ \left\{ \begin{aligned} \dot{x}&=v \\ \dot{v} &=\frac{F}{m} \end{aligned} \right. \]

where \(F\) is the resultant of external body forces. In our case, \(F\) will be given by the integral of fluid forces over the boundary of the swimmer. The Euler equations describe the rotational dynamics of the rigid body, by describing the dynamics of the attached reference system. There are several ways to write these equations, but we choose to use the version which involves unit quaternions, since it does not present any singular configuration. The equations can be written as follows:

\[ \left\{ \begin{aligned} \dot{q} &=\frac{1}{2}[0, \omega] * q \\ \dot{(I\omega)} &= T \end{aligned} \right. \]

where \(I\) is the inertia tensor of the rigid body, written in the moving principal frame, \(T\) is the external torque and \(a* b\) represents the product between two quaternions.

Another viable approach is the usage of Lagrange multipliers to impose rigid translations and rotations on the elastic body.

We aim to reproduce the laboratory experiment. In order to do so we need to input in our simulation the same fluid and solid parameters that characterized the experimental setting.

The microswimmer is made of two parts: the head is a magnetic cylinder made of Neodymium-Iron-Boron (NdFeB) N-45 alloy and the tail is realised by means of a 3D printer using a bi-component (base and catalyst) poly-additive RTV (Room-Temperature-Vulcanizing) silicon elastomer, which solidifies at ambient temperature. The parameters that we need to insert in the model, regarding the solid part, are the Young modulus (E), the Poisson coefficient (\(\nu\)) and the density (\(\rho_s\)). For what concerns the magnet, we have \(\rho_s=7500-7800\) \(kg/m ^ 3\), Young modulus \(E=150-160\) \(GPa\) and Poisson coefficient \(\nu=0.24\).

The fluid, instead is a mixture of water and glycerine, whose parameters are tuned in order to get a low Reynolds number \(Re\approx 10 ^ {-4}\) so that, despite the macroscopic dimensions of the swimmer, the conditions are those of a micro-swimmer. The fluid density is \(\rho_f=1260\) \(kg/m ^ 3\) while its dynamic viscosity is \(\mu_f=1.524\) \(kg/(m\cdot s)\).

The magnetization is aligned with the cylinder’s head, and it is also aligned with the swimmer’s rotational symmetry axis.

The swimmer is globally \(8.5 mm\) long and its tail is \(8 mm\) long and conic shaped.