- Email: [email protected]
Calculation of medium-frequency vibrations over a wide frequency range
Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191 www.elsevier.com/locate/cma
Calculation of medium-frequency vibrations over a wide frequency range P. Ladeve`ze *, H. Riou LMT Cachan (ENS Cachan/CNRS/Paris 6 University), 61 Avenue du President Wilson, 94235 Cachan Cedex, France Received 17 March 2004; received in revised form 27 August 2004; accepted 30 August 2004
Abstract An enhancement of the ‘‘variational theory of complex rays’’ (VTCR) is being developed in order to calculate the vibrations of slightly damped elastic structures in the medium-frequency range. Here, the emphasis is put on the extension of this theory to the analysis over a wide frequency range. The solution is sought in the form of an original decomposition: the mean value and its complementary part over the frequency range. Numerical examples show the capability of the VTCR to predict the vibrational response of a structure in a frequency range. 2004 Elsevier B.V. All rights reserved. Keywords: Medium-frequency vibrations; Frequency range; Multiscale numerical strategy
1. Introduction The calculation of medium-frequency vibrations is a key element in the design of car chassis or satellites and it represents a real challenge. The difficulty resides in the fact that the length of variation of the phenomena being studied is very small in relation to the characteristic dimension of the structure. Therefore, the use of an approximation based on continuous, piecewise polynomial shape functions leads to models with huge numbers of degrees of freedom (DOFs). Furthermore, pollution effects occur and affect the robustness of the finite element method (FEM) with respect to the wave number k [7,18]. Different approaches have been proposed in order to circumvent this problem. These approaches include predefined reduced bases [32,38], the Galerkin/least-squares FEM [13], the partition of unity method [30], the generalized finite element method [3,39], various multiscale finite element methods [16,27] such as the
*
Corresponding author. Tel.: +33 1 4740 2253; fax: +33 1 4740 2785. E-mail addresses: [email protected] (P. Ladeve`ze), [email protected] (H. Riou).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.08.009
3168
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
residual-free bubbles method [12], the ultra weak variational method [4], the wave envelope method [5], the discontinuous enrichment method [11], the Trefftz methods [8,14,15,31], and the boundary element method [6,33]. However, all these methods are essentially limited to low-frequency problems. Indeed, the use of element-based methods in the medium-frequency range requires that the number of interpolation functions be increased in order to represent the dynamic behaviour correctly. For medium frequencies, these methods are either inaccurate or costly. Statistical energy analysis (SEA, [28]) involves the description of the energy exchanged among various systems and yields the dynamic response of each system averaged over time and space. However, the SEA, which is characterized by a single energy level, cannot provide the spatial variation of the response within each system. Extensions of the SEA have been developed [2,17,25,26], but these extensions still require additional information (e.g. coupling loss factor, energy reflection coefficient, energy transmission coefficient . . .) in order to yield predictive results; therefore, they have been developed only for specific geometries, such as bars or beams. Many improvements still have to be made in order for these methods to become applicable to the medium-frequency range. Moreover, regardless of the methods chosen for discretizing the equations, the resolution of a vibrational problem with multiple frequencies involves the resolution of a series of linear systems of equations with different left- and right-hand sides. The resolution of such systems by a direct method requires the inversion of a number of matrices equal to the number of specified frequencies, which leads to very CPU-intensive computations. Therefore, some techniques have been developed in order to maximize the computational efficiency of the calculations over frequency bands [9,29,36]. These works showed that TaylorÕs, Pade´Õs or WynnÕs approximations of the solution with respect to the frequency seem to be efficient. Nevertheless, these approaches are based on a finite element model of the problem and, therefore, they involve large matrices in the medium-frequency range. The alternative approach discussed here is called the ‘‘variational theory of complex rays’’ (VTCR) [20,22]. It is a predictive tool specifically intended for dealing with these kinds of problems. Its fundamental aspects are described in [21]. The first characteristic feature of the VTCR is that it uses a new variational formulation of the problem being considered which was developed in order to allow the approximations within the substructures to be a priori independent of one another: in other terms, it is not necessary for these approximations to verify a priori the boundary conditions and the continuity conditions at the interfaces between substructures. Instead, these conditions are included in the variational formulation. The second feature which characterizes the VTCR is the introduction of two-scale approximations with a strong mechanical meaning: the solution is assumed to be properly described locally as the superposition of an infinite number of local vibration modes over each substructures, which are supposed to be homogeneous. Each basic mode (which can be an interior ray, an edge ray or a corner ray) associated with a substructure verifies the governing equation and the constitutive law over the substructureÕs domain. All wave directions are taken into account. The unknowns are discretized amplitudes with relatively long wavelengths. The third and last feature characterizing the VTCR is that one retains from the calculated discretized amplitudes only effective quantities related to the elastic energy, the kinetic energy, the dissipation work, the effective displacements, etc. Numerical examples of complex structures consisting of plate and shell assemblies have already demonstrated the capabilities of this formulation at a fixed frequency [34,35]. Recently, Ladeve`ze et al. [24] proposed a new algorithm for the calculation of multiple-frequency solutions using the VTCR. The proposed method uses a set of parameters q to derive a discrete approximation of the rapidly-varying quantities (exponential terms) inside the matrix of the VTCR. These parameters tell how the matrix equation of the VTCR varies in terms of the frequency. Then, three strategies are available in order to obtain the solution: the first is based on averaging over the frequency band; the second and third are based on Newton and
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3169
quasi-Newton algorithms. These strategies do not require the calculation of the solution at each frequency of the frequency band, which is a marked improvement over the standard VTCR. Nevertheless, as pointed out in [24], they require a large number of parameters q to represent the rapidly-varying exponential terms in the matrix properly and, therefore, they lead to large numerical models. In this paper, we present an improvement of the VTCR which enables one to calculate the response in a frequency range using only a small number of parameters to represent the rapidly-varying quantities of the VTCR matrix equations. In the technique presented, the matrix and the right-hand side of the system to be solved are expanded into Taylor series with respect to the frequency. These series are calculated at the center of the frequency range. Then, an ad hoc strategy uses the mean value and the complementary part of the solution over the frequency range to predict the behaviour of effective quantities, such as spatial averages, which leads to the resolution of a system of equations with multiple right-hand sides; this resolution is computationally efficient because the matrix of the system needs to be factored only once. The remainder of this paper is organized as follows: the reference problem to be solved is presented in Section 2; the VTCR as a general approach for solving fixed-frequency vibration problems is reviewed in Section 3; the principle of the approximate formulation, the discretized problem and the numerical strategy are outlined in Section 4; the basic features of the proposed calculation method in a wide frequency range using the VTCR are described in Section 5; in Section 6, numerical illustrations are given and a comparison with the strategy proposed in [36] is made using the example of a simply-supported plate; finally, the conclusion is presented in Section 7.
2. The reference problem In order to simplify the presentation, the problem will be formulated for an assembly of two substructures, but it is easy to generalize to an assembly of n substructures. Given two substructures S and S 0 , let oS and oS 0 represent the boundaries of S and S 0 respectively. We are studying the harmonic vibration of these two structures at a fixed frequency x. All quantities can be defined in the complex domain. An amplitude Q(X) is associated with Q(X)eixt. The excitations applied to S and represented in Fig. 1 are 1. a displacement field Ud on a portion oUS of the boundary oS, 2. a density of forces Fd on a portion oFS, 3. a density of forces f d over the whole domain S.
Fig. 1. The reference problem.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3170
Similar quantities are defined on S 0 . For a structure S, let Sad denote the field of the pairs (U, r) such that: 3
U 2U
ðfinite-energy displacement set ½H 1 ðSÞ Þ;
r2S
ðfinite-energy stress set ½L2 ðSÞ Þ;
divr þ f d ¼ x2 qð1 ihÞU r ¼ ð1 þ igÞKðU Þ
on S; on S;
3
ð1Þ
where K is HookeÕs operator, q the mass density, and h and g the damping coefficients (which depend on the frequency). The subspace Sad associated with homogenized conditions (f d = 0) is designated by Sad,0. In the same manner, we introduce the spaces S0ad and S0ad;0 . The reference problem can be formulated as follows: find (U(X), r(X), X 2 S) and (U 0 (X), r 0 (X), X 2 S 0 ) such that s ¼ ðU ; rÞ 2 Sad ; U ¼ U d on oU S; rn ¼ F d U ¼ U0
on oF S;
s0 ¼ ðU 0 ; r0 Þ 2 S0ad ; U 0 ¼ U 0d on oU S 0 ; 0 0
rn ¼
F 0d
ð2Þ
0
on oF S ;
on C;
rn þ r0 n0 ¼ 0 on C; where n (resp. n 0 ) is the outward normal vector of the boundary of S (resp. S 0 ).
3. Variational formulation associated with the VTCR This section reviews the main characteristics of the VTCR. For more details, the reader is referred to [20]. The VTCR is a global formulation of the boundary conditions (2) in terms of both displacements and forces. It is based on a priori independent approximations within the substructures: find s = (U, r) 2 Sad and s0 ¼ ðU 0 ; r0 Þ 2 S0ad such that: Z Z Z drnðU U d Þ dl þ ðrn F d ÞdU dl þ dr0 n0 ðU 0 U 0
Re ix d Þ dl 0 oU S oF S oU S Z Z 1 0 0 0 0
0 0
0
0 0
0
þ ðr n F d ÞdU dl þ ðdrn dr n ÞðU U Þ þ ðrn þ r n ÞðdU þ dU Þ dl ¼ 0 2 C oF S 0 8ðdU ; drÞ 2 Sad;0
8ðdU 0 ; dr0 Þ 2 S0ad;0 ;
ð3Þ
where Re[A] denotes the real part of quantity A and A* the conjugate of A. It is easy [20] to prove that the variational form is equivalent to the reference problem provided that: 1. the reference problem (2) has a solution, 2. HookeÕs operators K and K 0 are positive definite, 3. one of the damping factors h, g (resp. h 0 , g 0 ) is greater than zero on S (resp. S 0 ). Contrary to the usual formulations, displacements and forces may be applied simultaneously on the same boundary oLS. This enables one, for example, to take into account all the data of the problem available from an experiment in which both the forces and the displacements may be known at the same boundary. In this case, the following term should be considered:
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
1 Re ix 2
Z
drnðU U d Þ þ ðrn F d ÞdU dl
3171
:
oL S
It is possible [20] to express the formulation as: find s = (U, r) 2 Sad and s0 ¼ ðU 0 ; r0 Þ 2 S0ad such that: s s s 0 0 xdðED ðU Þ þ ED ðU ÞÞ þ 8ds0 2 S0ad;0 ; ;d 0 ¼ LD ; d 0 8ds 2 Sad;0 0 s s s where ED is the dissipated energy, LD a linear form, and h. , .i a bilinear form defined at the boundaries which verifies hu, vi = hv, ui*.
4. Approximate formulations 4.1. Principle: the general case All that one needs to do in order to derive approximate formulations from the VTCR is to define the 0h 0 0 subspaces Shad and Shad;0 (resp. S0h ad and Sad;0 ) from Sad and Sad,0 (resp. Sad and Sad;0 ). The approximate forh 0h h 0h h h 0h 0h mulation can be expressed as: find s ¼ ðU ; r Þ 2 Sad and s ¼ ðU ; r Þ 2 Sad such that: h h h s s s xdðED ðU h Þ þ E0D ðU 0h ÞÞ þ ; d 8ds0h 2 S0h ; d ¼ L 8dsh 2 Shad;0 D ad;0 : s0h s0h s0h The VTCR uses two approximation scales (Uh, rh), which have a strong mechanical meaning, by distinguishing three zones: the interior zone, the edge zone and the corner zone. For example, in the neighborhood of a point X of the interior zone, the solution is assumed to be properly described locally as the superposition of an infinite number of complex rays (propagative and evanescent waves) which can be written as follows: U h ðX ; Y ; P Þ ¼ W h ðX ; P ÞeixP Y ; rh ðX ; Y ; P Þ ¼ Ch ðX ; P ÞeixP Y ; where both X and Y represent the position vector, X being associated with slow variations and Y with rapid variations. More precisely, the terms related to the position vector X vary slowly when X moves along the structure, whereas the terms related to the position vector Y vary rapidly when Y moves along the structure. P is a vector which characterizes the complex ray. In order for these local modes (Uh, rh) to be admissible, they must lie in Shad and verify (1). Thus, we obtain some properties on P. These properties enable one to classify the rays into interior rays (propagative waves), and edge and corner rays (evanescent waves). 4.2. Admissible complex rays for a plate Let us consider the out-of-plane bending motion of a thin, flat, homogeneous and isotropic plate. According to KirchhoffÕs thin plate theory, the steady-state displacement u of the plateÕs mid-surface perpendicularly to the plate is governed by the dynamic equation: Eh3 ð1 þ igÞMMu ¼ qhx2 ð1 ihÞu on S; 12ð1 m2 Þ where n is the Laplacian operator, E the YoungÕs modulus, h the thickness of the plate, m the PoissonÕs ratio, q the mass density, x the frequency, and h and g the damping factors. A complex ray for the interior zone is
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3172
uhi ðX ; Y ; P Þ
¼
whi ðX ; P Þ exp
pffiffiffi h þ g pffiffiffiffi xP X ei xP Y : 4
This complex ray corresponds to a plane bending wave which propagates through the plate in direction P. This ray is admissible (i.e. it belongs to Sad) only if: Eh3 ð1 þ igÞMMuhi ¼ qhx2 ð1 ihÞuhi 12ð1 m2 Þ
on S:
Therefore, the properties on P are 12qð1 m2 Þ 2 ðP P Þ ¼ r4 with r4 ¼ : h2 E
ð4Þ
Eq. (4) shows that P lies on a circle Ci with a radius r defined by the material properties. By following this circular path, one takes into account all the directions of the plate. Examples of such interior rays are shown in Fig. 2. A similar approach can be used for the edge and corner zones. The properties on P for an edge complex ray of the form pffiffiffi pffiffiffi h þ g pffiffiffiffi r x n X e xP n nY þi xP t tY ; uhe ðX ; Y ; P Þ ¼ whe ðX ; P Þ exp 4 Pn which verifies Eh3 ð1 þ igÞMMuhe ¼ qhx2 ð1 ihÞuhe 12ð1 m2 Þ
on S
and does not propagate faster than the interior rays can be written as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P t ¼ r sinðuÞ P n ¼ r 1 þ sin2 ðuÞ;
ð5Þ
(where n and t denote respectively the normal and tangent vectors to the edge). Eq. (5) shows that P lies on a curve Cb, defined by the material properties, spanned by u. By following this curve, one takes into account all edge complex rays. Examples of such modes are shown in Fig. 3 (left and center). The properties on P for a corner complex ray of the form pffiffiffi h þ g pffiffiffiffi h h xP X e xP Y ; uc ðX ; Y ; P Þ ¼ wc ðX ; P Þ exp i 4 which verifies Eh3 ð1 þ igÞMMuhc ¼ qhx2 ð1 ihÞuhc 12ð1 m2 Þ
on S
Fig. 2. Admissible interior rays for a plate.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3173
Fig. 3. Admissible edge modes (left and center) and corner modes (right) for a plate.
can be written as 2
ðP P Þ ¼ r4
with r4 ¼
12qð1 m2 Þ : h2 E
ð6Þ
Eq. (6) shows that P lies on the same circle, denoted Cc, as in the case of interior complex rays. Examples of such modes are shown in Fig. 3 (right). One must note that the admissible complex rays are exact solutions of the governing equation of the problem. Therefore, unlike conventional finite element-based methods, the VTCR yields no numerical dispersion [7,18]. 4.3. The discretized problem The displacement of each point in the substructure is generated by a basis of admissible complex rays. The unknowns are the generalized amplitudes whi ðX ; P Þ; whb ðX ; P Þ and whc ðX ; P Þ of the basis (nth-order polynomials in X and large-wavelength quantities). Taking into account all the directions u in Ci, Cb and Cc, one ends up with an integral which takes the form: Z pffiffiffi hþgpffiffiffi whi ðX ; P ðuÞÞe 4 xP ðuÞX ei xP ðuÞY dC uh ðX ; Y Þ ¼ u2C i
þ þ
Z
Z
whb ðX ; P ðuÞÞexp
ffiffiffi
hþgp r xP n nX 4
e
pffiffiffi pffiffiffi xP n nY þi xP t tY
dC
u2C b
whc ðX ; P ðuÞÞei
ffiffiffi
hþgp xP ðuÞX 4
e
pffiffiffi xP ðuÞY
dC:
ð7Þ
u2C c
In order to obtain a finite-dimension problem, one can discretize this integral (7) and assume that the amplitudes whi ðX ; P ðuÞÞ; whb ðX ; P ðuÞÞ and whc ðX ; P ðuÞÞ are constant over angular sectors: whi ðX ; P ðEP ÞÞ; whb ðX ; P ðEP ÞÞ and whc ðX ; P ðEP ÞÞ (for example, see Fig. 4 for the interior rays). The angular distribution of the plane waves for any point in the substructure is assumed to be properly described by such a discontinuous angular distribution. Ci, Cb and Cc is discretized into N elements, which can be of different sizes (on Fig. 4, Ci is discretized into eight elements of the same size). Once the discretization has been chosen for each plate, the VTCR leads to a system of linear equations in the complex domain: K U ¼ F;
ð8Þ
where K ¼ xEhD þ Z h and F ¼ LhD . EhD is the symmetric, positive definite damping matrix associated with ED. Zh is the matrix associated with the bilinear form h. , .i (see Section 3), defined such that Z Th ¼ Z h . LhD is the vector associated with the linear form LD. U is the vector corresponding to the unknown amplitudes associated with the complex polynomial wh(X, EP). As a consequence of the last properties, Eq. (8) has a unique solution.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3174
Fig. 4. Discretized amplitudes.
4.4. The numerical strategy The system of linear equations leads to the solution of (8). K corresponds to the bilinear terms in the h 0h variational formulation (3) projected onto ðShad [ S0h ad Þ ðSad;0 [ Sad;0 Þ, U is the vector of the unknown genh h h eralized amplitudes wi ðX ; P ðEP ÞÞ, wb ðX ; P ðEP ÞÞ and wc ðX ; P ðEP ÞÞ, and F corresponds to the projection onto ðShad [ S0h ad Þ of the load terms of the variational formulation (3). A typical example of a term of K would be the coupling of an interior sector EP with an edge sector E0P on a plate S, as shown in Fig. 5. The coupling of these two angular sectors leads to the following integral: Z Z Z I¼ aðP ; P 0 ; X Þ exp½ixðP P 0 Þ Y du0 du dS; ð9Þ S
EP
E0P
Fig. 5. Typical term of the symmetric part of the coupling matrix.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3175
where a is a polynomial, and P and P 0 are the wave vectors belonging to the two angular sectors EP and E0P . The structure S is subdivided into elements for integration purposes. The size of these elements is chosen to be approximately one wavelength, so that slowly varying quantities can be considered to be constant. Fig. 6 shows the spatial decomposition used for the integration. Then, Eq. (9) becomes: Z Z Z X I¼ aðP ; P 0 ; X Þ exp½ixðP P 0 Þ Y du0 du dD: ð10Þ subdomains D
D
EP
E0P
Assuming the polynomial a to be constant over each subdomain D and denoting XE the center of the cell D, one can rewrite Eq. (10) as Z Z Z X I¼ aðP ; P 0 ; X E Þ exp½ixðP P 0 Þ Y du0 du dD: ð11Þ subdomains D
EP
E0P
D
With a change of variable transforming D into a symmetric domain in N1 and N2, Eq. (11) becomes: Z Z X I¼ aðP ; P 0 ; X E Þ exp½ixðP P 0 Þ X E du0 duJ ð12Þ subdomains D
EP
E0P
with J¼
Z
h1 2 h
21
exp½ixðP 1
¼ h1 Sinc xðP 1
P 01 ÞX 1 dX 1
h1 P 01 Þ 2
Z
h2 2 h
22
exp½ixðP 2 P 02 ÞX 2 dX 2
h1 Sinc xðP 2
h2 P 02 Þ 2
;
where Sinc ½x ¼
sin½x : x
Now, one must integrate over the angular domains EP and E0P . Using the same technique, one decomposes these domains into subdomains within which the polynomial can be considered to be constant. Fig. 7 shows this angular decomposition.
Fig. 6. Spatial decomposition for the integration.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3176
Fig. 7. Angular decomposition for the integration.
Denoting ehP and eh0P the angular subdomains, I becomes: I¼
X subdomains D
with
2
X X Z
4
ehP 2EP eh0 2E0 P P
ehP
3
Z 0
ehP
H du du0 5
ð13Þ
0 h1 0 h2 H ¼ h1 h2 aðP ; P ; X E Þ exp½ixðP P Þ X E Sinc xðP 1 P 1 Þ Sinc xðP 2 P 2 Þ : 2 2 0
0
Once again, one calculates the quantities subjected to rapid variations. These quantities are integrated analytically, whereas the others are calculated numerically. Because of the domains of variation of the two functions Sinc, these can be treated as slow functions and, by linearizing the angular dependency of the wave vectors P and P 0 as shown in Fig. 8.
Fig. 8. Linearization of the wave vector P.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
One can write Eq. (13) as Z Z Z 0 0 G¼ exp½ixðP P Þ X E du du ¼ ehP
0
ehP
¼ exp½ixðP h P 0h Þ X E
s1
exp½ixP X E ds
s1
Z
s1
exp½ixsT h X E ds
s1
Z
s01 s01
Z
s01 s01
3177
exp½ ixP 0 X E ds0
exp½ ixsT 0h X E ds0
¼ 4s1 s01 exp½ixðP h P 0h Þ X E Sinc ½xs1 T h X E Sinc ½xs01 T 0h X E : Since the angular dependencies of the edge and corner modes are constructed in the same way as those of the interior modes (see Section 4.1), their linearization yields exactly the same results as in Fig. 8. Finally, Eq. (13) becomes: X X X I¼ 4h1 h2 s1 s01 aðP h ; P 0h ; X E Þ subdomains D eh 2EP eh0 2E0 P P P
h1 h2 Sinc xðP h1 P 0h1 Þ Sinc xðP h2 P 0h2 Þ Sinc ½xs1 T h X E 2 2 Sinc ½xs01 T 0h X E exp½ixðP h P 0h Þ X E ;
which, because of the slow evolution of the functions Sinc in their respective domains, can be written as X X X ½AðP h ; P 0h ; X E Þ exp½ixðP h P 0h Þ X E ; ð14Þ I¼ subdomains D eh 2EP eh0 2E0 P P P
where A represents a function which contains all the slow-variation terms. All the rapid-variation terms are included in the exponential functions. This particular decomposition will be used to calculate the solution within a frequency range. 4.5. Effective quantities From a mechanical point of view, the spatial distribution of the solution in the medium-frequency range has no physical meaning. Indeed, the response at a specific location and at a specific frequency is extremely sensitive to perturbations. Spatially averaged quantities are more useful for interpreting the response. Therefore, it is necessary to extract from the solution some effective quantities defined over a domain D greater than one wavelength. These quantities extracted from the local data A are called effective quantities and are defined as Z 1 Aeff;D ¼ jAðX Þj dX : D D From here on, we will use the effective displacement over a domain D: Z 1 U eff;D ¼ jU ðX Þj dX : D D 5. Calculation in a frequency range The objective of this section is to present a variation of the VTCR capable of calculating the solution in a frequency range B defined by its central frequency x0 and its bandwidth Dx: find U which verifies KðxÞUðxÞ ¼ FðxÞ 8x 2 B ¼ ½x0 Dx; x0 þ Dx:
ð15Þ
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3178
This approach, first introduced in [19], uses the decomposition presented in Eq. (14). Considering the nature of the frequency range, let us first assume that the function A is constant on B. Consequently, one has AðP h ; P 0h ; X E ; xÞ ’ AðP h ; P 0h ; X E ; x0 Þ
for x 2 B:
If necessary, one could build a better approximation (for example, a linear approximation): AðxÞ ’ Aðx0 Þ þ
Aðxþ Þ Aðx Þ ðx x0 Þ xþ x
for x 2 B;
where x+ and x represent the limits of the frequency range. Let us now study the behaviour of the exponential argument. This quantity can be expanded into a Taylor series up to order k: " # k X ðiðP h P 0h Þ X E Þk k 0 0 ðx x0 Þ : exp½ixðP h P h Þ X E ¼ exp½ix0 ðP h P h Þ X E k! k¼0 Then, according to Section 4, the matrix K and the load F can be approximated as KðxÞ ¼
k X
Kk ðx0 Þðx x0 Þk ;
k¼0
FðxÞ ¼
k X
ð16Þ k
Fk ðx0 Þðx x0 Þ :
k¼0
Let us define the mean value over the frequency range B: Z 1 dx: hi ¼ 2Dx B The matrix K and the load F can be expressed as KðxÞ ¼ hKi þ DKðxÞ;
ð17Þ
FðxÞ ¼ hFi þ DFðxÞ: According to the Taylor expansion (16), one has DKðxÞ ¼
k X
k
DKk ðx0 Þðx x0 Þ ;
k¼0
DFðxÞ ¼
k X
k
DFk ðx0 Þðx x0 Þ :
k¼0
One can also define U using the same approach: UðxÞ ¼ hUi þ DUðxÞ:
ð18Þ
By combining (15), (17) and (18), one gets ½hKi þ DK ½hUi þ DU ¼ ½hFi þ DF:
ð19Þ
Eq. (19) involves terms with very different magnitudes: some are very large while others are very small. Using the techniques of perturbation methods, let us now identify terms of different orders and rewrite Eq. (19) for order 0 and order 1:
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
hKihUi ¼ hFi
order 0;
hKiDU ¼ DF DKhUi
order 1:
3179
ð20Þ
Then, DU is equal to DU ¼
k X
k
DUk ðx x0 Þ :
ð21Þ
k¼0
Finally, the displacement over the frequency range B which we retain is UðxÞ ¼ khUi þ lDU: where k and l minimize the error defined through the approximations in (20): Z 1 T E2 ¼ ½F K U ½KS D ½F K U dx; 2Dx B where [KS]D is the diagonal of the symmetric part of K, which is a slowly-varying quantity [20]. Remarks: (1) One can write U = hKi 1L(x x0), where L denotes an operator. The strategy we described can be viewed as the resolution of a finite-dimension linear problem with multiple right-hand sides. (2) One must note that the only parameter which must be chosen in order to ensure good accuracy is k, which corresponds to the degree of the Taylor expansion. (3) The mean dissipated energy over the frequency range is Z k2 l2 e ¼ hUihKS ihUi þ DUKS DU dx: 2 4Dx B (4) In [9], Pade´Õs or WynnÕs approximations seem to improve the interval of convergence of the Taylor series. In this paper, Taylor series were chosen for four reasons: (a) in [29] and [9], the whole vibrational solution was discretized, whereas in the VTCR only the slowly-varying part, i.e., the amplitude of the waves, a quantity less sensitive to a variation of frequency, is discretized (in the examples in the next section, only five terms of the Taylor series suffice to accurately represent the amplitude of the waves across a range of frequencies), (b) as pointed out in [9], if the number of approximation terms is large, Pade´Õs approximation leads to an ill-conditioned matrix and WynnÕs algorithm becomes computationally expensive, (c) Taylor series are easier to implement, and (d) the proposed strategy is not limited to a Taylor series expansion, as k and l can be x-varying functions.
6. Numerical examples 6.1. A simple example to test the method over the frequency range To illustrate the approach discussed above, let us consider the dynamic response of a beam subjected to an external excitation (see Fig. 9). The purpose of this simple example is to test the method without introducing any approximation, since for one-dimensional structures the VTCR yields the exact solution at a fixed frequency. We used the method to find the reference solution over the frequency range [800; 4000 Hz]. The reference solution is the analytical solution. The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz.
3180
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
F
Fig. 9. Model geometry for a beam in bending; E = 200 GPa; m = 0.3; g = 0.05; h = 0; q = 7800 kg/m3; S = 0.0001 m2; I = 8.3 · 10 10 m4; L = 2.5 m; F = 1 N/m.
In Figs. 10–12, the energy density at the center of the beam is plotted against the frequency for different orders of Taylor series expansion (see Section 5, Eq. (16)). One can see that the solution obtained with the frequency range calculation is very good, even on the portions of the curves with highly dynamic behaviour. As the order of the Taylor series expansion increases, the precision improves. At order 4, the approximate solution matches the reference solution very closely. This can also be observed in Fig. 13, which shows the error as a function of the expansion order: as expected, as the order increases, the error decreases. Figs. 14–16 show the solution of the problem at 2610 Hz (which is not the central frequency of a frequency range) for different orders of Taylor series expansion. These results show that the accuracy increases as the expansion order increases. Again, one can see that order 4 provides sufficient accuracy. In the following examples, all the approximate frequency range solutions were calculated to order 4. 6.2. Plate example This second example shows the effectiveness of the VTCR in the case of a plate structure whose geometry is shown in Fig. 17. The out-of-plane bending motion is induced by an external force applied in the middle
0. 2 0.18
Computation order 0 Exact theory
0.16
0.12 0. 1
u
E (L/2) (J/m)
0.14
0.08 0.06 0.04 0.02 0 500
1000
1500
2000
2500
3000
3500
4000
Frequency ω (Hz)
Fig. 10. Result for the beam problem (see Fig. 9). The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz and a Taylor series expansion up to order 0 (see Section 5, Eq. (16)).
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3181
0.2 0.18
Computation order 2 Exact theory
0.16
u
E (L/2) (J/m)
0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 500
1000
1500
2000
2500
3000
3500
4000
Frequency ω (Hz)
Fig. 11. Result for the beam problem (see Fig. 9). The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz and a Taylor series expansion up to order 2 (see Section 5, Eq. (16)).
0.2 0.18
Computation order 4 Exact theory
0.16
0.12 0.1
u
E (L/2) (J/m)
0.14
0.08 0.06 0.04 0.02 0 500
1000
1500
2000
2500
3000
3500
4000
Frequency ω (Hz)
Fig. 12. Result for the beam problem (see Fig. 9). The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz and a Taylor series expansion up to order 4 (see Section 5, Eq. (16)).
of the plate, perpendicular to the plateÕs mid-surface. The exact analytical solution for this problem is given in [10].
3182
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191 1
0.9
Energy relative error (x=L/2)
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0 1
1.5
2
2.5
3
3.5
4
4.5
5
k
Fig. 13. Relative energy error for the beam problem (see Fig. 9).
0.04
0.035
Computation order 0 Exact theory
0.025
0.02
u
E (L/2) (J/m)
0.03
0.015
0.01
0.005
0 0
0.5
1
1.5
2
2.5
x (m)
Fig. 14. Result for the beam problem (see Fig. 9) at x = 2p · 2610 Hz. The approximate solution was calculated with a Taylor series expansion up to order 0 (see Section 5, Eq. (16)).
This example was chosen so that we could compare our approach to that developed by Soize (1982) in [36], based on a low-frequency equation associated with the medium-frequency problem (LFEA). This method was used, in particular, in [37]. In [38], a numerical method for the analysis of medium-frequency
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3183
0.04
0.035
0.03
Computation order 2 Exact theory
E (L/2) (J/m)
0.025
u
0.02
0.015
0.01
0.005
0 0
0.5
1
1.5
2
2.5
x (m)
Fig. 15. Result for the beam problem (see Fig. 9) at x = 2p · 2610 Hz. The approximate solution was calculated with a Taylor series expansion up to order 2 (see Section 5, Eq. (16)).
0.04
0.035
Computation order 4 Exact theory
0.025
0.02
u
E (L/2) (J/m)
0.03
0.015
0.01
0.005
0 0
0.5
1
1.5
2
2.5
x (m)
Fig. 16. Result for the beam problem (see Fig. 9) at x = 2p · 2610 Hz. The approximate solution was calculated with a Taylor series expansion up to order 4 (see Section 5, Eq. (16)).
linear vibrations of elastic structures is proposed. This technique is used to characterized the vibratory state of a structure in frequency ranges. The finite element approximation of the solution sought in the frequency domain for an n-dimensional subspace is written as
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3184
F
Fig. 17. Model geometry for the plate; E = 210 GPa; m = 0.3; g = 0.001; h = 0; q = 7800 kg/m3; h = 0.005 m; dimensions 1 m · 1 m; F = 1 N.
Uðx; xÞ ¼
n X
U i ðxÞN i ðxÞ;
i¼1
where Ni denotes the finite element shape functions and Ui the ith nodal response quantity. The governing equation of motion is expressed as ð x2 M þ ixD þ KÞU ¼ F
8x 2 B ¼ ½x0 Dx; x0 þ Dx;
ð22Þ
where M, D, K and F are the mass, damping and stiffness matrices, and the force vector respectively. Since this medium-frequency equation (22) cannot be solved numerically using a time integration method (because the high frequencies involved would lead to too small a time step), the method consists in constructing a low-frequency equation using the new variable U0 ðxÞ ¼ Uðx þ x0 Þ;
3
x 10
–6
Exact solution
2.5
Effective displacement
VTCR 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 18. Result for the simply-supported plate problem (see Fig. 17) using the VTCR; Dx = 2p · 100 Hz and x0 2 2p · {500; 700; 900} Hz.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3185
–6
3
x 10
Exact solution
2.5
Effective displacement
LFEA 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 19. Result for the simply-supported plate problem (see Fig. 17) using the LFEA; Dx = 2p · 100 Hz and x0 2 2p · {500; 700; 900} Hz.
3
x 10
–6
Exact solution
2.5
Effective displacement
VTCR 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 20. Result for the simply-supported plate problem (see Fig. 17) using the VTCR; Dx = 2p · 300 Hz and x0 = 2p · 700 Hz.
3186
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3
x 10
–6
Exact solution
2.5
Effective displacement
LFEA 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 21. Result for the simply-supported plate problem (see Fig. 17) using the LFEA; Dx = 2p · 300 Hz and x0 = 2p · 700 Hz.
where U0 verifies the equations ð x2 M þ ixDn þ Kn ÞU0 ¼ F0 Dn ¼ D þ 2ix0 M; Kn ¼ K þ ix0 D x20 M;
8x 2 B0 ¼ ½ Dx; Dx; ð23Þ
F0 ðxÞ ¼ Fðx0 þ xÞ: The low-frequency Eq. (23) can be resolved using step-by-step time integration, which yields the desired solution of Eq. (22): X UðxÞ ¼ sL U0 ðmsL Þe imsL x 8x 2 B; m2Z ð24Þ p : sl ¼ Dx The parameters which must be chosen to ensure good accuracy are mI and mS, which correspond to the initial and final time integration in (24). To obtain the following results, we chose mI = 10 and mS = 40. This choice yielded good enough results with a reasonable calculation time. If we had taken larger parameters, the calculation time would have become really excessive. For the simply-supported plate example defined in Fig. 17, the two strategies (VTCR and LFEA) were applied in two cases. For the first case, one has Dx = 2p · 100 Hz and x0 2 2p · {500; 700; 900} Hz (Figs. 18 and 19). For the second case, one has Dx = 2p · 300 Hz and x0 = 2p · 700 Hz (Figs. 20 and 21). For the VTCR, we used 20 interior rays, 20 edge rays and no corner ray, i.e., 40 of freedom (DOFs). The Taylor series expansion (see Section 5) was taken up to order 4. For the LFEA, we calculated the solution according to the classical rule which says that one must use between 8 and 10 elements per wavelength (see [1]), which leads to 5766 DOFs.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3187
1
F Fig. 22. Model geometry for a 3D plate assembly; E = 200 GPa; m = 0.3; g = 0.001; h = 0; q = 7800 kg/m3; h = 0.004 m; F = 1 N/m.
In Figs. 18 and 19, one can see that both approaches match the exact analytical solution with good enough accuracy provided the bandwidth Dx is not too large (here, Dx = 2p · 100 Hz). In Figs. 20 and 21, one can see that neither approach leads to the exact analytical solution easily if the bandwidth becomes too large (here, Dx = 2p · 300 Hz). In the case of the VTCR, some peaks are shifted (see Fig. 20 at 900 Hz), whereas with the LFEA the levels of the peaks are not good enough (see Fig. 21). On these examples, we can observe that both strategies operate on the same frequency range. By comparing the numbers of DOFs required (40 DOFs for the VTCR needs, 5766 DOFs for the LFEA), one notes that the VTCR is less memory-consuming. 6.3. 3D plate assembly This third example shows the effectiveness of the VTCR in the case of a 3D structure intended to represent a car, whose geometry is shown in Fig. 22. The calculation using the VTCR involved 360 DOFs. The displacement obtained at x = 2p · 100 Hz is shown in Fig. 23 and compared with results from NASTRAN. The NASTRAN solution used QUAD4 elements and a mesh size of 10 elements per wavelength, yielding 33,230 DOFs. For this example, on a single-processor PC Linux system (INTEL 1.1 Ghz, 512 Mb RAM), it took 13 s to compute and inverse the matrix with the VTCR vs. 1 min and 13 s with NASTRAN. The two results are similar. The evolution of the effective displacement of Plate 1 of this assembly (the roof of the car; see Fig. 22) with the frequency of the excitation is shown in Figs. 24–26. The reference curve (dotted line) was calculated frequency by frequency. The approximate solution was calculated using the approach presented above (see Section 5) with a Taylor series expansion up to order 4. Again, these good results confirm the capability of the proposed method to find the solution over a frequency range. Fig. 27 shows the same results as in Figs. 24–26, this time collected next to one another. Thus, the frequency range extends from x = 2p · 50 Hz to
Fig. 23. Comparison of the displacements of the solutions for the 3D plate assembly (see Fig. 22) at x = 2p · 300 Hz: VTCR (left) vs. NASTRAN (right).
3188
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191 x 10
2
–5
1.8
Reference solution Frequency band analysis
1.6
Effective displacement
1.4 1.2 1 0.8 0.6 0.4 0.2 0 50
60
70
80
90
100
110
120
130
140
150
Frequency (Hz)
Fig. 24. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 = 2p · 100 Hz, Dx = 2p · 50 Hz.
–5
2
x 10
1.8
Reference solution Frequency band analysis
1.6
Effective displacement
1.4 1.2 1
0.8 0.6 0.4 0.2 0 150
160
170
180
190
200
210
220
230
240
250
Frequency (Hz)
Fig. 25. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 = 2p · 200 Hz, Dx = 2p · 50 Hz.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3189
–6
3
x 10
2.5
Reference solution
Effective displacement
Frequency band analysis 2
1.5
1
0.5
0 250
260
270
280
290
300
310
320
330
340
350
Frequency (Hz)
Fig. 26. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 = 2p · 300 Hz, Dx = 2p · 50 Hz.
–5
2
x 10
1.8 1.6
Effective displacement
1.4 1.2 1 0.8 0.6 0.4 0.2 0 50
100
150
200
250
300
350
Frequency (Hz)
Fig. 27. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 2 2p · {100, 200, 300} Hz, Dx = 2p · 50 Hz.
3190
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
x = 2p · 350 Hz. This figure shows that the method works on two different types of structural behaviour: highly dynamic behaviour (frequency range [50; 250 Hz]) and relatively smooth behaviour (frequency range [250; 350 Hz]). 7. Conclusion The approach proposed here, called the ‘‘variational theory of complex rays’’, was originally introduced in order to calculate the vibrations of slightly damped elastic structures in the medium-frequency range. It is a general multiscale approach with a strong mechanical basis. This paper emphasizes the extension of this theory to the analysis across a range of frequencies. The method used to find the solution within a frequency range consists in expanding the matrix and the loading of the system of equations to be solved into a Taylor series. This can be done by separating rapidly-varying quantities and slowly-varying quantities with respect to the frequency. Then, by using an original decomposition of the solution into its mean value over the frequency range and the complementary part, the method can derive an approximate solution for all frequencies in the frequency range. The numerical examples showed the capability of the VTCR to predict the vibrational response of a structure within a frequency range. In this work, the structure is supposed to be an assembly of homogeneous substructures, but this hypothesis can easily be removed [23]. References [1] P.E. Barbone, J.M. Montgommery, O. Michael, I. Harari, Scattering by a hybrid asymptotic/finite element method, Comput. Methods Appl. Mech. Engrg. 164 (1998) 141–156. [2] V.D. Belov, S.A. Ryback, Applicability of the transport equation in the one-dimensional wave propagation problem, Akust. Zh. Sov. Phys. Acout. 21 (1975) 173–180. [3] P. Bouillard, V. Lacroix, E. De Bel, A wave-oriented meshless formulation for acoustical d and vibro-acoustical applications, Wave Motion 39 (4) (2004) 295–305. [4] O. Cessenat, B. Despres, Application of an ultra weak variational formulation of elliptic pdes to the two-dimensional Helmholtz problem, SIAM J. Numer. Anal. 35 (1) (1998) 255–299. [5] E.A. Chadwick, P. Bettess, Modelling of progressive short waves using wave envelopes, Int. J. Numer. Methods Engrg. 14 (1997) 3229–3246. [6] R.D. Ciskowski, C.A. Brebbia, Boundary Element Methods in Acoustics, Elsevier Applied Science, London, 1991. [7] A. Deraemaeker, I. Babuˇska, P. Bouillard, Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions, Int. J. Numer. Methods Engrg. 46 (1999) 471–499. [8] W. Desmet, B. Van Hal, P. Sas, D. Vandepitte, A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems, Adv. Engrg. Software 33 (2002) 527–540. [9] R. Djellouli, C. Farhat, R. Tezaur, A fast method for solving acoustic scattering problems in frequency bands, J. Comput. Phys. 168 (2) (2001) 412–432. [10] J.F. Doyle, Wave Propagation in Structures, Springer Verlag, 1997. [11] C. Farhat, I. Harari, U. Hetmaniuk, A discontinuous Galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the mid-frequency regime, Comput. Methods Appl. Mech. Engrg. 192 (2003) 1389–1419. [12] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residual-free bubbles for the Helmholtz equation, Int. J. Numer. Methods Engrg. 40 (1997) 4003–4009. [13] I. Harari, T.J.R. Hughes, Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains, Comput. Methods Appl. Mech. Engrg. 98 (3) (1992) 411–454. [14] I. Harari, S. Haham, Improved finite element method for elastic waves, Comput. Methods Appl. Mech. Engrg. 166 (1998) 143– 164. [15] I. Harari, P. Barai, P.E. Barbone, Numerical and spectral investigations of Trefftz infinite elements, Int. J. Numer. Methods Engrg. 46 (4) (1999) 553–577. [16] T.J.R. Hughes, Multiscale phenomena: GreenÕs functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1–4) (1995) 387–401. [17] M.N. Ichchou, A. Le Bot, L. Je´ze´quel, Energy model of one-dimensional, multipropative systems, J. Sound Vib. 201 (1997) 535– 554.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3191
[18] F. Ihlenburg, I. Babuˇska, Dispersion analysis and error estimation of Galerkin finite element methods for Helmoltz equation, Int. J. Numer. Methods Engrg. 38 (1995) 3745–3774. [19] P. Ladeve´ze, A new computational approach for medium-frequency vibrations (in french), NT Ae´rospatiale YX/SA (1995) 119– 139. [20] P. Ladeve`ze, A new computational approach for structure vibrations in the medium frequency range, C.R. Acad. Sci. Paris 322, Se´rie II b No. 12 (1996) 849–856. [21] P. Ladeve`ze, L. Arnaud, P. Rouch, C. Blanze´, The variational theory of complex rays for the calculation of medium-frequency vibrations, Engrg. Comp. 18 (1999) 193–214. [22] P. Ladeve`ze, L. Arnaud, A new computational method for structural vibrations in the medium frequency range, Comput. Assisted Mech. Engrg. Sci. 7 (2000) 219–226. [23] P. Ladeve`ze, L. Blanc, P. Rouch, C. Blanze´, A multiscale computational method for medium-frequency vibrations of assemblies of heterogeneous plates, Comput. Struct. 81 (12) (2003) 1267–1276. [24] P. Ladeve`ze, P. Rouch, H. Riou, X. Bohineust, Analysis of medium-frequency vibrations in a frequency range, J. Comput. Acoust. 11 (2) (2003) 255–284. [25] R.S. Langley, A wave intensity technique for the analysis of high frequency vibrations, J. Sound Vib. 159 (1992) 483–502. [26] Y. Lase, M.N. Ichchou, L. Je´ze´quel, Energy flow analysis of bars and beams: theoretical formulation, J. Sound Vib. 192 (1996) 281–305. [27] W.K. Liu, Y.F. Zhang, M.R. Ramirez, Multiple scale finite element methods, Int. J. Numer. Methods Engrg. 32 (1991) 969–990. [28] R.H. Lyon, G. Maidanik, Power flow between linearly coupled oscillators, J. Acoust. Soc. Am. 34 (1962) 623–639. [29] M. Malhotra, P.M. Pinsky, Efficient computation of multi-frequency far-field solutions of the Helmholtz equation using pad approximation, J. Comput. Acoust. 8 (1) (2000) 223–240. [30] J.M. Melenk, I. Babusˇka, The partition of unity finite element method: basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289–314. [31] P. Monk, D.Q. Wang, A least-squares method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 175 (1999) 121–136. [32] J.P.H. Morand, A modal hybridization method for the reduction of dynamic models in the medium frequency range, in: Proc. European Conference on New Advances in Computational Structural Mechanics, 1991, pp. 347–365. [33] E. Perrey-Debain, J. Trevelyan, P. Bettess, Plane wave interpolation in direct collocation boundary element method for radiation and wave scattering: numerical aspects and applications, J. Sound Vib. 261 (2003) 839–858. [34] H. Riou, P. Ladeve´ze, P. Rouch, Extension of the variational theory of complex rays to shells for medium-frequency vibrations, J. Sound Vib. 272 (1–2) (2004) 341–360. [35] P. Rouch, P. Ladeve´ze, The variational theory of complex rays: a predictive tool for medium-frequency vibrations, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3301–3315. [36] C. Soize, Medium frequency linear vibrations of anisotropic elastic structures, Recherche Ae´rospatiale 5 (1982) 65–87. [37] C. Soize, P.M. Hutin, A. Desanti, J.M. David, F. Chabas, Linear dynamic analysis of mechanical systems in the medium frequency range, J. Comput. Struct. 23 (1986) 605–637. [38] C. Soize, Reduced models in the medium frequency range for the general dissipative structural dynamic systems, Eur. J. Mech., A/ Solids 17 (1998) 657–685. [39] T. Strouboulis, K. Copps, I. Babusˇka, The generalized finite element method, Comput. Methods Appl. Mech. Engrg. 190 (32–33) (2001) 4081–4193.
Calculation of medium-frequency vibrations over a wide frequency range P. Ladeve`ze *, H. Riou LMT Cachan (ENS Cachan/CNRS/Paris 6 University), 61 Avenue du President Wilson, 94235 Cachan Cedex, France Received 17 March 2004; received in revised form 27 August 2004; accepted 30 August 2004
Abstract An enhancement of the ‘‘variational theory of complex rays’’ (VTCR) is being developed in order to calculate the vibrations of slightly damped elastic structures in the medium-frequency range. Here, the emphasis is put on the extension of this theory to the analysis over a wide frequency range. The solution is sought in the form of an original decomposition: the mean value and its complementary part over the frequency range. Numerical examples show the capability of the VTCR to predict the vibrational response of a structure in a frequency range. 2004 Elsevier B.V. All rights reserved. Keywords: Medium-frequency vibrations; Frequency range; Multiscale numerical strategy
1. Introduction The calculation of medium-frequency vibrations is a key element in the design of car chassis or satellites and it represents a real challenge. The difficulty resides in the fact that the length of variation of the phenomena being studied is very small in relation to the characteristic dimension of the structure. Therefore, the use of an approximation based on continuous, piecewise polynomial shape functions leads to models with huge numbers of degrees of freedom (DOFs). Furthermore, pollution effects occur and affect the robustness of the finite element method (FEM) with respect to the wave number k [7,18]. Different approaches have been proposed in order to circumvent this problem. These approaches include predefined reduced bases [32,38], the Galerkin/least-squares FEM [13], the partition of unity method [30], the generalized finite element method [3,39], various multiscale finite element methods [16,27] such as the
*
Corresponding author. Tel.: +33 1 4740 2253; fax: +33 1 4740 2785. E-mail addresses: [email protected] (P. Ladeve`ze), [email protected] (H. Riou).
0045-7825/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2004.08.009
3168
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
residual-free bubbles method [12], the ultra weak variational method [4], the wave envelope method [5], the discontinuous enrichment method [11], the Trefftz methods [8,14,15,31], and the boundary element method [6,33]. However, all these methods are essentially limited to low-frequency problems. Indeed, the use of element-based methods in the medium-frequency range requires that the number of interpolation functions be increased in order to represent the dynamic behaviour correctly. For medium frequencies, these methods are either inaccurate or costly. Statistical energy analysis (SEA, [28]) involves the description of the energy exchanged among various systems and yields the dynamic response of each system averaged over time and space. However, the SEA, which is characterized by a single energy level, cannot provide the spatial variation of the response within each system. Extensions of the SEA have been developed [2,17,25,26], but these extensions still require additional information (e.g. coupling loss factor, energy reflection coefficient, energy transmission coefficient . . .) in order to yield predictive results; therefore, they have been developed only for specific geometries, such as bars or beams. Many improvements still have to be made in order for these methods to become applicable to the medium-frequency range. Moreover, regardless of the methods chosen for discretizing the equations, the resolution of a vibrational problem with multiple frequencies involves the resolution of a series of linear systems of equations with different left- and right-hand sides. The resolution of such systems by a direct method requires the inversion of a number of matrices equal to the number of specified frequencies, which leads to very CPU-intensive computations. Therefore, some techniques have been developed in order to maximize the computational efficiency of the calculations over frequency bands [9,29,36]. These works showed that TaylorÕs, Pade´Õs or WynnÕs approximations of the solution with respect to the frequency seem to be efficient. Nevertheless, these approaches are based on a finite element model of the problem and, therefore, they involve large matrices in the medium-frequency range. The alternative approach discussed here is called the ‘‘variational theory of complex rays’’ (VTCR) [20,22]. It is a predictive tool specifically intended for dealing with these kinds of problems. Its fundamental aspects are described in [21]. The first characteristic feature of the VTCR is that it uses a new variational formulation of the problem being considered which was developed in order to allow the approximations within the substructures to be a priori independent of one another: in other terms, it is not necessary for these approximations to verify a priori the boundary conditions and the continuity conditions at the interfaces between substructures. Instead, these conditions are included in the variational formulation. The second feature which characterizes the VTCR is the introduction of two-scale approximations with a strong mechanical meaning: the solution is assumed to be properly described locally as the superposition of an infinite number of local vibration modes over each substructures, which are supposed to be homogeneous. Each basic mode (which can be an interior ray, an edge ray or a corner ray) associated with a substructure verifies the governing equation and the constitutive law over the substructureÕs domain. All wave directions are taken into account. The unknowns are discretized amplitudes with relatively long wavelengths. The third and last feature characterizing the VTCR is that one retains from the calculated discretized amplitudes only effective quantities related to the elastic energy, the kinetic energy, the dissipation work, the effective displacements, etc. Numerical examples of complex structures consisting of plate and shell assemblies have already demonstrated the capabilities of this formulation at a fixed frequency [34,35]. Recently, Ladeve`ze et al. [24] proposed a new algorithm for the calculation of multiple-frequency solutions using the VTCR. The proposed method uses a set of parameters q to derive a discrete approximation of the rapidly-varying quantities (exponential terms) inside the matrix of the VTCR. These parameters tell how the matrix equation of the VTCR varies in terms of the frequency. Then, three strategies are available in order to obtain the solution: the first is based on averaging over the frequency band; the second and third are based on Newton and
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3169
quasi-Newton algorithms. These strategies do not require the calculation of the solution at each frequency of the frequency band, which is a marked improvement over the standard VTCR. Nevertheless, as pointed out in [24], they require a large number of parameters q to represent the rapidly-varying exponential terms in the matrix properly and, therefore, they lead to large numerical models. In this paper, we present an improvement of the VTCR which enables one to calculate the response in a frequency range using only a small number of parameters to represent the rapidly-varying quantities of the VTCR matrix equations. In the technique presented, the matrix and the right-hand side of the system to be solved are expanded into Taylor series with respect to the frequency. These series are calculated at the center of the frequency range. Then, an ad hoc strategy uses the mean value and the complementary part of the solution over the frequency range to predict the behaviour of effective quantities, such as spatial averages, which leads to the resolution of a system of equations with multiple right-hand sides; this resolution is computationally efficient because the matrix of the system needs to be factored only once. The remainder of this paper is organized as follows: the reference problem to be solved is presented in Section 2; the VTCR as a general approach for solving fixed-frequency vibration problems is reviewed in Section 3; the principle of the approximate formulation, the discretized problem and the numerical strategy are outlined in Section 4; the basic features of the proposed calculation method in a wide frequency range using the VTCR are described in Section 5; in Section 6, numerical illustrations are given and a comparison with the strategy proposed in [36] is made using the example of a simply-supported plate; finally, the conclusion is presented in Section 7.
2. The reference problem In order to simplify the presentation, the problem will be formulated for an assembly of two substructures, but it is easy to generalize to an assembly of n substructures. Given two substructures S and S 0 , let oS and oS 0 represent the boundaries of S and S 0 respectively. We are studying the harmonic vibration of these two structures at a fixed frequency x. All quantities can be defined in the complex domain. An amplitude Q(X) is associated with Q(X)eixt. The excitations applied to S and represented in Fig. 1 are 1. a displacement field Ud on a portion oUS of the boundary oS, 2. a density of forces Fd on a portion oFS, 3. a density of forces f d over the whole domain S.
Fig. 1. The reference problem.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3170
Similar quantities are defined on S 0 . For a structure S, let Sad denote the field of the pairs (U, r) such that: 3
U 2U
ðfinite-energy displacement set ½H 1 ðSÞ Þ;
r2S
ðfinite-energy stress set ½L2 ðSÞ Þ;
divr þ f d ¼ x2 qð1 ihÞU r ¼ ð1 þ igÞKðU Þ
on S; on S;
3
ð1Þ
where K is HookeÕs operator, q the mass density, and h and g the damping coefficients (which depend on the frequency). The subspace Sad associated with homogenized conditions (f d = 0) is designated by Sad,0. In the same manner, we introduce the spaces S0ad and S0ad;0 . The reference problem can be formulated as follows: find (U(X), r(X), X 2 S) and (U 0 (X), r 0 (X), X 2 S 0 ) such that s ¼ ðU ; rÞ 2 Sad ; U ¼ U d on oU S; rn ¼ F d U ¼ U0
on oF S;
s0 ¼ ðU 0 ; r0 Þ 2 S0ad ; U 0 ¼ U 0d on oU S 0 ; 0 0
rn ¼
F 0d
ð2Þ
0
on oF S ;
on C;
rn þ r0 n0 ¼ 0 on C; where n (resp. n 0 ) is the outward normal vector of the boundary of S (resp. S 0 ).
3. Variational formulation associated with the VTCR This section reviews the main characteristics of the VTCR. For more details, the reader is referred to [20]. The VTCR is a global formulation of the boundary conditions (2) in terms of both displacements and forces. It is based on a priori independent approximations within the substructures: find s = (U, r) 2 Sad and s0 ¼ ðU 0 ; r0 Þ 2 S0ad such that: Z Z Z drnðU U d Þ dl þ ðrn F d ÞdU dl þ dr0 n0 ðU 0 U 0
Re ix d Þ dl 0 oU S oF S oU S Z Z 1 0 0 0 0
0 0
0
0 0
0
þ ðr n F d ÞdU dl þ ðdrn dr n ÞðU U Þ þ ðrn þ r n ÞðdU þ dU Þ dl ¼ 0 2 C oF S 0 8ðdU ; drÞ 2 Sad;0
8ðdU 0 ; dr0 Þ 2 S0ad;0 ;
ð3Þ
where Re[A] denotes the real part of quantity A and A* the conjugate of A. It is easy [20] to prove that the variational form is equivalent to the reference problem provided that: 1. the reference problem (2) has a solution, 2. HookeÕs operators K and K 0 are positive definite, 3. one of the damping factors h, g (resp. h 0 , g 0 ) is greater than zero on S (resp. S 0 ). Contrary to the usual formulations, displacements and forces may be applied simultaneously on the same boundary oLS. This enables one, for example, to take into account all the data of the problem available from an experiment in which both the forces and the displacements may be known at the same boundary. In this case, the following term should be considered:
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
1 Re ix 2
Z
drnðU U d Þ þ ðrn F d ÞdU dl
3171
:
oL S
It is possible [20] to express the formulation as: find s = (U, r) 2 Sad and s0 ¼ ðU 0 ; r0 Þ 2 S0ad such that: s s s 0 0 xdðED ðU Þ þ ED ðU ÞÞ þ 8ds0 2 S0ad;0 ; ;d 0 ¼ LD ; d 0 8ds 2 Sad;0 0 s s s where ED is the dissipated energy, LD a linear form, and h. , .i a bilinear form defined at the boundaries which verifies hu, vi = hv, ui*.
4. Approximate formulations 4.1. Principle: the general case All that one needs to do in order to derive approximate formulations from the VTCR is to define the 0h 0 0 subspaces Shad and Shad;0 (resp. S0h ad and Sad;0 ) from Sad and Sad,0 (resp. Sad and Sad;0 ). The approximate forh 0h h 0h h h 0h 0h mulation can be expressed as: find s ¼ ðU ; r Þ 2 Sad and s ¼ ðU ; r Þ 2 Sad such that: h h h s s s xdðED ðU h Þ þ E0D ðU 0h ÞÞ þ ; d 8ds0h 2 S0h ; d ¼ L 8dsh 2 Shad;0 D ad;0 : s0h s0h s0h The VTCR uses two approximation scales (Uh, rh), which have a strong mechanical meaning, by distinguishing three zones: the interior zone, the edge zone and the corner zone. For example, in the neighborhood of a point X of the interior zone, the solution is assumed to be properly described locally as the superposition of an infinite number of complex rays (propagative and evanescent waves) which can be written as follows: U h ðX ; Y ; P Þ ¼ W h ðX ; P ÞeixP Y ; rh ðX ; Y ; P Þ ¼ Ch ðX ; P ÞeixP Y ; where both X and Y represent the position vector, X being associated with slow variations and Y with rapid variations. More precisely, the terms related to the position vector X vary slowly when X moves along the structure, whereas the terms related to the position vector Y vary rapidly when Y moves along the structure. P is a vector which characterizes the complex ray. In order for these local modes (Uh, rh) to be admissible, they must lie in Shad and verify (1). Thus, we obtain some properties on P. These properties enable one to classify the rays into interior rays (propagative waves), and edge and corner rays (evanescent waves). 4.2. Admissible complex rays for a plate Let us consider the out-of-plane bending motion of a thin, flat, homogeneous and isotropic plate. According to KirchhoffÕs thin plate theory, the steady-state displacement u of the plateÕs mid-surface perpendicularly to the plate is governed by the dynamic equation: Eh3 ð1 þ igÞMMu ¼ qhx2 ð1 ihÞu on S; 12ð1 m2 Þ where n is the Laplacian operator, E the YoungÕs modulus, h the thickness of the plate, m the PoissonÕs ratio, q the mass density, x the frequency, and h and g the damping factors. A complex ray for the interior zone is
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3172
uhi ðX ; Y ; P Þ
¼
whi ðX ; P Þ exp
pffiffiffi h þ g pffiffiffiffi xP X ei xP Y : 4
This complex ray corresponds to a plane bending wave which propagates through the plate in direction P. This ray is admissible (i.e. it belongs to Sad) only if: Eh3 ð1 þ igÞMMuhi ¼ qhx2 ð1 ihÞuhi 12ð1 m2 Þ
on S:
Therefore, the properties on P are 12qð1 m2 Þ 2 ðP P Þ ¼ r4 with r4 ¼ : h2 E
ð4Þ
Eq. (4) shows that P lies on a circle Ci with a radius r defined by the material properties. By following this circular path, one takes into account all the directions of the plate. Examples of such interior rays are shown in Fig. 2. A similar approach can be used for the edge and corner zones. The properties on P for an edge complex ray of the form pffiffiffi pffiffiffi h þ g pffiffiffiffi r x n X e xP n nY þi xP t tY ; uhe ðX ; Y ; P Þ ¼ whe ðX ; P Þ exp 4 Pn which verifies Eh3 ð1 þ igÞMMuhe ¼ qhx2 ð1 ihÞuhe 12ð1 m2 Þ
on S
and does not propagate faster than the interior rays can be written as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P t ¼ r sinðuÞ P n ¼ r 1 þ sin2 ðuÞ;
ð5Þ
(where n and t denote respectively the normal and tangent vectors to the edge). Eq. (5) shows that P lies on a curve Cb, defined by the material properties, spanned by u. By following this curve, one takes into account all edge complex rays. Examples of such modes are shown in Fig. 3 (left and center). The properties on P for a corner complex ray of the form pffiffiffi h þ g pffiffiffiffi h h xP X e xP Y ; uc ðX ; Y ; P Þ ¼ wc ðX ; P Þ exp i 4 which verifies Eh3 ð1 þ igÞMMuhc ¼ qhx2 ð1 ihÞuhc 12ð1 m2 Þ
on S
Fig. 2. Admissible interior rays for a plate.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3173
Fig. 3. Admissible edge modes (left and center) and corner modes (right) for a plate.
can be written as 2
ðP P Þ ¼ r4
with r4 ¼
12qð1 m2 Þ : h2 E
ð6Þ
Eq. (6) shows that P lies on the same circle, denoted Cc, as in the case of interior complex rays. Examples of such modes are shown in Fig. 3 (right). One must note that the admissible complex rays are exact solutions of the governing equation of the problem. Therefore, unlike conventional finite element-based methods, the VTCR yields no numerical dispersion [7,18]. 4.3. The discretized problem The displacement of each point in the substructure is generated by a basis of admissible complex rays. The unknowns are the generalized amplitudes whi ðX ; P Þ; whb ðX ; P Þ and whc ðX ; P Þ of the basis (nth-order polynomials in X and large-wavelength quantities). Taking into account all the directions u in Ci, Cb and Cc, one ends up with an integral which takes the form: Z pffiffiffi hþgpffiffiffi whi ðX ; P ðuÞÞe 4 xP ðuÞX ei xP ðuÞY dC uh ðX ; Y Þ ¼ u2C i
þ þ
Z
Z
whb ðX ; P ðuÞÞexp
ffiffiffi
hþgp r xP n nX 4
e
pffiffiffi pffiffiffi xP n nY þi xP t tY
dC
u2C b
whc ðX ; P ðuÞÞei
ffiffiffi
hþgp xP ðuÞX 4
e
pffiffiffi xP ðuÞY
dC:
ð7Þ
u2C c
In order to obtain a finite-dimension problem, one can discretize this integral (7) and assume that the amplitudes whi ðX ; P ðuÞÞ; whb ðX ; P ðuÞÞ and whc ðX ; P ðuÞÞ are constant over angular sectors: whi ðX ; P ðEP ÞÞ; whb ðX ; P ðEP ÞÞ and whc ðX ; P ðEP ÞÞ (for example, see Fig. 4 for the interior rays). The angular distribution of the plane waves for any point in the substructure is assumed to be properly described by such a discontinuous angular distribution. Ci, Cb and Cc is discretized into N elements, which can be of different sizes (on Fig. 4, Ci is discretized into eight elements of the same size). Once the discretization has been chosen for each plate, the VTCR leads to a system of linear equations in the complex domain: K U ¼ F;
ð8Þ
where K ¼ xEhD þ Z h and F ¼ LhD . EhD is the symmetric, positive definite damping matrix associated with ED. Zh is the matrix associated with the bilinear form h. , .i (see Section 3), defined such that Z Th ¼ Z h . LhD is the vector associated with the linear form LD. U is the vector corresponding to the unknown amplitudes associated with the complex polynomial wh(X, EP). As a consequence of the last properties, Eq. (8) has a unique solution.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3174
Fig. 4. Discretized amplitudes.
4.4. The numerical strategy The system of linear equations leads to the solution of (8). K corresponds to the bilinear terms in the h 0h variational formulation (3) projected onto ðShad [ S0h ad Þ ðSad;0 [ Sad;0 Þ, U is the vector of the unknown genh h h eralized amplitudes wi ðX ; P ðEP ÞÞ, wb ðX ; P ðEP ÞÞ and wc ðX ; P ðEP ÞÞ, and F corresponds to the projection onto ðShad [ S0h ad Þ of the load terms of the variational formulation (3). A typical example of a term of K would be the coupling of an interior sector EP with an edge sector E0P on a plate S, as shown in Fig. 5. The coupling of these two angular sectors leads to the following integral: Z Z Z I¼ aðP ; P 0 ; X Þ exp½ixðP P 0 Þ Y du0 du dS; ð9Þ S
EP
E0P
Fig. 5. Typical term of the symmetric part of the coupling matrix.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3175
where a is a polynomial, and P and P 0 are the wave vectors belonging to the two angular sectors EP and E0P . The structure S is subdivided into elements for integration purposes. The size of these elements is chosen to be approximately one wavelength, so that slowly varying quantities can be considered to be constant. Fig. 6 shows the spatial decomposition used for the integration. Then, Eq. (9) becomes: Z Z Z X I¼ aðP ; P 0 ; X Þ exp½ixðP P 0 Þ Y du0 du dD: ð10Þ subdomains D
D
EP
E0P
Assuming the polynomial a to be constant over each subdomain D and denoting XE the center of the cell D, one can rewrite Eq. (10) as Z Z Z X I¼ aðP ; P 0 ; X E Þ exp½ixðP P 0 Þ Y du0 du dD: ð11Þ subdomains D
EP
E0P
D
With a change of variable transforming D into a symmetric domain in N1 and N2, Eq. (11) becomes: Z Z X I¼ aðP ; P 0 ; X E Þ exp½ixðP P 0 Þ X E du0 duJ ð12Þ subdomains D
EP
E0P
with J¼
Z
h1 2 h
21
exp½ixðP 1
¼ h1 Sinc xðP 1
P 01 ÞX 1 dX 1
h1 P 01 Þ 2
Z
h2 2 h
22
exp½ixðP 2 P 02 ÞX 2 dX 2
h1 Sinc xðP 2
h2 P 02 Þ 2
;
where Sinc ½x ¼
sin½x : x
Now, one must integrate over the angular domains EP and E0P . Using the same technique, one decomposes these domains into subdomains within which the polynomial can be considered to be constant. Fig. 7 shows this angular decomposition.
Fig. 6. Spatial decomposition for the integration.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3176
Fig. 7. Angular decomposition for the integration.
Denoting ehP and eh0P the angular subdomains, I becomes: I¼
X subdomains D
with
2
X X Z
4
ehP 2EP eh0 2E0 P P
ehP
3
Z 0
ehP
H du du0 5
ð13Þ
0 h1 0 h2 H ¼ h1 h2 aðP ; P ; X E Þ exp½ixðP P Þ X E Sinc xðP 1 P 1 Þ Sinc xðP 2 P 2 Þ : 2 2 0
0
Once again, one calculates the quantities subjected to rapid variations. These quantities are integrated analytically, whereas the others are calculated numerically. Because of the domains of variation of the two functions Sinc, these can be treated as slow functions and, by linearizing the angular dependency of the wave vectors P and P 0 as shown in Fig. 8.
Fig. 8. Linearization of the wave vector P.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
One can write Eq. (13) as Z Z Z 0 0 G¼ exp½ixðP P Þ X E du du ¼ ehP
0
ehP
¼ exp½ixðP h P 0h Þ X E
s1
exp½ixP X E ds
s1
Z
s1
exp½ixsT h X E ds
s1
Z
s01 s01
Z
s01 s01
3177
exp½ ixP 0 X E ds0
exp½ ixsT 0h X E ds0
¼ 4s1 s01 exp½ixðP h P 0h Þ X E Sinc ½xs1 T h X E Sinc ½xs01 T 0h X E : Since the angular dependencies of the edge and corner modes are constructed in the same way as those of the interior modes (see Section 4.1), their linearization yields exactly the same results as in Fig. 8. Finally, Eq. (13) becomes: X X X I¼ 4h1 h2 s1 s01 aðP h ; P 0h ; X E Þ subdomains D eh 2EP eh0 2E0 P P P
h1 h2 Sinc xðP h1 P 0h1 Þ Sinc xðP h2 P 0h2 Þ Sinc ½xs1 T h X E 2 2 Sinc ½xs01 T 0h X E exp½ixðP h P 0h Þ X E ;
which, because of the slow evolution of the functions Sinc in their respective domains, can be written as X X X ½AðP h ; P 0h ; X E Þ exp½ixðP h P 0h Þ X E ; ð14Þ I¼ subdomains D eh 2EP eh0 2E0 P P P
where A represents a function which contains all the slow-variation terms. All the rapid-variation terms are included in the exponential functions. This particular decomposition will be used to calculate the solution within a frequency range. 4.5. Effective quantities From a mechanical point of view, the spatial distribution of the solution in the medium-frequency range has no physical meaning. Indeed, the response at a specific location and at a specific frequency is extremely sensitive to perturbations. Spatially averaged quantities are more useful for interpreting the response. Therefore, it is necessary to extract from the solution some effective quantities defined over a domain D greater than one wavelength. These quantities extracted from the local data A are called effective quantities and are defined as Z 1 Aeff;D ¼ jAðX Þj dX : D D From here on, we will use the effective displacement over a domain D: Z 1 U eff;D ¼ jU ðX Þj dX : D D 5. Calculation in a frequency range The objective of this section is to present a variation of the VTCR capable of calculating the solution in a frequency range B defined by its central frequency x0 and its bandwidth Dx: find U which verifies KðxÞUðxÞ ¼ FðxÞ 8x 2 B ¼ ½x0 Dx; x0 þ Dx:
ð15Þ
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3178
This approach, first introduced in [19], uses the decomposition presented in Eq. (14). Considering the nature of the frequency range, let us first assume that the function A is constant on B. Consequently, one has AðP h ; P 0h ; X E ; xÞ ’ AðP h ; P 0h ; X E ; x0 Þ
for x 2 B:
If necessary, one could build a better approximation (for example, a linear approximation): AðxÞ ’ Aðx0 Þ þ
Aðxþ Þ Aðx Þ ðx x0 Þ xþ x
for x 2 B;
where x+ and x represent the limits of the frequency range. Let us now study the behaviour of the exponential argument. This quantity can be expanded into a Taylor series up to order k: " # k X ðiðP h P 0h Þ X E Þk k 0 0 ðx x0 Þ : exp½ixðP h P h Þ X E ¼ exp½ix0 ðP h P h Þ X E k! k¼0 Then, according to Section 4, the matrix K and the load F can be approximated as KðxÞ ¼
k X
Kk ðx0 Þðx x0 Þk ;
k¼0
FðxÞ ¼
k X
ð16Þ k
Fk ðx0 Þðx x0 Þ :
k¼0
Let us define the mean value over the frequency range B: Z 1 dx: hi ¼ 2Dx B The matrix K and the load F can be expressed as KðxÞ ¼ hKi þ DKðxÞ;
ð17Þ
FðxÞ ¼ hFi þ DFðxÞ: According to the Taylor expansion (16), one has DKðxÞ ¼
k X
k
DKk ðx0 Þðx x0 Þ ;
k¼0
DFðxÞ ¼
k X
k
DFk ðx0 Þðx x0 Þ :
k¼0
One can also define U using the same approach: UðxÞ ¼ hUi þ DUðxÞ:
ð18Þ
By combining (15), (17) and (18), one gets ½hKi þ DK ½hUi þ DU ¼ ½hFi þ DF:
ð19Þ
Eq. (19) involves terms with very different magnitudes: some are very large while others are very small. Using the techniques of perturbation methods, let us now identify terms of different orders and rewrite Eq. (19) for order 0 and order 1:
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
hKihUi ¼ hFi
order 0;
hKiDU ¼ DF DKhUi
order 1:
3179
ð20Þ
Then, DU is equal to DU ¼
k X
k
DUk ðx x0 Þ :
ð21Þ
k¼0
Finally, the displacement over the frequency range B which we retain is UðxÞ ¼ khUi þ lDU: where k and l minimize the error defined through the approximations in (20): Z 1 T E2 ¼ ½F K U ½KS D ½F K U dx; 2Dx B where [KS]D is the diagonal of the symmetric part of K, which is a slowly-varying quantity [20]. Remarks: (1) One can write U = hKi 1L(x x0), where L denotes an operator. The strategy we described can be viewed as the resolution of a finite-dimension linear problem with multiple right-hand sides. (2) One must note that the only parameter which must be chosen in order to ensure good accuracy is k, which corresponds to the degree of the Taylor expansion. (3) The mean dissipated energy over the frequency range is Z k2 l2 e ¼ hUihKS ihUi þ DUKS DU dx: 2 4Dx B (4) In [9], Pade´Õs or WynnÕs approximations seem to improve the interval of convergence of the Taylor series. In this paper, Taylor series were chosen for four reasons: (a) in [29] and [9], the whole vibrational solution was discretized, whereas in the VTCR only the slowly-varying part, i.e., the amplitude of the waves, a quantity less sensitive to a variation of frequency, is discretized (in the examples in the next section, only five terms of the Taylor series suffice to accurately represent the amplitude of the waves across a range of frequencies), (b) as pointed out in [9], if the number of approximation terms is large, Pade´Õs approximation leads to an ill-conditioned matrix and WynnÕs algorithm becomes computationally expensive, (c) Taylor series are easier to implement, and (d) the proposed strategy is not limited to a Taylor series expansion, as k and l can be x-varying functions.
6. Numerical examples 6.1. A simple example to test the method over the frequency range To illustrate the approach discussed above, let us consider the dynamic response of a beam subjected to an external excitation (see Fig. 9). The purpose of this simple example is to test the method without introducing any approximation, since for one-dimensional structures the VTCR yields the exact solution at a fixed frequency. We used the method to find the reference solution over the frequency range [800; 4000 Hz]. The reference solution is the analytical solution. The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz.
3180
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
F
Fig. 9. Model geometry for a beam in bending; E = 200 GPa; m = 0.3; g = 0.05; h = 0; q = 7800 kg/m3; S = 0.0001 m2; I = 8.3 · 10 10 m4; L = 2.5 m; F = 1 N/m.
In Figs. 10–12, the energy density at the center of the beam is plotted against the frequency for different orders of Taylor series expansion (see Section 5, Eq. (16)). One can see that the solution obtained with the frequency range calculation is very good, even on the portions of the curves with highly dynamic behaviour. As the order of the Taylor series expansion increases, the precision improves. At order 4, the approximate solution matches the reference solution very closely. This can also be observed in Fig. 13, which shows the error as a function of the expansion order: as expected, as the order increases, the error decreases. Figs. 14–16 show the solution of the problem at 2610 Hz (which is not the central frequency of a frequency range) for different orders of Taylor series expansion. These results show that the accuracy increases as the expansion order increases. Again, one can see that order 4 provides sufficient accuracy. In the following examples, all the approximate frequency range solutions were calculated to order 4. 6.2. Plate example This second example shows the effectiveness of the VTCR in the case of a plate structure whose geometry is shown in Fig. 17. The out-of-plane bending motion is induced by an external force applied in the middle
0. 2 0.18
Computation order 0 Exact theory
0.16
0.12 0. 1
u
E (L/2) (J/m)
0.14
0.08 0.06 0.04 0.02 0 500
1000
1500
2000
2500
3000
3500
4000
Frequency ω (Hz)
Fig. 10. Result for the beam problem (see Fig. 9). The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz and a Taylor series expansion up to order 0 (see Section 5, Eq. (16)).
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3181
0.2 0.18
Computation order 2 Exact theory
0.16
u
E (L/2) (J/m)
0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 500
1000
1500
2000
2500
3000
3500
4000
Frequency ω (Hz)
Fig. 11. Result for the beam problem (see Fig. 9). The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz and a Taylor series expansion up to order 2 (see Section 5, Eq. (16)).
0.2 0.18
Computation order 4 Exact theory
0.16
0.12 0.1
u
E (L/2) (J/m)
0.14
0.08 0.06 0.04 0.02 0 500
1000
1500
2000
2500
3000
3500
4000
Frequency ω (Hz)
Fig. 12. Result for the beam problem (see Fig. 9). The approximate solution was calculated in successive frequency ranges with bandwidth Dx = 2p · 50 Hz and a Taylor series expansion up to order 4 (see Section 5, Eq. (16)).
of the plate, perpendicular to the plateÕs mid-surface. The exact analytical solution for this problem is given in [10].
3182
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191 1
0.9
Energy relative error (x=L/2)
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0 1
1.5
2
2.5
3
3.5
4
4.5
5
k
Fig. 13. Relative energy error for the beam problem (see Fig. 9).
0.04
0.035
Computation order 0 Exact theory
0.025
0.02
u
E (L/2) (J/m)
0.03
0.015
0.01
0.005
0 0
0.5
1
1.5
2
2.5
x (m)
Fig. 14. Result for the beam problem (see Fig. 9) at x = 2p · 2610 Hz. The approximate solution was calculated with a Taylor series expansion up to order 0 (see Section 5, Eq. (16)).
This example was chosen so that we could compare our approach to that developed by Soize (1982) in [36], based on a low-frequency equation associated with the medium-frequency problem (LFEA). This method was used, in particular, in [37]. In [38], a numerical method for the analysis of medium-frequency
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3183
0.04
0.035
0.03
Computation order 2 Exact theory
E (L/2) (J/m)
0.025
u
0.02
0.015
0.01
0.005
0 0
0.5
1
1.5
2
2.5
x (m)
Fig. 15. Result for the beam problem (see Fig. 9) at x = 2p · 2610 Hz. The approximate solution was calculated with a Taylor series expansion up to order 2 (see Section 5, Eq. (16)).
0.04
0.035
Computation order 4 Exact theory
0.025
0.02
u
E (L/2) (J/m)
0.03
0.015
0.01
0.005
0 0
0.5
1
1.5
2
2.5
x (m)
Fig. 16. Result for the beam problem (see Fig. 9) at x = 2p · 2610 Hz. The approximate solution was calculated with a Taylor series expansion up to order 4 (see Section 5, Eq. (16)).
linear vibrations of elastic structures is proposed. This technique is used to characterized the vibratory state of a structure in frequency ranges. The finite element approximation of the solution sought in the frequency domain for an n-dimensional subspace is written as
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3184
F
Fig. 17. Model geometry for the plate; E = 210 GPa; m = 0.3; g = 0.001; h = 0; q = 7800 kg/m3; h = 0.005 m; dimensions 1 m · 1 m; F = 1 N.
Uðx; xÞ ¼
n X
U i ðxÞN i ðxÞ;
i¼1
where Ni denotes the finite element shape functions and Ui the ith nodal response quantity. The governing equation of motion is expressed as ð x2 M þ ixD þ KÞU ¼ F
8x 2 B ¼ ½x0 Dx; x0 þ Dx;
ð22Þ
where M, D, K and F are the mass, damping and stiffness matrices, and the force vector respectively. Since this medium-frequency equation (22) cannot be solved numerically using a time integration method (because the high frequencies involved would lead to too small a time step), the method consists in constructing a low-frequency equation using the new variable U0 ðxÞ ¼ Uðx þ x0 Þ;
3
x 10
–6
Exact solution
2.5
Effective displacement
VTCR 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 18. Result for the simply-supported plate problem (see Fig. 17) using the VTCR; Dx = 2p · 100 Hz and x0 2 2p · {500; 700; 900} Hz.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3185
–6
3
x 10
Exact solution
2.5
Effective displacement
LFEA 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 19. Result for the simply-supported plate problem (see Fig. 17) using the LFEA; Dx = 2p · 100 Hz and x0 2 2p · {500; 700; 900} Hz.
3
x 10
–6
Exact solution
2.5
Effective displacement
VTCR 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 20. Result for the simply-supported plate problem (see Fig. 17) using the VTCR; Dx = 2p · 300 Hz and x0 = 2p · 700 Hz.
3186
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3
x 10
–6
Exact solution
2.5
Effective displacement
LFEA 2
1.5
1
0.5
0 400
500
600
700
800
900
1000
Frequency (Hz)
Fig. 21. Result for the simply-supported plate problem (see Fig. 17) using the LFEA; Dx = 2p · 300 Hz and x0 = 2p · 700 Hz.
where U0 verifies the equations ð x2 M þ ixDn þ Kn ÞU0 ¼ F0 Dn ¼ D þ 2ix0 M; Kn ¼ K þ ix0 D x20 M;
8x 2 B0 ¼ ½ Dx; Dx; ð23Þ
F0 ðxÞ ¼ Fðx0 þ xÞ: The low-frequency Eq. (23) can be resolved using step-by-step time integration, which yields the desired solution of Eq. (22): X UðxÞ ¼ sL U0 ðmsL Þe imsL x 8x 2 B; m2Z ð24Þ p : sl ¼ Dx The parameters which must be chosen to ensure good accuracy are mI and mS, which correspond to the initial and final time integration in (24). To obtain the following results, we chose mI = 10 and mS = 40. This choice yielded good enough results with a reasonable calculation time. If we had taken larger parameters, the calculation time would have become really excessive. For the simply-supported plate example defined in Fig. 17, the two strategies (VTCR and LFEA) were applied in two cases. For the first case, one has Dx = 2p · 100 Hz and x0 2 2p · {500; 700; 900} Hz (Figs. 18 and 19). For the second case, one has Dx = 2p · 300 Hz and x0 = 2p · 700 Hz (Figs. 20 and 21). For the VTCR, we used 20 interior rays, 20 edge rays and no corner ray, i.e., 40 of freedom (DOFs). The Taylor series expansion (see Section 5) was taken up to order 4. For the LFEA, we calculated the solution according to the classical rule which says that one must use between 8 and 10 elements per wavelength (see [1]), which leads to 5766 DOFs.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3187
1
F Fig. 22. Model geometry for a 3D plate assembly; E = 200 GPa; m = 0.3; g = 0.001; h = 0; q = 7800 kg/m3; h = 0.004 m; F = 1 N/m.
In Figs. 18 and 19, one can see that both approaches match the exact analytical solution with good enough accuracy provided the bandwidth Dx is not too large (here, Dx = 2p · 100 Hz). In Figs. 20 and 21, one can see that neither approach leads to the exact analytical solution easily if the bandwidth becomes too large (here, Dx = 2p · 300 Hz). In the case of the VTCR, some peaks are shifted (see Fig. 20 at 900 Hz), whereas with the LFEA the levels of the peaks are not good enough (see Fig. 21). On these examples, we can observe that both strategies operate on the same frequency range. By comparing the numbers of DOFs required (40 DOFs for the VTCR needs, 5766 DOFs for the LFEA), one notes that the VTCR is less memory-consuming. 6.3. 3D plate assembly This third example shows the effectiveness of the VTCR in the case of a 3D structure intended to represent a car, whose geometry is shown in Fig. 22. The calculation using the VTCR involved 360 DOFs. The displacement obtained at x = 2p · 100 Hz is shown in Fig. 23 and compared with results from NASTRAN. The NASTRAN solution used QUAD4 elements and a mesh size of 10 elements per wavelength, yielding 33,230 DOFs. For this example, on a single-processor PC Linux system (INTEL 1.1 Ghz, 512 Mb RAM), it took 13 s to compute and inverse the matrix with the VTCR vs. 1 min and 13 s with NASTRAN. The two results are similar. The evolution of the effective displacement of Plate 1 of this assembly (the roof of the car; see Fig. 22) with the frequency of the excitation is shown in Figs. 24–26. The reference curve (dotted line) was calculated frequency by frequency. The approximate solution was calculated using the approach presented above (see Section 5) with a Taylor series expansion up to order 4. Again, these good results confirm the capability of the proposed method to find the solution over a frequency range. Fig. 27 shows the same results as in Figs. 24–26, this time collected next to one another. Thus, the frequency range extends from x = 2p · 50 Hz to
Fig. 23. Comparison of the displacements of the solutions for the 3D plate assembly (see Fig. 22) at x = 2p · 300 Hz: VTCR (left) vs. NASTRAN (right).
3188
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191 x 10
2
–5
1.8
Reference solution Frequency band analysis
1.6
Effective displacement
1.4 1.2 1 0.8 0.6 0.4 0.2 0 50
60
70
80
90
100
110
120
130
140
150
Frequency (Hz)
Fig. 24. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 = 2p · 100 Hz, Dx = 2p · 50 Hz.
–5
2
x 10
1.8
Reference solution Frequency band analysis
1.6
Effective displacement
1.4 1.2 1
0.8 0.6 0.4 0.2 0 150
160
170
180
190
200
210
220
230
240
250
Frequency (Hz)
Fig. 25. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 = 2p · 200 Hz, Dx = 2p · 50 Hz.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3189
–6
3
x 10
2.5
Reference solution
Effective displacement
Frequency band analysis 2
1.5
1
0.5
0 250
260
270
280
290
300
310
320
330
340
350
Frequency (Hz)
Fig. 26. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 = 2p · 300 Hz, Dx = 2p · 50 Hz.
–5
2
x 10
1.8 1.6
Effective displacement
1.4 1.2 1 0.8 0.6 0.4 0.2 0 50
100
150
200
250
300
350
Frequency (Hz)
Fig. 27. Effective displacement of Plate 1 of the 3D plate assembly example (see Fig. 22): reference solution (dotted line) and approximate solution (continuous line). x0 2 2p · {100, 200, 300} Hz, Dx = 2p · 50 Hz.
3190
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
x = 2p · 350 Hz. This figure shows that the method works on two different types of structural behaviour: highly dynamic behaviour (frequency range [50; 250 Hz]) and relatively smooth behaviour (frequency range [250; 350 Hz]). 7. Conclusion The approach proposed here, called the ‘‘variational theory of complex rays’’, was originally introduced in order to calculate the vibrations of slightly damped elastic structures in the medium-frequency range. It is a general multiscale approach with a strong mechanical basis. This paper emphasizes the extension of this theory to the analysis across a range of frequencies. The method used to find the solution within a frequency range consists in expanding the matrix and the loading of the system of equations to be solved into a Taylor series. This can be done by separating rapidly-varying quantities and slowly-varying quantities with respect to the frequency. Then, by using an original decomposition of the solution into its mean value over the frequency range and the complementary part, the method can derive an approximate solution for all frequencies in the frequency range. The numerical examples showed the capability of the VTCR to predict the vibrational response of a structure within a frequency range. In this work, the structure is supposed to be an assembly of homogeneous substructures, but this hypothesis can easily be removed [23]. References [1] P.E. Barbone, J.M. Montgommery, O. Michael, I. Harari, Scattering by a hybrid asymptotic/finite element method, Comput. Methods Appl. Mech. Engrg. 164 (1998) 141–156. [2] V.D. Belov, S.A. Ryback, Applicability of the transport equation in the one-dimensional wave propagation problem, Akust. Zh. Sov. Phys. Acout. 21 (1975) 173–180. [3] P. Bouillard, V. Lacroix, E. De Bel, A wave-oriented meshless formulation for acoustical d and vibro-acoustical applications, Wave Motion 39 (4) (2004) 295–305. [4] O. Cessenat, B. Despres, Application of an ultra weak variational formulation of elliptic pdes to the two-dimensional Helmholtz problem, SIAM J. Numer. Anal. 35 (1) (1998) 255–299. [5] E.A. Chadwick, P. Bettess, Modelling of progressive short waves using wave envelopes, Int. J. Numer. Methods Engrg. 14 (1997) 3229–3246. [6] R.D. Ciskowski, C.A. Brebbia, Boundary Element Methods in Acoustics, Elsevier Applied Science, London, 1991. [7] A. Deraemaeker, I. Babuˇska, P. Bouillard, Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions, Int. J. Numer. Methods Engrg. 46 (1999) 471–499. [8] W. Desmet, B. Van Hal, P. Sas, D. Vandepitte, A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems, Adv. Engrg. Software 33 (2002) 527–540. [9] R. Djellouli, C. Farhat, R. Tezaur, A fast method for solving acoustic scattering problems in frequency bands, J. Comput. Phys. 168 (2) (2001) 412–432. [10] J.F. Doyle, Wave Propagation in Structures, Springer Verlag, 1997. [11] C. Farhat, I. Harari, U. Hetmaniuk, A discontinuous Galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the mid-frequency regime, Comput. Methods Appl. Mech. Engrg. 192 (2003) 1389–1419. [12] L.P. Franca, C. Farhat, A.P. Macedo, M. Lesoinne, Residual-free bubbles for the Helmholtz equation, Int. J. Numer. Methods Engrg. 40 (1997) 4003–4009. [13] I. Harari, T.J.R. Hughes, Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains, Comput. Methods Appl. Mech. Engrg. 98 (3) (1992) 411–454. [14] I. Harari, S. Haham, Improved finite element method for elastic waves, Comput. Methods Appl. Mech. Engrg. 166 (1998) 143– 164. [15] I. Harari, P. Barai, P.E. Barbone, Numerical and spectral investigations of Trefftz infinite elements, Int. J. Numer. Methods Engrg. 46 (4) (1999) 553–577. [16] T.J.R. Hughes, Multiscale phenomena: GreenÕs functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Methods Appl. Mech. Engrg. 127 (1–4) (1995) 387–401. [17] M.N. Ichchou, A. Le Bot, L. Je´ze´quel, Energy model of one-dimensional, multipropative systems, J. Sound Vib. 201 (1997) 535– 554.
P. Ladeve`ze, H. Riou / Comput. Methods Appl. Mech. Engrg. 194 (2005) 3167–3191
3191
[18] F. Ihlenburg, I. Babuˇska, Dispersion analysis and error estimation of Galerkin finite element methods for Helmoltz equation, Int. J. Numer. Methods Engrg. 38 (1995) 3745–3774. [19] P. Ladeve´ze, A new computational approach for medium-frequency vibrations (in french), NT Ae´rospatiale YX/SA (1995) 119– 139. [20] P. Ladeve`ze, A new computational approach for structure vibrations in the medium frequency range, C.R. Acad. Sci. Paris 322, Se´rie II b No. 12 (1996) 849–856. [21] P. Ladeve`ze, L. Arnaud, P. Rouch, C. Blanze´, The variational theory of complex rays for the calculation of medium-frequency vibrations, Engrg. Comp. 18 (1999) 193–214. [22] P. Ladeve`ze, L. Arnaud, A new computational method for structural vibrations in the medium frequency range, Comput. Assisted Mech. Engrg. Sci. 7 (2000) 219–226. [23] P. Ladeve`ze, L. Blanc, P. Rouch, C. Blanze´, A multiscale computational method for medium-frequency vibrations of assemblies of heterogeneous plates, Comput. Struct. 81 (12) (2003) 1267–1276. [24] P. Ladeve`ze, P. Rouch, H. Riou, X. Bohineust, Analysis of medium-frequency vibrations in a frequency range, J. Comput. Acoust. 11 (2) (2003) 255–284. [25] R.S. Langley, A wave intensity technique for the analysis of high frequency vibrations, J. Sound Vib. 159 (1992) 483–502. [26] Y. Lase, M.N. Ichchou, L. Je´ze´quel, Energy flow analysis of bars and beams: theoretical formulation, J. Sound Vib. 192 (1996) 281–305. [27] W.K. Liu, Y.F. Zhang, M.R. Ramirez, Multiple scale finite element methods, Int. J. Numer. Methods Engrg. 32 (1991) 969–990. [28] R.H. Lyon, G. Maidanik, Power flow between linearly coupled oscillators, J. Acoust. Soc. Am. 34 (1962) 623–639. [29] M. Malhotra, P.M. Pinsky, Efficient computation of multi-frequency far-field solutions of the Helmholtz equation using pad approximation, J. Comput. Acoust. 8 (1) (2000) 223–240. [30] J.M. Melenk, I. Babusˇka, The partition of unity finite element method: basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289–314. [31] P. Monk, D.Q. Wang, A least-squares method for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 175 (1999) 121–136. [32] J.P.H. Morand, A modal hybridization method for the reduction of dynamic models in the medium frequency range, in: Proc. European Conference on New Advances in Computational Structural Mechanics, 1991, pp. 347–365. [33] E. Perrey-Debain, J. Trevelyan, P. Bettess, Plane wave interpolation in direct collocation boundary element method for radiation and wave scattering: numerical aspects and applications, J. Sound Vib. 261 (2003) 839–858. [34] H. Riou, P. Ladeve´ze, P. Rouch, Extension of the variational theory of complex rays to shells for medium-frequency vibrations, J. Sound Vib. 272 (1–2) (2004) 341–360. [35] P. Rouch, P. Ladeve´ze, The variational theory of complex rays: a predictive tool for medium-frequency vibrations, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3301–3315. [36] C. Soize, Medium frequency linear vibrations of anisotropic elastic structures, Recherche Ae´rospatiale 5 (1982) 65–87. [37] C. Soize, P.M. Hutin, A. Desanti, J.M. David, F. Chabas, Linear dynamic analysis of mechanical systems in the medium frequency range, J. Comput. Struct. 23 (1986) 605–637. [38] C. Soize, Reduced models in the medium frequency range for the general dissipative structural dynamic systems, Eur. J. Mech., A/ Solids 17 (1998) 657–685. [39] T. Strouboulis, K. Copps, I. Babusˇka, The generalized finite element method, Comput. Methods Appl. Mech. Engrg. 190 (32–33) (2001) 4081–4193.