Efficient split-step schemes for fluid-structure interaction involving incompressible generalised Newtonian flows

Blood flow, dam or ship construction and numerous other problems in biomedical and general engineering involve incompressible flows interacting with elastic structures. Such interactions heavily influence the deformation and stress states which, in turn, affect the engineering design process. Therefore, any reliable model of such physical processes must consider the coupling of fluids and solids. However, complexity increases for non-Newtonian fluid models, as used, e.g., for blood or polymer flows. In these fluids, subtle differences in the local shear rate can have a drastic impact on the flow and hence on the coupled problem. There, existing (semi-)implicit solution strategies based on split-step or projection schemes for Newtonian fluids are not applicable, while extensions to non-Newtonian fluids can lead to substantial numerical overhead depending on the chosen fluid solver. To address these shortcomings, we present here a higher-order accurate, added-mass-stable fluid-structure interaction scheme centered around a split-step fluid solver. We compare several implicit and semi-implicit variants of the algorithm and verify convergence in space and time. Numerical examples show good performance in both benchmarks and an idealised setting of blood flow through an abdominal aortic aneurysm considering physiological parameters.


Introduction
Fluid-structure interaction (FSI) problems are characterised by a strong mutual dependence of fluid flow and structural deformation, exchanging momentum at the interface. Fluid forces acting on the solid cause deformation and induce strains, thereby influencing the stress state in the solid phase. A moving or deforming solid, in turn, alters the fluid flow domain and thus has a large impact on the flow quantities. Countless applications of FSI are found in science, engineering and biomedicine, ranging from airfoils or whole wind-turbines [1,2], bridge-decks [3], offshore engineering [4,5], insect flight [6] to blood flow through the circulatory system [7][8][9][10][11][12], human phonation [13] or respiration [14]. Consequently, the development of suitable models and solution procedures has been an active area of research over the past 40 years, which led to great advances in the field. Among the most popular numerical techniques to handle flows in moving domains -a central part of any FSI scheme -are arbitrary Lagrangian-Eulerian (ALE) [15][16][17][18][19][20][21][22], immersed boundary [23][24][25][26][27][28] and fictitious-domain methods [29][30][31][32], either tracking or capturing the motion of the fluid-structure interface.
Partitioned coupling schemes, on the other hand, alleviate the development of efficient preconditioners by iteratively enforcing the interface conditions. They allow a "reuse" of already well-advanced solution algorithms tailored to specific applications and ease the inclusion of complex physics or discretisation methods in the individual fields [56][57][58][59]. Using separate solvers and exchanging updates to the interface variables allows simpler software design, but shifts some of the intricacy to an outer coupling procedure. The added-mass effect [60,61] is notorious for severely impairing performance for certain parameter combinations using standard partitioned schemes such as the serial staggered method [62] and even implicit partitioned methods as demonstrated, e.g., in [63][64][65]. To counteract the hampered convergence behaviour in such cases, several remedies have been presented ranging from simple and effective Aitken relaxation [65], interface (quasi-)Newton or Newton-Krylov solvers [66][67][68][69][70], using Robin interface conditions [71][72][73] or artificial compressibility [74][75][76].
Nonetheless, performing several coupling iterations between fluid and solid phases remains costly. This was first alleviated by Fernández et al. [77], introducing the concept of semi-implicit FSI. Therein, the fluid flow problem is solved via a projection or split-step scheme, decoupling velocity and pressure unknowns. Additionally, the fluid domain deformation is extrapolated from previous time steps, which then allows for coupling only the fluid pressure and solid deformation implicitly. Based on this rationale of avoiding the implicit coupling of all components, a rapid development was seen in the following years (see, e.g., [20,[78][79][80][81][82][83][84][85]). In numerous challenging settings, these techniques were demonstrated to yield accurate and stable results while substantially increasing performance. Fully explicit treatment of interface conditions in FSI for problems with large added-mass effect was presented for shells (see, e.g., [85][86][87][88][89][90]) and also for three-dimensional continua [91][92][93][94][95]. Unfortunately, those fully explicit coupling schemes for bulk solids are as of now either limited to simple constitutive behaviour or only first-order accurate in time unless multiple correction steps are performed.
Although FSI has been a major research topic in the past decades, approaches specifically targeting generalised Newtonian fluid flow are still scarcely found in literature [88,[96][97][98][99]. And while partitioned FSI approaches and performant acceleration schemes do allow a straight-forward incorporation of complex constitutive laws, the promising concept of semi-implicit FSI builds upon methods decoupling fluid velocity and pressure, which then need to be capable of considering non-Newtonian models. Within such fluid solvers segregating velocity and pressure, necessary projection methods need to be free of nonphysical pressure boundary layers, which was not the case in the first generation of projection methods (see the work of Guermond et al. [100,101] for an excellent discussion on such schemes).
An alternative solution was presented in [102,103], replacing the continuity equation by a pressure Poisson equation (PPE) equipped with fully consistent boundary conditions. Further using extrapolated pressure and convective velocities in the fluid's balance of linear momentum, as often done also in mixed velocity-pressure formulations [104][105][106][107][108], allows to decouple velocity components and pressure unknowns completely. The PPE framework was recently extended to the generalised Newtonian case [109,110], enabling an efficient parallel solution using open-source finite element libraries and linear algebra packages as black boxes.
In this context, we present a novel split-step framework for partitioned FSI with incompressible, generalised Newtonian fluid flows and three-dimensional continua. Fluid velocity and pressure are completely decoupled using higher-order and possibly adaptive time-stepping and extrapolation formulae in a split-step scheme. This allows for equal-order, standard C 0 -continuous interpolation, which is a major advantage compared to similar projectionbased methods. Additionally, we derive fully consistent boundary and coupling conditions to preserve accuracy on the fluid-structure interface. Semi-implicit variants of the scheme are designed in an added-mass-stable way by implicitly coupling merely the solid displacement and fluid pressure. The remaining subproblems are treated in an explicit fashion, which enhances performance without degrading accuracy or stability. Moreover, we improve mass conservation through so-called divergence damping and achieve great flexibility concerning the rheological law, such that modifying the fluid material behaviour is as simple as exchanging the right-hand side of the viscosity projection step. All resulting linear systems are easily tackled using off-the-shelf black-box preconditioning techniques available as open-source scientific software, making the scheme an attractive alternative to available methods.
The remainder of this paper is organised as follows: Individual field equations are presented in Section 2, where time integration is carried out using higher-order extrapolation and backward-differentiation formulae of order m (BDF-m) for the fluid phase and the generalised-α scheme for the solid phase. In Section 3, fully implicit and semiimplicit variants of the coupling scheme are presented, enforcing consistent Robin conditions at the fluid-structure interface. The computational performance of the schemes is assessed in Section 4 with i) numerical tests of temporal and spatial convergence in two cases with analytic solutions [94,95], ii) the classical pressure pulse benchmark in three space dimensions (see, e.g., [34,96,111]), and iii) a final numerical experiment in the context of aortic blood flow, highlighting the potential of the presented approach for practical application.

Fluid and structure models in an ALE framework
The computational domain at time t, denoted by Ω t ⊂ R d with d = 2 or 3, is composed of the fluid and solid subdomains Ω t f and Ω t s with the moving interface Σ t = ∂Ω t f ∩ ∂Ω t s . Further, let us introduce for any function g(x, t) with x ∈ Ω t its counterpartĝ(x, t) = g(x, t) living in the reference configurationΩ. The transformations from reference to current domains are defined as with deformations d f and d s , following the arbitrary Lagrangian-Eulerian [16,112] and total Lagrangian approaches [113,114]. This gives rise to the deformation gradients F f and F s with respective Jacobians J f and J s :

Mesh update
The Lagrangian transformation L t is naturally defined by the deformation itself, but in the fluid domain, the mapping A t is constructed, e.g., via harmonic extension: with a possibly nonlinear stiffening parameter c = c 0 J f + J −1 f and suitably chosen c 0 [115,116]. This is only the simplest choice out of a wide variety of existing methods (cf. [115][116][117][118]), and may be further decomposed into equations in individual components of d f . Having introduced the basic setting, let us proceed with the formulation of balance equations in the individual subdomains.

Structure models
The balance of linear momentum in the reference configuration of the solid domainΩ s is expressed in terms of the solid displacement d s as with the solid's density ρ s and the first Piola-Kirchhoff stress tensor P, omitting body forces for brevity. Initial conditions and boundary conditions on the non-overlapping Dirichlet, Neumann and Robin boundary sections, denoted bŷ Γ D,s ,Γ N,s andΓ R,s , are given by withn s denoting the unit outward normal in the reference configuration and the Robin parameter η R s > 0. Constitutive equations linking stress and strain measures in the solid are herein formulated in terms of the second Piola-Kirchhoff stress tensor S, additionally using the relation P = F s S. Assuming isotropic, linear elastic material behavior leads to the St. Venant-Kirchhoff model with Lamé parameters expressed in terms of Young's modulus E s and Poisson ratio ν s . Further assuming ||∇d|| 1, one obtains for the case of linear elasticity Other constitutive relations of particular interest in biomedical engineering are those describing rubber-like materials such as arterial tissue (see, e.g., [113,[119][120][121][122]). Herein, a quasi-incompressible neo-Hookean model [121] P NH = µ s J −2/3 4 with the invariant I 1 = tr (C) and bulk modulus κ b = E/ [3(1 − 2ν)] is considered. This hyperelastic model is either used as is, choosing P = P NH , or together with contributions from dispersed collagen fibers. In the latter case, P is decomposed into P NH as defined in Equation (12) and an additional term, such that P is given as [119] Here, we introduce fiber parameters k 1 , k 2 and κ c , the tensors A i =m i ⊗m i with mean fiber directionsm i , the corresponding invariants I i = C : A i and G i for ease of notation [123]. Usually, in biomechanical applications, one defines the mean fiber directionsm i relative to circumferential (ê 1 ) and longitudinal (ê 2 ) directions of the vessel viâ with || · || denoting the Euclidean norm, leading to ±α c describing the deviation from circumferential vessel direction. The vector fieldsê 1 andê 2 are herein constructed with a two step procedure: First, two auxiliary scalar Laplace equations are solved with boundary conditions on the in-and outlets and the fluid-structure interface, such that the gradients of the resulting fields approximate the longitudinal (ê 2 ) and radial vessel orientation. Afterwards, the circumferential vessel direction (ê 1 ) is constructed using the radial and longitudinal vectors. Then, (14) gives mean fiber directions rotated by ±α c from the circumferential vessel direction into longitudinal direction. With these constitutive relations defined, differentiating between the individual stress tensors is not necessary, simply denoting the first Piola-Kirchhoff stress tensor by P for all of the above material laws.
The structural equations are discretised in time for t ∈ (0, T ], decomposing the interval into N t steps with size ∆t = t n+1 − t n , n = 0, ..., N t and using Newmark formulae [124] together with the generalised-α method [125], to obtain the time-discrete form of momentum balance in terms of where we introduce the shorthand-notations α m = 1 − α m and α f = 1 − α f . The nonlinear terms in Equation (17) are integrated via the generalised trapezoidal rule. Setting the parameters γ, β, α m and α f in the above time integration scheme determines accuracy and stability properties of the resulting scheme. Following [125], we achieve secondorder accuracy and unconditional stability for linear problems by choosing

5
The analysis in [125] was extended by Erlicher et al. [126], proving second-order accuracy and energy stability in the high-frequency range depending on the algorithmic parameter ρ ∞ even for nonlinear problems. The user-specified spectral radius in the high frequency limit ρ ∞ is used to specify α m and α f according to Table 1, resulting in the Newmark-β (N-β) [124], HHT-α [127], WBZ-α [128] or CH-α [125] schemes. Given the time-discrete form of the Table 1: Algorithmic parameters α m and α f in the generalised-α time integration scheme. momentum balance residual (17), one proceeds by employing Newton's method (see, e.g., [113,114,123] is fulfilled. Therein, the increment ∆d k s is the solution of the Jacobian system written in matrix-vector form which is based on the standard problem of finding δd k s ∈ [H 1 (Ω s )] d with δd k s |Γ D,s = 0, such that for all ϕ ∈ [H 1 (Ω s )] d , with ϕ|Γ D,s = 0 and ·, · Ω s or ·, · Γ denoting the L 2 (Ω s ) or L 2 (Γ) inner products. Moreover, ∂ ∂d s P denotes the directional derivative of P with respect to d s (see, e.g., [113,114]), which is omitted here for brevity. Known Neumann (t s ) and Robin (h s ) boundary data are plugged into the boundary terms arising from integrating the stress-divergence terms by parts. This completes the solution procedure for the nonlinear elastodynamics equations with Neumann and Robin boundary terms, various material models and generalised-α time integration.

Fluid models
The Navier-Stokes equations for incompressible flow, comprised of the linear momentum balance and continuity equations, reads in ALE form (assuming zero body forces): with the fluid velocity u f , stress tensor σ f , unit-outward normal in the current configuration n f , Robin parameter η R f > 0, ALE time-derivative ∂ ∂t u f A t and mesh velocity u m := ∂ ∂t d f , the latter two of which are connected via The fluid stress tensor σ f of a generalised Newtonian (or, as often called, quasi-Newtonian) fluid is given by being the symmetric gradient of u f , the fluid's pressure p f and dynamic viscosity µ f , which is most commonly expressed through a nonlinear map µ f = η(γ), η : R + → R * + , dependent on the fluid's shear rate defined aṡ Rheological models of great interest in biomedical and industrial applications may describe shear-thickening or shearthinning behaviour, as of particular importance when considering flows of polymer melts or blood, and may be put into the general form [129] Therein, η 0 and η ∞ denote viscosity limits and the fitting parameters κ f , λ f , a and n can be used to retrieve the power- or Carreau-Yasuda (κ f = 1) models, but also the standard Newtonian In the following, we extend the split-step scheme from [110] to a moving grid by replacing the standard velocitypressure form (22)-(27) by the following set of equations to advance the fluid velocity and pressure in time: with additional consistent boundary and initial conditions given as Altogether, this set of equations allow in their final form allow using C 0 -continuous finite element discretisations and decouple the balance of linear momentum and continuity equations, yielding standard discrete problems, for which off-the-shelf black-box preconditioners can be employed. This is a straight-forward extension of [110] to moving domains, which itself considers fluid flows on fixed grids and is based on [103] for the Newtonian case with open/traction boundary conditions and [102] for pure Dirichlet problems. Proof. See appendix.
This split-step scheme is not plagued by spurious pressure boundary layers, since the pressure is recovered from a fully consistent PPE rather than updated as in classical pressure-correction methods [100,102,103,110]. However, Liu [130] observed that stability can be improved significantly by performing a Leray projection, which is of particular importance when considering nonsmooth solutions. So, we aim to improve stability of the overall scheme and suppress accumulation of errors in mass conservation by solving the simple Poisson problem and updatingǔ f := u f − ∇ψ via projection. It is then easily verified that with any tangential vector τ f . Both u f and its (weakly) divergence-free counterpartǔ f converge at the same rates [101], making them equally valuable options from an accuracy point of view. However, it is also clear that one either settles for improved mass conservation, consideringǔ f , or chooses u f , fulfilling the boundary conditions exactly. There is an intermediate alternative available, though: As done by Liu [103] in the original scheme, we apply Leray projection on the past velocities only, effectively skipping the L 2 -projection step to obtainǔ f , and thereby fulfill the Dirichlet boundary conditions on the velocity u f exactly. This technique is also referred to as divergence damping [131,132] and has been shown to effectively reduce mass conservation errors and improve overall stability, while being cheaper than standard Leray projection and preserving boundary conditions on the velocity. Then, we construct the time-discrete weak form of the split-step scheme using BDF-m schemes of the form and higher-order accurate extrapolation exemplarily shown for u f with coefficients α m j and β m j according to Table 2, to effectively decouple balance of linear momentum and pressure Poisson equations. Since for generalised Newtonian fluids the viscosity depends on the shear rate, which itself is a function of the velocity gradient, standard Lagrangian finite elements cannot be applied in a straight-forward way due to increased regularity requirements on the velocity interpolant. Therefore, the viscosity µ f is introduced as an additional unknown and recovered through a simple to the scheme presented in [110], momentum balance and PPE steps are executed in reversed order, which is motivated by observations made, showing that updating the velocity with an implicitly coupled pressure increases stability in semi-implicit schemes of FSI as discussed later. So, we first project the pressure Dirichlet boundary conditions (33) and (34) on the respective boundary segments via for all ϕ ∈ H 1 (Ω t f ), with ϕ| Γ t p = 0 to apply divergence suppression on the current time step's velocity to be used in the next time step's momentum balance equation.
In summary, it is thus possible to construct a weak form containing only first-order derivatives, rendering our beloved C 0 -continuous, standard Lagrangian finite elements applicable to the problem at hand. In fact, we might even employ equal-order finite element pairs for velocity and pressure as already pointed out. The presented weak forms contain a generous set of boundary conditions, which will in the coupled FSI-problem (partly) depend on the solid subproblem's solution as shall be seen next.

The coupled FSI problem
The strong form of the FSI problem incorporating mesh, fluid and solid subproblems as discussed in Section 2 including only interface conditions for brevity reads where (51)-(53) enforce the continuity of displacements, velocities and tractions on the fluid-structure interface.
Note that the continuity of tractions in Equation (53) is formulated on Σ t , using the Cauchy stresses and the current configuration's normal vectors n f , but is easily rewritten as to enforce balance of tractions in the reference configuration. The Robin-Robin (RR) coupling conditions (see, e.g., [72,86] or [79] in a projection-based semi-implicit scheme), linearly combine the interface conditions enforcing continuity of velocities and normal tractions, i.e., (52) and (53) or (54), yielding in the respective configurations with Robin parameters η R f , η R s > 0. This type of interface condition leads to a coupling algorithm with good convergence properties even in the case of high added-mass effects (cf. [72,73,79,134]), which is of particular importance in biomedical applications [10,60,61]. As the basic algorithm, we perform an implicit single-loop coupling scheme (see, e.g., [10,[135][136][137]), which is executed until convergence criteria of the form are fulfilled. In the following, we will denote the last iterates by a superscript k, and the newly computed iterate by a superscript k + 1. Moreover, we directly present the RR scheme, which is obtained inserting interface conditions into the Robin terms of the respective subproblems. So, at each time step n, given the solutions from previous time steps 2. Extrapolation/initial guess: Compute d s , µ f , u f and p f based on old time step solutions via (40) and set d k s = d s , u k f = u f and p k f = p f as initial guess. 3. Implicit coupling loop: WHILE not converged according to Equation (57) DO and linearised stiffening parameter c(d n f ). (b) Mesh velocity update: Compute u k+1 m via the BDF-m formula (39).
(d) Pressure boundary projection: Update the pressure Dirichlet condition by projecting ζ k+1 on Γ t N, f using (e) Pressure Poisson step: holds for all ϕ ∈ H 1 (Ω t f ) with ϕ| Γ t N, f = 0, using the last solid iterate d k s to computed Note here that the treatment of Robin boundary conditions in the momentum balance and PPE steps is not simply assigning Robin boundary conditions to the fluid subproblem, but rather enforcing Robin conditions on Σ t in the fluid momentum equation and treating the interface as a Dirichlet boundary for all steps related to the fluid pressure. This combination is equivalent to the strategy adopted by [79] and is herein solely based on numerical observations. The sequence of fluid steps and viscosity projection turned out to be the most stable choice when confronted with large time steps and sudden jumps in fluid pressure boundary conditions as present in the pressure pulse benchmark in Section 4.3.
The RR coupling algorithm as introduced above includes the standard Dirichlet-Neumann (DN) coupling scheme in the asymptotic limit, when η R f → ∞ and η R s = 0. In the discrete setting, however, η R f has to be assigned a bounded value, which motivates including the interface Dirichlet condition on the fluid velocity in a more direct way. Thus, in the DN case, the interface is treated as part of the fluid Dirichlet boundary Σ t ⊂ Γ t D, f together with setting η R s = 0, effectively leading to small changes in the function space definitions only, but not introducing any additional terms.
To counteract decreased convergence for high added-mass effects, Aitken's acceleration is applied to relax the discrete solution vectorx k+1 in iteration k with a recursively defined ω k [65] x
Another option to increase efficiency drastically is to give up on fully implicit coupling of all the involved subproblems, but rather settling for a semi-implicit variant of the scheme. This option is directly accessible having formulated the fully implicit algorithms by simply moving the mesh and fluid momentum subproblems and the viscosity projection out of the coupling loop, similar to the methods proposed by [20,77,[79][80][81][82]84] for Newtonian fluids and [99] considering a non-Newtonian, viscoelastic fluid. Interestingly, the sequence of substeps in the fluid phase had to be changed in order to yield satisfying results. First, the fluid pressure and solid displacement are implicitly coupled, before explicitly treating the fluid balance of linear momentum and viscosity projection rather than performing the update on u f before the implicit loop. This way, the dependence on an extrapolated, non-coupled pressure is eliminated, which results in improved robustness of the scheme. An additional tuning possibility is available via the convergence criterion in the solid's Newton scheme, where one may tweak the tolerance N or even exit after a fixed number of steps in the nonlinear solver in the spirit of [137,138]. Using suitable higher-order extrapolation schemes, temporal accuracy is preserved, while the fully implicit coupling of pressure and structural displacements is sufficient to obtain a stable method as numerically observed and proven for simplified model problems (cf. [77,81,136,139]). In a nutshell, the following distinct features and benefits arise in the proposed scheme when compared to related methods: (i) Higher-order and possibly adaptive time-stepping schemes are available based on standard time integration and extrapolation formulae.
(ii) Equal-order finite element pairs can be employed, which would be unstable in the classical coupled velocitypressure formulation of the Navier-Stokes equations.
(iii) Exchanging the rheological model of the generalised Newtonian fluid is as simple as changing the right-hand side of the projection step.
(iv) The semi-implicit design reduces computing times tremendously, while preserving stability properties and accuracy. Only structural displacements and fluid pressure are iteratively coupled in each time step.
(v) Robin interface conditions might improve convergence even for high added-mass effects when suitable parameters are available and standard acceleration methods are directly applicable.
(vi) Divergence suppression avoids the accumulation of errors in mass conservation and neither spoils interface conditions nor requires a velocity projection step.
(vii) All linear systems can be effectively tackled using off-the-shelf black-box preconditioning techniques available as open-source scientific software.

Computational results
This section is devoted to the thorough testing of the presented schemes in terms of accuracy and robustness as well as critically comparing their individual performance. All of the showcased results were obtained with the finite element toolbox deal.II [140], solving each of the arising linear systems involved in the FSI-algorithm (58)- (63) iteratively. We employ algebraic multigrid methods provided by Trilinos' ML package [141] for preconditioning each linear solve. A preconditioned conjugate gradient method is used for the mass matrices in the viscosity projection step (60), the pressure Dirichlet data projection (61) and also for the Poisson problems, i.e., the PPE (62) and the mesh motion equation (59). The linear systems corresponding to fluid and solid momentum balance equations are solved adopting a flexible generalised minimal residual method.
The studied test cases are two analytical solutions taken from [94,95] to demonstrate convergence rates numerically, a classical benchmark of a pressure pulse travelling a straight pipe in three spatial dimensions (see, e.g., [34,96,111]) and, finally, we study the flow through an idealised abdominal aortic aneurysm to showcase performance in a practically relevant setting.

Analytical solution: rectangular piston
An analytical solution is taken from [94,95], which oscillates vertically with amplitude a and frequency ω. The structural displacements are thus given bŷ The fluid's vertical velocity resulting from the continuity equation is solely dependent on time, u f,2 (t) = ∂ ∂td Σ,2 (t), and the pressure is given by We refer to the original publications [94,95] for details on the derivation and proceed in defining problem parameters.
The fluid density and dynamic viscosity are set to ρ f = 1 kg/m 3  When refining the time step, expected convergence rates are observed in all primary variables when using the Newmark-β scheme, as can be seen in Figure 2. However, when introducing numerical high-frequency dissipation for the HHT-α and WBZ-α, whereas the so-called asymptotic annihilation case, i.e., ρ ∞ = 0, leads to a clearly linear convergence rate in p f for WBZ-α and CH-α methods as can be seen in Figure 4 where H s = 1, β = 0 and γµ max = 0.01 were chosen, the latter of which was found uncritical in this example. Considering the substantial increase in efficiency, the semi-implicit schemes are very attractive for practical applications as indicated at various places in the literature [20,77,[79][80][81][82]84].

Analytical solution: circular piston
Another analytical solution is employed to confirm the expected convergence rates in space. A circular piston as depicted in Figure 1b, where the structure occupies the region fromr = 0 to the interface atr Σ = r Σ (t = 0) in the  reference configuration, pulsates in radial direction, driving the fluid. The fluid adheres to the piston's periodic motion and since the time-dependent fluid domain ranging from r ≥ r Σ (t) to r = R changes in volume, the fluid freely exits and enters the computational domain due to incompressibility. For a detailed derivation see [94,95], whereas herein we conveniently express the solution in terms of the radial component of structure displacement d s,r (r, t) := βJ 1 ωr c p sin(ωt) , such that r Σ (t) = r 0 Σ +d s,r (r 0 Σ , t) with frequency ω, a parameter β scaling the amplitude, initial piston radius r 0 Σ , c p as defined in Equation (66) and J 1 denoting the Bessel function of the first kind and order one, which gives the fluid's radial velocity component and the fluid pressure where J 1 denotes the first derivative of J 1 . The piston with initial radius r 0 Σ = 0.5 pulsates with frequency ω = π and β = 0.1 is set. Fluid's density and dynamic viscosity are considered as ρ f = 1 g/m 3  The experimental orders of convergence are reported in Figure 6, allowing for a direct comparison of the errors obtained with Q 2 /Q 1 and Q 1 /Q 1 elements and various coupling schemes, where the Robin parameters are again set according to (67). The deliberately chosen parameters result in directly recovering the exact solution of d s up to the specified tolerance/time integration error, but allow easily measuring experimental convergence rates of fluid velocity and pressure. The fluid velocities converge at a rate of eoc = 2 for the fully implicit and semi-implicit schemes with almost identical errors obtained. This is exactly as expected, even for the Q 2 /Q 1 finite element pairing (cf. [102,103,110]). Pressure rates of order 1 are observed for linear pressure interpolation and a slightly higher rate of eoc ≈ 1.5 when using Q 2 /Q 1 interpolation. These results are optimal for the considered split-step scheme in the fluid pressure, where one might expect convergence rates of order 1 using (bi-)linear elements in the observed norm.
However, the measured rates of ≈ 1.5 for the Q 2 /Q 1 pair are lower than what one would hope for, judging from the basic split-step scheme, which gives rates of 2 in the pressure norm observed here. These higher rates are found initially, but decrease under refinement. Additionally, the solid subproblem yielded optimal convergence rates when tested using different parameter settings, which is omitted for the sake of brevity. The sole reason for not choosing those different parameter settings altogether is the inherent difficulty in demonstrating all convergence rates at the same time, which was not found possible here due to the problem setup and relative sizes of physical quantities and 19 tolerance choices.    [111,[142][143][144][145][146][147]. Concerning the material parameters, we set the fluid and solid densities to ρ s = 1200 kg/m 3 and ρ f = 1060 kg/m 3 , which triggers strong added-mass effects. Further, we consider a Newtonian fluid with a viscosity of µ f = 3.5 mPa s, which needs to be suitably chosen (see, e.g., [96,148,149]) for comparison to the more general Carreau model with η 0 = 56 mPa s, η ∞ = 3.45 mPa s, λ f = 3.313 s and n = 0.3568 taken from [150]. For the solid phase, linear elasticity (11) Figure 8), but negligible effects on displacements and pressure. Additionally, we observe "backflow" instabilities as a consequence of fluid entering over the inlet Neumann boundary due to mass conservation. Despite those oscillations close to the inlet, the solver still converges. Possible remedies to stabilise such effects are backflow stabilisation or similar techniques [111,142,146,[154][155][156], which are not considered at this point.
Let us note here that the presented solutions can merely indicate that suitable parameters and material laws themselves are central aspects in the description of the system behavior. This, however, does not lie within the scope of this contribution, since the aim here is simply to demonstrate versatility by easily switching between constitutive equations of particular interest depending on the available data or application. Within this work, we present this problem setup to assess computational performance of the various coupling schemes presented. From now on, we only consider linear elasticity with a Newtonian fluid and the HGO model together with the Carreau law.
Inspecting the accumulated FSI iterations depicted in Figure 10, we see that the total iteration counts are decreasing with the semi-implicit schemes. Interestingly, treating the fluid mesh motion and momentum balance equations     in [72,79]. In those papers, the Robin condition on the fluid phase was found to substantially decrease the number of needed coupling steps, which does not seem to transfer to the present split-step approach. To diminish a possible influence of the Robin parameter, tests are conducted in the linear elasticity/Newtonian fluid case varying the Robin parameters in the fluid momentum balance and PPE. The latter option modifies Equations (61) and (62) to include Robin interface conditions as well, but was found diverging for any η R f . Scaling the Robin parameter η R f computed via Equation (67) in the fluid momentum balance by some α R , iteration counts stay almost identical for α R large enough, as can be seen in Figure 11. Choosing α R → 0 and α R → ∞ correspond to Neumann-Neumann and the standard Dirichlet-Neumann coupling schemes, first of which needs more than the preset maximum of 150 steps to converge leading to divergence in the fifth time step.
The Robin parameter η R s in the solid's balance of linear momentum has two distinct interpretations depending on the coupling scheme applied: in implicit schemes, a standard Robin-Robin coupling is recovered, possibly mildly increasing efficiency as reported by Badia et al. [72]. In our semi-implicit scheme, however, it is recalled that the interface Robin term is responsible for information transfer from fluid to solid phase (with n indicating the time step, k the iterate in the FSI coupling scheme and l the iterate in Newton's method) as which as η R s → 0 reduces to a Neumann condition on the interface, but as η R s → ∞ smoothly transitions to an explicit coupling scheme based on the extrapolated velocity u k f . Despite the latter greatly decreasing the number of coupling iterations depending on the specific choice for η R s , one also observes decreased temporal stability. Therefore, η R f = 10 2 ρ s ∆t and η R s = 10 −4 ρ f ∆t are chosen for all of the RR results presented in Figure 10. Before moving on to a final, more challenging example, let us summarise the insights gained from this simple pressure pulse benchmark: As expected, semi-implicit treatment of both the fluid mesh motion and fluid balance of linear momentum is found sufficiently accurate, while the constitutive equations considered for fluid and solid properties, only mild improvements are seen in some cases, but instabilities arise from unsuitable choices. The semiimplicit Robin-Robin or Dirichlet-Robin schemes transition to fully explicit schemes as η R s → ∞, which are found rather unstable in first numerical tests. Consequently, the method of choice is the SIDN algorithm, which delivered low iteration counts, stable solutions and requires less data transfer than the Robin variants.

Flow through an idealised aneurysm
In this final numerical example, we consider the flow through an idealised abdominal aortic aneurysm (AAA) using physiological parameter sets to challenge the devised framework. The setup chosen is similar to [157], combining a prototypical geometry as recommended by surgeons [158,159] with flow data from [160]. Comparable settings in a biomedical context with benchmark character were presented by Turek et al. [161] and Balzani et al. [162], while patient-specific geometries were considered, e.g., in [10,53,163].
The AAA geometry with a length of l = 20 cm and with inner/lumen radius of r i = 1 cm at the in-and outlet is discretised using ≈ 1.3 × 10 5 trilinear elements as depicted in Figure 12, resulting in ≈ 1.3 × 10 5 nodes. These radii are expanded to a maximum of ≈ 6.5 cm in lateral direction and ≈ 5.    [122,143,145].
The fluid inlet condition is computed from periodic mean inlet velocityū in as shown in Figure 13b, prescribing the normal inlet velocity component u 1 in terms of the distance from the circular inlet center r: and a factor of 2 to match the volumetric flow rate computed from the mean velocityū in with a parabolic velocity profile. The outlet pressure is approximated via a three-element Windkessel model (see, e.g., [122,142,146,147]) which represents the characteristics of the excluded downstream vasculature. Therein, the flow over the outlet is denoted as Q(u f ) and parameters for the capacitance C = 1.25 × 10 −9 m 4 s 2 /kg and the proximal and distal resistances are specified as R p = 266.66 × 10 5 kg/m 4 s and R d = 6.8 × 10 8 kg/m 4 s, respectively. The outlet pressure is then weakly enforced through the standard traction boundary integral term, setting also including backflow stabilisation according to [164], where (u f · n f ) − = −u f · n f for u f · n f ≤ 0 and 0 otherwise.
Moreover, Galerkin least-squares stabilisation [165] is added to the fluid momentum balance equation to counteract dominant convective terms.
Concerning material properties, we keep the physiological parameters as set in the pressure pulse benchmark, i.e., adventitia layers of the aorta [151,152]. The prestress present in the tissue in a geometry reconstructed from medical image data is not considered in this setup. However, we can interpret the present mesh as the zero-stress geometry being the result of the prestress strategy "deflating" the model as in [166,167]. This slight discrepancy is simply ignored, since we do not investigate prestress effects. Altogether, three pulses are considered, i.e., t ∈ (0, 3], using a uniform time step of ∆t = 1 ms in the second-order accurate scheme, i.e., BDF-2 and linear extrapolation for linearisation. The WBZ−α and CH−α time integrators with ρ ∞ = 0 are selected, since they were again found to be more robust when using the Windkessel model. Aitken's relaxation is initiated with ω 0 = 0.001, coupling the pressure Poisson (PPE) and solid momentum subproblems in the semi-implicit Dirichlet-Neumann (SIDN) approach until reaching abs = 10 −6 or rel = 5 × 10 −4 , while a relative tolerance in the Newton solver (19) of N = 10 −5 was selected conservatively.
Looking at solution snapshots at three distinct time instances at t = 2.2, 2.4 and 2.6 s, we observe strong recirculations in the AAA, especially after the rapid drop of the inflow velocity. In Figure 14, selected streamlines indicate areas of recirculatory flow, leading to large gradients in the velocity fields and hence to a large shear rate, which ultimately results in strong gradients in the viscosity field. Further, the maximal and minimal viscosities observed at any point in time span almost the entire range from η ∞ to η 0 . The pressure in the lumen is spatially rather uniform, but rapidly changes in time due to the Windkessel model setting the pressure level, which dominates the deformation rather than the velocity acting on the vessel. Inspecting the time evolution of the displacement norm and fluid pressure at the apex point atx = (0.1, 0.039, 0) T depicted in Figure 15, we see that the periodic state is not yet reached. This is due to the Windkessel model applied on the outflow boundary, which only gradually increases the pressure level in the AAA. Observing the displacement norm and fluid pressure only, effects of the more advanced solid constitutive models cannot be investigated, especially since the displacement remains rather small. The parameter choice for the tissue, however, leads to a difference easily observable with the naked eye.
Regarding the difference when applying Newtonian or Carreau rheological models, Figure 15 is not sufficient, since the relevant quantities are neither the displacement nor the pressure, but rather the time averaged shear rate and shear stress. After defining the time averagef (x) of a function f (x, t) over a period T p bȳ calculateγ setting i = T p = 1, while the shear stress is computed following John et al. [168]: Inspecting these time averaged quantities in the third cycle, as shown in Figure 16, a striking difference is observed as expected. Nonetheless, within this work we focus on the FSI solver rather than the phenomenological influence of different model decisions, but still want to demonstrate the versatility of the framework.
Having thoroughly inspected the solution itself, let us now turn our attention to the numerical aspects of this test case. Motivated by the previous example, the semi-implicit Dirichlet-Neumann (SIDN) coupling scheme is considered. Throughout the entire simulation time, less than 30 steps coupling solid momentum balance and PPE are needed, as shown in Figure 17. Only a slight dependence on the inflow profile and pressure level in the geometry are observed. Interestingly, accumulated FSI iteration counts show that NH and HGO models lead to fewer FSI iterations, which is most likely triggered by higher pressure levels due to a stiffer material response caused by the selected (higher) Young's modulus. For comparison, the implicit Dirichlet-Neumann (IDN) scheme and a scheme treating the mesh motion explicitly, leading to a geometry explicit (GEDN) approach are also included in Figure 17b, showing a decrease of up to 31% in FSI coupling steps needed when using the semi-implicit variant.
Regarding the iteration counts in the subproblems, consider first systems solved only once per time step. Solving              pay off using the SIDN scheme, which is why we refer to [110] for further details.
For the overall performance, the steps within the semi-implicit coupling loop, namely the PPE and solid momentum balance solves are dominant in terms of computational effort. Thus, the good performance of the linear solvers as depicted in Figure 19 is crucial. We observe little dependence on the fluid velocity and pressure level, given the relatively small time step size. Additionally, nonlinear convergence in the Newton solver is reached in any time step within 3 iterations, which is also a result of the small time step size and the quadratic extrapolation used as an initial guess. Comparing the constitutive models, the linear solvers show equal behaviour independent of the fluid model employed. The solid model, however, heavily influences the iteration count in the linear system solve, resulting from the more complex terms and matrices when considering NH or HGO models.   compare the absolute and relative time spent per time step using different constitutive models. As can be seen from Table 3, the number of FSI coupling steps needed for convergence are slightly higher for stiffer material parameters, while the viscosity projection step is negligible in terms of computational effort. The solid constitutive models, however, do have a strong influence on the time spent per time step. This additional effort is caused by both the increased complexity of element integration and ≈ 2.5× more iterations in the linear system solve (see also Figure 19).
The SIDN variant outperforms the GEDN and IDN schemes being at least twice as fast overall, needing less FSI iterations and less time per coupling step. Interestingly, the number of coupling steps increases when using the GEDN scheme, but individual steps and the overall time interval are completed faster compared to the IDN scheme, since the mesh motion equation is only solved once per time step.
Timings for the computationally relevant steps within the SIDN coupling scheme are further listed in Table 4, where we notice the increased (relative) time spent in those steps. Also, this demonstrates how little time is spent on the semi-implicit steps, being the mesh motion equation (d f ), fluid momentum balance (u f ), Leray projection (ψ) and viscosity projection, the last one of which is not even listed since the time spent is less than 0.5%. The pressure boundary projection (ζ) is performed in each coupling step, but still does not lead to a major computational load.
Solving the PPE (p f ) and solid momentum balance (d s ) has the largest influence on the overall time spent in the SIDN scheme, where one may further optimise element integration routines and linear solvers. Nonetheless, iteration counts are nicely bounded, indicating good performance of the chosen setup. Parallel scalability of the overall algorithm is not investigated at this point, but aspect of future investigations.

Concluding remarks
Within this work, we presented a family of coupling schemes involving incompressible viscous flows interacting with three-dimensional solids. The standard Navier-Stokes equations were solved by decoupling velocity and pressure via an equivalent time-splitting or split-step scheme in arbitrary Lagrangian-Eulerian formulation, based on a pressure Poisson equation with consistent boundary conditions. Consequently, substantial parts of the algorithm can be treated semi-implicitly in an added-mass-stable way without degrading accuracy. The mesh motion equation, fluid momentum balance equation, Leray projection (included for improved mass conservation) and viscosity projection are only solved once per time step. Generalised Newtonian rheological laws are exchanged effortlessly, while also allowing for equal order, C 0 -continuous interpolation with Lagrangian finite elements. Moreover, we considered a general setup including various constitutive models applied for the solid, demonstrating the flexibility of the framework and robustness of the coupling algorithm. The accuracy of the scheme is assessed comparing to analytical solutions in space, applying Q 2 /Q 1 and Q 1 /Q 1 finite element pairs, and in time, combining BDF-2 with generalised-α time integration schemes. The framework is then further tested in a benchmark example in the context of arterial blood flow and an idealised abdominal aortic aneurysm, highlighting its great performance even with black-box solvers and preconditioners.
At this point, some open questions remain to be further analysed in the future: First, we note that owing to the extrapolation of the fluid pressure in time needed in the split-step scheme for the fluid problem, the restrictions on the time step size grow stronger with increasing the temporal order of accuracy. Thus, using high-order time stepping schemes beyond the presented second-order time integrators might reduce the maximum allowed step size drastically.
Moreover, when using even higher polynomial orders, the splitting error of the scheme itself might dominate spatial higher-order accuracy as is observed, e.g., for projection-schemes. Second, we want to mention that the Robin variants of the coupling scheme did not perform as expected in our studies. Depending on the parameters, the considered Robin variants can substantially decrease the number of coupling steps needed, but the correct settings remained elusive, demanding further investigation to not spoil overall stability. That being said, the acceleration scheme used herein is Aitken's relaxation for simplicity, but much more advanced schemes can be applied straight-forwardly. In this regard, first tests with the interface quasi-Newton inverse least-squares method [56,69] are very promising.
Future and ongoing work is centered around scalability tests and algorithmic optimisation to further reduce runtimes. From a modelling perspective, other constitutive equations and mixed/hybrid formulations for the solid phase will be included, harnessing the partitioned design. Moreover, we are testing semi-implicit coupling to thrombus formation models, which is a central element in various cardiovascular conditions.
where we can further use relations similar to (69) and (68) to get where we use (69) and (70), but without inserting ∇ · u f = 0, to finally arrive at being a heat equation in the new variable Φ := ∇ · u f in ALE form. Neumann conditions for (71) are obtained by taking the scalar product of (30) and n f and adding the result to (35), which gives Further using once again (69) and (68), we end up with Dirichlet conditions for (71) are constructed from (33) and (34), which after inserting the definition of t f and h f , respectively, directly result in Given these zero Neumann and Dirichlet conditions for Φ, and assuming the geometric conservation law is fulfilled, i.e., a constant quantity is preserved on the moving grid (cf. [169][170][171][172]) we obtain Φ ≡ 0 as the only admissible solution to (71)- (73). Consequently, the modified system also inherently enforces incompressibility and we conclude the proof by quickly noting that equivalence of (30) and (22) is easily shown using (68) again.