Adjacency-based, non-intrusive model reduction for Vortex-Induced Vibrations

Vortex-induced vibrations (VIV) pose computationally expensive problems of high practical interest to several engineering fields. In this work we develop a non-intrusive, reduced-order modelling methodology for two-dimensional (2D) VIV simulations. We consider an elliptical, non-deformable solid mounted on springs, subject to a laminar, incompressible flow. Approximating the Arbitrary Lagrangian-Eulerian (ALE) incompressible Navier-Stokes (NS) formulation, a discrete-time, quadratic-bilinear data-driven model structure is assigned for the velocity flowfield prediction. The full-order, data-driven model operators have a predefined sparsity pattern, following the adjacency-based sparsity of the discretized ALE-NS operators. Thus, the data-driven operators inference requires solving many low-dimensional least squares problems, isolating the contribution of"nearest neighbours"for each degree of freedom. Numerical aspects such as data centering and regularization are extensively discussed. With this approach, Dirichlet boundary conditions at the inlet and the fluid/solid interface can be enforced on the full-order, non-intrusive level. Consequently, a non-intrusive reduced-order model (ROM) for the velocity flowfield is obtained after linear projection. The resulting ROM is coupled with the first-principle, 2D solid oscillation equations and is simulated using an implicit time integration scheme. The coupled solution is then mapped on the deformed grid through the inverse ALE map. This methodology is showcased for two testcases, at different Reynolds numbers (Re = 90, 180). Numerical results indicate a successful coupling between the data-driven velocity flowfield and the solid oscillation, with prediction errors of less than 3% for both the flowfield and the solid oscillation. A comparative study with respect to the ROM dimension indicates the robustness and potential of the approach.


Introduction
Vortex-Induced vibrations (VIV) comprise a class of Fluid-Structure Interaction (FSI) problems with high practical interest to numerous engineering fields, among which are wind, offshore and aerospace engineering [41,46].Vortex-induced vibrations concern the two-way coupled system response of a non-deformable solid body mounted on elastic supports, subject to a fluid flow.A multitude of complex dynamical phenomena arise from the interplay between solid and fluid dynamics, through the mechanism of asymmetric vortex shedding in the wake of the solid body [3].
Vortex shedding patterns, as well as oscillation coupling properties can greatly vary, depending on the flow conditions and solid parameters of a given configuration [16,22].
In practice, the dynamically varying loads resulting from vortex shedding can lead to strong vibrations and thus a significant decrease in the fatigue life of structures.Such vibrations could be proven particularly destructive in applications such as underwater pipelines [21], turbine poles [9] or structural cables [20] and are thus carefully accounted for during engineering design.Apart from VIV suppression, intensive research is being performed on efficient harnessing of the kinetic energy of vortex-induced vibrations for energy production [17,49].Both these fields indicate the strong motivation in developing efficient VIV simulation tools for engineering design optimization and control.
Over the past decades, a number of methods have been developed for the direct numerical simulation of FSI and VIV problems [2,14,45], confronting issues such as numerical stability and computational overhead [27,38].Concise reviews of such numerical methods for FSI problems are given in [12,19].Nonetheless, even for regimes of incompressible, laminar flows, the two-way coupling of two highdimensional subsystems (i.e. a CFD and a structural finite element model) remains a challenging task.Depending on the problem at hand and the employed numerical solution, stability is not always guaranteed [26], while the stiff coupling and high dimensionality introduce a significant computational cost [38].
The growing field of model order reduction (MOR) comprises a potential alternative to bridge the gap between accuracy and computational overhead.MOR methods have proven to be efficient for several demanding, high-dimensional fluid and solid dynamics applications, leveraging computational cost and accuracy [8,15,33,42].A number of research works have also been presented for different FSI aspects and applications.Notably, Liberge et al. employed the method of Proper Orthogonal Decomposition (POD) to the flow around an oscillating cylinder, transferring the FSI problem to a global, fixed domain [24].Using experimental Particle Image Velocimetry (PIV) data, Riches et al. also employed POD, to investigate the intrinsic mechanisms of VIV wake-dynamics.A reduced basis method differentiating between solid, fluid and deformation modes for intrusive model reduction was developed by Nonino et al., applied to both partitioned and monolithic FSI numerical models [32].Focusing on parameter-dependent problems, Benner et al. proposed a low-rank method to increase the Newton iteration's computational efficiency for FSI computations [6], while Lieu et al. produced a POD-based model for the aeroelastic response of a complete F-16 aircraft, interpolating the model response over varying free-stream Mach numbers [25].
In cases where only simulation or experimental data are available, non-intrusive model reduction methods are being employed for modeling and prediction purposes.Working in this direction, Poussot-Vassal et al. [34] used the Loewner framework for the prediction of aircraft response to gust forcing.Focusing on the aeroelastic phenomenon of stall flutter, Dai et al. developed a Recurrent Neural Network architecture for the non-intrusive prediction of pitching airfoil dynamics [11].In a more general framework, Xiao et al. [47] presented a non-intrusive reduced-order model for FSI problems, based on radial basis function interpolation over time.In a slightly different direction, Yao et al. [48] focused on VIV problems at low Reynolds numbers, proposing a linear, Eigenvector Realization Algorithm (ERA) approach to explore the mechanism of transition to vortex shedding for different VIV configurations.In this study, we focus on laminar, low mass ratio, 2-dimensional VIV testcases [30,35], as a prototype example for the proposed non-intrusive, ROM methodology.
In further detail, we present a methodology for datadriven model reduction, focusing on the prediction of 2D laminar, vortex-induced vibration dynamics.Assuming state access, a non-intrusive model for the incompressible fluid flow is constructed on an ALE reference grid.To account for the grid deformation, the velocity data being interpolated to a reference domain through an ALE map, given by the solution of a Laplace equation.We show that the structure of the ALE-NS formulation for the velocity flowfield is approximately quadratic-bilinear.Motivated by the adjacency-based sparsity of the discretized ALE-NS operators, we propose a novel method to construct a full-order, data-driven model for the velocity flowfield, with an a priori adjacency-based, sparsity pattern.As a result, we locally infer the data-driven operators by solving one least squares problem for each internal grid node with an L 2 regularization term.In parallel, a proper treatment of the boundary nodes allows for a direct enforcement of the Dirichlet boundary conditions at the fluid inlet and the fluid/solid interface on the full-order level.The data-driven fluid subsystem is then coupled with the first-principle solid motion oscillations through a time-implicit scheme.The map to the deformed configuration is also considered as a post-processing step of the VIV predictions.The proposed methodology is showcased for two VIV cases of an oscillating solid along the streamwise and transverse directions, subject to a Re = 90 and Re = 180 flow.Through these testcases, several properties as well as the predictive capabilities of the method are highlighted.
The rest of this work is structured in the following way: A review of the theoretical background on vortex-induced vibrations is given in Section 2, motivating the non-intrusive model formulation for the velocity flowfield.The developed methodology is analyzed in Section 3, for the construction of a non-intrusive fluid dynamics ROM and its coupling with the solid oscillations along the transverse and streamwise directions.Finally, the model is tested for two testcases with Re = 90 and Re = 180 flow past an ellipse-shaped body, with corresponding results given in Section 4. Conclusions and potential future work is discussed in Section 5.

Theoretical Background
In this section, we present the physical modeling of laminar, vortex-induced vibration problems, which motivates the structure of the developed non-intrusive, data-driven model for the fluid velocity field.

Vortex-Induced Vibrations
As schematically depicted in Figure 1, the problem consists of an incompressible fluid flow over a non-deformable body that can oscillate along the streamwise and transverse (x and y) directions.The equations of motion for the solid (s) displacement d s are where n is the unit normal vector on ∂S pointing from the solid to the fluid and K is a diagonal matrix with the spring constants k x , k y .The cross-sectional area of the solid is A s , the solid and fluid densities are ρ s and ρ f respectively and σ is the 2-dimensional fluid stress tensor.
The integral term of (1) encodes the dynamic coupling condition between the fluid and solid subsystems, such that the equations of motion are essentially two uncoupled, externally forced oscillations.Internal damping has been set to zero to encourage the onset of vortex-induced vibrations.
For the fluid flow, we focus on the regime of 2D incompressible, laminar flows with velocity u, pressure p, kinematic viscosity ν f , gravity acceleration g and density ρ f .The stress tensor σ for a Newtonian fluid is: where I is an appropriate identity operator.The imposed boundary conditions to the flow are the kinematic coupling condition at the fluid/solid interface, the velocity at the inlet and the "do-nothing" condition at the outlet [38]  From the Eulerian perspective, the fluid domain Ω(t) is changing over time (see Figure 1).Through the ALE framework [12], the Navier Stokes are mapped to a reference domain Ω, which could be selected to be the domain configuration at t = 0.The map from the reference domain to the current configuration is then: x → x = x + d(t).The deformation field of the current configuration with respect to the reference one is given through function d(t) (ALE map).For an analytical discussion of possible definitions of the ALE map, the reader is directed to Chapter 2.5 of [38].Considering modest gird deformations, the ALE map d(t) is computed through the solution of a Laplace problem with boundary conditions on the interface Ω ∩ S, and on the domain external boundaries: We notice that using definition (5), the only time dependence of the ALE map originates from the boundary conditions (6) and that d depends solely on the solid deformation d s (t).This observation will be proven useful in the ROM level, since the ALE map can be inverted at a very low added computational cost.Denoting F as the Jacobian of the ALE map, and J = det F , we get the ALE formulation for the Navier-Stokes equations [32,38]: From the model reduction point of view, it is of particular interest to construct the fluid dynamics model on domain Ω, following a different approach from [24].In this way, we would isolate the flowfield snapshots from the solid motion, thus avoiding typical issues of transport-dominated equations projection [40] close to the FSI interface.The ALE formulation (7) of the Navier-Stokes is suitable for this task.However, the corresponding structure of the equations is no longer quadratic, as in typical, fixed domain Navier-Stokes [5].Instead, the structure of ( 7) is highly-nonlinear.
If we assume that the grid deformation is small, we can approximate F ≈ I and J ≈ 1.This simplifying assumption is made to reveal an approximate model structure on the data-driven level, however it could constrain the datadriven model towards modest grid deformations.Such a limitation is either way dictated by the usage of the ALE map (5).Simplifying (7) based on the above, we get while the pressure can be computed by taking the divergence of (8): The structure of (8), as well as the algebraic link between velocity and pressure (9) will comprise the foundation towards the fluid data-driven model.

Data-Driven Fluid Model Structure
Discretizing equation (8) of the previous section with given time-implicit and spatial schemes, results in a quadraticbilinear structure for the evolution of the velocity flowfield where A u encodes information from the viscous term and the time derivative discretization in (8).Through the Kronecker product, quadratic (H u ) and bilinear (K u ) terms are included, corresponding to the advection term of (8).Q u encodes pressure dependence, while terms B u , L u are introduced to enforce the kinematic boundary conditions on the solid-fluid interface and the inlet.From here on, we refer to the velocity flowfield on the reference configuration and thus notation ˆis dropped in the following.All involved operators are sparse, originating from the discretization of the operators in (8).The model is written in an implicit way; apart from the linear terms, all right-hand-side velocity terms are taken at the (k + 1)st timestep.This is related to the coupling of fluid and solid subsystems, later discussed in this work.Furthermore, discretizing (9) in a similar way reveals the algebraic relation of pressure with velocity, extensively discussed in [5] where A p originates from the outflow "do nothing" boundary condition and C p is the reference pressure.Matrices E p , A p , H p , K p are banded and sparse, based on the employed spatial discretization stencil.
The inverse of E p can be shown to be band-dominated [7], while the corresponding band can be reduced through several re-ordering methods, e.g.[36].We thus assume that after substituting equation ( 11) in (10), the resulting operators can be approximated by sparse matrices, with non-zero elements only in the positions corresponding to adjacent nodes of any given internal grid node and thus, any matrix row.
Finally, the structure of the data-driven model for the flowfield velocities on a reference domain takes the following, quadratic-bilinear form (12) with sparse operators A, H, K and a priori known operators L, B for the enforcement of the Dirichlet boundary conditions at the flow inlet and fluid-solid interface.

Average Flowfield Removal
We hereby make a note on the decomposition of the flow into a dynamical and a stationary, average flowfield.For the flow around an oscillating cylinder, the mean flowfield ū over the reference domain Ω is omnipresent and includes a significant part of the flow energy.This can be viewed in Figure 2, where the singular values of a 3000 × 500 flowfield snapshot matrix is plotted for velocity components u x and u y , before and after removing the corresponding average of the time series.It is observed that the first u x singular value corresponds to the average ūx , indicated by the observed singular values shift for ũx (t) = u x (t) − ūx .This effect is less evident for the u y component, since the average ūy is almost zero in this case.As a result, we write the velocity field as: Excluding the mean ū is known to have a beneficial effect [29] on the numerical manipulation of the data.As indicated by Figure 2, removing the mean reduces the condition number of the corresponding data matrix.Thus, assuming we have state access over some training time [0, T 1 ], we subtract the time-average velocity field over this interval from the data and build a ROM only for the ũ(t) component.From this point and on, we refer to the prediction of ũ(t), dropping the ˜notation for simplicity.It can be observed that subtracting ū in (12) by substituting (13) leaves the prescribed sparse, quadratic-bilinear structure unchanged.

Grid Deformation
Model ( 12) encodes the physics-based structure of the moderatelydeformed domain, incompressible Navier-Stokes ALE formulation in (8).Thus, it corresponds to the flow solution on a reference configuration.Before examining the datadriven solution of ( 12) for the reference configuration, we need to complement the model with a map from the reference domain to the current, deformed one.
For the non-intrusive modeling procedure, the flow velocity data over some training time t = [0, T 1 ] is imported, on a given grid Ω(t).The grid adjacency information is needed for the application of the aforementioned adjacency-based sparsity pattern.
The deformation of the grid over [0, T 1 ] is often also provided in case of proprietary software.However, different solvers employ variations of ( 5) to derive the grid deformation [32,38], often not traceable by the end-user.Also, in the case of experimental data (e.g.PIV campaigns [29]) there is no notion of grid deformation.
Thus, it is necessary to construct a new fluid mesh and compute its deformation based on a predefined ALE map (in our case ( 5)).Hence, the imported velocity data need to be interpolated at each training timestep from the solver grid positions (or e.g.PIV domain in experimental campaigns), to the positions of a re-computed grid, based on the selected ALE map.For the constructed mesh, grid adjacency information is stored, allowing for the application of an adjacency-based sparsity pattern in the non-intrusive model.
The motion of each grid node is given by discretizing and solving (5).As in standard practice, we denote Λ as the discretized Laplace operator, augmented by identity matrix rows for the Dirichlet boundary conditions in (6).I denotes the set of nodes on the interface Ω(t) ∩ S, which should follow the motion of the solid (see eq. ( 6)).The displacement of the grid nodes is then given by where the reference configuration is taken as the one at t = 0. We observe that i∈I Λ −1 :,i should be computed once and then the grid deformation can be computed directly for a given time t, based on the solid motion ((d s (t)) prediction.In essence, (14) gives an analytical expression of a basis for the grid displacement d with respect to the solid displacement.

Data-driven, Adjacency-based Sparsity
The task in the case of non-intrusive modeling is to infer the unknown operators of (12), i.e. to construct a datadriven model for the 2D velocity field.Typically, the model operators are inferred by first projecting the data to the leading SVD modes of the data matrix [33,42].However, we hereby follow a different approach, motivated by the discussed adjacency-based sparsity of the discretized operators which we want to infer.
Following this approach is also accommodated by the application-based need for grid generation, already discussed in Section 3.3.In particular, using the grid adjacency information we can identify adjacent nodes and extract corresponding velocity data.By examining each degree of freedom, we define a sparsity pattern for the full-order, inferred operators.This idea, briefly analyzed in [1], is here applied to all internal nodes of the 2D grid, for the prediction of u x and u y velocity components.
Enforcing sparsity enables computing the full-order, datadriven operators in (12), since only neighbouring node products of terms K ∂ t d s ⊗ u k+1 and Hu k ⊗ u k are non-zero.This also allows enforcing uniqueness of the data-driven model with respect to the states (velocities), by avoiding commutative velocity products.Such a property cannot be enforced when constructing the non-intrusive model in the reduced space, since the reduced vector elements consist of linear combinations of the states.Nonetheless, it is evident that the offline cost of this approach scales with the number of nodes on the interpolated grid and thus is expected to be significantly elevated, compared to that of projection-first methods.
In practice, to apply the aforementioned sparsity pattern based on grid adjacency we examine one row of (12) corresponding to some internal node.Dirichlet boundary conditions for nodes at x = 0 (i in ) and nodes on F ∩ S (i s ) can be a priori satisfied by setting A iin,. = H iin,. = K iin,.= 0, C iin = 0 and A is,. = H is,. = K is,.= 0, C is = 0, while entries of 1 are registered in respective positions of B, L. Thus, matrices B, L are grid-dependent, a priori known, permutation matrices.
We split the quadratic terms into two contributions h A and h B ; the first one includes velocity products including the velocity of the examined node, while the second one includes products between adjacent node velocities.Based on this, we now look more closely into a row of (12).The equation for the u x or u y velocity of a node i with a set of adjacent nodes q(i) is then where we set z = 1 for u x and z = −1 for u y respectively.In this general form, the terms of h A , h B are not unique.In practice, we eliminate the terms of h B involving the velocity components of node i so that h A and h B share no common terms.We also eliminate non-unique, commutative terms in h B .As a result, if node i has m neighbouring DOFs, then h Bi includes m 2 terms [33].Equation (15) for node i can be written for every timestep in [0, T 1 ].This leads to the formulation of a least squares problem for the entries of operators A, H, K min ai,.,hA i,. ,h B i,. ,ki,.
where, following the notation of [5], The solution of ( 16) for an internal DOF i is a vector with the entries of line i of the operators A, H, K. Based on grid adjacency information, these entries can then be stored in corresponding positions of the operators, thus leading to an adjacency-based sparsity pattern of the full-order, nonintrusive model.

Regularization
Adding regularization is required for the practical solution of (16), to avoid numerical errors due to small singular values of D. We employ Tikhonov regularization [28,44] to penalize solutions with high L 2 norm by modifying each least squares problem as follows: For the current application, setting λ 1 = λ 2 yields satisfactory results.Hence, for each degree of freedom, we solve problem (18), for different values of λ = λ 1 = λ 2 .The optimal regularization parameter is chosen based on the L-curve criterion [18].In particular, we aim to leverage between the least squares solution error and the solution norm.Denoting b 2 as the least squares solution error and x 2 as the solution norm, we seek for where ˆindicates normalization of the b origin is chosen.
In the full-order setting, we can reduce the computational cost by spatially interpolating the optimal regularization parameter.In essence, since velocity varies smoothly over the domain (dealing with an incompressible flow), matrix D in (17) and thus its singular values are also expected to vary smoothly over the domain.We can thus compute the optimal regularization parameter at randomly sampled nodes and approximate λ at the remaining nodes by interpolation.
By summing the contribution of x 2 and b 2 for each λ value over all the internal DOFs, we can plot a global L-curve.The corresponding sums can be calculated for the optimal λ (at each DOF), i.e. the training error of the model and norm of the inferred operators.The inferred model training error and operators norm will not necessarily lie on the global L-curve, since the optimal λ is independently computed for each DOF.
Figure 4 indicates the trade-off between the non-intrusive model training error and inferred operators norm, for a Re = 180 flowfield, further analyzed in Section 4. By spatial interpolation of the regularization parameter, the least squares computational cost is significantly reduced.For example, if 20 λ values are considered and the optimal λ is computed for 10% of the DOFs, then only one least squares problem should be solved for the remaining 90% of the DOFs.This translates to a decrease in the involved least squares computational cost by a factor of ≈ 6.5 (without considering the cost of λ interpolation).A zoomed-in view of the global L-curve is depicted in Figure 4a.A minor increase in the model training error when interpolating λ is recorded, compared to the model computed from solving (18) for all degrees of freedom.Figure 4b illustrates the optimal regularization parameter for each grid node based on (19), for the u x velocity component.We observe that λ is higher in the wake region of the flow, where the flow exhibits more complex dynamics.Figure 4c corresponds to the interpolated λ case, computing the optimal parameter for 10% of the internal grid nodes.Figure 4b and Figure 4c exhibit differences on the selected λ, however the resulting difference in training error is only minor.It is noted that parallelization of the least squares problems could also significantly reduce the offline clock time required for inferring the operators through the employed methodology, however, this has not been considered in this work.
At this step, the sparse, full-order, data-driven model of ( 12) has been computed.Next, the force induced from the flow to the solid should be predicted, as to couple the fluid and solid dynamics.

Surface Forces Modeling
As shown in (1), the forcing from the fluid flow to the solid motion originates from the integral of stresses on the body surface (eq.( 3)).Focusing on non-deformable solid motion along the two axes, a model for the forces F x , F y (Figure 1) is required.The mean velocity flowfield over the training time [0, T 1 ] results in a constant force F = Fx , Fy T .As a result, we can similarly write the force as By consideration of the examined application, F should tend to zero as training time increases, since the solid motion is oscillatory.We aim to predict the dynamically evolving component F, in line with the focus on the dynamically varying velocity flowfield component.In the following we drop ˜for simplicity.
We assume solid motion parameters (mass, spring constants) to be known.Mass and (linear) spring constants can be easily determined both numerically, but also experimentally through the impulse response for the solid oscillation.Therefore, only kinematic data is necessary for constructing the non-intrusive model.The velocity training data for the solid, as well as the known mass and spring constants are substituted in (1) to obtain numerical training data for the integral forcing terms we aim to predict.Based on the typical decomposition of flow-induced forces to a shear term (linear) and a pressure term (quadratic), we hereby employ a non-intrusive, quadratic model for the force, in the following form: For this part, a second adjacency matrix is formulated.Equations ( 1) and ( 3) indicate that only the velocity of nodes sufficiently close to the body will contribute to the force acting on it.As a result, the full-order output model ( 21) has 2n F DOFs (velocities u x and u y ), where n F << n is the number of nodes at the proximity of the interface Ω ∩ S.An illustration of this proximity grid is given in Figure 5.  terms of u ⊗ u were included in the model, but only the ones corresponding to adjacent pairs of DOFs.In fact, using only quadratic terms of the sort u i u q(i) , where q(i) denotes the adjacent DOFs for each DOF i and u can be u x or u y , proved to be sufficiently accurate for body forces prediction, while restraining the involved least square problems dimension.
As presented for the velocity flowfield in (3.4) and (3.5), we formulate the full-order matrices A F , H F by solving an L 2 -regularized least squares problem min where A F , H F are of dimension 2 × 2n F .

Coupled ROM-Oscillation Dynamics
Up to this point, the sparse, full-order, non-intrusive models for the velocity flowfield and the force on the solid have been computed.The non-intrusive ROM is then computed via projection, using a POD basis.The flowfield model is projected to the leading singular modes of data matrix U .The columns of U are vectorized velocity data over the new fluid mesh (see (3.3)), stacked over training time [0, T 1 ].Following the logic of Proper Orthogonal Decomposition, we compute the SVD U = ΦΣΨ T and truncate the leading r singular values.Denoting Φ r = Φ .,1:r, we project the grid velocities u k to the leading r modes, such that Applying (23) to model (12), we obtain the reduced matrices as follows: We observe that due to the sparsity of matrix H, we can efficiently perform the projection to H.In particular, we split matrix H into H A with entries including the "self node" terms (see eq. ( 15)) and H B containing unique products of neighbouring node velocities.We store only the non-trivial entries of these two matrices; as a result, we need to store matrices of dimensions 2n × 2 max i q(i) and 2n × 2(maxi q(i)−1)

2
, respectively for H A and H B , where n is the total number of nodes (velocities u x and u y ).This allows using a different computational stencil for each of the two terms H A and H B .In this application, we use the second-degree adjacent nodes for H A and the first-degree adjacent nodes for H B , limiting the number of required Kronecker products.Since the two matrices share no common non-zero entries, we can then write Knowing the sparsity pattern of H A and H B , we compute only the relevant Kronecker products of Φ r ⊗ Φ r for each case.Once matrix HΦ r ⊗Φ r is computed by the above procedure, it is multiplied on the left with Φ T r .This procedure significantly reduces the number of operations required for the projection of the quadratic matrix.
For the force acting from the flow to the solid, we follow a similar procedure: We identify the lines of Φ r corresponding to the nodes close to the body from the grid adjacency information (as shown in Figure 5).We denote Φ s the 2n F × r matrix by retaining the necessary lines of Φ r .Then, the projection step for ( 21) is Having obtained a non-intrusive model for both the velocity field and the resulting forces on the solid along the streamwise and transverse direction, we are ready to couple the fluid (data-driven) model with the solid (first-principle) one.The data-driven ROM writes as while the solid oscillation equations ( 1) are integrated using a Crank-Nicholson scheme, with timestep ∆t: An interesting observation during simulation was that the known stiff coupling between the fluid and solid subsystems in FSI numerical simulation [10,39] transfers also to this data-driven/first principle two-way coupling, rendering an explicit model unstable.Thus, the model is constructed with an implicit formulation, requiring convergence at each timestep.In detail, the force F k resulting from the datadriven flowfield ROM ũk is driving the solid oscillation.The resulting solid displacement d k s is in turn affecting the surrounding flow ũk .A combined convergence criterion of ũk j+1 − ũk j 2 < 10 −12 and d s < 10 −12 turned out to provide good results, where j denotes the index of iterations within timestep k.For each implicit iteration, a typical successive relaxation method was used [26], with ũk * = (1 − ω)ũ k j + ω ũk j+1 .During simulation, a maximum of 4 − 5 iterations per timestep where required to reach convergence.

Simulation Results
The above methodology was coded in MATLAB and applied to two numerical testcases for the VIV of an ellipseshaped solid.The details for the two testcases are given in Table 1.The solid has streamwise and transverse eigenfrequencies (in vacuum) of f s = 4.38 Hz and is subject to a laminar, incompressible channel flow.A parabolic velocity profile is prescribed at the channel inlet, with a maximum velocity of u max = 1.5ū in , while the domain size is 4 × 1 m 2 .
The reference length is twice the average of the semiaxes D = R x + R y and the reference velocity u ∞ is the average inlet velocity.The two testcases then correspond to Re = 90 and Re = 180.A low mass ratio is also selected (ρ s /ρ = 1.2), to promote a wide VIV lock-in velocity range and considerable peak amplitude [46].The added mass is m A = 0.011 kg, which leads to eigenfrequencies of the submerged solid of f N ≈ 3.24 Hz.The reduced velocity U * = u ∞ /(f N D) in the two testcases is then U * = 1.93 and U * = 3.86.This indicates that they respectively belong to the right and left parts of the initial VIV branch [37].All numerical simulations were performed with the GASCOIGNE open-source software [4] on a computational grid of 10256 nodes and simulation time of 15 and 10 seconds respectively for Re = 90, 180.

Table 1: VIV simulation properties: Two testcases with
Re = 90, 180 were examined, at a low mass ratio The grid interpolation step described in (3.3) was performed focusing on a region close to the body.An unstructured mesh of 2085 nodes is constructed on a truncated domain of [0.2, 2.4] × [0.1, 0.9], with the use of the MESH2D MATLAB tool [13].The data from the last 5.9 seconds of simulation were used in both cases, including a transient response.63% of the time series was used as training data, based on which the reduced-order model is constructed.The coupled system of equations ( 27) and ( 28) was then simulated and VIV predictions beyond the training time were examined.

Velocity flowfield POD Basis
Before presenting the ROM predictions, it is interesting to examine the POD basis onto which the data-driven model for the velocity flowfield is projected (see Section 3.6.1).We will hereby focus on selected modes Φ i for the two examined cases, which correspond to three significant phenomena in VIV dynamics.This analysis will help to interpret the sim-ulation results obtained with our approach in Sections 4.2 and 4.3, as well as provide insight into the projection step of Section 3.6.1 and the main mechanisms of VIV dynamics.
For a detailed study focused on the POD of experimental VIV data, the reader is directed to [37].
After removing the data average (see Section 3.2), the first POD mode pair corresponds to convective vortex shedding for both Re = 90 and Re = 180.This is similar to the case of the flow over a cylinder [31], with this mode pair containing more than 50 % of the total kinetic energy of the flow.We hereby show only the first mode for Re = 90, u x and u y (Figure 6a,b), since the second mode is just shifted by one quarter of the spatial wavelength.For a Strouhal number of St = 0.21 and choosing a reference velocity as u max = 1.13 m/s, the expected vortex shedding frequency is f v ≈ 1.97 Hz.This is in good accordance with the peak at 1.72 Hz of the Ψ 1 Fourier spectrum, given in Figure 6c.
A mode indicatory of the slow-drift of the mean flow is also detected in both datasets.This mode (Φ 5 for Re = 90) captures the energy of flow transition to a limit-cycle behaviour and is illustrated in Figure 7 for both velocity components.It is clear that the significance of such a mode depends on the presence of transient dynamics in the VIV simulation dataset and in many cases is individually treated [37].The involved low-frequency dynamics are indicated by the corresponding Fourier spectrum in Figure 7c, while two minor peaks are noticed close to the vortex shedding and solid natural frequencies.
In the case of Re = 90, other POD modes correspond mainly to multiples of the vortex-shedding frequency and corresponding wavelengths, since the solid forcing frequency is sufficiently below f N = 3.24 Hz.However, in the case of Re = 180, the two frequencies are very close.This is specifically evident in Φ 5 of the corresponding basis, shown in Figure 8.The Fourier transform of Ψ 5 (Figure 8c) shows a strong peak at 3.12 Hz, indicating the presence of the solid natural frequency in the system dynamics.This is expected theoretically [46] for Re = 180 and can be captured by the ROM with more than 5 modes.It is noted that a departure from the exact value of f v /f N = 1 can be expected in cases of low mass ratio and damping [23], as the one considered here.

Solid Oscillations
The oscillation velocity time series and corresponding ROM predictions are given in Figure 9 for both testcases.The vertical dashed line indicates the end of the training time.From that point and on, the model predictions are tested against unseen simulation data.In both cases, the streamwise oscillation amplitude (along x) is more than one order of magnitude lower compared to the amplitude of the trans-   verse oscillation.Furthermore, the transverse oscillations dominating frequency (f ≈ 1.72 Hz for Re = 90) is half than that of streamwise motion, as theoretically expected.
For the case of Re = 90, the oscillations exhibit a very low-amplitude (Figure 9b), hinting a link to the lower end of the initial VIV branch [30].On the contrary, the Re = 180 case shows strong transverse oscillations (Figure 9d), almost one order of magnitude higher than that for Re = 90, often exhibited in the higher end of the initial VIV branch.Due to the vicinity of f N and f v , beating phenomena are observed for Re = 180 [43].The coupled data-driven-oscillation model captures the solid dynamics with high accuracy in both cases.
The qualitative differences in the transverse and streamwise oscillations are reflected on the ROM predictions error.The truncation of the SVD modes ( 23) leads to a corresponding truncation of the flow energy.If the VIV coupling is strong (as in the streamwise direction), even leading flow modes interact with the solid oscillation (see e.g. Figure 8).However, in cases of weaker VIV coupling (as in Figure 9a,c), it is likely that low energy modes that interact with the low-amplitude solid motion are truncated.Hence, the relative error for the r = 30 ROM prediction e.g. in Figure 9c might seem significant at specific times, however the absolute ROM prediction error is negligible.

Flowfield Predictions
For each simulation timestep, we can reconstruct the predicted velocity field in the reference configuration, by u r = Φ T r Φ r ũ.However, it is of interest to also map the reconstructed flowfield to the predicted current configuration.This can be done by using the dimensionless deformation field i∈I Λ −1 :,i from ( 14) and multiplying with the predicted displacement of the solid body along both axes.The solution is then mapped to node positions of the predicted deformed grid.This last step entrails a negligible added cost, since i∈I Λ −1 :,i has been precomputed.The CFD and predicted ROM flowfields for the final testing time are given in Figure 10 and Figure 11, for both velocity components of Re = 90.A corresponding comparison is given for Re = 180 in Figure 12 and Figure 13.The CFD solution indicates a 2S mode for both testcases, where two single vortices appear per oscillation cycle.This matches well-known results for laminar VIV with free body oscillations [30,46].ROM predictions in both cases match the vortex dynamics qualitatively and quantitatively, employing r = 30 DOFs.It is noted that the full-order numerical simulation (e.g. for Re = 180) requires approximately 40 minutes on a personal laptop.In contrast, the offline ROM cost corresponds to approximately 3 minutes, while the ROM numerical simulation requires just several seconds   also compute an average, relative ROM error over the u x and u y flowfield.This reads (e.g. for the u x component) as Error e is recorded in Figure 14 over time, for both testcases and both velocity components.For Re = 180, the error exhibits oscillations, possibly linked to both modeling and projection errors of the stronger vortex dynamics, compared to Re = 90.We observe that in the transient part of the data, the error reaches its peak value, approximately 1.5%.It is noted that the flowfield error originates from the combined effect of modelling and projection errors, as well as an error from the coupling of ( 27) with (28).In the testing time interval, the errors for Re = 180 (both u x and u y ) oscillate for values below 1%.For Re = 90, the error shows a noticeable increase after the training time, reaching up to 2% at the end of the testing time.This is primarily linked with a phase shift, slightly noticeable in Figure 9a.Transient dynamics induced by the slow, mean flow drift   are not fully resolved during the training time (see the lowfrequency peak in Figure 8c).Thus, a gradual phase shift of the ROM predictions is observed during the testing time interval, leading to a respective, almost linear increase of the error in Figure 14.Finally, we perform a comparative study by averaging the error (29) over the testing time, with respect to the dimension r of the ROM.As presented in Section 3.4 and Section 3.6.1, the full-order data-driven flow model is independent of the projection basis Φ r (eq.( 23)).After truncating the projection basis to 10 different r values, we simulate the resulting coupled system and record the average of (29) from the end of training time T 1 to the final simulation time T .Figure 15 illustrates the obtained results for both testcases and both velocity components.In the case of Re = 180, the resulting u y error is slightly higher than u x , especially for r ≤ 10.For small values of r, the average error is higher for the Re = 180 case.That can be explained by the more rich dynamics exhibited for this case, as shown in Figure 8 and Figure 9c,d.Thus, comparably more POD modes are required to sufficiently capture the system dynamics.In both cases, there appears a sharp corner for a ROM dimension close to 20, meaning that adding further POD modes does not provide any further information on the model dynamics [5].However, Figure 15 also indicates a convergence of the ROM error to a certain value with increasing ROM dimension r.This property is typically acquainted in intrusive MOR (e.g.[25,32]), here inherited by the presented non-intrusive ROM methodology.A similar error convergence behaviour with the one observed for the non-intrusive model in Figure 15 was reported for an intrusive FSI ROM in [32].This indicates that the independence of the inferred sparse, full-order operators from the projec-tion basis Φ r could result to an increased robustness of the non-intrusive ROM.[15,20] is sufficient for the ROM, with an observed error convergence for a further r increase.

Outlook and Future Work
In this work, a non-intrusive model order reduction methodology was presented, applied to vortex-induced vibration problems.The dynamical model structure for the fluid velocity was motivated by the problem physics under the ALE formulation, while a physics-based sparsity pattern was enforced through the grid adjacency information.An L 2 regularization term was added to the corresponding least squares problem for each grid node velocity component and computational cost reduction strategies for the optimization of the regularization parameter were discussed.The data-driven velocity field ROM was subsequently coupled with the first-principle solid oscillations and a method for mapping the solution from the reference to the deformed configuration was presented.Finally, the methodology was applied to two transient VIV testcases with Re = 90, Re = 180, which exhibit different VIV phenomena.Results on the prediction of both the solid oscillation and the surrounding flowfield with only 30 DOFs indicated the potential of the followed methodology; the solid oscillation was accurately predicted over the testing time interval, while the average velocity error over the domain was found to be less than 3%.Finally, a parametric study with respect to the ROM dimension showcased the increased ROM robustness offered by the inference of sparse, full-order operators.The current work could comprise a first step towards a complete non-intrusive MOR framework for FSI problems.Hence, future work could incorporate several different aspects; in particular, the offline computational efficiency of the approach could be significantly enhanced through parallelizing least-squares problems and domain segregation.Moreover, testing the developed non-intrusive ROMs under varying inputs and solid oscillation parameters would provide further insight on the attributes of the proposed methodology.In a similar direction, the showcased sparse full-order model independence from the SVD basis renders parametric MOR specifically interesting, aiming towards non-smooth parameter dependencies.Such studies could be performed for VIV problems with respect to different Reynolds numbers, given the recorded rich coupled system dynamics and strong parametric dependence.Finally, the extension to deformable solids and consequently the formulation of a predictive, non-intrusive ROM methodology for FSI problems lies in the future scope of this work.

Figure 1 :
Figure 1: Problem schematic representation:A nondeformable body with translational degrees of freedom along x, y, subject to an incompressible channel flow.

Figure 2 :
Figure 2: Singular value decomposition of 3000 × 500 data matrices of the u x , u y velocity components of a 2D flowfield for a VIV case.The first singular value of u x corresponds to the mean flow.

2 and x 2 Figure 3 : 2 , ˆ x 2 2
Figure 3: Example of L 2 regularization L-curve criterion: The model (and thus the λ value) which gives the smallest normalized distance from the ˆ b 2 L-curve for the DOF-wise and interpolated λ models: Slightly increased error for the interpolated λ model, at a significantly lower cost.Optimal λ computed at each grid node (ux velocity).Optimal λ computed at 10% of grid nodes, interpolated for the rest 90%.

Figure 4 :
Figure 4: Regularization parameter selection for fluid nodes: Indicatory example for a Re = 180 dataset (discussed in detail in Section 4).Increased λ values in the flow wake for both the exact and the interpolated optimal regularization strategies.

Figure 5 :
Figure 5: Body-proximity nodes: Given a distance threshold, only 2n F DOFs are included in the body forces non-intrusive model k x (N/m) k y (N/m) ūin (m/s

Figure 6 :
Figure 6: First POD mode for Re = 90: The first POD mode pair contains more than 50 % of the total kinetic energy of the flow (here only the first mode is displayed) and is linked with the 2S vortex shedding frequency (peak at f = 1.72 Hz).The values of Φ 1 are normalized to a [−1, 1] scale for illustration purposes.

Figure 7 :
Figure 7: Fifth POD mode for Re = 90: A slow-drift of the mean flow is captured, linked to the transient wake dynamics.Minor peaks at f v and f N , possibly also due to transient phenomena.

Figure 8 :
Figure 8: Fifth POD mode for Re = 180: Solid natural frequency affects the wake dynamics, since the case lies in the right side of the initial VIV branch.A significant peak at 3.12 Hz, very close to the theoretical f N value is observed in the Fourier spectrum of Ψ 5 .

Figure 9 :
Figure 9: Solid oscillations predictions: The coupled dynamics of the data-driven flow ROM (27) and first-principle oscillation equations (28) yield highly accurate predictions for the solid oscillations, capturing well the low-amplitude motion (Re = 90) as well as beating phenomena (Re = 180), even past the training time.

Figure 10 :
Figure 10: Contour plots at t = 5.9s for the u x flowfield at Re = 90: The non-intrusive model captures the 2S vortical mode as well as the flow features close to the oscillating body.

Figure 11 :
Figure 11: Contour plots at t = 5.9s for the u y flowfield at Re = 90: Slight phase mismatch, also indicated by the solid transverse oscillation prediction in Figure 9b.

Figure 12 :
Figure 12: Contour plots at t = 5.9s for the u x flowfield at Re = 180: Accurate ROM prediction, with a double vortex shedding frequency compared to Re = 90.

Figure 13 :
Figure 13: Contour plots at t = 5.9s for the u y flowfield at Re = 180: The non-intrusive model is again in good accordance with the corresponding CFD results at the end of testing time.

Figure 14 :
Figure 14: Average error over the flowfield, for both u x and u y relative to max u x , max u y over time.

Figure 15 :
Figure 15: Average error (29) over testing time, for different ROM dimensions r: A projection subspace dimension of[15,20] is sufficient for the ROM, with an observed error convergence for a further r increase.