Adaptive Reduced-Order Modeling for Non-Linear Fluid-Structure Interaction

We present an adaptive reduced-order model for the efficient time-resolved simulation of fluid-structure interaction problems with complex and non-linear deformations. The model is based on repeated linearizations of the structural balance equations. Upon each linearization step, the number of unknowns is strongly decreased by using modal reduction, which leads to a substantial gain in computational efficiency. Through adaptive re-calibration and truncation augmentation whenever a non-dimensional deformation threshold is exceeded, we ensure that the reduced modal basis maintains arbitrary accuracy for small and large deformations. Our novel model is embedded into a partitioned, loosely coupled finite volume - finite element framework, in which the structural interface motion within the Eulerian fluid solver is accounted for by a conservative cut-element immersed-boundary method. Applications to the aeroelastic instability of a flat plate at supersonic speeds, to an elastic panel placed within a shock tube, and to the shock induced buckling of an inflated thin semi-sphere demonstrate the efficiency and accuracy of the method.


Introduction
Fluid-Structure Interaction (FSI) occurs in a broad range of applications, such as blood flow though heart valves (Peskin, 1972), flutter of aircraft wings (Farhat and Lesoinne, 2000) and shock-induced deformations of rocket nozzles and panels (Garelli et al., 2010;Pasquariello et al., 2015). FSI simulations involve two different branches of computational physics: Computational Fluid Dynamics (CFD), which is often based on an Eulerian finite-volume representation, and Computational Solid Mechanics (CSM), for which a finite-element discretization is frequently chosen.
FSI algorithms can be divided into monolithic and partitioned methods. The monolithic approach is characterized by solving a single set of discrete equations describing entire coupled system (Klöppel et al., 2011). While this procedure may be time-consuming, it is robust, accurate and stable. On the other hand, a partitioned approach is frequently used to conveniently couple off-the-shelf CFD and CSM codes, which are extensively validated and employ the most efficient numerical schemes for each discipline. Partitioned methods can be further classified as strongly or loosely coupled, where the distinction is based upon whether or not the complete set of coupling conditions at the conjoined FSI interface is satisfied. For similar mass density of fluid and solid, loosely coupled methods may suffer from the artificial added-mass effect, possibly leading to computational instabilities (Farhat and Lesoinne, 2000;Causin et al., 2005). Stability can be recovered by introducing sub-iterations (Küttler and Wall, 2008), which however increase the computational cost significantly (Farhat et al., 2006). Badia et al. (2008) obtained very promising results by employing a Robin-type boundary condition at the FSI interface. Similarly, Banks et al. (2013Banks et al. ( , 2014 introduced so-called Added-Mass Partitioned algorithms to overcome the added-mass effect for incompressible flow as well as for very light rigid bodies in compressible flow. For the majority of FSI applications, such as compressible aeroelasticity, however, a loosely coupled method is sufficient (Farhat and Lesoinne, 2000).
In practice, one of the main challenges is to keep the computational cost of time-resolved FSI simulations at a reasonable level without sacrificing the required accuracy. The present paper addresses the main performance bottleneck of high-fidelity turbulence resolving FSI simulations with two loosely coupled domain-specific codes, where the time step size is restricted by the resolution requirements of the fluid flow: CFD solvers for Large Eddy Simulations (LES) and Direct Numerical Simulations (DNS) usually employ high-order finite-volume or finite-difference space discretizations and explicit time marching methods, which efficiently satisfy the high resolution requirements of the fluid flow, and can be parallelized easily with excellent scalability on massively-parallel supercomputers. The parallelization of finite-element CSM methods is less straight-forward; in our experience, most parallelized CSM solvers require a shared memory architecture and usually show limited parallel scalability. This leads to the curious situation that one time step for advancing the orders of magnitude more expensive (in terms of CPU time) CFD equations requires much less run time (wall-clock time) than one time step of the CSM problem (Pasquariello et al., , 2016. This does not mean that finite-element CSM methods are per se inefficient; they are designed for much larger time steps, which would be effective in a stand-alone setting. Piperno et al. (1995) addressed this bottleneck with a sub-cycling approach, where the fluid solver is advanced multiple times before the structural solution is advanced in one large time step. The efficiency and stability of various sub-cycling methods are discussed by Farhat et al. (1997).
As an alternative, Reduced-Order Modeling (ROM) can be used to improve the efficiency of the CSM solver (Dowell, 1996;Hall et al., 1999;Rixen, 2001). One of the first reported model reductions was presented by Rayleigh (1887), who employed the Mode Superposition Method (MSM) to approximate the displacement field with a low number of free vibration modes. The method truncates the vibration modes of the structure at low number, i.e. N << n, where n is the number of degrees of freedom and N is the number of dynamically important modes. Several improvements for truncated modal superpositions have been proposed since then, one being the Static Correction Method (SCM) used by Rixen (2001), Wilson et al. (1982) and Besselink et al. (2013), among others. This method accounts for the omitted modes by including the truncated modes statically, which leads to a more adequate representation of the modal loads. However, this method is only effective if the structure has very low natural frequencies. Dickens et al. (1995) introduced the Modal Truncation Augmentation (MTA) method for the dynamic correction of the load representation.
MTA improves the overall accuracy compared to the SCM and is also effective for a broader frequency range (Dickens et al., 1995). Another popular reduction method is the Craig-Bampton Method (CBM) presented in Bampton and Craig (1968) and the family of Component Mode Synthesis (CMS), see, e.g., Qu (2010). These methods divide the global finite-element structure into several sub-structures connected with an appropriate interface description. The CMS method is also known as Super Element Method in the sense that each substructure can be considered as a single finite element.
Linear ROM approaches generally fail for non-linear problems. Some of the earliest reported work on non-linear ROM was presented by Morris (1977). It follows a modal superposition that incorporates the stiffness matrix at a deformed state and was applied for frame structures. Similarly, Remseth (1979) employed a Taylor linearization of the reduced-order finite-element equations with respect to a deformed state to account for non-linear effects on frame structures. The projection on a new eigenmode basis (around a new deformed state) generates truncation errors, which can accumulate. Nickell (1976) used a Rayleigh-Ritz analysis to derive a reduction model that includes non-linear effects through modal derivates along with the eigenmodes and exploits those derivates along with tangent modes to limit the need of eigenspace updates. Alternative methods that minimize the need to update the reduction basis have been presented by Idelsohn and Cardona (1985) & Tiso and Rixen (2011). Mignolet et al. (2013) reviewed ROMs for non-linear geometric structures based on indirect methodologies, where the non-linear stiffness terms are approximated by cubic polynomials. A key aspect of the ROM effort is to properly select the basis functions. The authors present a strategy that enriches the basis of free vibration modes by dual modes for capturing non-linear structural behavior. A set of non-linear static simulations with representative loads is needed to determine the dual modes, which are calculated based on the proper orthogonal decomposition of the series of non-linear displacement fields. Recently, Wu and Tiso (2016) presented a ROM for flexible multi-body systems with large non-linear deflections, where the ROM is based on a combination of the CBM and cubic polynomials of the configuration dependent terms.
In this paper we derive an adaptive ROM (AROM) based on adaptive re-calibration of the reduced modal basis with MTA correction, which allows us to maintain arbitrary accuracy also in the case of large and non-linear structural deformations. The AROM is imbedded into the loosely coupled partitioned FSI algorithm proposed by Pasquariello et al. (2016). We employ an established high-fidelity turbulence resolving finite-volume method for solving the threedimensional compressible Navier-Stokes equations on block-structured adaptive Cartesian grids (INCA, https:// inca-cfd.org) and an unstructured finite-element method for the discretization of the structural domain (CalculiX, http://www.calculix.de). The time-varying fluid-structure interface is accounted for by the cut-element based 3 Immersed Boundary Method (IBM) that was introduced byÖrley et al. (2015) and then extended to deformable structures and compressible FSI applications by Pasquariello et al. (2016). The CFD solver is designed to provide (close to) ideal parallel scalability on thousands of cores (Egerer et al., 2015). Due the small time step size imposed by the resolution requirements of compressible fluid flows and the necessarily synchronous nature of alternating between the CFD and CSM at each time step, the CSM solver constitutes the critical run-time bottleneck in this framework. The essential original contribution of this work is the development and demonstration of a novel FSI-AROM algorithm, which is capable of handling structures with large, non-linear deformations accurately with high computational efficiency.
The paper is structured as follows: The governing equations of the fluid and structure and the numerical formulation are presented in Sections 2 and 3. The coupling algorithm for non-matching time-varying interfaces is presented in Section 4. The novel FSI-AROM method is derived in Section 5. In Section 6 we validate the FSI-AROM method for linear and non-linear problems, and demonstrate the prediction capabilities for flow-induced buckling of a threedimensional thin semi-spherical shell. A final discussion and concluding remarks are given in Section 7.

Governing Equations
The domain of interest Ω = Ω F ∪ Ω S is divided into non-overlapping fluid Ω F and solid Ω S subdomains with a conjoined interface Γ = Ω F ∩ Ω S . The interface normal vector n Γ is assumed to point into the fluid domain. In the following a brief description of the mathematical models required for both subdomains is given. Unless specified otherwise, we use the Einstein summation convention.

Fluid
The fluid flow within the domain Ω F is governed by the three-dimensional compressible Navier-Stokes equations which describe the conservation of mass, linear momentum and total energy. We use Cartesian coordinates where The state vector w and flux tensor H (w) = H (1) , H (2) , H (3) are given as where u i is the velocity, ρ F the fluid density, and ρ F e t is the total energy density. The viscous stress tensor for a Newtonian fluid is where the first Lamé parameter is related to the dynamic viscosity µ F according to Stoke's hypothesis: The heat flux is evaluated according to Fourier's law, with the coefficient of thermal conductivity k. We consider air as a perfect gas with γ = 1.4 and specific gas constant of R = 287.058 J kg·K . The pressure p and temperature T are calculated from the definition of total energy and the ideal-gas equation of state

Solid
The governing equations for the solid are based on the local form of the balance of linear momentum which describes an equilibrium between the work done by inertia, internal and external forces expressed in the underformed configuration. The vector of displacements is denoted by d, ρ S ;0 is the material density of the solid, ∇ 0 · ( ) is the material divergence operator, P = F · S is the first Piola-Kirchhoff stress tensor, where F is the deformation gradient tensor, and external material body forces are represented byb 0 . The second Piola-Kirchhoff stress tensor is In this work, a hyperelastic Saint Venant-Kirchhoff material model is chosen. Its associated strain energy density function Ψ is given as where λ S and µ S are the first and second Lamé parameter. The Green-Lagrange strain tensor is defined as and D is the displacement gradient tensor. The Cauchy stress tensor σ S , also called true stress tensor, is defined as where J denotes the Jacobian determinant.
Boundary conditions need to be specified on ∂Ω S = Γ S ,D ∪ Γ S ,N ∪ Γ to make the system (7) well posed. Two different types are considered in this work, namely Dirichlet Γ S ,D and Neumann Γ S ,N boundaries for which we either prescribe displacementsd or tractionst d =d on Γ S ,D and P · n 0 =t on Γ S ,N .
Here n 0 denotes the unit normal vector in material configuration. Further, initial conditions for displacements and velocities must be specified:

Fluid-structure interface conditions
The interface between fluid and structure requires coupling conditions. Tractions on Γ have to be in equilibrium, that is, Herein, σ S is the Cauchy stress tensor given by Eq. (11) and denotes the fluid stress tensor comprising an inviscid and viscous contribution. In addition, the kinematic no-slip boundary condition must be satisfied, which in case of an inviscid flow reduces to matching normal velocities on Γ 3. Numerical models

Fluid
The compressible Navier-Stokes equations, Eq. (1), are discretized with a finite volume method based on the integral form where Gauss's theorem has been applied. The integral is taken over Ω i, j,k ∩ Ω F , i.e., the part of a Cartesian computational cell Ω i, j,k that belongs to the fluid domain Ω F , and over the time step ∆t = t n+1 − t n . In order to account for the FSI interface within the fluid solver, which operates on Cartesian grids, we employ the cut-element IBM ofÖrley et al.
S Ω (e) S Figure 1: Schematic triangulation of a structural interface element Γ (e) S . Resulting interface triangles Γ tri are used as an input for the cut-algorithm to compute individual cut-elements Γ ele and cut-cell related geometric quantities.
(2015) and Pasquariello et al. (2016). The discrete FSI interface is composed of several structural interface elements S . Each structural interface element Γ (e) S is triangulated as shown exemplarily in Fig. 1 for a quadratic hexahedral element. The resulting interface triangles Γ tri are used as an input for the IBM algorithm. A fluid cell that is cut by at least one interface triangle Γ tri is referred to as a cut-cell. The fluid-solid interface within a cut-cell is composed of one or several cut-elements Γ ele = Γ tri ∩ Ω i, j,k , each representing one or a part of one interface triangle, see Fig. 1.
Applying a volume average of the conserved variables and considering (for demonstration purposes) a simple forward Euler time integration scheme leads to the following discrete form of Eq. (18) Herein α i, j,k is the fluid volume fraction of the cut-cell, V i, j,k = ∆x i ∆y j ∆z k is the total volume of cell Ω i, j,k and A is the effective fluid wetted cell-face aperture, see also For the spatial reconstruction and numerical flux functions we either use the Adaptive Local Deconvolution Method (ALDM) by Hickel et al. (2006Hickel et al. ( , 2014, or the 5 th -order WENO (Weighted Essentially Non-Oscillatory) scheme by Liu et al. (1994) with the HLLC flux (Toro et al., 1994). In order to avoid modified interpolation stencils in the FV reconstruction near the interface, we assign special ghost fluid states that depend on the interface boundary conditions to non-cut fluid cells within the solid part of the domain (Mittal et al., 2008;Pasquariello et al., 2016).
Finally, time integration is performed with a conditionally stable, explicit third-order Runge-Kutta scheme.

Solid
We cast the structural equations, Eq. (7), into their weak form by applying the principle of virtual work with virtual displacements δd and subsequently integrating the balance equation over the structural subdomain Ω S . Following this procedure and applying Gauss's theorem yields where dA 0 and dV 0 are infinitesimal surface and volume elements, respectively, and δE is a result of the variation of the strain expression in Eq. (10), The weak form, Eq. (21), represents the balance of virtual work δW, namely where the work at the FSI interface is δW Γ S . We use the FEM to discretize the integral equation (21) in space. The structural domain Ω S is composed of n e solid elements Ω (e) S with consistent basis functions for representing the displacement field. The semi-discrete form of Eq. (21) is then obtained by assembling the contributions of all elements with the mass matrix M, the discrete acceleration vectord and the discrete displacement vector d. The forces are divided into internal forces f S ;int , external forces f S ;ext and interface forces f Γ S resulting from the fluid. In contrast to Pasquariello et al. (2016), who used linear FE together with element technology based on Enhanced Assumed Strains to avoid locking phenomena, we use quadratic shape functions for interpolating the displacements on Ω (e) S unless stated otherwise. The final step is to discretize Eq. (24) in time. We employ the Hilber-Hughes-Taylor α-method (Hilber et al., 1977) for time integration. Due to its implicit character a coupled set of non-linear equations needs to be solved, which is done by the Newton-Raphson method.

Load and motion transfer
The cut-cell discretization inevitably leads to non-matching grids at the conjoined interface Γ and requires interpolation methods for the load transfer between both subdomains. Specifically, we search for the discrete force vector f Γ S that results from the fluid tractions acting on the wetted structure. We follow the approach suggested by Farhat et al. (1998) and use the shape functions of Ω (e) S for interpolating the fluid loads on the adjacent structural nodes. Consider a single cut-element Γ ele as shown in Fig. 1. The fluid forces f ele F follow directly from the pressure and viscous contributions to the IBM interface flux χ = ele χ ele in Eq. (20). Since the structural interface, i.e., the triangulation Γ tri , directly serves as an input for the IBM, there is no extra need for a pairing algorithm to associate a single face centroid x c ele to the closest wet structural interface segment Γ (e) S . However, we have to determine the natural coordinates ξ c ele (x c ele ) of this fluid point. This inverse mapping problem is solved iteratively with a Newton-Raphson method. Finally, the force contribution from a single cut-element to an individual node of the paired structural interface segment where N k denotes the shape function of the k−th structural node on Γ (e) S . Summing up the contributions of all cutelements in Ω F leads to the interface force vector f Γ S . It is easy to verify that this interpolation guarantees a global conservation of loads over the interface by recalling that all shape functions at one specific location sum up to unity.
The cut-element IBM requires the velocity at the face centroid x c ele for evaluating the work done at the interface. We use the same interpolation strategy based on the shape functions of the structural domain whereḋ k is the velocity of the k−th structural node on Γ (e) S . The motion of the structure within the fluid domain is accounted for by updating the cut-elements (and cut-cells) after each time step based on the triangulated interface Γ tri (Pasquariello et al., 2016). Consequently the compatibility between the displacement fields of the structure and the fluid at the FSI interface is implicitly fulfilled in a discrete sense for all structural nodes k ∈ Γ S and no further interpolation is required.

Summary of the coupling procedure
We use an explicit, first-order in time accurate, loosely coupled FSI algorithm to advance the system from time level t n to t n+1 = t n + ∆t n . The main steps are summarized below: 1. At time level t n the structural displacements d Γ;n and velocitiesḋ Γ;n at the interface are used to update the cut-cell geometry. The triangulated interface Γ tri , see Fig. 1, is used as an input for the cut-algorithm.

Linearization and modal truncation
In this section, we propose a numerical framework for switching between a full FEM description and a more efficient Adaptive Reduced-Order Model (AROM) that maintains accuracy also when a structure undergoes large non-linear deformations. The algorithm is based on Taylor expansion around a reference state d re f . Linearizing

Eq. (24) around this reference leads to
Rearranging Eq. (27) leads to where the tangent stiffness matrix K (d ref ) represents the Jacobian of the internal forces The initial conditions for δd are where the superscript n denotes the (last) results obtained with the full FEM model, Eq. (24), before switching to AROM. Since this initial condition is also considered as the reference state, i.e., d n = d ref , the initial condition for the deflections δd 0 is zero.
Equation (29) is expressed in the physical space. For reduced-order modeling we shrink the system of equations by the mode superposition method (Rayleigh, 1887;Dowell and Hall, 2001). In a first step, the eigenmodes of the structure are obtained by the following general eigenvalue problem of order m where the columns of Φ = φ 1 , . . . , φ m are the orthonormalized (with respect to M) eigenvectors (natural vibration modes) and Ω = diag (ω 1 , . . . , ω m ) is a diagonal matrix listing associated eigenvalues (natural vibration frequencies).
Note that Eq. (34) is only exact when the sizes of K and Φ are equal. We define the following transformation from modal to physical space where δq denotes the vector of perturbations expressed in generalized coordinates, i.e. modal amplitudes. Substituting the latter expression into Eq. (29) and left-multiplying all terms with Φ T leads to The size of the generalized matrices M G and K G directly depends on the number of eigenvectors considered, i.e., including the first N eig eigenmodes reduces the system to rank N eig . Furthermore, the principle of orthogonality implies that Eq. (36) can be written for the i−th mode as recalling that M G is a unit matrix and K G is a diagonal matrix with eigenfrequencies squared on the diagonal (Vahdati et al., 1999). An unconditionally stable Newmark scheme is used for time integration of the modal equations with the following initial conditions prescribed in modal space Equations (40) and (41) are derived using the orthogonality principle, i.e. Φ T MΦ = I, and the relation in Eq. (35).

Modal truncation augmentation
The Modal Truncation Augmentation (MTA) method was derived by Dickens et al. (1995) in order to improve the representation of the load vector in modal space. The generalized loads can be computed as where f S ;tot = f S ;ext + f Γ S is the total load vector including external and interface loads. We transform the generalized forces back to the physical domain byf which consequently results in a projection error that can be summarized in a residual The MTA method attempts to correct for the projection error by appending a pseudo eigenvectorφ to the original modal basis Φ. Note that the pseudo eigenvector does not satisfy the eigenvalue problem defined in Eq. (34), but it satisfies the orthogonality principle (Dickens et al., 1995). In a first step we solve for the displacements d cor due to the residual force vector Following this, we compute where K cor and M cor are the stiffness and mass matrices projected with respect to the displacement vector d cor . Note that in the special case of a single right-hand-side vector the matrices, K cor and M cor reduce to simple scalars and the following trivial eigenvalue problem can be formulated where φ cor can be arbitrarily scaled. Following the work by Dickens et al. (1995), the pseudo eigenvector is calculated throughφ = φ cor d cor , which in our case reduces toφ = d cor . The final step is to append the pseudo eigenvector to the original eigenvector matrix Φ as follows and subsequently solve for the balance equation in modal space defined in Eq. (36).

Model re-calibration
Linear ROM generally fail for problems that involve large deformations because the structural properties (stiffness matrix and internal forces) used for constructing the ROM are valid only for small δd. around the reference configuration d ref .
We solve this problem by updating the FEM discretization once the solution deviates significantly from the expansion point d re f used for linearization. This implies that also new augmented eigenmodes need to be computed.
Constructing and updating the ROM is expensive (due to the eigenvalue problem which needs to be solved) while applying it is very cheap. Efficiency for the proposed FSI method is achieved by re-using the reduced-order model as long as possible. We define a non-dimensional parameter based on the maximum absolute deflection δd max with respect to the reference frame d re f , i.e. the most recent linearization state, normalized by a characteristic length L of the structure. The ROM is adapted whenever exceeds a prescribed threshold. The efficiency and accuracy of the resulting Adaptive ROM (AROM) depends on this threshold.
Note that the limit case = ∞ corresponds to using the same ROM throughout the simulation, which minimizes the computational cost but will give inaccurate results if non-linear effects are significant, while = 0 corresponds to updating the ROM at each time step, which essentially yields the same accuracy as non-linear FEM at slightly increased computational cost.

Validation of the FSI-AROM algorithm
We analyze the accuracy and efficiency of the algorithm for three application examples. The first problem considers a purely linear structure and hence the update threshold is set to = ∞, which implies that the ROM model is built only once at the beginning of the simulation. The second and third example include large deformations and we search for a suitable problem independent threshold . For all cases, we report typical computation run times on a workstation (Intel XEON E5-2650).

Supersonic panel flutter
The first example is the aeroelastic instability of a thin plate exposed to a supersonic inviscid flow. This FSI test problem is often considered in literature (Teixeira and Awruch, 2005;Sanches and Coda, 2014;Pasquariello et al., 2016). Dowell (1974) has derived the critical flutter speed using linear stability theory and found that limit cycle the panel with 196 quadratic hexahedral elements in the streamwise direction and two elements along its thickness.
Since we are dealing with a two-dimensional problem, we use one element across the span. The plate has a Young's modulus of E S = 77.28 GPa, a Poission's ratio of ν S = 0.33 and a density of ρ S ;0 = 2710 kg/m 3 . The pressure of the free-stream is set to p ∞ = 28 kPa and the fluid density is ρ F;∞ = 0.339 kg/m 3 . For the fluid domain a grid-converged resolution with a total number of 16, 500 cells is used (Pasquariello et al., 2016). The grid is uniform with a cell size of ∆x = 4.25 × 10 −3 m and ∆y = 4.8 × 10 −4 m in proximity to the panel, see Fig. 2b. A cavity with a height of h = 2.2 × 10 −2 m is defined below the panel to account for its motion within the IBM framework. Slip-wall boundary conditions apply except for the inflow and outflow patch. As the flow is supersonic, we prescribe all flow variables at the inflow and use linear extrapolation at the outflow boundary. We use ALDM for the flux discretization and a CFL number of 0.6 for the Runge-Kutta time-integration method. The upper panel surface is coupled to the fluid while a constant pressure of p ∞ is applied at the bottom side within the cavity. The cavity pressure is reduced by 0.1% the first 4 ms to provide an initial perturbation.
Main results for the flutter analysis are presented in Fig. 3  eigenmodes for a Mach number range of 1.9 ≤ Ma ∞ ≤ 2.1. Flutter onset is predicted to occur at a critical Mach number of Ma ∞;crit = 2.08, with an error of 4.0 % with respect to linear stability theory (Dowell, 1974). Almost identical results can be found in the work of Pasquariello et al. (2016) who found a critical speed of Ma ∞;crit = 2.09, and in the work of Sanches and Coda (2014) who predicted flutter onset at Ma ∞;crit = 2.05. Figure 3b shows the influence of the number of eigenmodes, N eig , used in the modal database on the flutter prediction at Ma ∞ = 2.09. We observe monotonic convergence; 7 to 10 eigenmodes lead to an identical structural response as the reference FSI-FEM  The reduced-order model significantly improves computational performance of the FSI simulation. Solving the structural problem with the classical FEM approach costs about 77% of the total simulation time (wall-clock time).
The ROM reduces the time required for the structural problem to almost zero and leads to a more than fourfold speedup in this case, see table 1. We note that as a side-effect also the performance of the CFD solver improves when the ROM is used. This is due to the much smaller memory requirements of the ROM, which leads to a more efficient use of the CPU cache for the memory bandwidth limited CFD solver.

Elastic panel in a shock tube
Next, we study the impact of a shock wave on an elastic panel. This case is based on an experiment of Giordano et al. (2005), and was later numerically investigated by Sanches and Coda (2014) and Pasquariello et al. (2016). The setup is shown in Fig. 4 Poisson's ratio of ν S = 0.33. It is discretized using 55 × 2 quadratic hexahedral elements. The air flow is considered inviscid and compressible. We use ALDM for the flux discretization and a CFL number of 0.6 for time integration.
The fluid domain is discretized with 123, 400 cells with grid refinement around the panel, see Fig. 4. The inflow condition is based on Riemann invariants (Poinsot and Lele, 1992), and the remaining boundary patches mimic a slip-wall condition. The motion of the panel is mostly affected by its 1st bending mode, but in the following analyses we enrich the reduced model with the first 10 eigenmodes to ensure convergence.
We start our analysis with results obtained by the non-linear FEM approach. In Fig. 5 we show contours of the density gradient magnitude |∇ρ| at different times. Note that at t = 0 µs the shock wave has already hit the panel, which is the same definition as used by Giordano et al. (2005). At t = 140 µs the shock has passed through the small gap between the tip and the upper wall. A reflected shock due to the collision with the panel is also seen. Subsequently, the vortex generated at the panel's tip grows and moves downstream, followed by a shedding of small-scale vortices after t = 560 µs. The initial shock wave is reflected at the shock tube's end and then interacts with the main vortex, which results in a complex flow field at t > 840 µs.
The panel-tip displacement history is plotted in Fig. 6a and Fig. 6b for the 0.04 m and 0.05 m panel length case, respectively. We compare our results to experimental data of Giordano et al. (2005) and with numerical data from Sanches and Coda (2014) and Pasquariello et al. (2016) Larger deviations between the non-linear and linear models can be observed for the long panel. We will therefore only consider the case with l = 0.05 m for the AROM simulations. Figure 9 shows   Table 2 shows the devision of the wall-clock computation time between the CFD solver and the CSM solver for both the non-linear FEM (FSI-NLFEM) and the adaptive ROM (FSI-AROM) with an adaptation threshold of = 1.70 × 10 −3 . The FEM solver is responsible for about 31 % of the total wall-clock time, whereas the proposed AROM with the highest update frequency (lowest adaptation threshold of = 1.70 × 10 −3 ) consumes only a negligible amount (0.14%) of the computation time. In Fig. 10, the computational cost of the AROM with different thresholds is compared with the cost of a non-linear FEM simulation. As expected, the performance gain can be even larger if the threshold is relaxed.

Buckling of a shock-loaded thin semi-spherical membrane
The final application example is a three-dimensional FSI simulation of a thin shock-loaded membrane undergoing buckling (Pasquariello et al., 2016). It can be seen as an extension of the previous FSI cases to three-dimensional problems with complex structural behavior. Dynamic buckling is a non-linear structural phenomenon and highly sensitive with respect to any kind of imperfections, including grid resolution and modeling parameters (Ramm and Wall, 2004). This implies that tiny spatial variations in the loading of the structure may excite different buckling modes, which becomes even more evident for FSI problems, where the loads themselves are sensitive to the shape of the deformed body. Pasquariello et al. (2016) found that the occurring buckling mode can be affected by the structural resolution, while the sensitivity with respect to the fluid grid plays a minor role for the present test case.
The geometry and other setup details are shown in Fig. 11. The thin semi-spherical structure is hit by a right-  pressure signal has two distinct jumps, which indicate when the shock wave passes the sensor the first time (t = 0.12 ms) and the second time (t = 1.22 ms) after reflection at the end wall. The pressure signal is in excellent agreement with the data provided by Pasquariello et al. (2016). The sensor location is above the membrane and thus the pressure signal is not very sensitive to the motion of the structural interface. Contrary to the previous FSI examples, which were two-dimensional cases where a few number of eigenmodes sufficed, the current three-dimensional case is expected to require many more eigenmodes for capturing the local buckling of the structure. This is better understood by considering Fig. 14, where we show selected eigenmodes of the semi-sphere. Low-frequency modes represent a global motion of the structure and higher-frequency modes involve local deformations that are equally important for the present case.  Fig. 15a. We observe deviations from our non-linear FEM solution especially after the membrane collapses, i.e., after t ≥ 1.2 ms, which is not unexpected as multi-mode buckling is highly sensitive to numerical details (Ramm and Wall, 2004).
Next we study the effect of the threshold value for updating the AROM. In Fig. 16, the time-evolution of the RMS displacements are shown for various tolerances in the left column, and relative errors with respect to the nonlinear FEM reference solution are shown in the right column. For the error plots we blank the initial part, where very small reference displacements values would lead to ambiguously high relative errors. The cases with = 17 × 10 −3 and = 12.9 × 10 −3 have maximum errors above 25 % and 10 %, respectively, with the largest errors with N eig = 50 eigenmodes. For the highest update frequency for AROM, i.e., the smallest threshold of = 1.70 × 10 −3 , N eig = 50 modes lead to a maximum error of 5.5 % occurring at approximately t = 1.5 ms. Extending the modal base to N eig = 75 and N eig = 100 eigenmodes, while keeping the same threshold, further reduces the maximum error down to 1.6 % (at t = 1.5 ms) and 0.7 % (at t = 1.2 ms), respectively. In general, a threshold of = 8.60 × 10 −3 results in acceptable errors of less than 5.0 % when considering N eig = 75 or N eig = 100 eigenmodes. eigenmodes. The top row shows the deformation at t = 1.2 ms, just before the shock wave hits the structure for the second time. We clearly identify a compression of the windward side initiated by the initial shock passage. At t = 1.5 ms (bottom row), the shock has passed the sphere a second time and high-order (local) buckling becomes significant. We observe excellent agreement between the non-linear FEM and AROM results.    1.70 × 10 −3 ) considered in this example saves approximately 50 % (with N eig = 100 eigenmodes), 70 % (with N eig = 75 eigenmodes) and 80 % (with N eig = 50 eigenmodes) with respect to non-linear FEM. The symbols in Fig. 18b exemplarily depict the update instances of AROM for N eig = 75 eigenmodes when using different thresholds. The threshold and eigenmode-number dependence of the computational cost stems mostly from the eigenvalue solver. We solve the eigenvalue problem using a shift-invert method, which is very efficient for finding the lowest eigenvalues, while it results in strongly increased computational cost when searching for relatively high eigenvalues (Lehoucq et al., 1997). The computational costs are provided in table 3 for the case with N eig = 75 at selected thresholds. Again, it is evident that the AROM is very efficient compared to the conventional non-linear FEM method. AROM can reduce the total simulation time by a factor of 2 without noticeably compromising the accuracy of the FSI simulation ( = 4.30 × 10 −3 ).

Conclusions
We The accuracy of the AROM is controlled by a non-dimensional displacement-based parameter that triggers a recalibration step. Constructing and updating the modal basis is expensive due to the eigenvalue problem that needs to be solved. Efficiency is achieved by re-using the modal basis as long as possible. With very small threshold values, the AROM is adapted very frequently and the computational results as well as the computational cost converge to a space and time resolved non-linear finite-element simulation. A too large threshold, on the other hand, leads to an essentially linear model and possibly inaccurate results. We performed sensitivity studies for several test cases and found that a non-dimensional threshold value of about 4 × 10 −3 leads to the best balance between computational efficiency and accuracy for all cases. For future applications, we recommend to use an update threshold criterion based on the maximum strain, which can be determined a-priori for a given material.
The proposed method can be used with any partitioned Fluid-Structure Interaction (FSI) solver framework; the algorithm is independent of the baseline discretizations of fluid and structure. Our loosely coupled FSI implementation employs an unstructured finite-element discretization of the structural domain and a finite-volume method for solving the three-dimensional compressible Navier-Stokes equations on block-Cartesian grids with a cut-element immersedboundary method for representing the moving interface between fluid and solid. Using this FSI solver framework, the accuracy and efficiency of the AROM have been demonstrated and quantified for two-and three-dimensional FSI problems: We have shown that the model is accurate and very efficient for predicting the onset of flutter in supersonic flows. The AROM approach can by construction yield predictions with any required accuracy for shock-loaded structures undergoing large deformations, for which classical linear ROM would fail. The AROM also correctly reproduced the multi-modal buckling of a thin semi-spherical membrane with the same accuracy as the non-linear finite-element method and at significantly reduced computational cost. FSI simulations with the AROM maintain the excellent parallel scalability of the CFD solver by reducing the run-time requirements of the structural problem to a minimum, without noticeably compromising accuracy.