FSI by OOMPH-lib

Comput Mech DOI 10.1007/s00466-008-0270-6 ORIGINAL PAPER Solvers for large-displacement fluid–structure interaction problems: segregated versus monolithic approaches Matthias Heil · Andrew L. Hazel · Jonathan Boyle Received: 31 October 2007 / Accepted: 2 March 2008 © Springer-Verlag 2008 Abstract We compare the relative performance of monolithic and segregated (partitioned) solvers for large- displacement fluid–structure interaction (FSI) problems within the framework of oomph- lib, the object-oriented multi-physics finite-element library, available as open-source software at http://www.oomph-lib.org. Monolithic solvers are widely acknowledged to be more robust than their seg- regated counterparts, but are believed to be too expensive for use in large-scale problems. We demonstrate that mono- lithic solvers are competitive even for problems in which the fluid–solid coupling is weak and, hence, the segregated solvers converge within a moderate number of iterations. The efficient monolithic solution of large-scale FSI problems requires the development of preconditioners for the iterative solution of the linear systems that arise during the solution of the monolithically coupled fluid and solid equations by New- ton’s method. We demonstrate that recent improvements to oomph- lib’s FSI preconditioner result in mesh-independent convergence rates under uniform and non-uniform (adaptive) mesh refinement, and explore its performance in a number of two- and three-dimensional test problems involving the interaction of finite-Reynolds-number flows with shell and beam structures, as well as finite-thickness solids. Keywords Fluid-structure interaction · Monolithic solvers · Preconditioning 1 Introduction Numerical methods for the solution of large-displacement fluid–structure interaction (FSI) problems can be classified M. Heil (B) · A. L. Hazel · J. Boyle School of Mathematics, University of Manchester, Oxford Road, Manchester M13 9PL, UK e-mail: mheil@maths.manchester.ac.uk as either segregated (partitioned) or monolithic. In a mono- lithic approach, the complete system of nonlinear algebraic equations that arises from the coupled discretisation of the equations of motion in the fluid and solid domains is solved as a whole, typically using a variant of Newton’s method. This approach is generally acknowledged to be more robust than a segregated approach in which separate fluid and solid problems are coupled via a Picard (fixed point) iteration. It is widely believed, however, that monolithic solvers (i) are too computationally expensive and (ii) cannot take advantage of software modularity to the same extent as segregated solvers. The alleged superiority of the latter approach is generally attributed to the fact that in a segregated scheme “smaller and better conditioned subsystems are solved instead of one overall problem” [1]. It is also believed to be “difficult to devise efficient global preconditioners and to maintain state- of-the-art schemes in each solver” [2] when a monolithic solver is used. As a result, the monolithic approach is often regarded as unsuitable for large-scale problems. oomph- lib, the object-oriented multi-physics finite- element library, available as open-source software at http:// www.oomph-lib.org, was developed to address these con- cerns. A main design goal is to provide an environment that facilitates the monolithic discretisation and solution of multi- physics problems. The overall design of the library has been discussed in Ref. [3] and further details can be found in the online documentation. Here, we provide a brief overview of the library’s overall structure and an outline of the imple- mentation for FSI problems (Sect. 2). The main objective of the paper is to assess the performance of different solution strategies for large-displacement FSI problems. In Sect. 3 we compare the performance of monolithic and segregated solvers and demonstrate that the former are competitive even in cases where the fluid–solid coupling is relatively weak and, hence, segregated solvers converge in a moderate num- 123 Comput Mech ber of iterations. In Sect. 4 we discuss recent improvements to oomph- lib’s FSI preconditioner and demonstrate its excellent performance in a number of test problems solved using the library’s monolithic Newton solver. 2 Problem formulation and solution 2.1 oomph- lib’s overall design oomph- lib’s design is based on a (finite-)element-like framework in which the system of nonlinear algebraic equa- tions arising from the fully coupled discretisation of multi- physics problems is generated using an element-by-element assembly procedure. Typically, time-derivatives are treated implicitly and for the unsteady simulations described in this paper a second-order, backward-difference method (BDF2) is used for the unsteady terms in the fluid equations, and a Newmark method is used for those in the solid equations. The library provides a large number of single-physics ele- ments that are easily (re-)used in multi-physics problems, via a combination of inheritance and template (generic) pro- gramming. In addition, the library provides a variety of node- update strategies to accommodate moving boundaries. In all the test problems below, an algebraic node-update strategy was used to adjust the nodal positions within the fluid domain in response to changes in its (solid) boundaries. The key feature of this strategy is that the position of each node in the fluid domain depends only on the positions of a small number of solid nodes, resulting in fast mesh updates during the assembly of the fully coupled system. Furthermore, the approach leads to sparse shape-derivative matrices whose entries are computed automatically by finite-differencing. Further details can be found in Ref. [3] and on the oomph- lib webpages. Also present in the library are numerous high- level helper functions that facilitate the specification of multi- physics interactions; for instance, functions that automati- cally determine the (fluid mechanics) degrees of freedom that affect the fluid load on the solid. It is, therefore, straight- forward to combine two (or more) single-physics problems to create a monolithically coupled multi-physics problem. Once a monolithic discretisation has been specified, the coupled problem can be solved by oomph- lib’s Newton solver, acting on the fully coupled system of nonlinear alge- braic equations. A variety of direct and iterative linear solvers, together with appropriate preconditioners, are provided for the solution of the linear systems that arise during the Newton iteration. Alternatively, a segregated solution strategy may be employed: oomph- lib’s segregated FSI solver starts by “pin- ning” the degrees of freedom associated with the solid mechanics problem and modifies the assembly procedure to include only those elements associated with the fluids prob- lem. Thus, the Newton solver will solve the equations gov- erning the fluid motion for the current, “frozen”, wall shape. Next, the original boundary conditions for the solid problem are re-assigned, the fluid degrees of freedom are “pinned” and the assembly procedure is restricted to the solid ele- ments. After these changes, the Newton solver computes a new wall shape for the given flow field. If desired, differ- ent linear solvers/preconditioners can be employed for the fluid and solid solves, allowing the (re-)use of optimal solu- tion methodologies for the solution of the sub-problems. The basic fixed-point iteration can be augmented by constant or adaptive under-relaxation (the latter based on Irons & Tuck’s convergence acceleration technique for vector sequences [4]), and predictors for the wall displacement at the next timestep can be employed in time-dependent simulations. Numeri- cal experiments showed that for problems with strong FSI the convergence of the solid sub-solves within the segre- gated solution strategy were dramatically improved by updat- ing the fluid mesh (and hence the applied viscous stresses) after each linear solve. The algebraic node-update procedure allows mesh updates to be performed very quickly and so this step was performed by default as it had very little effect on the timings but significantly improved the robustness of the fixed-point iteration. Because oomph- lib’s monolithic and segregated solvers are implemented within the same overall framework, it is possible to perform a direct comparison between the two approaches. 2.2 oomph- lib’s fluid and solid elements and their interaction Within oomph- lib it is generally assumed that all lengths and coordinates have been non-dimensionalised on a problem-specific lengthscale L, while time is non- dimensionalised on some reference timescale T . Assuming that the fluid velocities are non-dimensionalised on a repre- sentative velocity U , the dimensionless Navier–Stokes equa- tions, which govern the flow of an incompressible Newtonian fluid with density ρ f and viscosity µ, are then given by1 Re ( St ∂ui ∂t + u j ∂ui ∂x j ) = − ∂p ∂xi + ∂ ∂x j ( ∂ui ∂x j + ∂u j ∂xi ) and ∂u j ∂x j = 0, (1) where Re=ρ f UL/µ is the Reynolds number, St =L/(UT ) is the Strouhal number, and the pressure is non- dimensionalised on the viscous scale µU/L. Equation (1), implemented in the Arbitrary–Lagrangian–Eulerian (ALE) 1 Throughout this paper, Latin indices take the values i = 1, 2, 3, Greek indices take the values α = 1, 2, and we use the summation convention that repeated indices are summed over all possible values of the index. 123 Comput Mech form, are the basis of all Navier–Stokes elements in the library. oomph- lib’s geometrically nonlinear Kirchhoff–Love beam and shell elements, used in the test cases in Sects. 3.1, 4.2.2 and 4.2.4, are based on the variational principle ∫∫ [ (σ γ δ 0 + Eαβγ δγαβ) δγγ δ + 1 12 ( h L )2 Eαβγ δκαβ δκγ δ ]√ a dξ1 dξ2 = ∫∫ [(L h ) √ A a f − Λ2 ∂ 2Rs ∂t2 ] ×δRs√a dξ1 dξ2, (2) which describes the large-displacements of an incrementally linearly elastic shell of thickness h and density ρs , Young’s modulus E and Poisson ratio ν. Here, Rs is the position vector to the shell’s deformed midplane, parameterised by the Lagrangian coordinates ξ1 and ξ2; a beam is understood to be the (obvious) restriction to the case when the mid- plane is parameterised by a single Lagrangian coordinate. The stresses, the fourth-order elasticity tensor Eαβγ δ , and the load vector f are non-dimensionalised by the structure’s effective Young’s modulus, Eef f = E/(1 − ν2). The (sec- ond Piola–Kirchhoff) pre-stress is represented by σαβ0 , and γαβ and καβ are the midplane strain and bending tensors.√ A dξ1dξ2 and √ a dξ1dξ2 represent infinitesimal area ele- ments of the shell’s deformed and undeformed midplanes, respectively. The parameter Λ=(L/T )√ρs/Eef f is the ratio of the structure’s natural timescale (for free in-plane vibra- tions) to the timescale T used in the non-dimensionalisation of the equations. oomph- lib’s general large-displacement elasticity ele- ments (available in displacement and pressure-displacement forms for compressible and (near-)incompressible behav- iour) are based on the variational principle ∫ σ i j δγi j dv= ∫ ( b−Λ2 ∂ 2Rs ∂t2 ) · δRs dv+ ∮ f · δRs d A, (3) where the two integrals are performed over the undeformed reference volume and over the deformed surface of the body, respectively. The second Piola–Kirchhoff stress tensor, σ i j , is determined as a function of the Green strain tensor γi j via a user-specified constitutive equation. σ i j is assumed to be non-dimensionalised on some characteristic stiffness para- meter, S, such as Young’s modulus E , which is also used for the consistent non-dimensionalisations of the body force b and the surface traction f ; and the timescale ratio is now given by Λ = (L/T )√ρs/S . The solid displacements affect the fluid via the induced changes in the domain geometry and via the no-slip condition u = St ∂Rs ∂t on fluid–solid interfaces. (4) The Cartesian components of the traction that the Newtonian fluid exerts onto the solid (on the solid stress scale) is given by f [FSI]i = Q ( −pδi j + ( ∂ui ∂x j + ∂u j ∂xi )) N j , (5) where the N j are the Cartesian components of the outer unit normal on the deformed solid (pointing into the fluid) and Q is the ratio of the stress scales used in the non- dimensionalisation of the solid and fluid equations. The para- meter Q indicates the strength of the FSI: as Q → 0, the fluid stresses acting on the structure become negligible, effectively decoupling the fluid and solid problems. We stress, however, that this statement is to be understood in an asymptotic sense. A “small” but finite value of Q does not necessarily imply that FSI can be neglected, particularly if the stiffness para- meter used to non-dimensionalise the solid stresses does not provide a good indication of the structure’s stiffness. For instance, a small value of Q in a thin-shell problem indicates that the fluid stresses are small relative to the shell’s exten- sional stiffness. If the shell deforms in a bending mode (in which it is much more flexible) the fluid stresses may still induce large deformations at small Q. 3 Comparing monolithic and segregated solvers 3.1 The test problem: flow in a collapsible channel We shall explore the relative performance of segregated and monolithic solvers using the well-studied FSI problem of flow in a collapsible channel. Incompressible, Newtonian fluid is driven through a 2D channel whose width is used as the lengthscale L. A Poiseuille (parabolic) velocity pro- file of mean velocity U is imposed at the inflow boundary and the outflow velocity is set to be parallel and axially traction free. A section of the upper channel wall is elastic and mod- elled as a pre-stressed Kirchhoff–Love beam loaded by an external pressure, pext, and the fluid traction (5). The system is known to develop large-displacement, self-excited oscil- lations (see, e.g., Ref. [5] and Fig. 2), provided that Re and Q are sufficiently large. We employed oomph- lib’s QTaylorHood (Q2 Q1) Navier–Stokes elements to discretise the ALE form of the unsteady Navier–Stokes equations and used the Kirchhoff– Love HermiteBeamElements to discretise the flexible part of the channel wall. Displacement control (prescribing 123 Comput Mech Fig. 1 a Load–displacement diagram: the vertical position, x [ctrl]2 , of a control point on the elastic wall (located at 50, 50, 60 and 70% of its length for Q = 0, 10−4, 10−3 and 10−2, respectively) as a function of the external pressure, pext . b, c Steady flow fields (streamlines and pressure contours in the vicinity of the elastic segment) at Re = 500, for b Q = 10−4 and c Q = 10−2 the vertical displacement of a control point near the antici- pated point of strongest collapse and solving for the unknown external pressure required to achieve that deformation) was employed to handle the system’s possible snap-through behaviour in steady simulations; see Fig. 1a. 3.2 Results of the simulations Simulations were performed for a channel in which a massless (Λ = 0) flexible wall segment of thickness h/L = 1/20 and length L = 5 was subjected to an axial prestress of σ0 = 103 and mounted on two rigid channels of length Lup = 1 and Ldown = 10. In Fig. 1a we illustrate the system’s load–displacement characteristics by plotting the vertical position of a control point on the elastic section of the wall, x [ctrl]2 , as a function of the non-dimensional external pressure, pext. For small values of the FSI parameter Q, the fluid traction is a small fraction of the load on the wall and an increase in pext leads to an approximately proportional increase in the wall deflection; the fluid reacts passively to the changes in the wall geome- try. For instance, in Fig. 1b the wall deformation is virtually symmetric, with the point of strongest collapse located at the centre of the elastic segment, indicating that the load on the wall is dominated by the spatially constant external pressure. As Q increases, the fluid traction noticeably affects Fig. 2 a Time-trace of the vertical position of the control point on the wall (located at 70% of its length) during a self-excited oscillation at Re = 500, St = 1 and Q = 10−2. b–d Snapshots of the unsteady flow fields (instantaneous streamlines and pressure contours) at b t = 32.4, b t = 34.9 and c t = 37.4 the load–displacement curve. At larger values of Q, a larger external pressure is required to keep the wall in its (approxi- mately) undeformed position because pext must balance the increasingly large fluid pressure in the elastic section, and so the load–displacement curves shift to the right. The increase in fluid pressure is a consequence of the fact that the outflow boundary conditions impose a zero fluid pressure at the end of the rigid downstream section, but the viscous pressure drop through that section increases with Q. Another consequence of the increased viscous pressure drop is that the point of strongest collapse moves downstream and so the control point is varied accordingly. In addition, the Bernoulli effect causes a reduction in fluid pressure in the most strongly collapsed part of the channel, locally increasing the compressive load on the wall. Hence, at larger values of Q a smaller increase in external pressure is required to collapse the channel to a given degree. At Q = 10−2 this destabilising effect is so strong that limit points develop in the load-displacement curve, indicat- ing that, at sufficiently large pext, the wall “snaps through” into a strongly deformed equilibrium configuration. The cor- responding plot of the wall shape in Fig. 1c illustrates that the wall deformation is now strongly affected by fluid pressure distribution. Figure 2a shows the time-trace of the wall displacement in an unsteady simulation at Re = 500, St = 1 and Q = 10−2. For this simulation the steady solution for pext = 1.68 was used as the initial condition for a time-dependent simulation in which pext was set to 2.51 at t = 0. The time-trace and 123 Comput Mech Fig. 3 Convergence histories of the segregated (dashed lines) and monolithic (solid lines) solvers for a steady computation at Re = 500 and various values of the interaction parameter Q. In a-c the convergence is characterised by the maximum residual of the coupled equations; in d convergence is assessed in terms of the maximum change in the solid variables between two successive iterations the snapshots of the flow field show that, following the initial perturbation, the system performs sustained, large-amplitude self-excited oscillations during which complex vortical flow structures develop downstream of the oscillating wall seg- ment. 3.3 The relative performance of the segregated and monolithic solvers Figures 3a–c show the convergence histories (maximum residual of the coupled equations vs. CPU time on a 3.60 GHz Intel Xeon processor with 2 GB of memory) of the mono- lithic and segregated solvers for a steady computation at Re = 500, and for three different values of the interaction parameter Q. The six curves in each figure represent the convergence histories during six nonlinear solves performed in the course of a parameter study during which the chan- nel’s maximum collapse was increased from 0 to 35% of its width. For simplicity, all linear systems were solved with oomph- lib’s default linear solver, SuperLU [6], a sparse direct solver which is efficient for this moderately sized 2D problem (20,390 unknowns); see Tables 1 and 3. For all cases considered, both nonlinear solvers converge in a moderate number of iterations, with roughly comparable Table 1 Average CPU times (in seconds) for the monolithic Newton solver applied to the steady collapsible channel problem with Re = 500 and Q = 10−2, for different linear solvers/preconditioners NDOF SuperLU GMRES & PPP