Galerkin Generalized Least Squares finite element methods for time-harmonic structural acoustics

Computer methods in applied mechaaics and engineering Comput.

ELSEVIER

Methods Appl. Mech. Engrg.

154

( 1998) 299-3 18

Gale&in Generalized Least Squares finite element methods for time-harmonic structural acoustics Karl Grosh”‘“, Peter M. Pinskyb “Department of Mechanical Engineering and Applied Mechanics, University of Michigan, Ann Arbor, MI 48109-2125, bDepartnzent of Civil Engineering, Stanford University, Stanford, CA 94305-4020, USA Received

24 March

1996; revised

1 February

USA

1997

Abstract Galerkin Generalized Least Squares (GGLS) finite element methods are developed and implemented for modeling the response of fluid-loaded structures. A general formalism for the application of this methodology to problems in structural acoustics (the vibration of an elastic structure in contact with an acoustic medium) is presented. A consistent approach for the development and implementation of the design parameters inherent to the GGLS methods for coupled problems is given. This general methodology is then applied to an exterior problem of a fluid-loaded Reissner-Mindlin plate. Complex wave-number dispersion analysis and numerical experiments demonstrate clearly the reduction of storage requirements and solution time of the new GGLS methods over standard discretization techniques. 0 1998 Elsevier Science S.A.

1. Introduction This work is motivated by the need to reduce the computational requirements of mid-frequency structural acoustics calculations. Structural acoustics involve the prediction of the vibration response of an elastic structure loaded by an acoustic fluid. In this paper we consider the time harmonic vibration problem. As the frequency increases, discretization requirements for accurate solutions also increase, limiting the size and complexity of models that may be studied. In this paper, Gale&in Generalized Least Squares (GGLS) finite element methods are developed for the coupled problem which result in a significant decrease in the real costs of computation (storage size and computational time) over standard finite element methods. The benefits of these methods increase with frequency. The consistent framework for the application of GGLS methods to structural acoustics problems is presented. Inherent to GGLS methods are design parameters which are selected to enhance the accuracy and stability of finite element formulations. An approach to appropriate selection of these design parameters for coupled problems is given which results in a formulation combining ease of implementation and general applicability to the important classes of structural acoustics problems (namely, interior and exterior problems). Least-squares methodologies has been successfully applied to time harmonic acoustics and vibration problems to enhance the accuracy of the approximation while reducing the computational burden. Harari and Hughes [ 1,2] developed Galerkin least squares (GLS) methods for the Helmholtz equation improving the accuracy of these approximations for the single and multi-dimensional settings. The formulation and accuracy of applying GLS to the acoustic pressure alone on structural acoustics has been studied using complex wave-number dispersion analysis [3]. In part, the dispersion analysis work motivated the development of GGLS methods for the

* Corresponding

author.

0045-7825/98/$19.00 0 1998 Elsevier PII SOO45-7825(97)00131-X

Science S.A. All rights reserved

300

K. Grosh, P.M. Pinsky

I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

dynamics of structures. It was found that the accuracy of To, selective reduced integrated Timoshenko beam elements deteriorated rapidly with increasing frequency. GGLS methods were applied with success to the vector Timoshenko beam equations in [4,5]. This work is extended to fluid-loaded structures in the current paper. Gradient least-squares methods were first developed by Franca and do Carma [6) and convergence of the method is proved for the scalar advection diffusion equation. They also suggest the use of multiple least squares operators in that paper. Harari and Hughes develop a multiple least squares method (combining GLS with a gradient least-squares method) for the steady scalar advection diffusion equation [7]. Finite element methods with residual based improvements for the static response of Timoshenko beams and Reissner-Mindlin plates are presented in [8,9]. Selection of the design parameters for GGLS formulations has followed the concept of matching the dispersion relations for model problems, put forth, for instance, by Franca and Dutra do Carmo. This approach was applied in [ 1,2] for GLS and Galerkin gradient least squares applied to the Helmholtz equation, and for vector problems in [4,5] for the Timoshenko beam. In this paper, an approach to the selection of the design parameters for coupled problems (with fluid variables and structural variables) and accuracy assessment techniques are presented. The organization of this paper follows. First, the strong and weak form of the boundary value problem for the coupled system are given. The general formulation for the GGLS finite element method for structural acoustics is next presented. In the following section, the general framework is made specific, as the methodology is applied to the problem of a fluid-loaded plate. The accuracy and computational requirements of the new method are studied using complex wave-number dispersion analysis and numerical experiments. Finally, conclusions are made and future work is discussed.

2. Boundary

value problem

for structural

acoustics

2.1. Strong form of the exterior structural acoustics problem In this section of the paper, the general elasto-dynamic coupled fluid-structure interaction problem is considered. For a specific application of the method, see Section 4 where a fluid-loaded plate is considered under a much less abstract and general setting. Consider an elastic body immersed in an infinite acoustic medium. The elastic body occupies the domain Q C [Wd, where d = 2 or 3 is the number of spatial dimensions (Fig. 1). If the dynamics of the elastic body are to be approximated by a reduced shell or plate theory, as will be studied in this paper, then the elastic domain is parameterized by a surface if d = 3 or a curve if d = 2. The interior surface of the structure is denoted by r-. The domain interior to Q is assumed to be evacuated; therefore only mechanical forces act on r-. The surface exterior to OS is denoted by r’. This surface (or curve) is the fluid-structure interface (sometimes called the wet surface) where the acoustic pressure loads the structure. The unit normal outward from the structure (inward toward the fluid) on r is ff E Rd. The infinite

Fig. 1. Geometry

for an elastic structure

immersed

in an acoustic

medium.

K. Grosh, P.M. Pinsky

I Comput. Methods Appl. Mech. Engrg.

154 (1998) 299-318

301

domain exterior to the surface r is denoted by $28C Rd. In Fig. 1, a finite fluid domain, defined by the region fif C 9 is pictured. 4 occupies the region between the surfaces r and %2&, where aaR is the truncation boundary of the fluid. partition The boundary f is assumed to admit the nonoverlapping

r=r; ur:,. The boundary

(1)

of the structural

domain,

admits the partition

u r" h, ’ ans =I-” 8,

(2)

where i = 1, . . . , ndof with ndof being the number degrees of freedom of the structural model. All variables are assumed to undergo time harmonic vibrations of the form e-‘“’ where i = fi and w is the angular frequency. With the geometrical definitions in hand we move to the formal statement of the strong form and the definitions involved therein. The strong form of the exterior structural acoustics problem is: find the displacement vector u : iIs + Cdof and the acoustic pressure p : !i%+ @ satisfying the governing differential equations: %P =f

in fir and 22 ,

(3)

2$=q-AspI,

in OS,

(4)

essential

and natural boundary

conditions:

P=g

on ri,

(5)

P,, = h

on ri,

(6)

Ui = gi

on ri,

i=l,...,ndof,

(7)

.LBju= hi

on ri.

i = 1, . . . , ndof ,

(8)

plus a consistency continuity

condition

relating

of normal displacement

n^~vp=&w2n^“~u

and radiation

conditions

f;

and r:norma,, see below:

at the fluid-solid

on r; expressed

(9) either as

lim r(d-2)‘2( p,, - ik,p) = 0

(10)

r-_)m

or through the Dirichlet-to-Neumann p,,

= -Mp

on &%R.

interface:

map as (11)

The pressure, p, satisfies the Helmholtz equation, (3), where -rip : = -V2 - kt, Vz is the Laplacian, the acoustic wave number is defined as k, = w/c, where c0 is the speed of sound in the fluid medium, p,, is the density of the fluid, and f: 0, + @, the acoustic sources are restricted to &?r. The structural matrix differential operator, ZS, is left as generic; symbolically if represents the governing differential equations of the structure and all attached substructures. ZS represents the differential equations arising from a reduced shell or plate theory. The structural response may be modeled using the equations of linear elasticity. Modeling the structure via elasticity theory is nicely presented in the thesis of Abboud [lo]. The applied mechanical forces are given by q : OS+ C” dof.As the displacement vector may contain rotational and translational degrees of freedom, we define n^’ E R” dof, the vector which when dotted with the displacement vector yields the normal displacement. The effect of the acoustic pressure acting in the normal direction on f is denoted by n^‘PI,-. r is in flS and represents the surface of the structure that is in contact with the fluid (r= flS Uaf) such that the acoustic pressure, p, is non-zero in OS only on IY The given functions which define the boundary conditions are g: ri + @, h: rl+ @, gi: ri, + @ and

302

K. Grosh, P.M. Pinsky

I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

h, : ri, + @. The normal derivative is defined as p,, : = Vp* II where n is the unit outward pointing normal from the domain 4, V is the gradient operator and the inferior comma denotes differentiation. On the fluid-structure interface, there is a consistency condition relating natural boundary conditions on the pressure (on r;) to essential boundary conditions on the normal structural displacement (define this domain as Tin,,,,,). For the remainder of this work, we enforce natural boundary conditions on the pressure through the essential boundary conditions on the normal component of the displacement. The components of u are ui. 9Ji is vector differential operator which operates on the displacement to yield the natural boundary condition. However, the exact form of the natural boundary conditions in (8) cannot be explicitly written down until we fix the governing differential equations for the structure (this is just a consequence of the abstract and formal nature of the presentation given in this section). The radiation condition (lo), is stated in terms of a radial variable r, where r is the distance from the origin of a spherical or polar coordinate system (see e.g. [ll]). This condition requires that the scattered pressure is outgoing at infinity. Application of this form of the radiation condition, however, is not suitable for domain based calculations. Therefore, an artificial boundary, a91R, is introduced. The artificial boundary is obtained by holding constant the radial variable of a Helmholtz separable [l l] coordinate system (e.g. the radius for a spherical system). On this special class of boundaries, one may write the exact impedance boundary condition, (1 l), imposition of which results in a problem equivalent to the original one posed on the infinite domain 9?_ The integral operator which relates the Dirichlet pressure data to the Neumann normal derivative data on the separable surface is the Dirichlet-to-Neumann (or DtN) map, M. The use of the DtN map to truncate finite element calculations for exterior acoustics was probably first described by Feng [12], but the methodology was crystallized in the work of Givoli and Keller [ 13,141. The application of the DtN for exterior acoustics problems is discussed in detail [ 13-151. The DtN formulation is valid only for acoustic sourcesf E of. Boundary integral techniques based on the free space Green’s function may also be used to achieve the truncation of the acoustic domain. REMARK 1. The above formulation assumes that p is purely outgoing. The inclusion of incoming acoustic excitation of the structure ,and fluid is straightforward. First one sets p = p”““’+ pi”’ where p”‘“’ is the unknown scattered pressure and p’“” is the known incident wave field. By using the fact that pi”’ satisfies the terms due homogeneous Helmholtz equation, a formulation for p”““’and u which now includes inhomogeneous to the incident pressure is obtained. 2.2. Weak form of an exterior structural

acoustics problem

The weak form of the exterior structural acoustics problem is given next. Denote the space of square integrable functions on 0 as L’(0) and the Sobolev space H”(0) as the space of functions with their first m derivatives in L*(n). First, we define the variational spaces necessary for the structural variables. Let the space of solution vectors for the displacement be denoted by 9” where YU = {U 1u : f& + Cndof, ui E H”~(L&), ui = gi on 4, where i=l,.. . , ndof}. The theory used to describe the structural response determines the value of m, for the ith component of 0. The space of weighting functions, YU, to which U belongs, differs from 9, only in that the components of U. The space of weighting functions, ‘I$, to which U belongs, differs from YU only in that the components of U, defined by Ui, satisfy homogeneous essential boundary conditions on r,,. Let the solution space for the acoustic pressure be denoted by Y;, where YP = {p 1: f& + @, p E H’(L&), p = g on G}. The space of weighting functions for the pressure, ‘VP, is the same as YP except that homogeneous boundary conditions are indicated for the functions in VP on 4. The weak form of the problem is: Given q, f, g, h, gi and hi for i = 1, . . . , ndof find u E 9, and p E 9, such that for all ii E VU and p E VP a’@, PI + (jTn^“, po02u)r

+ 6,

a’(& 24)+ (Z, n^“p), = L”(i) , where

MP),,~

= L’(F),

(12) (13)

K. Grosh, P.M. Pinsky

I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

303

(14) (15)

and

64)n =

I

nUHqdOr

(17)

the superscript H indicates the Hermitian of the vector. The structural operator, a‘(., .) represents an integral of the potential and kinetic energy densities for the structure which cannot be explicitly defined until governing equations for the structure are given. Note that all of the operators are sesquilinear (i.e. conjugate linear in the first slot and linear in the second).

3. Gale&in

Generalized

Least Squares methods

for structural

acoustics

In this section, Gale&in Generalized Least Squares (GGLS) methods are presented for the fluid-loaded structure. A particular case of these methods is the classical Galerkin method, which is obtained by setting the least-squares modifications to zero. The notation fl is used to denote a union of finite element partitions covering the domain R. Let finite element spaces .Y’i C Y,, %$ C Y,, 9: C 9” and Yi C “y; be defined by the restriction of YP, VP, YU and YU to the space piecewise polynomials of selected order for each space for p and U. The discrete approximation for the pressure and displacement are denoted as ph and uh [16]. The GGLS finite element method consists of the addition of a consistent residual proportional or least-squares term to the standard Galerkin finite element method. The GGLS finite method can be stated as: Find uh E Y”f: and ph E Lf’i such that for all iih E Yt and p” E 7:: A’c&Fht

Ph) + @“is,

po@2”h)r

+ ($3

Mph),%,

= &&jh)

3

(18)

nts fGGLS(Uh, z2) + (iih, rispqr

+ z

(H;~szh,

H;n^sph)r

= L&&?)

i=l

)

(19)

where 4s 4X,,

(Fh, ph) := afGh,

ph) + c

(Hf‘Tf;ph,

H;zfph)n,,

(H;Tsiih,

H;Lfsuh)o,

(20)

i=l

A;,,,(Uh,uh):= as(ih,uh)+ 2

,

(21)

i=l

4s

6”) := LftFh)+ tF,tHf%F*,Hff)n,, GGLS

(22)

%_s

(Uh):= L."t$") + c GGLS i=l

(H”2$ih,

H’q)n,,

(23)

where rr& and ne, are the number of least-squares operators added to the standard Gale&in equations for the fluid and structural equations, respectively. These methods are defined by the operator Hf and the ndof X ndof matrix differential operator H: each of which contain the design parameters enabling the development of new methods with enhanced accuracy characteristics. If Hf and Hr are set to zero, the Gale&in finite element method is recovered.

304

K. Grosh, P.M. Pinsky I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

For the Helmholtz equation, two methods that have been studied previously (GLS) method which is defined by & = 1 and

are the Gale&in

Least Squares

H’ =fif,

(24)

where rr is the GLS design parameter. defined by & = 1 and

The second is Galerkin

Gradient Least Squares (GVLS) method which is

Hf =&V,

(25)

where rfv is the GVLS design parameter. Harari and Hughes have also studied a combined GLS-GVLS approach (e.g. & = 2) for the advection equation [7]. Optimal selection of the GLS and GVLS parameter for scalar problems is discussed in [16,17]. The structural least-squares operator depends on the particular theory used. Suppose for definiteness that ndof = 2 then the GLS method is defined by n[s = 1 and

The GVLS method would include a gradient operator on the diagonal of the matrix in (26). The development of the optimal r’s for the GVLS method applied to the transverse vibrations of Timoshenko beams is presented in r4s1. REMARK 2. The essence of the GGLS methods is to add to the standard Galerkin finite element equations, terms proportional to the residual of the original governing differential equations integrated over element interiors. If we fix n& = nts = 1 and drop the subscript, the GGLS equations (18) and (19) may then be written as

HfGfph -f))n,+af(ph,ph)+(phd,p/g2~h)y+(ph,MPh)~~~=Lf(ph),

(HfY$,

(HSLQih, H”(.L@

-q

+ tisph))n, + a”@, Uh) + (2, fi “p”), =

L”@) ,

where the residual proportional term is on the first line in each of these two equations and standard Gale&in equations are on the second line. Note that p is non-zero in Q only on r, the wet surface. These methods method (18) and design of methods Hy and Hf must frame).

4. Application

are consistent by construction (i.e. the exact solution, when substituted into the numerical (19), satisfies the equations identically). The addition of the least-squares terms allows the that enhance the accuracy of the approximate solution. Finally, the definition of the operators also satisfy conditions of rotational invariance (method independence of global coordinate

of GGLS for fluid-loaded

plates

Next, we make definite the formalism of Section 3 through application of the method to the exterior problem of a fluid-loaded Reissner-Mindlin plate. The new method combines the optimal GGLS parameters for the in vacua dynamic beam response [4,5], with a GLS method in the fluid [2] and for the extensional response. 4.l.Coupled

Fluid-plate

difSerentia1 equations

Let 4 = 0 E Iw* represent the region which the acoustic fluid occupies and the curve r denote the fluid-structure interface and the parameterization of the mid-surface of the plate. The domain interior to r, i.e. W*\Q is taken to be evacuated. We consider a plate which is infinite in the third space dimension with uniform loading in that coordinate. Thus, the acoustic problem becomes two-dimensional and the plate problem is equivalent to a Timoshenko beam formulation. Finally, for purposes of this study, fis = f, i.e. the entire structure is ‘wet’.

K. Grosh, P.M. Pinsky I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

p: f2 + @ and the displacement

The coupled equations for the fluid-pressure, effects of fluid loading are

u : f + C3 including

305

the

,iip,u =q-n^“p

on r,

(27)

Lztp =o

in 0 ,

(28)

where 5’ is the outward normal from the plate accounting for the rotational degree of freedom, Zf is the Helmholtz operator. The definition of the elements of the matrix differential operator ss for Timoshenko beam and classical rod theory is (considering only interconnected straight beams): -KGA

Z$

a2 2 - pAw2

-KGA~

=

KGA z

a

0

2

a

KGA-E+-PIo~

0

(29)

2

0

0

;

+A-$-pAw2

I1

where s is the coordinate along the midsurface of the beam/rod, uT = {w 8 u} where w is transverse displacement (i.e. normal to the midsurface), 8 is rotation and u is longitudinal displacement (along the midsurface). The shear correction factor is K, the shear modulus is G, the Young’s modulus is E and the mass of the beam is p. The cross-sectional area of the beam is A and the moment of inertia is I. The force vector is 0 q,}. The displacement vector is also subject to either natural or essential boundary conditions. qT = (4, The natural boundary conditions operators can now be defined (see (8)): B,=~KGA{~ B2=T

{

-1

0

B3=2{0

EI; 0

0

O},

>

,

(30)

EA;},

where the + arises from the direction of the outward normal of the beam at the boundary (see e.g. [18]). The homogeneous Helmholtz equation is considered with the extension of the inhomogeneous case being straightforward. The fluid pressure is coupled to the structural displacement through the continuity of normal displacement boundary condition (9) on r The fluid pressure must also satisfy radiation boundary conditions which we express in terms of the DtN condition on the truncation boundary, a,%&, (11). As mentioned in Section 2.1, this is the exact boundary condition producing a problem that is equivalent to the original infinite domain problem. This boundary condition involves an infinite sum of angular harmonics. When implemented in practice, the sum is truncated and an approximate boundary condition results. The effect of the number of terms taken in the series on the accuracy and uniqueness of the formulation is discussed by Givoli and Keller [13], Keller and Givoli [ 141 and Harari and Hughes [ 151. 4.2. Variational The variational

equations for the coupled system equations

are

a”@, 24)+ (U, rip), = L”(iE) , a’@, PI +
+
(31) = 0,

(32)

where @ and UT = {W U 8) are the weighting functions. As mentioned in Section 2.2, any natural boundary conditions for p on r are enforced via essential boundary conditions on the normal displacement. The operators for the coupled beam/rod and fluid equations are

306

K. Grosh, P.M. Pinsky I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

(33) and

(34) with

(35) The operator af(*, *) is

af(j5,p)

:=

I

a [Vj *VP - I&p]

do.

(36)

Let the variational function spaces Ya and “yb, be in H’(T). The solution space, Ya, is the set of all such functions satisfying the essential boundary conditions for (Y.The weighting function space, YU, is made up of functions whose value is zero where essential boundary conditions are specified. The function spaces are then defined as u E Y,, := Y,+,X YO X Yu and U E “v; := VWX 7fo X Tu. The weak form of the problem is: given q, f -and appropriate boundary conditions for u and p, find (u, p) E Yu X .Y, such that for all (u,p) E “Yu X ‘I(, the variational equations, (3 1) and (32), are satisfied.

4.3. Galerkin

Generalized

Least Squares jinite element formulations

In this section, the GGLS method for the fluid-loaded plate is presented. In the notation of Section 3, we will study methods with nLs = n& = 1, i.e. one least-squares operator each for the fluid and structural equations. Further, the approach involves a generalized least-squares method as the GVLS method is applied to the transverse structural vibrations and the GLS method is applied to the Helmholtz equation and to the longitudinal beam motion.

4.3.1. Finite element formulation The finite element interpolations of the displacement and pressure are uh and ph, respectively. Linear interpolation is used for the plate variables and bilinear interpolation is used for the acoustic pressure. Hence, the finite element subspaces are the restriction of the original subspaces (defined in Section 4.2) to the space of piecewise linear functions. The GGLS equations may be written as A&,,

Hn sph)r = L”@) )

(Uh, Uh) + (Uh, li sph>r + (H~jih,

A;Ls@h, ph) + ($6

po~2uh)r

+ ($3 Mph),%, = 0,

(37) (38)

where ‘CIS = af@,

ph> + (Tf2’,

r$fPh)n

9

(39) (40)

where H is the GGLS matrix differential

operator including

the design parameters.

K. Grosh, P.M. Pinsky I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

307

(41)

Next, the motivation described.

for the selection of the above form of the GGLS operators in the fluid and structure will be

4.3.2. Selection of design parameters and operators for jkid-loaded structures Our approach for developing a GGLS method for the coupled setting was to combine methods which were successful in the uncoupled cases. This was motivated from our experience in [3] as we found that refinement of only one component of a coupled problem (i.e. either the structure or the fluid representation) did not result in an improvement of the coupled solution. However, if the representation of both media can be improved, an increase in the accuracy of a particular level of discretization is obtained. The GGLS approach presented in this paper, least-squares improvements are used for the transverse and longitudinal beam response and for the fluid pressure. The design parameters for each of these problems has been selected such that there is zero dispersion error in uncoupled model problems. That is, the allowed free waves in the finite element solution for the Timoshenko beam, longitudinal rod and one-dimensional Helmboltz equation match their respective analytic counterparts. The GVLS method described in [5] for the transverse plate vibrations is used. This gradient method consistently removes shear locking and provides a frequency dependent correction to ensure zero dispersion error in both the propagating and evanescent waves. The development and presentation of results for this method applied to in vacua beam vibrations are given in Grosh and Pinsky [4,51 along with a detailed description of the design parameters. The gradient approach was found to most directly achieve the desired alleviation of locking in the zero frequency limit along with enabling a frequency dependent correction for zero dispersion error. It is emphasized that these elements are fully integrated and are not based on any form of reduced integration. This gradient least-squares method is coupled to a least-squares approach applied to the longitudinal plate and the acoustic fluid elements. The use of GLS for the one-dimensional and two-dimensional Helmholtz equation is described in [ 1,2]. The value of the GLS parameter that yields no dispersion error for the one-dimensional acoustic (non-coupled) problem using linear interpolation functions can be written as 1 - cos(wh/c,) 2 + cos(wh/c,)

’

(42)

where c0 is the speed of sound for the longitudinal waves in the rod or in the fluid. For the rod element, h is the length of the element. For two-dimensional problems, a characteristic element length dimension is used for h (see e.g. [2,19]). Hence, rf and r3 have similar form, with different material properties and element lengths. Note that the 7, terms inherent to the GGLS method for the transverse beam response equations result in a new fluid-pressure coupling term, (HZsiih, Hfi’~~)~ This term adds an unsymmetric component to the formulation. This non-symmetry has implications when evaluating computational efficiency and storage requirements as will be addressed in the later results sections. The symmetric formulation is given in Appendix A. Of central importance is the approach determining the design parameters for the GGLS method. In choosing the design parameters for the coupled problem, the r’s determined optimal for the uncoupled model problems are used. These methods are consistently combined to form the complete method. Alternative approaches, such as developing design parameters optimal specifically for the coupled problem were not pursued. These methods would lack generality, being specialized to only one particular problem type (e.g. exterior problems). The specialization arises from the fact that the dispersion relations for an interior (wave guide) problem is different than the dispersion relation for an exterior problem. New methods would need to be developed for other problem configurations (e.g. the interior problem). Further, it is anticipated that design parameters developed specifically for the coupled problem would be quite cumbersome to calculate. Our approach amounts to improving the accuracy of the uncoupled problems separately, then synthesizing the coupled method from the consistent combination of the improved methods. This approach combines the simplicity of design parameters

308

K. Grosh, P.M. Pinsky I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

and the generality of application (as no one problem has been targeted). of this approach can be found in the results, Section 5.1.

Further discussion

of the effectiveness

5. Results for GGLS method Two methods for the solution of the structure and fluid response variables are studied. One method is denoted as the standard Galerkin approach which consists of SRI elements for the transverse vibration of the plate (therefore this method is not a strictly displacement based method as SRI is equivalent to a mixed method) combined with Gale&in elements for the longitudinal plate vibrations and acoustic fluid elements. The other is the GGLS approach described in Section 4.3.2. Complex wave-number dispersion analysis and numerical experiments are used to analyze and compare the accuracy of the two methods. 5.1. Accuracy

of dispersion

relations for the &id-touded

plate

To study the accuracy of the new formulation, the model problem of the infinite fluid-loaded plate is used. The derivation of the dispersion relations for the discrete representation is found in [3] and with some discussion in [20]. In this paper, we outline the manner in which these dispersion relations are derived. The finite element dispersion relations provide a general framework for the a priori error measure of the new design. Derivation of the analytic dispersion relations for fluid-loaded plates can be found in [11,21] and are summarized in [3]. To obtain the finite element dispersion relations for the fluid-loaded Reissner-Mindlin plate, a uniform mesh in the fluid and structure is used, see Fig. 2. Just as for uncoupled systems, the assembly of the finite element matrix equations results in repetitive stencils for the fluid and structural variables. Exponential solutions in the plate and fluid are assumed as follows: w,

=

w.

e

ik,n

Ax

n

= WOPU,1

0, = 0, eikxnAx = sop: , Pn,m

=

POe

ik,n

Ax

eikYmAy = pop:py

(43) ,

where k, and k, are the wave numbers of the coupled system (the k, component of the fluid wave vector that coincides with the plate direction while the k,, component is normal to the plate, into the fluid), n and m are the node indices (see Fig. 2); AX and Ay are the interelement spacing in the x and y directions; we study meshes with AX = Ay = h in this paper. Substitution of the assumed forms into the stencils yields a non-standard condition for the existence of allowed waves in the coupled system. By seeking non-zero solutions for the pressure amplitude, the non-standard condition is reduced to a standard determinantal condition from which the non-trivial free wave solutions are calculated (see [3]).

WC-1 gec-l PLATE MESH Fig. 2.

Mesh for the fluid-loaded plate.

K. Grosh, P.M. Pinsky

I Comput. Methods Appl. Mech. Engrg.

309

154 (1998) 299-318

The allowed waves in the finite element mesh are compared to the analytic dispersion relations for the fluid-loaded Reissner-Mindlin plate which are quintic in k,. For a description of the analytic dispersion relations see [21-231 for a detailed comparison of the analytic and finite element dispersion relations for the coupled system (see [3]). We recount some of the important points regarding the k, dependence of the wave numbers satisfying the dispersion relation. There are five branches of the dispersion relation. One is a purely real branch which is denoted as the subsonic locus (the wave speed is less than the acoustic wave speed). This wave number is most important as it corresponds to a flexural mode which dominates the surface displacement. The next most important branches are the two evanescent wave-number loci. The wave numbers on these loci are complex conjugates of each other over the entire frequency range and consist of a small real part and dominant decaying or imaginary part. The final two wave-number loci, whose relation to one another changes with frequency (although they exist as complex conjugates over much of the frequency range) are called the leaky wave-number loci. The k, dependence of the leaky wave numbers, whose contribution is considered least important over the frequencies studied, is that of a predominantly real part with a smaller imaginary part. The discussion of the importance of the various wave-number loci involves a discussion of branch cuts in the complex plane (see ]3,21,221). Results in the section are for a steel plate (E = 210 dynes/cm*, p = 7.8 gr/cm3) in water (pO = 1 gr/cm3, c0 = 148 100 cm/s). The element spacing to plate thickness ratio is 1 (results for thinner plates are similar). The qualitative similarity of the finite element dispersion relations to the analytic dispersion relations for the Galerkin and GGLS formulation is excellent (see [3]). Each of the five wave-number loci, with their respective complex components, are replicated in the finite element dispersion relations. We do not present these curves but rather move directly to quantitative error plots which enable more direct accuracy comparisons. The phase and amplitude error of the kx component wave-number loci are presented. The error plots for the SRI plate element coupled to the Gale&in fluid approximation are shown in Figs. 3 and 4. In Fig. 3 the error in the real part of the three wave-number branches with non-negligible real parts is plotted; this represents the phase error for a wave propagating along the plate. In Fig. 4, the error in the imaginary part of the wave-number branches with non-negligible imaginary components is plotted; this represents an amplitude error of the wave as it propagates down the plate. For the new GGLS method, the phase error is shown in Fig. 5 and the amplitude error in Fig. 6. The error of the important subsonic wave-number locus remains below two percent over the entire frequency range (compare with Fig. 3 where the error in the Galerkin method subsonic wave-number locus rises above 15%). The evanescent wave-number branch is perfectly replicated up to the cut-off frequency for this branch in the k, direction (the k,” wave numbers for this branch have a dominant real component which reaches a value of k,h = T around oh/c, = 2.7, hence the mesh is too coarse to represent this wave (see Grosh and Pinsky [3])). This represents a major improvement over the Gale&in method with the significant amplitude error of this branch as seen in Fig. 4. The evanescent branch is an important contributor around discontinuities and drive points; hence, accurate representation of this component is critical to accurate solutions for realistic (finite) 1.20

1.20

1.15

_

1.15

2 z-

1.10 1.05

z* $w

.

.

42 v

1.00 0.95

,” 6

xp:

0.90

d

1.05 1.00

x d

iii/,

0.85 0.80

1.10

0

0.5

1.0

1.5

2.5

2.0

, 0

3.0

0.5

1.0

1.5

,

,

2.0

2.5

,I 3.0

whlco

Fig. 3. Phase error for the real components the fluid. Steel in water, h/t = 1.0.

of k, wave numbers for a fluid-loaded

Fig. 4. Amplitude error for the real components Steel plate in water, h/t = 1.0.

of k, for a fluid-loaded

plate. Finite elements: SRI for the plate and Galerkin for

plate. Finite elements:

SRI for the plate and Gale&n

in the fluid.

310

K. Grosh, P.M. Pinsky

/ Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

Fig. 5. Phase error for the real components the fluid. Steel in water, h/t = 1.0.

of k, wave numbers for a fluid-loaded

pIate. Finite elements: GGLS,,,

for the plate and GLS in

Fig. 6. Phase error for the real components the fluid. Steel in water, h/t = 1.0.

of kx wave numbers for a fluid-loaded

plate. Finite elements: GGLS,,,

for the plate and GLS in

structures. If we consider vertical branch cuts from kk, in the complex plane for the square roots in the analytic dispersion relation (see e.g. [21,22]), the leaky wave-number root only enters the first Riemann sheet for w/r/c, > 1.0; free waves exist only when the roots of the dispersion relation are on the first Riemann sheet [21] and [3]. At these frequencies, the error of the GGLS leaky wave-number branches, in both phase and amplitude, remains below three percent. The standard Galerkin approach has an error which increases with frequency and is greater than five percent over this frequency range. Therefore, the GGLS approach has a lower error than the standard Gale&in approach for every branch of the dispersion relation at frequencies where free waves are possible. This superior accuracy was attained by a consistent combination of optimal methods for the uncoupled systems. The complex wave-number dispersion analysis methodology developed in [3] enables an error estimate for the new method to be obtained in a general setting. An alternative approach to developing new methods for the coupled problem would have been to match all or part of the fifth order, fluid-loaded plate dispersion relations. Although possible in principle, a method matching certain of the loci would result in design parameters without a closed form solution and in a cumbersome method to implement. Further, this design would be specialized to the case of exterior fluid-loading and perhaps not be accurate for application to the interior problem. As mentioned previously, the approach taken here is not biased toward either the interior or exterior problems as the design parameters are selected from model uncoupled problems. In the next section, the significant reduction in dispersion error achieved by the GGLS method will be seen to translate into a significant reduction in actual computational requirements for the coupled problem. 5.2. Numerical

experiment

In this section, a two dimensional fluid-loaded plate problem is solved numerically. The results of the new formulation are compared to that of the standard Galerkin approach. In Fig. 7, a schematic of four interconnected plates immersed in an infinite fluid is shown. The long side of the rectangular configuration is 6.0 cm while the short side is 1.O cm. The plates are steel with a thickness of 0.15 cm. The plate material parameters are for a steel plate in water, given in the previous section. The DtN boundary is at R = 6.5 cm. The forcing function for this problem is a unit point load (one dyne) in the vertical direction at x = 3.0 cm and y = 0.0 cm. The plate admits longitudinal motion as well as transverse displacement and rotation. The plates are rigidly connected with continuity of displacement and rotation conditions at the joints. Hence, the longitudinal force will excite both extensional and transverse waves in the plate system. Due to the symmetry of the loading, the solutions exhibit amplitude symmetry about the horizontal axis. Results for two frequencies, k, = 1.0 cm-’ and k, = 2.0 cm-‘, are presented. Different meshes were used for comparing convergence and accuracy. The number of elements used in each mesh is given in Table 1. The meshes are used at both frequencies, and a typical mesh is shown in Fig. 8. These meshes are generated in the following manner. Uniform spacing is applied along each plate section and equi-angular circumferential spacing

K. Grosh, P.M. Pin&y I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

Fig. 7. Schematic

311

of four plate problem.

Fig. 8. Intermediate

mesh.

Table 1 Meshes for plate problem Mesh name

Circumferential

Coarse Intermediate Fine Finest

56 140 224 420

elements

Total elements 1904 4760 9408 17640

is used on the DtN boundary. Each radial line connects the inner plate node to the DtN node. The radial spacing of the nodes along each line is uniform, except for the first four nodes in the fluid, whose spacing is half that of the other nodes. The amplitude of the fluid pressure and normal plate displacement are presented. Contour plots of the amplitude of the pressure in the fluid domain are studied at each frequency to examine the global character of the solution and the general accuracy of the solution schemes. The normal displacement on the fluid-plate interface is then presented to obtain a quantitative estimate of the accuracy of both methods. The normal displacement gives rise to the acoustic field and serves as the initial amplitude for the outgoing waves. The presence on the DtN boundary is also used as quantitative measure of the error of the approximate solution. Amplitude errors in the normal plate displacement are seen to propagate out to the DtN boundary. The DtN boundary is at R = 6.5 cm, hence k,R > 6.5 for both frequencies studied. The values of this boundary yield the far field directivity pattern of the solution. 5.2.1. Results for k, = 1.0 cm-’ The amplitude of the fluid pressure and normal plate displacement for k, = 1.0 cm-i are given next. The pressure amplitude contours for the GGLS solution on the intermediate mesh, shown in Fig. 9, is visually the same as for the converged solution. For this frequency we consider the converged solution achieved when the maximum change in nodal values is less than a one percent change in the nodal amplitudes between consecutive finite element meshes. This occurs as the refinement is made from the intermediate to the fine mesh GGLS solution. In addition to convergence of the numerical solution, the spatial variation of the pressure and displacement on the fluid-plate interface is examined. For the converged solution, the wavelength with which the amplitude varies is the same as the subsonic wavelength (see [3] and Section 5.1). The subsonic flexural wave dominates the response and the correspondence between the analytic and numeric results gives further confidence as to the accuracy of the approximation. The GGLS solution on the fine mesh is used to provide the converged solution. The GGLS coarse mesh solution is shown in Fig. 10. Here, excellent coarse mesh accuracy is seen as the pressure amplitudes and directivity pattern are well represented. The Gale&in solution on the coarse mesh is shown in Fig. 11 and on the fine mesh in Fig. 12. While the

K. Grosh, Phf. Pinsky I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

312

Y

Fig. 9. Pressure

amplitude

Fig. 10. Pressure

Fig.

contours

amplitude

for GGLS finite element method, intermediate

contours

mesh, k, = 1.0 cm-‘.

for GGLS finite element method, coarse mesh, k, = l.Ocm-‘.

11. Pressure amplitude contours for Galerkin finite element method, coarse mesh, k, = 1.O cm-‘.

Fig. 12. Pressure

amplitude

contours

for Galerkin

finite element method, coarse mesh, k, = l.Ocm-‘.

overall characteristics and directivity pattern are well approximated by the coarse mesh Galerkin solution, there is a noticeable amplitude error. Errors in the solution on the Gale&in fine mesh solution can be seen in the contour plots (compare Figs. 9 and 12). These qualitative remarks are made more definite by examining the normal displacement on the plate surface and the pressure on the DtN boundary. The normal displacement on the upper plate is presented in Fig. 13 for the GGLS method and in Fig. 14 for the Galerkin method for various meshes. The normal displacement on the surface of the plate predicted by the Gale&in method on the finest mesh (not shown) overlays the converged solution. The accelerated convergence for the new method is clearly seen as the intermediate mesh provides sufficient resolution for a converged GGLS solution while a fine mesh is required for the Gale&in approach. Similar convergence occurs for the pressure amplitude on the DtN boundary shown in Figs. 15 and 16. amplitude errors present in the normal plate displacement are seen to propagate out to the DtN boundary. This effect is present in both methods. In Fig. 15, the convergence of the pressure amplitude at the DtN boundary is shown for the GGLS method and in Fig. 16 the results for the Galerkin approximation are given. The GGLS approach achieves the converged result on the intermediate mesh while the fine mesh is required for the Gale&in method to attain a similar but less accurate solution. Next, storage and computational requirements of the application of the two methods on these two meshes are compared. Let nzqns and nzqns denote the number of equations required for satisfactory accuracy to be achieved by the symmetric Gale&in and the unsymmetric GGLS approaches respectively. Noting that each plate node has three degrees of freedom, we have for the fine mesh nEqns = 9856 and for the intermediate mesh nzqns = 5040. Let b, and b, denote the mean half bandwidth of the meshes for the symmetric Gale&n and unsymmetric GGLS solution schemes, respectively. Note also that for the unsymmetric method, the upper and lower bandwidths are the same. For this frequency, b, = 217 (line mesh) and 6, = 136 (intermediate mesh). The storage requirement for the symmetric method is bsn& and for the unsymmetric formulation 2bunrqns; the ratio of the storage requirement for the

313

K. Grosh, P.M. Pinsky I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318 GaIerkin - coarse mesh -A-

ggls - coarse mesh -ag&-intermediate

Gale&in

-B-

_ intermediate

Galerkin converged

soln. -

_ ibe

converged

5

+

mesh -os&L. -

5

lOE12

10E12

Y 1I

0

3I

24

40

60

54

0

1

2

distance

3

4

5

6

distance

Fig. 13. Upper plate normal displacement for GGLS finite element method, k, = 1.0 cm-‘. Fig. 14. Amplitude of the normal displacement on the upper plate for Galerkin finite element method, k,, = 1.O cm-‘.

0.06

0.06

0.01

0.01

Ol 3lM

200

100

0

I

0

0

200

100

R

300

e

Fig. 15. Pressure amplitude on DtN boundary for GGLS finite element method, k, = 1.O cm -I. Fig. 16. Pressure amplitude on DtN boundary for Gale&in finite element method, k, = l.Ocm~‘.

Galerkin

approach

operation

counts

to that for the GGLS approach

of floating point band solver is approximately nIqns(bt + 3b,) [24] and for the unsymmetric band solver 2n,“,,,bf. The ratio of the flop estimates for the direct solution of the equations resulting from the Galerkin approach to that for the GGLS approach is 2.5. Hence, even at this low frequency, significant storage and computational gains are achieved by the GGLS approach. Note, any added cost of calculating the design parameters for the GLS and GGLS methods has been neglected in the analysis; these would be proportional to the number of total elements if all elements were of different size. Reduction in the calculations needed for the design parameters is achieved when elements have the same characteristic size. (flops)

for direct

solution

is 1.56. Considering

of the matrix

equations,

asymptotic

the numbers

estimates

of flops for the symmetric

5.2.2. Results for k, = 2.0 cm-’ Next, the pressure and displacement amplitudes of the fluid-loaded plate at k, = 2.0 cm-’ are studied. At this frequency, the condition for convergence is relaxed to a less than five percent nodal amplitude change between meshes. The pressure amplitude contours for the fine mesh, intermediate mesh and coarse mesh GGLS approximations are shown in Fig. 17, 18 and 19, respectively. The intermediate mesh and fine mesh Gale&in solution are shown in Figs. 20 and 21. The intermediate Gale&in solution bears a qualitative resemblance to the

K. Grosh, P.M. Pinsky I Compur. Methods Appl. Mech. Engrg. 154 (1998) 299-318

Fig. 17. Pressure

amplitude

contours

for fine mesh using GGLS finite element method, tine mesh k, = 2.0 cm-‘.

Fig. 18. Pressure

amplitude

contours

results for GGLS finite element method, intermediate

Fig. 19. Pressure Fig. 20. Pressure

amplitude

amplitude

contours

contours

mesh k, = 2.0 cm-‘.

results for GGLS finite element method, coarse mesh k, = 2.0 cm-‘.

for results Galerkin

finite element method, intermediate

mesh k, = 2.0cm-‘.

3 2 lOEl1

1 0

1

2

4

3

5

distance

Fig. 21. Pressure Fig. 22. Amplitude

amplitude

contours

results for Galerkin

of the normal displacement

finite element method, finest mesh k, = 2.0 cm-‘.

on the upper plate for GGLS finite element method, k, = 2.0 cm-‘.

6

K. Grosh, P.M. Pin+

315

/ Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

converged solution. Note that the pressure amplitudes are four times greater than the converged solution. This large error is due to the poor phase accuracy of the SRI plate element giving rise to ‘numeric’ or false resonance behavior as seen in the in vacua plate finite element models [5]. The fine mesh solution has an amplitude error of about 50%, similar to the error in the coarse mesh GGLS solution. The results for the normal plate displacement obtained through the GGLS method are shown in Fig. 22. The rapid convergence of the GGLS approach is seen as the intermediate mesh provides sufficient resolution to obtain the converged solution. It is noted that the analytic subsonic wave number is again matched by the spatial variations on the surface of the plate-fluid interface for the converged solution. The results for the normal plate displacement obtained through the Galerkin approach are shown in Fig. 23. As the mesh is refined from the coarse mesh solution to the intermediate mesh solution, the poor phase accuracy of the SRI plate element clearly manifests itself, as an increase in the amplitude error occurs with the mesh refinement. From the intermediate mesh to the finest mesh, the finite element solution converges monotonically to the actual solution. Under the Galerkin approach, one must use the fine mesh (9408 elements) to exceed the accuracy of the coarse mesh GGLS solution (1904 elements). The finest mesh solution for the Galerkin scheme achieves the converged result with 17 640 elements. The pressure amplitude results on the DtN boundary for GGLS solution on the coarse and intermediate meshes are plotted in Fig. 24 and for the Galerkin solution on the coarse, intermediate, and fine meshes in Fig. 25. These results show the propagation of the amplitude errors of the normal displacement on the surface of the structure out to the far field. The application of the GGLS method results in an even greater advantage, as measured by the reduction in computational effort and storage, over the Gale&in approach as frequency is increased. Using the storage and asymptotic flop estimates of the previous section, we compare the requirements of the converged GGLS approach to the Gale&in approach. For the GGLS approach on the intermediate mesh we have b, = 136 and neUqns= 5040; for the Galerkin approach on the finest mesh we have b, = 404 and Al&“, = 18 480. The ratio of the storage requirements for the Galerkin approach to the GGLS approach is 5.5. The ratio of asymptotic flops for the direct solution of the equations resulting from the Gale&in approach to those required by the GGLS approach is 32.0. In addition to a reduction in the storage requirements, the amount of computational effort is also reduced using the GGLS approach. The benefits gained from the GGLS method increase with frequency, as the correction from the least-squares addition becomes more important. For the frequencies studied here, the GGLS

Gale&in

intermediate

Gale&in G&&in

&

- muse meah -A-

Gal&in

- finemesh - finest mesh

-c-

mesh -A-

ggls-intermediate

-m-

-oconverged

dconverged sol”. -

1

- coarse

4

3

2

5

e

of the normal displacement

Fig. 24. Pressure

amplitude

-

6

distance

Fig. 23. Amplitude

s&l.

on the upper plate for Galerkin

on DtN boundary

finite element method, k, = 2.0 cm-

for GGLS finite element method, k, = 2.0 cm-‘.

’

316

K. Grosh, P.M. Pinsky

I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

Gddin--mcrh* Galdin

- intermediate-D-

Gdbrkin-l%Umsb-OGakrkh-Bnest mesh-xcmlvagtd s&l. -

-

“0

100

290

300

3

Fig. 25. Pressure amplitude on DtN for Galerkin finite element method, k, = 2.0cm-‘.

method for the transverse vibrations of the beam was mainly responsible for the dramatic improvement in the accuracy and reduction in computational effort as the longitudinal wavelengths are several times longer than the length of plates. At higher frequencies we would expect the GLS correction for the longitudinal vibrations also to play a role. REMARK 3. The unsymmetric component of the GGLS method only enters on the fluid-structure interface, other than this the method is completely symmetric. Advantage may be made of this fact to easily reduce the storage requirements from 2b~r&~ to bunzsns + nfsibU, where nfsi is the number of nodes on the fluid-structure interface. REMARK 4. In the numerical study undertaken, a very regular mesh was used. The only refinement of the mesh to the anticipated solution was the use of smaller elements close to the structure, to capture the subsonic near field. Much of the circumferential refinement of the uniform mesh was needed on the fluid-structure interface and such refinement was not needed on the DtN boundary. As the subsonic waves propagate through the fluid, a filtering effect occurs, smoothing the solution greatly. A re-meshing strategy that truly adapts to the character of the solution would yield significant benefits, for it is the reduction of the number of elements on the DtN boundary where the greatest decrease in the profile could be achieved. A successful refinement scheme would locate a greater density of elements near the surface of the structure and fewer elements at the DtN boundary. Such a refinement scheme would yield great benefits in reduction of storage and computations for both the Galerkin and GGLS schemes [25].

6. Conclusion Gale&in generalized least-squares (GGLS) finite element methods for structural acoustics have been developed and implemented. The GGLS method and design parameter selection approach afforded a dramatic reduction in computational time and storage requirements for accurate solutions compared to the same performance metrics for the standard Gale&in approach. Further, the GGLS approach was more attractive as the frequency is increased. A general formulation was presented with the flexibility to include multiple least-squares operators for the structural response variables and the acoustic pressure. This flexibility was utilized in solving a model problem, consisting of a fluid-loaded Reissner-Mindlin plate. In the GGLS formulation implemented in this work, a truly generalized approach was undertaken as a Gale&in least squares (GLS) method for the fluid pressure was

K. Grosh, P.M. Pinsky I Comput.Methods Appl. Mech. Engrg. 154 (1998) 299-318

317

combined with a structural approach comprised of a Galerkin gradient least-squares method for the transverse plate motion and GLS for the longitudinal vibrations. The key element to the development of efficient GGLS methods is the appropriate selection of the design parameters. We have developed an approach which involves the consistent combination of optimally selected GGLS parameters from uncoupled model problems to form the complete coupled methodology. This approach results in a formulation which is straightforward to implement and possesses general applicability to structural acoustics problems with interior, exterior and simultaneous interior-exterior fluid loading. An important area of future work is the development of rigorous proofs for the convergence of the GGLS finite element method for the coupled problem. The success of this GGLS formulation for the two-dimensional structural acoustics problem holds the considerable promise for the application of GGLS methods to three-dimensional configurations. Central to the extension to higher dimension is the development of GGLS methods for plate elements with arbitrary geometries and axisymmetric and general shell elements. A combination of mixed finite element interpolations of the structural variables, to alleviate shear locking, with GGLS methods, to provide a frequency dependent correction to the formulation, should lead to a class of efficient structural elements. With the structural elements in hand, the design of methods for the coupled problem should follow directly from the approach advocated in this paper.

Acknowledgments We would like to thank Manish Malhotra for his help in coding and mesh generation aspects of this paper. We would also like to acknowledge the funders of this work, the Office of Naval Research.

Appendix

A. Symmetric

formulation

for residual modifications

The formulation presented in Section 3 is a non-symmetric version of the residual proportional modification. These elements were shown to provide a substantial improvement in computational cost as compared to the standard Gale&in approach. In this appendix, the symmetric version of the residual modification, a true ‘least-squares’, method is given. Implementation of this approach will be given in a future paper; however, the implication is clear, the symmetric version will result in a tremendous savings in both storage and computational cost over the standard Galerkin method (see the storage and flop considerations in Sections 5.2.1 and 5.2.2.). The symmetric of the GGLS finite method can be stated as: Find uh E 9’: and ph E 9: such that for all Z”EY,h andphEY:: ph) + Ghis,

&.,,s6h,

p,,mZuh)r + (ph, Mph),,

R

= L&J$)

nls Uh, u”) + (Uh, riSph), + c

f&LS(

nks (H;Lz$h,

H;iisph),

i=l 6,s

+

c

(H;_@h,

H;rv~h)r

,=I

“Fh, H,sn^

“p”>p = J!&&Zh)

+ c

@z;q, H;hs~h)r,

04.2)

i=l

where the operators are defined in Eqs. (20)-(23). If we specialize the GGLS operators as in Remark symmetric form of the GGLS equations are

(H”( 2$”

+ 2

4s (H$

i=l

(Hf.@?,

(A.1)

Hf(sfph

-f>)n,

2 (i.e. fix nts = &

= 1 and drop the subscripts)

+ afGh, ph) + @n^, pfW2Uh)r + Gh, l14ph),$?, = LfGh)

the

)

+ n^5” ), H”( 5$uh - q + n^sph))@2s(Uh, Uh) + (Uh, n^“p”), = .L”(tlh) ) L I * newterm

where the symmetrization

consists of adding the term ri “E” to the least-squares

operator. Note again that pi is in

318

K. Grosh, P.M. Pins.@ I Comput. Methods Appl. Mech. Engrg. 154 (1998) 299-318

4 and represents the surface of the structure that is in contact with the fluid (T = 4 U@) non-zero in OS only on r.

such that j and p are

References [l] I. Harari and T.J.R. Hughes, Finite element methods for the Hehnholtz equation in an exterior domain: model problems, Comput. Methods Appl. Mech. Engrg. 87 (1991) 59-96. [2] I. Harari and T.J.R. Hughes, Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains, Comput. Methods. Appl. Mech. Engrg. 98 (1992) 41 l-454. [3] K. Grosh and PM. Pinsky, Complex wave-number dispersion analysis of Gale&in and Galerkin least squares methods for fluid-loaded plates, Comput. Methods Appl. Mech. Engrg. 113 (1994) 67-98. [4] K. Grosh and P.M. Pinsky, New generalized Gale&in least squares finite element methods for wave propagation in Timoshenko beams, in: R. Kleinman, et al. eds., Proc. SIAM Second Intl. Conf. on Math. and Numer. Aspects of Wave Prop. Univ. of Delaware (June 1993) 237-245. [5] K. Grosh and PM. Pinsky, Galerkin Generalized least squares methods for Timoshenko beams, Comput. Methods Appl. Mech. Engrg. 132 (1996) 1-16. [6] L.P. Franca and E.G. Dutra do Carmo, The Gale&in gradient least-squares method, Comput. Methods Appl. Mech. Engrg. .74 (1989) 41-54. [7] I. Harari and T.J.R. Hughes, Stabilized finite element methods for steady advection-diffusion with production, Comput. Methods Appl. Mech. Engrg. 114 (1994) 165-191. [8] A.F.D. Loula, T.J.R. Hughes, L.P. Franca and I. Miranda, Mixed Petrov-Galerkin methods for the Timoshenko beam problem, Comput. Methods Appl. Mech. Engrg. 63 (1987) 133-154. [9] T.J.R. Hughes and L.P. Franca, A mixed finite element method formulation for Reissner-Mindlin plate theory: Uniform convergence of all higher-order spaces, Comput. Methods Appl. Mech. Engrg. 67 (1988) 223-240. [lo] N.N. Abboud, A mixed finite element formulation for the transient and harmonic exterior fluid-structure interaction problem, Ph.D. Dissertation, Department of Civil Engineering, Stanford University, 1990. [l l] M. Junger and D. Feit, Sound and Structures and their Interaction (MIT Press, Cambridge, MA, 1986). [12] K. Feng, Asymptotic radiation conditions for reduced wave equation, J. Comput. Math. 2(2) (1984) 130-138. [13] D. Givoli and J.B. Keller, A finite element method for large domains, Comput. Methods Appl. .Mech. Engrg. 76 (1989) 41-66. [14] J.B. Keller and D. Givoli, Exact non-reflecting boundary conditions, J. Comput. Phys. 82( 1) (1989) 172-192. 1151 I. Harari and T.J.R. Hughes, Analysis of continuous formulations underlying the computation of time harmonic acoustics in exterior domains, Comput. Methods Appl. Mech. Engrg. 97( 1) (1992) 103-124. [ 161 T.J.R. Hughes, 1987, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis (Prentice-Hall, Englewood Cliffs, NJ, 1987). [17] I. Harari, Computational methods for problems of acoustics with particular reference to exterior domains, Ph.D. Dissertation. Deparmtent of Mechanical Engineering, Division of Applied Mechanics, Stanford University, 1991. [18] E.B. Magrab, Vibrations of Elastic Structural Members (Sijithoff and Noordhoff, Germantown, MD, 1979). [19] L.L. Thompson and PM. Pinsky, A Galerkin least squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Methods Engrg. 38(l) (1994) 371-397. [20] R. Jasti, Mixed shell finite elements with applications in structural acoustics, Ph.D. Dissertation, Department of Civil Engineering, Stanford University, 1992. [21] D.G. Crighton, The free and forced waves on a fluid-loaded elastic plate, J. Sound Vib. 63(2) (1979) 225-235. [22] J.F.M. Scott, The free modes of propagation of an infinite fluid-loaded shell thin cylindrical shell, J. Sound Vib., 125(2) (1988) 241-290. [23] C. Seren and S.I. Hayek, Acoustic radiation

from an insonified elastic plate with a line discontinuity,

J. Acoust. SOC. Am. 86(l) (1989)

195-209. [24] G. Golub and C.F. Van Loan, Matrix Computations, 2nd edition (Johns Hopkins, Baltimore, MD, 1989). [25] I. Haran’, K. Grosh, T.J.R. Hughes, M. Malhotra, P.M. Pinsky, J.R. Stewart and L.L. Thompson, 1995, Recent developments element methods for structural acoustics, Arch. Comput. Methods Engrg., to appear.

in finite