Adjacency-based, Non-intrusive Reduced-order Modeling for Fluid-Structure Interactions

Non-intrusive model reduction is a promising solution to system dynamics prediction, especially in cases where data are collected from experimental campaigns or proprietary software simulations. In this work, we present a method for non-intrusive model reduction applied to Fluid-Structure Interaction (FSI) problems. The approach is based on the a priori known sparsity of the full-order system operators, which is dictated by grid adjacency information. In order to enforce this type of sparsity, we solve a local, regularized least-squares problem for each degree of freedom on a grid, considering only the training data from adjacent degrees of freedom, thus making computation and storage of the inferred full-order operators feasible. After constructing the non-intrusive, sparse full-order model, Proper Orthogonal Decomposition (POD) is used for its projection to a reduced dimension subspace and thus the construction of a reduced-order model (ROM). The methodology is applied to the challenging Hron-Turek benchmark FSI3, for Re = 200. A physics-informed, non-intrusive ROM is constructed to predict the two-way coupled dynamics of a solid with a deformable, slender tail, subject to an incompressible, laminar flow. Results considering the accuracy and predictive capabilities of the inferred reduced models are discussed.


Introduction
Fluid-Structure Interactions encompass coupled, strongly nonlinear fluid-solid dynamics systems, which arise in various engineering fields.For different regimes, materials and geometry configurations, the specific coupling conditions and corresponding FSI dynamical response can widely differ [1].In this work, we focus on incompressible fluid dynamics at low Reynolds numbers, coupled with deformable solids, at relatively low mass ratios.Under such conditions, the timescales and inertia of the two subsystems are comparable, leading to a strong, two-way coupling mechanism.
The development of non-intrusive, reduced-order models for FSI problems is motivated by needs in engineering design and control.Firstly, the high computational cost of direct FSI numerical simulation renders multi-query tasks such as design optimization and real-time control cumbersome.Such engineering tasks could benefit from model reduction techniques, which have proven to allow for considerable speed-ups of FSI simulation, without significant loss of accuracy [2].Secondly, in cases where experimental or FSI simulation data are available, without access to a numerical simulation source code, non-intrusive modeling methods are developed for system dynamics prediction [3].Exploiting knowledge on the structure of the underlying dynamical system has been shown to enhance the predictive capabilities of non-intrusive ROMs for specific FSI regimes such as aeroelasticity [4].
In this work, we implement an adjacency-based approach to FSI problems.The method was recently studied for linear advection and diffusion problems in [5] and [6] and has been applied to the inference of quadratic-bilinear, incompressible fluid dynamics operators in Vortex-Induced Vibrations problems by the authors [7].Two coupled, data-driven models are inferred for the fluid dynamics and solid dynamics subsystems, exploiting a physics-based structure of the governing equations under the Arbitrary Lagrangian-Eulerian (ALE) formulation.We present results for the FSI3 Hron-Turek benchmark [8], which has been previously studied for intrusive model reduction in [2].Through this testcase, we examine the predictive capabilities of the obtained, non-intrusive model.

Governing equations
For the fluid flow, we focus on the regime of 2D incompressible, laminar flows.In cases of moving fluid domains Ω(t), the equations can be expressed with respect to a reference configuration Ω using the ALE formulation.The map from the reference domain to the current configuration is given as: x → x = x + df (t), where df (t) is a deformation field.The incompressible ALE Navier-Stokes equations are then given by where û is the ALE velocity field, p is the pressure, σ is the ALE stress tensor given by F is the gradient of the ALE map and J = det F, with For solid dynamics modelling, we use the linear Navier-Lamé equations, assuming no material or geometrical nonlinearities.The boundary of the solid is split into Γ s 1 where the displacement is zero and Γ s 2 where the fluid-structure coupling occurs.The solid deformation d s is thus modeled by The boundary conditions of the coupled fluid-solid problem are given below, based on Figure 1.
The above formulation considers the coupled problem on a reference fluid domain Ω.We use a stiffened harmonic extension as an ALE map, to transfer the solution to the current, deformed configuration Ω(t).
with boundary conditions on the interface Ω ∩ S, and on the domain external boundaries: a(x) introduces a stiffening effect close to the FSI interface Ω ∩ S, which allows for larger mesh displacements without element degradation, compared to a purely harmonic extension (see Section 5.3.5 of [1]).t 1 is selected as a reference time, for which Ω = Ω(t 1 ).

Approximate model structure
By discretizing (1), ( 4) and ( 6), we aim to reveal an approximate structure of the dynamical system, which will be exploited for non-intrusive modeling.We assume that the gradient of the deformation field ∇d f is sufficiently small, such that F ≈ I (and thus J ≈ 1) is a valid approximation of (2).
In that case, the structure of ( 1) is approximately quadratic-bilinear.After discretizing in space and time, working similarly to [9], we cancel out the algebraic contribution of the pressure and come up with the following, approximate model structure for the fluid velocity field: where (u) 2 denotes the vector of 2n f 2 unique elements of the Kronecker product u ⊗ u.For the solid velocity field, we get a second-order, forced linear system by discretizing (4) in space and time The forcing f originates from the dynamic coupling condition on the interface (i.e. the projection of the fluid stress tensor on the normal to the interface).This can be approximated by a quadratic model, as follows.
where F S denotes nodes that are sufficiently close to the interface, given the sparsity of the discretized stress tensor (see (2)).This means that u F S = Ru, where R is an ordered selection matrix.Finally, the solid velocity field ∂ t d s should be integrated to update the solid displacement field d s .To this end, a Crank-Nicholson scheme is used We also need to discretize and solve (6), in order to compute the fluid mesh displacement in the current configuration, Ω(t), with respect to the reference configuration Ω.This leads to where A L is the inverse of the discretized ALE map and S is an ordered selection matrix, which picks only the solid DoFs on Γ s 2 .Equations ( 8) to (12) comprise a discretized, coupled dynamical system, which models ALE FSI problems under small deformations.In the following, we use this structure to fit the corresponding operators to numerical FSI data and thus construct a physics-based, non-intrusive model.

Physics-based adjacency
Having identified an approximate structure for the FSI dynamical system, we aim to infer the unknown operators of the fluid dynamics subsystem ((8)), the solid dynamics subsystem ((9)) and the solid forcing term ((10)), from simulation data.
The discretized operators in (9) are sparse, while the operators of ( 8) and (10) are only approximately sparse, due to canceling out the contribution from the pressure.The corresponding sparsity pattern originates from grid adjacency.Based on this observation, the employed method aims to infer fullorder, adjacency-based, sparse operators.
To proceed with the inference task, we should map the available data for the flowfield over time t = [0, T ] to a reference configuration Ω, since equation ( 8) was derived there.Knowing the motion of the FSI interface from the data, we define a reference time as t 1 = argmin( ds (t) 2 ).Based on the reference position of the FSI interface at t 1 , we can construct a fluid and a solid mesh (e.g. via Delaunay triangulation), with the constraint that the nodes of the two meshes on the interface should match.
For the fluid, we compute the deformed mesh Ω(t) from (12), given the motion of the FSI interface (see (7)) and the constructed reference Ω = Ω(t 1 ).Then, the fluid velocity data can be interpolated on the reference configuration Ω at each of the n T timesteps.For the solid, we interpolate the displacement field values for all n T timesteps from the solver mesh configuration at t 1 , to the constructed solid mesh.The difference between handling the solid and the fluid problem should be noted: Mesh construction is necessary only for the fluid side due to the ALE formulation.However, grid adjacency information for both the fluid and the solid mesh is essential for the following.Thus, the solid mesh of the simulation data can be readily used as long as adjacency information can be retrieved.

Least Squares formulation and regularization
The number of grid nodes is denoted by n f for the fluid and n s for the solid mesh.Then, there are 2n f and 2n s kinematic degrees of freedom (DoFs), respectively, for this 2D problem.For the i-th DoF (streamwise or transverse component of u or ∂ t d), we denote by q A (i) the set of geometrically adjacent DoFs in the fluid A ≡ f or solid A ≡ s mesh respectively.For example, if i is a fluid velocity component, q f (i) includes exclusively DoFs which correspond to fluid mesh nodes, adjacent to the node of i.Then, the local inference of the operators for DoF i corresponds to solving the following problem: min where we make a distinction between the fluid and solid dynamics problems: For the fluid dynamics subsystem, where T 1 corresponds to the training time and m = [1, ..., n T1 − 1].For fluid DoFs on Γ f 1 and Ω ∩ S we can directly impose the corresponding Dirichlet boundary conditions, by setting row i of A, H, K, C to zero, and assigning entries of value 1 in corresponding positions of L and B, accordingly.For the solid dynamics subsystem, we introduce sets q F S (i) and q I (i), to distinguish between internal DoFs and DoFs on Γ s 2 .As a result, q F S (i) = q f (i) for DoFs i on Γ s 2 and q F S (i) = ∅ otherwise, while q I (i) = q s (i) everywhere, except for Γ s 2 , where q I (i) = ∅.Then the variables of (13) for the local inference of the solid dynamics operators are In practice, (13) should be complemented with a regularization term, to bypass potential numerical errors due to small singular values of D i .To this end, we employ truncated SVD of D i , for different truncation thresholds η of the normalized singular values of D i .The optimal η is determined as the "corner" of an L-curve, with the solution norm β i 2 on the x axis and the error β T i D i − y m+1 i 2 on the y axis (see Figure 2).In essence, the SVD truncation corresponds to an L 2 regularization of (13).
Regularization plays a crucial role on the inferred model properties and thus determining the range of possible η is not trivial.The process of computing the optimal threshold η is performed for few, randomly sampled DoFs.An interpolation from the sampled, optimal η values is performed to the rest of the mesh DoFs, as shown in Figure 3, such that only one least-squares problem per i is solved for those, with a predefined truncation threshold.

Model reduction and coupling
After solving (13) for both the solid and the fluid nodes, we have inferred the sparse, full-order matrices of ( 8) and (9).In order to derive a reduced-order model, we use a POD projection.This is done by taking the SVD of the velocity flowfield u snapshot matrix, and of the solid velocity field ∂ t d s snapshot matrix, over some training interval [0, T 1 ].By retaining the first r f and r s singular modes accordingly, we obtain the orthonormal projection bases Φ f ∈ R 2n f ×r f and Φ s ∈ R 2ns×rs .Substituting the projected states u = Φ f ũ and ∂ t d s = Φ f ∂ t ds to (8) and ( 9) results to a coupled, non-intrusive, and the so- Figure 3: Interpolated, sub-optimal η value for the streamwise flowfield component, using the optimal η for 30 % of the DoFs.Regularization values are higher along the wake of the flow and close to the deformable solid tail.
reduced-order model The corresponding ROM matrices of ( 16) can be efficiently computed since the sparsity patterns of the full-order matrices are known.We then have for the fluid subsystem and for the solid It was observed that the stiff coupling exhibited in FSI numerical simulation [8] transfers also to the data-driven, reduced-order model.The implicit formulation of ( 8) is necessary for stability.The combined convergence criterion for each timestep ũk j+1 − ũk j 2 < res 1 and ds < res 2 was used during simulation, where j denotes the index of iterations within timestep k.For each implicit iteration, a successive under-relaxation method was used.

Numerical Results
The analyzed methodology was tested for the Hron-Turek FSI3 benchmark [8].Simulation data were obtained using the Gascoigne 3D open-source, Finite Element solver [10].A monolithic approach was used, with a moderate mesh of 4128 elements and a timestep of ∆t = 5 × 10 −3 s, for a total time of 5 s.We discard the first 1.5 s of data to avoid the influence of numerical solution initialization.The clock time required for the numerical simulation is approximately 1.5 hours on a laptop.
As presented in Section 3.1, we interpolate the data to a new mesh with 4060 nodes, for which we construct the ALE map (12).We use 65% of the overall available time series (T 1 = 2.3 s) to solve (13), computing the optimal SVD truncation value for 10% of the mesh nodes, using 10, logarithmically sorted values of η ∈ [10 −5 , 10 −1 ].We then construct the ROM through (17), (18), by selecting r f = 40, r s = 5, and solve the implicitly coupled, reduced-order system (16) by successive under-relaxation and res 1 = res 2 = 10 −6 .The "offline" step of mesh construction, interpolation and FOM inference requires approximately 5 minutes, while the simulation time of the ROM requires less than a minute.
The predicted, reconstructed ROM flowfield at T = 3.5 s, well beyond the training time T 1 , is given in Figure 4, in comparison to the CFD results for the same instance, for both velocity components.It is observed that the main dynamical features of vortex shedding as well as the motion of the deformable tail are well approximated, using only 45 reduced variables, from the initial 8120 DoFs.The vertical velocity over time for the upper right corner of the deformable tail is also presented in Figure 5.The predicted ROM solid motion is well in line with the numerical results and captures the limit cycle behaviour beyond training time T 1 .A minor phase shift can be attributed to the considerable presence of transient dynamics in the training data.Finally, it is interesting to investigate the stability of the derived FSI model with respect to the ROM dimensions r s and r f of the solid and fluid subsystems, respectively.In Figure 6, we illustrate the ROM streamwise velocity error with respect to the numerical simulation on the reference configuration Ω, averaged over the domain Ω and testing time t = [T 1 , T ].This is denoted with e x (t).Based on e x (t), we observe that r s ≥ 3 modes are required for the solid subsystem for an accurate FSI ROM, as anticipated by the deformable tail motion.For values r s ≥ 8, a coupling instability was exhibited, irrespective of the r f value.Similarly r f ≥ 20 SVD modes are necessary on the fluid side.However, we also observe that the combination of r f and r s can significantly affect the ROM predictive accuracy, due to the strong coupling on the FSI interface.For a significant range of ROM dimensions, the average error e x (t) is below 5%.

Conclusions
In this work, we employed an adjacency-based method for non-intrusive, reduced-order modelling of FSI problems.Accurate results were obtained for the Hron-Turek FSI3 benchmark testcase, with an average ROM prediction error of less than 5 %.Based on the potential of this framework, aspects of computational efficiency, parametric model reduction, as well as as a theoretical analysis of local non-intrusive inference comprise potential, future research directions.

Figure 1 :
Figure 1: Schematic representation (not scaled) of the Hron-Turek FSI benchmark testcase.A velocity profile is assigned at the inlet Γ f 1 and FSI coupling occurs on the surface Γ s 2 of a slender, deformable tail.

Figure 2 :
Figure 2: Left: Truncation of the singular values of D to threshold η, for a random DoF i. Right: Corresponding L-curve: The optimal η achieves a compromise between the leastsquares error β T i D i − y m+1 i 2

Figure 4 :
Figure 4: (a) CFD results at t = 3.5 s for streamwise velocity component (b) ROM results (r f = 40, r s = 5) at t = 3.5 s for streamwise velocity component (c) CFD results at t = 3.5 s for transverse velocity component (d) ROM results (r f = 40, r s = 5) at t = 3.5 s for transverse velocity component

Figure 5 :
Figure 5: Velocity ROM prediction (streamwise, transverse) for the upper right corner of the deformable solid.Training time T 1 = 2.3 s is denoted with a dashed vertical line.

Figure 6 :
Figure 6: Average error for streamwise flowfield over testing time e x (t), for different combinations of ROM dimensions (r f , r s ).