Branched structures are present in a diverse set of problems, from modelling branch pipe connections to simulating tree dynamics. Soft corals like the Bipinnate sea plume, have a branched geometry and are soft enough to bend under the waves. Due to their cylindrical cross section, a vortex street forms in the coral’s wake inducing vibrations of its branches. Despite extensive studies on VIV in straight geometries, the three-dimensional (3D) dynamics of flexible branched structures remains uninvestigated. In this numerical and experimental study, we develop a novel formulation for the accurate computation of in-line and cross-flow VIV of frame structures undergoing large defor- mation. The finite element approach is used to model arbitrarily complex geometries of branched frame structures. Our formulation allows to model complex geometries with forks or sharp angles. The consistent 3D corotational for- mulation for frame elements computes the internal, inertial and hydrodynamic forces. A wake-oscillator approach models the near wake dynamics with fluctuating fluid forces on the structure in the in-line and cross-flow directions. The drag and lift coefficients follow distributed Van der Pol oscillators. Moreover, we implement the numerical res- olution procedure in the open-source library ONSAS. The present formulation and numerical resolution procedure is validated by solving three examples, including comparisons with an analytical solution, a wake-oscillator, and experimental data from the literature. We also conduct experiments of a flexible and elastic cylinder clamped inside a water tunnel under a constant uniform flow. Amplitudes and power spectral density of the tip transverse displace- ments are compared with the model prediction. Finally, the proposed formulation is applied on a cylinder with two branches. The simulations demonstrate a multi-frequency response with higher amplitudes of displacements when additional branches are incorporated onto the cylinder, emphasizing the significance of considering VIV in nature and engineering applications for such geometries. Keywords: Vortex-induced vibrations, Wake-oscillator, Corotational framework, Biomechanics Preprint submitted to Journal of Fluids and Structures June 30, 2023 1. Introduction 1 Vortex-Induced Vibrations (VIV) are widely studied in engineering to mitigate offshore structure oscillations 2 [7, 50, 20, 47], for vibrations control [35, 13] and more recently as a potential source of clean energy [5, 71, 23]. 3 VIV are characterized by a coupling between the vibration s of a structure and the flow fluctuations in its wake. 4 An asymmetric vortex shedding exerts oscillating drag and lift forces, inducing in-line and cross-flow vibrations 5 [57, 34]. 6 Marine life has developed mechanisms through evolution to cope with VIV. These adaptations are observed 7 across various organisms and serve different purposes. For example, some plants benefit from flow-induced vibra- 8 tions, which enhance the dispersion rate of underwater plant canopies [58] or facilitate pollination [45]. A seal’s 9 whiskers structural design mitigates flow-induced vibrations to better chase preys in the dark [29, 4, 43]. On the 10 other hand, some animals use VIV to increase pheromone capture for moth males [61] or heat exchange for coral 11 tentacles [44]. In the context of soft corals, Gosselin first reported VIV on the flexible branches of Antillogorgia bipin- 12 nata. Two-dimensional (2D) simulations of a vibrating cross section conducted in Ref.[10] demonstrate an increase 13 in particle capture of up to 40 % more than a fixed branch during lock-in, suggesting that sea plumes may use VIV 14 as a feeding strategy. 15 Incorporating VIV into biomechanical models of plant is necessary to improve the accuracy and realism of 16 computational simulations. Numerous studies have highlighted the importance of accounting for flexibility when 17 modeling the behavior of plants submitted to hydrodynamic forces [42, 39, 63]. Experimental studies [11, 65] and 18 numerical simulations [70, 38] show that continuous flexible structures encounter multiple lock-ins due to the mul- 19 tiplicity of their natural frequencies. The intricate interplay between nonlinear wake dynamics and multi-modal 20 response presents significant challenges for modeling VIV in flexible structures and requires special numerical tools. 21 One approach to study VIV numerically is to solve a fully coupled system with the Navier Stokes equations 22 for the fluid and the elastodynamics equations for the structure. Several numerical methods in Computational 23 Fluid Dynamics (CFD) have been used to predict VIV, such as Quasi-three-dimensional Method [69], coupled three- 24 dimensional (3D) Finite Element Method (FEM) [31] or Finite Volume Method [46]. However, simulating all of the 25 surrounding fluid flow is extremely computationally demanding due to the vast range of length and time scales and 26 to the large amplitudes of motion that must be considered. 27 One alternative method is to empirically model the forcing of the wake on the structure. Bishop and Hassan and 28 Hartlen and Currie proposed the wake-oscillator model (WOM) to simulate the hydrodynamic forcing of the wake 29 on a rigid cylinder. Facchinetti et al. set a phenomenological coupling between the cross section acceleration and 30 the wake variable. Violette et al. showed that coupling the dynamics of a flexible cylinder with the nonlinear Van 31 der Pol oscillator exhibits good agreements with Direct numerical simulation (DNS) results. 32 Aquatic plants and soft corals are flexible and interact with the flow developing large deflections [27, 16, 17]. 33 Jain and Modarres-Sadeghi experimentally observed a vortex shedding parallel to the cylinder’s axis for inclination 34 2 angles up to 65 °. They validate the independence principle, stating that the inclined cylinders can be treated as the 35 normal-incidence ones by considering only the component of the free stream velocity normal to the cylinder axis. 36 With a finite element approach, the independence principle allows for the straightforward extension of the WOM 37 to inclined structures. Prior works implementing a FEM formulation of the WOM [38, 68, 40] as well as commercial 38 codes [41] have consistently defined the wake variables as additional degrees of freedom assigned to each node. 39 However, these studies primarily focused on analyzing structures with straight geometries, neglecting the presence 40 of forks, angles, or other complex geometric features. Consequently, the issue of having discontinuous tangent 41 vectors at specific nodes, which can result in ill-defined cross-flow directions, has not been addressed in these 42 prior investigations. Numerous wind-tree numerical models have been developed [59, 32, 55, 56, 18] and various 43 simulations of tree dynamics have been performed [1, 28, 53]. However, these studies often overlook VIV. The 44 assessment of VIV on 3D branched structures like corals or trees undergoing large displacements remains an open 45 problem. The motivation of this study is to develop a novel methodology for incorporating VIV in 3D simulations 46 of branched structures, using the WOM. 47 Structures undergoing VIV are typically slender and can exhibit large amplitude beam-like deformations. The 48 kinematics of a flexible frame element can be described by the total Lagrangian or the updated Lagrangian formu- 49 lations, however they require a high number of elements or a small time step in order to provide proper results 50 for nonlinear problems [19, 37]. Crisfield et al. presented a 3D corotational formulation in which the movement 51 is separated into two parts using an element-independent framework. The first part involves rigid body motion, 52 which includes translation and rotation from the global to the rigid system. The second part involves a small defor- 53 mation motion, which goes from the rigid to the deformed configuration. This formulation handles large rotations 54 and large displacements, obtaining accurate results with a smaller number of elements [37]. The corotational model 55 has been gradually extended to consider warping effects [3], dynamic problems [36], a consistent formulation for 56 frame elements [37] and hydrodynamic forces [66]. Corotational frame elements have been used in fluid-structure 57 interaction simulations to model living organisms subjected to fluid loads, such as insect wings during flapping [14], 58 or vibrating vocal folds [64]. 59 Here we present a unified corotational formulation for simulating branched frame structures undergoing large 60 displacements and rotations, coupled with a the wake-oscillator model to assess VIV. The consistent 3D corotational 61 frame element formulation computes the internal, inertial and hydrodynamic forces. Moreover, we implement the 62 formulation in the open-source structural analysis tool ONSAS1[51] and experiments are performed for validation. 63 As a result, a new free tool for VIV analysis of 3D forked structures is available. 64 The paper is structured as follows: in Section 2 we present the formulation of the in-line/cross-flow wake- 65 oscillator model in the 3D corotational framework together with the experimental setup. In Section 3 the in-line and 66 cross-flow VIV model is verified with one frame element (Example 1). The model is validated in large displacements 67 1www.onsas.org 3 with the wake-oscillator model from Leclercq and de Langre (Example 2) and with experimental data from Trim 68 et al. (Example 3). It is also compared with new experimental data in Example 4. Simulation results of VIV on a 69 forked structure are presented in Example 5. The paper ends with the conclusions in Section 4. 70 2. Methodology 71 In this section, we introduce our proposed corotational formulation for the WOM and added mass, along with 72 a numerical procedure for solving the governing equations. Additionally, we provide a detailed description of the 73 experimental setup used to validate our model. 74 The formulation is developed for a frame element with a length l0and a circular, uniform cross section of diameter 75 Dand second moment of inertia I. The material is considered to be elastic and isotropic, of density ρand elasticity 76 modulus E. In this section, extensions of the corotational formulation are developed and a numerical resolution 77 procedure is presented. 78 1. Corotational framework 79 We summarize here the key steps of the corotational formulation from [37]. The 3D corotational framework uses 80 three coordinate systems. The global reference system is defined by the triad ( ex,ey,ez). The local system attached 81 to the element after a rigid rotation without deformation is ( r1,r2,r3). The triad that is rigidly attached to each node 82 is (t1,t2,t3). These coordinates systems are illustrated in Fig. 1. 83 Let us consider a two-node frame element, where the nodal displacements are represented by the column vector 84 d= (u1)T,(w1)T,(u2)T,(w2)TT. Here, uiandwidenote the column vectors of linear displacements and rota- 85 tions, respectively, at node i. These displacements can be expressed in both the global coordinate system as dgand 86 the local coordinate system as dl. Similarly, a vector of nodal forces fcan also be represented in the corresponding 87 coordinate systems. 88 The internal force fintis derived based on the work of Ref.[3]. This forces vector is expressed in the local 89 coordinates system where the Euler-Bernoulli beam theory is assumed. Translations and rotations at the cross 90 section centroid are interpolated from the local nodal values, using linear shape functions for the axial displacement 91 and rotation, and cubic shape functions for transverse displacements. 92 The internal virtual work is the same in local and global coordinates, therefore: 93 δdT gfint,g=δdT lfint,l, (1) withδdgandδdlbeing the virtual displacements in global and local coordinates, respectively. We can express 94 δdl=Bδdg, (2) 4 t3 t1t2 t3 t1t2 r3 r1r2 ez exey˙ue 2˙ue 1Π23 UUUrel −˙u Upr tU ϕ Π23 Ucos(ϕ)t2t3 fl fd Upr tUtdtl Figure 1: 3D an planar illustration of corotational kinematics and hydrodynamics forces. 5 withBbeing a transformation matrix expressing the axial displacements and rotations from the local to the global 95 coordinate system. It follows 96 δdT gfint,g=δdT gBTfint,l. (3) This is valid for any virtual displacement δdg, so: 97 fint,g=BTfint,l, (4) where fint,l= faℓ(m1 ℓ)T(m2 ℓ)T is the vector of internal forces in local coordinates. It considers the axial force 98 faℓalong t1, bending moments and torsional moment mi ℓat node ialong t1andt2,t3, respectively. They are given 99 by a linear constitutive behavior. 100 The inertial force is obtained through the kinetic energy variations. The kinetic energy Kof the element is: 101 K=1 2Z l0ρ˙uTA˙u+ρ˙wTI˙wdl, (5) withu,wbeing the displacements and rotations vectors of the center of the cross section, Athe cross section area 102 andIthe geometric inertia tensor. The dot notation stands for time differentiation. Using the derivative chain rule, 103 the variations of Eq. (5) writes: 104 δK=−Z l0δuTρA¨u+δwT[ρI¨w+e˙wρI¨w]dl, (6) withe˙wbeing the skew operator associated to the vector ˙w. The inertial force vector of the element in global 105 coordinates is: 106 δK= (fine)Tδdg, (7) where 107 fine=Z l0n HT 1RT rρA¨u+HT 2RT r[ρI¨w+e˙wρI˙w]o dl. (8) The rigid rotation matrix Rrexpresses ( r1,r2,r3) in the global coordinate system. Matrices H1,H2depend on the 108 interpolation functions which are linear for axial displacements and hermitian for bending. Their exact expression 109 can be found in Ref.[37]. 110 2. Hydrodynamic forces 111 The structure is submitted to a uniform and steady flow of velocity Uand density ρf. The hydrodynamic forces 112 acting on the structure are the added mass, the steady drag, and the fluctuating drag and lift forces due to the shed 113 vortices. We define the instantaneous relative velocity Urelthat the cylinder perceives based on Uand the cylinder 114 velocity ˙u:Urel=U−˙u. In this formulation, the computation of the fluid forces is performed at every instant in 115 6 each section plane by considering only the normal component of the flow velocity, as stated by the independence 116 principle. Thus, the relative velocity is projected on the cross section plane Π23as shown in Fig. 1. This projection 117 is denoted Upr. We define the vector tUas the normalized projection of Uon the cross section plane in local 118 coordinates: 119 tU=L2 RrRTU/∥L2 RrRTU∥, (9) with RrRTbeing the rotation matrix used to transform a global vector into local coordinates. The operator L2 120 projects it onto the section plane Π23.tUprovides the flow direction seen by the cross section. The fluid velocity in 121 the cross section plane is U−U·t1=Ucos(ϕ)tU, with ϕbeing the angle between the fluid velocity UandΠ23. A 122 3D and planar representation of these vectors and angles is provided in Fig. 1. The added mass force per unit length 123 acting on the section is: 124 fa(s, t) =ma ˙U(s, t)−¨u(s, t) , (10) withma=CaρfDπ/4being the fluid added mass. The added mass coefficient Cais assumed to be 1for a circular 125 cross section [8]. Here, the flow is considered as steady: ˙U= 0. The corotational formulation of the added mass 126 force acting on a beam element of length l0is expressed as an inertial force of the same element having a density 127 Caρf: 128 fam=−Z l0n HT 1RT rCaρfA¨u+HT 2RT r[CaρfI¨w+e˙wCaρfI˙w]o dl. (11) At time tand curvilinear coordinate s, the instantaneous drag and lift forces per unit length are: 129 fdrag(s, t) =1 2CD(s, t)ρfDU2 pr(s, t)td(s, t), (12) flift(s, t) =1 2CL(s, t)ρfDU2 pr(s, t)tl(s, t), (13) withtd,tlbeing the drag and lift directions, respectively. These forces are derived from the empirical WOM that 130 models the fluctuating pressure field on the cylinder caused by vortices in its wake by oscillating drag and lift coef- 131 ficients CDandCL. In particular, the in-line vibrations are modeled by splitting the instantaneous drag coefficient 132 CD=C0 D+CDiinto a constant mean coefficient C0 Dand a fluctuating component induced by the wake C0 Di133 following [26, 68, 40, 25]. The instantaneous vortex-induced drag and lift coefficients are modeled by 134 CDi(s, t) =1 2C0 Dip(s, t), (14) CL(s, t) =1 2C0 Lq(s, t), (15) withC0 Di, C0 Lbeing the constant amplitudes of the fluctuating drag and lift coefficients on a stationary cylinder. The 135 variables pandqare the dimenionless in-line and cross-flow wake variables. They can be seen as the normalized drag 136 7 and lift coefficients. They both follow nonlinear Van der Pol oscillators and are coupled to the section acceleration: 137 ¨p(s, t) + 2ϵpΩf p2(s, t)−1 ˙p(s, t) + 4Ω2 fp(s, t) =Ap D¨up(s, t), (16) ¨q(s, t) +ϵqΩf q2(s, t)−1 ˙q(s, t) + Ω2 fq(s, t) =Aq D¨uq(s, t), (17) where Ωf= 2πStUcosϕ/D is the wake angular frequency, Ap, Aqare the coupling coefficients, ϵp, ϵqare the fluid 138 damping parameters, Stis the Strouhal number and ¨up,¨uqare the in-line and cross-flow accelerations, respectively. 139 For any cross-section in a frame element, ¨up,¨uqare defined as: 140 ¨up=¨u·tU, (18) ¨uq=¨u·(t1×tU). (19) At a bifurcation or branch node, both t1andtUare discontinuous. Therefore, ¨upand¨uqat this node are ill- 141 defined. To avoid this issue, we set the wake variables as uniform on each element where t1andtUremain continu- 142 ous. Therefore, the Van der Pol oscillator is solved once per element. This novel approach enables the application of 143 the WOM to interconnected frame elements, accommodating arbitrary connections such as forks and sharp angles. 144 For both in-line and cross-flow oscillators, the coupling is done by averaging the two nodal accelerations of the 145 element: 146 ¨p+ 2ϵpΩf p2−1 ˙p+ 4Ω2 fp=Ap D ¨u1 p+ ¨u2 p /2, (20) ¨q+ϵqΩf q2−1 ˙q+ Ω2 fq=Aq D ¨u1 q+ ¨u2 q /2, (21) where the subscripts ·pand·qdenote the in-line and cross-flow directions, respectively, and the superscripts ·1and 147 ·2refer to the element nodes. The cross product t1×tUprovides the cross-flow direction. ¨u·tUand¨u·(t1×tU)are 148 the accelerations of the cross section in the in-line and cross-flow directions. Eqs (16,17) are known to provide quasi- 149 harmonic oscillations of normalized amplitude P0= 2,Q0= 2 forΩf= 1. As the WOM is a phenomenological 150 model, these coefficients St, Ap, Aq, ϵp, ϵqare calibrated with experimental data. 151 The corotational formulation of the drag and lift force acting on one element is taken from Ref.[66]: 152 fd=1 2ρfDEZ l0n HT 1R∥L2 RrRTUrel∥CD(p)L2 RrRTUrelo dl, (22) fl=1 2ρfDEZ l0n HT 1R∥L2 RrRTUrel∥CL(q)L3L2 RrRTUrelo dl. (23) where the operator L3=exp([π/2,0,0]T)rotates the relative fluid velocity by π/2radians around the axis t1to 153 provide the local lift direction. Finally, RandEtransform the velocity from local to rigid and from rigid to global 154 8 coordinates, respectively [3, 66]. 155 Note that the fluid damping is included in the drag force. Indeed, as pointed out by Ogink and Metrikine, 156 a linearized form of the instantaneous drag force can lead to a fluid damping term proportional to the section 157 velocity, named stall term by Skop and Balasubramanian, Facchinetti et al.. In the 2D VIV, this hypothesis may not 158 be respected. Here, we consider the full nonlinear drag force following [48, 52]. 159 When attempting to solve for the dynamics of a branched structure, a system of differential equations needs to 160 be solved for each branch, and continuity conditions must be imposed at the bifurcation nodes, which results in a 161 significant increase in the complexity of the problem. In the proposed formulation, the FEM assembly approach is 162 considered. Therefore, no complexity is added by including a fork node in the structure. 163 3. Balance equations and numerical resolution procedure 164 The governing equations are derived by applying the principle of virtual work to all elements of the structure for 165 the forces described in Eqs.(8), (11), (23), and (22), while simultaneously solving the Van der Pol equations, Eqs.(20) 166 and (21), for each element. 167 We denote the structure assembled vectors with subscripts ·sand elements variables with superscripts ·e. The 168 solid and wake residuals are: 169   ru=fint,s(ds) +fine,s(ds,˙ds,¨ds)−fam,s(¨ds)−fd,s(ds,˙ds,ps)−fl,s(ds,˙ds,qs), re p= ¨pe+ 2ϵpΩf pe2−1 ˙pe+ 4Ω2 fpe−Ap D ¨ue,1 p+ ¨ue,2 p /2, re q= ¨qe+ϵpΩf qe2−1 ˙qe+ Ω2 fqe−Aq D ¨ue,1 q+ ¨ue,2 q /2,e= 1, . . . , N e(24) where ruis the residual of the Principle of Virtual Work equations, and the wake residuals re p, re qcorrespond to the 170 Van der Pol equations. They are both computed independently for each element. For each time step, the system 171 ru=0;rp=rq= 0 must be solved, obtaining the current configuration of the structure (with velocities and 172 accelerations) and wake variables for each element. 173 We define the following non-dimensional parameters: the reduced velocity, the Reynolds number, the aspect 174 ratio and the mass ratio, respectively, 175 UR=U Dℓ2r ms+ma EI, R e=UD νf,Γ =ℓ D,M=ρ ρf, (25) withmsbeing the structural mass per unit length and νfthe fluid kinematic viscosity. 176 The system of nonlinear governing Eq. (24) is solved numerically using iterative methods. The internal forces 177 are implemented according to the formulations described in [37, 3]. The inertial forces in Eq. (8) and hydrodynamic 178 forces in Eqs (11, 23, 22) are computed with 4 Gauss integration points. The tangent matrices of the hydrodynamic 179 forces are considered for large displacement examples as recommended in Ref. [66]. The trapezoidal Newmark 180 numerical method is used to solve the dynamic problem Eq. (24) with α= 1/4,δ= 1/2[2]. 181 9 To solve for rpandrqin Eq. (24), a fourth-order Runge-Kutta algorithm (implemented using the ode45 function 182 in Matlab) solves the nonlinear Van der Pol equations Eqs (16), 17) with constant nodal accelerations ¨u1,¨u2at each 183 Newton iteration. The wake variables (p, q)are defined as uniform in the element. In all of the examples, two 184 stopping criteria are considered based on either the norm of the relative displacement, ||∆u||/u<tolu, or the 185 norm of the residual force, ru<tolr. Both criteria are used to determine convergence of the numerical solution. 186 Algorithms 1 and 2 describe the iterative procedure for resolution. 187 For non-homogeneous initial configurations, the numerical method is initialized with a given deformed config- 188 uration and initial velocities and accelerations. To start from the steady state, the Newton-Raphson method solves 189 the static problem and obtain the initial configuration. Algorithm 1 Numerical procedure 1:Initialize the structure ut←u0,˙ut←˙u0,¨ut←¨u0 2:Initialize the wake (pk+1,˙pk+1, qk,˙qk+1)←(p0,˙p0, q0,˙q0) 3:while t < t fdo 4: Initial guess for t+ ∆t:(uk,˙uk,¨uk)←(ut,˙ut,¨ut) 5: Compute rk+1 uand assemble the tangent matrix 6: while ∥∆uk+1∥>tolu&∥rk+1 u∥>tolr&k <maxiterdo 7: Compute ∆uk+1 8: Update uk+1←uk+ ∆uk+1 9: Newmark method gives next candidates ˙uk+1,¨uk+1 10: Algorithm 2 computes wake variables: pk+1(¨uk+1),qk+1(¨uk+1) 11: Update fk+1 int,s(uk+1),fk+1 ine,s(uk+1,˙uk+1,¨uk+1),fk+1 flu,s(uk+1,˙uk+1, pk+1, qk+1) 12: Compute rk+1 u(uk+1,˙uk+1,¨uk+1)and assemble the tangent matrix 13: end while 14: (ut+∆t,˙ut+∆t,¨ut+∆t, pt+∆t, qt+∆t)←(uk+1,˙uk+1,¨uk+1, pk+1, qk+1) 15: t←t+ ∆t 16:end while Algorithm 2 Wake equations resolution 1:Returns Wake variables candidates pk+1,qk+1of the element for time t+ ∆t 2:Initial guess (pk+1,˙pk+1, qk,˙qk+1)←(pt,˙pt, qt,˙qt) 3:Compute averaged element in-line and cross-flow accelerations ¨uk+1 p,¨uk+1 qusing Eqs (18, 19) 4:Compute rp(pk+1,˙pk+1,¨uk+1 p)andrq(qk,˙qk+1,¨uk+1 q)using Eqs (20,21) 5:ode45 solves rp= 0 6:ode45 solves rq= 0 7:(pk+1,˙pk+1)←(pode,˙pode) 8:(qk,˙qk+1)←(qode,˙qode) 190 4. Experimental methodology 191 In order to validate the numerical model, we carry experimental measurements of cross-flow VIV. The elongated 192 cylinder is 3D printed with Selective Laser Sintering using a Sinterit Lisa Pro printer. The printed material is Ther- 193 moplastic Elastomer (TPE powder) known for its elasticity and flexibility. The cylinder is ℓ= 150 mm in height, 194 10 Figure 2: Experimental setup for U= 0.56m/s (UR= 20 ). Flow direction is from left to right. which is the maximum height achievable by the printer. We selected a circular cross section with a diameter of 195 D= 5mm. This choice strikes a balance between having a high slenderness ( Γ = 30 ) and achieving a satisfactory 196 surface resolution. Since the material is porous, its density ρ= 7.9×102kg/m3is calculated by dividing its wet 197 weight by its volume (given in the design software). A three-points bending test provides the bending rigidity of 198 the structure EI= 4×102N mm2. 199 The structure is placed in a closed water tunnel, which is operated by two pumps. The main pump allows for flow 200 velocities ranging from 0.13m/s to 9m/s, while the secondary pump is used to achieve flow velocities below 0.18 201 m/s. One flow meter measures the flow rate for each pump, which is then converted to a mean flow velocity in the 202 test section. The test section area is 250×250mm2. The bottom end of the cylinder is clamped by screwing an M5 203 hexagonal nut which is embedded to the support plate. The quantities of interest are the amplitude and frequency 204 of the transverse displacements at the tip of the cylinder as they vary with the reduced velocity. The flow velocity 205 Uis varied in steps, and for each step, we wait 30seconds for the steady-state regime to be established. Then, a Go 206 Pro Hero Black 12 records a 4K-quality movie for 20seconds ( ∼100periods) of the transverse displacements at a 207 frame rate of 120frames per second. The camera is positioned 340mm downstream from the structure (Fig. 2). This 208 distance is large enough not to interact with the wake which is a few diameters long. The quantity of interest is the 209 tip of the cylinder transverse displacements. A Matlab image processing script detects the cylinder tip. The frames 210 sequence is first converted to 8bit black and white images and binarized with a threshold, to separate the cylinder 211 from the background. A quadratic regression curve is then fitted on the cylinder pixels to follow its center line. The 212 tip stands on the end of this regression line. We estimate the detection error of the tip 1 pixel = 0.1mm. As the 213 displacements are normalized by the cylinder diameter D, the error is 2%. Each root mean square value of a 20- 214 11 Parameter Example 1 Example 2 Example 3 Example 4 Example 5 Dimensional ℓ(mm) 200 1000 37800 150 D(mm) 10 1 27 5 E(GPa) 0.5 50 36 .2 0 .014 Θ(kN) 0 0 5 0 ρ(kg/m3) 1000 1000 1600 792 ρf(kg/m3) 1000 U(m/s) 0.05 N/A 0.4 N/A Dimensionless M 1 1 1 .6 0 .79 Γ 20 1000 1400 30 St 0.2 C0 D 1.2 2 1 .2 N/A C0 Di 0.2 N/A 0.2 C0 L 0.3 Ap 96 N/A 96 Aq 12 ϵp N/A N/A 0.02 ϵq 0.04 0 .3 0 .04 0 .04 Table 1: Dimensional and dimensionless parameters values for each example. second recording counting more than 2000 values, the standard error becomes negligible compared to experiments 215 repeatability. 216 3. Numerical Results 217 This section presents comparisons to validate the implementation of the model. For all the examples, the 218 fluid considered is water. Gravitational and buoyancy forces are neglected. To ensure deterministic results, we 219 initialize the wake variables pandqwith a uniform value of 0.001. The stopping criteria have default values 220 tolr= 10−5,tolu= 10−8. 221 In the subcritical range 300<Re<1.5×105, the values of C0 D, C0 Di, C0 LandStfor a fixed cylinder are 222 commonly assumed constant [22, 25] with values given in Tab. 1. For all the examples, the fluid considered is water 223 with density ρf= 103kg/m3and kinematic viscosity νf= 10−6mm2/s. While Facchinetti et al. studied cross-flow 224 VIV on a rigid cylinder and recommended Aq= 12, ϵq= 0.3from experiments, so that Aq/ϵq≈40, calibrations 225 from Ref.[62] showed that this ratio for coupled in-line/cross-flow VIV strongly depends on the structural damping 226 and the mass ratio and can reach up to 2727 . In this work, unless otherwise mentioned, empirical parameters of 227 the wake-oscillator are set to Ap= 96, Aq= 12, ϵp= 0.02, ϵq= 0.04. This set of parameters is widely used for 228 in-line/cross-flow WOM and have been validated for slender pipes [68, 40]. The frequencies are normalized by the 229 12 Mode number i 1 2 3 λi 1.87510 4.69409 7.85476 Table 2: Modal parameters for the first three modes of a flexible cantilever beam with constant cross section first natural frequency f1. The ithmode natural frequency of a clamped-free cylinder in water is given by: 230 fi=λ2 i 2πℓ2r EI ms+ma, (26) where λiis the modal parameter for the ith mode of vibration [9]. It’s value for the three first modes are provided 231 in Tab. 2. The codes are made publicly available on GitHub on this link 232 1. Example 1 - Verification of in-line and cross-flow VIV in small displacements 233 This example verifies the implementation of the in-line and cross-flow VIV on a single linear beam element with 234 a semi-analytical solution. 235 The problem consists in a cantilever beam submitted to a transverse flow with a uniform and steady velocity 236 Uexas illustrated in Fig. 3. The material has a Young modulus of Eand a density ρ. The structure has a linear elastic 237 behavior for the internal forces and a lumped mass formulation for the inertial terms. Additionally, we consider 238 small displacements and no structural damping. 239 We only consider the transverse displacement in the yz-plane. The fluid has constant velocity Uand exerts 240 oscillating drag and lift forces. The parameter values are presented in Tab. 1. 241 Consider the distributed nonlinear drag and lift loads: 242 bd(p,˙uy,˙uy) =bd(p,˙uy,˙uy)ex=1 2 CD 0+pCD i 2 ρfDU2 rel, (27) bl(q,˙ux,˙uy) =bl(q,˙uy,˙uy)ey=1 2qCL 0 2ρfDU2 rel, (28) with the relative velocity U2 rel= U2−˙ux/22+ ( ˙uy/2)2. The nodal velocities are averaged to compute a relative 243 velocity constant along the element. For this specific validation case, the change in load direction is neglected as 244 small displacements are considered. The equivalent nodal forces and moments are given by the principle of virtual 245 work. Considering a single linear Bernoulli beam element with length ℓthe equivalent nodal force and moment at 246 the tip Aare1 2bdℓex+1 2blℓeyand−1 12bdℓ2ex+1 12blℓ2eyrespectively [49]. The equation of motion of node Ais: 247  m0 0 0 0 0 0 0 0 0 m0 0 0 0 0  ¨uy ¨θx ¨ux ¨θy +EI ℓ3 12−6ℓ0 0 −6ℓ4ℓ20 0 0 0 12 6 ℓ 0 0 6 ℓ4ℓ2  uy θx ux θy =1 12 6blℓ −blℓ2 6bdℓ bdℓ2 , (29) 13 exezeyA θxθy uyux ℓ D U Figure 3: 3D view of the cantilever beam with small displacements. where mcorresponds to the lumped solid mass and fluid added mass at the free end: 248 m=1 2 ms+Caρfℓπd 22! . (30) System (29) can be simplified to two equations: 249 m¨ux+ux12EI ℓ3+ −bdℓ 8−ux9EI ℓ3 =bdℓ 2, (31) m¨uy+uy12EI ℓ3+blℓ 8−uy9EI ℓ3 =blℓ 2. (32) 14 We set the equivalent stiffness ku= 3EI/ℓ3. The reference system is: 250   m¨ux+kuux=3bdℓ 8, m¨uy+kuuy=3blℓ 8, ¨p=−2ϵpΩ2 f(p2−1) ˙p−4Ω2 fp+Ap D¨ux 2, ¨q=−ϵqΩf(q2−1) ˙q−Ω2 fq+Aq D¨uy 2.(33) This nonlinear reference solution is used to validate the amplitude displacements of the coupled variables (ux, uy, p, q). 251 The Matlab function ode45 using a 4thorder Rung-Kutta algorithm solves this nonlinear system. 252 System (33) is solved numerically over 2000 time steps. The relative error between the two solutions is defined 253 by 254 ϵx=∥xmodel−xref∥2 ∥xref∥2, (34) where xis a vector with the magnitude values for all the time steps. Using this metric, the errors for this simulation 255 are displayed in Tab. 3. Fig. 4 compares temporal evolution of ux, uy, p, q from the model with the reference solution Error ϵux ϵuy ϵp ϵq Value 0.29% 0.23% 0.15% 0.000087% Table 3: Relative errors of ux, uy, p, q . 256 over one period Tof the motion. The model follows closely the reference solution. The doubled vortex shedding 257 pulsation in the in-line Van der Pol oscillator Eq. (16) causes a doubling of frequency of the oscillator p. The present 258 model is in agreement with the reference solution and the relative errors of the four investigated quantities are 259 below 0.3%. 260 15 012·10−3 ux/D −202 p −505·10−4 uy/D 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−202 t/Tq Figure 4: Verification of the in-line and cross-flow displacements ux, uyand wake variables p, qof a cantilever beam in a uniform flow with small displacements. Present model ( ), Reference solution ( ). 16 2. Example 2 - Validation of cross-flow VIV with large displacements 261 This example validates the cross-flow VIV model in large displacements with Leclercq and de Langre wake- 262 oscillator. 263 A cantilever beam of length ℓ, density ρand uniform circular cross section with diameter Dis clamped at one 264 end. It bends under a transverse uniform and steady flow of velocity U=Uexand density ρf(see Tab. 1 for exact 265 values). The beam bends in the xz-plane with large deformation and undergoes cross-flow vibrations along eyas 266 displayed on Fig. 3. In this specific validation example, in-line VIV forces are neglected: the drag coefficient takes the 267 constant value CD= 2. The wake-oscillator parameters for cross-flow VIV are Aq= 12, ϵq= 0.3as recommended 268 in Ref.[22] and used by Leclercq and de Langre. 269 The beam is discretized in 100 elements and the time using 106time steps. The static solution is first computed 270 with increasing reduced velocity. The dynamic problem is then solved using the steady solution as the initial con- 271 figuration. The analysis considers only the second half of the simulation. This cuts off the transient and allows the 272 vibrational dynamics to build up. The amplitudes are normalized by Dand the frequencies by the first structural fre- 273 quency f1obtained by Eq. (26). The root mean square amplitudes and frequencies of the transverse displacement of 274 the cylinder tip are compared in Figs (5,6), respectively. Note that the definition of URused in Ref.[38] incorporates 275 the Strouhal number St= 0.2, resulting in a stretching of the x-axis by a factor of 5. 0 20 40 60 80 100 120 140 160 180 200 220 24000.050.10.150.20.250.3 URuy,rms D Figure 5: Cross-flow RMS amplitude of the tip of the cylinder: present model ( ); Ref.[38] data ( ). 276 For small deflections ( UR<150), the cylinder is excited with single-mode lock-in: the peaks of amplitude at 277 UR= 2.8,UR= 18 ,UR= 56 andUR= 110 in Fig. 5 correspond to the lock-in between the wake frequency and 278 the first, second, third and fourth modes respectively. At those values, only one frequency is excited (Fig. 6). When 279 17 URf f1PSD Figure 6: Cross-flow displacements of the cylinder tip: dominating frequencies from the present model ( ∗) is plotted over the power spectral density (PSD) spectrogram from the Matlab code used in Ref.[38]. the deflection is large enough ( UR>150), the transverse amplitudes reach a lower constant value of 0.1(Fig. 5) and 280 several frequencies are excited at the same time (Fig. 6). The variations in the normal flow profile along the span of 281 the structure induces a lower and multi-frequency response compared with uniform flow, as described in Ref.[38]. 282 The present model is capable of reproducing the same amplitudes and frequencies as in Ref.[38] up to reduced 283 velocities as large as UR= 240 , validating the cross-flow VIV in large displacements and large rotations. 284 18 10110210−510−410−310−210−1 12 NeϵNe Figure 7: Mesh analysis with a reference solution of 1000 elements ( ). 3. Example 3 - Validation of in-line VIV 285 This example aims to validate the implementation of the in-line and cross-flow wake-oscillator model against 286 experimental results form Ref.[65] and DNS results from Ref.[31]. 287 We study a cantilever hollow pipe of length ℓ= 37.8m, external diameter D= 27 mm and internal diameter 288 Dint= 21 mm. It has pinned boundary conditions at both ends. Riser tension is constant, set to Θ = 5 kN (even 289 if experimental measurements show oscillations between 4kN and 6kN). The uniform and steady current has a 290 velocity Uand flows transverse to the pipe. The other parameters are summarized in Tab. 1. 291 A mesh analysis is performed on this example in order to estimate the number of elements needed to accurately 292 solve both in-line and cross-flow VIV. 293 We then compare the numerical displacements dNe(x, t)with a reference solution dref(x, t)overNt= 1000 294 time steps, equivalent to 20periods of motion. The comparison is performed for different numbers of elements, 295 Ne= 1000 being the reference value. The relative error between the displacement histories is given by: 296 ϵNe=Rℓ 0Rtf 0∥uNe(x, t)−uref(x, t)∥1dtdx Rℓ 0Rtf 0∥uref(x, t)∥1dtdx. (35) Figure 7 displays the evolution of the relative error with Ne. The error is inversely proportional to the square of 297 the number of elements: ϵNe∼1/N2 e. In this example, Ne= 100 is enough to have an error ϵNe<10−3. 298 Numerical results are presented with the reference solution Ne= 1000 in Figure 8. The first half of the simulation 299 is ignored in the analysis to cut out the transient. RMS values of both cross-flow and in-line displacements along the 300 19 0.05 0.1 0.15 0.200.20.40.60.81 ux,rms /Dz ℓ 0 0.10.20.30.40.50.60.7 uy,rms /D Figure 8: In-line and cross-flow RMS displacements along the pipe. Present model ( ), experimental data from Ref.[65] ( ), DNS data from [31] ( ). span are compared to experimental results from Ref.[65] as well as DNS simulation from Ref.[31]. The fifth mode 301 emerges as the dominant in-line mode in both the model and experimental data, while the third mode remains the 302 dominant cross-flow mode in all the presented data. The symmetry of the problem about the pipe midpoint leads to 303 a symmetric response of the pipe. This differs from the richer experimental response, underlying the complexity of 304 the VIV process. 305 The maximum RMS displacement is predicted with a relative error of 0.4% for in-line and 0.2% for cross-flow 306 compared with the experimental data. This is often a quantity of interest when assessing structural fatigue from 307 VIV [24], well estimated by this simple numerical model. We obtain general agreement with the experimental data 308 from Ref.[65]. 309 20 0 5 10 15 20 25 30 40 45 50 UR Figure 9: Cylinder center line static deflection in xz-plane with increasing UR. 4. Example 4 - Experimental validation 310 This example compares the model with lab experimental data obtained following the experimental methodology 311 described in Section 2. 312 A cylinder is submitted to a transverse flow with a constant and uniform velocity Uwith clamped and free ends. 313 The simulation and experimental parameters are displayed in Tab. 1. The cylinder has a length ℓ, diameter Dand 314 bending stiffness EI. Unlike previous examples, the mass ratio is lower than 1:M= 0.79. To account for the 315 Reynolds number effect on the drag, the constant drag coefficient C0 Dfollows the empirical model for a circular 316 cross section from Ref.[12]: 317 C0 D(Re) = 11 R−0.75 e + 0.9 1−e−1000/Re + 1.2 1−e−(Re/4500)0.7 . (36) Fig. 9 plots the cylinder center line in the xz-plane for increasing UR, providing a visual representation of the 318 deformations. 319 The experimental data are compared with the numerical results. Fig. 10 displays the transverse RMS amplitudes 320 of the tip uy,rms /Dover the range of interest UR= [0,50]. During the intermediate flow regimes between the two 321 lock-ins, the cylinder experiences very low oscillations, with the experimental RMS amplitude dropping below 0.1 322 forUR∈[8.5,14.3], and after the second lock-in for UR>39(Fig. 10). In Eq. (36), the drag coefficient decreases 323 with the reduced velocity, which might explain the higher amplitudes observed during the second lock-in. 324 Fig. 11 presents the spanwise localization of the RMS transverse amplitudes obtained from simulations for vary- 325 ing reduced velocities. The transverse beam shapes are plotted for two different flow velocities: UR= 4.5and 326 UR= 18 . In the first lock-in, the vibration is dominated by the first mode, while in the second lock-in, mode 2 327 becomes predominant. The ranges of URwhere the first and second lock-in occur in the simulation are coloured 328 with pink and purple rectangles, respectively. First and second lock-ins occur on the range UR= [1.7,7.8]and 329 UR= [11 ,50]respectively, with a gap in between the two lock-ins where the cylinder almost doesn’t vibrate. 330 Fig. 12 displays the frequencies PSD of the tip transverse displacements at different reduced velocities, as obtained 331 from both simulation and experimental data. Fig. 12 also displays the vortex shedding frequency fw=StU/D and 332 21 0 5 10 15 20 25 30 35 40 45 5000.10.20.30.40.50.60.7 URuy,rms /D Figure 10: Cross-flow RMS displacements of the cylinder tip: present model ( ), experimental data ( ). Pink and purple regions represent the mode 1 and 2 lock-in ranges, respectively, where uy,rms <0.2D. URz ℓuy,rms /D Figure 11: Spanwise localization of the RMS transverse displacement for varying UR. Beam shapes transverse displacements (—) appear for UR= 4.5andUR= 18 . Pink and purple regions represent the mode 1 and 2 lock-in ranges, respectively, where uy,rms <0.2D. 22 URff3 f2 f1fw a) b) URfw Figure 12: Spectrograms of the cross-flow displacements of the cylinder tip from the present model (a) and experimental data (b). All frequencies are normalized by f1. the three first natural frequencies of the cylinder f1, f2, f3on horizontal dotted lines, computed with Eq. (26). All 333 the frequencies are normalized by f1= 2.7Hz. 334 The main excited frequency in the simulation follows at first the vortex shedding frequency (Fig. 12), and a 335 strong lock-in on f2is observed for UR>12. The first two lock-in ranges in Fig. 12(a) are dominated with a single 336 frequency, while a multi-frequency response appears for UR>45. This multi-frequency response was observed 337 by Leclercq and de Langre for UR>150in Fig. 6. In a large deflection case, as the elements are closer to the tip 338 they perceive a decreasing normal flow velocity. Thus, when moving towards the tip, the structure is excited by the 339 wake at a decreasing frequency. This spreading of the wake excitation spectrum during large deformations is thus 340 attributed to the variations in the spanwise profile of the normal free-stream component [38]. This feature is also 341 captured in our model. 342 The proposed model effectively captures the overall physics of VIV measured experimentally. Figs 10, 12 demon- 343 strate that the model accurately represents the tip amplitudes, frequencies and the right mode number for the first 344 and second lock-in ranges. The observed decrease in vibration between them is also present. The relative errors on 345 the tip transverse maximum RMS amplitude in Fig. 10 for the first and second lock-in ranges are 19.6% and 0.9%, 346 respectively. However the second one is predicted on a wider range of URthan the experimental data. Note that 347 the WOM empirical parameters could be optimized for this specific study to improve the prediction, but it is beyond 348 the scope of this work. 349 23 U β ℓb ℓ ez exey Figure 13: Coral geometry with two branches. Arrows at the tips indicate for each branch the in-line direction (lighter color) and cross-flow direction (darker color) 5. Example 5 - Branching structure 350 This example presents an application of the model to structures with branching features and constitutes a novel 351 contribution as it is the first time the WOM has been applied to such structures. The planar branched geometry is 352 inspired by the gorgonian soft coral Antillogorgia bipinnata , which grows perpendicular to the ambient flow [54]. 353 The materials and method are the same as in Example 4. The only difference is the structure geometry: one pair 354 of symmetric branches with an angle of β= 60 °and length ℓb= 0.08m is added to the straight cylinder. Fig. 13 355 displays the problem geometry. The planar structure is placed perpendicular to the flow in the initial configuration 356 and the flow velocity is increased by steps. 357 The in-line displacement at all nodes is ux. We define the transverse displacements of the nodes located at the 358 tips of the trunk, the right and left branches respectively as: 359 utrunk t=uy, (37) ub1 t= cos ( β)uy−sin (β)uz, (38) ub2 t= cos ( β)uy+ sin ( β)uz. (39) We keep the definition of the reduced velocity: UR=Uℓ2D−1(m+ma)1/2(EI)−1/2using the initial height of the 360 structure ℓ. 361 Fig. 14 displays the coral displacements with two point of views over half a period of the motion at chosen reduced 362 velocities: UR,1= 2.0,UR,2= 6.5,UR,3= 13 ,UR,4= 19.6,UR,5= 29.3,UR,6= 39.9. The in-line and cross-flow 363 displacements RMS amplitude against URare displayed in Fig. 15. The two branches tip RMS amplitudes are similar, 364 24 UR,1= 2.0 UR,2= 6.5 UR,3= 13 UR,4= 19.6 UR,5= 29.3 UR,6= 39.9 Figure 14: Shape of the structure motion at the end of the simulation for different reduced velocities. The 3D point of view is pointing towards the vector [-1 -1 1] 25 0 5 10 15 20 25 30 35 40 45 5000.20.40.60.811.2 URut,rms/D Figure 15: In-line RMS displacements of the cylinder tips from simulations: trunk tip ( ), left branch tip ( ), right branch tip ( ). Cross-flow RMS displacements of the cylinder tips from simulations: trunk tip ( ), left branch tip ( ), right branch tip ( ), previous example ( ). which is relevant with the problem symmetry. The two main peaks of the trunk tip transverse amplitude occur at 365 UR∼2andUR∼20, just as in the previous example for the cylinder. However, the two peaks are much higher: up 366 to0.6Dand1.2D. The magnitude of the in-line displacements is comparable to that of the transverse displacements 367 within the range UR∈[6−9]. In order to explain these observations, we study the excited frequencies. 368 The spectrograms presented in Fig. 16(a,b) depict the PSD of the in-line and cross-flow displacements at the 369 right branch tip, respectively, as a function of UR. Additionally, in-line and cross-flow vortex shedding frequencies 370 2fwandfw, respectively, are presented. The spectrogram from the left branch (not displayed here) is visually 371 indistinguishable from the right branch. 372 The in-line spectrogram in Fig. 16(a) presents several frequencies aligning with the 2fwline for UR∈[5−10], 373 exhibiting distinct jumps indicative of lock-in regions succession. The stronger lock-in regions on the doubled 374 vortex shedding frequency is consistent with the doubled frequency model in the in-line oscillator Eq. 20. The corals 375 shapes at these frequencies (Fig. 14, UR,1= 2.0,UR,2= 6.5) indeed suggest that vibration mode transitions occur 376 in this range. At low reduced velocities, the coral exhibits transverse vibrations without independent motion of the 377 branches (Fig. 14, UR,1= 2.0). This first mode of vibration explains the first amplitude peak of transverse vibrations, 378 while the in-line amplitudes remain small (Fig. 15). This is consistent with the high PSD region in the cross-flow 379 spectrogram, absent in the in-line spectrogram (Fig. 16). For low reduced velocities, lock-in occurs mainly on the 380 transverse displacements. The presence of the two branches contributes to increased lift in this vibration mode. As a 381 result, the transverse RMS amplitudes increase from 0.47Dfor the cylinder in the previous example to 0.65Dhere. 382 The second peak of amplitude occurs at UR= 20 where the trunk vibrates with a mode 2 (Fig. 14). At UR= 20 , 383 26 URffw 2fw a) b) URfw 2fw Figure 16: Spectrogram of the in-line (a) and cross-flow (b) displacements at the tip of the right branch. Vertical lines ( −·) are plotted for UR,1−6= 2.0,6.5,13,19.5,29.3,39.9. Frequencies are normalized by the first bending frequency of a cylinder in water f1. RMS transverse amplitudes reach up to 1.2D, compared to 0.7Dfor the straight cylinder. 384 In contrast to the straight cylinder case where the natural frequencies of the structure are excited one at a 385 time for UR<45, a notable difference is observed in the case of branched structures. Here, a multi-frequency 386 response emerges for UR,6<6.5, indicating simultaneous excitation of multiple natural frequencies. By adding 387 branches to the structure, the modal response of the structure is modified and several frequencies can be excited 388 at the same reduced velocity, causing a multi-frequency response (Fig. 16). Increased vibrations amplitudes are 389 observed, emphasizing the crucial role of VIV in branched structures. 390 27 4. Discussion and conclusion 391 This study aims to extend the capabilities of the WOM by developing a new formulation that accounts for cou- 392 pled in-line and cross-flow VIV in branching structures. We provide a simple yet accurate tool for simulating 3D 393 fluid-structure interaction problems. The proposed finite element corotational formulation allows to simulate ar- 394 bitrarily complex beam structures, including tree-like geometries. This constitutes the main novelty of this paper 395 since previous applications of the WOM are restricted to straight geometries such as pipes or cylinders. It also inte- 396 grates the WOM and added mass effects into the corotational framework. The numerical resolution of the proposed 397 formulation is implemented and made publicly available in the open-source FEM library ONSAS. 398 The present model is verified and validated through a comprehensive analysis involving four examples of a 399 clamped-free flexible cantilever cylinder in a uniform flow. 400 In Example 1, the in-line and cross-flow VIV formulation is verified against a semi-analytical solution in small 401 displacements. In Example 2, the cross-flow VIV model is validated in large displacements against existing wake- 402 oscillator simulation results. In Example 3, the in-line and cross-flow VIV model is validated against experimental 403 data from the literature. A mesh analysis is also conducted, revealing a second-order convergence. In Example 4, our 404 own experimental data are compared to the model. The maximum RMS transverse amplitudes are well predicted, 405 yet the velocity range of the second lock-in is overestimated. Eventually, two branches are added in Example 5 and 406 the amplitudes and frequency of the trunk and branch tips are presented. This study represents the first application 407 of the WOM on a structure with branches. 408 Our analysis reveals that compliant branched structures subjected to flow exhibit a multi-frequency response 409 with high amplitude of vibrations. The branching leads to a richer and denser spectrum of natural frequencies. We 410 could expect soft corals and other branched structures to exhibit a rich multi-frequency dynamics. 411 Future work can explore the effect of the wake variable interpolation order on element accuracy. It would be 412 valuable to validate the formulation with a fully coupled FSI study of a branched structure. Carrying a modal analysis 413 of the branched structure can enlighten the complex interaction between the structure and the wake. 414 Overall, this study highlights the potential of the wake-oscillator model combined with the corotational approach 415 for predicting VIV in complex geometries, and opens new avenues for the numerical study of the vortex-induced 416 dynamics of branched structures and tree-like organisms. 417 Acknowledgments 418 The authors gratefully acknowledge the funding from the G´enie par Simulation program, as well as NSERC 419 Discovery Grants No. RGPIN-2019-07072. The authors would like to thank Dr. Tristan Leclercq for providing his 420 Matlab code used for generating his data in [38]. 421 28