Projection Based Semi-Implicit Partitioned Reduced Basis Method for Fluid-Structure Interaction Problems

In this manuscript a POD-Galerkin based Reduced Order Model for unsteady Fluid-Structure Interaction problems is presented. The model is based on a partitioned algorithm, with semi-implicit treatment of the coupling conditions. A Chorin–Temam projection scheme is applied to the incompressible Navier–Stokes problem, and a Robin coupling condition is used for the coupling between the fluid and the solid. The coupled problem is based on an Arbitrary Lagrangian Eulerian formulation, and the Proper Orthogonal Decomposition procedure is used for the generation of the reduced basis. We extend existing works on a segregated Reduced Order Model for Fluid-Structure Interaction to unsteady problems that couple an incompressible, Newtonian fluid with a linear elastic solid, in two spatial dimensions. We consider three test cases to assess the overall capabilities of the method: an unsteady, non-parametrized problem, a problem that presents a geometrical parametrization of the solid domain, and finally, a problem where a parametrization of the solid’s shear modulus is taken into account.

The complex nature of these problems is reflected not only by their theoretical treatment, but also by the way they are treated numerically. There are two approaches that can be adopted in order to address a FSI problem: the first approach consists of a monolithic procedure [17][18][19], whereas the second approach consist of a partitioned (segregated), procedure [20][21][22].
In a monolithic algorithm the fluid and the solid problem are solved simultaneously: this results in algorithms that show good stability properties, with respect to time, independently of the density of the fluid and the solid, and independently of the geometrical properties of the physical domain; indeed the monolithic algorithm does not suffer from the so called added mass effect (see [23] for an analytical study of this phenomenon), which is very well known in the FSI community, and is responsible for the numerical instabilities in the design of partitioned algorithms. Stability in time is highly desirable in the framework of unsteady problems, especially if we wish to use large time-steps in the simulations; the main drawback of these monolithic algorithms is given by the fact that they deeply rely on the availability of legacy softwares that can be used to solve both the fluid problem and the solid problem: in this sense, monolithic algorithms are less flexible and more tailored to the particular problem at hand. In the literature there are many examples of works that are based on a monolithic approach: in [24] the author focuses on a monolithic algorithm to address a coupled problem, written within the Arbitrary Lagrangian Eulerian (ALE) formalism, which models the interaction between the blood flow and the arterial walls; another example of FSI problems related to the blood flow-arterial interaction can be found in [25]. In [26] the authors propose different preconditioners, to be used in a Newton-Krylov method for the nonlinear problem arising from solving in a monolithic fashion a coupled problem. For the reader interested in a general introduction to monolithic approaches to FSI problems, we refer to [27].
As an alternative to monolithic approaches, one can think of adopting a partitioned procedure; indeed, existing simulation tools for fluid dynamics and for structural dynamics are well developed and are used on a daily basis in industrial applications. It is therefore natural to try to combine these computational tools, to address coupled problems: this is exactly the rationale behind a partitioned algorithm. In a partitioned procedure, we solve separately the fluid and the solid problem, and then we couple the two physics with some iterative procedure, see for example [28]. Partitioned approaches are very flexible, as they allow to design the procedure in different ways, according to the problem under consideration. In [29], the authors propose a segregated algorithm to solve a FSI problem, where the coupling of the two physics at the fluid-structure interface is taken care of through a constrained optimization problem. In [30,31] the authors consider the problem of coupling an incompressible fluid with a thin structure; in [30] the authors propose a Robin-Neumann type of coupling at the fluid-structure interface, whereas in [31] the authors propose and explain different couplings techniques at the fluid-structure interface, within an explicit coupling setting. On the contrary, in [32] the authors deal with a problem that has a strong added mass effect, which is typically the case for the blood in the vessels: here, an implicit coupling is the way to control the stability issues due to the added mass effect. Nevertheless, it is clear that a fully implicit treatment of the coupling conditions leads to prohibitive computational costs; for this reason, in [32] the authors propose a semi-implicit coupling technique, which is the approach that will be adopted in this manuscript.
Addressing a coupled problem by means of a partitioned procedure is advantageous in terms of computational efficiency, also from the Reduced Order Model (ROM) point of view: indeed, in the online phase of the Reduced Basis Method [33][34][35][36][37][38][39], we have to solve, separately, smaller systems. Moreover, with some minor changes such as change of variables and appropriate choices for the couplings, it is possible to further reduce the dimension of the online systems, as we will see in the following. In the model order reduction framework, there is also a fair amount of work that is being carried out and which focuses on ROM-ROM and ROM-FOM coupling, see for example [40][41][42]. All these works represent an extremely interesting approach from which many FSI applications of interest could benefit; for this reason, the authors believe that this direction represents a future line of work within partitioned algorithms.
In this manuscript we design a segregated procedure, combined with a Reduced Order Model based on a Proper Orthogonal Decomposition. The goal is to extend the work done in [20], moving to the treatment of a two dimensional structure within an Arbitrary Lagrangian Eulerian formalism, and the work done in [43,44], adapting the computation and the treatment of the Robin coupling condition, also to the case of a thick, two dimensional structure. The present manuscript represents also an extension of the work done in [45], where the problem under consideration was only unsteady, but no geometric or physical parametrization has been considered.
This manuscript is structured as follows: in Sect. 2 we briefly introduce the Arbitrary Lagrangian Formulation, and we set the notation that will be used throughout the manuscript. In Sect. 3 we introduce the first test case, namely a time dependent, non parametrized FSI problem that models the interaction of a fluid with a thick, two dimensional, structure; in Sect. 3.2 we introduce the partitioned procedure at the high order level. In Sect. 3.3 we derive the partitioned procedure at the reduced order level, and in Sect. 3.4 we present the numerical results. In Sect. 4 we consider the same problem of interest, with the addition of a shape parametrization: in Sect. 4.1 we present the ALE formalism in the presence of a geometrical parametrization of the domain; in Sect. 4.2 we give the strong formulation of the problem of interest, and in Sect. 4.3 we describe the algorithm at the high order level. In Sect. 4.4 we introduce the reduced order model, and then we present some numerical results in Sect. 4.5, for the geometrical parametrization only. Then, in Sect. 4.6 we show some numerical results also in the presence of a physical parameter. Conclusions and considerations on future possible lines of work are presented in Sect. 5.

Configurations, Definitions and Notation
In this section we are going to introduce briefly the Arbitrary Lagrangian Eulerian (ALE) formalism, in order to set the notation that will be used throughout the rest of this manuscript.
In FSI problems the fluid domain is a moving domain (except for those situations in which the displacement of the solid is very small, and thus the whole physical domain can be considered as fixed). In solid mechanics, on the other hand, it is common to deal with deforming domains, and the deformation itself is the unknown of the problem; for fluid dynamics instead one usually considers fixed domains. This different point of view is the motivation behind a formalism, very known and widely used in the community, which is called the Arbitray Lagrangian Eulerian formulation [27,[46][47][48].
Let (t) ⊂ R 2 be the physical domain over which the FSI problem is formulated, where f (t) ⊂ R 2 and s (t) ⊂ R 2 are the fluid and the solid domain at time t, respectively; we assume that the two domains do not overlap, i.e. f (t)∩ s (t) = ∅, and finally, the fluid-structure interface F SI (t) is defined as To describe the behavior of a solid it is common practice to use the so called Lagrangian formalism: all the quantities and the conservation laws are formulated on the reference configurationˆ s = s (t = 0). On the contrary, when describing the behavior of a fluid, the Eulerian formalism is used instead: all the quantities and the conservation laws are formulated on the configuration f (t) at the current time t. In order to be able to describe both the fluid and the solid, a mixed formulation (the ALE formulation indeed) is used: the underlying idea is that of pulling back the fluid equations to an arbitrary time-independent configurationˆ f : one possible choice forˆ f isˆ f = f (t = 0), the domain at initial time.
In Fig. 1 we can see an example of a reference configuration and the configuration of the domain at the current time t. Let us see in the next paragraph how to introduce the ALE formalism; for a more detailed discussion about different approaches to describe coupled systems we refer to [27,49].
Let [0, T ] be a time interval, and letˆ f be a reference configuration for the fluid.

Definition 1
The ALE mapping A f (t), for every t ∈ [0, T ] is defined as follows: There are different possibilities for the definition of the mesh displacement: in this manuscript, we decide to defined f as an harmonic extension of the solid displacementd s on the whole fluid domainˆ f : and homogeneous Dirichlet boundary conditions on the remaining portion of the boundary. Hereˆ F SI is the fluid-structure interface in the reference configuration. A great deal of attention has to be paid to the definition of the mesh displacement, as different choices ford f lead to different degrees of regularity: if we lose regularity due to the mesh displacement, we consequently lose regularity at the FSI interface, which is exactly where the coupling between the two physics takes place. It is beyond the scope of this work to discuss the regularity of different definitions of the mesh displacement; nonetheless we refer the interested reader to Chapter 5.3.5 of [27]. Remark 1d f represents the displacement of the grid points, therefore it is not a quantity with a real physical meaning, but rather a geometrical quantity that describes the deformation of the mesh, according to the deformation of the physical domain. It is also important to underline that ∂ td f =û f : in fact, whileû f represents the velocity of the fluid, ∂ td f is again a geometrical quantity, that can be interpreted as the velocity with which the mesh moves.
Let us now define the gradient F of the ALE map and its determinant J , respectively: With these quantities we are ready to present the strong formulation of the FSI problem of interest, within an ALE formalism.

First Test Case: Time Dependent FSI Problem
We present the first FSI problem of interest: a time-dependent, non parametrized, nonlinear multiphysics test case. The goal is to simulate the behavior of an incompressible, Newtonian fluid interacting with a linear elastic solid, in the time interval [0, T ]; Fig. 2 shows the physical domain in its reference configuration.

Strong Formulation
The coupled FSI problem, formulated over the original configuration, reads as follows: find In system (1), thed iv denotes the fact that the divergence is computed with respect tox, the space variable in the reference configuration. ρ f and ρ s are the fluid and the solid density, while σ f is the fluid Cauchy stress tensor for an incompressible Newtonian fluid, andP is the solid first Piola-Kirchoff tensor (compressible, linear elastic solid). System (1) is then completed by some initial conditions (we assume the system to be at rest at the starting time of the simulation), and by the following boundary conditions: where p in (t) is a time-dependent pressure pulse, n is the outward unit normal to the boundary being considered, in and out represent the inlet and outlet boundaries depicted in Fig. 2, andˆ D s is the portion of the leaflets' boundary that is attached to the top and bottom walls top and bott . It only remains to state the coupling conditions that take place at the fluid-structure interface: where d s is the solid deformation in the current configuration s (t), F s :=∇d s + I, I the 2 × 2 identity matrix, J s := detF s , n f and n s are the unit normals to the FSI interface F SI (t) outgoing the fluid and the solid domain, respectively. In system (2) we can interpret the first condition as a kinematic coupling, which requires the continuity of the velocities at the interface (the fluid sticks to the moving boundary), whereas the second equation is a dynamic coupling corresponding to an action-reaction principle, which is simply stating that the two stresses have to balance out at the interface.
With the formalism introduced in Sect. 2, we are able to perform a pull-back of the fluid equations onto the fluid reference configurationˆ f ; the FSI problem on the reference configu-rationˆ :=ˆ f ∪ˆ s now reads: for every t ∈ [0, T ], find the fluid velocityû f (t) :ˆ f → R 2 , the fluid pressurep f (t) :ˆ f → R, the fluid displacementd f (t) :ˆ f → R 2 and the solid deformationd s (t) :ˆ s → R 2 such that The fluid tensorσ f is the representation in the reference configuration of the Cauchy stress tensor:σ where λ s and μ s are the first and second Lamé constant of the solid, respectively (μ s is also referred to as shear modulus). System (3) is completed by the same initial conditions, by the boundary conditions (4) and by the following coupling conditions: In the previous equations the vectorn represents the normal vector to the inlet (or outlet) boundary in the reference configuration, whereasn f andn s are the unit normals to the FSI interfaceˆ F SI , outgoing the fluid and the solid domain, respectively. In system (5), the first equation is a geometric condition, which states that the fluid and the solid domain do not overlap. The second condition is the kinematic condition, expressed in the reference configuration, and similarly the third condition is the dynamic condition in the reference configuration.

Remark 2
The gradient and the divergence in Eq. (3) are computed with respect to the spatial coordinates in the reference configuration, namelyx. Nevertheless, from now on, since everything will be formulated and computed on the reference configuration, in order to ease the exposition, we will drop theˆnotation.

Offline Computational Phase
We are now going to describe the offline phase of the partitioned procedure that we use to solve the FSI problem of this section. The algorithm is based on a Chorin-Temam projection scheme for the incompressible Navier-Stokes equations [50,51], and we choose to treat the coupling conditions (5) in a semi-implicit way (see also [20,32,52]). We first apply a time stepping procedure to design the algorithm, and then we show the space discretization of the whole procedure.

High Fidelity Semi-implicit Scheme
We present the offline phase of the partitioned procedure: we use an operator splitting approach, based on a Chorin-Temam projection scheme with pressure Poisson formulation. Let T be a time-step: we discretize the time interval [0, T ] with an equispaced sampling {t 0 , . . . , t N T }, where t i = i T , for i = 0, . . . , N T and N T = T T . We discretize the partial derivative of a function f with a first backward difference BDF1: where f i+1 = f (t i+1 ). Hereafter, we will make use of the BDF1 time discretization for both the fluid and the solid problem. We consider the following semi-implicit time discretization of (3): for i = 0, . . . , N T − 1

Extrapolation of the Mesh Displacement d f
find d i+1 f : f → R 2 such that: In this step we are imposing the first of our three coupling conditions, namely the continuity of the displacements at the fluid-structure interface. This condition is imposed in an explicit way, with respect to time, because we are taking into account the solid displacement d i s at the previous time iteration i, and the mesh displacement unknown d i+1 f at the current time iteration i + 1. The choice of imposing the continuity of the displacements in an explicit way is inspired by many works present in the literature of FSI, see for example [20].

Remark 3
Here we treat the first coupling condition in an explicit way: this approach is less strong than a monolithic one, where this coupling condition would be imposed weakly for the fluid and solid displacement at the same timestep t i+1 (see for example [45]). Nevertheless this choice will allow us to build the fluid displacement d f in a cheap way in the online phase.

Fluid Explicit
Step with ε(u i+1 Here, we are imposing the dynamic condition (continuity of the velocities at the FSI interface), again in an explicit fashion with respect to time, since now the fluid displacement d i+1 f is already known.

Implicit
Step subject to the boundary conditions: • Structure projection substep find d i+1 s : s → R 2 such that: subject to the boundary condition d i+1

Remark 4
We remark that the time-stepping schemes for the fluid problem and for the solid problem are implicit. The denomination "semi-implicit" comes from the fact that the coupling conditions are treated differently. Indeed in system (7), the geometrical coupling condition (the second equation) is treated explicitly; on the other hand, the coupling on the fluid and solid velocity (second equation in system (8)), as well as the coupling of the stresses at the fluid-structure interface (second equation in system (10)), are treated implicitly.
The implicit step couples pressure stresses to the structure, and it is iterated until convergence is reached. We would like to remark that, throughout this manuscript, the BDF1 time stepping scheme is used also for the solid problem: there are other alternatives, such as a Newmark-Beta or HHT-alpha methods, that represent a standard choice in solid mechanics, as they have some desirable properties regarding stability and dissipation. In our work we did not encounter any problem with the stability in time of the algorithm, and therefore chose to use a BDF1 scheme for its easiness of implementation.

Remark 5
In the implicit step (8) we have chosen a pressure Poisson formulation; an alternative is to use a Darcy formulation, which is defined as follows: find p i+1 Throughout this manuscript we choose to employ a Poisson formulation, for the sake of a more efficient reduced order model, since the Darcy formulation requires the introduction of an additional unknownũ f , which translates in a larger system, comprised of both velocity and pressure, at the implicit step.
Let us now have a look at Eqs. (7)-(10): the fluid problem is solved using Dirichlet boundary conditions (the displacement computed at the previous timestep), and the solid problem is solved using Neumann boundary conditions (the fluid normal stress just computed). However, as it is mentioned in [44] and references therein, these kind of partitioned schemes (Dirichlet-Neumann couplings) usually require a large amount of sub-iterations of the implicit step, before a convergence between fluid and solid problem is reached, especially in those situations where the added mass effect is particularly heavy (e.g. blood flow simulations). Motivated by this, in order to have a better control on the number of sub-iterations needed in our algorithm, we decide to replace the Neumann condition is system (8) with a Robin coupling condition, as suggested by [20,43,44]. In [20] the authors propose a Robin coupling condition that is based on a coefficient α RO B that has been computed for the one dimensional structure. For our problem however, the solid is two dimensional and elastic: we therefore rely on the work presented in [53], where the authors compute the constant α RO B in the case of an elastic solid. We therefore just have to incorporate the expression of α RO B found in [53] into the Robin coupling condition presented in [20], and remember to pull back the condition onto the reference fluid-structure interface, using the ALE map. The final expression of the Robin coupling condition is: In Eq. (11), p i+1, and d i+1, s are suitable extrapolations of the fluid pressure and the solid displacement, respectively; we show in the next paragraph which kind of extrapolation we Condition (11) is imposed only on the fluid side: Robin conditions are indeed nonstandard in solid mechanics, therefore a lot of already existing codes would not allow to impose such a condition for the solid problem.

Space Discretization of the Semi-implicit Procedure
We now present the semi-discretized version of the algorithm introduced. We define the following function spaces for the fluid: )and the L 2 norm respectively, and the function space for the solid: norm. We discretize in space the FSI problem, using second order Lagrange Finite Elements for the fluid velocity, resulting in the discrete space while for the fluid pressure, the fluid displacement and the solid displacement we use first order Lagrange Finite Elements, resulting in the discrete space we further assume here that the fluid and the solid discretizations match at the FSI interface. The non-homogeneous boundary condition in system (9) can be easily treated by introducing, at timestep t i+1 , a lifting function i+1 such that i+1 = p in (t i+1 ) on in and i+1 = 0 on out ; we refer, for example, to [45,54] for more details concerning the use of a lifting function within model order reduction. By introducing the homogenized pressure p 0,i+1 The discretized version of the semi-implicit procedure reads as follows: for i = 0, . . . , N T ,

Fluid Explicit
Step This step results in a nonlinear system, which, at the computational level, is being solved with a Newton method. We remark that in the weak formulation, the boundary terms for the velocity vanish: this is in part due to the homogeneous Dirichlet boundary condition on top ∪ bott , and in part due to the fact that, in the Chorin-Temam scheme, the inlet boundary condition (first equation in system (4)) is split between the velocity and the pressure in the following way: The imposition of Neumann boundary conditions within a Chorin-Temam scheme is not at all trivial, and we refer the interested reader to a more detailed discussion presented in [55,56].

Implicit
Step for any j = 0, . . . until convergence: We iterate between the two implicit substeps, until a convergence criteria is satisfied; we choose as stopping criteria a relative error on the increments of the pressure and the solid displacement, namely: where ε is a fixed tolerance.
In the pressure Poisson formulation, to impose the Robin coupling condition, we have chosen the pressure at the previous implicit iteration, namely p i+1, j f , as an extrapolation for the fluid pressure, and the same goes for the extrapolation of the structure displacement.

POD and Reduced Basis Generation
For the generation of the reduced basis for the fluid velocity u f and the fluid displacement d f we pursue here the idea that was first proposed in [20]. For the homogenized fluid pressure p 0 f and for the solid displacement d s we employ a standard POD, giving rise to the reduced spaces Q 0 N and E s N respectively, though the authors would like to mention the fact that, as an alternative to the POD modes for the solid problem, the so-called vibrational modes can be used: these are obtained solving a generalized eigenvalue problem involving the mass and the stiffness matrix of the solid problem. Vibrational modes show very good results, especially for linear problems, and the authors refer the interested reader to [38,57].

Change of Variable for the Fluid Velocity
The main idea here is to introduce a change of variable in the fluid problem, in order to transform the non homogeneous Dirichlet condition at the FSI interface in system (13) into a homogeneous boundary condition. The motivation of this choice is that, to impose the second condition in system (13), we could use a Lagrange multiplier λ, thus increasing the dimension of the system to be solved in the online phase. In order to avoid this and in order to design a more efficient reduced method, we choose to transform the non-homogeneous coupling condition into a homogeneous one: we refer to [58] for a detailed discussion on the treatment of non-homogeneous Dirichlet boundary conditions within a model order reduction framework. We begin by defining a new variable z i+1 f ,h : With this change of variable, the second equation in (13) is now equivalent to the homogeneous boundary condition for the new variable: for which no imposition by means of Lagrange multiplier is needed. Therefore, during the offline phase of the algorithm, at every iteration i + 1, after we have computed the velocity u i+1 f ,h , we compute the change of variable z i+1 f ,h . We then consider the following snapshots matrix: We then apply a POD to the snapshots matrix S z and we retain the first N z POD modes 1 z , . . . , N z z . We therefore have the reduced space: and now it is clear that, since every k z satisfies the condition k z = 0 on F SI , then also every element of V N will satisfy the same condition.

Harmonic Extension of the Fluid Displacement
In order to generate the reduced basis for the fluid displacement d f , we pursue again the idea presented in [20]. We start by generating the snapshots matrix related to the solid displacement: We then apply a POD to the snapshots matrix and retain the first N d POD modes 1 d s , . . . , N d d s , thus defining the reduced space for the solid problem: We then employ an harmonic extension of each one of the reduced basis k d s to the fluid domain, thus obtaining the functions k d f such that: We impose homogeneous Dirichlet boundary conditions on the remaining part of the boundaries. We can then define the reduced space for the fluid displacement: The reason for for this choice, instead of employing a standard POD on the set of snapshots of d f , is given by the fact that we can avoid the introduction of another Lagrange multiplier to impose the non-homogeneous boundary condition present in system (22). We avoid to solve the reduced system related to (22): instead of solving an harmonic extension problem at every time-step in the online phase, we solve once and for all N d harmonic extension problems in the expensive offline phase. Then, during the online phase, the reduced fluid displacement will be computed just as a linear combination of the basis i d f , with coefficients that are the coefficients of the reduced solid displacement at the previous time-step. We will see in the next section the final formulation of the online phase of the algorithm. Before moving on, we summarize the offline computational phase, with the aim of helping the reader to better understand the whole procedure so far. Let i + 1 be the index of the current time iteration: f ,h , using the previously computed snapshot d i s,h ; 2. Solve the fluid explicit part, and find u i+1 f ,h in the snapshot matrix S z ; 6. Iterate until tolerance ε is reached: • Solve the pressure Poisson problem, using the solid displacement at the previous subiteration, and find p i+1, j+1 f ,h ; • Solve the solid problem, using the fluid stress tensor

Online Computational Phase
We are now ready to present the online formulation of the partitioned procedure, which is obtained by means of a Galerkin projection over the reduced spaces Then the online phase of the partitioned procedure reads as follows:

Mesh Displacement
let d i+1 f ,N be defined by the reduced solid displacement at the previous time-step: we then recover the reduced fluid pressure p

Numerical Results
We now present some numerical results obtained with the semi-implicit scheme. The reference physical configuration of the problem of interest is the one represented in Fig. 2: the geometrical properties of the domain are reported in Table 1; the leaflets are situated 1 cm downstream the inlet boundary. For our simulation we used a time-step T = 10 −4 , and a final time T = 0.05 s, for a total of N T = 500 iterations. Figure 3 shows the mesh used for the spatial discretization of the original problem.  The values of the physical constants used in the simulation are reported in Table 1. A pressure impulse p in (t) is applied at the inlet boundary, and after some time this impulse becomes constant: where T in = 0.1 s. We fix a tolerance of ε = 10 −6 as a stopping criterion for the subiterations between the pressure Poisson problem and the solid problem.
Since we do not consider the top and the bottom walls of the fluid domain to be deformable, we impose a homogeneous boundary condition for the fluid velocity on these walls. Figure 4a shows the rate of decay of the first 100 eigenvalues associated with three unknowns of the problem, namely the change of variable for the fluid velocity change of variable z f , the pressure p f and the solid displacement d s . It can be noticed that the rate of decay of the eigenvalues for the pressure and for the fluid change of variable is slower Fig. 6 Reduced solid displacement d s,N at time-step t = 0.035s (left), t = 0.04s (center) and at time-step t = 0.05s (right). The displacement has been obtained with N d = 10 reduced basis. The displacement has been magnified for visualization purposes   Fig. 4b we can notice that the first mode of the solid displacement retains 2% more energy compared to the first mode of the pressure, and 8% more energy with respect to the first mode of z f , which is the one that retains less energy. Fig. 4a is also important to have a first insight on the dimension of the reduced spaces that we are going to take during the online phase: indeed the rate of decay of the eigenvalues returned by the POD gives us an idea of the behavior of the approximation error that we commit by approximating the FE solution with the RB one.
The following relationship holds true (we state it for the fluid pressure, but the same holds also for the other components): where N p is the orthogonal projector onto the POD space of dimension N p , and λ i p are the eigenvalues returned by the POD. Figures 5 and 6 show two reproductive reduced order solutions: the fluid velocity and the solid displacement, respectively; as we can see from Fig. 6, the reduced order model shows a good capability also in reproducing very small deformations: these results were obtained using N z = 15, N p = 10 and N d = 10 modes, respectively. Figure 7 shows that, with N z = 15, N p = 10 and N d = 10 basis functions for each component of the solution, we have a good relative approximation error behavior over time: as we can see from the figure, at the final timesteps of the simulation the error increases, and we think this is due to some error accumulation phenomenon. The error has been computed as the L 2 error for the fluid pressure, and as the H 1 error for the fluid velocity and the solid displacement: We were interested in seeing how the average approximation error and the average number of internal iterations changes, by changing the number of reduced basis N p used for the fluid pressure in the reduced order model: results are reported in Tables 2 and 3 respectively. As we can see from Table 2, the average approximation error decreases up to when we use N p = 25 modes for the pressure, then we observe an increment in the approximation error: we read this result as the fact that with 40 modes for the fluid pressure, we are just adding noise to the online system. It is also interesting to see from Table 3 that the average number of internal iterations required from the algorithm, in order to reach a coupling tolerance of ε = 10 −6 is relatively higher for the reduced order model, when compared to the full order one: this is due to the reduction of the two problems, but in any case we can see that this number stabilizes around 26 − 27 subiterations. Finally, Fig. 8 depicts the behavior with respect to time of the number of internal iterations: a comparison is drawn between the full order model and the reduced order one, where we used N p = 10 reduced basis functions. We can see that the number of internal iterations, both for the offline and for the online part stabilizes towards the end of the simulation. We would like to make the following remark: all these results are computed by varying the number of modes used for the approximation of the fluid pressure, while keeping fixed both N z and N d . The motivation behind our choice is the fact that we want to see how the number of modes directly impacts the performance of the method, and, more precisely, of the implicit step, where the coupling between the two physics is imposed, by coupling the pressure Poisson problem with the solid problem. The authors are aware that these results are by no means exhaustive, and this is a further testimony to the capability that such a partitioned procedure offers: many more tests are possible, where for example N d s is varied, and N z is kept fixed, or both can vary. For this reason, we now present some additional results, that have been instead obtained by using the same number of modes for all the components of the FSI solution, namely N z = N p = N d = N . Figure 9 shows the relative approximation error in time, for the three components of the FSI solution: as we can see, by increasing the number of modes, we get a better approximation error. In particular, from Fig. 9, we observe an oscillating behavior of the pressure relative error towards the last timesteps of the simulation: this behavior can be seen also in Fig. 7, and can be interpreted as the result of an accumulation phenomenon, where the error accumulates and starts to oscillate. Another important observation is that the online pressure, with this partitioned approach, has been obtained without the supremizer enrichment technique [54]; as already remarked in [20], this may lead to non optimal error convergence. Motivated by this, we think it can be a very interesting point for a future work to see if this has some implications in the error oscillation that we observe in Fig. 9. Figure 10 shows the number of internal iterations for each timestep of the simulation: the online computations are performed using N = 25 modes for all the components, and are compared to the high fidelity computations. Also in this case the average number is higher for the online simulation, due to model order reduction, and also in this case the number stabilizes around 17 iterations towards the end of the simulation. Finally, in Table 4 we report   the average number of internal iterations, for an increasing number N of modes used in the online phase. As we can see, by increasing N from 10 to 25 there is almost no improvement in the number of internal iterations: it stabilizes around 25. This number then drops to 16 (which is very close to the FOM results) for N = 25: we did not increase N further, because N d = 25 is the total number of modes retained by the POD on the solid displacement, before hitting a very small magnitude (less than 10 −9 ) for the corresponding eigenvalues.

Shape Parametrization of the Leaflets
In this section we are going to address a slightly different situation, the difference being now the presence of a geometrical parameter μ g , that represents the length of the leaflets; we also admit the possibility of a further physical parameter μ p , so that, to summarize, we consider a parameter μ ∈ P ⊂ R d , where d = 1 if just a geometrical parametrization is considered (and thus μ = μ g ), or d = 2 (and thus μ = (μ g , μ p )).

FSI in the Presence of Shape Parametrization
Let us denote by (t; μ g ) := f (t; μ g ) ∪ s (t; μ g ) the current physical domain: we now have a time dependence and a parameter dependence. We introduce the time-independent intermediate configuration˜ (μ g ) :=˜ f (μ g ) ∪˜ s (μ g ), where we are considering the reference configuration of both physics, still taking into account the parameter dependence. Finally, we have the time-independent, parameter-independent reference configurationˆ := We call T the shape parametrization map; for every μ g we have a map T μ g defined as follows: We then have the ALE map A f (t; μ), already introduced in Sect. 2, which is now a map from the current parametrized fluid configurationˆ f (t; μ g ) and the intermediate fluid con-figuration˜ f (μ g ): whered f is the mesh displacement already defined in Sect. 2.
Let us define the gradients and the determinants of the deformation maps: We can pull-back the gradientF(x; μ) to the reference domainˆ f , and we obtain F(x; μ) = Id +∇d f (μ)G −1 (x, μ g ). With this notation, we can conclude that the gradient of the deformation map from the reference configuration to the current configuration is given by F(x, μ)G(x, μ g ); let us denote byF μ , F μ and G μ g the gradientsF(x, μ), F(x, μ) and G(x, μ g ) respectively, and by J μ and K μ g the determinants of F μ and G μ g . We are now ready to state the strong form of the problem of interest.

Strong Formulation
The strong form of the parametrized FSI problem reads as follows: for every t ∈ [0, T ] and for every μ ∈ P, find the fluid velocity u f (t; μ) : f (t; μ g ) → R 2 , the fluid pressure p f (t; μ) : f (t; μ g ) → R, the mesh displacementd f (t; μ) :˜ f (μ g ) → R 2 and the solid displacementd s (t; μ) :˜ s (μ g ) → R 2 such that: Here we notice that, again, the fluid problem is formulated in the current parametrized configuration f (t; μ g ), whereas the solid problem is formulated in the parametrized intermediate configuration˜ s (μ g ). The quantity ∂ t u f |x represents the ALE time derivative: t; μ g ). Again, σ f is the fluid Cauchy stress tensor, andP is the first Piola-Kirchoff stress tensor: their definition has been given in Sect. 3. The previous system is completed by some suitable initial conditions, by the same boundary conditions prescribed in (4) and by the following coupling conditions: beingσ f the Cauchy stress tensor in the parametrized intermediate fluid domain˜ f (μ g ): Thanks to the introduction of the pull-back maps, we can reformulate our problem in the reference configurationˆ : for every t ∈ [0, T ] and for every μ ∈ P, find the fluid velocitŷ where:σ We have the coupling conditions and the following boundary conditions: Again,n represents the normal vector to the relative part of the boundary of the domain.

Remark 6
In this section we stressed the difference between entities on the current configuration, the parametrized intermediate configuration and the reference configuration, by using the superscriptsˆand˜. However, since from now on everything will be cast in the reference configuration, and in order to make the notation as light as possible, we will drop all the superscripts.

Offline Computational Phase
Hereafter we present the offline phase of the partitioned procedure in the presence of a parameter. We employ again a Chorin-Temam projection scheme for the Navier-Stokes equation; we define a time-stepping procedure by sampling the time interval [0, T ] with an equispaced sampling {t 0 , . . . , t N T }, where t i = i T , for i = 0, . . . , N T and N T = T T . We discretize the time derivative of a function f with a first backward difference: . We then discretize the parameter space P ⊂ R d , d = 1, 2 with an equispaced sampling, and we obtain P train = [μ 1 g , . . . , μ We define here N train to be the cardinality of the training set P train .
In the following, we use the same function spaces that we have introduced in Sect. 3.2.2: and E s ( s )) and the L 2 norm respectively. We remark that in the previous definitions, the domains f and s are the reference configurations (both parameter-and time-independent). Again we discretize in space the FSI problem, using second order Lagrange Finite Elements for fluid velocity resulting in the discrete space V h ⊂ V , while the fluid pressure, the fluid displacement and the solid displacement are discretized with first order Lagrange Finite Elements, resulting in the discrete spaces Q h ⊂ Q, E f h ⊂ E f and E s h ⊂ E s ; we make again use of a lifting function for the fluid pressure, thus we introduce also the discrete space Q 0 h , which is defined exactly as in Sect. 3.2.2.

Remark 7
In this case the lifting function (t) does not depend on the parameter μ g , as we can deduce from the fact that the quantity p in (t) is parameter-independend. Therefore we can compute the lifting function, during the offline phase, once and for all for every timestep t i . We would also like to refer the interested reader to the work presented in [58] in the case p in is parameter dependent: indeed, in [58] the authors present a detailed description of the work that has to be done in order to implement reduced order models in the presence of multiple parameters in the boundary data.
The space discretized version of the partitioned procedure now reads as follows: for i = 0, . . . , N T , for μ = (μ g , μ p ) ∈ P train :

Fluid Explicit
Step

Implicit
Step for any j = 0, . . . until convergence: subject to the boundary conditions (9). We then retrieve the original fluid pressure p = 0 on s D . In the fluid projection step, in order to enhance the stability of the method we have employed again a Robin boundary condition, which in the case of shape parametrization reads as follows: previous two test cases; in this way we obtain a set { k

Implicit
Step for any j = 0, . . . until convergence: we then recover the reduced fluid pressure p

Numerical Results: Geometrical Parametrization Only
We now present some numerical results concerning the parametrized version of the two dimensional FSI test case presented in Sect. 3. The original domain is shown in Fig. 11, together with the reference configuration, and the parametrized reference configuration. The fluid domain is represented in blue, while the solid (the leaflets) is depicted in red. The geometrical constants defining the physical domain are reported in Table 5. Only one geometrical parameter is considered here: the length μ g of the leaflets, where we have chosen μ g ∈ P = [0.8, 1.0]. An affine mapping T μ g is chosen to deform the reference domainˆ , obtained for μ g = 1.0 cm, to the parametrized configuration˜ (μ g ): such a map is computed analytically. Top and bottom walls of the blue domain are rigid, thus both the displacement d f and the fluid velocity u f are set to zero. Homogeneous Neumann condition is imposed on u f on the outlet; a pressure profile p in (t) is described at the inlet, where: for t ≤ 0.025s 5 for t > 0.025s, and T in = 0.1 s. Also in this case we set a tolerance of ε = 10 −6 as a stopping criterion for the subiterations between the pressure problem and the solid problem.
For the simulation, we use the same mesh used for the previous test case, we set t = 10 −4 , for a maximum number of time-steps N T = 500, thus T = 0.05s. Table 5 summarizes the details of the offline stage and of the FE discretization. The number of parameter samples used during the offline phase to train the algorithm varies between N g = 8 (for a total of 4000 snapshots generated) to N g = 16 (for a total of 8000 snapshots generated). Figure 12a shows the rate of decay of the eigenvalues returned by the final run of the POD on z f , p f and d s , respectively. As we can see, now the eigenvalues display and overall slower decay, showing that the complexity of the problem is higher with respect to the non parametrized case: this slower decay will then reflect into a higher number of modes that need to be used in the online phase. Figure 12b shows the energy retained by the modes returned by the final run of the POD: the first mode of the pressure is indeed the most energetic one, while the first mode of the velocity is the least energetic one. All these results were obtained by using N g = 16 sampling parameters. In Fig. 13 we depict the behavior in time of the number of iterations required in the implicit step, at the ROM level, to reach a tolerance of ε = 10 −6 , according to the number of modes used: for this test, the modes were generated using N g = 16 sampling parameters, for a total of 8000 snapshots generated, and the leaflet length is μ g = 0.84. It is very interesting to see that, except for N = 15, refining the number     This result highlights very well the deep influence that the geometrical parametrization has; indeed, the length μ g = 1.0 corresponds to the reference leaflets length that we chose: in this case the geometrical deformation map T μ g defined in Sect. 4.1 is the identity. We would like to stress also the following: the test case with μ g = 0.84 represents a prediction test case; indeed this value not only corresponds to a significant geometrical deformation, but it also corresponds to a value that has not been selected to train our algorithm at the offline level. For this reason, this test case is a stress test for our algorithm. In the framework of predictive problems, we believe that the result presented in Fig. 14 can be improved, by refining the geometrical parameter sampling set used to train our method, thus by generating, during the offline phase, more snapshots: these simulations are very expensive, and for this reason we did not proceed any further with the sampling set refinement. Figure 15 shows the reduced displacement, for the same values of the geometrical parameter: μ g = 0.84 and μ g = 1.0; again, the influence of μ g is clear: the longer the leaflets, the bigger their deformation is going to be, under the same physical parameters. Table 6 represent the average approximation error for the fluid velocity, the fluid pressure and solid displacement, with refinement of the training sample P train : all the reduced solutions have been obtained using N = 30 reduced basis functions for all the components u f , p f and d s ; it is interesting to see that, for the highest number of training samples, namely N g = 16 (which corresponds to a total of 8000 snapshots generated) we observe a slight increase in the average approximation error of the velocity and the solid displacement. We think this could be due to the fact that the online model could benefit from more reduced basis for u f and  Table 7, where we show the average approximation error for N g = 16, with basis refinement: as we can see, the model benefits from the increment in the number of modes used in the online simulations. Also here the average error is intended as average over time, and it is computed as the H 1 relative error for the velocity, the L 2 error for the pressure and the L 2 error for the displacement. Finally, in Fig. 16 we present the behavior of the relative approximation error in time: we consider the reduced solution to be obtained with N = 15, 30, 40, 50 modes for each component (see color legend in the figure). As we observe, the behavior in time of the approximation error seems to confirm the results reported in Table 7: our model benefits from an increase in the number of modes used. It is interesting to observe the behavior of the pressure relative error: the error accumulates over time, thus steadily increases. This represents a starting point for future studies and future work development; indeed, this steady growth in time of the approximation error is related only to the pressure component of the FSI solution. We ask ourselves if this is somehow related to the fact that we are not using a supremizer enrichment of the velocity space at the online level; it would be therefore interesting to study if this has some effect in the online approximation of the pressure, especially for long time simulations.

Numerical Results: Physical and Geometrical Parametrization
We present some numerical results for the test case with a geometrical and physical parametrization: now μ = (μ g , μ p ), where the physical parameter μ p represents the shear modulus of the leaflets, and thus μ p = μ s , the second Lamé constant. The original configuration, the intermediate configuration and the reference configuration are the ones depicted in Fig. 11. The geometrical constants of the problem are the same ones of the previous test case, and are reported in Table 5. The physical parameter μ s varies within the range [10 5 , 8 × 10 5 ]. We use the same boundary conditions, the same inlet pressure profile, time step and tolerance used in the previous test case. Figure 17 shows two different examples of the behavior of the leaflets, according to the change of the physical and/or of the geometrical parameters: all the pictures represent the online displacement of the solid at the final timestep of the simulation, namely for t = 0.05s.     The results have been obtained using N = 30 reduced basis for all the components; as we can see from Fig. 17a, b, an increase in the shear modulus leads to a material that is much more hard to deform. On the contrary, for a fixed value of the properties of the material under consideration, and increase of the length of the leaflets leads to an increase in the displacement. Figure 18 shows the behavior in time of the number of iterations of the implicit step, for different number of modes N used in the online phase, compared against the FOM. Figure 19 represent the relative error approximation for u f (top left), p f (top right) and d s (bottom center) as a function of time. From Fig. 19 we can see that the relative error for the fluid velocity stabilizes after some iterations, reaching a magnitude of 2 × 10 −1 when using N = 20, 30 and 40 modes. Also the relative error for d s shows a plateau around 9 × 10 −2 , except for when N = 20, in which case the error decreases in time: we read this result as the consequence of using too many modes in the online model, when N = 30, 50. For the fluid pressure we observe the same accumulation phenomenon that we observe in the previous parametrized test case. To conclude, even though these relative errors seem high, we are again testing our algorithm with a prediction problem, since the value (μ g , μ s ) = (0.9, 10 5 ) has not been used to generate modes during the training phase of the algorithm. In addition to this, we would like to remark that, by increasing the number of sampling parameters used at the FOM level, we should be able to drive the error down, and this seems to be confirmed by the results in Table 10: we did not proceed with a sampling refinement because of time constraints, since the generation of 40, 000 snapshots is very demanding. In Tables 8, 9 and 10 we report the average approximation error for the fluid velocity u f , the fluid pressure p f and the solid displacement d s , when μ g = 0.9 cm and μ s = 10 5 : the average has been taken over time, and we used the L 2 norm for the fluid pressure and the solid displacement, and the H 1 norm for the fluid velocity. In Table 8 we computed the relative approximation error, by using N g = 8 and N s = 10 training samples, and refining the number of basis functions used online, from N = 15 to N = 50. In this case, we used the same number of basis functions for all the components of the FSI solution: again, we remark that other tests are possible, for the interested reader, for example keeping the number of modes fixed for the fluid velocity and the solid displacement, and varying the number of modes used for the fluid pressure. Also here, like for the previous test case, the online model seems to benefit from a higher number of modes used. In Table 9, we used N = 30 modes for the online solution, N s = 10 training samples for the physical parameter, and we refined the geometrical training set, from N g = 5 to N g = 8. On the contrary, in Table 10, we kept N g = 8 training samples for the geometrical parameter, and refined the physical training set from N s = 6 to N s = 10. We remark that analyzing the training samples with N g = 8 and N s = 10 corresponds to generating 40.000 snapshots. Indeed, in the framework of the POD, we have to compute all the snapshots for each value of the training parameter, in an unsteady framework: this procedure is extremely expensive, and it required a total of 5 days, by running the simulations on two different computers: a computer Intel(R) Core(TM) i5-4670 S CPU with 3.10 GHz and 7.7 G of RAM, and a supercomputer Intel(R) Xeon(R) CPU E5-2687W v4 with 3.00 GHz and 540 G of RAM. For this reason, we did not proceed further with the refinement of the parameter sampling. As the results show, however, by refining the parameter space, thus by using more training samples, we are able to improve the average approximation error. We remark also that, in case of Tables 9 and 10, the online computations have been made using N = 30 modes for each component of the solution. We think that in this case the average approximation errors are higher (compared to the ones obtained for the non parametrized test case) expecially due to the presence of the geometrical parametrization. Finally, in Table 11 we show some CPU times: we choose μ g = 0.9 and μ s = 8 × 10 5 , and fix N = 30 for all the components (these modes were obtained fixing N g = 8 and N s = 10). The simulation is run on a computer Intel(R) Core(TM) i5-4670 S CPU with 3.10 GHz and 7.7 G of RAM; the CPU execution time is reported in Table 11 (times measured in seconds): the FOM requires 24763.47 seconds for the same geometrical and physical parameters (the computation comprises also the computation of the change of variable z f and the homogenization of the fluid pressure through the lifting function). As we can see, and as expected, by strengthening the tolerance, the computational time grows, as more iterations are needed at the implicit step to reach the convergence. What is interesting to notice, however, is that, by strengthening the tolerance we do not see an important decrease in the average approximation error: this result seems to therefore suggest that online computations could be carried out also using a coarser tolerance (for example 10 −4 instead of 10 −6 ).

Conclusions
In this manuscript we presented a Reduced Order Model algorithm designed to address FSI problems, in the unsteady case, and possibly in the presence of a parameter dependence. The ROM is based on a partitioned procedure: the main advantage of this is given by the fact that, by solving separately the fluid and the solid problem, we are not only able to lower the dimension of the systems to be solved in the online phase, but we also have a better control on the number of variables that are needed. The reduced basis functions are generated through a Proper Orthogonal Decomposition on the set of snapshots, with the introduction of a change of variable in the fluid problem formulation. The procedure that we have proposed aims at extending the work presented in [20,45] to the case of the coupling between an incompressible fluid and a thick, two dimensional structure, also in the presence of geometrical parametrization. The results that we have obtained confirm the following aspects: • Introducing a change of variable in fluid explicit step allows us to avoid the introduction of a further unknown in the system, namely a Lagrange multiplier, in order to impose non homogeneous boundary conditions at the fluid-structure interface; • The choice of not performing a POD on the snapshots d f , but rather performing an harmonic extension of some modes, allows us to build the online d N , f in a cheap way. Moreover, the coupling condition that imposes the continuity of the displacements at the fluid-structure interface is automatically satisfied, thanks to the way we have defined the reduced basis for d f .
In addition to the list of remarks presented, another very important detail of the procedure presented in this manuscript is the following: we did not rely a supremizer enrichment of the fluid velocity space, as it is the usual case for reduction methods, in order to obtain a stable approximation of the fluid pressure in the online phase. Our choice is motivated by the fact that, even at the Finite Element level, the Chorin-Temam projection scheme with the pressure Poisson formulation can be applied succesfully also to velocity-pressure FE spaces that do not satisfy the inf-sup condition, see [60]. This allows us to limit the dimension of the system to be solved online, and it is a big motivation for the choice of a Chorin-Temam projection scheme within our partitioned approach. While testing this algorithm, we have seen that a partitioned procedure is demanding from the computational time point of view: this drawback is represented by the fact that, in the imposition of the coupling conditions through a Robin boundary condition, the constant α RO B that makes the procedure more stable depends on the time-step used. If we choose a time-step that is too big, then our Robin coefficient α RO B becomes very small, and we recover the original Dirichlet-Neumann coupling, which is known to have stability problems, i.e. the implicit step may not converge. Always in the direction of the computational effort of the offline phase, we remark that in this manuscript we considered a physical parameter for the solid, but no parametrization of the fluid has been taken into account: considering a fluid parameter, instead of a solid one, does not change the design of the algorithm. However having a test case which comprises both a parameter for the solid and a parameter for the fluid does highlight the boundaries of the POD in this framework; indeed, generating the snapshots for the chosen sampling set would be extremely time demanding, even though it would be carried out at an offline level: for this reason, future perspectives also include the design of an error estimator and the integration of a Greedy algorithm based on such estimator. It is also important to mention the fact an efficient online-offline decoupling is very important in terms of model order reduction efficiency: for the fluid part, this can be recovered thanks to the Empirical Interpolation Method (EIM), see for example [17,61,62], whereas for the solid mechanics part, one may think about using hyper-reduction procedures that preserve the Hamiltonian structure, such as, for example, the ECSW [63]. These techniques have not been used in this work, as we wanted to focus on the development and test of a reduced order segregated procedure for FSI problems which involve the coupling of an incompressible fluid with an elastic structure; this can be seen as a natural future perspective for this work.
As it was mentioned in the Introduction, there is currently a fair amount of interest in approaches that are able to couple high fidelity models with reduced order models: one may think about using a FOM for the structure, and a ROM for the fluid. The authors believe that there is a possible extension (with some required modifications) of the algorithm proposed, to this kind of situations: again, this represents another interesting future line of research. Another remark that we would like to make concerns the timestep used in our numerical simulations: in the results that we have showed, a timestep of T = 10 −4 has been used both for the fluid and for the solid problem. However, one may be interested in using different timesteps for the two physics, as the solid model is an hyperbolic equation, and the fluid is a parabolic one: we do believe that this is possible, with the procedure presented. Indeed, given the time interval I n := (t n−1 , t n ], given a solid timestep T s and a fluid timestep T f , assume that the resulting time discretization I s n of the interval I n is finer than the time discretization I f n of I n . Then, one just needs to be able to evaluate fluid quantities on the times t s i ∈ I s n and solid quantities on the times t f i ∈ I f n : with the definition of suitable interpolation operators (interpolation in time), we should be able to implement a partitioned scheme with different timesteps for the two physics. This idea is presented for example in [64] for a monolithic scheme, but it represents an interesting starting point for a future work within partitioned schemes. In this case, maybe manually tuning the Robin parameter α RO B to achieve optimal convergence can be a better idea, in order to drop the dependence on T : further research in this direction needs to be carried out.