Methodology for Comparing Coupling Algorithms for Fluid-Structure Interaction Problems

The multi-physics simulation of coupled fluid-structure interaction problems, with disjoint fluid and solid domains, requires one to choose a method for enforcing the fluid-structure coupling at the interface between solid and fluid. While it is common knowledge that the choice of coupling technique can be very problem dependent, there exists no satisfactory coupling comparison methodology that allows for conclusions to be drawn with respect to the comparison of computational cost and solution accuracy for a given scenario. In this work, we develop a computational framework where all aspects of the computation can be held constant, save for the method in which the coupled nature of the fluid-structure equations is enforced. To enable a fair comparison of coupling methods, all simulations presented in this work are implemented within a single numerical framework within the deal.ii [1] finite element library. We have chosen the two-dimensional benchmark test problem of Turek and Hron [2] as an example to examine the relative accuracy of the coupling methods studied; however, the comparison technique is equally applicable to more complex problems. We show that for the specific case considered herein the monolithic approach outperforms partitioned and quasi-direct methods; however, this result is problem dependent and we discuss computational and modeling aspects which may affect other comparison studies.


Introduction
In this work, we focus on constructing a numerical framework in which algorithms for solving the class of engineering problems commonly referred to as Fluid-Structure Interaction (FSI) problems can be benchmarked and compared in a consistent and quantitative manner.The primary challenge in finding a numerical solution to a FSI problem is not in solving the field equations over either the solid of fluid components of the domain, but rather handling the compatibility conditions that must be enforced on the interface between these components. 1  Recently, two primary strategies for enforcing this coupl-ing in the Arbitrary Lagrangian-Eulerian (ALE) framework have emerged: solving the solid and fluid components of the domain separately, while iteratively passing information between them, and a simultaneous solution of all field equations at once.The former method is referred to in this work as "partitioned", and the latter as "monolithic".
The choice of coupling strategy is an important decision to make from the start of any endeavor to implement a numerical ALE framework for the simulation of fluidstructure interaction, as it places inherent constraints on the entire effort moving forward.On one hand, an author must be sure to choose an approach that allows for all of the relevant physics to be captured, but must also consider the performance of the resulting computer code, the ease of implementation, and its ability to address problems of interest.If one turns to the work of other authors for guidance, conflicting statements can be found regarding the performance of various coupling approaches.In addition to this, there is seemingly conventional wisdom in the community today that partitioned approaches are faster than their monolithic counterparts due to the optimized nature of the individual component solvers and smaller systems inverted by the linear solvers [4][5][6][7]; however, these notions have been challenged in [8].There is agreement within the FSI community that, in general, the choice of coupling type is problem dependent, although no specific guidance exists on exactly when a coupling strategy should be preferred.
Michler et al. [9] state that in a simple one dimensional FSI problem "the computational cost of the monolithic procedure is three to four times the one of the partitioned procedure;" however, he does not hold each approach to the same convergence tolerance.In fact, he requires more strict convergence of the monolithic solver, which contributes to its longer run time.In contrast, he later hypothesizes that for more complex problems the monolithic approach would be superior in a comparison of computational cost for a desired accuracy against his partitioned method.
Degroote et al. [6] utilize ADINA [10], a commercial finite element software package, to compare partitioned and monolithic FSI.They claim that ADINA converges both methods to the same tolerances, for the norms of displacement and force on the interface; however, it is unclear exactly what convergence criterion is used as no equations or methods are presented.They find that for certain cases their partitioned method requires half the computational time of their monolithic method, while for other cases it requires as much as four times the computational time of the monolithic method, and in some cases, the partitioned method does not converge at all.
On the other hand, Heil et al. [8] use the open source finite element library OOMPH-LIB [11] to compare partitioned and monolithic methods.In their tests, they found that, in terms of computation time required, monolithic methods outperform partitioned methods when FSI effects were significant and that this behavior was amplified with proper pre-conditioners.Again, it is unclear exactly what parameters they use to signify convergence.They claim that the convergence criterion is either the maximum global residual, or the maximum change in position of nodes between iterations; however, in all of their cases they do not require the partitioned method to converge to the same tolerances as the monolithic method.
As can be seen by these examples, although there has been past work in comparing coupling techniques, there is no clear agreement on which method is optimal for a given problem.To aid in providing clarity to this situation, we have developed an ALE-based framework that can hold constant many of the confounding factors and can focus solely on the fluid-structure coupling methodology.This includes performing all calculations on the same mesh, using the same discretization scheme for all approaches, and using the same benchmark problem for all approaches.Most importantly, we use the same convergence criterion and linear solvers for all coupling methods, ensuring that we are eliminating as many confounding factors as possible.We believe that this type of framework is necessary in order to provide a fair comparison between FSI coupling types, and has clearly been lacking in prior work.
We note that we are interested in algorithms appropriate for fully-coupled FSI problems, i.e., where the structural and fluid responses are tightly intertwined and require a truly simultaneous solution.Problems requiring only one-way or weak coupling are outside the scope of this study.We begin this presentation by deriving the governing equations necessary for FSI in Sections 2 -3.The discretization procedure is presented in Section 4, the three coupling strategies we are comparing are discussed in Section 5, and formal verification of the implemented algorithm is demonstrated in Section 6.Finally, results of the study and our conclusions are discussed in Sections 7 -8.

Kinematics
Fluid-structure interactions are modeled as initial boundary value problems (IBVPs).This section presents an overview of the mathematical background necessary to describe these problems.This material is presented from a continuum mechanics perspective, see Gurtin [12], Bowen [13], Chadwick [14], or Spencer [15] for a more detailed discussion.
Kinematics is the geometrical description of motion.It provides descriptions for the movement of bodies (a collection of material points) and changes of reference frames.We denote a body 3 β ∈  in the reference configuration 3   κ as κ β , where the reference confi- guration is chosen to be the configuration of the body initially before deformation.The same body β can be defined in subsequent configurations through a deformation function χ .The position of a material point in κ is ( ) ,t X X , while the position of a material point in the deformed configuration, known as the spatial point, is x .The deformation function ( ) ,t χ X maps a material point X to the spatial point x at the instant of time t , as shown in Figure 1.
Considering a smooth sequence of configurations ordered in time, a motion is defined as where u is the displacement vector, given by ( ) The domain Ω of an FSI problem, as shown in Figure 2, is defined as the union of the disjoint solid and fluid subdomains: The displacement vector is piecewise defined over the solid and fluid subregions as where s u is the solid displacement and m u is the (arbitrary) displacement of the fluid subdomain.Continuity of displacement on the fluid-solid interface requires The deformation gradient is the spatial gradient of the motion, defined as ( ) where I is the identity tensor.In the next section, we will use F for push-forward and pull-back operations in the Arbitrary-Lagrangian-Eulerian (ALE) description of the fluid region.Operators with capitalized first letters (e.g.Grad) denote differential operators with respect to X , while those with lower-case first letters (e.g.grad) denote differential operators with respect to x .The determinant of F is ( )

Governing Equations
Governing equations are generated by complementing balance laws with constitutive relations for specific materials and kinematic relations.We consider solid materials that are modeled as St. Venant-Kirchhoff hyperelastic isotropic solids.Fluids considered are all incompressible and linearly viscous (Newtonian).We choose to use the ALE description for the fluid to address the disparity on the fluid-solid interface between the reference frames naturally used for the fluid (Eulerian) and solid (Lagrangian).Other approaches, such as space-time techniques [16,17], or immersed boundary/body methods [18,19], could also have been implemented, but we chose to use the ALE description for the purposes of our study.The ALE fluid deformation is described with the equations of linearized elasticity.

Elasticity
The governing differential equations for an elastic solid are cast in the Lagrangian (referential) frame.The balance of linear momentum is where s v is the solid's velocity field, ( ) is the body force per unit mass, and is the first Piola-Kirchhoff stress.The first Piola-Kirchhoff stress is related to the second Piola-Kirchhoff stress ( ) (8) For an isotropic material, The separate displacement and velocity fields necessitate the use of the kinematic compatibility condition Strong enforcement of (11) yields the familiar governing differential equation for elastic solids solely in terms of s u :

Div
. t In this work, ( 7) and ( 11) are referred to as the twofield , u v formulation of elastodynamics, while (12) is referred to as the one-field ( ) s u formulation of elastodynamics.We choose to use the two-field formulation in our FSI formulation, as it is very convenient to have the velocity as a primitive variable on the fluidsolid interface (see Section 3.4).

Incompressible Navier-Stokes
The Eulerian frame is the natural setting for the mathematical description of fluid behavior.In the context of FSI descriptions, however, we need to unify the mathematical description of the fluid and solid domains (on the interface, at least).The ALE description of fluid behavior utilizes a smooth map between the referential and spatial domains.We begin by providing the incompressible Navier-Stokes equations in the Eulerian frame, and then transform those equations into the ALE description.

Eulerian Description
The balance of linear momentum in the Eulerian frame is ( ) where ( ) is the current deformed fluid domain, ρ is the mass density, and T is the Cauchy stress tensor.We denote the velocity field in the fluid domain as f = v v .Mass balance for an incompressible flow reduces to the divergence free velocity condition ( ) The Cauchy stress for a linearly viscous, incompressible fluid, is given as where p is the fluid pressure and µ is the fluid dynamic viscosity.Equations ( 13), (14), and (15) are collectively known as the Navier-Stokes equations for an incompressible flow.

Arbitrary Lagrangian-Eulerian Description
The ALE formulation is used to provide a unified reference frame for fluid-structure interaction problems.This requires transforming the fluid governing differential Equations ( 13), (14), and (15), into a form defined on the reference configuration ( ) The resulting equations are termed the 'incompressible Navier-Stokes equations in ALE form.' The following relations are necessary to construct the transformation: (16) and where m v is the fluid domain velocity, defined as . The κ subscript is used to denote quan- tities transformed from spatial coordinates into the ALE reference frame.Using these relations, the ALE form of the incompressible Navier-Stokes equations 2 are ( )

Mesh Motion
The computational mesh is a non-physical object which can be constructed in any manner we see as beneficial; this is the "arbitrary" component of the ALE formulation.We choose to treat the mesh as a linearly elastic material governed by an elastostatic equation (derived from the balance of linear momentum (7)): For linearly elastic (infinitesimal deformation) materials, which require that m Grad 1 u  , the first Piola- 2 Hron and Turek [20], Richter and Wick [21], and perhaps others, present a different relation for the balance of mass in an ALE frame than the one we present in (19).It can be shown that ( 19) is a simplified but identical, version of these through application of equation (2.2.30) from Bowen [13].
Kirchhoff stress, when linearized with respect to m Gradu , is given as where m λ and m µ are arbitrarily defined mesh material properties and e is the infinitesimal strain tensor.The infinitesimal strain tensor ( ) ( ) results from linearizing E with respect to m Gradu .

Boundary and Initial Conditions
In addition to the governing differential equations, there are also boundary conditions, fluid-structure interfacial conditions, and initial conditions necessary to fully define a FSI problem.A boundary condition that prescribes a value for a field variable is a Dirichlet boundary condition, while a boundary condition involving spatial derivatives of a field variable is a Neumann boundary condition.Initial conditions prescribe a value for a field variable over the entire domain (or a subdomain) at the initial time.Specific boundary and initial conditions are problem dependent, but the general form can be defined.
The primitive variables can be can be prescribed on the Dirichlet boundaries as ( ) ( ) where time differentiation is denoted with a superposed dot as t . The tractions can be prescribed on the Neumann boundaries as For the mesh, a constraint is placed on the entire boundary that there be no normal mesh displacement, to maintain a conforming geometry: The initial conditions for the primitive variables can be prescribed over their respective domains at FSI coupling imposes three additional conditions on the fluid-solid interface.The solid displacement governs the mesh displacement on the interface: ( ) the velocity fields are continuous across the interface: ( ) and the tractions are continuous across the interface: where the orientation of the interface requires s κ = − n n .

Discretization via Finite Elements
The finite element method (FEM) is used to solve an IBVP by transforming its governing partial differential equations (PDEs) into a system of algebraic equations that approximately solve the original problem.The PDEs are weighted by an arbitrary function, and integrated over their respective domains.Integration by parts yields the weak form of the equations.In this work, the notation is used that ( ) . The ALE form of the fluid equations requires us to compute the spatial integrals over a reference configuration; as such, the following transformation from the Eulerian to the ALE frame is necessary: In all future equations ( ) for brevity.A sufficiently regular computational mesh is used to discretize the domain Ω , and we ensure that the fluidstructure interface FS Γ is coincident with element edges.The approximate solutions are denoted by the superscript h .We define the following standard finite-dimensional trial solution spaces on our mesh: Our weighting function spaces are With these definitions, it is clear that our finite element method is of the standard (Bubnov-) Galerkin weak formulation, in which weighting function spaces and trial solution spaces are taken to be the same except on the Dirichlet boundaries [22].
The temporal discretization is handled with the implicit, second-order accurate Crank-Nicolson method.The choice of time integration, while important for accuracy and efficiency of a method, is not the focus of our current work.We are interested in investigating different strategies for solving the fully-coupled FSI system.As such, in the remaining mathematical formulation of the FSI problem, we shall omit the time discretization aspect of the problem for the sake of brevity.The weak formulation of the problem is understood to be integrated over a time step with a suitable approximation of the time derivatives.

Weak Formulation
An FSI problem consists of three distinct sub-problems which are coupled through boundary conditions on the fluid-structure interface FS Γ .We introduce the Galerkin formulation for each sub-problem separately, followed by the statement of the global FSI problem.
Problem 1 (Solid sub-problem) where k Ω is a constant over the solid domain with dimensions of mass per length cubed per time cubed to give our equations dimensional consistency.For the solids considered in this work, the first Piola-Kirchhoff stress R T uses the St. Venant-Kirchhoff strain tensor given in (10).
Problem 2 (Fluid sub-problem) where F and J are as given in ( 5) and ( 6), and, for the linearly viscous, incompressible fluids considered in this work, the Cauchy stress κ T is as given in (20).λ and m µ can can result in poor quality meshes if the solid displacement is large enough3 .To prevent this, an artificial stiffness is imposed on the mesh, using ( ) , in a manner such that it is stiffer near the interface and less stiff away from the interface.This propagates the distortion caused by the solid displacement throughout the entire mesh, instead of it being localized at the interface.
The global FSI problem is a direct sum of Problems 1, 2, and 3.
The individual sub-problems are coupled solely via the boundary conditions on the FSI interface FS Γ .The globally defined system in Problem 4 is a fully coupled nonlinear system of equations.We have employed the Newton-Raphson method to solve these equations, including all non-linear sub-problems that arise.We present some of the details in Appendix 9 since the calculations can be daunting, especially for the monolithic method.

FSI Algorithms
The numerical solution of Problem 4 is the focus of our study.We examine three different types of solution procedures: partitioned, monolithic, and quasi-direct.These are often referred to as FSI coupling methods, although they can also be thought of as a family of algebraic solution methods to solve the fully coupled Problem 4.
One of the key features of our study is that all three solution algorithms are implemented within the same software package (deal.ii).This feature allows us to use the same finite element assembly routines for all FSI coupling types, ensuring that computational cost differences between the methods studied is not due to different implementations.Furthermore, the convergence criterion chosen are applied in exactly the same manner for each FSI algorithm, as described in Section 7.2.
The global solution procedure of the three coupling types is the same.After a setup procedure and application of initial conditions, a time loop advances the solution with a prescribed time step size until the final simulation time is reached.For each time step, the three FSI coupling algorithms solve Problem 4 in different manners.We refer to the non-linear solution procedure within a single time step as FSI iterations.

Partitioned Coupling
The partitioned coupling method couples Problems 1, 2, and 3 in a sequential fashion to solve the global FSI problem.This coupling method is also referred to as Dirichlet-Neumann coupling, due to the way constraints are imposed on the FSI interface.There are other types of partitioned coupling, but we choose Dirchilet-Neumann coupling as it is by far the most common type used in practice: used by Turek et al. [2], Degroote et al. [6], Küttler et al. [23], and Campbell et al. [24], amongst many others.The FSI iterations for the partitioned method involve solving the fluid, solving the solid, and solving the mesh; this process is iterated until convergence is attained.The partitioned algorithm is detailed in Figure 3(a).
We solve (48) with prescribed velocities on , , respectively.The strong enforcement of the velocity coupling is achieved in our implementation through the use of deal.ii'sconstraint tools [1].Newton-Raphson iterations are used to solve the non-linear fluid subproblem.
Once a fluid solution is obtained, Problem 1 is solved over the solid domain using the fluid traction as the driving force on the FSI interface.That is,  ( ) (52) which comes from ( 35) and (20).The function space definitions on S Ω do not need modification for the partitioned method.Due to the geometric non-linearities in Problem 1, a Newton-Raphson method is used to solve for the solid displacements and velocities.
If the change in position of FS Γ is non-trivial (see Section 7.2 for discussion of this criterion), we deem it necessary to update the mesh position by solving Problem 3. Before solving for the mesh displacement, it is well known [23] that the interfacial displacements should be under-relaxed to ensure convergence.The solid displacements and velocities are relaxed as ( ) where the subscript i denotes the iteration count and the superposed tilde denotes the non-relaxed solution from solving Problem 1.We find it beneficial (or, at a minimum, non-detrimental) to relax the entire solid solution, rather than just the degrees of freedom whose support is on FS Γ .The scalar relaxation parameter 1 i ω + can be chosen to be constant, but it has been shown [23] that allowing it to change between iterations can accelerate convergence.Aitken's 2 δ convergence acceleration method seems to be the most popular algorithm for choosing ω , and we have used it in this work.The position of the fluid-solid interface is viewed as a set of fixed points for the partitioned FSI iterations, and the secant method is used to extrapolate the expected position of the interface.The result is the recursive formula for computing the relaxation factor ( ) .The Euclidean inner product and norm are used for the discrete vector.For the first FSI iteration of each timestep a relaxation factor of 0.75 is used.
Once the solid solution has been relaxed, we update the mesh displacements by solving Problem 3. The function space definitions (41) and ( 46) are augmented as respectively, to account for enforcing a Dirichlet condition on the fluid-structure interface.Unlike the fluid and solid sub-problems, the mesh motion is governed by a linear system of equations and only requires a single algebraic solve.

Quasi-Direct Coupling
The second algorithm of our study for quasi-direct coupling, as first performed by Tezduyar [25], is detailed in Figure 3(b).The algorithm solves Problems 1 and 2 monolithically, and an iterative procedure is used to couple with Problem 3 to fully solve the complete FSI Problem 4. For each time step, we solve the coupled fluid-solid system, update the mesh displacements, and iterate on these two solves until convergence is reached.The quasi-direct algorithm does not use the Dirichlet-Neumann coupling scheme that the partitioned method uses.Rather, we construct a conforming velocity field over the entire reference domain Ω as The consequence of this choice is that (34) is strongly satisfied, which implies that after enforcing (35) and recognizing that the velocity weighting functions are chosen from the same function space.Another way to recognize that the tractions on = ∅ by construction of a conforming (continuous) velocity space.Therefore, the traction continuity is implicitly satisfied and does not enter the quasi-direct (or monolithic) equations [26].
The solution from the fluid-solid does not need to be relaxed as it was for the partitioned algorithm.The relaxation procedure is only used for iterative Dirichlet-Neumann coupling schemes for stability.Our only requirements on the mesh motion are the same as for the partitioned case.Namely, m h u is computed using the function spaces given in ( 56) and (57).

Monolithic Coupling
The full monolithically coupled solution procedure is the third algorithm of our study and is detailed in Figure 3(c).Problem 4 is solved as a single non-linear system, with no FSI iterations necessary.We use the conforming velocity space given in ( 58) and (59) to strongly couple the fluid and solid domains.
The displacement fields s h u and m h u can be combined into a single conforming field, but doing so raises some computational challenges.Unlike the fluid-solid domains which are two-way coupled (they affect each other mutually), the solid-mesh displacement fields are one-way coupled.That is, solving Problem 3 for the mesh displacements should not directly influence the displacement field of the solid, while the solution of Problem 1 provides the driving boundary condition for the mesh solution.Thus, constructing a typical conforming finite element space over the entire domain has the detrimental affect of allowing the artificial mesh 'stiffness' to limit the motion of the solid structure.
We have used a simple and efficient method to strongly enforce (33) while avoiding the artificial stiffness phenomena.We solve Problem 4 using function spaces defined by (37), ( 40), ( 41), ( 42), ( 45), ( 46), ( 58), (59).The Newton-Raphson method is used to compute the solution, and we provide details of the linearization in the Appendix.Within each Newton-Raphson iteration, we weakly enforce the displacement continuity constraint within Problem 3. At the end of the iteration, the nodal mesh displacements on the interface are set to be equal to the solid displacement.In our numerical computations, we have found that weakly enforcing continuity yields mesh values that are relatively close to those computed for the solid.Also, modifying the mesh solution at the end of each Newton-Raphson iteration does not seem to affect the convergence characteristics of the method.

Code Verification
Verification of the discretized equations of motion, and the software in which it is implemented, is performed via the method of manufactured solutions (see, e.g., [27,28]).We perform the following sequence of analysis for the solid and fluid solvers separately: 1) choose an exact solution for the governing equations, 2) substitute this solution in the equations to obtain analytical expressions for the corresponding forcing functions and boundary conditions, 3) implement these into the code, and 4) compute an approximate solution.The error between the exact and approximate solutions is then computed according to some norm, and the process is repeated for a sequence of uniformly refined meshes.
Figure 4 shows the results of this process for several iterations of mesh refinement.The Figure 4(a) shows the individual components of the elastic displacement and velocity error declining as mesh size decreases, with a linear rate of 1 3 p + = under the 2  L norm when using bi-quadratic 2  Q tensor-product elements for all fields.The Q tensor product elements used for the velocity and pressure, respectively.
Overall, this simple verification exercise demonstrates that the individual components of the FSI solver are functioning correctly.We will evaluate the correctness of the coupling implementation in the next section; more details regarding verification of the solvers can be found in Sheldon [29].

Model Validation
An FSI benchmark was proposed by Turek and Hron [2] for the purpose of facilitating comparison between different FSI modeling approaches and algorithms, regardless of coupling types or discretization schemes.This  benchmark consists of a numerical experiment involving laminar incompressible channel flow around a bluff body, which subsequently sheds vortical structures that excite a flexible tail until it reaches a state-state flow induced oscillation.Data were provided by Turek and Hron for the tip displacement of the elastic structure, hereafter called the flag (see Figure 5), and the drag and lift about the combined boundary of the flag and the cylinder to which the flag is attached.We used this benchmark as a means to validate the FSI model developed in this work, via direct comparison of our model's predictions to that of Turek and Hron, and as a baseline for comparing the coupling strategies to one another.

Numerical Experiment Definition
The problem domain is depicted in Figure 5, where the channel, rigid cylinder, and elastic flag are noted.Flow is from left to right, with the cylinder slightly offset in the channel.Relevant dimensions are reproduced for the reader in Table 1.We will present our analysis of repro- ducing the "FSI2" case of Turek and Hron, as it involves the largest structural deformations, and is therefore the most intense test for our solver.The material properties for this simulation are also reproduced for the reader in Table 1, where ρ is the density, ν is the Possion's ratio, µ is the viscosity, E is the Young's modulus, and U is the mean inflow velocity.As before, the s and f subscripts indicate solid and fluid properties, respectively.Figure 6 shows an example of the flow field and solid deformation, using results from FSI2 after 16.5 simulation seconds, to help in visualizing the benchmark.
The boundaries of the computational domain are labeled in Figure 7, which also shows the mesh we used for our calculations.The boundary conditions are as follows: a parabolic fluid velocity inflow on 1 Γ : ( ) a no-slip fluid velocity on 2 Γ , 4 Γ , and 6 Γ : a constant fluid pressure on 3 Γ , which is prescribed to be zero: a fixed solid displacement and velocity on 5 Γ : 5 5 s s 0, the interface conditions given in (33) and (34) on FS Γ , and the mesh displacement condition given in (29) Γ , and 6 Γ .All fields are initially zero.For the first two seconds the fluid inlet velocity is subject to a smooth increase defined by ( ) was used in all cases to satisfy the Courant-Friedrichs-Lewy (CFL) condition [30], computed from a maximum flow velocity and minimum cell diameter of approximately Q elements for the fluid velocity/pressure, 2  Q elements for the mesh displacement, and the Crank-Nicholson method for the time discretization.For consistency in comparison of the timing results, all simulations were performed on the same computer.Linear algebraic solves were done with the sparse direct solver UMFPACK [31].

Convergence Criterion
The FSI coupling strategies that we have employed in this work all require some sort of sub-iteration at each time step in the simulation process.For the partitioned schemes, we make fixed-point iterations between the solver components to converge the solution, while for the monolithic approaches we use Newton-Raphson iterations.Consistency in defining convergence criteria for the global coupled FSI solution, that is applicable to both the partitioned and monolithic approaches, is a key point of this work.If this metric is not chosen consistently, we will not be able to say that the solutions are similarly converged to the same level of error, and therefore be unable to make comparisons regarding the wall clock time.
For the global FSI solution, the convergence criterion [23,32] is given by is some convergence tolerance, which we typically defined to be 1 E−10, for all coupling methods.For the partitioned method, we perform this solution check after every instance where the solid displacement is solved, while for the monolithic approaches, we perform this check after each Newton Iteration.In all coupling approaches, once this criterion is satisfied, we exit the sub-iteration loop and move on to the next time step.

Computations
In Turek and Hron's FSI2 simulation benchmark, the system starts at rest and the speed of the inflow grows until it levels off to some steady value.After this initial transient period, which is approximately the first twelve seconds of simulation time, the solution reaches a steady state of oscillation.All results over a given periodic interval after this time are interchangeable, and as such we have normalized the comparison time interval on all forthcoming images.
We were able to reproduce the results given by Turek and Hron, for the x-and y-displacement of point A, at the mid-point of the tip of the tail, and the lift and drag generated about the combined surface of the flag and cylinder.Our calculations are nearly indistinguishable whether using the partitioned, quasi-direct, or monolithic coupling, as seen in Figures 8 and 9. Table 2 provides numerical values for the maximum and minimum of the displacement, drag, and lift from both Turek and Hron's work and our own, along with the percent error between our solutions and the validation data.The maximum deviation from the validation data for any quantity is 4%, with the three coupling methods each being the most accurate for some quantities and the least for others.Overall, the visual and numerical results demonstrate the coupling methods are all within an acceptable deviation from the validation data, and can be confirmed to functioning correctly.Even more so, we have demonstrated that using the same convergence criterion for the global FSI problem has resulted in three solutions that are quite similar to each other.oscillations began to be significant, thus requiring more sub-iterations of the FSI algorithm.As can be seen, for this benchmark case the partitioned scheme requires over twice the amount of CPU time to solve any given timestep compared to the quasi-direct scheme, and the monolithic scheme requires even less CPU time than that.Given that we have held each solver to the same convergence criterion, and given that all other aspects of the simulation are held constant, we can state that for this problem, the monolithic scheme is certainly more efficient.

Summary and Conclusions
A methodology has been presented for comparing FSI coupling algorithms in the most fair and unbiased manner we could conceive.We developed a numerical ALE framework that allows for different FSI coupling algo-

OPEN ACCESS WJM
rithms to be evaluated for accuracy and efficiency in a quantitative manner.By holding constant many of the usually confounding factors, such as the grid, mesh motion, discretization approach, and solver tolerances, we were able to compare three different coupling strategies.These include a Dirichlet-Neumann partitioned coupling method, a monolithic coupling method, and quasi-direct coupling method.Each of these is tightly-coupled and uses the same convergence criterion to define when the global FSI solution is converged during a given time step.
All three coupling methods use the same software infrastructure, system of equations assembly code, algebraic solver, and parameters.Other than their inherent differences in methodology, they are identical.
The benchmark proposed by Turek and Hron [2] was used to test the FSI coupling strategies developed in this work.We were able to achieve results which accurately reproduce those provided by Turek and Hron to within four percent at worst, generally within two percent.For this case, our results show that monolithic coupling and quasi-direct coupling are slightly more accurate, with respect to the validation data, than partitioned coupling and slightly better at capturing fine details along the interface.We also found that monolithic coupling and quasi-direct coupling are significantly faster than partitioned coupling when held to tight coupling tolerances.This substantiates Michler's hypothesis discussed earlier, that monolithic schemes can outperform partitioned schemes when held to the same tolerances [9], at least for the specific problem consider herein.
While this work aims to present a methodology for proper comparison between ALE FSI coupling algorithms, we acknowledge that the benchmark example considered is neither exhaustive nor addresses many of the other practical simulation issues.In the weak-coupling limit, we might expect that a partitioned or quasidirect approach would be more efficient than a monolithic approach.Large three-dimensional problems may also change the comparison results: as the number of degrees of freedom grows quickly, the cost of solving a large linear system becomes the main computational cost and scalability needs to be addressed.The choice of linear algebraic solvers could also greatly sway the study, as extremely efficient methods for solving each subproblem exist, but cannot be applied to the monolithic system of equations.The governing equations and the choice of discretization scheme will also undoubtedly affect any sort of computational cost comparisons.Finally, there exist numerous other ALE partitioned coupling methods (e.g.Robin-Robin coupling rather than Dirichlet-Neumann) that may be more accurate and efficient than the methods we investigated.These aspects of the ALE FSI solution can be studied within our numerical framework, and we leave them as future research efforts.ˆ, ,  The Jacobian matrix is formed by taking the the derivative of the residual vector with respect to the nodal unknowns: In the most general form, the discrete equations for the monolithic system are written as Although this system seems daunting, application of the residual vector components, from (69), (70), (71), (72), and (73), to the Jacobian definition in (74) shows that many of the terms are zero for all coupling algorithms, reducing the system to OPEN ACCESS WJM u 0 0 0 0 0 0 , 0 0 0 0 0 0 where the terms marked by C arise due to the con- forming velocity space in (58).

Discrete Quasi-Direct Equations
The quasi-direct discretization of ( 76) is similar to the monolithic discretization in that the tractions do not contribute on the interface.This results in the same terms being eliminated as were in the monolithic case.Furthermore, as the quasi-direct approach only solves Problems 1 and 2 monolithically, and iterates their solution with that of Problem 3, (76) can be separated into two smaller systems which reflect the coupled fluid-solid and seperate mesh: s s s s s s s s s s f s

Discrete Partitioned Equations
With partitioned coupling, Problems 1, 2, and 3 are all solved sequentially.This allows for (76) to be simplified into three independent smaller systems for the solid, fluid, and mesh:

Jacobian Component Calculations
The actual Jacobian calculations are simply derivatives of the residuals presented previously, and so they have been ommited.When evaluating the Jacobians, the following tensor calculus identities from Gurtin [12] are useful: ( ) In our calculations, we approximate the mesh (fluid domain) velocity within each time step as

Figure 1 .
Figure 1.Configurations of a body before and after deformation.
is the solid subdomain and F Ω is the fluid subdomain.Similarly, the domains's boundary ∂Ω is split into four parts: the Dirichlet boundaries for the fluid anddaries can be subdivided further to suit a given problem.The interface between the solid and fluid subdomains is Dirichlet and Neumann parts with the subscripts.The usefulness of this arises when defining proper finite-dimensional function spaces and setting boundary conditions for our discrete problems.

Figure 2 .
Figure 2. Domain Ω of an arbitrary two-dimensional fluidstructure interaction problem.The solid and fluid subdomains, Dirichlet and Neumann boundaries, and the fluidsolid interface are labeled.
using spatially constant material properties m The corresponding change to our finite dimensional vector spaces requires us to replace (39) and (44) with

Figure 4 ( 2 Q − 1
b) shows the analogous situation for the velocity and pressure fields in the fluid solver, where the difference in convergence rate is due to the Taylor-Hood

Figure 4 .
Figure 4. Error convergence plots for the two-field elasticity formulation (a) and Navier-Stokes formulation (b).In each case, the expected rates of error convergence are achieved, with respect to the L2-norm.This demonstrates that each component of the FSI solver is functioning correctly.

Figure 5 .
Figure 5. Turek and Hron [2] FSI benchmark domain with expanded view of cylinder and flag.Table 1. Turek and Hron [2] FSI benchmark domain dimensions and material properties.Parameter Symbol Value [m] Property Value

Figure 6 .
Figure 6.Visualization of the velocity magnitude and flag deformation at t = 16:5 [s].Solid displacement and fluid velocity are shown with their own color scales.The flag is shown as a wireframe to distinguish it from the fluid and the rigid cylinder.

2 Q 2 Q − 1
6, all subsequent results in the section were computed with elements for the solid displacement and velocity, Taylor-Hood

u
refer to the solid displacement at the current and previous FSI iteration, whether that iteration is fixed-point (partitioned) or Newton-Raphson (monolithic).Here,

Finally, Figure 10 Figure 8 .
Figure 8.The x-and y-displacement of point A on the flag tip.Our monolithic and quasi-direct schemes show nearly identical results to Turek and Hron [2], while our partitioned scheme is slightly different.

Figure 9 .
Figure 9.The drag and lift generated about the combined surface of the flag and cylinder.All three of our schemes show larger minimum values than Turek and Hron [2] do when calculating the drag, with otherwise good agreement.Our monolithic and quasi-direct schemes capture the complex shape of of Turek and Hron's lift results [2], while our partitioned scheme matches the global maximum and minimum values but not all local the minimums and maximums in-between.

Figure 10 .
Figure 10.Time study for partitioned, quasi-direct, and monolithic coupling over 20 simulation seconds (8000 timesteps).Our partitioned coupling scheme required over twice the amount of CPU time per time step compared to our quasi-direct coupling scheme, which in turn required approximately twice that of our monolithic coupling scheme. }