A theoretical and experimental investigation of a family of immersed finite element methods

TL;DR

This paper provides a comprehensive theoretical and experimental analysis of immersed finite element methods, comparing explicit, implicit, and one-field fictitious domain approaches, highlighting their relationships, robustness, and efficiency in simulating fluid-structure interactions.

Contribution

It demonstrates that the one-field FDM is more robust and efficient than explicit IFEM, and closely approximates implicit IFEM results without iterative steps.

Findings

One-field FDM is more robust than explicit IFEM.

One-field FDM can simulate wider parameter ranges with larger time steps.

Results of one-field FDM closely match implicit IFEM outcomes.

Abstract

In this article we consider the widely used immersed finite element method (IFEM), in both explicit and implicit form, and its relationship to our more recent one-field fictitious domain method (FDM). We review and extend the formulation of these methods, based upon an operator splitting scheme, in order to demonstrate that both the explicit IFEM and the one-field FDM can be regarded as particular linearizations of the fully implicit IFEM. However, the one-field FDM can be shown to be more robust than the explicit IFEM and can simulate a wider range of solid parameters with a relatively large time step. In addition, it can produce results almost identical to the implicit IFEM but without iteration inside each time step. We study the effect on these methods of variations in viscosity and density of fluid and solid materials. The advantages of the one-field FDM within the IFEM framework are illustrated through a selection of parameter sets for two benchmark cases.

Tables

Table 1. Table 1: Parameter sets for lid-driven cavity flow with a solid disc (the first parameter set is used in [ 31 , 8 ] ).

Parameter sets	$Re$	${\tilde{c}}_1$	${\rho }^r$	${\mu }^r$
Parameter 1	$100$	$0.1$	$1$	$1$
Parameter 2	$100$	$1$	$1$	$1$
Parameter 3	$100$	$1$	$2$	$1$
Parameter 4	$100$	$1$	$0.5$	$1$
Parameter 5	$100$	$1$	$1$	$1.5$
Parameter 6	$100$	$1$	$1$	$2$
Parameter 7	$500$	$0.5$	$2$	$2$

Table 2. Table 2: Parameter sets for the oscillating leaflet in a channel.

Parameter sets	$Re$	${\tilde{c}}_1$	${\rho }^r$	${\mu }^r$
Parameter 1	$100$	$1000$	$1$	$1$
Parameter 2	$100$	$1000$	$1.2$	$1$
Parameter 3	$300$	$10000$	$2$	$2$

Equations

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

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

$$\mathbf{\sigma }={\mathbf{\sigma }}^f={\mathbf{\tau }}^f-p^f\mathbf{I},$$

$${\mathbf{\tau }}^f={\mu }^f\mathrm{D}{\mathbf{u}}^f$$

$$\mathbf{\sigma }={\mathbf{\sigma }}^s={\mathbf{\tau }}^s-p^s\mathbf{I},$$

$${\mathbf{\tau }}^s=c_1\left(\mathbf{F}{\mathbf{F}}^T-\mathbf{I}\right)+{\mu }^s\mathrm{D}{\mathbf{u}}^s$$

$$\mathbf{F}=\frac{\partial \mathbf{x}}{\partial \mathbf{X}}=\frac{\partial \mathbf{x}}{\partial {\mathbf{x}}_0}\frac{\partial {\mathbf{x}}_0}{\partial \mathbf{X}}={\nabla }_0\mathbf{x}{\nabla }_{\mathbf{X}}{\mathbf{x}}_0$$

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

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

$${\mathbf{u}}^f=\bar{\mathbf{u}}on{\Gamma }_D$$

$${\mathbf{\sigma }}^f{\mathbf{n}}^s=\bar{\mathbf{h}}on{\Gamma }_N,$$

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

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

$$\begin{array}{rl}{\int }_{\Omega }Eq.(\text{ref}momentu{m}_{e}quation)\left(\mathbf{\sigma }\right)\cdot \mathbf{v}& \equiv {\int }_{{\Omega }_t^f}Eq.(\text{ref}momentu{m}_{e}quation)\left({\mathbf{\sigma }}^f\right)\cdot \mathbf{v}+{\int }_{{\Omega }_t^s}Eq.(\text{ref}momentu{m}_{e}quation)\left({\mathbf{\sigma }}^s\right)\cdot \mathbf{v}\\ & \equiv {\int }_{\Omega }Eq.(\text{ref}momentu{m}_{e}quation)\left({\mathbf{\sigma }}^f\right)\cdot \mathbf{v}+{\int }_{{\Omega }_t^s}\left(Eq.(\text{ref}momentu{m}_{e}quation)\left({\mathbf{\sigma }}^s\right)-Eq.(\text{ref}momentu{m}_{e}quation)\left({\mathbf{\sigma }}^f\right)\right)\cdot \mathbf{v}.\end{array}$$

$$\begin{array}{rl}& {\rho }^f{\int }_{\Omega }\frac{d\mathbf{u}}{dt}\cdot \mathbf{v}+{\int }_{\Omega }{\mathbf{\tau }}^f:\nabla \mathbf{v}-{\int }_{\Omega }p\nabla \cdot \mathbf{v}\\ +& \left({\rho }^s-{\rho }^f\right){\int }_{{\Omega }_t^s}\frac{d\mathbf{u}}{dt}\cdot \mathbf{v}+{\int }_{{\Omega }_t^s}\left({\mathbf{\tau }}^s-{\mathbf{\tau }}^f\right):\nabla \mathbf{v}\\ =& {\int }_{\Omega }{\rho }^f\mathbf{g}\cdot \mathbf{v}+{\int }_{{\Omega }_t^s}\left({\rho }^s-{\rho }^f\right)\mathbf{g}\cdot \mathbf{v}+{\int }_{{\Gamma }_N}\bar{\mathbf{h}}\cdot \mathbf{v}.\end{array}$$

$$-{\int }_{{\Omega }_t^f}Eq.(\text{ref}continuit{y}_{e}quation)q-{\int }_{{\Omega }_t^s}Eq.(\text{ref}continuit{y}_{e}quation)q\equiv -{\int }_{\Omega }Eq.(\text{ref}continuit{y}_{e}quation)q,$$

$$\begin{array}{rl}& {\rho }^f{\int }_{\Omega }\frac{\partial \mathbf{u}}{\partial t}\cdot \mathbf{v}+{\rho }^f{\int }_{\Omega }\left(\mathbf{u}\cdot \nabla \right)\mathbf{u}\cdot \mathbf{v}+\frac{{\mu }^f}{2}{\int }_{\Omega }\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{v}-{\int }_{\Omega }p\nabla \cdot \mathbf{v}\\ +& {\rho }^{\delta }{\int }_{{\Omega }_t^s}\frac{\mathfrak{d}\mathbf{u}}{\mathfrak{d}t}\cdot \mathbf{v}+\frac{{\mu }^{\delta }}{2}{\int }_{{\Omega }_t^s}\mathrm{D}\mathbf{u}:\mathrm{D}\mathbf{v}+c_1{\int }_{{\Omega }_t^s}\left(\mathbf{F}{\mathbf{F}}^T-\mathbf{I}\right):\nabla \mathbf{v}\\ =& {\int }_{\Omega }{\rho }^f\mathbf{g}\cdot \mathbf{v}+{\int }_{{\Omega }_t^s}{\rho }^{\delta }\mathbf{g}\cdot \mathbf{v}+{\int }_{{\Gamma }_N}\bar{\mathbf{h}}\cdot \mathbf{v},\end{array}$$

$$-{\int }_{\Omega }q\nabla \cdot \mathbf{u}=0.$$

$$\begin{array}{rl}& {\int }_{\Omega }\frac{\partial \tilde{\mathbf{u}}}{\partial \tilde{t}}\cdot \tilde{\mathbf{v}}+{\int }_{\Omega }\left(\tilde{\mathbf{u}}\cdot \nabla \right)\tilde{\mathbf{u}}\cdot \tilde{\mathbf{v}}+\frac{1}{2Re}{\int }_{\Omega }\mathrm{D}\tilde{\mathbf{u}}:\mathrm{D}\tilde{\mathbf{v}}-{\int }_{\Omega }\tilde{p}\nabla \cdot \tilde{\mathbf{v}}\\ +& \left({\rho }^r-1\right){\int }_{{\Omega }_t^s}\frac{\mathfrak{d}\tilde{\mathbf{u}}}{\mathfrak{d}\tilde{t}}\cdot \tilde{\mathbf{v}}+\frac{{\mu }^r-1}{2Re}{\int }_{{\Omega }_t^s}\mathrm{D}\tilde{\mathbf{u}}:\mathrm{D}\tilde{\mathbf{v}}+{\tilde{c}}_1{\int }_{{\Omega }_t^s}\left(\mathbf{F}{\mathbf{F}}^T-\mathbf{I}\right):\nabla \mathbf{v}\\ =& Fr{\int }_{\Omega }\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \tilde{\mathbf{v}}+Fr\left({\rho }^r-1\right){\int }_{{\Omega }_t^s}\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \tilde{\mathbf{v}}+{\int }_{{\Gamma }_N}\tilde{\mathbf{h}}\cdot \tilde{\mathbf{v}},\end{array}$$

$$Re={\rho }^f\tilde{U}\tilde{L}/{\mu }^f.$$

$${\rho }^r={\rho }^s/{\rho }^f.$$

$${\mu }^r={\mu }^s/{\mu }^f.$$

$${\tilde{c}}_1=c_1/{\rho }^f{\tilde{U}}^2.$$

$$Fr=|\mathbf{g}|\tilde{L}/{\tilde{U}}^2.$$

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

$$\begin{array}{rl}& {\int }_{\Omega }\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}\cdot \mathbf{v}+{\int }_{\Omega }\left({\mathbf{u}}_{n+1}\cdot \nabla \right){\mathbf{u}}_{n+1}\cdot \mathbf{v}\\ +& \frac{1}{2Re}{\int }_{\Omega }\mathrm{D}{\mathbf{u}}_{n+1}:\mathrm{D}\mathbf{v}-{\int }_{\Omega }{p}_{n+1}\nabla \cdot \mathbf{v}\\ +& \left({\rho }^r-1\right){\int }_{{\Omega }_{n+1}^s}\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}\cdot \mathbf{v}+\frac{{\mu }^r-1}{2Re}{\int }_{{\Omega }_{n+1}^s}\mathrm{D}{\mathbf{u}}_{n+1}:\mathrm{D}\mathbf{v}\\ +& {\tilde{c}}_1{\int }_{{\Omega }_{n+1}^s}\left(\mathbf{F}{\mathbf{F}}^T-\mathbf{I}\right):\nabla \mathbf{v}\\ =& Fr{\int }_{\Omega }\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \mathbf{v}+Fr\left({\rho }^r-1\right){\int }_{{\Omega }_{n+1}^s}\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \mathbf{v}+{\int }_{{\Gamma }_N}\tilde{\mathbf{h}}\cdot \mathbf{v},\end{array}$$

$$-{\int }_{\Omega }q\nabla \cdot {\mathbf{u}}_{n+1}=0.$$

$${\Omega }_{n+1}^s=\left\{\mathbf{x}:\mathbf{x}={\mathbf{x}}_n+\Delta t{\mathbf{u}}_{n+1},{\mathbf{x}}_n\in {\Omega }_n^s\right\}.$$

$$\begin{array}{rl}& {\int }_{\Omega }\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}\cdot \mathbf{v}+{\int }_{\Omega }\left({\mathbf{u}}_{n+1}\cdot \nabla \right){\mathbf{u}}_{n+1}\cdot \mathbf{v}\\ +& \frac{1}{2Re}{\int }_{\Omega }\mathrm{D}{\mathbf{u}}_{n+1}:\mathrm{D}\mathbf{v}-{\int }_{\Omega }{p}_{n+1}\nabla \cdot \mathbf{v}\\ =& \left(1-{\rho }^r\right){\int }_{{\Omega }_n^s}\frac{{\mathbf{u}}_n-{\mathbf{u}}_{n-1}}{\Delta t}\cdot \mathbf{v}-\boxed{\frac{{\mu }^r-1}{2Re}{\int }_{{\Omega }_n^s}\mathrm{D}{\mathbf{u}}_n:\mathrm{D}\mathbf{v}}\\ -& {\tilde{c}}_1{\int }_{{\Omega }_n^s}\left({\mathbf{F}}_n{\mathbf{F}}_n^T-\mathbf{I}\right):\nabla \mathbf{v}\\ +& Fr{\int }_{\Omega }\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \mathbf{v}+Fr\left({\rho }^r-1\right){\int }_{{\Omega }_n^s}\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \mathbf{v}+{\int }_{{\Gamma }_N}\tilde{\mathbf{h}}\cdot \mathbf{v}.\end{array}$$

$$\begin{array}{rl}& \left(1+\left({\rho }^r-1\right)I\left(\mathbf{x}\right)\right){\int }_{\Omega }\frac{{\mathbf{u}}_{n+1}-{\mathbf{u}}_n}{\Delta t}\cdot \mathbf{v}+{\int }_{\Omega }\left({\mathbf{u}}_{n+1}\cdot \nabla \right){\mathbf{u}}_{n+1}\cdot \mathbf{v}\\ +& \frac{1}{2Re}{\int }_{\Omega }\mathrm{D}{\mathbf{u}}_{n+1}:\mathrm{D}\mathbf{v}-{\int }_{\Omega }{p}_{n+1}\nabla \cdot \mathbf{v}\\ =& \frac{1-{\mu }^r}{2Re}{\int }_{{\Omega }_n^s}\mathrm{D}{\mathbf{u}}_n:\mathrm{D}\mathbf{v}-{\tilde{c}}_1{\int }_{{\Omega }_n^s}\left({\mathbf{F}}_n{\mathbf{F}}_n^T-\mathbf{I}\right):\nabla \mathbf{v}\\ +& Fr{\int }_{\Omega }\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \mathbf{v}+Fr\left({\rho }^r-1\right){\int }_{{\Omega }_n^s}\frac{\mathbf{g}}{|\mathbf{g}|}\cdot \mathbf{v}+{\int }_{{\Gamma }_N}\tilde{\mathbf{h}}\cdot \mathbf{v}.\end{array}$$

Full text

††journal: Journal of Fluids and Structures

1 Introduction

Three major questions arise when considering a finite element method for the problem 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 to these questions lead to different types of numerical method. For example, [1, 2] solve for fluid velocity and solid displacement sequentially (partitioned/segregated) using an Arbitrary Lagrangian-Eulerian (ALE) fitted mesh, whereas [3, 4, 5] 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) [6, 7, 8, 9, 10, 11] and the Fictitious Domain Method (FDM) [12, 13, 14, 15, 16, 17] use two meshes to represent the fluid and solid separately. Although IFEM could be monolithic [18], the classical IFEM only solves for velocity, while the solid information is arranged on the right-hand side of the fluid equation as a prescribed force term. Conversely, although the FDM may be partitioned [17], 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, [19] introduces a 1D model using a one-field FDM formulation based on two meshes, and [20, 21] introduces an energy stable monolithic method (in 2D) based on one Eulerian mesh and discrete remeshing.

In a previous study [22], we present a one-field monolithic fictitious domain method (subsequently referred to as the one-field FDM) which has the following main features: (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. The primary purpose of this paper is to highlight the relationship between the one-field FDM and the IFEM family of methods: demonstrating that the former is in fact a new variant of the latter which possesses a number of practical advantages. Before describing these in detail however we briefly provide further context for the one-field FDM based on comparing its features with those of existing monolithic schemes.

FDM/DLM methods [12, 13, 14, 15, 16, 17] solve the solid equation, but for a displacement field, and couple this displacement with the velocity of the fictitious fluid via a Lagrange multiplier. This leads to a large discrete linear algebra system. The one-field FDM solves a smaller discrete system since it rewrites the solid equation in terms of a velocity variable and couples the fictitious fluid through a finite element interpolation. Monolithic Eulerian methods [20, 21] also express the solid equation in terms of velocity, and the fluid and solid are coupled naturally on an interface-fitted mesh. The one-field FDM also uses two meshes to represent the fluid and solid respectively. Consequently, before discretization in space, these two methods have many similarities, the advantage of the one-field FDM being that interface fitting is not required.

As discussed above, the primary purpose of this paper is to demonstrate that the proposed one-field FDM is a particular linearization of the fully implicit IFEM. It is more robust than the classical explicit IFEM, and presents almost the same accuracy as the implicit IFEM. We will show:

• (1)

The one-field FDM is based upon the implicit expression of the FSI force. This FSI force term is linearized using the velocity at the current configuration (instead of displacement at the reference configuration), which is an approximation of the fully implicit IFEM [23, 24, 25] but without requiring a nonlinear iteration at each time step.

• (2)

In the simple case of equal density and viscosity for both fluid and solid, the only difference between the one-field FDM and the explicit IFEM (explicitly expressing the FSI force term) [7, 8, 9, 10, 11] is that there are some terms of order $O(\Delta t)$ and $O(\Delta t^{2})$ retained in the former. However these terms have a helpful stabilizing effect, which can allow the one-field FDM to use a time step that is significantly larger than the explicit IFEM, with almost the same accuracy.

• (3)

The one-field FDM can naturally deal with the case of different densities and different viscosities in the fluid and solid.

The paper is organized as follows. The control equations and a general finite element weak formulation are introduced in Section 2 and 3 respectively, followed by a dimensionless weak formulation in Section 4 and time discretization in Section 5. Different types of IFEM approaches are then discussed in Section 6, and their relationship with the weak form of Section 4 is highlighted. An operator splitting scheme is introduced in Section 7 followed by a comparison between the IFEM and one-field FDM approaches in Section 8 and 9. Numerical examples are given in Section 10, and conclusions are presented in Section 11.

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 Figure 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}$ the initial coordinates of the solid.

Let $\rho,\mu,{\bf u},{\bm{\sigma}}$ and ${\bf g}$ denote the density, viscosity, velocity, stress tensor and acceleration due to gravity respectively. We assume both an incompressible fluid in $\Omega_{t}^{f}$ and incompressible solid in $\Omega_{t}^{s}$. The conservation of momentum and conservation of mass take the same form in the fluid and solid (just differing in the specific expressions of ${\bm{\sigma}}$):

Momentum equation:

[TABLE]

and continuity equation:

[TABLE]

An incompressible Newtonian constitutive equation in $\Omega_{t}^{f}$ can be expressed as:

[TABLE]

where ${\rm D}{\bf u}=\nabla{\bf u}+\nabla^{\scriptsize T}{\bf u}$, and

[TABLE]

is the deviatoric part of stress ${\bm{\sigma}}^{f}$.

In $\Omega_{t}^{s}$ we use an incompressible neo-Hookean solid with viscosity $\mu^{s}$ [13, 17, 12] (see appendices for a brief discussion of other solid material models). The constitutive equation may be expressed as:

[TABLE]

with

[TABLE]

being the deviatoric part of stress ${\bm{\sigma}}^{s}$, where

[TABLE]

is the deformation tensor of the solid, and $c_{1}$ a solid material parameter.

Finally the system is completed with the following boundary and initial conditions.

Interface continuity conditions:

[TABLE]

and

[TABLE]

Dirichlet and Neumann boundary conditions:

[TABLE]

and

[TABLE]

where $\Gamma=\Gamma_{D}\cup\Gamma_{N}$. Initial conditions:

[TABLE]

and

[TABLE]

3 Weak formulation

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 Dirichlet boundary of $\omega$.

Let ${p}=\left\{\begin{matrix}{{p}^{f}\quad in\quad\Omega_{t}^{f}}\\ {{p}^{s}\quad in\quad\Omega_{t}^{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 (3) and (5) and interface condition (9), give:

[TABLE]

Note that the integrals on the interface $\Gamma_{t}$ are cancelled out using boundary condition (9), because they are internal forces for the whole FSI system. Combining with the following symbolic operations for $q\in L^{2}(\Omega)$,

[TABLE]

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

Given ${\bf u}_{0}$ and $\Omega_{0}^{s}$, for each $t>0$ find ${\bf u}(t)\in H^{1}(\Omega)^{d}$, $p(t)\in L^{2}(\Omega)$ and $\Omega_{t}^{s}$, such that $\forall{\bf v}\in H_{0}^{1}(\Omega)^{d}$, $\forall q\in L^{2}(\Omega)$, the following two equations hold:

[TABLE]

(where $\rho^{\delta}=\rho^{s}-\rho^{f}$ and $\mu^{\delta}=\mu^{s}-\mu^{f}$) and

[TABLE]

In equation (15), $\frac{\mathfrak{d}}{\mathfrak{d}t}$ is the time derivative with respect to a frame moving with the solid velocity ${\bf u}^{s}=\left.{\bf u}\right|_{\Omega_{t}^{s}}$.

4 Dimensionless weak formulation

We may rewrite the weak form (15) and (16) in a dimensionless form by introducing the following scales [17]: $\tilde{L}$ for length, $\tilde{U}$ for velocity, $\tilde{L}/\tilde{U}$ for time and $\rho^{f}\tilde{U}^{2}$ for pressure $p$ and parameter $c_{1}$. Dividing by $\rho^{f}\tilde{U}^{3}/\tilde{L}$ on both sides of equation (15), we have:

[TABLE]

where $\tilde{\bf u}={\bf u}/\tilde{U}$, $\tilde{\bf v}={\bf v}/\tilde{U}$, $\tilde{p}=p/\rho^{f}\tilde{U}^{2}$, $\tilde{t}=t\tilde{U}/\tilde{L}$, $\tilde{\bf h}=\bar{\bf h}\tilde{L}/\rho^{f}\tilde{U}^{2}$ and the following parameters.

Reynolds number:

[TABLE]

• 2.

density ratio:

[TABLE]

• 3.

viscosity ratio:

[TABLE]

• 4.

material parameter:

[TABLE]

• 5.

Froude number:

[TABLE]

The dimensionless form of equation (16) is straightforward by dividing by $\rho^{f}\tilde{U}^{3}/\tilde{L}$:

[TABLE]

with $\tilde{q}=q/\rho^{f}\tilde{U}^{2}$. For the sake of notation convenience, we shall still use ${\bf u}$, ${\bf v}$, $p$ and $q$ instead of $\tilde{\bf u}$, $\tilde{\bf v}$, $\tilde{p}$ and $\tilde{q}$ in equation (17) and (23) in the following context.

5 Discretization in time

Using the backward Euler method to discretize in time, equations (17) and (23) may be discretized as follows:

Given ${\bf u}_{n}$, $p_{n}$ and $\Omega_{n}^{s}$, find ${\bf u}_{n+1}\in H^{1}(\Omega)^{d}$, $p_{n+1}\in L^{2}(\Omega)$ and $\Omega_{n+1}^{s}$, such that $\forall{\bf v}\in H_{0}^{1}\left(\Omega\right)^{d}$, $\forall q\in L^{2}(\Omega)$, the following two equations hold:

[TABLE]

and

[TABLE]

Remark 1**.**

$\Omega_{n+1}^{s}$* is updated from $\Omega_{n}^{s}$ by the following formula:*

[TABLE]

6 Implementations of different IFEM approaches

Having introduced the weak formulation and time discretization in the previous sections, we now consider implementation details, and demonstrate that different choices lead to methods that correspond to existing IFEM schemes, as well as the proposed one-field FDM [22]. We can see from (24) that the integrals are carried out over two different domains: the whole domain $\Omega$ and the moving solid domain $\Omega_{n+1}^{s}$. The IFEM methods compute these two types of integrals based on two different meshes, and use projection/distribution functions to transfer information between the two meshes [12, 13, 14, 15, 16, 17, 7, 8, 9, 10, 11]. The one-field FDM also adopts two meshes, and uses the FEM interpolation function (also used in [8]) to transfer information between the two meshes. In the remainder of this section we focus on how different IFEM approaches fit within the framework of this weak formulation.

The classical IFEM is introduced in [11, 10], in which all the solid integrals (in $\Omega_{n+1}^{s}$) are moved to the right-hand side of equation (24) as a force term and evaluated at the previous time step as follows:

[TABLE]

The above formulation differs from [11, 10] only in the following respects:

1. (1)

The boxed term in (27) vanishes in [11, 10] because the fluid stress within the solid domain is neglected (which is equivalent to setting $\mu^{r}=1$). 2. (2)

[11, 10] use the stabilized equal-order finite element method to treat convection and pressure after discretization in space, while we shall use a splitting scheme to deal with convection and a stable mixed-order finite element space for the velocity and pressure. 3. (3)

The above formulation is expressed in a dimensionless form.

However these differences are not the distinguishing features of IFEM, and do not influence any comparisons (in Section 10 we show that our implementation of IFEM produces the same results as in the literature). The distinguishing feature of IFEM is its development from the Immersed Boundary Method (IBM) [6]: the solid information is based on the previous time step and arranged on the right-hand side of the fluid equation as a force term, which is computed on the solid mesh, distributed to the fluid mesh and then added to the fluid equation as an extra term.

The IFEM formulation (27) represents an explicit forcing strategy, which approximates the time derivative in the solid using values from the previous two time steps. Errors may be accumulated as time evolves in this case, and this may be observed in numerical tests (see, for example, Figures 6 and 7). There is a semi-implicit formulation which introduces an indicator function $I({\bf x})$ to smear out the densities across the fluid-solid interface [7]. Based on this indicator function $I({\bf x})$ (see [7] for the definition), the formulation (27) may be modified as follows:

[TABLE]

Furthermore, a fully implicit forcing strategy may also be considered as follows:

[TABLE]

It can be seen that the force term on the right-hand side of above equation is computed using the velocity at the current time step, which then needs to be iteratively constructed. This idea of implicit forcing was first utilized in the immersed boundary method (IBM) [26, 27], and then also adopted within IFEM in [23, 24, 25], where a Newton-Krylov iterative procedure is used to solve this nonlinear system. In our implementation, for simplicity, we use fixed point iteration to demonstrate the implicit IFEM. The fixed point iteration generally converges more slowly than the Newton-Krylov method, however it is not our purpose to compare the efficiency of these two implicit forcing strategies. Instead we shall demonstrate that both the implicit IFEM (iterating at each time step) and the one-field FDM (which needs no iteration inside the time step) produce almost identical results.

Based upon the above discussion, we next introduce an operator splitting scheme in Section 7. In Section 8 we then present an explicit and an implicit forcing strategy for IFEM (corresponding to (27) and (29) respectively), and in Section 9 the one-field FDM formulation is illustrated in detail with the proposed operator splitting scheme.

7 An operator splitting scheme

The operator splitting scheme is introduced here in order to treat the non-linear convection term in the Navier-Stokes equation, and simplify the saddle-point problem, which is widely adopted to solve pure fluid equations [28, 29]. The fluid-structure coupling process can still be designed either explicitly or implicitly inside the diffusion step (as discussed in Section 8 and 9). In this article, we focus on studying the FSI coupling strategies rather than different methods to deal with the convection or saddle-point problem.

• (1)

Convection step:

[TABLE]

• (2)

Diffusion step:

[TABLE]

In the above, $\nabla_{n}(\cdot)=\frac{\partial(\cdot)}{\partial{\bf x}_{n}}$ and ${\rm D}_{n}=\nabla_{n}+\nabla_{n}^{T}$.

• (3)

Pressure step:

[TABLE]

and

[TABLE]

Remark 2**.**

Notice that the variables ${\bf u}_{n+1/3}$ or ${\bf u}_{n+2/3}$ are just intermediate values, not the velocity at time $t=t_{n}+\frac{\Delta t}{3}$ or $t=t_{n}+\frac{2\Delta t}{3}$. The notation ${\bf F}_{n+1/3}$ or ${\bf F}_{n+2/3}$ is interpreted as follows:

[TABLE]

with $t=n+1/3$ or $n+2/3$.

Remark 3**.**

Compared with the 2-step splitting scheme in our original paper [22], this 3-step splitting scheme decouples the FSI problem and Stokes/saddle-point problem into two separate steps. The fluid and solid are coupled in the diffusion step, which may be effectively solved by the preconditioned Conjugate Gradient algorithm. The pressure step becomes a “degenerate” Stokes problem [28], which can also be efficiently solved (readers may refer to [28, Section 34] for more discussion). There are a variety of numerical methods to treat the convection equation (30), such as wave-like methods [28], characteristic based schemes [28, 30, 20], upwind schemes (including the Streamline Upwind Petrov Galerkin (SUPG) method) [28, 30] or the Least-squares method [30]. In our implementations we primarily use this latter approach.

It can be seen that the fluid-structure interaction only occurs in the diffusion step (31) based upon this operator splitting scheme. In order to solve equation (31) the one-field FDM and IFEM use different strategies. In the following two sections we focus on this diffusion step, and present the differences between the one-field FDM, and IFEM strategies.

8 Explicit and implicit forcing

For notational convenience let us define the following force term:

[TABLE]

where $t=n+1/3$ or $n+2/3$ as in Remark 2.

Based upon the splitting scheme adopted here, we use ${\bf u}_{n+1/3}$, obtained from the previous convection step, to evaluate $\mathcal{F}_{t}^{FSI}$, and solve equation (31) as follows.

Explicit forcing:

[TABLE]

As noted above, expression (36) corresponds to a formulation of the classical explicit IFEM. The implicit IFEM may be expressed in a similar form, based upon the splitting scheme, but using the current value ${\bf u}_{n+2/3}$ to construct $\mathcal{F}_{n+2/3}^{FSI}$.

Implicit forcing:

[TABLE]

9 One-field FDM

It can be seen that the solid velocity is hidden in the nonlinear term ${\bf F}{\bf F}^{T}-{\bf I}$ in equation (24) or (31). The difference between the one-field FDM and the explicit IFEM is how to treat this nonlinear term: the former extracts this hidden velocity by linearizing ${\bf F}{\bf F}^{T}-{\bf I}$ in the current configuration, while the latter evaluates this term as a force term on the right-hand side of the equation. In this section, we shall demonstrate how the nonlinear term ${\bf F}{\bf F}^{T}-{\bf I}$ is linearized in the one-field FDM, and expressed in terms of velocity in the current configuration. Also notice that this idea is not limited to the operator splitting. The splitting is just a specific implementation that allows us to express the IFEM and the one-field FDM in a similar form so as to facilitate direct comparison with each other.

Denoting ${\bf F}_{t}{\bf F}_{t}^{T}-{\bf I}$ by ${\bf s}_{t}$, and according to the definition of F (7), ${{\bf s}_{t}}$ may be computed as follows:

[TABLE]

Using the chain rule, this last equation can also be expressed as:

[TABLE]

or

[TABLE]

Then ${{\bf s}_{t}}$ can be expressed based on the previous coordinate ${\bf x}_{n}$ as follows:

[TABLE]

Using ${\bf x}_{t}={\bf x}_{n}+\Delta t{\bf u}_{t}$ (see Remark 2 (34)), the last equation can finally be expressed as:

[TABLE]

There are two nonlinear terms in the last equation, which may be linearized as

[TABLE]

and

[TABLE]

Substituting ${\bf s}_{n+2/3}={{\bf F}_{n+2/3}{\bf F}_{n+2/3}^{T}-{\bf I}}$, using expression (42), (43) and (44), into diffusion step (31), we finally get the one-field FDM formulation as follows:

[TABLE]

In the above, the linear operators ${\rm D}_{n}^{1}$, ${\rm D}_{n}^{2}$ and ${\rm D}_{n}^{3}$ are defined as:

[TABLE]

[TABLE]

and

[TABLE]

Remark 4**.**

When $\rho^{r}=\mu^{r}=1$, comparing equations (36) and (45), we see that the only difference between the one-field FDM and explicit IFEM is the additional boxed terms in equation (45) of $O(\Delta t)$ or $O(\Delta t^{2})$ respectively. We shall demonstrate, with numerical tests, that these terms are not trivial at all: in fact they have a very positive stabilizing effect, such that a significantly larger time step may be adopted.

10 Numerical experiments

Having analyzed the proposed one-field FDM, and illustrated its close relationship with the IFEM family of methods, in this section we focus on validation of the three advantages, as claimed in Section 1, of the one-field FDM compared with IFEM. We shall use the Least-squares method to approximate the convection step [22]. For the diffusion step in which FSI interaction happens, we use the finite element isoparametric interpolation functions to transfer information between the solid mesh and fluid mesh for both the one-field FDM and IFEM. The finite element interpolation function is suggested to be capable of producing sharper interfaces than the traditional discretized Dirac delta function or the reproducing kernel function in [8]. The pressure step is a “degenerate” Stokes equation, and we shall use the standard $Taylor$-$Hood$ element to discretize this step. Gravity will not be considered in this paper, so the Froude number (22) will be zero ($Fr=0$) in each of the following tests.

10.1 Lid-driven cavity flow with a solid disc

This example is taken from papers [31, 8], in which IFEM approaches are adopted. A sketch of the problem and boundary conditions are shown in Figure 2. We consider the parameter sets displayed in Table 1, and all the simulations use a sufficiently small time step to ensure stability: parameter set 1 (very soft solid) uses $\Delta t=1.0\times 10^{-3}$ and all other tests have a time step of $\Delta t=5.0\times 10^{-3}$. To illustrate the meshes that we use, velocity norms on the fluid mesh ($40\times 40$ biquadratic quadrilaterals) and solid mesh (1373 linear triangles with 771 nodes) for Parameter set 1 are presented in Figure 3.

For all these tests, we aim to run to $t=10$. However for the Parameter sets 3 to 5 ($\rho^{r}\neq 1$), our implementation of the implicit IFEM scheme cannot reach $t=10$. Therefore, in these cases, we show a comparison shortly before the IFEM breaks down. The following criterion is used to check whether the implicit IFEM converges.

[TABLE]

where ${\bf u}_{k}$ and ${\bf u}_{k+1}$ are the iterative values of the last and current step respectively, and $tol=10^{-6}$ is the convergence tolerance used in our tests.

The first two parameter sets are simple cases because $\rho^{r}=\mu^{r}=1$. We can see from Figure 4 that the one-field FDM and the explicit IFEM present almost identical results in the case of a very soft solid ($\tilde{c}_{1}=0.1$), both of which are themselves indistinguishable from the published results in the literature [8]. Figure 5 shows the solid deformation for a slightly harder disc ($\tilde{c}_{1}=1$). Although the explicit IFEM, implicit IFEM and the one-field FDM present very similar results, a close look at the interface shape shows that the solution of the one-field FDM is almost identical to that of the implicit IFEM, and different from the explicit IFEM.

The case of a larger solid density ($\rho^{r}=2$) and a smaller solid density ($\rho^{r}=0.5$) are tested by Parameter sets 3 and 4 respectively. Both results (Figures 6 and 7) show that the one-field FDM and the fully converged implicit IFEM have almost the same accuracy. Furthermore, in neither IFEM case do the results converge up to $t=10$ when using the same time step as the one-field FDM: $\Delta t=5.0\times 10^{-3}$. The explicit IFEM uses velocities from the previous two time steps to compute the effect of the solid: it can be seen from Figure 6 that the disc moves more slowly using this explicit IFEM. We also note that reducing the time step cannot solve the problem in this case, because the temporal term in the FSI force $\mathcal{F}_{t}^{FSI}$ (35) becomes huge and has a negative effect on the stability. Figure 7 demonstrates similar problems for the explicit IFEM, but the disc using the explicit IFEM moves faster than the one-field FDM or the implicit IFEM.

Parameter sets 5 ($\mu^{r}=1.5$) and 6 ($\mu^{r}=2$) are included to show the case of different viscosities between fluid and solid. It can be seen from Figure 8 that the results of the one-field FDM and the implicit IFEM match very well. Using the selected time step, our IFEM implementations break down after the first few time steps when testing Parameter set 6, therefore we only show the results of the one-field FDM in Figures 9 and 10 (as future test cases for others to validate against).

We purposely choose Parameter set 7 to be difficult, with large Reynolds number and differences in viscosity and density between fluid and solid. The one-field FDM can stably run up to $t=10$. We show the result in Figure 11 when the solid disc arrives at the top of the cavity.

Remark 5**.**

Notice that we have not considered the case of $\mu^{r}<1$, because we find that all the three methods (one-field FDM, explicit IFEM and implicit IFEM) may be unstable when $\Delta t\to 0$. However we shall not discuss this stability issue in more detail here as it is not the primary purpose of this paper. Please refer to [32, 33] for stability analysis.

10.2 Oscillating leaflet in a channel

This numerical example is taken from [17, 12, 16]. The computational domain is a channel of dimension $L\times H$, with a leaflet of dimension $h\times w$ initially located across the channel at its midpoint as shown in Figure 12. A periodic flow condition is prescribed on the inlet and outlet boundaries, given by

[TABLE]

with $T$ being the dimensionless period of the input flow and equal to 10. In this test, $L=4$, $H=1$, $h=0.8$ and $w=0.0212$.

The leaflet is approximated with 152 linear triangles with 116 nodes, and the fluid mesh is made up of $189\times 47$ biquadratic squares with 36005 nodes. We extend parameters ($\rho^{r}=1$) used in the above three publications to two more general cases as shown in Table 2. Using the first group of parameters, we demonstrate that the one-field FDM can use a time step of $\Delta t=5.0\times 10^{-3}$ while the explicit IFEM has to use a time step of $\Delta t=1.0\times 10^{-4}$ in order to remain stable (we reduce the time by a half over consecutive tests to check the convergence until finding a stable time step $\Delta t=7.8125\times 10^{-5}$, and then slightly increasing it we find $\Delta t=1.0\times 10^{-4}$ is also stable). However both simulations lead to almost identical results, as shown in Figure 13, which match the results of Fig. 3 in [17]. The reason for the difference in $\Delta t$ is due to the additional stabilizing terms added in the one-field FDM, as discussed in Remark 4. Also notice that we have to use the same time step for both our implicit and explicit IFEM schemes in order to converge for this example, because of the huge forcing term on the right-hand side in equation (36) and (37).

We then test a case with different density between fluid and solid: $\rho^{r}=1.2$. We use the same time step $\Delta t=1.0\times 10^{-4}$ for both the one-field FDM and the implicit IFEM, and their results, observed from Figure 14, are very close. The case of $\rho^{r}=2$ has also been tested, but our nonlinear implicit IFEM solver cannot converge at the first time step for any time step size. For completeness, we present the results using the proposed one-field FDM in Figures 15 and 16.

11 Conclusions

In this article we have illustrated the relationship between the recently proposed one-field FDM and both the explicit (explicitly expressing the force term) and implicit (implicitly expressing the force term) IFEM. This is facilitated through the use of a particular operator splitting scheme. Furthermore, we show that the one-field FDM produces the same accuracy of results as the implicit IFEM, but requires no iteration, whilst it significantly improves upon the classical IFEM at very little additional computational complexity. The one-field FDM is shown to naturally deal with the case of different densities and/or viscosities of the fluid and solid. Therefore, whilst we may view the scheme of [22] as a fictitious domain method, it is also legitimate to consider it to be an alternative, highly efficient and robust, approximate solution strategy for the fully implicit IFEM methods of [23, 24, 25]. Even in the simple case of the same density and viscosity, where the explicit IFEM is known to be successful, we find that the additional terms added in the formulation of the one-field FDM have a helpful stabilizing effect, such that a larger time step can be adopted compared with the explicit IFEM.

Appendix A Extension to the compressible neo-Hookean solid model

In this section, we extend the incompressible neo-Hookean solid model to a compressible case, in which case the constitutive equation can be expressed as [34]:

[TABLE]

where $\nu$ is the Poisson’s ratio. For a compressible solid, the continuity equation can simply be expressed as:

[TABLE]

where $\rho_{0}^{s}$ is the initial solid density. Then the corresponding weak forms (15) and (16) can be expressed as [34]:

[TABLE]

with $\rho^{\delta}=\rho_{0}^{s}/J-\rho^{f}$, and

[TABLE]

The boxed terms in the above equations indicate the differences compared with equation (15) and (16). For a compressible solid model, the incompressibility equation (2) ($\nabla\cdot{\bf u}=0$) only holds in the fluid domain $\Omega_{t}^{f}$. Therefore we cannot solve it in the whole domain $\Omega$ using a fictitious domain method, because this never matches the velocity of a compressible solid ($\nabla\cdot{\bf u}\neq 0$). The pressure computed in the solid domain $\left.p\right|_{\Omega_{t}^{s}}$ is meaningless, which is weakly imposed to be zero in (54) with $\kappa$ playing the role of a bulk modulus [34].

In order to implement the one-field FDM approach, after time discretization one could update the solid stress as described in Section 9. Alternatively, one could also update the deformation tensor as follows.

[TABLE]

Using (55), equation (53) and (54) may be solved implicitly, which can also use the operator spitting scheme introduced in Section 7. See [33] for more details about the implicit solver and different types of explicit splitting schemes.

Appendix B Extension to the compressible Saint Venant-Kirchhoff solid model

The constitutive equation of the Saint Venant-Kirchhoff solid model can be expressed as [35]:

[TABLE]

where

[TABLE]

is the Lagrangian Green strain, $\mu$ and $\lambda$ are the Lam$\acute{\rm e}$ constants. Then the corresponding weak forms (15) and (16) can be expressed as:

[TABLE]

with $\delta{\bf E}={\bf F}^{T}\left(\nabla_{\bf X}{\bf v}\right)+\left(\nabla_{\bf X}^{T}{\bf v}\right){\bf F}$ and $\rho^{\delta}=\rho_{0}^{s}/J-\rho^{f}$, and

[TABLE]

The term ${\bf S}:\delta{\bf E}$ may be linearized at a given displacement $\tilde{\bf d}$ as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

As with the previous implementation of the one-field FDM, one may update the solid displacement ${\bf d}$ after time discretizaiton as follows:

[TABLE]

One still can use the operator splitting scheme by choosing $\tilde{\bf d}={\bf d}_{n}$, although we omit the full details here.