Validation of a non‐conforming monolithic fluid‐structure interaction method using phase‐contrast MRI

Abstract This paper details the validation of a non‐conforming arbitrary Lagrangian‐Eulerian fluid‐structure interaction technique using a recently developed experimental 3D fluid‐structure interaction benchmark problem. Numerical experiments for steady and transient test cases of the benchmark were conducted employing an inf‐sup stable and a general Galerkin scheme. The performance of both schemes is assessed. Spatial refinement with three mesh refinement levels and fluid domain truncation with two fluid domain lengths are studied as well as employing a sequence of increasing time step sizes for steady‐state cases. How quickly an approximate steady‐state or periodic steady‐state is reached is investigated and quantified based on error norm computations. Comparison of numerical results with experimental phase‐contrast magnetic resonance imaging data shows very good overall agreement including governing of flow patterns observed in the experiment.


INTRODUCTION
Many industrial and engineering problems, for example, in aeronautics, 1 power generation, 2 defense, 3 and biomedical engineering, 4,5 involve complex multiphysics phenomena, such as the interaction between fluids and solids. In the field of biomedical engineering, collaborative work of researchers, modelers, physicians, medical imaging technicians, and others becomes increasingly important to ultimately provide patient-specific models [6][7][8][9] for therapy planning.
Innumerable mathematical modeling techniques and numerical solution methods for fluid-structure interaction (FSI) problems have been proposed to date and can be classified based on the nature of the underlying algorithmic approach and solution strategy. Classification based on the nature of underlying meshes distinguishes immersed methods and moving domain methods. Immersed methods employ a fixed background mesh and a Eulerian description for the motion of the fluid over the flow domain, whereas the solid motion and deformation are described within a Lagrangian coordinate frame on an embedded mesh. The presence of the solid is accounted for either via adding a body force term to the fluid equations (immersed boundary method 10 ) to constrain local flow with a similar effect as the no-slip condition at fluid-solid boundaries, or by coupling the fluid and solid equations by introducing a Lagrange multiplier at the fluid-solid interface (ficticious domain method 11,12 ). Immersed methods are particularly well suited for FSI problems involving large structural deformation or solids moving through flow domains. On the other hand, moving domain methods enable inherent interface-tracking by using an arbitrary Lagrangian-Eulerian (ALE) coordinate frame for the fluid domain 13 and constraining the fluid and solid equations by requiring equal but opposite tractions and the no-slip condition at the fluid-solid interface. Although, the ability of this method to deal with large structural deformation is often noted as a limitation, and remeshing steps might become necessary to maintain mesh quality. On the other hand, Lagrangian meshless methods, such as the study of Idelsohn et al, 14 avoid the requirement for remeshing. Fluid-structure interaction methods can also be classified based on the solution strategy, for example, monolithic/partitioned approach (global assembly and solution of a single matrix system vs solution of fluid and solid subsystems with exchange of boundary values) and explicit/implicit discretization and time integration, where a given choice may impact computational cost, accuracy, and stability.
In this paper, we consider a monolithic ALE FSI technique that is able to use non-conforming meshes at the interface, 5,[15][16][17] such that meshes can be designed based on the requirements of the physics of the coupled subsystems leading to improved accuracy and to decreased computational cost (by avoiding underrefinement and overrefinement, respectively). Coupling of subdomain equations is achieved via introduction of an additional coupling domain and enforcing interface constraints by means of a Lagrange multiplier variable. The method has been studied regarding stability and convergence 15,18 and successfully applied to various biomedical engineering problems, such as the simulation of whole-heart and left ventricular mechanics. [15][16][17]19,20 Previously, the method has been used for coupling the non-conservative ALE Navier-Stokes equations and the governing equations for quasi-static/transient finite elasticity. It has been extended recently to enable modeling of turbulent flow phenomena by a stabilized cG(1)cG(1) scheme 21 to extend the use of the method over a larger range of Reynolds numbers. Besides the various biomedical engineering applications, 5,17,20 the method has been assessed and verified using test problems; however, it was not validated in any previous work. Thus, validation of the method will be the focus of this work as well as validation of using the cG(1)cG(1) scheme within the Lagrange multiplier-based coupling method. In this work, we focus on the validation of the method as well as the comparative performance of an inf-sup stable scheme and the cG(1)cG(1) approach.
Verification and validation are important to confirm fidelity and assess the capabilities of established and newly proposed mathematical models and numerical algorithms. Thus, standard numerical benchmark problems [22][23][24][25][26] and FSI experiments [27][28][29][30][31][32][33][34][35] have been developed over the last decades and found widespread use. In this tradition, a recently developed 3D FSI experiment 36,37 introduces two new challenging benchmark test cases that involve steady and periodic interaction between a moderately viscous incompressible fluid and an incompressible nonlinear solid in a 3D setting. 36,37 With key aspects (for example, flow regime, material parameters, and mechanical properties) of the experiment being in line with those given in typical translational biomedical engineering applications (such as simulation of left ventricular mechanics under support of a left ventricle assist device 5 ), the benchmark is considered in this paper for validation of the non-conforming monolithic FSI method. An inf-sup stable (iss) and the cG(1)cG(1) scheme are considered for the fluid model. Numerical predictions of both methods are compared and contrasted with experimental phase-contrast magnetic resonance imaging (MRI) data. Spatial refinement, fluid domain truncation, and convergence to (periodic) steady-state are studied in this paper. Further, the employed coupling technique is assessed, and performance of both schemes is investigated regarding computational cost and prediction of steady and dynamic behavior (flow patterns, deformation of solid, summed forces at fluid / solid boundary, and others).
In the following, the numerical solution procedure is based on the use of a non-conforming monolithic ALE FSI technique. The details of this method are outlined in Section 2. Further, in Section 3, we present results obtained for both benchmark test cases and validate our results via comparison with experimental results. Further, aspects such as spatial refinement and early truncation of the fluid domain are assessed. Finally, we discuss the quality of our numerical results and conclude with possible future improvements in Section 4.

METHODOLOGY
In this section, a brief overview of the 3D FSI experiment is given. The model problem and the respective model equations are detailed. The weak form for the FSI problem within a finite element formulation using Lagrange multipliers is given using an inf-sup stable and stabilized scheme. The section concludes with details about the numerical solution and the definition of error norms for comparison of numerical and experimental data.

3D FSI experiment
The flow domain in the relevant section (that is, where FSI phenomena occur) of the experiment features two inlets (diameter ⊘21.9 mm) that merge smoothly into a single outlet (diameter ⊘76.2 mm), with a solid (volume 11 × 2 × 65 mm 3 ) attached to the wall in the merging section ( Figure 1). A right-handed Cartesian coordinate system is used 37 with the origin chosen to be at the center of the attachment point of the solid to the wall. The selected pump rates create steady (Phase I) and periodic (Phase II) inflow. This yields steady and time-dependent periodic interaction between a moderately viscous incompressible fluid (aqueous glycerol solution) and an incompressible nonlinear solid (silicone material).
Experimental data were acquired using MRI techniques and are available for comparison with numerical results. The data include the geometry under zero inflow conditions (for example, the deformed state of the structure with maximum deflection of 29.50 mm and 25.65 mm for Phases I and II, respectively), inflow boundary condition data and time-resolved flow and deformation fields.
For Phase I, parabolic inflow profilesṽ I f were defined on Γ I f with peak value [0, 0,ṽ z ] T  and v x =ṽ y = 0 ∀ t ∈ [0, T]. B, Phase II: recorded peak velocityṽ i (i∈{x,y,z}) for parabolic inflow and data fit. We note, thatṽ y | y<0 = 0 37 with a smooth increase over 0.5 s (Figure 2A) for the upper (y > 0) and lower (y < 0) inlet. For Phase II, parabolic profilesṽ I f were defined on Γ I f with peak values [ṽ x ,ṽ y ,ṽ z ] T fit to the experimental data, see Figure 2B. Fluid material parameters and solid density are defined in the study of Hessenthaler et al 37 and collected in Table 1, which also contains the solid material parameter for an isotropic incompressible Neo-Hookean material law. It was shown in the study of Hessenthaler et al 37 that this solid constitutive model can be used to model the mechanical response of the silicone material and that it is able to represent test data obtained from a uniaxial tensile load-displacement test.
However, the silicone material used in the FSI experiment undergoes a continuous curing process, 37 such that the solid material model was calibrated using the zero inflow displacement data recorded for both test cases.

Reference frames
The 3D FSI domain arising from the definition of the 3D FSI experiment consists of the fluid and solid subdomains, Ω f ⊂ R 3 × I and Ω s ⊂ R 3 × I (with I = [0, T]), and its respective subdomain boundaries, f and s . The fluid domain boundary Figure 1. Similarly, the solid domain boundary Γ s = Γ W s ∪ Γ C s is partitioned into subdomains for the wall Γ W s and the coupling surface Γ C s . A third coordinate frame on the coupling domain Ω = Γ C f = Γ C s is introduced to enforce coupling of the fluid and solid equations. Further, Ω 0 i denotes reference domain i (with respective boundary Γ 0 i ), and n i denotes the outward boundary normal.
To relate reference domains, Ω 0 i , to moving domains, Ω i (t), bijective mappings are introduced, 38,39 where w f is the fluid domain velocity (referred to as ALE velocity), u i is the deformation of the domain (generally, u i (·,0) ≠ 0 such that Ω 0 i ≠ Ω i (0) and Γ 0 i ≠ Γ i (0)), and X is the spatial coordinate of a point in the reference domain. More precisely, Equations 2 and 3 denote an ALE mapping for the fluid domain and a Lagrangian mapping for the solid domain, repectively, such that conservation laws can be equivalently formulated on the reference or moving domain. Using the ALE and Lagrangian mapping, we can link a function given on reference domain Ω 0 f to its counterpart on the moving domain Ω f , ie, where x is the spatial coordinate of a point in the moving domain and likewise for functions given on Ω 0 s and Ω s . Further, the Jacobian mapping of the domain displacement is given as is the deformation gradient tensor. In the following, the hat notation in Equation 4 will be omitted; t denotes the temporal derivative with respect to a fixed point in the reference domain (eg, the studies of Nordsletten et al, 16,19 and Formaggia and Nobile 40 ), and ∇ x and ∇ X are the Eulerian and Lagrangian gradient operators.

FSI model problem
The non-conservative ALE Navier-Stokes equations 19 where f and s are the fluid and solid density (Table 1), f and s are the Cauchy stress tensors of the fluid and solid, = (X) is a diffusion coefficient tensor, and and t n f are Dirichlet and Neumann boundary condition (BC) data.
with given inflowṽ I f (Section 2.1). To be able to deal with potential reflow on the outflow boundary Γ O f and prevent backflow divergence, we further introduce outflow stabilization, 45 with parameter = 0.2 as suggested in the study of Moghadam et al. 46 The Cauchy stresses in Equations 5 and 6 are written as where f is the dynamic viscosity of the fluid and the silicone material was modeled as an isotropic incompressible Neo-Hookean material with material parameter s (Table 1). Further, a change-of-variables (COV) for the pressure variables, p f and p s , was introduced, where f and s are the substituted fluid and solid pressure variables and P o is the unknown mean outlet pressure. The COV is introduced to account for the contribution of a gravitational field (with g = [0, − 9.80665,0] T m/s 2 ) to the momentum balance in Equation 5 (and likewise for Equation 6 to ensure compatibility of stresses). Without COV and under zero inflow conditions, a linear pressure gradient along the y-axis would yield flow in the outlet region due to violation of a zero traction assumption on Γ O f . On the other hand, the same conditions do not cause flow if a COV is employed.
The real pressure is recovered by reversing the COV with In general, the shape of the fluid domain changes because of the deforming solid. To extend the deformation of the fluid-solid boundary (with fixed to the interior of the fluid domain, a mesh velocity w f is introduced that satisfies Equation 7. Here, the boundary motion is extended to the interior using the typical harmonic extension augmented with a coefficient tensor, , to preferentially weight motion in the boundary normal direction. Specifically, the tensor was defined as . As the variable yields a decaying diffusion field from the coupling interface, this introduces anisotropy in preferentially distributing motion faster in the direction of ∇ X . The value set for was selected based on simple test problems in three dimensions, illustrating preservation of the mesh quality over time.

Finite element formulation using Lagrange multipliers
To couple the fluid and solid models, dynamic and kinematic interface constraints are enforced with equal but opposite tractions, t f := f ·n f and t s := s ·n s , and the no-slip condition, as detailed in Equation 8. We note that the interface constraints are not influenced by the introduced COV because yielding an equivalent formulation. On the coupling domain Ω , a Lagrange multiplier variable is defined as = t f =− t s , maintaining the interface conditions. Within the continuous weak form for both fluid and solid mechanical subsystems, the resulting boundary terms are substituted with , constraining the multiplier through the no-slip interface condition. Finally, the respective fluid and solid models (given in Equations 5 and 6) and the interface constraints (Equation 8) are coupled monolithically. Here, we consider an inf-sup stable scheme and the cG(1)cG(1) scheme 21 as a variant of a General Galerkin method. We note, that Equation 7 is partitioned from the main FSI problem and solved in turn.

Spatiotemporal discretization
A first-order backward Euler scheme is employed for the time The fluid domain Ω f was discretized using tetrahedral elements, whereas hexahedral elements were selected to discretize the solid domain Ω s . The coupling domain Ω was discretized using triangular elements that conform with the fluid domain surface elements on Γ C f . In the following, we use to denote meshes consisting of N e i many, non-overlapping elements Ω e i,h and Ω n f ,h =  and to denote current fluid and solid meshes.

Weak formulation I: inf-sup stable scheme
Finite element discretizations were constructed using inf-sup stable P 2 − P 1 Taylor-Hood elements for fluid velocity and pressure, P 2 for fluid domain velocity and inf-sup stable Q 2 − Q 1 Taylor-Hood elements for solid displacement and pressure. Further, P 2 elements were employed for the Lagrange multiplier variable that was nested into the trace of the richest space on Ω , which was on the fluid side. The discrete solution at each step n in time can then be written as follows: stress tensor, u n+ = u n + Δ n t v n+1 , and the notation Y n + = Y n + 1 + (1 − )Y n .
Following the study of Nordsletten et al, 16 we introduce, defining all k th -order piecewise continuous polynomials on respective discretized domains (i = f,s), where P k are the polynomials of degree k, k ⩾ 1. Finally, from the spaces, we select only those functions satisfying the Dirichlet BC or homogeneous BC, and define the Lagrange multiplier space as, is the trace operator on Γ C f ,h .

Weak formulation II: cG(1)cG(1) scheme
Here, the cG(1)cG(1) scheme as given for the Navier-Stokes equations in the study of Hoffman et al 21 is considered to be able to employ equal-order interpolation and to govern potential turbulent effects. Finite element discretizations were constructed using P 1 − P 1 Taylor-Hood elements for fluid velocity and pressure, P 1 elements for fluid domain velocity and inf-sup stable Q 2 −Q 1 Taylor-Hood elements for solid displacement and pressure. Further, P 1 elements were employed for the Lagrange multiplier variable, that was again nested into the trace of the richest space on Ω , which was on the fluid side. The discrete solution at each step n in time can then be written as follows: with the stabilization term,

SD n+1
( and stabilization parameters, where h is the mesh size and v max is the expected peak velocity on Ω  Further, the Shamanskii-Newton-Raphson (SNR) method 49 was employed to reduce computational cost by re-using the Jacobian matrix (and its inverse) 5,18 as long as a sufficient decrease in the residual norm was observed. A pseudocode description of the employed algorithm is given in Algorithm 1, where c = 100 is used for scaling the initial residual, = 3/4 requires a sufficient residual decrease, and the general notation, is used to define the current approximate solution and test functions (with basis functions and ). Further, the residuals and Jacobians are denoted by Because of the iterative nature of the algorithm (Algorithm 1) and the partitioned approach, w d f can be updated from the current solution of v f on Γ C f . *http://cheart.co.uk

Error norms
Numerical results are compared with experimental data at equispaced data point positions at cross sections along the z-direction with z ∈ {3.5, 13.5, … , 93.5}mm and 5 mm spacing along the x-and y-direction, as specified in the studies of Gaddum et al 36  We note that a mask is employed such that only points are considered that exist inside both fluid domains in the current configuration. For example, consider discretized domains Ω i,h (t p ) and Ω i,h (t q ) at times t p and t q . Then, the mask is given as, For comparison of results at sampled data points x i , consider an approximation w(x i , t p ) ∈ R d and a reference v(x i , t q ) ∈ R d at times t p and t q with d ∈ {1, 3}. Then, a scalar field is defined at each point x i by computing the Euclidean distance, ||·|| 2 , between approximation w and reference v. To detect variation in the fluid and solid variables, we compute the maximum distance, (14) and the mean distance, between an approximation and a given reference.

3D FSI BENCHMARK RESULTS AND COMPARISON WITH EXPERIMENTAL DATA
Numerical results for both phases of the 3D FSI benchmark were obtained using the inf-sup stable (Section 2.4.2) and cG(1)cG(1) (Section 2.4.3) schemes. In this section, we present and compare the numerical results with the experimental data, while studying spatial refinement and studying how quickly steady-state and periodic steady-state are reached.

Mesh construction
For Phase I, meshes were constructed based on undeformed geometries corresponding to a stress-free state of the coupled FSI system and reference domains Ω 0 f , Ω 0 s , and Ω 0 were selected accordingly (Figure 3).
For Phase II, the Phase I meshes were deformed based on the coupled FSI system being exposed to gravity loading of 70% strength. Then, fluid and interface domains were remeshed, to counteract deteriorating mesh quality due to large mesh deformations, and the fluid and reference domains redefined to be in this deformed state ( Figure 3). Finally, the coupled FSI system was subjected to 100% gravity loading to obtain a deformed state that matches the state during the calibration step of Phase II of the FSI experiments. 37 In preliminary numerical experiments, the static deflection of the solid in its hydrostatic equilibrium (with material parameters for Phase II; Table 1) was studied regarding mesh resolution. Buoyancy forces in a surrounding resting fluid were mimicked by applying the net gravity load, specifically ( s − f )/ s g. Equation 6 was used for the solid model, but the inertial term was neglected (because predictions of a static and transient model are expected to tend to the same asymptotic limit subject to static loading conditions). As expected, because of the high-aspect ratio of the solid side lengths, the observed maximum deflection of the tip of the solid is predominantly influenced by the mesh size in z-direction (Figure 4). A similar result is anticipated for dynamic loading cases, such that a fine mesh resolution of 11 × 2 × 65 hexahedral elements seems a good choice for both Phases I and II.
For the model given in Section 2.4.2, the fluid domain boundary f was discretized using a fine triangular surface mesh. On the coupling boundary Γ C f , the triangular mesh was constructed, such that a quadrilateral surface element on Γ C s would have two triangular surface elements as  Table 2 FIGURE 4 The graph depicts the position of the center of the solid tip, u y + Y, under static gravity loading depending on mesh resolution, where the number of elements in one coordinate direction is kept fixed and varied in the other coordinate directions counterparts on Γ C f with matching nodes. Aim of selecting a fixed but fine mesh resolution on f was to avoid dominant boundary effects while studying spatial refinement in the interior of the domain, where three different mesh refinement levels were considered (referred to as rl-1, rl-2, rl-3), see Figure 3 and Table 2. On the other hand, only one refinement level rl-1 was considered for the model given in Section 2.4.3 ( Figure 3G), as we expect similar behavior with respect to spatial refinement. Here, the fluid coupling boundary Γ C f was discretized using triangle  elements, such that a quadrilateral surface element on Γ C s would have two triangular surface elements as counterparts on Γ C f . We note that the outlet of the computational domain was extended by 250 mm to regularize flow downstream. To study the necessity of such an elongation, a truncated domain ( Figure 3A) in combination with the model given in Section 2.4.2 was also considered with two refinement levels (referred to as rl-1-trunc and rl-2-trunc) for Phases I and II, respectively.
The triangular surface mesh on the fluid domain's coupling boundary was extracted from the mesh discretizing Ω 0 f to define the mesh on Ω .

Phase I experiment
Initially, the FSI system is at rest and subject to no body forces (no flow, no displacement). Then, parabolic inflow is prescribed at both inlets as given in Equation 1 (Figure 2A). Simultaneously, gravitational forces are increased according to (24t 3 − 8t 2 ) · g over 0.5 s and kept constant at 1 · g for t ⩽0.5 s. Because a steady-state is expected, 37 an increasing sequence of time step sizes was employed; eg, we simulated 2 s with Δ t = 1 ms, 60 s with Δt ′ = 10 ms, and 138 s with Δt ′′ = 100 ms (such that the total simulated time was T = 200 s). Because of increased inflow and gravitational forces, the silicone filament moves into the way of the incoming flow jet.
The impact causes strong disturbance in the flow field, which is then advected toward the outlet causing reflow regions at the outflow boundary. Here, backflow divergence is avoided effectively by employing outflow stabilization (Equation 9).
After a short transition phase, the silicone filament reaches a steady position and flow stabilizes. Figure 5 illustrates that the final position of the silicone filament along its centerline is predicted consistently by the inf-sup stable scheme with a slightly overestimated deflection. On the other hand, the cG(1)cG(1) scheme slightly underestimates the deflection of the silicone filament near the tip. However, agreement with the experimental data is very good considering finite voxel sizes in MRI. 37 To investigate how quickly the silicone filament reaches its steady-state position (where the position predicted by the inf-sup stable scheme with fluid domain refinement level 3 is assumed as the steady-state position), we employ Equation 14 to quantify the maximum Euclidean distance to the steady-state position. Figure 6 shows that the solid reaches its final position quickly irrespective of the fluid domain mesh resolution or employed scheme, such that the fluid domain deformation is negligible after the initial transition phase.
Once flow has stabilized and fluctuations become smaller, flow in the v x and v y components is less pronounced (Figure 7) than during the transition phase. Flow mainly occurs in the v z component, and the presence of the silicone filament does not seem to cause significant disturbance or deflec-tion of the flow jets entering from the inlets. In fact, the tip of the silicone filament is positioned just below the upper flow jet. Between the two flow jets, a large recirculation zone is observed, where approximately one third of the outflow boundary is covered by fluid flowing back into the computational domain, see Figure 7C,D. Investigating differences in the flow predictions obtained for various fluid domain refinement levels using both schemes, it becomes clear that an early truncation of the domain can change local flow in the considered region significantly, see Figure 6. Further, flow predictions computed for the elongated fluid domain are consistent for all meshes indicating that rl-1 provides a sufficient refinement. Similarly to the solid deflection, a steady-state flow field is obtained quickly (Figure 6), such that the simulated time can be reduced drastically without corrupting the approximation accuracy of the predicted steady-state. Finally, comparing flow predictions and experimental results at plane z ≈ 30 mm shows good qualitative agreement for all velocity components, see Figure 8A-C.

Phase II experiment
Initially, the FSI system is at rest with the initial configuration corresponding to the hydrostatic equilibrium (that is, solid deformation due to gravitational forces). Parabolic inflow is prescribed at both inlets with peak values fit to experimental data ( Figure 2B). A total of K = 10 cycles of the periodic inflow pulse seen in Figure 2B were simulated with time step size 1 ms. In the following, we refer to each cycle individually by using the notation t n k ∈ [0, 6] s (with k = 1,2, … ,10 and n = 0,1, … ,6000) to simplify comparison of corresponding time steps of two different cycles, eg, cycles k and K with t n k and t n K . Starting from its initial position at the beginning of each cycle (referred to as resting position), the silicone filament is pushed downwards because of the impacting flow jet entering from the upper inlet (the time-dependent position of the solid's centerline is given in Figure 9). Once flow decelerates because of decreased inflow, gravitational forces become more dominant, and the silicone filament moves back upwards performing a swing before reaching its resting position. The silicone filament follows the same repeatable deflection pattern as observed in the experiment; however, its resting position coincides with the deflection under zero inflow conditions, which was not found in the experiment. The relative displacement of the silicone filament ( Figure 9A) is predicted well. Even though the maximum relative displacement at t ≈ 1.15 s is overestimated by approximately 1.8 mm (on average), the instant when the maximum relative displacement is reached is predicted precisely and the swing at t ≈ 1.59 s governed (however, a short delay is observed). While the deflection pattern is consistently predicted by the inf-sup stable scheme on all refinement levels (irrespective of fluid domain length), the cG(1)cG(1) scheme yields a slightly different relative displacement pattern with a lagging solid motion, thus the silicone filament reaching its resting position later ( Figure 9A). Further, the swing is less pronounced. We note that the numerical results in Figure 9B are Lagrangian positions. On the other hand, the position of the solid under flow conditions was obtained from intersecting the solid's centerline with MRI image planes at Z s ≈ 3 mm + s·10 mm, s = 0, … , 5. Thus, the position of the solid is not available for its entire length. Similarly, the position of the solid under zero inflow conditions was extracted from an MRI image plane at x ≈ 0 mm.
Flow patterns show strong similarities with those observed in the experiment. For example, a double Ω-shape is observed at flow planes that drifts in negative x-direction along the v z -direction ( Figure 10A) because of the non-zero v x -component at the inflow boundary. Further, the flow dynamics and local flow phenomena are governed well and compare well with the experimental data † (compare Figure 8D-L and the study of Hessenthaler et al 37 ). We note that minor differences in the local flow velocities are observed if the flow domain is truncated as exemplified for the v y -component in Figure 11B. For example, finer resolution yields a more symmetric prediction of the flow field, and stronger flow in the upward direction is observed at the truncation and near the tip of the silicone filament. However, it does not affect the position of the silicone filament ( Figure 9A).
Forces exerted onto the solid were computed by integrating the Lagrange multiplier variable (which represents the † http://cheart.co.uk/other-projects/fsi-benchmark/  Velocity v z at t ≈ 2.02 s. B, Phase II: Velocity v y at t = 3.96 s (volume rendered with opacity ranging from 0 to 1) obtained with the inf-sup stable scheme and the elongated (left; rl-3) and the truncated computational domain (right; rl-1). Here, both domains are truncated at z = 160 mm to simplify comparison traction on the coupling boundary Ω ), where n , t , and s are linearly independent normal and tangent unit vectors on Ω . As Figure 9C illustrates, the f y component is negative for downward motion of the silicone filament and positive when it starts moving upwards after the largest relative displacement. For the f x and f y components, the truncation of the domain does not induce significant differences; however, the maximum force in f z -direction during cycle 10 is observed earlier. Forces predicted by the cG(1)cG(1) scheme are significantly different; however, match with observations of relative solid motion (for example, the smaller f y force after the largest relative displacement is in good agreement with the silicone filament moving toward its resting position more slowly).
To investigate how quickly the FSI system reaches a periodic steady-state, we employ Equations 14 and 15 to detect fluctuations in the solid motion and flow (assuming cycle 10 with the inf-sup stable scheme and rl-3 as the periodic steady-state). Here, we consider the inf-sup stable scheme and refinement levels 1-3. As Figure 12A illustrates, the solid deformation reaches a periodic steady-state after a small number of cycles on all refinement levels, and the same holds for the flow field, see Figure 12B,C. Although Figure 12B indicates that a different periodic steady-state flow field was found on refinement levels 1 and 2, the quite large error is related to insufficient fluid domain refinement in regions of large velocity gradients (for example, note the more diffused appearance of the double Ω-shape in Figure 10B) such that the error norm given in Equation 15 is a much more suitable measure of error in the predicted flow field ( Figure 12C).
An increased average runtime for Phase II was observed for the cG(1)cG(1) scheme (Table 2), which arises from more Newton iterations per solve step. This stems from a general delay in rebuilding the Jacobian matrix (Algorithm 1). However, we note that the SNR solver parameters were tailored for the inf-sup stable approach and subsequently used for both methods. Runtime would likely improve significantly for the cG(1)cG(1) approach with adaption of SNR solver.

DISCUSSION AND CONCLUSION
An inf-sup stable FSI scheme has been validated and compared with, and a stabilized cG(1)cG(1) scheme (presented in Sections 2.4.2 and 2.4.3) has been validated using benchmark data obtained from a recently developed 3D FSI experiment. 37 For the Phase I experiment, numerical results for predicted solid displacement and flow velocity were found to agree exceptionally well with available experimental data with solid deflection slightly overestimated/underestimated by the inf-sup stable / cG(1)cG(1) scheme. On the other hand, numerical results for the Phase II experiment yield a good prediction of flow patterns (eg, Figure 10A), but it was found that a fine mesh resolution is required to sufficiently resolve flow regions with large velocity gradients ( Figure 10B).
Overall, the inf-sup stable scheme performed better in the considered cases regarding prediction of relative solid displacement or shape of flow patterns. Although, employed mesh resolution for the cG(1)cG(1) can be considered coarse.
It is recommended to simulate Phase I for about 30 seconds and Phase II for at least 5 cycles to reach an approximate (periodic) steady-state with sufficient accuracy (with respect to initial / boundary conditions and error norms employed in this paper, eg, no significant fluctuations in solid motion and flow are observed; Figures 6 and 12). Results show that an elongation of the computational domain by 250 mm yields significantly better accuracy than an elongation by 50 mm. Further, mesh resolution has to be tuned depending on the phase of the benchmark because dynamic effects need to be governed in Phase II, thus requiring a finer fluid domain discretization as compared to Phase I. However, computational cost depends on the selected spatial mesh refinement of the fluid domain, as can be appreciated in Table 2. We note that the increased average runtime for Phase II in the case of the cG(1)cG(1) scheme stems from an observed increase in the number of Newton iterations. For the cG(1)cG(1) scheme, parameters to trigger rebuilding the Jacobian (see Algorithm 1) were chosen based on the inf-sup stable scheme, however, delay the rebuild step in the transient test case and have thus to be tuned for better performance. From numerical experiments presented in this paper, we conclude that solid motion is not significantly altered by fluid mesh resolution; however, solid motion is influenced by the selected scheme.
Future investigations on the considered 3D FSI benchmark will include, for example, studying temporal convergence, sensitivity to the selected solid material law, spatial mesh refinement for the cG(1)cG(1) scheme, a second-order time-integration scheme, and sensitivity to outflow stabilization (preliminary results with = 0.02 indicate no dependency on the outflow stabilization; however, final observations are likely to depend on the selected fluid domain length). Further, computational cost can potentially be decreased by optimizing the length of the computational domain, using a coarser solid mesh resolution and making use of the symmetry in the Phase I setup.