Model Studies of Fluid-Structure Interaction Problems

In this work, we employ fluid-structure interaction (FSI) systems with immersed flexible structures with or without free surfaces to explore both Singular Value Decomposition (SVD)-based model reduction methods and mode superposition methods. For acoustoelastic FSI systems, we adopt a three-field mixed finite element formulation with displacement, pressure, and vorticity moment unknowns to effectively enforce the irrotationality constraint. We also propose in this paper a new Inf-Sup test based on the lowest non-zero singular value of the coupling matrix for the selection of reliable sets of finite element discretizations for displacement and pressure as well as vorticity moment. Our numerical examples demonstrate that mixed finite element formulations can be effectively used to predict resonance frequencies of fully coupled FSI systems within different ranges of respective physical motions, namely, acoustic, structural, and slosh motions, without the contamination of spurious (non-physical) modes with nonzero frequencies. Our numerical results also confirm that SVD-based model reduction methods can be effectively used to reconstruct from a few snapshots of transient solutions the dominant principal components with moderate level of signal to noise ratio, which may eventually open doors for simulation of long-term behaviors of both linear and nonlinear FSI systems.


Introduction
Over the past two decades, there have been a surge of research in the modeling of fluidsolid interaction (FSI) systems, many of which are related to the combination of nonlinear and meshless finite element methods with the immersed boundary method concepts [Boffi, Gastaldi, Heltai et al. (2008); Liu, Zhang, Wang et al. (2004); Wang (2006); Wang, Zhang and Liu (2009); Wang and Liu (2004) ;Liu, Liu, Farrell et al. (2006); Zhang, Gerstenberger, Wang et al. (2004)]. In this paper, we explore a further improvement of monolithic FSI modeling techniques by adopting model reduction techniques. More specifically, we to the traditional u/p formulation (as used for solids with zero shear modulus) [Wang and Bathe (1997a)], an effective three-field mixed finite element formulation, namely, the displacement/pressure/vorticity moment (u − p − Λ) formulation was presented in Wang [Wang (2008)]. In this paper, we briefly review the u/p and u − p − Λ formulations for acoustic fluid media and the related fluid-solid interaction (FSI) models. Moreover, we adopt the same mixed finite element formulations for both fluid and solid domains with mixed finite elements with continuous pressure. Hence, the coupling of fluids and solids in our monolithic fluid-solid is identical to standard finite element assemblage processes. Our numerical examples demonstrate that mixed finite element formulations are feasible for the prediction of coupled frequencies and mode shapes even if primary slosh, structural, and acoustic modes are within separate frequency ranges. A similar yet more straightforward numerical Inf-Sup test based on the lowest non-zero singular value of the coupled matrix. The focus of this paper however is on issues related to model reduction methods with mixed finite element formulations. One of the most popular model reduction methods is the mode superposition method [Daniel (1980); Bathe (1996)], the other one is the so-called Singular Value Decomposition (SVD)-based method or Principal Component Analysis (PCA) [Berrar, Dubitzaky and Granzow (2004); Jolliffe (2002)]. Recent studies include the overview of modeling of large-scale dynamical systems [Antoulas and Sorsensen (2001); Antoulas (2005)] and interpolating method based on diffeomorphism or differentiable manifold instead of tangent plane [Amsallem, Cortial, Carlberg et al. (2009)]. In general, the mode superposition method is limited to linear systems, whereas various forms of PCA can be extended to nonlinear systems with frequent updates of the principal directions. In order to illustrate the procedures of these two approaches, in this paper, we employ a fairly comprehensive two-dimensional acoustoelastic FSI model with both submerged twodimensional elastic structure and free surface. If we do not know the material properties of the system, with SVD-based principal component analysis, combinations of eigen modes can still be constructed by taking the first few dominant principal components of snapshots of available transient data which could also be derived from experiments or other black box simulation tools. In general, excellent agreements are confirmed between the original transient solutions and the data reconstructed with a few dominant principal components. The figures of energy are also plotted in order to verify the realization of this objective, namely, recovering the transient data with a few principal components without losing dominant characteristics. Finally, some temporal coarse-grained techniques are also proposed for the study of long-term behaviors of FSI systems.
2 Mixed formulations for acoustic continua

Formulations for fluids
In linear analysis, for isentropic and inviscid fluids undergoing small vibrations, the momentum balance and constitutive equations show where ρ, u, p, f B , and β represent the density, displacement vector, pressure, body force vector, and fluid bulk modulus, respectively. Assume the inertia force −ρü is included in the body force term and combine Eqs.
(1) and (2), we have In addition, it is known that the inviscid fluid must satisfy the following irrotationality constraint: The variational form of Eq.
(3) can then be expressed as where V f and S f stand for the fluid domain and Neumann boundary, respectively; u s n is the displacement normal to S f ; andp is the given pressure on S f corresponding to u s n with the following variational indicator In practice, it has been observed that with such a displacement-based formulation circulation modes with non-zero frequencies may emerge with the zero frequency circulation modes [Bathe (1996)]. In engineering practice, it is difficult to delineate such non-zero frequencies, often named as spurious or non-physical, from non-zero frequency physical modes. Furthermore, when the fluid acts as almost incompressible, the usual constraint of near incompressibility is shown as Eq.
(2), with the bulk modulus β → ∞. However, the displacement-based formulation with irrotationality and (almost) incompressibility conditions often lead to a much too stiff response, i.e., 'locking' behavior. In order to eliminate these deficiencies, the penalty-form variational indicator based on constraint (2) is introduced where λ p is the Lagrange multiplier. For this indicator, the first term corresponds to the strain energy with respect to pressure, the second term stands for the external body force potential, which also includes the inertia force, and the last term is the potential due to the boundary pressure. Eq. (6) is often called u/p mixed finite element formulation or simply u/p formulation. Furthermore, because the irrotationality constraint is not imposed in the u/p formulation, therefore exist too many exact zero frequencies. In reliable finite element analysis, these zero frequency modes are often eliminated with the following modification of Eq. (4) where Λ is the so-called vorticity moment and α stands for a very large number. With this additional constraint, another variational indicator Π is defined by the Lagrange multiplier method on the bases of the u/p formulation where λ Λ is a Lagrange multiplier vector corresponding to the vorticity moment. The fourth term is introduced to statically condense out the degrees of freedom of the vorticity moment Λ. From the numerical tests, as discussed in Refs. [Wang and Bathe (1997a,b)], it is found that α can be any numerically reasonable value larger than β, for example, 100β ≤ α ≤ 10 6 β. In this work, α is defined as 1000β. Invoke the stationarity of Π with respect to u, p, Λ and identify λ p as p and λ Λ as Λ, we have the three field equations with the boundary conditions where S u indicates the Dirichlet boundary with S = S u S f and S u S f = .
Notice that the pressure unknown p is a physical quantity in the formulation as governed by Eq. (10) as well as the boundary condition (13). However, the vorticity moment Λ is an artificial vector introduced to enforce the irrotationality constraint (4). Moreover, the governing equation (11) implies a small value for the vorticity. It is therefore natural to simplify the problem by imposing a zero value for the vorticity moment on all boundaries as described in (14). Detailed descriptions of this point and the relative issues are available in Refs. [Bao, Wang and Bathe (2001); Wang and Bathe (1997b)]. Applying the standard Galerkin finite element discretization, and considering a typical element, we have where U, P, and Λ list all the nodal point displacements, pressure, and vorticity moment, respectively in global coordinates and H, H p , and H Λ refer to the displacement, pressure, and vorticity moment interpolation matrix, respectively. Finally, the field equations can be expressed in matrix form with respect to nodal unknown vectors U, P, and Λ with Combining the coefficients of different variables, we can systematically rewrite the entries of the mass matrix and stiff matrix as follows with subscripts indicating their corresponding unknowns.

Formulations for solids
Consider a structure with isotropic material where τ , f B , n, and f Sf represent the stress tensor, body force vector, outward normal vector, and the applied surface traction vector on the Neumann boundary S f , respectively. Based on the governing equation and boundary conditions, the principle of virtual displacements can be derived: where δ is the virtue strain tensor; R k c is the k th concentrated load vector; δu, δu Sf , and δU are the virtue displacement over the volume, the Neumann boundary, and isolated points, respectively. If the material is almost incompressible, the constitutive relations are derived as where κ is the bulk modulus with κ = E/3(1 − 2ν); G is the shear modulus with G = E/2(1 + ν); v is the volumetric strain with v = kk in Cartesian coordinates; and ij are the deviatoric strain components with respect to Kronecker delta, ij = ij − v δ ij /3. As the Poisson's ratio ν gradually approaches to 0.5, the bulk modulus κ approaches to infinity and the volumetric strain v approaches to zero. Naturally, it is difficult to evaluate the pressure from the relationship p = −κ v , even though the pressure is finite. Therefore, the governing equation should involve another unknown, pressure p, which is independent of the displacement unknown u. Consequently, with the shear modulus G, the relation can be achieved as Now the principle of virtual work is introduced in terms of independent variables u and p with S = C , S = τ + pδ, and = − v δ/3, where S and are the deviatoric stress and strain tensors, respectively; the material matrix C is defined by Eq. (24) and linearly proportional to the shear modulus G; and δR corresponds to the net external virtual work. Implement the standard Galerkin discretization we derive the following governing equation in matrix form with respect to discretized unknowns U and P where B D and B v are the deviatoric and volumetric strain interpolation matrix, respectively; R represents the external force corresponding to the net external virtual workR; and the components of stiffness matrix are defined as In transient analysis, the inertia forces need to be considered. Using d'Alembert's principle, the element inertia forces can be directly included as part of the body forces where R B is the body force vector; ρ is the mass density; andÜ lists the nodal point accelerations, i.e. the second time derivative of the nodal point displacements U.
At this point, f B no longer includes the inertia forces and the dynamic equilibrium equation becomes where the matrix M = V ρH T HdV is the so-called consistent mass matrix of the solid.

Coupled formulation and mixed elements
In this section, we will explicitly illustrate the coupling procedures for mixed formulations.
Notice that there is no need to use the so-called mixed formulations with additional pressure and vorticity moment unknowns for the solid domain. Therefore, the additional pressure and vorticity moment unknowns will only show up for the fluid domain. Let's first employ the pure displacement-based formulation for solid and the u/p formulation for fluid respectively, where U s , U f , and U c refer to the nodal displacement unknowns in the structure interior, fluid interior, and fluid-structure interface, respectively; and R s , R f , R s c , and R f c stand for the nodal forces in the structure interior, fluid interior, structure, and fluid side of the fluid-structure interface,respectively. It is important to notice that in this case the kinematic matching of the coupling of fluid and solid domains share only the displacement unknown U c between the fluid and solid domains, whereas the dynamic matching will cancel the effects of R s c and R f c according to the Newton's third law. In Ref. [Wang (2008)], wang provides a more comprehensive explanation for the coupled fluid-structure matrices. In the practical analysis, U c might only consist of normal displacement components of fluid elements. In the assemblage process, Eqs. (33) and (34) are combined and yield the final governing equation for FSI systems with M cc = M s cc + M f cc and K cc = K s cc + K f cc . In order to choose appropriate mixed elements, we must consider the Inf-Sup condition [Chapelle and Bathe (1993)]. It has been proven that 9 − 4c − 4c and 9 − 3 − 3 mixed elements are appropriate choices in two-dimensional analysis [Bao, Wang and Bathe (2001); Bathe (1996)]. Figure 1 depicts the element configurations for these two types of mixed elements schematically. Take the second one as an example, there are 9 continuous nodal displacements and 4 continuous pressure and vorticity moment nodes in each element. When elements are adjacent to each other. The same assemblage procedures are utilized for the displacement, pressure, and vorticity moment unknowns. Letter c indicates that these quantities are continuous among elements. In order to obtain the continuity of the pressure and vorticity moment unknowns, we can also employ 9 − 4c − 4c elements for our test problems. In the mesh generation figure, it shows that each pressure node, except the nodes on the boundary, is shared by four neighboring elements. Based on the configuration of 9−4c−4c element, to describe the pressure of a generic element, a bilinear interpolation function is introduced where r and s are the natural coordinates for a single element. Comparing Eq.
(2) with (4), it can be observed that these two constraints are very similar in nature and are not coupled in the displacement-based formulation. Therefore, it can be assumed that they should be imposed with the same interpolation functions in generic element discretization. That is to say, for the inviscid acoustic fluid formulation, the interpolation for the vorticity moment is

Mesh generation and stability analysis
In the following discussion, the stability analysis is based on the u − p − Λ formulation for acoustic fluid and the u/p formulation for structure in a typical acoustoelastic fluid-solid interaction system. The basic observations and conclusions are also directly applicable to the u/p formulation for fluid and displacement-based formulation for structure or solid. Assign therefore, the coupled formulation with either u − p − Λ or u − p formulation for acoustic fluid can be expressed as follows where U h contains all the variables of the nodal point displacements and S h lists all the pressure and vorticity moment unknowns. For step-by-step transient analysis, after transforming, the following equation is considered for every time step: where R * h is an effective load vector. Moreover, it is known that the shear modulus G = 0 for invisid fluid, according to Eq. (24), we have where C is a constant depending on the time integration method, for example, C = 1/α t 2 , with α ≥ 0.25 for Newmark method.
Since (M uu ) h is positive definite, (K * uu ) h is always a positive-definite matrix. The stability and effectiveness of the u/p or u − p − Λ formulation and corresponding mixed elements which will be employed in the test problems are verified through ellipticity and Inf-Sup conditions, which are necessary and sufficient conditions for well-posedness: (a) Ellipticity Condition: with where the constant c 2 is independent of the mesh size h, the constant α in the irrotationality constraint and the bulk modulus β.
Generally, whether the Inf-Sup condition is satisfied depends on the specific chose of finite element discretization, the boundary condition, and the macromesh generation of the mathematical models [Bao, Wang and Bathe (2001)]. The Inf-Sup condition or often called Brezzi-Babuska condition, was popularized by Babuska andBrezzi in early 1970s [Babuska (1973); Brezzi (1974)] and stated in an earlier theoretical work by Necas in 1962.
References [Brezzi and Fortin (1991); Ern and Guermond (2003)] introduce this condition from the viewpoint of Banach space theorems. Although the Inf-Sup condition for mixed formulations was proposed some time ago, the analytical proof of whether the Inf-Sup condition is satisfied by a specific set of elements can be very difficult [Chapelle and Bathe (1993)]. Reference [Bathe (1996)] gives its very first numerical proof of Inf-Sup condition derived from the matrix equations and feasible for engineering applications. As a part of the novelties presented in this paper, we will explore a similar yet different type of numerical proof of Inf-Sup condition based on the singular values of the coupling matrix (K us ) h , which is a direct extension or benefit of SVD-based model reduction method and does not introduce additional numerical challenges.

Model reduction methods
Since the finite element procedure is realized on the basis of elements, in engineering practice, we often encounter a large number of elements and the degrees of freedom which lead to high-dimensional matrices. To accommodate the challenges in both spatial and temporal resolutions, we must develop efficient methods to construct equivalent lowdimensional dynamic descriptions which can capture the essence of the full fledged finite element model with high fidelity. In this paper, two different model reduction methods are introduced, namely, mode superposition methods and SVD-based model reduction methods.

Mode superposition method
The mode superposition method depends on the explicit construction of mass and stiffness matrices. If the coupled system is well known along with all material properties, this method can be efficient and useful for all linear systems. Suppose we have the following governing dynamic equilibrium equation without damping effects We propose to transform this equation into a more effective form for direct integration by using the following transformation x(t) = Pu(t) = u i P i , where P is a square modal matrix, each column P i of which represents a mode shape, and u is a time-dependent general coordinate vector. Pre-multiply Eq. (45) by P T , we derivẽ Mü +Ku =r(t), withM = P T MP,K = P T KP, andr = P T r.
The transformation matrix is established by solving the generalized eigenvalue problem The homogeneous solution of Eq. (46) can be written as where c i is the scaling factor corresponding to the i-th eigenvector p i and eigenvalue ω i , t i is the corresponding time constant related to the phase shifts. Note that both the scaling factor c i and the time shift t i depend on the initial conditions. To further simplify the formulation, we can also scale the eigenvectors according to the so-called M-orthonormal condition, i.e., p T i Mp j = δ ij , we can finally define the transformation matrix as P = [p 1 , p 2 , . . . , p n ], Λ = diag{ω 2 1 , ω 2 2 , . . . , ω 2 n } and Eq. (47) can be rewritten as with P T KP = Λ and P T MP = I, where I is an n × n dimensional identity matrix. Note that in engineering practice, due to the restraints of physical conditions, it is safe to assume the eigenvalue problem in Eq. (47) is simple, namely, the matrices are diagonalizable and no Jordan form is needed [Strang (1988)].

Singular value decomposition
In comparison with the mode superposition method, SVD-based model reduction methods do not depend on the explicit construction of mass and stiffness matrices. If a few snapshots of the transient solution are available, it is possible to identify the dominant principal directions and the proportionality of the energy content. In fact, this type of model reduction method is also applicable to nonlinear systems in which the principal directions can be altered or updated through time [Gunzburger (2003)]. It is also important to recognize that for linear systems, the principal directions are often a linear combination of the mode shapes which depends very much on the transient conditions or snapshots around that particular time.
Consider an n × m matrix A with the column number m much smaller than the dimension of the problem n. We can easily construct two symmetric matrices A T A (m × m) and AA T (n × n). Suppose the rank of the matrix A is r (r ≤ m n), there will be r non-zero eigenvalues of these two symmetric matrices. Based on the linear algebra theories [Jolliffe (2002)], the null space along with the row space span the entire R m whereas the left null space along with the column space span the entire R n . As a consequence, the eigenvectors of the matrix A T A denoted as v i in R m , span the row space of the matrix A whereas the eigenvectors of the matrix AA T represented with u i in R n comprise the column space of the matrix A. In fact, u i and v i are similar to the right and left eigenvectors of the square matrix provided we have m = n. Denote a matrix U, the first r columns of which are the eigenvectors of the symmetric matrix AA T , and the last m − r columns are constructed with the left null space vectors. In addition, introduce a matrix V, the first r columns of which are the eigenvectors of the symmetric matrix A T A, and the last n − r columns are made of the null space vectors. The singular value decomposition of the matrix A is then depicted as where A i = u i v T i are the basis matrices and the first r diagonal entities of the n × m matrix Σ are the square root of the non-zero eigenvalues of the matrix A T A or AA T , namely, singular values of the matrix A. Another popular and efficient modal reduction procedure based on singular value decomposition is the so-called principal component analysis. Instead of dealing with a much larger n-dimensional problem in the time domain, it is often more practical to use the model reduction method to construct generalized coordinates in a much smaller m-dimensional Krylov subspace [Dunteman (1989)]. For example, let us have m time snapshots within the time interval [t 0 , t 1 ]. Notice here that such a time interval should be shifted to explore the landscapes of interest to the research subjects in nonlinear problems; whereas such a shift is not necessary for linear problems. For a typical channel or snapshot x k , denotex k = n j=1 x jk /n, the variance of this channel can be expressed as n j=1 (x jk −x k ) 2 /n. Therefore, we have a matrix x = (x 1 , x 2 , · · · , x m ), the column of which represents each time snapshot. Moreover, we can also denote a new shifted matrix y = (y 1 , y 2 , · · · , y m ), the column of which represents each time snapshot shifted by its own average, namely, y ik = x ik −x k . The procedure to maximize the variance in R n leads to an optimization problem: find φ ∈ R m such that φ T y T yφ is maximized subject to the constraint φ T φ = I, where I is an m × m dimensional identity matrix. It is then clear that the principal component φ is in fact the eigenvector of the covariance matrix C, the ij th entity can be expressed as C ij = n k=1 (x ki −x i )(x kj −x j )/(n − 1), and the number of principal components r depends on the rank of the covariant matrix C. The eigenvalues, which are also called principal values, represent the variances of each principal component. Then the singular values are sorted in descending order with the eigenvectors (principal components) arranged accordingly. Furthermore, in engineering practice, the directions in which the data set has the most significant amounts of energy can be found. Once the first l principal components are selected, a m × l matrix Ψ = [φ 1 , φ 2 , · · · , φ l ] can be built up. This matrix is called feature vectors. Using the feature vectors, the projected data set can be obtained by xΨ. PCA involves a procedure that transforms a number of correlated variables into a much smaller set of variables and respective principal components or independent base vectors. The objective of PCA is to discover and to reduce the dimensionality of the data in the most effective way with respect to the variance of the data. Besides, it provides the most efficient way of capturing the dominant components of many processes with only few instead. Another method closely related to PCA and SVD is called proper orthogonal decomposition (POD) or Karhunen-Loève decomposition (KLD). The mathematical proof of the equivalence of these methods is presented in Ref. [Liang, Lin, Lee et al. (2002)].

Ito process and Brownian motion
In physical world, random perturbations often present themselves in dynamical modeling of fluid-solid interaction systems. A familiar way to quantify such random perturbations has been elegantly characterized as Einstein's Random Walk [Haw (2005)]. In fact, the coupled governing equation can adopt the following general form in which the external forces include a perturbation vector R r , characterized as the white noise V t , formal derivative of a Wiener process W t [Ostoja-Starzewsk and Wang (1999)]. Let y =ẋ, we can express the governing equation into a typical dynamical system form In this standard dynamical system form, with X = (x T y T ) T , so-called Langevin equation [Lemons and Gythiel (1997)], the governing equation has been abbreviated as in which the drift term a(t, X) could be simplified as AX with A a constant matrix for linear dynamical systems, and the diffusion term b(t, X) is often characterized as BV t with a constant matrix B multiplied by a vector V t , every entity of which is represented by a white noise random process. Assume the constant matrix A is simple, namely, diagonalizable, with A = φΛφ −1 , the generalized coordinate ξ = φ −1 X is introduced. Following the so-called Ornstein-Uhlenbeck process [Higham (2001)], we can then employ the Ito process for all the generalized coordinate ξ i satisfying the stochastic differential equation(SDE) with socalled additive noise c i dW t : where for the i th mode the drift term is represented by λ i ξ i and the diffusive term increment c i dW t yields from φ −1 BV t dt along with the definition of the white noise as the formal derivative of a Wiener process W t [Ostoja-Starzewsk and Wang (1999)].
For the i th mode, if λ i is zero, the general solution ξ i for that mode will be Martingale. Consequently, we can express the explicit solution as or an implicit approximation as W n represents a normal distribution with zero mean and variance proportional to t n . Assume N is the number of realization, the sample variance can be expressed as wherex N is the sample mean. Moreover, based on the Ito calculus, the mean of the solutions for SDE (53) for i th mode, can be expressed as For any given initial value, repeat N different simulations of sample path by implicit Euler method. Then use M different initial values, i.e., generating a series of trajectories with distinct starting points, the mean error of the j th batch by the statistic iŝ for j = 1, 2, . . . , M , and their averagê where ξ i (t) k,j represents the k-th numerical simulation at time t of the i th mode with the j th initial value. The variance of the batch averageμ j can be estimated bŷ For the student t-distribution with M − 1 degree of freedom, a 90% confidence interval (α = 0.1) for the mean error µ is

Numerical examples
The mathematical models considered in this paper are not restrictive, although for simplicity, we choose to work in two-dimensional solutions, the application can be directly extended to similar three-dimensional problems. To demonstrate the capabilities of proposed mixed formulations for both fluid and structure, we present the numerical results of two typical acoustoelastic fluid-solid interaction problems.

Immersed deformable solid block
The first model is a typical FSI system as shown in Fig. 2. It depicts a square deformable solid immersed in an open square cavity filled with an acoustic fluid. The immersed solid is connected to the bottom with two springs to compensate the effects of buoyancy and introduce some kind of resistance to rotational motion. In this problem, there exist three distinctive sets of modes, namely, surface wave sloshing modes, structural modes, and acoustoelastic modes. Table 1 lists the first few frequencies of sloshing/structural/acoustic modes based on different formulations and different mesh density. Convergence issues for mixed formulations will be addressed in a separate paper. In general, mixed formulations do   Figure 3: Displacement of sloshing/structural/acoustic modes of FSI Model I not yield monotonic convergence. However, it is clear in Tab. 1 that there exists a clear separation of sloshing, structure, and acoustic modes and the precise number of zero modes from our implemented program matches exactly with the theoretical prediction based on the same mesh. Furthermore, the first four sloshing, structure, and acoustic modes predicted by both displacement and displacement/pressure mixed formulations with refined meshes match very well with each other. It is also important to note that displacement/pressure formulation is not as accurate as displacement formulation for compressible solid. Displacement/pressure formulation is only needed when the solid is also incompressible.
To ensure there exists no non-zero spurious (not physical) frequency, we compare the number of zero frequencies with the mathematical prediction [Wang and Bathe (1997a)]. Assume that we have n displacement unknowns, m pressure and vorticity moment unknowns, and k free surface nodes, then the total number of zero frequencies is n − m − (k − 1). Here, we assume that the physical constant pressure mode arisen from the boundary condition v · n = 0 on S has been eliminated. From the table, we can see the numerical results are identical to the predictions. It is also important to point out that the u − p − Λ formulation for acoustic fluid works equally well with both the displacement formulation or the u/p formulation for immersed solid. Figure 3 describes the first four structural modes in addition to the first two sloshing and acoustic modes. Note here for mode shapes, the absolute scale is not relevant since there always exists a scaling constant. The upward direction is defined as positive. The top two sub-figures depict the sloshing modes, which describe the fluid motion on the free surface. The surface displacements are involved and the free surface is no longer flat.
Since the gravitational phenomena depend very much on the gravity, the corresponding coupled frequencies would be within low frequency ranges. The next four sub-figures are structural modes, which involve significant motions of the immersed solid. Naturally,  K up K ul ]. N represents the square root of the number of the elements used for any local element such modes will have higher ranges of frequencies than sloshing modes. Based on the configurations of the immersed solid, it can be recognized that the first and fourth structural modes are flexural or bending modes; the second one is a stretching mode and the third is a shearing mode. The last two sub-figures are acoustic modes, which focus on the coupling of acoustic modes between both immersed solid and the surrounding fluid. Of course, the key characteristics for the acoustic mode is reflected in the constant pressure condition on the free surface and much higher ranges of coupled frequencies.
As for the stability conditions, the ellipticity condition is automatically satisfied, since (K * uu ) h is always a positive-definite matrix. Therefore, the problem becomes to verify that 9 − 4c − 4c elements with corresponding mixed formulation satisfies the Inf-Sup condition. Fig. 4(a) shows the Inf-Sup value vs. the number of elements N in log − log coordinates. It is clear that the Inf-Sup value approaches some value without decaying to zero, which indicates that the 9−4c−4c elements pass the test. If an appropriate basis for displacement and pressure unknowns is selected, the K up matrix can be composed by a diagonal matrix followed by its kernel. Therefore, the minimum non-zero singular value is attempted to generate consistent information with respect to numerical Inf-Sup tests. As shown in Fig. 4(b), the lowest (non-zero) singular value does not decay as the number of elements employed increases. Since we are using the SVD in our model reduction methods, it is straightforward to employ the singular value of the coupling matrix to replace the standard numerical Inf-Sup test. With the given mesh densities which are self-accommodating, Fig. 4(a) demonstrates a possibility of the plateau in the standard numerical Inf-Sup test as described in Bathe [Bathe (1996)]; whereas utilizing the lowest singular value of the coupling matrix [K up K ul ], as shown in Fig. 4(b), our proposed new numerical Inf-Sup test yield the same if not more visible conclusion. dimensional acoustoelastic FSI model. In the following example, an elastic slender solid, similar to a structure, is immersed in acoustic fluid with the left hand side attached with the cavity boundary, a build-in support. Fig. 6 depicts the pressure distributions of the first two sloshing, structural, and acoustic modes, respectively. According to the convention, the pressure is assumed positive in compression. As shown in Fig. 6, for acoustic modes, the free surface essentially is a surface with a constant pressure which is consistent with traditional assumptions for underwater acoustics. Likewise, for sloshing modes, the pressure distributions match directly with the surface elevations in tank sloshing problems. Since the coupled structural natural frequencies are much higher than the coupled sloshing natural frequencies, it is not difficult to conclude that in this case, surface sloshing waves do not significantly influence the structural dynamical behaviors. Similarly, the structural modes hardly trigger the sloshing motion. Finally, for a test mesh, we have the following theoretical prediction of the zero frequencies as well as the ranges of the sloshing, structural, and acoustic modes in Tab. 2. In general, if we arbitrarily choose a primary unknown vector u, the Rayleigh-Ritz quotient α calculated according to the following formula could be a combination of all natural frequencies ω i . However, if the vector u represents a particular eigen mode such as sloshing or structural modes, the Rayleigh-Ritz quotient α will be approaching to the respective natural frequencies as we increase the number of snapshots, which is confirmed by our numerical test shown in Fig. 7. Assume that a white noise in Gaussian distribution N (0; t) is added to all the fluid nodes. Due to the disturbance of the white noise, the first singular vector, on which the transient data projected, could lose the characteristics of the corresponding eigen modes of the FSI system. Therefore, noise level is detrimental to the SVD-based model reduction method.
Raleigh-Ritz quotient: In Fig. 8, we illustrate a procedure to use the dominant principal directions constructed with the small time snapshot interval dt to extrapolate at a much larger time interval dT . Suppose we start from the initial time and take 5 snapshots at time step dt and we have n overall degrees of freedom, a n × 5 matrix A can be assembled. We can recover the transient data A by a selected number of principal components. Essentially, at every time instance, the coefficients or general coordinates related to these principal components are calculated by solving a normal equation.
Using the recovered transient information, we can then extrapolate the transient system based on a much larger time step dT . At the new time instance t + dT , after relaxation within a prescribed time duration which is often smaller than dT . We can then repeat the same process with this type of combination of fine and coarse time step sizes. Fig. 9 shows the Root Mean Square Deviation (RMSD) error between the original and extrapolated results with respect to dT . From the figure, we find that the slopes at the point 40dt and 60dt are very small (less than 0.03). Moreover, Fig. 10 suggests that as soon as the number of snapshots are sufficiently large, there is no need to further increase the number   Similar studies can also be carried out for fine time step size dt. Fig. 10 shows some results based on different time steps. In this case, considering the memory and numerical flops, dt is set to be 0.0005. Then we go back to the former procedure and relax the extrapolated solution for ever time step dT . Further discussions are also available in the conference publication of earlier results of this work [Yang, Wu and Wang (2010)]. In order to verify the effect of PCA, we also employ the first order Sobolev norm to represent the stored energy in the FSI system denoted with the total volume V It is shown in Fig. 11 that the sample variance as defined by Eq. (56) begins to converge when the number of samples N approaches 50. Figure 11 also shows the confidence intervals of displacement and velocity for increasing number of batches M , respectively. It can be seen from the figure that the length of the interval decreases as the number of batches increases. In other words, 90% of the times that upper and lower thresholds are calculated, the true mean lies below the upper threshold and above the lower threshold. As the number of batches M increases, the range between upper and lower thresholds gets smaller and smaller. Fig. 12 shows the energy of original transient displacements and recovered data by the first few basis matrices. In this case, the white noise is handled with Ito integration. The solid curve depicts the energy at discretized time snapshots while the dashed curve describes the energy of recovered data. This test confirms that even with a limited number of time snapshots the recovered data is fairly close to the original data for the aforementioned acoustoelastic FSI systems. Traditional mode superposition methods can only be applied to linear systems. Furthermore, the entire system behaviors have to be well established. Therefore, it is an effective tool to handle different external loads with different frequencies, for example, the spectrum analysis in civil engineering disciplines. However, SVD based model can be used without the prior knowledge of systems, the selection of dominant response is based on the magnitude of singular values which is a direct presentation of the proportionalities of energies within different response modes. Finally, we must point out the Krylov space constructed with the snapshots of the transient solutions can only be valid for the time vicinity at which the snapshots are recorded. The relationship between the linear Krylov space and the actual nonlinear transient behavior is very much the same as the tangent plane at a particular point of a curved surface. Therefore, for nonlinear problems, such Krylov space must be frequently updated depending on the actual nonlinear behaviors.

Conclusion
The acoustoelastic FSI problems are very challenging, in particular when free surface modes, structural modes, and acoustic modes are fully coupled together. In this paper, utilizing the monolithic u/p and u − p − Λ mixed finite element formulations for fluid and solid domains, we have demonstrated that coupled frequencies separated in different frequency ranges, can be effectively and efficiently captured in the same implementation. The crucial advantage of u−p−Λ formulation is that it can solve the acoustic FSI problems without introducing the spurious non-zero frequencies and many zero frequencies. In addition, we have confirmed that similar to the traditional numerical Inf-Sup test the singular values of the coupled matrix can also be used to identify mixed finite elements which satisfy the Inf-Sup condition.
In comparison with the mode superposition methods, we have shown through our numerical experiments that SVD-based model reduction methods do not necessarily depend on the explicit construction of mass and stiffness matrices and can be used to identify the dominant principal components with the proportionality of the energy content. In SVD-based model reduction methods, a few snapshots of the transient solution can be derived from linear or nonlinear FSI formulations as well as experiments. It is also important to recognize that for linear systems, the principal directions are often a linear combination of the mode shapes which depends very much on the initial and transient conditions. Again, unlike the mode superposition method, SVD-based model reduction methods are applicable to nonlinear systems in which the principal directions can be altered through time.