A numerical study of a spectral problem in solid-fluid type structures

A Numerical Study of a Spectral Problem in Solid-Fluid Type Structures Carlos Conca Deparfamento de Ingenieria Matematica, Facultad de Ciencias Fisicas y Matematicas, Universidad de Chile, Casilla 1 70/3-C0rreo 3, Santiago, Chile Mario Duran* Centre de Mathematiques Appliquees, Ecole Polytechnique, 9 1 128 Palaiseau, France and Deparfamento de Ingenieria Matematica, Facultad de Ciencias Fisicas y Matematicas. Universidad de Chile, Casilla 170/3-C0rreo 3, Santiago, Chile. Received 30 September 1991; revised manuscript received I0 May 1994

This article presents a numerical study of a spectral problem that models the vibrations of a solid-fluid structure. It is a quadratic eigenvalue problem involving incompressible Stokes equations. In its numerical approximation we use Lagrange finite elements. To approximate the velocity, degree 2 polynomials on triangles are used, and for the pressure, degree 1 polynomials. The numerical results obtained confirm the theory, as they show in particular that the known theoretical bound for the maximum number of nonreal eigenvalues admitted by such a system is optimal. The results also take account of the dependence of vibration frequencies with respect to determined physical parameters, which have a bearing on the model. 0 1995 John Wiley & Sons, Inc.

1. INTRODUCTION

This article aims to solve numerically a mathematical model that describes the vibrations of a solid-fluid type system. The problem of determining the vibratory eigenfrequencies and eigenmotions of this type of system has considerable importance in engineering, as this occurs naturally in the design and simulation of various sorts of industrial equipment. In recent decades much effort has been devoted to experimental and theoretical research into this type of problem, examples of which include the articles of R. Blevins [1,2], H. J. Connors [3], S. S. Chen [4], [5], D. J. Gorman [6], M. J. Pettigrew [ 7 ] , M. Paidoussis [8], and J. Planchard [9, 10, 11,121. In particular, the model that we study describes the vibrations of a bundle of K metallic tubes immersed in an incompressible viscous fluid. Among the publications that study this problem using similar mathematical techniques *Please address all correspondence to Dr. M. DurBn, Centre de MathCmatiques AppliquCes, Ecole Polytechnique, 91 128 Palaiseau Cedex, France. Numerical Methods for Partial Differential Equations, 11, 423-444 (1995) 0 1995 John Wiley & Sons, Inc.

CCC 0749- 159W95IO40423-22

424

CONCA AND DURAN

and models, one might cite the joint works of F. Aguirre and C. Conca [13], C. Conca, J. Planchard, and M. Vanninathan [14,15, 161. From the mathematical point of view, the model in question is a nonstandard differential eigenvalue problem, which involves the stationary Stokes equations. The variational formulation of this model is a quadratic eigenvalue problem whose coefficients are self-adjoint linear operators, acting on an infinite dimensional Hilbert space. In C. Conca, M. DurBn, and J. Planchard [17], a result is proved of the existence of eigenvalues and eigenvectors of the model and it is shown that the spectrum of the problem consists of a countable infinite quantity of complex eigenvalues, which converge in modulus to infinity. Furthermore, it is shown that the number of eigenvalues with a nonzero imaginary part is finite. An estimate of the number of non-real eigenvalues is provided in [18], where it is shown that the spectrum of the solid-fluid system allows at most 4 K imaginary eigenvalues. The mathematical techniques used to prove the above consist, in the first case, of identifying the unknown eigenvalues with the characteristic values of a non-self-adjoint compact operator, and using classical spectral theorems of functional analysis. In the second case, to obtain the estimate, we use abstract results concerning the variation of the spectrum of a bounded linear operator perturbed by a non-negative operator, of finite range. More details can be found in the above cited works and the references contained therein. In the Section I1 we present a summary of the above. The problem is then solved numerically. To this end a mixed variational formulation of the differential problem is introduced, in which the condition of incompressibility of the fluid, or of zero divergence, is dealt with implicitly. Then, the problem is approximated by the discretization of the Hilbert spaces involved. We approximate the velocity and pressure through continuous functions, which on each element of the triangulation of the domain are degree 2 and degree 1 polynomials, respectively. The basis chosen for the finite dimensional spaces are those that are usually used in the finite elements method. The family of triangulations of the domain is constructed in such a way as to satisfy a pair of conditions of nondegeneracy and regularity. The discrete problem becomes, equivalently, a quadratic eigenvalue matrix problem. Using classical techniques, we show that this can be reduced to a generalized eigenvalue problem of the Sylvester type, which can be dealt with numerically using standard computer software packages for numerical calculations. The numerical experiments were carried out in a test geometry with a single tube, using three domain triangulations. The numerical evidence collected confirms the previous theoretical results, in the sense that a countable spectrum was obtained, which converges to infinity and which admits nonreal eigenvalues, and which can in fact be found in the region of the complex plane indicated theoretically. Furthermore, this numerical experience provides valuable information about the spectral variation of the problem. In particular, it clearly shows the dependence that exists between the imbedded rigidity of the bundle of tubes (which we represent by k) and the nonreal part of the spectrum of the problem. In this context, the numerical experiments take account of the bifurcation to the complex plane of a set of real eigenvalues, which collapse together for certain threshold values of k. The insensitivity with respect to k of the higher eigenvalues should also be noted. The above can be seen from simple inspection of the result tables. Finally, the graphs are presented of velocity fields and isobaric lines of some of the structure’s vibratory eigenmotions.

. . . A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

425

II. MATHEMATICAL FORMULATION OF THE PROBLEM

In this section we present a summary of the deduction and theoretical study of the mathematical model, which describes the physical problem in question. Firstly, the physical hypotheses assumed for the system are established, and the mathematical model that this represents is deduced. Then, the differential problem is rigorously formulated, and finally, the most important theorems of existence and location are summarized. A. The Physical Problem

The physical problem that interests us is the study of the vibrations of a solid-fluid type system. To be more exact, the problem consists of determining the vibratory eigenfrequencies and eigenmotions of a bundle of metallic tubes immersed in an incompressible viscous fluid. The fluid is assumed to be contained in a three-dimensional cavity with rigid walls. It is assumed that the walls are parallel to each other, that they are perfectly rigid (they do not allow deformations) and that they are elastically mounted in such a way that they can only vibrate on a transverse plane, perpendicular to the bundle. Furthermore, axial effects are not considered and it is assumed that the tubes are of infinite length. The problem is then studied in two dimensions, restricting it to any of the sections of the cavity that are perpendicular to the tubes. With respect to dynamics, it is assumed that the solid-fluid system undergoes small vibrations around a state of equilibrium. 6. Formulation of the Eigenvalue Problem

Let Ro be an open bounded subset of R2,with a locally Lipschitz continuous boundary To (see J. NeEas [19] Chapter I) and let { @ ; } ; = ] , K be a family of K open subsets of Ro, which has the following properties:

V i

=

1,. . . ,K ,

V i

=

1, . . . , K ,

0; is a nonempty connected open subset of Ro . 0;C 0,.

(la) (1b)

v i z j , 0; n Gj = 4 .

(lc)

Each 0; has a locally Lipschitz boundary T i .

(14

Using the above notation we can define R as follows: K

LR

=

no\ U O , . i=l

+

It should be observed that the boundary of R has ( K 1) connected components, which are To,rl ,..., rK. The mathematical model, which describes the solid-fluid interaction, is a differential eigenvalue problem with nonlocal boundary conditions on the velocity. In the model, the and , the domain R represents tube sections are represented by the perforations { O i } i = l , K the area occupied by the fluid. If the velocity of the fluid is denoted by uo = uo(x,r ) and the pressure by PO = po(x, t ) , then (UO,PO) satisfies:

426

CONCA AND DURAN a U0 -

at

2v div e(u0)

+ Vpo = 0

in R ,

div uo = 0 in R ,

(2a) (2b)

uo = 0 on To, (2c) ds; Vi = 1,..., K , V t E R, uo = - on Ti, (24 dt where, in (2a), v represents the kinematic viscosity of the fluid ( v is a strictly positive given constant) and e(uo) is the linear strain tensor, defined by 2e(uo) = Vuo

+ (Vuo)'.

In Eq. (2d), s , is the transverse displacement vector of the ith tube, which, due to the physical assumptions made, only depends on t . If it is furthermore assumed that there is no interaction between the tubes, and given that small oscillations are being considered around a state of equilibrium, the movement of the ith tube obeys a simple harmonic oscillation with a forced term, implied by its interaction with the fluid. Thus s, satisfies the equation: d2sj mi+ klsi = dt2 where m, is the mass per length unit of tube i and k , is a strictly positive real constant, which represents the stiffness constant of the spring system supporting the ith tube (see [9]). The term a(uo,p o ) represents the stress tensor of the system. This satisfies Stokes's law:

d u o , P O ) = -POI + 2ve(u0), (3) where I is the identity matrix. Finally, in (2e), n represents the outward unit normal on the boundary of R . The Eqs. (c, d, e) describe the interactions between the fluid and the tubes. In particular, they model the fact that the fluid, being viscous, adheres to the rigid walls. As is usual in vibration problems, we additionally assume a periodic time dependence. We, therefore, want to find (UO, PO) by: uo(x, t) = u ( x ) e"', (44 po(x,t) = p(x)e"', (4b) where w is the unknown vibratory pulsation of the system. Replacing (4)on the right-hand side of (2e) one can explicitly calculate the unique solution of (2e). We have:

s;(t)

=

k;

+ m;w2

Now, combining Eqs. (4), ( 5 ) with (2a-d), it follows that the triplet ( w , u, p ) ought to be the solution of the following spectral problem on R :

- 2v div e(u) + V p div u = 0 in R , u = O on To,

+ wu

=

0 in

a,

. . . A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

427

To obtain the variational formulation of problem (6), the following Sobolev space is introduced:

H

=

{v E H ’ ( f l ) ’ I div v

=

0 in fl, v

=

0 on To and

v is a constant vector on

Ti,V i

=

1,. . . ,K } .

Clearly H is a closed vector subspace of H ’ ( f 2 ) 2 , and, therefore, a Hilbert space with the induced norm. If one now considers the semi-norm: r

by virtue of Korn’s inequality (see P.-A. Raviart and J.-M. Thomas [20] Chapter 2), it is the case that le(.)lo,n is a norm in H , equivalent to the standard norm induced by H ’ ( f l ) 2 . From now, we will consider H to be equipped with this norm. Multiplying (6a) by V in H , integrating by parts in f2, and using (6c,d) and (3), it follows that if the triplet ( w , u, p ) is a solution to (6), then the pair ( w , u ) is a solution of the following variational eigenvalue problem: Find w E C

u E H , u # 0 such that Vv E H

where, in (7c), yi(u)denotes the trace of u on Ti.Conversely, a standard application of the De Rham Theorem (see J.-L. Lions’ book 1211 Chapter 1 or R. Temam [22] Propositions 1.1 and 1.2) shows that if the pair ( w , u ) is a solution of (7), then there is a function p E L2(f2) such that ( w , u , p ) is a solution of (6). Below, we deal with problem (7) with a view to solving it theoretically and numerically.

C. Theorem of Existence and Location of Eigenvalues

In this section we expound the existence and location theorem of eigenfrequencies of problem (7), which is proved in [17]. The technique used to prove existence consists of showing that the spectrum of (7) coincides with the characteristic values of a quadratic eigenvalue problem involving three linear operators T I ,T2, and Q, which are defined as follows:

-

T2:C2K

T2s = and

Q:H QU =

-

402

H,

Vs E C2K,

C2K,

(YI(u),...,~K(u))V U E H ,

( 104 ( 1Ob)

where functions cpl, 402 are the unique solutions of the following variational problems:

428

CONCA AND DURAN

Find pl E H such that

and Find p2 E H such that K

r

Given that w

=

0 is not a solution of (7), then w

is correctly defined. As is shown in [17] problem (7) is equivalent to the following eigenvalue problem: Find A E C, u E H , u # 0

+ A2T2ddiag(ki)Qu = 0 ,

( T I + T2ddiag(m;)Q)u - Au

(14a) ( 14b)

where ddiag(m,) and ddiag(k,) are diagonal matrices of dimension 2K X 2K, whose diagonals entries are { m l , m i , .. . , m K , m K }and { k l ,k l ,. . . ,k K ,k K } ,respectively. Problem ( 14) is a quadratic eigenvalue problem, which brings into play self-adjoint compact linear operators. The result expounded below summarizes the principal characteristics of the spectrum of (14); more details are available in the work referred to above.

Theorem 1. The spectrum of (14) consists of a countable infinite quantity of complex numbers which converge to Zero: A I , . . . ,A t , . . .

0.

(15)

Moreover, the eigenvalues have the following properties: (i) R e ( & ) < 0 V i 2 1 . (ii) y Im(Ai) # 0, then 1

k/m9

where k / m denotes the quantity max { k j / m j } . ISjSK

Proof. Let us define 1

+ T2ddiag(mi)Q)u. (17) A It is easy to show that if (A, u) is a solution of (14), then v belongs to the finite dimensional space F , defined by v

=

u

F

- --(TI

=

Im(T2) = T2(C2K)C H ,

and (A, ( u ,v ) ) is a solution of the following generalized eigenvalue problem:

.. . A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

429

where IH and I F are the identity operators of spaces H and F, respectively. On the right-hand side of (1 S), i F denotes the canonic imbedding of F into H . Conversely, it is straightforward to see that if ( A , ( u , v ) )E C X ( H X F ) is a solution of (18), then the couple ( A , u ) is a solution of (14) and v verifies (17). Thus, (14) and (1 8) are equivalent. It is appropriate to note that the right-hand side of (18) is a non-self-adjoint linear operator in H X F . Thus, the spectrum of (IS), as well as the spectrum of (14), may admit nonreal eigenvalues. Let us next define A and B by the following rules: A:H X F

and

B:H X F B ( u , v )=

-

-

H X F,

H

X

-

-

F,

[ T2ddiz(k;)Q

Therefore, (IS), [and also (14)] is nothing but the generalized eigenvalue problem associated with the operators A and B , i.e., find A E C and ( u , v ) E H X F , ( u , v ) # (0,O) such that A ( u , v ) = AB(u, v ) .

(21) We can easily prove, using classical theorems of Functional Analysis, that the operator B has a linear continuous inverse. Then, it follows that the problems (14) and (21) are in fact equivalent to the following eigenvalue problem: find A E C and ( u ,v ) E H X F , ( u ,v ) # (0,O) such that T ( u ,v ) = A(u,v ) ,

(22)

where the operator T is defined by T

=

B-'A

Since T I and T2ddiag(mi)Q are compact operators, the operator A is compact. Thus, T is also compact. Applying the spectral theory of compact operators (see T. Kato [23] Chapter HI), it follows that the spectrum of T consists of a countable sequence of complex numbers whose only possible accumulation point is zero. We denote the eigenvalues of T by A I , . . . ,A t , . . . , where they are always assumed to be numbered so that the same eigenvalue is repeated according to its geometric multiplicity. Now, we prove the inequalities. Let ( A , u ) be a solution of (14). Taking the scalar product of Eq. (14b) by u, we obtain [ ( T I + T d d i a g ( m ; ) Q ) uu, ] - Ale(u)lt,, Taking v

=

+ A2[T2ddiag(k,)Qu,u ] = 0 .

u in (7) and using (12), this identity can be rewritten as follows:

430

CONCA AND DURAN

which implies

If Irn(A) = 0, (24a) clearly implies that Re(A) 5 0. On the other hand, if Im(A) # 0, (24b) reduces to K

+ Re(A)EkiIyi(u)I’ = 0 ,

vk(u)l&

i= I

which clearly implies that Re(A) < 0. Combining these results, and using the fact that A = 0 is not eigenvalue of (14), we conclude that Re(A) < 0. Let us next prove (ii). To this end, we assume that Im(A) # 0 and we introduce the vector functions 4 , ,42,which are defined by

41= Tiu, 42 In terms of

41,42,(14b) can

=

T2ddiag(A2ki

+ mi)Qu.

be rewritten as follows:

41 + 42

=

Au.

(26)

Using the definitions of the operators T I , T (cf. ~ (8), ( 1 1) and (9), (12), respectively), it follows that 4 I , 42 verify ~ V [ ~ I , V= ] -

I,

u . Vdx

Vv EH ,

(274

K

2 ~ [ 4 2 ,= ~ ]-

2 (A’k;

+ W Z ; ) ~ ; ( U .) y;(V)

i= I

Multiplying (26) by

we obtain le(4l)l;,il

Taking v

=

+ [42,411= A[u, 413.

u as test function in (27a), we have

and replacing this identity in (28), we can write equivalently

On the other hand, multiplying (26) by

42,

we obtain

Vv EH.

(27b)

. . . A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES But, taking v

=

431

u in (27b), we have

0

and so, replacing [ u , + * ]in (30), we obtain

[417421=

-1442)li.n -

A

2y

Combining (29) with (31), we deduce

hence.‘ we infer that

which implies the following inequality,

and so, we conclude that

k

IN2

>o.

This shows (16) and completes the proof. It is interesting to observe in (ii) that the spectrum of (14) can contain only a finite quantity of eigenvalues whose imaginary part is different from zero. It should be noted that Theorem 1 does not provide an estimate of the number of imaginary (nonreal) eigenvalues; this is provided by the following result, the proof of which can be referred to in [17].

Theorem 2. The spectrum of (14) admits a maximum of 4K imaginary eigenvalues. Moreover, if we denote by 0 < ...

5

p m 5 ... 5 PI

the additive inverses of the real eigenvalues of (14), and by

0<

*.*

the eigenvalues of the operator - ( T I inequalities hold

5J. 5 t j

5

trn 5 ... I

I

+ T 2 d d i a g ( m i ) Q ) ,then

+ 1,

p, It j - 2 K

Vj

2

2K

pj

Vj

=

1 , ..., 2 K .

the following interlacing (32d (32b)

432

CONCA AND DURAN

111. NUMERICAL ANALYSIS OF THE PROBLEM

In this section a numerical solution of problem (6) is undertaken, by means of Lagrange type finite elements on triangles. To this end, certain simplifications are made in the geometry of the problem, which allow a more appropriate numerical approach. First, the discretization of the problem is presented, and then its matrix formulation, which turns out to be a generalized eigenvalue problem of the Sylvester type. Finally, the numerical results that were obtained are presented. A. Discretization of the Problem. The first geometrical hypothesis consists in assuming that R C R2 is a bounded polygonal domain with which we associate a regular family of triangulations { T h } h > 0 (in the sense of Ph. Ciarlet [24], Chapter 2, Section 2.1), such that

where h is defined by h = maxTErhh,; hT being the diameter of triangle T. Another simplification adopted is to assume that all the tubes are identical and that they are supported by the same system of springs. Due to these considerations the constants ki and m iVi = 1 , . . . ,K , have a single value, say k and m, respectively. That is,

k = kl

=

... = kK

m = ml

=

...

=

mK.

(33)

To carry out the numerical analysis of (6), we shall begin by reformulating the problem variationally. While the theoretical study of (6) requires the use of only a primal variational formulation (see (7)), its numerical study needs a mixed variational formulation, in which the incompressibility condition is dealt with implicitly. To this end, we introduce the Sobolev space V , defined by:

V

=

2

{v E H 1 ( R ) I v

=

0 on To, v is a constant vector on Ti,V i

=

1,. . . ,K } . (34)

Clearly, V is a closed vector subspace of H'(R)2(thus, a Hilbert space by itself with the standard induced norm). Multiplying (6a) by V E V and (6b) by 7 E L2(R), it is easy to prove that a mixed variational formulation of (6) is: Find ( w , u , p )E Q= X V X L2(R),* such that

a(u,v)- b(v,p) + wd(u,v)+

(

f m m w 2 ) ~ ( ~ , v )= 0 V v

(35a) E V ,

b(u, q) = 0 v q E L 2 ( R ) , where, in (35b, c), the bilinear forms a(., .), b(., c ( - , d ( . , are defined by: e),

a ( u ,v )

=

b(v,p) =

2v

/I,

I,

a ) ,

e(u) : e ( v ) d x ,

p div V d x ,

0

)

(35b) (35c)

... A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

d(u,v) =

J

n

433

u * Vdx

To approximate V and L 2 ( f l ) ,we introduce the finite dimensional spaces:

o on roand

vh = {vh E c0(IT)* IvhlT E P ~ ( T ) 'v T

E

7 h , v / ~=

is a constant vector on Ti, V i

=

1,. . . ,K},

vh

L'h = { q h E c

" m l q h l T E

P,(T) V T E

Th},

where, in the above definitions, P , ( T )( P k ( T ) k, 2 0 in the general case) denotes the space of polynomials of degree less than or equal to 1 ( 5 k, respectively) and C"( is the space of continuous functions on IT. The positive integers 2Nh and M h are defined as the dimensions of spaces Vh and L i , respectively. Unless it would lead to confusion, we shall omit the subindex h from now on, with the aim of simplifying notation. Let { V I , . . . ,VZN} and { + I , . . . , +,+,} be bases of VIZand L i , respectively. In what follows we shall assume that the family of triangulations { ~ ~ } satisfies h , ~ the following hypothesis of nondegeneracy, concerning the degrees of freedom associated with the velocity and pressure:

a)

2N - M > O .

(36)

In these conditions, we can now approximate problem (35) by: Find

(Wh,Uh,ph)

E

~ i { d ( ~vh) h ,+ mc(uh,vh)}

X

v/,x L;,*

such that V

+ w / ~ ~ ( w , ,+ v ~kC(Uh, ) vh)

b ( u / , , q / ,= ) 0 V qh E L;.

vh

E vh

= ~ h b ( ~p hh ), ,

(374 (37b) (37c)

The rate of convergence estimates can be derived for the approximation of eigenvalues and eigenvectors of this mixed method. We refer to C . Conca, M. Durhn, and J. Rappaz [25]. Briefly, we prove that the problem (30) (discretized problem (37), respectively) can be put in the standard form of a Brezzi-BabuSka saddle-point problem satisfying some adequate hypotheses, which allow the use of several closely related abstract results on eigenvalue approximation. These abstract results can be found in the works of J.H. Bramble and J. E. Osborn [26], I. BabuSka and A. Aziz [27], J. E. Osborn [28], J. Descloux, N. Nassif, and J. Rappaz [29], [30], W.G. Kolata [31], and B. Mercier, J.E. Osborn, J. Rappaz, and P.-A. Raviart [32]. 6. Matrix Formulation of (37)

In this section, problem (37) is reformulated from a matrix point of view. This will allow the approximated solutions ( w h U, h , p h ) to be calculated explicitly. We begin by finding (uh,ph) of the form:

434

CONCA AND DURAN

x x 2N

uh

=

ffrvr

(38a)

9

r=l

M

Ph

=

PI419

(38b)

r=l

where a , , P , E C. By introducing (38a) and (38b) into (37b), we obtain: 2N 6).;

2N

xar{d(vr,Vh) + mc(vr,Vh)} + wh ~ a r a ( v r . v h ) r=l

r=l

2N

M

kxQ'ic(vr.Vh)

f

r=l

This expression is valid for all vh E follows that

vh,

= W / ,x P l b ( V h , 4 , ) . r=1

in particular for every v,, j

=

1,. . . , 2 N . It

Thus, defining the matrices A, B , C, D by: A

=

[ai,] = a ( v i , v j ) i , j

B

=

[b;,] = b(vi,q5,) i

C

=

[c;,] = c ( v ; , v , ) i , j

=

D

=

[dij] = d(vi,vj) i , j

=

=

=

1,. . . , 2 N ,

1 , . .. ,2N,

j

(40a) =

1 , . .., M ,

.,, 2 N , 1,. . . , 2 N ,

1,.

(40b) (40~) (404

and the vectors a,P by: (Y

=

((YI,

. . . ,( Y 2 N ) '

(4 1 a)

p

=

(PI , - . * P d '

(4 1b)

it follows that (39), expressed in matrix terms, is equivalent to

w ~ { D-k mC}a

+ wh A(Y f kCa

=

w,,BP.

(42)

In the same way, introducing (38a) in (37c), we obtain 2N

x a ; b ( v i ,4 j )= 0

vj

=

1,.

. ., M

(43)

i= I

Using the above notation, and conjugating (43), it follows that B* a = 0 , where B* is the transposed conjugate matrix of B. The eqs. (42), (44) reduce problem (37) to: Find (wh,(Y,

P) E C

X

C2NX CM,* such that

(44)

(454

. .. A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

+ mC}a

,i{D

Oh

A(Y 4- k C a =

435

(45b)

WhBP,

B*(Y = 0 ,

(45c)

which is a quadratic eigenvalue problem of finite dimension. Problem (45) can still be transformed if the family of triangulations (Th}h>O satisfies certain conditions of nondegeneracy and regularity. Let us thus observe that if (Y satisfies (49, then (Y

E Ker(B*)

and (Y can be written as a linear combination of elements of a basis of Ker(B*). The following result provides information about this subspace.

Theorem 3. I f the family of triangulations where

1 ~ h denotes l

the cardinality of

Th,

{Th)h>o

satisfies (36) and if

then

dim K e r ( B * ) = 2N - M where 2N and M are the dimensions of

+ 1,

vh and L i , respectively.

Proof. Let x E K e r ( B ) be any element. By definition, the vector x satisfies: Fx,h,

=

0 Vi

=

= ( X I , ..., X M ) '

1, ..., 2 N .

J=I

By introducing (40b) into the above expression, we obtain: Fxj/,

d j div V , d x = 0 V i

=

1,..., 2 N .

(474

j= 1

M

Thus, if we define cp

x j d j , then (47a) is equivalent to:

= j= I

J,

cp div V d x = 0

Vv E

Vh.

Integrating (47b) by parts, we obtain:

Now, as c p l T E P I ( T ) V T E T of T h . Let us denote

Th,

then Vcp is a constant vector of

CT =

Using this notation, and the fact that

VPIT V T E

c2,in every triangle

71,.

re(v) is constant on reV e = 1,. . . ,K , we have:

436

CONCA AND DURAN

Furthermore, (47d) can be rewritten as follows:

where AY = {T E 7 h I T having at least one side on re}. For every T E Ae, we define LT = n T. If it can be proved that a V * E v h exists that fulfills the following properties:

re

we can then conclude the proof. Effectively, taking v CT =

VplT

=

0 VT E

which implies that p is a constant function on

in (47e), if follows that

= v *

Th,

a , which in turn implies that

K e r ( B ) = ((1, ..., 1)'). Thus, range ( B ) = M - 1. However, dim Ker(B*) = dim I m ( B ) L ,

and it can then be concluded that dim K e r ( B n )= 2N - M

+

1.

Let us now prove the existence of such v * E v h . To do this, let us begin by introducing a basis { z I , . . ,Z 2 N } of v h , with the following properties: yp(Z;) = 0

VC

=

1,. . . ,K ,

V i

=

e = 1 , . . . ,K , Ve 1,. . . ,K,

ye(z2;-1) = (1,O)Se; V y e ( z z i )= (0, l)Sei

=

2K

+ 1,. . . , 2 N ,

Vi V i

= =

(494

1,... , K ,

(49b)

. . ,K .

(49c)

1,.

The existence of a basis that fulfills these conditions is evident, if (36) is satisfied. We look for V * in the form 2N

Vf = ~ p ; z ; . i= I

The Eqs. (48a) impose 2 1 7 h l linearly independent restrictions on the coefficients which are:

On the other hand, since the functions { z ; } ~ !satisfy ~ (49), v* will satisfy (48b) if the coefficients {p;}?:~verify the following equations:

. . . A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

437

Thus, on the 2N coefficients p i , we impose 2 1 7 h l + K linearly independent restrictions. As triangulation r h fulfills (46), one can conclude the existence of at least one V * E v h , which satisfies (48). This completes the proof of Theorem 3. Let us denote n = dim Ker(B*),and let ( P I , .. . ,'pn} be a basis of K e r ( B * ) . Let us decompose (Y in this basis, as follows:

where r], E C, V j = 1 , . . . ,n, and where @ is the 2 N by n matrix whose columns are the vectors 'p,. Multiplying (45b) by @*, the following homogeneous equation is obtained: wi@*(D + mC)@r] + w h @*A @r ] + k@*C@r]= 0 .

(504

Let us now fix the following notation:

E

=

@*(D+ mC)@,

F

=

@*A@,

G

=

k@*C@,

and the matrix Eq. (50a) can then be rewritten as follows:

wiEq

+ w,,Fr] + Gr] = 0 ,

(50b)

which is clearly a quadratic eigenvalue problem. The above procedure consists in first solving Eq. (50b), that is, calculating the velocity, and then obtaining the coefficients of the pressure from (45b). That is, first, the following problem is solved: Find (wh, r ] ) E C X d=",

~ i E r+ ] whF7

r]

# 0

such that

+ Gr] = 0 ,

(

5

(5 1b)

and next the coefficients of the pressure are computed by solving the linear system: k BP = W h ( D + mC)@r] 4- A@r] -car]. (52)

+

Oh

Problem (51) is solved by reducing it to a generalized eigenvalue problem. To be more exact, by introducing 6, defined by:

6

=

(F

+

-G)r]. 1 wh

From (51b), it follows that

6

satisfies the equation:

By writing (53) as a matrix system, we obtain:

Thus, if (wh, r ] ) is a solution to (51), then ( O h , r ] , 6) is a solution of (54). Conversely, every solution ( w h , rl.6) is such that 6 fulfills (53a), and accordingly ( w h ,r ] ) is a solution

~

438

CONCA AND

DURAN

of (51). If we use P and Q, respectively, to denote the left-hand and right-hand matrices of (54),we find that (51) is equivalent to: Find ( w h ,7 , l ) E C X C2',(7,&)f 0 such that

(553)

The following theorem summarizes the above and describes the spectrum of,'(37).

Theorem 4. The following statements hold true: (i) Problem (55) admits at most (n + range(G))different solutions. (ii) Problem ( 5 5 ) and system ( 5 2 ) provide all the solutions of the discretized problem (37). Proof. Matrix E is invertible, therefore P is also invertible. Thus, defining R

=

P-'Q,

it is sufficient to note that the degree of the polynomial that provides the nonzero eigenvalues of R is ( n + range(G)), which proves (i). Claim (ii) follows from what was done to deduce (55).

It is interesting to note that the spectrum of the discretized problem (37) might contain nonreal eigenvalues, as seen in problem ( 5 5 ) , in which matrices P and Q are not Hermi tian. C. Numerical Results

The numerical experiments were carried out in a test geometry that consists of a square cavity (transverse section) with a single perforation, also square ( K = l), located in the center of the cavity. This geometry corresponds to the case of having only one tube immersed in the fluid. For the effective numerical solution of (37), the selected bases of spaces Vh and L i are those usually used in the finite elements method. To be more precise, we used finite Lagrange elements on triangles, with the characteristic that the velocity is approximated by continuous functions, which, when restricted to each triangle, are degree 2 polynomials, and that the pressure is approximated by continuous functions, which, on each triangle, are degree 1 polynomials. Thus, associated with each traingle not bearing boundary conditions, there are six degrees of freedom for the velocity (vertices and middle points of the triangles TABLE I. Main features of the triangulations. Feature Number of traingles Total number of vertices and middle points Number of degrees of freedom for the velocity (2N) Number of degrees of freedom for the pressure (M) Number of eigenvalues of (14) (2N - M + 1 + r a n g e ( G ) ) a "See Theorem 2.2.

Th I

Th2

7h3

32

64 160 194 48 149

128 288 450 80

80 98 24 77

373

. .. A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

439

of 7 h ) , and with every triangle of 7 h , three degrees of freedom for the pressure (all vertices of r h )are associated. The boundary condition on TI [see (6d)l was handled explicitly, taking the whole rl as a single (independent) node of the triangulation. The generalized eigenvalue problem of the Sylvester type [see (54)]was solved by using a standard computer software package for numerical calculations. Table I shows Three triangulations were considered: T h , , T h 2 , and Th3 of the domain a summary of the principal characteristics of each of them. It might be observed that all these triangulations satisfy the hypothesis (36) of nondegeneracy, and hypothesis (46) of Theorem 3. Tables 11-IV provide the computed values of the eigenfrequencies for each triangulations of 0. The vicosity of the fluid was fixed on the value Y = 1. In all three tables the calculations were made with five different values of the constant k / m : 0.01, 0.1, 1.O, 10.0, and 100.0. The graphs of some eigenfunctions and their pressure fields (isobaric lines) are included (see Figs. 1 and 2, where we illustrate by means of some examples the different types of symmetries that the solutions can present). As can be seen in Tables 11-IV, the numerical experiments confirm the existing abstract results for this mathematical model. To be more precise, it can be observed

a.

TABLE 11. Characteristic values of problem (14) (Triangulation 7,,,). k/m

1

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 etc.

76 77

w,

0.01 1/A,

k/m

= =

1.967 X 1.967 X -2.686 -3.942 - 3.942 -7.132 -8.937 -9.164 -9.164 - 11.841 -12.173 - 12.173 - 12.596 - 14.285 - 14.637 -16.110 - 16.110 - 16.643 - 16.669 - 17.073 - 17.073 - 17.912 - 18.750 - 20.004 -

etc.

-88.292 -88.292

w,

0.1 1/A,

=

=

1.967 X lo-' -1.967 X lo--' -2.686 -3.941 -3.941 -7.132 -8.937 -9.164 -9.164 - 1 1.841 - 12.173 - 12.173 - 12.596 - 14.285 - 14.637 -16.110 - 16.1 10 - 16.643 - 16.669 - 17.073 - 17.073 - 17.912 - 18.750 -20.004 -

etc.

-88.292 -88.292

k/m

=

1.0

w , = 1/A,

1.975 X -1.975 X -2.686 -3.926 -3.926 -7.132 -8.937 -9.163 -9.163 - 11.841 - 12.173 - 12.173 - 12.596 - 14.285 - 14.637 -16.110 -16.110 - 16.643 - 16.669 - 17.073 - 17.073 - 17.912 - 18.750 - 20.004 -

etc.

-88.292 -88.292

k / m = 10.0 w , = 1/A,

k/m

=

w,

=

100.0 1/A,

-0.206 -2.170 + -2.170 -0.206 -2.170 + -2.686 -2.170 -3.777 -2.689 -3.777 -7.132 -7.132 - 8.937 - 8.937 -9.153 -9.060 -9.153 -9.060 -11.841 - 11.841 -12.172 - 12.155 - 12.172 - 12.156 - 12.596 - 12.596 - 14.285 - 14.285 - 14.637 - 14.637 - 16.1 10 -16.106 - 16.110 - 16.106 - 16.643 - 16.643 - 16.669 - 16.669 - 17.073 - 17.073 - 17.073 - 17.073 -17.912 - 17.912 - 18.750 - 18.750 -20.004 -20.004 etc.

-88.291 -88.291

etc.

-88.283 -88.283

1.798 i 1.798 i 1.798 i 1.798 i

440

CONCA AND DURAN

TABLE 111. Characteristic values of problem (14) (Triangulation T h , ) . k/m

1

w,

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

-1.861 X -1.861 X -2.650 -3.907 -3.907 -8.51 1 - 8.705 -8.901 -8.901 -9.442 - 10.633 - 10.633 - 10.839 - 10.921 - 12.877 - 12.877 - 14.466 - 14.491 - 16.281 - 16.704 - 16.704 - 19.408 -20.378 -20.378

etc.

148 149

0.01 1/A,

k/m

= =

-1.862 X lo-’ -1.862 X lo-’ -2.650 - 3.905 -3.905 -8.511 -8.705 -8.901 -8.901 - 9.442 - 10.633 - 10.633 - 10.839 - 10.921 - 12.877 - 12.877 - 14.466 - 14.491 - 16.281 16.704 - 16.704 - 19.408 -20.378 -20.378 ~

etc.

etc.

198.154 - 198.776

198.154 - 198.776

-

0.1 1/A,

k/m

=

w, =

=

1.0

w , = ]/A,

-1.869 X lo-’ -1.869 X lo-’ -2.650 -3.891 -3.891 -8.51 1 -8.705 -8.900 -8.900 -9.442 - 10.632 - 10.632 - 10.839 - 10.921 - 12.877 - 12.877 - 14.466 - 14.491 - 16.281 - 16.703 - 16.703 - 19.408 - 20.377 -20.377

k/m

198.154 - 198.776

k/m

100.0

=

w , = 1/A,

-0.195 -0.195 -2.650 -3.739 -3.739 -8.511 - 8.705 -8.898 -8.898 -9.442 - 10.627 - 10.627 - 10.839 - 10.921 - 12.874 - 12.874 - 14.466 - 14.491 - 16.281 - 16.701 - 16.701 - 19.408 -20.376 -20.376

-2.085 + -2.085 -2.085 + -2.085 -2.650 -8.51 1 -8.705 -8.879 -8.879 -9.442 - 10.572 - 10.572 - 10.839 - 10.921 - 12.845 - 12.845 - 14.466 - 14.491 - 16.28 1 - 16.675 - 16.675 - 19.408 -20.361 - 20.36 1

etc.

etc.

etc.

-

10.0 1/A,

=

w, =

1.743 i 1.743 i 1.743 i 1.743 i

-198.154 -198.154 - 198.776 - 198.776

-

TABLE IV. Characteristic values of problem (14) (Triangulation q 3 ) . 1

k/m

=

w, =

1 2 3 4 5 6 7 8 9 10 11 12 13 14

0.01 1/A,

-1.868 X -1.868 X -2.625 -3.863 -3.863 -7.909 -7.934 -8.167 -8.167 -9.223 -9.554 -9.844 -9.844 - 10.318

k/m w,

0.1 1/A,

=

=

-1.869 X lo-’ - 1.869 X lo-’

-2.625 -3.862 -3.862 -7.909 -7.934 -8.167 -8.167 -9.223 -9.554 -9.844 -9.844 -10.318

k/m

=

1.0

w , = 1/A,

-1.876 X lo-* -1.876 X lo-‘ -2.625 -3.847 -3.847 -7.909 -7.934 -8.167 -8.167 -9.223 -9.554 -9.844 -9.844 -10.318

k/m

10.0 1/A,

=

w, =

-0.196 -0.196 -2.625 -3.696 - 3.696 -7.909 -7.934 -8.164 -8.164 -9.223 -9.554 -9.841 -9.841 -10.318

k/m

=

w, =

100.0 1/A,

-2.067 + -2.067 -2.067 + -2.067 -2.625 -7.909 -7.934 -8.138 -8.138 -9.223 -9.554 -9.818 -9.818 - 10.318

1.752 i 1.752 i 1.752 i 1.752 i

. . . A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

441

TABLE IV (Continued)

k/m

=

0.01

k/m

=

0.1

k/m

=

1.0

w , = 1/A,

w , = 1/A,

w , = 1/A,

- 11.574 - 11.574

1.574 - 1 1.574 - 1 1.924 -12.173 - 13.724 - 13.724 - 14.096 - 17.256 - 17.758 - 17.758

-11.573 - 11.573 - 11.924 -12.173 - 13.724 - 13.724 - 14.096 - 17.256 - 17.758 - 17.758

k/m

=

10.0

k/m

w , = 1/A,

=

100.0

w , = 1/A, ~~~~~

15 16 17 18 19 20 21 22 23 24 etc.

372 373

-

11.924

- 12.173 - 13.724 - 13.724 - 14.096 - 17.256 -

17.758

- 17.758

etc.

-43 1.5 18 -431.5 18

-1

etc.

-431.5 18 -431.518

etc.

-43 1.518 -431.5 18

Eigenduc

- 11.567

1 1.567 - 11.924 - 12.173 - 13.722 - 13.722 - 14.096 - 17.256 - 17.757 - 17.757 -

etc.

-43 1.5 18 -431.5 18

-11.504 - 1 1.504 - 11.924 -12.173 - 13.701 - 13.701 - 14.096 - 17.256 - 17.748 - 17.748 etc.

-431.518 -431.518

= -2.625

Eigenvalue w = -7 934

FIG. 1. Graphical view of velocity field and isobaric lines (triangulation T , , ~ )case k/m

=

1.O.

442

CONCA AND

DURAN Eigendue a = -8.167

Eigenvalue wlo = -9.223

Eigenvalue wL1 = -9.554

FIG. 2. Graphical view of velocity field and isobaric lines (Triangulation

T ~ , case )

k/m

=

1.0.

that the countable spectrum of the problem converges in modulus to infinity and it admits nonreal eigenvalues, which can be found in the region of the complex plane indicated theoretically (see Theorem 1). It is interesting to note that the quantity of nonreal eigenvalues is no more than 4K (see Theorem 2). In addition, the numerical results provide valuable information on the behavior of the spectrum when k varies (considering that m = 1.0, fixed). In particular, it can be observed that for small values of the parameter k, the spectrum is made up only of real numbers, and that as its value increases, imaginary eigenvalues appear, due to the bifurcation on the complex plane of pairs of real eigenvalues, which collapse together. It is also interesting to note that the eigenvalues of simple multiplicity and large magnitude are insensitive to k's variation, as well as all the graphs of velocity and pressure fields.

... A SPECTRAL PROBLEM IN SOLID-FLUID TYPE STRUCTURES

443

Carlos Conca was partially supported by Fondecyt through grant 1201-91 and D. T. I. under grant E3099-9335. Mario Durhn was partially supported by Fondecyt through grant 004690 and by Commission of the European Communities under The International Scientific Cooperation Scheme B/CI 1*-9232 17.

References 1. R. Blevins, Flow-Induced Vibrations, Van Nostrand Reinhold, New York, 1977. 2. R. Blevins, “Sound induced by vortex shedding from cylinders,” J. Sound Vibration 92, 455 ( 1984). 3. H. J. Connors, “Fluid elastic vibrations of tube arrays excited by cross-flow,’’ in Flow Induced Vibration in Heat Exchangers, R. D. Reiff Ed., A.S.M.E., New York, 1950. 4. S . S . Chen, “Vibrations of nuclear fuel-bundles,’’ Nuclear Engineering Design 35, 399 (1975). 5. S . S . Chen, “Vibration of a row of circular cylinders in a liquid,” J. Engrg. Industry Trans. A.S.M.E. Ser. B 97, 1212 (1975). 6. D. J. Gorman, “Experimental study of peripheral problems to liquid flow-induced vibration in heat exchangers and steam generators,” Paper F 4 6 , Trans. SMIRT, San Francisco, 1977. 7. M. J. Pettigrew, “Flow induced vibration phenomena in nuclear power station components,” Power Industry Research 1, 97 (1981). 8. M. Paidoussis, “A review of flow-induced vibrations in reactors and reactor components,” Nuclear Engineering Design 74, 3 1 (1982). 9. J. Planchard, “Computation of the acoustic eigenfrequencies of cavities containing a tube bundle,” Cornput. Methods Appl. Mechanics Engrg. 24, 125 (1980). 10. J. Planchard, “Eigenfrequencies of a tube-bundle placed in a confined fluid,” Cornput. Methods Appl. Mech. Engrg. 30, 75 (1982). 11. J. Planchard, “Vibrations of nuclear fuel assemblies: a simplified model,” Nuclear Engineering and Design 86, 383 (1985). 12. J. Planchard, “Global behaviour of large elastic tube bundles immersed in a fluid,” Cornputational Mechanics 2, 105 (1987). 13. F. Aguirre and C. Conca, “Eigenfrequencies of a tube bundle immersed in a fluid,” Appl. Math. Optim. 18, 1 (1988). 14. C. Conca and M. Vanninathan, “A spectral problem arising in fluid-solid structures,” Cornput. Methods Appl. Mech. Engrg. 69, 215 (1988). 15. C. Conca, J. Planchard and M. Vanninathan, “Existence and location of eigenvalues for fluid-solid structures,” Cornput. Methods Appl. Mech. Engrg. 77, 253 (1989). 16. C. Conca, J. Planchard and M. Vanninathan, “Limits of the resonance spectrum of tube arrays immersed in a fluid,” J. Fluids Structures 4, 541 (1990). 17. C. Conca, M. Durln and J. Planchard, “A quadratic eigenvalue problem involving Stokes equations,” Cornput. Methods Appl. Mech. Engrg. 100,295 (1992). 18. C. Conca, M. Durhn, and J. Planchard, “A bound for the number of nonreal solutions of a quadratic eigenvalue problem,” Adv. Math. Sci. Appl. 1, 229 (1992). 19. J. NeEas, Les Me‘thodes Directes en The‘orie des Equations Elliptiques, Masson, Paris, 1967. 20. P.-A. Raviart and J.-M. Thomas, Introduction a I’Analyse Nurne‘rique des Equations aux De‘rive‘es Partielles, Masson, Paris, 1983. 21. J.-L. Lions, Quelques Me‘thodes de Re‘solutiondes P r o b l h e s aux Limites Non-Line‘aires, Dunod, Paris, 1969. 22. R. Temam, Navier-Stokes Equations, North-Holland, Amsterdam, 1977. 23. T. Kato, Perturbation Theory for Linear Operators, Springer-Verlag, Berlin, 1976. 24. P. H. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1979. 25. C. Conca, M. Duran and J. Rappaz, “Rate of convergence estimates for the spectral approximation of a generalized eigenvalue problem,” Rapport Interne, Centre de Mathkmatiques AppliquCes, Ecole Polytechnique de Paris, Palaiseau, France (to appear). 26. J. Bramble and J. E. Osbom, “Rate of convergence estimates for non-selfadjoint eigenvalue approximations,” Math. Cornp. 27, 525 (1973).

444

CONCA AND DURAN

27. I. BabuSka and A. Aziz, “Survey lectures on the mathematical foundations of the finite element method,” in The Mathematical Foundations of the Finite Element Method with Applications t o Partial Differential Equations, A. K. Aziz Ed., Academic Press, New York, 1973. 28. J. E. Osborn, “Spectral approximation for compact operators,” Math. Comp. 29, 712 (1975). 29. J. Descloux, N. Nassif and J. Rappaz “Various results on spectral approximation,” Rapport Interne. Department de MathCmatiques, Ecole Polytechnique Fkdkrale de Lausanne, Lausanne, Suisse, 1977. 30. J. Descloux, N. Nassif and J. Rappaz, “On spectral approximation, Part 2: Error estimates for the Galerkin method,” RAIRO Model. Math. Anal. Nume‘r. 7 , 113 (1978). 3 1. W. G. Kolata, “Approximation of variationally posed eigenvalue problems,” Numer. Math. 29, 159 (1978). 32. B. Mercier, J.E. Osborn, J. Rappaz and P.-A. Raviart, “Eigenvalue approximation by mixed and hybrid methods,” Math. Comp. 36, 427 (1981).