On the Choice of Interface Parameters in Robin–Robin Loosely Coupled Schemes for Fluid–Structure Interaction

: We consider two loosely coupled schemes for the solution of the ﬂuid–structure interaction problem in the presence of large added mass effect. In particular, we introduce the Robin–Robin and Robin–Neumann explicit schemes where suitable interface conditions of Robin type are used. For the estimate of interface Robin parameters which guarantee stability of the numerical solution, we propose a new strategy based on the optimization of the reduction factor of the corresponding strongly coupled (implicit) scheme, by means of the optimized Schwarz method. To check the suitability of our proposals, we show numerical results both in an ideal cylindrical domain and in a real human carotid. Our results showed the effectiveness of our proposal for the calibration of interface parameters, which leads to stable results and shows how the explicit solution tends to the implicit one for decreasing values of the time discretization parameter.


Introduction
The numerical solution of fluid-structure interaction (FSI) problems is very challenging and many different strategies have been considered so far. Among them, we mention monolithic strategies where the whole space and time discretized problem (e.g., due to Finite Elements) is solved by means of efficient linear solvers and ad-hoc preconditioners that should accelerate convergence which in general is quite problematic due to the elevated condition number of the FSI problem, see e.g., [1][2][3][4][5][6]. Another family of strategies relies on partitioned or segregated schemes which introduce the separate solution of the fluid and structure (and possibly fluid geometry) subproblems. Among them, strongly coupled or implicit strategies solve such problems until convergence of the interface conditions at each time step, see, e.g., [7][8][9][10][11][12][13][14][15][16].
Within partitioned schemes, loosely coupled or explicit schemes for the numerical solution of the FSI problem are based on an overall explicit time discretization which leads to the solution of just one fluid and one structure problem per time step. This makes this family of methods very attractive from the computational and implementation point of view and for these reasons they have been widely used in many engineering applications such as aeroelasticity [17][18][19]. However, loosely coupled schemes suffer from a lack of stability when the added mass effect is relevant (i.e., when the densities of fluid and structure are comparable). This happens, for example, in hemodynamics [20]. In this respect, it is known that the explicit Dirichlet-Neumann (DN) scheme is unconditionally unstable in the hemodynamic regime, see [7,11,21].
In such works different proposals for the interface parameters were addressed with the aim of improving the stability properties when the added mass effect is relevant with respect to the explicit DN scheme. In our recent study [29], we have provided for a model problem a stability analysis of the explicit Robin-Neumann (RN) scheme. In particular, we have found sufficient conditions for the interface Robin parameter, guaranteeing both unconditional instability and conditional stability.
In this paper, we address the issue of selecting suitable and easily computable interface parameters both for the explicit RN and the explicit Robin-Robin (RR) schemes, able to guarantee the stability of such loosely coupled schemes. In particular, we discuss the case of cylindrical-like geometries as happens, for example, in vascular hemodynamics and we propose a new and effective way to select such parameters. We start from the analysis based on the optimized Schwarz method [30] provided for the implicit (i.e., strongly coupled) RR scheme in the FSI context in [31]. This allowed us to determine effective values for the Robin interface parameters which guarantee excellent converge property even in presence of large added mass effect (for strongly coupled scheme, a large added mass effect yields a very slow convergence [7]). Here we provide also a new way to easily estimate an effective Robin interface parameter for the RN strongly coupled scheme. The idea of the present work is to use such estimates in the corresponding loosely coupled RN and RR schemes. In particular, we verify the stability of the corresponding numerical solution in 3D FSI numerical experiments.

The Continuous Problem
We consider the coupling between the Navier-Stokes equations for an incompressible fluid solved in the Arbitrary Lagrangian-Eulerian (ALE) formulation [32] and the linear infinitesimal elasticity [33]. Let Ω f and Ω s be the fluid and structure domains, Σ the fluidstructure interface, Σ out the external structure surface, n = n f the unit normal outgoing the fluid domain, and n s the unit normal outgoing the structure domain. We have for each t [33]: where T f (u, p) = −pI + µ(∇u + (∇u) T ) is the Cauchy stress tensor for the fluid and with µ the dynamic viscosity. ∂ A t represents the ALE time derivative, i.e., with respect to the ALE framework, and ω is the velocity of the fluid domain obtained by solving an harmonic extension of the interface velocity with homogeneous Dirichlet or Neumann boundary conditions on ∂Ω f \ Σ. Notice that, accordingly, Ω f changes in time. Instead, the structure problem (1e) is solved in a Lagrangian framework and for this reason we have indicated with the corresponding quantities. For the sake of notation, in what follows will be understood. Moreover, T s is the structure Cauchy stress tensor given by where λ 1 and λ 2 are the Lamé constants that can be defined in terms of the Young modulus E and the Poisson ratio ν as follows Finally, we observe that condition (1f) represents a Robin condition at the external surface to account for the effect of an elastic surrounding tissue with elasticity modulus γ ST [34]. The previous problem needs to be completed with other boundary conditions and initial conditions for both fluid and structure. Examples of applications where the FSI problem (1) has been considered are hemodynamics, where blood interacts with the vascular wall [33], the respiratory system [35], and aeroelasticity for the study of airfoils [17].

Robin Robin Loosely Coupled Scheme
In order to write a suitable algorithm for the numerical solution of the FSI problem (1), we consider the following linear combinations of the interface conditions (1c)-(1d), for given scalars α f = α s : The coupled FSI problem (1) where (1c)-(1d) are substituted by (2a)-(2b) is still equivelent to (1), thus we can introduce suitable numerical strategies based on the exchange of conditions (2a)-(2b). To this aim, we first need to detail the time discretizaton and how we manage the geometric coupling, i.e., the fact that the fluid domain movement depends on the structure displacement. Regarding the time discretization, we used a first order implicit method for both fluid and structure, with a semi-implicit treatment of the fluid convective term, relying on a CFL-like bound for the time discretization ∆t. We also consider an explicit treatment of the no-slip condition, allowing in fact to split the two subproblems. Regarding the geometric coupling, it has been shown that in the hemodynamic regime an explicit treatment is enough to provide stable and accurate results [12,14,22,36,37]. This means that the harmonic extension problem for the fluid domain displacement and velocity is solved with structure data that comes from previous time steps. This in fact decouples the geometric and FSI problems, thus at each time step the time discretization of problem (1b)-(1c) is in fact solved in a known domain Ω f . Let t n = n∆t, n = 0, . . . , the discrete time instants and v n v(t n ) the approximation at time t n of a function of time v(t). Thus, for the numerical solution of problem (1b), (1c), (2a), (2b), (1e), and (1f), we introduce in Algorithm 1 the Explicit Robin-Robin loosely coupled scheme, obtained after time discretization and by prescribing condition (2a) as boundary condition for the fluid problem, with structure quantities taken from the previous time step, and condition (2b) to the structure problem.

Remark 1.
As observed, at the continuous level conditions (2a) and (2b) are perfectly equivalent to (1c) and (1d). After the numerical discretization and the selection of an explicit treatment, we obtain conditions (3c) and (4b) which satisfy the original interface conditions up to an error of the order of ∆t.

Remark 2.
The implicit (strongly coupled) Robin-Robin scheme is obtained from Algorithm 1 by replacing the right hand side of (3c) with α f η n+1 −η n ∆t + T s (η n+1 ) and then subiterating with the structure problem (4).

Algorithm 1 Explicit Robin-Robin scheme
Given two scalars α f = α s and quantities at previous time steps, at time step t n+1 solve in sequence: 1: A fluid problem with a Robin condition at the fluid-structure interface: 2: A structure problem with a Robin condition at the fluid-structure interface: In what follows we discuss the choice of the interface parameters α f and α s in the case of cylindrical-like geometries and interface, a situation which occurs in many applications with large added mass effect, e.g., in hemodynamics.

Convergence Analysis of the Implicit Robin-Robin Scheme
Our starting point is the optimization procedure to properly select the interface parameters in the implicit Robin-Robin scheme for a simplified FSI problem in the case of cylindrical geometries in [31], whose main results are here reviewed for the sake of completeness.
We consider the problem arising from the interaction between an incompressible, inviscid, and linear fluid occupying the fixed domain Ω f = {(x 1 , x 2 , y) ∈ R 3 : x 2 1 + x 2 2 < R 2 }, and a linear elastic structure modeled with the wave equation occupying the domain In Algorithm 2 we report the implicit Robin-Robin scheme at time t n+1 for the solution of this simplified FSI problem. Actual temporal index n + 1 is understood. Notice that the coupling occurs only in the radial direction r since the fluid is inviscid. We have indicated with u r and η r the radial fluid velocity and structure displacement, respectively, and with F 1 and F 2 terms coming from the previous time step.
Notice that in [31] we considered general operators S f and S s to build the interface linear combinations. Here for the sake of exposition, we limit ourselves to the scalar constant case since the forthcoming optimization is performed over the subset of the scalars. Algorithm 2 Implicit Robin-Robin scheme for the simplified FSI problem Given two scalars α f = α s and quantities at previous time steps, solve for k ≥ 1 until convergence: 1: A fluid problem with a Robin condition at the fluid-structure interface: 2: A structure (wave) problem with a Robin condition at the fluid-structure interface:

Selection of Effective Interface Parameter Values for the Explicit Robin-Robin Scheme
Following [31], set where k ≥ 0 and m = 0, 1, 2, . . . are the frequencies related to the axial and circumferential coordinates, respectively, and which belong to the set K, I m and K m are the modified Bessel functions [38]. Then, in [31] it has been proven, through the optimized Schwarz method, that the reduction factor related to Algorithm 2 is given by In particular, it has been proposed to look for parameter values along the straight line α f = p and α s = −p + 2M for varying p ∈ R. With this specific choice, the reduction factor (6) becomes and it has been proved that it satisfies The previous result provided an easy way to compute a range of values of p (and thus of α f and α s ) which guarantees convergence for any frequencies (m, k) ∈ K. Moreover, this range contains the optimal value p * of p which minimizes the reduction factor for the choice α f = p and α s = −p + 2M, which could be easily found manually. The efficiency of such procedure has been shown in [31] both in ideal and in realistic carotid geometries in the context of hemodynamics. An extension to the case of spherical geometries and interfaces has been provided in [39].
The idea proposed in this paper is to use the range p ∈ [p − , p + ] given by (8) to properly select the interface parameters in the explicit Robin-Robin Algorithm 1, still with the specific choice α f = p and α s = −p + 2M. In particular, we propose here to use the optimal value p * also for the explicit RR scheme. Indeed, we expect that the interface Robin parameters that guarantee a fast convergence in the implicit case should guarantee stability and accuracy for the explicit case. Although there is not yet a proof of this, we provide here an experimental analysis to support our choices, see Section 4.

An Alternative Way to Select the Interface Parameter in the Explicit Robin-Neumann Scheme
We consider now the implicit Robin-Neumann scheme, i.e., Algorithm 2 with α s = 0. We propose here a new way to efficiently select the parameter α f for this algorithm. In particular, we still look for α f = p with p a scalar independent of the frequencies.
We have the following result.
Theorem 1. Suppose to have for a given iterative algorithm a reduction factor of the form for suitable scalar functions A(k) and B(k) defined on K, where k is a general scalar or vector variable, and with A(k)B(k) = 0 for all k ∈ K. Then, by setting for any k ∈ K provided that Proof. Notice that for any θ ∈ [0, 1) we have By imposing a + θb i.e., for Now, setting we can write the reduction factor (9) as follows One only has to observe now that inequality (10), that is which makes sense for any k ∈ K owing to (12). Thus the thesis follows.
We can now apply the above result to our case, i.e., the implicit Robin-Neumann scheme (Algorithm 2 with α s = 0), interpreting k in Theorem 1 as the couple of frequencies (m, k) and by taking A and B as in (5a)-(5d). Indeed, from (6), by setting α f = p and α s = 0, we obtain We now only need to prove that, as requested by Theorem 1, b < a. To this aim, we prove the following result.
with β and χ given by (5c) and(5d). When k = 0, the expression for b(m, k) has to be intended as Since K m /I m is decreasing on (0, +∞), it suffices to show that for all x > 0, and this follows immediately since K m (x) < 0 and I m (x) > 0. Similarly, the denominator of a(m, k) is positive if and only if is decreasing in (0, +∞), it suffices to show that for all x > 0, and this follows immediately since K m (x) > 0 and I m (x) > 0.
Thus, all hypotheses of Theorem 1 are satisfied by our application. The idea we propose in this paper is to use also for the explicit RN scheme (Algorithm 1 with α s = 0) an optimized value p * found in the range of convergence of its implicit counterpart (11).

Generalities
In this section we report numerical results aiming at showing the effectiveness of our proposal for the interface Robin parameters in the explicit RR and RN schemes. In particular, we want to understand if the values of the parameters derived by the simplified FSI problem (see Algorithm 2) work well for the complete three-dimensional FSI problem (1) solved by means of Algorithm 1. The use of simplified FSI models to perform analyses whose results are then used for more complex problems is a standard procedure due to the difficulty in analyzing directly such problems. For example, simplified FSI problems have been used to study the convergence of strongly coupled partitioned procedures, then successfully tested over 3D general problems, in [7,9,16,29,40], and to derive the well-known result about the instability of the explicit Dirichlet-Neumann scheme for large added-mass effect in [7,21].
All the simulations are run in the hemodynamic regime, characterized by a large added mass effect and where the stability of loosely coupled methods is a challenging issue.
We consider problem (1) and for its time discretization we used the BDF schemes of order 1 for both fluid and structure, with a semi-implicit treatment of the fluid convective term. For the space discretization we used P1bubble − P1 Finite Elements for the fluid and P1 Finite Elements for the structure. The fluid domain at each time step is obtained by extrapolation of the previous time step (semi-implicit approach [12,14,37]). We also used the following data: fluid density ρ f = 1 g/cm 3 , fluid viscosity µ = 0.035 g/(cm s), structure density ρ s = 1.1 g/cm 3 , Young modulus E = 3 × 10 6 dyne/cm 2 , Poisson ratio ν = 0.49.
Notice that to compute A(m, k) given by (5a) which is needed for the calibration of the interface parameters, we need the value of λ in the wave equation representing the structure problem in the simplified FSI problem. To do this, we assumed that the value of λ could be approximated by Gλ 1 , with G = π 2 /12 the Timoshenko correction factor.

Test in the Cylinder-Test I
In the first numerical test (test I), the fluid domain is a cylinder with length L = 5 cm and radius R = 0.5 cm, whereas the structure domain is the external cylindrical crown with thickness H s = 0.1 cm. The meshes are composed by 4680 tetrahedra and 1050 vertices for the fluid and 1260 vertices for the structure.
At the inlet we prescribed a Neumann condition with the following pressure function with absorbing resistance conditions at the outlets [22,41].
If not otherwise specified, we used the following parameters (referred to as "basic"): P = 500, ∆t = 0.0005 s, γ ST = 1.5 × 10 6 dyne/cm 3 . All the numerical results have been obtained with the parallel Finite Element library LIFEV [42].
We will refer to RR-explicit simulation when using α f and α s selected in the range (8) described in Section 3.2, whereas to RN-explicit simulation when using α f in the range (11) as in Section 3.3 and α s = 0. We reported also the numerical solution obtained by using an implicit method, in particular the Robin-Neumann scheme, with an absolute tolerance of 10 −7 on the convergence of the interface conditions.
In Table 1 we reported the values of the optimized Robin interface parameters a priori estimated via an empirical procedure. In particular, we take K = [m min , m max ] × [k min , k max ] with m min = 0 and m max = 10, k min = π/L = 0.6 and k max = π/h = 12.5 (remember that m and k are the angular and longitudinal frequencies) and we empirically look for p that minimizes max (m,k)∈K ρ(m, k) when either α f = p, α s = 2M − p and p varies in the range (8) (RR-explicit), or α f = p, α s = 0 and p satisfies (11) (RN-Explicit). Notice from the analyses reported in Sections 3.2 and 3.3 that any change of ∆t and γ ST influences the estimates of the interface parameters, whereas the choice of P does not.
In Figure 1 we report the time behavior of the mean pressure at the section located at half of the pipe (z = 2.5 cm) for different values of ∆t. In Figure 2, instead we show the same quantity in the case of an increased Reynolds number ( P = 5000) for two values of the surrounding tissue parameter (γ ST = 3 × 10 6 dyne/cm 3 and γ ST = 3 × 10 6 dyne/cm 3 ) and for a reduced values of the time step (∆t = 1.25 × 10 −4 ). From these results, we observe that the RR-Explicit and RN-Explicit solutions are in any case stable and feature a behavior which is reasonable if compared with the implicit solution, also depicted in the figures. As expected, decreasing ∆t the two solutions tend to coincide with the implicit one, see Figures 1 and 2, top and middle. In addition, from Figure 2 we notice that the performances of the proposed explicit schemes seem to be robust with respect to the Reynolds number and the value of the surrounding tissue. In any case, the two explicit solutions seem to be very similar, thus the new strategy proposed in Section 3.3 to build an explicit RN scheme could be an effective way to obtain a loosely coupled scheme characterized by only one parameter. Remark 3. Sometimes in the literature it is asserted that for FSI problems, among the solutions obtained with a strongly coupled (SC) scheme and with a loosely coupled (LC) scheme, the former is the reference one and the latter should be close to it to be considered accurate. However, we believe that a SC scheme features in general only better stability properties (being associated to an implicit time discretization) than a LC one (which is associated to an explicit time discretization). For values of ∆t which guarantee stability of LC schemes, the accuracy of SC and LC methods is the same if discretization methods of the same order have been used. For example, in our case both methods are first order in time. The fact that the explicit solution seems to be more diffusive than the implicit one for not enough small values of ∆t (see Figure 1, top and middle) could be ascribed to the fact the our Robin based LC schemes could be seen as stabilized methods where α f and α s play the role of stabilization parameters, whose influence on the solution increases with ∆t.

Test in a Human Aortic Abdominal Aneurysm-Test II
In the second test (test II), we consider a human aortic abdominal aneurysm (AAA) reconstructed from CT images. At the inlet we prescribed the physiological Neumann condition (13) with P in = 53, 320 sin(5πt) dyne/cm 2 t ≤ 0.2 s, 0 dyne/cm 2 0.2 s < t ≤ 0.8, with absorbing resistance conditions at the outlets. For the structure we imposed the Robin condition (1f) with γ ST = 3 × 10 6 dyne/cm 3 at the external surface, and fixed inlet and outlets. We set ∆t = 0.001 s. We considered the explicit RN and RR scheme with parameters estimated as described in Sections 3.2 and 3.3, respectively, by using the radius of the aneurysm (1.9 cm) as representative value for R. Moreover, we have used H = 0.17 cm for the structure thickness. The fluid and structure meshes were formed by 153 k and 74 k tetraedra, respectively, leading to the following ranges for the frequencies: k ∈ [0.3, 19.6] cm −1 and m ∈ [0, 38]. In particular, the optimized interface parameters found by our analyses were α f = 2174 g/(cm 2 s), α s = −85 g/(cm 2 s) for the explicit RR scheme and α f = 2229 g/(cm 2 s) for the explicit RN scheme. In Figures 3-5 we report the pressure field, the fluid velocity magnitude, and the vessel displacement magnitude at t = 0.05 s and t = 0.1 s (the latter being the systolic peak) for the explicit RR and RN schemes, together with an implicit solution. From these results we can observe stability of the proposed explicit schemes, reached for a value of ∆t which is commonly used in hemodynamic applications to guarantee accuracy. Moreover, we observe very similar results between the explicit and implicit solutions, much more similar than those reported in Figure 1, upper, for the same value of ∆t. This is probably due to the lower frequency of the physiological input signal used in test II with respect to that used for test I. All these observations are very encouraging in view of using the proposed loosely coupled schemes for real hemodynamic applications, a field where the strong added mass effect made the use of loosely coupled schemes a very challenging issue.

Final Remarks
We have considered a loosely coupled (explicit) partitioned algorithm based on Robin interface conditions for the numerical solution of FSI problem in the presence of large added mass effect. In particular, we focused on the selection of the interface parameters in the coupling conditions, which plays a crucial role in the stability of the method. The novelty of the paper relies on the description of an effective way to select the parameters in the interface conditions. This choice is crucial to improve stability of the proposed schemes, which, in general, is very poor when added mass effect is large. All the reported results (obtained in the hemodynamic regime when the added mass effect is very large) showed the stability of the numerical solution obtained by using the interface parameters proposed in this work.
The interest in using a partitioned scheme relies in its modularity, i.e., in the fact that we can use pre-existing separate fluid and structure codes and no ad-hoc implementation of a FSI solver is needed. Moreover, regarding the choice of an explicit/loosely coupled method, this is very attractive since, if stable, it allows for the reduction of computational costs with respect to implicit/fully coupled methods, if the need of using a smaller time step is compensated by the reduction of number of the subproblems we have to solve for each time step. In [29] we show for a test case that the CPU time of the explicit RN scheme was reduced by about three times with respect to the implicit RN case. Regarding the validation of our results, we observe that the results obtained by the library LifeV have been compared with analytical solutions and with clinical or experimental data. For the fluid solution, CFD results in the context of hemodynamics were successfully compared with experimental results in [43], with ECD measures in [44] and with PC-MRI data in [45]. Regarding FSI studies, in [41] the authors validated the numerical solution of implicit RR methods by comparing the results with an analytical solution.
Regarding test II in a human AAA, it provides the first results obtained with the proposed loosely coupled schemes in a real geometry in the context of hemodynamics and using data all in the hemodynamic regime. The results are stable and very similar to the implicit ones and highlight the ability of the method of representing the wave propagation and the deformation of the domains. FSI results could provide important clinical indications about the evolution of a disease or the improvement of a therapy. We mention, for example, studies that quantified the stresses in carotid atherosclerotic plaques [46,47] or compared different surgical strategies for coronary bypasses [48,49]. However, often the application of FSI methods in such contexts has been limited due to the very high computational effort needed to numerically solve such problem. Our test in human AAA gives an important preliminary answer toward the use of explicit Robin-Robin and Robin-Neumann methods for clinical applications, with a dramatic savings in computational times with respect to strongly coupled (implicit) schemes. Author Contributions: Conceptualization, G.G. and C.V.; methodology, G.G. and C.V.; software, G.G. and C.V.; validation, G.G. and C.V.; formal analysis, G.G. and C.V.; investigation, G.G. and C.V.; data curation, G.G. and C.V.; writing-original draft preparation, G.G. and C.V.; writing-review and editing, G.G. and C.V.; visualization, G.G. and C.V.; supervision, G.G. and C.V.; project administration, G.G. and C.V. Both 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.
Data Availability Statement: Data about geometries and meshes will not be available.
Acknowledgments: C. Vergara has been partially supported by the H2020-MSCA-ITN-2017, EU project 765374 "ROMSOC-Reduced Order Modelling, Simulation and Optimization of Coupled systems" and by the Italian research project MIUR PRIN17 2017AXL54F. "Modeling the heart across the scales: from cardiac cells to the whole organ".

Conflicts of Interest:
The authors declare no conflict of interest.