Let \(\Omega \subseteq \mathbb{R}^d\), where \(d=2\) or \(d=3\), be the domain of the fluid, \(u : \Omega \mapsto \mathbb{R}^d\)
and \(p : \Omega \mapsto \mathbb{R}\) respectively its velocity
and its pressure. We note \(\mu \in \mathbb{R} ^ {+}\) the dynamic viscosity of the fluid, \(\rho_f \in \mathbb{R}^{+}\) its density
and \(\sigma = -p I_d + \mu ( \nabla u + \nabla u^T)\) the stress tensor.
Let \(B_i\) be the \(i\)-th rigid body and \(\partial B_i\) its boundary.
We define \(d^{*} = 1\) if \(d=2\) and \(d^{*} = 3\) if \(d=3\).
We note \(U_i \in \mathbb{R}^d\) and \(\omega_i \in \mathbb{R}^{d^{*}}\) the translational and
angular velocities of \(B_i\). \(F^e_i\) and \(M^e_i\) are the external force and momentum
acting on the \(i\) body, \(x^{CM}_i \in \mathbb{R}^d\) its center of mass, \(m_i \in \mathbb{R} ^ {+}\) its mass and \(J_i \in \mathbb{R}^{d^{*}\times d^{*}}\), positive definite, its momentum of inertia. Let \(R\) be the rotation matrix that expresses the change in reference frame from the body’s to the laboratory’s frame. The dynamics of the bodies, falling under the effect of gravity, is described by the Newton equations:
\[\begin{align*}\label{p:2}
\left\{\begin{array}{rcl}
m_i \frac{d U_i}{d t} &=& (m_i - \int_{B_i} \rho_f )g + F^e_i \mbox{,}\\
\frac{d [ R J_i R^T \omega_i ]}{d t} &=& M^e_i \mbox{.}
\end{array}\right.\;
\end{align*}\]
where \(g \in \mathbb{R}^d\) is the gravity vector.
The hydrodynamical forces \(F^{e}_i\) and torque \(M^{e}_i\) acting on the body \(B_i\) are defined by:
\[ F^{e}_i = - \int_{\partial B_i} \sigma \cdot \overrightarrow{n} \quad \mbox{and} \quad M^{e}_i = - \int_{\partial B_i} \sigma \cdot \overrightarrow{n}\times (x_i - x^{CM}_i) \mbox{,}\]
where \(\overrightarrow{n}\) is the external unit vector on \(\partial B_i\).
The fluid-body interaction is defined by a coupling condition on the velocities on the boundary of the bodies \(\partial B_i\):
\[u = U_i + \omega_i \times (x - x^{CM}_i)\mbox{.}\]
The fluid-body problem links these equations by the
following system:
\[\begin{equation*}\label{p:3}
\left\{\begin{array}{rcl}
\rho_f \frac{\partial u}{\partial t} + \rho_f (u \cdot \nabla)u - \mu \Delta u + \nabla p &=& 0_{\mathbb{R}^d} \mbox{,} \quad \mbox{in } \Omega \texttt{,} \\
\nabla \cdot u &=& 0 \mbox{,} \quad \quad \mbox{in } \Omega \texttt{,}\\
u &=& U_i + \omega_i \times (x - x^{CM}_i) \mbox{, } \quad \mbox{on } \partial B_i \texttt{,}\\
m_i \frac{dU_i}{dt} &=& (m_i - \int_{B_i} \rho_f)g + F^{e}_i \texttt{,}\\
\frac{d [ R J_i R^T \omega_i ]}{d t} &=& M^{e}_i \texttt{.}
\end{array}\right.\;
\end{equation*}\]
Let \(X = H^1(\Omega)^d, M = L^2_0(\Omega), T = \{U | U(t) \in \mathbb{R}^d\}, W = \{\omega | \omega(t) \in \mathbb{R}^{d^{*}}\} \). The variational formulation is given by :
Find \((u,p,U,\omega) \in X \times M \times T\times W \) such that:
\[\begin{align*}
\int_{\Omega} \rho_f \frac{\partial u}{\partial t} \cdot v + \int_{\Omega} \rho_f (u \cdot \nabla)u \cdot v + 2 \mu \int_{\Omega} D(u):D(v) - \int_{\Omega} p \nabla \cdot v &= \int_{\partial B_i} \sigma \cdot \overrightarrow{n} \cdot v + \int_{\Omega} f \cdot v \;, \quad \forall v \in X \mbox{,} \\
\int_{\Omega} q \nabla \cdot u &= 0 \;, \quad \forall q \in M \mbox{,}\\
m_i \frac{d U_i}{d t} \cdot \tilde{U}_i &= (m_i - \int_{B_i} \rho_f)g \cdot \tilde{U}_i + F^{e}_i \cdot \tilde{U}_i \;, \quad \forall \tilde{U}_i \in T \mbox{,} \\
\frac{d [ R J_i R^T \omega_i ]}{d t} \cdot \tilde{\omega}_i &= M^{e}_i \cdot \tilde{\omega}_i \;, \quad \forall \tilde{\omega}_i \in W \mbox{.}
\end{align*}\]
The test functions satisfy the constraint on the velocities on the boundary of the body. Thus, the boundary conditions are rewritten as:
\[\begin{align*}
v &= \tilde{U}_i + \tilde{\omega}_i \times (x - x^{CM}_i)\\
\Rightarrow \int_{\partial B_i} \sigma \cdot \overrightarrow{n} \cdot v &= \int_{\partial B_i} \sigma \cdot \overrightarrow{n} \cdot \tilde{U}_i + \int_{\partial B_i} \sigma \cdot \overrightarrow{n} \times (x - x^{CM}_i) \cdot \tilde{\omega}_i \\
\Rightarrow \int_{\partial B_i} \sigma \cdot \overrightarrow{n} \cdot v &= F^{e}_i \cdot \tilde{U}_i + M^{e}_i \cdot \tilde{\omega}_i \\
\Rightarrow \int_{\partial B_i} \sigma \cdot \overrightarrow{n} \cdot v &= - m_i \frac{dU_i}{dt} \cdot \tilde{U}_i - (m_i - \int_{B_i} \rho_f)g \cdot \tilde{U}_i - \frac{d [ R J_i R^T \omega_i ]}{d t} \cdot \tilde{\omega}_i \mbox{.}
\end{align*}\]
One can insert this new expression for the boundary conditions into the variation formulation and one obtains:
Find \((u,p,U,\omega) \in X \times M \times T \times W\) such that:
\[\begin{align*}
\int_{\Omega} \rho_f \frac{\partial u}{\partial t} \cdot v + \int_{\Omega} \rho_f (u \cdot \nabla)u \cdot v + 2 \mu \int_{\Omega} D(u):D(v) - \int_{\Omega} p \nabla \cdot v &= - m_i \frac{dU_i}{dt} \cdot \tilde{U}_i - (m_i - \int_{B_i} \rho_f)g \cdot \tilde{U}_i - \frac{d [ R J_i R^T \omega_i ]}{d t} \cdot \tilde{\omega}_i + \int_{\Omega} f \cdot v \;, \quad \forall v \in X \mbox{, } \tilde{U} \in T \mbox{, } \tilde{\omega} \in W \mbox{,} \\
\int_{\Omega} q \nabla \cdot u &= 0 \;, \quad \forall q \in M \mbox{.}
\end{align*}\]
\[\begin{align*}
a(u,v) &= \int_{\Omega} \rho_f \frac{\partial u}{\partial t} \cdot v + \int_{\Omega} \rho_f (u \cdot \nabla)u \cdot v + 2 \mu \int_{\Omega} D(u):D(v) \mbox{,} \\
b(v,p) &= - \int_{\Omega} p \nabla \cdot v \mbox{,} \\
G(v) &= \int_{\Omega} f \cdot v \mbox{,} \\
\end{align*}\]
The variational formulation can be rewritten as:
Find \((u,p,U,\omega) \in X \times M \times T \times W\) such that:
\[\begin{align*}
a(u,v) + b(v,p) &= - m_i \frac{dU_i}{dt} \cdot \tilde{U}_i - (m_i - \int_{B_i} \rho_f)g \cdot \tilde{U}_i - \frac{d [ R J_i R^T \omega_i ]}{d t} \cdot \tilde{\omega}_i + G(v) \;, \quad \forall v \in X \mbox{, } \tilde{U} \in T \mbox{, } \tilde{\omega} \in W \mbox{,} \\
b(u,q) &= 0 \;, \quad \forall q \in M \mbox{.}
\end{align*}\]
We approximate this solution by the Galerkin method using
a \(\mathbb{P}_2/\mathbb{P}_1/\mathbb{P}_0\) finite element discretization.
We consider in a first step a single body, so \(i=1\). To consider the conditions on the boundary of the body \(\partial B\) without having
to compute the integral, we distinguish the degrees of freedom of \(u\) on this boundary,
denoted by \(u_{\partial B}\) and on the remaining domain, \({u_{\Omega}}\).
Using the same approach as for the Stokes case, one find the following coupling-free matrix formulation:
\[\begin{bmatrix}
A_{\Omega, \Omega} & A_{\Omega, \partial B} & 0 & 0 & B_{\Omega}^T \\
A_{\partial B, \Omega} & A_{\partial B, \partial B} & 0 & 0 & B_{\partial B}^T \\
0 & 0 & m I_d & 0 & 0 \\
0 & 0 & 0 & R J R^T & 0 \\
B_{\Omega} & B_{\partial B} & 0 & 0 & 0 \\
\end{bmatrix} \begin{bmatrix}
u_{\Omega}\\
u_{\partial B}\\
U \\
\omega \\
p
\end{bmatrix} = \begin{bmatrix}
G_{\Omega}\\
G_{\partial B}\\
(m - \int_{B} \rho_f)g\\
0\\
0
\end{bmatrix} \mbox{,}\]
\[A_{J,K} = (a( \phi_{J_i} , \phi_{K_j} ))_{i,j} \;, \quad (\phi_J)_i \mbox{ basis of } (X_h)_J \mbox{ and } (\phi_K)_j \mbox{ basis of } (X_h)_K \\
B_{J} = (b( \phi_{J_i} , \varphi_{K_j} ))_{i,j} \;, \quad (\phi_J)_i \mbox{ basis of } (X_h)_J \mbox{ and } (\varphi_K)_j \mbox{ basis of } (M_h)_K \\
G_{J} = (G(\phi_{J_i}))_{i} \;, \quad (\phi_J)_i \mbox{ basis of } (X_h)_J \mbox{.}\]
The fluid-body coupling is then performed using an operator \(\mathcal{P}\).
This operator satisfies the following equality:
\[\mathcal{P} \begin{bmatrix}
u_{\Omega}\\
U\\
\omega\\
p
\end{bmatrix} = \begin{bmatrix}
u_{\Omega}\\
u_{\partial B} = U + \omega \times (x - x^{CM}) = U \tilde{P}_U + \omega \tilde{P}_\omega \\
U\\
\omega\\
p
\end{bmatrix}\]
Thus, the definition of \(\mathcal{P}\) is given by:
\[\mathcal{P} = \begin{bmatrix}
I & 0 & 0 & 0 \\
0 & \tilde{P}_U & \tilde{P}_\omega & 0 \\
0 & I & 0 & 0 \\
0 & 0 & I & 0 \\
0 & 0 & 0 & I
\end{bmatrix} \mbox{,}\]
where \(I\) is the identity matrix.
The final matrix formulation is then defined by:
\[\mathcal{P}^T \begin{bmatrix}
A_{\Omega, \Omega} & A_{\Omega, \partial B} & 0 & 0 & B_{\Omega}^T \\
A_{\partial B, \Omega} & A_{\partial B, \partial B} & 0 & 0 & B_{\partial B}^T \\
0 & 0 & m I_d & 0 & 0 \\
0 & 0 & 0 & R J R^T & 0 \\
B_{\Omega} & B_{\partial B} & 0 & 0 & 0 \\
\end{bmatrix} \mathcal{P} \begin{bmatrix}
u_{\Omega}\\
U \\
\omega \\
p
\end{bmatrix} = \mathcal{P}^T \begin{bmatrix}
G_{\Omega}\\
G_{\partial B}\\
(m - \int_{B} \rho_f)g\\
0\\
0
\end{bmatrix}\]
.pdf