New iterative and staggered solution schemes for incompressible fluid‐structure interaction based on Dirichlet‐Neumann coupling

In the presence of strong added mass effects, partitioned solution strategies for incompressible fluid‐structure interaction are known to lack robustness and computational efficiency. A number of strategies have been proposed to address this challenge. However, these strategies are often complicated or restricted to certain problem classes and generally require intrusive modifications of existing software. In this work, the well‐known Dirichlet‐Neumann coupling is revisited and a new combined two‐field relaxation strategy is proposed. A family of efficient staggered schemes based crucially on a force predictor is formulated alongside the classical iterative approach. Both methodologies are rigorously analyzed on the basis of a linear model problem derived from a simplified fluid‐conveying elastic tube. The investigation suggests that both the robustness and the efficiency of a partitioned Dirichlet‐Neumann coupling scheme can be improved by a relatively small nonintrusive modification of a standard implementation. The relevance of the model problem analysis for finite element‐based computational fluid‐structure interaction is demonstrated in detail for a submerged cylinder subject to an external force.


Summary
In the presence of strong added mass effects, partitioned solution strategies for incompressible fluid-structure interaction are known to lack robustness and computational efficiency. A number of strategies have been proposed to address this challenge. However, these strategies are often complicated or restricted to certain problem classes and generally require intrusive modifications of existing software. In this work, the well-known Dirichlet-Neumann coupling is revisited and a new combined two-field relaxation strategy is proposed. A family of efficient staggered schemes based crucially on a force predictor is formulated alongside the classical iterative approach. Both methodologies are rigorously analyzed on the basis of a linear model problem derived from a simplified fluid-conveying elastic tube. The investigation suggests that both the robustness and the efficiency of a partitioned Dirichlet-Neumann coupling scheme can be improved by a relatively small nonintrusive modification of a standard implementation. The relevance of the model problem analysis for finite element-based computational fluid-structure interaction is demonstrated in detail for a submerged cylinder subject to an external force. or even prohibitive. In many cases, the iteration is inherently divergent. This phenomenon has been investigated by Causin et al, 2 Förster et al, 3 Joosten et al, 4 van Brummelen, 5,6 and others. It has been attributed to the fact that the inertia term in the solid solver does not account for the added fluid mass. The convergence of the subiteration can be improved to some extent, if relaxation is employed. It can be shown that the beneficial effect of relaxation is directly linked to the added mass phenomenon. The problem of added mass related instabilities is particularly notable for incompressible fluid flows. The stabilizing effect of fluid compressibility has been investigated by van Brummelen 5,6 and exploited by Farhat and coworkers. [7][8][9] The latter have proposed efficient staggered solution schemes for compressible fluid-structure interaction, which require the execution of the fluid and solid subsolvers only once per time step. Such schemes are often referred to as semiexplicit or weakly coupled. Pioneering work in this area of research was presented by Felippa and Park. [10][11][12][13] In recent years, a number of iterative or staggered schemes suitable for incompressible fluid-structure interaction with strong added mass effects have been proposed. Significant contributions to this field of research were presented by Gerbeau, Fernandez, Burman, and their coworkers [14][15][16][17][18][19][20][21][22] and also by Jaiman, 23,24 Wüchner 25 and others. 26 Some strategies are based on alternative types of boundary conditions, namely, Robin or Nitsche boundary conditions, while others are restricted to thin-walled structures or integrate the structural solver with certain steps of a segregated fluid solver. However, these methodologies are generally far more complex than Dirichlet-Neumann coupling and, with few exceptions, do not allow for the utilization of existing subsolver software. Some of them depend on the appropriate choice of artificial physical control parameters, which may be problematic for complex real-world applications.
In an earlier article, 27 the authors of this work have presented a provably second-order accurate staggered scheme for incompressible fluid-structure interaction that exploits the benefit of relaxation normally used in Dirichlet-Neumann subiteration. The computer implementation is nonintrusive and requires only a small modification of the Dirichlet-Neumann iterative scheme. The integration of existing subsolver software into the strategy is straightforward. It is shown that the critical amount of added mass is identical to the corresponding iterative solution scheme. However, while the performance of the subiteration deteriorates for increasing added mass effects, the proposed staggered scheme remains unconditionally stable and equally efficient up to the critical amount of added mass. In the article, 27 a basic one-dimensional linear model problem is used to support the proposed scheme with an analytical investigation. The authors have demonstrated the good performance of this staggered scheme in the context of arbitrary Lagrangian-Eulerian strategies 27 as well as immersed boundary methods 28,29 for fluid-structure interaction. The strategy has been applied to hydraulic valves and to the interaction between wind and a thin-walled solar collector. 30 For the sake of completeness, it is pointed out that the added mass related instabilities of partitioned solvers for fluid-structure interaction can be avoided by the employment of monolithic solution schemes based on Newton-Raphson-type strategies. The development of such solvers for fluid-structure interaction has been advanced by Tezduyar,31,32 Bazilevs,32,33 and their coworkers, who have applied them to large-scale problems including parachute and wind turbine simulations, but also by Heil, 34 Dettmer and Perić, [35][36][37][38] Wall, 39 and others. 40 The drawbacks of such strategies consist in the nonstandard derivatives required to ensure robustness and convergence, in the large and often prohibitive computational times as well as in the difficulties associated with the integration of existing subsolver software. In this context, the effect of neglecting certain cross-derivatives has been investigated by Dettmer and Perić. 41 They have also presented a methodology based on Schur complements, 35,36 which, for small number of interface degrees of freedom and linear fluid solver, competes with staggered schemes in terms of computational cost, but at the same time is insensitive to added mass. This methodology is ideally suited for fluid-rigid body interaction. 35

Contribution of this work
As mentioned above, Dirichlet-Neumann coupling can be implemented without intrusive software modification. While this is an outstanding advantage, the instabilities generally experienced in the presence of strong added mass effects render this type of coupling unfeasible for a wide range of practical applications. In a sense, the contribution of this work is to extract improved performance out of basic Dirichlet-Neumann coupling without resorting to intrusive software modification. This is done on a theoretical level on the basis of a linear model problem. The focus is clearly on the proposed staggered schemes. For the sake of completeness and comparative considerations, the analysis of the associated iterative strategies is presented in similar depth. The different schemes evolve from different application of numerical relaxation, namely, single-field and sequential two-field relaxation. One of the described staggered strategies, based on single-field relaxation, was presented in a second-order accurate setting in an earlier article by the authors 27 and has been successfully applied to a number of numerical examples. [27][28][29][30] In the following, the article's contribution is summarized under four headings: Model problem: It is shown that the linear one-dimensional model problem used in a number of the authors' previous publications 4,27,42 is equivalent to the model problem used in the classic article by Causin et al. 2 Essentially, the linearized problem of a fluid-conveying thin-walled elastic tube reduces to the superposition of fundamental solutions that are equivalent to one-dimensional linear oscillators with split inertia terms and potentially substantially different amounts of added mass. One aspect of this work is the investigation of the effect of the presence of both low and high added mass modes in a generic fluid-structure interaction problem on the performance of iterative solution schemes.
More insight into existing schemes: The model problem is used to assess the performance of iterative and staggered Dirichlet-Neumann solution schemes with velocity or force relaxation. The results obtained for the iterative strategies confirm those reported by other researchers, while those obtained for the staggered scheme provide more insight into the good performance of the algorithm proposed and tested in the authors' earlier article. 27 The effect of fluid damping on the solution schemes is also demonstrated.
Improved performance and stability with two-field relaxation: Iterative and staggered schemes with combined two-field relaxation are proposed. Namely, interface velocity and traction forces are relaxed in a sequential manner with the same relaxation factor. The analysis of the model problem shows that this significantly increases the critical amount of added mass for both staggered and iterative schemes. Furthermore, it accelerates the convergence of the iterative scheme in the presence of multiple added mass modes and increases the accuracy of the staggered scheme. It is suggested that the combined relaxation presents a nonintrusive modification of basic Dirichlet-Neumann coupling that significantly increases the efficiency and range of applicability of traditional iterative schemes and of the earlier staggered scheme. 27 Relevance of the model problem analysis: The presented solution schemes are implemented for a two-dimensional finite element-based example problem of incompressible fluid-structure interaction. A detailed investigation shows the close qualitative and quantitative agreement of the simulation with the conclusions drawn from the model problem. In order to allow for a well-founded conclusion, the example problem is run with a number of finite element formulations that employ different strategies to satisfy the incompressibility constraint.
All developments shown in this work are based on backward Euler time integration in the fluid and solid phases. The iterative and staggered coupling schemes discussed in this article are therefore restricted to first-order accuracy in time. A particular version of the staggered strategy covered in this work was presented in a second-order accurate setting in an earlier article by the authors. 27 More results on second-order accurate strategies will be presented in forthcoming publications.

Layout of this article
This article is organized as follows. In Section 2, the model problem is derived. In Section 3, Dirichlet-Neumann iteration and the staggered scheme are applied to the model problem and the tools for the analysis of the solution schemes are introduced. The effect of fluid damping on the solution schemes is investigated. Section 4 presents different relaxation strategies and their effect on the iterative and staggered schemes. In particular, velocity, traction force, and combined relaxation are considered. The results of the model problem analysis are summarized in Section 5. In Section 6, the proposed schemes are implemented for different finite element-based fluid solvers. A detailed comparison between the model problem analysis and the numerical examples is presented on the basis of a submerged rigid cylinder that is subject to an external force. Conclusions are drawn in Section 7.

Summary of the model problem presented by Causin et al 2
Consider a thin-walled elastic tube with circular cross-section that contains fluid. The fluid flow is governed by the Navier-Stokes equations and the tube wall is linear elastic. The radius R, the length L, the interface Σ, the fluid and solid domains Ω F and Ω S are defined as shown in Figure 1. • The fluid is incompressible and inviscid, • convection is negligible, • structural displacements are small, • the stiffness of the tube wall in axial direction is negligible.
Hence, the conservation of momentum and mass of the fluid reduce to where u and p denote, respectively, the fluid velocity and pressure fields, while f is the fluid density. The radial displacement of the tube wall is governed by where s , h s , a, and f represent, respectively, the wall density, thickness, stiffness, and the external traction force acting at the interface with the fluid. The interaction between the fluid and the solid phases must satisfy the requirements of kinematic consistency and equilibrium of the traction forces. Thus, the following equations hold where n={1,0} T denotes the outward normal unit vector on Σ. The boundary conditions at the orifices of the tube determine particular rather than fundamental solutions and are not relevant for the following analysis. The solution to the coupled problem given by Equations (1), (2), and (3) can be expressed as the superposition of fundamental solutions, that is, with The modified Bessel functions I 0 (•) and I 1 (•) are given by Using Equations (4) to (9) and defining k ∶= the problem given by Equations (1) to (3) reduces to a set of differential equations for the scalar amplitudes of the wall displacement̂k One observes that the acceleration of the solid mass s h s requires the simultaneous acceleration of the fluid mass f k . Hence, the expression f k is identified as the added mass associated with the kth fundamental solution. The study of Equation (11) reveals that the added mass effect is weak for large values of k (short wave lengths ), while it is strong for small values of k (long wave lengths ). The schematic representation of a number of different mode shapes in Figure 2 illustrates that longer wave length responses are associated with more fluid inertia than shorter ones. It follows that max = 1 . Furthermore, in the limit of slender tubes with L > >R, one obtains max ≈ By combining Equations (12) and (13) one observes that high added mass modes are associated with low frequencies in the time domain and vice versa. Hence, in the context of time discretization, a time step size Δt that is large for a low added mass/high frequency mode is small for a high added mass/low frequency mode.
Equations (12) and (13) can be rewritten as Introducing a viscous term with a damping coefficient in Equation (17) renders the model problem used by the authors in their earlier article, 27 that is, It is noted that the parameter is equal to the ratio of the solid mass over the sum of the solid mass and the added mass of fluid. Thus, it follows that 0 ≤ ≤ 1. Problems involving strong added mass effects are characterized by values of close to zero, while small added mass effects are associated with values of close to one. In the limit of slender tubes characterized by Equation (14), the smallest value of for a tube with L = 100 mm, R = h s = 1 mm, and s = f is obtained as 4.93 × 10 −4 . Hence, it is clear that many applications including, for instance the modeling of components of the vascular system, are associated with very small values of .

2.3
The challenge arising from multiple modes The representation of the response of the tube to mechanical excitation such as, for instance, an initially deformed configuration or time variable inflow boundary conditions, generally requires the superposition of a large number of fundamental solutions with a range of different wave lengths . Hence, the fluid-structure interaction problem of the tube consists of a number of systems described by Equations (18) and (19) that span a wide range of values of . As is generally the case with large complex and nonlinear mechanical systems, the modal decomposition of a realistic fluid-structure interaction problem is impossible or, at best, impractical. Any viable solution strategy for the coupled system must work simultaneously on all modes, since it is impossible to tune any solution control parameters to each individual mode. For a viable methodology one set of such solution parameters must work for the whole system. This is important for the investigation in Section 3.

Time discretization of the model problem
The investigations in Sections 3 and 4 are based on the standard implicit first-order accurate and unconditionally stable backward Euler discretization of the time derivatives in Equations (18) and (19). Hence, the following approximations are employedu where the subscript n identifies all quantities associated with the time instance t n and Δt = t n+1 − t n . Using this in Equations (18) and (19) renders, after some manipulation, In the following, Equation (21) is referred to as the solid solver and Equation (23) is termed the fluid solver.
As presented in Equation (21), the solid solver returns the interface velocities rather than the displacements. It is important to note that this choice has been made exclusively to allow for a shorter presentation of the analysis in Sections 3 and 4. The results of the analysis are not affected by this choice. In practice, any standard finite element software for solid dynamics returns both displacements and velocities.

Dirichlet-Neumann iteration
The most intuitive strategy for the resolution of the strong coupling between the fluid and the solid phases is commonly referred to as iterative Dirichlet-Neumann coupling. 1 Based on an initial guess of the interface configuration the fluid solver is executed and the traction field exerted by the fluid on the interface is passed to the solid solver where it is used as external loading, that is, Neumann boundary condition. Next, the interface displacement and velocity fields obtained from the solid solver are transferred to the fluid solver and used as Dirichlet boundary conditions. This process is repeated iteratively. In terms of the model problem given by Equations (18) and (19) this procedure can be expressed as or, shorter,u where the superscript i identifies the iteration step. Recalling Equations (21) and (23) renderṡ The factor A is known as the amplification factor. Since A multiplies the numerical error in each iteration step, convergence requires that |A| < 1. Convergence is poor if |A| is close to one and fast if |A| is close to zero. To investigate the performance of the scheme in Equation (25), it is convenient to consider the small and large time step limits of A. One obtains It is shown in Section 3.3.1 that the study of the limit cases is sufficient for the assessment of the convergence unless the damping in the fluid is excessive. The condition |A 0 | < 1 requires while |A ∞ | < 1 is satisfied in any case.
Recalling the nature of the parameter as described in Section 2.2, it is observed that, in the small time step regime, the scheme only converges if the added fluid mass is smaller than the structural mass. This condition must be met for each mode of the model response. Furthermore, convergence is slow if there exists at least one mode with a value of close to 1/2. On the other hand, if all values of are close to 1, convergence is fast. For larger time steps the dependency of convergence on the added mass phenomenon weakens and vanishes as Δt → ∞.

Staggered solution strategy
Based on the backward Euler scheme in the fluid and solid subsolvers (see Section 2.4), the following staggered solution procedure is considered: This procedure corresponds closely to the first step of the Dirichlet-Neumann iteration described in Section 3.1. While the latter can be started with an initial guess of the traction forces or the interface configuration without affecting the performance and the critical condition (29), it is essential for stability that this staggered scheme is based on a traction force predictor rather than a predictor of the interface configuration. The solid solver requires the displacement u n as an input argument. Hence, Step 2 of the staggered scheme includes the calculation of the displacement u n + 1 , based on the backward Euler formula, in readiness for the next time step.
A key feature of the scheme (30) is the explicit coupling of the implicit solid and fluid subsolvers. The solid and the fluid phases share the interface kinematics, whereas the equilibrium of the traction forces at the interface is satisfied only approximately. The solid interface tractions are equal to the predictor f n , while the fluid interface tractions are equal to f n + 1 .
The staggered scheme (30) can be essentially reduced to the relationṡ Due to the linearity of the expressions, the above system of equations can be written in matrix form as with For stability in the time domain, the spectral radius (A) = max(| j |, j = 1, 2, 3) of the amplification matrix A, where j are the eigenvalues of A, must not exceed the value of 1. Similarly to the analysis in Section 3.1 the small and large time step limits are considered. One obtains F I G U R E 3 Spectral radius for the staggered scheme for = 1, = 0.01, and different values of . The spectral radius of the strongly coupled scheme is obtained from Equation (43). It is independent of and shown as a reference It is argued in Section 3.3.2 that the analysis of the limit cases is sufficient unless the numerical damping in the fluid is excessive. The requirement for stability | 0 | ≤ 1 renders while | ∞ | ≤ 1 is met in all cases. Therefore, if the condition in (38) is met, the Scheme (30) is unconditionally stable. Except for the inclusion of the critical value = 1∕2 into the range of admissible values of , the condition in (38) is identical to the one obtained for the iterative Dirichlet-Neumann iteration in (29). Thus, with regard to the added mass phenomenon, the two strategies have the same range of applicability. However, while a value of close to 1/2 jeopardizes the convergence of the Dirichlet-Neumann iteration and hence drastically increases the computational cost, the same value of does not affect the performance of the staggered scheme. This is a crucial advantage. The stability of the staggered scheme for any time step size is further illustrated in Figure 3 in terms of the spectral radius of the amplification matrix. It is noted that, despite its semiexplicit nature, the staggered scheme in (30) behaves very similarly to the strongly coupled solution. Interestingly, it maintains the high frequency limit ∞ = 0 of the backward Euler schemes in the subsolvers.
Next, the accuracy of the staggered scheme is investigated. The procedure adopted for this purpose is the same as the one used by Chung and Hulbert 43 and in earlier work by the authors. 27,44 First, the exact solution of the system of Equations (18) and (19) is determined and written in the format of an amplification matrix The eigenvalues of A exact are obtained as Inserting them into the characteristic polynomial of the numerical amplification matrix A from Equation (35) and expanding the resulting expression in terms of powers of Δt renders If all terms involving factors of , which is generally small, are neglected, this reduces to The order of accuracy of the scheme is obtained by reducing the lowest exponent of Δt by 2. Hence the weakly coupled staggered scheme is first-order accurate and therefore maintains the accuracy of the backward Euler scheme in the subsolvers on the level of the overall coupled problem.
For a detailed comparison, the temporal accuracy of the solution obtained with the Dirichlet-Neumann iteration described in Section 3.1 can be found by performing a similar analysis. This requires the derivation of the two-dimensional numerical amplification matrix associated with the strongly coupled problem. One obtains and or, neglecting higher orders of , The comparison of Equations (42) and (45) shows that, depending on the value of with 1 2 ≤ ≤ 1, the approximation error in time of the strongly coupled solution is between 1 and 3 times smaller than for the staggered one. However, by reducing the time step size by a small factor, the staggered scheme will outperform the accuracy of the Dirichlet-Neumann iteration and still be significantly more efficient.

The effect of the fluid damping
The analyses in Sections 3.1 and 3.2 of the performance of the Dirichlet-Neumann iteration and of the staggered scheme is based on the time step limits Δt → 0 and Δt → ∞. In the following paragraphs, it is shown that this is indeed sufficient.

Dirichlet-Neumann iteration
Consider the case without damping in the fluid, that is, = 0. The amplification factor in Equation (26) reduces to Clearly, the function A(Δt) transitions smoothly and monotonously from Δt = 0 to Δt → ∞ without any stationary points. It follows that, for = 0, the consideration of the small and large time step limits in the convergence analysis in Section 3.1 is sufficient. Hence, provided that the critical condition in (29) is met, any values of |A(Δt)| > 1 must be due to the effect of the fluid viscosity.
Reintroducing the damping > 0 in the amplification factor, differentiating it with respect to Δt and equating the result to zero renders the position of the minimum value A min . Solving A min = −1 for gives Hence, if there is no added mass in the system, that is, = 1, even excessive fluid damping with values of close to 1 does not jeopardize the convergence of the scheme. If is close to the critical value of 0.5, where convergence is slow, values of > 0 can cause divergence. However, a numerical study shows that, for moderate values of , this occurs only in the very close vicinity of = 0.5. For instance, for = 0.51, one obtains max = 0.1414, which is still a very substantial and unlikely amount of damping. Examples of A(Δt) for different values of and are visualized in Figure 4.
Therefore, one concludes that the effect of fluid damping on the convergence of the Dirichlet-Neumann iteration is negligible in most cases and the consideration of the small and large time step limits in the convergence analysis is sufficient. However, close to the critical amount of added mass, high physical or artificial fluid viscosities may trigger convergence issues.

Staggered scheme
An analysis similar to Section 3.3.1 requires the analytical expressions for the spectral radius of the amplification matrix in Equation (35). This gives rise to lengthy expressions and is omitted. Instead, analogously to Figure 4, the spectral radius of the amplification matrix is shown in Figure 5 as a function of Δt and for different values of and . The figure suggests that the effect of fluid damping on the stability of the staggered scheme is similar to its effect on the convergence of the Dirichlet-Neumann iteration: For = 0, the spectral radius transitions monotonously from Δt = 0 to Δt → ∞. For larger values of , a relative maximum develops. Numerical evidence suggests the same critical value for as given by Equation (47). Hence, very large fluid viscosities may cause stability issues close to the critical amount of added mass. This effect is negligible and the analysis of the staggered scheme can be based on the small and large time step limits.

Summary
It has been established that the schemes presented in Sections 3.1 and 3.2 are unsuitable for the simulation of problems with strong added mass effects where at least one of the modes features < 0.5: The Dirichlet-Neumann iteration fails to converge and the staggered scheme is unstable. It has also been shown that, in the regime of smaller added mass effects with ≥ 0.5, the staggered scheme is substantially more efficient than the Dirichlet-Neumann iteration without jeopardizing the stability or accuracy in time. This is especially the case for problems involving one or more modes with values of close to 0.5. In this regime the Dirichlet-Neumann iteration requires a large number of iteration steps since the amplification factor is close to 1.
In the next section, it is investigated to what extent basic modifications of the two schemes can remove the limitations posed by the added mass phenomenon.

MODIFICATION OF SOLUTION SCHEMES
The following subsections present modifications of the Dirichlet-Neumann iteration and of the staggered scheme presented in Sections 3.1 and 3.2, respectively. In particular, the effects of velocity and of traction force relaxation are investigated.

Dirichlet-Neumann iteration
Introducing relaxation of the interface velocity in the Dirichlet-Neumann iteration as described by Equation (24) renders or, in short,u where is the relaxation factor and can be chosen freely to optimize the rate of convergence. It is noted that, by using Equation (20), this can be rephrased as relaxation of the displacement u n + 1 .
Recalling Equations (21) and (23) the iteration can be expressed aṡ where the factor A introduced in Equation (26) is reused. Convergence requires that the modified amplification factor A ′ satisfies |A ′ | < 1. Similarly to Section 3.1, the small and large time step limits are considered and one obtains This result has also, in different notation, been reported by Causin et al 2 and others 3,4,6 . It is noted that the condition in (53) implies (54). Setting = 1 correctly recovers the inequality in (29). Furthermore, there exists an optimal relaxation factor * that renders A ′ = 0 and therefore requires only one iteration step to obtain the exact solution. In the small time step limit * = and, in the large time step limit, * = 1.

F I G U R E 6
Dirichlet-Neumann iteration; visualization of the overall amplification factorÂ( ) as given by Equations (51)  Considering the regime of small time steps and recalling the nature of the parameter , one observes that problems with large added mass effects require the choice of a value of close to zero, while for small added mass effects the interval of admissible values for is much larger. However, for all with 0 < ≤ 1, can be chosen such that convergence is ensured and there exists an optimal value * which renders the solution after a single iteration step.
As described in Section 2, realistic fluid-structure interaction problems involve a large spectrum of values . For convergence, the relaxation parameter for the overall problem must be chosen such thatÂ Furthermore, the optimal relaxation factor is found from If min and max denote the smallest and the largest values of present in a particular fluid-structure interaction problem, then one can deduce that convergence is conditional on and, in the small time step limit, the optimal relaxation factor and amplification factor are, respectively, * = 2 max min max + min This is illustrated in Figure 6. As described in Section 2.1, a time step size which is small for the low frequency high added mass mode with min is large for the high frequency low added mass mode with max . Thus, by using A ′ 0 for min and A ′ ∞ for max , one obtains, instead of Equation (59), * = 2 min 1 + min andÂ( * ) = 1 − min 1 + min .
As an example, consider two modes with 1 = 0.99 and 2 = 0.01. Convergence requires ∈]0, 0.02[. Equation (59) renders * = 0.0198 andÂ( * ) = 0.98. It follows that even the optimal amplification factor requires 228 iterations in order to reduce the solution error by a factor of 100. Also note the close vicinity of * to the critical value of 0.02. The values obtained from Equation (60) are slightly more disadvantageous. Hence, it is clear that a large spectrum of added mass modes substantially jeopardizes the convergence of the Dirichlet-Neumann iteration and potentially makes the strategy unviable.

Staggered solution strategy
The integration of interface velocity relaxation into the staggered scheme described in Section 3.2 requires the distinction of solid and fluid interface kinematics. Therefore, the superscripts s and f are reintroduced to denote the quantities associated with the solid and the fluid phases, respectively. The staggered scheme (30) This procedure can be rewritten in matrix form as where, for brevity, = + 2 Δt 2 and = 1 − + 2 Δt. Analysis similar to Section 3.2 renders It follows that the scheme is unconditionally stable provided In the context of multiple solution modes with a range of values of , the relaxation parameter must satisfy which is identical to the condition for convergence of the Dirichlet-Neumann iteration (58), except for the critical value lying inside the admissible interval of . The minimization of the term |1 − ∕ | which renders the optimal performance of the Dirichlet-Neumann iteration given in Equation (59) is not relevant in the context of the staggered scheme. Inserting the eigenvalues (40) of the exact amplification matrix (39) into the characteristic polynomial of the amplification matrix (63) renders If all terms involving factors of , which is generally small, are neglected this reduces to Clearly, for = 1 the expression in Equation (42) is recovered. The presence of the term 1∕ in Equation (69) indicates that the approximation error grows as added mass effects increase, that is, the results of the simulation deviate more strongly from the strongly coupled computation as → 0.

4.2
Force relaxation

Dirichlet-Neumann iteration
Introducing relaxation of the traction force in the Dirichlet-Neumann iteration as described by Equation (24) renders Formulating the second equation at (i) instead of (i + 1) and inserting it into the first yields a shorter representation of the scheme, Recalling Equations (21) and (23) the iteration can be expressed as where the factor A introduced in Equation (26) is reused. The comparison of Equations (72) and (50) shows that force relaxation yields the same amplification factor A ′ as displacement relaxation. Hence, the performance of the iterative schemes described by Equations (49) and (71) is identical. It has been investigated in detail in Section 4.1.1.

Staggered solution strategy
The integration of interface displacement relaxation into the staggered scheme described in Section 3.2 is straightforward and leads to the following scheme: 1. set traction force predictor f P n+1 = f n 2. execute solid solver 5. proceed to next time step n ← n + 1 .
This procedure can be rewritten in matrix form as (75) Analysis of the small and large time step limits renders Hence, 0 and ∞ are identical to the staggered scheme based on displacement relaxation and the results of the stability analysis described in Section 4.1.2 apply. The accuracy analysis also renders exactly the same error terms as in Equations (68) and (69). A second-order accurate version of this staggered scheme based on the generalized-method has been published and thoroughly tested by the authors. [27][28][29]

Dirichlet-Neumann iteration
Using relaxation for both the traction force and the velocity in the iteration described by Equation (24) renders The same relaxation factor is applied to the force and to the velocity. The employment of two different values does not render any benefit and is not considered here. The above system of equations can be written in matrix form as The analysis of the small and large time step limits renders , It follows from 0 < 1 that F I G U R E 7 Single-field relaxation: Contour plot of the amplification factor A 0 for the Dirichlet-Neumann iteration and of the spectral radius 0 of the amplification matrix for the staggered scheme in the small time step limit Δt → 0

F I G U R E 8 Combined
relaxation: Contour plot of the spectral radii 0 of the amplification matrices for the Dirichlet-Neumann iteration and the staggered scheme in the small time step limit Δt → 0 and from ∞ < 1 that The latter represents the same large time step convergence criterion that applies to the single-field relaxation of velocity or traction force as shown in Sections 4.1.1 and 4.2.1. The small time step regime requires more attention. Figure 8 shows the contour plot of the spectral radius 0 . The comparison with Figure 7 illustrates that convergence is achieved in a larger region of the -plane. It follows from the critical values of the relaxation factor as given in the conditions in (53) and (83), that combined relaxation allows for substantially larger values of in the high added mass regime, where is close to zero. In fact, for combined relaxation, the ratio of the critical value of over tends to infinity for → 0, whereas it is equal to 2 for single-field relaxation. As an example, = 0.01 requires < 0.02 or < 0.1827 for, respectively, single-field or combined relaxation.
The optimal value of the relaxation factor in the small time step limit is obtained as Hence, cannot be chosen such that 0 disappears. This is illustrated in Figure 8. Yet, this is disadvantageous only in the case of a single-solution mode. As discussed in Section 2.3, physically relevant fluid-structure interaction problems typically feature multiple solution modes that range from low to high added mass effects and, therefore, do not share the same optimal value of . Considering the example of min = 0.01 and max = 0.99, one obtains * = 0.1818 and 0 ( * ) = 0.8182. It follows that 23 iteration steps are required to reduce the approximation error by a factor of 100. This compares very favorably to 228 steps that are required by single-field relaxation as shown in Section 4.1.1.

Staggered solution strategy
Combining the relaxation of the interface velocity and traction forces, that is, merging the algorithms in (61) and (73), one obtains the following staggered scheme: 1. set traction force predictor f P n+1 = f n 2. execute solid solver 6. proceed to next time step n ← n + 1 .
This procedure can be rewritten in matrix form as where, for brevity, = + 2 Δt 2 and = 1 − + 2 Δt. The analysis of the small and large time step limits renders It follows from 0 ≤ 1 that and from ∞ ≤ 1 that Hence, the critical values of the relaxation factor correspond exactly to those obtained for the iterative scheme in Section 4.3.1 and are larger than their counterparts associated with single-field relaxation.
Inserting the eigenvalues (40) of the exact amplification matrix (39) into the characteristic polynomial of the amplification matrix (88) renders If all terms involving factors of , which is generally small, are neglected, this reduces to The error term includes the expression 2 ∕ . Even though the critical is much larger for this scheme than it is for single-field relaxation, the condition in (91) states that smaller values of still require the choice of a smaller value of . Thus, the term 2 ∕ does not cause the loss of accuracy for small values of that is associated with the expression 1∕ in the leading order error term for the staggered scheme with single-field relaxation in Equation (69).

SUMMARY OF MODEL PROBLEM ANALYSIS
The results of the analysis presented in Sections 3 and 4 are summarized in Table 1. The table shows the critical conditions required to ensure convergence or stability of, respectively, the iterative or the staggered scheme. It also shows the expressions that have been derived for the amplification factors and spectral radii in the small and large time step limits. For the staggered schemes, the leading order error terms are provided. The correspondence of the critical amounts of added mass between the iterative and the proposed staggered schemes is evident. It must be recalled, however, that the performance of the iterative scheme deteriorates significantly as the critical added mass is approached, while the efficiency of the staggered schemes is not affected until the critical added mass is exceeded.
For both the iterative and the staggered schemes, velocity and traction force relaxation render identical results. The model problem analysis suggests that combined relaxation allows for the resolution of substantially larger added mass effects.
In the context of the staggered schemes, combined relaxation has two advantages over single-field relaxation: It features larger critical values of , especially in the regime of strong added mass effects, and the approximation error (the deviation from the strongly coupled solution) does not grow proportionally to the inverse of .

NUMERICAL EXAMPLE
In this section, the agreement of the model problem analysis of Sections 2 to 5 with different finite element models of a basic two-dimensional fluid-structure interaction system is investigated. In order to allow for meaningful conclusions, different finite element-based fluid subsolvers for the incompressible Navier-Stokes equations are considered and presented in Section 6.1. In Section 6.2, these fluid subsolvers are integrated in the staggered schemes and used to simulate the problem of a rigid cylinder submerged in two-dimensional fluid flow and pulled by a prescribed external force. In Section 6.3 the iterative schemes are applied to the same problem.
It is important to note that, despite the simple setup described in Section 6.2, the problem of the submerged rigid cylinder features strongly coupled fluid-structure interaction and fully resolved incompressible fluid flow and is not based on any simplifying assumptions. Yet, it closely resembles the single-mode model problem analyzed in the preceding sections. Hence, if the model problem analysis is indeed relevant for real applications of fluid-structure interaction, it must predict the performance of the numerical strategies for this test case with a good degree of accuracy. This is verified in the following subsections.
The implementation of the iterative and staggered schemes can be based on the algorithms shown in Equations (48), (70), (78) and (61), (73), (86), respectively, and is straightforward. The fluid and solid subsolvers used in this section employ the same backward Euler time discretization scheme as the model problem. Note that the convective velocity is evaluated at the previous time instant, thereby avoiding nonlinearity and allowing for fully linear fluid solvers.

6.1
Fluid flow solver

Governing equations
Consider a domain Ω ⊂ R d (d ≤ 3) with boundary Γ which is separable into Dirichlet and Neumann subsets, Γ D and Γ N . The velocity field u and pressure field p are described by the incompressible Navier-Stokes equations where f is the body force vector, u D is the velocity prescribed on Γ D , and n is the outward normal to Γ. The velocity of the reference frame is denoted by the vector field v. The symmetric gradient is expressed as ∇ s (•) ∶= 1 2 (∇(•) + ∇ T (•)).

Finite element method: Mixed Taylor-Hood elements
In order to satisfy the LBB stability condition, a Taylor-Hood element is chosen for the spatial discretization (see Figure 9). Thus, the velocity interpolation is piecewise quadratic, while the pressure interpolation is piecewise linear.

F I G U R E 9 Two dimensional Taylor-Hood interpolations
The weak form reads: where  h ,  h , and  h are the appropriate finite element spaces of piecewise continuous quadratic and linear basis functions.

Finite element method: Stabilized SUPG/PSPG elements
The velocity and the pressure fields are both interpolated with the same piecewise continuous linear basis functions. Following the original work by Tezduyar et al, 45 the SUPG/PSPG stabilized weak form reads: The scalar stabilization parameters are defined as where h denotes the characteristic element size and is a scaling factor which allows for different weighting of the pressure stabilization. Note that, in order to avoid nonlinearity, the norm of the convective velocity in the expression for au u is taken from the previous time instant. Much more sophisticated stabilization expressions have been proposed by Tezduyar and coworkers. However, for the purpose of the present work, the above terms, which work well in the regime of laminar flow, are sufficient.

Finite element method: Projection scheme
For the projection scheme considered in this work, the velocity and the pressure fields are interpolated, respectively, with piecewise quadratic and piecewise linear basis functions as shown in Figure 9. Following a decoupling strategy-based on the original nonincremental projection method by Chorin, 46 the semidiscrete form of the Navier-Stokes equations is partitioned into three steps: Step 1: This is essentially the momentum equation after the application of the backward Euler time integration scheme, where the pressure gradient has been removed. It is solved for the so-called intermediate velocityũ n+1 . Note that u n+1 satisfies the Dirichlet boundary condition in Equation (97). The weak form reads: Findũ h n+1 ∈ h , such that for where h and h represent the appropriate finite element spaces based on piecewise quadratic basis functions.
Step 2: Taking the divergence of the terms in Equation (105) and making use of Equation (106) renders the so-called pressure Poisson problem from which the pressure p n + 1 is obtained. The weak form of step 2 is written as follows: where  h represents the appropriate finite element space based on piecewise linear basis functions. In Step 3, the physical velocity is recovered from Equation (105) as Step 3: It is noted that the evaluation of the right-hand side is not trivial since the velocities are nodal quantities, while the pressure gradient is available only in Gauß points. In the present work, Equation (109) is solved weakly, which requires a global matrix system. Finally, the traction forces are computed by substituting u n + 1 and p n + 1 into the weak form of the momentum equation (95) and evaluating the reaction forces at the mesh boundary.
It is generally possible to eliminate the end-of-step velocity u n + 1 and thereby obtain a more efficient procedure. However, the traction forces obtained from Equation (104) are less accurate than those obtained by the procedure described above.
More information on segregated solution schemes for the incompressible Navier-Stokes equations is provided in a recent article 47 by the authors and many references therein. In their article, the authors employ a model problem similar to the one used in the present work to provide insight into such methodologies.

Submerged cylinder subject to external force: Staggered schemes
Consider a rigid solid cylinder of density s submerged in a square domain of fluid with density f = 1 and viscosity = 0.01 as shown in Figure 10. The cylinder is free to move in x-direction and subject to an external force F ext for t < 10 0.6 for t ≥ 10.
Hence, the external force accelerates the cylinder until it is balanced by the increasing drag force. The motion of the cylinder at terminal velocity corresponds to flow around a fixed cylinder at Reynolds number Re ≈ 100. The finite element F I G U R E 10 Submerged cylinder: Geometry and mesh with 12 258 elements mesh used in all simulations is also shown in Figure 10. The number of nodes and of degrees of freedom depends on which of the fluid solvers described in Section 6.1 is employed. This problem is chosen due to its simplicity and due to the fact that it is well known that the amount of added fluid mass corresponds to the mass of the fluid displaced by the cylinder. 48 Arguing that the solution to the problem can be described by a single mode in the sense of Sections 2 to 5, it follows that a representative value of can be obtained from It can be varied by using different values for the solid density s . The solid solver resolves only the inertia of the cylinder and, therefore, reduces toḋ whereḋ s and F represent, respectively, the cylinder velocity and the interface traction force in x-direction, while R is the radius of the cylinder. The fluid solver takes as input the velocity of the cylinder and applies it as a Dirichlet boundary condition at the interface boundary and also as reference frame velocity throughout the fluid domain. Hence, the fluid mesh moves with the same velocity as the cylinder and does not deform. Similarly to the model problem of Sections 2 to 5, the solid displacement at the interface does not need to be transferred to the fluid domain. The fluid solver returns the interface traction force which is obtained by summing the nodal reaction forces on the interface boundary. Periodic vortex shedding develops in all simulations shortly after the terminal velocity has nearly been reached. This is, however, not relevant for the current investigation.

Robustness
In order to establish critical values for the relaxation parameter , the simulations are performed for different combinations of and . Figure 11 shows  2. The accuracy of the agreement with the analytical expressions varies with the type of fluid solver employed. In particular, one observes that the largest deviation from the analytical expression amounts to approximately 10% and is obtained for single-field relaxation with Taylor-Hood elements. On the other hand, the critical values of obtained for the projection scheme agree very accurately with analytical values for all values of . The stabilized finite element method with = 1 recovers the analytical expressions quite accurately, but is inherently unstable for = 0.01. and combined relaxation. Reducing to 0.5 allows for the stable simulation of = 0.01, but causes a slightly larger deviation of the critical values of for all values of considered. Experiments with smaller time steps also require that is reduced. Small spurious oscillations observed in the response of the Taylor-Hood elements with combined relaxation disappear if the pressure boundary condition at the outflow boundary is removed. The methodology is then found to be stable in the range of small values of and Δt with critical values of close to the values predicted by the model problem.
Crucially, Observation 1 confirms the relevance of the model problem and the model problem analysis for multidimensional computational fluid-structure interaction.
Observation 2 suggests that the accuracy of the agreement of the critical value of with the model problem analysis depends on the mechanism used in the fluid solver to satisfy the incompressibility of the fluid. The effect of changing the factor in the stabilized finite element method shows this very clearly. This is consistent with the well known sensitivity of the robustness of iterative partitioned solvers to the treatment of the fluid incompressibility.

Accuracy
The diagrams in Figures 13 to 15 show the evolution of the interface traction force and the cylinder velocity for three sets of simulations which have all been performed with stabilized finite elements. The results obtained from Taylor-Hood elements and from the projection scheme deviate very little and are not shown in order to avoid repetition. The cases investigated are summarized in Table 2.
The following observations are made from Figures 13 to 15: 1. For = 0.2 (Case 1, Figure 13), velocity relaxation, traction force relaxation and combined relaxation render almost identical results. 2. For = 0.01 and Δt = 0.1 (Case 2, Figure 14), velocity relaxation renders accurate traction forces and fluid interface velocity, while the solid interface velocity features a strong artificial maximum before it approaches the terminal velocity. Traction force relaxation, on the other hand, yields an accurate solid traction force, but results in a poor fluid traction force and an artificial velocity maximum prior to approaching the terminal velocity. Combined relaxation yields accurate results for the solid force and the fluid velocity, while the results obtained for the fluid force and the solid velocity are not as accurate but of good quality. Altogether it is evident that the accuracy of the staggered scheme with combined relaxation is superior to single-field relaxation. 3. For = 0.01 and Δt = 0.01 (Case 3, Figure 15), the artificial effects in the solutions obtained for Δt = 0.1 (Case 2) from single-field relaxation are substantially reduced. For the combined relaxation method , the fluid and solid velocities and forces coincide.
Observation 1 confirms that, in the low added mass regime with large values of , all three variants of the staggered scheme are viable. Observation 2 suggests that, in the high added mass regime with small values of , the combined  Table 2). The fluid velocities and solid forces computed with combined relaxation are shown in all diagrams as reference solutions relaxation outperforms single-field relaxation. Observation 3 confirms that the poor accuracy of single-field iteration in the large added mass regime is effectively reduced by using smaller time steps. All observations are consistent with the model problem analysis in Section 4.

Submerged cylinder subject to external force: Iterative schemes
The same submerged cylinder is considered as in Section 6.2. The fluid solver is based on stabilized SUPG/PSPG elements with = 1 and Δt = 0.1. The parameters and are varied. The simulation is run until t = 5. The fluid-structure interaction is resolved by Dirichlet-Neumann iteration. The iterations are stopped once the absolute change in the interface velocities and in the interface forces, respectively, is smaller than 10 −3 . For a small problem such as the one under consideration, this is a relatively large tolerance. Figures 16 to 18 show the number of iterations required by the different relaxation strategies. The optimal relaxation factors that render minimum numbers of iteration steps are extracted from Figures 16 to 18 and compared to the analytical expressions derived from the model problem in Figure 19.  Table 2). The fluid velocities and solid forces computed with combined relaxation are shown in all diagrams as reference solutions The following observations are made: 1. Figure 19 shows that the optimal values for the relaxation factor agree well with those obtained from the model problem analysis. 2. The qualitative behavior of the iteration numbers shown in Figures 16 to 18 agrees well with the model problem analysis: For single-field relaxation and small values of , the number of iterations rises sharply on either side of the optimal relaxation factor * . For larger values of iteration numbers increase more slowly as deviates from * . This corresponds to the narrowing interval of admissible values of as approaches zero and is visualized in Figure 7. The performance for = * is similar for all values of , approximately two iteration steps per time step. This is due to the fact that for any value the optimal relaxation factor results in an amplification factor of zero in the small time step limit. For combined relaxation, the iteration numbers rise sharply as exceeds * while they rise more slowly and follow the same path as is reduced. The minima of the iteration numbers increase as becomes smaller. This behavior is in perfect agreement with the variation of the spectral radius shown in Figure 8. 3. For an application involving the combination of low and high added mass modes, the combined relaxation results in fewer iteration steps than single-field iteration, confirming the claim in Section 4.3.1: Assume a problem involving  Table 2). The fluid velocities and solid forces computed with combined relaxation are shown in all diagrams as reference solutions two cylinders, one corresponding to 1 = 0.05 and the other corresponding to 2 = 0.5. The relaxation factor has to be chosen such that the iteration converges for both values of . It is observed from Figures 16 to 18 that the minimum numbers of total iteration steps are obtained as 730 ( ≈ 0.07), 600 ( ≈ 0.07), and 500 ( ≈ 0.32) for, respectively, velocity relaxation, traction force relaxation and combined relaxation. In order to avoid excessive numbers of iterations the problem has not been run for values < 0.05. However, it can be deduced from the Figures 16 to 18, that for smaller minimum values of the difference between the numbers of iteration steps required is much more significant.

CONCLUSIONS
The achievements and conclusions of this work may be summarized as follows: Model problem derivation and relevance: The equivalence of the linear one-dimensional model problem used earlier by the authors 4,27,42 with the work by Causin et al 2 has been shown and a model problem based setting for the analysis of computational strategies for fluid-structure interaction has been described. The accurate agreement of the findings from the model problem analysis with the finite element simulation of a submerged rigid cylinder accelerated by an external  (85) force has been shown in detail, including the consideration of different finite element-based solvers for incompressible fluid flow, namely, stabilized, LBB-stable, and segregated methods.
More insight into existing schemes: The inadequacy of standard iterative relaxation techniques to support convergence in the simultaneous presence of low and high added mass modes has been demonstrated. A detailed investigation of a backward Euler based version of the staggered scheme proposed earlier by the authors 27 has been presented, including the derivation of the critical relaxation factor and the leading order error term. The effect of fluid viscosity on convergence has been investigated.
Improved performance and stability with two-field relaxation: Combined velocity and traction force relaxation has been shown to allow for substantially larger relaxation factors and to improve convergence of the iterative scheme in the case of multiple added mass modes. For the staggered scheme with standard relaxation, the leading-order term of the approximation error increases for stronger added mass effects. This is not the case for combined relaxation, which allows for the simulation of stronger added mass effects with larger values of the relaxation parameter.
Staggered schemes: A family of staggered schemes with force predictor and Dirichlet-Neumann coupling has been proposed and proven to be unconditionally stable below a critical added mass limit. This limit has been shown to be identical to the corresponding iterative schemes. While the performance of iterative schemes deteriorates in the vicinity of the critical relaxation factor, the performance of the staggered scheme is not affected. This is especially beneficial in the presence of multiple added mass modes, typically encountered in fluid-flexible structure interaction. Hence, the study presented in this article suggests that the staggered schemes are more robust in the added mass limit and significantly more efficient than their iterative counterparts. The strategy presented earlier by the authors 27 is a second-order accurate member of this family of staggered schemes and has been applied to a number of numerical examples and applications, [27][28][29][30] including hydraulic valves and a solar collector. Future work will include the numerical testing of the staggered scheme with combined relaxation and its formulation in a second-order accurate setting.