A Monolithic and a Partitioned, Reduced Basis Method for Fluid–Structure Interaction Problems

Abstract

The aim of this work is to present an overview about the combination of the Reduced Basis Method (RBM) with two different approaches for Fluid–Structure Interaction (FSI) problems, namely a monolithic and a partitioned approach. We provide the details of implementation of two reduction procedures, and we then apply them to the same test case of interest. We first implement a reduction technique that is based on a monolithic procedure where we solve the fluid and the solid problems all at once. We then present another reduction technique that is based on a partitioned (or segregated) procedure: the fluid and the solid problems are solved separately and then coupled using a fixed point strategy. The toy problem that we consider is based on the Turek–Hron benchmark test case, with a fluid Reynolds number $$Re=100$$.

1. Introduction

The bridging between approximation techniques and high-performance computing finds numerous fields of applications in the industry as well as in academia: it is sufficient to think about heat transfer problems, electromagnetic problems, structural mechanics problems (linear/nonlinear elasticity), fluid problems, and acoustic problems. In all of these examples, the models are described using a system of partial differential equations (PDE) that usually depends on a given number of parameters that describe the geometrical configuration of the physical domain over which the problem is formulated or that describe some physical quantities (e.g., the Reynolds number for a fluid or the Lamè constants for a solid) or some boundary conditions. For all of these models, we usually focus on a particular quantity of interest, also called an output of interest, such as the maximum temperature of a system, a pressure drop, or a channel flowrate. Unfortunately, computing such an output for each new value of the parameter is a difficult task that is expensive both in terms of time computation and in terms of computer memory, even on modern HPC systems. With these premises in mind, it is clear why the Reduced Basis Method [1,2,3,4,5,6,7] (RBM) comes into play and shows a wide range of advantages: the idea at the core of the method is to simulate the behavior of the solution of our system of interest for some chosen values of the parameters in the PDE. This is usually performed using some well-established discretization technique, such as the Finite Element Method (FEM); another discretization method used, for example, in the compressible framework in computational fluid dynamics is the Finite Volume Method (FVM), and another possibility is the Cut Finite Element Method (CutFEM) (see for example [8,9,10]). Once we compute these solutions, in an expensive offline phase, we can use them to build some other basis functions: with these new basis functions, in the inexpensive online phase, we can approximate the solution of the system for a new value of the parameter.

Among the numerous applications of the RBM, Fluid–Structure Interaction (FSI) problems definitely represent a great challenge as well as an extremely interesting topic (see for example [4,11,12,13,14,15,16,17], just to cite a few). Indeed, despite their instrinsic complicated nature (see [18,19]), FSI problems are frequently used in everyday life: in naval engineering, they are used to study interactions between the water and the hull of a ship [20]; in biomedical applications, FSI problems are used to model the interaction between the blood flow and the deformable walls of a vessel [21,22,23,24,25,26,27]. Finally, in aeronautical engineering, FSI describes the way the air interacts with a plane or with (parts of) a shuttle; see [15,28,29,30].

The goal of this work is to present an extensive overview on the formulation of two model order reduction procedures that are applied to sets of snapshots obtained using two different approaches: a monolithic approach and a partitioned (segregated) approach. We present the entire formulation of the reduction procedures as well as some numerical results that were obtained using the two reduced order model techniques. The same test case is considered: the problem was inspired by the Turek–Hron benchmark test case FSI2, for which well-established results and analyses have already been presented in the literature; see for example [31,32]. The main difference in our work is that, for the structure, we consider a linear problem (linearized strain tensor), different from what that considered by Turek and Hron.

The rest of the work is structured as follows: in Section 2, we define the mathematical formalism behind coupled systems, and in particular, we introduce the Arbitrary Lagrangian Eulerian (ALE) formulation, which is used throughout the rest of the manuscript. In Section 3, we briefly introduce the two main approaches that can be adopted when dealing with FSI, namely monolithic and partitioned approaches. In Section 4, we present a monolithic reduced order model: in Section 4.1, we define the time discretization; in Section 4.2, we introduce the space discretization of the problem; and in Section 4.2, we discuss the imposition of coupling conditions through Lagrangian multipliers. In Section 4.3, we introduce the supremizer enrichment technique, and in Section 4.4, we present more in detail the implementation of the Proper Orthogonal Decomposition. In Section 4.5, we formulate the online reduced order coupled system, and finally in Section 5, we present some numerical results. Section 6 is devoted to a partitioned algorithm: in Section 6.1 and Section 6.2, we introduce the time discretization and the space discretization, respectively. In Section 6.3, we present the Proper Orthogonal Decomposition technique used in this case, and finally, in Section 6.4 and Section 7, we present the online reduced order problem and some numerical results. Finally, Section 8 is devoted to a discussion on the two algorithms.

2. Fluid–Structure Interaction Problems

In the following, we introduce the mathematical formalism for FSI problems: we assume a two dimensional setting. Therefore, let $\Omega \left(t\right)\subset {\mathbb{R}}^2$ be the physical domain of interest, at time $t\in [0,T]$. The physical domain can be naturally divided into two subdomains: $\Omega \left(t\right)={\Omega }_f\left(t\right)\cup {\Omega }_s\left(t\right)$, where ${\Omega }_f\left(t\right)$ is the fluid domain, and ${\Omega }_s\left(t\right)$, which is the solid domain; we further assume that ${\Omega }_f\left(t\right)\cap {\Omega }_s\left(t\right)=\varnothing$ and that ${\overline{\Omega }}_f\left(t\right)\cap {\overline{\Omega }}_s\left(t\right)={\Gamma }_{FSI}\left(t\right)$ is the fluid–structure interface: the left side of Figure 1 shows the fluid domain in blue and the solid domain in red.

The fluid is assumed to be Newtonian and incompressible, and therefore, its behavior can be modelled using incompressible Navier–Stokes equation: for every $t\in [0,T]$, find ${\mathit{u}}_f\left(t\right):{\Omega }_f\left(t\right)\mapsto {\mathbb{R}}^2$ and $p_f\left(t\right):{\Omega }_f\left(t\right)\mapsto \mathbb{R}$ such that

(1)$\left\{\begin{array}{c}{\rho }_f({\partial }_t{\mathit{u}}_f+(u_f·\nabla ){\mathit{u}}_f)-\mathrm{div}{\sigma }_f({\mathit{u}}_f,p_f)={\mathit{b}}_f\mathrm{in}{\Omega }_f\left(t\right)\times (0,T],\hfill \\ \mathrm{div}{\mathit{u}}_f=0\mathrm{in}{\Omega }_f\left(t\right)\times (0,T],\hfill \end{array}\right.$

where ${\mathit{b}}_f$ is the fluid volume external force and ${\sigma }_f({\mathit{u}}_f,p_f)$ is the Cauchy stress tensor that, given the fluid is Newtonian, can be expressed in the following way:

${\sigma }_f({\mathit{u}}_f,p_f)={\rho }_f{\nu }_f(\nabla {\mathit{u}}_f+{\nabla }^T{\mathit{u}}_f)-p_f\mathit{I}.$

where $\mathit{I}$ is the $2\times 2$ identity matrix, ${\rho }_f$ is the fluid density, and ${\nu }_f$ is the kinematic viscosity. We remark that, in this case, the differential operator ∇ is intended to be the differentiation with respect to the time-dependent spatial coordiantes and that the same is valid for the divergence operator $div$, which also is considered with respect to the time-dependent coordinates.

Throughout this manuscript, for the sake of simplicity of the exposition, the solid behavior is described using a linear elasticity equation; nevertheless, we remark that everything we say can also be applied to nonlinear models for the solid. The structure problem reads as follows: for every $t\in [0,T]$, find the solid displacement ${\widehat{\mathit{d}}}_s\left(t\right):{\widehat{\Omega }}_s\mapsto {\mathbb{R}}^2$ such that

(2)${\partial }_{tt}{\widehat{\mathit{d}}}_s-\widehat{\mathrm{div}}\widehat{P}\left({\widehat{\mathit{d}}}_s\right)={\widehat{\mathit{b}}}_s\mathrm{in}{\widehat{\Omega }}_s\times (0,T],$

where ${\widehat{\mathit{b}}}_s$ is the external force acting on the solid. We assume small deformations here: this means that the solid Piola–Kirchhoff stress tensor $\widehat{P}\left({\widehat{\mathit{d}}}_s\right)$ can be expressed as

$\widehat{P}\left({\widehat{\mathit{d}}}_s\right)=2{\mu }_s\widehat{\epsilon }\left({\widehat{\mathit{d}}}_s\right)+{\lambda }_s\mathrm{tr}\widehat{\epsilon }\left({\widehat{\mathit{d}}}_s\right)\mathit{I},$

with ${\mu }_s$ and ${\lambda }_s$ being the Lamé constants of the material and with $\widehat{\epsilon }\left({\widehat{\mathit{d}}}_s\right)$ being the linearized strain operator, which is defined as follows:

$\widehat{\epsilon }\left({\widehat{\mathit{d}}}_s\right):=\frac{1}{2}(\widehat{\nabla }{\widehat{\mathit{d}}}_s+{\widehat{\nabla }}^T{\widehat{\mathit{d}}}_s).$

In the previous equations, $\widehat{\nabla }$ is the gradient with respect to the coordinates in the reference configuration and $\widehat{\mathrm{div}}$ is the divergence operator, always in the reference configuration.

The first thing we notice from Equations (1) and (2) is that the two problems are formulated over two “different kind” of domains: indeed, the Navier–Stokes equation is formulated over a time-dependent domain ${\Omega }_f\left(t\right)$, whereas the linear elasticity equation is formulated over a time-independent domain ${\widehat{\Omega }}_s$. In this section, in order to clearly introduce the Arbitrary Lagrangian Eulerian formulation, we make use of the following notation: all fields that are defined on the reference, the time-independent configuration, are denoted with a hat$\widehat{}$; on the contrary, all of the fields that are defined on a time-dependent configuration are denoted without the hat. This distinction between the time-dependent configuration and the time-independent configuration is a rather natural peculiarity of FSI problems that arises from the different natures of the two problems involved: we refer the reader to Figure 1 for a graphical representation of both configurations. Solving a FSI problem in a time-dependent domain can be extremely costly from a computational point of view, as it requires remeshing at every time-step in order to update the entire configuration. One possible alternative to avoid remeshing is to try instead to solve the problems in a reference configuration. For the solid problem, the definition of the reference configuration is quite easy and natural: indeed, we can say that the solid reference configuration ${\widehat{\Omega }}_s$ is exactly the solid domain, undeformed: ${\widehat{\Omega }}_s={\Omega }_s(t=0)$. For the fluid problem, defining a reference configuration is instead rather complicated. To this aim, we introduce the ALE map

(3)$\begin{array}{cc}\hfill {\mathcal{A}}_f\left(t\right):{\widehat{\Omega }}_f& \mapsto {\Omega }_f\left(t\right)\hfill \\ \hfill \widehat{\mathit{x}}& \mapsto \mathit{x}=\widehat{\mathit{x}}+{\widehat{\mathit{d}}}_f\left(t\right),\hfill \end{array}$

which maps an arbitrary time-independent fluid domain to the fluid current configuration. In Equation (3), the mesh displacement ${\widehat{\mathit{d}}}_f\left(t\right)$ at time t is defined as an extension to the whole domain ${\widehat{\Omega }}_f$ of the solid deformation ${\widehat{\mathit{d}}}_s\left(t\right)$. This extension can be carried out in different ways. Here, we adopt an harmonic extension:

$\left\{\begin{array}{c}-\widehat{\mathrm{div}}\left(\frac{1}{J}\widehat{\nabla }{\widehat{\mathit{d}}}_f\right)=0,\mathrm{in}{\widehat{\Omega }}_f\times [0,T]\hfill \\ {\widehat{\mathit{d}}}_f={\widehat{\mathit{d}}}_s,\mathrm{on}{\widehat{\Gamma }}_{FSI}\times [0,T].\hfill \end{array}\right.$

where J is defined as $J\left(t\right):=\det F\left(t\right)$, where $F\left(t\right)$ is the gradient of the ALE map ${\mathcal{A}}_f\left(t\right)$. From now on, for ease of notation, we drop the time dependence of J and F. We complete the previous system with homogeneous Dirichlet boundary conditions on the whole $\partial {\widehat{\Omega }}_f$. In the previous system, we used the scaling factor $\frac{1}{J}$: the reason for doing this is that the mesh displacement obtained is sligthly more regular than the one we would obtain without the scaling factor, meaning that the deformation of the triangles of the mesh is more regular. We refer the reader to Chapter 5 of [33] for some interesting comparison results regarding the quality of deformation of the triangles of the mesh as well as for some alternative ways to define the mesh displacement (linear elasticity and biharmonic extension).

Thanks to the introduction of the ALE map, we can decide to take the domain at time $t=0$ as the arbitrary time-independent fluid domain; hence, ${\widehat{\Omega }}_f:={\Omega }_f(t=0)$. After some computations (see for example [33]), we can reformulate the Navier–Stokes problem in the reference configuration: for every $t\in [0,T]$, find the fluid velocity ${\widehat{\mathit{u}}}_f\left(t\right):{\widehat{\Omega }}_f\mapsto {\mathbb{R}}^2$, the fluid pressure ${\widehat{p}}_f\left(t\right):{\widehat{\Omega }}_f\mapsto \mathbb{R}$, and the fluid displacement ${\widehat{\mathit{d}}}_f\left(t\right):{\widehat{\Omega }}_f\mapsto {\mathbb{R}}^2$ such that

(4)$\left\{\begin{array}{c}{\rho }_{f}J({\partial }_t{\widehat{\mathit{u}}}_f+\widehat{\nabla }{\widehat{\mathit{u}}}_f{F}^{-1}({\widehat{\mathit{u}}}_f-{\partial }_t{\widehat{\mathit{d}}}_f))-\widehat{\mathrm{div}}\left(J{\widehat{\sigma }}_f({\widehat{\mathit{u}}}_f,{\widehat{p}}_f)F^{-T}\right)=J{\widehat{\mathit{b}}}_f\mathrm{in}{\widehat{\Omega }}_f\times (0,T],\hfill \\ \widehat{\mathrm{div}}\left(J{F}^{-1}{\widehat{\mathit{u}}}_f\right)=0\mathrm{in}{\widehat{\Omega }}_f\times (0,T],\hfill \\ -\widehat{\mathrm{div}}\left(\frac{1}{J}\widehat{\nabla }{\widehat{\mathit{d}}}_f\right)=0\mathrm{in}{\widehat{\Omega }}_f\times (0,T],\hfill \end{array}\right.$

where ${\widehat{b}}_f$ is now the external force acting on the fluid, expressed with respect to the coordinates in the reference configuration, and ${\widehat{\sigma }}_f({\widehat{\mathit{u}}}_f,{\widehat{p}}_f)$ is the fluid Cauchy stress tensor in the arbitrary reference configuration ${\widehat{\Omega }}_f$:

${\widehat{\sigma }}_f({\widehat{\mathit{u}}}_f,{\widehat{p}}_f):={\rho }_f{\nu }_f(\widehat{\nabla }{\widehat{\mathit{u}}}_f{F}^{-1}+F^{-T}{\widehat{\nabla }}^T{\widehat{\mathit{u}}}_f).$

Together, Equations (2) and (4) form the coupled FSI problem, expressed in a time-independent configuration. In order to have a coupled system, we still need some coupling conditions that describe the interaction between the two physics:

(5)$\left\{\begin{array}{c}{\widehat{\mathit{d}}}_f={\widehat{\mathit{d}}}_s\mathrm{on}{\widehat{\Gamma }}_{FSI}\times (0,T],\hfill \\ {\widehat{\mathit{u}}}_f={\partial }_t{\widehat{\mathit{d}}}_s\mathrm{on}{\widehat{\Gamma }}_{FSI}\times (0,T],\hfill \\ J{\widehat{\sigma }}_f{F}^{-T}{\widehat{\mathit{n}}}_f={\widehat{P}}_s{\widehat{\mathit{n}}}_s\mathrm{on}{\widehat{\Gamma }}_{FSI}\times (0,T],\hfill \end{array}\right.$

where ${\widehat{\mathit{n}}}_f$ is the normal to the fluid–structure interface (reference configuration) from the fluid domain and ${\widehat{\mathit{n}}}_s$ is the normal to the fluid–structure interface from the solid domain.

The entire coupled system is then completed by some initial conditions and by some boundary conditions. In the following, we work with these boundary conditions:

$\left\{\begin{array}{c}{\widehat{\mathit{u}}}_f=\overline{\mathit{u}}\left(t\right)\mathrm{on}{\widehat{\Gamma }}_D^f\times (0,T],\hfill \\ {\widehat{\mathit{d}}}_s=0\mathrm{on}{\widehat{\Gamma }}_D^s\times (0,T],\hfill \\ J{\widehat{\sigma }}_f({\widehat{\mathit{u}}}_f,{\widehat{p}}_f)F^{-T}\widehat{\mathit{n}}=0\mathrm{on}{\widehat{\Gamma }}_N^f\times (0,T].\hfill \end{array}\right.$

where ${\widehat{\Gamma }}_D^f:={\widehat{\Gamma }}_{in}\cup {\widehat{\Gamma }}_{walls}$ is the Dirichlet boundary and is made by the inlet boundary and by ${\widehat{\Gamma }}_{walls}$, which is the union of the top boundary, the bottom boundary, and the boundary of the cylinder immersed in the fluid. ${\widehat{\Gamma }}_N^f$ is the fluid Neumann boundary (outflow boundary), and $\widehat{\mathit{n}}$ is the normal to ${\widehat{\Gamma }}_N^f$ from the fluid domain; in addition $\overline{\mathit{u}}\left(t\right)$ is the Dirichlet data, and we have

(6)$\overline{\mathit{u}}\left(t\right)=\left\{\begin{array}{c}{\mathit{u}}_{in}\left(t\right)\mathrm{on}{\widehat{\Gamma }}_{in}\times (0,T],\hfill \\ 0\mathrm{on}{\widehat{\Gamma }}_{walls}\times (0,T].\hfill \end{array}\right.$

An example of the reference configuration together with an illustration of ${\widehat{\Gamma }}_{FSI}$, ${\widehat{\Gamma }}_D^f$, ${\widehat{\Gamma }}_D^s$, ${\widehat{\Gamma }}_N^f$, ${\widehat{\Gamma }}_{in}$, and ${\widehat{\Gamma }}_{walls}$ is represented in Figure 1.

3. Approaches to Fluid–Structure Interaction Problems

As we previously mentioned, FSI problems are characterized by the presence of two different physics interacting with one another: we have the fluid problem, represented by the Stokes or Navier–Stokes equation, and we have the structure problem, represented by a string equation, by linear elasticity, or by nonlinear elasticity. It is therefore almost natural to conclude that there are two main different routes one can take to address such problems. Indeed, we can decide to solve the fluid and the solid problems separately and then take care of the coupling between the two physics: this gives rise to the so-called partitioned (segregated) approach. On the other hand, we can decide to solve the two problems together, and this gives rise to a monolithic procedure instead.

3.1. Monolithic Approach

In a monolithic algorithm, the fluid and the solid problems are solved simultaneously. These kind of algorithms are more stable, and they usually allow for the use of bigger time-steps during the discretization of the problem in time. The main drawback is that one really relies on available ad hoc softwares that are capable of handling a computational fluid dynamics problem as well as a solid mechanics problem at the same time. In addition to this, in order to pursue a Galerkin discretization of the original problem, one needs to introduce Lagrange multipliers to impose the coupling conditions at the fluid–structure interface, thus increasing the number of unknowns in the coupled problem. For the reader interested in looking into more detail about monolithic algorithms, we refer to [11,17,33,34,35,36,37], even though this list is by no means complete.

3.2. Segregated Approach

The rationale behind a partitioned approach is to deal with the two physics separately: this indeed allows us to better exploit existing simulation tools for fluid dynamics and for structural dynamics, which are well developed nowadays and are used on a daily basis in industrial applications. A partitioned approach has many advantages: indeed, we have the possibility to combine different discretization tools for the two physics (e.g., finite volumes for the fluid and finite elements for the structure), we have the possibility to refine them for one of the two physics in time, as required by the situation. Unfortunately, in this case, there are also some drawbacks, as it turns out that, under some physical and geometrical conditions, partitioned algorithms are unstable; this situation may occur when the physical domain has a slender shape or when the fluid density ${\rho }_f$ is close to the solid density ${\rho }_s$ (this is almost always the case in hemodynamics applications, where the density of the blood is quite close to the density of the walls of the vessel). This instability occurs because of the well-known added mass effect: the fluid acts similar to an added mass to the solid, thus changing its natural behavior. The reader interested in the analysis behind the added mass effect and in its derivation is referred to [38]. With a partitioned procedure, we can give rise to a variety of different algorithms according to the strategy used to impose the coupling conditions at the fluid–structure interface:

• Explicit algorithms: after time discretization, the coupling conditions are treated explicitly at every time-step. These algorithms, also known as weakly or loosely coupled algorithms [39], are successfully applied in aerodynamics applications (see [40,41]), but some studies (see [38,42,43]) showed that they are unstable under some physical and geometrical conditions due to the added mass effect, as we previously mentioned.

• Implicit algorithms: in these algorithms, also known as strongly coupled algorithms, the coupling conditions are treated implicitly at every time-step; see for example [44,45]. This implicit coupling represents a way to circumvent the instability problems due to the added mass effect; nevertheless, an implicit treatment of the coupling conditions leads to algorithms that are more expensive in terms of computational time.

• Semi-implicit algorithms: in these algorithms (see [46,47,48]), the continuity of the displacement is treated explicitly whereas the other coupling conditions are treated implicitly. This alternative represents a tradeoff between the computational cost of the algorithm and its stability in relation to the physical and geometrical properties of the problem. In Section 6.4, we present a reduced order method that is based on this kind of partitioned approach.

4. Monolithic Approach

In the following, we propose a monolithic approach for coupled problems. As already mentioned in the Introduction, a monolithic approach means that we solve the fluid and the solid problem all together; as we see in Section 4.2, this leads to a more subtle treatment of the coupling conditions at the fluid–structure interface; at the same time, adopting a monolithic procedure allows us to have better control of the global behavior of the coupled system. Before going any further, we now remark that, from now on, unless otherwise stated, we assume that everything is formulated over the reference configuration. Therefore, for ease of notation, we drop the hat symbol over the variables.

4.1. Time Discretization

We discretize the time interval $[0,T]$ with equispaced time sample points, thus obtaining $\{0=t_0,\cdots ,t_{N_T}=T\}$, where $t_i=i\Delta T$ and where $\Delta T$ is the time-step chosen. In the following, we make use of the notation $f^i:=f\left(t^i\right)$ for any given function f.

We discretize the time derivatives in the fluid problem with a second-order backward difference formula (BDF2) (see [37]):

$D_t{\mathit{u}}_f^{i+1}=\frac{3}{2\Delta T}{\mathit{u}}_f^{i+1}-\frac{4}{2\Delta T}{\mathit{u}}_f^i+\frac{1}{2\Delta T}{\mathit{u}}_f^{i-1}.$

Using a Newmark scheme for the solid problem as suggested in [49], we can define the structure first- and second-order time derivatives:

$\begin{array}{cc}\hfill D_{tt}{\mathit{d}}_s^{i+1}& =\frac{1}{\beta \Delta T}({\mathit{d}}_s^{i+1}-{\mathit{d}}_s^i)-\frac{1}{\beta \Delta T}{D}_t{\mathit{d}}_s^i-\left(\frac{1}{2\beta }-1\right)D_{tt}{\mathit{d}}_s^i,\hfill \\ \hfill D_t{\mathit{d}}_s^{i+1}& =\frac{\gamma }{\beta \Delta T}({\mathit{d}}_s^{i+1}-{\mathit{d}}_s^i)-\left(\frac{\gamma }{\beta }-1\right)D_t{\mathit{d}}_s^i-\Delta T\left(\frac{\gamma }{2\beta }-1\right)D_{tt}{\mathit{d}}_s^i,\hfill \end{array}$

where the constants $\gamma$ and $\beta$ are chosen in order to ensure unconditional stability of the Newmark scheme, as suggested in [49].

4.2. Space Discretization

We now introduce the discrete version of the original problems (2)–(4) from a monolithic point of view. Let us define the following function spaces:

$\begin{array}{cc}\hfill V^f& :=\{{\mathit{u}}_f\in {\left[H^1\left({\Omega }_f\right)\right]}^2\mathrm{s}.\mathrm{t}.{\mathit{u}}_f=\overline{\mathit{u}}\mathrm{on}{\Gamma }_D^f\times (0,T]\},\hfill \\ \hfill V_0^f& :=\{{\mathit{v}}_f\in {\left[H^1\left({\Omega }_f\right)\right]}^2\mathrm{s}.\mathrm{t}.{\mathit{v}}_f=0\mathrm{on}{\Gamma }_D^f\times (0,T])\},\hfill \\ \hfill Q& :=\{q\in L^2\left({\Omega }_f\right)\},\hfill \\ \hfill E^f& :=\{{\mathit{d}}_f\in {\left[H^1\left({\Omega }_f\right)\right]}^2\mathrm{s}.\mathrm{t}.{\mathit{d}}_f=0\mathrm{on}\partial {\Omega }_f\times (0,T])\},\hfill \\ \hfill E^s& :=\{{\mathit{d}}_s\in {\left[H^1\left({\Omega }_s\right)\right]}^2\mathrm{s}.\mathrm{t}.{\mathit{d}}_s=0\mathrm{on}{\Gamma }_D^s\times (0,T]\}.\hfill \end{array}$

We consider the spaces $V^f$, $V_0^f$, and $E^f$ to be endowed with the $H^1$ seminorm ${(\nabla (·),\nabla (·))}_{{\Omega }_f}$; Q to be endowed with the $L^2$ norm; and the space $E^s$ to be endowed with the $H^1$ seminorm ${(\nabla (·),\nabla (·))}_{{\Omega }^s}$.

We discretize the FSI problem in space, using the inf–sup stable Taylor–Hood pair $(V_h^f,Q_h)=({\mathbb{P}}_2,{\mathbb{P}}_1)$ for the fluid velocity and the fluid pressure and, similarly, for the space $V_0^f$. We use second-order Lagrange finite elements obtaining the discrete space $E_h^f\subset E^f$ for the mesh displacement and, similarly, for the solid displacement, resulting in the discretized space $E_h^s\subset E^s$. The weak formulation of problems (2) and (4) now reads: find $({\mathit{u}}_{f,h},p_{f,h},{\mathit{d}}_{f,h},{\mathit{d}}_{s,h})\in V_h^f\times Q_h\times E_h^f\times E_h^s$ such that, for every $({\mathit{v}}_{f,h},q_{f,h},{\mathit{e}}_{f,h},{\mathit{e}}_{s,h})\in V_h^0\times Q_h\times E_h^f\times E_h^s$,

(7)$\left\{\begin{array}{c}{({\partial }_{tt}{\mathit{d}}_{s,h},{\mathit{e}}_{s,h})}_{{\Omega }_s}+{(P\left({\mathit{d}}_{s,h}\right),\nabla {\mathit{e}}_{s,h})}_{{\Omega }_s}-{(P\left({\mathit{d}}_{s,h}\right){\mathit{n}}_s,{\mathit{e}}_{s,h})}_{{\Gamma }_{FSI}}={({\mathit{b}}_s,{\mathit{e}}_{s,h})}_{{\Omega }_s},\hfill \\ {\rho }_{f}J{({\partial }_t{\mathit{u}}_{f,h},{\mathit{v}}_{f,h})}_{{\Omega }_f}+{\rho }_{f}J{(\nabla {\mathit{u}}_{f,h}{F}^{-1}{\mathit{u}}_{f,h},{\mathit{v}}_{f,h})}_{{\Omega }_f}-{\rho }_{f}J{(\nabla {\mathit{u}}_{f,h}{F}^{-1}{\partial }_t{\mathit{d}}_{f,h},{\mathit{v}}_{f,h})}_{{\Omega }_f}+\hfill \\ +{(J{\sigma }_f({\mathit{u}}_{f,h},p_{f,h})F^{-T},\nabla {\mathit{v}}_{f,h})}_{{\Omega }_f}-{(J{\sigma }_f({\mathit{u}}_{f,h},p_{f,h})F^{-T}{\mathit{n}}_f,{\mathit{v}}_{f,h})}_{{\Gamma }_{FSI}}={(J{\mathit{b}}_f,{\mathit{v}}_{f,h})}_{{\Omega }_f}\hfill \\ -{(\mathrm{div}\left(J{F}^{-1}{\mathit{u}}_{f,h}\right),q_{f,h})}_{{\Omega }_f}=0\hfill \\ {(\frac{1}{J}\nabla {\mathit{d}}_{f,h},\nabla {\mathit{e}}_{f,h})}_{{\Omega }_f}={(\frac{1}{J}\nabla {\mathit{d}}_{f,h}{\mathit{n}}_f,{\mathit{e}}_{f,h})}_{{\Gamma }_{FSI}}.\hfill \end{array}\right.$

In the previous system, ${\mathit{n}}_f$ and ${\mathit{n}}_s$ are the normals to the fluid–structure interface, from the fluid domain and the solid domain, respectively.

Coupling Conditions through Lagrange Multipliers

Let us analyze the boundary integrals appearing in system (7) a little bit more in detail. From condition (5), it is clear that we have to satisfy the following:

${(P\left({\mathit{d}}_{s,h}\right){\mathit{n}}_s,{\mathit{e}}_{s,h})}_{{\Gamma }_{FSI}}=-{(J{\sigma }_f({\mathit{u}}_{f,h},p_{f,h})F^{-T}{\mathit{n}}_f,{\mathit{v}}_{f,h})}_{{\Gamma }_{FSI}}.$

Now, let us have a look at condition (5). It is easy to see that, in the time-continuous regime, thanks to (5), the following are equivalent:

${\mathit{d}}_{f,h}={\mathit{d}}_{s,h}\mathrm{on}{\Gamma }_{FSI},{\partial }_t{\mathit{d}}_{s,h}={\mathit{u}}_{f,h}\mathrm{on}{\Gamma }_{FSI},{\partial }_t{\mathit{d}}_{s,h}={\partial }_t{\mathit{d}}_{f,h}\mathrm{on}{\Gamma }_{FSI}.$

This equivalence does not hold anymore in general after time discretization, as one may choose different time integration schemes for the fluid and for the solid. In this work, we chose to weakly enforce the two following conditions:

(8)${\mathit{d}}_{f,h}={\mathit{d}}_{s,h}\mathrm{on}{\Gamma }_{FSI},{\partial }_t{\mathit{d}}_{s,h}={\mathit{u}}_{f,h}\mathrm{on}{\Gamma }_{FSI}.$

In order to do so, let us introduce the space $L:={\left[H^{-\frac{1}{2}}\left({\Gamma }_{FSI}\right)\right]}^2$ first and let us consider a finite dimensional subspace $L_h\subset L$ (we use first-order Lagrange finite elements). Let us take a discretized Lagrangian multiplier field ${\mathit{\lambda }}_{u,h}\in L_h$ and its corresponding test function ${\mathit{\mu }}_{u,h}\in L_h$. We can then write

${({\mathit{\mu }}_{u,h},{\mathit{u}}_{f,h}-{\partial }_t{\mathit{d}}_{s,h})}_{{\Gamma }_{FSI}}=0,$

and, by identifying the Lagrange multiplier with the unknown surface traction, we can finally rewrite the surface integrals in (7) as follows:

$\begin{array}{cc}\hfill -{(P\left({\mathit{d}}_{s,h}\right){\mathit{n}}_s,{\mathit{e}}_{s,h})}_{{\Gamma }_{FSI}}& =-{({\mathit{\lambda }}_{u,h},{\mathit{e}}_{s,h})}_{{\Gamma }_{FSI}},\hfill \\ \hfill -{(J{\sigma }_f({\mathit{u}}_{f,h},p_{f,h})F^{-T}{\mathit{n}}_f,{\mathit{v}}_{f,h})}_{{\Gamma }_{FSI}}& =+{({\mathit{\lambda }}_{u,h},{\mathit{e}}_{s,h})}_{{\Gamma }_{FSI}}.\hfill \end{array}$

We treat the continuity of the displacements at the interface similarly: we take another Lagrangian multiplier field ${\mathit{\lambda }}_{d,h}\in L_h$ and its corresponding test function ${\mathit{\mu }}_{d,h}\in L_h$. We can then write

${({\mathit{\mu }}_{d,h},{\mathit{d}}_{f,h}-{\mathit{d}}_{s,h})}_{{\Gamma }_{FSI}}=0,$

and now we identify the Lagrange multiplier with the surface traction caused by the mesh displacement; thus, we can rewrite the surface integral in (7)${}_4$ as follows:

${(\frac{1}{J}\nabla {\mathit{d}}_{f,h}{\mathit{n}}_f,{\mathit{e}}_{f,h})}_{{\Gamma }_{FSI}}={({\mathit{\lambda }}_{d,h},{\mathit{e}}_{f,h})}_{{\Gamma }_{FSI}}.$

Finally, after the space and the time discretization, our monolithic coupled system reads as follows: for every $t^i$, $i=0,\cdots ,N_T$, find ${\mathit{u}}_{f,h}^{i+1}\in V_h^f$, $p_{f,h}^{i+1}\in Q_h$, ${\mathit{d}}_{f,h}^{i+1}\in E_h^f$, ${\mathit{d}}_{s,h}^{i+1}\in E_h^s$, ${\mathit{\lambda }}_{u,h}^{i+1}\in L_h$, and ${\mathit{\lambda }}_{d,h}^{i+1}\in L_h$ such that, for every ${\mathit{v}}_{f,h}\in V_h^0$, $q_{f,h}\in Q_h$, ${\mathit{e}}_{f,h}\in E_h^f$, ${\mathit{e}}_{s,h}\in E_h^s$, ${\mathit{\mu }}_{u,h}\in L_h$, and ${\mathit{\mu }}_{d,h}\in L_h$,

(9)$\left\{\begin{array}{c}{(D_{tt}{\mathit{d}}_{s,h}^{i+1},{\mathit{e}}_{s,h})}_{{\Omega }_s}+(P{({\mathit{d}}_{s,h}^{i+1},\nabla {\mathit{e}}_{s,h})}_{{\Omega }_s}+{({\mathit{\lambda }}_{u,h}^{i+1},{\mathit{e}}_{s,h})}_{{\Gamma }_{FSI}}={({\mathit{b}}_s,{\mathit{e}}_{s,h})}_{{\Omega }_s},\hfill \\ {\rho }_{f}J{(D_t{\mathit{u}}_{f,h}^{i+1},{\mathit{v}}_{f,h})}_{{\Omega }_f}+{\rho }_{f}J{(\nabla {\mathit{u}}_{f,h}^{i+1}{F}^{-1}{\mathit{u}}_{f,h}^{i+1},{\mathit{v}}_{f,h})}_{{\Omega }_f}-{\rho }_{f}J{(\nabla {\mathit{u}}_{f,h}{F}^{-1}{D}_t{\mathit{d}}_{f,h}^{i+1},{\mathit{v}}_{f,h})}_{{\Omega }_f}+\hfill \\ +{(J{\sigma }_f({\mathit{u}}_{f,h}^{i+1},p_{f,h}^{i+1})F^{-T},\nabla {\mathit{v}}_{f,h})}_{{\Omega }_f}-{({\mathit{\lambda }}_{u,h}^{i+1},{\mathit{v}}_{f,h})}_{{\Gamma }_{FSI}}={(J{\mathit{b}}_f,{\mathit{v}}_{f,h})}_{{\Omega }_f}\hfill \\ -{(\mathrm{div}\left(J{F}^{-1}{\mathit{u}}_{f,h}^{i+1}\right),q_{f,h})}_{{\Omega }_f}=0\hfill \\ {(\frac{1}{J}\nabla {\mathit{d}}_{f,h}^{i+1},\nabla {\mathit{e}}_{f,h})}_{{\Omega }_f}={({\mathit{\lambda }}_{d,h}^{i+1},{\mathit{e}}_{f,h})}_{{\Gamma }_{FSI}},\hfill \\ {({\mathit{u}}_{f,h}^{i+1}-D_t{\mathit{d}}_{s,h}^{i+1},{\mathit{\mu }}_{u,h})}_{{\Gamma }_{FSI}}=0,\hfill \\ {({\mathit{d}}_{f,h}^{i+1}-{\mathit{d}}_{s,h}^{i+1},{\mathit{\mu }}_{d,h})}_{{\Gamma }_{FSI}}=0.\hfill \end{array}\right.$

4.3. Lifting Function and Supremizer Enrichment

We now present the technique used to perform a compression on the set of snapshots obtained with the monolithic FE discretization in order to create a set of reduced basis functions. The first detail that we explain is the introduction of a lifting function for the fluid velocity ${\mathit{u}}_{f,h}$. The use of a lifting function is quite common in the RBM approach; see for example [1,50]: the advantage of this technique is represented by the fact that, sometimes, in the problem of interest, we have to deal with non-homogeneous Dirichlet boundary conditions, as in our case, where we remind the reader that we ask for ${\mathit{u}}_{f,h}\left(t\right)={\mathit{u}}_{in}\left(t\right)$ at the inlet boundary for every $t\in (0,T]$. From the implementation point of view, this condition does not present any problems during the offline phase, where we solve the FE discretization of the original problem. However, this non-homogeneous condition represents a problem during the solution of the online system. Indeed, if we perform a POD on a set of snapshots that satisfy the same inlet condition, we obtain a set of reduced basis functions for the fluid velocity that has a given value at the inlet boundary. Now, let us assume that we have solved the online system and, hence, have found the coefficients that allow us to represent the reduced velocity approximation as a linear combination of the velocity reduced basis functions. It is easy to imagine that a linear combination of these basis functions does not automatically satisfy the original inlet condition in general. A solution to this problem is represented by the introduction of a lifting function ${\ell }_h\left(t\right)\in V_h^f$ during the offline phase such that ${\ell }_h\left(t\right)={\mathit{u}}_{in}\left(t\right)$ on ${\Gamma }_D^f$ for every $t\in (0,T]$. By subtracting the lifting function to the fluid velocity snapshots ${\mathit{u}}_{f,h}$ before performing a POD, we obtain a new variable ${\mathit{u}}_{0,h}:={\mathit{u}}_{f,h}-{\ell }_h\in V_0^f$ that satisfies a homogeneous Dirichlet boundary condition also at the fluid inlet boundary ${\Gamma }_{in}$: these are the snapshots on which we perform a POD, thus obtaining basis functions that are all zero at the inlet boundary. We point out that the definition of the lifting function is not unique; for our problem, the lifting function ${\ell }_h$ is defined by solving a steady Stokes problem over ${\Omega }_f$ at every time-step, with a prescribed inlet velocity ${\mathit{u}}_{in}$ and with homogeneous Dirichlet boundary conditions on ${\Gamma }_{walls}\cup {\Gamma }_{FSI}$, but other choices are possible as well.

Another important technique that has been used in this manuscript is the supremizer enrichment technique, which is necessary in order to obtain a stable approximation of the fluid pressure also at the reduced order level. The main reason for the introduction of the supremizer enrichment is that, even if the FE spaces $(V_h^f,Q_h)$ satisfy the inf–sup condition (which guarantees that the Navier–Stokes problem is uniquely solvable with respect to the pressure; see for example [51,52]), this may not hold true anymore once we move to the reduced spaces generated by the reduced basis functions. With the introduction of the supremizer enrichment, therefore, we aim to construct a pair of reduced function spaces for the fluid velocity and the fluid pressure, that also satisfies the inf–sup condition. The supremizer variable ${\mathit{s}}_h\in V_h^0$, is defined by solving the following problem: find ${\mathit{s}}_h\in V_h^0$ such that

$-{(\mathrm{div}{\mathit{v}}_{f,h},p_{f,h})}_{{\Omega }_f}={(\nabla {\mathit{s}}_h,\nabla {\mathit{v}}_{f,h})}_{{\Omega }_f}\forall {\mathit{v}}_{f,h}\in V_h^0.$

In the previous equation, $p_{f,h}$ is the FE pressure solution of the Navier–Stokes problem whereas the right-hand side is the scalar product that defines the $H^1$ seminorm, which we consider for the velocity function space $V_h^0$. For the reader interested in the details about this technique as well as the computations that lead to the formulation of the previous problem, we refer to [50,53].

4.4. Reduced Basis Generation

Once we obtain the FE supremizer snapshots ${\mathit{s}}_h^i$, $i=0,\cdots ,N_T$, and once we homogenize the fluid velocity snapshots thanks to the lifting function, we are ready to generate a set of reduced basis by performing a compression by Proper Orthogonal Decomposition.

In order to perform a POD, we need two main ingredients: the matrices of the inner products and the snapshots matrices. First, we need to introduce the basis functions for the FE spaces that we consider. We define the following:

$\begin{array}{c}\hfill \{{\mathit{\phi }}_1^u,\cdots ,{\mathit{\phi }}_{{\mathcal{N}}_h^u}^u\}\mathrm{the}\mathrm{FE}\mathrm{basis}\mathrm{of}\mathrm{the}\mathrm{discretized}\mathrm{space}{V}_h^0,\\ \hfill \{{\phi }_1^p,\cdots ,{\phi }_{{\mathcal{N}}_h^p}^p\}\mathrm{the}\mathrm{FE}\mathrm{basis}\mathrm{of}\mathrm{the}\mathrm{discretized}\mathrm{space}{Q}_h,\\ \hfill \{{\mathit{\phi }}_1^{d_f},\cdots ,{\mathit{\phi }}_{{\mathcal{N}}_h^{d_f}}^{d_f}\}\mathrm{the}\mathrm{FE}\mathrm{basis}\mathrm{of}\mathrm{the}\mathrm{discretized}\mathrm{space}{E}_h^f,\\ \hfill \{{\mathit{\phi }}_1^{d_s},\cdots ,{\mathit{\phi }}_{{\mathcal{N}}_h^{d_s}}^{d_s}\}\mathrm{the}\mathrm{FE}\mathrm{basis}\mathrm{of}\mathrm{the}\mathrm{discretized}\mathrm{space}{E}_h^s,\\ \hfill \{{\mathit{\phi }}_1^{\lambda },\cdots ,{\mathit{\phi }}_{{\mathcal{N}}_h^{\lambda }}^{\lambda }\}\mathrm{the}\mathrm{FE}\mathrm{basis}\mathrm{of}\mathrm{the}\mathrm{discretized}\mathrm{space}{L}_h,\end{array}$

where ${\mathcal{N}}_h^u$ is the dimension of the FE space $V_h^0$, ${\mathcal{N}}_h^p$ is the dimension of the FE space $Q_h$, ${\mathcal{N}}_h^{d_f}$ is the dimension of the FE space $E_h^f$, ${\mathcal{N}}_h^{d_s}$ is the dimension of the FE space $E_h^s$, and ${\mathcal{N}}_h^{\lambda }$ is the dimension of the FE space $L_h$ (which we remember we used to approximate both the Lagrange multiplier ${\mathit{\lambda }}_u$ and the Lagrange multiplier ${\mathit{\lambda }}_d$). We begin by constructing the snapshots matrices ${\mathcal{S}}_u\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}$, ${\mathcal{S}}_s\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}$,${\mathcal{S}}_p\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}$, ${\mathcal{S}}_{d_f}\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}$, ${\mathcal{S}}_{d_s}\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}$, ${\mathcal{S}}_{{\lambda }_u}\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}$, and ${\mathcal{S}}_{{\lambda }_d}\in {\mathbb{R}}^{{\mathcal{N}}_h\times M}$ defined as follows:

$\begin{array}{ccc}\hfill {\mathcal{S}}_u& & =[{\mathit{u}}_{0,h}^1,\cdots ,{\mathit{u}}_{0,h}^{N_T},0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,0\cdots ,0,0,\cdots ,0],\hfill \\ \hfill {\mathcal{S}}_s& & =[{\mathit{s}}_h^1,\cdots ,{\mathit{s}}_h^{N_T},0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,0\cdots ,0,0,\cdots ,0],\hfill \\ \hfill {\mathcal{S}}_p& & =[0,\cdots ,0,p_{f,h}^1,\cdots ,p_{f,h}^{N_T},0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,0\cdots ,0],\hfill \\ \hfill {\mathcal{S}}_{d_f}& & =[0,\cdots ,0,0,\cdots ,0,{\mathit{d}}_{f,h}^1,\cdots ,{\mathit{d}}_{f,h}^{N_T},0,\cdots ,0,0,\cdots ,0,0,\cdots ,0],\hfill \\ \hfill {\mathcal{S}}_{d_s}& & =[0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,{\mathit{d}}_{s,h}^1,\cdots ,{\mathit{d}}_{s,h}^{N_T},0,\cdots ,0,0,\cdots ,0],\hfill \\ \hfill {\mathcal{S}}_{{\lambda }_u}& & =[0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,{\mathit{\lambda }}_{u,h}^1,\cdots ,{\mathit{\lambda }}_{u,h}^{N_T},0,\cdots ,0],\hfill \\ \hfill {\mathcal{S}}_{{\lambda }_d}& & =[0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,0,\cdots ,0,0,\cdots ,0{\mathit{\lambda }}_{d,h}^1,\cdots ,{\mathit{\lambda }}_{d,h}^{N_T}].\hfill \end{array}$

In the previous definition, we have that ${\mathcal{N}}_h={\mathcal{N}}_h^u+{\mathcal{N}}_h^p+{\mathcal{N}}_h^{d_f}+{\mathcal{N}}_h^{d_s}+2{\mathcal{N}}_h^{\lambda }$ is the sum of the dimensions of all FE spaces that we used for the FE approximation of each component of the solution of the FSI system and that $M=6N_T$.

Next, we need to define the inner product matrices $X_u$, $X_p$, $X_{d_s}$, $X_{d_f}$, $X_{{\lambda }_u}$, and $X_{{\lambda }_d}$, all belonging to ${\mathbb{R}}^{{\mathcal{N}}_h\times {\mathcal{N}}_h}$. These matrices are block diagonal matrices and have the following form:

$\begin{array}{c}\hfill X_u=\mathrm{diag}({\mathit{x}}_u,{\mathbf{0}}_p,{\mathbf{0}}_{d_f},{\mathbf{0}}_{d_s},{\mathbf{0}}_{{\lambda }_u},{\mathbf{0}}_{{\lambda }_d})\\ \hfill X_p=\mathrm{diag}({\mathbf{0}}_u,{\mathit{x}}_p,{\mathbf{0}}_{d_f},{\mathbf{0}}_{d_s},{\mathbf{0}}_{{\lambda }_u},{\mathbf{0}}_{{\lambda }_d})\\ \hfill X_{d_f}=\mathrm{diag}({\mathbf{0}}_u,{\mathbf{0}}_p,{\mathit{x}}_{d_f},{\mathbf{0}}_{d_s},{\mathbf{0}}_{{\lambda }_u},{\mathbf{0}}_{{\lambda }_d})\\ \hfill X_{d_s}=\mathrm{diag}({\mathbf{0}}_u,{\mathbf{0}}_p,{\mathbf{0}}_{d_f},{\mathit{x}}_{d_s},{\mathbf{0}}_{{\lambda }_u},{\mathbf{0}}_{{\lambda }_d})\\ \hfill X_{{\lambda }_u}=\mathrm{diag}({\mathbf{0}}_u,{\mathbf{0}}_p,{\mathbf{0}}_{d_f},{\mathbf{0}}_{d_s},{\mathit{x}}_{{\lambda }_u},{\mathbf{0}}_{{\lambda }_d})\\ \hfill X_{{\lambda }_d}=\mathrm{diag}({\mathbf{0}}_u,{\mathbf{0}}_p,{\mathbf{0}}_{d_f},{\mathbf{0}}_{d_s},{\mathbf{0}}_{{\lambda }_u},{\mathit{x}}_{{\lambda }_d}).\end{array}$

In the previous definition, for simplicity, we used the following notation: ${\mathbf{0}}_{*}\in {\mathbb{R}}^{{\mathcal{N}}_h^{*}\times {\mathcal{N}}_h^{*}}$ is a zero block of dimension ${\mathcal{N}}_h^{*}\times {\mathcal{N}}_h^{*}$, where $\ast \in \{u_f,p,d_f,d_s,{\lambda }_u,{\lambda }_d\}$. In addition to this, we have the following nonzero blocks:

$\begin{array}{ccc}& & {\left({\mathit{x}}_u\right)}_{i,j}={(\nabla {\mathit{\phi }}_i^u,\nabla {\mathit{\phi }}_j^u)}_{{\Omega }_f}\mathrm{for}i,j=1,\cdots ,{\mathcal{N}}_h^u,\hfill \\ & & {\left({\mathit{x}}_p\right)}_{i,j}={({\phi }_i^p,{\phi }_j^p)}_{{\Omega }_f}\mathrm{for}i,j=1,\cdots ,{\mathcal{N}}_h^p,\hfill \\ & & {\left({\mathit{x}}_{d_f}\right)}_{i,j}={(\nabla {\mathit{\phi }}_i^{d_f},\nabla {\mathit{\phi }}_j^{d_f})}_{{\Omega }_f}\mathrm{for}i,j=1,\cdots ,{\mathcal{N}}_h^{d_f},\hfill \\ & & {\left({\mathit{x}}_{d_s}\right)}_{i,j}={(\nabla {\mathit{\phi }}_i^{d_s},\nabla {\mathit{\phi }}_j^{d_s})}_{{\Omega }_s}\mathrm{for}i,j=1,\cdots ,{\mathcal{N}}_h^{d_s},\hfill \\ & & {\left({\mathit{x}}_{{\lambda }_u}\right)}_{i,j}={({\mathit{\phi }}_i^{\lambda },{\mathit{\phi }}_j^{\lambda })}_{{\Gamma }_{FSI}}\mathrm{for}i,j=1,\cdots ,{\mathcal{N}}_h^{\lambda },\hfill \\ & & {\left({\mathit{x}}_{{\lambda }_d}\right)}_{i,j}={({\mathit{\phi }}_i^{\lambda },{\mathit{\phi }}_j^{\lambda })}_{{\Gamma }_{FSI}}\mathrm{for}i,j=1,\cdots ,{\mathcal{N}}_h^{\lambda }.\hfill \end{array}$

As the reader may notice, the inner product matrices are very big, given the fact that ${\mathcal{N}}_h^u,\cdots ,{\mathcal{N}}_h^{\lambda }\gg 1$: this is because the structure of the Proper Orthogonal Decomposition itself reflects the fact that we use a monolithic approach to solve the FSI problem.

We are now able to define the correlation matrices ${\mathcal{C}}_u$, ${\mathcal{C}}_s$, ${\mathcal{C}}_p$, ${\mathcal{C}}_{d_f}$, ${\mathcal{C}}_{d_s}$, ${\mathcal{C}}_{{\lambda }_u}$, and ${\mathcal{C}}_{{\lambda }_d}$, all belonging to the space ${\mathbb{R}}^{M\times M}$:

$\begin{array}{cc}\hfill {\mathcal{C}}_u& :={\mathcal{S}}_u^T{X}_u{\mathcal{S}}_u\\ \hfill {\mathcal{C}}_s& :={\mathcal{S}}_s^T{X}_u{\mathcal{S}}_s\\ \hfill {\mathcal{C}}_p& :={\mathcal{S}}_p^T{X}_p{\mathcal{S}}_p\\ \hfill {\mathcal{C}}_{d_f}& :={\mathcal{S}}_{d_f}^T{X}_{d_f}{\mathcal{S}}_{d_f}\\ \hfill {\mathcal{C}}_{d_s}& :={\mathcal{S}}_{d_s}^T{X}_{d_s}{\mathcal{S}}_{d_s}\\ \hfill {\mathcal{C}}_{{\lambda }_u}& :={\mathcal{S}}_{{\lambda }_u}^T{X}_{{\lambda }_u}{\mathcal{S}}_{{\lambda }_u}\\ \hfill {\mathcal{C}}_{{\lambda }_d}& :={\mathcal{S}}_{{\lambda }_d}^T{X}_{{\lambda }_d}{\mathcal{S}}_{{\lambda }_d}\end{array}$

Remark: in the correlation matrices, all of the snapshots are defined on a common mesh. Indeed, as we mentioned at the beginning of the section, we dropped the hat notation$\widehat{}$, with the understanding that all of the quantities are defined on the common reference configuration. This aspect is extremely important, as mapping everything back onto a reference configuration greatly simplifies the implementation of the Proper Orthogonal Decomposition.

Once we built the correlation matrices, we carried out a POD compression on the set of snapshots, following for example [54]. We did this by solving the following (seven) eigenvalue problems:

(10)${\mathcal{C}}_{*}{\mathcal{Q}}_{*}={\mathcal{Q}}_{*}{\Lambda }_{*},$

where $*\in \{u_f$, s, p, $d_f$, $d_s$, ${\lambda }_u$, ${\lambda }_d\}$, ${\mathcal{Q}}_{*}$ is the eigenvectors matrix, and ${\Lambda }_{*}$ is the diagonal eigenvalues matrix. The kth reduced basis function related to problem (10) is obtained by applying the snapshots matrix ${\mathcal{S}}_{*}$ to the kth column of the matrix ${\mathcal{Q}}_{*}$; we therefore end up with the following basis functions (assume for simplicity $\ast =u_f$):

${\mathsf{\Phi }}_k^u:=\frac{1}{\sqrt{{\lambda }_k^u}}{\mathcal{S}}_u{\underline{\mathit{v}}}_k^u,$

where ${\lambda }_k^u$ is the eigenvalue corresponding to the eigenvector ${\mathit{v}}_k^u$. Similar definitions hold true for the other components of the solution, namely $p_f$, ${\mathit{d}}_s$, ${\mathit{d}}_f$, ${\mathit{\lambda }}_u$, and ${\mathit{\lambda }}_d$, as well as for the supremizer $\mathit{s}$. We refer the reader interested in more details about the POD to [1,50].

We therefore end up with the following set of reduced basis: $\{{\mathsf{\Phi }}_1^{s,u},\cdots ,{\mathsf{\Phi }}_{N_u}^{s,u},\cdots ,{\mathsf{\Phi }}_1^{{\lambda }_d},\cdots ,$${\mathsf{\Phi }}_{N_{{\lambda }_d}}^{{\lambda }_d}\}$, where each basis function is a block function of six components (one for each variable in the FSI problem):

${\mathsf{\Phi }}_k^{s,u}=\left(\begin{array}{c}{\Phi }_k^{s,u}\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right),\cdots ,{\mathsf{\Phi }}_k^{{\lambda }_d}=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ {\Phi }_k^{{\lambda }_d}\end{array}\right).$

where we have used the following notation: let $\{{\mathsf{\Phi }}_1^u,\cdots ,{\mathsf{\Phi }}_{N1}^u\}$ be the basis functions obtained by a compression by POD on the fluid velocity snapshots and let $\{{\mathsf{\Phi }}_1^s,\cdots ,{\mathsf{\Phi }}_{N2}^s\}$ be the basis functions obtained by running a POD on the supremizer snapshots. In order to carry on the supremizer enrichment technique, we consider the union of the two sets of basis functions $\{{\mathsf{\Phi }}_1^u,\cdots ,{\mathsf{\Phi }}_{N1}^u,{\mathsf{\Phi }}_1^s,\cdots ,{\mathsf{\Phi }}_{N2}^s\}$ and then denote by ${\mathsf{\Phi }}_k^{s,u}$ a generic element of the last set. We indicated ${\mathsf{\Phi }}_k^{s,u}$ in order to remark that the reduced basis functions for the fluid velocity consist of the reduced basis generated by the fluid velocity snapshots and of the reduced basis generated by the supremizer snapshots.

Finally, we introduce the reduced order finite dimensional space $V_N=\mathrm{span}\{{\mathsf{\Phi }}_1^{s,u},\cdots ,$${\mathsf{\Phi }}_{N_u}^{s,u},\cdots ,{\mathsf{\Phi }}_1^{{\lambda }_d},\cdots ,{\mathsf{\Phi }}_{N_{{\lambda }_d}}^{{\lambda }_d}\}$, with $N_u=N1+N2$ and $N=N_u+\cdots +N_{{\lambda }_d}$.

4.5. Online Phase

Once we have the reduced basis functions, we can define the reduced solution ${\mathit{u}}_N^{i+1}:=({\mathit{u}}_{0,N_u}^{i+1},p_{N_p}^{i+1},{\mathit{d}}_{N_{d_f}}^{i+1},{\mathit{d}}_{N_{d_s}}^{i+1},{\mathit{\lambda }}_{N_{{\lambda }_u}}^{i+1},{\mathit{\lambda }}_{N_{{\lambda }_d}}^{i+1})$ of our FSI problem:

(11)$\begin{array}{cc}\hfill {\mathit{u}}_{0,N_u}^{i+1}& :=\sum _{k=0}^{N_u}{\underline{u}}_{0,k}^{i+1}{\Phi }_k^{s,u},\hfill \end{array}$

(12)$\begin{array}{cc}\hfill p_{N_p}^{i+1}& :=\sum _{k=0}^{N_p}{\underline{p}}_k^{i+1}{\Phi }_k^p,\hfill \end{array}$

(13)$\begin{array}{cc}\hfill {\mathit{d}}_{N_{d_f}}^{i+1}& :=\sum _{k=0}^{N_{d_f}}{\underline{d}}_{f,k}^{i+1}{\Phi }_k^{d_f},\hfill \end{array}$

(14)$\begin{array}{cc}\hfill {\mathit{d}}_{N_{d_s}}^{i+1}& :=\sum _{k=0}^{N_{d_s}}{\underline{d}}_{s,k}^{i+1}{\Phi }_k^{d_s},\hfill \end{array}$

(15)$\begin{array}{cc}\hfill {\mathit{\lambda }}_{N_{{\lambda }_u}}^{i+1}& :=\sum _{k=0}^{N_{{\lambda }_u}}{\underline{\lambda }}_{u,k}^{i+1}{\Phi }_k^{{\lambda }_u},\hfill \end{array}$

(16)$\begin{array}{cc}\hfill {\mathit{\lambda }}_{N_{{\lambda }_d}}^{i+1}& :=\sum _{k=0}^{N_{{\lambda }_d}}{\underline{\lambda }}_{d,k}^{i+1}{\Phi }_k^{{\lambda }_d}.\hfill \end{array}$

In the previous equations, the underline bar indicates the vector of coefficients of the reduced solution; therefore, it indicates an element of $\mathbb{R}$ (for scalar components of the reduced solution) or ${\mathbb{R}}^2$ for vectorial components of the reduced solution (such as the fluid velocity for example). The online monolithic reduced order system reads as follows: for every $t^{i+1}$, $i=0,\cdots ,N_T-1$, find $u_N^{i+1}\in V_N$, $u_N^{i+1}=({\mathit{u}}_{0,N_u}^{i+1},p_{N_p}^{i+1},{\mathit{d}}_{N_{d_f}}^{i+1},{\mathit{d}}_{N_{d_s}}^{i+1},{\mathit{\lambda }}_{N_{{\lambda }_u}}^{i+1},{\mathit{\lambda }}_{N_{{\lambda }_d}}^{i+1})$ such that, for all $v_N\in V_N$, $v_N=({\mathit{v}}_{N_u},q_{N_p},{\mathit{e}}_{N_{d_f}},{\mathit{e}}_{N_{d_s}},{\mathit{\mu }}_{N_{{\lambda }_u}},{\mathit{\mu }}_{N_{{\lambda }_d}})$, the following holds

(17)$\left\{\begin{array}{c}{(D_{tt}{\mathit{d}}_{N_{d_s}}^{i+1},{\mathit{e}}_{N_{d_s}})}_{{\Omega }_s}+{(P\left({\mathit{d}}_{N_{d_s}}^{i+1}\right),\nabla {\mathit{e}}_{N_{d_s}})}_{{\Omega }_s}+{({\mathit{\lambda }}_{N_{{\lambda }_u}}^{i+1},{\mathit{e}}_{N_{d_s}})}_{{\Gamma }_{FSI}}={({\mathit{b}}_s,{\mathit{e}}_{N_{d_s}})}_{{\Omega }_s},\hfill \\ {\rho }_{f}J{(D_t({\mathit{u}}_{0,N_u}^{i+1}+{\ell }_{N_u}^{i+1}),{\mathit{v}}_{N_u})}_{{\Omega }_f}+{\rho }_{f}J{(\nabla ({\mathit{u}}_{0,N}^{i+1}+{\ell }_{N_u}^{i+1})F^{-1}({\mathit{u}}_{0,N}^{i+1}+{\ell }_{N_u}^{i+1}-D_t{\mathit{d}}_{N_{d_f}}^{i+1}),{\mathit{v}}_{N_u})}_{{\Omega }_f}+\hfill \\ +{(J{\sigma }_f({\mathit{u}}_{0,N}^{i+1}+{\ell }_{N_u}^{i+1},p_{N_p}^{i+1})F^{-T},\nabla {\mathit{v}}_{N_u})}_{{\Omega }_f}-{({\mathit{\lambda }}_{N_{{\lambda }_u}}^{i+1},{\mathit{v}}_{N_u})}_{{\Gamma }_{FSI}}={(J{\mathit{b}}_f,{\mathit{v}}_{N_u})}_{{\Omega }_f},\hfill \\ -{(\mathrm{div}(J{F}^{-1}{\mathit{u}}_{0,N}^{i+1}+{\ell }_{N_u}^{i+1}),q_{N_p})}_{{\Omega }_f}=0\hfill \\ {(\frac{1}{J}\nabla {\mathit{d}}_{N_{d_f}}^{i+1},\nabla {\mathit{e}}_{N_{d_f}})}_{{\Omega }_f}={({\mathit{\lambda }}_{N_{{\mathit{\lambda }}_d}}^{i+1},{\mathit{e}}_{N_{d_f}})}_{{\Gamma }_{FSI}},\hfill \\ {({\mathit{u}}_{0,N}^{i+1}+{\ell }_{N_u}^{i+1},{\mathit{\mu }}_{N_{{\lambda }_u}})}_{{\Gamma }_{FSI}}-{(D_t{\mathit{d}}_{N_{d_s}}^{i+1},{\mathit{\mu }}_{N_{{\lambda }_u}})}_{{\Gamma }_{FSI}}=0,\hfill \\ {({\mathit{d}}_{N_{d_f}}^{i+1}-{\mathit{d}}_{N_{d_s}}^{i+1},{\mathit{\mu }}_{N_{{\lambda }_d}})}_{{\Gamma }_{FSI}}=0.\hfill \end{array}\right.$

In the previous system, ${\ell }_{N_u}^{i+1}$ is the projection of the finite element lifting function ${\ell }_h^{i+1}$ on the finite dimensional space generated by the velocity reduced basis functions. Once we solve the online system, we can restore the reduced fluid velocity ${\mathit{u}}_{N_u}^{i+1}$ that satisfies the inlet condition ${\mathit{u}}_{N_u}^{i+1}={\mathit{u}}_{in}\left(t^{i+1}\right)$ on ${\Gamma }_{in}$ by using the relation ${\mathit{u}}_{N_u}^{i+1}={\mathit{u}}_{0,N_u}^{i+1}+{\ell }_{N_u}^{i+1}$.

5. Results

We now present some numerical results that were obtained by adopting a monolithic approach for the toy problem inspired by the Turek–Hron benchmark test case [31,32]; in particular, we refer to the test case FSI2 therein, which corresponds to a fluid with Reynolds number $Re=100$. We remark again before going any further that, while in the original benchmark problem the Green strain tensor was used for the solid, here, we used the linearized strain tensor, with the reason being that we are not interested in modelling large deformations; however, the values of the parameters are taken from the benchmark FSI2 presented in [31,32].

All of the numerical simulations for the offline phase were obtained with the use of multiphenics [55], whereas the online simulations were implemented with RBniCS [56].

Figure 2 represents the physical domain of the problem of interest. The channel has a length $L_f=2.5$ cm and a height of $h_f=0.41$ cm. The cylinder, which is assumed to be at rest and therefore is not considered as part of the solid domain, has center $C=(0.2,0.2)$ and a radius of $r=0.05$ cm. The deformable bar is $0.35$ cm long and $0.02$ cm thick.

In Table 1, we summarized the values of the physical parameters identifying the fluid and the solid behavior: as we can see, the test case corresponds to a fluid with Reynolds number $Re=100$, where $Re=\frac{\overline{U}2r}{{\nu }_f}$, and r is the radius of the immersed cylinder. We impose homogeneous Dirichlet boundary conditions around the cylinder and on the top and bottom walls of the cavity ${\Gamma }_{walls}$ for the fluid velocity and impose no conditions on ${\Gamma }_N^f$. We impose a velocity profile at the inlet boundary:

${\mathit{u}}_{in}(t,y):=\left\{\begin{array}{c}\overline{\mathit{u}}(0,y)\frac{1-\cos \left(\frac{\pi }{2}t\right)}{2}\mathrm{if}t<2 \mathrm{s}\hfill \\ \overline{\mathit{u}}(0,y)\mathrm{otherwise},\hfill \end{array}\right.$

where

$\overline{\mathit{u}}(0,y)=1.5\overline{U}\times \frac{4}{0.1681}y(0.41-y),$

and the value of $\overline{U}$ is reported in Table 1. We also require the bar to be attached to the cylinder; therefore, ${\mathit{d}}_s=0$ on ${\Gamma }_D^s$.

We use a time-step $\Delta t={10}^{-2}$ for the discretization in time, and the two constants $\gamma$ and $\beta$ used for the discretization of the structure time derivatives have the following values: $\gamma =0.25$ and $\beta =0.5$. The total number of iterations of the simulation is $N_T={10}^3$: all of these values are summarized in Table 2.

Since we are mostly interested in investigating the performance and the ability of the monolithic approach to reproducing the behavior of the coupled system, we adopted a mixed approach: we use a standard FE method, until the elastic bar starts to oscillate because of the action of the fluid; then, we run the reduced method. The oscillating behavior of the system takes approximatively $i=800$ iterations to occur: we therefore run the monolithic reduced order method for the remaining 200 iterations.

Figure 3a represents the behavior of the first 100 eigenvalues obtained with the POD on the snapshots of the monolithic system. The eigenvalues of the solid displacement have a faster decay with respect to the others, and this can be justified by the fact that the solid displacement behavior is periodic (the bar oscillates up and down) and that therefore the reduced order model is able to capture this periodic behavior with few modes. On the other hand the fluid velocity eigenvalues present a slower decay with respect to the other components: we expect this to be caused by the fact that, due to the periodic oscillation of the solid, we have the formation of some small vortices in the fluid that propagate into the domain, and this is a more complex phenomenon to reproduce with just a few modes. In Figure 3b, we can see the behavior of the energy $E_N$ retained by the first N modes for different components of the solution. Here, we give the definition of the retained energy for the fluid velocity component $u_f$, with the understanding that the energy retained for the modes of the other components of the FSI solution is defined in the exact same way:

$E_N^u=\frac{{\sum }_{k=1}^N\left|{\lambda }_k^u\right|}{{\sum }_{k=1}^{N_u}\left|{\lambda }_k^u\right|}.$

As we can see from the definition, the energy retained gives us some idea on the amount of information about the physical phenomenon that the first modes carry within: as it has to be expected by looking at the behavior of the eigenvalues, the first modes for the solid retain almost all information about the solid behavior, whereas for the fluid velocity, the energy retained slowly increases to 1.

Figure 4, Figure 5, Figure 6 and Figure 7 are intended to help the reader visualize a few outputs of the Proper Orthogonal Decomposition on the snapshots obtained with a monolithic approach. Figure 4 represents the first three modes for the fluid velocity and the first three modes for the supremizer enrichment: we observe how the fluid velocity modes are all zero at the inlet boundary thanks to the lifting function. Figure 5 represents the first three modes for the fluid pressure: as we can see, the modes present a highly oscillating behavior, thus suggesting that the supremizer may be needed in the online phase in order to get rid of any instability in the pressure approximation. Figure 6 and Figure 7 represent the modes for the mesh displacement and the solid displacement: we remark that, by looking at the two figures, we can see how the continuity of the displacement along the FSI interface and hence the need for the Lagrange multiplier also in the online phase of the algorithm are not automatically satisfied (as in the partitioned approach on the other hand) by the two sets of modes.

Figure 8 represents the deformation of the elastic bar at time $t=0.9$ s. The online solution was obtained using $N_{d_s}=21$ basis functions for the component $d_s$ of the solution of the FSI system. As we can see, with this number of reduced bases, the method is able to reconstruct the behavior of the solid component of the coupled system with good accuracy: indeed, the absolute value of the approximation error over the solid domain has a magnitude of ${10}^{-3}$, which, considered the magnitude of the solid deformation, represents a percentage error of $0.07\%$. In Figure 9 we can see the behavior of the fluid pressure again at time $t=0.9$ s. The online approximation of the fluid pressure was obtained by employing the supremizer enrichment technique, which requires the introduction of further basis functions in the fluid velocity reduced space: without this technique, the approximation of the fluid pressure becomes unstable and the whole algorithm diverges after a few time–steps, as we can see in Figure 10. For the fluid pressure, as we can see from Figure 9 bottom, the approximation error is good, and it represents a percentage error of $0.17\%$. Figure 11 represents the fluid velocity: as we can see, with a Reynolds number of 100, after some time, we have the developement of some Karman vortices that propagate into the fluid domain: the solid bar starts to oscillate and the whole system aquires a periodical behavior. The monolithic algorithm is capable of capturing and reproducing these complex phenomena, such as the Karman vortices in the fluid. Figure 11 shows the behavior of the approximation error: we observe that the error is mostly localized in the region of the domain that is close to the two main vortices that are detached from the solid bar, as it has to be expected, since this represents, from the fluid point of view, the most difficult physical aspect to reproduce.

Finally, Figure 12 represents the behavior of the average relative error of approximation for the different components of interest of the FSI problem: the average is taken over the number of time-steps. The relative error of approximation is computed in the norm considered in this manuscript for each component of the solution; hence, the $H^1$ seminorm over ${\Omega }_f$ for ${\mathit{u}}_f$ and ${\mathit{d}}_f$, the $L^2$ norm over ${\Omega }_f$ for $p_f$, and the $H^1$ seminorm over ${\Omega }^s$ for ${\mathit{d}}_s$. We would like to remark that, for this error analysis, we kept the number of reduced basis functions for the Lagrange multipliers fixed, in this specific case, to $N_{{\lambda }_u}=N_{{\lambda }_d}=5$. Indeed, we observed an increased presence of stability issues of the online algorithm as we increased the number of modes for the approximation of the Lagrange multipliers. One possible explanation to this is the following: in the online system (see Equation (17)), the reduced Lagrange multiplier ${\mathit{\lambda }}_{N_{{\lambda }_u}}$ represents the surface traction created by the homogenized fluid velocity $u_{0,N}$ and not by the fluid velocity $u_N$. For this reason, the FE Lagrange multiplier and the reduced order Lagrange multiplier have a different physical interpretation. We suspect this may be the source of instabilities arising by increasing the number of modes for the multipliers: further investigations on this subject will be carried out as future steps.

It is interesting to see that the average error of approximation decreases up to a certain number of basis function: in this specific case, it decreases up until, more or less, 20 reduced bases are used. After this threshold, the average relative approximation error increases and reaches a plateau, which indicates the fact that we have no gain in adding more reduced basis functions in the simulation. The increase in the average approximation error indicates that we reached a number of modes after which, if we add more basis functions, we just add noise to the online system. The same behavior is observed in Figure 13, which represents the average approximation error of the solid stress at the fluid–structure interface, where the average is taken over the number of time–steps, and the error is computed in the $L^2$ norm: again, the approximation error increases if we use a number of reduced basis that is greater than 20, confirming the fact that we are just adding noise to the system.

6. Partitioned Approach

In this section, we propose an alternative approach that is instead based on a segregated procedure: the idea is now to solve the fluid and the solid problems separately and to couple the two physics through some iterative procedure. As we will see, this idea leads to some advantages from the reduction point of view, and it allows us to work without the employment of additional Lagrange multipliers for the coupling conditions at the fluid–structure interface. The procedure that we propose is based on a Chorin–Temam projection scheme for the incompressible Navier–Stokes equations [57,58]; in addition, we employ a semi-implicit treatment of the coupling conditions [12,47,59].

6.1. Time Discretization

We first present the partitioned scheme after time discretization. $\Delta T$ is the time-step: we also employ an equispaced discretization of the time interval $[0,T]$ in this case. For discretization of the time derivative of a function f, we now use first backward difference BDF1:

$D_t{f}^{i+1}=\frac{f^{i+1}-f^i}{\Delta T},D_{tt}{f}^{i+1}=D_t\left(D_t{f}^{i+1}\right).$

The partitioned algorithm reads as follows for $i=0,\cdots ,N_T$:

• Extrapolation of the mesh displacement: find ${\mathit{d}}_f^{i+1}:{\Omega }_f\mapsto {\mathbb{R}}^2$ such that

  (18)$\left\{\begin{array}{c}-\Delta {\mathit{d}}_f^{i+1}=0\mathrm{in}{\Omega }_f,\hfill \\ {\mathit{d}}_f^{i+1}={\mathit{d}}_s^i\mathrm{on}{\Gamma }_{FSI}.\hfill \end{array}\right.$

• Fluid explicit step: find ${\mathit{u}}_f^{i+1}:{\Omega }_f\mapsto {\mathbb{R}}^2$ such that

  (19)$\left\{\begin{array}{c}J{\rho }_f\left(\frac{{\mathit{u}}_f^{i+1}-{\mathit{u}}_f^i}{\Delta T}+\nabla {\mathit{u}}_f^{i+1}{F}^{-1}({\mathit{u}}_f^{i+1}-D_t{\mathit{d}}_f^{i+1})\right)-\hfill \\ -{\rho }_f{\nu }_f\mathrm{div}\left(J\epsilon \left({\mathit{u}}_f^{i+1}\right)F^{-T}\right)+J{F}^{-T}\nabla p_f^i=J{\mathit{b}}_f\mathrm{in}{\Omega }_f,\hfill \\ \\ {\mathit{u}}_f^{i+1}=D_t{\mathit{d}}_f^{i+1}\mathrm{on}{\Gamma }_{FSI},\hfill \end{array}\right.$

  subject to the Dirichlet boundary condition ${\mathit{u}}_f^{i+1}=\overline{\mathit{u}}\left(t^{i+1}\right)$ on ${\Gamma }_D^f$, where $\overline{\mathit{u}}\left(t\right)$ is defined as in Equation (6). In the above system, $\epsilon \left({\mathit{u}}_f^{i+1}\right)$ is defined as follows:

  $\epsilon \left({\mathit{u}}_f^{i+1}\right):=\nabla {\mathit{u}}_f^{i+1}{F}^{-1}+F^{-T}{\nabla }^T{\mathit{u}}_f^{i+1}.$

• Implicit step:

  • Fluid projection substep (pressure Poisson formulation): find $p_f^{i+1}:{\Omega }_f\mapsto {\mathbb{R}}^2$ such that

    (20)$\left\{\begin{array}{c}-\mathrm{div}(J{F}^{-1}{F}^{-T}\nabla p_f^{i+1})=-\frac{{\rho }_f}{\Delta t}\mathrm{div}\left(J{F}^{-1}{\mathit{u}}_f^{i+1}\right)\mathrm{in}{\Omega }_f,\hfill \\ -F^{-T}\nabla p_f^{i+1}·J{F}^{-T}{\mathit{n}}_f={\rho }_f{D}_{tt}{\mathit{d}}_s^{i+1}·J{F}^{-T}{\mathit{n}}_f\mathrm{on}{\Gamma }_{FSI},\hfill \end{array}\right.$

    subject to the boundary conditions:

    (21)$p_f^{i+1}=\overline{p}\mathrm{on}{\Gamma }_{in},$

    where $\overline{p}$ is a prescribed pressure that we computed: a more detailed discussion about this aspect is presented after the structure problem substep.

  • Structure projection substep: find ${\mathit{d}}_s^{i+1}:{\Omega }_s\mapsto {\mathbb{R}}^2$ such that

    (22)$\left\{\begin{array}{c}{\rho }_s{D}_{tt}{\mathit{d}}_s^{i+1}-\mathrm{div}P\left({\mathit{d}}_s^{i+1}\right)={\mathit{b}}_s\mathrm{in}{\Omega }_s,\hfill \\ -P\left({\mathit{d}}_s^{i+1}\right){\mathit{n}}_s=J{\sigma }_f({\mathit{u}}_f^{i+1},p_f^{i+1})F^{-T}{\mathit{n}}_f\mathrm{on}{\Gamma }_{FSI},\hfill \end{array}\right.$

    subject to the boundary condition ${\mathit{d}}_s^{i+1}=0$ on ${\Gamma }_D^s$.

Before going any further, we briefly summarize some important remarks of this formulation of the original FSI problem:

• The original Navier–Stokes problem was divided into two subproblems, namely the fluid explicit step and the fluid projection step. In the explicit step, we take care of the momentum balance of the fluid problem, whereas in the projection step, we take care of the divergence free condition: this subdivision is a peculiarity of the Chorin–Temam projection scheme; we refer the interested reader to [57,58]. The advantage of adopting such a numerical scheme for the fluid problem is given by the fact that, in this case, we can also use pairs of discrete spaces ($V_h$ and $Q_h$) for the fluid velocity and pressure that do not necessarily satisfy the inf–sup condition; this represents a great advantage for the forthcoming online phase of the method, since we are able to obtain a stable approximation of the fluid pressure also without employing the supremizer enrichment technique, unlike in the monolithic approach.

• The treatment of the boundary conditions (in our specific case, the inlet boundary condition) is a delicate aspect of partitioned schemes. If the original problem of interest is provided with a boundary condition for the fluid Cauchy stress tensor ${\sigma }^f({\mathit{u}}_f,p_f){\mathit{n}}_{in}={\mathit{g}}_{in}$, where ${\mathit{n}}_{in}$ is the normal outgoing the inlet boundary, then during the Chorin–Temam projection scheme, this condition can be splitted into two: a natural condition for the fluid velocity explicit step $\epsilon \left({\mathit{u}}_f\right){\mathit{n}}_{in}=\mathbf{0}$ and a Dirichlet condition for the pressure $p_f{\mathit{n}}_{in}={\mathit{g}}_{in}$ (see for example [58]). However, in our particular toy problem, we have a Dirichlet inlet condition for the velocity: we therefore need a Dirichlet boundary condition also for the pressure in order to obtain uniqueness of solution of the pressure Poisson problem. In this case, we have no specific indication of what value for the pressure to choose at the inlet boundary: for our test case, we decided to compute the inlet pressure value by computing the quantity ${\sigma }^f({\mathit{u}}_f,p_f){\mathit{n}}_{in}$ on the inlet boundary and use the fact that, there, we have ${\mathit{u}}_f={\mathit{u}}_{in}$.

• In the projection step (20), we chose a pressure Poisson formulation; it is possible to use a Darcy formulation instead: find $p_f^{i+1}$ and ${\tilde{\mathit{u}}}_f^{i+1}$ such that

  $\left\{\begin{array}{c}{\rho }_{f}J\frac{{\tilde{\mathit{u}}}_f^{i+1}-u_f^{i+1}}{\Delta T}+J{F}^{-T}\nabla p_f^{i+1}=0\mathrm{in}{\Omega }_f,\hfill \\ \mathrm{div}\left(J{F}^{-1}{\tilde{\mathit{u}}}_f^{i+1}\right)=0\mathrm{in}{\Omega }_f.\hfill \end{array}\right.$

  However, in view of an efficient model order reduction, we chose to employ a Poisson formulation, since the Darcy formulation requires the introduction of an additional unknown ${\tilde{\mathit{u}}}_f$, which translates in a larger system, comprised of both velocity and pressure, at the implicit step.

In order to enhance the stability of the projection scheme, we employ Robin–Neumann coupling, as proposed in [12,60]; for other references on this kind of coupling, we refer to [61,62]. We thus replace condition (20) with the following:

(23)${\alpha }_{ROB}{p}^{i+1}+F^{-T}\nabla p^{i+1}·J{F}^{-T}{\mathit{n}}_f={\alpha }_{ROB}{p}^{i+1,☆}-{\rho }_f{D}_{tt}{\mathit{d}}_s^{i+1,☆}·J{F}^{-T}{\mathit{n}}_f.$

In Equation (23), $p^{i+1,☆}$ and ${\mathit{d}}_s^{i+1,☆}$ are suitable extrapolations of the fluid pressure and the solid displacement, respectively; we show in the next paragraph which kind of extrapolation we use. The constant ${\alpha }_{ROB}$ is defined as ${\alpha }_{ROB}=\frac{{\rho }_f}{z_p\Delta T}$, where $z_p$ is called the solid impedance:

$\begin{array}{cc}\hfill z_p& ={\rho }_s{c}_p,\hfill \\ \hfill c_p& =\sqrt{\frac{{\lambda }_s+2{\mu }_s}{{\rho }_s}}.\hfill \end{array}$

6.2. Space Discretization

We now aim to provide the final formulation of the partitioned problem; in order to do so, we need to discretize in space the original problem. We consider hereafter the same function spaces $V^f$, $V_0^f$, $E^f$, and $E^s$ that have been defined in Section 4.2. As far as the fluid pressure is concerned, with the Chorin–Temam projection scheme, the solution of the Poisson problem is now in $H^1$, and therefore, we introduce the following pressure function spaces:

$\begin{array}{cc}\hfill Q& :=\{q\in H^1\left({\Omega }_f\right)\mathrm{st}.q=\overline{p}\mathrm{on}{\Gamma }_{in}\},\hfill \\ \hfill Q_0& :=\{q\in H^1\left({\Omega }_f\right)\mathrm{st}.q=0\mathrm{on}{\Gamma }_{in}\}.\hfill \end{array}$

Additionally, in this case, we discretize the FSI problem in space using second-order Lagrange finite elements for the fluid velocity, the fluid displacement, and the solid displacement, resulting in the discrete spaces $V_h^f\subset V^f$, $V_h^0\subset V_0^f$, $E_h^f\subset E^f$, and $E_h^s\subset E^s$, while the fluid pressure is discretized with first-order Lagrange finite elements, resulting in the discrete spaces $Q_h\subset Q$, $Q_h^0\subset Q_0$. We are now ready to present the weak formulation of the original problem for every $i=0,\cdots ,N_T$:

• Extrapolation of the mesh displacement: find ${\mathit{d}}_{f,h}^{i+1}\in E_h^f$ such that $\forall {\mathit{e}}_{f,h}\in E_h^f$:

  (24)$\left\{\begin{array}{c}{(\nabla {\mathit{d}}_{f,h}^{i+1},\nabla {\mathit{e}}_{f,h})}_{{\Omega }_f}=0\hfill \\ {\mathit{d}}_{f,h}^{i+1}={\mathit{d}}_{s,h}^i\mathrm{on}{\Gamma }_{FSI}.\hfill \end{array}\right.$

• Fluid explicit step: find ${\mathit{u}}_{f,h}^{i+1}\in V_h^f$ such that $\forall {\mathit{v}}_{f,h}\in V_h^0$:

  (25)$\left\{\begin{array}{c}\begin{array}{cc}\hfill & {\rho }_f{(J\left(\frac{{\mathit{u}}_{f,h}^{i+1}-{\mathit{u}}_{f,h}^i}{\Delta T}\right),{\mathit{v}}_{f,h})}_{{\Omega }_f}+{\rho }_f{(J(\nabla {\mathit{u}}_{f,h}^{i+1}{F}^{-1}({\mathit{u}}_{f,h}^{i+1}-D_t{\mathit{d}}_{f,h}^{i+1})),{\mathit{v}}_{f,h})}_{{\Omega }_f}\hfill \\ \hfill & +{\rho }_f{\nu }_f{(J\epsilon \left({\mathit{u}}_{f,h}^{i+1}\right)F^{-T},\nabla {\mathit{v}}_{f,h})}_{{\Omega }_f}+{(J{F}^{-T}\nabla p_{f,h}^i,{\mathit{v}}_{f,h})}_{{\Omega }_f}=(J{\mathit{b}}_f,{\mathit{v}}_{f,h})\hfill \end{array}\hfill \\ {\mathit{u}}_{f,h}^{i+1}=D_t{\mathit{d}}_{f,h}^{i+1}\mathrm{on}{\Gamma }_{FSI},\hfill \end{array}\right.$

• Implicit step: for any $j=0,\dots$ until convergence:

  • Fluid projection substep (pressure Poisson formulation): find $p_{f,h}^{i+1,j+1}\in Q_h$ such that $\forall q_{f,h}\in Q_h^0$:

    (26)$\begin{array}{cc}& -\frac{{\rho }_f}{\Delta T}{(\mathrm{div}\left(J{F}^{-1}{\mathit{u}}_{f,h}^{i+1}\right),q_{f,h})}_{{\Omega }_f}-{\rho }_f{(\left(D_{tt}{\mathit{d}}_{s,h}^{i+1,j}\right),J{F}^{-T}{\mathit{n}}_f{q}_{f,h})}_{{\Gamma }_{FSI}}+\hfill \\ & +{\alpha }_{ROB}{(p_{f,h}^{i+1,j},q_{f,h})}_{{\Gamma }_{FSI}}={\alpha }_{ROB}{(p_{f,h}^{i+1,j+1},q_{f,h})}_{{\Gamma }_{FSI}}+{(J{F}^{-T}\nabla p_{f,h}^{i+1,j+1},F^{-T}\nabla q_{f,h})}_{{\Omega }_f}.\hfill \end{array}$

  • Structure projection substep: find ${\mathit{d}}_{s,h}^{i+1,j+1}\in E_h^s$ such that $\forall {\mathit{e}}_{s,h}\in E_h^s$:

    (27)${\rho }_s{(D_{tt}{\mathit{d}}_{s,h}^{i+1,j+1},{\mathit{e}}_{s,h})}_{{\Omega }_f}+{(P\left({\mathit{d}}_{s,h}^{i+1,j+1}\right),\nabla {\mathit{e}}_{s,h})}_{{\Omega }_s}=-{(J{\sigma }_f({\mathit{u}}_{f,h}^{i+1},p_{f,h}^{i+1,j+1})F^{-T}{\mathit{n}}_f,{\mathit{e}}_{s,h})}_{{\Gamma }_{FSI}}+{({\mathit{b}}_s,{\mathit{e}}_{s,h})}_{{\Omega }_s}$

    subject to the boundary condition ${\mathit{d}}_{s,h}^{i+1,j+1}=0$ on ${\Gamma }_D^s$.

We iterate between the two implicit substeps using a fixed point strategy:

(28)$\max \left(\frac{\left|\right|p_{f,h}^{i+1,j+1}-p_{f,h}^{i+1,j}{\left|\right|}_{Q_h}}{\left|\right|p_{f,h}^{i+1,j+1}{\left|\right|}_{Q_h}};\frac{\left|\right|{\mathit{d}}_{s,h}^{i+1,j+1}-{\mathit{d}}_{s,h}^{i+1,j}{\left|\right|}_{E_h^s}}{\left|\right|{\mathit{d}}_{s,h}^{i+1,j+1}{\left|\right|}_{E_h^s}}\right)<\epsilon ,$

where $\epsilon$ is a fixed tolerance.

In the pressure Poisson formulation, to impose the Robin coupling condition, we chose the pressure at the previous implicit iteration, namely $p_f^{i+1,j}$, as an extrapolation for the fluid pressure, and the same goes for the extrapolation of the structure displacement.

6.3. Reduced Basis Generation

For generation of the reduced basis for the fluid velocity ${\mathit{u}}_{f,h}$ and the fluid displacement ${\mathit{d}}_{f,h}$ we pursue the idea that was first proposed in [12]. For the solid displacement ${\mathit{d}}_{s,h}$, we employ a standard POD; for fluid pressure $p_{f,h}$, we first introduce a lifting function ${\ell }_p\left(t\right)$ and obtain the homogenized pressure $p_{0,h}\left(t\right)=p_{f,h}\left(t\right)-{\ell }_p\left(t\right)$ such that $p_{0,h}\left(t\right)=0$ on ${\Gamma }_{in}\times (0,T]$ and then we perform a standard POD. We therefore define the reduced pressure space $Q_N^0:=\mathrm{span}\{{\Phi }_1^p,\cdots ,{\Phi }_{N_p}^p\}$.

6.3.1. Change of Variable for the Fluid Velocity

The main idea here is to introduce a change of variable for ${\mathit{u}}_{f,h}$ in the fluid problem in order to transform condition (25) into a homogeneous boundary condition. The reason for this is similar to the reason why we introduced the lifting function in the first place: it is more convinient for the sake of the online system to work with homogeneous boundary conditions. Indeed, even if the reduced basis functions satisfy condition (25), it is not guaranteed, in general, that an element of the linear space generated by these reduced basis will also satisfy condition (25). For this reason, we need to introduce a Lagrange multiplier in order to make sure that Equation (25) is satisfied also at the reduced order level. Therefore, in order to avoid this and in order to design a more efficient reduced method, we chose to transform the non-homogeneous coupling condition into a homogeneous one. First, we transform the non–homogeneous inlet boundary condition ${\mathit{u}}_f^{i+1}={\mathit{u}}_{in}\left(t^{i+1}\right)$ on ${\Gamma }_{in}$, by introducing a lifting function ${\ell }_u$, similar to that performed for the monolithic approach in Section 4.3: we therefore obtain a homogenized fluid velocity ${\mathit{u}}_{0,h}^{i+1}$ such that ${\mathit{u}}_{0,h}^{i+1}=0$ on ${\Gamma }_{in}$. Afterwards, we define a new variable ${\mathit{z}}_{f,h}^{i+1}$:

${\mathit{z}}_{f,h}^{i+1}:={\mathit{u}}_{0,h}^{i+1}-D_t{\mathit{d}}_{f,h}^{i+1}={\mathit{u}}_{f,h}^{i+1}-{\ell }_u^{i+1}-D_t{\mathit{d}}_{f,h}^{i+1}$

With this change in variable, Equation (25) is equivalent to the homogeneous boundary condition for the new variable:

${\mathit{z}}_{f,h}^{i+1}=0\mathrm{on}{\Gamma }_{FSI},$

for which no imposition by means of Lagrange multiplier is needed. Therefore, during the offline phase of the scheme, at every iteration $i+1$, after we compute the homogenized velocity $u_{0,h}^{i+1}$, we compute the change of variable ${\mathit{z}}_{f,h}^{i+1}$. We then consider the following snapshots matrix:

${\mathcal{S}}_z=[{\mathit{z}}_{f,h}^1,\cdots ,{\mathit{z}}_{f,h}^{N_T}]\in {\mathbb{R}}^{{\mathcal{N}}_h^u\times N_T},$

where ${\mathcal{N}}_h^u=\dim V_h$. We then apply a POD to the snapshots matrix ${\mathcal{S}}_z$, and we retain the first $N_z$ POD modes ${\mathsf{\Phi }}_1^z,\cdots ,{\mathsf{\Phi }}_{N_z}^z$. We therefore have the reduced space:

$V^N:=\mathrm{span}{\left\{{\mathsf{\Phi }}_k^z\right\}}_{k=1}^{N_z},$

and now it is clear that, since every ${\mathsf{\Phi }}_k^z$ satisfies the condition ${\mathsf{\Phi }}_k^z=0$ on ${\Gamma }_{FSI}$, then also every element of $V^N$ satisfies the same condition.

6.3.2. Harmonic Extension of the Fluid Displacement

In order to generate a reduced basis for the fluid displacement ${\mathit{d}}_f$, we pursue the idea presented in [12]. Therefore, we start by generating the snapshots matrix related to the solid displacement:

${\mathcal{S}}_{d_s}=[{\mathit{d}}_{s,h}^1,\cdots ,{\mathit{d}}_{s,h}^{N_T}]\in {\mathbb{R}}^{{\mathcal{N}}_h^{d_s}\times N_T},$

where ${\mathcal{N}}_h^{d_s}=\dim E_h^s$, and again, the underline notation denotes the vector of the FE degrees of freedom corresponding to each solution of the solid displacement. We then apply a POD to the snapshots matrix and retain the first $N_{d_s}$ POD modes ${\mathsf{\Phi }}_1^{d_s},\cdots ,{\mathsf{\Phi }}_{N_{d_s}}^{d_s}$, thus defining the reduced space for the solid problem:

$E_N^s:=\mathrm{span}{\left\{{\mathsf{\Phi }}_k^{d_s}\right\}}_{k=1}^{N_{d_s}}.$

We then employ a harmonic extension of each one of the reduced basis ${\mathsf{\Phi }}_k^{d_s}$ to the fluid domain, thus obtaining the functions ${\mathsf{\Phi }}_k^{d_f}$ such that

$\left\{\begin{array}{c}-\Delta {\mathsf{\Phi }}_k^{d_f}=0\mathrm{in}{\Omega }_f,\hfill \\ {\mathsf{\Phi }}_k^{d_f}={\mathsf{\Phi }}_k^{d_s}\mathrm{on}{\Gamma }_{FSI}.\hfill \end{array}\right.$

We can then define the reduced space for the fluid displacement:

$E_N^f:=\mathrm{span}{\left\{{\mathsf{\Phi }}_k^{d_f}\right\}}_{k=1}^{N_{d_s}}.$

The reason for defining the basis functions for ${\mathit{d}}_f$ in such a way instead of employing a standard POD on the set of snapshots for the fluid displacement computed in the offline phase lies in the fact that we want to avoid the introduction of another Lagrange multiplier to impose the non-homogeneous boundary condition (18) With our method, we avoid solving the reduced system related to Equation (18): indeed, instead of solving an harmonic extension problem at every time-step in the online phase, we solve once and for all $N_{d_s}$ harmonic extension problems in the expensive offline phase. Then, during the online phase, the reduced fluid displacement is computed just as a linear combination of the basis ${\mathsf{\Phi }}_k^{d_f}$, with coefficients that are the coefficients of the reduced solid displacement at the previous time-step. We see the final formulation of the online phase of the algorithm in the next section.

6.4. Online Computational Phase

We are now ready to present the online formulation of the partitioned procedure. For every $i=0,\dots ,N_T$, we introduce the reduced functions $z_{f,N}^{i+1}$, $p_{f,N}^{0,i+1}$, $d_{s,N}^{i+1}$ of the following form:

(29)$\begin{array}{c}{z}_{f,N}^{i+1}=\sum _{k=1}^{N_{z_f}}{\underline{\mathit{z}}}_k^{i+1}{\mathsf{\Phi }}_k^{z_f},\end{array}$

(30)$\begin{array}{c}{p}_{f,N}^{0,i+1}=\sum _{k=1}^{N_p}{\underline{p}}_k^{0,i+1}{\mathsf{\Phi }}_k^p,\end{array}$

(31)$\begin{array}{c}{d}_{s,N}^{i+1}=\sum _{k=1}^{N_{d_s}}{\underline{\mathit{d}}}_k^{i+1}{\mathsf{\Phi }}_k^{d_s}.\end{array}$

Then, the online phase of the partitioned procedure reads as follows:

Mesh displacement:

let ${\mathit{d}}_{f,N}^{i+1}$ be defined by the reduced solid displacement at the previous time-step:

(32)${\mathit{d}}_{f,N}^{i+1}=\sum _{k=1}^{N_{d_s}}{\underline{\mathit{d}}}_k^i{\mathsf{\Phi }}_k^{d_f};$

Fluid explicit step (with change of variable):

find ${\mathit{z}}_{f,N}^{i+1}\in V_N$ such that $\forall {\mathit{v}}_{f,N}\in V_N$:

(33)$\begin{array}{cc}& {\rho }_f\Delta T^{-1}{(J\left({\mathit{z}}_{f,N}^{i+1}-{\mathit{u}}_{f,N}^i\right),{\mathit{v}}_{f,N})}_{{\Omega }_f}+{\rho }_f{(J(\nabla ({\mathit{z}}_{f,N}^{i+1}+D_t{\mathit{d}}_{f,N}^{i+1}+{\ell }_{u,N}^{i+1})F^{-1}{\mathit{z}}_{f,N}^{i+1}),{\mathit{v}}_{f,N})}_{{\Omega }_f}\hfill \\ & +{\rho }_f{\nu }_f{(J\epsilon \left({\mathit{z}}_{f,N}^{i+1}\right)F^{-T},\nabla {\mathit{v}}_{f,N})}_{{\Omega }_f}+{(J{F}^{-T}\nabla p_{f,N}^i,{\mathit{v}}_{f,N})}_{{\Omega }_f}+{\rho }_f{(J(\nabla {\mathit{z}}_{f,N}^{i+1}{F}^{-1}{\ell }_{f,N}^{i+1}),{\mathit{v}}_{f,N})}_{{\Omega }_f}=\hfill \\ & -{\rho }_f\Delta T^{-1}{(J\left(D_t{\mathit{d}}_{f,N}^{i+1}\right),{\mathit{v}}_{f,N})}_{{\Omega }_f}-{\rho }_f\Delta T^{-1}{(J{\ell }_{u,N}^{i+1},{\mathit{v}}_{f,N})}_{{\Omega }_f}-{\rho }_f{\nu }_f{(J\epsilon \left(D_t{\mathit{d}}_{f,N}^{i+1}\right)F^{-T},\nabla {\mathit{v}}_{f,N})}_{{\Omega }_f}\hfill \\ & -{\rho }_f{\nu }_f{(J\epsilon \left({\ell }_{f,N}^{i+1}\right)F^{-T},\nabla {\mathit{v}}_{f,N})}_{{\Omega }_f}-{\rho }_f{(J(\nabla (D_t{\mathit{d}}_{f,N}^{i+1}+{\ell }_{u,N}^{i+1})F^{-1}{\ell }_{f,N}^{i+1}),{\mathit{v}}_{f,N})}_{{\Omega }_f}+{(J{\mathit{b}}_N^f,{\mathit{v}}_{f,N})}_{{\Omega }_f},\hfill \end{array}$

where again ${\mathit{b}}_N^f$ is the projection of the fluid volume external force on the space $V_N$. We then restore the reduced fluid velocity: ${\mathit{u}}_{f,N}^{i+1}={\mathit{z}}_{f,N}^{i+1}+D_t{\mathit{d}}_{f,N}^{i+1}+{\ell }_{u,N}^{i+1}$. Here ${\ell }_{u,N}^{i+1}$ is the projection of the FE lifting function ${\ell }_{u,h}^{i+1}$ over the reduced basis space $V_N$.

Implicit step:

for any $j=0,\dots$ until convergence:

• Fluid projection substep: find $p_{f,N}^{0,i+1,j+1}\in Q_N^0$ such that $\forall q_{f,N}\in Q_N^0$:

  $\begin{array}{cc}\hfill & -\frac{{\rho }_f}{\Delta T}{(\mathrm{div}\left(J{F}^{-1}{\mathit{u}}_{f,N}^{i+1}\right),q_{f,N})}_{{\Omega }_f}-{\rho }_f{(\left(D_{tt}{\mathit{d}}_{s,N}^{i+1,j}\right),J{F}^{-T}{\mathit{n}}_f{q}_{f,N})}_{{\Gamma }_{FSI}}\hfill \\ \hfill & +{\alpha }_{ROB}{(p_{f,N}^{i+1,j},q_{f,N})}_{{\Gamma }_{FSI}}-{\alpha }_{ROB}{({\ell }_p^{i+1},q_{f,N})}_{{\Gamma }_{FSI}}\hfill \\ \hfill & -J{F}^{-T}{(\nabla {\ell }_p^{i+1},F^{-T}\nabla q_{f,N})}_{{\Omega }_f}={\alpha }_{ROB}{(p_{f,N}^{i+1,j+1},q_{f,N})}_{{\Gamma }_{FSI}}\hfill \\ \hfill & +{(J{F}^{-T}\nabla p_{f,N}^{i+1,j+1},F^{-T}\nabla q_{f,N})}_{{\Omega }_f};\hfill \end{array}$

  we then recover the reduced fluid pressure $p_{f,N}^{i+1,j+1}=p_{f,N}^{0,i+1,j+1}+{\ell }_p^{i+1}$.

• Structure projection substep: find ${\mathit{d}}_{s,N}^{i+1,j+1}\in E_N^s$ such that $\forall {\mathit{e}}_{s,N}\in E_N^s$:

  (34)${\rho }_s{(D_{tt}{\mathit{d}}_{s,N}^{i+1,j+1},{\mathit{e}}_{s,N})}_{{\Omega }_s}+{(P\left({\mathit{d}}_{s,N}^{i+1,j+1}\right),\nabla {\mathit{e}}_{s,N})}_{{\Omega }_s}=-{(J{\sigma }_f({\mathit{u}}_{f,N}^{i+1},p_{f,N}^{i+1,j+1})F^{-T}{\mathit{n}}_f,{\mathit{e}}_{s,N})}_{{\Gamma }_{FSI}}+{({\mathit{b}}_N^s,{\mathit{e}}_{s,N})}_{{\Omega }_s},$

  where ${\mathit{b}}_N^s$ is the projection of the solid volume external force on the space $E_N^s$.

Additionally, at the reduced order level, the implicit steps are iterated between one another until the stopping criteria Equation (28) is satisfied.

7. Results

We now present some results obtained by implementing the partitioned algorithm previously described for the same FSI test case, namely the toy problem inspired by the Turek–Hron benchmark test case: we remark here also that there is a difference with respect to the original test case presented by Turek and Hron, since in [31,32], the authors consider the Green strain tensor for the solid, whereas we restored a linear elasticity problem for the structure. The physical constants describing the fluid and the solid properties are the same, and we therefore refer to the values of Table 1.

The time-step used for the segregated approach is $\Delta t={10}^{-3}$, for a total of $N_T={10}^4$ iterations. Additionally, in this case, we adopt a standard FE approach for the first 8000 iterations and we employ the partitioned reduced order model for the remaining 2000 iterations. All of the values are reported in Table 3.

Figure 14a show the rate of decay of the eigenvalues for the fluid pressure $p_f$, the solid displacement ${\mathit{d}}_s$, and the fluid velocity change of variable ${\mathit{z}}_f$. As we can see, in this case, the eigenvalues for the fluid variable ${\mathit{z}}_f$ and for the pressure $p_f$ show almost the same rate of decay, whereas the eigenvalues for the solid displacement decay much faster. In Figure 14b, we can see the energy retained by the first one hundred modes for ${\mathit{z}}_f$, $p_f$ and ${\mathit{d}}_s$. As we can see, the first modes for the fluid pressure are the most energetic ones: they retain almost $10\%$ more energy with respect to the first modes of the solid displacement; additionally, the first modes of the change of variable retain a larger amount of energy with respect to the energy retained by the first modes of the solid displacement.

Figure 15 represents the first three modes for the fluid change of variable ${\mathit{z}}_f$, and it is interesting to notice how the modes are zero not only on the inlet boundary but also along the FSI interface thanks to the implementation of the change of variable. Figure 16 shows the first three modes for the fluid pressure: as we can see, in this case, the reduced basis also show a highly oscillatory behavior; nevertheless, thanks to the choice of the Chorin–Temam projection scheme, we are now able to obtain stable approximations of the fluid pressure even without the supremizer enrichment. Finally, Figure 17 shows the first three modes for ${\mathit{d}}_s$ (left column); on the right column, we have portraied the corresponding first three basis functions for ${\mathit{d}}_f$, obtained with the harmonic extensions (and not with a POD on the mesh displacement snapshots!).

In Figure 18, we can see the deformation of the elastic beam at the last time-step of the simulation: in this case, we used 13 basis functions for the reduced order approach. We can see that there is some difference between the FE solution and the reduced solution, as the pointwise relative error tells us, with the highest approximation error again at the point where the bar bends downwards: we suspect that this is because, in order to speed up the computations in the online phase, we relaxed the tolerance for the implicit iterations from ${10}^{-8}$ (FE discretization) to ${10}^{-5}$. Future investigations will be performed on how to improve the approximation of the solid behavior. Figure 19 shows the fluid pressure that was obtained in the partitioned scheme by solving a Poisson problem: this clearly leads to numerical results that are slightly different from the ones in the monolithic approach, but this has to be expected, since we introduced a “fictious” Dirichlet boundary condition for the fluid pressure in order to guarantee uniqueness of the solution of the Poisson problem. As we can see, the reduced order fluid pressure accurately represents the FE snapshot even without the use of the supremizer enrichment technique. Finally, in Figure 20, we can see the behavior of the fluid velocity: additionally, in this case, we recognize the Karman vortices that develop after a while, and in this case, we can see that most of the error is localized in the regions where these vortices detach from the bar and start to propagate into the fluid domain.

In Figure 21, we depicted the behavior of the average relative approximation error (average with respect to time) of the various components of the problem, namely ${\mathit{u}}_f$, $p_f$, ${\mathit{d}}_f$ and ${\mathit{d}}_s$. As we can see, the fluid velocity has the best approximation error compared to the other components, and adding reduced basis functions does not seem to provide a significant improvement in the approximation quality; on the contrary, adding basis functions for $p_f$ does seem to affect the quality of the approximation. As we can see, from $N=19$, it looks like we just add noise to the online system. Finally, we can see that, as expected also from the previous representation of the pointwise approximation error (Figure 18, the solid displacement is the one for which we have the higher approximation error: again, this could depend on the relaxation of the convergence tolerance of the implicit steps in the online system, but future investigations will be carried out on how to improve the approximation of the solid behavior. Figure 22 represents the average approximation error of the solid stress at the interface: as we can see, we have a good approximation for the first 19 basis functions, and then again from $N=19$, we get the worst approximation error: again, we suspect this is because we add noise to the system and because, from $N=19$, the average approximation error of all of the components becomes worse.

As a final result, we present in Figure 23 the average number of implicit iterations for the partitioned online phase, plotted against the number of modes used. It is interesting to note that we experience an increase in the average number of iterations in the implicit step (and thus an increase in the computational cost of the online phase) as we increase the number of basis functions; we remember that these results were obtained for a chosen tolerance of $\epsilon ={10}^{-5}$ for the fixed point iterations (recall Equation (28)).

8. Discussion

The aim of the work was to provide the reader with an extensive overview about the combination of the RBM with the two main different approaches that are used to address a Fluid–Structure Interaction problem, namely the monolithic and the partitioned approach. We provided insights on the two different kind of algorithms that differ from one another in various aspects: starting from the way the finite element discretization of the original problem is carried out, we then have a significant difference also in the Proper Orthogonal Decomposition. Indeed, even if the underlying idea is to perform a compression of the set of snapshots by means of a POD, we observe that, in the monolithic POD, we work with huge snapshots matrices and with huge inner product matrices because we keep the monolithic structure of the whole approach; in the partitioned approach instead, the snapshots matrices are much smaller. This difference reflects the time that the two PODs take to perform, with the monolithic one being slower than the partitioned one (because we have bigger correlation matrices for which we solve the eigenvalue problem). In addition, the difference between the two approaches also reflects the formulation of the online system: we have a big block system for the monolithic RBM and a series of much smaller systems for the partitioned approach.

We applied the two different model order reduction approaches to a benchmark test case of interest, namely a test case inspired by the Turek–Hron FSI test case FSI2 in order to provide the reader with some additional numerical results that provide a better insight on the two procedures.

Figure 24 represents a comparison between the amount of energy retained by the first one hundred modes for the fluid velocity ${\mathit{u}}_f$ (modes obtained with the monolithic approach) and the first one hundred modes for the change of variable ${\mathit{z}}_f$ (modes obtained with a partitioned approach). It is interesting to see that the first mode for ${\mathit{z}}_f$ retains almost $10\%$ more energy with respect to the first mode for ${\mathit{u}}_f$, and the general trend is that the first modes for ${\mathit{z}}_f$ retain more energy, thus leading us to believe, at a first glance, that we need fewer reduced basis functions for the fluid momentum equation in the partitioned approach.

In Table 4, we summarize some important aspects of both the reduction procedures in order to highlight the differences of the two approaches. As we can see, the Newton method hwa used in the monolithic approach to solve the huge system and, in the partitioned approach, to solve the fluid explicit step; we used a tolerance for the max. norm of the residual between two consecutive iterations of the solver equal to $6·{10}^{-6}$. We remark that it is further possible, if necessary, to speed up the time required for solving the big online system in the monolithic RBM by using some preconditioners: plenty of results on the improvement of the performance of the RBM with preconditioners exist, we refer for example to [16]. For the implicit steps of the partitioned online system (pressure Poisson recovery and linear elasticity), we solve the two linear systems with the numpy.linalg solver for linear systems, which is based on a LAPACK routine that performs the LU factorization of the matrix on the left hand side. All of the numerical simulations were performed on a computer with $3.50$ GHz per CPU, and the mesh used in the two approaches is the same, as we can see from the mesh resolution reported in Table 2 and Table 3. For the numerical simulation of the offline phase of both approaches, we relied on multiphenics [55], which is a Python-based library that helps with the implementation of simulations on conformal meshes of multiphysic problemsM; for the numerical simulation of the online phase, we relied on rbnics [56], which is a Python implementation of several reduced order modelling techniques.

The monolithic algorithm, as we have seen, brings along an increase in the number of unknowns to be used in the coupled system: this is because, in order to impose the coupling conditions, we used two Lagrange multipliers, and this leads to an increase in the dimension of the algebraic system to be solved during the online phase. In addition to this, in order to obtain a stable approximation of the reduced order fluid pressure, we adopted a supremizer enrichment of the reduced fluid velocity space: this also leads to an increase of the number of basis functions to be used in the online phase, with a further increase of the dimension of the algebraic system. On the other hand, a monolithic approach is more stable and, therefore, allows for bigger time-steps in the numerical simulations: this turns out to be extremely useful, especially in the case of physical phenomena that take some time to develop, such as the Karman vortices in the Turek–Hron benchmark that we considered. We also make the following remark: for the monolithic approach, one could have used globally continuous spaces for the velocity and the displacement. This approach at the FE level is the one originally adopted by Turek and Hron (see [31]), and within the RBM, there are results present in the literature (see [11]). At the finite element level, globally continuous spaces for velocity and displacement result in many advantages; however, from the RBM point of view, there are some aspects that we feel are much more easily handled with the approach that we have proposed. First, having globally defined reduced basis functions does not automatically guarantee the balance of the stresses at the interface: this coupling condition would in any case require a weak imposition by means of Lagrange multipliers. A second important aspect is the implementation of the supremizer enrichment technique: if we were to use globally continuous spaces, the supremizer problem would be solved not only in the fluid domain (which is exactly where the supremizer is needed, for stabilization purposes) but also in the entire domain and, hence, in the solid subdomain, adding unnecessary DOFs to the supremizer enrichment problem. Finally, as it is remarked in [11], in the online system, it would not be possible to highlight the contributions corresponding to fluid and solid DOFs, as the unknowns of the online system would be coefficients of a modal expansion (global in both fluid and solid subdomains) and thus not related to a spatial location in the domain. Due to all of these details, we restored to a block formulation from the FE stage: the implementation of the code with this block structure is also much easier thanks to the use of multiphenics.

The partitioned algorithm is very useful, and it gives a lot of control on the systems to be solved during the online phase of the reduction method. In this work we chose a semi-implicit treatment of the coupling conditions: this is, in our opinion, the best choice for the test case considered. Indeed, an implicit treatment of the coupling conditions would be too expensive and would drastically increase the number of sub-iterations at each time-step, both in the offline phase and in the online phase. We also tried to adopt an explicit treatment of the coupling conditions, as suggested for example in [40]: this approach is unstable for the benchmark considered here, because, due to the slender shape of the domain, the added mass effect plays an important role and leads to an algorithm that diverges after a few time-steps. Therefore, the implicit coupling is the best tradeoff between stability and computational cost. The Chorin–Temam projection scheme has allowed us to work without the need of a supremizer enrichment technique: this gives good control in the number of basis functions to be used in the online phase of the method. In addition to this, the choice of a pressure Poisson formulation allows us to discard the so-called end of step velocity and to work just with the intermediate velocity, leading again to a decrease in the dimension of the online system. Finally, the harmonic extension of the fluid displacement basis functions gives the possibility to efficiently compute the mesh displacement in the online phase of the method without the need to solve an additional system. The drawback of a partitioned reduced order model, as we have seen in the Numerical Results section, is that the time-step required in order to have a stable algorithm is, in general, much smaller with respect to the time-step that can be used in a monolithic approach, thus resulting in a larger number of snapshots to be processed in the offline phase. In addition to this, the treatment of the boundary conditions is rather delicate with partitioned approaches and needs to be tailored to the problem at hand, whereas with a monolithic approach, the imposition of these conditions is simpler, either by incorporating them in the weak formulation or by using Lagrange multipliers.

In this work, we presented numerical results concerning the behavior of the FSI system up to a time of $0.9$ s; however, inspired by the many results present in the literature for the Turek–Hron benchmark test case FSI2, we expect that, given the periodic oscillatory behavior of the structure, for a longer period of observation, vortex shedding phenomena may occur also in our test case. For the Reynolds number considered here, we expect, for both the monolithic and the partitioned approaches, that a much larger number of fluid (pressure and velocity) reduced basis are needed in order to obtain a good approximation of the fluid behavior. If we were to consider a higher Reynolds number for the fluid, at the reduced order level, the supremizer enirchment technique may not be sufficient to obtain a stable approximation and other stabilization techniques may be required, such as SUPG (see for example [53,63,64]). If the Reynolds number increases significantly and the behavior of the fluid becomes turbulent, then a finite volume based approximation may be the right choice, and in this case, we see an advantage in the partitioned RBM with respect to the monolithic RBM: the partitioned approach indeed potentially allows us to use two different spatial discretization techniques for the fluid and for the solid problem; we refer the interested reader to [6,65,66,67] for some results on the implementation of finite volume discretization within the RBM.

In this work, we also considered a linearized strain tensor for the solid domain because we were interested in developing a procedure for the small deformations range; nonetheless, everything we said can be applied to a nonlinear solid problem as well. If the structure in the system undergoes large deformations, then an ALE formulation may not be the right formalism within which to study the behavior of the coupled system: indeed, it is known that, for large deformation, the ALE formalism may lead to some complications. In this case, we expect a Cut Finite Element approach to be more suited; see for example [68]: the investigation of the performance of a Cut Finite Element based RBM for FSI within a monolithic approach is currently under investigation (we refer to [8,69] for some preliminary results on computational fluid dynamics problems); to the best of our knowledge, there are no results in the literature concerning the application of Cut Finite Element discretization to the RBM within a partitioned approach instead.

We conclude the discussion with some reflections and perspectives for future works and investigations. For the monolithic reduced order model, even though we introduce a few new variables in the problem formulation (the Lagrange multipliers), we can still make the entire approach computationally feasible: for example, a suggestion could be to perform a rather sparse Proper Orthogonal Decomposition on the first snapshots (when the important physical phenomenon, such as the Karman vortices in our test case, is not yet developed) and then to refine the sampling. We are aware that this requires some sort of “a priori” knowledge of the physical phenomenon that we are simulating: for the benchmark test case considered, plenty of numerical results already exist that can provide an insight on when to refine the sampling procedure. For other applications, this opens up the possibility of bridging a monolithic reduction procedure with some machine learning algorithm that can be used to investigate, in an efficient way, the overall behavior of the solution to be approximated.

As long as the partitioned approach is concerned, its application is very well indicated for those problems that do not require long time simulations and for industrial applications since the idea of a segregated algorithm is to combine already existing state-of-the-art softwares for computational fluid dynamics and computational solid mechanics. Nonetheless, until alternative stopping criteria and/or alternative treatments of the coupling conditions are further investigated, their application in long time simulations results in an increase in the computational time during the online phase.

9. Conclusions

In this work, we presented an overview on two possible reduced order models for FSI problems that are based on two different approaches: a monolithic or a partitioned approach. We provided the details of the implementation of the two reduction procedures, both at the FE level and at the reduced order level, analysing the different aspects of the two algorithms, such as the change in variable for the fluid velocity in the partitioned procedure, the creation of mesh displacement basis functions thanks to an harmonic problem, the block structure of the matrices in the monolithic POD, and the different treatment of the coupling conditions at the fluid–structure interface. Finally, we implemented the aforementioned algorithms for a toy problem of interest, which was inspired by the Turek–Hron benchmark test case FSI2. We provided numerical results for the monolithic and for the partitioned RBM, showing, among other things, the behavior of the average approximation error as a function of the number of modes used and the average error between the solid stress at the FSI interface with the FE discretization and with the reduced order discretization. We have seen how the RBM can be modified and adapted in order to be tailored to a monolithic or to a partitioned approach, and the work presented represents an interesting overview, expecially for what concerns segregated reduced order procedures, for which not many results in the literature exist, to the best of our knowledge.

Author Contributions

Conceptualization, methodology, resources, and supervision: G.R. and F.B. Conceptualization, methodology, software, validation, formal analysis, investigation, data curation, and writing—original draft preparation: M.N.; writing—review and editing: G.R., F.B., and M.N.; project administration: G.R.; funding acquisition: G.R. and M.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We acknowledge the support by European Union Funding for Research and Innovation—Horizon 2020 Program—in the framework of European Research Council Executive Agency: Consolidator grant H2020 ERC CoG 2015 AROMA-CFD project 681447 “Advanced Reduced Order Methods with Applications in Computational Fluid Dynamics” (Gianluigi Rozza). We also acknowledge the INDAM-GNCS project “Tecniche Numeriche Avanzate per Applicazioni Industriali”. Numerical simulations were obtained with the extension multiphenics of FEniCS [2,34] for the high fidelity part and from RBniCS [1] for the reduced order part. We acknowledge the developers and contributors of each of the aforementioned libraries. Nonino also acknowledges the support of the Austrian Science Fund (FWF) project F65 “Taming complexity in Partial Differential Systems” and the Austrian Science Fund (FWF) project P 33477.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Fluid–structure interaction domains. Time-dependent configuration (left): fluid domain Ωf(t) (blue) and solid deformed domain Ωs(t) (red). Reference configuration (right): the inlet boundary Γ^in (magenta), the wall boundaries for the fluid Γ^walls (light blue), and the Neumann (outlet) boundary Γ^Nf (orange). The fluid arbitrary reference configuration Ω^f is depicted in blue, the solid reference configuration Ω^s is depicted in red, and Γ^Ds is the solid Dirichlet boundary. The fluid–structure interface in the reference configuration Γ^FSI is highlighted in green.

Figure 2. Reference configuration of the benchmark test case. The solid is depicted in red, while the fluid domain is in blue.

Figure 3. The outcomes of the monolithic POD: eigenvalue decay (a) and retained energy (b) for the first 100 modes.

Figure 4. The first POD modes for the fluid velocity uf and for the supremizer s (magnitude).

Figure 5. The first three POD modes for the fluid pressure pf.

Figure 6. The first three POD modes for the fluid displacement df (magnitude).

Figure 7. The first three POD modes for the solid displacement ds (magnitude).

Figure 8. The deformation of the bar at time t=0.9 s: the FE solution (top left) and the reduced order solution (top right). Bottom: approximation error |ds,h−ds,N|, represented over the solid reference configuration (undeformed state). Nds=21 basis functions were used for the solid displacement. The deformation was magnified by a factor 5 for visualization purposes.

Figure 9. Fluid pressure at time t=0.9 s: the FE solution (top left) and the reduced order solution (top right). Bottom: approximation error |pf,h−pf,N|. Np=21 basis functions were used for the fluid pressure, with the supremizer enrichment technique.

Figure 10. Fluid pressure approximation without the implementation of the supremizer enrichment: solution before the code diverges.

Figure 11. Fluid velocity at time t=0.9 s: the FE solution (top left) and the reduced order solution (top right). Bottom: approximation error |uf,h−uf,N|. Nu=42 basis functions were used for the fluid velocity.

Figure 12. Average relative error as a function of the number of basis functions used in the online system.

Figure 13. Average approximation error for the solid stress at the FSI interface.

Figure 14. The outcomes of the partitioned POD: eigenvalue decay (a) and energy retained (b) for the first 100 modes.

Figure 15. The first three POD modes for the fluid velocity change in variable zf (magnitude): notice how, for the partitioned approach, the magnitude of the modes is zero not only on the inlet boundary (thanks to the lifting function) but also on the fluid–structure interface (thanks to the implementation of the change of variable).

Figure 16. The first three POD modes for the Poisson recovery of the fluid pressure pf.

Figure 17. The first three POD modes for the solid displacement ds (left column) and the corresponding mesh displacement modes (right column) obtained with an harmonic extension of the basis functions on the left column.

Figure 18. Deformation of the structure: FE solution (top left), reduced order solution (top right), and local approximation error (bottom). The approximation was obtained with Nds=13 basis functions. The deformation was magnified by a factor 10 for visualization purposes.

Figure 19. Pressure Poisson recovery: FE solution (top left), reduced order solution (top right), and local approximation error (bottom). The approximation was obtained with Npf=13 basis functions.

Figure 20. Fluid velocity: FE solution (top left), reduced order solution (top right), and local approximation error (bottom). The approximation was obtained with Nuf=13 basis functions.

Figure 21. Average relative approximation error as a function of the number N of modes used in the online phase.

Figure 22. Average approximation error for the solid stress at the interface.

Figure 23. Average number of iterations needed in the implicit step in order to reach convergence, as a function of the number N of reduced basis for uf, pf, ds.

Figure 24. A comparison between the energy retained by the first modes for zf and for uf.

Table 1

Values of the physical constants of the fluid and of the solid.

Parameter	Value	Parameter	Value
Fluid density ${\rho }_f$$\left[{10}^3\frac{\mathrm{kg}}{{\mathrm{m}}^3}\right]$	1	Solid density ${\rho }_s$$\left[{10}^3\frac{\mathrm{k}\mathrm{g}}{{\mathrm{m}}^3}\right]$	1
Fluid kinematic viscosity ${\nu }_f$$\left[{10}^{-3}\frac{{\mathrm{m}}^2}{\mathrm{s}}\right]$	1	2nd Lamé constant (solid shear modulus) ${\mu }_s$$\left[{10}^6\frac{\mathrm{k}\mathrm{g}}{\mathrm{m}{\mathrm{s}}^2}\right]$	$0.5$
Mean inflow velocity $\overline{U}$ $\left[\frac{\mathrm{m}}{\mathrm{s}}\right]$	1	1st Lamé constant ${\lambda }_s$	$0.4$
Fluid external force ${\mathit{b}}_f$	$(0,0)$	Solid external force ${\mathit{b}}_s$	$(0,0)$

Table 2

Values of the parameters in the time discretization and in the spatial discretization.

Time Discretization Parameters	Value
Timestep $\Delta T$	$0.01$ s
Total number of iterations $N_T$	${10}^3$
$\gamma$	$0.25$
$\beta$	$0.5$
Space Discretization Parameters	Value
FE velocity order	2
FE pressure order	2
FE displacement order (${\mathit{d}}_f$ and ${\mathit{d}}_s$)	2
FE multiplier order (${\mathit{\lambda }}_u$ and ${\mathit{\lambda }}_d$)	1
mesh resolution using mshr mesh generator	128

Table 3

Values of the parameters in the time discretization and in the spatial discretization.

Time Discretization Parameters	Value
$\Delta T$	$0.0001$ s
total number of iterations $N_T$	${10}^4$
Space Discretization Parameters	Value
FE velocity order	2
FE pressure order	1
FE displacement order	2
mesh resolution using mshr generator	128
tolerance $\epsilon$ for the implicit iterations	${10}^{-5}$

Table 4

A comparison between the two approaches.

Monolithic Approach	Value
$\Delta T$	$0.01$ s
number of time iterations i of the RB solver	200
solver for the system	Newton method
average iterations of Newton method	4
absolute tolerance for Newton method	$\left|\right|·{\left|\right|}_{\infty }<6·{10}^{-6}$
computational time to solve the online system for one time iteration	$224.9$ s
Partitioned Approach	Value
$\Delta T$	$0.001$ s
number of time iterations i of the RB solver	2000
solver for explicit fluid step	Newton method
average number of iterations of Newton method	2
absolute tolerance for Newton method	$\left|\right|·{\left|\right|}_{\infty }<6·{10}^{-6}$
computational time to solve the online systems (explicit + implicit) for one time iteration	$162.68$ s