2. Formulation of the problem

Consider a physical system, comprised of a rigid body ${\mathscr{B}}$ (closure of a bounded domain of $\mathbb{R}^3$ of class $C^2$ ) moving in a Navier–Stokes liquid ${\mathscr{L}}$ that fills the whole space outside ${\mathscr{B}}$ , under the action of a force $\boldsymbol{\mathsf{{f}}}$ . We assume that the motion of ${\mathscr{B}}$ is translatory, which can be achieved by applying on ${\mathscr{B}}$ a suitable torque preventing rotation. Concerning $\boldsymbol{\mathsf{{f}}}$ , we suppose that it is time periodic of period $T$ ( $T$ -periodic) and that its average, $\overline {\boldsymbol{\mathsf{{f}}}}$ , over a period is zero

(2.1) \begin{align} \overline {\boldsymbol{\mathsf{{f}}}}:=\frac 1T\int _0^T\boldsymbol{\mathsf{{f}}}(t)\,\textrm{ d}t=\boldsymbol{0}. \end{align}

A remarkable special case that motivated our research is when $\boldsymbol{\mathsf{{f}}}$ is the result of the periodic displacement of a mass $m$ in the interior of ${\mathscr{B}}$ . More specifically, ${\mathscr{B}}$ houses a mechanism consisting of a cavity in which $m$ slides along a prescribed constant (relative to ${\mathscr{B}}$ ) direction $\widehat {\boldsymbol{b}}$ in a time-periodic manner (see figure 1). As such, if $y=y(t)$ is the displacement of $m$ from some reference point in ${\mathscr{B}}$ , by Newton’s law of motion, this mechanism induces on ${\mathscr{B}}$ the force

(2.2) \begin{align} \boldsymbol{\mathsf{{f}}}(t)=-m\ddot {y}(t)\widehat {\boldsymbol{b}}. \end{align}

Thus, if $y(t)$ is $T$ -periodic and sufficiently smooth, we deduce $\overline {\boldsymbol{\mathsf{{f}}}}=\boldsymbol{0}$ .

Our analysis will be devoted to the general case, where $\boldsymbol{\mathsf{{f}}}$ is any (sufficiently smooth) $T$ -periodic force satisfying (2.1). The equations of motion for the coupled body–liquid system $(\mathscr{B},\mathscr{L})$ , referred to a frame attached to ${\mathscr{B}}$ , are then given by (Galdi Reference Galdi2002, Section 1)

(*) \begin{align} \left.\begin{array}{r@{\,}l} \dfrac {\partial \boldsymbol{v}}{\partial t}+(\boldsymbol{v}-\boldsymbol{\gamma })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v} & = \textrm{div}\,\boldsymbol{T}(\boldsymbol{v},p) \\ \textrm{div}\,\boldsymbol{v} &=0 \\ \end{array} \right\} &\quad \text{in}\kern2.5pt\quad \varOmega \times \mathbb{R} \nonumber \\[5pt]\boldsymbol{v}=\boldsymbol{\gamma }&\quad \text{on}\quad \partial \varOmega \times \mathbb{R} \\[5pt]\lim _{|\boldsymbol{x}|\rightarrow \infty }\boldsymbol{v}(\boldsymbol{x},\boldsymbol{\cdot })=\boldsymbol{0}&\quad \text{in}\kern2.5pt\quad \mathbb{R}. \nonumber\\[5pt]M\dot {\boldsymbol{\gamma }}=\boldsymbol{f}-\int _{\partial \varOmega }\boldsymbol{T}(\boldsymbol{v},p)\boldsymbol{\cdot }\boldsymbol{n}\,\textrm{d}S &\quad \text{in}\kern2.5pt\quad \mathbb{R}, \nonumber \end{align}

where, $\boldsymbol{v}$ and $\rho p$ are the velocity and pressure fields of ${\mathscr{L}}$ , respectively, with $\rho$ its density, whereas $\boldsymbol{T}(\boldsymbol{v},p):=-p\boldsymbol{1}+2\nu \boldsymbol{D}(\boldsymbol{v})$ , with $\boldsymbol{1}$ identity tensor, $\nu$ kinematic viscosity and where

(2.3) \begin{align} \boldsymbol{D}(\boldsymbol{v}):=\frac {1}{2}\big (\boldsymbol{\nabla }\boldsymbol{v}+\left (\boldsymbol{\nabla }\boldsymbol{v}\right )^{\top }\big ), \end{align}

is the Cauchy stress tensor. Moreover, $\rho M$ represents the mass and $\boldsymbol{\gamma }$ the (translational) velocity of ${\mathscr{B}}$ , and $\boldsymbol{n}$ stands for the outer unit normal to $\partial \varOmega$ . We also set

(2.4) \begin{align} \boldsymbol{\mathsf{{f}}}=:\rho \boldsymbol{f}. \end{align}

Thus, in view of (2.1), we have $\overline {\boldsymbol{f}}=\boldsymbol{0}$ .

As mentioned earlier, our goal is to find out when and how $\boldsymbol{f}$ induces a propulsion on ${\mathscr{B}}$ . Since $\boldsymbol{f}$ is $T$ -periodic, it is expected – and shown in Theorem 4.5 – that problem (*) admits $T$ -periodic solutions. For such solutions, denote by $\boldsymbol{\eta }=\boldsymbol{\eta }(t)$ the position vector of the centre of mass $G$ of ${\mathscr{B}}$ counted, say, from time $t=0$ . We then have

(2.5) \begin{align} \boldsymbol{\eta }(t)=\int _0^t\boldsymbol{\gamma }(s)\,\textrm{d}s+\boldsymbol{\eta }(0). \end{align}

Because $\boldsymbol{\gamma }$ is $T$ -periodic, it follows that the distance $d_T$ covered by $G$ in any interval of length $T$ is given by

(2.6) \begin{align} d_T=\left |\int _0^T\boldsymbol{\gamma }(s)\,\textrm{d}s\right | . \end{align}

As a result, ensuring propulsion occurs is equivalent to showing (1.1); namely, that the time average of $\boldsymbol{\gamma }$ over a period is non-zero.

Remark 2.1. Our propulsion problem is strictly nonlinear. In fact, suppose $(\boldsymbol{v},p,\boldsymbol{\gamma })$ is a $T$ -periodic solution to (*). Denote by $\delta \gt 0$ the magnitude of $\boldsymbol{f}$ , so that $\boldsymbol{f}=\delta \boldsymbol{F}$ . If we write $\boldsymbol{v}=\delta \boldsymbol{v}_0$ , $p=\delta p_0$ , $\boldsymbol{\gamma }=\delta \boldsymbol{\gamma }_0$ , and disregard terms in $\delta ^2$ and higher, we find that $(\boldsymbol{v}_0,p_0,\boldsymbol{\gamma }_0)$ obeys the following problem:

(2.7) \begin{align} \left.\begin{array}{r} \dfrac {\partial \boldsymbol{v}_0}{\partial t}={\rm{div}}\,\boldsymbol{T}(\boldsymbol{v}_0,p_0) \\ {\rm{div}}\,\boldsymbol{v}_0=0 \end{array}\right \}& \quad {\rm{in}\kern2.5pt\quad \varOmega \times \mathbb{R}} \nonumber \\[5pt]\boldsymbol{v}_0=\boldsymbol{\gamma }_0& \quad\rm{on}\quad \partial \varOmega \times \mathbb{R} \\[5pt]\mathop {\rm{lim}} \limits _{|\boldsymbol{x}|\rightarrow \infty }\boldsymbol{v}_0(\boldsymbol{x},\boldsymbol{\cdot })=\boldsymbol{0}& \quad \rm{in}\kern2.5pt\quad \mathbb{R}. \nonumber\\[5pt] M\dot {\boldsymbol{\gamma }}_0=\boldsymbol{F}-\int _{\partial \varOmega }\boldsymbol{T}(\boldsymbol{v}_0,p_0)\boldsymbol{\cdot }\boldsymbol{n}\,\rm{d}S & \quad {\rm{in}\kern2.5pt\quad \mathbb{R}}.\end{align}

Therefore, the averaged fields $(\overline {\boldsymbol{v}}_0,\overline {p}_0,\overline {\boldsymbol{\gamma }}_0)$ and $\overline {\boldsymbol{F}}$ obey the boundary-value problem

(2.8) \begin{align} \left .\begin{array}{r} {\rm{div}}\,\boldsymbol{T}(\overline {\boldsymbol{v}}_0,\overline {p}_0)=\boldsymbol{0} \\[3pt] \rm{div}\,\overline {\boldsymbol{v}}_0=0 \\ \end{array}\right \}& \quad{\rm{in}\kern2.5pt\quad \varOmega } \nonumber\\[5pt]\overline {\boldsymbol{v}}_0=\overline {\boldsymbol{\gamma }}_0& \quad\rm{on}\quad \partial \varOmega \\[5pt]\mathop {\rm{lim}} \limits _{|\boldsymbol{x}|\rightarrow \infty }\overline {\boldsymbol{v}}_0(\boldsymbol{x})=\boldsymbol{0}& \nonumber \\[5pt]\int _{\partial \varOmega }\boldsymbol{T}(\overline {\boldsymbol{v}}_0,\overline {p}_0)\boldsymbol{\cdot }\boldsymbol{n}\,\rm{d}S=\overline {\boldsymbol{F}} \nonumber \end{align}

From classical results on the Stokes problem (Galdi Reference Galdi2011, Section V.7), we infer that, in the very large class of weak solutions, we may have $(\overline {\boldsymbol{v}_0},\boldsymbol{\nabla }\overline {p_0},\overline {\boldsymbol{\gamma }_0})\neq (\boldsymbol{0},\boldsymbol{0},\boldsymbol{0})$ if and only if $\overline {\boldsymbol{F}}\neq \boldsymbol{0}$ ; namely, $\overline {\boldsymbol{f}}\neq \boldsymbol{0}$ . Since, in our case $\overline {\boldsymbol{f}}=\boldsymbol{0}$ , we thus conclude that propulsion can only occur at the order of $\delta ^2$ or higher; that is, only when nonlinear effects are taken into account.

3. Notation and functional setting

We begin to recall some basic notation. By $\varOmega \subset \mathbb{R}^3$ we denote the complement in $\mathbb{R}^3$ of a compact domain ${\mathscr{B}}$ of class $C^2$ (the ‘body’); see Remark 4.2. By $B_r$ we indicate the ball of radius $r\gt 0$ in $\mathbb{R}^3$ . For a domain $A\subseteq \mathbb{R}^3$ , $L^q(A)$ denotes the usual Lebesgue space endowed with the norm $\|\boldsymbol{\cdot }\|_{L^q(A)}$ and, for $m\in \mathbb{N}$ and $q\in [1,\infty ]$ , $W^{m,q}(A)$ stands for the Sobolev space with norm $\|\boldsymbol{\cdot }\|_{W^{m,q}(A)}$ . We use $(\boldsymbol{\cdot },\boldsymbol{\cdot })_{L^2(A)}$ to denote the standard inner product in $L^2(A)$ . Moreover, $D^{m,q}(A)$ will denote the homogeneous Sobolev space with semi-norm $|u|_{D^{m,q}(A)}:=\sum _{m=|\alpha |}\|D^{\alpha }u\|_{L^q(A)}$ . Where the context is clear, we shall simply write $\|\boldsymbol{\cdot }\|_q\equiv \|\boldsymbol{\cdot }\|_{L^q(A)}$ , $\|\boldsymbol{\cdot }\|_{m,q}\equiv \|\boldsymbol{\cdot }\|_{W^{m,q}(A)}$ , $|\boldsymbol{\cdot }|_{m,q}\equiv |\boldsymbol{\cdot }|_{D^{m,q}(A)}$ and $(\boldsymbol{\cdot },\boldsymbol{\cdot })\equiv (\boldsymbol{\cdot },\boldsymbol{\cdot })_{L^2(A)}$ . When any of the above function spaces are used with the subscript ‘per’, we shall mean that a function $u$ from this space has the additional property of being $T$ -periodic; namely, $u(t+T)=u(t)$ , for almost all $t\in \mathbb{R}$ . Finally, given a function $w=w(t)$ defined in the interval $(0,T)$ , we define its average

(3.1) \begin{align} \overline {w}:=\frac 1T\int _0^Tw(t)\,\textrm {d}t. \end{align}

We introduce the set

(3.2) \begin{align} \begin{aligned} \mathcal{C}(\varOmega )&:= \left \{\boldsymbol{\varphi }\in C_0^{\infty }(\overline {\varOmega }): \begin{array}{l} \text{$\textrm{div}\,\boldsymbol{\varphi }=0$ in $\varOmega $;} \\[5pt]\text{$\boldsymbol{\varphi }=\boldsymbol{\gamma }_{\boldsymbol{\varphi }}$ in a neighbourhood of $\partial \varOmega $, for some $\boldsymbol{\gamma }_{\boldsymbol{\varphi }} \in \mathbb{R}^3$} \end{array} \right \}, \end{aligned} \end{align}

together with the inner products

(3.3) \begin{align} \left (\boldsymbol{u},\boldsymbol{w}\right )_{\mathcal{K}(\varOmega )}:=\left ( \boldsymbol{u},\boldsymbol{w}\right )_{L^2(\varOmega )}+M\boldsymbol{\gamma }_{\boldsymbol{u}}\boldsymbol{\cdot }\boldsymbol{\gamma }_{\boldsymbol{w}},\quad \text{for all} \quad \boldsymbol{u},\boldsymbol{w}\in \mathcal{C}(\varOmega ), \end{align}

and (see (2.3))

(3.4) \begin{align} \left (\boldsymbol{u},\boldsymbol{w}\right )_{\mathcal{H}(\varOmega )}:=\int _{\varOmega }\boldsymbol{D}(\boldsymbol{u}):\boldsymbol{D}(\boldsymbol{w}),\quad \text{for all} \quad \boldsymbol{u},\boldsymbol{w}\in \mathcal{C}(\varOmega ), \end{align}

and associated norms,

(3.5) \begin{align} \|\boldsymbol{u}\|_{\mathcal{K}(\varOmega )}:=\left (\|\boldsymbol{u}\|^2_{L^2(\varOmega )}+M|\boldsymbol{\gamma }_{\boldsymbol{u}}|^2\right )^{1/2},\quad \text{for all} \quad \boldsymbol{u}\in \mathcal{C}(\varOmega ), \end{align}

and

(3.6) \begin{align} {\|\boldsymbol{u}\|}_{\mathcal{H}(\varOmega )}:={\|\boldsymbol{D}(\boldsymbol{u})\|}_{L^2(\varOmega )},\quad \text{for all} \quad \boldsymbol{u}\in \mathcal{C}(\varOmega ), \end{align}

respectively. We then define the completions

(3.7) \begin{align} \mathcal{K}(\varOmega ):=\overline {\mathcal{C}(\varOmega )}^{\|\boldsymbol{\cdot }\|_{\mathcal{K}(\varOmega )}}, \quad \mathcal{H}(\varOmega ):=\overline {\mathcal{C}(\varOmega )}^{\|\boldsymbol{\cdot }\|_{\mathcal{H}(\varOmega )}}. \end{align}

We shall often simply use the symbols $\mathcal{K}$ and $\mathcal{H}$ in place of $\mathcal{K}(\varOmega )$ and $\mathcal{H}(\varOmega )$ , respectively.

It can be shown that both $\mathcal{K}(\varOmega )$ and $\mathcal{H}(\varOmega )$ are Hilbert spaces when endowed with the inner products defined above. Moreover, (see (Galdi Reference Galdi2002, Lemma 4.11)),

(3.8) \begin{align} \begin{aligned} \mathcal{H}(\varOmega )&:= \left \{\boldsymbol{v}\in W^{1,2}_{\textit{loc}}(\overline {\varOmega })\cap L^6(\varOmega ): \begin{array}{l} \text{$\boldsymbol{D}(\boldsymbol{v})\in L^2(\varOmega )$, $\textrm{div}\,\boldsymbol{v}=0$, and} \\[5pt]\text{$\boldsymbol{v}=\boldsymbol{\gamma }_{\boldsymbol{v}}$, for some $\boldsymbol{\gamma }_{\boldsymbol{\boldsymbol{v}}} \in \mathbb{R}^3$ on $\partial \varOmega $} \end{array} \right \}. \end{aligned} \end{align}

For $m\in \mathbb{N}\cup \{\infty \}$ and fixed $T\gt 0$ , we introduce the test function spaces

(3.9)

where we use $\mathcal{C}^m_{\textit{per}}(\varOmega \times [0,T])$ to denote the functions of $\mathcal{C}^m_{\textit{per}}(\varOmega \times \mathbb{R})$ restricted to $[0,T]$ . Similarly, we will use $C^m_{\textit{per}}([0,T])$ to denote the functions of $C^m_{\textit{per}}(\mathbb{R})$ restricted to $[0,T]$ .

Given $q\in [1,\infty ]$ , we establish the conventions

(3.10) \begin{align} \widehat {L}^q_{\textit{per}}(\mathbb{R}):= \left\{\boldsymbol{\xi }\in L^q_{\textit{per}}(\mathbb{R}):\overline {\boldsymbol{\xi }}=\boldsymbol{0}\right\} \quad \text{and}\quad \widehat {W}_{\textit{per}}^{1,q}(\mathbb{R}):=\left\{\boldsymbol{\xi }\in \widehat {L}^q_{\textit{per}}(\mathbb{R}):\dot {\boldsymbol{\xi }}\in L^q_{\textit{per}}(\mathbb{R})\right\} , \end{align}

and, using the shorthand $X(Y)\equiv X(\mathbb{R};Y)$ for function spaces $X$ and $Y$ (e.g. $W^{m,q}_{\textit{per}}(L^r)\equiv W^{m,q}_{\textit{per}}(\mathbb{R};L^r(\varOmega ))$ , we define the space

(3.11) \begin{align} \mathcal{W}:=\left [\widehat {W}^{1,2}_{\textit{per}}\big(L^2\big)\cap \widehat {L}^2_{\textit{per}}\big(W^{2,2}\cap \mathcal{H}\big)\right ] , \end{align}

with norm

(3.12) \begin{align} \|\boldsymbol{u}\|_{\mathcal{W}}:=\|\boldsymbol{u}\|_{W^{1,2}(L^2)}+\|\boldsymbol{u}\|_{L^2(W^{2,2})}. \end{align}

4. On the existence of a class of solutions to problem (*)

The purpose of this section is to introduce a certain class of solutions to problem (*), suitable for the investigation of propulsion. To this end, set

(4.1) \begin{align} \delta := \|\boldsymbol{f}\|_{L^2(0,T)}, \quad \boldsymbol{F}:=\boldsymbol{f}/\|\boldsymbol{f}\|_{L^2(0,T)} \end{align}

and consider the following linear problem:

(4.2) \begin{align} \left .\begin{array}{r} \dfrac {\partial \boldsymbol{V}_0}{\partial t} =\textrm{div}\,\boldsymbol{T}(\boldsymbol{V}_0,p_0) \\\textrm{div}\,\boldsymbol{V}_0 =0 \end{array}\right \}& \quad{\text{in}\quad \varOmega \times \mathbb{R}} \nonumber\\[5pt] \boldsymbol{V}_0(\boldsymbol{x},t)=\boldsymbol{\xi }_0(t) &\quad \text{on}\quad \partial \varOmega \times \mathbb{R} \\[5pt] M\dot {\boldsymbol{\xi }}_0=\boldsymbol{F}-\int _{\partial \varOmega } \boldsymbol{T}(\boldsymbol{V}_0,p_0)\boldsymbol{\cdot }\boldsymbol{n}\;\textrm{d}S &\quad \text{in}\quad \mathbb{R}, \nonumber \end{align}

which is known to have a solution. In fact, we have the following (Galdi Reference Galdi2020, Lemma 5.2).

Theorem 4.1. For any $\boldsymbol{F}\in \widehat {L}_{\textit{per}}^2(\mathbb{R})$ , there exists a unique solution

(4.3) \begin{align} (\boldsymbol{V}_0,p_0,\boldsymbol{\xi }_0)\in \mathcal{W}\times \widehat {L}_{\textit{per}}^2(D^{1,2})\times \widehat {W}_{\textit{per}}^{1,2}(\mathbb{R}) \end{align}

to problem (4.2). In addition, this solution satisfies the estimate (observe that $\|\boldsymbol{F}\|_{L^2(0,T)}=1$ )

(4.4) \begin{align} \|\boldsymbol{V}_0\|_{\mathcal{W}}+\|p_0\|_{L^2(D^{1,2})}+\|\boldsymbol{\xi }_0\|_{W^{1,2}(0,T)}\leqslant C_0\,, \end{align}

with $C_0\gt 0$ depending only on $\varOmega$ and $T$ .

We shall then look for a solution to (*) ‘around’ the solution $(\boldsymbol{V}_0,p_0,\boldsymbol{\xi }_0)$ . Precisely, for each $\delta \gt 0$ , we introduce the formal decompositions of the fields $(\boldsymbol{v},p,\boldsymbol{\gamma })$ in (*)

(4.5) \begin{align} \boldsymbol{v}:=\boldsymbol{u}+\delta \boldsymbol{V}_0,\;\;\;\;\;\;\;\boldsymbol{\gamma }:=\boldsymbol{\xi }+\delta \boldsymbol{\xi }_0,\;\;\;\;\;\;\;p:=\pi +\delta p_0. \end{align}

Thus, substituting (4.5) into (*), we obtain the following nonlinear problem in the unknowns $(\boldsymbol{u},\boldsymbol{\xi },\pi )$ :

(4.6) \begin{align} \left .\begin{array}{r} \dfrac {\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u}-\boldsymbol{\xi })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} +\delta [(\boldsymbol{u}-\boldsymbol{\xi })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V}_0+ (\boldsymbol{V}_0 -\boldsymbol{\xi }_0)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} ]\\ =\textrm{div}\,\boldsymbol{T}(\boldsymbol{u},\pi )+\delta ^2\boldsymbol{F}_0 \\\textrm{div}\,\boldsymbol{u} =0 \end{array}\right \} &\quad {\text{in}\kern2.5pt\quad \varOmega \times \mathbb{R}} \nonumber\\ \boldsymbol{u}(\boldsymbol{x},t)=\boldsymbol{\xi }(t) &\quad \text{on}\quad \partial \varOmega \times \mathbb{R} \\ M\dot {\boldsymbol{\xi }}=-\int _{\partial \varOmega } \boldsymbol{T}(\boldsymbol{u},\pi )\boldsymbol{\cdot }\boldsymbol{n}\;\textrm{d}S &\quad \text{in}\kern2.5pt\quad \mathbb{R}, \nonumber \end{align}

where

(4.7) \begin{align} \boldsymbol{F}_0:=-(\boldsymbol{V}_0-\boldsymbol{\xi }_0)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V}_0. \end{align}

We shall prove existence of (4.6) in a class of weak solutions. To this end, we first dot-multiply (4.6) $_1$ by arbitrary $\boldsymbol{\varphi }\in \mathcal{C}_{\textit{per}}^1(\varOmega \times \mathbb{R})$ and integrate by parts using (4.6) $_{2,3,4}$ and also periodicity to get

(4.8) \begin{align} \int _0^T \Bigg [\left (\boldsymbol{u},\frac {\partial \boldsymbol{\varphi }}{\partial t}\right )_{\mathcal{K}} & -\delta \left [ \left ((\boldsymbol{u}-\boldsymbol{\xi })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V}_0,\boldsymbol{\varphi }\right )_{2}+ \left ((\boldsymbol{V}_0-\boldsymbol{\xi }_0)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u},\boldsymbol{\varphi }\right )_{2}\right ] \nonumber\\& - \left ((\boldsymbol{u}-\boldsymbol{\xi })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u},\boldsymbol{\varphi }\right )_{2}-2\nu \left ( \boldsymbol{u},\boldsymbol{\varphi }\right )_{\mathcal{H}}+ \delta ^2(\boldsymbol{F}_0,\boldsymbol{\varphi })_{2}\Bigg ]\textrm{d}t=0. \end{align}

With this in hand, we now can introduce the following definition.

Remarking on (i) of the above definition, we wish to note that, in fact, due to (B3), the condition $\boldsymbol{\xi }\in L^2_{\textit{per}}(\mathbb{R})$ is automatic if we take $\boldsymbol{\xi } \equiv \boldsymbol{\xi}_{u}$ . Now with this definition established, existence of weak solutions may then given by the following theorem, whose proof we provide in Appendix C.

Theorem 4.5. Let $(\boldsymbol{V}_0,\boldsymbol{\xi }_0)\in \mathcal{W}\times \widehat {W}_{\textit{per}}^{1,2}(\mathbb{R})$ . Then, there exists $\delta _0\gt 0$ such that, for all $\delta \in (0,\delta _0)$ , there is at least one corresponding $T$ -periodic weak solution $(\boldsymbol{u},\boldsymbol{\xi })$ to problem (4.6) satisfying the estimate

(4.9) \begin{align} \|\boldsymbol{u}\|_{L^2(\mathcal{H})}+\|\boldsymbol{\xi }\|_{L^2(0,T)}\leqslant c\,\delta ^2. \end{align}

5. Sufficient conditions for propulsion

We are now ready to address the problem of propulsion. Continuing the procedure we started in § 4, we now formally substitute the secondary scaling

(5.1) \begin{align} \boldsymbol{u}:=\delta ^2\boldsymbol{w},\;\;\;\;\;\;\;\pi =\delta ^2\textit{p},\;\;\;\;\;\;\;\boldsymbol{\xi }:=\delta ^2\boldsymbol{\chi }, \end{align}

into (4.6), using (4.2) and take the average over $(0,T)$ to obtain

(5.2)

Formally taking the limit in (5.2) as $\delta \rightarrow 0$ , we finally arrive at the following steady Stokes problem:

(5.3)

which is exactly problem (5.2) at second order in $\delta$ . Dot-multiplying (again, formally) equation (5.3 $)_1$ by arbitrary $\boldsymbol{\psi }\in \mathcal{H}(\varOmega )$ and integrating by parts over $\varOmega$ , as was done to obtain (4.8), we are lead to a weak formulation of (5.3), made precise by the following definition.

Now, for each $\delta \gt 0$ , thanks to Theorem 4.5, there exists a corresponding weak solution $(\boldsymbol{u}_{\delta },\boldsymbol{\xi }_{\delta })$ to the nonlinear problem (4.6). We claim that, as $\delta \rightarrow 0$ , their scaled time-averaged parts $(\overline {\boldsymbol{w}}_{\delta },\overline {\boldsymbol{\chi }}_{\delta })$ converge to the weak solution $(\boldsymbol{w}_0,\boldsymbol{\chi }_0)$ of (5.3), whose existence and uniqueness must be verified first. To this end, we begin by recalling the following result (Galdi Reference Galdi2011, Section V.4), (Happel & Brenner Reference Happel and Brenner1983, Section 5.2–5.4).

Lemma 5.3. Let $s\in (1,\infty )$ , $q\in ({3}/{2},\infty )$ and $r\in (3,\infty )$ . For each $i=1,2,3$ , there exists a unique solution

(5.5) \begin{align} \big(\boldsymbol{h}^{(i)},p^{(i)}\big)\in \!\big[D^{2,s}(\varOmega )\cap D^{1,q}(\varOmega )\cap L^r(\varOmega )\cap C^\infty (\varOmega )\!\big]\times \!\big[D^{1,s}(\varOmega )\cap L^q(\varOmega )\cap C^\infty (\varOmega )\!\big] \end{align}

to the Stokes problem

(5.6) \begin{align} \begin{aligned} \left .\begin{aligned} \displaystyle \mathrm{div}\,\boldsymbol{T}\big(\boldsymbol{h}^{(i)},p^{(i)}\big)&=\boldsymbol{0} \\ \mathrm{div}\,\boldsymbol{h}^{(i)}&=0 \\ \end{aligned}\right \} &\quad \mathrm{in} \quad \varOmega \\ \boldsymbol{h}^{(i)}=\boldsymbol{e}_i &\quad \mathrm{on}\quad \partial \varOmega . \end{aligned} \end{align}

Moreover, for $i,k\in \{1,2,3\}$ , we have that the matrix $\mathbb{K}$ , defined, component-wise, by

(5.7) \begin{align} (\mathbb{K})_{ki}:=\boldsymbol{e}_k\boldsymbol{\cdot }\int _{\partial \varOmega }\boldsymbol{T}\big(\boldsymbol{h}^{(i)},p^{(i)}\big)\boldsymbol{\cdot }\boldsymbol{n}\, \mathrm{d}S \end{align}

is both symmetric and invertible.

Before we move on to the existence of weak solutions to (5.3), we wish to make the following remark about the physical meaning of (5.6). Indeed, this system describes the flow of a viscous liquid around a body with the prescribed motion of pure translation along basis vector $\boldsymbol{e}_i$ . In turn, $(\mathbb{K})_{ki}$ is the resistance matrix and represents the $k^{\text{}}$ component of the hydrodynamic force exerted on $\partial \varOmega$ due to translational motion along the direction $\boldsymbol{e}_i$ .

Lemma 5.4. Problem (5.3) admits one and only one weak solution $(\boldsymbol{w}_0,\boldsymbol{\chi }_0)$ given by

(5.8) \begin{align} \boldsymbol{w}_0:=\sum _{i=1}^3\chi _{0i}\boldsymbol{h}^{(i)},\quad \boldsymbol{\chi }_0=\mathbb{K}^{-1}\boldsymbol{\cdot }\sum _{i=1}^3\big(\overline {\boldsymbol{F}}_0,\boldsymbol{h}^{(i)}\big)\boldsymbol{e}_i, \end{align}

where $\boldsymbol{\chi }_0=\chi _{0i}\boldsymbol{e}_i$ . Furthermore, letting also

(5.9) \begin{align} \textit{p}_0:=\sum _{i=1}^3\chi _{0i}p^{(i)}, \end{align}

we have that $(\boldsymbol{w}_0,\textit{p}_0,\boldsymbol{\chi }_0)$ solves (5.3) in the ordinary sense.

Proof. Since $\boldsymbol{h}^{(i)}\in \mathcal H(\varOmega )$ , $i=1,2,3$ , it follows that $\boldsymbol{w}_0$ and $\boldsymbol{\chi }_0$ as given in (5.8) are well defined; see Remark 5.2. Moreover, multiplying (5.6 $)_{\text{1-3}}$ by $\chi _{0i}$ and summing over $i$ , we immediately see that the choice of $(\boldsymbol{w}_0,\textit{p}_0,\boldsymbol{\chi }_0)$ , given in the theorem statement, satisfies (5.3 $)_{\text{1-3}}$ . Next, applying $\mathbb{K}$ to both sides of (5.8) $_2$ , from the definition (5.7), we also see the validity of (5.3 $)_4$ , completing the proof of existence. This solution is also clearly the unique one in its class, since after setting $\boldsymbol{\psi }\equiv \boldsymbol{w}_0$ in the homogeneous problem corresponding to (5.4) (i.e. for $\overline {\boldsymbol{F}}_0$ ), we have $\|\boldsymbol{w}_0\|_{\mathcal{H}}=0$ .

We are now equipped to prove the convergences claimed above.

Lemma 5.5. Let $\delta \in (0,\delta _0)$ , where $\delta _0$ is the positive number given in Theorem 4.5, and let $(\boldsymbol{u}_{\delta },\boldsymbol{\xi }_{\delta })$ be a weak solution to problem (4.6) corresponding to the unique pair $(\boldsymbol{V}_0,\boldsymbol{\xi }_0)$ given in Theorem 4.1. Apply the same scaling in (5.1) to $\boldsymbol{u}_{\delta }$ and $\boldsymbol{\xi }_{\delta }$

(5.10) \begin{align} \boldsymbol{u}_{\delta } :=\delta ^2\boldsymbol{w}_{\delta },\;\;\;\;\;\;\;\;\;\;\;\boldsymbol{\xi }_{\delta }:=\delta ^2\boldsymbol{\chi }_{\delta }. \end{align}

Then, as $\delta \rightarrow 0$ , we have

(5.11) \begin{align} \overline {\boldsymbol{w}}_{\delta } \rightharpoonup \boldsymbol{w}_0 & \quad \mathrm{in}\quad \mathcal{H}(\varOmega ) , \nonumber \\\overline {\boldsymbol{\chi }}_{\delta }\longrightarrow \boldsymbol{\chi }_0 & \quad \mathrm{in}\quad \mathbb{R}^3 ,\end{align}

where $(\boldsymbol{w}_0,\boldsymbol{\chi }_0)$ is the weak solution to problem (5.3) furnished by Lemma 5.4.

Proof. By the uniqueness property afforded by Lemma 5.4, it suffices to show the properties (5.11) along a subsequence, say $\{\delta _n\}_{n\in \mathbb{N}}$ , of positive numbers with $\lim _{n\rightarrow \infty }\delta _n=0$ . Given such a subsequence, for each $n\in \mathbb{N}$ , write

(5.12) \begin{align} \boldsymbol{u}_n=\delta _n^2\boldsymbol{w}_n,\;\;\;\;\;\boldsymbol{\xi }_n=\delta _n^2\boldsymbol{\chi }_n. \end{align}

Then, from (4.9) and (5.12), we easily obtain the estimate

(5.13) \begin{align} \|\overline {\boldsymbol{w}}_n\|_{\mathcal{H}}+|\overline {\boldsymbol{\chi }}_n|\leqslant \kappa , \end{align}

where, from now on, by $\kappa$ we denote a generic positive constant depending, at most, on $T$ , $\boldsymbol{V}_0$ and $\boldsymbol{\xi }_0$ . Then, by standard compactness theorems, one finds $\widetilde {\boldsymbol{w}}\in \mathcal{H}(\varOmega )$ and $\widetilde {\boldsymbol{\chi }}\in \mathbb{R}^3$ , such that (up to subsequence)

(5.14) \begin{align} \overline {\boldsymbol{w}}_n \rightharpoonup {\;\;\;\;\,} \widetilde {\boldsymbol{w}} \quad \text{in} \quad {\mathcal{H}(\varOmega )},\;\;\;\;\;\;\;\overline {\boldsymbol{\chi }}_n\longrightarrow \widetilde {\boldsymbol{\chi }}\quad \text{in} \quad \mathbb{R}^3 \end{align}

as $n\rightarrow \infty$ and also

(5.15) \begin{align} \widetilde {\boldsymbol{w}}=\widetilde {\boldsymbol{\chi }} \quad \mathrm{on} \quad \partial \varOmega . \end{align}

Next, in (4.8) we take, in particular, $\boldsymbol{\varphi }\in \mathcal{C}(\varOmega )$ and $\boldsymbol{u}\equiv \boldsymbol{u}_n$ , $\boldsymbol{\xi }\equiv \boldsymbol{\xi }_n$ as in (5.12). This yields

(5.16) \begin{align} 2\nu (\overline {\boldsymbol{w}}_n,\boldsymbol{\varphi })_{\mathcal{H}(\varOmega )}=\delta _n^2 A_n+\delta _nB_n+\big(\overline {\boldsymbol{F}}_0,\boldsymbol{\varphi }\big), \end{align}

where

(5.17) \begin{align} \begin{aligned} A_n &:=-\big(\overline {(\boldsymbol{w}_n-\boldsymbol{\chi }_n)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{w}_n},\boldsymbol{\varphi }\big), \\ B_n &:=-\big(\overline {(\boldsymbol{w}_n-\boldsymbol{\chi }_n)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V}_0},\boldsymbol{\varphi }\big)- \big(\overline {(\boldsymbol{V}_0-\boldsymbol{\xi }_0)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{w}_n},\boldsymbol{\varphi }\big). \end{aligned} \end{align}

Then, comparing (5.16) with (5.4), one immediately sees from (5.14 $)_1$ and (5.11 $)_1$ that the lemma is proved once we show that $A_n$ and $B_n$ are bounded uniformly in $n$ ; indeed, then we can pass to the limit in (5.16) as $n\rightarrow \infty$ and, with (5.15), use the uniqueness property of Lemma 5.4. But this is easily accomplished by first observing, from (4.9) and (5.12), that

(5.18) \begin{align} \|\boldsymbol{w}_n\|_{L^2(\mathcal{H})}+\|\boldsymbol{\chi }_n\|_{L^2(0,T)}\leqslant \kappa , \end{align}

and then employing in (5.17) this uniform bound (5.18) in combination with Lemma B.2 and Hölder inequality. The proof is thereby complete.

As a corollary to Lemma 5.5, we have now the main result of this section. We recall that

(5.19) \begin{align} \boldsymbol{F}_0=-(\boldsymbol{V}_0-\boldsymbol{\xi }_0)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V}_0. \end{align}

Theorem 5.6. Let $(\boldsymbol{v},\boldsymbol{\gamma })$ be a weak solution to problem (*) corresponding to the force $\boldsymbol{f}\in L_{\textit{per}}^{\infty }(\mathbb{R})$ , where $\overline {\boldsymbol{f}}=:\delta \, \overline {\boldsymbol{F}}=0$ . Then, if

(5.20) \begin{align} \boldsymbol{G}:=-\sum _{i=1}^3\left(\overline {(\boldsymbol{V}_0-\boldsymbol{\xi }_0)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V}_0},\boldsymbol{h}^{(i)}\right)\boldsymbol{e}_i \neq \boldsymbol{0}\,, \end{align}

we necessarily have $\overline {\boldsymbol{\gamma }}\neq \boldsymbol{0}$ ; that is, ${\mathscr{B}}$ experiences propulsion. Precisely, there is $\delta _0\gt 0$ such that

(5.21) \begin{align} \overline {\boldsymbol{\gamma }} = \delta ^2\mathbb{K}^{-1}\boldsymbol{\cdot }\boldsymbol{G}+\boldsymbol{{R}}(\delta ), \quad \mathrm{for\ all} \quad \delta \in (0,\delta _0) , \end{align}

where

(5.22) \begin{align} \lim _{\delta \to 0}\frac {1}{\delta ^2}\,\boldsymbol{{R}}(\delta )=0. \end{align}

Proof. Dot-multiplying both sides of (5.6) $_1$ by the difference $\overline {\boldsymbol{w}}_{\delta }-\boldsymbol{w}_0$ and integrating by parts over $\varOmega$ , using (5.3 $)_3$ , (5.10), property (i) of Definition 4.4, and (5.7), we get

(5.23) \begin{align} 2\nu \sum _{i=1}^3 \big(\overline {\boldsymbol{w}}_{\delta }-\boldsymbol{w}_0,\boldsymbol{h}^{(i)}\big)_{\mathcal{H}}\boldsymbol{e}_i=\mathbb{K}\boldsymbol{\cdot }(\overline {\boldsymbol{\chi }}_{\delta }-\boldsymbol{\chi }_0). \end{align}

Taking

(5.24) \begin{align} \boldsymbol{{R}}(\delta ):=2\nu \delta ^2\mathbb{K}^{-1}\boldsymbol{\cdot }\sum _{i=1}^3 \big(\overline {\boldsymbol{w}}_{\delta }-\boldsymbol{w}_0,\boldsymbol{h}^{(i)}\big)_{\mathcal{H}}\boldsymbol{e}_i, \end{align}

this is

(5.25) \begin{align} \boldsymbol{{R}}(\delta )=\delta ^2\big(\overline {\boldsymbol{\chi }}_{\delta }-\boldsymbol{\chi }_0\big), \end{align}

from which, thanks to (5.11) $_2$ , property (5.22) immediately follows. Then (5.21) is a consequence of applying (5.10) $_2$ , (5.8) $_2$ , (4.5) $_2$ , and the fact that $\overline {\boldsymbol{\xi }}_0=\boldsymbol{0}$ , to (5.25).

6. Numerical results

In the following section, we numerically investigate the following: if the periodic force $\boldsymbol{f}$ acting on the body ${\mathscr{B}}$ is due to an internal oscillating mass system (see figure 1), for which forces $\boldsymbol{f}$ and which shapes of ${\mathscr{B}}$ will $\boldsymbol{f}$ result in a non-zero net motion of ${\mathscr{B}}$ ? We restrict ourselves to bodies having rotational symmetry around the $z$ -axis. We also assume that $\boldsymbol{f}$ and $\boldsymbol{\gamma }$ are parallel to the $z$ -axis. Under these conditions, we seek for solutions possessing the above symmetry, so that the complete numerical set-up can be reduced to a two-dimensional setting.

We begin by analysing the functional value of $\boldsymbol{G}$ in (5.20) and the linearised periodic system (4.2) for which we must also discretise the stationary Stokes problem (5.6). Successively, we will consider the fully nonlinear coupled system (*) to explore whether a non-zero motion can be obtained even when the functional value of $\boldsymbol{G}$ is identically zero.

To simplify the numerical analysis, we shall also non-dimensionalise problems (*) and (4.2). However, before doing so, recalling the expressions (2.2) and (2.4) relating to the force, we first set

(6.1) \begin{align} \boldsymbol{f}(t):={\mathfrak{m}}\ddot {\boldsymbol{y}}(t), \end{align}

where we have ‘normalised’ the mass $m:=\rho {\mathfrak{m}}$ and set $\boldsymbol{y}:= -y\widehat {\boldsymbol{b}}$ , for convenience. Then, letting $a$ and $\omega$ denote the characteristic length of ${\mathscr{B}}$ and the frequency of oscillations, respectively, and setting $\boldsymbol{F}:={\mathfrak{m}}\ddot {\boldsymbol{Y}}$ , for $\boldsymbol{y}=\delta \boldsymbol{Y}$ , we proceed to scale

(6.2) \begin{align} \boldsymbol{x}, \boldsymbol{y} \ \, \text{and} \ \, \boldsymbol{Y} & \quad \text{by} \quad a, \nonumber \\[0pt]t & \quad \text{by} \quad 1/\omega , \nonumber \\[0pt]\boldsymbol{v}, \boldsymbol{\gamma }, \boldsymbol{\xi }_0, \ \text{and} \ \, \boldsymbol{V}_0 & \quad \text{by} \quad \omega a, \\[0pt]p \ \, \text{and} \ \, p_0 &\quad \text{by} \quad \nu \omega , \nonumber \\[0pt]M \ \, \text{and} \ \, {\mathfrak{m}} &\quad \text{by} \quad a^3, \nonumber \end{align}

to obtain the resulting dimensionless forms

(6.3) \begin{align} \left. \begin{array}{r} \displaystyle 2h^2\left [\frac {\partial \boldsymbol{v}}{\partial t}+(\boldsymbol{v}-\boldsymbol{\gamma })\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}\right ] =\textrm{div} \,\boldsymbol{T}(\boldsymbol{v}, p) \\ \textrm{div}\,\boldsymbol{v} =0 \end{array}\right \}& \quad{\text{in}\kern2.5pt\quad \varOmega \times \mathbb{R}} \nonumber \\[5pt] \boldsymbol{v}=\boldsymbol{\gamma } &\quad\text{on}\quad \partial \varOmega \times \mathbb{R} \\[5pt] \lim _{|\boldsymbol{x}|\rightarrow \infty }\boldsymbol{v}(\boldsymbol{x},t)=\boldsymbol{0} &\quad \text{in}\kern2.5pt\quad \mathbb{R}. \nonumber \\[5pt] \displaystyle 2h^2\dot {\boldsymbol{\gamma }}=2h^2\ddot {\boldsymbol{y}}-\int _{\partial \varOmega }\boldsymbol{T}(\boldsymbol{v}, p)\boldsymbol{\cdot }\boldsymbol{n}\,\textrm{d}S & \quad {\text{in}\kern2.5pt\quad \mathbb{R}}, \nonumber \end{align}

and

(6.4) \begin{align} \left. \begin{array}{r} \displaystyle 2h^2\frac {\partial \boldsymbol{V}_0}{\partial t} =\textrm{div}\,\boldsymbol{T}(\boldsymbol{V}_0, p_0) \\ \textrm{div}\boldsymbol{V}_0 =0 \end{array}\right \}& \quad {\text{in}\kern2.5pt\quad \varOmega \times \mathbb{R}} \nonumber \\[5pt] \boldsymbol{V}_0=\boldsymbol{\xi }_0 &\quad \text{on}\quad \partial \varOmega \times \mathbb{R} \\[5pt] \lim _{|\boldsymbol{x}|\rightarrow \infty }\boldsymbol{V}_0(\boldsymbol{x},t)=\boldsymbol{0}& \nonumber \\[5pt]\displaystyle 2h^2\dot {\boldsymbol{\xi }}_0=2h^2\ddot {\boldsymbol{Y}}-\int _{\partial \varOmega } \boldsymbol{T}(\boldsymbol{V}_0, p_0)\boldsymbol{\cdot }\boldsymbol{n}\,\textrm{d}S &\quad {\text{in}\kern2.5pt\quad \mathbb{R}}, \nonumber \end{align}

of (*) and (4.2), respectively, where we have set the dimensionless masses to 1 and defined the dimensionless parameter

(6.5) \begin{align} h:=a\sqrt {\frac {\omega }{2\nu }}, \end{align}

representing the, so-called, ‘Stokes number’. Finally, we similarly non-dimensionalise the auxiliary system (5.6) by scaling $\boldsymbol{h}^{(i)}$ with 1 $(=|\boldsymbol{e}_i|)$ and $p^{(i)}$ with $\nu /a$ to obtain the following dimensionless form of the vector $\boldsymbol{G}$ , defined in (5.20):

(6.6) \begin{align} \boldsymbol{G}=-2h^2\sum _{i=1}^3\left (\int _{\varOmega }\overline {(\boldsymbol{V}_0-\boldsymbol{\xi }_0)\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{V}_0}\boldsymbol{\cdot }\boldsymbol{h}^{{(i)}}\,\textrm{d}V\right )\boldsymbol{e}_i. \end{align}

6.1. Definition of the test cases

Recall that we construct the force acting on ${\mathscr{B}}$ via the internal motion of a mass. As such, we consider three test forces which we identify with the following (dimensionless) accelerations (see (6.3) $_5$ and (6.4) $_5$ ):

(6.7) \begin{align} \begin{aligned} \ddot {y}_1(t)=\ddot {Y}_1(t) &= \frac {\textrm{d}^2}{\textrm{d}t^2}\big [\sin ^2(\pi t)\big ] ,\\\ddot {y}_2(t)=\ddot {Y}_2(t) &= \frac {\textrm{d}^2}{\textrm{d}t^2}\big [\sin ^2\big(\pi t^2\big)\big ] ,\\\ddot {y}_3(t)=\ddot {Y}_3(t) &= \frac {\textrm{d}^2}{\textrm{d}t^2}\big [\sin \big(\pi t^2\big)-t\left (1-t\right )\pi \big ]. \end{aligned} \end{align}

As such, $\ddot {y}_1$ represents a ‘fully symmetric’ force with high regularity, whereas $\ddot {y}_2$ and $\ddot {y}_3$ are ‘non-symmetric’ (see figure 2). Each force is periodically extended from $[0,T)$ to $\mathbb{R}$ . As such, $\ddot {y}_1$ is $C^\infty (\mathbb{R})$ and $\ddot {y}_2$ and $\ddot {y}_3$ have anti-derivatives $\dot {y}_2$ and $\dot {y}_3$ that are continuous on $[0,T)$ (see also figure 2).

Figure 2. Three forces used in the simulation. (a) internal motion of the rigid body $y(t)$ , (b) relative velocity of the rigid body $\dot {y}(t)$ and (c) its acceleration $\ddot {y}(t)$ . From top to bottom: symmetric smooth force $\ddot {y}(t)$ , (b) non-symmetric forces $\ddot {y}_2(t)$ and (c) $\ddot {y}_3(t)$ .

Although all all forces $\ddot {y}_1,\ddot {y}_2$ and $\ddot {y}_3$ have average zero, numerical quadrature on the time grid $0=t_0\lt \ldots \lt t_N=T$ introduces an error

(6.8) \begin{align} 0 \,=\, \int _0^T \ddot {y}(t)\,\textrm{d}t\,\approx \, \sum _{i=1}^N\Delta t\boldsymbol{\cdot }\ddot {y}(t_i)\,\neq \, 0. \end{align}

For numerical simulation we therefore correct the force such that its discrete average is exactly zero independent of the step size $\Delta t\gt 0$ . Instead of $\ddot {y}(t)$ we use

(6.9) \begin{align} \ddot {y}(t)-\frac {1}{T}\sum _{i=1}^N\Delta t\boldsymbol{\cdot }\ddot {y}(t_i). \end{align}

We study three different bodies, defined as follows:

(6.10) \begin{align} \begin{aligned} \mathscr{B}_e &= \Big\{ \boldsymbol{x}\in \mathbb{R}^3,\; 1.5\big(x_1^2+x_2^2\big)+0.7 x_3^2 \lt 1 \Big\},\\[5pt]\mathscr{B}_d &= \left \{ \boldsymbol{x}\in \mathbb{R}^3,\; 1.5\frac {x_1^2+x_2^2}{(1+0.3 x_3)^2}+0.7 x_3^2 \lt 1 \right \},\\[5pt]\mathscr{B}_{fd} &= \left \{ \boldsymbol{x}\in \mathbb{R}^3,\; 1.5\frac {x_1^2+x_2^2}{(1-0.3 x_3)^2}+0.7 x_3^2 \lt 1 \right \}. \end{aligned} \end{align}

Here, ${\mathscr{B}}_e$ is an ellipsoid, ${\mathscr{B}}_d$ is drop-like shape, and ${\mathscr{B}}_{\textit{fd}}$ is the body ${\mathscr{B}}_d$ in the opposite orientation with respect to the axis of symmetry (see figure 3).

Figure 3. Ellipsoid with symmetry in the direction of flow (a), drop-like shape (b) and flipped drop (c).

6.2. Discretisation

The discretisation of (6.3), (6.4) and (5.6) is fairly standard, but some care is required, as we are looking for periodic solutions of a system that is expected to undergo oscillations of substantial amplitude with a mean velocity that is zero in the linear and close to zero (compared with its amplitude) in the nonlinear case. Even small absolute errors in the amplitude of the oscillatory solution might have a large impact on the resulting average velocity which then gives rise to a substantial relative error.

To discretise in time, we use the second-order backward difference formula BDF2 for both the rigid body system and for the Navier–Stokes equation. It is sufficiently accurate, strongly A-stable but it introduces little numerical damping. Special care is taken in the last time step of each interval, e.g. $(0,T)$ , $(T,2T)$ , etc., as the forces $\ddot y_2$ and $\dot y_3$ are not globally continuous. Here, we make certain that this force is always evaluated from within the periodic interval.

We split the period $T$ into $N$ discrete time steps of size $\Delta t=T/N$ and solve the coupled system of Stokes and rigid body motion (6.4) in a semi-explicit iteration. We start with

(6.11) \begin{align} \boldsymbol{\gamma }_n^{(0)} = \boldsymbol{\gamma }_{n-1},\quad \boldsymbol{v}_n^{(0)} = \boldsymbol{v}_{n-1} \end{align}

and iterate for $l=1,2,\ldots$ with a relaxation parameter $\omega \in (0,1]$ which is usually set to $\omega =0.8$

\begin{align*} \dfrac {\dfrac {3}{2}\widetilde {\boldsymbol{\gamma }}_n^{(l)}-2\boldsymbol{\gamma }_{n-1}+\dfrac {1}{2}\boldsymbol{\gamma }_{n-2}}{\Delta t}&=\boldsymbol{f}(t_n) -\boldsymbol{e}_z\int _{\partial \varOmega }\boldsymbol{T}\big (\boldsymbol{v}_{n}^{(l-1)},p_n^{(l-1)}\big )\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_z\,\textrm{d}s, \nonumber\\\boldsymbol{\gamma }_n^{(l)} &= \omega \widetilde {\boldsymbol{\gamma }}_n^{(l)}+(1-\omega )\boldsymbol{\gamma }_n^{(l-1)},\end{align*}

(6.12) \begin{align}\dfrac {\dfrac {3}{2}\boldsymbol{v}_n^{(l)}-2\boldsymbol{v}_{n-1}+\dfrac {1}{2}\boldsymbol{v}_{n-2}}{\Delta t}&= \textrm{div}\, \boldsymbol{T}\big (\boldsymbol{v}_n^{(l)},p_n^{(l)}\big )\qquad {\boldsymbol{v}}_n^{(l)}=\boldsymbol{\gamma }_n^{(l)}\text{ on}\quad \partial \mathscr{B}, \nonumber \\\textrm{div}\, \boldsymbol{v}_n^{(l)} &=0. \end{align}

By $\widetilde {\boldsymbol{\gamma }}_n^{(l)}$ we denote an intermediate solution that is damped by $\omega$ for better robustness and convergence. Since we only consider set-ups with cylindrical symmetric and motion only along the axis of symmetry, the Stokes problem is reformulated in cylindrical coordinates reduced to the $(r,z)$ -plane. For numerical simulation, we have to restrict the domain $\varOmega =\mathbb{R}^3\setminus \overline {\mathscr{B}}$ by introducing, for each $R\gt \text{diam}\,\mathscr{B}$ , the numerical domain $D_{R}$ as

(6.13) \begin{align} D_{R} = \Big \{ (r,z)\in \mathbb{R}^2\,:\, r\in (0,R),\; z\in (-R,R)\} \setminus \overline {\mathscr{B}}\Big \}, \end{align}

where the rigid body ${\mathscr{B}}$ is centred at $(0,0)$ . We take $R=16$ . This choice has been shown to have minimal impact on the artificial domain restriction.

Figure 4. Computational domain $D_R$ with indication of the boundaries (a) and part of the mesh enclosing the flipped drop ${\mathscr{B}}_{\textit{fd}}$ (b). The largest edges have the length $\Delta x=0.5$ whereas $\Delta x=0.0073$ at $\partial \mathscr{B}_{fd}$ .

On the symmetry boundary at $r=0$ we impose the Dirichlet condition $v_r=0$ for the radial component. On the outer boundary at $r=R$ a no-slip condition $\boldsymbol{v}=\boldsymbol{0}$ is set. On the surface of the body $\partial \mathscr{B}$ the fluid’s velocity matches that of the solid (i.e. $\boldsymbol{v}=\boldsymbol{\gamma }$ ).

For the linear test cases aiming at the functional values, the do-nothing outflow condition $\partial _n \boldsymbol{v}-p\boldsymbol{n}=\boldsymbol{0}$ holds, see Heywood, Rannacher & Turek (Reference Heywood, Rannacher and Turek1992) and is imposed at $z=R$ and $z=-R$ . This condition, however, is not compatible with the transformed body-centred reference system of the nonlinear problem as $\boldsymbol{v}$ is not zero on the numerical boundary.

The surface force $\boldsymbol{T}(\boldsymbol{v},p)\boldsymbol{\cdot }\boldsymbol{n}$ on the body is evaluated as a volume integral

(6.14) \begin{align} \int _{\partial \varOmega }\boldsymbol{T}(\boldsymbol{v},p)\boldsymbol{\cdot }\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_z\,\textrm{d}s = \int _\varOmega \frac {d\boldsymbol{v}}{dt}\boldsymbol{\cdot }\boldsymbol{\overline {z}}+\boldsymbol{T}(\boldsymbol{v},p):\boldsymbol{\nabla }\boldsymbol{\overline {z}}\,\textrm{d}x, \end{align}

where $\boldsymbol{\overline {z}}\in H^1(\varOmega )$ with $\boldsymbol{\overline {z}} = \boldsymbol{e}_z$ on $\partial \mathscr{B}$ and $\boldsymbol{\overline {z}}=\boldsymbol{0}$ on the other boundaries. In comparison with the evaluation of the surface integral, this approach gives higher-order convergence, see Braack & Richter (Reference Braack and Richter2006).

In space we discretise with quadratic finite elements for velocity and pressure using the local projection stabilisation technique to cope with the missing inf-sup stability, see Becker & Braack (Reference Becker and Braack2001). The nonlinear algebraic problems are solved by Newton’s method and the linear systems by Generalized minimal residual method, using a parallel geometric multigrid iteration as preconditioner. Details on the discretisation are given in (Richter Reference Richter2017, Chapter 4) and the nonlinear and linear solver is detailed in Failer & Richter (Reference Failer and Richter2020). Figure 4 shows a sketch of the computational domain as well as the relevant section of the mesh enclosing the body. The coarse mesh is constructed by hand and then successively refined. New nodes on $\partial \mathscr{B}$ are pulled onto the curved surface. The mesh has an average mesh size $\Delta x=0.25$ and this size varies from $\Delta x\approx 0.0073$ at the body to $\Delta x\approx 0.5$ on the outer boundaries.

Table 1. Computed values of $\boldsymbol{G}$ for Stokes number $h=4$ depending on the shape of the domain (elliptic ${\mathscr{B}}_e$ or drop-like ${\mathscr{B}}_d$ and ${\mathscr{B}}_{\textit{fd}}$ ) and the type of force (symmetric $\ddot {y}_1$ , non-symmetric $\ddot {y}_2$ and $\ddot {y}_3$ ). As expected, the values of $\boldsymbol G$ for ${\mathscr{B}}_{d}$ and ${\mathscr{B}}_{\textit{fd}}$ are opposite to each other.

In dealing with the linear Stokes problem, we know that the resulting periodic solution of (6.4) will result in zero average velocities $\boldsymbol{\overline {v}}=\boldsymbol{0}$ and $\boldsymbol{\overline {\gamma }}=\boldsymbol{0}$ . Hence, to accelerate convergence towards the periodic solution we use the following simple algorithm that aims at projecting the average to the expected one (see Richter (Reference Richter2021) for a detailed description of the algorithm and an extension to the time-periodic Navier–Stokes problem). Starting with $\boldsymbol{v}_0:=\boldsymbol{0}$ and $\boldsymbol{\gamma }_0:=\boldsymbol{0}$ we iterate for $l=1,2,\ldots$ .

1. (i) Approximate the coupled Stokes problem on $I=[0,T]$ using the initial values

   (6.15) \begin{align} \boldsymbol{v}^{(l)}(0) = \boldsymbol{v}^{(l-1)}_0,\quad \boldsymbol{\gamma }^{(l)}(0) = \boldsymbol{\gamma }^{(l-1)}_0. \end{align}

2. (ii) If periodicity is reached up a given tolerance $\epsilon _P\gt 0$ (we choose $\epsilon _P=10^{-7}$ if not otherwise stated), namely

   (6.16) \begin{align} \|\boldsymbol{v}^{(l)}(T)-\boldsymbol{v}^{(l)}(0)\|+ \|\boldsymbol{\gamma }^{(l)}(T)-\boldsymbol{\gamma }^{(l)}(0)\| \lt \epsilon _P, \end{align}
   stop and accept $\boldsymbol{v}^{(l)},\boldsymbol{\gamma }^{(l)}$ .

3. (iii) Correct the average and predict new initial value

   (6.17) \begin{align} \boldsymbol{v}^{(l)}_0:=\boldsymbol{v}^{(l)}(T) - \frac {1}{T}\int _0^T \boldsymbol{v}^{(l)}(t)\,\textrm{d}t,\quad \boldsymbol{\gamma }^{(l)}_0:=\boldsymbol{\gamma }^{(l)}(T) - \frac {1}{T}\int _0^T \boldsymbol{\gamma }^{(l)}(t)\,\textrm{d}t. \end{align}

This iteration quickly converges to the periodic solution and the periodic tolerance is usually reached within a few cycles. For the nonlinear coupled problem, we apply this algorithm in the very first step to correct the initial values.

6.3. Computation of the thrust vector $\boldsymbol{G}$

We study the linearised problem (6.4). In table 1, we indicate the functional values of $\boldsymbol{G}$ , as given in (6.6), for the three types of forces and three shapes. The functional value is non-zero only in the case of the non-symmetric bodies ${\mathscr{B}}_d$ and ${\mathscr{B}}_{\textit{fd}}$ . The combination of the symmetric force together with a non-symmetric body does generate a non-zero functional value. The value of $\boldsymbol{G}$ changes its sign for the mirrored body ${\mathscr{B}}_{\textit{fd}}$ . It is noteworthy that the three forces $\ddot {y}_1$ , $\ddot {y}_2$ and $\ddot {y}_3$ do not give substantially different values.

Table 2. Top: dependency of $\overline {\boldsymbol{\gamma }}_0\equiv {\mathbb{K}}^{-1}\boldsymbol{\cdot }\boldsymbol{G}$ on the Stokes number for the flipped drop $B_{fd}$ with $\ddot {y}_2$ . Bottom: visualisation of the dependency compared with the average velocity $\overline {\boldsymbol{\gamma }}$ of the corresponding full nonlinear Navier–Stokes problem.

Table 3. Average body velocity $\boldsymbol{\overline {\gamma }}$ computed from the full nonlinear problem for Stokes number $h=4$ and the three forces and shapes. It is interesting to observe that the absolute values of $\overline {\boldsymbol{\gamma }}$ for ${\mathscr{B}}_{d}$ are larger than those for ${\mathscr{B}}_{\textit{fd}}$ , indicating, as expected, that the viscous drag on ${\mathscr{B}}_{\textit{fd}}$ is larger than the one on ${\mathscr{B}}_{d}$ .

Next, we study the dependence on the Stokes number $h$ of the average body velocity $\overline {\boldsymbol{\gamma }}_0:=\mathbb{K}^{-1}\boldsymbol{\cdot }\boldsymbol{G}$ determined in Theorem 5.6. We limit the study to the case of the flipped drop-shaped body ${\mathscr{B}}_{\textit{fd}}$ forced by $\ddot {y}_2$ . Numerical values and their corresponding plot are reported at the top and left of table 2. These results present interesting features. In the first place, $\overline {\boldsymbol{\gamma }}_0$ has a quadratic behaviour for sufficiently large values of $h$ , which implies that, for such $h$ , the speed of the body is proportional to the frequency of oscillations. Moreover, for small $h$ , there is a ‘critical’ value of $h$ , between $h=3$ and $h=4$ , around which $\overline {\boldsymbol{\gamma }}_0$ increases slower before growing quadratically. For comparison we show the resulting average velocity $\overline {\boldsymbol{\gamma }}$ of the full nonlinear problem. These values are larger than the corresponding results for $\overline {\boldsymbol{\gamma }}_0$ but they match its trajectory for larger values of $h$ .

Table 4. Dependency on the Stokes number $h$ of the average body velocity $\boldsymbol{\overline {\gamma }}$ computed from the full nonlinear problem, for the three different shapes and forces.

6.4. Numerical integration of the nonlinear coupled problem

The second-order analysis shows that the non-symmetry of the body yields a non-zero functional value of thrust vector (6.6) implying a non-zero net motion of the body. For the (round) ellipsoid, however, the functional output is zero, regardless of the choice of force. In order to gain insight into the role of symmetry breaking of body and force, we provide numerical results for the full nonlinear coupled problem.

All tests are run at Stokes number $h=4$ . We consider the same nine combinations tested previously for $\boldsymbol{G}$ in table 1 and report the values of the average velocity of the body $\boldsymbol{\overline {\gamma }}$ resulting, this time, directly from the numerical integration of the coupled nonlinear problem (*); the results are given in table 3. They show, in particular, that, in contrast to the second-order prediction, a non-zero net motion occurs in the case of the ellipsoid whenever the force is not symmetric. In the case of the symmetric force $\ddot {y}_1$ , the resulting average velocity for the drop and the flipped drop have opposite sign and the same magnitude up to the numerical discretisation error. This is expected due to the symmetry of the problem set-up. Studying the two non-symmetric forces $\ddot {y}_2$ and $\ddot {y}_3$ reveals the distinct influence of the symmetry of force and shape on the average velocity. In both cases, the mean of the velocities of drop and flipped drop is close to the velocity of the ellipsoid which can be considered as the shape resulting from overlaying ${\mathscr{B}}_d$ and ${\mathscr{B}}_{\textit{fd}}$ .

In table 4, we report the resulting average body velocity $\boldsymbol{\overline {\gamma }}$ depending on the Stokes number $h$ . We find approximately quadratic growth of the average velocity with increasing Stokes number. This is the same behaviour observed for the thrust vector $\boldsymbol{G}$ .

Finally, in figure 5 we report, in the case of the flipped drop and for different forces, the variation with time of the position of the centre of mass along with its net distance travelled, computed from the full nonlinear problem.

Figure 5. Position of the centre of mass of the flipped drop ${\mathscr{B}}_{\textit{fd}}$ ( $\int _0^t\boldsymbol{\gamma } (s)\,\textrm{d}s$ ) and net distance covered $(\overline {\boldsymbol{\gamma }}\,t$ ) vs. time, computed from the full nonlinear problem.

6.5. Robustness of the discretisation

To validate the robustness of the numerical discretisation we show in table 5 a convergence analysis. We consider the flipped drop ${\mathscr{B}}_{\textit{fd}}$ forced with $\ddot y_1$ . We start from the coarsest discretisation with $\Delta t\approx 0.067$ and $\Delta x =0.25$ (in average) which is the one that has been used for all studies of this work. For comparison we further consider all combinations of two further uniform refinements in space and time. The finest discretisation uses the time step size $\Delta t= ({1}/{60})\approx 0.0167$ and the average mesh size $\Delta x=0.0625$ with elements ranging from $\Delta x=0.125$ to $\Delta x\approx 0.0018$ at the body.

Table 5. Resulting average velocity $\bar {\boldsymbol{\gamma }}$ of the flipped frop ${\mathscr{B}}_{\textit{fd}}$ for the nonlinear problem forced with $\ddot {y}_1$ . We also indicate the experimentelly observed order of convergence for $\Delta x\to 0$ and $\Delta t\to 0$ . Here, EOC-estimated order of convergence.

The results show that the original discretisation has been chosen adequately, as the numerical errors are small and do not dominate the results. The experimental orders of convergence show a slightly suboptimal performance in space (1.39 vs. the expected second order) and a better than expected result in time (2.45 vs. the expected second order). However, it appears that the convergence is ‘from below’ in space and ‘from above’ in time such that cancellations take place which might distort the identified order of convergence.