A One-Field Energy-conserving Monolithic Fictitious Domain Method for Fluid-Structure Interactions

TL;DR

This paper introduces a novel one-field, energy-conserving fictitious domain method for fluid-structure interactions that simplifies computations by solving a single velocity field while maintaining stability and accuracy.

Contribution

The paper presents a new implicit fictitious domain method that conserves energy, decouples fluid and solid interactions, and provides stability and error analysis.

Findings

Energy conservation is proven for the method.

The method demonstrates stability and accuracy in 2D and 3D simulations.

Numerical examples validate the conservation properties.

Abstract

In this article, we analyze and numerically assess a new fictitious domain method for fluid-structure interactions in two and three dimensions. The distinguishing feature of the proposed method is that it only solves for one velocity field for the whole fluid-structure domain; the interactions remain decoupled until solving the final linear algebraic equations. To achieve this the finite element procedures are carried out separately on two different meshes for the fluid and solid respectively, and the assembly of the final linear system brings the fluid and solid parts together via an isoparametric interpolation matrix between the two meshes. In this article, an implicit version of this approach is introduced. The property of energy conservation is proved, which is a strong indication of stability. The solvability and error estimate for the corresponding stationary problem (one time step of the transient problem) are analyzed. Finally, 2D and 3D numerical examples are presented to validate the conservation properties.

Figures

Tables

Table 1. Table 1: Parameter sets for test problem of the oscillating disc.

Parameter sets	${\rho }^f$	${\rho }^s$	${\nu }^f$	${\nu }^s$	${\mu }^s$
Parameter 1	$1$	$1$	$0.01$	$0.01$	$1$
Parameter 2	$1$	$10$	$0.01$	$0$	$10$

Equations

$$\rho \frac{d\mathbf{u}}{dt}=\nabla \cdot \mathbf{\sigma },$$

$$\nabla \cdot \mathbf{u}=0.$$

$$\mathbf{\sigma }={\mathbf{\sigma }}^f={\nu }^f\mathrm{D}{\mathbf{u}}^f-p^f\mathbf{I},$$

$$\mathbf{\sigma }={\mathbf{\sigma }}^s={\mathbf{\tau }}^s+{\nu }^s\mathrm{D}{\mathbf{u}}^s-p^s\mathbf{I},$$

$${\mathbf{\tau }}^s=J^{-1}\mathbf{P}{\mathbf{F}}^{\mathrm{T}},$$

$${\mathbf{u}}^f={\mathbf{u}}^{s}on{\Gamma }_t,$$

$${\mathbf{n}}^s{\mathbf{\sigma }}^f={\mathbf{n}}^s{\mathbf{\sigma }}^{s}on{\Gamma }_t,$$

$${\mathbf{u}}^f=0on\Gamma ,$$

$${{\mathbf{u}}^f|}_{t=0}={\mathbf{u}}_0^f,$$

$${{\mathbf{u}}^s|}_{t=0}={\mathbf{u}}_0^s.$$

$${\int }_{{\Omega }_t^f}Eq.\text{lx@cref}creftyperefnummomentu{m}_{e}quation\cdot \mathbf{v}d\mathbf{x}+{\int }_{{\Omega }_t^s}Eq.\text{lx@cref}creftyperefnummomentu{m}_{e}quation\cdot \mathbf{v}d\mathbf{x}.$$

$$\begin{array}{rl}& {\rho }^f{\int }_{\Omega }\frac{d\mathbf{u}}{dt}\cdot \mathbf{v}d\mathbf{x}+\frac{{\nu }^f}{2}{\int }_{\Omega }\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{v}d\mathbf{x}-{\int }_{\Omega }p\nabla \cdot \mathbf{v}d\mathbf{x}\\ & +{\rho }^{\delta }{\int }_{{\Omega }_t^s}\frac{d\mathbf{u}}{dt}\cdot \mathbf{v}d\mathbf{x}+\frac{{\nu }^{\delta }}{2}{\int }_{{\Omega }_t^s}\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{v}d\mathbf{x}+{\int }_{{\Omega }_t^s}{\mathbf{\tau }}^s:\nabla \mathbf{v}d\mathbf{x}=0,\end{array}$$

$$-{\int }_{{\Omega }_t^f}Eq.\text{lx@cref}creftyperefnumcontinuit{y}_{e}quationqd\mathbf{x}-{\int }_{{\Omega }_t^s}Eq.\text{lx@cref}creftyperefnumcontinuit{y}_{e}quationqd\mathbf{x},$$

$$\begin{array}{rl}& {\rho }^f{\int }_{\Omega }\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{v}d\mathbf{x}+{\rho }^f{\int }_{\Omega }\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}\cdot \mathbf{v}d\mathbf{x}+\frac{{\nu }^f}{2}{\int }_{\Omega }\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{v}d\mathbf{x}-{\int }_{\Omega }p\nabla \cdot \mathbf{v}d\mathbf{x}\\ & +{\rho }^{\delta }{\int }_{{\Omega }_t^s}\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{v}d\mathbf{x}+\frac{{\nu }^{\delta }}{2}{\int }_{{\Omega }_t^s}\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{v}d\mathbf{x}+{\int }_{{\Omega }_{\mathbf{X}}^s}\mathbf{P}(\mathbf{F}):{\nabla }_{\mathbf{X}}\mathbf{v}d\mathbf{X}=0,\end{array}$$

$$-{\int }_{\Omega }q\nabla \cdot \mathbf{u}d\mathbf{x}=0,$$

$${\int }_0^t{\int }_{{\Omega }_{\mathbf{X}}^s}\mathbf{P}:{\nabla }_{\mathbf{X}}\mathbf{u}d\mathbf{X}={\int }_{{\Omega }_{\mathbf{X}}^s}\Psi (\mathbf{F})\mathbf{d}\mathbf{X}$$

$$\begin{array}{r}\frac{d}{dt}{\int }_{{\Omega }_{\mathbf{X}}^s}\Psi (\mathbf{F})d\mathbf{X}={\int }_{{\Omega }_{\mathbf{X}}^s}\frac{\partial \Psi }{\partial \mathbf{F}}:\frac{d\mathbf{F}}{dt}d\mathbf{X}={\int }_{{\Omega }_{\mathbf{X}}^s}\mathbf{P}:\frac{d}{dt}\left(\mathbf{I}+{\nabla }_{\mathbf{X}}\mathbf{d}\right)d\mathbf{X}={\int }_{{\Omega }_{\mathbf{X}}^s}\mathbf{P}:{\nabla }_{\mathbf{X}}\mathbf{u}d\mathbf{X},\end{array}$$

$${\int }_{\Omega }\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}\cdot \mathbf{u}d\mathbf{x}=0.$$

$${\int }_{\Omega }\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}\cdot \mathbf{u}d\mathbf{x}={\int }_{\Omega }\nabla \left(\mathbf{u}\otimes \mathbf{u}\right)\cdot \mathbf{u}d\mathbf{x}-{\int }_{\Omega }{\left|\mathbf{u}\right|}^2\nabla \cdot \mathbf{u}.$$

$${\int }_{\Omega }\nabla \left(\mathbf{u}\otimes \mathbf{u}\right)\cdot \mathbf{u}d\mathbf{x}={\int }_{\Gamma }{\left|\mathbf{u}\right|}^2\mathbf{u}\cdot \mathbf{n}-{\int }_{\Omega }\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}\cdot \mathbf{u}d\mathbf{x}.$$

$${\int }_{\Gamma }{\left|\mathbf{u}\right|}^2\mathbf{u}\cdot \mathbf{n}={\int }_{\Omega }{\left|\mathbf{u}\right|}^2\nabla \cdot \mathbf{u}=0,$$

$$\begin{array}{rl}& \frac{{\rho }^f}{2}{\int }_{\Omega }|\mathbf{u}{|}^{2}d\mathbf{x}+\frac{{\nu }^f}{2}{\int }_0^t{\int }_{\Omega }\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{u}d\mathbf{x}\\ +& \frac{{\rho }^{\delta }}{2}{\int }_{{\Omega }_t^s}|\mathbf{u}{|}^{2}d\mathbf{x}+\frac{{\nu }^{\delta }}{2}{\int }_0^t{\int }_{{\Omega }_t^s}\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{u}d\mathbf{x}+{\int }_{{\Omega }_{\mathbf{X}}^s}\Psi (\mathbf{F})\mathbf{d}\mathbf{X}=0.\end{array}$$

$$\begin{array}{rl}& {\rho }^f{\int }_{\Omega }\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}\cdot \mathbf{v}d\mathbf{x}+{\rho }^f{\int }_{\Omega }\left({\mathbf{u}}_{\ast }\cdot \nabla \right){\mathbf{u}}_{\ast }\cdot \mathbf{v}d\mathbf{x}+\frac{{\nu }^f}{2}{\int }_{\Omega }\mathrm{D}{\mathbf{u}}_{\ast }:\mathrm{D}\mathbf{v}d\mathbf{x}\\ & -{\int }_{\Omega }{p}_{n+1}\nabla \cdot \mathbf{v}d\mathbf{x}+{\rho }^{\delta }{\int }_{{\Omega }_{n+1}^s}\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}\cdot \mathbf{v}d\mathbf{x}+\frac{{\nu }^{\delta }}{2}{\int }_{{\Omega }_{n+1}^s}\mathrm{D}{\mathbf{u}}_{\ast }:\mathrm{D}\mathbf{v}d\mathbf{x}\\ & +{\int }_{{\Omega }_{\mathbf{X}}^s}\mathbf{P}\left({\mathbf{F}}_{\ast }\right):{\nabla }_{\mathbf{X}}\mathbf{v}d\mathbf{X}=O\left(\Delta t^2\right),\end{array}$$

$$-{\int }_{\Omega }q\nabla \cdot {\mathbf{u}}_{\ast }d\mathbf{x}=0,$$

$$\Delta t\mathbf{P}\left({\mathbf{F}}_{\ast }\right):{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\ast }+O\left(\Delta t^3\right)=\Psi ({\mathbf{F}}_{n+1})-\Psi ({\mathbf{F}}_n).$$

$$\Delta t{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\ast }={\nabla }_{\mathbf{X}}{\mathbf{x}}_{n+1}-{\nabla }_{\mathbf{X}}{\mathbf{x}}_n={\mathbf{F}}_{n+1}-{\mathbf{F}}_n,$$

$$\begin{array}{rl}& \Psi ({\mathbf{F}}_{n+1})-\Psi ({\mathbf{F}}_n)=w\left(\Delta t\right)-w(0)=\Delta t{w}^{\prime }\left(\frac{\Delta t}{2}\right)+O\left(\Delta t^3\right).\end{array}$$

$$w^{\prime }\left(\frac{\Delta t}{2}\right)={\frac{\partial \Psi }{\partial \mathbf{F}}|}_{\xi =\frac{\Delta t}{2}}:{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\ast }=\mathbf{P}\left({\mathbf{F}}_{\ast }\right):{\nabla }_{\mathbf{X}}{\mathbf{u}}_{\ast }.$$

$${\int }_{\Omega }\left({\mathbf{u}}_{\ast }\cdot \nabla \right){\mathbf{u}}_{\ast }\cdot {\mathbf{u}}_{\ast }d\mathbf{x}=0.$$

$$\begin{array}{rl}& \frac{{\rho }^f}{2}{\int }_{\Omega }{\left|{\mathbf{u}}_{n+1}\right|}^{2}d\mathbf{x}+\frac{{\rho }^{\delta }}{2}{\int }_{{\Omega }_{n+1}^s}{\left|{\mathbf{u}}_{n+1}\right|}^{2}d\mathbf{x}+{\int }_{{\Omega }_{\mathbf{X}}^s}\Psi \left({\mathbf{F}}_{n+1}\right)d\mathbf{X}\\ & +\frac{\Delta t{\nu }^f}{2}{\int }_{\Omega }\mathrm{D}{\mathbf{u}}_{\ast }:\mathrm{D}{\mathbf{u}}_{\ast }d\mathbf{x}+\frac{\Delta t{\nu }^{\delta }}{2}{\int }_{{\Omega }_{n+1}^s}\mathrm{D}{\mathbf{u}}_{\ast }:\mathrm{D}{\mathbf{u}}_{\ast }d\mathbf{x}\\ & =\frac{{\rho }^f}{2}{\int }_{\Omega }{\left|{\mathbf{u}}_n\right|}^{2}d\mathbf{x}+\frac{{\rho }^{\delta }}{2}{\int }_{{\Omega }_n^s}{\left|{\mathbf{u}}_n\right|}^{2}d\mathbf{x}+{\int }_{{\Omega }_{\mathbf{X}}^s}\Psi \left({\mathbf{F}}_n\right)d\mathbf{X}+O\left(\Delta t^3\right).\end{array}$$

Full text

{AMS}

97N40, 74S05, 76M10

1 Introduction

Three major questions arise when considering a finite element method for the problems of Fluid-Structure Interactions (FSI): (1) what kind of meshes are used (interface fitted or unfitted); (2) how to couple the fluid-structure interactions (monolithic/fully-coupled or partitioned/segregated); (3) what variables are solved (velocity and/or displacement). Combinations of the answers of these questions lead to different types of numerical methods. For examples, [8, 16] solve for fluid velocity and solid displacement sequentially (partitioned/segregated) using an Arbitrary Lagrangian-Eulerian (ALE) fitted mesh. Whereas [12, 13, 17] use an ALE fitted mesh to solve for fluid velocity and solid displacement simultaneously (monolithic/fully-coupled) with a Lagrange Multiplier to enforce the continuity of velocity/displacement on the interface. The Immersed Finite Element Method (IFEM) [5, 18, 21, 22, 23, 26, 27] and the Fictitious Domain Method (FDM) [3, 6, 10, 14, 15, 25] use two meshes to represent the fluid and solid separately. Although IFEM could be monolithic [5], the classical IFEM only solves for velocity, while the solid information is arranged on the right-hand side of the fluid equation as a known force term. Although the FDM may be partitioned [25], usually the FDM approach solves for both velocity in the whole domain (fluid plus solid) and displacement of the solid simultaneously via a distributed Lagrange multiplier (DLM) to enforce the consistency of velocity/displacement in the overlapped solid domain.

In the case of one-field and monolithic numerical methods for FSI problems, [2] introduces a 1D model using a one-field FD formulation based on two meshes, and [11, 19] introduces an energy stable monolithic method (in 2D) based on one Eulerian mesh and discrete remeshing. In a previous study [24], we present a one-field monolithic fictitious domain method which is straightforward to implement in both 2D and 3D. The main features of this method are: (1) only one velocity field is solved in the whole domain, based upon the use of an appropriate $L^{2}$ projection; (2) the fluid and solid equations are solved monolithically. In this paper, we construct an implicit version of this one-field fictitious domain and, developing the ideas from [2, 19], we analyze properties of energy stability and error estimate of the proposed approach.

The paper is organized as follows. Control equations are introduced in Section 2. Sections 3 and 4 discuss the weak form and the energy estimate in the continuous case respectively. Sections 5 and 6 discuss the weak form and the energy estimate after time discretization respectively. In Section 7, the corresponding stationary problem (one time step of the transient problem) is analyzed after time discretization and further after space discretization. In Section 8, implementation details are presented. Numerical examples are given in Section 9, and conclusions are presented in Section 10.

2 Control equations

In the following context, $\Omega_{t}^{f}\subset\mathbb{R}^{d}$ and $\Omega_{t}^{s}\subset\mathbb{R}^{d}$ with $d=2,3$ denote the fluid and solid domain respectively which are time dependent regions as shown in Fig. 1. $\Omega=\Omega_{t}^{f}\cup\Omega_{t}^{s}$ is a fixed domain (with outer boundary $\Gamma$) and $\Gamma_{t}=\partial\Omega_{t}^{f}\cap\partial\Omega_{t}^{s}$ is the moving interface between fluid and solid. We denote by ${\bf X}$ the reference (material) coordinates of the solid, by ${\bf x}={\bf x}(\cdot,t)$ the current coordinates of the solid, and by ${\bf x}_{0}$ is the initial coordinates of the solid. We assume ${{\bf x}(t):\Omega_{\bf X}^{s}\to\Omega_{t}^{s}}$ is one-to-one and invertible with Lipschitz inverse, i.e.: for all $t\in\left[0,T\right]$ and ${\bf s}_{1},{\bf s}_{2}\in\Omega_{\bf X}^{s}$, $\exists$ $c_{0}>0$ such that $\|{{\bf x}({\bf s}_{1},t)}-{{\bf x}({\bf s}_{2},t)}\|_{2}\geq c_{0}\|{\bf s}_{1}-{\bf s}_{2}\|_{2}$.

Let $\rho,{\bf u},{\bm{\sigma}}$ denote the density, velocity and stress tensor respectively. We assume both an incompressible fluid and incompressible solid, then the conservation of momentum and conservation of mass take the same form as follows:

Momentum equation:

[TABLE]

Continuity equation:

[TABLE]

Let superscripts $f$ and $s$ refer to the fluid and solid respectively, and ${\rm D}{\bf u}=\nabla{\bf u}+\nabla^{\rm T}{\bf u}$, then the constitutive equations may be expressed as follows:

Incompressible Newtonian fluid in $\Omega^{f}$:

[TABLE]

Incompressible viscous-hyperelastic solid in $\Omega^{s}$ [6]:

[TABLE]

where

[TABLE]

with ${\bf P}=\frac{\partial\Psi\left({\bf F}\right)}{\partial{\bf F}}$ (${\bf P}_{ij}=\frac{\partial\Psi\left({\bf F}\right)}{\partial{F}_{ij}}$) being the first Piola-Kirchhoff stress tensor and $\Psi\left({\bf F}\right)$ being the energy function for a hyperelastic solid material. $J=det_{\bf F}$ is the determinant of ${\bf F}$, where ${\bf F}=\frac{\partial{\bf x}}{\partial{\bf X}}=\frac{\partial{\bf x}}{\partial{\bf x}_{0}}\frac{\partial{\bf x}_{0}}{\partial{\bf X}}$=$\nabla_{0}{\bf x}\nabla_{\bf X}{\bf x}_{0}$ is deformation tensor of the solid. In the above, $\nu^{f}$ and $\nu^{s}$ are viscosity of the fluid and solid respectively, $p^{f}$ and $p^{s}$ are pressure of the fluid and solid respectively.

The system is complemented with the following boundary and initial conditions.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Other boundary conditions are possible on $\Gamma$ but Eq. 8 are used here for simplicity.

3 Weak form on the continuous level

In the following context, let $L^{2}(\omega)$ be the square integrable functions in domain $\omega$, endowed with norm $\left\|u\right\|_{0,\omega}^{2}=\int_{\omega}\left|u\right|^{2}$ ($u\in L^{2}(\omega)$). Let $H^{1}(\omega)=\left\{u:u,\nabla u\in L^{2}(\omega)\right\}$ with the norm denoted by $\left\|u\right\|_{1,\omega}^{2}=\left\|u\right\|_{0,\omega}^{2}+\left\|\nabla u\right\|_{0,\omega}^{2}$. We also denote by $H_{0}^{1}(\omega)$ the subspace of $H^{1}(\omega)$ whose functions have zero values on the boundary of $\omega$, and denote by $L_{0}^{2}(\omega)$ the subspace of $L^{2}(\omega)$ whose functions have zero mean value.

Let ${p}=\left\{\begin{matrix}{{p}^{f}\quad in\quad\Omega^{f}}\\ {{p}^{s}\quad in\quad\Omega^{s}}\\ \end{matrix}\right.$. Given ${\bf v}\in H_{0}^{1}(\Omega)^{d}$, we perform the following symbolic operations:

[TABLE]

Integrating the stress terms by parts, the above operations, using constitutive equation Eq. 3 and Eq. 4 and boundary condition Eq. 7, gives:

[TABLE]

where $\rho^{\delta}=\rho^{s}-\rho^{f}$, $\nu^{\delta}=\nu^{s}-\nu^{f}$. Note that the integrals on the interface $\Gamma_{t}$ are cancelled out using boundary condition Eq. 7. This is not surprising because they are internal forces for the whole FSI system considered here.

Substitute the expression of ${{\bm{\tau}}^{s}}$ Eq. 5 into Eq. 11 and transfer the integral of the last term to the reference coordinate system. Then, the following symbolic operations for $q\in L^{2}(\Omega)$,

[TABLE]

lead to the weak form of the FSI system as follows.

Problem 1**.**

Find ${\bf u}(t)\in H_{0}^{1}(\Omega)^{d}$ and $p(t)\in L_{0}^{2}(\Omega)$, such that

[TABLE]

$\forall{\bf v}\in H_{0}^{1}(\Omega)^{d}$* and*

[TABLE]

$\forall q\in L^{2}(\Omega)$*. *

Remark 1**.**

*Because domain $\Omega$ is stationary (the Eulerian description will be used) and $\Omega_{t}^{s}$ is transient which will be updated by its own velocity (the updated Lagrangian description), there is a convection term from the total derivative of time in $\Omega$, but there is no convection term in $\Omega_{t}^{s}$. *

4 Energy conservation on the continuous level

A property of energy conservation for the weak forms Eq. 12 and Eq. 13 will be proved in this section.

Lemma 4.1**.**

Assume the solid energy function $\Psi(\cdot)\in C^{1}$ over the set of second order tensors, then

[TABLE]

Proof 4.2**.**

Using the fact ${\bf A:C}=tr_{{\bf AC}^{\rm T}}$ (${\bf A}$ and ${\bf C}$ are arbitrary matrices), we have:

[TABLE]

*where ${\bf d}$ is displacement of the solid at time $t$ in the above. *

Lemma 4.3**.**

If $\left({\bf u},p\right)$ is the solution pair of 1, then

[TABLE]

Proof 4.4**.**

First,

[TABLE]

Integrate by parts:

[TABLE]

According to the boundary condition Eq. 8 and equation Eq. 13, we have

[TABLE]

*which further indicates equation Eq. 15. *

Proposition 4.5** (Energy Conservation on the Continuous Level).**

Let $\left({\bf u},p\right)$ be the solution pair of 1, then

[TABLE]

Proof 4.6**.**

*first choose ${\bf v}={\bf u}$ in Eq. 12 and integrate from time [math] to $t$, then let $q=p$ in Eq. 13 and substitute into Eq. 12. Finally we can construct the above equation of energy balance due to Lemma 4.1 and Lemma 4.3. *

5 Weak form after discretization in time

Using the Crank-Nicolson scheme to discretize equation Eq. 12 in time, 1 becomes:

Problem 5.1**.**

For each time step, find ${\bf u}_{n+1}\in H_{0}^{1}(\Omega)^{d}$ and $p_{n+1}\in L_{0}^{2}(\Omega)$, such that

[TABLE]

$\forall{\bf v}\in H_{0}^{1}\left(\Omega\right)^{d}$* and*

[TABLE]

$\forall q\in L^{2}(\Omega)$*, where ${\bf u}_{*}=\frac{{\bf u}_{n+1}+{\bf u}_{n}}{2}$ and ${\bf F}_{*}=\frac{{\bf F}_{n+1}+{\bf F}_{n}}{2}$. *

Remark 5.2**.**

*Notice that the subscript $n+1$ indicates that $\Omega_{n+1}^{s}$ is transient but $\Omega$ is stationary. $\Omega_{n+1}^{s}$ is updated from $\Omega_{n}^{s}$ by $\Omega_{n+1}^{s}=\left\{{\bf x}_{n+1}:{\bf x}_{n+1}={\bf x}_{n}+{\bf u}_{*}\Delta t\right\}$, for all ${\bf x}_{n}\in\Omega_{n}^{s}$. *

6 Energy conservation after discretization in time

Energy conservation for the weak form Eq. 17 after time discretization will be proved in this section.

Lemma 6.1**.**

Assume $\Psi(\cdot)\in C^{1}$ over the set of second order tnesors, then

[TABLE]

Proof 6.2**.**

Let $w(\xi)=\Psi\left({\bf F}_{n}+\xi\nabla_{\bf X}{\bf u}_{*}\right)$, and notice that

[TABLE]

then

[TABLE]

Using the chain rule, Eq. 19 holds thanks to

[TABLE]

Similarly to Lemma 4.3, we have:

Lemma 6.3**.**

If $\left({\bf u}_{n+1},p_{n+1}\right)$ is the solution pair of 5.1, then

[TABLE]

Proposition 6.4** (Local Energy Conservation).**

Let $\left({\bf u}_{n+1},p_{n+1}\right)$ be the solution pair of 5.1, then

[TABLE]

Proof 6.5**.**

First, let ${\bf v}={\bf u}_{n+1}$ in Eq. 17 and $q=p_{n+1}$ in Eq. 18, to get the quation:

[TABLE]

then let ${\bf v}={\bf u}_{n}$ in Eq. 17 and $q=p_{n+1}$ in Eq. 18, to get the further equation:

[TABLE]

*Add the above two equations, we finally have Proposition 6.4 due to Lemma 6.1 and Lemma 6.3. *

Corollary 6.6** (Total Energy Conservation).**

*Let $N+1=t/\Delta t$, where $t$ is the computational time, and use the following notation for the different contributions to the total energy at the time t:

• •

Kinetic energy in $\Omega$:

[TABLE]

• •

Kinetic energy in $\Omega_{t}^{s}$:

[TABLE]

• •

Viscous dissipation in $\Omega$:

[TABLE]

• •

Viscous dissipation in $\Omega_{t}^{s}$:

[TABLE]

• •

Potential energy of solid:

[TABLE]

Denote the total energy at time $t$ as:

[TABLE]

and the error of total energy as:

[TABLE]

then,

[TABLE]

Proof 6.7**.**

*First let $O(\Delta t^{3})=\sum_{k\geq 3}C_{n+1}^{k}\Delta t^{k}$ in Eq. 21, where $C_{n+1}^{k}$ is a constant dependent on the specific time step. Then add equation Eq. 21 from $n=0$ to $n=N$, we have $Err=\sum_{k\geq 3}\Delta t^{k}\sum_{n=0}^{N}C_{n+1}^{k}$. Eq. 30 holds due to $\sum_{n=0}^{N}\left|C_{n+1}^{k}\right|\leq C_{max}^{k}t/\Delta t$, where $C_{max}^{k}=max\{\left|C_{n+1}^{k}\right|,n=0,\cdots,N\}$. *

Remark 6.8**.**

*Note that, in subsequent sections, the viscous term in Eq. 26 and Eq. 27 will introduce an error $\Delta tO(h^{m})$ ($m>0$) after discretization in space, where $h$ represents the mesh size. In that case it will be seen that the error of total energy will be reduced to $O(\Delta t)$ for a fixed mesh. *

7 Analysis of the stationary problem corresponding to 5.1

We now focus on the stationary problem corresponding to one step of the 5.1 and consider its well-posedness and discretization in space. In order to make the analysis tractable we make the following simplifying assumptions, equivalent to those presented in [6].

• •

Neglect the convection term.

• •

Assume $\rho^{\delta}\geq 0$ and $\nu^{\delta}\geq 0$.

• •

Assume a linear model for ${\bf P}$, i.e., ${\bf P}=\mu^{s}{\bf F}$.

Remark 7.1**.**

*We have implemented numerical examples which include convection and the cases of $\rho^{\delta}<0$ or $\nu^{\delta}<0$ (see Section 9, and [24] for more numerical tests). The linear assumption for ${\bf P}$ is important in order to define the following bilinear form Eq. 33, and we have only implemented an incompressible neo-Hookean model in this linear case. For non-linear cases, it may be that linearization and/or modification of 5.1 are required. This will be a topic of our further investigation. *

Under the assumption of a linear model for ${\bf P}$, we have

[TABLE]

Define the following bilinear forms:

[TABLE]

[TABLE]

and

[TABLE]

where $\alpha=2\rho^{f}/\Delta t$, $\beta=2\rho^{\delta}/\Delta t$ and $\gamma=\mu^{s}\Delta t/2$. We also define the following linear forms:

[TABLE]

and

[TABLE]

Let

[TABLE]

and

[TABLE]

then the weak form corresponding to one step of 5.1 with ${\bf u}={\bf u}_{*}$ and $p=p_{n+1}$, using the above notation and assumptions, can be stated as:

Problem 7.2**.**

Find ${\bf u}\in\mathbb{V}$, $p\in\mathbb{P}$, such that

[TABLE]

*where $\mathbb{V}=H_{0}^{1}\left(\Omega\right)^{d}$ and $\mathbb{P}=L_{0}^{2}(\Omega)$. *

We shall use the following norms for vector functions and matrix functions respectively: $\left\|{\bf u}\right\|_{0,\omega}^{2}=\sum_{i=1}^{d}\left\|u_{i}\right\|_{0,\omega}^{2}$ and $\left\|{\bf A}\right\|_{0,\omega}^{2}=\sum_{i=1}^{d}\sum_{j=1}^{d}\left\|A_{ij}\right\|_{0,\omega}^{2}$, and denote $\|{\bf u}\|_{\mathbb{V}}=\|{\bf u}\|_{1,\Omega}$.

7.1 Well-posedness of 7.2

Lemma 7.3**.**

$a({\bf u},{\bf v})$* is coercive, i.e., there exists a positive real number $c_{a}$, such that for $\forall{\bf u}\in\mathbb{V}$,*

[TABLE]

Proof 7.4**.**

First, according to definitions of $\left\|{\rm D}{\bf u}\right\|_{0,\Omega}$ and $\left\|{\bf u}\right\|_{1,\Omega}$,

[TABLE]

and

[TABLE]

the above estimate is due to the Poincar$\rm\acute{e}$ inequality $\left\|{\bf u}\right\|_{0,\Omega}\leq C_{\Omega}\left\|\nabla{\bf u}\right\|_{0,\Omega}$, where $C_{\Omega}$ is constant. Then,

[TABLE]

*using our assumption $\beta=2\rho^{\delta}/\Delta t\geq 0$ and $\nu^{\delta}\geq 0$. *

Lemma 7.5**.**

$a({\bf u},{\bf v})$* is bounded, i.e., there exists a positive real number $C_{a}$, such that for $\forall{\bf u},{\bf v}\in{\mathbb{V}}$,*

[TABLE]

Proof 7.6**.**

First

[TABLE]

Noticing that

[TABLE]

and

[TABLE]

*we can observe Eq. 40 holds from LABEL:proof_1_1, with $C_{a}=\alpha+\beta+2\nu^{f}+2\nu^{\delta}+\gamma$. *

Lemma 7.7**.**

$\ell({\bf v})$* is bounded, i.e., there exists a positive real number $C_{\ell}$, such that for $\forall{\bf v}\in{\mathbb{V}}$,*

[TABLE]

Proof 7.8**.**

[TABLE]

*Notice that $\left\|\nabla_{\bf X}{\bf x}_{n}\right\|_{0,\Omega_{\bf X}^{s}}$ is bounded owing to the assumption that ${\bf x}={\bf x}({\bf s},t)$ (${\bf s}\in\Omega_{\bf X}^{s}$) is Lipschitz invertible. We further get Eq. 43 using relation Eq. 42. *

Proposition 7.9** (Inf-Sup Condition).**

There exist $c_{b}>0$ such that

[TABLE]

Proof 7.10**.**

*can be found in [7, Section 12.2]. *

The above three lemmas and the inf-sup condition imply [7, Lemma 12.2.12]:

Theorem 7.11**.**

7.2* has a unique solution $({\bf u},p)\in\mathbb{V}\times\mathbb{P}$. *

7.2 Discretization in space

We shall use a fixed Eulerian mesh for $\Omega$ and an updated Lagrangian mesh for $\Omega_{n+1}^{s}$ to discretize 7.2. First, we discretize $\Omega$ as $\Omega^{h}$ with the corresponding finite element spaces as

[TABLE]

and

[TABLE]

The approximated solution ${\bf u}^{h}$ and $p^{h}$ can be expressed in terms of these basis functions as

[TABLE]

We further discretize $\Omega_{n+1}^{s}$ as $\Omega_{n+1}^{sh}$ with the corresponding finite element spaces as:

[TABLE]

and approximate $\left.{\bf u}^{h}({\bf x})\right|_{{\bf x}\in\Omega_{n+1}^{sh}}$ as:

[TABLE]

where ${\bf x}_{i}^{s}$ is the nodal coordinate of the solid mesh. Notice that the above approximation defines an $L^{2}$ projection $P_{n+1}$ from $V^{h}\left(\Omega^{h}\right)^{d}$ to $V^{sh}\left(\Omega_{n+1}^{sh}\right)^{d}$: $P_{n+1}\left({\bf u}^{h}({\bf x})\right)={\bf u}^{sh}\left({\bf x}\right)$.

Remark 7.12**.**

*The solid domain $\Omega_{t}^{s}$ is actually discretized once on $\Omega_{0}^{s}$, and then $\Omega_{n}^{s}(n\geq 1)$ is updated from the mesh at the previous time step as illustrated in Remark 5.2. *

We then discretize 7.2 as follows.

Problem 7.13**.**

Find ${\bf u}^{h}\in\mathbb{V}^{h}$, $p^{h}\in\mathbb{P}^{h}$, such that

[TABLE]

where $\mathbb{V}^{h}=V^{h}(\Omega^{h})^{d}$, $\mathbb{P}^{h}=L^{h}(\Omega^{h})$, $a^{h}$ can be expressed as:

[TABLE]

and $\ell^{h}$ can be expressed as:

[TABLE]

*where ${\bf u}^{sh}=P_{n+1}\left({\bf u}^{h}\right)$ and ${\bf v}^{sh}=P_{n+1}\left({\bf v}^{h}\right)$. *

From the error estimate for interpolation [7, Chapter 4, Corollary 4.4.24], the approximation Eq. 46 can be bounded by:

[TABLE]

where $C_{h}$ is a constant depending on the element of $\Omega^{sh}$. In order to prove 7.13 is well-posed, we still need to prove the continuity of $a^{h}({\bf u}^{h},{\bf v}^{h})$ and $\ell^{h}({\bf v}^{h})$, and the coercivity of $a^{h}({\bf u}^{h},{\bf v}^{h})$. The main idea is to bound ${\bf u}^{sh}$, $\nabla_{\bf X}{\bf u}^{sh}$ and ${\rm D}{\bf u}^{sh}$ in $\mathbb{V}$, which can be achieved by the above interpolation error estimate.

Lemma 7.14**.**

There exists a positive real number $c_{1}$ such that

[TABLE]

Proof 7.15**.**

From Eq. 49, we have:

[TABLE]

*We have Eq. 50 due to $\left\|{\bf u}^{sh}({\bf x})\right\|_{0,\Omega_{n+1}^{s}}\leq\left\|{\bf u}^{sh}({\bf x})\right\|_{1,\Omega_{n+1}^{s}}$. *

Lemma 7.16**.**

There exists a positive real number $c_{2}$ such that

[TABLE]

Proof 7.17**.**

Using Eq. 51 we have:

[TABLE]

*In the above, $\left\|{\bf F}\right\|_{0,\Omega_{\bf X}^{s}}=\left\|\nabla_{\bf X}{\bf x}_{n+1}\right\|_{0,\Omega_{\bf X}^{s}}\leq C_{F}$ thanks to the Lipschitz inversible assumption of ${\bf x}={\bf x}\left({\bf s},t\right)$ (${\bf s}\in\Omega_{\bf X}^{s}$). *

Lemma 7.18**.**

There exists a positive real number $c_{3}$ such that

[TABLE]

Proof 7.19**.**

Using Eq. 51 we have:

[TABLE]

From the definition of $a^{s}({\bf u},{\bf v})$ Eq. 33 and $\ell^{s}({\bf v})$ Eq. 36, using the above three lemmas, we have the following two corollaries.

Corollary 7.20**.**

$a^{s}({\bf u}^{sh},{\bf v}^{sh})$* is bounded in $\mathbb{V}^{h}$, i.e., there exists a positive real number $C_{s}^{h}$, such that for $\forall{\bf u}^{h},{\bf v}^{h}\in{\mathbb{V}^{h}}$,*

[TABLE]

Corollary 7.21**.**

$\ell^{s}({\bf v}^{sh})$* is bounded in $\mathbb{V}^{h}$, i.e., there exists a positive real number $C_{\ell}^{h}$, such that for $\forall{\bf v}^{h}\in{\mathbb{V}^{h}}$,*

[TABLE]

Lemma 7.22**.**

$a^{h}({\bf u}^{h},{\bf v}^{h})$* is coercive, i.e., there exists a positive real number $c_{a}^{h}$, such that for $\forall{\bf u}^{h}\in\mathbb{V}^{h}$,*

[TABLE]

Proof 7.23**.**

*This proof follows the same procedure as the proof for Lemma 7.3 by changing $a({\bf u},{\bf u})$ to $a^{h}({\bf u}^{h},{\bf u}^{h})$. *

Lemma 7.24**.**

$a^{h}({\bf u}^{h},{\bf v}^{h})$* is bounded, i.e., there exists a positive real number $C_{a}^{h}$, such that for $\forall{\bf u}^{h},{\bf v}^{h}\in{\mathbb{V}^{h}}$,*

[TABLE]

Proof 7.25**.**

*By the Cauchy-Schwarz inequality, $a^{f}({\bf u}^{h},{\bf v}^{h})$ is a bounded bilinear functional. Using the definition Eq. 47, $a^{h}({\bf u}^{h},{\bf v}^{h})$ is bounded due to Eq. 54. *

Lemma 7.26**.**

$\ell^{h}({\bf v}^{h})$* is bounded, i.e., there exists a positive real number $C_{\ell}^{h}$, such that for $\forall{\bf v}^{h}\in{\mathbb{V}^{h}}$,*

[TABLE]

Proof 7.27**.**

*By the Cauchy-Schwarz inequality, $\ell^{f}({\bf v}^{h})$ is a bounded linear functional. According to the definition Eq. 48, $\ell^{h}({\bf v}^{h})$ is bounded due to Eq. 55. *

We choose the $P_{2}/P_{1}$ or $P_{2}/(P_{1}+P_{0})$ element which satisfies the following discrete inf-sup condition [4][7, Section 12.6]:

Proposition 7.28** (Discrete Inf-Sup Condition).**

There exist $c_{b}^{h}>0$ such that

[TABLE]

Using the above three lemmas (Lemma 7.22 - 7.26) and the discrete inf-sup condition, we further have the following optimal error estimate result [7, Corollary 12.5.18].

Theorem 7.29**.**

Let $({\bf u},p)$ and $({\bf u}^{h},p^{h})$ be the solution pairs of 7.2 and 7.13 respectively, then there is a constant $c$ depending on $C_{a}^{h}$, $c_{a}^{h}$ and $c_{b}^{h}$ such that

[TABLE]

8 Implementation

We shall use an incompressible neo-Hookean model as described in section 8.1. The treatment of convection is presented in section 8.2, and the preconditioned iterative solver is introduced in section 8.3.

8.1 Neo-Hookean hyperelastic solid

For an incompressible neo-Hookean hyperelastic material, the energy function is described as:

[TABLE]

Notice that $tr_{{\bf F}^{\rm T}{\bf F}}=tr_{{F}_{ki}{F}_{kj}}={F}_{ki}^{2}$, therefore $\frac{\partial tr_{{\bf F}^{\rm T}{\bf F}}}{\partial{\bf F}}=\frac{\partial}{\partial{F}_{mn}}\left({F}_{ki}^{2}\right)=2{\bf F}$ and ${\bf P}=\mu^{s}{\bf F}$, hence according to (5),

[TABLE]

where ${\bf B}={\bf F}{\bf F}^{\rm T}$ is the left Cauchy-Green deformation tensor.

Remark 8.1**.**

It is convenient to choose the reference configuration $({\bf B}={\bf I})$ the same as the initial configuration if the solid is initially stress free. In this case, $J\equiv 1$, and if both the fluid and solid are initially stationary, it can be seen from (4) and (62) that the solid is subjected to a hydrostatic stress $\mu^{s}-p^{s}$, while the fluid is subject to a hydrostatic stress $-p^{f}$ due to (3). From the boundary condition (7), we can observe that $\mu^{s}-p^{s}=-p^{f}$ or $p^{s}-p^{f}=\mu^{s}$. In order to avoid this unnecessary jump of pressure across the interface, we shall replace (62) as:

[TABLE]

*This actually does not influence the solid constitutive equation (4), because $J^{-1}\mu^{s}{\bf I}$ can be integrated into the pressure term. However, this is important for a numerical scheme based on unfitted meshes as used in this article, where difficult to accurately capture the jump of pressure across the fluid-solid interface. This also does not change analysis of the energy-conservation and well-posedness for the proposed scheme. *

Expressing ${\bf F}$ by velocity as shown in Eq. 31, then the weak form Eq. 17 becomes:

[TABLE]

Remark 8.2**.**

*Like noted in Remark 6.8, term $\frac{\mu^{s}\Delta t}{2}\int_{\Omega_{\bf X}^{s}}\nabla_{\bf X}{\bf u}_{*}:\nabla_{\bf X}{\bf v}d{\bf X}$ in above equation will introduce an error $\Delta tO(h^{m})$ ($m>0$) after discretization in space, which then reduce the order to $O(\Delta t)$ for a fixed mesh. *

8.2 Treatment of convection

Without considering the non-linear convection term, equation Eq. 64 is almost a linear equation except the moving domain $\Omega_{n+1}^{s}$. So we have to iteratively construct $\Omega_{n+1}^{s}$ and take derivative on it. For the low Reynolds number ($Re<50$) that we consider in this article, the convection term can also be arranged on the right-hand side of the equation and included in one iteration loop, i.e., a fixed point iteration is adopted to solve the nonlinear system at each time step. For other methods to treat convection, readers may refer to [20, 28]. Notice that no artificial diffusion terms are added in this case, and we can measure the system energy exactly using the energy functions defined from Eq. 24 to Eq. 28.

Remark 8.3**.**

$J_{n+1}^{-1}=1/det_{{\bf F}_{n+1}}$* makes the term $\mu^{s}\int_{\Omega_{n+1}^{s}}J_{n+1}^{-1}\nabla\cdot{\bf v}d{\bf x}$ non-linear in equation Eq. 64. This term is also iteratively solved together with the convection term in one loop. *

8.3 Iterative linear algebra solver

7.13 leads to the following linear equation system [24]:

[TABLE]

where

[TABLE]

and

[TABLE]

In the above, matrix ${\bf D}$ is the isoparametric interpolation matrix derived from equation (46) which can be expressed as

[TABLE]

All the other matrices and vectors arise from standard FEM discretization: ${\bf M}$ and ${\bf M}^{s}$ are mass matrices from discretization of the first term of $a^{f}({\bf u}^{h},{\bf v}^{h})$ and the first term of $a^{s}({\bf u}^{sh},{\bf v}^{sh})$ respectively, and similarly the stiffness matrices ${\bf K}$ and ${\bf K}^{s}$ are from the last term of $a^{f}({\bf u}^{h},{\bf v}^{h})$ and the last two terms of $a^{s}({\bf u}^{sh},{\bf v}^{sh})$ respectively. ${\bf B}$ is from discretization of the linear functional $b\left({\bf v}^{h},p^{h}\right)$. The force vectors ${\bf f}$ and ${\bf f}^{s}$ come from discretization of $\ell^{f}({\bf u}^{h})$ and $\ell^{s}({\bf u}^{sh})$ respectively.

We use the block matrix $\begin{bmatrix}{{\bf M/}\Delta t+{\bf K}}&{\bf 0}\\ {\bf 0}^{\rm T}&{{\bf M/}\Delta t}\end{bmatrix}$ as a preconditioner, computed by an incomplete Cholesky decomposition, and MinRes algorithm [9] to solve equation Eq. 65. A stable convergence performance can be observed although this is not the topic of this article.

9 Numerical experiments

In this section, we focus on the validation of energy conservation of the proposed numerical method in two and three dimensions. For more two-dimensional numerical examples and validation see [24]. The $P_{2}/(P_{1}+P_{0})$ elements will be used, i.e., the standard $Taylor$-$Hood$ elements $P_{2}P_{1}$ is enriched by a constant $P_{0}$ for approximation of the pressure. This element has the property of local mass conservation and the constant $P_{0}$ may better capture the element-based jump of pressure [1, 4]. We shall demonstrate the improvement of mass conservation and energy conservation by using the $P_{2}/(P_{1}+P_{0})$ elements compared to the $P_{2}P_{1}$ elements. We shall also compare the Crank-Nicolson scheme and the backward Euler scheme (see Appendix A) in this section.

9.1 Oscillating disc driven by an initial kinetic energy

In this test, we consider an enclosed flow (${\bf n}\cdot{\bf u}=0$) in $\Omega=[0,1]\times[0,1]$ with a periodic boundary condition. A solid disc is initially located in the middle of the square $\Omega$ and has a radius of $0.2$. The initial velocity of the fluid and solid are prescribed by the following stream function

[TABLE]

where $\Psi_{0}=5.0\times 10^{-2}$ and $a=b=2\pi$. Two parameter sets (Table 1) are used to test the energy and mass conservation based on four different uniform meshes on $\Omega$: (1) $50\times 50$ ($h=0.02$), (2) $66\times 66$ ($h=0.015$), (3) $100\times 100$ ($h=0.01$) and (4) $133\times 133$ ($h=0.0075$). The solid mesh is constructed to have a similar node density to the fluid mesh. In order to visualise the flow a snapshot of the velocity and deformation fields for the first set of parameters is presented in Fig. 2. For a comparison between the flows energy plots are presented in Fig. 3, from which it can be seen that the solid corresponding to the second parameter set is harder and has a larger frequency than the first one. In addition, $E_{d}(\Omega_{t}^{s})<0$ and both the viscosity and density jump across the interface.

We commence by comparing the Crank-Nicolson (CN) scheme and backward Euler (BE) scheme, and comparing $P_{2}/P_{1}$ elements and $P_{2}/(P_{1}+P_{0})$ elements. The evolution of mass variation and energy error are demonstrated in Fig. 4 and Fig. 5 respectively, from which it is clearly apparent that the Crank-Nicolson scheme has an advantage for both the mass and energy conservation. It can be seen from Fig. 4 that the enrichment of the pressure field by a constant $P_{0}$ dramatically improves the mass conservation, although the effect for energy conservation appears to be very slightly negative from Fig. 5.

We then use element $P_{2}/(P_{1}+P_{0})$ and investigate time convergence of the proposed method. It can be observed, from Fig. 6 (a) and Fig. 7 (a), that both Crank-Nicolson and backward Euler scheme have a first order time convergence (see Remark 6.8 and Remark 8.2 for the reason that Crank-Nicolson scheme is also first order), however the former is more accurate than the latter using the same time step. It can also be observed, from Fig. 6 and Fig. 7, that Crank-Nicolson scheme introduces more oscillation than backward Euler scheme in order to gain this accuracy, however, the oscillation can be reduced by mesh refinement (see Fig. 8 and Fig. 9).

Finally, the mesh convergence is clearly demonstrated in Fig. 8 and Fig. 9 via mass and energy evolution respectively.

9.2 Oscillating disc driven by an initial potential energy

In the previous example, the disc oscillates because a kinetic energy (the 2nd and 3rd terms in Eq. 30) is prescribed for the FSI system at the beginning. In this test, we shall stretch the disc and create a potential energy in the solid (the last term in Eq. 30), then release it causing the disc to oscillate due to this potential solid energy. The computational domain is a square $\Omega=[0,1]\times[0,1]$. One quarter of a solid disc is located in the left-bottom corner of the square, and initially stretched as an ellipse as shown in Fig. 10. Notice the equation of an ellipse $\frac{x^{2}}{a^{2}}+\frac{y^{2}}{b^{2}}=1$ and its area $\pi ab$, hence we ensure that this stretch does not change mass of the solid.

We choose $\rho^{f}=1$, $\rho^{s}=2$, $\mu^{s}=2$ and $\nu^{f}=\nu^{s}=0.01$. The fluid adopts the same meshes defined in section 9.1, and the solid has similar node density as the fluid. A snapshot of pressure on the fluid mesh and corresponding solid deformation with its velocity norm are displayed in Fig. 11. Time and mesh convergence are shown in Fig. 13 and Fig. 13 respectively.

9.3 Oscillating ball driven by an initial kinetic energy

In this section, we consider a 3D oscillating ball, which is an extension of the example in section 9.1. The ball is initially located at the center of $\Omega=[0,1]\times[0,1]\times[0,0.6]$ with a radius of $0.2$. Using the property of symmetry this computation is carried out on $1/8$ of domain $\Omega$: $[0,0.5]\times[0,0.5]\times[0,0.3]$. The initial velocities of $x$ and $y$ components are the same as that used in section 9.1 and the $z$ component is set to be 0 at the beginning. We adopt the Parameter set 1 and Mesh (1) defined in section 9.1 (with the same mesh size in z direction). A snapshot of the $1/8$ solid ball and the corresponding fluid velocity norm are presented in Fig. 14, and the result of time convergence is presented in Fig. 16 and Fig. 16, from which a first order convergence can be observed.

10 Conclusions

In this article, the energy conservation of a new one-field fictitious domain method for fluid-structure interactions is proved after discretization in time. Furthermore, for a special case of the Stokes fluid equation and the neo-Hookean solid model, the well-posedness of the proposed scheme is analyzed and demonstrated.

We then present a selection of numerical tests that demonstrate the theoretical energy estimate in both $2$ and $3$ dimensions. Moreover, we show that the Crank-Nicolson scheme is more accurate than the backward Euler scheme in terms of mass and energy conservation, although both exhibit first order convergence in time (see Remark 6.8 and Remark 8.2 for the reason that Crank-Nicolson scheme is also first order).

Appendix A Energy estimate of the backward Euler scheme

A.1 Weak form

Using the backward Euler method, the discretized weak form corresponding to 1 is:

Problem A.1**.**

For each time step, find ${\bf u}_{n+1}\in H_{0}^{1}(\Omega)^{d}$ and $p_{n+1}\in L_{0}^{2}(\Omega)$, such that

[TABLE]

$\forall{\bf v}\in H_{0}^{1}\left(\Omega\right)^{d}$* and*

[TABLE]

$\forall q\in L^{2}(\Omega)$*. *

Remark A.2**.**

*For backward Euler scheme, $\Omega_{n+1}^{s}$ is updated from $\Omega_{n}^{s}$ by $\Omega_{n+1}^{s}=\left\{{\bf x}_{n+1}:{\bf x}_{n+1}={\bf x}_{n}+{\bf u}_{n+1}\Delta t\right\}$, for all ${\bf x}_{n}\in\Omega_{n}^{s}$. *

A.2 Energy conservation

Lemma A.3**.**

Assume the solid energy function $\Psi(\cdot)\in C^{1}$ over the set of second order tensors, then

[TABLE]

Proof A.4**.**

Let $w(\xi)=\Psi\left({\bf F}_{n}+\xi\nabla_{\bf X}{\bf u}_{n+1}\right)$, and notice that

[TABLE]

then,

[TABLE]

Using the chain rule, Eq. 70 holds thanks to

[TABLE]

Similarly to Lemma 4.3, we have:

Lemma A.5**.**

If $\left({\bf u}_{n+1},p_{n+1}\right)$ is the solution pair of A.1, then

[TABLE]

Proposition A.6** (Local Energy Conservation).**

Let $\left({\bf u}_{n+1},p_{n+1}\right)$ be the solution pair of A.1, then

[TABLE]

Proof A.7**.**

Let ${\bf v}={\bf u}_{n+1}$ in Eq. 68, $q=p_{n+1}$ in Eq. 69 and use Lemma A.5, we have:

[TABLE]

Using Cauchy-Schwarz inequality and the fact $ab\leq\frac{a^{2}+b^{2}}{2}$, we have:

[TABLE]

*where $\omega=\Omega$ or $\Omega_{n+1}^{s}$. Substituting the above equation into Eq. 73, we get Proposition A.6 due to Lemma A.3. *

Corollary A.8** (Total Energy Conservation).**

Let $N+1=t/\Delta t$, where $t$ is the computational time, and denote the error of total energy as:

[TABLE]

then,

[TABLE]

Proof A.9**.**

*First let $O(\Delta t^{2})=\sum_{k\geq 2}C_{n+1}^{k}\Delta t^{k}$ in Eq. 72, where $C_{n+1}^{k}$ is a constant dependent on the specific time step. Then add equation Eq. 72 from $n=0$ to $n=N$, we have $Err\leq\sum_{k\geq 2}\Delta t^{k}\sum_{n=0}^{N}C_{n+1}^{k}$. Eq. 74 holds due to $\sum_{n=0}^{N}C_{n+1}^{k}\leq C_{max}^{k}t/\Delta t$, where $C_{max}^{k}=max\{C_{n+1}^{k},n=0,\cdots,N\}$. *