This paper presents a staggered scheme with second-order temporal accuracy for ﬂuid-structure interaction problems involving ultra-lightweight rigid bodies. The staggered scheme is based on the Dirichlet-Neumann coupling and is non-intrusive. First, the spectral properties of the staggered scheme are studied and also com- pared against the monolithic scheme using a linear model problem. Later, the suitability and effectiveness of the staggered scheme for problems involving incompressible ﬂows and lightweight rigid solids are illustrated by using the examples of galloping of a square cylinder and lock-in of a circular cylinder for mass-ratio values as low as 0.01. This is the ﬁrst time in the literature ﬂow-induced vibrations of rigid bodies with such low mass ratio values are successfully simulated using a staggered scheme. Two different ﬂuid solvers are considered to illustrate the non-intrusive nature of the proposed scheme. Guidelines for choosing the relaxation parameter are also provided. With its iteration-free nature and with a single (relaxation) parameter, the proposed stag- gered scheme renders itself as an accurate and computationally efﬁcient scheme for ﬂuid-rigid body interaction problems, including those involving lightweight structures. Keywords: Fluid-structure interaction; partitioned approach; staggered scheme; added-mass; lightweight structures 1. Introduction Fluid-structure interaction (FSI) is a multiphysics phenomenon that is common in mechanical, civil, aerospace, biomedical, ocean and naval engineering. V ortex-induced vibrations (VIV) of chimneys, marine risers and sub- marine pipelines; galloping of transmission cables and bridge decks; wave-structure interaction of ﬂoating and ﬁxed offshore structures; the aeroelastic ﬂutter of aeroplane wings; and blood ﬂow through heart valves and arteries, are a few examples of FSI problems. Computer simulation of FSI is a complex and expensive task, and it poses numerous difﬁculties depend- ing upon the complexity of the problem (in terms of problem geometry, structural deformations, topological changes and the model for ﬂuid ﬂow) and the type of numerical scheme used for the simulation [1–4]. The frequently encountered issues in simulating an FSI problem are those associated with (i) capturing large struc- tural deformations and topological changes in the ﬂuid domain and (ii) overcoming added-mass instabilities. A vast amount of research work has gone into addressing many of the issues associated with large structural deformations and topological changes by adapting the arbitrary Lagrangian-Eulerian (ALE) formulation with or without remeshing [5–9] and immersed/embedded/unﬁtted/ﬁctitious-domain methods [10–27]. However, ro- bust and accurate partitioned FSI algorithms that can successfully overcome signiﬁcant added-mass instabilities are still lacking. Depending upon whether the governing equations of the FSI problems are solved together or separately, the FSI schemes are termed as monolithic or partitioned schemes. In a monolithic FSI scheme [7, 20, 27– 34], the coupled nonlinear governing equations are solved together using an iterative technique, for example, the Newton-Raphson scheme. Granted that monolithic FSI schemes are devoid of added-mass instabilities, their disadvantages are due to the unavoidable requirement for computing off-diagonal (coupling) matrices and their high computational cost. The need to compute the coupling terms in the monolithic schemes becomes Corresponding author Email address: c.kadapa@swansea.ac.uk (Chennakesava Kadapa) Accepted for publication in Ocean Engineering August 10, 2020 a serious bottleneck in combining various ﬂuid and solid solvers in a non-intrusive manner; thus, diminishing their applicability to a limited set of solvers that satisfy certain requirements. On the other hand, in the partitioned/staggered approaches for FSI, the solution for the coupled equations is obtained by solving the ﬂuid and solid sub-problems separately, with an appropriate algorithm for ﬁeld transfer between individual solvers. While the partitioned approaches involve sub-iterations at every time step, the staggered schemes do not. Therefore, partitioned/staggered schemes are not only computationally appealing because of their ﬂexibility in combining various ﬂuid and solid solvers, but also are signiﬁcantly less expensive when compared to the monolithic FSI schemes. However, the main drawback of the partitioned/staggered schemes is their instabilities due to the added-mass effect. The major cause of instabilities in the partitioned FSI algorithms is the added mass [35] which becomes predominant in the FSI problems involving lightweight structures, see Causin et al. [36], Conca et al. [37] and F ¨orster et al. [38], van Brummelen [39] and Dettmer and Peri ´c [42]. The added-mass instabilities can be characterised by the mass ratio ( m=solid mass/displaced ﬂuid mass) for rigid structures, and density ratio (=solid density/ﬂuid density) for ﬂexible structures. The lower the value of m(or), the higher the added mass. The stability of a partitioned/staggered scheme is dictated by the value of m(or). A wide variety of techniques have been proposed towards addressing the added-mass instabilities in the partitioned FSI algorithms. While some techniques solve the FSI problem without using sub-iterations, for ex- ample, the conventional sequential staggered (CSS) scheme [40], others use sub-iterations at each time step until a pre-deﬁned convergence criterion is satisﬁed. Some noteworthy contributions on partitioned FSI schemes are due to Florentie et al. [41], Dettmer and Peri ´c [42], Fern ´andez and Moubachir [43], Burman and Fern ´andez [44, 45], dos Santos et al. [46], Breue et al. [47], Badia et al. [48], Wood et al. [49], Matthies and Steindorf [50], Degroote et al. [29], Farhat et al. [51], Heil et al. [32], Guidoboni et al. [52], Jaiman et al. [53], Michler et al. [54], Causin et al. [36], van Brummelen [39, 55], F ¨orster et al. [38], He [56, 57], He and Zhang [58], He et al. [59, 60], Idelsohn et al. [62], Banks et al. [63, 64], Bukac et al. [65, 66], and Winterstein et al. [67]. In spite of numerous contributions towards addressing the added-mass issues in partitioned/staggered FSI schemes, the accuracy and robustness of these algorithms are still open problems. Note that the majority of the research work is focussed on the FSI problems with ﬂexible lightweight structures, and, of course, it is for all the correct reasons. Nevertheless, the literature on partitioned/staggered approaches that are accurate and robust for the simulation of FSI problems involving lightweight rigid solids is limited. To the best of author’s knowledge, only Jaiman et al. [68], He [56, 57] and Banks et al. [63, 64] successfully demonstrated partitioned approaches to FSI problems involving lightweight rigid solids for which m<1. With the novel nonlinear interface force correction (NIFC) scheme, Jaiman et al. [68, 105] successfully demonstrated vortex-induced vibrations of a circular cylinder for mvalues as low as 0.26. He [69] successfully adapted a partitioned implicit coupling using a displacement-predictor to successfully simulate vortex-induced vibrations of a circular cylinder for m0:298. Banks et al. [63, 64] proposed the added-mass partition (AMP) algorithm that can cope with even zero mass-ratio value. Despite its second-order accuracy, the NIFC scheme [68] requires sub-iterations, and the overall accuracy of the partitioned approach followed by He [69] is ﬁrst-order accurate because the accuracy of the mesh update procedure is only of order one. While the AMP scheme can accommodate even zero mass-ratio value, the need to evaluate the surface integrals for comput- ing the added-mass and added-damping tensors somewhat limits the applicability of the scheme to simplistic geometries. Moreover, the added-mass and added-damping tensors proposed in [63, 64] modify the dynamic characteristics of the original solid sub-system. Due to the lack of examples involving vortex-induced vibra- tions of rigid solids in [63, 64], the effect of the added-mass and added-damping tensors on the amplitude and frequency of oscillations is unknown. To perform accurate and computationally efﬁcient simulations of the FSI problems involving lightweight rigid solids by overcoming the issues mentioned above, we present a new staggered scheme for ﬂuid-rigid body interaction that (i) is second-order accurate in time, (ii) involves no sub-iterations at every time step, (iii) requires only Dirichlet boundary conditions for the ﬂuid problem and Neumann boundary conditions for the solid problem and (v) can be used to develop an FSI framework by coupling ﬂuid and solid solvers in a non- intrusive manner. The proposed staggered scheme is an extension of the scheme proposed by Dettmer and Peri ´c [42] to FSI problems with lightweight rigid solids. For the sake of simplicity, and for future extensions to ﬂuid- structure-control interaction problems, the staggered scheme is reformulated using a single time integration scheme for both the ﬂuid and solid problems. The accuracy of the staggered scheme is demonstrated using the 2 examples of ﬂow-induced vibrations involving rigid solids with mass ratio values as low as 0:01. To the best of author’s knowledge, this is the ﬁrst time in the literature that successful numerical solutions are computed using a staggered scheme for ﬂuid-rigid body interaction problems with such low mvalues. The non-intrusive nature of the proposed scheme is illustrated using two different ﬂuid solvers. The rest of the paper is organised as follows. The governing equations for the ﬂuid-rigid body interac- tion problems are presented in Section 2, followed by the discussion of monolithic and partitioned solution approaches to FSI problems in Section 3. The accuracy and stability of the staggered scheme are assessed by studying the single-degree-of-freedom linear spring-mass-damper system in Section 4. Numerical damping characteristics and condition for the stability of the staggered scheme are also established in Section 4. The remarkable ability of the staggered scheme in computing accurate numerical solutions of FSI problems involv- ing lightweight rigid solids and its ability to combine different solid and ﬂuid solvers in a non-intrusive manner are demonstrated using the example of vortex-induced vibrations of a circular cylinder, galloping of a square body and lightweight rectangular body rising in a liquid, in Section 5. Guidelines for choosing the relaxation parameter are outlined in Section 6. The paper is concluded with Section 7 by providing a summary of key observations made and conclusions drawn. 2. Governing equations of the ﬂuid-rigid body interaction problem This section provides a brief summary of the governing equations for the ﬂuid-structure interaction prob- lems involving isothermal, laminar, viscous, incompressible ﬂows and rigid body dynamics. The domains of the ﬂuid and solid problems are denoted as fand s, respectively, with fandsas their respective boundaries. The interface between the ﬂuid and solid domains is denoted as fs. 2.1. Governing equations for the uid ow problem For isothermal, laminar, viscous and incompressible ﬂows, the initial-boundary value problem is stated as: Find velocity, vf: f!R3; and pressure, pf: f!R, for the given ff: f!R3;vf: f D!R3; and tf: f N!R3, such that: f@vf @t+f(vfr)vffvf+rpf=ffin f; (1a) rvf= 0 in f; (1b) vf=vfon f D; (1c) fnf=tfin f N; (1d) vf(t= 0) = vf 0 in f; (1e) where,tis the time variable, ris the gradient operator, fis the dynamic viscosity of the ﬂuid, fis the density of the ﬂuid, ffis the body force, nfis the unit outward normal on the boundary f,vf 0is the initial velocity of the ﬂuid in the domain, vfis the prescribed velocity on Dirichlet boundary f D, and tfis applied traction on the Neumann boundary f N. The conditions f= f D[f Nandf D\f N=;, apply for the Dirichlet and Neumann boundaries. The pseudo-stress fis deﬁned as, f:=frvfpfI; (2) where, Iis the second-order identity tensor. 2.2. Governing equations for rigid body dynamics The equations governing the dynamics of rigid solids can be written as M as+C vs+K ds=Fs; (3a) ds(t= 0) = ds 0; (3b) vs(t= 0) = vs 0; (3c) 3 where, Mis the mass matrix, Cis the damping matrix, Kis the stiffness matrix, dsis the displacement vector of the rigid body, vs(= dds=dt)is the velocity vector and as(= d2ds=dt2)is the acceleration vector, ds 0 is the initial displacement vector, vs 0is the initial velocity vector and Fsis the force vector. Note that only translational degrees of freedom (DOFs) are considered in this work for the sake of simplicity. The inclusion of rotational DOFs is straightforward, of course, with appropriate modiﬁcations. 2.3. Governing equations at the uid-solid interface At the interface fsbetween the ﬂuid and solid domains, the kinematic constraint and the equilibrium of tractions need to be satisﬁed by the solution variables in the ﬂuid and solid domains. These two conditions can be written as, Kinematic constraint: vf=vs; (4) Traction equilibrium: fnf+sns=0; (5) where, nsis the unit outward normal to the boundary of the solid domain. At the ﬂuid-solid interface, we have, ns=nf. The traction equilibrium (5) can be transformed into the following force equilibrium condition for the rigid-body interaction problem, Ff+Fs=0; (6) where, Ff=Z fsfnfd: (7) 3. Formulation Several numerical formulations exist in the literature for solving the ﬂuid ﬂow problem over moving solids. Such schemes can be broadly classiﬁed into two categories: (i) Arbitrary Lagrangian-Eulerian (ALE) based methods using body-ﬁtted meshes [5–9, 70, 71] and (ii) immersed/embedded/unﬁtted/CutFEM/ﬁctitious- domain methods [10–27]. The speciﬁc details of the formulation used for the ﬂow problem are not particularly important in the context of the present work since the motivation behind the proposed staggered scheme is to be able to develop FSI solvers by coupling different ﬂuid and solid sub-solvers in a non-intrusive fashion. The only requirement is that the ﬂuid solver treats the interface velocity as the Dirichlet boundary condition and the solid solver treats the interface traction as the Neumann boundary condition , which is quite easily satisﬁed by the majority of opensource as well as commercial computational ﬂuid dynamics (CFD) and rigid body dynamics solvers. After selecting the appropriate ﬂuid and solid solvers, the coupled FSI problem given by the equations (1)- (5) can be solved using either a monolithic or a partitioned approach, as discussed in the following subsections. 3.1. Monolithic approach In the monolithic approach, the coupled governing equations (1)-(5) are solved together for the solution ﬁeld, typicallyfvf;pf;vsg, using an iterative procedure such as the Newton-Raphson scheme. After adapting a spatial discretisation technique and appropriate time integration schemes for the ﬂuid and solid sub-problems, the coupled FSI problem can be written in the matrix-vector form, KffKfs KsfKssVf Vs =Rf Rs (8) where, Vfis the vector of DOFs, Kffis the effective stiffness matrix and Rfis the residual vector, for the ﬂuid sub-problem. Vs,KssandRfare the corresponding entries for the solid sub-problem. KfsandKsf are the coupling matrices which are the result of linearisation of the interface terms. The speciﬁc details of the 4 individual matrices in the coupled matrix system (8) are not important, except the fact that the computation of coupling terms, KfsandKsf, is essential for the success of the monolithic scheme. The unique advantage of the monolithic scheme is in its ability to accommodate any mass-ratio ( m) values, includingm= 0. However, in addition to being computationally expensive, monolithic schemes are known to pose serious challenges in their computer implementation, especially in computing the coupling terms as well as achieving good scalability on high-performance computing architectures. Albeit the coupling terms can be computed in the in-house codes for CFD solvers based on body-ﬁtted meshes, this task is not without its programming difﬁculties for complex geometries. Moreover, CFD solvers based on cut-cell methods [22] pose additional challenges in computing the coupling terms due to activation/deactivation of ﬂuid cells from one iteration to the other. Furthermore, it is almost impossible to compute the coupling terms for combining various ﬂuid and solid solvers that adapt different modelling techniques, for example, ﬁnite volume method (FVM) for the ﬂuid problem and ﬁnite element method (FEM) for the solid problem. Therefore, to achieve ﬂexibility in combining various ﬂuid and solid solvers and solve the two-way FSI problem in a computationally efﬁcient manner, partitioned approaches are the preferred choice. 3.2. Partitioned or staggered approach In the partitioned approaches for FSI problems, the solid and ﬂuid sub-systems are solved separately in an iterative manner until a predeﬁned convergence criterion is satisﬁed. In contrast, staggered FSI schemes do not involve sub-iterations [40]. Several possibilities exist in the partitioned approaches for FSI problems; some schemes approximate the coupled matrices or use quasi-Newton methods with ﬁxed-point iterations to achieve convergence [29], and others ignore coupled matrices altogether. Schemes based on inexact Jacobians or quasi-Newton methods are still computationally expensive and pose several challenges in parallelisation for the high-performance computing architectures [72, 73]. In the explicit coupling approaches, for example, the conventional sequentially staggered scheme (CSS) [40] and other similar variations, the ﬂuid and solid problems are solved separately without considering the effects of coupling matrices. This feature is what makes the explicit coupling approaches computationally very appealing because such an approach is relatively easy to implement by developing wrappers for combining various ﬂuid and solid solvers. Nevertheless, the partitioned approaches for FSI problems are not without their drawbacks, among which the instability due to added-mass is the most prominent one. A partitioned or staggered scheme can be designed using different approaches. The majority of the parti- tioned schemes available in the literature are based on displacement predictors for the solid problem [38, 51, 69, 74]. Staggered schemes based on force-predictors for the solid problem also exist [42, 61, 75, 106]. Other variations in partitioned schemes include added-mass partitioned (AMP) scheme by Bank et al. [63, 64], com- bined interface boundary condition by Jaiman et al. [53], and nonlinear interface force correction(NIFC) by Jaiman et al. [68]. The staggered scheme considered in the present work is based on the force predictor for the solid problem. The coupling type employed is the Dirichlet-Neumann coupling in which the interface traction from the ﬂuid solver is applied as a Neumann boundary condition for the solid domain, and the interface velocity is imposed as a Dirichlet boundary condition for the ﬂuid problem. Dirichlet-Neumann coupling is the most intuitive type of coupling because of the nature of applied boundary conditions on the ﬂuid and solid problems. Following the terminology of Felippa et al. [40], the temporal ﬂow diagram for the proposed staggered scheme is depicted in Fig. 1, and the corresponding pseudocode is presented in Fig. 2. 5 F F FS S S P P P PS S S tn−1 tn tn+1∆t ∆tFigure 1: The temporal ow diagram for proposed staggered scheme: Sis the solid problem, Fis the uid problem, P is the predictor step, and S is the substitution step. tn1,tnandtn+1are the time instants, and  tis the time step. Update time: tn+1=tn+ ∆ t Predict the force on the solid: FsP n+1= 2Fs n−Fs n−1 Solve the solid problem Update the ﬂuid mesh Solve the ﬂuid problem Compute the force F∗ n+1from the ﬂuid Correct the force on the solid: Fs n+1=−βF∗ n+1+ [1−β]FsP n+1 Figure 2: Key steps in the staggered scheme. Note that the scheme involves no sub-iterations at each time step. 6 4. Model problem Analysis of a staggered scheme for the FSI problems is quite a challenging task under the general settings of incompressible Navier-Stokes equations. Therefore, in this work, the accuracy and stability characteristics of the staggered scheme are assessed using the single-degree-of-freedom linear spring-mass-damper system and the observations are extended to the FSI problems involving actual CFD solvers. Despite its simplicity, the linear spring-mass-damper system provides a reasonable framework for conducting the spectral analysis of the staggered scheme. The original linear problem consisting of a spring-mass-damper system is shown in Fig. 3(a), and it is modiﬁed into the one with two masses, as shown in Fig. 3(b), by distributing the mass mto the solid (spring) and the ﬂuid (damper) sub-systems such that ms= m;mf= [1 ]m; andms+mf=m; (9) where,msthe mass of the solid sub-system and mfis the mass of the ﬂuid sub-system. The parameter gives an indication of the solid mass with respect to the total mass, and it is a crucial parameter in the study of the staggered scheme. The parameter is related to the mass-ratio parameter ( m=ms=mf) which is well-known in the literature on ﬂuid-rigid solid interaction [76] as, =ms m=m 1 +m: (10) The parameter lies in the range [0;1]for any positive value of mand it approaches zero as mapproaches zero, as shown in Fig. 4. dk cm (a) Standard spring-mass-damper model dk cms mf (b) Modi ed spring-mass-damper model Figure 3: Linear model problem. 10−310−210−1100101102103 m∗0.00.20.40.60.81.0α Figure 4: Variation of with respect to m. 7 4.1. Governing equations Under the absence of external force, the governing equations for the original spring-mass-damper model are given by, md+c_d+kd= 0; (11) d(t= 0) =d0; (12) _d(t= 0) = _d0; (13) where,mis the mass, cis the damping coefﬁcient, kis the stiffness, dis the displacement, _d( =dd dt)is the velocity and d( =d2d dt2)is the acceleration, d0is the initial displacement, and _d0is the initial velocity. Other useful quantities, namely, the undamped natural frequency, !, the damping ratio, , and time period, Tare deﬁned as !:=r k m;:=c 2p km;T:=2 !: (14) The governing equations for the modiﬁed model problem can be written as, ds+!2ds=fs; (15a) (1 ) _vf+ 2!vf=ff; (15b) _ds=vf; (15c) fs+ff= 0; (15d) where,ds,_dsanddsare displacement, velocity and acceleration of the solid sub-system, respectively; vf(= _df),_vf(=df)are the velocity and acceleration of the ﬂuid sub-system, respectively. fsandffare the forces on the solid and ﬂuid sub-systems, respectively. The solid and ﬂuid sub-system satisfy the kinematic constraint (15c) and the force equilibrium (15d) at the ﬂuid-solid interface. Remark 1: The staggered scheme originally proposed in [42] is based on the generalised-alpha scheme proposed by Chung and Hulbert [77] for second-order ordinary differential equations (ODEs) for structural dynamics. The present work uses the generalised-alpha scheme originally proposed by Jansen et al. [78] for the ﬂuid problem (ﬁrst-order ODE) and then later extended to structural dynamics in Kadapa [79]. The decision to use the time integration scheme for the ﬁrst-order equations is for the ease of future extensions to ﬂuid-structure-control interaction problems. Nevertheless, the staggered scheme can be used with any other second-order accurate time integration scheme, for example, the Newmark-beta method [80] or the HHT-alpha method [81]. Now, introducing an auxiliary variable vs=_ds, the system of equations (15) for the modiﬁed model prob- lem, including the initial conditions, can be re-written as the following system of ﬁrst-order equations, _ds=vs; (16a) _vs+!2ds=fs; (16b) [1 ] _vf+ 2!vf=ff; (16c) vs=vf; (16d) fs+ff= 0; (16e) ds(t= 0) =ds 0; (16f) _ds(t= 0) =vs(t= 0) =vf(t= 0) = _ds 0: (16g) (16h) 8 4.2. Generalised-alpha time integration scheme for the rst-order ODE For the ﬁrst-order ODE, _u#u=f; (17) where,uis the solution variable, #is a constant and fis the external force, the governing equation for the generalised-alpha scheme of Jansen et al. [78] is given as, _un+ m#un+ f=fn+ f; (18) with, un+ f= [1 f]un+ fun+1; (19) fn+ f= [1 f]fn+ ffn+1; (20) _un+ m= [1 m] _un+ m_un+1; (21) un+1=un+ th [1 ] _un+ _un+1i ; (22) _un+1=1 th un+1uni + 1 _un: (23) Here, tis the time step size, and m, fand are the parameters of the time integration scheme. The generalised-alpha scheme is second-order accurate if the parameters are chosen such that m=1 231 1 +1; f=1 1 +1; =1 2+ m f;for011; (24) where,1is the spectral radius at inﬁnite time step size, and it controls the amount of numerical damping. For1= 0, this scheme is spectrally equivalent to the second-order accurate backward differentiation formula (BDF2). The interested reader is referred to Jansen et al. [78] for the comprehensive details on the spectral characteristics of the generalised-alpha time integration scheme, and to Kadapa et al. [79] for its extension to structural dynamics problems. 4.3. Monolithic scheme In a monolithic approach for FSI problems, the coupled equations (16) are all solved together for the solu- tion variables in the solid and ﬂuid sub-systems. Using the generalised-alpha scheme, the system of equations for the monolithic scheme for the solution of equations (16) can be written as, _ds n+ smvs n+ s f= 0; (25a) _vs n+ sm+!2ds n+ s f=fs n+ s f; (25b) [1 ] _vf n+ f m+ 2!vf n+ f f=ff n+ f f; (25c) vs n+1=vf n+1; (25d) fs n+1+ff n+1= 0: (25e) Withvs n=vf nandfs n+1=ff n+1=fn+1, the system of equations (25) can be written in the matrix-vector form as, 9 8 >>>>>>< >>>>>>:ds n+1 vs n+1t _ds n+1t _vs n+1t2 _vf n+1t2 fn+1t29 >>>>>>= >>>>>>;=2 6666664Amono3 77777758 >>>>>>< >>>>>>:ds n vs nt _ds nt _vs nt2 _vf nt2 fnt2;9 >>>>>>= >>>>>>;(26) where, Amono is the ampliﬁcation matrix for the monolithic scheme. 4.4. Staggered scheme Using the generalised-alpha scheme, the governing equations for solution of equations (16) with the stag- gered approach are given as fsP n+1= 2fs nfs n1; (27a) fs n+ s f= s ffsP n+1+ [1 s f]fs n; (27b) _ds n+ smvs n+ s f= 0; (27c) _vs n+ sm+!2ds n+ s f=fs n+ s f; (27d) [1 ] _vf n+ f m+ 2!vf n+ f f=ff n+ f f; (27e) ff n+ f f= f fff n+1+ [1 f f]ff n; (27f) fs n+1= ff n+1+ [1 ]fsP n+1; (27g) vs n+1=vf n+1; (27h) fs n+1+ff n+1= 0: (27i) The above system of equations (27) can be written in the matrix-vector form as, 8 >>>>>>>>>< >>>>>>>>>:ds n+1_ds n+1t vs n+1t _vs n+1t2 _vf n+1t2 fn+1t2 fnt29 >>>>>>>>>= >>>>>>>>>;=2 666666664Astag3 7777777758 >>>>>>>>>< >>>>>>>>>:ds n_ds nt vs nt _vs nt2 _vf nt2 fnt2 fn1t29 >>>>>>>>>= >>>>>>>>>;(28) where, Astagis the ampliﬁcation matrix for the staggered scheme. Using the ﬁeld variables for the model problem, the pseudocode for the staggered scheme is shown in Algorithm. 1. Algorithm 1 Key steps in the staggered scheme for the model problem 1:fort=t0tot=tfdo 2: Predict the force on the solid: fsP n+1= 2fs nfs n1 3: Solve the solid problem: fds n+1;vs n+1g=S(fsP n+1) 4: Solve the uid problem: fvf n+1;f n+1g=F(ds n+1;vs n+1) 5: Correct the force on the solid: fs n+1= f n+1+ [1 ]fsP 6:end for 10 Remark 2: The staggered scheme has also been experimented with third-order ( fsP n+1= 3fs n3fs n1+fs n2) and fourth-order ( fsP n+1= 4fs n6fs n1+ 4fs n2fs n3) force predictors. Although these higher-order force predictors improve the accuracy of the staggered scheme for problems with signiﬁcant added-mass, i.e., smaller values ofmor , it has been observed that their stability limits are much more restrictive (with respect to ) than the second-order predictor. Therefore, they are not pursued in this work. Remark 3: Our experience indicates that the instability condition for the partitioned FSI schemes based on the displacement predictors (for the solid problem) is much more restrictive with respect to mwhen compared with the one based on the force predictors used in the present work. Since there is no unique way of designing a staggered scheme using displacement predictors, the interested reader can compute the stability limits for the staggered scheme of their choice by using the linear model problem considered in this work. Remark 4: For adaptive time stepping, the force predictor given in Eq. (27a) is to be modiﬁed as fsP n+1= [1 + ]fs nfs n1;with  =tn+1 tn; (29) where, tn+1andtnare the time step sizes, respectively, for the current and previously converged time steps. 4.5. Stability and accuracy analysis Analytical proof of stability and accuracy of systems (26) and (28) is quite challenging. In this work, we assess the spectral characteristics of the monolithic and staggered schemes by evaluating the spectral radii of the corresponding ampliﬁcation matrices for several combinations of system parameters. The spectral radius of a matrix Ais deﬁned as the largest value of its eigenvalues [82, 83]. A system is said to be stable for the chosen parameter values if the spectral radius of its ampliﬁcation matrix is less than or equal to one; otherwise, the system is unstable. Plots of spectral radii for four different values of m, 10, 2, 0.2 and 0.02, are shown in Fig. 5 for = 0and = 0:1. The numerical solutions computed with the staggered scheme for the same mvalues are shown in Fig. 6. The important characteristics of the staggered scheme are: (i)The staggered scheme is second-order accurate in time . (ii)The leading error term for = 0ands 1=f 1= 0is obtained as, " O(t2) =4h 3 + 2 [1 +m]i 27m: (30) Although the accuracy of the staggered scheme deteriorates as mapproaches zero, accurate numerical results can be obtained for the values of mthat are encountered in real-world applications, as demon- strated with the numerical examples of galloping of a square and vortex-induced vibrations of a circular cylinder, see Section 5. Moreover, since the scheme is second-order accurate, the time step size required for obtaining accurate numerical results is not prohibitively small. (iii) For the unconditional stability of the staggered scheme for the undamped case, i.e., = 0, the condition on the relaxation parameter is given by 4 3 +1 ; fors 1=f 1=1: (31) This condition implies that smaller values of are required for smaller values of , which is the case for problems with signiﬁcant added-mass. The condition for unconditional stability becomes very com- plicated for the case with damping. Therefore, the relaxation parameter is chosen based on numerical experiments. (iv) The amount of numerical damping for the monolithic scheme remains the user-speciﬁed value, and is independent of the relaxation parameter , as shown in Fig. 5. However, for the staggered scheme, the amount of numerical damping in the overall algorithm is a complicated function of s 1,f 1and . Using the limit case analysis as t!1 , the relation for the numerical damping is obtained as, lim t!1= max s 1;f 1;p 1  : (32) 11 This means that, for small values of , the spectral radius approaches one, sincep1 1for 0:0. This is the case irrespective of the physical damping and the user-speciﬁed numerical damping, as shown in Fig. 5. 12 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5)(a)m= 10:0,= 0:0 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5) (b)m= 10:0,= 0:1 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5) (c)m= 2:0,= 0:0 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5) (d)m= 2:0,= 0:1 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5) (e)m= 0:2,= 0:0 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5) (f)m= 0:2,= 0:1 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5) (g)m= 0:02,= 0:0 −3 −2 −1 0 1 2 3 log(∆ t/T)0.00.20.40.60.81.01.2Spectral radiusMonolithic( ρ∞= 0.0) Staggered ( ρ∞= 0.0) Monolithic( ρ∞= 0.5) Staggered ( ρ∞= 0.5) (h)m= 0:02,= 0:1 Figure 5: Model problem: spectral radii for != 1:0,mf= 1:0,s 1=f 1=1and = . 13 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered(a)m= 10:0,= 0:0 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered (b)m= 10:0,= 0:02 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered (c)m= 2:0,= 0:0 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered (d)m= 2:0,= 0:02 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered (e)m= 0:2,= 0:0 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered (f)m= 0:2,= 0:02 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered (g)m= 0:02,= 0:0 0 10 20 30 40 50 Time−1.5−1.0−0.50.00.51.01.5DisplacementMonolithic Staggered (h)m= 0:02,= 0:02 Figure 6: Model problem: spectral radii for != 1:0,mf= 1:0,s 1=f 1= 0:0, = and t=T=200. 14 5. Numerical examples The suitability of the staggered scheme for applications in FSI problems, especially those involving lightweight rigid solids often encountered in the interaction of ﬂoating bodies in Ocean Engineering and marine structures [84–87], is demonstrated using the examples of galloping of a square, vortex-induced vibrations of a circular cylinder and an extremely lightweight rectangular block rising in a ﬂuid chamber. In this work, the ﬂow is assumed to be isothermal, laminar, viscous and incompressible. The proposed staggered scheme for resolving ﬂuid-solid interactions is based on the Dirichlet-Neumann coupling; therefore, the input to the ﬂuid solver is the position and velocity of the solid body and the input to the solid solver is the force exerted by the ﬂuid. Thus, the proposed scheme is non-intrusive, meaning that it can be used for combining different ﬂuid and solid solvers by treating them as black-box solvers. To demonstrate this non-intrusive nature of the staggered scheme, two different ﬂuid solvers are considered in this work. The ﬁrst ﬂuid solver is based on the ﬁctitious-domain/distributed Lagrangian multiplier methodology proposed in Kadapa et al. [27] and the second one is based on the cut-cell based approach proposed in Kadapa et al. [14, 15, 104]. Both ﬂuid solvers are based on the concept of unﬁtted meshes. The reader is referred to [11, 14, 15, 27, 104] for the comprehensive details of the formulation for the ﬂuid solvers. In all the numerical simulations, the spectral radius is considered as s 1=f 1= 0. The guidelines for the selection of the relaxation parameter are outlined in Section 6. 5.1. Galloping of a square Galloping is a self-excited instability which results in high-amplitude and low-frequency oscillations unlike lock-in which is caused by resonance [88, 89]. The geometry and boundary conditions of the problem are as shown in Fig. 7(a). The length of the side of the square is D= 1. The properties of the ﬂuid are: density, f= 1, and dynamic viscosity, f= 0:01. The value of v1at the inlet is adjusted to match the Reynolds number ( Re) of interest, where, Re=fDv1=f. The simulation framework used for this example is the ﬁctitious-domain/distributed Lagrange multiplier methodology [27]. The main reason behind using this framework is its ability to compute numerical solutions of FSI problems using a monolithic scheme which serve as the reference solutions for assessing the accuracy of the results obtained with the proposed staggered scheme. Before presenting the results for galloping, spatial and temporal convergence of the formulation are assessed by simulating the unsteady ﬂow past a ﬁxed square. 5.1.1. Spatial and temporal convergence study - ow past a xed square For the convergence study, unsteady ﬂow past a ﬁxed square for Re= 100 andRe= 200 is considered. The ﬁnite element mesh with multiple levels of hierarchical reﬁnement, as shown in Fig. 7, is used for the analysis. The number of points per side of the square is 10 and 20, respectively, for the level-2 (L2) and level-3 (L3) meshes. The number of DOFs, including Lagrange multipliers, for the ﬂuid problem is 32441 and 46867, respectively, for the level-2 and level-3 meshes. Two different time step sizes, t= 0:1andt= 0:05 are considered for each mesh. The graphs of the evolution of the lift coefﬁcient, as shown in Figs. 8 and 9, respectively, for Re= 100 andRe= 200 , indicate convergence of the numerical solutions with respect to spatial and temporal discretisation. The computed values of mean drag coefﬁcient ( CDmean), RMS lift coefﬁcient (CLrms) and Strouhal number ( St) are tabulated in Tables. 1 and 2, respectively, for Re= 100 andRe= 200 . As shown, the present numerical results match well with the corresponding values from the published literature [90–98]. 15 10 201515 1vf x=v∞ vf y= 0 vf y= 0 , tf x= 0vf y= 0 , tf x= 0 tfx=tfy= 0 XY(a) (b) Figure 7: Flow past a square: (a) geometry and boundary conditions and (b) the nite element mesh used for the simulation. 0 25 50 75 100 125 150 175 200 Time−0.4−0.20.00.20.4Lift coeﬃcientL2(∆ t= 0 .1) L2(∆ t= 0 .05) L3(∆ t= 0 .1) L3(∆ t= 0 .05) Figure 8: Flow past a xed square: evolution of the lift coecient ( CL) forRe= 100. 0 10 20 30 40 50 60 70 80 Time−1.0−0.50.00.51.0Lift coeﬃcientL2(∆ t= 0 .1) L2(∆ t= 0 .05) L3(∆ t= 0 .1) L3(∆ t= 0 .05) Figure 9: Flow past a xed square: evolution of the lift coecient ( CL) forRe= 200. 16 Study CDmeanCLrmsSt Sharma and Eswaran [90] 1.4936 0.1922 0.1488 Sen et al. [91] 1.5287 0.1928 0.1452 Sahu et al. [92] 1.4878 0.1880 0.1486 Sohankar et al. [93] 1.4770 0.1560 0.1460 Singh et al. [94] 1.5100 0.1600 0.1470 Darekar and Sherwin (3D) [95] 1.4860 0.1860 0.1460 Zhao et al. [96] 1.4520 0.1908 0.1447 Present (Level-2,  t= 0:1) 1.6136 0.1876 0.1440 Present (Level-2,  t= 0:05) 1.6136 0.1877 0.1420 Present (Level-3,  t= 0:1) 1.5523 0.1843 0.1460 Present (Level-3,  t= 0:05) 1.5525 0.1842 0.1460 Table 1: Flow past a xed square: mean drag coecient ( CDmean), RMS lift coecient ( CLrms) and Strouhal number (St) forRe= 100. Study CDmeanCLrmsSt Sohankar et al. [97] 1.4600 0.3200 0.1700 Singh et al. (3D) [94] 1.5400 0.5500 0.1420 Joly et al. [98] 1.4400 0.4200 0.1510 Present (Level-2,  t= 0:1) 1.7830 0.5519 0.1509 Present (Level-2,  t= 0:05) 1.7934 0.5658 0.1530 Present (Level-3,  t= 0:1) 1.6320 0.5376 0.1470 Present (Level-3,  t= 0:05) 1.6385 0.5501 0.1470 Table 2: Flow past a xed square: mean drag coecient ( CDmean), RMS lift coecient ( CLrms) and Strouhal number (St) forRe= 200. 5.1.2. Galloping with mass ratio m= 20 First, we consider the simulation case for m= 20 that has been originally proposed in Robertson et al. [99] and later studied in [7, 14, 28, 58, 100, 101]. For this case, the simulation parameters are: m= 20 ,Re= 250 , reduced velocity, Ur=v1 fnD= 40 , and= 0:0037 . The natural frequency for this system is, fn= 0:0625 Hz. Level-3 mesh with t= 0:05is used for performing simulations with both monolithic and staggered schemes. The relaxation parameter for the staggered scheme is = 0:9 ( ). The displacement response of the square body obtained with the staggered scheme, as shown in Fig. 10, is in excellent agreement with that of the monolithic scheme. The values of displacement amplitude ( Y=D ), frequency of rigid body oscillations ( fo), and the frequency of vortex-shedding ( fv) match well with the reference values, as shown in Table. 3. Thus, the numerical results obtained with this case establish the accuracy of the staggered scheme against a benchmark example. 17 0 50 100 150 200 250 300 350 400 450 500 Time−1.5−1.0−0.50.00.51.01.52.0Y-displacementMonolithic StaggeredFigure 10: Galloping of a square: displacement response for m= 20. For this case, = 0:9. Study Y=Dfo=fnfv=fn Robertson et al. [99] 1.150 0.938 - Dettmer and Peri c [28] 1.117 0.943 6.33 He et al. [56] 1.279 0.952 6.367 He et al. [57] 1.159 0.939 6.525 Kadapa et al. [14] 1.197 0.944 - Present - Monolithic 1.174 0.939 6.488 Present - Staggered 1.180 0.958 6.489 Table 3: Galloping of a square: displacement amplitude and frequency of oscillations for m= 20,Re= 250,Ur= 40 and= 0:0037. 5.1.3. Galloping with mass ratio 0:01m5 Simulations are performed for different values of mass ratio in the range 0:01m5by keeping the rest of the parameters constant. The displacement response of the square for m=f5;2;0:2;0:02g, are presented in Fig. 11. These graphs illustrate that the displacement responses obtained with the proposed scheme are in excellent agreement with the respective solutions obtained with the monolithic scheme. The time step size of t= 0:02used form=f2;0:2;0:02gcorresponds to only 100 time steps per one period of sustained oscil- lation. Successful simulations for even lower values of mcan also be performed by reducing the relaxation parameter and time step size appropriately, as demonstrated in Fig. 12 for m= 0:01, and also in the following examples. The presented results illustrate the exceptional accuracy of numerical results that can be obtained with the proposed scheme. It is remarkable that such accurate numerical results can be obtained with the proposed staggered scheme, especially considering that the time step size is large and that the scheme does not require any sub-iterations at each time step. The ability of the staggered scheme to compute accurate results using large time steps and without sub-iterations reaps signiﬁcant computational beneﬁts in performing large-scale 3D FSI simulations. The graphs in Fig. 11 show that the amplitude of periodic oscillations of the square in galloping reduces substantially for low values of m. This behaviour of reduced amplitude of oscillations of a galloping body for low values of mis consistent with the earlier observations reported in Joly et al. [98]. 18 0 40 80 120 160 200 240 Time−0.75−0.50−0.250.000.250.500.751.00Y-displacementMonolithic Staggered(a)m= 5 ( = 0:8;t= 0:05) 0 10 20 30 40 50 60 70 80 90 100 Time−0.75−0.50−0.250.000.250.500.75Y-displacementMonolithic Staggered (b)m= 2 ( = 0:5;t= 0:02) 0 10 20 30 40 50 60 70 80 90 100 Time−1.25−1.00−0.75−0.50−0.250.000.250.500.75Y-displacementMonolithic Staggered (c)m= 0:2 ( = 0:1;t= 0:02) 0 10 20 30 40 50 60 70 80 90 100 Time−1.25−1.00−0.75−0.50−0.250.000.250.500.75Y-displacementMonolithic Staggered (d)m= 0:02 ( = 0:01;t= 0:02) Figure 11: Galloping of a square: displacement response for four di erent values of m. The corresponding values of the relaxation parameter and time step size  tare indicated in the captions of sub gures. 19 0 10 20 30 40 50 Time−1.25−1.00−0.75−0.50−0.250.000.250.500.75Y-displacementMonolithic, ∆ t= 0 .02 Staggered, ∆ t= 0 .02 Staggered, ∆ t= 0 .01Figure 12: Galloping of a square: displacement response for m= 0:01. Relaxation parameter, = 0:005. 5.2. Vortex-induced vibrations of a circular cylinder This example is concerned with the problem of vortex-induced oscillations of a circular cylinder at Re= 100 . To demonstrate the non-intrusive nature of the staggered scheme, a different ﬂuid solver from the one used for the previous example is considered. The ﬂuid problem in this example is solved using the stabilised ﬁnite element formulation proposed in Kadapa et al. [11, 14, 15, 104]. The key difference between the formulation for the present example and previous example is in the technique used for imposing the kinematic constraint between the ﬂuid and solid interface. While interface condition is enforced using distributed Lagrange multi- pliers in the ﬁctitious domain method used in the previous example, it is enforced using the Nitsche method in the cut-cell-based method used in this example. For comprehensive details regarding the ﬁnite element for- mulations, we refer the reader to [11, 14, 15, 104]. The spatial and temporal convergences of the formulation for the example of ﬂow past a ﬁxed circular cylinder and the benchmark examples of transverse galloping of a square and rotational galloping a rectangle are already established in [14]. Therefore, convergence studies are not repeated in this work. The geometry and boundary conditions of the ﬂuid domain are shown in Fig. 13(a). Important parameters of the problem are: cylinder diameter, D= 1; ﬂuid density, f= 1; ﬂuid viscosity, f= 0:01; inlet velocity, v1= 1; natural frequency, fn= 0:2Hz; and damping ratio, = 0. For these parameters, the reduced velocity is, Ur=v1 fnD= 5. The displaced mass of the ﬂuid, mf, is kept constant, and the solid mass and spring stiffness are adjusted based on the mass ratio, m. To accommodate the signiﬁcant streamwise movement of the cylinder for low mass ratio values, the ﬁnite element mesh with extended reﬁnement in the streamwise direction, as shown in Fig. 13(b), is used. The total number of DOFs for the ﬂuid problem is approximately 125600 at each time step. The circumference of the cylinder is discretised with 80 straight edges. Simulations are performed for different values of min the range 0:01m10. The time step for all the simulations is t= 0:05, which corresponds to only 100 time steps per one natural time period of the spring-mass system. 20 10 201515 1vf x=v∞ vf y= 0 vf y= 0 , tf x= 0vf y= 0 , tf x= 0 tfx=tfy= 0 XY(a) (b) Figure 13: VIV of a circular cylinder: (a) geometry and boundary conditions and (b) the nite element mesh with three levels of hierarchical re nement. 5.2.1. 2-DOF motion with mass ratio m= 10 The relaxation parameter considered for this case is = 0:45 (0:5 ). The evolution of force coefﬁcients and displacements of the circular cylinder obtained with the staggered scheme is shown in Fig. 14. The amplitudes of Y-displacement and lift coefﬁcient obtained in the present work are in good agreement with the results presented in Jaiman et al. [68], see Table 4. 0 50 100 150 200 250 300 Time−1.2−0.60.00.61.21.82.43.03.6Force coeﬃcientCD CL 0 50 100 150 200 250 300 Time−0.6−0.4−0.20.00.20.40.6Displacementds x ds y Figure 14: VIV of a circular cylinder: (a) force coecients and (b) displacement response, for m= 10. The relaxation parameter is = 0:45. Study Displacement amplitude ( Y=D)Lift coecient ( CL) Jaiman et al. [68] 0.56 0.88 Present 0.59 0.83 Table 4: VIV of a circular cylinder: displacement amplitude ( Y=D) and lift coecient ( CL) form= 10. 5.2.2. 2-DOF motion with mass ratio 0:1m2 Analogous to the studies performed in Jaiman et al. [68] and He [69], a number of simulations are per- formed with different mass ratio values in the range 0:1m2for the vortex-induced vibrations of a circular cylinder with 2-DOF motion. Note that the lowest mvalue reported in Jaiman et al. [68] is 0.26, and in He 21 [69] it is 0.298, while the lowest value in the present work for the circular cylinder with 2-DOF motion is 0.1. Simulations with even lower values for mare performed using the 1-DOF model in the next subsection. The time step is t= 0:01fort5andt= 0:05for the rest of the simulation. The Lissajous orbits of the centre of mass of the cylinder for different values of m, as shown in Fig. 15, demonstrate the sustained periodic oscillations resulting in the well-known 8-proﬁle for the vortex-induced vibrations of a circular cylinder at low Reynolds numbers [68, 69, 76, 102, 103]. The evolution of force coefﬁcients and displacements of the cylinder for m= 0:408, as presented in Fig. 16, agree well with the results presented in He [69]. The contour plots of pressure at two different time instants are shown in Fig. 17. To the best of the author’s knowledge, this is the ﬁrst time in the literature vortex-induced vibrations of a circular with 2-DOF motion are reported for such a low value of 0.1 for the mass ratio. 1 2 3 4 5 6 7 8 X-displacement−0.8−0.6−0.4−0.20.00.20.40.60.8Y-displacementm∗=0.5 m∗=0.408m∗=0.4m∗=0.3 m∗=0.2 m∗=0.1 Figure 15: VIV of a circular cylinder: trajectories of the centre of mass of the cylinder for di erent mvalues. 0 50 100 150 200 250 300 Time−0.50.00.51.01.52.02.5Force coeﬃcientCD CL 0 50 100 150 200 250 300 Time−1.0−0.50.00.51.01.52.02.53.0Displacementds x ds y Figure 16: VIV of a circular cylinder: (a) force coecients and (b) displacement response, for m= 0:408. The relaxation parameter is = 0:14. 22 -0.36 -0.22 -0.08 0.06 -0.50 0.20pressure -0.36 -0.22 -0.08 0.06 -0.50 0.20pressureFigure 17: VIV of a circular cylinder: instantaneous contour plots of pressure at two di erent time instants for m= 0:408. 5.2.3. 1-DOF motion with mass ratio 0:01m<0:1 The streamwise displacement of the cylinder becomes excessive for m<0:1, which requires an exten- sion of the ﬂuid domain further downstream. To avoid the associated computational cost, cylinder with only transverse DOF is considered for 0:01m<0:1. Simulations are performed for different values of min the range 0:01m<0:1using a uniform time step of t= 0:05and the responses of the cylinder for four different values, m=f0:04;0:03;0:02;0:01g, are presented in Fig. 18. The respective relaxation parameters used for these simulations are 0.02, 0.015, 0.01 and 0.005, which correspond to 0:5 . As shown in Fig. 18, the response of the cylinder is smooth, and there is no signiﬁcant difference in the amplitude and frequency of oscillations. The ever so slight increase in the time period of oscillations with decreasing mis consistent with the expected behaviour: as the mass ratio decreases the added-mass increases, thereby increasing the effective mass. This increased effective mass reﬂects in the increased time period (or reduced frequency) of oscillations. The results presented for 0:01m<0:1illustrate the remarkable ability of the proposed iteration-free staggered scheme in performing successful simulations of FSI problems involving extremely lightweight rigid solids. 23 0 25 50 75 100 125 150 175 200 Time−0.50−0.250.000.250.500.75Y-displacementm∗= 0.04 m∗= 0.03m∗= 0.02 m∗= 0.01(a) 20 30 40 50 60 70 80 90 100 Time−0.50−0.250.000.250.500.75Y-displacementm∗= 0.04 m∗= 0.03m∗= 0.02 m∗= 0.01 (b) Figure 18: VIV of a circular cylinder: displacement response of the cylinder for four di erent values of m. (a) time range [0;200] and (b) time range [20 ;100]. 5.3. Rectangular body rising to the top of a uid chamber As the last example, we consider the problem of rising of a lightweight rectangular block with rounded corners that was originally studied in Banks et al. [64]. The geometry and boundary conditions of the problem are shown schematically in Fig. 19. The density of the ﬂuid is f= 1, and its dynamic viscosity is f= 0:025. The density of the solid block is b= 0:001, which is 1=1000thof that of the ﬂuid. Due to the differences in the densities, the solid rises in the ﬂuid under the inﬂuence of gravity. Following [64], the buoyancy on the solid is given as, Fbuoyancy =Ab[bf]R(t); (33) where,Abis the area of the solid and R(t)is a ramp function, given by R(t) =8 < :[35 + [84 + [7020t]t]t]t4if 0t1; 1; ift>1:(34) 24 For this example, the mass-ratio parameter is m= 0:001. Such a signiﬁcantly low value of mposes serious difﬁculties in simulating this example using the conventional staggered scheme. In fact, the conventional staggered scheme without any relaxation parameters fails to simulate this particular problem. Banks et al. [64] report that the traditional partitioning algorithm for rigid bodies (TP-RB), see [63], requires an average of 85 sub-time-step iterations for each time step on the coarsest grid used in [64]. Banks et al. [64] successfully simulated this problem using their added-mass partition algorithm for rigid bodies (AMP-RB). For the simulations, we consider two uniform meshes, as shown in Fig. 20. With each mesh, we perform the simulations with three different time steps, t=f0:02;0:01;0:005g. The ﬂuid solver is the same as the one used for the circular cylinder example. The relaxation parameter for all the simulations is = 0:0002 . The motion of the square body during its rise is plotted in terms of displacement and velocity response in Fig. 21 for all the six simulation cases, along with the reference solution from Banks et al. [64]. As shown in Fig. 21, there is no noticeable difference in the numerical results obtained with different meshes, and more importantly, the present results converge to the reference solution as the time step is reduced. Except for the initial spikes, the displacement response of the solid body is indistinguishable from the reference solution for t3. Typical contour plots of ﬂuid pressure and plots of streamlines for time instants t= 2,t= 4andt= 10 are presented in Fig. 22. The present example demonstrates the excellent capability of this seemingly simple staggered scheme in computing accurate numerical solutions of very challenging FSI problems involving lightweight rigid solids. We reiterate that the ﬂuid and solid solvers do not require any additional information other than the standard Dirichlet and Neumann boundary conditions, respectively, for the ﬂuid and solid sub-problems, from their counterparts. 0.75 0.25 0.25 0.75 0.5 0.5 0.5 0.5XYvf=0 vf=0vf=0 vf=0Ωf Ωsr=0.1 Figure 19: Rising rectangle: geometry and boundary conditions. 25 (a) Mesh 1 (5151 elements) (b) Mesh 2 (101101 elements) Figure 20: Rising rectangle: nite element meshes with the original position of the rectangle. 0 2 4 6 8 10 12 Time−0.10.00.10.20.30.40.50.60.7Y-displacementReference Mesh 1, ∆ t= 0 .02 Mesh 1, ∆ t= 0 .01 Mesh 1, ∆ t= 0 .005 Mesh 2, ∆ t= 0 .02 Mesh 2, ∆ t= 0 .01 Mesh 2, ∆ t= 0 .005 (a) Displacement response 0 2 4 6 8 10 12 Time−0.10.00.10.20.3Y-velocityReference Mesh 1, ∆ t= 0 .02 Mesh 1, ∆ t= 0 .01 Mesh 1, ∆ t= 0 .005 Mesh 2, ∆ t= 0 .02 Mesh 2, ∆ t= 0 .01 Mesh 2, ∆ t= 0 .005 (b) Velocity response Figure 21: Rising rectangle: displacement and velocity response of the solid body with M1 and M2 meshes and di erent time steps. The reference solution is from Banks et al. [64]. 26 0.09 0.18 0.27 0.36 -0.01 0.45pressure(a)t= 2:0 0.040 0.080 0.120 0.160 0.000 0.200abs(vel) (b)t= 2:0 -0.83 -0.71 -0.60 -0.48 -0.94 -0.35pressure (c)t= 4:0 0.060 0.120 0.180 0.240 0.000 0.300abs(vel) (d)t= 4:0 -0.15 0.00 0.15 0.31 -0.34 0.42pressure (e)t= 10:0 0.006 0.012 0.018 0.024 0.000 0.030abs(vel) (f)t= 10:0 Figure 22: Rising rectangle: contour plots of pressure and streamlines at three di erent time instants obtained with Mesh 2 and  t= 0:01. 27 6. Selection criterion for the relaxation parameter ( ) Almost all of the staggered schemes rely on relaxation parameter(s) for their success. Such a parameter in the proposed scheme that needs to be supplied by the user is the relaxation parameter . The minor in- convenience in selecting a suitable value for is negligible when compared with the remarkable accuracy of the computed numerical results as well as the computational efﬁciency of the staggered scheme. Moreover, the selection of is not completely arbitrary, as discussed below. We offer the following insights to help in choosing the value of for the given value of m. (i)Considering the stability condition (31) obtained from the analysis of the linear model problem in Section 4.5 into account, the valid range of for any non-zero value of mbecomes, 0<  ; (35) for the worst-case scenario of zero numerical damping, i.e. s 1=f 1= 1. While this stability condition might vary due to the spatial and temporal discretisation schemes for the ﬂuid and solid solvers, the inclusion of physical damping and the manner in which boundary conditions are applied, our experience indicates that the variation is not by orders of magnitude for a particular value of m. Therefore, the above condition (35) serves as the starting point for choosing an appropriate value of for the given m. (ii)The task of selecting is also eased by the fact that the simulation crashes quickly if the chosen value of is larger than that dictated by the stability limit; the farther the value of away from the stability limit the quicker the simulation crashes . It has been observed that the simulation crashes within the ﬁrst 10-100 time steps, as illustrated in Fig. 23 for the example of galloping of a square for m= 0:2for different values of . Thus, a value for that works can be selected relatively quickly, considering especially that the present staggered scheme is free from sub-iterations at every time step. 0 20 40 60 80 100 120 Time step number−1.0−0.50.00.51.0Y-displacementβ= 1.25α β= 1.00αβ= 0.75α β= 0.50α Figure 23: Sensitivity of the simulation crash with respect to for the example of galloping of square for m= 0:2. 28 120 140 160 180 200 220 240 260 Time step number−0.75−0.50−0.250.000.250.500.751.00Y-displacementβ= 0.9(≈α) β= 0.45(≈0.5α)Figure 24: Displacement response of circular cylinder for m= 10 for di erent values of . Based on our extensive experience in performing successful FSI simulations using the staggered scheme in the present work as well as the previous works [14, 15, 104], the recommended value of form0:01is = 0:5 . The above recommendation is demonstrated to work with the example of galloping of a square for m= 0:2and m= 0:02, and VIV of a circular cylinder for which m0:01. Further reduction of was found necessary only whenmis extremely small, as seen in the example of rising of the rectangular block. Since such extreme cases are uncommon in real-world applications, = 0:5 is the recommended choice. This choice is further justiﬁed by the fact that it has a negligible effect on the accuracy of the results for the cases with large values of m, as demonstrated in Fig. 24 using the example of VIV of circular cylinder for m= 10 with two different values of . Note: To further aid the curious reader in the numerical experiments, the following procedure is suggested for choosing a value of that works: Start with = . If the simulation crashes, then reduce by 10% of and repeat the simulation. With decrements of 10%, the value of at the sixth in- crement is reduced to 0:5 , which is the recommended value. 7. Summary and Conclusions We have presented a second-order accurate staggered scheme for the simulation of ﬂuid-structure interac- tion problems involving ultra-lightweight rigid solids. The novelty of the present work lies in the demonstrated abilities of the proposed scheme in performing accurate numerical simulations of ﬂuid-rigid body interaction problems with signiﬁcant added-mass using large time steps and without sub-iterations. The capability of the staggered scheme to compute accurate solutions of FSI problems with signiﬁcant added-mass is demonstrated with the example of transverse galloping of a square for 0:01m20. The results obtained with the staggered scheme are shown to match remarkably well with those computed with the monolithic scheme. Later, the staggered scheme is applied to study the vortex-induced vibrations of a circular cylinder with two translational DOFs for 0:1m10and with only transverse DOF for 0:01m0:1. Finally, the excellent accuracy of the staggered is illustrated with an even more challenging example of a rectangular block with m= 0:001rising in a liquid chamber. The recommended value of the relaxation parameter is = 0:5 form0:01. Guidelines for choosing the relaxation parameter ( ) are also provided. The amount of effort required in selecting a value of that works is negligible when compared with the computational beneﬁts that can be realised with the proposed scheme. In addition to its excellent accuracy when using large time steps, especially for problems with signiﬁcant added-mass, the standout feature of the proposed staggered scheme is its non-intrusive nature. The scheme does not require any additional information from the ﬂuid and solid solvers except the interface velocity from 29 the solid solver and the interface traction from the ﬂuid solver. Therefore, it can be adapted to develop an FSI framework by combining different ﬂuid and solid solvers in a non-intrusive manner by treating the solvers as black boxes. This aspect of the staggered scheme is illustrated using two different CFD solvers in this work. The ongoing work is focussed at the extension of the proposed scheme for 3D problems as well as for problems with high Reynolds number ﬂows. Towards this, wrappers for coupling ANSYS and OpenFOAM CFD solvers with an in-house rigid-body dynamics solver are being explored. ACKNOWLEDGEMENTS The author acknowledges the support of the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government.