The Lagrange-Galerkin method for fluid-structure interaction problems

Abstract

In this paper, we consider a Lagrange-Galerkin scheme to approximate a two-dimensional fluid-structure interaction problem. The equations of the system are the Navier-Stokes equations in the fluid part, coupled with ordinary differential equations for the dynamics of the solid. We are interested in studying numerical schemes based on the use of the characteristics method for rigid and deformable solids. The schemes are based on a global weak formulation involving only terms defined on the whole fluid-solid domain. Convergence results are stated for both semi and fully discrete schemes. This article reviews known results for rigid solid along with some new results on deformable structure yet to be published.

1 Introduction

In this article, we present a modified characteristics method for the discretization of the equations modelling the motion of a solid immersed in a cavity filled by a viscous incompressible fluid. We are interested in rigid and deformable solids modelling some particulate flows in the case of rigid solid and the swimming of slender, neutrally buoyant fish, for the deformable structure (see [1]). The presented methods are generalizations of the numerical scheme introduced in [2], where the solid immersed in the fluid is rigid and has the same density with the fluid.

The fluid-structure interaction problem that we study is characterized by the strong coupling between the nonlinear equations of the fluid and those of the structure, as well as the fact that the equations of the fluid are written in a variable domain in time, which depends on the displacement of the structure. From the numerical point of view, in this kind of problems it is necessary to solve equations on moving domains. For this reason, in recent years various authors have proposed a number of different techniques [3–9].

For the numerical treatment of convection term in the Navier-Stokes equations, we discretize the material derivative along trajectories (see [10]) combined with the Lagrange-Galerkin mixed finite element approximation of Navier-Stokes equations in a velocity/pressure formulation studied in [11]. In [12], the convergence analysis of a finite element projection/Lagrange-Galerkin method for the incompressible Navier-Stokes equations is done.

The numerical analysis of some time decoupling algorithms in the case, where the deformation of the structure induces an evolution in the fluid domain has been developed in [13] (one-dimensional problem). For the ALE method applied to interaction problems, we may cite [14] in the case of the unsteady Stokes equations in a time dependent domain and [15] for a two-dimensional problem describing the motion of a rigid body in a viscous fluid. In [2, 16], the authors have introduced a convergent numerical method based on finite elements with a fixed mesh for a two-dimensional fluid-rigid body problem, where the densities of the fluid and the solid are equal. In [17, 18], we have introduced crucial modifications on the characteristic function, and we have proposed a convergent numerical scheme for a two dimensional fluid-rigid body problem where the densities of the fluid and the solid are different. In this paper, we go further, and we present a new characteristic function which gives us convergent algorithms for the simulation of aquatic organisms (for the existence and regularity of the solution in this kind of interactions, see [19]).

2 Setting of the problem

2.1 Notation and hypothesis

Let us now introduce some notation following paper [19], where the existence and uniqueness for the solution of similar problem are treated. We denote by ${\mathcal{S}}_0$ the domain occupied by the solid in a reference configuration. We assume that ${\mathcal{S}}_0$ is an open connected set with $C^{\mathrm{\infty }}$ boundary, and we choose a system of coordinates with the origin at the mass center of ${\mathcal{S}}_0$. For the deformable solid, we suppose that the motion is given by a smooth mapping

$$\mathbf{X}:{\mathcal{S}}_0\times [0,\mathrm{\infty })\to {\mathbb{R}}^2,$$

which satisfies

$$\mathbf{X}(\mathbf{y},t)=\mathit{\xi }(t)+{\mathcal{R}}_{\theta (t)}{\mathbf{X}}^{\ast }(\mathbf{y},t)\mathrm{\forall }\mathbf{y}\in {\mathcal{S}}_0,t\ge 0,$$

(1)

where for every $t\ge 0$, $\mathit{\xi }(t)$ is the trajectory of the mass center, and $\theta (t)$ represents the angle giving the orientation of the solid. ${\mathcal{R}}_{\theta }$ denotes the matrix associated to the rotation of angle θ. In (1), ${\mathbf{X}}^{\ast }$ denotes an appropriate smooth mapping, representing the undulatory deformation of the creature. The rigid solid is obtained by considering the map

$${\mathbf{X}}^{\ast }=\mathrm{I},$$

where I denotes the identity map.

Throughout this paper, the deformable body will be called the creature or sometimes just the body in particular, when considering the rigid body case.

In the remaining part of this work, the functions ξ, θ are unknowns to be determined from the governing equations below, whereas the undulatory motion ${\mathbf{X}}^{\ast }$ will be supposed to be known and to satisfy several assumptions as in [19], which will be recalled in the sequel.

(H1) For every $t\ge 0$, the mapping $\mathbf{y}\mapsto {\mathbf{X}}^{\ast }(\mathbf{y},t)$ is a $C^{\mathrm{\infty }}$ diffeomorphism from $\overline{{\mathcal{S}}_0}$ onto $\overline{{\mathcal{S}}^{\ast }}(t)$, where ${\mathcal{S}}^{\ast }(t)={\mathbf{X}}^{\ast }({\mathcal{S}}_0,t)$. Moreover, for every $\mathbf{y}\in {\mathcal{S}}_0$, the mapping $t\mapsto {\mathbf{X}}^{\ast }(\mathbf{y},t)$ is of class $C^{\mathrm{\infty }}$ and ${\mathbf{X}}^{\ast }(\mathbf{y},0)=\mathbf{y}$.

For every $t\ge 0$ we denote by ${\mathbf{Y}}^{\ast }$ the inverse of ${\mathbf{X}}^{\ast }$, i.e., the diffeomorphism satisfying

$${\mathbf{X}}^{\ast }({\mathbf{Y}}^{\ast }({\mathbf{x}}^{\ast },t),t)={\mathbf{x}}^{\ast },{\mathbf{Y}}^{\ast }({\mathbf{X}}^{\ast }(\mathbf{y},t),t)=\mathbf{y}$$

for every ${\mathbf{x}}^{\ast }\in {\mathcal{S}}^{\ast }(t)$, $t\ge 0$ and $\mathbf{y}\in {\mathcal{S}}_0$.

(H2) The total volume of the creature is preserved, i.e.,

$${\int }_{{\mathcal{S}}^{\ast }(t)}d{\mathbf{x}}^{\ast }={\int }_{{\mathcal{S}}_0}d\mathbf{y}\mathrm{\forall }t\ge 0.$$

(2)

Denote by ${\mathbf{w}}^{\ast }$ the undulatory velocity of creature, written as a vector field on ${\mathcal{S}}^{\ast }(t)$, i.e.,

$${\mathbf{w}}^{\ast }({\mathbf{x}}^{\ast },t)={\frac{\partial {\mathbf{X}}^{\ast }}{\partial t}(\mathbf{y},t)|}_{\mathbf{y}={\mathbf{Y}}^{\ast }({\mathbf{x}}^{\ast },t)}\mathrm{\forall }{\mathbf{x}}^{\ast }\in {\mathcal{S}}^{\ast }(t),t\ge 0.$$

(3)

Let ${\rho }_{S,0}$ be the density field of the solid in the reference configuration ${\mathcal{S}}_0$, and let ${\rho }_S^{\ast }(t)$ be the density field of ${\mathcal{S}}^{\ast }(t)$. The mass conservation principle applied to the whole body gives

$${\int }_{{\mathcal{S}}^{\ast }(t)}{\rho }_S^{\ast }({\mathbf{x}}^{\ast },t)d{\mathbf{x}}^{\ast }={\int }_{{\mathcal{S}}_0}{\rho }_{S,0}(\mathbf{y})d\mathbf{y}\mathrm{\forall }t\ge 0,$$

(4)

whereas the local form of the conservation of mass yields

$${\rho }_S^{\ast }({\mathbf{x}}^{\ast },t)=\frac{{\rho }_{S,0}({\mathbf{Y}}^{\ast }({\mathbf{x}}^{\ast },t))}{det(\mathrm{\nabla }{\mathbf{X}}^{\ast })({\mathbf{Y}}^{\ast }({\mathbf{x}}^{\ast },t))}\mathrm{\forall }t\ge 0,\mathrm{\forall }{\mathbf{x}}^{\ast }\in {\mathcal{S}}^{\ast }(t),$$

(5)

where $\mathrm{\nabla }{\mathbf{X}}^{\ast }$ stands for the Jacobian matrix of ${\mathbf{X}}^{\ast }(\cdot ,t)$.

(H3) ${\int }_{{\mathcal{S}}^{\ast }(t)}{\rho }_S^{\ast }({\mathbf{x}}^{\ast },t){\mathbf{w}}^{\ast }({\mathbf{x}}^{\ast },t)d{\mathbf{x}}^{\ast }=0$ for all $t\ge 0$.

(H4) ${\int }_{{\mathcal{S}}^{\ast }(t)}{\rho }_S^{\ast }({\mathbf{x}}^{\ast },t){\mathbf{x}}^{\ast \perp }\cdot {\mathbf{w}}^{\ast }({\mathbf{x}}^{\ast },t)d{\mathbf{x}}^{\ast }=0$ for all $t\ge 0$, where we denote by ${\mathbf{x}}^{\perp }$ the vector ${\mathbf{x}}^{\perp }=\left(\begin{array}{c}-x_2\\ x_1\end{array}\right)$ for $\mathbf{x}=\left(\begin{array}{c}{x}_1\\ x_2\end{array}\right)$.

Conditions (H3), (H4) correspond to the so-called self-propelling conditions which are natural requirements for understanding swimming viewed as a self-propelled phenomena.

In particular, hypotheses (H1) and (H3) imply that the position of the center of mass of the creature is not affected by the undulatory motion, that is,

$${\int }_{{\mathcal{S}}_0}{\rho }_{S,0}(\mathbf{y}){\mathbf{X}}^{\ast }(\mathbf{y},t)d\mathbf{y}=0\mathrm{\forall }t\ge 0.$$

(6)

From (1), it follows that the region occupied by the creature at time t is given by

$$\mathcal{S}(\mathit{\xi }(t),\theta (t),t)={\mathcal{R}}_{\theta (t)}{\mathcal{S}}^{\ast }(t)+\mathit{\xi }(t)\mathrm{\forall }t\ge 0.$$

(7)

Moreover, by differentiating equation (1) with respect to t, it follows that the Eulerian velocity field of the solid is given for every $t\ge 0$ by

$${\mathbf{u}}_S(\mathbf{x},t)={\mathit{\xi }}^{\prime }(t)+{\theta }^{\prime }(t){(\mathbf{x}-\mathit{\xi }(t))}^{\perp }+\mathbf{w}(\mathbf{x},t)\mathrm{\forall }\mathbf{x}\in \mathcal{S}(\mathit{\xi }(t),\theta (t),t),$$

(8)

where

$$\mathbf{w}(\mathbf{x},t)={\mathcal{R}}_{\theta (t)}{\mathbf{w}}^{\ast }({\mathcal{R}}_{-\theta (t)}(\mathbf{x}-\mathit{\xi }(t)),t)\mathrm{\forall }\mathbf{x}\in \mathcal{S}(\mathit{\xi }(t),\theta (t),t).$$

(9)

The Eulerian density field of the body is given by

$${\rho }_S(\mathbf{x},t)={\rho }_S^{\ast }({\mathcal{R}}_{-\theta (t)}(\mathbf{x}-\mathit{\xi }(t)),t)\mathrm{\forall }t\ge 0,\mathbf{x}\in \mathcal{S}(\mathit{\xi }(t),\theta (t),t),$$

(10)

with ${\rho }_S^{\ast }$ given by (5). The mass M of the body and its moment of inertia with respect to an axis orthogonal to the plane of the motion and passing by the mass center of $\mathcal{S}(\mathit{\xi }(t),\theta (t),t)$, are as usually given by

$$M={\int }_{\mathcal{S}(\mathit{\xi }(t),\theta (t),t)}{\rho }_S(\mathbf{x},t)d\mathbf{x}.$$

(11)

$$I(t)={\int }_{\mathcal{S}(\mathit{\xi }(t),\theta (t),t)}{\rho }_S(\mathbf{x},t){|\mathbf{x}-\mathit{\xi }(t)|}^{2}d\mathbf{x}.$$

(12)

Let us notice that from (4), (10) and (11), we have that

$$M={\int }_{{\mathcal{S}}_0}{\rho }_{S,0}(\mathbf{y})d\mathbf{y}.$$

Remark 2.1 In the case ${\mathbf{X}}^{\ast }=\mathrm{I}$ (rigid solid), all hypotheses (H1)-(H4) are satisfied, and the undulatory velocity field ${\mathbf{w}}^{\ast }$ is equal to zero.

2.2 Equations

Let Ω be an open bounded set in ${\mathbb{R}}^2$ representing the domain occupied by the solid-fluid system. Recalling that $\mathcal{S}(\mathit{\xi }(t),\theta (t),t)$ is the domain occupied by the solid at instant t, we have that the fluid fills, at instant t, the domain $\mathcal{F}(\mathit{\xi }(t),\theta (t),t)=\mathrm{\Omega }\setminus \overline{\mathcal{S}(\mathit{\xi }(t),\theta (t),t)}$.

With the notation above, the full system describing the self-propelled motion of the creature can be written as

$${\rho }_F(\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \mathrm{\nabla })\mathbf{u})-\mu \mathrm{\Delta }\mathbf{u}+\mathrm{\nabla }p={\rho }_F\mathbf{f},\mathbf{x}\in \mathcal{F}(\mathit{\xi }(t),\theta (t),t),t\in (0,T),$$

(13)

$$div\mathbf{u}=0,\mathbf{x}\in \mathcal{F}(\mathit{\xi }(t),\theta (t),t),t\in (0,T),$$

(14)

$$\mathbf{u}=0,\mathbf{x}\in \partial \mathrm{\Omega },t\in (0,T),$$

(15)

$$\mathbf{u}(\mathbf{x},t)={\mathit{\xi }}^{\prime }(t)+{\theta }^{\prime }(t){(\mathbf{x}-\mathit{\xi }(t))}^{\perp }+\mathbf{w}(\mathbf{x},t),\mathbf{x}\in \partial \mathcal{S}(\mathit{\xi }(t),\theta (t),t),t\in (0,T),$$

(16)

$$M{\mathit{\xi }}^{″}(t)=-{\int }_{\partial \mathcal{S}(\mathit{\xi }(t),\theta (t),t)}\mathit{\sigma }\mathbf{n}d\mathrm{\Gamma }+{\int }_{\mathcal{S}(\mathit{\xi }(t),\theta (t),t)}{\rho }_S(\mathbf{x},t)\mathbf{f}(\mathbf{x},t)d\mathbf{x},t\in (0,T),$$

(17)

$$\begin{array}{c}{\left(I{\theta }^{\prime }\right)}^{\prime }(t)=-{\int }_{\partial \mathcal{S}(\mathit{\xi }(t),\theta (t),t)}{(\mathbf{x}-\mathit{\xi }(t))}^{\perp }\cdot \mathit{\sigma }\mathbf{n}d\mathrm{\Gamma }\hfill \\ \phantom{{\left(I{\theta }^{\prime }\right)}^{\prime }(t)=}+{\int }_{\mathcal{S}(\mathit{\xi }(t),\theta (t),t)}{\rho }_S(\mathbf{x},t){(\mathbf{x}-\mathit{\xi }(t))}^{\perp }\cdot \mathbf{f}(\mathbf{x},t)d\mathbf{x},t\in (0,T).\hfill \end{array}$$

(18)

In the system above, ${\rho }_F>0$ and $\mu >0$ stand for the density and the viscosity of the fluid, which are supposed to be constant, u is the Eulerian velocity field of the fluid, and p denotes the pressure field of the fluid. A prime stands for the derivation operator with respect to time. By using the classical notation

$$\mathbf{D}(\mathbf{u})=\frac{1}{2}((\mathrm{\nabla }\mathbf{u})+{(\mathrm{\nabla }\mathbf{u})}^T),$$

(19)

the stress tensor field σ is defined by

$$\sigma (\mathbf{u},p)=2\mu \mathbf{D}(\mathbf{u})-p\mathrm{Id},$$

(20)

where Id is the identity matrix in $M_2(\mathbb{R})$. Moreover, for $t\in [0,T]$ and $\mathbf{x}\in \partial \mathcal{S}(\mathit{\xi }(t),\theta (t),t)$ we denote by $n(\mathbf{x},t)$ the unit normal to $\mathcal{S}(\mathit{\xi }(t),\theta (t),t)$ oriented towards the solid. Recall that the mass M and the moment of inertia $I(t)$ of the solid at instant t are defined by (11) and (12).

System (13)-(20) is completed by the initial conditions

$$\mathbf{u}(\mathbf{x},0)={\mathbf{u}}_0(\mathbf{x}),\mathbf{x}\in \mathcal{F}(\mathit{\xi }(0),\theta (0),0),$$

(21)

$$\mathit{\xi }(0)={\mathit{\xi }}_0,\theta (0)={\theta }_0,{\mathit{\xi }}^{\prime }(0)={\mathit{\xi }}_1,{\theta }^{\prime }(0)={\omega }_0.$$

(22)

Remark 2.2 In the case of rigid solid, equation (16) becomes

$$\mathbf{u}(\mathbf{x},t)={\mathit{\xi }}^{\prime }(t)+{\theta }^{\prime }(t){(\mathbf{x}-\mathit{\xi }(t))}^{\perp },\mathbf{x}\in \partial \mathcal{S}(\mathit{\xi }(t),\theta (t),t),t\in [0,T],$$

because the undulatory velocity field w is equal to zero.

2.3 Weak formulation

Let $\mathit{\xi }\in H^2((0,T);{\mathbb{R}}^2)$, $\theta \in H^2((0,T);\mathbb{R})$ be two functions such that $\overline{\mathcal{S}(\mathit{\xi }(t),\theta (t),t)}\subset \mathrm{\Omega }$ for all $t\in [0,T]$. In the sequel, we define $\mathcal{F}=\mathcal{F}(\mathit{\xi }(0),\theta (0),0)$ and $\mathcal{S}=\mathcal{S}(\mathit{\xi }(0),\theta (0),0)$. Moreover, if no confusion is possible, we define

$$\mathcal{S}(t)=\mathcal{S}(\mathit{\xi }(t),\theta (t),t),\mathcal{F}(t)=\mathcal{F}(\mathit{\xi }(t),\theta (t),t)=\mathrm{\Omega }\setminus \overline{\mathcal{S}(t)}.$$

Let $\mathrm{\Psi }:{\mathbb{R}}^2\times [0,T]\to {\mathbb{R}}^2$ be a mapping such that for every $t\in [0,T]$, the function $\mathrm{\Psi }(\cdot ,t)$ is a $C^{\mathrm{\infty }}$-diffeomorphism from ℱ onto $\mathcal{F}(\mathit{\xi }(t),\theta (t),t)$ and such that the derivatives

$$\frac{{\partial }^{i+k_1+k_2}\mathrm{\Psi }}{\partial t^i\partial y_1^{k_1}\partial y_2^{k_2}},i\le 1,k_1\ge 0,k_2\ge 0$$

exist and are continuous. The existence of such a function is due, in particular, to the fact that $dist(\mathcal{S}(t),\partial \mathrm{\Omega })>0$ for all t (see [19]). We can now define the following functions spaces:

$$\begin{array}{c}{L}^2(0,T;H^2{(\mathcal{F}(t))}^2)=\{\mathbf{u}\mid {\mathbf{u}}_{\mathrm{\Psi }}\in L^2(0,T;H^2{(\mathcal{F})}^2)\},\hfill \\ H^1(0,T;L^2{(\mathcal{F}(t))}^2)=\{\mathbf{u}\mid {\mathbf{u}}_{\mathrm{\Psi }}\in H^1(0,T;L^2{(\mathcal{F})}^2)\},\hfill \\ C([0,T];H^1{(\mathcal{F}(t))}^2)=\{\mathbf{u}\mid {\mathbf{u}}_{\mathrm{\Psi }}\in C([0,T];H^1{(\mathcal{F})}^2)\},\hfill \\ L^2(0,T;H^1(\mathcal{F}(t)))=\{p\mid p_{\mathrm{\Psi }}\in L^2(0,T;H^1(\mathcal{F}))\},\hfill \end{array}$$

where $v_{\mathrm{\Psi }}$ denotes the function defined by $v_{\mathrm{\Psi }}(\mathbf{y},t)=v(\mathrm{\Psi }(\mathbf{y},t),t)$ for $(\mathbf{y},t)\in \mathcal{F}\times (0,T)$.

In order to introduce the weak formulation, we first define some additional functions spaces. For every $t\ge 0$, let $(\mathit{\xi },\theta )$ be an arbitrary position of the creature at time t, such that $\overline{\mathcal{S}(\mathit{\xi },\theta ,t)}\subset \mathrm{\Omega }$. We denote

$$\mathcal{K}(\mathit{\xi },\theta ,t)=\{\mathbf{u}\in H_0^1{(\mathrm{\Omega })}^2:\mathbf{D}(\mathbf{u})=0\text{ in }\mathcal{S}(\mathit{\xi },\theta ,t)\},$$

(23)

$$\mathcal{M}(\mathit{\xi },\theta ,t)=\{p\in L^2(\mathrm{\Omega }):{\int }_{\mathrm{\Omega }}pd\mathbf{x}=0\text{ and }p=0\text{ in }\mathcal{S}(\mathit{\xi },\theta ,t)\},$$

(24)

where $\mathbf{D}(\mathbf{u})$ is the strain rate tensor defined by (19).

Let $(\mathbf{u},p,\mathit{\xi },\theta )$ be a solution of (13)-(22). The vector velocity field u and the pressure p can be extended to Ω by setting

$$\mathbf{u}(\mathbf{x},t)={\mathit{\xi }}^{\prime }(t)+{\theta }^{\prime }(t){(\mathbf{x}-\mathit{\xi }(t))}^{\perp }+\mathbf{w}(\mathbf{x},t)\text{if }\mathbf{x}\in \mathcal{S}(\mathit{\xi }(t),\theta (t),t),$$

(25)

$$p(\mathbf{x},t)=0\text{if }\mathbf{x}\in \mathcal{S}(\mathit{\xi }(t),\theta (t),t).$$

(26)

The extended vector $\mathbf{u}(\cdot ,t)$ belongs to $H_0^1{(\mathrm{\Omega })}^2$. In the remaining part of this paper, the solution u and p of (13)-(22) will be extended as above.

We also need to extend the density field ${\rho }_S$ of the creature (defined in (10)) to the whole domain Ω by setting

$$\rho (\mathbf{x},t)=\{\begin{array}{cc}{\rho }_F\hfill & \text{for all }\mathbf{x}\in \mathcal{F}(\mathit{\xi }(t),\theta (t),t),\hfill \\ {\rho }_S(\mathbf{x},t)\hfill & \text{for all }\mathbf{x}\in \mathcal{S}(\mathit{\xi }(t),\theta (t),t),\hfill \end{array}t\ge 0.$$

(27)

By a slight variation of the argument in Ladyzhenskaya [[20], p.27], it can be shown that for every $\delta >0$, there exists a continuous function $({\mathbf{x}}^{\ast },t)\mapsto {\mathrm{\Lambda }}^{\ast }({\mathbf{x}}^{\ast },t)$ such that, for every $t\ge 0$, the map ${\mathbf{x}}^{\ast }\mapsto {\mathrm{\Lambda }}^{\ast }({\mathbf{x}}^{\ast },t)$ is $C^{\mathrm{\infty }}$ on ${\mathbb{R}}^2\setminus {\mathcal{S}}^{\ast }(t)$ and such that the function $t\mapsto {\mathrm{\Lambda }}^{\ast }({\mathbf{x}}^{\ast },t)$ is of class $C^{\mathrm{\infty }}$ for every ${\mathbf{x}}^{\ast }\in {\mathbb{R}}^2\setminus \overline{{\mathcal{S}}^{\ast }(t)}$ and

$$\{\begin{array}{cc}div{\mathrm{\Lambda }}^{\ast }=0\hfill & \text{in }{\mathbb{R}}^2\setminus \overline{{\mathcal{S}}^{\ast }(t)},t\in (0,T),\hfill \\ {\mathrm{\Lambda }}^{\ast }({\mathbf{x}}^{\ast },t)=0\hfill & \text{if }dist({\mathbf{x}}^{\ast },{\mathcal{S}}^{\ast }(t))\ge \delta >0,t\in (0,T),\hfill \\ {\mathrm{\Lambda }}^{\ast }({\mathbf{x}}^{\ast },t)={\mathbf{w}}^{\ast }({\mathbf{x}}^{\ast },t)\hfill & \text{if }{\mathbf{x}}^{\ast }\in {\mathcal{S}}^{\ast }(t),t\in (0,T).\hfill \end{array}$$

(28)

For every $t\ge 0$, let $(\mathit{\xi },\theta )$ be an arbitrary position of the creature at time t, such that $\overline{\mathcal{S}(\mathit{\xi },\theta ,t)}\subset \mathrm{\Omega }$ . We then define $\mathrm{\Lambda }(\cdot ,t;\mathit{\xi },\theta )$ by

$$\mathrm{\Lambda }(\mathbf{x},t;\mathit{\xi },\theta )={\mathcal{R}}_{\theta }{\mathrm{\Lambda }}^{\ast }({\mathcal{R}}_{-\theta }(\mathbf{x}-\mathit{\xi }),t)\mathrm{\forall }\mathbf{x}\in {\mathbb{R}}^2.$$

(29)

Then the function Λ satisfies

$$\{\begin{array}{cc}div\mathrm{\Lambda }=0\hfill & \text{in }{\mathbb{R}}^2\setminus \overline{\mathcal{S}(\mathit{\xi },\theta ,t)},\hfill \\ \mathrm{\Lambda }=0\hfill & \text{on }\partial \mathrm{\Omega },\hfill \\ \mathrm{\Lambda }=\mathbf{w}\hfill & \text{in }\mathcal{S}(\mathit{\xi },\theta ,t).\hfill \end{array}$$

(30)

An important ingredient of the numerical method we use is given by the characteristic functions whose level lines are the integral curves of the velocity field. More precisely (see, for instance, [10, 11]) the characteristic function $\mathit{\psi }:{[0,T]}^2\times \mathrm{\Omega }\to \mathrm{\Omega }$ is defined as the solution of the initial value problem

$$\{\begin{array}{c}\frac{d}{dt}\mathit{\psi }(t;s,\mathbf{x})=\mathbf{u}(\mathit{\psi }(t;s,\mathbf{x}),t)\mathrm{\forall }t\in [0,T],\hfill \\ \mathit{\psi }(s;s,\mathbf{x})=\mathbf{x}.\hfill \end{array}$$

(31)

It is well known that the material derivative $D_t\mathbf{u}=\partial \mathbf{u}/\partial t+(\mathbf{u}\cdot \mathrm{\nabla })\mathbf{u}$ of the velocity field u at instant $t_0$ satisfies:

$$D_t\mathbf{u}(\mathbf{x},t_0)=\frac{d}{dt}{\left[\mathbf{u}(\mathit{\psi }(t;t_0,\mathbf{x}),t)\right]}_{{|}_{t=t_0}}.$$

(32)

Remark 2.3 By using a classical result of Liouville (see, for instance, [[21], p.251]), if

$$\mathit{\xi }\in H^2{(0,T)}^2,\theta \in H^2(0,T),\mathbf{u}\in C([0,T];H_0^1{(\mathrm{\Omega })}^2)$$

are such that for any $t\in [0,T]$, we have $\mathcal{S}(\mathit{\xi }(t),\theta (t),t)\subset \mathrm{\Omega }$ and

$$div\mathbf{u}=0\text{in }\mathcal{F}(\mathit{\xi }(t),\theta (t),t),$$

then we get

$$det{\mathbf{J}}_{\mathit{\psi }}(t,s,\mathbf{x})=1\mathrm{\forall }\mathbf{x}\in \mathcal{F}(\mathit{\xi }(t),\theta (t),t),$$

(33)

where we have denoted by

$${\mathbf{J}}_{\mathit{\psi }}={\left(\frac{\partial {\mathit{\psi }}_i}{\partial y_j}\right)}_{i,j}$$

the Jacobian matrix of the transformation $\mathbf{y}\mapsto \mathit{\psi }(\mathbf{y})$.

In order to give the global weak formulation of our problem, we need to introduce the bilinear forms

$$a:H_0^1{(\mathrm{\Omega })}^2\times H_0^1{(\mathrm{\Omega })}^2\to \mathbb{R},b:H_0^1{(\mathrm{\Omega })}^2\times L^2(\mathrm{\Omega })\to \mathbb{R},$$

defined by

$$\begin{array}{c}a(\mathbf{u},\mathbf{v})=2\mu {\int }_{\mathrm{\Omega }}\mathbf{D}(\mathbf{u}):\mathbf{D}(\mathbf{v})d\mathbf{x}\mathrm{\forall }\mathbf{u},\mathbf{v}\in H_0^1{(\mathrm{\Omega })}^2,\hfill \\ b(\mathbf{u},q)=-{\int }_{\mathrm{\Omega }}(div\mathbf{u})qd\mathbf{x}\mathrm{\forall }\mathbf{u}\in H_0^1{(\mathrm{\Omega })}^2,q\in L^2(\mathrm{\Omega }).\hfill \end{array}$$

Proposition 2.4 Assume that

$$\begin{array}{c}\mathbf{u}\in L^2(0,T;H^2{(\mathcal{F}(t))}^2)\cap H^1(0,T;L^2{(\mathcal{F}(t))}^2)\cap C([0,T];H^1{(\mathcal{F}(t))}^2),\hfill \\ p\in L^2(0,T;H^1(\mathcal{F}(t))),\mathit{\xi }\in H^2{(0,T)}^2,\theta \in H^2(0,T),\hfill \end{array}$$

and that u and p are extended to Ω as above.

Then $(\mathbf{u},p,\mathit{\xi },\theta )$ is the solution of (13)-(22) if and only if for all $t\in [0,T]$, $\mathbf{u}(\cdot ,t)-\mathrm{\Lambda }(\cdot ,t;\mathit{\xi }(t),\theta (t))\in \mathcal{K}(\mathit{\xi }(t),\theta (t),t)$, $p(\cdot ,t)\in \mathcal{M}(\mathit{\xi }(t),\theta (t),t)$, and $(\mathbf{u},p)$ satisfies

$$(\rho D_t\mathbf{u},\mathbf{v})+a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)=(\rho \mathbf{f}(t),\mathbf{v})\mathrm{\forall }\mathbf{v}\in \mathcal{K}(\mathit{\xi }(t),\theta (t),t),$$

(34)

$$b(\mathbf{u},q)=0\mathrm{\forall }q\in \mathcal{M}(\mathit{\xi }(t),\theta (t),t),$$

(35)

for a.e. $t\in (0,T)$.

More details on the existence and uniqueness of the solution and the complete proof of this result could be found in [19].

In the remainder of the paper, we suppose that f and ${\mathbf{u}}_0$ satisfy

$$\begin{array}{r}\mathbf{f}\in C([0,T];H^1{(\mathrm{\Omega })}^2),{\mathbf{u}}_0\in H^2{(\mathcal{F})}^2,div({\mathbf{u}}_0)=0\text{in }\mathcal{F},\\ {\mathbf{u}}_0=0\text{on }\partial \mathrm{\Omega },\\ {\mathbf{u}}_0(\mathbf{x})={\mathit{\xi }}_1+{\omega }_0{(\mathbf{x}-{\mathit{\xi }}_0)}^{\perp }+{\mathcal{R}}_{{\theta }_0}{\mathbf{w}}^{\ast }({\mathcal{R}}_{-{\theta }_0}(\mathbf{x}-{\xi }_0),0)\text{for }\mathbf{x}\in \partial \mathcal{S},\end{array}$$

(36)

where ${\mathit{\xi }}_0,{\mathit{\xi }}_1\in {\mathbb{R}}^2$ and ${\theta }_0$, ${\omega }_0\in \mathbb{R}$ are given as initial data in (22). Let us also assume that the corresponding solution $(\mathbf{u},p,\mathit{\xi },\theta )$ of problem (13)-(22) satisfies the following regularity properties:

$$\{\begin{array}{c}\mathbf{u}\in C([0,T];H^2{(\mathcal{F}(t))}^2)\cap H^1(0,T;L^2{(\mathcal{F}(t))}^2),\hfill \\ D_t^2\mathbf{u}\in L^2(0,T;L^2{(\mathcal{F}(t))}^2),\mathbf{u}\in C([0,T];C^{0,1}{(\overline{\mathrm{\Omega }})}^2),\hfill \\ p\in C([0,T];H^1(\mathcal{F}(t))),\mathit{\xi }\in H^3{(0,T)}^2,\theta \in H^3(0,T).\hfill \end{array}$$

(37)

Moreover, we suppose that there exists a nonempty open connected subset ${\mathrm{\Omega }}_0$ of Ω such that for any ${\mathit{\xi }}_1,{\mathit{\xi }}_2\in {\mathrm{\Omega }}_0$, we have

$$\mathcal{S}({\mathit{\xi }}_1+\lambda ({\mathit{\xi }}_2-{\mathit{\xi }}_1),\theta ,t)\subset \mathrm{\Omega }\mathrm{\forall }\lambda \in [0,1],\mathrm{\forall }\theta \in [0,2\pi ],\mathrm{\forall }t\in [0,T].$$

(38)

Using this notation, we assume that

$$\xi (t)\in {\mathrm{\Omega }}_0,dist(\mathcal{S}(t),\partial \mathrm{\Omega })>0\mathrm{\forall }t\in [0,T].$$

(39)

Remark 2.5 The hypotheses (37) and (39) imply the existence of $\eta >0$ such that

$$dist(\xi (t),\partial {\mathrm{\Omega }}_0)>3\eta ,dist(\mathcal{S}(t),\partial \mathrm{\Omega })>3\eta \mathrm{\forall }t\in [0,T].$$

(40)

3 Time discretization and first main result

In this section, based on a weak form of the governing equations, we describe a method for the time discretization of (13)-(22).

Let us first divide the time interval $[0,T]$ into subintervals $[t^k,t^{k+1}]$ with $t^{k+1}-t^k=\mathrm{\Delta }t=\frac{T}{N}$, where N is a positive integer and $k\in \{1,\dots ,N\}$. Let $({\mathbf{u}}^k,p^k,{\mathit{\xi }}^k,{\theta }^k)$ be the approximation of the solution of (34)-(35) at time $t=t^k$ (remark that ${\mathbf{u}}^k$ and $p^k$ are functions defined on the whole domain Ω). We denote

$$\begin{array}{r}{\mathcal{K}}^k=\mathcal{K}({\mathit{\xi }}^k,{\theta }^k,t^k),{\mathcal{M}}^k=\mathcal{M}({\mathit{\xi }}^k,{\theta }^k,t^k),\\ {\mathcal{S}}^k=\mathcal{S}({\mathit{\xi }}^k,{\theta }^k,t^k),{\mathcal{F}}^k=\mathcal{F}({\mathit{\xi }}^k,{\theta }^k,t^k),\end{array}$$

(41)

and we consider the functions

$$\begin{array}{r}{\mathrm{\Lambda }}^k(\mathbf{x})={\mathcal{R}}_{{\theta }^k}{\mathrm{\Lambda }}^{\ast }({\mathcal{R}}_{-{\theta }^k}(\mathbf{x}-{\mathit{\xi }}^k),t^k)\mathrm{\forall }\mathbf{x}\in {\mathbb{R}}^2,\\ {\rho }^k(\mathbf{x})=\{\begin{array}{cc}{\rho }_S^{\ast }({\mathcal{R}}_{-{\theta }^k}(\mathbf{x}-{\mathit{\xi }}^k),t^k)\hfill & \text{if }\mathbf{x}\in {\mathcal{S}}^k,\hfill \\ {\rho }_F\hfill & \text{if }\mathbf{x}\in {\mathcal{F}}^k.\hfill \end{array}\end{array}$$

(42)

Remark 3.1 Combining the regularity properties of ${\mathrm{\Lambda }}^{\ast }$ and ${\mathbf{w}}^{\ast }$, it follows that

$${\mathrm{\Lambda }}^k\in {\mathcal{C}}^{0,1}{(\mathrm{\Omega })}^2.$$

(43)

Moreover, taking $\delta <\eta$ in definition (28) of ${\mathrm{\Lambda }}^{\ast }$, where η is defined in (40), we have that ${\mathrm{\Lambda }}^k+{\mathcal{K}}^k\subseteq H_0^1{(\mathrm{\Omega })}^2$.

Now, let us describe the numerical scheme for approximating the solutions of (13)-(22). This procedure is based on the weak form derived in Proposition 2.4.

The first step of our scheme consists in computing the new position of the mass center and the new orientation of the creature by setting

$${\mathit{\xi }}^{k+1}={\mathit{\xi }}^k+\mathrm{\Delta }t{\mathbf{u}}^k\left({\mathit{\xi }}^k\right),$$

(44)

$${\theta }^{k+1}={\theta }^k+\frac{\mathrm{\Delta }t}{I(t^k)}{\int }_{{\mathcal{S}}^k}{\rho }^k({\mathbf{u}}^k(\mathbf{x})-{\mathbf{u}}^k\left({\mathit{\xi }}^k\right))\cdot {(\mathbf{x}-{\mathit{\xi }}^k)}^{\perp }d\mathbf{x}.$$

(45)

The second step consists in computing the global velocity field ${\mathbf{u}}^{k+1}$ and the global pressure field $p^{k+1}$. To this end, we look for ${\mathbf{u}}^{k+1}\in {\mathrm{\Lambda }}^{k+1}+{\mathcal{K}}^{k+1}$ and $p^{k+1}\in {\mathcal{M}}^{k+1}$ such that for all $\mathbf{v}\in {\mathcal{K}}^{k+1}$, and for all $q\in {\mathcal{M}}^{k+1}$, we have

$$({\rho }^{k+1}\frac{{\mathbf{u}}^{k+1}-{\mathbf{u}}^k\circ {\overline{\mathbf{X}}}_{{\mathbf{u}}^k}^k}{\mathrm{\Delta }t},\mathbf{v})+a({\mathbf{u}}^{k+1},\mathbf{v})+b(\mathbf{v},p^{k+1})=({\rho }^{k+1}{\mathbf{f}}^{k+1},\mathbf{v}),$$

(46)

$$b({\mathbf{u}}^{k+1},q)=0,$$

(47)

where ${\mathbf{f}}^{k+1}(\mathbf{x})=\mathbf{f}(\mathbf{x},t^{k+1})$ for any $\mathbf{x}\in \mathrm{\Omega }$.

In the equations above, the approximate characteristic is given by

$${\overline{\mathbf{X}}}_{{\mathbf{u}}^k}^k(\mathbf{x})={\mathit{\chi }}^k(t^k;t^{k+1},\mathbf{x}),$$

(48)

for all $\mathbf{x}\in \mathrm{\Omega }$, where ${\mathit{\chi }}^k$ is the solution of the problem

$$\{\begin{array}{c}\frac{d}{dt}{\mathit{\chi }}^k(t;t^{k+1},\mathbf{x})=\widetilde{{\mathbf{u}}^k}({\mathit{\chi }}^k(t;t^{k+1},\mathbf{x})),t\in [t^k,t^{k+1}],\hfill \\ {\mathit{\chi }}^k(t^{k+1};t^{k+1},\mathbf{x})={\mathit{\xi }}^k+{\mathcal{R}}_{{\theta }^k}\mathrm{\Pi }(t^k;t^{k+1},{\mathcal{R}}_{-{\theta }^{k+1}}(\mathbf{x}-{\mathit{\xi }}^{k+1})),\hfill \end{array}$$

(49)

with

$$\widetilde{{\mathbf{u}}^k}(\mathbf{z})={\mathbf{u}}^k(\mathbf{z})-{\mathbf{u}}^k\left({\mathit{\xi }}^k\right)-\frac{{\theta }^{k+1}-{\theta }^k}{\mathrm{\Delta }t}{(\mathbf{z}-{\mathit{\xi }}^k)}^{\perp }-{\mathrm{\Lambda }}^k(\mathbf{z})\mathrm{\forall }\mathbf{z}\in {\mathbb{R}}^2,$$

where ${\mathbf{u}}^k(\mathbf{z})$ is extended by zero outside of Ω.

For all $s\in [0,T]$, and $\mathbf{z}\in {\mathbb{R}}^2$, function $\mathrm{\Pi }(\cdot ;s,\mathbf{z})$ corresponds to the characteristic function of the extended undulatory velocity ${\mathrm{\Lambda }}^{\ast }$, defined by

$$\{\begin{array}{c}\frac{d}{dt}\mathrm{\Pi }(t;s,\mathbf{z})={\mathrm{\Lambda }}^{\ast }(\mathrm{\Pi }(t;s,\mathbf{z}),t),t\in [0,T],\hfill \\ \mathrm{\Pi }(s;s,\mathbf{z})=\mathbf{z}.\hfill \end{array}$$

(50)

Remark 3.2 Let us note that for any $\mathbf{z}\in {\mathcal{S}}^{\ast }(t^{k+1})$, equation (50) has the explicit solution

$$\mathrm{\Pi }(t;t^{k+1},\mathbf{z})={\mathbf{X}}^{\ast }({\mathbf{Y}}^{\ast }(\mathbf{z},t^{k+1}),t)\in {\mathcal{S}}^{\ast }(t).$$

(51)

Then for any $\mathbf{x}\in \mathcal{S}({\mathit{\xi }}^{k+1},{\theta }^{k+1},t^{k+1})$, since

$${\mathcal{R}}_{-{\theta }^{k+1}}(\mathbf{x}-{\mathit{\xi }}^{k+1})\in {\mathcal{S}}^{\ast }\left(t^{k+1}\right),$$

we obtain that the initial condition in (49) is

$${\mathit{\chi }}^k(t^{k+1};t^{k+1},\mathbf{x})={\mathit{\xi }}^k+{\mathcal{R}}_{{\theta }^k}{\mathbf{X}}^{\ast }({\mathbf{Y}}^{\ast }({\mathcal{R}}_{-{\theta }^{k+1}}(\mathbf{x}-{\mathit{\xi }}^{k+1}),t^{k+1}),t^k).$$

Moreover, since $\mathbf{z}={\mathit{\chi }}^k(t^{k+1};t^{k+1},\mathbf{x})\in \mathcal{S}({\mathit{\xi }}^k,{\theta }^k,t^k)$, we have $\widetilde{{\mathbf{u}}^k}(\mathbf{z})=0$ and

$${\mathit{\chi }}^k(t;t^{k+1},\mathbf{x})={\mathit{\xi }}^k+{\mathcal{R}}_{{\theta }^k}{\mathbf{X}}^{\ast }({\mathbf{Y}}^{\ast }({\mathcal{R}}_{-{\theta }^{k+1}}(\mathbf{x}-{\mathit{\xi }}^{k+1}),t^{k+1}),t^k)\mathrm{\forall }t\in [t^k,t^{k+1}].$$

In particular, for $t=t^k$ we have that

$${\overline{\mathbf{X}}}_{{\mathbf{u}}^k}^k(\mathbf{x})={\mathit{\xi }}^k+{\mathcal{R}}_{{\theta }^k}{\mathbf{X}}^{\ast }({\mathbf{Y}}^{\ast }({\mathcal{R}}_{-{\theta }^{k+1}}(\mathbf{x}-{\mathit{\xi }}^{k+1}),t^{k+1}),t^k)$$

(52)

for all $\mathbf{x}\in \mathcal{S}({\mathit{\xi }}^{k+1},{\theta }^{k+1},t^{k+1})$.

It is easy to see that for any $k\in \{1,\dots ,N\}$, equations (46)-(47) represent a mixed formulation of a well-posed Navier-Stokes-type system, so that our scheme is well defined.

Let us now state our first main result concerning the convergence of the semi-discrete scheme (46)-(47).

Theorem 3.3 Suppose that Ω is an open smooth bounded domain in ${\mathbb{R}}^2$, f and ${\mathbf{u}}_0$ satisfy (36), and the exact solution $(\mathbf{u},p,\mathit{\xi },\theta )$ of problem (13)-(22) satisfies hypotheses (37)-(39). Then there exist a constant $K>0$ depending on T and a constant $\delta >0$ independent of T such that for all $0<\mathrm{\Delta }t\le \delta$, the solution ${({\mathbf{u}}^k,p^k,{\mathit{\xi }}^k,{\theta }^k)}_{k\in \{1,\dots ,N\}}$ of the time-discretization problem (46)-(47) satisfies

$$\underset{1\le k\le N}{sup}({\parallel \mathbf{u}\left(t^k\right)-{\mathbf{u}}^k\parallel }_{L^2{(\mathrm{\Omega })}^2}+|\mathit{\xi }\left(t^k\right)-{\mathit{\xi }}^k|+|\theta \left(t^k\right)-{\theta }^k|)\le K\mathrm{\Delta }t.$$

(53)

The complete proof of this result could be found in [18] for the case of rigid solid and in the forthcoming paper [22] for the deformable structure.

4 Fully discrete formulation and second main result

In order to discretize problem (46)-(47) with respect to the space variable, we introduce two families of finite element spaces which approximate spaces ${\mathcal{K}}^k$ and ${\mathcal{M}}^k$ defined in (41), (23) and (24). To this end, for any discretization parameter $h\in (0,1)$, we consider a quasi-uniform triangulation ${\mathcal{T}}_h$ of the domain Ω. Suppose that Ω is a bounded convex domain with a polygonal boundary. We denote by ${\mathcal{W}}_h$ the ${\mathbb{P}}_1+\mathrm{bubble}$ finite elements space associated with ${\mathcal{T}}_h$ for the velocity field and by $E_h$ the ${\mathbb{P}}_1$-finite elements space for the pressure, that is,

$$\begin{array}{r}{\mathcal{W}}_h=\{\mathbf{v}\in \mathcal{C}{(\overline{\mathrm{\Omega }})}^2:\mathrm{\forall }T\in {\mathcal{T}}_h,{\mathbf{v}}_{{|}_T}\in {\mathbb{P}}_1+\mathrm{bubble}(T)\},\\ E_h=\{q\in \mathcal{C}(\overline{\mathrm{\Omega }}):\mathrm{\forall }T\in {\mathcal{T}}_h,q_{{|}_T}\in {\mathbb{P}}_1\}.\end{array}$$

Then, we define the following finite elements spaces for a conform approximation of the fluid-solid system:

$$\begin{array}{c}{\mathcal{K}}_h(\mathit{\xi },\theta ,t)={\mathcal{W}}_h\cap \mathcal{K}(\mathit{\xi },\theta ,t)\mathrm{\forall }\mathit{\xi }\in {\mathrm{\Omega }}_0,\mathrm{\forall }\theta \in [0,2\pi ],\mathrm{\forall }t\in [0,T],\hfill \\ {\mathcal{M}}_h(\mathit{\xi },\theta ,t)=E_h\cap \mathcal{M}(\mathit{\xi },\theta ,t)\mathrm{\forall }\mathit{\xi }\in {\mathrm{\Omega }}_0,\mathrm{\forall }\theta \in [0,2\pi ],\mathrm{\forall }t\in [0,T].\hfill \end{array}$$

Let us recall an approximation property of the projection on $K_h(\mathit{\xi },\theta ,t)\times {\mathcal{M}}_h(\mathit{\xi },\theta ,t)$ (see [2]).

Lemma 4.1 Suppose that $\mathbf{V}\in \mathcal{K}(\mathit{\xi },\theta ,t)$ and $P\in \mathcal{M}(\mathit{\xi },\theta ,t)$. Then there exists a unique couple $({\overline{\mathbf{V}}}_h,{\overline{P}}_h)$ in ${\mathcal{K}}_h(\mathit{\xi },\theta ,t)\times {\mathcal{M}}_h(\mathit{\xi },\theta ,t)$ such that

$$\{\begin{array}{c}a(\mathbf{V}-{\overline{\mathbf{V}}}_h,\mathbf{v})+b(\mathbf{v},P-{\overline{P}}_h)=0\mathrm{\forall }\mathbf{v}\in {\mathcal{K}}_h(\mathit{\xi },\theta ,t),\hfill \\ b(\mathbf{V}-{\overline{\mathbf{V}}}_h,q)=0\mathrm{\forall }q\in {\mathcal{M}}_h(\mathit{\xi },\theta ,t).\hfill \end{array}$$

(54)

In addition, if we suppose that ${\mathbf{V}}_{|\mathrm{\Omega }\setminus \mathcal{S}(\mathit{\xi },\theta ,t)}\in H^2{(\mathrm{\Omega }\setminus \mathcal{S}(\mathit{\xi },\theta ,t))}^2$ and $P_{|\mathrm{\Omega }\setminus \mathcal{S}(\mathit{\xi },\theta ,t)}\in H^1(\mathrm{\Omega }\setminus \mathcal{S}(\mathit{\xi },\theta ,t))$, then there exists a positive constant C, independent of h, such that

$${\parallel \mathbf{V}-{\overline{\mathbf{V}}}_h\parallel }_{L^2{(\mathrm{\Omega })}^2}\le Ch.$$

In order to define the approximate characteristics, let us denote by $F_h$ the ${\mathbb{P}}_2$-finite elements space associated with the triangulation ${\mathcal{T}}_h$, and we introduce the space

$${\mathcal{R}}_h(\mathit{\xi },\theta ,t)=\{{\mathrm{\nabla }}^{\perp }{\varphi }_h:{\varphi }_h\in F_h,{\varphi }_h=0\text{ on }\partial \mathrm{\Omega }\}\cap \mathcal{K}(\mathit{\xi },\theta ,t),$$

$\mathrm{\forall }\mathit{\xi }\in {\mathrm{\Omega }}_0$, $\mathrm{\forall }\theta \in [0,2\pi ]$, $\mathrm{\forall }t\in [0,T]$, where ${\mathrm{\nabla }}^{\perp }{\varphi }_h=\left(\begin{array}{c}-\frac{\partial {\varphi }_h}{\partial y}\\ \frac{\partial {\varphi }_h}{\partial x}\end{array}\right)$.

We denote by $\mathbf{P}(\mathit{\xi },\theta ,t)$ the orthogonal projection from $L^2{(\mathrm{\Omega })}^2$ onto ${\mathcal{R}}_h(\mathit{\xi },\theta ,t)$, i.e., for any $\mathbf{u}\in L^2{(\mathrm{\Omega })}^2$, then the projection $\mathbf{P}(\mathit{\xi },\theta ,t)\mathbf{u}\in {\mathcal{R}}_h(\mathit{\xi },\theta ,t)$ is such that $(\mathbf{u}-\mathbf{P}(\mathit{\xi },\theta ,t)\mathbf{u},{\mathbf{r}}_h)=0$ for all ${\mathbf{r}}_h\in {\mathcal{R}}_h(\mathit{\xi },\theta ,t)$.

Let N be a positive integer. We denote $\mathrm{\Delta }t=T/N$ and $t^k=k\mathrm{\Delta }t$ for all $k\in \{0,\dots ,N\}$. For $k=0$, we define

$${\mathbf{u}}_h^0(\cdot )=\overline{\mathbf{u}(\cdot ,0)},{\mathit{\xi }}_h^0={\mathit{\xi }}_0\text{and}{\theta }_h^0={\theta }_0,$$

(55)

where $(\overline{\mathbf{u}(\cdot ,0)},\overline{p(\cdot ,0)})\in {\mathcal{K}}_h({\mathit{\xi }}_0,{\theta }_0,0)\times {\mathcal{M}}_h({\mathit{\xi }}_0,{\theta }_0,0)$ is the projection of the initial condition $(\mathbf{u}(\cdot ,0),p(\cdot ,0))$ on ${\mathcal{K}}_h({\mathit{\xi }}_0,{\theta }_0,0)\times {\mathcal{M}}_h({\mathit{\xi }}_0,{\theta }_0,0)$ defined in (54).

Assume that the approximate solution $({\mathbf{u}}_h^k,p_h^k,{\mathit{\xi }}_h^k,{\theta }_h^k)$ of (13)-(18) at time $t=t^k$ is known. We describe below the numerical scheme allowing to determinate the approximate solution $({\mathbf{u}}_h^{k+1},p_h^{k+1},{\mathit{\xi }}_h^{k+1},{\theta }_h^{k+1})$ at $t=t^{k+1}$. First, we compute ${\mathit{\xi }}_h^{k+1}\in {\mathbb{R}}^2$ and ${\theta }_h^{k+1}\in \mathbb{R}$ by

$${\mathit{\xi }}_h^{k+1}={\mathit{\xi }}_h^k+{\mathbf{u}}_h^k\left({\mathit{\xi }}_h^k\right)\mathrm{\Delta }t,$$

(56)

$${\theta }_h^{k+1}={\theta }_h^k+\frac{\mathrm{\Delta }t}{I(t^k)}{\int }_{{\mathcal{S}}^k}{\rho }_h^k({\mathbf{u}}_h^k(\mathbf{x})-{\mathbf{u}}_h^k\left({\mathit{\xi }}_h^k\right))\cdot {(\mathbf{x}-{\mathit{\xi }}_h^k)}^{\perp }d\mathbf{x},$$

(57)

where ${\rho }_h^k$ is defined by the identity (60) below.

We consider the approximated characteristic function ${\mathit{\chi }}_h^k$ defined as the solution of the following ordinary differential equation:

$$\{\begin{array}{c}\frac{d}{dt}{\mathit{\chi }}_h^k(t;t^{k+1},\mathbf{x})=\mathbf{P}({\mathit{\xi }}_h^k,{\theta }_h^k,t^k)\widetilde{{\mathbf{u}}_h^k}({\mathit{\chi }}_h^k(t;t^{k+1},\mathbf{x})),t\in [t^k,t^{k+1}],\hfill \\ {\mathit{\chi }}_h^k(t^{k+1};t^{k+1},\mathbf{x})={\mathit{\xi }}_h^k+{\mathcal{R}}_{{\theta }_h^k}\mathrm{\Pi }(t^k;t^{k+1},{\mathcal{R}}_{-{\theta }_h^{k+1}}(\mathbf{x}-{\mathit{\xi }}_h^{k+1})),\hfill \end{array}$$

(58)

with

$$\widetilde{{\mathbf{u}}_h^k}(\mathbf{z})={\mathbf{u}}_h^k(\mathbf{z})-{\mathbf{u}}_h^k\left({\mathit{\xi }}_h^k\right)-\frac{{\theta }_h^{k+1}-{\theta }_h^k}{\mathrm{\Delta }t}{(\mathbf{z}-{\mathit{\xi }}_h^k)}^{\perp }-{\mathrm{\Lambda }}_h^k(\mathbf{z})\mathrm{\forall }\mathbf{z}\in {\mathbb{R}}^2,$$

where ${\mathbf{u}}_h^k(\mathbf{z})$ is extended by zero outside of Ω and

$${\mathrm{\Lambda }}_h^k(\mathbf{x})={\mathcal{R}}_{{\theta }_h^k}{\mathrm{\Lambda }}^{\ast }({\mathcal{R}}_{-{\theta }_h^k}(\mathbf{x}-{\mathit{\xi }}_h^k),t^k)\mathrm{\forall }\mathbf{x}\in {\mathbb{R}}^2.$$

The characteristic function Π is defined by (50).

Finally, we define

$${\overline{\mathbf{X}}}_h^k(\mathbf{x})={\mathit{\chi }}_h^k(t^k;t^{k+1},\mathbf{x})\mathrm{\forall }\mathbf{x}\in \mathrm{\Omega }.$$

(59)

In the sequel, we shall split the mesh into the union of 4 different types of triangle subsets. We first introduce $A_h$ as the union of all triangles intersecting the solid $\mathcal{S}({\mathit{\xi }}_h^k,{\theta }_h^k,t^k)$, i.e.,

$$A_h=\underset{\stackrel{\circ }{T}\cap \stackrel{\circ }{\mathcal{S}}({\mathit{\xi }}_h^k,{\theta }_h^k,t^k)\ne \mathrm{\varnothing }}{\bigcup _{T\in {\mathcal{T}}_h}}T.$$

We also denote by $Q_h$ the union of all triangles such that all their vertices are contained in $\overline{A_h}$. The triangles of ${\mathcal{T}}_h$ are then split into the following four categories (see Figure 1):

• ${\mathcal{F}}_1$ is the subset of ${\mathcal{T}}_h$ formed by all triangles $T\in {\mathcal{T}}_h$ such that $\overline{T}\subset \mathcal{S}({\mathit{\xi }}_h^k,{\theta }_h^k,t^k)$.

• ${\mathcal{F}}_2$ is the subset formed by all triangles $T\in {\mathcal{T}}_h\setminus {\mathcal{F}}_1$ such that $\overline{T}\subset \overline{Q_h}$.

• ${\mathcal{F}}_3$ is the subset formed by all triangles $T\in {\mathcal{T}}_h$ such that $\overline{T}\cap \overline{Q_h}\ne \mathrm{\varnothing }$ and $T\not\subset \overline{Q_h}$.

• ${\mathcal{F}}_4={\mathcal{T}}_h\setminus ({\mathcal{F}}_1\cup {\mathcal{F}}_2\cup {\mathcal{F}}_3)$.

Figure 1

In this figure, we see the splitting of the fixed triangulation related to the position of the solid at time t .

We introduce the approximated density function ${\rho }_h^k$ as follows:

$${\rho }_h^k(\mathbf{x})=\{\begin{array}{cc}{\rho }_S^{\ast }({\mathcal{R}}_{-{\theta }_h^k}(\mathbf{x}-{\xi }_h^k),t^k)\hfill & \text{if }\mathbf{x}\in \mathcal{S}({\mathit{\xi }}_h^k,{\theta }_h^k,t^k),\hfill \\ {\rho }_F\hfill & \text{if }\mathbf{x}\in \mathrm{\Omega }\setminus \mathcal{S}({\mathit{\xi }}_h^k,{\theta }_h^k,t^k).\hfill \end{array}$$

(60)

With these notations, we introduce the following mixed variational fully discrete formulation: Find ${\mathbf{u}}_h^{k+1}\in {\mathrm{\Lambda }}_h^{k+1}+{\mathcal{K}}_h({\mathit{\xi }}_h^{k+1},{\theta }_h^{k+1},t^{k+1})$, $p_h^{k+1}\in {\mathcal{M}}_h({\mathit{\xi }}_h^{k+1},{\theta }_h^{k+1},t^{k+1})$ such that

$$\begin{array}{c}({\rho }_h^{k+1}\frac{{\mathbf{u}}_h^{k+1}-{\mathbf{u}}_h^k\circ {\overline{\mathbf{X}}}_h^k}{\mathrm{\Delta }t},\mathbf{v})+a({\mathbf{u}}_h^{k+1},\mathbf{v})+b(\mathbf{v},p_h^{k+1})\hfill \\ =({\rho }_h^{k+1}{\mathbf{f}}_h^{k+1},\mathbf{v})\mathrm{\forall }\mathbf{v}\in {\mathcal{K}}_h({\mathit{\xi }}_h^{k+1},{\theta }_h^{k+1},t^{k+1}),\hfill \end{array}$$

(61)

$$b({\mathbf{u}}_h^{k+1},q)=0\mathrm{\forall }q\in {\mathcal{M}}_h({\mathit{\xi }}_h^{k+1},{\theta }_h^{k+1},t^{k+1}),$$

(62)

where ${\mathbf{f}}_h^{k+1}$ is the $L^2{(\mathrm{\Omega })}^2$-projection of ${\mathbf{f}}^{k+1}=\mathbf{f}(t^{k+1})$ on ${(E_h)}^2$.

Let us now state the second main result of this paper, which asserts the convergence of the fully-discrete scheme (61)-(62). The complete proof of this result could be found in [18] for the case of rigid body and in the forthcoming paper [22] if the structure is deformable.

Theorem 4.2 Let Ω be a convex domain with a polygonal boundary. Suppose that f and ${\mathbf{u}}_0$ satisfy the conditions from (36), and that $(\mathbf{u},p,\mathit{\xi },\theta )$ is a solution of (13)-(18) satisfying regularity properties (37)-(39). Let $C_0>0$ and $0<\alpha \le 1$ be two fixed constants. Then there exist two positive constants K and ${\tau }^{\ast }$, independent of h and Δt such that for all $0<\mathrm{\Delta }t\le {\tau }^{\ast }$ and for all $h\le C_0\mathrm{\Delta }{t}^{1+\alpha }$, we have

$$\underset{1\le k\le N}{sup}({\parallel \mathbf{u}\left(t^k\right)-{\mathbf{u}}_h^k\parallel }_{L^2{(\mathrm{\Omega })}^2}+|\mathit{\xi }\left(t^k\right)-{\mathit{\xi }}_h^k|+|\theta \left(t^k\right)-{\theta }_h^k|)\le K\mathrm{\Delta }{t}^{\alpha }.$$

Let us mention that in order to get an approximation of first order in time (i.e., $O(\mathrm{\Delta }t)$ in Theorem 4.2), we have to choose $\alpha =1$. In this case, the corresponding condition on h becomes $h\le C_0\mathrm{\Delta }{t}^2$ which is similar to the one obtained in [[2], Theorem 3.2], where the densities of the fluid ${\rho }_F$ and of the solid ${\rho }_S$ are equal.

Remark 4.3 Let us give some comments on the condition of h and Δt required for the convergence result in Theorem 4.2. First, we emphasize that the same type of condition appears in several works for approximation in a Lagrangian framework of the Navier-Stokes equations without any rigid body. We may cite [10], where convergence is obtained under condition $h\le C_0\mathrm{\Delta }t$ and [11], where h and Δt are chosen such that $h^2\le C\mathrm{\Delta }t\le C_1{h}^{\sigma }$ and $\sigma >1/2$ (with h and Δt small enough). We also mention [14] for an ALE scheme applied to Stokes equations in a time-dependent domain, where the authors obtain an error estimate of order $\mathcal{O}(\mathrm{\Delta }t)$ under condition $h\le C\mathrm{\Delta }{t}^{3/4}$.