Three-field partitioned analysis of fluid–structure interaction problems with a consistent interface model

Available online at www.sciencedirect.com ScienceDirect Comput. Methods Appl. Mech. Engrg. 414 (2023) 116134 www.elsevier.com/locate/cma Three-field partitioned analysis of fluid–structure interaction problems with a consistent interface model José A. Gonzáleza,∗, K.C. Parkb aEscuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, Sevilla 41092, Spain bAnn and H. J. Smead Aerospace Engineering Sciences, University of Colorado, Boulder, CO 80309-429, USA Received 8 March 2023; received in revised form 16 May 2023; accepted 16 May 2023 Available online 6 June 2023 Abstract This paper proposes a new Fluid–Structure Interaction computational framework. The coupling between the solid and an incompressible fluid is formulated by means of the method of localized Lagrange multipliers (LLM). Instead of applying a direct coupling between the fluid and the structure, which is the traditional approach, LLM introduces an intermediate surface with its own degrees of freedom that is connected to the fluid and structure sides using independent fields of localized Lagrange multipliers. This approach allows the connection of non-matching meshes with mortar or classical localized methods and provides consistent dynamic equations of motion for the interface that can be integrated in parallel. Interface multipliers are later eliminated and the interface motion is used to update the fluid and structure states. This way, dedicated stand-alone software modules for the fluid and the structure are connected to a third interface system treating their interaction, thus preserving the modularity of the single-discipline software modules. Different numerical examples are solved with the proposed methodology to prove its efficiency and accuracy by running a series of classical dynamic FSI benchmark problems. © 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Keywords: FSI; Partitioned analysis; Coupled systems; Localized Lagrange multipliers 1. Introduction The simulation of fluid–structure interaction (FSI) problems with finite elements is today a mature numerical technique with many applications in science and engineering. A combined Lagrangian description of motion for the structure with an arbitrary Lagrangian–Eulerian (ALE) form for the incompressible fluid is the preferred approach when the interface motions are not excessive and many different solution schemes are readily available for its solution. However, there is still room for improvement in some aspects related with the accuracy, stability and computational efficiency of these methods. In particular, we believe that a broader view of the partitioning concept can still bring important advantages to the formulation and implementation of more efficient FSI algorithms. Historically, staggered partitioned solution procedures represent the first and most elementary computational strategy for the simulation of fluid–structure interaction problems. These methods usually involve the use of a predictor followed by the separate execution of each subdomain solver maintaining stability and accuracy. The idea ∗Corresponding author. E-mail addresses: [email protected] (J.A. Gonz´ alez), [email protected] (K.C. Park). https://doi.org/10.1016/j.cma.2023.116134 0045-7825/© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http: //creativecommons.org/licenses/by-nc-nd/4.0/). J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 goes back to the early works of Park and Felippa in the field of elasto-acoustic FSI problems [1–3] where the first conditionally stable staggered solution procedures were proposed. Since then, much work has been dedicated to the improvement of this concept [4,5]. For a general review of numerical methods for coupled FSI problems, see [6–8]. A coupled problem generally consists of two or more physical systems interacting at common boundaries by exchanging information affecting their internal state. The physical equations governing the different systems are usually discretized by using nonconformal space discretizations and often approximated with different time integration schemes for the different domains. Three basic approaches may be distinguished for the solution of multi-physics coupled problems: loosely-coupled schemes,partitioned (staggered) schemes and monolithic schemes. Particularly in FSI applications, loosely-coupled schemes solve the problem sequentially in an uncoupled fashion. First the increment in the fluid solution is found on a fixed spatial domain, then the fluid forces on the structure are collected, and finally the structural-solution increment is computed with these forces, which is followed by an update of the mesh position. This enables the use of existing fluid and structural solvers, a significant motivation for adopting this approach. If the fluid dynamical problem and the structural problem are solved in a loosely coupled manner, without the feedback of subiterations, the coupling conditions at the interface will not be exactly satisfied at each time step. This solution process works well particularly in aeroelastic applications or in applications involving compressible viscous fluids where added mass effects are small [9,10]. For cases of moderate to high fluid–structure density ratios, typically when the mass of the fluid subsystem is greater than the mass of the structure, the addedmass effect of the surrounding fluid on the structure plays a dominant role and loose coupling schemes tend to be unstable, especially in the context of incompressible flows [11]. In this case, the methods are conditionally stable, inaccurate and very often the time-step limitations are so severe that the use of partitioned or monolithic schemes is preferred. In monolithic FSI schemes, the strongly-coupled equations of the fluid, the solid, and the mesh deformation are solved simultaneously. The system of equations is obtained for each Newton iteration and solved directly or iteratively. Monolithic schemes in FSI have been employed in isogeometrical analysis IGA [12,13] with the effect of the structural and the mesh motion on the fluid equations included on the left-hand-side matrix for robustness and the system solved with preconditioned GMRES methods. Monolithic approaches have been applied on multiple FSI methods, like multigrid solvers [14,15], single field predictors [16], unified field formulations [17], Lagrangian meshfree methods [18], hybrid Eulerian-ALE formulations [19] and applications in free-surface contact [20] or coupled thermo-FSI problems [21]. Their main advantage is that monolithic FSI solvers are more robust. Many of the convergence problems encountered with the loosely coupled approaches are completely avoided. However, the monolithic approach necessitates writing a dedicated coupled fluid–structure solver, thus precluding the reuse of existing fluid and structure solvers. A compromise between the loosely and monolithic schemes is given by the partitioned FSI methods. Partitioned coupling schemes try to compute the solution of the monolithic system by using an iterative process, involving the repeated solution of the fluid and structure partitioned systems. The strength of this approach originates from its modularity and the capability of using for each physics the most effective numerical method. Mathematically, they are equivalent to an iterative solution of coupled equations by blocks using algebraic methods like Jacobi, Gauss– Seidel or block Newton methods. Three different types of partitioned strategies exist. First, fixed-point iteration methods [22] that algebraically resemble the block Gauss–Seidel or multiplicative Schwarz iteration for the solution of systems of equations [23]. Block iterative Gauss–Seidel solves one system first, then transfers the interface newly computed values to the other system, solves the second system, and this is repeated iteratively until convergence. This process may not converge, even for starting values close to the solution. The stability of Block Gauss–Seidel partitioned solution of the coupled problem was studied by [24] exchanging interface displacements and forces until convergence. Among these methods are the fixed-point iterations with fixed or dynamic under-relaxation [25,26]. Block Gauss–Seidel partitioned solution methods have also been successfully employed for coupling finite element and finite volume solvers [27,28]. Sencond, vector extrapolation methods constitute a generalization of the fixedpoint schemes that incorporate many history values in one relaxation step. This interpolation of the improves the convergence [29]. Finally, Block-Newton methods solve the linearized equations by blocks using independent solvers. For example, Block-Gauss elimination has been applied on the monolithic equations [30], four steps of solutions exchanges are required per iteration and quadratic convergence is attained. Partitioned schemes represent the predominant approach for FSI solvers with mixed Lagrangian-ALE formulations. In their most basic form, involve a close loop between the fluid, the structure and the dynamic mesh, see 2 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 1. Representation of a classical fixed-point iteration approach for the three-field coupled FSI problem. In a sequence of Dirichlet– Neumann boundary solutions exchange, the Structure is solved first with the boundary forces λFexerted by the fluid to find the structural displacements uS. Then, the Dynamic Mesh of the fluid is informed of the new structure boundary displacements BSTuSand its position is updated to uM. Finally, the Fluid is solved using the mesh velocities ˙ uMand its variables, velocity vFand pressure pF, are recalculated for the next iteration. Fig. 1. For example in a fixed-point iteration method, inside a Dirichlet–Neumann iteration, the ALE fluid is defined as the Dirichlet system, driven by the structural interface motion, while the structure system reacts to the Neumann conditions applied by the fluid on the interface. Clearly, this scheme is more flexible than a monolithic approach. It allows for the use of independent software modules, it is possible to employ different numerical methods for each system, implement locally different time-integrators, adopt asynchronous time-integration. Another important aspect of FSI modeling is to define how the state of one physical system is transferred across the fluid–structure interface to the other. Note that normally, due to different space resolution requirements or simply due to the existence of complex interface geometries between the solid and the structure, we have to deal with nonmatching meshes at the interface and there is not a one-to-one mapping between boundary variables. For this task, there exists two basic groups of coupling techniques [31]: •Direct coupling methods: Finding direct operators relating the fluid and structure interface motions. The classical techniques can be categorized into: point-to-point mapping based on interpolation with scattered data points using for example radial basis functions, point-to-element projection where kinematical constraints are enforced at the fluid nodes or quadrature points and common-refinement projection introducing a common refinement surface for an accurate load transferring between fluid and structure. •Lagrange multiplier methods: The interface constraints are enforced in a variationally consistent way by using Lagrange multipliers to couple non-matching fluid and structure meshes. The most common method of this group is the mortar method, which was originally proposed in the context of non-overlapping domain decomposition and is referred in the literature as a coupling method with desirable mathematical and numerical properties. Another member of this group is the method of localized Lagrange multipliers that uses independent fields of Lagrange multipliers to connect the coupled systems with the shared interface. The integration of the dual mortar method within a general ALE-FSI framework was first proposed in [16,32] where a constrained coupled four-field system with a saddle point type structure is obtained, that is later reduced to a classical three-field system by condensation of the Lagrange multiplier unknowns. The method adopts a monolithic solution approach with Newton–Raphson iterations and attains optimal temporal convergence. Mortar method provides higher accuracy but also requires, for the construction of the interface mapping operator, an elaborate process of projection of slave and master nodes onto the common interface, perform polygon clipping, divide clip polygon into triangular integration cells and perform Gaussian integration. On the other hand, the LLM method introduces a dedicated interface system endowed with its particular displacement field and couples separately the fluid and the structure to the interface with independent fields of Lagrange multipliers. Specifically in the field of FSI, the method of localized Lagrange multipliers has been applied to study the interaction of flexible structures with acoustic fluids [33,34] and incompressible fluids adopting a pure Lagrangian description of motion [35,36]. The extension of the method to FSI problems with incompressible fluids described in ALE form is treated in this paper. The potential advantages of separately coupling fluid and structure to an independent interface with its own discretization were also highlighted by Dettmer et al. [37–39] noting that, this way, the interface degrees-of-freedom 3 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 2. Simple one-dimensional FSI partitioned linear model coupled with classical Lagrange multipliers (left) or using localized Lagrange multipliers (right). Introduction of localized Lagrange multipliers allows to solve for the interface accelerations that are used to update the location of fluid and structure boundaries. (DOFs) can then be connected to the fluid and the solid boundaries using independent interpolations. This allows a more modular computer implementation and thus facilitates the maintenance and extension of the fluid and solid solvers. It is also possible to employ connection techniques like the method of Lagrange multipliers for the data transfer between the coupled systems. This work proposes a new computational framework based on the use of LLM for the solution of FSI dynamic coupled problems. Instead of applying a direct coupling between the fluid and structure meshes, which is the traditional approach, the LLM introduces an intermediate system at the interface with its own degrees of freedom that is connected to the fluid and structure sides using independent fields of localized Lagrange multipliers. This approach brings a more modular computer implementation, allows the connection of non-matching meshes using different methods and produces consistent dynamic equations of motion for the interface that can be efficiently solved in parallel. The paper is organized as follows. The following Section 2introduces the idea of LLM in FSI illustrating with a simple example the relevance of the dynamic interface equations of motion in a partitioned solution scheme of FSI problems. Section 3defines the governing equations of the four systems involved, dynamic mesh, fluid, structure and interface. Section 4derives the variational formulation of the coupled problem. The Galerkin finite element approximation of these variational forms is performed in Section 5including the stabilization of the fluid equations to allow equal order velocity–pressure interpolations. The interface model presents an independent discretization that needs to be designed in order for maximum accuracy. Its construction process is described in Section 7. Time integration of the obtained nonlinear coupled equations using the generalized α-method is described in Section 8. A monolithic version of the solution algorithm applying Newton method and subiterations is described in Section 9. A culminating stage of the paper is the development of the partitioned algorithm presented in Section 10. We believe that it is simple to understand and easy to implement. In essence, utilizes a time integration of the dynamic interface equations for enhancing stability and parallelism of the partitioned scheme. In Section 11, a series of numerical examples are compared with analytical solutions to demonstrate the accuracy and robustness of the proposed methodology. Finally, Section 12 presents some concluding remarks. 2. Localized Lagrange multipliers in FSI The method of Lagrange multipliers first entered into the FSI arena as a tool for improving the accuracy at the interface between non-matching meshes. In practical applications, the multipliers are eliminated and the coupled system solved monolithically. In this Section it is demonstrated how Lagrange multipliers through LLM also provides a theoretically consistent dynamic equilibrium equation for the interface that can be treated as a truly separated system under the influence of fluid and structure interactions. The stability of FSI solution methods with Lagrange multipliers can be analyzed by using a simple onedimensional linear example proposed in [40,41]. The model is represented in Fig. 2 and consists of two linear one-DOF models of the fluid and the structure, in our case connected by classical Lagrange multipliers (left) and localized Lagrange multipliers (right). The equations of the FSI model treated with classical Lagrange multipliers 4 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 3. Convergence of the interface force for a partitioned 1D FSI model with ω=1, ξ=0.01 solved using classical Gauss–Seidel with relaxation (GS-θ), classical (LM) and localized (LMM) Lagrange multipliers. It is shown the iterative solution for the first time step, obtained using Newmark time integration method with average constant acceleration (γ=1 2, β =1 4)and time step ∆t=0.01. Convergence of the interface forces for structure masses α=0.25 (left) and α=0.75 (right). Gauss–Seidel without relaxation (θ=1) is unstable for α < 0.75 while implicit Lagrange multiplier methods are always stable solving the problem in only one iteration. λ, see Fig. 2(left), are expressed: ⎧ ⎪ ⎨ ⎪ ⎩ rF(vF, λ)=(1 −α)˙vF+2ξωvF−fF−λ=0 rS(uS, λ)=α¨uS+ω2uS−fS+λ=0 rI(vF,uS)= −˙vF+ ¨uS=0 (1) where ω=√k/mis the natural frequency of the coupled system and ξ=c/(2 mω) is the damping ratio associated to the fluid dissipation. Also the mass parameter 0 ≤α≤1 where 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. Next, the stability of a staggered scheme is investigated. The Gauss–Seidel with under-relaxation solution process involves iterating on the interface multipliers in (1) as follows: knowing λi, first the structure equation rS(uS i+1, λi) is solved for the increment of structural displacements ∆uS i+1and then, using the interface condition rI(vF i+1,uS i+1) a fraction θof this structural displacement is applied to the fluid interface, i.e., ∆vF i+1=θ∆uS i+1. Then, the fluid equation rF(vF i+1, λi+1) can be solved for the updated interface forces λi+1and the process is repeated until the convergence criteria |λi+1−λi|< ε is satisfied. The result of this procedure is shown in Fig. 3 where it is represented the convergence of the interface force for increasing structural masses, α=0.25 (left) and α=0.75 (right), combined with different levels of underrelaxation. It is observed that the unrelaxed solution (θ=1) is only stable for large values of the structural mass (α > 0.75), hence some relaxation is needed and significant structural deformations can be taken into account by an under-relaxation of the boundary geometry. We can see that small values of mass-ratio parameter αrequire lower values of the under-relaxation factor θto achieve unconditional stability [40]. Moreover, the required level of under-relaxation for stability is not constant and depends on the fluid–structure mass ratio. Hence, this solution strategy is parameter dependent and also difficult to parallelize. On the other hand, a monolithic solution permits to eliminate the Lagrange multipliers and reduce the structure interface displacements finally needing to solve a full system with the rest of unknowns [32]. The result of this solution method, indicated as LM in Fig. 3, is always stable and requires a single iteration for convergence in a linear case. However, this strategy is not commonly employed in practical applications because it is computationally very expensive and difficult to parallelize. Also, iterating on the Lagrange multipliers would not be very effective when combined with an ALE description of the fluid motion, because direct updating of the fluid background mesh requires first the determination of the interface displacements. These difficulties can be avoided by using the method of localized Lagrange multiplies (LLM) exemplified in Fig. 2(right). The method introduces an explicit definition of the interface with displacement uIand now fluid and 5 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 structure are connected to this interface with independent Lagrange multipliers (λF, λS). The semidiscrete equations of motion in this case become: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ rF(vF, λF)=(1 −α)˙vF+2ξωvF−fF+λF=0 rS(uS, λS)=α¨uS+ω2uS−fS+λS=0 rF I(vF,uI)= ˙vF− ¨uI=0 rS I(uS,uI)= ¨uS− ¨uI=0 rI λ(λF, λS)= −λF−λS=0 (2) with three coupling equations, instead of one, corresponding to the two-side kinematic constraints plus the interface equilibrium condition, given by the last equation rI λ(λF, λS)=0. Suppose now that these LLM equations are being solved for a particular time-step by adopting the Newton– Raphson iteration scheme. In this case, to solve for a new subiteration (i+1), the incremental form of equations is written: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (1 −α)(∆˙vF i+1+ ˙vF i)+2ξωvF i−fF+λF i+1=0 α(∆¨uS i+1+ ¨uS i)+ω2uS i−fS+λS i+1=0 ∆˙vF i+1−∆¨uI i+1=0 ∆¨uS i+1−∆¨uI i+1=0 −λF i+1−λS i+1=0 (3) a linear system to be solved for the acceleration increments and the new multipliers. However, after isolating the fluid and structure acceleration increments from the two first equations, substituting into the third and fourth boundaryinterface acceleration constraints and using the last interface equilibrium condition, the following interface update equation is found: ∆¨uI i+1= −[(1−α)˙vF i+2ξωvF i−fF   Fluid residual ]−[α¨uS i+ω2uS i−fS   Structure residual ] (4) where the right hand side represents the addition of fluid and structure boundary residuals at the previous iteration, which is equal to the preceding unbalance of interface forces. These fluid and structure reaction forces can be computed independently and in parallel from each physics, then the interface position updated and finally the fluid and structure systems recalculated for the next iteration. It is noted that, when the problem is linear and the iteration solution process made implicit, the described method converges in only one iteration as shown in Fig. 3 for Lagrange multiplier methods. Therefore, the basic idea of the proposed LLM method is simple, localized Lagrange multipliers can be eliminated while the independent residuals of the fluid and the structure are used to update the position of the interface. After this update, both systems can be actualized independently in parallel using the new position of the interface. In the following, we will describe how this technique of LLMs can be applied to complex FSI interaction problems presenting non-matching meshes. 3. Governing equations In this section we introduce the governing equations of the FSI problem, which consists of a structural domain ΩS(t) interacting with a deformable fluid domain field ΩF(t) in Rd, where d∈{2,3}is the number of space dimensions and t∈[0,T]is the considered time interval, see Fig. 4. Both fields interact through a fluid–structure common interface Γ=ΓF∩ΓS, where equilibrium and kinematic compatibility conditions need to be fulfilled. The structure uses a Lagrangian description of motion and the fluid is described in Arbitrary Lagrangian Eulerian (ALE) form. In the following, the fluid quantities are denoted by the superscript (·)F, the structure by (·)Sand the ALE fields by (·)M. Classical FSI methods are formulated in terms of interaction between these three systems. Here we introduce an additional explicit representation of the interface, denoted as Γ, which follows the fluid and structure interfaces and reformulate the FSI system in terms of interaction of the previous three systems with the interface. An extra superscript (·)Iis introduced to indicate that the quantity is located at this FSI interface Γ. 6 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 4. Representation of a typical FSI scenario where a new interface Γ=ΓF∩ΓSis defined between the two systems and connected using localized Lagrange multiplier fields (λF,λS)to the fluid and the structure. 3.1. Fluid equations in ALE form The fluid is considered incompressible and governed by the Navier–Stokes equations in a deformable domain ΩF(t) kinematically described in ALE form by through mapping x:ΩM 0→ΩM(t): x(X,t)=X+uM(X,t),∀X∈ΩM 0,t≥0 (5) where X∈ΩM 0is a location in the reference configuration that through the ALE mapping locates the point x∈ΩF in the actual fluid domain configuration and uM(X,t)is the displacement of the reference configuration at time t∈[t0,T]. Then, the velocity of points in the reference configuration is given as: vM=∂x(X,t) ∂t⏐⏐⏐X(6) where (·)|Xindicates ALE time derivative. The incompressible Navier–Stokes equations in ALE form are expressed in terms of the fluid velocity vF(x,t) and the fluid pressure pF, satisfying conservation of momentum and the incompressibility condition in the domain of interest together with the boundary conditions: ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ρF{∂vF ∂t⏐⏐⏐X+(vF−vM)·∇vF−bF}−∇·σF(vF,pF)=0in ΩF(t) ∇·vF=0 in ΩF(t) vF=¯ vFon ΓF v σF·nF=¯ tFon ΓF t (7) where ρFis the constant fluid density, bFis a given body force per unit mass and (vF−vM)is the ALE convective velocity. The second order tensor σFis the Cauchy stress tensor, which for a Newtonian fluid is given by: σF(vF,pF)= −pFI+µFε(vF)(8) with µFthe dynamic viscosity of the fluid, Ithe identity tensor and ε(vF) the symmetric strain-rate tensor calculated as: ε(vF)=1 2(∇vF+∇vFT)(9) a symmetric second order tensor function of the fluid velocity gradients. Finally, the fluid velocity ¯ vFon Dirichlet boundaries ΓF vand the fluid traction ¯ tFon Neumann boundaries ΓF t need to be prescribed, together with the initial conditions: vF=vF 0in ΩF 0(10) where the initial velocity field vF 0should also fulfill the incompressibility condition. 7 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 3.2. Structure equations The structure is assumed to present a nonlinear elastic behavior and its equations motion given by the equations of elastodynamics expressed in strong form as: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ρS(∂2uS ∂t2−bS)−∇·PS(uS)=0in ΩS(t) uS=¯ uSon ΓS u σS·nS=¯ tSon ΓS t (11) where uS(x,t)is the displacement vector field, ρSthe structural density, bSare the applied body forces per unit mass and PSis the first Piola–Kirchhoff stress tensor with σS=JS−1PSFSTthe Cauchy stress tensor, FSthe deformation gradient tensor and JS=det FSthe Jacobian of the transformation. The interface tractions ¯ tSare acting on the Neumann boundaries ΓS twith normal vector nSand ¯ uSdefined on ΓS urepresent the Dirichlet boundary conditions. The material of the structure is assumed to exhibit a compressible Neo-Hookean behavior. This material presents characteristics that can be identified with the familiar material parameters found in linear elastic analysis. The strain energy density function and stress tensor of such a material is defined as: ΨS(CS)=µS 2(CS:I−d)−µSln JS+λS 2(ln JS)2,PS=∂ΨS ∂FS(12) where λSand µSare the material Lam´ e parameters and CS=FSTFSis the right Cauchy–Green tensor. Note that in the absence of deformation, that is, when CS=Iand JS=1, the stored energy function vanishes as expected. Furthermore, adequate initial conditions with prescribed displacement and velocity fields, uS 0and vS 0, are defined as: {uS=uS 0in ΩS 0 vS=vS 0in ΩS 0 (13) for the complete structure at the initial time t0. 3.3. Interface equations Instead of performing a direct coupling between the fluid and the structure as in common FSI formulations, we include an additional surface at the fluid–structure interface Γ(t) and describe its motion in time introducing an independent displacement field uI(x,t). Then, see Fig. 4(right), fluid and structure are connected to this explicit representation of the interface by using two independent fields of localized Lagrange multipliers, named λFon the fluid side and λSon the structure side. These two fields of Lagrange multipliers are identified with the interface tractions on the fluid tFand the structure tS, appearing in Eqs. (7) and (11) as boundary tractions acting on Γ(t)=ΓF(t)∩ΓS(t), respectively. Finally, the coupling conditions on the fluid–structure interface are expressed in strong form: ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ vF−∂uI ∂t=0on Γ(t) uS−uI=0on Γ(t) λF+λS=0on Γ(t) (14) where the two first equations enforce the kinematical compatibility constraints of the fluid and structure with the interface, while the last equation represents the equilibrium condition of the interface. 4. Weak formulation In this section we derive the weak formulation of the different systems involved; the incompressible Navier– Stokes equations expressed in ALE form for the fluid, the dynamic equilibrium equation of a nonlinear elastic structure and the interface constraints coupling both systems. 8 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 4.1. Weak form of fluid-mesh motion The deformation of the ALE background domain ΩMis produced by the motion of the fluid–structure interface Γand assumed to satisfy the equations of linear elastostatics. In this case, the elastic material properties are purely artificial parameters that are selected in order to minimize the fluid mesh distortion during the deformation process. It is also needed to include a compatibility condition on Γbetween the fluid mesh and the interface motion that is enforced using an additional field of Lagrange multipliers. To formulate this part of the fluid problem, adequate solution spaces UMand LMare defined respectively for the mesh displacement fields uMand Lagrange multipliers λMwith corresponding test function spaces WMand WM I. We can then derive the discrete weak form of the fluid-mesh motion, which consists in finding (uM,λM)∈UM×LM such that ∀(wM,wM I)∈WM×WM I, it is fulfilled: ∫ΩM ε(wM):DM:ε(uM)dΩ+∫Γ wM·λMdΓ+∫Γ wM I·(uM−uI)dΓ=0 (15) where DMis the elastic fourth order constitutive tensor with material parameters (λM, µM)selected to minimize the fluid mesh distortion. Appropriated spaces for the interface displacements uIare introduced later in the definition of the interface problem. 4.2. Weak form of the fluid problem Let us now introduce admissible spaces VFand PFof trial functions on ΩFfor the fluid velocities vFand pressure pF, with their corresponding spaces of test functions WF mand WF c, respectively. The incompressible Navier–Stokes balance equations expressed in ALE form (7) are multiplied by corresponding test functions (wF, wF c)from their respective spaces and integrated on the current configuration. Additionally, partial integration of the pressure and viscous stresses is performed in the resulting balance of linear momentum. This reduces the occurring order of spatial derivatives to one and leads to additional boundary integrals in the weak form of the fluid problem. After these steps, the weak form of the fluid problem is derived, which consists in finding the fields (vF,pF,λF)∈ VF×PF×LFsuch that: ∫ΩF ρFwF·{∂vF ∂t⏐⏐⏐X+(vF−˙ uM)·∇vF−bF}dΩ+∫ΩF ε(wF):σFdΩ− ∫ΓF t wF·tFdΓ+∫Γ wF·λFdΓ+∫ΩF wF c∇·vFdΩ=0 (16) is fulfilled ∀(wF, wF c)∈WF m×WF c, where the unknown boundary tractions λF= −tFcan be interpreted as a Lagrange multiplier vector appearing to compatibilize the fluid interface motion with the structure. The complete approximation of these interface Lagrange multiplier fields will be specified later in combination with the interface problem. 4.3. Weak form of the structure problem For the structure, admissible spaces USof trial functions and WSof weighting functions are defined on ΩSfor the displacement field uSand its corresponding test functions wS. Then, the nonlinear structural equations (11) multiplied by arbitrary test function wS∈WSand integrated in the domain ΩSwith subsequent partial integration leads to the weak form that is expressed in the following terms: find (uS,λS)∈US×LSsuch that: ∫ΩS ρSwS·{∂2uS ∂t2−bS}dΩ+∫ΩS ∇wS:PSdΩ−∫ΓS t wS·tSdΓ+∫Γ wS·λSdΓ=0 (17) is fulfilled ∀wS∈WS. Here, the interface tractions on the structural side λS= −tSarise from the needed compatibility condition between the structure and the fluid interface. These multipliers are chosen from a corresponding solution space defined in the following subsection. 9 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 6. Classical LLM interface mesh construction from a defined position of the fluid boundary nodes xFand structure boundary nodes xS. The interface nodes are located at xI, the zero-moment roots of the bending moment diagram M(x) produced by the multiplier forces acting on the interface. where M0represents the roots of the bending moment function for certain position of the fluid and the structure. There exists a number of different works in the literature treating the problem of constructing LLM interfaces for more complex cases. For example, the coupling of different numerical methods [51], treating pathologies of nonmatching meshes with gap elements [52–54] and automatic methods for the construction of LLM interfaces for nonmatching 3D meshes [55]. 8. Time integration Generalized α-method is a general time-integrator with second-order accuracy and very small dissipation in the low frequency range, initially proposed by Chung and Hulbert [56] as an extension of the HHT α-method [57] for the time integration of the second-order equations of structural dynamics. The method was later extended to solve the first-order Navier–Stokes equations by Jansen [58] and analyzed for FSI problems in [59]. To describe its incorporation into the formulation, let us first divide the time interval of interest [0,T] into a sequence of subintervals (tn,tn+1) with fixed time-step size ∆t=tn+1−tn. The classical Newmark formulas are then applied to approximate velocities and displacements at next time step (n+1) using information from the previous time step (n) in the following form, for the fluid: vF n+1=vF n+∆t{(1 −γ)˙ vF n+γ˙ vF n+1}=vF p+δv˙ vF n+1(52) and equivalently for the structure: {˙ uS n+1=˙ uS n+∆t{(1 −γ)¨ uS n+γ¨ uS n+1}=˙ uS p+δv¨ uS n+1 uS n+1=uS n+∆t˙ uS n+∆t2{(1 2−β)¨ uS n+β¨ uS n+1}=uS p+δu¨ uS n+1 (53) where (γ, β)∈[0,1] are the Newmark time integration parameters; δv=γ∆tand δu=β∆t2are contribution factors of the next-time acceleration to the velocity and displacement vectors; and the velocity predictions ˙ uF pand ˙ uS p, together with the structural displacement prediction uS pthat are introduced to simplify upcoming expressions. Next, the different terms of the semi-discrete equations are evaluated by separately averaging, between both time instants (tn,tn+1), the inertia forces and the internal/external forces. To do that, two time averaging parameters (αm, α f)∈[0,1] are introduced to define separately the time instants for the evaluation of inertia forces and internal/external forces respectively. This average is defined for the fluid: {˙ vF n+αm=˙ vF n+αm(˙ vF n+1−˙ vF n) vF n+αf=vF n+αf(vF n+1−vF n)(54) and similarly for the structure: {¨ uS n+αm=¨ uS n+αm(¨ uS n+1−¨ uS n) uS n+αf=uS n+αf(uS n+1−uS n)(55) 16 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 where inertia terms are computed at instant (n+αm) and internal/external terms at instant (n+αf). The specific value of the time averaging parameters (αm, α f) is selected in order to maximize the accuracy and stability of the algorithm. Note also that different time integration methods can be reproduced by a proper selection of these parameters. In particular, αm=αf=0 corresponds to the Newmark method, αm=αf=1 2to the Midpoint algorithm and αm=1 with αf=αto the HHT α-method. A final important aspect to discuss is the selection of the time-integration parameters. The optimal value of these parameters is not the same for structural dynamics than for fluid dynamics and here we give priority to the stability of the fluid solution. It is known that for the first-order linearized Navier–Stokes equations [58], the generalized α-method provides unconditional stability if the time averaging parameters satisfy the condition: 1 2≤αf≤αm(56) while second-order accuracy is obtained with an specific choice of the time integration parameters: γ=1 2+αm−αf, β =1 4(1 +αm−αf)2(57) introducing into the solution a small numerical damping of intensity controlled by the values of the two averaging parameters αmand αf. These time averaging parameters can be selected from a desired spectral radius at infinity condition by using the expressions: αm=1 2(3−ρ∞ 1+ρ∞), α f=1 1+ρ∞ (58) where ρ∞∈[0,1] is the spectral radius at infinity and its value can be used as a single parameter to control the high-frequency dissipation with the guarantee that second-order accuracy and unconditional stability are maintained. 9. Monolithic quasi-direct FSI solver with LLMs The finite element formulations described in Section 5lead to a group of coupled nonlinear equations that need to be solved for each time step of the simulation. In particular, note that the discrete ALE fluid momentum residual equation in (28) presents an implicit dependence on the mesh velocity ˙ uM. This means that a full linearization of the fluid residual would need to include the derivatives with respect to the mesh deformation, i.e., the so-called shape derivatives. However, it is a common practice to disengage the mesh motion from the coupled FSI problem by assuming a small influence of the fluid-mesh variation. This approach is known as quasi-direct coupling [60–62]. By adopting a quasi-direct approach, the formulation of the coupled FSI problem starts by collecting the individual residual contributions coming from the fluid conserving equations (28), the structure equations (41) and the interface equations (48), in a semi-discrete residual with the form: r(z)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ rF m(˙ vF,vF,pF,λF) rF c(˙ vF,vF,pF) rS(¨ uS,uS,λS) rF I(vF,˙ uI) rS I(uS,uI) rI(λF,λS) ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˆ rF m(˙ vF,vF,pF)+BFλF rF c(˙ vF,vF,pF) ˆ rS(¨ uS,uS)+BSλS BFTvF−LF˙ uI BSTuS−LSuI −LFTλF−LSTλS ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ =0(59) where z=(˙ vF,pF,¨ uS,λF,λS,¨ uI)is the global vector with the problem unknowns. Next step is to discretize in the time the residual (59) by using the generalized α-method described in Section 8. This way, the time collocation equations are given by: r(zn+1)= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˆ rF m(˙ vF n+αm,vF n+αf,pF n+1)+BFλF n+1 rF c(˙ vF n+αm,vF n+αf,pF n+1) ˆ rS(¨ uS n+αm,uS n+αf)+BSλS n+1 BFTvF n+1−LF˙ uI n+1 BSTuS n+1−LSuI n+1 −LFTλF n+1−LSTλS n+1 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ =0(60) 17 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 satisfying dynamic equilibrium and compatibility of the three systems at the interface. Inertia forces are evaluated at intermediate step (n+αm), internal/external forces at a different middle step (n+αf), while the constraints are fulfilled at the end of the time-step (n+1). The system of simultaneous equations (60) is nonlinear, in general, and can be solved using Newton–Raphson method. Let subscript (i) now denote the Newton method iteration step and zithe solution estimate in the ith iteration step. At the beginning of the iteration loop, a first guess z0is required, which can usually be chosen as the solution at previous time step tn, i.e., z0=zn. Then, a linearization of the residual is used to find the next solution estimate zi+1=zi+∆ziby solving the linear system: J(zi)∆zi= −r(zi)(61) where the full Jacobian matrix defined as: J(z)=∂r ∂z(62) and the iteration process is repeated until the norm of the residual falls below a certain epsilon ∥r(zi+1)∥< ε. In our case, the Jacobian of system (61) is constructed by collecting contributions from the partial derivatives of the different residuals; corresponding to the fluid momentum conservation (31), the fluid incompressibility constraint (36), the structure dynamic equilibrium (43) and the linear interface constraints (48), ordered in the linear system: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ JF mm JF mc 0 BF0 0 JF cm JF cc 0 0 0 0 0 0 JS0 BS0 BFT0 0 0 0 −LF 0 0 BST0 0 −LS 0 0 0 −LFT−LST0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∆˙ vF ∆pF ∆¨ uS λF λS ∆¨ uI ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = − ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˆ rF m rF c ˆ rS 0 0 0 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (63) where the contributions from the fluid subsystem are given by the expressions: JF mm =αmMF+αfδvCF,JF mc =GF(64) JF cm =αmTF+αfδvDF,JF cc =SF(65) and for the structure: JS=αmMS+αfδuKS,(66) with some secondary terms, related with the stabilization residuals, that have been discarded. Strictly speaking, the solution scheme should be considered a quasi-direct Newton method because the position of the fluid ALE mesh should be updated between iterations. The different steps of a monolithic solution algorithm are described in Algorithm 1. Algorithm 1: Monolithic quasi-direct FSI solver with LLMs Assume that the solution process has progressed until time level n. Apply Newton–Raphson predictor–corrector method for the new solution at time n+1 using the following steps. Step 1: Predictor step. At this stage, given the solution at time level tn, we predict the solution at time level tn+1using previous information: ⎧ ⎪ ⎨ ⎪ ⎩ ˙ vF 0=(γ−1) γ˙ vF n vF 0=vF n pF 0=pF n ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uS 0=(γ−1) γ¨ uS n ˙ uS 0=˙ uS n uS 0=uS p+δu¨ uS 0 ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uI 0=(γ−1) γ¨ uI n ˙ uI 0=˙ uI n uI 0=uI p+δu¨ uI 0 18 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 where the subscript n+1 represents value at the current iteration. This choice corresponds to the same velocity predictor, meaning that the initial guess for the fluid velocity at time level tn+1is the velocity from time level tn. Step 2: Multicorrector step. At this stage, we iterate on the solution until the residual equations are satisfied. For this, we repeat the following steps for i=0,1,...,imax , where iis the iteration counter and imax is the maximum number of nonlinear iterations specified for the current time step. Step 3: Evaluate at the intermediate time levels: {˙ vF n+αm=˙ vF n+αm(˙ vF i−˙ vF n) vF n+αf=vF n+αf(vF i−vF n){¨ uS n+αm=¨ uS n+αm(¨ uS i−¨ uS n) uS n+αf=uS n+αf(uS i−uS n) Step 4: Use the intermediate values to assemble the linear system of equations (63) corresponding to the linearization with respect to the nodal unknowns: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ JF mm JF mc 0 BF0 0 JF cm JF cc 0 0 0 0 0 0 JS0 BS0 BFT0 0 0 0 −LF 0 0 BST0 0 −LS 0 0 0 −LFT−LST0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∆˙ vF i ∆pF i ∆¨ uM i λF i λS i ∆¨ uI i ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ = − ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ˆ rF m rF c ˆ rS 0 0 0 ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ Step 5: Update the solution: ⎧ ⎪ ⎨ ⎪ ⎩ ˙ vF i+1=˙ vF i+∆˙ vF i vF i+1=vF i+δv∆˙ vF i pF i+1=pF i+∆pF i ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uS i+1=¨ uS i+∆¨ uS i ˙ uS i+1=˙ uS i+δv∆¨ uS i uS i+1=uS i+δu∆¨ uS i ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uI i+1=¨ uI i+∆¨ uI i ˙ uI i+1=˙ uI i+δv∆¨ uI i uI i+1=uI i+δu∆¨ uI i Step 6: Calculate new position of the fluid mesh using the interface advance: [KMBM BMT0]{uM i+1 λM i+1}={0 LMuI i+1} Step 7: Evaluate the global residual r(zi+1)and if ∥r(zi+1)∥> ϵ, increment iteration count i=i+1 and return to Step 2 . The dimension of the linear system appearing in Step 4 makes difficult the application of this monolithic solution algorithm to very large problems. Therefore, a partitioned method is preferred. This is the solution approach discussed in the next Section. 10. Partitioned FSI solver with LLMs With the aim of decreasing the computational costs entailed by the monolithic scheme, favoring software modularity and parallelization, we introduce now the so-called interface driven partitioned scheme. We proceed simplifying system (63) by retaining only the fluid, structure and interface accelerations, together with the Lagrange 19 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 multipliers, to yield the following system: ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ JF0 BF0 0 0 JS0 BS0 BFT0 0 0 −LF 0 BST0 0 −LS 0 0 −LFT−LST0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ∆˙ vF ∆¨ uS λF λS ∆¨ uI ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ = −⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ ˆ rF m ˆ rS 0 0 0 ⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭ (67) eliminating the fluid and structure kinematic variables: ⎡ ⎢ ⎣ MF b−10 LF 0 MS b−1LS LFTLST0 ⎤ ⎥ ⎦⎧ ⎨ ⎩ λF λS ∆¨ uI⎫ ⎬ ⎭=⎧ ⎨ ⎩ ˆ aF b ˆ aS b 0⎫ ⎬ ⎭ (68) with subsystem equivalent mass-matrices: MF b=[BFTJF−1BF]−1 ,MS b=[BSTJS−1BS]−1 (69) and corresponding boundary residual accelerations: ˆ aF b= −BFTJF−1ˆ rF m,ˆ aS b= −BSTJS−1ˆ rS(70) for the fluid and the structure, respectively. A solution of the AFETI problem (68) provides the unknown interface displacement increments needed to perform a completely uncoupled NR subiteration solving in parallel the fluid and structure domains. Very efficient iterative parallel algorithms exist to solve this interface problem exactly, see [63] and references therein, however only a sufficiently good approximation of its solution is required to continue with the NR iterations. This distinct aspect of the partitioned algorithm is described next. 10.1. Approximate solution of the interface problem We can go one step further in (68) by eliminating the localized Lagrange multipliers to finally obtain the dynamic equations of motion for the interface: [MF I+MS I]∆¨ uI=LFTMF bˆ aF b+LSTMS bˆ aS b(71) where the interface equivalent mass-matrices emerge from the fluid and structure side defined as: MF I=LFTMF bLF,MS I=LSTMS bLS(72) that are obtained as projections of the individual system equivalent mass-matrices on the non-matching interface using the coupling operators. The assembly and solution of the interface system (71) depends significantly on the number of interface degrees of freedom. Normally, this number is closer to the coarser boundary discretization allowing to obtain a good computational model of the fluid–structure interaction at relatively low computational cost. Remark 2. The interface dynamic equilibrium equation (71) reveals that the interface acceleration increment is calculated from an equivalent mass-matrix average of the fluid and structure interface residual accelerations (70). These residuals can be computed independently and in parallel. However, to proceed with the solution and update of the interface, we also need to calculate the boundary equivalent mass-matrices (69) and these terms are computationally expensive. We have seen that the main difficulty for an efficient solution the interface problem lies in the computation of the boundary equivalent mass-matrices. Nevertheless, according to the interface dynamic equilibrium equation, these matrices are only needed to calculate the interface acceleration by performing a mass-average of the system interface accelerations and then continue with the NR iterations. This means that the mass-average part of the serial 20 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 7. Proposed interface-centered coupled four-field approach for the solution of FSI problems with LLMs. First, fluid and structure interface residuals are received by the Interface system to calculate its displacements uI. Then, the interface boundary advances are sent to the Structure and Dynamic Mesh systems. Finally, the state of the Structure and the Fluid systems are recalculated independently and in parallel. computation could be approximated, while the interface accelerations are calculated independently and in parallel. For the former task, we construct the following approximations of the equivalent boundary mass matrices: M⊛ b=[B⊛TJ⊛−1B⊛]−1 ≈{B⊛TJ⊛B⊛lumped approximation B⊛TS⊛B⊛Dirichlet approximation (73) with superscript ⊛indicating either fluid or structure domains and where S⊛is the Schur complement of the Jacobian matrix reduced to the interface DOFs. Note that these are the two classical preconditioners used in AFETI methods for the iterative solution of coupled interface problems. Although the Dirichlet approximation is more accurate, a computation of the Schur complement of the Jacobian for every NR iteration is very costly. Due to its economy, the lumped approximation of the equivalent boundary mass-matrix is preferred and adopted in this work. The proposed partitioned algorithm is schematically represented in Fig. 7 and described in Algorithm 2. Algorithm 2: Partitioned FSI solver with LLMs Assume that the solution process has progressed until time level n. Apply Newton–Raphson predictor–corrector method for the new solution at time n+1 using the following steps. Step 1: Predictor step. At this stage, given the solution at time level tn, we predict the solution at time level tn+1using previous information: ⎧ ⎪ ⎨ ⎪ ⎩ ˙ vF 0=(γ−1) γ˙ vF n vF 0=vF n pF 0=pF n ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uS 0=(γ−1) γ¨ uS n ˙ uS 0=˙ uS n uS 0=uS p+δu¨ uS n+1 ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uI 0=(γ−1) γ¨ uI n ˙ uI 0=˙ uI n uI 0=uI p+δu¨ uI n+1 where the subscript n+1 represents value at the current iteration. This choice corresponds to the same velocity predictor, meaning that the initial guess for the fluid velocity at time level tn+1is the velocity from time level tn. Step 2: Multicorrector step. At this stage, we iterate on the solution until the residual equations are satisfied. For this, we repeat the following steps for i=0,1,...,imax , where iis the iteration counter and imax is the maximum number of nonlinear iterations specified for the current time step. Step 3: Parallel stage. Compute interface residual accelerations. Step 3.1: Evaluate the iterates at the intermediate time levels: 21 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 {˙ vF n+αm=˙ vF n+αm(˙ vF i−˙ vF n) vF n+αf=vF n+αf(vF i−vF n){¨ uS n+αm=¨ uS n+αm(¨ uS i−¨ uS n) uS n+αf=uS n+αf(uS i−uS n) Step 3.2: Solve in parallel fluid and structure to find the boundary accelerations: [JF mm JF mc JF cm JF cc ]{ˆ aF ˆ pF}={−ˆ rF m −rF c}→ˆ aF b=BFTˆ aF JSˆ aS= −ˆ rS b→ˆ aS b=BSTˆ aS Step 4: Serial stage. Solve interface problem to update interface location. Step 4.1: Solve the interface dynamic equilibrium equation for the interface accelerations update: [MF I+MS I]∆¨ uI i=LFTMF bˆ aF b+LSTMS bˆ aS b Step 4.2: Actualize the position of the interface at time tn+1: ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uI i+1=¨ uI i+∆¨ uI i ˙ uI i+1=˙ uI i+δv∆¨ uI i uI i+1=uI i+δu∆¨ uI i Step 5: Parallel stage. Update fluid and structure with new interface location. Step 5.1: Calculate new position of the fluid mesh using the interface advance: [KMBM BMT0]{uM i+1 λM i+1}={0 LMuI i+1} Step 5.2: Solve in parallel the fluid and the structure with the updated interface accelerations: ⎡ ⎢ ⎣ JF mm JF mc BF JF cm JF cc 0 BFT0 0 ⎤ ⎥ ⎦⎧ ⎪ ⎨ ⎪ ⎩ ∆˙ vF i ∆pF i λF i ⎫ ⎪ ⎬ ⎪ ⎭=⎧ ⎨ ⎩ −ˆ rF m −rF c LF∆¨ uI i ⎫ ⎬ ⎭[JSBS BST0]{∆¨ uS i λS i}={−ˆ rS LS∆¨ uI i} Step 5.3: Update the tn+1solution of the fluid and the structure: ⎧ ⎪ ⎨ ⎪ ⎩ ˙ vF i+1=˙ vF i+∆˙ vF i vF i+1=vF i+δv∆˙ vF i pF i+1=pF i+∆pF i ⎧ ⎪ ⎨ ⎪ ⎩ ¨ uS i+1=¨ uS i+∆¨ uS i ˙ uS i+1=˙ uS i+δv∆¨ uS i uS i+1=uS i+δu∆¨ uS i Step 6: Evaluate the global residual r(zi+1)and if ∥r(zi+1)∥> ϵ, increment iteration count i=i+1 and return to Step 2. The proposed partitioned is composed of three important stages. In Step 3.2 fluid and structure are solved independently without boundary conditions to find the free interface boundary accelerations. Then, Step 4.1 solves 22 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 8. Submerged cylinder. Geometry configuration with dimensions and boundary conditions (left) and finite element mesh used for the computations with 4996 linear triangular elements and 2548 nodes (right). A very refined mesh is needed around the cylinder to capture the drag force accurately. the dynamic interface equation of motion using previously calculated boundary forces coming from the fluid and the structure. Finally, once it is known the new interface position, Step 5.2 solves separately fluid and structure with known interface boundary accelerations. Remark 3. Note that in the proposed partitioned Algorithm 2, the fluid and structure localized Lagrange multipliers have been eliminated and they do not participate in the solution process. They only appear separately in Step 5.2 as the boundary forces needed to enforce the boundary constraints dictated by the known interface motion 11. Numerical examples The following examples serve to illustrate the capabilities verify the results of the proposed partitioned formulation for modeling and simulation of FSI problems. The tests are presented in ascending order of complexity, including FSI cases of rigid bodies combined with Eulerian and ALE descriptions of the fluid, followed by more complicated FSI problems with flexible structures. In the simulations, ALE-VMS formulation is applied for the fluid with triangular finite elements and equal order for velocity and pressure using the stabilization parameter described in (27). For the fluid mesh updating, artificial elastic properties EM e=1 and νM=0.3 have been used as defined in (23) with linear variation of the stiffness based on the element size. The structure displacements are discretized using four-node quadrilateral finite elements and the solid material model in the flexible cases is compressible Neo-Hookean (12). Time integration is always made with the Generalized-αmethod using an spectral radius ρ∞=0.5 which provides an adequate balance between accuracy and robustness. 11.1. Cylinder accelerating on still fluid This first example is taken from [41] and considers the case of a rigid cylinder of diameter D=1 m submerged in still fluid that accelerates under the action of an increasing external force f(t), see Fig. 8(left). The problem is solved using a reference system moving with the cylinder and adopting an Eulerian description of the fluid motion around the structure. This means that the fluid mesh is fixed with uM=0and no dynamic mesh solver is required. For the fluid properties, density ρF=1 kg/m3and viscosity µF=0.01 N m/s are considered. On the structure side, it is assumed that the cylinder has density ρS=1 kg/m3, hence its total mass is mS=ρSVS where VS=1 4πD2is the volume of the cylinder per unit length. The structure is modeled as a free point mass mS 23 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 9. Cylinder accelerating on still fluid. Velocity of the cylinder compared with the analytical Morison’s formula computed with added mass coefficient Ca=1 and drag coefficient Cd=1 (left). Contour of velocity magnitudes with fully developed flow around the cylinder at time t=20 s (right). with an applied external force fSand a reaction from the fluid λS. Therefore, the structure residual can be written: rS(uS, λS)=mS¨uS−fS+λS=0 (74) where the external force applied on the structure is considered to vary in time according to the law: fS(t)={0.3(1−cos (πt/t1)) for t≤t1 0.6 for t>t1 (75) increasing monotonically up to time t1=10 s and then maintained constant until a stationary equilibrium is reached. Then, the fluid–structure interface is modeled as single point with velocity ˙uI= ˙uSthat is connected through the coupling matrix LF=[0 1 . . . 0 1]Tto all the vertical components of the fluid boundary velocities by a group of localized Lagrange multipliers λF. A theoretical solution of this problem can be approximated by adopting the Morison semi-empirical equation to model the fluid reaction acting on the cylinder [64]. The equation of motion is then obtained by equilibrating the external force with the Morison force to yield the expression: fS(t)=fi+fa+fd=ρSVS¨uS+ρFCaVS¨uS+1 2ρFCdAS˙uS|˙uS|(76) where fi=ρSVS¨uSis the cylinder inertia force, fa=ρFCaVS¨uSis the hydrodynamic mass force with Cathe added mass coefficient, fd=1 2ρFCdAS˙uS|˙uS|is the drag force with the drag coefficient Cdand AS=Dis the frontal area per unit length, e.g. the cross-sectional area of the body perpendicular to the flow direction. The solution of this ODE for constant added mass coefficient Ca=1 and drag coefficient Cd=1 is presented as Morison in Fig. 9(left). It is observed that the velocity gained by the cylinder is approximately ˙uS=1 m/s when the external force reaches its maximum at instant t1and it rapidly decelerates until the constant external force is equilibrated with the fluid drag force. This equilibrium condition occurs with a final steady velocity of the cylinder ˙uS=√fS(t1)/(1 2ρFCdAS)=1.095 m/s corresponding to a Reynolds number Re =109.5. Note that Morison solution is only an approximation because the drag coefficient Cdvaries with the Reynolds number, an effect that can be taken into account by including a finite element model for the fluid. For the finite element solution of the problem, see Fig. 8(right), a square domain of side H/D=30 is considered around the cylinder that is discretized using a mesh of 2548 nodes and 4996 stabilized triangular elements with linear approximation of velocities and pressure. The fluid mesh presents a strong refinement near the cylinder surface in order to capture correctly the boundary layer and obtain a good approximations of the drag force during the cylinder acceleration. The maximum–minimum ratio of mesh element sizes is 180 and no slip boundary conditions 24 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 are enforced on the interface. For the external fluid boundaries, zero vertical velocity is assumed on the top and bottom sides of the fluid domain together with outflow boundary condition on the left side. The coupled FSI equations are successfully integrated with time-step ∆t=0.1 s by performing only N=3 NR iterations per time-step. The structure velocity is compared with Morison in Fig. 9(left) where it is observed that the obtained drag force is perfectly stable. It is also presented in Fig. 9(right) a snapshot at time t=20 s of the fluid velocity magnitudes around the cylinder. At this particular instant, with a steady state velocity ˙uS=1.21 m/s, the obtained numerical drag coefficient is Cd≈0.82. The numerical solution approximates well the Morison solution during the cylinder acceleration until this instant. Note that the obtained value of drag coefficient can be influenced by the finite size of the cavity and also affected by a poor resolution of the wake behind the cylinder at higher Reynolds. 11.2. Vortex-induced vibrations of an elastically supported cylinder In this second example, we analyze the stability of the proposed scheme by studying the case of a rigid cylinder in cross flow. The problem considers a rigid cylinder of diameter D=0.1 m, attached to an elastic spring of stiffness kSand immersed in a two-dimensional channel of length L/D=20 and height H/D=10 as sketched in Fig. 10(top). The cylinder is placed in the middle of the channel section with the center located at coordinate x/D=5 downstream of the inflow section and allowed to move with one degree of translational freedom in the vertical direction. The analytical solution for a fully developed laminar channel flow is used as inflow condition, with a parabolic velocity profile given by vx(y)/vmax =(y/H)(2−y/H) with maximum value vmax =15 m/s and mean ¯vx=2 3vmax . A relatively coarse unstructured mesh of 1470 triangular linear elements and 803 nodes is used to model the fluid, see Fig. 10(bottom). No-slip boundary conditions are enforced at the channel walls and outlet conditions are prescribed on the right side, allowing vortices to leave the integration domain without significant disturbances. The fluid is assumed incompressible with density ρF=1 kg/m3and viscosity µF=0.01 N m/s. Based on the mean inflow velocity, the Reynolds number of the flow is Re =100. Experimentally, for fixed cylinders at Reynolds number above 40, small oscillations in the wake induce asymmetry and finally the wake starts shedding vortices into the stream. This situation is termed as onset of periodicity and the wake keeps on undulating up to a Reynolds number of 90 when eddies are shed alternately from top and bottom sides of the cylinder and a regular pattern of alternately shed clockwise and counterclockwise vortices form the Von Karman vortex street. In non-dimensional form, the vortex shedding frequency is expressed as St =ωD/¯vxknown as the Strouhal number. The Strouhal number for fixed cylinders shows a slight but continuous variation with Reynolds number [65], reaching St =0.165 for Re =100 which represents an alternating vortex detachment with frequency 21 Hz. Our cylinder with density ρS=2 kg/m3and mass mS=ρSVS=0.0157 kg is attached to a linear spring of stiffness kS=100 N/m and restricted to move in the vertical direction. This means that the in-vacuo natural frequency of the structure is ωS n=12.7 Hz. Based on this frequency and in order to obtain a resolution of about 15 points per oscillation, a time-step integration ∆t=1/200 s is selected. Then, the problem is solved for N=(2,5,10) using different number of NR iterations per time-step. The results presented in Fig. 11(left) show the history of the vertical displacement obtained for the cylinder, where it is observed that the coupled system behaves as an undamped harmonic oscillator. After a short initial transient of about 0.3 s, a stable harmonic oscillation is achieved with a frequency of approximated frequency ωS=15.2 Hz so ωS/ωS n=1.2. It is also shown in Fig. 11(right) the last interface displacement increment performed by the cylinder after N={2,4,8}NR iterations of the partitioned solver. The problem is solved for three different numbers N={2,4,8}of NR iterations and it is observed that the obtained results are very close for the three cases. In Fig. 11(left) it is represented the cylinder displacement history obtained for N=4, where it is observed that the coupled system behaves as an undamped harmonic oscillator. After a short initial transient of 0.3 s it is achieved a stable harmonic oscillation of frequency ωS=15.2 Hz which represents a numerical Strouhal number St =0.152 and a forcing frequency ratio ωS/ωS n=1.2 over the structure natural frequency. The increments of interface displacement obtained after NNR iterations are represented with time in Fig. 11(right). Approximately, one order of magnitude is reduced per iteration until N=4 on the interface displacement and we are moving the interface at the scale of 10−4. It is noted that the convergence is not quadratic because the interface problem solved in each NR iteration is not exact. 25 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 Fig. 18. Flow-induced vibration of a flexible strip for Reynolds number Re =150. Snapshots with detail of the structure deformation and contours of fluid velocity magnitude. Severe deformation of the structure causes intense fluid mesh distortion that requires special treatment. Declaration of competing interest The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Jose A. Gonzalez reports financial support was provided by Consejer´ ıa de Transformaci´ on Econ´ omica, Conocimiento, Empresas y Universidades de la Junta de Andaluc´ ıa. Acknowledgment The work of Jos´ e A. Gonz´ alez was supported by Consejer´ ıa de Transformaci´ on Econ´ omica, Conocimiento, Empresas y Universidades de la Junta de Andaluc´ ıa (Spain) through the research project P18-RT-3128. References [1] K.C. Park, C. Felippa, J. Deruntz, Stabilization of staggered solution procedures for fluid–structure interaction analysis, Am. Soc. Mech. Eng. Appl. Mech. Div. AMD 26 (1977). [2] C. Felippa, K. Park, Staggered transient analysis procedures for coupled mechanical systems: Formulation, Comput. Methods Appl. Mech. Engrg. 24 (1980) 61–111, http://dx.doi.org/10.1016/0045-7825(80)90040-7. [3] C.A. Felippa, T.L. Geers, Partitioned analysis for coupled mechanical systems, Eng. Comput. 5 (1988) 123–133, http://dx.doi.org/10. 1108/eb023730. [4] S. Piperno, C. Farhat, B. Larrouturou, Partitioned procedures for the transient solution of coupled aroelastic problems Part I: Model problem, theory and two-dimensional application, Comput. Methods Appl. Mech. Engrg. 124 (1995) 79–112, http://dx.doi.org/10.1016/ 0045-7825(95)92707-9. 32 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 [5] S. Piperno, C. Farhat, Partitioned procedures for the transient solution of coupled aeroelastic problems – Part II: energy transfer analysis and three-dimensional applications, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3147–3170, http://dx.doi.org/10.1016/S0045- 7825(00)00386-8. [6] C.A. Felippa, K.C. Park, C. Farhat, Partitioned analysis of coupled mechanical systems, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3247–3270, http://dx.doi.org/10.1016/S0045-7825(00)00391-1. [7] G. Hou, J. Wang, A. Layton, Numerical methods for fluid–structure interaction - a review, Commun. Comput. Phys. 12 (2012) 337–377, http://dx.doi.org/10.4208/cicp.291210.290411s. [8] T.E. Tezduyar, K. Takizawa, Y. Bazilevs, Fluid–Structure Interaction and Flows with Moving Boundaries and Interfaces, John Wiley & Sons, Ltd, 2017, pp. 1–53, http://dx.doi.org/10.1002/9781119176817.ecm2069. [9] C. Farhat, K.G. van der Zee, P. Geuzaine, Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity, Comput. Methods Appl. Mech. Engrg. 195 (2006) 1973–2001, http://dx.doi.org/10.1016/j.cma.2004.11.031, fluid-Structure Interaction. [10] C. Farhat, CFD-Based Nonlinear Computational Aeroelasticity, John Wiley & Sons, Ltd, 2017, pp. 1–21, http://dx.doi.org/10.1002/ 9781119176817.ecm2063. [11] C. Förster, W.A. Wall, E. Ramm, Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows, Comput. Methods Appl. Mech. Engrg. 196 (2007) 1278–1293, http://dx.doi.org/10.1016/j.cma.2006.09.002. [12] Y. Bazilevs, V.M. Calo, Y. Zhang, T.J.R. Hughes, Isogeometric fluid–structure interaction analysis with applications to arterial blood flow, Comput. Mech. 38 (2006) 310–322, http://dx.doi.org/10.1007/s00466-006-0084-3. [13] Y. Bazilevs, V.M. Calo, T.J.R. Hughes, Y. Zhang, Isogeometric fluid–structure interaction: theory, algorithms, and computations, Comput. Mech. 43 (2008) 3–37, http://dx.doi.org/10.1007/s00466-008-0315-x. [14] M.W. Gee, U. Küttler, W.A. Wall, Truly monolithic algebraic multigrid for fluid–structure interaction, Internat. J. Numer. Methods Engrg. 85 (2011) 987–1016, http://dx.doi.org/10.1002/nme.3001. [15] T. Richter, A monolithic geometric multigrid solver for fluid–structure interactions in ale formulation, Internat. J. Numer. Methods Engrg. 104 (2015) 372–390, http://dx.doi.org/10.1002/nme.4943. [16] M. Mayr, T. Klöppel, W.A. Wall, M.W. Gee, A temporal consistent monolithic approach to fluid–structure interaction enabling single field predictors, SIAM J. Sci. Comput. 37 (2015) B30–B59, http://dx.doi.org/10.1137/140953253. [17] S. Bordère, J.-P. Caltagirone, A unifying model for fluid flow and elastic solid deformation: A novel approach for fluid–structure interaction, J. Fluids Struct. 51 (2014) 344–353, http://dx.doi.org/10.1016/j.jfluidstructs.2014.09.010. [18] J. Fan, H. Liao, R. Ke, E. Kucukal, U.A. Gurkan, X. Shen, J. Lu, B. Li, A monolithic lagrangian meshfree scheme for fluid–structure interaction problems within the otm framework, Comput. Methods Appl. Mech. Engrg. 337 (2018) 198–219, http://dx.doi.org/10.1016/ j.cma.2018.03.031. [19] B. Schott, C. Ager, W. Wall, A monolithic approach to fluid–structure interaction based on a hybrid eulerian-ale fluid domain decomposition involving cut elements, Internat. J. Numer. Methods Engrg. 119 (2019) 208–237, http://dx.doi.org/10.1002/nme.6047. [20] R.A. Sauer, T. Luginsland, A monolithic fluid–structure interaction formulation for solid and liquid membranes including free-surface contact, Comput. Methods Appl. Mech. Engrg. 341 (2018) 1–31, http://dx.doi.org/10.1016/j.cma.2018.06.024. [21] V. Gravemeier, S.M. Civaner, W.A. Wall, A partitioned-monolithic finite element method for thermo-fluid–structure interaction, Comput. Methods Appl. Mech. Engrg. 401 (2022) 115596, http://dx.doi.org/10.1016/j.cma.2022.115596. [22] P. Le Tallec, J. Mouro, Fluid structure interaction with large structural displacements, Comput. Methods Appl. Mech. Engrg. 190 (2001) 3039–3067, http://dx.doi.org/10.1016/S0045-7825(00)00381-9. [23] S. Badia, F. Nobile, C. Vergara, Fluid–structure partitioned procedures based on robin transmission conditions, J. Comput. Phys. 227 (2008) 7027–7051, http://dx.doi.org/10.1016/j.jcp.2008.04.006. [24] M.M. Joosten, W.G. Dettmer, D. Peri´ c, Analysis of the block gauss–seidel solution procedure for a strongly coupled model problem with reference to fluid–structure interaction, Internat. J. Numer. Methods Engrg. 78 (2009) 757–778, http://dx.doi.org/10.1002/nme.2503. [25] U. Küttler, W.A. Wall, Fixed-point fluid–structure interaction solvers with dynamic relaxation, Comput. Mech. 43 (2008) 61–72, http://dx.doi.org/10.1007/s00466-008-0255-5. [26] J. Degroote, Partitioned simulation of fluid–structure interaction, Arch. Comput. Methods Eng. 20 (2013) 185–238, http://dx.doi.org/ 10.1007/s11831-013-9085-5. [27] C. Kassiotis, A. Ibrahimbegovic, R. Niekamp, H.G. Matthies, Nonlinear fluid–structure interaction problem. Part I: implicit partitioned algorithm, nonlinear stability proof and validation examples, Comput. Mech. 47 (2011) 305–323, http://dx.doi.org/10.1007/s00466-010- 0545-6. [28] C. Kassiotis, A. Ibrahimbegovic, R. Niekamp, H.G. Matthies, Nonlinear fluid–structure interaction problem. Part II: space discretization, implementation aspects, nested parallelization and application examples, Comput. Mech. 47 (2011) 335–357, http://dx.doi.org/10.1007/ s00466-010-0544-7. [29] U. Küttler, M. Gee, C. Förster, A. Comerford, W.A. Wall, Coupling strategies for biomedical fluid–structure interaction problems, Int. J. Numer. Methods Biomed. Eng. 26 (2010) 305–321, http://dx.doi.org/10.1002/cnm.1281. [30] H.G. Matthies, J. Steindorf, Partitioned strong coupling algorithms for fluid–structure interaction, Comput. Struct. 81 (2003) 805–812, http://dx.doi.org/10.1016/S0045-7949(02)00409-1. [31] R.K. Jaiman, V. Joshi, Computational Mechanics of Fluid-Structure Interaction: Computational Methods for Coupled Fluid-Structure Analysis, Springer, 2022, http://dx.doi.org/10.1007/978-981-16-5355-1. [32] T. Klöppel, A. Popp, U. Küttler, W.A. Wall, Fluid–structure interaction for non-conforming interfaces based on a dual mortar formulation, Comput. Methods Appl. Mech. Engrg. 200 (2011) 3111–3126, http://dx.doi.org/10.1016/j.cma.2011.06.006. [33] M.R. Ross, C.A. Felippa, K.C. Park, M.A. Sprague, Treatment of acoustic fluid–structure interaction by localized Lagrange multipliers: formulation, Comput. Methods Appl. Mech. Engrg. 197 (2008) 3057–3079, http://dx.doi.org/10.1016/j.cma.2008.02.017. 33 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 [34] L. Rodríguez, J.A. González, A. Cerrato, Partitioned solution strategies for coupled BEM-FEM acoustic fluid–structure interaction problems, Comput. Struct. 152 (2015) 45–58, http://dx.doi.org/10.1016/j.compstruc.2015.02.018. [35] K.C. Park, R. Ohayon, C.A. Felippa, J.A. González Pérez, Partitioned formulation of internal and gravity waves interacting with flexible structures, Comput. Methods Appl. Mech. Engrg. 199 (2010) 723–733, http://dx.doi.org/10.1016/j.cma.2009.11.005. [36] J.A. González, K.C. Park, I. Lee, C.A. Felippa, R. Ohayon, Partitioned vibration analysis of internal fluid–structure interaction problems, Internat. J. Numer. Methods Engrg. 92 (2012) 268–300, http://dx.doi.org/10.1002/nme.4336. [37] W. Dettmer, D. Peri´ c, A computational framework for fluid–structure interaction: Finite element formulation and applications, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5754–5779, http://dx.doi.org/10.1016/j.cma.2005.10.019. [38] W.G. Dettmer, D. Peri´ c, A fully implicit computational strategy for strongly coupled fluid–solid interaction, Arch. Comput. Methods Eng. 14 (2007) 205–247, http://dx.doi.org/10.1007/s11831-007-9006-6. [39] P. Saksono, W. Dettmer, D. Peri´ c, An adaptive remeshing strategy for flows with moving boundaries and fluid–structure interaction, Internat. J. Numer. Methods Engrg. 71 (2007) 1009–1050, http://dx.doi.org/10.1002/nme.1971. [40] W.G. Dettmer, D. Peri´ c, A new staggered scheme for fluid–structure interaction, Internat. J. Numer. Methods Engrg. 93 (2013) 1–22, http://dx.doi.org/10.1002/nme.4370. [41] W.G. Dettmer, A. Lovri´ c, C. Kadapa, D. Peri´ c, New iterative and staggered solution schemes for incompressible fluid–structure interaction based on Dirichlet-Neumann coupling, Internat. J. Numer. Methods Engrg. (2020) http://dx.doi.org/10.1002/nme.6494. [42] P. Tonon, R.A. Sanches, K. Takizawa, T.E. Tezduyar, A linear-elasticity-based mesh moving method with no cycle-to-cycle accumulated distortion, Comput. Mech. 67 (2021) 413–434, http://dx.doi.org/10.1007/s00466-020-01941-y. [43] A. Shamanskiy, B. Simeon, Mesh moving techniques in fluid–structure interaction: robustness, accumulated distortion and computational efficiency, Comput. Mech. 67 (2021) 583–600, http://dx.doi.org/10.1007/s00466-020-01950-x. [44] T.J. Hughes, G.R. Feijóo, L. Mazzei, J.-B. Quincy, The variational multiscale method-a paradigm for computational mechanics, Comput. Methods Appl. Mech. Engrg. 166 (1998) 3–24, http://dx.doi.org/10.1016/S0045-7825(98)00079-6. [45] Y. Bazilevs, K. Takizawa, T.E. Tezduyar, Computational Fluid-Structure Interaction: Methods and Applications, John Wiley & Sons, Ltd., 2013, http://dx.doi.org/10.1002/9781118483565. [46] T.E. Tezduyar, Stabilized finite element formulations for incompressible flow computations, Adv. Appl. Mech. 28 (1991) 1–44, http://dx.doi.org/10.1016/S0065-2156(08)70153-4. [47] K.C. Park, C.A. Felippa, G. Rebel, A simple algorithm for localized construction of non-matching structural interfaces, Internat. J. Numer. Methods Engrg. 53 (2002) 2117–2142, http://dx.doi.org/10.1002/nme.374. [48] M.R. Ross, M.A. Sprague, C.A. Felippa, K. Park, Treatment of acoustic fluid–structure interaction by localized lagrange multipliers and comparison to alternative interface-coupling methods, Comput. Methods Appl. Mech. Engrg. 198 (2009) 986–1005, http: //dx.doi.org/10.1016/j.cma.2008.11.006. [49] M.A. Crisfield, Re-visiting the contact patch test, Internat. J. Numer. Methods Engrg. 48 (2000) 435–449, http://dx.doi.org/10.1002/ (SICI)1097-0207(20000530)48:3%3C435::AID-NME891%3E3.0.CO;2-V. [50] G. Rebel, K.C. Park, C.A. Felippa, A contact formulation based on localized lagrange multipliers: formulation and application to two-dimensional problems, Internat. J. Numer. Methods Engrg. 54 (2002) 263–297, http://dx.doi.org/10.1002/nme.426. [51] J.A. González, K.C. Park, FEM and BEM coupling in elastostatics using localized Lagrange multipliers, Internat. J. Numer. Methods Engrg. (2007) http://dx.doi.org/10.1002/nme.1833. [52] Y.U. Song, S.K. Youn, K.C. Park, A gap element for treating non-matching discrete interfaces, Comput. Mech. 56 (2015) 551–563, http://dx.doi.org/10.1007/s00466-015-1186-6. [53] Y.U. Song, S.K. Youn, K.C. Park, Virtual gap element approach for the treatment of non-matching interface using three-dimensional solid elements, Comput. Mech. 60 (2017) 585–594, http://dx.doi.org/10.1007/s00466-017-1423-2. [54] Y.U. Song, G.E. Jeong, S.K. Youn, K. Park, Virtual tetrahedral gap element to connect three-dimensional non-coincident interfaces, Finite Elem. Anal. Des. 152 (2018) 18–26, http://dx.doi.org/10.1016/j.finel.2018.08.005. [55] G.E. Jeong, Y.U. Song, S.K. Youn, K. Park, A new approach for nonmatching interface construction by the method of localized lagrange multipliers, Comput. Methods Appl. Mech. Engrg. 361 (2020) 112728, http://dx.doi.org/10.1016/j.cma.2019.112728. [56] J. Chung, G.M. Hulbert, A time integration algorithm for structural dynamics with improved numerical dissipation: The generalized-α method, J. Appl. Mech. 60 (1993) 371–375, http://dx.doi.org/10.1115/1.2900803. [57] H.M. Hilber, T. Hughes, R.L. Taylor, Improved numerical dissipation for time integration algorithms in structural dynamics, Earthq. Eng. Struct. Dyn. 5 (1977) 283–392, http://dx.doi.org/10.1002/eqe.4290050306. [58] K.E. Jansen, C.H. Whiting, G.M. Hulbert, A generalized-αmethod for integrating the filtered Navier–Stokes equations with a stabilized finite element method, Comput. Methods Appl. Mech. Engrg. 190 (2000) 305–319, http://dx.doi.org/10.1016/S0045-7825(00)00203-6. [59] M.M. Joosten, W.G. Dettmer, D. Peri´ c, On the temporal stability and accuracy of coupled problems with reference to fluid–structure interaction, Internat. J. Numer. Methods Fluids 64 (2010) 1363–1378, http://dx.doi.org/10.1002/fld.2333. [60] T.E. Tezduyar, S. Sathe, R. Keedy, K. Stein, Space–time finite element techniques for computation of fluid–structure interactions, Comput. Methods Appl. Mech. Engrg. 195 (2006) 2002–2027, http://dx.doi.org/10.1016/j.cma.2004.09.014. [61] T.E. Tezduyar, S. Sathe, K. Stein, Solution techniques for the fully-discretized equations in computation of fluid–structure interactions with the space–time formulations, Comput. Methods Appl. Mech. Engrg. 195 (2006) 5743–5753, http://dx.doi.org/10.1016/j.cma.2005. 08.023. [62] T.E. Tezduyar, S. Sathe, Modeling of fluid–structure interactions with the space–time finite elements: Solution techniques, Internat. J. Numer. Methods Fluids 54 (2007) 855–900, http://dx.doi.org/10.1002/fld.1430. [63] J.A. González, K. Park, Accelerating the convergence of AFETI partitioned analysis of heterogeneous structural dynamical systems, Comput. Methods Appl. Mech. Engrg. 360 (2020) 112726, http://dx.doi.org/10.1016/j.cma.2019.112726. 34 J.A. Gonz´ alez and K.C. Park Computer Methods in Applied Mechanics and Engineering 414 (2023) 116134 [64] J. Morison, J. Johnson, S. Schaaf, The force exerted by surface waves on piles, J. Pet. Technol. 2 (1950) 149–154, http://dx.doi.org/ 10.2118/950149-G. [65] C.H.K. Williamson, The natural and forced formation of spot-like ’vortex dislocations’ in the transition of a wake, J. Fluid Mech. 243 (1992) 393–441, http://dx.doi.org/10.1017/S0022112092002763. [66] W.A. Wall, Fluid-Struktur-Interaktion mit Stabilisierten Finiten Elementen (Ph.D. thesis), Universität Stuttgart, 1999, http://dx.doi.org/ 10.18419/opus-127. 35