Stable symmetric finite element - boundary integral coupling methods for fluid-structure interface problems

Engineering Analysis with BoundaryElements 15 (1995) 79-91 Copyright © 1995 Elsevier Science Limited 0955-7997(95)00021-6

ELSEVIER

Printed in Great Britain. All rights reserved 0955-7997/95/$09.50

Stable symmetric finite element - boundary integral coupling methods for fluid-structure interface problems* Xiaogang Zeng Algor Inc., 150 Beta Drive, Pittsburgh, Pennsylvania 15238, USA

&

Jaeobo Bielak:~ Computational Mechanics J"aboratory, Department of Civil and Environmental Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA In a recent paper (Zeng, et al., Numerical Meth. Partial Differential Equations, 1992, 8 451-467) a stable domain and boundary integral method was proposed for the solution of scattering and radiation problems of elastic structures submerged in an acoustic medium. No numerical examples, however, were included in this paper. This method is based on a unified energy-based variational functional which couples the domain representation for the interior deformable sWucture with the boundary integral equation representation for the exterior acoustic medium. Two important features of this finite element/boundary element representation are that (1) it is valid for all frequencies of excitation, and (2) upon discretization, it results in a symmetric algebraic system of equations. We present results of a numerical study conducted in order to test the validity and accuracy of the proposed methodology, for a variety of two- and three-dimensional prototype problems involving two-dimensional thin and thick shells and three-dimensional rigid scatterers. One important observation is that thick shell theory is inadequate for modeling the scattering problem for an annular infinitely long elastic cylinder even of moderate thickness. Full two-dimensional elastic theory is needed to represent this structure. The accurate results and excellent convergence demonstrated from the various numerical experiments suggest that the proposed coupled finite element - boundary element method can be used for solving a wide range of engineering interface problems in a systematic and robust way.

Key words: Structural acoustics, boundary integral equations, finite elements, symmetric coupling, boundary element method. fluid or solid, whereas domain FE methods are used to represent the regions which may include, for example, inhomogeneous or nonlinear media. The main purpose of introducing the BIE method to represent the homogeneous regions of a structure is to reduce by one the dimension of the subproblems associated with these regions. BIE methods are particularly effective for problems defined on infinite regions since the corresponding numerical discretizations are performed on finite boundary curves or surfaces. There is a vast literature on this subject. See, e.g. Refs 1 and 2 for theoretical developments, and Refs 3-10 for a number of different applications. The majority of available F E and BIE coupling procedures to date are based on various collocation methods.

1 INTRODUCTION In this paper we are concerned with the problem of scattering by deformable elastic structures in an acoustic medium, and its solution via a coupled finite element (FE) and boundary inlLegral equation (BIE) method. This is intended as a prototype for a class of more general interface problems that occur in such varied disciplines as elastostatics, elastodynamics, electromagnetics, etc. In coupled FE and BIE methods, the BIE methods are generally used to represent homogneous regions, either by L. J. Gray & J. H. Kane. ~Communicated To whom correspondence should be addressed. 79

80

X. Zeng, J. Bielak

Drawbacks of such methods are that the coupling destroys the symmetry of the stiffness matrix derived from the regions represented by FE methods, which usually is responsible for the majority of the unknowns in the final system of equations, and that the transition conditions on the interfaces of different regions are not natural. One well known difficulty with most of the existing BIE formulations is that they suffer from the nonuniqueness or nonexistence of solutions at a discrete set of frequencies, called critical or characteristic frequencies. Near these frequencies the discretized equations, while not singular, become ill-conditioned. Consequently, at high frequencies for which critical frequencies become closely spaced, the numerical results are unreliable. This problem is due to the failure of the boundary integral representation of the domain governing conditions. In the coupled FE-BIE method, the governing equations of a BIE subdomain is represented by BIE defined on the subdomain's boundary. However, when the boundary itself is a nodal line (or surface) of an eigenfunction of the integral operator, the operator becomes singular and thus fails to represent the associated domain governing conditions. This failure of the BIE at critical frequencies was first noted in Ref. 11, and later in Ref. 12. To circumvent this deficiency, various procedures have been developed in the past 30 years. A discussion of these developments is given in Ref. 13, and a comparison of numerical results corresponding to several of these methods is presented in Ref. 14. In Ref. 15 a new energy-based F E - B I E methodology which is both valid at all frequencies and leads to a symmetric formulation was proposed for the analysis of structures submerged in an acoustic medium and subjected to steady-state harmonic excitation, either in the form of applied loads (radiation) or incident waves (scattering). The main objective of the present work is to assess the validity and accuracy of this methodology via numerical examples. In the next two sections we present a summary of the theory and of the discretization procedure. Additional details may be found in Refs 15 and 16. Various numerical examples involving twodimensional shells of circular and elliptical shape and three-dimensional spherical and ellipsoidal rigid scatterers are considered subsequently.

2 STATEMENT OF PROBLEM AND VARIATIONAL FORMULATION Consider the geometry of Fig. 1. The bounded region fte in R 3 is occupied by a linear, isotropic, inhomogeneous elastic solid with density Pe and Lam6 material parameters #e, Ae. The exterior region f~f is an infinite space filled with a compressible, inviscid homogeneous fluid with constant material properties cf, the speed of sound in the fluid, and pf, the density of the fluid, r represents the interface between the two media. We assume that there is an incident steady-state harmonic fluid motion

J

Fig. 1. Model of elastic obstacle immersed in a fluid fullspace. given by a pressure P°(x, t) = Re[p°(x)exp(iwt)], where w is the frequency of excitation. We denote the total pressure in the fluid by Q(x, t) = Re[q(x) exp(iwt)], and the scattered pressure by e(x, t) = Q(x, t) - P°(x, t) = Re[p(x)exp(i~t)], and use U ( x , t ) = Re[u(x)exp(iwt)] to represent the resulting displacement field within the elastic solid. The linearized fluid-structure interaction problem then consists in finding p and u such that:

Ap + k2p --_ O,

k2

0";2

= c-T,

V. o'[u] q- peW2U= O,

in

in

u-.n =

on

1 (p+ +p0), pf ~d2

(la)

f~e

(o'[u-].n)-n = - ( p + +p0), (tr[u-].n) x n = 0 ,

f~f

(lb)

on

r

F

on

(lc) (ld)

F

p satisfies a radiation condition in ~ f

(le) (If)

In these equations k~ is the wave number, and tr[u] ----#e[VU + (Vu) r] ÷ Ae(V" u)I

(lg)

is the stress tensor within the solid, n is the outward unit normal to F, the superscript + denotes limit from f~f, and the subscript n denotes normal derivative. Equation (la) is the Helmholtz equation that governs the pressure in the fluid and (lb) is the reduced elastodynamic equation governing the motion of the inclusion. Equations (lc) and (ld) are the two traction transition conditions in the normal and tangential direction, respectively, while (le) represents the normal displacement transition condition on the interface F. To derive a variational principle corresponding to eqn (1) following the procedure described in Ref 15 we first construct an appropriate functional and then proceed to show that the vanishing of its variation yields (1). We start with the generalized lagrangian functional: nl[u,p]

=

5

[¢[u] : v u - po

2ua]dne

e

If 1 ,2 2 + ~J,¢ ~ tKYq - (Vq)2]df~f -

-

(2)

Fluid-structure interface problems The first two integrals o:a the right side of eqn (2) represent the total potential and kinetic energy of the interior elastic solid and the exterior fluid media, respectively, while the last integral represents the total background energy, i.e. the potentia]i and kinetic energy of the fluid due to the incident wave p0. Thus, II 1 denotes the difference in energy due to scattering. It is easy to show that by setting the first variation of the functional II 1[u,p] to zero for arbitrary 6p and 6u (subject to sufficient regularity) will lead to the original problem defined by eqn (1). Functional r[ l involves integrals over the infinite region f~f, which is impractical for numerical approximations. To reduce this problem to a finite domain we proceed as follows. If one integrates by parts the second and third integrals on the right of eqn (2) and makes use of the divergence theorem the domain integrals over f~f can be eliminated. There results:

II2[u,p) =

We recall the following jump relations for smooth ~b and F:

= S[ l(x)

(7a)

= :r ½ (x) +

(7b)

0

anx 5g[~b]+(x) = 0

Onx

(7c)

+lq~(x) + N[q~](x)

=

(70)

Here S, D, N, and M are integral operators on F, which satisfy the following symmetry relationships: Ir S[~b](x)¢(x)dF =

Ir S[¢](x)¢(x)dr

(8a)

Iv D[~b](x)~b(x)dF =

JrN[¢](x)¢(x)dr

(8b)

~1 1ft. [~r[u] : Vu - petO2u2]dQe

+ Ir(p ° + p + ) u -

• ndF

- Ir 2-~p+(p+n + 2p°)dF

(3)

in which the total pressure q has been replaced by p0 + p. In writing eqn (3) we have used the fact that p0 is a solution of Helmholtz equation (la), and have assumed that p satisfies the governing equation (la) plus the radiation condition (1 f). To actually ensure that (la) and (If) are satisfied, we use an indirect integral approach and express p in f~f as a linear combination of a single- and double-layer with density ~b, i.e. p = 6e[~b] + ag[~b]

(4)

ila af

where 5P[~b](x) - Ir ¢(y)G(lx - YI)dF

(5a)

...~[~](X)----fr, J t~(Y)--~-0: un, G(lx- yl)dF

(5b)

a is a constant with nonzero imaginary part, and G(z) is the fundamental solution, or Green's function, for eqns (la) and (lg),

G(z) =

81

1 - ~-~zeXp(ikfz)

in

R3

(6a)

That is, S and M are self-adjoint and D and N are adjoints. We remark that S has a weak singularity, D is continuous in R 2 and has an integrable singularity in R 3, and M is a hypersingular operator. From eqns (4) and (7) it follows that p+ = S[~b] + a(-½~b + D[~b])

(9a)

aM[~]

pn+ = 14> + N[~b] +

(9b)

An important characteristic of eqns (9a) and (9b) is that either equation can be solved for ~b uniquely for arbitrary p+ and Pn, + respectively, for all wave numbers kf, provided that Im(a) ¢ 0. The proof of the first statement is given in Ref. 18 and that of the second in Ref. 1. At this point, the functional associated with the indirect approach (4) could be obtained simply by inserting eqns (9a) and (9b) into eqn (3). This, however, would lead to an undesirable triple integral because of the term P+P+nin eqn (3). A remedy to circumvent this difficulty is obtained by first replacing pn+ with u and p0 using the transition condition (le). This leads to

p+p+=p+[pfcv2(u-.n)-p°],

on

F

(10)

We then consider eqn (le) as an additional constraint for u and Pn, + and introduce a Lagrange multiplier ~b to incorporate this constraint into the functional• This results in a modified functional

lIo

II 3 [u,p, g)] = ~

[tr[u] : Vu - pe~v2]dae

e

G(z)=4H(o2)(ikfz)

in

R2

(6b) + Ir(p ° + p + ) u - . nd£

Thus, p in eqn (4) auto:matically satisfies eqns (la) and (If) for arbitrary ~b. The reason we choose to represent p in f~f by eqn (4) rather than by either a single- or a double-layer by itself, is to preclude the occurrence of critical frequencies, thus ensuring that our procedure will be valid for all frequencies. The representation (4) was proposed independently in Refs 17 and 18.

l__}___p+ - Ir2p

'I

+-2 r

+p° ldr

' +p °)] Cdr -~(p+n (11)

X. Zeng, J. Bielak

82

Replacing p+ and p+ in eqn (11) by (9a) and (9b) then yields

lIo

H 4 [u, ¢, ¢] = ~

normal and tangential displacement components of the interior solid, i.e. tnn = (tY-. n). n, t;r = (tr-. n)- "r, Un = U- .n, u~- = u- .'r.

[tr[u] : Vn - peW2U2ldf~e

e

By setting 6114 to zero for arbitrary 6u, 6¢ and 6¢ there results:

4- Ir{p ° -4- [S[¢1 + a ( - ½ ¢ -4-D[¢])]} × (u-. n)dr

V.tr[u]+peW2U 2 = 0 , tnr = 0,

- Jr 2 p + {S[¢1 ÷ c~(- ½¢ ÷ D[¢])} o

(13b, c,d,e)

in f~e

(14a) (14b)

on F

2tnn -4- S[¢] -4- a ( - ½ ¢ -4-D[¢]) ÷pO _ ¢ = 0,

on F (14c)

1

x [~2(u- • n) +pn]dF + ~ Jr {u- "n 1 --(½¢

pfw2

o2 (½~ + N[¢] + c~g[¢] +pO) Cdr

6II4(u, ¢, ¢ ) = -- [~ ( V . o'[n] -4- pew2n) • 6udQ¢

J

o

(½pO+D[pO] +

aM[pOD=

0,

on F (14e)

Now, since p in f~f is specified by the ansatz (4), eqn (9) gives p+ and p+, with which eqns (14c) and (14d) reduce on 1-' on

(15)

F

(16)

Inserting eqn (16) into (14e) yields !2 (¢ +p0) -4-D[¢ ÷p0] -4-aM[¢ ÷p0] -4- S[p+]

1f

1 [½¢ + D[¢] + aM[e]]

+S[un + a [ - ~ (Un

+~2

1

1 (p++p0)_Un=0,

+ ~[_1¢ + D[¢]I +pO _ ¢}6uXdr

Jr(

1

pfw2

)-j[½¢+u[¢]

1

1 (I¢+D[¢I+aM[¢])+S[Un-p-~P° 1

2tnn +p+ +p0 _ ¢ = 0,

+ f [t;n6Un-4-t;~6u;]dF F

+

F

to:

d~t e

-4-p°n] - .x

on

(lad)

(12) This functional will serve as the basis for our variational principle. Its main attractive feature is that it satisfies the transition conditions automatically, and that it leads upon discretization to a symmetric formulation. This may be seen by examining the first variation of H4. After making use of the divergence theorem and the adjointness of D and N, and the self-adjointness of S and M, eqn (12) may be rewritten as:

-4- N[¢] -4- aM[el ÷pO) _ u; = 0,

p__~ pO)

+ pfofl [ pO + D[pO] + aM[pO]] 6¢dr (13a) where tnn and tn~ are, respectively, the normal and tangential components of the surface traction vector t n - tr-. n, whereas un and u~ represent the surface

.÷O~[_ 2PB 1 + -4- N[p+]] = 0

(17)

Thus far, the vanishing of the variation 6174 ensures that eqns (14a), (14b), and (16), that is, the second, fourth and fifth equations in (1), are satisfied. In addition, eqn (la) and the radiation condition (If) are automatically satisfied for any density ¢ in view ofeqn (4). It only remains to show that the normal stress transition condition (lc) is also satisfied, which may be done by using relations (17) and (15) as follows. Since p+ and p+, as given in eqn (9), are the limits from the solution p in ~f given by eqn (4), they satisfy automatically the integral representation formulae: ½p+ -4-D[p +] - S[p+.l = o, 1 + 2pn-N[p+n]+M[p+]=O,

on F

(18a, b)

With these relationships, eqn (17) reduces to

½A÷D[A]÷aM[A]=O, on

F

A-p+ ÷p°÷¢, (19a, b)

Fluid-structure interface problems It is shown in Refs 1 and 16 that eqn (19a) has only the zero solution provided I m a ~ 0. Thus eqn (19b) leads to p+ + p0 + ~b = 0

(f)

(20)

It immediately follows firom eqns (20) and (15) that

p+ +p°+tnn=O,

~b=t~n

(21a, b)

Therefore, (lc) is satisfied. In addition, it is found that the Lagrange multiplier ~, physically represents the normal traction of the interior material on the interface 1-'. Thus, it has been shown that a triad [u, ~b,~b] for which ~ I I 4 = 0 is a solution of the problem defined by eqn (1). That the converse is also true may be shown from eqn (1) using (13) and (18). Then, one has the following result:

Variation Principle VP. Given the functional II4 defined by eqn (12), the triad [u, 0,~b] is a solution of the problem defined in eqn (1) if and only if the variation 61-I4 vanishes for arbitrary variations 6u, The normal traction t~n on 1-' is given by ~b, and the scattered pressure p in f~f may be obtained from eqn (4) once the density, ~b, has been determined. Remarks: (a) The transition conditions are natural in this formulation. Consequently, the approximating functions for the interior solid and the exterior fluid need not satisfy any constraints on the interface I'. (b) The variational principle may be used to obtain finite and boundary element approximations for the coupled problem defined in eqn (1). The resulting discretized equations will be symmetric provided real-valtled basis functions are used. (c) The functional 1-114contains the duality pairing J'r~PM[~b]dI' that contains the hypersingular operator M. For performing computations this pairing may be evaluated readily by using the following identity: 19

r ¢M[¢]dF = '~ Ir S[•n] • (•n)dF - Ir S[n x V~]. (n x r e ) d r (22) which involves tangential derivatives of ~b and ¢ and has only a weak singularity. This makes any special treatment completely unnecessary. (d) A rigorous analysis of the above method, showing existence, uniqueness, and optimal convergence, is provided in Ret: 20. (e) The variational formulation is valid for all wave

(g)

(h)

(i)

83

numbers for which the problem defined in eqn (1) has a unique solution. I" has been taken to be smooth in this derivation. If r has edges or comers the limiting values on I" for the integral representation for p in f~f can be modified readily and the resulting equations will still be symmetric. In fact, the current formulation can still be used provided any edges or comers occur only at element interfaces. 16'21 The variational formulation VP remains valid if the material in the elastic body is anisotropic. Only the constitutive relations (lg) need to be modified. The exterior region f~f has been taken as a homogeneous full-space in the present formulation. With a simple modification, ~2f may also be a halfspace or a homogeneously multi-layered space. The only require-ment for the exterior configuration is the availability of the corresponding fundamental solution. An alternative variational procedure based on the direct boundary integral approach may be obtained by first considering a linear combination of the two equations given in eqn (18) as an additional constraint and then incorporating this constraint into functional II2 through a Lagrange multiplier ~b. This leads to a new functional

++

xj.

IIs(u,p ,p~, 4)) = ~

{~[u] : VII -- peW2u2}d~e e

+ Ir(p ° + p + ) u - .ndF

- Ir 2-~p+(p+ + 2p°)dI"

-

2llrp-~{(½P++D[P+] S[pn+]) + a(½p+~- U[p+~]

+ M[p+])}¢dr

'

(23)

in which p+ is considered as an independent unknown and a is a constant with nonzero imaginary part. Again, the reason for choosing a linear combination of the two conditions in eqn (18) as the additional constraint rather than either of them is to preclude the occurrence of critical frequencies. 22 It can be easily shown, by following the same steps as those in the indirect approach, that the set [u,p+,p +, ~b] is a solution of the problem defined in eqn (1) if and only if the variation 6II 5 vanishes for arbitrary variations Au, 6p+, 6p+, and 64~. u is the solution within f~e; P+ and p+ are the exterior limiting value of p and Pn on F; and ~b, the Lagrange

X. Zeng, J. Bielak

84

multiplier, turns out to be the density of the layer potentials in eqn (4) which again yields the solution for p in 9tf. One advantage of the direct formulation is that with a minor change, the functional H5 may be easily extended to multi-domain problems to obtain a stable domain decomposition functional which upon discretization leads to a sparse system of algebraic equations. 23

may be approximated as: M~

N~

= Z

=

k=l

N~ erh (x) = Z ¢ ~ ( X ) ,

(25)

k=l

in which

and

are u n k n o w n nodal point values.

Next substitute eqn (25) into (24c) and perform the integrations shown. One arrives at the corresponding approximation,

3 FINITE ELEMENT DISCRETIZATION We consider the numerical implementation of the indirect variational problem. The corresponding implementation for the direct variational problem is entirely analogous. The indirect functional (12) may be divided into two parts:

r • 1 T 1] 4 [U , ¢, ~b] = ~ Vr Krvr

+ fTvr

(24a)

where

lIo

1"[4~° [U] ---- ~

(26a)

where v r and fr represent, respectively, the total unknown and the effective load vectors on F, Vr - [Urh, Crh, Crh]T,

H 4 [u, q~,/~] = H4~ [u] + 1-I4F [ u - , ¢, ~b]

k=l

fr = [fur,f*, f~]T

(26b, c)

and Kr is the effective stiffness matrix of the coupling integral over F, K~r¢

(O'[U] : VU -- p e W 2 U 2 ) d a e

(24b)

Kr -

e

represents the total energy of the interior elastic solid, and

[s[¢] + a(-10 + D[¢])]}u~dr Hr[u-, ¢, ¢] - I{p0+ r

× [p~2u~ + p°]dr

• 1 T * [u]=~u~K~u~ -

(24c)

is the energy corresponding to the exterior fluid and the interaction at the interface between the two regions. The functional Hr [u-, ¢, ¢] will be subsequently referred to as the indirect coupling integral, or simply the coupling integral. The numerical implementation for the interior energy functional H4~°[u] may he performed by a standard finite element discretization procedure which is completely independent of the numerical methods for the coupling integral. Our discussion will focus on the coupling integral H4r [u-, ¢, ¢]. To derive the matrix form of the coupling integral II4r[u-, ¢, ¢], we introduce two finite dimensional subspaces U~ and S ~ in F depending on the mesh size h. U~ has the real basis of the form a h, a~,..., a~,, and S h has the real basis of/31~,/32~,...,/3Nhh, in w ~ h /3k~ are at least continuous, piecewise linear prolynomiats due to eqn (22) while a~ must be selected such that u is continuous on F. With these basis functions, u-, ¢ and ¢

(26d)

(27a)

1 2 u~w~fiu~ T-. -1 T ~w = ~ufiKfiu~

where K~ and M~ are, respectively, the usual stiffness and mass matrices of the interior structure, and ufi = [uf2,Ur]T,

1 ~ (½¢+N[0] + ~1Jr [Un - p~o--+ aM[e] + p0)] e d r

0 LKo~ur K¢¢

The matrix form approximation for the interior energy functional can be obtained from standard FE procedures, H4

- Jr 2-~[S[¢] + a(-½¢ + D[¢])]

IKcu°r

Kfi =

?uo 1 Kur,

Kuar Kurr

(27b, c)

are, respectively, the nodal displacement vector of the discretized region (~ = f~ U F, and the effective partitioned dynamic stiffness or impedance matrix of the interior structure. To obtain the resulting system of linear equations, substitute eqns (27) and (26) into (24), take the differential of the resulting expression for lI4, and set to zero the coefficients of each dUk h, d~b~ and d¢~. This will lead to a symmetric system of linear algebraic equations,

Kuur o o

//o/

Kuru, Kurur Kur¢ Ku~¢ [~

KCur K¢~r

0 K¢¢

K0~¢ J

Ur ~b ¢

=

fu~ f¢ f, (28)

The submatrices Ku, u~, Kunur, Kuru~, and Kurur in eqn (28) are sparse as they represent partitions of the standard impedance symmetric matrix K ~ - w2Mfi. Kur¢ and K~r are small and also sparse as they involve only products of the approximations for the normal

Fluid-structure interface problems displacements and tractions on F. Kur¢, KCur, K0¢, and K ~ , on the other hand, are full but small since ~b and ¢ are defined only on the boundary F. These properties can be used for the efficient numerical solution of eqn (28). The vectors fur, and re, f~ are effective generalized forces and displacements applied on F, and are functions of the prescribed driving field p0 and its normal derivative/'~, on F. It is quite easy to verify that the system (28) is symmetric, that is, K T = K. For instance, the relationship, K ~ = K ~ , may be easily verified by inserting the expressions ~b(x) and ~(x) given in eqn (25) into the coupling integral expression (24c) and in view of relations (8) to obtain

= Ir{½ /(x) +

+

og.ML~i ]x) }

x ~(x)dr(x)

Jr{½ (x) + x ~ ( x ) d r ( x ) = (K~)ji

85

possible, also, to represent the density ~ of the layer potentials by discrete point sources of unknown intensity, in which case all the double integrals reduce to single integrals of the same form as those that occur in exterior problems. In this case, special techniques such as those described in Refs 19, 24 and 25 are needed to evaluate the hypersingular integrals. (k) Various iterative solution schemes may also be effective for solving the coupled system of algebraic equations. It has been shown in Ref. 21 through numerical examples for rigid scatterers that the resulting system of linear equations are well-conditioned provided I m ( a ) ~ 0. This may result in rapid convergence of iterative solution methods when applied to the coupled fluid-structure problem. 4 NUMERICAL EXPERIMENTS

(29)

Remarks:

In order to assess the stability and accuracy of our procedure, we first consider several two-dimensional

(J) Several double integrals must evaluated for con-

r/a= 1

structing the submatrices appearing in eqn (28). The added computational effort is not large as only a few Gaussian integration points are sufficient for the second integral even when the source point and the observer involved in the integral operators are in the same element. It is

r/a= 10

2.0

(a) 1.5

Ca)

Exact Results Shell Results Elastic Results

1.0

0.5

0.0

(c)

(d)

(e)

(0

1.5

1.0 a_

0.5

I

0.0

I

I

I

13

1.0 0.5 O.C

Fig. 2. T y p i c a l finite e l e m e n t m e s h e s f o r c i r c u l a r a n d elliptic (a/b = 2) shells.

0

2

4

6

~ 8 0 2 4 Wave Number, ka

6

' 8

10

Fig. 3. N o r m a l i z e d a m p l i t u d e o f s c a t t e r e d p r e s s u r e a t v a r i o u s l o c a t i o n s o n a c i r c u l a r cylindrical shell a n d w i t h i n the fluid as a f u n c t i o n o f w a v e n u m b e r : (a) r/a = 1, 0 = 270°; (b) r/a --- 10, O = 270°; (c) r/a= 1, 0 = 9 0 ° ; (d) r/a= 10, O = 90°; (e) r/a = 1, 0 = 0 °, 180°; (f) r/a = 10, 0 = 0 °, 180°; h/a = 0"01.

X. Zeng, J. Bielak

86 2.0

2.0

1.5

1.5

(a) ~

1.0

i

[i

; '.*

i

~,

"~

1.0

"~""*%

0.5 II 0.0

£ i

Lo i

::i

~

--

Shell Results

3.0 ~

2.0

~-.

ii

1.0

~

°

s*'~'~'*" ~

.~¢,,,~

t~

*

* Elastic Retmlts

o

¢J r,t}

:

0.0

ii

i

i

i

%"

I

t~

.o o ¢J O~

0.0

'

I

I

(c) 1.5 D

1.0

~%",

1.0

0.5

~

0.5

0.0

0

- - - - - - - - - ~ ~ - - " " " " ~ " - - 2 4 6 8

10

Wave Number, ka

0.0

0

2

4

6

8

10

Wave Number, ka

Fig. 4. Normalized amplitude of scattered pressure at various locations on a circular cylindrical shell as a function of wave number: (a) 0 = 270°, (b) 0 = 90°, (c) 0 = 0°, 180°; h/a = 0.1.

Fig. 5. Normalized amplitude of scattered pressure at various locations within the fluid as a function of wave number: (a) 0 = 270°, (b) 0 = 90°, (c) 0 = 0°, 180°, h/a = O"1.

fluid-cylindrical shell interface problems. The configurations investigated here include circular thin and thick shells. We also consider the case of a shell with elliptical cross-section. In all the examples the axis of the cylinder is perpendicular to the incident wave front. We consider first circular thin and thick elastic shells immersed in a fluid and subjected to an incident plane wave p°(x) =piexp[-ikf(xsin• + y c o s 0 i ) ] . pi and 0i represent the amplitude and the angle of inclination of the incident wave, respectively. The shell is modeled in two different ways. In one approach, it is represented by the general shell theory for both thin and thick shells described in Ref. 26. In the second approach, the shell is modeled as an annular continuous isotropic elastic body in plane strain. Both problems are solved by the coupled FE-BIE procedure described in the previous sections. (Details for the coupling treatment between the shell approximation and the acoustic medium can be found in Ref. 16.) The exact analytical solution of the problem when the shell is considered as a hollow cylinder is given in Refs 27 and 28. This solution is presented here also, for comparison. In all our calculations the shell interior is taken to be

the general shell theory is given in Fig. 2. One layer of eight-noded quadratic isoparametric elements is used for the annular formulation while three noded general shell isoparametric elements are used for the shell formulation. In all the calculations, three noded quadratic isoparametric elements are used to interpolate the interface F, the Lagrange multiplier ~p, and the density of the layer potential ¢~. Between 16 and 64 elements are used along the periphery depending on the value of the wavenumber. The numerical results for a circular thin shell are presented in Fig. 3, The computations are performed using a radius-to-thickness ratio a/h = 100. The ratios of mass densities and wave speeds of the shell and fluid are given by Pe/Pf = 7"65 and ce/cf = 5"97. The dotted line, the thin solid line, and the thick solid line correspond, respectively, to the analytical results, the results from two-dimensional general shell theory, and the hollow elastic model. The three plots on the left column show the total pressures at various points directly on the shell, for a range of values of the normalized wave number kfa q (0, 10] and with a step-size of 0.1. The plots on the right solumn give the scattered pressures at the same angular locations as those of the corresponding left column, but on a circle of radius r = 10a within the fluid.

in vacuo. A typical finite element mesh for the circular shell structure represented either as a hollow cylinder or by

Fluid-structure interface problems

87

3

(a)

o

O r~ I.

u

Co)

e~

eL t~

2

"O

1

~

~

~

t~

"I3

1

r/)

"O e¢

0

J

- - TotalPressure 2 1 ~ 1 .......... Scattered....~Plt:ssure

I

L~ r.q

0

;

"O

(¢)

J

4

6

(e)

- - r/b= i00 .......... r,/b= 2(~

.-"

b,,

J

Z O

'.. ............

0" 0

2

'"

4

.. . . . .

6

$

10

Wave Number, ka

0

0

2

8

10

Wave Number, ka

Fig. 6. Normalized amplitude of total and scattered pressure at various locations on an elliptic cylindrical shell as a function of wave number: (a) 0 = 2 7 0 °, (b) 0=.90 °, (c) 0 = 0 °, 180°, a/b = 2, h/b = 0"1, ¢ = 0°.

Fig. 7. Normalized amplitude of scattered pressure at various locations within the fluid as a function of wave number: (a) 0 = 2 7 0 °, (b) 0 = 9 0 °, (c) O=0°, 180°, a/b=2, h/b=O'l, 0' = 0 °.

The results in Fig. 3 serve to demonstrate the robustness and accuracy o f the proposed coupling technique, as the numerical results are practically indistinguishable from the exact ones for all the frequencies considered herein. In addition, it is clear that shell theory provides an excellent approximation to the more accurate, but more expensive, hoUow continuum model when the shell thickness is small compared to its radius. To assess how the two models work for thick shell structures, Figs 4 and 5 depict, respectively, the amplitude of the normalized back-, forward- and sidescattered pressure, directly on the shell and at a distance of ten radii from the shell, respectively. The material properties in this case are the same as those for the thin shell, but the radius-to-thickness ratio is taken as a/h = 10. Figures 4 and 5 show that unlike the thin shell case, the results from the shell model depart signficantly from the analytical results, while the results from the elastic model remain practically identical to the exact ones. There are several possible reasons for the failure of the shell model in the thick shell case. First, in the shell model, the interaction between the shell and the fluid is approximated with the displacement of the mid-plane of the shell; this causes errors as the thickness of the

shell increases. Secondly, as the thickness of the shell increases, the assumption from shell theory that straight lines remain straight is no longer valid. Finally, high modes may also contribute to errors. In fact, even low frequencies excite higher modes both along the periphery and in the radial direction; the former has the effect of making the wavelength in the shell shorter, thus relatively thicker, while the latter produces displacement of the shell along the radial direction, in contradiction to the linear displacement assumption of shell theory. The small deviation at higher frequencies between the hollow cylinder model and the corresponding analytical results in Figs 4 and 5 may be due to the fact that as the frequency becomes high, one layer of eight-noded quadratic elements along the radial direction is insufficient for representing the radial displacement. Also, additional elements may be needed in the angular direction. By comparing the results from the thin and thick shell cases, it can be seen that while the thin shell pressure responses are quite smooth with respect to the wavenumber, the responses from the thick shell structures show increased fluctuations. Each peak corresponds to a physical resonance of the shell-fluid system. In fact,

X. Zeng, J. Bielak

88

2.0

(a)

1.5 -

--

TotalPressure .......... Scau~redPressure

1.0 J

0.5 / 0.0

.... ,.....-, ......................... ',..,..........

2.0 -

1.5

.,:

..

.

.

.

.

.

(b)

'..

..,.

l.O

~

13.5-

,/

:"

/

~ Fig. 8. A typical finite element mesh for a three-dimensional

i

0,0"

•~.

1.5 -

Z-

~.o

(c) --------

sphere. it has been pointed out in Ref. 29 that even for the thin shell cases, there is a large number of resonant wave numbers in the range of (0,10). However, the peaks are so highly localized that an extremely small stepsize (of the order of 10 -8 ) is needed for the peaks to become manifest. To investigate further how the shell response is affected by the complexity of the shell geometry,

0.50.0 ,......

0

i 2

4

i 6

r 8

10

Wave Numberka Fig. 9. Normalized amplitude of total and scattered pressure at

p/pi

Table 1. Total pressure generated at various surface locations by a rigid sphere scatterer, for three critical wavenumbers, and various numbers of elements

various locations on a sphere surface as a function of wave number, 7 i = 0°: ( a ) back-scatter, (b) forward-scatter, (c) sidescatter.

ka

Figs 6 and 7 show the total and normalized scattered pressures at several points located both directly on the surface of an elliptic shell and in the far field within the fluid. The results are calculated from the fluid hollow cylinder formulation. Except for the finite element mesh for the interior elastic solid, which is shown in Fig. 2, all the material properties, interpolation shape functions, etc., are the same as those for the circular shell cases. The incident wave in this example travels in the direction parallel to the semiminor axis of the ellipse (0i = 0°). The ellipse has an aspect ratio of a / b = 2 and a thickness to radius ratio a / h = 0.1. The ellipse is a rather simple geometry, yet mathematically, even for this simple shape it is already quite difficult to find the corresponding analytical solution. Such a solution does not seem to be available in the literature. Notice from Fig. 6 the various peaks corresponding to resonant frequencies of the shellfluid system. As seen in Fig. 7, these frequencies appear as well in the far field. It is of interest that the scattered wave attenuates as O((r/b)l/2), since in Fig. 7 the normalized scattered pressures at two different distances away from the scatterer practically coincide. Next we consider the scattering problem of a threedimensional rigid sphere immersed in a fluid full-space.

~

7(

No. of Elements 16(96) 24(288) 32(384) Exact

Back scatter

Forward scatter

Real

Imaginary

Real

1.7507 1.7544 1-7553 1.7553

-0-3986 -0.3979 -0.3978 -0.3977

-0-3897 -0.3818 -0.3804 -0-3792

Imaginary 1.0882 1.0876 1.0881 1'0884

ka = 27r No. of Elements

Back scatter Real

16(96) 24(288) 32(384) Exact

-1"8466 -1-8962 -1-8999 -1.9026

Imaginary 0.2908 0'2660 0'2656 0.2664

Forward scatter Real 0.8661 0.8211 0.8172 0.8157

Imaginary 0.8491 0.7814 0"7809 0.7744

ka = 3re No. of Elements 16(96) 24(288) 32(384) 40(600) Exact

Back scatter Real

Imaginary

1.9081 1.9578 1.9496 1.9481 1.9471

-0'1157 -0.1929 -0-1953 -0.1958 -0.1948

Forward scatter Real 1.2723 1.0347 0-9887 0"9739 0-9622

Imaginary -0"6658 -0.4732 -0.4484 -0.4428 -0.4407

Fluid-structure interface problems

89

1.0

--r/a = 20 ......... r/a = 100

0.8

(a)

0.6 0.4 0.2 t~ 0.0 4.0

,.........'"'"'

a.

..................... (b)

3.0

o

2.0 1.0 0.0

'o

(c)

0.8 0 Z

0.6 0,

¥

0

2

4

6

8

I0

Wave N u m b e r k a

Fig. 10. Normalized amplitude of scattered pressure at various locations on two far field sphere surfaces as a function of wave number, 7i = 0°: (a) back-sratter z = -20, -100; (b) forwardscatter z = 20, 100; (c) side-scatterer x, y = 4-20, 4-100. The goal of the investigation for this simple case is, again, to assess the robustness and accuracy of the proposed method for three-dimensional problems using available exact solutions. Later we discuss the case of a rigid ellipsoid. A typical finite element mesh on the surface of the sphere is depicted in Fig. 8. In order to compare and assess the convergence through different mesh grids, the mesh is generated in such way that most of the elements have approximately the same shape and aspect ratio. For this numerical experiment, nine-node quadratic isoparametric elements are used to represent the boundary 1-', the Lagrange multiplier ~p, and the density of the single-layer ~b. Three by three Gauss quadrature points are used for each regular element. The integration scheme for singular elements (in which the source point coincides with the field point) is described in detail in Ref. 16. The incident wave is, again, taken as a plane wave: p0(x )

=

pi e x p

[-ikf(x c o s

oLi -{- y c o s / ~ i

+ z c o s ,yi)]

(29) where t~i,/3 i and ,,/i are, respectively, the three direction angles between the incident wave and the x, y and z axes.

Fig. 11. Three projections of a typical boundary element mesh of an ellipsoid: a : b : c = 1 : 0.5 : 2, (from top: yz-projection, zx-projection, xy-projection). The total pressure at several points directly on the scatterer corresponding to three critical wave numbers for which standard boundary integral methods are expected to fail, using different mesh sizes is shown in Table 1. The incident wave travels along the positive z axis, i.e. (ai,Bi,Ti) = (90°,90°,0°). The first column of Table 1 gives the average number of elements along any equator and, in parentheses, the corresponding total number of elements on the surface of the sphere. Results shown in this table indicate excellent accuracy and convergence, even for the antipole position, opposite to the point where the incident wave first impinges upon the scatterer, where the largest errors normally occur. As expected, the agreement with the exact solution is closest at the lower frequency for a fixed number of elements. The total and scattered pressures of the back-, forward- and side-scattering points directly on the scatterer as a function of wave number are shown in Fig. 9. The corresponding analytical results are not shown here because graphically, they are indistinguishable from the numerical ones. The normalized scattered pressure on the back-, forward- and side-scattering points on two spheres within the fluid with radii of, respectively, 20 and 100 is shown in Fig. 10. Note that the scale factor in this

90

X. Zeng, J. Bielak 2.5

6.0 ---

2.0

- - -

-- -

1.5

Total Pressure (back) Total Pressure (forward) S c a t m r e , d P r e s s u r e (back) Scattered Pressure (forward)

- - ---

4.0

(x,y,z) = (-20,-20,-20) (x,y,z) = (-100,-100,-100) (x,y,z) = (20,20,20) (x,y,z)=(100,100,100) ~ ' ~

~ "~" ~

1.0

2.0 (a) z-sides

0.5 ~

.0

,~

2.0

~

-"" P

I

i

0.0

(b) x - s i d e s

--- - -

1.0

.5,

"~

(a) b a c k - a n d f o r w a r d - s c a t t e r i n g

(x,y,z) = (-20,0,0) (x,y,z) = (-100,0,0) (x,y,z) = (20,0,0)

(b) x - s i d e s

1.0 0.5 0.5

.~--_0.0

~

....................... ;

~.o

I

[

-

./,,

I

r~

1.5

~

1.1)

--- - - - -

o

z 0.5

0.3

0.5 0.0 0

". ,,"

'\\./

(c) y - s i d e s 0.8

~

"~

0.0

i

i

l

i

2

4

6

8

10

(x,y,z) (x,y,z) (x,y,z) (x,y~)

/

0.0

(c) z - s i d e s

= (0,0,-20) = (0,0,-100) = (0,0,20) -- (0,0,100)

\

/

\

z

i

0

2

4

6

8

10

Wave Number ka Wave Number ka

Fig. 12. Normalized amplitude of total and scattered pressure at various 19catigns on an ellipsoid surface as a function of wave number, a ' = ff = 7 ' = 54"73°: (a) back (x,y,z) = (0,0,-2), forward (x,y,z) = (0,0,2); (b) back (x,y,z) = (-1,0,0), forward (x,y,z)= (1,0,0); (c) back (x,y,z)= (0,-0.5,0) forward (x, y, z) = 0, 0.5, 0). three-dimensional case is r/a, instead of (r/a) 1/2 as for the two-dimensional one. The agreement of the two sets of scattered curves verifies that, as expected, the far-field scattered pressure decays as (r/a) -l. As a final example of the application of the BIE methodology to three-dimensional scatterers we consider the case of a rigid ellipsoid with semi-axes ratios a : b : c = l : 0 - 5 : 2 (Fig. 11). Obviously, this geometrical configuration cannot be reduced to a twodimensional problem. The surface and far-field pressures at several points on the scatterer and in the fluid for the case in which the incident plane wave forms equal angles with the three coordinate axes are shown in Figs 12 and 13. The trends of the total and scattered pressures directly on the scatterer are in accordance with intuitive behavior. It is of interest that in the far field the forward-scattered pressure is much greater than the back-scattered one. Also, the attenuation of the response varies generally with the inverse of the distance from the center of the ellipsoid. Remarks: (1) During the course of obtaining the numerical results presented in this study we observed that

Fig. 13. Normalized amplitude of amplitude of scattered pressure at various locations in the fluid as a function of wave number, a i = fli = ? i = 54.730: (a) back- and forwardscattering, (x,y, z) = (+20, +20, ±20), (x,y, z) = (±100, ±100, ±100); (b) x-sides, (x,y,z)= ±20,0,0), (x,y,z)= (±100,0, 0); (c)z-sides (x,y,z) = (0,0,±20), (x,y,z) = (0,0,±100). the accuracy of the far-field results is significantly higher than that directly on the scatterer, i.e. far fewer elements were required to achieve a given accuracy in the far-field than directly on the scatterer. This, if confirmed for more complex structures may result in efficient numerical schemes for the analysis of the far-field response. Thus, in engineering applications where the main concern is in the far-field effects, relatively coarse meshes may adequately fulfill the practical requirements for accuracy. This property, however, may pose problems if the interest is focused on inverse problems. (m) Another important property revealed through the numerical experiments is that, unlike the rigid scattering problems where the accuracy of the results can be predicted by the number of nodes per wavelength in the fluid, the interface problems are more complicated. Other factors, such as the complexity of geometrical configuration and the wave speed in the structure is greatly responsible for the accuracy of the results.

Fluid-structure interface problems

5 C O N C L U D I N G REMARKS In this paper we have presented a set of numerical examples designed to illustrate the applicability of a recently developed coupled FE--BIE method for solving radiation and scattering problems that involve elastic bodies immersed in an infinite acoustic medium. Numerical results confirm that this procedure is robust and exact for all frequencies within discretization error. The greatest challenge for coupled FE-_BIE methods still lies ahead, however, due to the dense nature of the equations on the structure boundary. By combining integral methods with domain decomposition ideas or other sparsification techniques and domain finite elements, 23'3°-32 it seems likely that large threedimensional interface problems with complex geometry, which have eluded investigators for some time, now can become more tractable.

REFERENCES 1. Hsiao, G. C., The coupling of BEM and FEM - - a brief review. In Boundary Elements X, Vol. 1, ed. C. A. Brebbia. Springer, Berlin, 1988, pp. 431-445. 2. Johnson, C. & Nedelec, J. C., On the coupling of boundary integral and finite eleraent methods. Math. Comput., 1980 35, 1063-1079. 3. Bielak, J. MacCamy, R. C. & McGhee, D. S. On the coupling of finite element and boundary integral methods. In AMD Vol. 60, ed. S. K. Datta. ASME 1984, pp. 115-132. 4. Costabel, M. Symmetric methods for the coupling of finite elements and boundary elements. In Boundary Elements IX Vol. 1, ed. C. A. Brebbia et al. Springer, Berlin, 1987, pp. 411-420. 5. Han, H. A new class of variational formulations for the coupling of finite artd boundary element methods. J. Comput. Math, 1988. 6. Bielak, J. & MacCarny, R. C., Symmetric finite element and boundary integral coupling methods for fluid-solid interaction. Q. Appl. Math. 1991, 49, 107-119. 7. Hsiao, G. C., The coupling of boundary element and finite element methods. Z. Angew. Math. Mech., 1990, 70, 493503. 8. Jeans, R. A. & Mathews, I. C., Solution of fluid-structure interaction problems using a coupled finite element and variational boundary element technique. J. Acoust. Soc. Am., 1990, 88, 2459-2466. 9. Everstine, G. C. & Henderson, F. M., Coupled finite element/boundary element approach for fluid-structure interaction. J. Acoust. Soc. Am., 1990, 87, 1983-1947. 10. Ginsberg, J. H., Chen, P. T. & Pierce, A. D., Analysis using variational principles of the surface pressures and displacements along axisymmetrically excited disk in a baffle. J. Acoust. Soc. Am., 1990, 88, 548-559. 11. Lamb, H. Hydrodynamics, 6th edn, Article 290. Dover, New York, 1945. 12. Schenck, H. A., Iraproved integral formulation for acoustic radiation problems. J. Acoust. Soc. Am., 1968, 44, 41-58.

91

13. Colton, D. & Kress, R., Integral equation methods in scattering theory. John Wiley, New York, 1983. 14. Amini, S. & Harris, P. J., A comparison between various boundary integral formulations of the exterior acoustic problem. Comput. Methods Appl. Mech. Engng., 1990, 84, 59-75. 15. Zeng, X., Bielak, J. & MacCamy, R. C., Unified symmetric finite element and boundary integral variational coupling methods for fluid-structure interaction. Numerical Meth. Partial Differential Equations, 1992, 8, 451-467. 16. Zeng, X., Stable symmetric FEM-BEM coupling methods for fluid-structure interface problems with applications. PhD thesis, Department of Civil Engineering, Carnegie Mellon University, 1992. 17. Leis, R., Zur Dirichletschen Randwertaufgbe des Aussenraums der Schwingungsgleichung. Math. Z., 1965, 90, 205-211. 18. Brakhage, H. & Werner, P. Uber das Dirichletsche Aussenraumproblem ffir die Helmholtzsche Schwingungsgleichung. Arch. Math., 1965, 16, 325-329. 19. Maue, A. W. Zfir Formulierung eines allgemeinen Beugungsproblems durch eine Integralgleichung. Z. Phys., 1949, 126, 601-618. 20. Bielak, J., MacCamy, R. C. & Zeng, X., Stable coupling method for interface scattering problems, Research Report No. R-91-199. Department of Civil Engineering, Carnegie Mellon University, 1991. 21. Zeng, X. & Bielak, J. Stability assessment of a unified variational boundary integral method applicable to thin scatterers and scatterers with comers. Comput. Methods Appl. Meth. Engng., 1994, 111, 305-321. 22. Burton, A. J. & Miller, G. F. The application of integral equation methods to the numerical solution of some exterior boundary-value problems. Proc. Roy. Soc. Lond., A., 1971, 323, 201-210. 23. Zeng, X. & Bielak, J., Exterior stable domain segmentation integral equation method for scattering problems. Int. J. Num. Meth. Engng, 1994, 37, 777-792. 24. Chien, C. C., Rajiyah, H. & Atluri, S. H., An effective method for solving the hypersingular integral equations in 3-D acoustics. J. Acoust. Soc. Am., 1990, 88, 918-937. 25. Meyer, W. L., Bell, W. A., Stallybrass, M. P. & Zinn, B. T., Prediction of the sound field radiated from axisymmetric surface. J. Acoust. Soc. Am., 1979, 65, 631-638. 26. Bathe, K.-J., Finite element procedures in engineering analysis. Prentice-Hall, Englewood Cliffs, NJ, 1982. 27. Doolittle, R. D. & Uberall, H., Sound scattering by elastic cylindrical shells. J. Acoust. Soc. Am., 1965, 272-275. 28. Gaunaurd, G. C. & Brill, D., Acoustic spectrogram and complex-frequency poles of a resonantly excited elastic tube. J. Acoust. Soc. Am., 1984, 75, 1680-1693. 29. Kallivokas, L. F., Zeng, X. & Bielak, J., Stable localized BEM-FEM for fluid-structure interaction (in preparation). 30. Canning, F. X., The impedance matrix localization (IML) method for moment-method calculations. IEEE Antennas Propagation Mag., 1990, 32, 18-30. 31. Rokhlin, V., Rapid solution of integral equations of classical potential theory. J. Comput. Phys., 1985, 60, 187-207. 32. Zeng, X., Kallivokas, L. F. & Bielak, J., Stable localized symmetric integral equaion method for acoustic scattering problems. J. Acoust. Soc. Am., 1992, 91, 2510-2518.