Model Studies of Fluid-Structure Interaction Problems

X.Sheldon Wang,Ye Yang and Tao Wu

1 McCoy School of Engineering,Midwestern State University,Wichita Falls,TX 76308,USA.

2 Former Graduate Students,New Jersey Institute of Technology,Newark,NJ 07102,USA.

Abstract: 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.

Keywords: Model reduction,fluid-structure interaction,mixed finite element,singular value decomposition,acoustic.

1 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 employ traditional mixed finite element formulations for very complicated acoustoelastic FSI problems in which there exist three different ranges of coupled modes,namely,sloshing,structural,and acoustic modes.It is important to notice that although we name these three frequency bands as sloshing,structural,and acoustic ranges,the coupled effects are fully incorporated in the formulation.For example,the structural frequencies do depend on added mass effects of the surrounding fluid.In engineering practice,it is essential to focus on particular problems which involve different frequency ranges and to develop corresponding efficient spectrum methods.In this paper,in addition to a comprehensive review of mixed formulations for acoustoelastic FSI systems,we also explore two different types of model reduction procedures,namely,mode superposition methods and singular value decomposition (SVD)-based model reduction methods.For mode superposition methods,the advantage is the direct calculation of mode shapes and natural frequencies,whereas the disadvantage is the requirement of an explicit formulation for the FSI system at that particular configuration,which is essential for nonlinear FSI problems such as vocal track acoustic resonance.For SVD-based model reduction methods,snap shots of transient solutions derived from linear or nonlinear FSI formulations or experiments can be used to formulate the dominant principal modes surrounding the physical configuration as illustrated in Schmid[Schmid(2010)].

A number of finite element formulations have been proposed to model an acoustic fluid for the analysis of acoustoelastic FSI problems,namely,the displacement formulation(Belytschko [Belytschko (1980)];Bathe et al.[Bathe and Hahn (1979)];Hamdi et al.[Hamdi,Ousset and Verchery (1978)];Belytschko et al.[Belytschko and Kennedy(1976)];Olson et al.[Olson and Bathe (1983)]),the displacement potential and pressure formulation,and the velocity potential formulation (Morand et al.[Morand and Ohayon(1995)];Everstine[Everstine(1981)];Olson et al.[Olson and Bathe(1985)];Felippa et al.[Felippa and Ohayon(1990)]).

Since its inception,primitive variable formulation has received considerable attention for simple and direct interface conditions no different from element assemblage processes[Kiefling and Feng (1976);Olson and Bathe (1983)].The key reason is of course the monolithic formulations for both media around the interfaces,which are useful in frequency calculations and response spectrum analysis.Various approaches have been introduced to eliminate spuriousnon-zerofrequency circulation modes and to obtain improved formulations [Wilson and Khalvati (1983);Chen and Taylor (1990);Bermúdez and Rodríguez (1994)].It has been proposed that the spurious non-zero frequencies are caused by the irrotationality constraint [Wilson and Khalvati (1983);Hamdi,Ousset and Verchery (1978)].Only recently,the true origins of the spurious non-zero frequency rotational modes have been identified [Wang (2008);Wang and Bathe (1997a)].It was concluded that the reported non-zero frequency spurious modes are caused either by the pure displacement formulation,the use of mixed elements which do not satisfy the Inf-Sup condition,or the improper treatment of the boundary conditions.Therefore,a proper way of eliminating the non-zero frequency spurious modes is to use the displacement/pressure(u/p)based mixed formulation with elements satisfying the Inf-Sup condition.In addition to the traditional u/pformulation (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/pand 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

2.1 Formulations for fluids

In linear analysis,for isentropic and inviscid fluids undergoing small vibrations,the momentum balance and constitutive equations show

whereρ,u,p,fB,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

whereVfandSfstand for the fluid domain and Neumann boundary,respectively;usnis the displacement normal toSf;andis the given pressure onSfcorresponding tousnwith 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λpis 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/pmixed finite element formulation or simply u/pformulation.Furthermore,because the irrotationality constraint is not imposed in the u/pformulation,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/pformulation

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β ≤α ≤106β.In this work,αis defined as 1000β.Invoke the stationarity of Π with respect to u,p,Λ and identifyλpaspandλΛas Λ,we have the three field equations

with the boundary conditions

whereSuindicates the Dirichlet boundary with

Notice that the pressure unknownpis 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,Hp,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.

2.2 Formulations for solids

Consider a structure with isotropic material

whereτ,fB,n,and fSfrepresent the stress tensor,body force vector,outward normal vector,and the applied surface traction vector on the Neumann boundarySf,respectively.Based on the governing equation and boundary conditions,the principle of virtual displacements can be derived:

whereis the virtue strain tensor;Rkcis thekthconcentrated load vector;δu,δuSf,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ν);Gis the shear modulus withG=E/2(1+ν);is the volumetric strain within Cartesian coordinates;andare the deviatoric strain components with respect to Kronecker delta,

As the Poisson’s ratiogradually approaches to 0.5,the bulk modulusκapproaches to infinity and the volumetric strainapproaches to zero.Naturally,it is difficult to evaluate the pressure from the relationshipeven though the pressure is finite.Therefore,the governing equation should involve another unknown,pressurep,which is independent of the displacement unknown u.Consequently,with the shear modulusG,the relation can be achieved as

Now the principle of virtual work is introduced in terms of independent variables u andp

we derive the following governing equation in matrix form with respect to discretized unknowns U and P

where BDand Bvare the deviatoric and volumetric strain interpolation matrix,respectively;R represents the external force corresponding to the net external virtual work ˉR;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 RBis the body force vector;ρis the mass density;andlists the nodal point accelerations,i.e.the second time derivative of the nodal point displacements U.At this point,fBno longer includes the inertia forces and the dynamic equilibrium equation becomes

2.3 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/pformulation for fluid

respectively,where Us,Uf,and Ucrefer to the nodal displacement unknowns in the structure interior,fluid interior,and fluid-structure interface,respectively;and Rs,Rf,Rsc,and Rfcstand 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 Ucbetween the fluid and solid domains,whereas the dynamic matching will cancel the effects of Rscand Rfcaccording 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,Ucmight 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 Mcc=Mscc+Mfccand Kcc=Kscc+Kfcc.

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?4cand 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.Lettercindicates 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?4celements for our test problems.

Figure 1:Two mixed elements for u-p-Λ formulation.Full numerical integration is used(i.e.,3×3 Gauss integration)

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?4celement,to describe the pressure of a generic element,a bilinear interpolation function is introduced

whererandsare 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

2.4 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/pformulation for structure in a typical acoustoelastic fluid-solid interaction system.The basic observations and conclusions are also directly applicable to the u/pformulation for fluid and displacement-based formulation for structure or solid.Assign

therefore,the coupled formulation with either u?p?Λ or u?pformulation for acoustic fluid can be expressed as follows

where Uhcontains all the variables of the nodal point displacements and Shlists 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?his an effective load vector.

Moreover,it is known that the shear modulusG=0 for invisid fluid,according to Eq.(24),we have

whereCis a constant depending on the time integration method,for example,withα ≥0.25 for Newmark method.

Since(Muu)his positive definite,(K?uu)his always a positive-definite matrix.The stability and effectiveness of the u/por 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(b)Inf-Sup Condition:

where the constantc2is independent of the mesh sizeh,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-Babuka condition,was popularized by Babuka and Brezzi in early 1970s[Babuka(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 (Kus)h,which is a direct extension or benefit of SVD-based model reduction method and does not introduce additional numerical challenges.

3 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.

3.1 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) =uiPi,where P is a square modal matrix,each column Piof which represents a mode shape,and u is a time-dependent general coordinate vector.Pre-multiply Eq.(45)by PT,we derive

The transformation matrix is established by solving the generalized eigenvalue problem

The homogeneous solution of Eq.(46) can be written asui(t) =cisinωi(t?ti),whereciis the scaling factor corresponding to thei-th eigenvector piand eigenvalueωi,tiis the corresponding time constant related to the phase shifts.Note that both the scaling factorciand the time shifttidepend 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.,pTiMpj=δij,we can finally define the transformation matrix as

can be rewritten as

with PTKP=Λ and PTMP=I,where I is ann×ndimensional 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)].

3.2 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 ann×mmatrix A with the column numbermmuch smaller than the dimension of the problemn.We can easily construct two symmetric matrices ATA(m×m)and AAT(n×n).Suppose the rank of the matrix A isr(),there will bernon-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 entireRmwhereas the left null space along with the column space span the entireRn.As a consequence,the eigenvectors of the matrix ATA denoted as viinRm,span the row space of the matrix A whereas the eigenvectors of the matrix AATrepresented with uiinRncomprise the column space of the matrix A.In fact,uiand viare similar to the right and left eigenvectors of the square matrix provided we havem=n.

Denote a matrix U,the firstrcolumns of which are the eigenvectors of the symmetric matrix AAT,and the lastm?rcolumns are constructed with the left null space vectors.In addition,introduce a matrix V,the firstrcolumns of which are the eigenvectors of the symmetric matrix ATA,and the lastn?rcolumns are made of the null space vectors.The singular value decomposition of the matrix A is then depicted as

where Ai=uivTiare the basis matrices and the firstrdiagonal entities of then×mmatrix Σ are the square root of the non-zero eigenvalues of the matrix ATA or AAT,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 largern-dimensional problem in the time domain,it is often more practical to use the model reduction method to construct generalized coordinates in a much smallerm-dimensional Krylov subspace [Dunteman (1989)].For example,let us havemtime snapshots within the time interval[t0,t1].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 xk,denotethe variance of this channel can be expressed asTherefore,we have a matrix x = (x1,x2,··· ,xm),the column of which represents each time snapshot.Moreover,we can also denote a new shifted matrix y = (y1,y2,··· ,ym),the column of which represents each time snapshot shifted by its own average,namely,yik=xik?.The procedure to maximize the variance inRnleads to an optimization problem:findφ ∈Rmsuch thatφTyTyφis maximized subject to the constraintφTφ=I,where I is anm×mdimensional identity matrix.It is then clear that the principal componentφis in fact the eigenvector of the covariance matrix C,theijthentity can be expressed asand the number of principal componentsrdepends 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 firstlprincipal components are selected,am×lmatrix Ψ=[φ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`eve decomposition(KLD).The mathematical proof of the equivalence of these methods is presented in Ref.[Liang,Lin,Lee et al.(2002)].

4 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 Rr,characterized as the white noise Vt,formal derivative of a Wiener processWt[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 = (xTyT)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 BVtwith a constant matrix B multiplied by a vector Vt,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ξ=φ?1X is introduced.Following the so-called Ornstein-Uhlenbeck process [Higham (2001)],we can then employ the Ito process for all the generalized coordinateξisatisfying the stochastic differential equation(SDE) with socalled additive noisecidWt:

where for theithmode the drift term is represented byλiξiand the diffusive term incrementcidWtyields fromφ?1BVtdtalong with the definition of the white noise as the formal derivative of a Wiener processWt[Ostoja-Starzewsk and Wang(1999)].

For theithmode,ifλiis zero,the general solutionξifor that mode will be Martingale.Consequently,we can express the explicit solution as

or an implicit approximation as

represents a normal distribution with zero mean and variance proportional toAssumeNis the number of realization,the sample variance can be expressed as

Moreover,based on the Ito calculus,the mean of the solutions for SDE(53)forithmode,can be expressed as

For any given initial value,repeatNdifferent simulations of sample path by implicit Euler method.Then useMdifferent initial values,i.e.,generating a series of trajectories with distinct starting points,the mean error of thejthbatch by the statistic is

forj=1,2,...,M,and their average

whereξi(t)k,jrepresents thek-th numerical simulation at timetof theithmode with thejthinitial value.

The variance of the batch averagecan be estimated by

For the student t-distribution withM?1 degree of freedom,a 90% confidence interval(α=0.1)for the mean errorμis

with

5 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.

5.1 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 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.

Figure 2:Model I:Acoustoelastic FSI system with a square acoustic solid vibrating in an acoustic fluid

Table 1:Frequencies of Sloshing/Structural/Acoustic Modes of FSI Model I

Figure 3:Displacement of sloshing/structural/acoustic modes of FSI Model I

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 havendisplacement unknowns,mpressure and vorticity moment unknowns,andkfree surface nodes,then the total number of zero frequencies isn?m?(k?1).Here,we assume that the physical constant pressure mode arisen from the boundary condition v·n = 0 onShas 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/pformulation 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,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.

Figure 4:(a)Numerical Inf-Sup test of 9?4c?4c elements for u?p?Λ formulation.(b) Lowest non-zero singular value of matrix [Kup Kul]. N represents the square root of the number of the elements used for any local element

As for the stability conditions,the ellipticity condition is automatically satisfied,since(K?uu)his always a positive-definite matrix.Therefore,the problem becomes to verify that 9?4c?4celements with corresponding mixed formulation satisfies the Inf-Sup condition.Fig.4(a)shows the Inf-Sup valuevs.the number of elementsNin log?log coordinates.It is clear that the Inf-Sup value approaches some value without decaying to zero,which indicates that the 9?4c?4celements pass the test.If an appropriate basis for displacement and pressure unknowns is selected,the Kupmatrix 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[KupKul],as shown in Fig.4(b),our proposed new numerical Inf-Sup test yield the same if not more visible conclusion.

Figure 5:Model II:Another acoustoelastic FSI model with a free surface and an immersed structure

5.2 Immersed slender solid/structure

Fig.5 describes the paradigm of another two-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.

Table 2:Frequency range of FSI Model II

Figure 6:Pressure distributions of sloshing,structural,and acoustic modes of FSI Model II

Figure 7:Rayleigh-Ritz quotient vs. nstep for sloshing and structural modes of FSI Model II

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 distributionis 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.

In Fig.8,we illustrate a procedure to use the dominant principal directions constructed with the small time snapshot intervaldtto extrapolate at a much larger time intervaldT.Suppose we start from the initial time and take 5 snapshots at time stepdtand we havenoverall degrees of freedom,an×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.

Figure 8:An illustration of coarse graining in time

Figure 9:RMSD error vs.dT

Figure 10:RMSD error vs.nstep.

Figure 11:Sample variance vs. sample number and confidence intervals of displacement for increasing number of batches based on mean error

Figure 12:Energy of transient combined modes and recovered data by SVD and relaxation

Using the recovered transient information,we can then extrapolate the transient system based on a much larger time stepdT.At the new time instancet+dT,after relaxation within a prescribed time duration which is often smaller thandT.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 todT.From the figure,we find that the slopes at the point 40dtand 60dtare 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 of snapshots.Therefore,we setdT= 60dtwith the number of time snapshotsnstepequal to 2000.Similar studies can also be carried out for fine time step sizedt.Fig.10 shows some results based on different time steps.In this case,considering the memory and numerical flops,dtis set to be 0.0005.Then we go back to the former procedure and relax the extrapolated solution for ever time stepdT.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 volumeV

It is shown in Fig.11 that the sample variance as defined by Eq.(56) begins to converge when the number of samplesNapproaches 50.Figure 11 also shows the confidence intervals of displacement and velocity for increasing number of batchesM,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 batchesMincreases,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.

6 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/pand 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.

Acknowledgement:The supports from NSF (DMI-0503652 and BES-0503649) and ASEE Summer Faculty Fellowship in 2008 and 2009 through AFRL as well as the Department of Mathematical Sciences at New Jersey Institute of Technology and the McCoy College of Science,Mathematics,and Engineering at Midwestern State University are greatly acknowledged.